Hello,

I've been asked to analyze a small set of illumina beadarray microarray
data. This is my first time dealing with beadarray data and I'm trying to
use lumi for the analysis. I've asked them to send me the data exported
from BeadStudio as explained in the lumi manual. here's how it looks:

TargetID PROBE_ID OP9Strom_DMSO_A.AVG_Signal OP9Strom_DMSO_A.BEAD_STDERR
OP9Strom_DMSO_A.Avg_NBEADS OP9Strom_DMSO_A.Detection Pval
OP9Strom_DMSO_B.AVG_Signal OP9Strom_DMSO_B.BEAD_STDERR
OP9Strom_DMSO_B.Avg_NBEADS OP9Strom_DMSO_B.Detection Pval
OP9Strom_DMSO_C.AVG_Signal OP9Strom_DMSO_C.BEAD_STDERR
OP9Strom_DMSO_C.Avg_NBEADS OP9Strom_DMSO_C.Detection Pval
OP9Strom_SB216_A.AVG_Signal OP9Strom_SB216_A.BEAD_STDERR
OP9Strom_SB216_A.Avg_NBEADS OP9Strom_SB216_A.Detection Pval
OP9Strom_SB216_B.AVG_Signal OP9Strom_SB216_B.BEAD_STDERR
OP9Strom_SB216_B.Avg_NBEADS OP9Strom_SB216_B.Detection Pval
OP9Strom_SB216_C.AVG_Signal OP9Strom_SB216_C.BEAD_STDERR
OP9Strom_SB216_C.Avg_NBEADS OP9Strom_SB216_C.Detection Pval
0610006I08RIK ILMN_2896528 168.3699 10.09909 29 0.5114213 212.3592 8.833445
42 0.2941177 281.087 12.216 43 0 161.7814 7.293227 38 0.5191327 210.6965
9.079482 42 0.2701665 326.766 13.8865 49 0.006393862
0610006I08RIK ILMN_2721178 189.1175 8.836789 40 0.1256345 225.4542 9.498487
38 0.09590793 223.0708 10.58419 43 0.003826531 177.2381 7.986325 40
0.1607143 212.4594 10.95498 37 0.2381562 306.3548 15.03577 51 0.03836317
0610007C21RIK ILMN_3033922 196.6734 9.167346 30 0.05456853 283.8329 16.99586
35 0 501.1899 28.64095 45 0 183.8619 7.665568 40 0.08035714 263.0788
12.38318 32 0 674.6611 29.22499 48 0
0610007C21RIK ILMN_3092673 271.1504 10.24042 45 0 735.3365 31.2681 45 0
1696.855 89.37159 45 0 273.522 8.028311 50 0 504.6748 14.43287 46 0 2348.904
98.40449 48 0
0610007P08RIK ILMN_2816356 179.2056 7.086736 44 0.2982233 188.1702 8.057949
26 0.7378517 177.0092 7.181882 36 0.5484694 168.3153 6.504077 27 0.3686225
216.793 8.503628 39 0.1677337 226.0183 10.03401 24 0.9104859
0610007P14RIK ILMN_2808939 140.6596 5.374153 36 0.9048223 235.4146 9.731593
49 0.03324808 417.4424 15.76849 33 0 146.7121 6.569923 44 0.8061224 197.1432
10.51995 37 0.5416133 678.2877 37.60738 35 0
0610007P22RIK ILMN_2634564 174.2713 10.35142 25 0.3908629 197.6312 9.78494
32 0.5997443 196.7326 16.95064 24 0.1364796 172.9638 9.338004 28 0.255102
202.1019 9.070997 37 0.4379001 287.9155 10.94573 42 0.1393862
0610008C08RIK ILMN_2737647 189.4498 7.198838 47 0.1218274 217.6156 9.236488
34 0.1969309 199.2474 9.024244 45 0.1020408 160.0271 6.664541 31 0.557398
210.8939 11.33168 29 0.2676056 283.2556 12.53874 37 0.1930946
0610009B22RIK ILMN_2734484 175.8993 10.71337 36 0.35533 214.3315 9.845508 42
0.2557545 209.3384 9.369405 35 0.03188775 150.3579 5.357093 41 0.755102
236.2385 11.68408 45 0.02816901 326.6299 14.28444 44 0.006393862


there are 6 datasets, 2 conditions in triplicate. This is what i've done on
lumi:

> filename <- "kimberly_mouse_ref8_onlySB_for_lumi.txt" #define the file
name
> allisa.raw <- lumiR(filename) #import the data as a lumi object
> alissa.quant <- lumiExpresso(allisa.raw)
Background Correction: bgAdjust
Variance Stabilizing Transform method: vst
Normalization method: quantile


Background correction ...
Perform bgAdjust background correction ...
There is no control probe information in the LumiBatch object!
 No background adjustment will be performed.
done.

Variance stabilizing ...
Perform vst transformation ...
2012-04-13 10:15:31 , processing array  1
2012-04-13 10:15:31 , processing array  2
2012-04-13 10:15:31 , processing array  3
2012-04-13 10:15:31 , processing array  4
2012-04-13 10:15:31 , processing array  5
2012-04-13 10:15:31 , processing array  6
done.

Normalizing ...
Perform quantile normalization ...
done.

