################################################### ### chunk number 1: setup ################################################### #line 50 "vignettes/VanillaICE/inst/doc/VanillaICE.Rnw" options(width=70) ################################################### ### chunk number 2: data ################################################### #line 111 "vignettes/VanillaICE/inst/doc/VanillaICE.Rnw" library(VanillaICE) data(locusLevelData) ################################################### ### chunk number 3: integerCopynumber ################################################### #line 121 "vignettes/VanillaICE/inst/doc/VanillaICE.Rnw" par(las=1) plot(locusLevelData[["copynumber"]][, 1]/100, pch=".", ylab="copy number", log="y") abline(h=1:3, col="grey70") ################################################### ### chunk number 4: createLocusSet ################################################### #line 131 "vignettes/VanillaICE/inst/doc/VanillaICE.Rnw" oligoSet <- new("oligoSnpSet", copyNumber=locusLevelData[["copynumber"]]/100, call=locusLevelData[["genotypes"]], callProbability=locusLevelData[["crlmmConfidence"]], annotation=locusLevelData[["platform"]]) oligoSet <- oligoSet[!is.na(chromosome(oligoSet)), ] ################################################### ### chunk number 5: eval=FALSE ################################################### ## #line 149 "vignettes/VanillaICE/inst/doc/VanillaICE.Rnw" ## sds <- robustSds(log2(locusLevelData[["copynumber"]]/100)) ################################################### ### chunk number 6: logscale ################################################### #line 171 "vignettes/VanillaICE/inst/doc/VanillaICE.Rnw" copyNumber(oligoSet) <- log2(copyNumber(oligoSet)) oligoSet <- oligoSet[order(chromosome(oligoSet), position(oligoSet)), ] ##trace(hmm.setup, browser) hmmOpts <- hmm.setup(oligoSet, copynumberStates=log2(c(1, 2, 2, 3)), states=c("hem-del", "ROH", "normal", "amp"), normalIndex=3, log.initialP=rep(log(1/4), 4), prGenotypeHomozygous=c(0.99, 0.99, 0.7, 0.7)) ################################################### ### chunk number 7: emission.gt ################################################### #line 184 "vignettes/VanillaICE/inst/doc/VanillaICE.Rnw" dim(hmmOpts$log.emission) ################################################### ### chunk number 8: fit_van ################################################### #line 196 "vignettes/VanillaICE/inst/doc/VanillaICE.Rnw" fit.van <- hmm(oligoSet, hmmOpts) ################################################### ### chunk number 9: fig2 ################################################### #line 202 "vignettes/VanillaICE/inst/doc/VanillaICE.Rnw" 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") ################################################### ### chunk number 10: ice ################################################### #line 235 "vignettes/VanillaICE/inst/doc/VanillaICE.Rnw" hmmOpts <- hmm.setup(oligoSet, ICE=TRUE, copynumberStates=log2(c(1, 2, 2, 3)), states=c("hem-del", "ROH", "normal", "amp"), normalIndex=3, log.initialP=rep(log(1/4), 4), rohStates=c(TRUE, TRUE, FALSE, FALSE)) fit.ice <- hmm(oligoSet, hmmOpts) ################################################### ### chunk number 11: fig3 ################################################### #line 247 "vignettes/VanillaICE/inst/doc/VanillaICE.Rnw" fit.chr1 <- fit.ice[fit.ice$chrom==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[fit.chr1$state[i]], border=cols[fit.chr1$state[i]]) } legend("topleft", fill=cols, legend=hmmOpts$states, bty="n") ################################################### ### chunk number 12: copyNumberOnly ################################################### #line 274 "vignettes/VanillaICE/inst/doc/VanillaICE.Rnw" cnSet <- new("CopyNumberSet", copyNumber=log2(locusLevelData[["copynumber"]]/100), annotation=locusLevelData[["platform"]]) cnSet <- cnSet[order(chromosome(cnSet), position(cnSet)), ] cnSet <- cnSet[!is.na(chromosome(cnSet)), ] hmmOpts <- hmm.setup(cnSet, copynumberStates=log2(c(0, 1, 2, 3)), states=c("hom-del", "hem-del", "normal", "amp"), normalIndex=3, log.initialP=rep(log(1/4), 4)) fit.cn <- hmm(cnSet, hmmOpts) ################################################### ### chunk number 13: genotypesOnly ################################################### #line 316 "vignettes/VanillaICE/inst/doc/VanillaICE.Rnw" snpSet <- new("SnpSet", call=locusLevelData[["genotypes"]], callProbability=locusLevelData[["crlmmConfidence"]], annotation=locusLevelData[["platform"]]) featureData(snpSet) <- addFeatureAnnotation(snpSet) fvarLabels(snpSet) snpSet <- snpSet[order(chromosome(snpSet), position(snpSet)), ] snpSet <- snpSet[!is.na(chromosome(snpSet)), ] hmmOpts <- hmm.setup(snpSet, states=c("ROH", "normal"), normalIndex=2, log.initialP=rep(log(1/2), 2), prGenotypeHomozygous=c(0.99,0.7), TAUP=5e7) fit.gt <- hmm(snpSet, hmmOpts) ################################################### ### chunk number 14: fig5 ################################################### #line 336 "vignettes/VanillaICE/inst/doc/VanillaICE.Rnw" fit.chr1 <- fit.gt[fit.gt$chrom==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") ################################################### ### chunk number 15: smoothing ################################################### #line 369 "vignettes/VanillaICE/inst/doc/VanillaICE.Rnw" 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] } ################################################### ### chunk number 16: ################################################### #line 409 "vignettes/VanillaICE/inst/doc/VanillaICE.Rnw" toLatex(sessionInfo())