How to interpret RseQC result for -r/--mate-inner-dist in Tophat2?

Question

[CandiceChuDVM]

Hi all,

I have been playing with options of Tophat2 with my paired-read RNA-seq data. When I was wondering how to use the Tophat2 option -r/--mate-inner-dist to optimize my alignment result, I saw the Q&A from the official website:

"The SAM output of Bowtie2 for paired reads is especially helpful as the 9th field in the SAM alignment lines should show the estimated fragment length, from which you should subtract twice the read length to get the value of the "inner distance" that can be used with the -r parameter"

So I did a Bowtie2 test run for a subset of data and got my test.sam back. However, the number in column 9 is not consistent. Here is a portion of my data after filtered out unmapped reads and unmapped mates:

0 -338 338 310 -310 150 -150 -159 159 -3786 3786 -379 379 -173 173 248 -248 260 -260 0 -128 128 164 -164 -246 246 0

I was wondering how can I deduce the value for -r/--mate-inner-dist from these wide numbers? Could anybody shed some lights on this output?

Thanks!!!

Answer 1

[Devon Ryan]

Have a look at either picard's CollectInsertSizeMetrics or RSeQC. Both can give you a mean and standard deviation that you can then plug into tophat2. I expect that RSeQC is simpler, since it's output can directly be used, as opposed to Picard's, where you'll need to subtract 2*(average read length).

Answer 2

[CandiceChuDVM]

Thanks!

I did inner_distance.py —k 35069833 -i test.sam -o output -r Genome/CanFam3.1.bed and got this back:

 Total read pairs  used 1000000    
Name    Mean               Median     sd    
output  -67.108743152891   -83  55.5160777011065
null device 
      1

Should the result with a negative value ring alarm bells? Or that simply means I should put 67 as the input for -r/--mate-inner-dist ?

[Devon Ryan]

Try -67, if it complains about the negative value then use 0.