Unexpectedly high FDR using contrasts in Limma?
0
0
Entering edit mode
@jdelasherasedacuk-1189
Last seen 8.7 years ago
United Kingdom
Hi everyone, I'm analysing (limma 2.9.10) a set of twelve 2-colour cDNA arrays (4 experiments of 3 slides each) using a common reference design. I find that I get very high adjusted P values (BH) using contrasts to compare some of these experiments. No adjusted P value is lower than 0.9, which would indicate there's not enough evidence that any gene behaves different in either experiment, and I find that surprising. Yes, all four experiments are quite similar, but there are some genes that do stand out enough (I'd think!) and I have confirmation by RT that this is the case. So I am wondering whether I don't fully understand where these adjusted p values come from, or whether I'm not asking the question I think I am asking, or maybe an error on my code? Here's a summary of what I did: - 12 files read with read.maimages > targets #H226 is the common reference SlideNumber FileName Cy3 Cy5 1 25 140025-3.gpr H226 MBD1 2 26 140026-3.gpr MBD1 H226 3 23 140023-3.gpr H226 MBD1 4 19 140019-3.gpr MBD2 H226 5 17 140017-3.gpr H226 MBD2 6 15 140015-3.gpr MBD2 H226 7 20 140020-3.gpr H226 MeCP2 8 18 140018-3.gpr MeCP2 H226 9 16 140016-3.gpr MeCP2 H226 10 22 140022-3.gpr H226 DBL 11 24 140024-3.gpr DBL H226 12 14 140014-3.gpr DBL H226 - background corrected using a custom method: RGList object 'RG.s' > design # H226 as common reference DBL MBD1 MBD2 MeCP2 1 0 1 0 0 2 0 -1 0 0 3 0 1 0 0 4 0 0 -1 0 5 0 0 1 0 6 0 0 -1 0 7 0 0 0 1 8 0 0 0 -1 9 0 0 0 -1 10 1 0 0 0 11 -1 0 0 0 12 -1 0 0 0 > MA.sw<-normalizeWithinArrays(RG.s, layout, bc.method="none", > method="printtiploess") > fit.sw<-lmFit(MA.sw, design, method="ls") > eb.sw<-eBayes(fit.sw, proportion=0.01) - I get the adjusted P values from 'topTable': > tt.sw.MBD1<-topTable(eb.sw, coef=2, n=number, genelist=gal, > adjust.method="BH", sort.by="P") > tt.sw.MBD2<-topTable(eb.sw, coef=3, n=number, genelist=gal, > adjust.method="BH", sort.by="P") > tt.sw.MeCP2<-topTable(eb.sw, coef=4, n=number, genelist=gal, > adjust.method="BH", sort.by="P") > tt.sw.DBL<-topTable(eb.sw, coef=1, n=number, genelist=gal, > adjust.method="BH", sort.by="P") The values I get here look alright and make sense. From the 4 experiments, the 4th one (DBL) is a control experiment. I would like to compare the other three to it, to see what genes are differentially expressed between each of the first three experiments and the control. I wanted to use contrasts for this. This is how I did it: > levels<-c("DBL", "MBD1", "MBD2", "MeCP2") > contrasts<-makeContrasts(DBL,MBD1,MBD2,MeCP2,DBL-MBD1,DBL-MBD2,DBL- MeCP2,levels=levels) > contrasts Contrasts Levels DBL MBD1 MBD2 MeCP2 DBL - MBD1 DBL - MBD2 DBL - MeCP2 DBL 1 0 0 0 1 1 1 MBD1 0 1 0 0 -1 0 0 MBD2 0 0 1 0 0 -1 0 MeCP2 0 0 0 1 0 0 -1 > fitc<-contrasts.fit(fit.sw,contrasts) > ebfitc<-eBayes(fitc) then I use 'topTable' again to get the adjusted p values etc: > tt.cMBD1<-topTable(ebfitc, coef=5, n=number, genelist=gal, > adjust.method="BH", sort.by="P") > tt.cMBD2<-topTable(ebfitc, coef=6, n=number, genelist=gal, > adjust.method="BH", sort.by="P") > tt.cMeCP2<-topTable(ebfitc, coef=7, n=number, genelist=gal, > adjust.method="BH", sort.by="P") the top of the list makes sense, I picked up the genes I was expecting, however the adjusted p values are terrible, so I wonder if this is real or I am making a mistake somewhere. As an example, here's a gene I know to be differentially expressed: CDKN1C This gene has negligible expression in the common reference, it has negligible expression after the control experiment (DBL) is performed, but it is clearly expressed after the MBD1 experiment is performed (verified by RT). This is what I get: A.mean = 10.48 (it goes up to 11.6 in the individual slides of the MBD1 experiment, it stays low in all others) M (MBD1) = 2.98 P (MBD1) = 6.07e-05 adj p val (MBD1) = 0.00044 M (DBL) = -0.44 P (DBL) = 0.38 adj p val (DBL) = 0.56 M (DBL-MBD1) = -3.43 P (DBL-MBD1) = 0.00036 adj p val (DBL-MBD1) = 0.95 I am surprised this gene (and several others like this one) give me such poor adjusted p values... why? Jose -- Dr. Jose I. de las Heras Email: J.delasHeras at ed.ac.uk The Wellcome Trust Centre for Cell Biology Phone: +44 (0)131 6513374 Institute for Cell & Molecular Biology Fax: +44 (0)131 6507360 Swann Building, Mayfield Road University of Edinburgh Edinburgh EH9 3JR UK
• 731 views
ADD COMMENT

Login before adding your answer.

Traffic: 918 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