September 25, 2010
setwd (“/ home/wanguan2000/array/Original Data”) library (limma) library (“affy”) targets <- readTargets ("targetpairs.txt") targets 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) exprs (eset) -> myexpre # do not normalizeWithinArrays # M = log2R-log2G M = as.matrix (log2 (myexpre [, sampleNames (eset) [i]] / my expre [, sampleNames (eset) [ i 1]])) # A = (log2R log2G) / 2 A = as.matrix (log2 (sqrt (myexpre [, sampleNames (eset) [i]] * myexpre [, sampleNames (eset) [i 1 ]]))) myMAList <- new ("MAList") for (i in seq (1,140,2)) {M = as.matrix (log2 (myexpre [, sampleNames (eset) [i]] / my expre [, sampleNames (eset) [i 1]])) A = as.matrix (log2 (sqrt (myexpre [, sampleNames (eset) [i]] * myexpre [, sampleNames (eset) [i 1 ]]))) colnames (M) <- colnames (A) <- sub (".*_","", sampleNames (eset) [i]) MA1 <- new ("MAList", list (M = M, A = A) ) cbind (myMAList, MA1) -> myMAList} myMAList $ genes <- matrix (featureNames (eset), ncol = 1, dimnames = list (c (), c ("affyID"))) # fit targets <- readTargets ( "targetpairs.txt") targets [seq (1,140,2), c ("FileName", "Type")] -> mytargets Type <- factor (mytargets $ Type) design <- model.matrix (~ 0 Type) colnames (design) <- levels (Type) rownames (design) <- sub (".*_","", sampleNames (eset) [seq (1,140,2)]) fit <- lmFit (myMAList, design) cont . matrix <- makeContrasts (NvsR = nonrecurrent-recurrence, levels = design) fit2 <- contrasts.fit (fit, cont.matrix) fit2 <- eBayes (fit2) topTable (fit2, coef = 1, adjust = "BH", n = 10) -> tab boxplot (myMAList $ M ~ col (myMAList $ M), names = colnames (myMAList $ M)) # use normalizeWithinArrays myRGList <- new ("RGList") for (i in seq (1,140,2 )) {R = as.matrix (myexpre [, sampleNames (eset) [i]]) rownames (R) <- NULL Rb = as.matrix (rep (0.001, times = length ((myexpre [, samp leNames (eset ) [i ]])))) G = as.matrix (myexpre [, sampleNames (eset) [i 1]]) rownames (G) <- NULL Gb = as.matrix (rep (0.001, times = length ( (myexpre [, samp leNames (eset) [i ]])))) colnames (Rb) <- colnames (Gb) <- colnames (R) <- colnames (G) <- sub (".*_"," ", sampleNames (eset) [i]) RG1 <- new (" RGList ", list (R = R, G = G, Rb = Rb, Gb = Gb)) cbind (myRGList, RG1) -> myRGList} myRGList $ genes <- matrix (featureNames (eset), ncol = 1, dimnames = list (c (), c ("affyID"))) MA <- normalizeWithinArrays (myRGList, method = "none") # fit targets <- readTargets ( "targetpairs.txt") targets [seq (1,140,2), c ("FileName", "Type")] -> mytargets Type <- factor (mytargets $ Type) design <- model.matrix (~ 0 Type) colnames (design) <- levels (Type) rownames (design) <- sub (".*_","", sampleNames (eset) [seq (1,140,2)]) fit <- lmFit (MA, design) cont . matrix <- makeContrasts (NvsR = nonrecurrent-recurrence, levels = design) fit2 <- contrasts.fit (fit, cont.matrix) fit2 <- eBayes (fit2) topTable (fit2, coef = 1, adjust = "BH", n = 100) -> tab # boxplot (myMAList $ M ~ col (myMAList $ M), names = colnames (myMAList $ M)) library (“annotate”) # annotation (eset2) library (“hgu133plus2.db”) genenames = as.character (tab $ affyID) 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”)