[BioC] Paired tests using limma and beadarray
Gordon K Smyth
smyth at wehi.EDU.AU
Tue Apr 19 00:37:49 CEST 2011
Dear Moritz,
Thanks for including your detailed code.
Yes, coef=19 stands for B-A. You could equally well write
topTable(ebFit,coef="treatmentB")
so there is no need to figure out that it is column 19.
There is something strange to me about your email, however, which is that
the code you give does not produce a contrast matrix, as there is no need
for it. So the contrast matrix you give must be from code that you do not
give here (and is not needed).
I personally recommend the preprocessing pipeline for Illumina data given
in Case Study 11.7 on page 91 of the limma User's Guide. The beadarray
code you give does not appear to be doing any background correction. You
could call "neqc" as an option in the normaliseIllumina() function in
beadarray, but that does not give the same results as calling neqc()
directly in limma.
Best wishes
Gordon
> Date: Sun, 17 Apr 2011 13:12:48 +0200
> From: Moritz Kebschull <endothel at gmail.com>
> To: bioconductor at r-project.org
> Subject: [BioC] Paired tests using limma and beadarray
> Message-ID: <BANLkTi=pYGOSU4OaD0js84xGgPmxx-udgA at mail.gmail.com>
> Content-Type: text/plain
>
> Dear list,
>
> I am looking at Illumina beadarrays of 18 patients sampled at baseline ("A")
> and after intervention ("B"), i.e. a total of 36 arrays. I think in this
> case, a paired assessment should be done. I am interested in diff. exp. at
> timepoint B over A.
>
> Thus far, I have read the data using the beadarray package, and tried to
> follow the limma manual's instructions for paired tests. I am, however,
> unsure about the contrasts, and my results look weird to me.
>
> I have done:
>
> # read and normalize data using beadarray
> library(beadarray)
> dataFile = "data.txt"
> qcFile = "control.txt"
> sampleSheet = "SampleSheet.csv"
> BSData = readBeadSummaryData(dataFile = dataFile, dec=",", qcFile = NULL,
> skip = 0, qc.skip = 0)
> BSData.quantile = normaliseIllumina(BSData, method = "quantile", transform =
> "log2")
> BSData.genes = BSData.quantile[which(fData(BSData)$Status == "Gene"), ]
> expressed = apply(Detection(BSData.genes) < 0.05, 1, any)
> BSData.filt = BSData.genes[expressed,]
>
> # set up paired tests using limma - what's wrong here?
> library(limma)
> patients = c(17, 2, 18, 5, 3, 8, 6, 14, 4, 13, 7, 15, 9, 10, 11, 12, 17, 18,
> 2, 5, 3, 8, 6, 14, 1, 12, 4, 13, 7, 15, 9, 16, 10, 1, 11, 16)
> patients = factor(patients)
> treatment = c(rep("B",7), "A", "B", "A", "B", "A", "B", rep("A",10), "B",
> "A", "B", "A", "B", "A", "B", "A", "A", rep("B", 4) )
> treatment = factor(treatment, levels = c("A", "B"))
> design = model.matrix(~patients+treatment)
> fit = lmFit(exprs(BSData.filt), design)
> ebFit=eBayes(fit)
>
> # left out some annotation steps
>
> # not sure what contrast is for what - does #19 stand for B-A??
> topTable(ebFit, coef=19)
>
>
> The primary issue is what contrast to choose - for normal, unpaired tests I
> would have set up a contrast matrix for "B-A". When I look at the different
> contrasts, I am missing patient #1.
>
>> ebFit$contrasts
> (Intercept) patients2 patients3 patients4 patients5 patients6
> (Intercept) 1 0 0 0 0 0
> patients2 0 1 0 0 0 0
> patients3 0 0 1 0 0 0
> patients4 0 0 0 1 0 0
> patients5 0 0 0 0 1 0
> patients6 0 0 0 0 0 1
> patients7 0 0 0 0 0 0
> patients8 0 0 0 0 0 0
> patients9 0 0 0 0 0 0
> patients10 0 0 0 0 0 0
> patients11 0 0 0 0 0 0
> patients12 0 0 0 0 0 0
> patients13 0 0 0 0 0 0
> patients14 0 0 0 0 0 0
> patients15 0 0 0 0 0 0
> patients16 0 0 0 0 0 0
> patients17 0 0 0 0 0 0
> patients18 0 0 0 0 0 0
> treatmentB 0 0 0 0 0 0
> patients7 patients8 patients9 patients10 patients11 patients12
> (Intercept) 0 0 0 0 0 0
> patients2 0 0 0 0 0 0
> patients3 0 0 0 0 0 0
> patients4 0 0 0 0 0 0
> patients5 0 0 0 0 0 0
> patients6 0 0 0 0 0 0
> patients7 1 0 0 0 0 0
> patients8 0 1 0 0 0 0
> patients9 0 0 1 0 0 0
> patients10 0 0 0 1 0 0
> patients11 0 0 0 0 1 0
> patients12 0 0 0 0 0 1
> patients13 0 0 0 0 0 0
> patients14 0 0 0 0 0 0
> patients15 0 0 0 0 0 0
> patients16 0 0 0 0 0 0
> patients17 0 0 0 0 0 0
> patients18 0 0 0 0 0 0
> treatmentB 0 0 0 0 0 0
> patients13 patients14 patients15 patients16 patients17
> patients18
> (Intercept) 0 0 0 0 0
> 0
> patients2 0 0 0 0 0
> 0
> patients3 0 0 0 0 0
> 0
> patients4 0 0 0 0 0
> 0
> patients5 0 0 0 0 0
> 0
> patients6 0 0 0 0 0
> 0
> patients7 0 0 0 0 0
> 0
> patients8 0 0 0 0 0
> 0
> patients9 0 0 0 0 0
> 0
> patients10 0 0 0 0 0
> 0
> patients11 0 0 0 0 0
> 0
> patients12 0 0 0 0 0
> 0
> patients13 1 0 0 0 0
> 0
> patients14 0 1 0 0 0
> 0
> patients15 0 0 1 0 0
> 0
> patients16 0 0 0 1 0
> 0
> patients17 0 0 0 0 1
> 0
> patients18 0 0 0 0 0
> 1
> treatmentB 0 0 0 0 0
> 0
> treatmentB
> (Intercept) 0
> patients2 0
> patients3 0
> patients4 0
> patients5 0
> patients6 0
> patients7 0
> patients8 0
> patients9 0
> patients10 0
> patients11 0
> patients12 0
> patients13 0
> patients14 0
> patients15 0
> patients16 0
> patients17 0
> patients18 0
> treatmentB 1
>
> Is contrast #19 the one I am looking for? Where is patient #1? The thing is
> that when looking at pt #2 - #18 seperately, I do find a couple of genes
> regulated pretty consistently, and these ones do not show up at all in the
> overall comparison (contrast 19). In fact, there's pretty much no difference
> using that contrast, which I find weird...
>
> Perhaps anyone has an idea what's wrong? Many thanks!
>
> Moritz (Bonn, Germany)
>
>
> sessionInfo()
> R version 2.12.1 (2010-12-16)
> Platform: i386-pc-mingw32/i386 (32-bit)
>
> locale:
> [1] LC_COLLATE=German_Germany.1252 LC_CTYPE=German_Germany.1252
> [3] LC_MONETARY=German_Germany.1252 LC_NUMERIC=C
> [5] LC_TIME=German_Germany.1252
>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods base
>
> other attached packages:
> [1] illuminaHumanv4.db_1.8.0 org.Hs.eg.db_2.4.6 RSQLite_0.9-4
>
> [4] DBI_0.2-5 AnnotationDbi_1.12.0 limma_3.6.9
>
> [7] beadarray_2.0.6 Biobase_2.10.0
______________________________________________________________________
The information in this email is confidential and intend...{{dropped:4}}
More information about the Bioconductor
mailing list