PICS.R
In PICS: Probabilistic inference of ChIP-seq

## ----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")