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


###################################################
### chunk number 2: data
###################################################
library(VanillaICE)
data(locusLevelData)


###################################################
### chunk number 3: integerCopynumber
###################################################
par(las=1)
plot(locusLevelData[["copynumber"]][, 1]/100, pch=".", ylab="copy number", log="y")
abline(h=1:3, col="grey70")


###################################################
### chunk number 4: createLocusSet
###################################################
oligoSet <- new("oligoSnpSet",
		copyNumber=locusLevelData[["copynumber"]]/100,
		call=locusLevelData[["genotypes"]],
		callProbability=locusLevelData[["crlmmConfidence"]],
		annotation=locusLevelData[["platform"]])


###################################################
### chunk number 5:  eval=FALSE
###################################################
## sds <- robustSds(log2(locusLevelData[["copynumber"]]/100))


###################################################
### chunk number 6: logscale
###################################################
copyNumber(oligoSet) <- log2(copyNumber(oligoSet))
oligoSet <- oligoSet[order(chromosome(oligoSet), position(oligoSet)), ]
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
###################################################
dim(hmmOpts$log.emission)


###################################################
### chunk number 8: fit_van
###################################################
fit.van <- hmm(oligoSet, hmmOpts)


###################################################
### chunk number 9: fig2
###################################################
library(RColorBrewer)
cols <- brewer.pal(5, "YlOrBr")[2:5]
chr1 <- oligoSet[chromosome(oligoSet)==1,]
fit.chr1 <- fit.van[space(fit.van)=="chr1", ]
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
###################################################
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
###################################################
fit.chr1 <- fit.ice[space(fit.ice)=="chr1", ]
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
###################################################
cnSet <- new("CopyNumberSet",
	     copyNumber=log2(locusLevelData[["copynumber"]]/100),
	     annotation=locusLevelData[["platform"]])
cnSet <- cnSet[order(chromosome(cnSet), position(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
###################################################
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)), ]
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
###################################################
fit.chr1 <- fit.gt[space(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")
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
###################################################
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[[3]], nrow(cnSet),1)
	copyNumber(cnSet)[50]
}


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