-r option - Does it make a difference to the alignment?

[arvid]

Can someone please address the issue I brought up here? How can the -r option that gives the best alignment results be way off? from the means I had estimated from the samples? In fact that option that gives the best result was chosen as an extreme in the + direction just to see what effect it would have on the alignment had one taken a very large value.

Help?

Define "best". If you say the more flagged as properly aligned, the better, then setting a really high allowed insert size will give you more properly paired reads - this seems quite obvious. However, you might risk more falsely mapped pairs than you would want to.
In my way of thinking, it would be better to try to estimate the real mean insert size first (e.g. like in post #7), and then adapt the variance to allow for a sensible inclusion of deviating fragments. Then you aren't lying to the tools that include that parameter in their alignment/assembly/quantification algorithms and heuristics.

[vkartha]

Still unsure why I see +50 giving the best result

So in terms of studying the effect of different -r options on aligning a given sample, I should look at the 'properly paired' statistic?

It is weird because I estimated the mean and std dev inner distance using picard tools. For the two groups (say diseased cases vs control) I got an average inner distance of -16 and -33 respectively - which means there is overlap of the two reads being sequenced, and I know Tophat accepts a negative -r value for the same reason.

On running samtools flagstat, I get the highest % for the 'properly paired' stat as +50 at 92% and it is lower for the above two -r options I tried! Why is that??

Thanks for all your help guys. Appreciate it

Can someone please address the issue I brought up here? How can the -r option that gives the best alignment results be way off? from the means I had estimated from the samples? In fact that option that gives the best result was chosen as an extreme in the + direction just to see what effect it would have on the alignment had one taken a very large value.

Help?

[cedance]

I usually estimate the inner distance by mapping reads using BWA and computing the mean inner distance and standard deviation. And for my samples, most of the time the reads overlap. The mean inner distance can be negative and I would say its important to choose a negative inner distance if your fragment size selection itself is low. If your chosen fragment size is 170 bases and your reads are 100 bases, of course your reads will have a lot of overlap and your mean inner distance is negative. So, I prefer to estimate it for each data every time. And it definitely is much better than choosing default parameters. The number of reads that map and the mapping efficiency is very reliable this way.

[rnaseek]

I think the properly paired is a good metric to use and I would use the one with higher "Properly Paired" reads.

Depending on the library prep, rna fragments can be shorter and result in overlapping paired end reads. I am not sure trying negative r value is a good idea.

Anyway with positive "r" values you can vary different variances and choose the one that gives the most number of properly paired reads.

[vkartha]

So in terms of studying the effect of different -r options on aligning a given sample, I should look at the 'properly paired' statistic?

It is weird because I estimated the mean and std dev inner distance using picard tools. For the two groups (say diseased cases vs control) I got an average inner distance of -16 and -33 respectively - which means there is overlap of the two reads being sequenced, and I know Tophat accepts a negative -r value for the same reason.

On running samtools flagstat, I get the highest % for the 'properly paired' stat as +50 at 92% and it is lower for the above two -r options I tried! Why is that??

Thanks for all your help guys. Appreciate it

[swbarnes2]

I tried using samtools flagstat on one of the output bam files and this is what I got :


43603331 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
43603331 + 0 mapped (100.00%:nan%)
43603331 + 0 paired in sequencing
21831004 + 0 read1
21772327 + 0 read2
38524104 + 0 properly paired (88.35%:nan%)
42357862 + 0 with itself and mate mapped
1245469 + 0 singletons (2.86%:nan%)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

Why am I seeing the +0's with every result???

Also, I'm sorry I'm new to this - how exactly do I interpret this?
Thanks a lot.

Just like it says...the +0 means that no reads in your .bam were flagged as QC-failures. Which makes sense, Bowtie wouldn't add that flag to the .bam it makes.

You have no unmapped reads, because bowtie by default doesn't include them, unlike, say, bwa, which does.

So basically, you want to compare the number of mapped reads, depending on the alignment settings.

[vkartha]

I tried using samtools flagstat on one of the output bam files and this is what I got :


43603331 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
43603331 + 0 mapped (100.00%:nan%)
43603331 + 0 paired in sequencing
21831004 + 0 read1
21772327 + 0 read2
38524104 + 0 properly paired (88.35%:nan%)
42357862 + 0 with itself and mate mapped
1245469 + 0 singletons (2.86%:nan%)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

Why am I seeing the +0's with every result???

Also, I'm sorry I'm new to this - how exactly do I interpret this?
Thanks a lot.

[rnaseek]

Try using samtools flagstat command. flagstat will give the simple alignement stats using the bam file. With different "r" values you can see how they vary.

[vkartha]

Hey guys,
I have paired end reads on RNASeq data and I wanted to test run Tophat2 to map these reads for a given sample using 5 different -r values(that is the mean inner distance between paired ends) just to see how it would affect the alignment.

I'm not exactly sure what I should look at once the alignment is over to figure out the so called 'alignment results'.

Tophat2 gives out a log folder along with the accepted_hits.bam, junction.bed files etc. in which I looked at two files namely:

bowtie.left_kept_reads.fixmap.log
bowtie.right_kept_reads.fixmap.log

These give summaries on the alignment from each end such as:

22630038 reads; of these:
22630038 (100.00%) were unpaired; of these:
2447457 (10.82%) aligned 0 times
15049769 (66.50%) aligned exactly 1 time
5132812 (22.68%) aligned >1 times
89.18% overall alignment rate

Now I notice there are also log files namely:

bowtie.right_kept_reads_seg1.fixmap.log
bowtie.right_kept_reads_seg2.fixmap.log
bowtie.right_kept_reads_seg3.fixmap.log
bowtie.right_kept_reads_seg4.fixmap.log

(And the same 4 for the left kept reads)

I don't understand what these files are or how they are different from the file I mentioned earlier (fixmap.log file)

I'm not sure if I'm looking at the correct log file, and if the fixmap.log file is indeed the right one to look at, I noticed that the alignment percentage is the SAME (i.e. 89%) for ALL -r options (I tried -50, -16, -33, 0, 50) using the same std deviation value. I just cannot understand why!!! - Does this mean that the -r option doesn't make a big difference in the alignment???

Hope to hear back from you guys soon. Thanks a lot.