### R code from vignette source 'vignettes/VanillaICE/inst/doc/VanillaICE.Rnw' ################################################### ### code chunk number 1: setup ################################################### options(width=70) ################################################### ### code chunk number 2: data ################################################### library(VanillaICE) data(locusLevelData) ################################################### ### code chunk number 3: integerCopynumber ################################################### par(las=1) plot(locusLevelData[["copynumber"]][, 1]/100, pch=".", ylab="copy number", log="y") abline(h=1:3, col="grey70") ################################################### ### code chunk number 4: createLocusSet ################################################### oligoSet <- new("oligoSnpSet", copyNumber=log2(locusLevelData[["copynumber"]]/100), call=locusLevelData[["genotypes"]], callProbability=locusLevelData[["crlmmConfidence"]], annotation=locusLevelData[["platform"]]) oligoSet <- oligoSet[!is.na(chromosome(oligoSet)), ] oligoSet <- oligoSet[chromosome(oligoSet) < 3, ] ################################################### ### code chunk number 5: sdMethod ################################################### sds <- sd(oligoSet) ################################################### ### code chunk number 6: assignSds ################################################### cnConfidence(oligoSet) <- 1/sds ################################################### ### code chunk number 7: logscale ################################################### oligoSet <- order(oligoSet) hmmOpts <- hmm.setup(oligoSet, is.log=TRUE) ################################################### ### code chunk number 8: fit_van ################################################### ## first 2 chromosomes oligoSet <- oligoSet[chromosome(oligoSet) <= 2, ] fit.van <- hmm(oligoSet, hmmOpts) ################################################### ### code chunk number 9: fig2 ################################################### library(RColorBrewer) cols <- brewer.pal(5, "YlOrBr")[2:5] chr1 <- oligoSet[chromosome(oligoSet)==1,] fit.chr1 <- fit.van[fit.van$chrom == 1, ] ##fit.chr1 <- fit.van[fit.van$chrom==1, ] isHet <- snpCall(chr1)==2 par(las=1) plot(position(chr1), copyNumber(chr1), pch=".", cex=2, col="royalblue", ylab="log2 copy number") points(position(chr1)[isHet], copyNumber(chr1)[isHet], col="red", pch=".", cex=2) abline(h=log2(1:3), col="grey70") sts <- start(fit.chr1); ends <- end(fit.chr1) xx <- range(c(sts,ends)) y <- c(-1,-1,-0.9,-0.9) polygon(x=c(xx, rev(xx)), y=y, col="white") for(i in 1:nrow(fit.chr1)){ polygon(x=c(sts[i], ends[i], ends[i], sts[i]), y=y, col=cols[fit.chr1$state[i]], border=cols[fit.chr1$state[i]]) } legend("topleft", fill=cols, legend=hmmOpts$states, bty="n") ################################################### ### code chunk number 10: range2 ################################################### fit.van[2, ] ################################################### ### code chunk number 11: findOverlaps ################################################### mm <- matchMatrix(findOverlaps(fit.van, oligoSet)) markersInRange <- featureNames(oligoSet)[mm[,1]==2] ################################################### ### code chunk number 12: multipanelSingeLocus (eval = FALSE) ################################################### ## xyplot(cn ~ x | id, oligoSet, range=RangedDataObject, ylim=c(-0.5,4), panel=xypanel) ################################################### ### code chunk number 13: multipanelSingeLocus (eval = FALSE) ################################################### ## xyplot(cn ~ x | id, oligoSet, range=RangedDataObject, ylim=c(-0.5,4), scales=list(x="free"), ## panel=xypanel) ################################################### ### code chunk number 14: xyplotExample ################################################### ranges.altered <- fit.van[state(fit.van) != 4, ] xy.example <- xyplot(cn ~ x | range, oligoSet, range=ranges.altered, frame=2e6, scales=list(x="free"), ylim=c(-0.5,4), panel=xypanel, cex=0.4, pch=21, border="grey", ylab=expression(log[2]("copy number"))) ################################################### ### code chunk number 15: xyplotPrint ################################################### pdf("VanillaICE-xyplot.pdf", width=8, height=7) print(xy.example) dev.off() ################################################### ### code chunk number 16: ice ################################################### hmmOpts <- hmm.setup(oligoSet, ICE=TRUE, rohStates=c(FALSE, TRUE, TRUE, FALSE, FALSE)) res <- tryCatch(hmm(oligoSet, hmmOpts), error=function(e) "platform not supported") print(res) ## supported platforms VanillaICE:::icePlatforms() ################################################### ### code chunk number 17: iceIllustration ################################################### ann <- annotation(oligoSet) annotation(oligoSet) <- "genomewidesnp6" fit.ice <- hmm(oligoSet, hmmOpts) fit.ice annotation(oligoSet) <- ann ################################################### ### code chunk number 18: fig3 ################################################### fit.chr1 <- fit.ice[chromosome(fit.ice)==1, ] widths <- width(fit.chr1) fit.chr1 <- fit.chr1[order(widths,decreasing=TRUE),] par(las=1) plot(position(chr1), copyNumber(chr1), pch=".", ylab="log2 copy number", xlab="physical position", cex=2, col="royalblue") points(position(chr1)[isHet], copyNumber(chr1)[isHet], col="red", pch=".", cex=2) abline(h=log2(1:3), col="grey70") sts <- start(fit.chr1); ends <- end(fit.chr1) xx <- range(c(sts,ends)) y <- c(-1,-1,-0.9,-0.9) polygon(x=c(xx, rev(xx)), y=y, col="white") for(i in 1:nrow(fit.chr1)){ polygon(x=c(sts[i], ends[i], ends[i], sts[i]), y=y, col=cols[state(fit.chr1)[i]], border=cols[state(fit.chr1)[i]]) } legend("topleft", fill=cols, legend=hmmOpts$states, bty="n") ################################################### ### code chunk number 19: copyNumberOnly ################################################### cnSet <- new("CopyNumberSet", copyNumber=log2(locusLevelData[["copynumber"]]/100), annotation=locusLevelData[["platform"]]) cnSet <- cnSet[chromosome(cnSet) <= 2 & !is.na(chromosome(cnSet)), ] hmmOpts <- hmm.setup(cnSet, is.log=TRUE) fit.cn <- hmm(cnSet, hmmOpts) ################################################### ### code chunk number 20: genotypesOnly ################################################### snpSet <- new("SnpSet", call=locusLevelData[["genotypes"]], callProbability=locusLevelData[["crlmmConfidence"]], annotation=locusLevelData[["platform"]]) featureData(snpSet) <- addFeatureAnnotation(snpSet) snpSet <- snpSet[chromosome(snpSet) < 3, ] hmmOpts <- hmm.setup(snpSet) fit.gt <- hmm(snpSet, hmmOpts) ################################################### ### code chunk number 21: fig5 ################################################### fit.chr1 <- fit.gt[chromosome(fit.gt)==1, ] widths <- width(fit.chr1) fit.chr1 <- fit.chr1[order(widths,decreasing=TRUE),] gt <- ifelse(snpCall(chr1) == 1 | snpCall(chr1) == 3, 1, 0) par(las=1) plot(position(chr1), jitter(gt, amount=0.05), pch=".", ylab="", xlab="physical position", ylim=c(-3, 1.2), yaxt="n") ##points(position(chr1)[isHet], copyNumber(chr1)[isHet,], pch=".", ylab="log2 copy number", xlab="physical position", cex=2, col="red") axis(side=2, at=c(0,1), labels=c("AB", "AA or BB"), cex.axis=0.7) sts <- start(fit.chr1); ends <- end(fit.chr1) xx <- range(c(sts,ends)) y <- c(-1,-1,-0.5,-0.5) polygon(x=c(xx, rev(xx)), y=y, col="white") for(i in 1:nrow(fit.chr1)){ polygon(x=c(sts[i], ends[i], ends[i], sts[i]), y=y, col=cols[fit.chr1$state[i]], border=cols[fit.chr1$state[i]]) } legend("bottomleft", fill=cols, legend=hmmOpts$states, bty="n") ################################################### ### code chunk number 22: smoothing ################################################### if(require("DNAcopy")){ ##create an outlier copyNumber(cnSet)[50] <- 10 copyNumber(cnSet)[45:55] cnaObj <- CNA(genomdat=copyNumber(cnSet), chrom=chromosome(cnSet), maploc=position(cnSet), data.type="logratio", sampleid=sampleNames(cnSet)) smoothed.cnaObj <- smooth.CNA(cnaObj) copyNumber(cnSet) <- matrix(smoothed.cnaObj[, "NA06993"], nrow(cnSet), 1) copyNumber(cnSet)[50] } ################################################### ### code chunk number 23: VanillaICE.Rnw:469-470 ################################################### toLatex(sessionInfo())