[BioC] Limma Script

Shawn Westaway westaway at ohsu.edu
Fri May 11 03:41:04 CEST 2007


Dear Listers,
I have a working Limma script below, but I don't know what the output
is comparing. Can anyone give me a clue? Also, you can see from my run
that i am having other syntax problems.
Thanks,
Shawn

Script:

dataMatrix <- exprs(lumi.N)
presentCount <- pData(featureData(lumi.N))$presentCount
selDataMatrix <- dataMatrix[presentCount > 0, ]
selProbe <- rownames(selDataMatrix)
sampleType <-
c('1st_litter','WT','2nd_litter','WT','1st_litter','2nd_litter','WT','1st_litter')
require(limma)
design <- model.matrix(~ factor(sampleType))
colnames(design) <- c('1st_litter', '2nd_litter', 'WT')
fit <- lmFit(selDataMatrix, design)
fit <- eBayes(fit)
if (require(lumiMouseV1) & require(annotate)) {
geneSymbol <- getSYMBOL(fit$genes$ID, 'lumiMouseV1')
fit$genes <- data.frame(fit$genes, geneSymbol=geneSymbol)
}
topTable(fit, coef='WT', adjust='fdr', number=10)
p.adj <- p.adjust(fit$p.value[,2])
sigGene.adj <- selProbe[ p.adj < 0.01]
sigGene <- selProbe[ fit$p.value[,2] < 0.001]
}


Transcript:

> dataMatrix <- exprs(lumi.N)
> presentCount <- pData(featureData(lumi.N))$presentCount
> selDataMatrix <- dataMatrix[presentCount > 0, ]
> selProbe <- rownames(selDataMatrix)
> sampleType <-
c('1st_litter','WT','2nd_litter','WT','1st_litter','2nd_litter','WT','1st_litter')
> require(limma)
[1] TRUE
> design <- model.matrix(~ factor(sampleType))
> colnames(design) <- c('1st_litter', '2nd_litter', 'WT')
> fit <- lmFit(selDataMatrix, design)
> fit <- eBayes(fit)
> if (require(lumiMouseV1) & require(annotate)) {
+ geneSymbol <- getSYMBOL(fit$genes$ID, 'lumiMouseV1')
+ fit$genes <- data.frame(fit$genes, geneSymbol=geneSymbol)
+ }
> topTable(fit, coef='WT', adjust='fdr', number=10)
                     ID geneSymbol      logFC         t      P.Value 
adj.P.Val
9576 ooK5UkVej.E737PSw0      Fkbp5  0.6698199  9.658213 6.216342e-06
0.06326372
5398 faEIr_9OZAUBsX0VLk       Per2  0.6880751  7.431683 4.825718e-05
0.24555665
9243 oyruyspXqhyohXojdA        Arc -1.0021833 -7.023328 7.395032e-05
0.25086413
3875 uefVK.bBeKqfv9XzuQ       Xbp1  0.4280023  6.613818 1.155933e-04
0.25790396
9581 HAyFSJK.n9pHhIQJS4      Dusp1 -1.0741752 -6.530455 1.269007e-04
0.25790396
6386 lKF2iniB8HnghGhHXk     Slc2a1  0.7937010  6.371030 1.520511e-04
0.25790396
9717 iUEl382FlILs4VFkFk      Tekt4  0.3838363  6.036845 2.244065e-04
0.32625504
8764 otfSbQ70_07x405Ptg        Fos -1.0311305 -5.863374 2.762056e-04
0.35136802
8467 r7fv79XdQPfX3_5FeU       Sox9 -0.6741969 -5.203197 6.317984e-04
0.71442364
7276 HbaB_NS3aeHQg7ogic     Polr3e  0.3346411  5.079054 7.431893e-04
0.75634378
              B
9576 -0.6906876
5398 -1.0774567
9243 -1.1777686
3875 -1.2909444
9581 -1.3156866
6386 -1.3647262
9717 -1.4753395
8764 -1.5372394
8467 -1.8044508
7276 -1.8608493
> p.adj <- p.adjust(fit$p.value[,2])
> sigGene.adj <- selProbe[ p.adj < 0.01]
> sigGene <- selProbe[ fit$p.value[,2] < 0.001]
> }
Error: syntax error
>



More information about the Bioconductor mailing list