搜索
楼主: Jimmy

生信编程直播第一题:人类基因组的外显子区域到底有多...

  [复制链接]

0

主题

6

帖子

67

积分

注册会员

Rank: 2

积分
67
发表于 2017-3-13 13:29:49 | 显示全部楼层
本帖最后由 highwind 于 2017-3-20 16:28 编辑

菜鸟一枚,以前就做过一些qPCR
生物数据的来源和格式特别多,解决问题的路径也多种多样,不过这里是来做作业的,就搬个小板凳吭哧吭哧开始算吧:
Q:数据来源是NCBI FTP里面的CCDS.current.txt(https://ftp.ncbi.nlm.nih.gov/pub/CCDS/current_human/),
目的是根据这个列表计算cds覆盖基因组的长度,大概推算Exome的大小


A:问题拆解1(看看数据)
用Excel先看下数据的样子:
前面7列中ccds_status有好几种状态(我后面计算以Public为准,顺便练习下怎么做数据过滤);
后面3列的数字就是今天计算的重点了,看这些数字位数这么大想必应该是指在染色体上的位置,随便挑一个SAMD11的第一个cds小区域925941-926012来看看:

咦?为啥感觉会错了一个碱基?(这里还请大大们来解答一下),但不管怎么说确定了这个数字的含义(虽然里面有一列特意标记了正义链、反义链,但是还好数字只是按正义链给出的),也知道了这里的数字是开头、结尾的闭合区间,所以之后计算碱基要(尾-头+1)或者计算间隔区的时候(头-尾-1)。
最后一列留下identical,好像连rstrip去换行符都不一定要用了(结果我好像用pandas去了

问题拆解2(想想算法)
因为会出现一个基因对应多个ORF,比如下面这个TNFRSF18:
或者不同基因覆盖同一区域的问题,比如下图APITD1-CORT和CENPS
所以这题的关键是染色体号和cds的头尾位置,感觉更像是一次reads mapping。
1.    如前所述导入数据后过滤一下
2.    按染色体号将每个小cds区域存起来,并给前后带[]符号的加个小标记,因为好像后面还有预估UTR部分
3.    按染色体将这些区域按起始碱基位置排序
4.    然后逐个判断有无覆盖
a)     无覆盖的直接长度累积
b)     有覆盖的将区域融合后在进入4判断
5.    最后返回结果

问题拆解3(算法实现)
大学旁听过C,python前天刚学,看之前答题的大大们也有用python的,而且不少Perl的正则和R的作图功能Python都有,以后我还是先从颜值不错的Python 入手吧。
[Python] 纯文本查看 复制代码
import pandas as pd
import re
import time
starttime = time.clock()

def import_file(path):
    f = pd.read_csv(path,sep = '\t') #用Pandas模块,以tab为间隔导入列表文件
    return f

def filter_DT(DT):#过滤在ccds_status列中值为Public的数据
    DT_P = DT[DT['ccds_status'].isin(['Public'])] #Public要以字符列表的形式出现,所以写成['Public']
    DT_I = DT_P[DT_P['match_type'].isin(['Identical'])]
    DT_U = DT_I.loc[:,['#chromosome','cds_locations']] #选择其中部分列用于下一步计算
    return(DT_U)

def count_DT(DT):
    N = len(DT.axes[0])#还有2种方法是: DT.shape[0] 和 len(DT.index)
    return(N)

def print_count_DT(D1,D2): #告知有多少条数据被采纳了
    NT = count_DT(D1)
    NP = count_DT(D2)
    print('表中共有{0}行数据, 本次统计使用了其中{1}行数据'.format(NT,NP))
    
def assign_DT(DT): 
    array = []
    l= []
    for i in range(0,len(DT.values)):
        cdss = DT.values[i][1].split(', ')  #这里的1对应了'cds_locations'列
        chro = DT.values[i][0]              #这里的0对应了'#chromosome'列
        for j in cdss: 
            k = re.split('\[|\]|-',j)       #调试可用 print(int(k[1]))
            if(len(k) == 4):              #[000001-0000002]这种形式
                l = [chro,1,1,int(k[1]),int(k[2])]
                array.append(l)
            elif(len(k) == 2):           #000001-0000002这种形式
                l = [chro,0,0,int(k[0]),int(k[1])]
                array.append(l)
            elif(k[0]):                       #000003-000004]这种形式
                l = [chro,0,1,int(k[0]),int(k[1])]