Quality control after preprocessing ...
Perform Quality Control assessment of the LumiBatch object ...
done.
Warning messages:
1: In function (x.lumi, method = c("vst", "log2", "cubicRoot"), ifPlot =
FALSE,  :
  Too few probes are detectable based on detection p-values!

Iteration method will be used for VST.
2: In function (x.lumi, method = c("vst", "log2", "cubicRoot"), ifPlot =
FALSE,  :
  Too few probes are detectable based on detection p-values!

Iteration method will be used for VST.

> presentCount <- detectionCall(alissa.quant)
> selDataMatrix <- dataMatrix[presentCount > 0,]
> probeList <- rownames(selDataMatrix)
> sampleType <- c("DMSO", "DMSO", "DMSO", "SB216", "SB216", "SB216")
> model.matrix(~ factor(sampleType))
  (Intercept) factor(sampleType)SB216
1           1                       0
2           1                       0
3           1                       0
4           1                       1
5           1                       1
6           1                       1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$`factor(sampleType)`
[1] "contr.treatment"

> colnames(design) <- c("DMSOcontrol","DMSOvsSB216")
> fit <- lmFit(selDataMatrix, design)
> fit <- eBayes(fit)
> library("annotate")
> library("AnnotationDbi")
> library("lumiMouseAll.db")
> library(lumiMouseIDMapping)
> nuIDList <- IlluminaID2nuID(probeList, species="Mouse")
> nuID <-  as.character(nuIDListDF$nuID)
> geneSymbol <- getSYMBOL(nuID, 'lumiMouseAll.db')
> geneName <- sapply(lookUp(nuID, 'lumiMouseAll.db', 'GENENAME'),
function(x) x[1])
> fit$genes <- data.frame(ID= probeList, geneSymbol=geneSymbol,
geneName=geneName, stringsAsFactors=FALSE)
> topTable(fit)
                       ID    geneSymbol
        geneName DMSOcontrol DMSOvsSB216  AveExpr        F      P.Value
 adj.P.Val
ILMN_2914418 ILMN_2914418         Ncstn
       nicastrin    8.472922  0.08387813 8.514861 581.8616 1.393301e-12
2.364993e-11
ILMN_2693594 ILMN_2693594         Abcf1 ATP-binding cassette, sub-family F
(GCN20), member 1    8.645798  0.04869082 8.670143 570.0970 1.571149e-12
2.364993e-11
ILMN_2503329 ILMN_2503329        Zcchc9                zinc finger, CCHC
domain containing 9    8.482354  0.08195255 8.523331 569.2882 1.584320e-12
2.364993e-11
ILMN_2648189 ILMN_2648189         Klf10
Kruppel-like factor 10    8.400388  0.01474584 8.407761 557.8118
1.785863e-12 2.364993e-11
ILMN_2670361 ILMN_2670361         Setd4                              SET
domain containing 4    8.339616  0.16911942 8.424176 553.2727 1.873740e-12
2.364993e-11
ILMN_2819311 ILMN_2819311        Tdpoz2                       TD and POZ
domain containing 2    8.220126  0.10358821 8.271920 552.3219 1.892782e-12
2.364993e-11
ILMN_3162692 ILMN_3162692        Gm5820
 predicted gene 5820    8.329354  0.16889402 8.413801 548.8573 1.964110e-12
2.364993e-11
ILMN_2743660 ILMN_2743660 2310061C15Rik                           RIKEN
cDNA 2310061C15 gene    7.932102  0.19539448 8.029799 520.1216 2.694079e-12
2.364993e-11
ILMN_2626950 ILMN_2626950       Spata21
 spermatogenesis associated 21    7.954388  0.18081143 8.044794 512.4081
2.941224e-12 2.364993e-11
ILMN_2656543 ILMN_2656543        Klraq1                             KLRAQ
motif containing 1    7.964677  0.05158519 7.990470 512.0854 2.952128e-12
2.364993e-11


While things seem to work fine - besides the warning I get from lumiExpress
command, the adjusted p values are either 1 or 2.34e-11 there's no middle
ground and I feel like this might be wrong. Is is the lack of background
correction?

As I said I've never tried to analyze beadarray data before and I'm
wondering I'm doing something wrong in the process. Can anyone pitch in?

Thanks a lot!

Cheers,

> sessionInfo()
R version 2.14.2 (2012-02-29)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

locale:
[1] C/en_US.UTF-8/C/C/C/C

attached base packages:
[1] grDevices datasets  splines   graphics  utils     grid      stats
methods   base

other attached packages:
 [1] lumiMouseIDMapping_1.10.0 lumiMouseAll.db_1.16.0    org.Mm.eg.db_2.6.4
       annotate_1.32.3           AnnotationDbi_1.16.19     limma_3.10.3

 [7] lumi_2.6.0                nleqslv_1.9.3             methylumi_2.0.13
       Biobase_2.14.0            plyr_1.7.1                reshape2_1.2.1

[13] survival_2.36-12          RSQLite_0.11.1            DBI_0.2-5
        knitr_0.4                 gplots_2.10.1
KernSmooth_2.23-7
[19] caTools_1.12              bitops_1.0-4.1            gdata_2.8.2
        gtools_2.6.2              RColorBrewer_1.0-5        ggplot2_0.9.0


loaded via a namespace (and not attached):
 [1] BiocInstaller_1.2.1   IRanges_1.12.6        MASS_7.3-17
Matrix_1.0-5          Rcpp_0.9.10           affy_1.32.1
affyio_1.22.0
 [8] codetools_0.2-8       colorspace_1.1-1      dichromat_1.2-4
digest_0.5.2          evaluate_0.4.1.1      formatR_0.3-4
hdrcde_2.15
[15] highlight_0.3.1       lattice_0.20-6        memoise_0.1
mgcv_1.7-13           munsell_0.3           nlme_3.1-103
 parser_0.0-14
[22] preprocessCore_1.16.0 proto_0.3-9.2         scales_0.2.0
 stringr_0.6           tools_2.14.2          xtable_1.7-0
 zlibbioc_1.0.1
-- 
*Ahmet Z.*

	[[alternative HTML version deleted]]

