Hi Khadeeja, On 10/7/2011 9:42 AM, khadeeja ismail wrote: >> In particular, I suspect >> that you are using something like: >> topTable(X) >> You might try: >> topTable(X,coef=4) > > It worked with the small dataset. When I tried the same with a larger set (11 pairs and and some 150k genes) I'm getting different results from the two approaches (even with using topTable(X, coef=12) in approach 1). > I'm finding no significant genes from the first approach, but using the > second approach I'm finding some. Is it possible? I thought I was doing > the same thing in two different ways. Have I overlooked something again? I don't think you have overlooked anything, but there may be some misunderstanding. In your dummy data below, you have very few genes, and fewer samples. In the example you cite above, you use 150k genes and more samples. This will increase your power to detect differences in two ways. First, by increasing the number of genes, you can more accurately estimate an overall variance, which is then used in the eBayes() step to 'shrink' your observed variance towards this overall variance. This is one of the reasons that limma is so popular - by using information from all genes, you can increase the power to detect differences in individual genes. Second, when you increase the number of pairs you are using, your power to detect differences increases as well. This has nothing to do with limma per se; it is just that the power of a t-test increases as you increase the number of observations. Best, Jim > > > Thanks, > Khadeeja > > > > > Targets<-readTargets(file="TargetsSelBMI.csv", sep="\t") > > Pair<- factor(Targets$Twin) > Classes<- factor(Targets$HL) > design<- model.matrix(~Pair + Classes) > > fit<-lmFit(data,design) > fit.pair<-eBayes(fit) > > topGenes<-topTable(fit.pair, coef=12,number=2000) > > > > > > > > > ________________________________ > From: Sean Davis<sdavis2@mail.nih.gov> > > Cc: "bioconductor@r-project.org"<bioconductor@r-project.org> > Sent: Thursday, October 6, 2011 2:25 PM > Subject: Re: [BioC] Regression analysis problem > > >> Hi everyone, >> >> I am wondering if anyone can help me understand the following output from limma. >> >> The problem goes like this: >> Say I have methylation values for four genes from three twin pairs. I created some dummy data so that values of gene1 is clearly different within twin pairs. My interest here is to find out the genes that are significantly different within the pairs (i.e. between T1a and T1b, T2a and T2b and T3a and T3b) >> >> >> T1a T1b T2a T2b T3a T3b >> >> "gene1" 16 4 18 4 19 5 >> "gene2" 3 1 4 5 2 6 >> "gene3" 2 1 4 3 4 5 >> "gene4" 3 1 2 3 3 4 >> >> >> I created a design matrix: >> >> > Pair<-c("TWIN1","TWIN1","TWIN2","TWIN2","TWIN3","TWIN3") >> >Classes<-c(2,1,2,1,2,1) >> >BMI<-c(45,20,34,12,25,19) >> >Pair<-factor(Pair) >> >Classes<-factor(Classes) >> >design<- model.matrix(~Pair+Classes) >> >>> design >> (Intercept) PairTWIN2 PairTWIN3 Classes2 >> 1 1 0 0 1 >> 2 1 0 0 0 >> 3 1 1 0 1 >> 4 1 1 0 0 >> 5 1 0 1 1 >> 6 1 0 1 0 >> >> >> ********** >> >> ...and fitted a linear model to the data to get the following results: >> ************************************************************ >> ID X.Intercept. PairTWIN2 PairTWIN3 Classes2 AveExpr F >> 1 gene1 3.333333 1.0 2.0 1.333333e+01 11.000000 106.298724 >> 2 gene2 2.500000 2.5 2.0 -1.000000e+00 3.500000 8.745648 >> 3 gene3 1.333333 2.0 3.0 3.333333e-01 3.166667 7.430245 >> 4 gene4 2.000000 0.5 1.5 9.064933e-16 2.666667 4.799441 >> P.Value adj.P.Val >> 1 9.993042e-91 3.997217e-90 >> 2 4.683758e-07 9.367515e-07 >> 3 5.578117e-06 7.437489e-06 >> 4 7.186523e-04 7.186523e-04 >> ************************************************************* > You didn't show us the rest of the code. In particular, I suspect > that you are using something like: > > topTable(X) > > You might try: > > topTable(X,coef=4) > > Sean > >> As you can see, all the genes have adjusted p values less than 0.05. >> >> But if I just compute the differnces between the pairs >> >> T1 T2 T3 >> gene1 12 14 14 >> gene2 2 -1 -4 >> gene3 1 1 -1 >> gene4 2 -1 -1 >> >> >> ....and used the design matrix (1,1,1), >> gene 1 comes up as significantly different within the pairs. >> >> ID logFC t P.Value adj.P.Val B >> 1 gene1 13.3333333 10.6660452 5.234557e-06 2.093823e-05 46.016217 >> 2 gene2 -1.0000000 -0.7999534 4.468388e-01 8.936777e-01 -5.513313 >> 3 gene3 0.3333333 0.2666511 7.964821e-01 1.000000e+00 -5.772418 >> 4 gene4 0.0000000 0.0000000 1.000000e+00 1.000000e+00 -5.804806 >> >> >> So, my questions are: >> Why am I not getting the same (expected) result from my first approach? Or am I misinterpreting the output? >> >> What does the p value in the 1st approach signify? >> >> And most importantly, am I applying the right model? >> >> Any help on this would be most appreciated. >> >> Thanks in advance, >> >> Khadeeja >> [[alternative HTML version deleted]] >> >> >> _______________________________________________ >> Bioconductor mailing list >> Bioconductor@r-project.org >> https://stat.ethz.ch/mailman/listinfo/bioconductor >> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor >> > [[alternative HTML version deleted]] > > > > _______________________________________________ > Bioconductor mailing list > Bioconductor@r-project.org > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor -- James W. MacDonald, M.S. Biostatistician Douglas Lab University of Michigan Department of Human Genetics 5912 Buhl 1241 E. Catherine St. Ann Arbor MI 48109-5618 734-615-7826 ********************************************************** Electronic Mail is not secure, may not be read every day, and should not be used for urgent or sensitive issues [[alternative HTML version deleted]]