搜索
查看: 8072|回复: 187

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

  [复制链接]

598

主题

1057

帖子

3525

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
3525
发表于 2016-8-27 21:20:39 | 显示全部楼层 |阅读模式
首先,我们要知道外显子的坐标文件,参考:http://www.bio-info-trainee.com/1101.html
可以看到根据NCBI的CCDS参考,统计到的外显子全才是34729283bp,也就是约35M,根很多文献里面说的都一样!同理,我统计了一下外显子的侧翼长度
54692160 外显子加上前后50bp
73066288  外显子加上前后100bp
90362533  外显子加上前后150bp

而hg19版本基因组里面有着entrez gene ID号的基因是23056个基因,所以我接下来探究一下这些基因的信息!
我们首先看看基因与基因之间的交叉情况
其中有12454266bp的位点,是多个外显子共有的,可能是一个基因的多个外显子,或者是不同基因的外显子
我们主要看看不同基因的交叉情况
不同基因交叉情况共516种,共涉及498种基因!
ZNF341   NFS1
ZNF397   ZSCAN30
我们来看看,如果最不擅长的R来做,会是怎么样的呢?
[Python] 纯文本查看 复制代码
rm(list=ls())

a=read.table(choose.files(),sep = '\t',stringsAsFactors = F,header = T)# 选择你下的CCDs文件

tmp <- apply(a[1:100,], 1, function(gene){# 取前100行数据分析调试
  # gene=a[3,]
  x=gene[10]# Column10 外显子坐标位置列
  if(grepl('\\]',x)){#判断x中是否存在有]这样的符号,如果有就利用正则替换掉。
    x=sub('\\[','',x)
    x=sub('\\]','',x)
    # 这个时候得到的对象还是像这样的“880073-880179, 880436-880525……”
    tmp <- strsplit(as.character(x),',')[[1]]# 我们先从逗号开始分割成小块
    start <- as.numeric(unlist(lapply(tmp,function(y){# 取开始位点
      strsplit(as.character(y),'-')[[1]][1]
    })))
    end <- as.numeric(unlist(lapply(tmp,function(y){#取结束位点
      strsplit(as.character(y),'-')[[1]][2]
    })))
    gene_d <- data.frame(gene=gene[3],#将基因名,染色体,开始、结束位点绑定为数据框
                      chr=gene[1],
                      start=start,
                      end=end
    )
    return (gene_d)#返回数据框
  }
}) 

tmp_pos=c() #构造一个空的向量
lapply(tmp[1:10], function(x){ #取前10个list文件计算调试
 # print(x)
  if ( !is.null(x)){
    apply(x, 1,function(y){
      #print(y)
      for ( i in as.numeric(y[3]):as.numeric(y[4]) )# y[3]为坐标起点,y[4]为终止坐标,历编
        tmp_pos <<- c(tmp_pos,paste0(y[2],"-",i))
    })
    
  }
})
 length(tmp_pos) # 计算exon的长度
 length(unique(tmp_pos)) #计算去重后的exon的长度


用perl我只需要两个单行命令,三五分钟就搞定了,可是在R里面,写这个代码花了我足足半个小时,因为我的R学的很差,而且,最后压根没有做出来!!
但是这个代码很有借鉴意义,你们好好听视频吧!
里面所有的R的apply系列代码,请务必拆开来看!!!

语言,本质上都是一样的,用R语言也可以编程制作一个QQ聊天软件。
实现一个需求,在R里面各种函数都可以,也可以用各种封装好的包,这就是我对R语言的实验,用包!







上一篇:人类hg19这个参考基因组里面各个染色体长度是多少?
下一篇:人类参考基因组中所有基因对应的转录本分布情况如何呢?
回复

使用道具 举报

0

主题

7

帖子

105

积分

注册会员

Rank: 2

积分
105
发表于 2017-1-8 00:03:06 | 显示全部楼层
028-Python-Ryu

1. 什么是外显子?什么是cds?
搜索后在biostars看到
CDS: "A contiguous sequence which begins with, and includes, a start codon and ends with, and includes, a stop codon."

