Question: Detecting highly variable genes within replicates
1
gravatar for Alternative
8 weeks ago by
Alternative230
Alternative230 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

ADD COMMENTlink modified 8 weeks ago by Michele Busby1.9k • written 8 weeks ago by Alternative230
3
gravatar for Friederike
8 weeks ago by
Friederike3.8k
United States
Friederike3.8k 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

ADD COMMENTlink modified 8 weeks ago • written 8 weeks ago by Friederike3.8k
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.

ADD REPLYlink written 8 weeks ago by ATpoint15k
2
gravatar for Michele Busby
8 weeks ago by
Michele Busby1.9k
United States
Michele Busby1.9k 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.

ADD COMMENTlink modified 8 weeks ago • written 8 weeks ago by Michele Busby1.9k

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

ADD REPLYlink written 7 weeks ago by mah190

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.

ADD REPLYlink written 7 weeks ago by Michele Busby1.9k
1
gravatar for JC
8 weeks ago by
JC7.7k
Mexico
JC7.7k 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.

ADD COMMENTlink written 8 weeks ago by JC7.7k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1450 users visited in the last hour