[/i][/i][i][i][i][i]                array.append(l)
[/i][/i]            else:                             #[000001-000002这种形式
                l = [chro,1,0,int(k[1]),int(k[2])]
                array.append(l)
    DT_A = pd.DataFrame(array)
    DT_A.columns = ['chr','ES','EE','CS','CE']
    return(DT_A)

def sort_DT(DT,UTR):
    DT['STA'] = DT['CS']-DT['ES']*UTR           #先给Exon前后添加估计是需要的UTR长度
    DT['END'] = DT['CE']+DT['EE']*UTR
    DT.sort_values(['chr','STA'],inplace=True) #这里是有顺序的,“chr”列(染色体号)比“STA”列(也就是cds其实碱基位)排序更优先
    DT.index = range(1,len(DT)+1)                #排序后重新写入索引
    print('预估使用UTR长度为单侧{}对碱基'.format(UTR))
    return(DT)

def cal_DT(DT):
    array = DT.loc[:,['chr','STA','END']].values.tolist()
    flag = 0
    for i in range(0,len(array)-1):
        j=i+1 # j代表i的下一行数据的索引值
        if(array[i][2] >= array[j][1]):      #判断为相同或叠加的cds区间
            if(array[i][0] == array[j][0]):  #判断在同一个染色体上出现
                array[j][2] = array[i][2]    #"吃进"叠加或重复区间,注意应保留j行而不是保留i行,否则多重覆盖就得做多轮循环了
                array[i][1],array[i][2] = 0,0#del array会导致array行列数变化不能完成for循环
                flag += 1
            else:
                continue
        else:
            continue
    print('输入cds区域{}个,合并重复后余下{}个区段'.format(len(array),len(array)-flag))
    sum = 0
    for i in range(0,len(array)):
        sum += array[i][2]-array[i][1]+1
    return(sum)

    
CCDS_file_path ='CCDS.current.txt' #文件路径
UTR = 0

DataRaw = import_file(CCDS_file_path) #导入数据
Data2use = filter_DT(DataRaw)              #按Public和identical过滤数据
print_count_DT(DataRaw,Data2use)      #告知采纳数据的大致统计

DataCDS = assign_DT(Data2use)           #摘取需要的数据并组成数据框
DataCDS_sorted = sort_DT(DataCDS,UTR) #数据框排序
print('人类外显子组覆盖约{}对碱基'.format(cal_DT(DataCDS))) #显示计算结果

endtime=time.clock()
print("本次计算用时{:.3f}秒".format(endtime - starttime))    #合计计算用时

计算结果:
[Bash shell] 纯文本查看 复制代码
表中共有34297行数据, 本次统计使用了其中32521行数据
预估使用UTR长度为单侧50对碱基
输入cds区域313228个,合并重复后余下174368个区段
人类外显子组覆盖约30053332对碱基
本次计算用时7.409秒

问题拆解4(算法优化)
这个。。我真的还不会啊。。。慢慢积累经验吧。

学习小结及疑问:
方法:
混合使用了pandas的数据框(用来完成过滤和排序)和数据列表(用来按行计算),不知道为啥如果用pandas数据框做for循环都快一个钟头了居然算不出来。。。
结果:因为考虑很多过滤和去重的问题,所有结果比老师和之前小伙伴们的要小,但至少没差出一个数量级去,算是按心愿完成了第一份作业了,开心。

当然还要请各位大大们指正啦。

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 2 反对 0

使用道具 举报

7

主题

22

帖子

126

积分

注册会员

Rank: 2

