### R code from vignette source 'VanillaICE.Rnw' ################################################### ### code chunk number 1: setup ################################################### options(width=70) ################################################### ### code chunk number 2: example_data ################################################### library(VanillaICE) library(foreach) registerDoSEQ() 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 ################################################### require(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 ################################################### ## ## A directory to keep parsed files ## parsedDir <- tempdir() views <- ArrayViews(rowRanges=fgr, sourcePaths=files, parsedPath=parsedDir) lrrFile(views) <- file.path(parsedDir, basename(fileName(views, "lrr"))) views@bafFiles <- file.path(parsedDir, basename(fileName(views, "baf"))) views@gtFiles <- file.path(parsedDir, basename(fileName(views, "gt"))) colnames(views) <- gsub(".csv", "", colnames(views)) 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: applyParseSourceFile ################################################### parseSourceFile(views, scan_params) ################################################### ### code chunk number 13: list_parsed_files ################################################### head(list.files(parsedPath(views)), n=3) ################################################### ### code chunk number 14: Views ################################################### lrr(views)[1:2, ] ## or lrr(views[1:2, ]) ## B allele frequencies baf(views[1:2, ]) ## Use :: to avoid masking by function of the same name in crlmm VanillaICE::genotypes(views)[1:2, ] ################################################### ### code chunk number 15: SnpArrayExperiment ################################################### snp_exp <- SnpExperiment(views[, 4:5]) show(snp_exp) ################################################### ### code chunk number 16: emission_param ################################################### param <- EmissionParam() ################################################### ### code chunk number 17: temper ################################################### param <- EmissionParam(temper=0.5) show(param) ################################################### ### code chunk number 18: hmm ################################################### fit <- hmm2(snp_exp, param) show(fit) length(fit) ## HMM object for the first sample fit[[1]] ################################################### ### code chunk number 19: HMMList_names ################################################### names(fit) ################################################### ### code chunk number 20: HMMList_subset ################################################### show(fit[[2]]) ################################################### ### code chunk number 21: HMMList_unlist ################################################### head(unlist(fit), n=3) ################################################### ### code chunk number 22: HMMList_grangeslist ################################################### grl <- split(unlist(fit), unlist(fit)$id) ################################################### ### code chunk number 23: filter_param ################################################### filter_param <- FilterParam() show(filter_param) ################################################### ### code chunk number 24: cnvFilter ################################################### cnvFilter(fit, filter_param) ################################################### ### code chunk number 25: cnvFilter_altered ################################################### select_cnv <- FilterParam(state=c("1", "2", "5", "6"), seqnames="chr22") cnvs <- cnvFilter(fit, select_cnv) cnvs ################################################### ### code chunk number 26: xyplotList ################################################### trellis_param <- HmmTrellisParam() cnvList <- split(cnvs, cnvs$id) figList <- xyplotList(cnvList, snp_exp, trellis_param) names(figList) ################################################### ### code chunk number 27: trellis ################################################### class(figList[["FinalReport6841.csv"]][[1]]) ################################################### ### code chunk number 28: viewports ################################################### fig1 <- figList[["FinalReport6841.csv"]][[1]] ################################################### ### code chunk number 29: fig1 ################################################### vps <- viewports() xygrid(fig1, vps, cnvList[[1]][1]) ################################################### ### code chunk number 30: fig2 ################################################### cnvs_sample2 <- cnvList[[2]] cnvs_sample2 xygrid(figList[[2]][[1]], vps, cnvs_sample2[1]) ################################################### ### code chunk number 31: reduce ################################################### cnvs_sample2r <- reduce(cnvs_sample2, min.gapwidth=500e3) fig2 <- xyplotList(cnvs_sample2r, snp_exp) ################################################### ### code chunk number 32: fig_reduced ################################################### invisible(print(fig2[[1]])) ################################################### ### code chunk number 33: parallelEnvironment (eval = FALSE) ################################################### ## library(foreach) ## library(snow) ## library(doSNOW) ## cl <- makeCluster(2, type="SOCK") ## registerDoSNOW(cl) ################################################### ### code chunk number 34: fitParallel (eval = FALSE) ################################################### ## results <- hmm2(snp_exp) ################################################### ### code chunk number 35: stopcl (eval = FALSE) ################################################### ## stopCluster(cl) ################################################### ### code chunk number 36: VanillaICE.Rnw:555-556 ################################################### toLatex(sessionInfo())