## ----m6aPeak, echo=FALSE, out.height = "80%", out.width = "80%", include=TRUE,fig.cap="Peaks from MeRIP-seq in a mouse brain study"----
knitr::include_graphics("m6APeak.png")

## ---- eval = FALSE, warning=FALSE, message=FALSE------------------------------
#  install.packages("devtools") # if you have not installed "devtools" package
#  library(devtools)
#  install_github("https://github.com/ZhenxingGuo0015/TRESS",
#                 build_vignettes = TRUE)

## ----eval=FALSE,warning=FALSE,message=FALSE-----------------------------------
#  library(TRESS)
#  vignette("TRESS")

## ---- eval= FALSE-------------------------------------------------------------
#  install_github("https://github.com/ZhenxingGuo0015/datasetTRES")

## ---- eval=FALSE, message=FALSE, warning=FALSE--------------------------------
#  ## Directly use "makeTxDbFromUCSC" function to create one
#  library(GenomicFeatures)
#  txdb = makeTxDbFromUCSC("mm9", "knownGene")
#  # saveDb(txdb, file = paste0("YourPATH", "/", "YourGenome.sqlite")
#  
#  ## or load a TxDb annotation package like
#  # library(TxDb.Mmusculus.UCSC.mm9.knownGene)
#  # txdb <- TxDb.Mmusculus.UCSC.mm9.knownGene
#  # saveDb(txdb, file = paste0("YourPATH", "/", "mm9_knownGene.sqlite")

## ---- eval=FALSE, message=FALSE, warning=FALSE--------------------------------
#  library(GenomicFeatures)
#  txdb = makeTxDbFromGFF("directory/to/your/xxx.gtf", format = "gtf")
#  # saveDb(txdb, file = paste0("YourPATH", "/", "YourGenome.sqlite")

## ----loadData,warning=FALSE,message=FALSE-------------------------------------
library(TRESS)
data("Basal")
Bins <- Basal$Bins$Bins
BinCounts <- Basal$Bins$Counts
dim(BinCounts)
Candidates <- Basal$Candidates$Regions
CandidatesCounts <- Basal$Candidates$Counts
sf <- Basal$Bins$sf

## ----BasalBins----------------------------------------------------------------
head(Bins, 3) 
head(Candidates, 3)

## ----BasalCount---------------------------------------------------------------
dim(BinCounts)
head(BinCounts, 3)
dim(CandidatesCounts)
head(CandidatesCounts,3)

## ----Basalsf------------------------------------------------------------------
sf

## ---- eval= FALSE,warning=FALSE, message=FALSE--------------------------------
#  ## Directly take BAM files in "datasetTRES" available on github
#  library(datasetTRES)
#  Input.file = c("cb_input_rep1_chr19.bam", "cb_input_rep2_chr19.bam")
#  IP.file = c("cb_ip_rep1_chr19.bam", "cb_ip_rep2_chr19.bam")
#  BamDir = file.path(system.file(package = "datasetTRES"), "extdata/")
#  annoDir = file.path(system.file(package = "datasetTRES"),
#                      "extdata/mm9_chr19_knownGene.sqlite")
#  OutDir = "/directory/to/output"
#  TRESS_peak(IP.file = IP.file,
#             Input.file = Input.file,
#             Path_To_AnnoSqlite = annoDir,
#             InputDir = BamDir,
#             OutputDir = OutDir, # specify a directory for output
#             experiment_name = "examplebyBam", # name your output
#             filetype = "bam")

## ---- eval= TRUE--------------------------------------------------------------
### example peaks
peaks = read.table(file.path(system.file(package = "TRESS"),
                           "extdata/examplebyBam_peaks.xls"),
                 sep = "\t", header = TRUE)
head(peaks[, -c(5, 14, 15)], 3)

## ---- eval=FALSE--------------------------------------------------------------
#  ## or, take BAM files from your path
#  Input.file = c("input_rep1.bam", "input_rep2.bam")
#  IP.file = c("ip_rep1.bam", "ip_rep2.bam")
#  BamDir = "/directory/to/BAMfile"
#  annoDir = "/path/to/xxx.sqlite"
#  OutDir = "/directory/to/output"
#  TRESS_peak(IP.file = IP.file,
#             Input.file = Input.file,
#             Path_To_AnnoSqlite = annoDir,
#             InputDir = BamDir,
#             OutputDir = OutDir,
#             experiment_name = "xxx",
#             filetype = "bam")
#  peaks = read.table(paste0(OutDir, "/",
#                            "xxx_peaks.xls"),
#                     sep = "\t", header = TRUE)
#  head(peaks, 3)

## ----eval=TRUE, message=FALSE, warning=FALSE----------------------------------
## load in first example dataset
data("Basal") 
Candidates = CallCandidates(
    Counts = Basal$Bins$Counts,
    bins = Basal$Bins$Bins
    )

## ---- eval=TRUE, message= TRUE, warning= FALSE--------------------------------
data("Basal") ### load candidate regions
Basal$Candidates$sf = Basal$Bins$sf
peaks1 = CallPeaks.multiRep(
  Candidates = Basal$Candidates,
  mu.cutoff = 0.5
  )
head(peaks1, 3)

## ---- eval=TRUE, message= TRUE, warning= FALSE--------------------------------
### use different threshold to filter peaks
peaks2 = CallPeaks.multiRep(
  Candidates = Basal$Candidates,
  mu.cutoff = 0.5, 
  fdr.cutoff = 0.01, 
  lfc.cutoff = 1.6  
  )

## ---- eval=TRUE, message= TRUE, warning= FALSE--------------------------------
### use different threshold to filter peaks
peaks3 = CallPeaks.multiRep(
  Candidates = Basal$Candidates,
  mu.cutoff = 0.5, 
  WhichThreshold = "fdr"
  )

## ---- eval = TRUE, message= FALSE, warning= FALSE-----------------------------
# A toy example
data("Basal")
bincounts = Basal$Bins$Counts[, 1:2]
sf0 = Basal$Bins$sf[1:2]
bins = Basal$Bins$Bins
peaks = CallPeaks.oneRep(Counts = bincounts, 
                         sf = sf0, bins = bins)
head(peaks, 3)

## ---- eval=FALSE, message= FALSE, warning= FALSE------------------------------
#  ShowOnePeak(onePeak, allBins, binCounts, ext = 500, ylim = c(0,1))

## ---- eval=TRUE, message= FALSE, warning= FALSE-------------------------------
peaks = read.table(file.path(system.file(package = "TRESS"),
                             "extdata/examplebyBam_peaks.xls"),
                   sep = "\t", header = TRUE) 
load(file.path(system.file(package = "TRESS"),
               "extdata/examplebyBam.rda"))
allBins = as.data.frame(bins$bins)
colnames(allBins)[1] = "chr"
allBins$strand = binStrand

head(peaks, 1)
ShowOnePeak(onePeak = peaks[1,], allBins = allBins, 
            binCounts = allCounts)

## ----sessionInfo, echo=FALSE--------------------------------------------------
sessionInfo()