刘小泽写于19.2.19
之前得到转录本的时候,我们使用TrinityStats.pl Trinity.fasta
这个命令对转录本进行了统计,结果包含了Contig Nx长度的统计,例如:
$TRINITY_HOME/util/TrinityStats.pl Trinity.fasta
################################
## Counts of transcripts, etc.
################################
Total trinity 'genes': 1388798
Total trinity transcripts: 1554055
Percent GC: 44.52
########################################
Stats based on ALL transcript contigs:
########################################
Contig N10: 5264
Contig N20: 3136
Contig N30: 1803
Contig N40: 989
Contig N50: 606
Median contig length: 288
Average contig: 511.61
Total assembled bases: 795066996
Contig N10: 5264
表示至少10%的转录本序列拥有至少5264个碱基长度。这个结果文件中我们一般关注N50这个指标,也就是说组装结果至少一半的转录本序列的长度。
当然,基于全部contigs的统计方式有一个问题:如果转录本比较长,并且拼接得到的isoform太多,那么N50会变大【因为Nx的计算方式是:先将长度从大到小排列,然后计算总长度,依次累加contig,看哪个contig先达到总长一半,就作为N50】,但这并不能真实反映拼接情况,因此还有一个选择:每个基因选择最长的一条isoform进行统计
#####################################################
## Stats based on ONLY LONGEST ISOFORM per 'GENE':
#####################################################
Contig N10: 3685
Contig N20: 1718
Contig N30: 909
Contig N40: 588
Contig N50: 439
Median contig length: 281
Average contig: 433.39
Total assembled bases: 601896081
可以看到,这种方式比一开始得到的N50要小,正是降低了”一条基因有多条长的isoform“这种影响,因此如果使用N50的话,推荐使用ONLY LONGEST ISOFORM per 'GENE'
这个指标中的
相比N50,还有比它更合适的ExN50 ,多了个Ex,意思就是Expressed,要看表达量
ExN50除了考虑上面N50的计算方法,还限制了必须是最高表达的转录本(代表了前x%的总表达量)
$TRINITY_HOME/util/misc/contig_ExN50_statistic.pl \
transcripts.TMM.EXPR.matrix Trinity.fasta >ExN50.stats
结果会生成这样的一个表格,第一列是Ex表达量(%),第二列是ExN50,第三列是num_transcripts
# 可以对统计结果可视化
$ TRINITY_HOME/util/misc/plot_ExN50_statistic.Rscript ExN50.stats
ExN50
之前的统计结果N50是606 bases,图中就是100这个位置(就是因为N50计算考虑的是全部的contigs),这里其实可以过滤掉一些低表达量的contigs
可以看到总体上N50统计结果是在2.4K 到3.5K,最高值再E90的位置,我们就可以使用E90N50这个统计值(大概是3.5K)作为转录本拼接质量的评价(因为拼接约完整越好,相比于之前的606 base,是不是结果好了许多?)
参考:https://github.com/trinityrnaseq/trinityrnaseq/wiki/Transcriptome-Contig-Nx-and-ExN50-stats
欢迎关注我们的公众号~_~
我们是两个农转生信的小硕,打造生信星球,想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到Bioplanet520@outlook.com
网友评论