Question: Detecting highly variable genes within replicates
1
11 months ago by
Alternative240
Alternative240 wrote:

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

modified 11 months ago by Michele Busby2.1k • written 11 months ago by Alternative240
4
11 months ago by
Michele Busby2.1k
United States
Michele Busby2.1k wrote:

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.

Hello Michele, If I understand your approach correctly, you are calculating the total variance on assumption that the data points( normalized read counts for a certain gene across replicates) follow a normal distribution and then calculating the shot noise on assumption that the same data points follow a Poisson distribution? On another note, are three replicates enough to make assumptions about the total variance for a given gene read counts ? Thank you

I don't think you need to assume a normal distribution but I haven't looked at the theory in a while. The main thing is that you can add uncorrelated variances to get the total variance. They are probably something like a lognormal distribution but the tail is constrained by biology so you can usually assume a normal, in my opinion, and that's what I did for the application we built. Three replicates isn't great but it's what you have.

3
11 months ago by
Friederike5.2k
United States
Friederike5.2k wrote:

You can calculate the variance (or any other measure of variability) fairly easily in R, e.g. with the `rowVars` or `rowWeightedVars`, e.g. :

``````# normalize and rlog-transform the data to get rid of most of the technical/uninteresting variation
rlog_ds <- rlog(DESeq.ds)
# calculate variance (or whatever measure you prefer) per gene
rv <- rowVars(assay(rlog_ds))
# sort, so that the most variable genes will be on top of the object
rv <- order(rv, decreasing = TRUE)
``````

You may want to do subset by condition to find genes that seem to be least reproducible between replicates.

EDIT: updated to include the rlog-transformation following ATPoint's suggestion

1

Adding on this (not sure if you did this), it might be advisable to normalize your data prior to calculating variance. See Calculating significant variance of gene expression across multiple samples for a code example pretty much the same as Friederike showed but with vst-transformed data.

Can a vst transformation be used insted of rlog ?

1

yes, you can (also see the documentation of either `?DESeq2::rlog` or `?DESeq2::vst`

In many cases you will be forced to use `vst` anyway as `rlog` is super slow for many samples. Once you have cohorts of samples, like dozens or hundreds or more, `rlog` is basically unusable as it would take ages. I think that `vst` is currently the recommendation by the `DESeq2` maintainer as a default choice unless there is a specific reason against it. Cannot find the website with the reference for this sentence though.

1
11 months ago by
JC9.4k
Mexico
JC9.4k wrote:

After you compute the dispersion in your samples, you can extract that value with `disp <- getDispersion(y)`, check the plotBCV code to see how is used to create dispersion plots.