对于双端测序RNA-seq数据，TopHat在运行时候，有两个参数-r/--mate-inner-dist和--mate-std-dev分别标识一对reads的间隔长度的期望平均值和标准差，其默认值分别为50bp和20bp。这两个参数本身是个估计值，用于TopHat在map过程中确定一对reads是否匹配到基因组正确位置。如果能够准确设定这两个数值，将会提升TopHat结果的准确性和完整性，参考一个例子。
有两种方法获得这对参数的准确值：
第一种：获取RNA-seq实验建库方法，之后按照以下网址说明计算，RNA-seq差异表达分析工作流程。
第二种：根据RNA-seq数据进行估算，具体步骤为：
使用TopHat默认参数先跑一遍。
使用MISO的pe_utils工具估算。
以下详细介绍pe_utils使用方法。
第一步， 下载对应物种的基因注释文件GTF或者GFF，比如USCS Table Browser（output format选择GTF）或者使用MISO提供的Ensembl版本。如果GTF文件，使用Cufflinks的gffread工具进行转换。
第二步，确定TopHat运行结果的bam文件与基因注释GFF文件，两者基因组命名方法一致。有的使用类似chr1命名，而另外一些使用1。如果不一致，建议修改GFF文件。
# 查看GFF文件中基因组命名
$ awk '{print $1}' hg19USCS.gff | sort -n | uniq -c
# 查看bam文件中基因组命名
samtools view accepted_hits.bam | head -1000 | awk '{print $3}' | sort -n | uniq -c
第三步，筛选较长外显子，比如长度大于1000bp。MISO提供了exon_utils工具用于提取长外显子，但是我们没有能够成功运行过。因此这里提供一个R版本的脚本，比如基因注释文件名为hg19USCS.gff，输出筛选的文件名为hg19USCS_selected.gff。
gffFile <- read.table('hg19USCS.gff', stringsAsFactors = FALSE)
gffExon <- gffFile[gffFile[, 3] == 'exon', ]
exonLen <- abs(gffExon[, 5] - gffExon[, 4])
gffExonSelect <- gffExon[exonLen >= 1000, ]
write.table(gffExonSelect, 'hg19USCS_selected.gff',
row.names = FALSE, col.names = FALSE,
quote = FALSE, sep = '\t')
第四步，使用pe_utils，实例如下：
## 输入bam文件和GFF文件
$ pe_utils --compute-insert-len accepted_hits.bam hg19USCS_selected.gff --output-dir insert-dist/
## 在insert-dist会出现类似accepted_filtered.bam.insert_len文件
## -r/--mate-inner-dist估计值为mean
## --mate-std-dev估计值为sdev
$ head -1 accepted_filtered.bam.insert_len
## mean=165.3,sdev=45.2,dispersion=3.5,num_pairs=5622239
2015年5月23日