积分
126
发表于 2017-3-18 16:39:28 | 显示全部楼层
这一题真的难,没有编程基础,好多步骤读不懂,对于基础差的人,可能需要每一步的具体函数都要略提一下功能(略提一下功能,我们再自己去进一步查询)这样可能好些,按照视频的R步骤,没有做出来,运行很久,依葫芦画瓢按照 “x2yline”的方法做出来了,很快,超级厉害,还要重复去看几遍。
[AppleScript] 纯文本查看 复制代码
setwd('F:/data/Rstudio') #设置工作路径
t1 <- Sys.time() #标记现在开始工作的时间
directory = 'CCDS.current.txt' #目录
# 读取数据并提取第1列和第10列
data <- read.table(directory, sep='\t',stringsAsFactors=F, header=T)[c(1,10)]

get_gene <-function(data_item){
  if (!data_item[2] =='-'){
    exon_ranges <- data_item[2]
    exon_ranges <- substr(exon_ranges, start=2, stop=nchar(exon_ranges)-1)
  }
}
# 用apply执行该函数
# 输入的数据为仅含原始数据第1列和第10列的dataframe
# 输出的数据为c('111-112, 115-135, 125-138', '254-258',...)

get_exon <- function(gene){
  exon <- unique(strsplit(gene,", ")[[1]])
}
# 输入的数据为c('111-112, 115-135, 125-138, 254-258,...')
# 输出的数据为c('111-112','115-135', '125-138', '254-258', ...)
# 注:strsplit的输出结果为列表


get_length <- function(exon){
  loc <- strsplit(exon,"-")[[1]]
  a <- c(as.numeric(loc[1]), as.numeric(loc[2])-as.numeric(loc[1]), as.numeric(loc[2]))
  a
}
# 输入的数据为lapply(c('111-112','115-135', '125-138', '254-258', ...),fun)
# 输出结果为两坐标值和左右两坐标之差
#if (a==0){
#print(loc)
#}


exon_length = 0
for (i in unique(data[,1])){
  gene_i <- paste(apply(data[which(data[1]==i & data[2] != '-'),], 1, get_gene),collapse=', ')
  exon_i <-  lapply(get_exon(gene_i), get_length)
  mat <- matrix(unlist(exon_i), ncol=3, byrow = T)
  base_loc <- matrix(unique(unlist(apply(mat, 1, function(x) c(x[1]:x[3])))))
  exon_length <- exon_length + dim(base_loc)[1] * dim(base_loc)[2]
}

# paste 函数把i号染色体的所有外显子的坐标合并为一个character对象
# gene_i的格式为'111-112, 115-135, 125-138, 254-258,...'
# exon_i的格式为c('111-112','115-135', '125-138', '254-258', ...)
#mat <- mat[order(mat[,2], decreasing = F),]
#mat <- mat[order(mat[,1], decreasing = F),]
# 使用matrix 是因为vector太长会报错
#R memory management / cannot allocate vector of size n MB


difftime(Sys.time(), t1, units = 'secs') # difftime函数用于计算时间差(这里指本次本次运算耗时长度)

print(paste('all exons length is',exon_length)) 


运行结果:"all exons length is 34889110"
感谢大神“x2yline”,都是复制你的运行的哦
与一群努力的人在一起,还有什么理由不努力
回复 支持 反对

使用道具 举报

2

主题

16

帖子

249

积分

中级会员

Rank: 3Rank: 3

积分
249
发表于 2017-4-5 15:50:07 | 显示全部楼层
本帖最后由 蔡卷子 于 2017-4-5 16:44 编辑

第一次做  居然拖到现在
很久不码了  突然一用发现好多不是不会就是忘记了 很方
请大家酌情看我的注释 这全是我的笔记
首先,我为了看清楚哪一步分别是干啥的,就只留了一个三行的test的文件
#chromosome        nc_accession        gene        gene_id        ccds_id        ccds_status        cds_strand        cds_from        cds_to        cds_locations        match_type
1        NC_000001.8        LINC00115        79854        CCDS1.1        Withdrawn        -        -        -        -        None
2        NC_000001.11        SAMD11        148398        CCDS2.2        Public        +        925941        944152        [925941-926012, 930154-930335, 931038-931088, 935771-935895, 939039-939128, 939274-939459, 941143-941305, 942135-942250, 942409-942487, 942558-943057, 943252-943376, 943697-943807, 943907-944152]        Identical

