Why total RNA-seq usually yields low mapping rate?

Question

Maybe this is a silly question, but I'm really wondering why we usually get low mapping rates if we map total RNA-seq, but not poly(A)-enriched (in particular, for human, mouse, and zebrafish datasets)?

Doesn't the genome fasta file contain ribosomal RNAs (which are abundantly expected in the total RNA-seq libraries) as well?

I don't have the numbers for human right now, but I recall previously having low mapping rates for human data as well. What triggered me at the moment is actually when I mapped zebrafish data. I got around 60% mapping to the genome and 45% to the transcriptome with one dataset and 36% - 35% with another dataset (note that these are from 2 different studies and both are total RNA-seq, mapping was done using STAR and Salmon, which have their own methods).

Answer 1

A likely explanation is that total RNA-Seq contains a high fraction of reads from ribosomal RNAs. Ribosomal RNAs are present in multiple copies across the genome, hence many reads map to multiple genomic locations and get discarded by the aligner. For example, STAR with default parameters considers a read as unmapped if it maps to more than 10 genomic loci (this behavior can be changed with the --outFilterMultimapNmax option). To confirm if this is actually the case, you could check the number of multi-mapping reads in the aligner's log files.

Answer 2

In particular for mouse, this can happen if the ribo-depletion wasn't terribly efficient, since there's no Rn45s sequence in the reference genome. That, combined with numerous copies of tRNAs, 5S rRNAs, etc. causing issues with multimappers (see the answer from @Tom) can heavily decrease alignment rates.

Note that the human reference genome contains a few copies of the 45S (e.g., on GL000220.1 and chr21), so these reads will all multimap (possibly too many times). I don't know about the zebrafish reference genome.

Answer 3

When sequencing RNA from freshly harvested tissues under proper conditions, you should generally expect > 50% mapped reads. In fact, everything < 80% would usually raise concerns.

From your description (in the question and comments) it sounds as if your samples are potentially degraded and therefore saturated with short RNA fragments, either because you’ve got old tissue, you’re looking at environmental samples, or from improper handling.

And since you’re not performing any enrichment, the sequencing data will then be highly saturated with degraded RNA, too. The reason for this is that RNA-seq is a (hopefully somewhat uniform) sampling of your sample aliquot: everything in it is proportionally represented in the RNA-seq data. That’s why enrichment steps are important, to increase the proportion of whatever fraction we’re interested in.

In particular,

1. Total RNA is mostly composed of genes with many copies: mostly rRNA and (much less) tRNA. Hence the importance of ribosome depletion or poly(A) selection.
2. RNA in general is unstable. Depending on the sample collection (non-tissue sample origin? age of the tissue before harvesting?), a large fraction of the RNA may have degraded into minuscule fragments, which are further shortened by the library preparation. Size selection would get rid of these undesired fragments.

In my experience, the second issue is exacerbated on certain sequencers (I’ve seen it on an Illumina HiSeq 1500) because a mixed library of short and long RNA doesn’t seem to cycle well through the whole fragment length. So even if you have a (small but still present) fraction of long RNA fragments, the sequencer may not be able to synthesise their whole length efficiently. The effect is that there are proportionally more short reads than there were short fragments in the sample.

You can verify this by calculating the insert size distribution of your reads after adapter trimming (e.g. using picard, or by simply tallying read lengths¹). Small reads (< ~14 nt) from degradation fragments are essentially unmappable, since they are too short to have a specific identity. Even if we allowed the mapper to align them their coordinates would be essentially random.

---

¹ This works:

awk 'NR % 4 == 2 {c[length($0)]++} END {for (i in c) print i, c[i]}' in.fastq

Answer 4

You can't assume that the rRNA genes are in the genome. There many copies of them and not all are placed on chromosomes. We've had this issue with human and mouse data for total RNA preps. You need to ensure you map against the whole genome not just the chromosomes.