用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 的全部信息。

具体解释如下

1. table genePred

2. "A gene prediction."

3.    (

4.    string  name;               "Name of gene"

5.    string  chrom;              "Chromosome name"

6.    char[1] strand;             "+ or - for strand"

7.    uint    txStart;            "Transcription start position"

8.    uint    txEnd;              "Transcription end position"

9.    uint    cdsStart;           "Coding region start"

10.    uint    cdsEnd;             "Coding region end"

11.    uint    exonCount;          "Number of exons"

12.    uint[exonCount] exonStarts; "Exon start positions"

13.    uint[exonCount] exonEnds;   "Exon end positions"

14.    )

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

1. #gpd

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

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

1. #gtf

2. Chr1    Araport11   5UTR    3631    3759    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"

3. Chr1    Araport11   exon    3631    3913    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"

4. Chr1    Araport11   start_codon 3760    3762    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"

5. Chr1    Araport11   CDS 3760    3913    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"

6. Chr1    Araport11   exon    3996    4276    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"

7. Chr1    Araport11   CDS 3996    4276    .   +   2   transcript_id "AT1G01010.1"; gene_id "AT1G01010"

8. Chr1    Araport11   exon    4486    4605    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"

9. Chr1    Araport11   CDS 4486    4605    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"

10. Chr1    Araport11   exon    4706    5095    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"

11. Chr1    Araport11   CDS 4706    5095    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"

12. Chr1    Araport11   exon    5174    5326    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"

13. Chr1    Araport11   CDS 5174    5326    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"

14. Chr1    Araport11   CDS 5439    5630    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"

15. Chr1    Araport11   exon    5439    5899    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"

16. Chr1    Araport11   stop_codon  5628    5630    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"

17. 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：

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

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

应用实例：

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

1. 取出所有gene的exon。

2. 对exon进行排序。

3. 对有overlap的exon进行merge。

4. 计算merge后的exon长度。

代码如下：

1. def cmpBed(x,y):

2.    return cmp(x[0],y[0]) or cmp(x[1],y[1])

3. def isOverlap(bed1,bed2):

4.    return not (bed1[2]!=bed2[2] or bed1[0] > bed2[1] or bed1[1]<bed2[0])

5. def mergeBed(beds):

6.    tbed = beds[0]

7.    for bed in beds[1:]:

8.        if isOverlap(bed,tbed):

9.            tbed = (min(bed[0],tbed[0]),max(bed[1],tbed[1]),tbed[2])

10.        else:

11.            yield tbed

12.            tbed = bed

13.    yield tbed

14. exons = []

15. for line in open("test.genepred"):

16.    gene = line.rstrip().split('\t')

17.    exons.extend([(int(s),int(e),gene[1]) for s,e in zip(gene[8].rstrip(',').split(','),gene[9].rstrip(',').split(','))])

18. exons = sorted(exons,cmp=cmpBed)

19. 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。

1. import pandas

2. df = pandas.read_table("/scratch/bcb/ywang52/TData/genomes/hg38/gencode.v22.annotation_GeneSymbol.gpd.gz",header=None)

3. # transcript count

4. df.ix[:,:1].groupby(1).count().plot.bar(legend=False)

5. ylabel('Transcript count')

6. title('Gencode V22')

7. # gene count

8. sdf = df.ix[:,[1,11]].drop_duplicates().groupby(1).count().plot.bar(legend=False)

9. ylabel('Gene count')

10. title('Gencode V22')

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

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

1. 把每个文件读成dataframe，用第一列做行名, 文件名做列名。

2. 按照行名merge dataframe

1. import os

2. import pandas

3. dfs = []

4. for f in os.listdir("."):

5.    if f.endswith(".readcount"):

6.        df = pandas.read_table(f,index_col=0,header=None)

7.        df.columns = [os.path.splitext(f)[0]] # file name as column name

8.        dfs.append(df)

9. 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，后画的部分会覆盖先画的部分。

1. import matplotlib.pyplot as pltcolors = ['red','blue','purple']  # color loops

2. gids = []

3. for cnt,line in enumerate(open("ANXA1.genepred")):

4.    gene = line.rstrip().split("\t")

5.    gids.append(gene[0])

6.    plt.plot([int(gene[3]),int(gene[4])],[cnt,cnt],linewidth=1,color='black') # draw intron

7.    txstart, txstop = int(gene[5]),int(gene[6])

8.    for start,stop in zip([int(s) for s in gene[8].rstrip(',').split(',')],[int(s) for s in gene[9].rstrip(',').split(',')]):

9.        if start > txstop or stop < txstart: # UTR: no overlap with ORF

10.            plt.plot([start,stop],[cnt,cnt],linewidth=3,color='grey')

11.        else: # draw ORF

12.            plt.plot([start,stop],[cnt,cnt],linewidth=3,color='grey')

13.            plt.plot([max(txstart,start),min(txstop,stop)],[cnt,cnt],linewidth=5,color=colors[cnt%len(colors)])                

14. plt.yticks(range(cnt+1),gids)

15. plt.ylim(-0.5,cnt+0.5)

16. plt.tight_layout()

17. 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)

本文编辑：思考问题的熊