对于双端测序RNA-seq数据,TopHat在运行时候,有两个参数-r/--mate-inner-dist--mate-std-dev分别标识一对reads的间隔长度的期望平均值和标准差,其默认值分别为50bp和20bp。这两个参数本身是个估计值,用于TopHat在map过程中确定一对reads是否匹配到基因组正确位置。如果能够准确设定这两个数值,将会提升TopHat结果的准确性和完整性,参考一个例子

有两种方法获得这对参数的准确值:

第一种:获取RNA-seq实验建库方法,之后按照以下网址说明计算,RNA-seq差异表达分析工作流程

0. 结论

在使用TopHat2匹配双端测序结果后,建议根据成对reads的map基因组位置唯一、方向正确和距离合适的标准,筛选得到的匹配结果。比如,TopHat2可能生成accepted_hits.bam文件,处理方法如下:

Filter TopHat Outputs
1
2
3
4
5
6
7
8
9
10
11
# 首先查看bam文件头部有多少行
$ samtools view -H accepted_hits.bam | wc -l
86

# 筛选成对且成功map到基因组唯一位置的reads,按照上一条输出结果,调整“NR <= 86”
$ samtools view -h accepted_hits.bam | \
    awk '{if (NR <= 86) print $0}; {if($5 == 50 && ($2 == 163 || $2 == 147 || $2 == 83 || $2 == 99)) print $0}' | \
    samtools view -b - > accepted_filtered.bam

$ samtools view accepted_filtered.bam | wc -l
79143942

1. TopHat输出sam文件的第五列

TopHat文档没有解释其输出bam文件(比如accepted_hits.bam)的第五列的意义。按照Bowtie2输出结果来看,是表示映射质量指标Mapping Quality scores(MAPQ),具体计算参考公式$\eqref{eq:1}$。