Hi everyone,

I’m analyzing a RNAseq dataset and I’ve come across something that is giving me some headaches. 

First some background: I have eight cyanobacterial samples from two different strains (four vs four) that show a (rather small) phenotypic difference under the same experimental conditions. My goal is to identify differentially expressed genes that can potentially explain this phenotypic difference. Though I wasn’t expecting to find a lot of DGEs since the strains are very similar.

Now the data; After aligning and counting, I found the total number of counts differs between samples (figure 1A), as it does too the ribosmal RNA depletion efficiency (figure 1B). All normal up to here, the rRNA depletion kits are not optimized for cyanobacteria and it is not the first time I see this discrepancy in efficiency between samples. A couple more observations, the percentage of rRNA correlates with the total number of counts (figure 1C) and after removing the rRNA counts, all samples have similar number of counts (figure 1D).

![figure_1][1]

I removed the rRNA counts and used DESeq2 to perform variance stabilization on the counts and perform PCA to see how the samples were clustering (figure 2). PC1 explains ~ 50% of the variance and, at first glance, it Iooked like it was separating the strains with the exception of S4. An outlier I thought.

![figure_2][2]

But then I looked again at the percentage of rRNA in each sample and realized the samples’ clustering in PC1 could also be explained by the rRNA percentage, I plotted this and the correlation is there (figure 3).

![figure_3][3]

So I’m wondering a few things.

1. Does it make any sense that most of the variation in the data is explained by the percentage of rRNA? 
2. Or is this just an artifact? 
3. If it is not, how can this be explained? 
4. Does it make the data unusable for my original question?

Thank you in advance!

Note: in case it matters, the ribosomal depletion was performed with MICROBExpress kit.

  [1]: /media/images/6fc7f854-0d9b-4b80-91c2-cf436c5d
  [2]: /media/images/f3b8b3e4-ba7a-43dc-82a9-8eab3bdb
  [3]: /media/images/b22c0ed1-3261-4c7c-beab-dcc07bbe
As a practical note: If you worry that rRNA proportion impacts the estimation of other genes, you can (and likely should) include it as a covariate in DESeq2 with a design `~ strain + rrna_proportion`.

As a technical note: it seems like you removed rRNA counts *prior* to normalization -- meaning that the library sizes associated with the counts are possibly not correct. The easiest way to adapt your current workflow would be to run `estimateSizeFactors` on the original counts; then remove ribosomes, and then use `sizeFactors(de.obj)` to pass in the originally-estimated size factors. This could result in a consistent over-estimation of expression of non-ribosomal genes; resulting in the observed PC1.