September 17, 2010
setwd (“/ home/wanguan2000/array/Original Data”)
library (limma)
library ( “affy”)
targets <- readTargets ("targets.txt")
rownames (targets) <-targets $ Id
targets < br />
celpath <- file.path ("/ home/wanguan2000/array/Original Data / CEL /")
abatch <- ReadAffy (filenames = targets $ FileName, celfile. path = c elpath)
sampleNames (abatch)
eset <- rma (abatch)
pData (eset)
# genefilter
library (genefilter)
featureFilter (eset, require.entrez = TRUE, require.GOBP = FALSE, require.GOCC = FALSE, require.GOMF = FALSE, require.CytoBand = FALSE, remove.dupEntrez = TRUE, feature.exclude = “^ AFFX”) -> eset2
# design
ct <- factor (targets $ Character)
design <- model.matrix (~ 0 ct)
colnames (design) <- levels (ct)
colnames (design)
fit <- lmFit (eset2, design)
# cont.matrix
cont.matrix <- makeContrasts (NRvsR_c = recurrence_c -nonrecurrent_c, NRvsR_p = recurrence_p-nonrecurrent_p, NRvsR = (recurrence_c recurrence_p) - (nonrecurrent_c nonrecurrent_p), levels = design)
fit2 <- contrasts.fit (fit, cont.matrix)
fit2 <- eBayes (fit2)
# annotation
library (“annotate”)
# annotation (eset2)
library (“hgu133plus2.db”)
topTable (fit2, coef = 1, adjust = “BH”, n = 10)
volcanoplot (fit2, coef = 1, highlight = 8, names = getSYMBO L (fit2 $ genes $ ID, “hgu133plus2″), main = “recurrence_c-n onrecurrent_c”)
topTable (fit2, coef = 2, adjust = “BH”, n = 10)
volcanoplot (fit2, coef = 2, highlight = 8, names = getSYMBO L (fit2 $ genes $ ID, “hgu133plus2″), main = “recurrence_p-n onrecurrent_p”)
topTable (fit2, coef = 3, adjust = “BH”, n = 10)
volcanoplot (fit2, coef = 3, highlight = 8, names = getSYMBO L (fit2 $ genes $ ID, “hgu133plus2″), main = “(recurrence_c recurrence_p) – (nonrecurrent_c nonrecurrent_p)”)
# html toptable
tab = topTable (fit2, coef = 3, adjust = “BH”, n = 100)
genenames = as.character (tab $ ID)
l1 = getEG (genenames, “hgu133plus2″)
sym = getSYMBOL (genenames, “hgu133plus2″)
tab = data.frame (sym, signif (tab [, - 1], 3))
htmlpage (list (l1), othernames = tab, filename = “mycha2. html”, title = “HTML reprot”, table.center = T, table.head = c (“Entrez ID”, colnames (tab)))
browseURL (“mycha2.html”)
# html table GO pathway
library (“annaffy”)
topTable (fit2, coef = 3, adjust = “BH”, n = 100) -> geneID
atab = aafTableAnn (geneID $ ID , “hgu133plus2.db”, aaf.ha ndler ())
saveHTML (atab, file = “lishi.html”)
browseURL (“lishi.html”)
# multest
library (“multtest”)
match (colnames (exprs (eset2)), targets $ FileName) -> mt < br />
targets $ Character [mt] -> bimatix
exprs (eset2) -> mydata
# recurrence_c-nonrecurrent_c

c (rep (0, times = sum (bimatix == “recurrence_c”)), rep (1, times = sum (bimatix == “nonrecurrent_c”))) -> cl
resTp = mt.maxT (cbind (mydata [, bimatix == "recurrence_c"], myd ata [, bimatix == "nonrecurrent_c"]), classlabel = cl, B = 1000)
head (resTp, n = 100) -> myresTp
l1 = getEG (rownames (myresTp), “hgu133plus2″)
symp = getSYMBOL (rownames (myresTp), “hgu133plus2″) < br />
tabp = data.frame (rownames (myresTp), symp, myresTp)
htmlpage (list (l1), othernames = tabp, filename = “mycha2 ss.html”, title = “HTML reprot”, table.center = T, table.head = c (“Entrez ID”, colnames (tabp)))
browseURL (“mycha2ss.html”)
# recurrence_p-nonrecurrent_p
c (rep (0, times = sum (bimatix == “recurrence_p”)), rep (1, times = sum (bimatix == “nonrecurrent_p”) )) -> cl
resTp = mt.maxT (cbind (mydata [, bimatix == "recurrence_p"], myd ata [, bimatix == "nonrecurrent_p"]), classlabel = cl, B = 1000)
head (resTp, n = 100) -> myresTp
l1 = getEG (rownames (myresTp), “hgu133plus2″)

symp = getSYMBOL (rownames (myresTp), “hgu133plus2″)
tabp = data.frame (rownames (myresTp), symp, myresTp)
htmlpage (list (l1 ), othernames = tabp, filename = “mycha2 ss.html”, title = “HTML reprot”, table.center = T, table.head = c (“Entrez ID”, colnames (tabp)))
< br /> browseURL (“mycha2ss.html”)