搜索
查看: 2411|回复: 13

建议:请放弃GTF文件格式,改用GenePred格式处理转录组注释!

[复制链接]

3

主题

21

帖子

128

积分

注册会员

Rank: 2

积分
128
发表于 2017-2-22 06:22:39 | 显示全部楼层 |阅读模式
本帖最后由 tsznxx 于 2017-2-22 06:22 编辑

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

GTF格式非常冗余。以人类转录组为例,Gencode V22的GTF文件为1.2G,压缩之后只有40M。大家知道压缩软件的压缩比是和软件的冗余程度。很少有文件能够压缩到1/30的大小。可见GTF格式多么冗余。GTF格式(及其早期版本GFF等)有很好的替代格式。从信息量上来讲:GTF 等价于 GenePred (Bed12) + Gene_Anno_table。GenePred是Jimmy Kent创建UCSC genome browser的时候建立的文件格式。UCSC的文件格式定义是非常smart的,包括之后我可能会讲到的2bit,bigwig格式。

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

GTF 转换成GenePred:
   
[Shell] 纯文本查看 复制代码
gtfToGenePred -genePredExt -ignoreGroupsWithoutExons -geneNameAsName2 test.gtf test.gpd


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

应用实例:
生信编程直播第一题:人类基因组的外显子区域到底有多长?
1. 取出所有gene的exon。
2. 对exon进行排序。
3. 对有overlap的exon进行merge。
4. 计算merge后的exon长度。
代码如下:

[Python] 纯文本查看 复制代码
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 tbed

exons = []
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。
[Python] 纯文本查看 复制代码
import pandas

df = pandas.read_table("/scratch/bcb/ywang52/TData/genomes/hg38/gencode.v22.annotation_GeneSymbol.gpd.gz",header=None)
# transcript count
df.ix[:,:1].groupby(1).count().plot.bar(legend=False)
ylabel('Transcript count')
title('Gencode V22')

# gene count
sdf = 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
[Python] 纯文本查看 复制代码
import os
import pandas
dfs = []
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,后画的部分会覆盖先画的部分。

[Python] 纯文本查看 复制代码
import matplotlib.pyplot as pltcolors = ['red','blue','purple']  # color loops
gids = []
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





本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:R语言ggplot2画图系列——Pathway富集分析气泡图
下一篇:对有临床信息的表达矩阵批量做生存分析
回复

使用道具 举报

1

主题

35

帖子

263

积分

中级会员

Rank: 3Rank: 3

积分
263
发表于 2017-2-22 10:33:11 | 显示全部楼层
这个看着挺好!楼主大神!
人生若只如初见!
回复 支持 反对

使用道具 举报

0

主题

4

帖子

146

积分

注册会员

Rank: 2

积分
146
发表于 2017-2-22 11:58:18 | 显示全部楼层
gtfToGenePred  的压缩率好高,
二进制使用前要给使用权限。
lz大神
回复 支持 反对

使用道具 举报

1

主题

41

帖子

266

积分

中级会员

Rank: 3Rank: 3

积分
266
发表于 2017-2-22 12:13:04 | 显示全部楼层
厉害!!,确实方便好多
回复 支持 反对

使用道具 举报

0

主题

6

帖子

49

积分

新手上路

Rank: 1

积分
49
发表于 2017-2-23 00:35:08 | 显示全部楼层
本帖最后由 snjiang12 于 2017-2-23 00:43 编辑

做第五题,读gtf文件多次R studio非常卡,觉得电脑要爆了,根本做不了,试试楼主的

可是用r如何转换啊?算了,自己去找吧
回复 支持 反对

使用道具 举报

3

主题

21

帖子

128

积分

注册会员

Rank: 2

积分
128
 楼主| 发表于 2017-2-23 02:34:54 | 显示全部楼层
snjiang12 发表于 2017-2-23 00:35
做第五题,读gtf文件多次R studio非常卡,觉得电脑要爆了,根本做不了,试试楼主的

可是用r如何转换啊?算 ...

我下载了gtf文件第一件事就是转换格式。 先转换了再用R往下做,不用在R里面做转换。
回复 支持 反对

使用道具 举报

6

主题

36

帖子

415

积分

中级会员

Rank: 3Rank: 3

积分
415
发表于 2017-2-23 22:46:43 | 显示全部楼层
不会pandas和numpy,不过看着真的强,简化好多。。。
回复 支持 反对

使用道具 举报

1

主题

7

帖子

110

积分

版主

Rank: 7Rank: 7Rank: 7

积分
110
发表于 2017-2-23 23:31:12 | 显示全部楼层
明天试一下
回复 支持 反对

使用道具 举报

10

主题

48

帖子

471

积分

版主

Rank: 7Rank: 7Rank: 7

积分
471
QQ
发表于 2017-3-5 21:43:59 | 显示全部楼层
厉害了,我的哥。这个格式再应对什么问题合适呢?它会保留gtf里所有重要信息吗?
回复 支持 反对

使用道具 举报

3

主题

21

帖子

128

积分

注册会员

Rank: 2

积分
128
 楼主| 发表于 2017-3-6 14:34:29 | 显示全部楼层
旭日早升 发表于 2017-3-5 21:43
厉害了,我的哥。这个格式再应对什么问题合适呢?它会保留gtf里所有重要信息吗? ...

genepred格式保留的是基因的所有结构信息,注释信息只保留了最主要的transcript ID 和gene name。GTF的注释信息都在第9列,只需要把所有的transcript行的第9列拿出来,做成一个表格就好了。

genepred格式是用于UCSC的网页呈现的,这些注释信息不属于呈现的内容,所以没必要放在这个表里面。如果你愿意,可以把这些注释信息统统放在genepred格式的最后。
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2017-10-20 04:03 , Processed in 0.117392 second(s), 30 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.