different results from limma using the same data set and two different methods
0
0
Entering edit mode
@gordon-smyth
Last seen 52 minutes ago
WEHI, Melbourne, Australia
At 11:40 AM 1/07/2004, xiaocui zhu wrote: >Hello all, I have a cDNA data set with 16 different conditions, each >with 3 pairs of dyeSwap measurements (so 6 repeats in total). I ranked >differentially expressed genes in each condition using two different >methods in Limma, and the two methods yielded different B values and >gene ranks. I do not know if this is normal or I did something wrong. It is normal and fully intentional. The two methods lead to the same estimated fold changes (fit$coef) and the same unscaled standard deviations (fit$stdev.unscaled) but to different residual standard deviations (fit$sigma) and different degrees of freedom (fit$df.residual). This is because, if you fit a model to all your arrays at once, the residuals from all the arrays are pooled to estimate the gene-wise standard deviations. Generally speaking, this is desirable because the extra arrays do provide extra information on the gene specific standard deviations and hence using all the arrays at once provides more reliable estimators of variability. Gordon >Here is how I did it: > >METHOD 1: Generate a MAList for each condition, and fit the linear model >with individual MAList separately. The code looks like the following: > >M.condition1<-read.delim("condition1-log2Ratio.txt", header=TRUE) >A.condition1<-read.delim("condition1-amplitude.txt", header=TRUE) >MA.condition1<-new("MAList", list(M=M.condition1, A=A.condition1)) >MA.condition1$genes<- read.delim("Genes.txt", header=TRUE) > > >fit.condition1<-lmFit(MA.condition1, design=c(1,-1,1,-1,1,-1)) >efit.condition1<-eBayes(fit.condition1) >output<-topTable(efit.condition1, number=16200, adjust="fdr") >write.table(output, file="condition1-ebayes-sep.txt", sep="\t") > > >METHOD 2: Generate a MAList for the entire data set (96 arrays in >total), and fit the linear model with the MAList. The code looks like >the following: > > >M.all<-read.delim("all-log2Ratio.txt", header=TRUE) >A.all<-read.delim("all-Amplitude.txt", header=TRUE) >MA.all<-new("MAList", list(M=M.all, A=A.all)) >MA.all$genes<-read.delim("Genes.txt", header=TRUE) > >treatments <- factor(c(1,1,1,1,1,1,2,2,2,2,2,2,3,3,3,3,3,3,4,4,4,4,4,4, > 5,5,5,5,5,5,6,6,6,6,6,6,7,7,7,7,7,7, > 8,8,8,8,8,8,9,9,9,9,9,9,10,10,10,10,10,10, > 11,11,11,11,11,11,12,12,12,12,12,12,13,13,13,13,13,13, > 14,14,14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16)) >design <- model.matrix(~ 0+treatments) >design<-edit(design) # change the table to reflect the dyeSwap design Multiplying the rows of the design matrix by a vector of 1s and -1s would be simpler and more reproducible than manual editing. >colnames(design) <- c("condtion1", " condtion2"," condtion3", > "condtion4", " condtion5"," condtion6", > "condtion7", " condtion8"," condtion9", > "condtion10", " condtion11"," condtion12", > > "condtion13", " condtion14"," condtion15", > "condtion16") > >fit <- lmFit(MA.lab, design) >efit<-eBayes(fit) ># To get the gene rank of condition1: >output<-topTable(efit, coef=1, number=16200,adjust="fdr") >write.table(output, file="condition1-ebayes-batch.txt", sep="\t") > > >The rank and the B value of the same gene in the two output files are >different. > >I may have done something incorrectly in METHOD 2, because some genes in >its output file have a different "A" value when compared to the average >of the six replicates. On the other hand, all genes in the output file >of METHOD 1 have the same "M" and "A" value when compared to the average >M and average A values of the six replicates. > >Any help would be greatly appreciated. Eventually I want to compare >across different conditions using METHOD 2. Since I am not sure if my >code in METHOD 2 is correct, I have to wait until I figure out what's >going on. > > >Thanks a lot! > >Xiaocui > >_______________________________________________ >Bioconductor mailing list >Bioconductor@stat.math.ethz.ch >https://www.stat.math.ethz.ch/mailman/listinfo/bioconductor
limma limma • 646 views
ADD COMMENT

Login before adding your answer.

Traffic: 657 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6