http://www.sequenceontology.org/ ... cvs/term/SO:0000316

Exon: "A region of the transcript sequence within a gene which is not removed from the primary RNA transcript by RNA splicing."

http://www.sequenceontology.org/ ... cvs/term/SO:0000147
那么个人观点是,外显子是形成mRNA后剪接剩下的部分,包括UTR区与CDS,而CDS则是真正编码蛋白的Coding Sequence;
2.此次作业虽然是统计人类基因组的外显子区域,但是还是统计CDS部分,不统计外显子;
3.去重不去除在染色体上overlap的部分,只去除完全一样的部分;

1.首先,根据老师的做法是,从NCBI上下载CCDS的最新版本文件,ftp://ftp.ncbi.nlm.nih.gov/pub/C ... an/CCDS.current.txt
代码如下
[Python] 纯文本查看 复制代码
import csv

ncbi_file="CCDS.20160908.txt"
with open(ncbi_file,'r')as f1:
    file=csv.reader(f1,delimiter="\t")
    next(file)
    sum=0
    exon_dict={} # dict is faster than list
    for record in file:
        if record[9] != "-":
            chr=record[0]
            exon_list=record[9].lstrip("[").rstrip("]").split(", ")
            # print(exon_list)
            for range in exon_list:
                exon=chr+":"+range
                # print(exon)
                if exon not in exon_dict:
                    exon_dict[exon]=""
                    exon_start=int(range.split("-")[0])
                    exon_end=int(range.split("-")[1])
                    # print(exon_start)
                    sum+=exon_end-exon_start
print(sum)
运行结果
36M,与老师的相符,虽然结果一致,但是如果仅仅是重复老师做的,那就没有意思了;其实在看老师的视频前,这道题我想到的方法是去UCSC或者Ensembl下载所有基因的区域再进行统计,老师从NCBI下载算是教会我NCBI还有这样的公共库,下面我就从自己的方法去验证结果:

2.从UCSC下载所有CDS的region来验证,下载方法是使用ucsc的table browser,链接在这里
下载后,文件打开,是这样的,统计“Exon_Starts”和“Exon_Ends”即为每个基因的cds起始和结束位置
代码如下
[Python] 纯文本查看 复制代码
import csv

ucsc_file="ucsc_cds.txt"
with open(ucsc_file,'r')as f1:
    exon_dict={}
    sum=0
    file=csv.reader(f1,delimiter="\t")
    next(file)
    for record in file:
        chr=record[1]
        exon_start_list=record[5].rstrip(",").split(",")
        exon_end_list=record[6].rstrip(",").split(",")
        pos=0
        while pos < len(exon_start_list):
            exon=chr+":"+exon_start_list[pos]+"-"+exon_end_list[pos]
            if exon not in exon_dict:
                exon_dict[exon]=""
                sum+=(int(exon_end_list[pos])-int(exon_start_list[pos]))
            pos+=1
print(sum)
运行结果
33M,基本相符

3.从Ensembl下载所有的CDS进行验证,使用的是Ensembl的Biomart,链接在这里
Biomart需要点选的数据库为
从上图可以看出,如果统计“Exon Chr Start”“Exon Chr End”会把一些非编码的区域也统计进去,需要统计的是“Genomic Coding Start”和“Genomic Coding End”,当然,取巧的统计“CDS Length”这一列应该也可以,但是要去重,我还没有尝试,我试着比较统计Exon和CDS的两个结果
代码如下
[Python] 纯文本查看 复制代码
import csv

ensembl_file="mart_export.txt"
with open(ensembl_file,'r')as f1:
    file=csv.reader(f1,delimiter="\t")
    next(file)

    cds_sum=0
    exon_sum=0

    cds_dict={} # dict is faster than list
    exon_dict={}

    for record in file:
        chr=record[3]
        exon_start=record[5]
        exon_end=record[6]
        exon=chr+":"+exon_start+"-"+exon_end
        if exon not in exon_dict:
            exon_dict[exon]=""
            exon_sum+=(int(exon_end)-int(exon_start))

        if record[7] != "":
            cds_start=record[7]
            cds_end=record[8]
            cds=chr+":"+cds_start+"-"+cds_end
            if cds not in cds_dict:
                cds_dict[cds]=""
                cds_sum+=(int(cds_end)-int(cds_start))

