[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