5
$\begingroup$

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).

$\endgroup$
  • 1
    $\begingroup$ It can also happen if you map to a transcriptome sequence instead of the genome sequence because there can be a lot of intron sequence from unprocessed mRNA and from unannotated non-coding RNAs $\endgroup$ – heathobrien Nov 21 '17 at 11:04
5
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ I tried to increase the number of allowed multimapping reads and indeed the number of reads mapped to multiple loci has increased, albeit slightly (6.92% to 7.50%). Most of the non-mapped are categorised into "too short" (in STAR). Do you have any ideas what are these (well, it refers back to my initial question, since this happen in Salmon as well)? $\endgroup$ – kaka01 Nov 22 '17 at 13:11
  • $\begingroup$ Are you trimming adapters from the reads prior to mapping? That might explain the many "too short" reads. Also, have a look at this issue on STAR's Github repo: github.com/alexdobin/STAR/issues/169 $\endgroup$ – Tom Nov 22 '17 at 13:36
  • 1
    $\begingroup$ @kaka01 If accounting for multi-mapping doesn’t solve your problem then there may simply be something wrong with your data: on high quality data sets, mapping total RNA to a genomic reference should typically yield >80% mapped reads. Many “too short” mappings may be indicative of RNA fragments. Did you perform size selection before sequencing? $\endgroup$ – Konrad Rudolph Nov 22 '17 at 16:33
  • $\begingroup$ Tom: No trimming was done prior to mapping. @Konrad The libraries weren't size selected. Your explanation makes sense, aligning the unmapped reads to the collection of rRNAs yielded high mapping rate (which is a sign that those are really rRNAs). Still not sure why STAR classified them as too short, maybe I just have to play around more with the parameters. $\endgroup$ – kaka01 Nov 23 '17 at 9:44
  • 1
    $\begingroup$ @kaka01 “Too short” here simply means that STAR is unable to align the read with high quality: either because the initial read (after trimming) is so short that it could match the reference virtually anywhere (< ~14 bases for something like the human genome) so we have low confidence in the correct origin. Or because STAR, when running with --alignEndsType Local (which is the default) is only able to match a small part of the read (which then leads back to the first problem). $\endgroup$ – Konrad Rudolph Nov 23 '17 at 10:06
2
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

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 lengths1). 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.


1 This works:

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

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.

$\endgroup$
  • $\begingroup$ I 'm mapping to the whole genome indeed. But anyways, assuming that you are mapping against the whole genome, did you get significantly higher mapping rate? $\endgroup$ – kaka01 Nov 22 '17 at 13:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.