print(cds_sum)
print(exon_sum)
运行结果,先输出CDS总长度再输出Exon总长度
可以看到,CDS结果是39M,基本符合预期。但是外显子统计结果则有117M,相差的有点大,估计是数据有些冗余,有基因的不同剪接和非编码的区域

欢迎各位同学指点,谢谢!


本帖子中包含更多资源

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

x
回复 支持 14 反对 0

使用道具 举报

0

主题

1

帖子

40

积分

新手上路

Rank: 1

积分
40
发表于 2017-1-17 21:31:12 | 显示全部楼层
105-R+Python
刚刚开始接触python ,此题以学习并理解老师的代码为主:
1.生物基本概念:
1)基因DNA分为编码区和非编码区,编码区包含内含子和外显子,一般非编码区具有基因表达的调控功能,如启动子在非编码区。编码区则转录为mRNA并最终翻译成蛋白质。
2)内含子和外显子都被转录到mRNA的前提hnRNA中,当hnRNA进行剪接变为成熟的mRNA时,内含子被切除,外显子被保留。实际上真正编码蛋白质的是外显子,而内含子则无编码功能。
3)CDS Sequencecoding for aminoacids in protein 蛋白质编码区 外显子并不一定编码蛋白质
2.Python
以研究学习别人的代码为主
1)Python模块
在Python中,一个.py文件就被称之为一个模块。
使用模块的最大好处提高了代码的可维护性,以及代码不必从零开始
模块的引用:模块名.函数名
2)代码中的3个模块
A.OS模块
OS 操作系统   与操作系统相关的功能 ,且不受平台限制
常用功能:
        • 显示当前使用平台 os.name()
        • 显示当前python 脚本工作路径 os.getcwd()
        • 删除一个文件  os.remove('filename')  
        • 重命名文件 os.rename("oldname","newname") 
B.CSV模块
csv模块是python的内置模块,直接import csv 即可
csv.reader()---读取csv文件数据
csv.writer()---写入csv文件数据
示例:1.improt csv  #读取csv
                 reader=csv.reader(open('test.csv','rb'))
                 for item in reader:
                        print line
               2.import csv  #写入csv
                  writer=csv.writer(open('test.csv','wb'))---此处必须加上参数'b',不然会出现各行现象
C.RE模块
Regular Expression 正则表达式
提供各种正则表达式的匹配操作,在文本解析、复杂字符串分析和信息提取时是一个非常有用的工具   
3)next()函数
 next()方法当一个文件被用作迭代器,典型例子是在一个循环中被使用,next()方法被反复调用。此方法返回下一个输入行,或引发StopIteration异常EOF时被命中。
for index in range(5):
  line = fo.next()
  print "Line No %d - %s" % (index, line)
4)with open(…) as myfile
对open的调用,会返回一个简单的文件对象,复制给myfile。可以使用一般的文件工具来使用myfile。此对象也支持with语句所使用的环境管理协议。在这个with语句执行后,环境管理机制保证由myfile所引起的文件对象会自动关闭,即使处理该文件时,for循环引发了异常也是如此。
                             


回复 支持 3 反对 0

使用道具 举报

1

主题

41

帖子

250

积分

中级会员

Rank: 3Rank: 3

积分
250
发表于 2017-1-5 22:02:17 | 显示全部楼层
本帖最后由 learnyoung 于 2017-1-8 00:49 编辑

