###################################################
### chunk number 1: setup
###################################################
options(width=80)


###################################################
### chunk number 2: createLocusSet
###################################################
library(VanillaICE)
library(Biobase)
library(oligoClasses)
data(locusLevelData)
copynumber <- locusLevelData[["copynumber"]]/100
chromosome <- locusLevelData[["annotation"]][, "chromosome"]
position <- locusLevelData[["annotation"]][, "position"]
names(position) <- names(chromosome) <- rownames(locusLevelData[["annotation"]])
locusAnnotation <- data.frame(list(chromosome=chromosome, position=position), row.names=names(chromosome))
featuredata <- new("AnnotatedDataFrame", data=locusAnnotation, varMetadata=data.frame(labelDescription=colnames(locusAnnotation)))
locusset <- new("oligoSnpSet",
		copyNumber=copynumber,
		calls=locusLevelData[["genotypes"]],
		callsConfidence=locusLevelData[["crlmmConfidence"]],
		phenoData=annotatedDataFrameFrom(copynumber, byrow=FALSE),		
		featureData=featuredata,
		annotation=locusLevelData[["platform"]])
stopifnot(all(c("chromosome", "position") %in% varLabels(featuredata)))
locusset <- locusset[order(chromosome(locusset), position(locusset)), ]


###################################################
### chunk number 3: joint.states
###################################################
joint.states <- c("hemDel", "normal", "ROH", "amplification")
copynumber.states <- log2(c(1, 2, 2, 3))
prob.genotype.is.homozygous <- c(0.999, 0.7, 0.999, 0.7)
names(prob.genotype.is.homozygous) <- joint.states
initialP <- (rep(1, length(joint.states)))/length(joint.states)
tau <- transitionProbability(chromosome=chromosome,
			     position=position,
			     TAUP=1e8)
log.sds <- VanillaICE:::robustSds(copyNumber(locusset))


###################################################
### chunk number 4: emissionProbabilities
###################################################
emission.cn <- copynumberEmission(copynumber=log2(copyNumber(locusset)),
				  states=joint.states,
				  mu=copynumber.states,
				  sds=log.sds,
				  takeLog=FALSE,
				  verbose=TRUE)
emission.cn[1:5, , ]
emission.gt <- genotypeEmission(genotypes=calls(locusset),
				states=joint.states,
				probHomCall=prob.genotype.is.homozygous)


###################################################
### chunk number 5: jointEmission
###################################################
emission.joint <- emission.gt + emission.cn


###################################################
### chunk number 6: fit1
###################################################
fit1 <- viterbi(initialStateProbs=log(initialP),
		emission=emission.joint,
		tau=tau[, "transitionPr"],
		arm=tau[, "arm"], 
		normalIndex=2)
(brks <- breaks(x=fit1, states=joint.states, position=tau[, "position"],
		chromosome=tau[, "chromosome"]))
rohBoundary <- as.integer(as.matrix(brks[brks$state == "ROH", ][c("start", "end")]))


###################################################
### chunk number 7: plotPredictions
###################################################
require(SNPchip)
gp <- plot(locusset[chromosome(locusset) == 1, ])
cns <- fit1
cns[cns == 3] <- 2
##show copy-neutral LOH by vertical lines
gp$abline.v <- TRUE ##plots vertical dashed lines
gp$col.predict[3] <- "white"
gp$hmm.ycoords <- c(0.7,0.9)
show(gp)
lines(position(locusset)[chromosome(locusset)==1], cns[chromosome(locusset) == 1, ])
abline(v=rohBoundary, col="royalblue", lty=2)
legend("topright", lty=c(2, 1), col=c("royalblue", "black"), legend=c("ROH boundary",
							     "copy number"), bty="n")


###################################################
### chunk number 8: supportedPlatforms
###################################################
path <- system.file("extdata", package="VanillaICE")
supportedPlatforms <- list.files(path)[grep("Conf", list.files(path))]
as.character(sapply(supportedPlatforms, function(x) strsplit(x, "Conf")[[1]][[1]]))


