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

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


###################################################
### code chunk number 2: example_data
###################################################
library(VanillaICE)
extdir <- system.file("extdata", package="VanillaICE", mustWork=TRUE)
files <- list.files(extdir, pattern="FinalReport")


###################################################
### code chunk number 3: annotation_file
###################################################
list.files(extdir, pattern="SNP_info")


###################################################
### code chunk number 4: FeatureAnnotation
###################################################
library(data.table)
features <- suppressWarnings(fread(file.path(extdir, "SNP_info.csv")))
fgr <- GRanges(paste0("chr", features$Chr), IRanges(features$Position, width=1),
               isSnp=features[["Intensity Only"]]==0)
fgr <- SnpGRanges(fgr)
names(fgr) <- features[["Name"]]


###################################################
### code chunk number 5: seqinfo
###################################################
library(BSgenome.Hsapiens.UCSC.hg18)
sl <- seqlevels(BSgenome.Hsapiens.UCSC.hg18)
seqlevels(fgr) <- sl[sl %in% seqlevels(fgr)]
seqinfo(fgr) <- seqinfo(BSgenome.Hsapiens.UCSC.hg18)[seqlevels(fgr),]
fgr <- sort(fgr)


###################################################
### code chunk number 6: sourceFiles
###################################################
files <- list.files(extdir, full.names=TRUE, recursive=TRUE, pattern="FinalReport")


###################################################
### code chunk number 7: ArrayViews
###################################################
##
## Where to keep parsed files
##
parsedDir <- "ParsedFiles"
if(!file.exists(parsedDir)) dir.create(parsedDir)
views <- ArrayViews(rowData=fgr, sourcePaths=files, parsedPath=parsedDir)
show(views)


###################################################
### code chunk number 8: fread
###################################################
## read the first file
dat <- fread(files[1])
head(dat,n=3)


###################################################
### code chunk number 9: select_columns
###################################################
## information to store on the markers
select_columns <- match(c("SNP Name", "Allele1 - AB", "Allele2 - AB",
                          "Log R Ratio", "B Allele Freq"), names(dat))


###################################################
### code chunk number 10: order_of_markers
###################################################
index_genome <- match(names(fgr), dat[["SNP Name"]])


###################################################
### code chunk number 11: scan_params
###################################################
scan_params <- CopyNumScanParams(index_genome=index_genome, select=select_columns,
                                 cnvar="Log R Ratio",
                                 bafvar="B Allele Freq",
                                 gtvar=c("Allele1 - AB", "Allele2 - AB"))


###################################################
### code chunk number 12: parseSourceFile
###################################################
parsedPath(views)
parseSourceFile(views[, 1], scan_params)


###################################################
### code chunk number 13: applyParseSourceFile
###################################################
invisible(sapply(views, parseSourceFile, param=scan_params))


###################################################
### code chunk number 14: list_parsed_files
###################################################
head(list.files(parsedPath(views)), n=3)


###################################################
### code chunk number 15: Views
###################################################
lrr(views)[1:2, 2:4]
## or
lrr(views[1:2, 2:4])
## B allele frequencies
baf(views[1:2, 2:4])
## potentially masked by function of the same name in crlmm
VanillaICE::genotypes(views)[1:2, 2:4]


###################################################
### code chunk number 16: SnpArrayExperiment
###################################################
snp_exp <- SnpExperiment(views[, 4:6])
show(snp_exp)


###################################################
### code chunk number 17: emission_param
###################################################
param <- EmissionParam()


###################################################
### code chunk number 18: temper
###################################################
param <- EmissionParam(temper=0.5)
show(param)


###################################################
### code chunk number 19: hmm
###################################################
if(!exists("fit"))
  fit <- hmm2(snp_exp, param)
show(fit)


###################################################
### code chunk number 20: HMMList_names
###################################################
names(fit)


###################################################
### code chunk number 21: HMMList_subset
###################################################
show(fit[[2]])


###################################################
### code chunk number 22: HMMList_unlist
###################################################
head(unlist(fit), n=3)


###################################################
### code chunk number 23: HMMList_grangeslist
###################################################
grl <- split(unlist(fit), unlist(fit)$id)


###################################################
### code chunk number 24: filter_param
###################################################
filter_param <- FilterParam()
show(filter_param)


###################################################
### code chunk number 25: cnvFilter
###################################################
cnvFilter(fit, filter_param)


###################################################
### code chunk number 26: cnvFilter_altered
###################################################
select_cnv <- FilterParam(state=c("1", "2", "5", "6"), seqnames="chr22")
cnvs <- cnvFilter(fit, select_cnv)
cnvs


###################################################
### code chunk number 27: xyplotList
###################################################
trellis_param <- HmmTrellisParam()
cnvList <- split(cnvs, cnvs$id)
figList <- xyplotList(cnvList, snp_exp)
names(figList)


###################################################
### code chunk number 28: trellis
###################################################
class(figList[["FinalReport6841.csv"]][[1]])


###################################################
### code chunk number 29: viewports
###################################################
fig1 <- figList[["FinalReport6841.csv"]][[1]]


###################################################
### code chunk number 30: fig1
###################################################
vps <- viewports()
xygrid(fig1, vps, cnvList[[1]][1])


###################################################
### code chunk number 31: fig2
###################################################
cnvs_sample2 <- cnvList[[2]]
cnvs_sample2
xygrid(figList[[2]][[1]], vps, cnvs_sample2[1])


###################################################
### code chunk number 32: reduce
###################################################
cnvs_sample2r <- reduce(cnvs_sample2, min.gapwidth=500e3)
fig2 <- xyplotList(cnvs_sample2r, snp_exp)


###################################################
### code chunk number 33: fig_reduced
###################################################
invisible(print(fig2[[1]]))


###################################################
### code chunk number 34: parallelEnvironment (eval = FALSE)
###################################################
## library(foreach)
## library(snow)
## library(doSNOW)
## cl <- makeCluster(2, type="SOCK")
## registerDoSNOW(cl)


###################################################
### code chunk number 35: fitParallel (eval = FALSE)
###################################################
## results <- hmm2(snp_exp)


###################################################
### code chunk number 36: stopcl (eval = FALSE)
###################################################
## stopCluster(cl)


###################################################
### code chunk number 37: VanillaICE.Rnw:556-557
###################################################
toLatex(sessionInfo())