059 python+R-冷洋
试了一些方法,一开始用存入字典的方法来去重,发现效率有点低,我的笔记本CORE_M的处理器跑了好久才出来结果(大概是去细胞房处理了一批细胞的时间)。所以想能不能提高效率一点,每一次通过判断exon在不在字典里,不在还要写进去来去重,太慢了(对于 我这个老爷本来说是这样,等的好焦急)。在python中set就是一个没有重复的集合,我把用正则 [0-9]+-[0-9]+匹配到外显子存入一个list中,然后将list 经过set化就去重了,验证了一下速度非常快,比写入字典的方法要快很多解题思路:
1,正则[0-9]+-[0-9]+直接逐行匹配CDS_location
2,将匹配到的CDS_location存入一个list中
3,set(list)去重复
4,首尾相减,算长度
方法一:
[Python] 纯文本查看 复制代码
#coding=utf-8
import re
f=open('CCDS.current.txt','r')
sum=0
list=[]
dit={}
for line in f:
    if line.startswith('#'):
        continue
    line=line.rstrip()
    exons=re.findall('[0-9]+-[0-9]+',line)#匹配像00000-00000这种格式的数据,检查了一下可以将所有的外显子匹配到
    for exon in exons:
        list.append(exon)#将每一个外显子存入list中
st=set(list)#将列表set去重
for i  in st:
    exon_one=i.split('-')#还可以用find函数找到-的位置然后用切片操作找到起始点和终止点
    exon_start=exon_one[0]
    exon_end = exon_one[1]
    sum+=int(exon_end)-int(exon_start)
print sum
。结果为36027064

方法二:
[Python] 纯文本查看 复制代码
#coding=utf-8
import re
f=open('CCDS.current.txt','r')
g=open('NCCDS.txt','w')
sum=0
dit={}
for line in f:
    if line.startswith('#'):
        continue
    line=line.rstrip()
    exons=re.findall('[0-9]+-[0-9]+',line)
    for exon in exons:
        exon_one=exon.split('-')
        exon_start=exon_one[0]
        exon_end=exon_one[1]
        if exon not in dit.keys():
            dit[exon]=1
        sum+=int(exon_end)-int(exon_start)
print sum

[Python] 纯文本查看 复制代码
#coding=utf-8
import re
f=open('CCDS.current.txt','rt')
sum=0
lt=[]
dic={}
for line in f:
    line=line.rstrip()#去掉每一行最后的换行符
    if line.startswith('#'):
        continue
    lst=line.split('\t')#以空格将
    if lst[-2]=='-':
        continue
    exons=re.sub("\[|\]"," ",lst[-2])#去除中括号
    exons=exons.split(',')
    for exon in exons:
        lt.append(exon)
st=set(lt)
for i in st:
    exon_one=i.split('-')
    exon_start = exon_one[0]
    exon_end = exon_one[1]
    sum += int(exon_end) - int(exon_start)
print sum







回复 支持 3 反对 0

使用道具 举报

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

使用道具 举报

1

主题

5

帖子

64

积分

注册会员

Rank: 2

积分
64
发表于 2017-2-8 11:07:44 | 显示全部楼层
006-perl
第一课就看了好几遍,理解了老师的思路,但是在os下的操作尚未掌握,未能成功运行。在学习共享的Linux 和perl语言的基础知识资料。寒假已经补了一些Linux基本命令。希望自己能把基础学习好后完成作业
回复 支持 2 反对 0

使用道具 举报

0

主题

30

帖子

331

积分

中级会员

Rank: 3Rank: 3

积分
331
发表于 2017-1-7 08:59:39 | 显示全部楼层
本帖最后由 x2yline 于 2017-2-18 10:29 编辑

077 R-python- shell  x2yline
1.我的解题思路:




