用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 条评论
登录 后参与评论

相关文章

来自专栏源码之家

图片类杂志域下载思路

983
来自专栏李成熙heyli

如何写一个webpack插件(一)

前言 最近由于用着html-webpack-plugin觉得很不爽,于是乎想自己动手写一个插件。原以为像gulp插件一样半天上手一天写完,但令人郁闷的是完全找不...

3618
来自专栏企鹅号快讯

phpjiami 数种解密方法

Pwnhub公开赛出了个简单的PHP代码审计题目,考点有两个: 如果说仅为了做出题目拿到flag,这个题目太简单,后台也有数十名选手提交了答案和writeup。...

2257
来自专栏葡萄城控件技术团队

【初学者指南】在ASP.NET MVC 5中创建GridView

介绍 在这篇文章中,我们将会学习如何在 ASP.NET MVC 中创建一个 gridview,就像 ASP.NET Web 表单中的 gridview 一样。服...

2059
来自专栏QQ音乐技术团队的专栏

打通Android Gradle编译过程的任督二脉

本文主要是基于自己在工作当中的一些Android Gradle实践经验,对gradle相关知识做的一个简单总结和分享,希望对大家有帮助。 首先会讲Gradle大...

7818
来自专栏吴老师移动开发

Flutter TextStyle参数解析关于学习

671
来自专栏ascii0x03的安全笔记

PySide——Python图形化界面入门教程(三)

PySide——Python图形化界面入门教程(三)          ——使用内建新号和槽               ——Using Built-In S...

2908
来自专栏炉边夜话

VC++ MFC 常用技巧 (一)

 <?xml:namespace prefix = o ns = "urn:schemas-microsoft-com:office:office" />

611
来自专栏Python入门

你还在为Python中文乱码而感到烦恼?今天老司机给你讲讲!

有没有遇到过这样的问题,读取文件被提示“UnicodeDecodeError”、爬取网页得到一堆乱码,其实这些都是编码惹的祸,如果不能真正理解编码的问题所在,就...

1023
来自专栏自动化测试实战

flask第二十二篇——模板【4】过滤器

2036

扫码关注云+社区