[Python] 纯文本查看 复制代码
#re 正则表达式,可以检查一个穿是否含有某种子串,并将其替换或提取
import re
#os: This module provides a portable way of using operating system dependent functionality.
#这个模块提供了一种方便的使用操作系统函数的方法。操作文件目录
import os
#ordereddict有序字典
#Python附带一个模块,它包含许多容器数据类型,名字叫作collections。
from collections import OrderedDict
#operator模块是python中内置的操作符函数接口,它定义了一些算术和比较内置操作的函数。
#operator模块是用c实现的,所以执行速度比python代码快。
from operator import itemgetter
#operator.itemgetter函数获取的不是值,而是定义了一个函数,通过该函数作用到对象上才能获取值。
exonLength = 0
overlapExons = OrderedDict()
#将C:\Users\Leona\Desktop\生信菜鸟团\第一讲的包导入到python中
os.chdir('D:\\myself\\study\\生信菜鸟团\\第一讲')

之后
[Python] 纯文本查看 复制代码
#检查前面的目录是否真的有改对
#sys: This module provides access to some variables used or maintained by the interpreter and to functions that interact strongly with the interpreter.
#这个模块可供访问由解释器使用或维护的变量和与解释器进行交互的函数。针对系统环境的交互


import sys

path = "/tmp"

## 查看当前工作目录
retval = os.getcwd()
print ("当前工作目录为 %s" % retval)


结果是这个
当前工作目录为 D:\myself\study\生信菜鸟团\第一讲已经检查过打开的地址是否正确了  就可以接下来的工作了
[Python] 纯文本查看 复制代码
#python可以行读
#with 所求值的对象必须要有一个enter(),as 对应的是exit(),将enter的值赋
#值给exit()
#在windouws平台用的是\r\n两个ASCII字符表示换行,python内部用\n表示换行;
#rt模式下,python在读文件的时候自动会把\r\n转换成\n
#这句话就是将text中的赋值到f中
with open ('test.txt', 'rt')as f:
    #将每一行f中的,逐行读入到line中
    for line in f:
        #因为第一行为’#‘是title行,所以跳过
        #print('evey line', line)
        if line.startswith('#'): 
            continue#‘/t’是一个tab character,删去每行最末尾的tab
        
        #每一行结尾都有一个看不见的tab其实,要把它去掉
        #r.strip()默认删除最后的字符
        #r.strip(x)删除字符串最后的的所有x
        #‘/t’是一个tab character,删去每行最末尾的tab
        #line.rstrip()删除string字符串末尾指定字符后生成的新字符串
        #line.rstrip('\t')为删除line末尾的‘\t’
        #?可是若python行读的话 为什么要去掉tab
        line = line.rstrip()
        print('\t')
        print("line = line.rstrip()", line)


        #对于每一行,按照制表符切割字符串,得到的结果构成一个数组,数组的每个元素代表一行中的一列。
        lst = line.split('\t')
        print('\t')
        print("lst = line.split('\t')", lst)
        
        #有时候倒数第二列为‘-’,则这一段为假基因
        #如果lst的倒数第二个为‘-’,跳过
        if lst[-2] == '-':
            continue
        print("lst[-2] == '-'",lst)   
        #提取cds_locations
        #去‘[’
        #re.sub 其中 re就是 regular expression,为正则表达式;sub为substitute表替换
        #inputstr.replace('a','b')将a换成b也可以替换,但是只能是一对一的
        #re.sub 对于输入的一个字符串,利用正则表达式(的强大的字符串处理功能),去实现(相对复杂的)字符串替换处理,然后返回被替换后的字符串
        #将倒数第二个字符的 换成‘’ 啥都没有
        #将中括号去掉,换成空格
        lst[-2] = re.sub('\[|\]', '', lst[-2])
        print('\t')
        print('let[-2]', lst)
        #cds_locations 的exon都是用‘,’隔开的,所以把逗号都去掉
        #cds_locations为lst[-2],split(',')为将其中元素用逗号隔开
        exons = lst[-2].split(',')
        print('exons',exons)
        #打印exons中的每一个元素
        for exon in exons:
            #int(x)将x转化成int
            start = int(exon.split('-')[0])
            end = int(exon.split('-')[1])
            coordinate = lst[0] + ':' + exon
            #我之前一直搞不懂lst[0]为什么会是1,因为我们始终没有改变原数据,只是再不断的创建赋值
            #加上染色体的位置是为了不会将不同的染色体混为一谈
            print('coordinate: is',coordinate)
            #可以看出exon中有overlap
            #所以重复的去掉
            #如果这个坐标以前用过现在就不要用了
            #将这些坐标存到dict中 这就用到了ordereddict
            #若coodinate没有再ordereddict中没出现过就可以继续往下走
            if coordinate not in overlapExons.keys():
                
                overlapExons[coordinate] = 1
                exonLength += end - start
                