2.python实现
(1)jupyter编辑器太强大了,非常好用,但是没有查看当前变量的功能,所以最终还是选择spyder作为python编写平台(有shift+enter键相当于Rstudiod的ctr+r键,也有查看当前已有变量数值的功能
(2)关于open(file, 'rt')的解释
w,r,wt,rt都是python里面文件操作的模式。w是写模式,r是读模式。
t是windows平台特有的所谓text mode(文本模式),区别在于会自动识别windows平台的换行符。
类Unix平台的换行符是\n,而windows平台用的是\r\n两个ASCII字符来表示换行,python内部采用的是\n来表示换行符。
rt模式下,python在读取文本时会自动把\r\n转换成\n.
wt模式下,Python写文件时会用\r\n来表示换行。

(3)代码实现
[Python] 纯文本查看 复制代码
import pandas as pd
import numpy as np
file = r'E:\r\biotrainee_demo1\CCDS.current.txt'


def calculate_exon(file):
  data = pd.read_csv(file, sep='\t',\
    usecols=[0,9])

#data.loc[1:10,:]
#  data[0:3]
#  data.iloc[1:3]
#  data.iloc[3]
  all_length = 0

  for i in data.iloc[:,0].unique():
    # get the data of chrosome i
    # iloc[row_vector,col_vect]
    # iloc[row_vector]
    data_i = data.loc[data.iloc[:,0] == i]
    type(data_i)
    type(data_i.iloc[:,1])
    # remove the '[]' in column2
    data_j = data_i.iloc[:,1].apply(lambda x: x[1:-1])
    data_p = data_j.apply(lambda x: x.split(', '))
    data_g = data_p.apply(lambda x: pd.Series(x))
    # 把nan填充为 0-0
    data_f = np.array(data_g.fillna('0-0'))
    # 去除重复的外显子
    data_f = np.unique(data_f.reshape((data_f.shape[0]*data_f.shape[1], 1)))
    data_f = pd.DataFrame(data_f)
    data_m = data_f.apply(lambda x: \
      x.apply(lambda y: (y.split('-')[0])))
    data_n = data_f.apply(lambda x: \
      x.apply(lambda y: (y.split('-')[-1])))
    # pd.to_numeric can only apply to a 1-d array
    data_mi = data_m.apply(lambda x: pd.to_numeric(x, downcast='float'))
    data_ni = data_n.apply(lambda x: pd.to_numeric(x, downcast='float'))
    all_length += (data_ni - data_mi).sum().sum()
  return(all_length)

length = calculate_exon(file)
print(length)

运算速度有点慢,因为是临时学的pandas和numpy,很多步骤还没有优化

未去重overlap结果为:36046283

3.R语言实现

第一次是未考虑外显子重叠情况写的R语言代码如下:
运行计算时间为:14.74084 secs
最后运行结果为:36048075


代码如下:
[AppleScript] 纯文本查看 复制代码
setwd('E:\\r\\biotrainee_demo1')
t1 <- Sys.time()
directory = 'CCDS.current.txt'

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)
  }
}

get_exon <- function(gene){
  exon <- unique(strsplit(gene,", ")[[1]])
}

get_length <- function(exon){
  loc <- strsplit(exon,"-")[[1]]
  a <- as.numeric(loc[2])-as.numeric(loc[1])
  if (a==0){
    print(loc)
  }
  a
}


exon_length = 0
exon_length_items = NULL
for (i in unique(data[,1])){
  
  gene_i <- paste(apply(data[which(data[1]==i & data[2] != '-'),], 1, get_gene),collapse=', ')
  exon_i <-  get_exon(gene_i)
  print('************************')
  print(paste('染色体', i, '的无效外显子区域为'))
  exon_i_length <- sapply(exon_i, get_length)
  exon_length <- exon_length + sum(exon_i_length)
  exon_length_items <- c(exon_i_length, exon_length_items)
  names(exon_length_items)[1:length(exon_i_length)] <- i
}

# 耗时长度
difftime(Sys.time(), t1, units = 'secs')

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

经过差不多一下午的思考和调试(大多数时间都在改写代码,因为很多种方法R语言都跑不动

最后终于也做出来了,速度也还可以:
运行计算时间为:18.71507 secs
最后运行结果为: 34889110
运行时间:2017-01-07 18:06:38 CST
系统配置:
R version 3.3.2 (2016-10-31)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1


代码如下:
[AppleScript] 纯文本查看 复制代码
setwd('E:\\r\\biotrainee_demo1')
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){
  # 用apply执行该函数
  # 输入的数据为仅含原始数据第1列和第10列的dataframe
  # 输出的数据为c('111-112, 115-135, 125-138', '254-258',...)
  if (!data_item[2] =='-'){
    exon_ranges <- data_item[2]
    exon_ranges <- substr(exon_ranges, start=2, stop=nchar(exon_ranges)-1)
  }
}

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

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

