yifanzhang — Apr 13, 2014, 10:22 PM
library(limma)
library(lattice)
prepareData <- function(myGenes, data, design) {
miniDat <- t(data[myGenes, ])
miniDat <- data.frame(gExp = as.vector(miniDat),
gene = factor(rep(colnames(miniDat), each =
nrow(miniDat)), levels = colnames(miniDat)))
miniDat <- suppressWarnings(data.frame(design, miniDat))
miniDat
}
stripplotIt <- function(myData, ...) {
stripplot(gExp ~ devStage | gene, myData,
jitter.data = TRUE, group = gType,
auto.key = TRUE, type = c('p', 'a'), grid = TRUE, ...)
}
prDat <- read.table("GSE4051_data.tsv")
prDes <- readRDS("GSE4051_design.rds")
x <- matrix(rnorm(1000 * 3), nrow = 1000)
obsVars <- apply(x, 1, var)
densityplot(~obsVars, n = 200)
wtDes <- subset(prDes, gType == "wt")
wtDat <- subset(prDat, select = prDes$gType == "wt")
wtDesMat <- model.matrix(~devStage, wtDes)
wtFit <- lmFit(wtDat, wtDesMat)
wtEbFit <- eBayes(wtFit)
(dsHits <- topTable(wtEbFit, coef = grep("devStage", colnames(coef(wtEbFit)))))
devStageP2 devStageP6 devStageP10 devStage4_weeks AveExpr
1440645_at 0.3990 0.1952 0.9200 3.961 6.528
1416041_at 0.1580 0.4797 0.3327 5.115 9.383
1425222_x_at 0.8820 0.7995 1.5488 5.532 7.028
1451635_at 1.3025 1.1900 2.0160 6.188 8.319
1429028_at -2.4433 -3.4073 -4.3105 -4.602 8.045
1422929_s_at -2.9117 -3.6182 -3.5472 -3.661 7.278
1424852_at 0.4575 0.2298 0.5740 3.979 7.454
1425171_at 0.9980 3.0530 5.2787 6.079 9.620
1451617_at 0.7255 2.5128 4.9837 6.685 8.817
1451618_at 0.6028 2.8903 5.0507 6.288 9.431
F P.Value adj.P.Val
1440645_at 425.4 1.588e-17 4.755e-13
1416041_at 195.5 1.522e-14 2.280e-10
1425222_x_at 173.4 4.348e-14 4.341e-10
1451635_at 157.3 1.013e-13 7.585e-10
1429028_at 148.8 1.646e-13 9.203e-10
1422929_s_at 146.9 1.844e-13 9.203e-10
1424852_at 143.2 2.290e-13 9.799e-10
1425171_at 138.8 3.002e-13 1.124e-09
1451617_at 136.5 3.485e-13 1.160e-09
1451618_at 134.2 4.032e-13 1.207e-09
hitdat <- prepareData(rownames(dsHits[c(3,6,9),]), wtDat, wtDes)
stripplotIt(hitdat)
hitdat1 <- prepareData(rownames(dsHits[3,]), wtDat, wtDes)
fit1 <- lm(gExp~devStage, hitdat1)
summary(fit1)
Call:
lm(formula = gExp ~ devStage, data = hitdat1)
Residuals:
Min 1Q Median 3Q Max
-0.4847 -0.2157 -0.0077 0.1781 0.8315
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.276 0.173 30.50 6.5e-15 ***
devStageP2 0.882 0.245 3.61 0.0026 **
devStageP6 0.800 0.245 3.27 0.0052 **
devStageP10 1.549 0.245 6.33 1.3e-05 ***
devStage4_weeks 5.532 0.245 22.61 5.3e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.346 on 15 degrees of freedom
Multiple R-squared: 0.977, Adjusted R-squared: 0.971
F-statistic: 159 on 4 and 15 DF, p-value: 4.31e-12
dsHits[3,]
devStageP2 devStageP6 devStageP10 devStage4_weeks AveExpr
1425222_x_at 0.882 0.7995 1.549 5.532 7.028
F P.Value adj.P.Val
1425222_x_at 173.4 4.348e-14 4.341e-10
#These two F stats and p-values are similar
cont.matrix <- makeContrasts(P10VsP6 = devStageP10 - devStageP6, fourweeksVsP10 = devStage4_weeks -
devStageP10, levels = wtDesMat)
Warning: Renaming (Intercept) to Intercept
wtFitCont <- contrasts.fit(wtFit, cont.matrix)
Warning: row names of contrasts don't match col names of coefficients
wtEbFitCont <- eBayes(wtFitCont)
cutoff <- 0.05
wtResCont <- decideTests(wtEbFitCont, p.value = cutoff, method = "global")
summary(wtResCont)
P10VsP6 fourweeksVsP10
-1 142 154
0 29776 29042
1 31 753
hits5 <- rownames(prDat)[which(wtResCont[, "P10VsP6"] == 0 & wtResCont[, "fourweeksVsP10"] ==
0)]
hits6 <- rownames(topTable(wtEbFit,
coef = grep("devStageP[26]", colnames(coef(wtEbFit))),
n = Inf, p.value=1e-04))
common <- intersect(hits5, hits6)
stripplotIt(prepareData(common[(length(common)-5):length(common)], wtDat, wtDes))
#These are probes that have some expression changes up to P6 and then hold steady all the way to 4_weeks