################################################### ### chunk number 1: setup ################################################### options(width=80) ################################################### ### chunk number 2: assaydata ################################################### library(VanillaICE) data(locusLevelData) copynumber <- locusLevelData[["copynumber"]]/100 genotypes <- locusLevelData[["genotypes"]] ################################################### ### chunk number 3: featuredata ################################################### chromosome <- locusLevelData[["annotation"]][, "chromosome"] position <- locusLevelData[["annotation"]][, "position"] names(position) <- names(chromosome) <- rownames(locusLevelData[["annotation"]]) ################################################### ### chunk number 4: chromosomeAnnotation ################################################### require(SNPchip) data(chromosomeAnnotation, package="SNPchip") ################################################### ### chunk number 5: checks ################################################### all(c(identical(names(position), rownames(copynumber)), identical(names(position), rownames(genotypes)), identical(colnames(genotypes), colnames(copynumber)))) ################################################### ### chunk number 6: order ################################################### ordering <- order(chromosome, position) chromosome <- chromosome[ordering] position <- position[ordering] genotypes <- genotypes[ordering, , drop=FALSE] copynumber <- copynumber[ordering, , drop=FALSE] ################################################### ### chunk number 7: TESTING eval=FALSE ################################################### ## ##For testing that the arm variable forces the HMM to run independently on each chromosomal arm ## data(locusLevelData, package="VanillaICE") ## cnConf <- cnConfidence(chromosome1) ## callsConf <- callsConfidence(chromosome1) ## callsConf <- matrix(as.integer(-1000*log(1-callsConf)), nrow(callsConf), ncol(callsConf)) ## callsConf <- rbind(callsConf, callsConf[1:100, , drop=FALSE]) ## dimnames(callsConf) <- dimnames(copynumber) ## locusLevelData[["crlmmConfidence"]] <- callsConf ## ## annotation <- locusLevelData[["annotation"]] ## copynumber <- locusLevelData[["copynumber"]]/100 ## genotypes <- locusLevelData[["genotypes"]] ## i <- order(annotation[, "chromosome"], annotation[, "position"]) ## genotypes <- genotypes[i, , drop=FALSE] ## copynumber <- copynumber[i, , drop=FALSE] ## annotation <- annotation[i, , drop=FALSE] ## copynumber <- copynumber[1:9165, , drop=FALSE] ## genotypes <- genotypes[1:9165, , drop=FALSE] ## annotation <- annotation[1:9165,, drop=FALSE] ## copynumber <- rbind(copynumber, copynumber) ## genotypes <- rbind(genotypes, genotypes) ## annotation2 <- annotation ## annotation2[, 1] <- 2 ## rownames(annotation2) <- paste("SNP_A-", 1:nrow(annotation2), sep="") ## annotation <- rbind(annotation, annotation2) ## rownames(genotypes) <- rownames(copynumber) <- rownames(annotation) ################################################### ### chunk number 8: annotatedDataFrame ################################################### require(oligoClasses) locusAnnotation <- data.frame(list(chromosome=chromosome, position=position), row.names=names(chromosome)) featuredata <- new("AnnotatedDataFrame", data=locusAnnotation, varMetadata=data.frame(labelDescription=colnames(locusAnnotation))) stopifnot(all(c("chromosome", "position") %in% varLabels(featuredata))) ################################################### ### chunk number 9: instantiateCallSet eval=FALSE ################################################### ## new("SnpCallSet", ## calls=genotypes, ## phenoData=annotatedDataFrameFrom(genotypes, byrow=FALSE), ## featureData=featuredata) ################################################### ### chunk number 10: instantiateCopyNumberSet eval=FALSE ################################################### ## new("SnpCopyNumberSet", ## copyNumber=copynumber, ## phenoData=annotatedDataFrameFrom(copynumber, byrow=FALSE), ## featureData=featuredata) ################################################### ### chunk number 11: instantiateOligoSnpSet ################################################### locusset <- new("oligoSnpSet", copyNumber=copynumber, calls=genotypes, phenoData=annotatedDataFrameFrom(copynumber, byrow=FALSE), featureData=featuredata) ################################################### ### chunk number 12: orderClass ################################################### locusset <- locusset[order(chromosome(locusset), position(locusset)), ] ################################################### ### chunk number 13: vanillaSegments ################################################### ##featuredata <- new("AnnotatedDataFrame", data=locusAnnotation, varMetadata=data.frame(labelDescription=colnames(locusAnnotation))) ##locusset <- new("oligoSnpSet", ## copyNumber=copynumber, ## calls=genotypes, ## phenoData=annotatedDataFrameFrom(copynumber, byrow=FALSE), ## featureData=featuredata) (brks <- hmm(object=locusset, states=c("hemDel", "normal", "ROH", "amplification"), mu=c(1, 2, 2, 3), probs=c(0.999, 0.7, 0.999, 0.7), takeLog=TRUE)) ################################################### ### chunk number 14: vanillaMatrix ################################################### fit2 <- hmm(object=locusset, states=c("hemDel", "normal", "ROH", "amplification"), mu=c(1, 2, 2, 3), probs=c(0.999, 0.7, 0.9999, 0.7), takeLog=TRUE, returnSegments=FALSE) ################################################### ### chunk number 15: vanillaPlot ################################################### show(plot(locusset[chromosome(locusset) == 1, ])) fit2[fit2 == 3] <- 2 ##So LOH line will not plot at 3 fit2[fit2 == 4] <- 3 lines(position(locusset)[chromosome(locusset) == 1], fit2[chromosome(locusset) == 1, 1]) rohRegion <- as.integer(brks[brks$state == "ROH", ][c("start", "end")]) abline(v=rohRegion, col="royalblue", lty=2) legend("topright", lty=c(1, 2), col=c("black", "royalblue"), legend=c("CN state", "ROH boundary"), bty="n") ################################################### ### chunk number 16: passingAnEnvironment ################################################### env <- new.env() fit2 <- hmm(object=locusset, states=c("hemDel", "normal", "ROH", "amplification"), mu=c(1, 2, 2, 3), probs=c(0.999, 0.7, 0.9999, 0.7), takeLog=TRUE, returnSegments=FALSE, envir=env) ls(env) ################################################### ### chunk number 17: plotSnp eval=FALSE ################################################### ## chr1 <- annotation[, "chromosome"] == 1 ## plot(annotation[chr1, "position"], copynumber[chr1, 1], pch=".", cex=2, col="grey60") ################################################### ### chunk number 18: addCalls eval=FALSE ################################################### ## plot(annotation[chr1, "position"], copynumber[chr1, 1], pch=".", cex=2, col="grey60") ## het <- genotypes[, 1] == 2 ## points(annotation[chr1 & !het, "position"], copynumber[chr1 & !het, 1], pch=".", cex=2, col="royalblue") ## points(annotation[chr1 & het, "position"], copynumber[chr1 & het, 1], pch=".", cex=2, col="red") ## abline(h=1:3) ################################################### ### chunk number 19: hiddenStates ################################################### states <- c("hemDel", "normal", "ROH", "amplification") copynumberStates <- c(1, 2, 2, 3) ################################################### ### chunk number 20: hiddenStateParameters-genotypes ################################################### probs <- c(0.999, 0.7, 0.999, 0.7) names(probs) <- states ################################################### ### chunk number 21: initialProb ################################################### initialP <- (rep(1, length(states)))/length(states) ################################################### ### chunk number 22: transitionProbability ################################################### tau <- transitionProbability(chromosome=chromosome, position=position, TAUP=1e8) table(tau[, "arm"]) str(tau) ################################################### ### chunk number 23: copynumberEmission ################################################### copynumberStates <- c(1, 2, 2, 3) emission.cn <- copynumberEmission(copynumber=copynumber, states=states, mu=copynumberStates, takeLog=TRUE, verbose=TRUE) ##or copynumberStates <- matrix(c(1, 2, 2, 3), nrow(copynumber), length(states), byrow=TRUE) emission.cn2 <- copynumberEmission(copynumber=copynumber, states=states, mu=copynumberStates, takeLog=TRUE) stopifnot(identical(emission.cn, emission.cn2)) dim(emission.cn) ################################################### ### chunk number 24: qq ################################################### par(pty="s", las=1) qqnorm(log2(copynumber), cex=0.6) qqline(log2(copynumber)) abline(v=c(-1.645, 1.645)) ################################################### ### chunk number 25: epsilon ################################################### emission.cn[emission.cn[, , 2] < -10, , 2] <- -10 ################################################### ### chunk number 26: robustSds ################################################### cn.sds <- VanillaICE:::robustSds(copynumber, takeLog=TRUE) dim(cn.sds) ################################################### ### chunk number 27: multipleSamples eval=FALSE ################################################### ## CT <- matrix(sample(copynumber, 100*200, replace=TRUE), 100, 200) ## sds <- VanillaICE:::robustSds(CT, takeLog=TRUE) ################################################### ### chunk number 28: genotypeEmission ################################################### emission.gt <- genotypeEmission(genotypes=genotypes, states=states, probHomCall=probs) ################################################### ### chunk number 29: emission ################################################### emission <- emission.gt + emission.cn ################################################### ### chunk number 30: fit ################################################### fit <- viterbi(initialStateProbs=log(initialP), emission=emission, tau=tau[, "transitionPr"], arm=tau[, "arm"]) table(fit) breaks(x=fit, states=states, position=tau[, "position"], chromosome=tau[, "chromosome"], sampleNames=colnames(copynumber)) ################################################### ### chunk number 31: oligoSnpSetCrlmm ################################################### ##load("../../../tmp/VanillaICE/data/chromosome1.RData") ##cnConf <- cnConfidence(chromosome1) copynumber <- locusLevelData[["copynumber"]]/100 crlmmConfidence <- locusLevelData[["crlmmConfidence"]] genotypes <- locusLevelData[["genotypes"]] fd <- locusLevelData[["annotation"]] featuredata <- new("AnnotatedDataFrame", data=data.frame(fd), varMetadata=data.frame(labelDescription=colnames(fd))) (locusset2 <- new("oligoSnpSet", copyNumber=copynumber, calls=genotypes, callsConfidence=crlmmConfidence, phenoData=annotatedDataFrameFrom(copynumber, byrow=FALSE), featureData=featuredata, annotation="genomewidesnp6")) locusset2 <- locusset2[order(chromosome(locusset2), position(locusset2)), ] ################################################### ### chunk number 32: ice ################################################### fit3 <- hmm(object=locusset2, states=c("hemDel", "normal", "ROH", "amplification"), mu=c(1, 2, 2, 3), takeLog=TRUE, returnSegments=FALSE, ice=TRUE) brks <- hmm(object=locusset2, states=c("hemDel", "normal", "ROH", "amplification"), probs=c(0.05, 0.99, 0.7, 0.999), mu=c(1, 2, 2, 3), takeLog=TRUE, ice=TRUE) ################################################### ### chunk number 33: iceFig ################################################### i <- which(position(locusset2) >= 45*1e6 & position(locusset2) <= 55*1e6 & chromosome(locusset2) == 1) gp <- plot(locusset2[i, ]) gp$pch <- 21; gp$col <- "black"; gp$cex <- 0.9; gp$ylim <- c(0.9, 6) gp$add.cytoband <- FALSE show(gp) fit3[fit3 == 3] <- 2 ##So LOH line will not plot at 3 fit3[fit3 == 4] <- 3 lines(position(locusset2)[i], fit3[i, 1]) rohRegion <- as.numeric(as.matrix(brks[brks$state == "ROH", ][c("start", "end")])) abline(v=rohRegion, col="royalblue", lty=2) legend("topright", lty=c(1, 2), col=c("black", "royalblue"), legend=c("CN state", "ROH boundary"), bty="o", bg="white") ################################################### ### chunk number 34: vanillaPlot eval=FALSE ################################################### ## gp <- plot(snpset[chromosome(snpset) == 1, ]) ## lines(position(snpset)[chromosome(snpset)==1], viterbiResults[chromosome(snpset) == 1, ]) ## ##show copy-neutral LOH by vertical lines ## gp$abline.v <- TRUE ##plots vertical dashed lines ## allParameters <- unlist(snpPar(gp)) ## gp$col.predict[3] <- "white" ## gp$hmm.ycoords <- c(0.7,0.9) ## show(gp) ################################################### ### chunk number 35: hiddenStates eval=FALSE ################################################### ## states <- c("homozygousDeletion", "hemizygousDeletion", ## "normal", "LOH", "3copyAmp", "4copyAmp") ## mu <- c(0.05, 1, 2, 2, 3, 4) ## probs <- c(0.99, 0.999, 0.99, 0.999, 0.99, 0.99) ## probMissing <- c(0.95, rep(0.5, 5)) ################################################### ### chunk number 36: ################################################### toLatex(sessionInfo())