exon_length = 0
for (i in unique(data[,1])){
  # paste 函数把i号染色体的所有外显子的坐标合并为一个character对象
  # gene_i的格式为'111-112, 115-135, 125-138, 254-258,...'
  gene_i <- paste(apply(data[which(data[1]==i & data[2] != '-'),], 1, get_gene),collapse=', ')
  # exon_i的格式为c('111-112','115-135', '125-138', '254-258', ...)
  exon_i <-  lapply(get_exon(gene_i), get_length)
  mat <- matrix(unlist(exon_i), ncol=3, byrow = T)
  #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
  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]
}

# 耗时长度
difftime(Sys.time(), t1, units = 'secs')

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












本帖子中包含更多资源

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

x
回复 支持 2 反对 0

使用道具 举报

7

主题

34

帖子

192

积分

注册会员

Rank: 2

积分
192
发表于 2017-1-6 12:05:17 | 显示全部楼层
本帖最后由 babytong 于 2017-5-2 16:00 编辑

025 perl +python

[Python] 纯文本查看 复制代码
import re
sum=0
list_all=[]

with open ('D:/sxjns/question1/CCDS.current.txt') as f:
    for line in f:
        if line.startswith('#'):continue
        lists=line.rstrip().split('\t')
        exons=re.findall('[0-9]+-[0-9]+',lists[-2])
        [list_all.append(exon) for exon in exons if not exon in list_all]  
    for i in list_all:
        exon_one=i.split('-')
        sum+=int(exon_one[1])-int(exon_one[0])
    print(sum)


1.with打开文件优于传统方法打开
2.事先判断存在性,再append 就不需要进行去重
3.以上代码运算结果为36027374 以上代码执行相当慢
所以做出如下修改
[Python] 纯文本查看 复制代码
import os
import re
os.chdir("D:")
sum=0
list_all=[]
with open('CCDS.current.txt','r') as f:
    for line in f :
        if line.startswith('#'):continue
        lists=line.rstrip().split('\t')
        exons=re.findall('[0-9]+-[0-9]+',lists[-2])
        for exon in exons:
            list_all.append(exon) 
        
    list_all=set(list_all)
    for i in list_all:
        exon_one=i.split('-')
        sum+=int(exon_one[1])-int(exon_one[0])
    print(sum)

025 perl +python


欢迎关注博客: 生信客部落
回复 支持 2 反对 0

使用道具 举报

1

主题

3

帖子

73

积分

注册会员

Rank: 2

积分
73
发表于 2017-1-7 13:54:16 | 显示全部楼层
我有一个问题,大家是不是没有考虑部分重叠的情况?

本帖子中包含更多资源

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

x
回复 支持 2 反对 0

使用道具 举报

598

主题

1057

帖子

3525

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
3525
 楼主| 发表于 2017-1-17 09:06:22 | 显示全部楼层
朝曦曦 发表于 2017-1-17 05:31
大家是怎么上传代码的?

http://www.biotrainee.com/forum-62-1.html 请自己修改一下
回复 支持 1 反对 0

使用道具 举报

0

主题

1

帖子

16

积分

新手上路

Rank: 1

积分
16
发表于 2017-1-10 09:59:46 | 显示全部楼层
本帖最后由 tony2015116 于 2017-1-16 19:34 编辑

这几周忙着复习考试,所以没做作业,抱歉了视频看完了,表示压力山大,很多地方不明白
回复 支持 1 反对 0

使用道具 举报

0

主题

12

帖子

157

积分

注册会员

Rank: 2

积分
157
发表于 2017-1-3 20:41:41 | 显示全部楼层
用perl写了一个,对于所注释的行做的不够好。

回复 支持 反对

使用道具 举报

0

主题

12

帖子

157

积分

注册会员

Rank: 2

积分
157
发表于 2017-1-3 20:42:40 | 显示全部楼层
本帖最后由 btrainee 于 2017-1-3 20:43 编辑

[Perl] 纯文本查看 复制代码
#!/usr/bin/perl
use strict;
use warnings;

