用GenePred注释文件进行数据分析
编者注：前几天在生信技能树我们发现了一个神奇的帖子（http://www.biotrainee.com/thread-928-1-1.html ）， 作者用一种并非特别常用的注释文件格式（GenePred table format）解决了多道生信编程直播练习题。小编今天首先简单介绍一下这种格式，随后为大家带来作者的文章。
小编预备知识
GFF/GTF
大多数生物信息学数据的分析和挖掘都十分依赖注释信息，注释文件的好坏对分析结果有着非常重要的影响。
目前，大家常用的有GFF和GTF两种文件。其中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.gpdGeneAnnotable: 其实就是把GTF的所有transcript行的第9列转换变成一个表格。
应用实例：
人类基因组的外显子区域到底有多长？
- 取出所有gene的exon。
- 对exon进行排序。
- 对有overlap的exon进行merge。
- 计算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每条染色体基因,转录本的分布
- 读取genepred格式文件为DataFrame。
- 按照chrom进行group，然后count，最后barplot。
- 按照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无关，但是还是做了吧。
- 把每个文件读成dataframe，用第一列做行名, 文件名做列名。
- 按照行名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基因为例：
- 按行读取genepred文件，第3，4列为转录本的区间，第4，5列为ORF的区间，第9和10列为exon起始和终止位置。
- 先画Intron，即基因全长：第3，4列。
- 再画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
本文编辑：思考问题的熊