## ----opts, echo = FALSE, message = FALSE, warning = FALSE--------------------- library(knitr) opts_chunk$set(echo = TRUE, fig.align = "center", message = FALSE, warning = F) ## ----lib---------------------------------------------------------------------- library(PICS) ## ----reading_data------------------------------------------------------------- path <- system.file("extdata", package = "PICS") dataIP <- read.table(file.path(path, "Treatment_tags_chr21_sort.bed"), header=TRUE, colClasses=c("factor", "integer", "integer", "factor")) dataIP <- as(dataIP, "GRanges") dataCont <- read.table(file.path(path, "Input_tags_chr21_sort.bed"), header = TRUE, colClasses = c("factor", "integer", "integer", "factor")) dataCont <- as(dataCont, "GRanges") ## ----mappability_prof--------------------------------------------------------- map <- read.table(file.path(path, "mapProfileShort"), header=TRUE, colClasses=c("factor", "integer", "integer", "NULL")) map <- as(map, "GRanges") ## ----genome_segmentation------------------------------------------------------ seg <- segmentPICS(dataIP, dataC = dataCont, map = map, minReads = 1) ## ----cluster_init------------------------------------------------------------- library(parallel) ## ----PICS--------------------------------------------------------------------- pics <- PICS(seg, nCores = 2) ## ----FDR---------------------------------------------------------------------- segC <- segmentPICS(dataCont, dataC = dataIP, map = map, minReads = 1) picsC <- PICS(segC) fdr <- picsFDR(pics, picsC, filter = list(delta = c(50, Inf), se = c(0, 50), sigmaSqF = c(0, 22500), sigmaSqR = c(0, 22500))) ## ----plotFDR------------------------------------------------------------------ plot(pics, picsC, xlim = c(2, 8), ylim = c(0, .2), filter = list(delta = c(50,300), se = c(0, 50), sigmaSqF = c(0, 22500), sigmaSqR = c(0, 22500)), type = "l") ## ----plotFDR2----------------------------------------------------------------- plot(fdr[, c(3, 1)]) ## ----------------------------------------------------------------------------- myFilter = list(delta = c(50, 300), se = c(0, 50), sigmaSqF = c(0, 22500), sigmaSqR = c(0, 22500)) rdBed <- makeGRangesOutput(pics, type="bed", filter = c(myFilter, list(score = c(1, Inf)))) ## ----bedfile, eval = F-------------------------------------------------------- # library(rtracklayer) # export(rdBed, "myfile.bed") ## ----wigfile, eval = F-------------------------------------------------------- # rdWig <- makeGRangesOutput(pics, type="wig", filter=c(myFilter, list(score=c(1,Inf)))) # export(rdWig, "myfile.wig")