################################################### ### code chunk number 5: compendium ################################################### library("ff") library("genefilter") library("IRanges") library("MASS") library("VanillaICE") library("crlmmCompendium") ################################################### ### code chunk number 15: loadObject-snr ################################################### if(!exists("SNR")) data("SNR") ################################################### ### code chunk number 12: snr ################################################### (snrfig <- histogram(~SNR, breaks=100)) ################################################### ### code chunk number 22: loadObject-genotypeSet ################################################### if(!exists("genotypeSet")) data("genotypeSet") ################################################### ### code chunk number 23: dataframeForClusterPlot ################################################### df <- prePredictPanel(genotypeSet) ################################################### ### code chunk number 24: genotypeColor ################################################### fill1 <- brewer.pal(3, "Set1")[df$gt] ################################################### ### code chunk number 25: confidenceColor ################################################### gt.conf <- df$gt.conf min.conf <- min(gt.conf) max.conf <- max(gt.conf) sc <- (gt.conf - min.conf)/(max.conf-min.conf) fill2 <- sapply(sc, grey) ################################################### ### code chunk number 36: loadObject-exampleData1 ################################################### if(!exists("exampleData1")) data("exampleData1") ################################################### ### code chunk number 37: dataCnFigs ################################################### a <- t(as.matrix(A(exampleData1))) gt <- t(as.matrix(calls(exampleData1))) nuA <- as.numeric(nu(exampleData1, "A")) phA <- as.numeric(phi(exampleData1, "A")) col <- brewer.pal(7, "Accent")[c(1, 4, 7)] NN <- Ns(exampleData1, i=1:16, j=1)[, , 1] fns <- featureNames(exampleData1) snpId <- matrix(fns, nrow(a), ncol(a), byrow=TRUE) snpId <- factor(snpId, levels=fns, ordered=TRUE) ##IMPORTANT ldat <- data.frame(A=as.integer(a), gt=as.factor(c("AA", "AB", "BB")[as.integer(gt)]), snp=snpId) boxplotfig <- bwplot(A~gt|snp, ldat, cex=0.6, panel=lmPanel, nu=nuA, ph=phA, fill="lightblue", Ns=NN, par.strip.text=list(lines=0.9, cex=0.6), ylab=expression(I[A]), xlab=expression(I[B]), ltext.y=2500, label.cex=0.6) ################################################### ### code chunk number 38: boxplotA ################################################### pars <- trellis.par.get() pars$axis.text$cex <- 0.7 pars$xlab.text$cex <- 0.8 pars$ylab.text$cex <- 0.8 trellis.par.set("axis.text", pars$axis.text) trellis.par.set("xlab.text", pars$xlab.text) trellis.par.set("ylab.text", pars$ylab.text) print(boxplotfig) ################################################### ### code chunk number 39: defineLatticeObjects ################################################### ldat <- prePredictPanel(exampleData1) shades <- makeTransparent(brewer.pal(6, "BrBG"), alpha=0.6)[c(1,2,3,5,6)] ##replace the middle color (white) with something else mykey <- simpleKey(as.character(0:4), points=FALSE, rectangles=TRUE, col="black", space="right", cex=0.7) mykey$rectangles[["col"]] <- shades (bvnfig <- xyplot(A~B|snp, ldat, cex=0.3, panel=cnPredictionPanel, object=exampleData1, x.axis="B", copynumber=0:4, line.col=shades, line.lwd=1.5, shades=shades, ylab=expression(log[2](I[A])), xlab=expression(log[2](I[B])), par.strip.text=list(lines=0.9, cex=0.6), xlim=c(6,12.5), ylim=c(6,12.5), key=mykey)) ################################################### ### code chunk number 40: ABscatterplots ################################################### pars <- trellis.par.get() pars$axis.text$cex <- 0.5 pars$xlab.text$cex <- 0.7 pars$ylab.text$cex <- 0.7 trellis.par.set("axis.text", pars$axis.text) trellis.par.set("xlab.text", pars$xlab.text) trellis.par.set("ylab.text", pars$ylab.text) print(bvnfig) ################################################### ### code chunk number 46: loadObject-redonSet ################################################### if(!exists("redonSet")) data("redonSet") ################################################### ### code chunk number 44: center ################################################### copyNumber(redonSet) <- copyNumber(redonSet) - median(copyNumber(redonSet), na.rm=TRUE) + 2 ################################################### ### code chunk number 47: HMM ################################################### hmmOpts <- hmm.setup(redonSet, c("hom-del", "hem-del", "normal", "amp"), copynumberStates=c(0:3), normalIndex=3, log.initialP=rep(log(1/4), 4), prGenotypeHomozygous=c(0.8, 0.99, 0.7, 0.75)) ################################################### ### code chunk number 48: viterbi ################################################### fit.cn <- hmm(redonSet, hmmOpts, verbose=FALSE, TAUP=1e10) hmm.df <- as.data.frame(fit.cn) print(hmm.df[, c(2:4,7:9)])