###################################################
### chunk number 1: setup
###################################################
#line 50 "VanillaICE.Rnw"
options(width=70)


###################################################
### chunk number 2: data
###################################################
#line 111 "VanillaICE.Rnw"
library(VanillaICE)
data(locusLevelData)


###################################################
### chunk number 3: integerCopynumber
###################################################
#line 121 "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 "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 "VanillaICE.Rnw"
## sds <- robustSds(log2(locusLevelData[["copynumber"]]/100))


###################################################
### chunk number 6: logscale
###################################################
#line 171 "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 "VanillaICE.Rnw"
dim(hmmOpts$log.emission)


###################################################
### chunk number 8: fit_van
###################################################
#line 196 "VanillaICE.Rnw"
fit.van <- hmm(oligoSet, hmmOpts)


###################################################
### chunk number 9: fig2
###################################################
#line 202 "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 "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 "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 "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 "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 "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 "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 "VanillaICE.Rnw"
toLatex(sessionInfo())