print(exonLength)

结果是


现在咱们就要把源文件CCDS导入 计算了

[Python] 纯文本查看 复制代码
import re
import sys
import os
from collections import OrderedDict
from operator import itemgetter

exonsLength = 0
overlapExons = OrderedDict()
#打开文件所在位置
os.chdir('D:\\myself\\study\\生信菜鸟团\\第一讲')
#正文
with open('CCDS current.txt','rt')as f:
    for line in f:
        
        if line.startswith('#'):
            continue
            
        line = line.rstrip()
        lst = line.split('\t')
        
        if lst[-2] == '-':
            continue
        lst[-2] = re.sub('\[|\]', '', lst[-2])
        
        exons = lst[-2].split(',')
        
        for exon in exons:
            start = int(exon.split('-')[0])
            end = int(exon.split('-')[1])
            coordinate = lst[0] + ':' + exon
            if coordinate not in overlapExons.keys():
                overlapExons[coordinate] = 1
                exonsLength += end - start
                
print(exonsLength)

结果
36055234


虽然程序码出来了  但是我还是有点不懂为什么
[Python] 纯文本查看 复制代码
if coordinate not in overlapExons.keys():
overlapExons[coordinate] = 1


中overlapExons不用定义就可以出现 还可以直接找是否not in 在其中
overlapExons[coordiante] = 1 是何解
虽然我知道这一步的目的是为了防止之前计算过的exon被再计算一遍
但是这两步我真心不懂为什么要这么写

谢谢各位大大了




本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 反对

使用道具 举报

0

主题

3

帖子

61

积分

注册会员

Rank: 2

积分
61
发表于 2017-5-5 23:13:44 | 显示全部楼层
351-R+python:
#通过集合进行去重

import re
import os

os.chdir(r'C:/Users/Administrator/Desktop/python_test/')

exonlength = 0
list = []

with open('CCDS.current.txt', 'rt') as f:
    for line in f:
        if line.startswith("#"):
            continue
        lists = line.rstrip().split('\t')
        if lists[-2] == '-':
            continue
        lists[-2] = re.sub('\[|\]', '', lists[-2] )
        exons = lists[-2].split(', ')
        for exon in exons:
            list.append(exon)
    st = set(list)
    for i in st:
        exon_one=i.split('-')
        exon_start = exon_one[0]
        exon_end = exon_one[1]
        exonlength += int(exon_end) - int(exon_start)

                 
print "exonlength= ", exonlength
回复 支持 反对

使用道具 举报

0

主题

16

帖子

149

积分

注册会员

Rank: 2

积分
149
发表于 2017-5-7 14:47:14 | 显示全部楼层
本帖最后由 xxnku 于 2017-8-8 10:17 编辑

[Python] 纯文本查看 复制代码
import sys
import re

exonlength = 0
overlapExons = {}
args = sys.argv

