# Analysis of the METH experiment library(limma) # Already run # RGnormexp50 <- backgroundCorrect(RG, method="normexp", offset=50) # use normal-exponential convolution model to background correct the data # Load data load("meth.Rdata") # Print-tip normalise the data MA <- normalizeWithinArrays(RGnormexp50) # Estimate array weights array.wts <- arrayWeights(MA, design, weights=NULL) round(array.wts,2) # [1] 0.29 0.96 1.31 0.36 1.10 1.02 2.04 1.31 1.53 1.67 # Linear model analysis # Make a matrix of array weights for lmFit array.wts <- matrix(rep(array.wts, each=dim(MA)[1]), dim(MA)[1], dim(MA)[2]) # Fit linear model with empirical array weights out.array.wts <- lmFit(MA, design, weights=array.wts) out.array.wts <- eBayes(out.array.wts) # Fit linear model after having removed suspect arrays (1st and 4th) # Set-up new design matrix design.indirect <- design[c(-1,-4),] out.remove2 <- lmFit(MA$M[,c(-1,-4)], design.indirect, weights=NULL) out.remove2 <- eBayes(out.remove2) # Set up contrasts of interest contrast.matrix <- cbind(d1vs0=c(1,0,-1),d3vs0=c(0,1,-1)) contrasts.array.wts <- contrasts.fit(out.array.wts, contrast.matrix) contrasts.array.wts <- eBayes(contrasts.array.wts) contrasts.remove2 <- contrasts.fit(out.remove2, contrast.matrix) contrasts.remove2 <- eBayes(contrasts.remove2) # Rank genes for contrasts of interest toptab.array.wts.1vs0 <- toptable(contrasts.array.wts, 1, 10368, genelist=RG$genes,adjust.method="fdr") toptab.array.wts.3vs0 <- toptable(contrasts.array.wts, 2, 10368, genelist=RG$genes,adjust.method="fdr") toptab.remove2.1vs0 <- toptable(contrasts.remove2, 1, 10368, genelist=RG$genes,adjust.method="fdr") toptab.remove2.3vs0 <- toptable(contrasts.remove2, 2, 10368, genelist=RG$genes,adjust.method="fdr") # Number of candidate DE genes by method sum(toptab.array.wts.1vs0$P.Value<0.05, na.rm=TRUE) # [1] 654 sum(toptab.array.wts.3vs0$P.Value<0.05, na.rm=TRUE) # [1] 1790 sum(toptab.remove2.1vs0$P.Value<0.05, na.rm=TRUE) # [1] 0 sum(toptab.remove2.3vs0$P.Value<0.05, na.rm=TRUE) # [1] 1263 length(intersect(toptab.array.wts.1vs0$Name[toptab.array.wts.1vs0$P.Value<0.05 & !is.na(toptab.array.wts.1vs0$P.Value)], toptab.array.wts.3vs0$Name[toptab.array.wts.3vs0$P.Value<0.05 & !is.na(toptab.array.wts.3vs0$P.Value)])) # [1] 413 length(intersect(toptab.remove2.1vs0$Name[toptab.remove2.1vs0$P.Value<0.05 & !is.na(toptab.remove2.1vs0$P.Value)], toptab.remove2.3vs0$Name[toptab.remove2.3vs0$P.Value<0.05 & !is.na(toptab.remove2.3vs0$P.Value)])) # [1] 0 # Number of genes responding in a consistent manner index <- toptab.array.wts.1vs0$P.Value<0.05 gene.names <- toptab.array.wts.1vs0$Name[index] gene.names <- gene.names[!is.na(gene.names)] index1vs0 <- match(gene.names, toptab.array.wts.1vs0$Name) index3vs0 <- match(gene.names, toptab.array.wts.3vs0$Name) sum(abs(toptab.array.wts.1vs0$M[index1vs0])