open  IN ,"CCDS.current.txt" or die $!;
my $loc;
my %hash;
while (<IN>){
>---chomp;
>---next if /^#/;
>---$_=~/\[(.*?)\]/;
>---$loc.=$1.",";
}
close IN;
my @location=split /,/,$loc;
foreach my $tt (@location){
#>--next if $tt=="-";
>---my @tmp=split /-/,$tt;
>--->---foreach my $i ($tmp[0]..$tmp[1]){
>--->--->---$hash{$i}=1;
>--->---}
>---}
my $hash_size=keys %hash;
print "$hash_size\n";
回复 支持 反对

使用道具 举报

5

主题

22

帖子

328

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
328
发表于 2017-1-4 10:57:15 | 显示全部楼层
回复

使用道具 举报

1

主题

7

帖子

90

积分

注册会员

Rank: 2

积分
90
发表于 2017-1-4 16:01:10 | 显示全部楼层
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w

use strict;
my @temp;
my @nextdata;
my @nextnextdata;
my $zuobiao;
my $site;
my %hash;
my $sum=1;
open my $data,"CCDS.current.txt" or die $!;

while(<$data>){
chomp;
next if /^#/;
 @temp=split/\t/,$_;
if ($temp[9]=~/\[(.*)\]/){
 @nextdata=split/,\s/,$1;
foreach $zuobiao (@nextdata){
 @nextnextdata=split/-/,$zuobiao;
foreach  $site ($nextnextdata[0]..$nextnextdata[1])
{$hash{$site}=1}
}
}
}
foreach (keys %hash){$sum++;}
print $sum;
close $data;


计算结果31691416
perl-024
回复 支持 反对

使用道具 举报

0

主题

12

帖子

157

积分

注册会员

Rank: 2

积分
157
发表于 2017-1-4 20:32:03 | 显示全部楼层
本帖最后由 btrainee 于 2017-1-4 20:35 编辑

048-perl+R+python
一开始写的版本有点问题,把所有重复的坐标全去除了,没有考虑到不同的染色体。重新写了个版本,使用最新下载的20170101CCDS.current.txt,如下:
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl
use strict;
use warnings;
open  IN ,"CCDS.current.txt" or die $!;
my $loc;
my %hash;
while (<IN>){
>---chomp;
>---next if /^#/;
>---my @F=split;
>---(my $match)=$_=~/\[(.*?)\]/;
>---$match=~s/(\d+)\-(\d+)/$F[0]\-$1\-$2/g;
>---$loc.=$match.",";
}
close IN;
my @location=split /,/,$loc;
foreach my $tt (@location){
#>--next if $tt=="-";
>---my @tmp=split /-/,$tt;
>--->---foreach my $i ($tmp[1]..$tmp[2]){
>--->--->---$hash{"$tmp[0]:$i"}=1;-
>--->---}
>---}
my $hash_size=keys %hash;
print "$hash_size\n";

运行结果为:35288998
然后又发现一个问题,把版主的perl,copy过来运行,发现结果不一致
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl
use strict;
use warnings;
#############格式化输出
my (@F,$match,@L,%hash,$tmp);
open IN,"CCDS.current.txt" or die $!;
open OUT,">","CCDS.20170101.exon.txt" or die $!;
while (<IN>){
>---next if /^#/;
>---@F=split;
>---($match)=$_=~/\[(.*?)\]/;
>---my @tmp=split/,/,$match;
>---foreach (@tmp){
>--->---$_=~s/-/\t/g;
>--->---print OUT "$F[2]\t$F[6]\t$F[0]\t$_\n";
>---}
}
close IN;
close OUT;
#############算外显子总长度,或者加上侧翼50/100/150后的总长度
open FH,"CCDS.20170101.exon.txt" or die $!;
while (<FH>){
>---chomp;
>---@L=split;
>---foreach (($L[3]-$ARGV[0])..($L[4]+$ARGV[0])){
>--->---$hash{"$L[2]:$_"}=1;
>---}
}
close FH;
$tmp++ foreach keys %hash;
print "$tmp\n";

