###################################################
### 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())


