Getting much higher coverage with bowtie2 than tophat2

[Studentlost]

Hello,


I ran an analysis on paired end reads through tophat2 using:
tophat2 -p12 -o <tophat_dir> --no-coverage-search <reference genome> R1.fq R2.fq
and the results gave 1.2% coverage.

I ran the same data through bowtie2
bowtie2 -x <index> -1 <R1.fq> -2 <R2.fq> -S <output.sam> and got a 42.75% overall alignment rate.

Why such a big discrepancy? I tried --coverage-search as well and got the same results.

I checked the tophat run.log and it's putting this into bowtie:
bowtie2 -k 20 -D 15 -R 2 -N 0 -L 20 -i S,1,1.25 --gbar 4 --mp 6,2 --np 1 --rdg 5,3 --rfg 5,3 --score-min C,-14,0 -p 12 --sam-no-hd -x

Any idea what's going on?

P.S. I've been getting a silent error in my tophat.log
bam2fastx: /usr/lib64/libz.so.1: no version information available
and
fix_map_ordering: /usr/lib64/libz.so.1: no version information available
I do have libz.so.1.2.3

The process still runs fine and the bam file output can be used for differential analysis... I just have terrible coverage. Any ideas?

[GenoMax]

libz.so.1 is a warning so your observation is unrelated (see #7 in http://seqanswers.com/forums/showthread.php?t=39873).

[Studentlost]

I'm sorry but to which #7 are you referring to? I didn't see any relevance in that thread?

I'm just curious why tophat is missing 40% of the alignment that bowtie2 is finding.

[GenoMax]

Since you had noted the silent errors in your post I was only pointing out that they are warnings and are not related to the difference you are seeing.

Have you tried to run bowtie2 with the same parameters as TopHat? Since bowtie2 in tophat is being run with different parameters it is not surprising that the result there is different.

You can try BBMap as an alternative to tophat.

[Studentlost]

The interesting thing is that I took a look at bowtie 2's defaults and tophat was pretty much on point with running them.

[GenoMax]

But the command options you included for bowtie2 example that you ran directly are not the same as what TopHat used.

The parameters used in TopHat are the defaults for bowtie2 is what you were saying. My apologies.

[GenoMax]

Have you done QC with these reads? Have they been trimmed in parallel (R1 and R2)?

[Studentlost]

Yeah I did QC before and after checking for adapter contamination and trimming using scythe and sickle, respectively.

This is really puzzling me. There's no reason for tophat to give me different results than bowtie2 would...

[colindaven]

Tophat can be quite sensitive to a few parameters especially insert size in my experience. That wouldn't account for the discrepancy here however.

I would try another aligner - and can highly recommend STAR for speed and accuracy.

[Studentlost]

There's no reason Tophat should be failing like this though. Any idea what parameters I can try to change in Tophat to fix the issue?

[sdriscoll]

To cut down on the possible complexity of what is going wrong, try aligning only one of the two mate files with Tophat (i.e. as a single-end alignment) and see if Tophat manages to align more data.

Also..how long are your reads and what are you aligning to?

[Studentlost]

My reads average around 1,000,000 base pairs. I aligned them to the Mmul_1 Rhesus build from Ensembl. I also tried the resMac3 build to compare. I tried with and without a transcriptome index and with and without a reference GTF file. Nothing made a difference. This is blowing my mind.

I can align the paired ends with bowtie2 and I get ~40-60% per sample, but with tophat I get between 0.5% - 4 % per sample.

I was told by another lab working on this that they were able to get the alignment I got with bowtie2 using gsnap. It doesn't make any sense to me why Tophat is the only tool doing this. That means it's unreliable for other alignments in my mind and that bothers me a lot.

[GenoMax]

You are not using original reads? You have reads/contigs that average a megabase each?

TopHat is designed for reads that are a kb or shorter.

[Studentlost]

Wait I'm sorry, I meant my total reads for each strand are at 1 MB. Each read per file is 251.

[sdriscoll]

Yeah, did I read that right? 1,000,000 base paired-end reads? No wonder bowtie2 returns a different result. Are you able to run the alignments with the original RNA-seq reads whatever they were (i.e. PE 100 or whatever)? Then you'll see Tophat actually function.