### 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"]],
		genome="hg19")
oligoSet <- oligoSet[!is.na(chromosome(oligoSet)), ]
oligoSet <- oligoSet[chromosome(oligoSet) < 3, ]
oligoSet <- chromosomePositionOrder(oligoSet)


###################################################
### code chunk number 4: fit_van
###################################################
fit.van <- hmm(oligoSet)
fit.van <- fit.van[[1]]


###################################################
### 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[chromosome(fit.van) == "chr1", ]
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 seq_len(length(fit.chr1))){
	polygon(x=c(sts[i], ends[i], ends[i], sts[i]),
		y=y, col=cols[elementMetadata(fit.chr1)$state[i]],
		border=cols[elementMetadata(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
###################################################
frange <- makeFeatureGRanges(oligoSet)
markersInterval2 <- subsetByOverlaps(frange, fit.van[2, ])
markersInterval2


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


###################################################
### code chunk number 9: multipanelSingleLocus (eval = FALSE)
###################################################
## xyplot2(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) & numberProbes(fit.van) > 5, ]
##elementMetadata(ranges.altered)$sampleId <- sampleNames(oligoSet)
xy.example <- xyplot2(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 <- fit.ice[[1]]
annotation(oligoSet) <- ann


###################################################
### code chunk number 14: fig3 (eval = FALSE)
###################################################
## fit.chr1 <- fit.ice[chromosome(fit.ice)=="chr1", ]
## widths <- width(fit.chr1)
## fit.chr1 <- fit.chr1[order(widths,decreasing=TRUE),]
## par(las=1)
## plot(position(chr1)/1e6, copyNumber(chr1)/100, pch=".", ylab="log2 copy number", xlab="physical position", cex=2, col="royalblue")
## points(position(chr1)[isHet]/1e6, 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 seq_len(length(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, "SnpSet2")
fit.gt <- hmm(snpSet, prGtHom=c(0.7, 0.99), normalIndex=1L, S=2L)
fit.gt <- fit.gt[[1]]


###################################################
### code chunk number 16: fig5
###################################################
fit.chr1 <- fit.gt[chromosome(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")
cols <- cols[c(1, length(cols))]
for(i in seq_len(length(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("bottomleft", fill=cols, legend=rev(c("region of homozygosity", "Normal")), bty="n")


###################################################
### code chunk number 17: readdata
###################################################
bset <- readRDS(file.path(system.file("extdata", package="VanillaICE"), "bset.rds"))


###################################################
### code chunk number 18: arm
###################################################
library(oligoClasses)
arm <- getArm(bset)
table(arm)


###################################################
### code chunk number 19: fit1p
###################################################
b1p <- bset[arm=="chr1p", ]
fit1p <- hmm(b1p, verbose=TRUE)


###################################################
### code chunk number 20: fig1p
###################################################
par(mfrow=c(2,1), las=1)
x <- position(b1p)/1e6
plot(x, lrr(b1p)/100, pch=".", col="grey")
abline(h=-0.07547)
plot(x, baf(b1p)/1000, pch=".", col="grey")
abline(h=0.542)


###################################################
### code chunk number 21: fit1q
###################################################
b1q <- bset[arm=="chr1q", ]
fit1q <- hmm(b1q, verbose=TRUE)


###################################################
### code chunk number 22: fig1q
###################################################
par(mfrow=c(2,1), las=1)
x <- position(b1q)/1e6
plot(x, lrr(b1q)/100, pch=".", col="grey")
abline(h=c(0, .359778))
plot(x, baf(b1q)/1000, pch=".", col="grey")
abline(h=c(0.235, 0.790))


###################################################
### code chunk number 23: VanillaICE.Rnw:463-464
###################################################
toLatex(sessionInfo())