Search
Question: How is the mean-variance plot calculated?
0
gravatar for neverstop
2.4 years ago by
neverstop0
neverstop0 wrote:

I'm trying to understand how the mean-variance plot is calculated with the meanVarPlot() of EDASeq package. I tried to reproduce it manually with the code below but the result is very different.

Given the expression matrix (row=genes, columns=samples), I used the logarithm to transform the counts before computing mean and variance. Since there are zero counts, I considered the  maximum between 0 and the log of the counts. After that I calculated the row means and the row variances (the mean and variance of each gene across all samples).

Here's the code I used:

rm(list=ls())
library(yeastRNASeq)
library(EDASeq)
data(geneLevelData)
head(geneLevelData)
Data=newSeqExpressionSet(as.matrix(geneLevelData))
meanVarPlot(Data,log=T)
geneLevelDataM=as.matrix(geneLevelData)

means=apply(pmax(log(geneLevelDataM),0),1,mean)
variances=apply(pmax(log(geneLevelDataM),0),1,var)
plot(means,variances)

means=apply(pmax(log2(geneLevelDataM),0),1,mean)
variances=apply(pmax(log2(geneLevelDataM),0),1,var)
plot(means,variances)

 

ADD COMMENTlink modified 2.4 years ago by Wolfgang Huber13k • written 2.4 years ago by neverstop0
2
gravatar for Wolfgang Huber
2.4 years ago by
EMBL European Molecular Biology Laboratory
Wolfgang Huber13k wrote:

Have a look at the source, e.g. https://github.com/Bioconductor-mirror/EDASeq/blob/master/R/methods-SeqExpressionSet.R (search for "meanVarPlot").

Other ways for accessing the source of a Bioconductor package are: download the .tar.gz from the package homepage and unpack; via the subversion server; via the R command line, you can list each function's body (it's a bit wordy for S4 methods).

It looks like the function does something else than what its manual page says, regarding the log argument.

 

ADD COMMENTlink modified 2.4 years ago • written 2.4 years ago by Wolfgang Huber13k

The function does indeed do something else than what the manual page says, i.e., it takes the log after computing mean and variance.

Thank you both for spotting this incongruence! I will fix it in devel.

Just to add to Wolfgang's answer re: source code, another way, within R, is to type:

getMethod("meanVarPlot", "SeqExpressionSet")

 

ADD REPLYlink written 2.4 years ago by davide risso750
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.2.0
Traffic: 226 users visited in the last hour