### R code from vignette source 'vignettes/VanillaICE/inst/doc/VanillaICE.Rnw'

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


###################################################
### code chunk number 2: data
###################################################
library(oligoClasses)
library(VanillaICE)
library2(SNPchip)
library2(Biobase)
library2(IRanges)
data(locusLevelData)


###################################################
### code chunk number 3: createLocusSet
###################################################
cn <- log2(locusLevelData[["copynumber"]]/100)
oligoSet <- new("oligoSnpSet",
		copyNumber=integerMatrix(cn),
		call=locusLevelData[["genotypes"]],
		callProbability=locusLevelData[["crlmmConfidence"]],
		annotation=locusLevelData[["platform"]])
oligoSet <- oligoSet[!is.na(chromosome(oligoSet)), ]
oligoSet <- oligoSet[chromosome(oligoSet) < 3, ]


###################################################
### code chunk number 4: fit_van
###################################################
oligoSet <- oligoSet[chromosome(oligoSet) <= 2, ]
fit.van <- hmm(oligoSet)


###################################################
### code chunk number 5: fig2
###################################################
if(require(RColorBrewer)){
	cols <- brewer.pal(6, "YlOrBr")
} else cols <- rainbow(n=6)
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)/100, pch=".", cex=2, col="royalblue",
     ylab="log2 copy number")
points(position(chr1)[isHet], copyNumber(chr1)[isHet]/100, 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=c("hom-del", "hem-del", "N", "N-noHets", "sDup", "dDup"), bty="n")


###################################################
### code chunk number 6: range2
###################################################
fit.van[2, ]


###################################################
### code chunk number 7: findOverlaps
###################################################
mm <- as.matrix(findOverlaps(fit.van, oligoSet))
markersInRange <- featureNames(oligoSet)[mm[,1]==2]


###################################################
### code chunk number 8: multipanelSingleLocus (eval = FALSE)
###################################################
## xyplot(cn ~ x | id, oligoSet, range=RangedDataObject, ylim=c(-0.5,4), panel=xypanel)


###################################################
### code chunk number 9: multipanelSingleLocus (eval = FALSE)
###################################################
## xyplot(cn ~ x | id, oligoSet, range=RangedDataObject, ylim=c(-0.5,4), scales=list(x="free"),
##        panel=xypanel)


###################################################
### code chunk number 10: xyplotExample
###################################################
ranges.altered <- fit.van[!state(fit.van) %in% c(3,4) & coverage2(fit.van) > 5, ]
xy.example <- xyplot(cn ~ x | range, oligoSet, range=ranges.altered, frame=2e6,
		     scales=list(x="free"), ylim=c(-1,3),
		     panel=xypanel, cex.pch=0.4, pch=21, border="grey",
		     ylab=expression(log[2]("copy number")))


###################################################
### code chunk number 11: xyplotPrint
###################################################
pdf("VanillaICE-xyplot.pdf", width=8, height=7)
print(xy.example)
dev.off()


###################################################
### code chunk number 12: ice
###################################################
VanillaICE:::icePlatforms()


###################################################
### code chunk number 13: iceIllustration
###################################################
ann <- annotation(oligoSet)
annotation(oligoSet) <- "genomewidesnp6"
fit.ice <- hmm(oligoSet, ICE=TRUE)
fit.ice
annotation(oligoSet) <- ann


###################################################
### code chunk number 14: fig3 (eval = FALSE)
###################################################
## fit.chr1 <- fit.ice[chromosome(fit.ice)==1, ]
## widths <- width(fit.chr1)
## fit.chr1 <- fit.chr1[order(widths,decreasing=TRUE),]
## par(las=1)
## plot(position(chr1), copyNumber(chr1)/100, pch=".", ylab="log2 copy number", xlab="physical position", cex=2, col="royalblue")
## points(position(chr1)[isHet], copyNumber(chr1)[isHet]/100, 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[state(fit.chr1)[i]],
## 		border=cols[state(fit.chr1)[i]])
## }
## legend("topleft", fill=cols, legend=c("hom-del","hem-del", "N", "N/no hets", "s-dup", "d-dup"), bty="n")


###################################################
### code chunk number 15: genotypesOnly
###################################################
snpSet <- as(oligoSet, "SnpSet")
fit.gt <- hmm(snpSet, prGtHom=c(0.7, 0.99), normalIndex=1L, S=2L)


###################################################
### code chunk number 16: fig5
###################################################
fit.chr1 <- fit.gt[chromosome(fit.gt)==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=c("region of homozygosity", "Normal"), bty="n")


###################################################
### code chunk number 17: VanillaICE.Rnw:446-447
###################################################
toLatex(sessionInfo())