Hi Mark, thank you for the reply. Would you include %rRNA as factor in the DE modelling even though the rRNA has been excluded from alignment?
Hello everyone,
Hoping to get some opinions from people with more experience in transcriptomics analyses
We have done total RNAseq sequencing on human whole blood samples using the Illumina Stranded Total RNA Prep, Ligation with Ribo-Zero Plus kit to deplete rRNA and globin RNA and prep the library; and analysing using DRAGEN (there is a rRNA filter step). The kit specifies that there should be less than 10% rRNA after depletion.
Out of the 300 samples we got 100 (30%) with > 11% rRNA; of these we have 25 > 20%, Of these, 11 >30% Of these, 3 > 40% Of these, 1 > 70%
DRAGEN runs the mapping and has a rRNA contig, so all rRNA maps there and is excluded from the alignment. So the final mapped reads count is excluding any rRNA reads.
For all of these samples we have generated plenty of mapped reads to still do DGE analysis (on average more than 100 million per sample, lowest is 70 million per sample)
BUT, for example we have a sample that has 500 million mapped reads and a sample that has 70 million. And the worst example is a sample that total reads was 700 million but 78% of those are rRNA. So technically we still have 100 million mapped reads there but are they affected in some way from having such a high rRNA reads profile?
What I'm worried about is - are these samples with higher rRNA still going to produce meaningful results if we do DGE or in your experience do they have an influence on the DGE analysis (even thought they are excluded and not aligned to the genome?). For example, the sample with 70 million reads vs one with 500 million reads will have much fewer "rare" transcripts detected than the 500 million reads one so these will not really be picked up as we will be comparing the samples so the transcripts that are detected in most/all samples? What other potential issues might be there?
Also any suggestions as to why such a variability of rRNA? We have used one sequencing lane and one flow cell, the samples were prepared on 96 well plates with the same reagents, there doesnt seem to be much difference ie 42% of high rRNA samples were on plate 3, 33% on plate 2, 24% on plate 1 and 18% on plate 4 but that also had fewer samples overall. We have also run a high and a low rRNA sample on DRAGEN with a different rRNA contig and got the same result. All samples have good RIN values that were checked before library prep.
Is it sample variation (ie leftover reagent from extraction)?
Would really appreciate any comments or suggestions! Thank you for reading
1 answer
Hi Anna, Due to human/robot error during sample prep, there can be variation in effectiveness of rRNA removal. This can be reduced by running smaller batch sizes. 96 samples at a time could be problematic. These rRNA reads generally aren't assigned to genes, so their main effect is to reduce the fraction of usable reads. This isn't a problem when rRNA proportion is less than 15%, but can be a problem with higher rates. This reduces resolution of expression levels of other genes. Normally small levels of rRNA isn't a major problem, but when it gets high, you may need to exclude those samples for QC reasons, or alternatively include %rRNA as a factor to be corrected in the DE modelling stage.
I would include %rRNA as factor in the DE modelling if there was evidence that it rRNA content was impacting PCA.
Log in to answer this question.
Just proceed with explroatory analysis and see in a PCA (or similar) whether these samples show any evidence for a bias by the rRNA content. I guess they won't if number of usable reads and total assigned genes counts are remotely comparable. There really is not much use in guessing impacts, just look at the data.
Thank you for the suggestion. I did look at the PCA with all samples and the positive (reference RNA sample) and negative controls (NTC). Posiitve controls and NTCs cluster nicely together and very separate from the sample cluster. 3 samples (one with the highest rRNA %) stand out as separate from the large cluster.
You say there is 100 samples with > 11%rRNA, yet only three overall outliers, so I conclude that rRNA has limited or no real impact here (at least in the PCs examined). Consider removing the outliers, or downweighting in the DE analysis with something like
voomWithQualityWeights()orvoomLmFit(y, sample.weights = TRUE).