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

相关文章

来自专栏Fish

CUDA C最佳实践-CUDA Best Practices(三)

10. 运行配置优化 10.1. 占用 10.1.1. 计算占用 10.2. 同步Kernel执行 10.3. 多上下文 10.4. 隐藏寄存器依赖 10.5....

18410
来自专栏数据结构与算法

11:潜伏者

个人感觉这道题在字符串里面算是一个比较综合比较难的题目了 对于选手的数据查找能力有很强的要求 我在做这道题时运用了桶排的思想 不多说了,看代码吧,全网最简 11...

2736
来自专栏小红豆的数据分析

小蛇学python(12)分析《今生今世》人物关系图谱

《今生今世》是渣男胡兰成所写的一部自传体小说。今天我们就来分析一下在他所写的自传中的人物关系图谱,分析一下胡兰成到底和多少女人有关系。

622
来自专栏生信技能树

初学者怎么样才能迅速学会一个软件呢

首先谷歌找到这个教程:http://nix-bio.blogspot.com/2013/10/installing-blat-and-blast.html

1164
来自专栏CDA数据分析师

不可不知的一点Python陷阱

于易于学习以及快速开发更大更复杂的应用,Python渐渐在计算环境中无处不在。尽管明显的语言清晰度和友好会麻痹软件工程师和系统管理员的警觉性 —— 诱使他们编码...

2028
来自专栏腾讯Bugly的专栏

手把手教你如何分析 iOS 系统栈 crash

先上栈,这个 crash 是我们目前开发产品的 top5 crash ? 第一步 对于死在 ojbc _ msgSend 的函数(不仅仅是 msgSend, o...

3668
来自专栏喔家ArchiSelf

MCU上的代码执行时间

在许多实时应用程序中,二八原则并不生效,CPU 可以花费95%(或更多)的时间在不到5% 的代码上。电动机控制、引擎控制、无线通信以及其他许多对时间敏感的应用程...

632
来自专栏木子昭的博客

Python3好用的原生api

对列表进行反序是一个很常见的操作, 但python反向切片的玩法实在是非常简洁, 让人无法拒绝, 其实对某一数据结构进行"反向"是一个很有意...

571
来自专栏微信公众号:Java团长

即将发布的 JDK 10 有 109 项新特性,你喜欢哪些?

按计划,JDK 10 将于 3 月 20 日正式发布。据前 Oracle 员工 Simon Ritter 的统计,JDK 10 总共包含 109 项新特性。当然...

502
来自专栏Ryan Miao

添加PMD插件扫描潜在的bug

上一节使用checkstyle来规范你的项目主要解决了代码编码规范问题,比如缩进换行等。这次继续代码健康工具类PMD。

653

扫描关注云+社区