用GenePred注释文件进行数据分析

编者注:前几天在生信技能树我们发现了一个神奇的帖子(http://www.biotrainee.com/thread-928-1-1.html ), 作者用一种并非特别常用的注释文件格式(GenePred table format)解决了多道生信编程直播练习题。小编今天首先简单介绍一下这种格式,随后为大家带来作者的文章。

小编预备知识

GFF/GTF

大多数生物信息学数据的分析和挖掘都十分依赖注释信息,注释文件的好坏对分析结果有着非常重要的影响。

目前,大家常用的有GFFGTF两种文件。其中GTF格式是对GFF格式文件的精炼和规范。

GFF文件要求每一行数据必须有由tab键分隔的九个字段,每一个字段代表的含义如下所示。

注:GTF文件前8列和GFF文件相同,第9列信息标签和值用空格分开,不同信息用分号分隔。

GenePred

那什么是GenePred table 格式呢?

这种格式主要用在基因浏览器中基因预测的track。如果有可变剪切的情况,那表格的每一行就是一个 transcript 的全部信息。

具体解释如下

table genePred"A gene prediction."    (    string  name;               "Name of gene"    string  chrom;              "Chromosome name"    char[1] strand;             "+ or - for strand"    uint    txStart;            "Transcription start position"    uint    txEnd;              "Transcription end position"    uint    cdsStart;           "Coding region start"    uint    cdsEnd;             "Coding region end"    uint    exonCount;          "Number of exons"    uint[exonCount] exonStarts; "Exon start positions"    uint[exonCount] exonEnds;   "Exon end positions"    )

如果觉得抽象,我们可以用示例来进行一下对比。小编在这里首先将模式植物拟南芥的gtf文件转化为gpd格式。 head -n 1 看一下gpd文件第一行的样式。

#gpdAT1G01010.1    Chr1    +   3630    5899    3759    5630    6   3630,3995,4485,4705,5173,5438,  3913,4276,4605,5095,5326,5899

再来 grep一下gtf文件中有AT1G01010.1信息的那些行是什么样

#gtfChr1    Araport11   5UTR    3631    3759    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"Chr1    Araport11   exon    3631    3913    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"Chr1    Araport11   start_codon 3760    3762    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"Chr1    Araport11   CDS 3760    3913    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"Chr1    Araport11   exon    3996    4276    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"Chr1    Araport11   CDS 3996    4276    .   +   2   transcript_id "AT1G01010.1"; gene_id "AT1G01010"Chr1    Araport11   exon    4486    4605    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"Chr1    Araport11   CDS 4486    4605    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"Chr1    Araport11   exon    4706    5095    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"Chr1    Araport11   CDS 4706    5095    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"Chr1    Araport11   exon    5174    5326    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"Chr1    Araport11   CDS 5174    5326    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"Chr1    Araport11   CDS 5439    5630    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"Chr1    Araport11   exon    5439    5899    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"Chr1    Araport11   stop_codon  5628    5630    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"Chr1    Araport11   3UTR    5631    5899    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"

我们可以看出,在GTF中,AT1G01010.1 这个 transcript 共有6个CDS,那么对应到相应gpd文件AT1G01010.1 这一行的第8列就是6,而第9列和第10列则是这6个CDS对应的起始和终止位置。

细心的话可能会发现,GTF文件中CDS起始位置在GenePred table中统统少了1,这其实就是两种文件的起始位置从1开始还是从0开始计数的区别。

对GTF和GendPred有了一个初步的了解,我们来看一下作者的原文。

作者正文

写在前面

GTF的产生和流行有其历史的原因。但是从技术角度来讲,这个文件格式是个非常差劲的格式。

GTF格式非常冗余。以人类转录组为例,Gencode V22的GTF文件为1.2G,压缩之后只有40M。大家知道压缩软件的压缩比是和软件的冗余程度。很少有文件能够压缩到1/30的大小。可见GTF格式多么冗余。GTF格式(及其早期版本GFF等)有很好的替代格式。从信息量上来讲:GTF 等价于 GenePred (Bed12) + GeneAnnotable。

GenePred是Jimmy Kent创建UCSC genome browser的时候建立的文件格式。UCSC的文件格式定义是非常smart的,包括之后我可能会讲到的2bit,bigwig格式。

GTF vs GenePred

  • 从文件大小上来看,压缩前:GTF(1.2G) >> Genepred(23M) + GeneAnnotable (2.8M)。压缩后:GTF(40M) >> GenePred(7.8M) +GeneAnnotable (580K)
  • 从可读性来讲,GTF是以gene interval 为单位(行),每行可以是gene,transcript,exon,intron,utr等各种信息,看起来什么都在里面,很全面,其实可读性非常差,而且容易产生各种错误。GenePred格式是以transcript为单位,每行就是一个transcript,简洁直观。
  • 从程序处理的角度来讲:以GTF文件作为输入的程序,如果换成以GenePred格式为输入,编程的难度会降低一个数量级,运算时间会快很多,代码的可读性强很多。

GTF 转换成GenePred:

gtfToGenePred -genePredExt -ignoreGroupsWithoutExons -geneNameAsName2 test.gtf test.gpd

GeneAnnotable: 其实就是把GTF的所有transcript行的第9列转换变成一个表格。

应用实例:

人类基因组的外显子区域到底有多长?

  1. 取出所有gene的exon。
  2. 对exon进行排序。
  3. 对有overlap的exon进行merge。
  4. 计算merge后的exon长度。

代码如下:

def cmpBed(x,y):    return cmp(x[0],y[0]) or cmp(x[1],y[1])def isOverlap(bed1,bed2):    return not (bed1[2]!=bed2[2] or bed1[0] > bed2[1] or bed1[1]<bed2[0])def mergeBed(beds):    tbed = beds[0]    for bed in beds[1:]:        if isOverlap(bed,tbed):            tbed = (min(bed[0],tbed[0]),max(bed[1],tbed[1]),tbed[2])        else:            yield tbed            tbed = bed    yield tbedexons = []for line in open("test.genepred"):    gene = line.rstrip().split('\t')    exons.extend([(int(s),int(e),gene[1]) for s,e in zip(gene[8].rstrip(',').split(','),gene[9].rstrip(',').split(','))])exons = sorted(exons,cmp=cmpBed)print sum([bed[1]-bed[0] for bed in mergeBed(exons)])

hg38每条染色体基因,转录本的分布

  1. 读取genepred格式文件为DataFrame。
  2. 按照chrom进行group,然后count,最后barplot。
  3. 按照gene symbol去重复,然后按照chrom进行group,然后count,最后barplot。
import pandasdf = pandas.read_table("/scratch/bcb/ywang52/TData/genomes/hg38/gencode.v22.annotation_GeneSymbol.gpd.gz",header=None)# transcript countdf.ix[:,:1].groupby(1).count().plot.bar(legend=False)ylabel('Transcript count')title('Gencode V22')# gene countsdf = df.ix[:,[1,11]].drop_duplicates().groupby(1).count().plot.bar(legend=False)ylabel('Gene count')title('Gencode V22')

多个同样的行列式文件合并起来

虽然这个题目与genepred无关,但是还是做了吧。

  1. 把每个文件读成dataframe,用第一列做行名, 文件名做列名。
  2. 按照行名merge dataframe
import osimport pandasdfs = []for f in os.listdir("."):    if f.endswith(".readcount"):        df = pandas.read_table(f,index_col=0,header=None)        df.columns = [os.path.splitext(f)[0]] # file name as column name        dfs.append(df)data = pandas.concat(dfs,axis=1) # merge dataframe by row names

根据GTF画基因的多个转录本结构

以ANXA1基因为例:

  1. 按行读取genepred文件,第3,4列为转录本的区间,第4,5列为ORF的区间,第9和10列为exon起始和终止位置。
  2. 先画Intron,即基因全长:第3,4列。
  3. 再画ORF 和UTR。把第9,10列分割成exon 的starts,stops。遍历每个exon: 如果exon和ORF没有overlap:画成UTR。 反之: 先画成UTR,再把overlap的部分画成ORF。 NOTE: 这里有个trick,后画的部分会覆盖先画的部分。
import matplotlib.pyplot as pltcolors = ['red','blue','purple']  # color loopsgids = []for cnt,line in enumerate(open("ANXA1.genepred")):    gene = line.rstrip().split("\t")    gids.append(gene[0])    plt.plot([int(gene[3]),int(gene[4])],[cnt,cnt],linewidth=1,color='black') # draw intron    txstart, txstop = int(gene[5]),int(gene[6])    for start,stop in zip([int(s) for s in gene[8].rstrip(',').split(',')],[int(s) for s in gene[9].rstrip(',').split(',')]):        if start > txstop or stop < txstart: # UTR: no overlap with ORF            plt.plot([start,stop],[cnt,cnt],linewidth=3,color='grey')        else: # draw ORF            plt.plot([start,stop],[cnt,cnt],linewidth=3,color='grey')            plt.plot([max(txstart,start),min(txstop,stop)],[cnt,cnt],linewidth=5,color=colors[cnt%len(colors)])                plt.yticks(range(cnt+1),gids)plt.ylim(-0.5,cnt+0.5)plt.tight_layout()plt.savefig("test.pdf")

总结

我没有数过以GTF文件作为输入程序解决上述问题究竟有多复杂,代码有多长。但是以GenePred格式为输入的程序,所有代码都不到20行,而且程序清晰简洁,可读性强(当然Python语言的支持也是功不可没),至少新手们有动力看下去,不会被动辄上百行的程序吓跑了。

UCSC 格式汇总:

https://genome.ucsc.edu/FAQ/FAQformat

gtfToGenePred 二进制文件:

http://hgdownload.cse.ucsc.edu/a ... 86_64/gtfToGenePred

本文作者:tsznxx(论坛ID)

本文编辑:思考问题的熊

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-03-01

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏圣杰的专栏

『设计模式』之小试牛刀

为了更好的学习设计模式,以及督促自己完成设计模式的学习,现提笔为记。 怎么的,每周至少也要学一个设计模式!!! 恳请大家的监督和不吝赐教,共同学习和进步!...

2006
来自专栏极客生活

真正好用的python库

这个库是我安装完python环境后第一个安装的库,装上这个库再开始写代码才有底气,作者 Kenneth Reitz 是公认python领域代码写的最好的两个人之...

983
来自专栏牛客网

C++后台腾讯WXG实习面经(已拿offer)

1865
来自专栏吉浦迅科技

DAY25: 阅读硬件的多线程

1124
来自专栏zhangdd.com

zabbix监控交换机出图断断续续问题解决

Zabbix通过snmp监控交换机流量时,之前经常出现出图很不稳定的情况。我先对zabbix以及snmp是怎么来监控流量的做个简单的说明:

2044
来自专栏程序员互动联盟

解密H264、AAC硬件解码的关键扩展数据处理

通过上一篇文章,我们用ffmpeg分离出一个多媒体容器中的音视频数据,但是很可能这些数据是不能被正确解码的。为什么呢?因为在解码这些数据之前,需要对解码器做一些...

3497
来自专栏机器学习和数学

[编程经验] 我是如何半自动抓取素材公社图片的

网络爬虫是一件比较繁琐的事情,特别考验人的耐心。但又是非常令人着迷的一件事,因为当你从网络上爬到了自己的想要的数据,满满的成就感油然而生。但是我对爬虫掌握的并不...

3405
来自专栏LET

CPU Cache简介

真空中光速为299,792,458米/秒,目前,Intel的i7频率可以达到4GHz,简单换算一下,可以得出结论:光(电流)在一个Cycle内移动的距离约为0....

1202
来自专栏编程

C语言怎么学习?以C语言模拟登录为例,学C语言该当如此方可成

好了不用看了,看直了也不会有妹子,看题了: C语言编写代码实现,模拟用户登录------看着是不是很简单?,是不是很轻松?,是不是就有这样的代码? (小编随便百...

3159
来自专栏张善友的专栏

基于DotNetOpenAuth实现OpenID 服务提供者

EverBox网盘(www.everbox.com)是由盛大创新院推出的一款网盘产品,提供了超大的免费存储空间(可升级到 10GB),支持文件同步、文件分享、在...

28910

扫码关注云+社区