### 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(oligoClasses) library(VanillaICE) library2(SNPchip) library2(Biobase) library2(IRanges) data(locusLevelData) ################################################### ### code chunk number 3: createLocusSet ################################################### cn <- log2(locusLevelData[["copynumber"]]/100) oligoSet <- new("oligoSnpSet", copyNumber=integerMatrix(cn), call=locusLevelData[["genotypes"]], callProbability=locusLevelData[["crlmmConfidence"]], annotation=locusLevelData[["platform"]], genome="hg19") oligoSet <- oligoSet[!is.na(chromosome(oligoSet)), ] oligoSet <- oligoSet[chromosome(oligoSet) < 3, ] oligoSet <- chromosomePositionOrder(oligoSet) ################################################### ### code chunk number 4: fit_van ################################################### fit.van <- hmm(oligoSet) fit.van <- fit.van[[1]] ################################################### ### code chunk number 5: fig2 ################################################### if(require(RColorBrewer)){ cols <- brewer.pal(6, "YlOrBr") } else cols <- rainbow(n=6) chr1 <- oligoSet[chromosome(oligoSet)==1,] fit.chr1 <- fit.van[chromosome(fit.van) == "chr1", ] isHet <- snpCall(chr1)==2 par(las=1) plot(position(chr1), copyNumber(chr1)/100, pch=".", cex=2, col="royalblue", ylab="log2 copy number") points(position(chr1)[isHet], copyNumber(chr1)[isHet]/100, 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 seq_len(length(fit.chr1))){ polygon(x=c(sts[i], ends[i], ends[i], sts[i]), y=y, col=cols[elementMetadata(fit.chr1)$state[i]], border=cols[elementMetadata(fit.chr1)$state[i]]) } legend("topleft", fill=cols, legend=c("hom-del", "hem-del", "N", "N-noHets", "sDup", "dDup"), bty="n") ################################################### ### code chunk number 6: range2 ################################################### fit.van[2, ] ################################################### ### code chunk number 7: findOverlaps ################################################### frange <- makeFeatureGRanges(oligoSet) markersInterval2 <- subsetByOverlaps(frange, fit.van[2, ]) markersInterval2 ################################################### ### code chunk number 8: multipanelSingleLocus (eval = FALSE) ################################################### ## xyplot2(cn ~ x | id, oligoSet, range=RangedDataObject, ylim=c(-0.5,4), panel=xypanel) ################################################### ### code chunk number 9: multipanelSingleLocus (eval = FALSE) ################################################### ## xyplot2(cn ~ x | id, oligoSet, range=RangedDataObject, ylim=c(-0.5,4), scales=list(x="free"), ## panel=xypanel) ################################################### ### code chunk number 10: xyplotExample ################################################### ranges.altered <- fit.van[!state(fit.van) %in% c(3,4) & numberProbes(fit.van) > 5, ] ##elementMetadata(ranges.altered)$sampleId <- sampleNames(oligoSet) xy.example <- xyplot2(cn ~ x | range, oligoSet, range=ranges.altered, frame=2e6, scales=list(x="free"), ylim=c(-1,3), panel=xypanel, cex.pch=0.4, pch=21, border="grey", ylab=expression(log[2]("copy number"))) ################################################### ### code chunk number 11: xyplotPrint ################################################### pdf("VanillaICE-xyplot.pdf", width=8, height=7) print(xy.example) dev.off() ################################################### ### code chunk number 12: ice ################################################### VanillaICE:::icePlatforms() ################################################### ### code chunk number 13: iceIllustration ################################################### ann <- annotation(oligoSet) annotation(oligoSet) <- "genomewidesnp6" fit.ice <- hmm(oligoSet, ICE=TRUE) fit.ice <- fit.ice[[1]] annotation(oligoSet) <- ann ################################################### ### code chunk number 14: fig3 (eval = FALSE) ################################################### ## fit.chr1 <- fit.ice[chromosome(fit.ice)=="chr1", ] ## widths <- width(fit.chr1) ## fit.chr1 <- fit.chr1[order(widths,decreasing=TRUE),] ## par(las=1) ## plot(position(chr1)/1e6, copyNumber(chr1)/100, pch=".", ylab="log2 copy number", xlab="physical position", cex=2, col="royalblue") ## points(position(chr1)[isHet]/1e6, copyNumber(chr1)[isHet]/100, 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 seq_len(length(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=c("hom-del","hem-del", "N", "N/no hets", "s-dup", "d-dup"), bty="n") ################################################### ### code chunk number 15: genotypesOnly ################################################### snpSet <- as(oligoSet, "SnpSet2") fit.gt <- hmm(snpSet, prGtHom=c(0.7, 0.99), normalIndex=1L, S=2L) fit.gt <- fit.gt[[1]] ################################################### ### code chunk number 16: fig5 ################################################### fit.chr1 <- fit.gt[chromosome(fit.gt)=="chr1", ] 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") cols <- cols[c(1, length(cols))] for(i in seq_len(length(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("bottomleft", fill=cols, legend=rev(c("region of homozygosity", "Normal")), bty="n") ################################################### ### code chunk number 17: readdata ################################################### bset <- readRDS(file.path(system.file("extdata", package="VanillaICE"), "bset.rds")) ################################################### ### code chunk number 18: arm ################################################### library(oligoClasses) arm <- getArm(bset) table(arm) ################################################### ### code chunk number 19: fit1p ################################################### b1p <- bset[arm=="chr1p", ] fit1p <- hmm(b1p, verbose=TRUE) ################################################### ### code chunk number 20: fig1p ################################################### par(mfrow=c(2,1), las=1) x <- position(b1p)/1e6 plot(x, lrr(b1p)/100, pch=".", col="grey") abline(h=-0.07547) plot(x, baf(b1p)/1000, pch=".", col="grey") abline(h=0.542) ################################################### ### code chunk number 21: fit1q ################################################### b1q <- bset[arm=="chr1q", ] fit1q <- hmm(b1q, verbose=TRUE) ################################################### ### code chunk number 22: fig1q ################################################### par(mfrow=c(2,1), las=1) x <- position(b1q)/1e6 plot(x, lrr(b1q)/100, pch=".", col="grey") abline(h=c(0, .359778)) plot(x, baf(b1q)/1000, pch=".", col="grey") abline(h=c(0.235, 0.790)) ################################################### ### code chunk number 23: VanillaICE.Rnw:463-464 ################################################### toLatex(sessionInfo())