StringTie在注释基因组时的注意事项

Posted 2023-05-09

参考技术A 在利用RNA-seq注释基因组时，有一个问题就是，我将不同组织来源的转录组数据和参考基因组比对之后，那下一步是1）先将这三个比对结果进行合并，然后用StringTie进行预测，还是2）用StringTie分别进行预测，然后用StringTie的merge模式进行合并？

这个问题的提出，是我采取第二种方式时，发现合并后的基因数减少，觉得哪里不太对劲，于是用IGV检查了不同分析策略的结果，

结论如下：

StringTie 的merge模式能够合并不同的来源的结果，但在合并的同时会根据FPKM，TPM和转录本长度过滤，最终结果可以认为是在所有样本里面都是有所表达的基因，因此最终的数目会少一些。同时由于某些基因表达量低，单个样本里由于read数少无法覆盖基因，因此最终的预测结构还完整。

因此，需要先将BAM合并后，然后进行预测，如果为了输出结果的可靠性，还可以根据FPKM和TPM做过滤。

Bedtools如何比较两个参考基因组注释版本的基因？

目录

• 问题
• 思路

问题

原问题来自：How to calculate overlapping genes between two genome annotation versions？

其实可分为两个问题：

• 一是我组装了一个新的基因组，做了多个注释版本，如何比较它们的feature？比如gene
• 二是我组装了一个新的参考基因组，并做了注释，想和其他已有的同物种参考基因组比较，如何寻找共有和特有的基因（或其他feature）？

思路

第一个问题是比较好解决的，使用bedtools即可。

bedtools比较gff、bed、bam的方法类似，具体可参考这篇教程：
bedtools求overlap
要比较gene，可先从gff中提取gene后再进行比较。或者比较所有feature后再筛选也行。

# 将所有overlap 区域成对输出
 bedtools intersect -a  A.gene.gff3 -b B.gene.gff3 -wa -wb >gene_wa_wb.out
#只要A中的这段区域与B中区域有交集，就输出，而且overlap几次，就输出几次
 bedtools intersect -a  A.gene.gff3 -b B.gene.gff3 -wa >gene_wa.out
#除了输出A中的overlap区域外，还会输出B中的整个区间
 bedtools intersect -a  A.gene.gff3 -b B.gene.gff3 -wb >gene_wb.out
#统计A中每个区域与B overlap的次数
 bedtools intersect -a  A.gene.gff3 -b B.gene.gff3 -c >gene_overlap.count
#只输出A中没有与B overlap的区域
 bedtools intersect -a  A.gene.gff3 -b B.gene.gff3 -v >gene_nonoverlap.count
 bedtools intersect  -a B.gene.gff3 -b  A.gene.gff3 -v >gene_msu_uniq.count

第二个问题需要用比对软件，如gmap进行比对，建立两个基因组的联系，得到gff文件。再利用bedtools比较。

/gmap/bin/gmap_build -D ./ -d A A.fa
/gmap/bin/gmap -D ./ -t 30 -d A -f gff3_gene ../B.cdna > B.gff3

最后的结果要注意，feature不是一一对应的，有一对多，多对一，unique等情况。