with open(args[1]) as f:
    for line in f:
        line = line.rstrip()
        if line not startswith('#'):

            line = line.split('\t')

        if line[-2] != '-':
            #exons = re.sub('\[|\]','',line[-2])
            exons = line[-2].replace('\[|\]','')
            
        for exon in exons:
            start = int(exon.split('-')[0])
            end = int(exon.split('-')[0])
            
            if exon not in overlapExons.keys():
                overlapExons[exon] = ""
                exonlength += end - start

print (exonlength)

或者
[Python] 纯文本查看 复制代码
import re
import sys

exonlength = 0
overlapExons = {}
args = sys.argv

with open(args[1]) as f:
    for line in f:
        if not line.startswith("#"):
            exons = re.findall('\d+-\d+',line)

            for exon in exons:
                start = int(exon.split('-')[0])
                end = int(exon.split('-')[1])

                if exon not in overlapExons:
                    overlapExons[exon] = ""
                    exonlength += end - start

print (exonlength)



运行结果为36027064
回复 支持 反对

使用道具 举报

16

主题

61

帖子

519

积分

高级会员

Rank: 4

积分
519
发表于 2017-6-13 14:35:53 | 显示全部楼层
btrainee 发表于 2017-1-4 20:32
048-perl+R+python
一开始写的版本有点问题,把所有重复的坐标全去除了,没有考虑到不同的染色体。重新写 ...

这个perl程序计算多久能算出来呢?
回复 支持 反对

使用道具 举报

1

主题

10

帖子

80

积分

注册会员

Rank: 2

积分
80
发表于 2018-5-19 19:47:14 | 显示全部楼层
本帖最后由 A_rrow 于 2018-5-19 19:52 编辑

觉得思路是没错的  不知道为什么算出来的是 32思路:
1.先提取cds_location这一段的序列
2.将该序列的起止坐标提取出来
3.去除存在并列关系的重复外显子
4.将去除后每个外显子的终止坐标减去起始坐标累加即为外显子的总长度


[Python] 纯文本查看 复制代码
import re
from collections import OrderedDict
exons_length = 0
exons_dict = OrderedDict()
with open(r'new.txt','r') as f:
    for line in f:
        line = line.rstrip()
        # 去掉第一行(第一行以 # 开头)
        if line.startswith('#'):
            continue
            '''
            接下来提取 cds_location 字符段
            观察序列可知,有些基因的 cds_location 字符段并没有数值
            '''
        # 每行的数据都是以 \t 分隔
        data = line.split('\t')
        # 去除 某些cds_location 无字符段的基因
        if data[-2] == '-':
            continue
        # 将 [] 换成空格
        data[-2] = re.sub('\[|\]',' ',data[-2])
        # 去掉 原来[] 内以 ,分隔符
        exons = data[-2].split(', ')
        for exon in exons:
            start = int(exon.split('-')[0])
            end = int(exon.split('-')[1])
            # 可能会有两个或同一个基因放在同一个 exons 中
            # 将exon存到字典,重读的的话就不再利用
            if exon not in exons_dict.keys():
                exons_dict[exon] = ''
                exons_length += end - start
    print(exons_length)                  


回复 支持 反对

使用道具 举报

0

主题

1

帖子

63

积分

注册会员

Rank: 2

积分
63
发表于 2019-7-5 12:09:29 | 显示全部楼层
Ryu 发表于 2017-1-8 00:03
028-Python-Ryu

1. 什么是外显子?什么是cds?

你好我想问下 关于第一个python脚本中 exon_dict是什么时候赋了值,最开始不是一个空字典吗
回复 支持 反对

使用道具 举报

0

主题

2

帖子

155

积分

注册会员

Rank: 2

积分
155
QQ
发表于 2019-7-19 17:13:38 | 显示全部楼层
awk -F"\t" '$10~/[0-9]/{OFS="\t";sub("\\]","");sub("\\[","");len=split($10,a,",");for(i=1;i<=len;i++){split(a[i],b,"-");print $3,$1,b[1]"\t"b[2]}}' CCDS.current.txt|awk '{OFS="\t";for(i=$3;i<$4;i++){print $2"_"i}}'|sort|uniq|wc -l
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-11-13 03:54 , Processed in 0.035296 second(s), 23 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.