### R code from vignette source 'cqn.Rnw' ################################################### ### code chunk number 1: cqn.Rnw:7-8 ################################################### options(width=70) ################################################### ### code chunk number 2: load ################################################### library(cqn) library(scales) ################################################### ### code chunk number 3: data1 ################################################### data(montgomery.subset) dim(montgomery.subset) montgomery.subset[1:4,1:4] colnames(montgomery.subset) ################################################### ### code chunk number 4: data2 ################################################### data(sizeFactors.subset) sizeFactors.subset[1:4] ################################################### ### code chunk number 5: data3 ################################################### data(uCovar) head(uCovar) ################################################### ### code chunk number 6: checkdata ################################################### stopifnot(all(rownames(montgomery.subset) == rownames(uCovar))) stopifnot(colnames(montgomery.subset) == names(sizeFactors.subset)) ################################################### ### code chunk number 7: cqncall ################################################### cqn.subset <- cqn(montgomery.subset, lengths = uCovar$length, x = uCovar$gccontent, sizeFactors = sizeFactors.subset, verbose = TRUE) cqn.subset ################################################### ### code chunk number 8: cqnplot1 ################################################### par(mfrow=c(1,2)) cqnplot(cqn.subset, n = 1, xlab = "GC content", lty = 1, ylim = c(1,7)) cqnplot(cqn.subset, n = 2, xlab = "length", lty = 1, ylim = c(1,7)) ################################################### ### code chunk number 9: normalizedvalues ################################################### RPKM.cqn <- cqn.subset$y + cqn.subset$offset RPKM.cqn[1:4,1:4] ################################################### ### code chunk number 10: rpkmvalues ################################################### RPM <- sweep(log2(montgomery.subset + 1), 2, log2(sizeFactors.subset/10^6)) RPKM.std <- sweep(RPM, 1, log2(uCovar$length / 10^3)) ################################################### ### code chunk number 11: groups ################################################### grp1 <- c("NA06985", "NA06994", "NA07037", "NA10847", "NA11920") grp2 <- c("NA11918", "NA11931", "NA12003", "NA12006", "NA12287") ################################################### ### code chunk number 12: whGenes ################################################### whGenes <- which(rowMeans(RPKM.std) >= 2 & uCovar$length >= 100) M.std <- rowMeans(RPKM.std[whGenes, grp1]) - rowMeans(RPKM.std[whGenes, grp2]) A.std <- rowMeans(RPKM.std[whGenes,]) M.cqn <- rowMeans(RPKM.cqn[whGenes, grp1]) - rowMeans(RPKM.cqn[whGenes, grp2]) A.cqn <- rowMeans(RPKM.cqn[whGenes,]) ################################################### ### code chunk number 13: maplots ################################################### par(mfrow = c(1,2)) plot(A.std, M.std, cex = 0.5, pch = 16, xlab = "A", ylab = "M", main = "Standard RPKM", ylim = c(-4,4), xlim = c(0,12), col = alpha("black", 0.25)) plot(A.cqn, M.cqn, cex = 0.5, pch = 16, xlab = "A", ylab = "M", main = "CQN normalized RPKM", ylim = c(-4,4), xlim = c(0,12), col = alpha("black", 0.25)) ################################################### ### code chunk number 14: gcmaplots ################################################### par(mfrow = c(1,2)) gccontent <- uCovar$gccontent[whGenes] whHigh <- which(gccontent > quantile(gccontent, 0.9)) whLow <- which(gccontent < quantile(gccontent, 0.1)) plot(A.std[whHigh], M.std[whHigh], cex = 0.2, pch = 16, xlab = "A", ylab = "M", main = "Standard RPKM", ylim = c(-4,4), xlim = c(0,12), col = "red") points(A.std[whLow], M.std[whLow], cex = 0.2, pch = 16, col = "blue") plot(A.cqn[whHigh], M.cqn[whHigh], cex = 0.2, pch = 16, xlab = "A", ylab = "M", main = "CQN normalized RPKM", ylim = c(-4,4), xlim = c(0,12), col = "red") points(A.cqn[whLow], M.cqn[whLow], cex = 0.2, pch = 16, col = "blue") ################################################### ### code chunk number 15: edgeRconstructor ################################################### library(edgeR) d.mont <- DGEList(counts = montgomery.subset, lib.size = sizeFactors.subset, group = rep(c("grp1", "grp2"), each = 5), genes = uCovar) ################################################### ### code chunk number 16: cqn.Rnw:225-227 (eval = FALSE) ################################################### ## ## Not run ## d.mont$offset <- cqn.subset$glm.offset ################################################### ### code chunk number 17: edgeRdisp ################################################### design <- model.matrix(~ d.mont$sample$group) d.mont$offset <- cqn.subset$glm.offset d.mont.cqn <- estimateGLMCommonDisp(d.mont, design = design) ################################################### ### code chunk number 18: edgeRfit ################################################### efit.cqn <- glmFit(d.mont.cqn, design = design) elrt.cqn <- glmLRT(efit.cqn, coef = 2) topTags(elrt.cqn, n = 2) ################################################### ### code chunk number 19: cqn.Rnw:264-265 ################################################### summary(decideTestsDGE(elrt.cqn)) ################################################### ### code chunk number 20: edgeRstd ################################################### d.mont.std <- estimateGLMCommonDisp(d.mont, design = design) efit.std <- glmFit(d.mont.std, design = design) elrt.std <- glmLRT(efit.std, coef = 2) summary(decideTestsDGE(elrt.std)) ################################################### ### code chunk number 21: cqn.Rnw:327-328 (eval = FALSE) ################################################### ## cqn$y + cqn$offset ################################################### ### code chunk number 22: sessionInfo ################################################### toLatex(sessionInfo())