###################################################
### chunk number 9: hetIndices
###################################################
hetIndices <- which(calls(locusset) == 2 & fit1 == 3)
confs <- round(1-exp(-callsConfidence(locusset)[hetIndices]/1000), 5)


###################################################
### chunk number 10: emissionCompare
###################################################
emission.crlmm <- genotypeEmissionCrlmm(genotypes=calls(locusset),
					conf=callsConfidence(locusset),
					annotation=annotation(locusset),
					pHetCalledHom=0.001,
					pHetCalledHet=0.995,
					pHomInNormal=0.8,
					pHomInRoh=0.999)
emission.crlmm[hetIndices, , ]


###################################################
### chunk number 11: updateEmission
###################################################
joint.emission2 <- emission.cn
joint.emission2[, , c("normal", "amplification")] <- joint.emission2[, , c("normal", "amplification")] + emission.crlmm[, , "norm"]
joint.emission2[, , c("hemDel", "ROH")] <- joint.emission2[, , c("hemDel", "ROH")] + emission.crlmm[, ,"ROH"]
fit2 <- viterbi(initialStateProbs=log(initialP),
		emission=joint.emission2,
		tau=tau[, "transitionPr"],
		arm=tau[, "arm"], 
		normalIndex=2)
table(fit2)
breaks(x=fit2, states=joint.states, position=tau[, "position"],
       chromosome=tau[, "chromosome"])


###################################################
### chunk number 12: copyNumberOnly
###################################################
states <- 0:5
mus <- log2(c(0.05, 1:5))
emission.cn <- copynumberEmission(copynumber=log2(copyNumber(locusset)),
				  states=states,
				  mu=mus,
				  sds=log.sds,
				  takeLog=FALSE,
				  verbose=TRUE)
initialP <- log(rep(1/length(states), length(states)))
fit3 <- viterbi(initialStateProbs=initialP,
		emission=emission.cn,
		tau=tau[, "transitionPr"],
		arm=tau[, "arm"], 
		normalIndex=3)
breaks(x=fit3, states=states, position=tau[, "position"],
       chromosome=tau[, "chromosome"])


###################################################
### chunk number 13: genotypesOnly
###################################################
states <- c("normal", "ROH")
prob.genotype.is.homozygous <- c(0.7, 0.999)
##without confidence scores
emission.gt <- genotypeEmission(genotypes=calls(locusset),
				states=states,
				probHomCall=prob.genotype.is.homozygous)
initialP <- log(rep(1/length(states), length(states)))
fit4 <- viterbi(initialStateProbs=initialP,
		emission=emission.gt,
		tau=tau[, "transitionPr"],
		arm=tau[, "arm"], 
		normalIndex=1)
breaks(x=fit4, states=states, position=tau[, "position"],
       chromosome=tau[, "chromosome"])


###################################################
### chunk number 14: qq
###################################################
par(pty="s", las=1)
qqnorm(log2(copyNumber(locusset)), cex=0.6)
qqline(log2(copyNumber(locusset)))
abline(v=c(-1.645, 1.645))


###################################################
### chunk number 15: epsilon eval=FALSE
###################################################
## emission.cn[emission.cn[, , "normal"] < -10, , "normal"] <- -10


###################################################
### chunk number 16: robustSds
###################################################
cn.sds <- VanillaICE:::robustSds(copyNumber(locusset), takeLog=TRUE)
dim(cn.sds)


###################################################
### chunk number 17: multipleSamples
###################################################
CT <- matrix(sample(copyNumber(locusset), 100*200, replace=TRUE), 100, 200)
sds <- VanillaICE:::robustSds(CT, takeLog=TRUE)


###################################################
### chunk number 18: fit eval=FALSE
###################################################
## fit <- viterbi(initialStateProbs=log(initialP),
## 	       emission=emission,
## 	       tau=tau[, "transitionPr"],
## 	       arm=tau[, "arm"], 
## 	       normalIndex=2,
## 	       normal2altered=0.005,
## 	       altered2normal=0.5,
## 	       altered2altered=0.0025)


###################################################
### chunk number 19: 
###################################################
toLatex(sessionInfo())