Question: Hello,

We are looking for a method to detect/extract genes that show variability within biological replicates (typically 3 replicates) of the same condition.

Differential expression algorithms (DESeq2, EdgeR) will consider only genes (or group of genes) that show similarity in their expression profiles between biological replicates. In our case, we are interested in genes that are significantly variable within replicates. Is this information available in DESeq2/EdgeR objects? If not, what would be the best approach for that?

Any though or hint will help,

Best and thanks

Answer: This is a good question and you were wise to ask it because I don't think the obvious approach is the correct one.

For RNA Seq, the observed variance that you get from calculating the variance in the replicates will have a mixture of biological + technical variance.

If the libraries are good, the total variance you observe in biological replicates will be as follows:

Total variance=Biological variance + library prep and sequencing variance + Poisson variance (shot noise, counting noise) from sampling the molecules

To isolate the biological variance I would:

1) Normalize the data

2) Calculate the total variance as the variance of the normalized counts across replicates

3) ignore the library prep and sequencing variance by assuming it is zero (which is usually OK, you can also estimate it from technical replicates if you have them)

4) Estimate the poisson variance as the mean of the expression. This is because the variance of a Poisson distribution is the mean of that distribution.

5)Therefore, the biological variance of the gene can be calculated as:

Total variance=Biological variance + [zero for technical]+the observed mean expression [for Poisson]

Biological variance = Total variance - the observed mean expression

6) Then you can sort by the isolated biological variance. If you don't do this you will see a lot more variance in low expression genes, and low expression genes will appear to be more variable*** just because the noise is higher (i.e. 1 vs 2 is more likely to vary as an artifact of sampling noise than 100 vs 200, which is certainly not counting noise).

This blog explains the ideas better:
http://michelebusby.tumblr.com/post/26913184737/thinking-about-designing-rna-seq-experiments-to

The supplement here goes into depth on all the formulas:
https://academic.oup.com/bioinformatics/article/29/5/656/252753

You don't need to read the paper (trust me, I wrote it), just the first section of the supplement. 
Supplement 1. Sources of variance in expression measurements

Here is a toy example in python. You can also do it in Excel (don't at me):

import numpy as np
measurements=[5, 10, 18]
v=np.var(measurements)
print("Observed variance:", v)

m=np.mean(measurements)
print("Mean expression:", m)

poissonVariance=m

biologicalVariance=v-poissonVariance
print("Isolated Biological Variance:", biologicalVariance)

Observed variance: 28.666666666666668

Mean expression: 11.0

Isolated Biological Variance: 17.666666666666668

***Itay Tirosh has some old papers in yeast looking at variability in expression and what it is correlated with but I can't remember all the details.