Low mapping TopHat

[emolinari]

Hi guys,

another problem to solve. I am running TopHat on some human cells samples I have submitted to our facility for RNAseq. The bioanalyzer values were great (all above RIN 8.5) and thus a library based exclusively on ribo depletion has been created.
I have multiplexed 4 samples in a lane (75bp, single-end) and obtained roughly 50 mil reads for each sample.
This is my TopHat summary:

Reads kept by tophat2 for mapping: 48777778
Reads discarded by tophat2: 321249

Number of unique reads that mapped: 5922988
Percentage of unique reads that mapped: 12.06%

I have used transcriptome for Tophat2...but I cannot understand why such a low mapping percentage (not the same across the samples, ranges across 12-40%).
Do you have any idea???
Other groups in our lab (doing polyA selection, on same cells, using same RNAseq pipeline) map almost 90% of the reads.

Help me please!!!
Thanks
Manu

[kopi-o]

Perhaps a lot of rRNA left in the samples, leading to many multimapping reads due to repetitiveness, or "duplicates" which are really rRNA genes with very high coverage. It could also be adapter contamination but I have seen such low mapping rates often in e g RiboMinus samples due to poor rRNA removal.

[sdriscoll]

It might be useful to grab some random read sequences that didn't align (tophat makes a file called unmapped.bam) and try aligning them with BLAST and see what it comes up with. not very efficient but if you repeatedly see poor / no alignments from BLAST you can be pretty sure that Tophat isn't messing up.

You can pull some reads from it like this:

samtools view unmapped.bam | \
  perl -slane 'if(rand() < 0.000001) { print ">$F[0]\n$F[9]"; }' > random_reads.fa

You can tweak that value in the if test to control how much gets extracted. What I did here should pull on average 1 read per million.

[emolinari]

Thanks for both suggestions!
I have run the script and I BLASTed the unmapped reads of one of my samples...
apparently most of my reads match with a pre-rRNA (or rRNA in general) and also incredibly to some exotic fish, which I assume is a mere coincidence and an artifact.
So my question is: shouldn't rRNA map anyhow to the transcriptome?
More: is there any strategy to limit such a huge contamination? (I really need to do ribo depletion rather than polyA selection)

Thanks a lot!!!
Manu

[Gonza]

Hi,

I have a quick follow up question......I am using that perl script to see why i am getting sooo low % of mapped reads (60%). When i run it i only pull ~ 11 samples of ~ 120 bp (not really a million reads). Am i missing something?

[GenoMax]

Hi,

I have a quick follow up question......I am using that perl script to see why i am getting sooo low % of mapped reads (60%). When i run it i only pull ~ 11 samples of ~ 120 bp (not really a million reads). Am i missing something?

@sdriscoll's code pulls one read from a million sampled so you must have about 11 million reads that are in unmapped.bam. You can change the "0.000001" value to sample more reads.

While it is possible that your samples are contaminated with something (either because of something you did or the sequencing facility did) you will want to take a more careful look before concluding anything.

Brian asked a set of questions in your other post so let us see what you have to say about those.