$ARGV[0]参数给为0,运行结果为:34889110
由于python还不太会,所以直接copy了视频的代码,自己运行了一下,神奇了,结果也不同,代码如下:
[Python] 纯文本查看 复制代码
import re
import os
from collections import OrderedDict
from operator import itemgetter

os.chdir('~/bob/btraine/1_exon_length')
exonlength = 0
overlapExons = OrderedDict()
with open ('CCDS.current.txt','rt') as f:
>---for line in f:
>--->---if line.startswith('#'):
>--->--->---continue
>--->---line = line.rstrip() #去掉右边的换行符空格TAB等,lstrip strip
>--->---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
>--->--->--->---exonlength += end-start
print (exonlength)

运算结果:36048075
一时有点晕,用的同一个文件,怎么计算结果都不一样,有点费解。
回复 支持 反对

使用道具 举报

3

主题

8

帖子

74

积分

注册会员

Rank: 2

积分
74
发表于 2017-1-4 20:52:27 | 显示全部楼层
本帖最后由 猪兔子2号 于 2017-1-4 21:25 编辑

摩羯座:007

首先下载文件,linux 命令行下输入 wget -bc ftp://ftp.ncbi.nlm.nih.gov/pub/C ... n/CCDS.20160908.txt

然后写脚本,前面跟群主的一样

[Perl] 纯文本查看 复制代码
#命令行输入perl test.pl CCDS.20160908.txt >test.txt
while(<>){
               chomp;
               next if /^#/;
               @F=split("\t",$_);
                /\[(.*?)\]/;
                @tmp=split(",",$1);
          foreach (@tmp){ 
         @region=split("-",$_);
          print "$F[2]\t$F[6]\t$F[0]\t$region[0]\t$region[1]\n";                                   } 
              }#命令行输入perl test1.pl test.txt >test2.txt
$i=0
while(<>){
             @tmp=split("\t",$_);
             foreach ($tmp[3]..$tmp[4]) {
                  if($hash{"$tmp[2]:$_"}){print $hash{"$tmp[2]:$_"}."\t"$tmp[0]}
                  else{
                            $hash{"$tmp[2]:$_"}=$tmp[2]; $i++;}                                                                                       
         }                         
   }print $i;


这里不一样的地方时,命名一个变量,起始为零,判断该位点是否存在,不存在加一,输出所有外显子位点之和,并且输出那些基因具有重复位点两端的基因有哪些,因为有些基因不知道为什么在这个文件里会重复出现,但是起始和终止位置稍有不同,可能因为编码出来的蛋白不同吧,所以同一基因也会被记录多次,所以存在部分同一位点交叉基因是同一基因的情况。
cat test1.txt | sort |uniq >test2.txt
当然也可以查看没有去除重复是有多少位点cat test.txt |awk -F "\t" '{sum+=$1-$2}END{print sum}'






回复 支持 反对

使用道具 举报

0

主题

4

帖子

68

积分

注册会员

Rank: 2

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

您运行的结果和视频是一致的呀
回复 支持 反对

使用道具 举报

0

主题

12

帖子

157

积分

注册会员

Rank: 2

积分
157
发表于 2017-1-4 22:23:39 | 显示全部楼层
003-python 发表于 2017-1-4 22:03
您运行的结果和视频是一致的呀

是我自己使用的同一个文件,但是用的两个不同的perl程序和那个python程序,3个程序运行出来的结果不同
回复 支持 反对

使用道具 举报

0

主题

4

帖子

68

积分

注册会员

Rank: 2

积分
68
发表于 2017-1-4 22:35:48 | 显示全部楼层
btrainee 发表于 2017-1-4 22:23
是我自己使用的同一个文件,但是用的两个不同的perl程序和那个python程序,3个程序运行出来的结果不同 ...

哦哦。。不好意思,只关注Python了
回复 支持 反对

使用道具 举报

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

本版积分规则

QQ|关于我们|手机版|小黑屋|生信技能树    

GMT+8, 2017-6-24 09:57 , Processed in 0.060280 second(s), 33 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.