## ----setup, echo=FALSE--------------------------------------------------------
knitr::opts_chunk$set(
  message = FALSE, warning = FALSE, fig.height = 8, fig.width = 10
)

## ----install, eval = FALSE----------------------------------------------------
#  if (!"BiocManager" %in% rownames(installed.packages()))
#    install.packages("BiocManager")
#  pkg <- c("tidyverse", "Rsamtools", "csaw", "BiocParallel", "rtracklayer")
#  BiocManager::install(pkg)
#  BiocManager::install("extraChIPs")

## ----wincounts----------------------------------------------------------------
library(tidyverse)
library(Rsamtools)
library(csaw)
library(BiocParallel)
library(rtracklayer)
bfl <- system.file(
  "extdata", "bam", c("ex1.bam", "ex2.bam", "input.bam"), package = "extraChIPs"
) %>% 
  BamFileList()
names(bfl) <- c("ex1", "ex2", "input")
rp <- readParam(
  pe = "none",
  dedup = TRUE,
  restrict = "chr10"
)
wincounts <- windowCounts(
  bam.files = bfl,
  spacing = 60,
  width = 180,
  ext = 200,
  filter = 1,
  param = rp
)

## ----add-totals---------------------------------------------------------------
wincounts$totals <- c(964076L, 989543L, 1172179L)

## ----update-coldata-----------------------------------------------------------
wincounts$sample <- colnames(wincounts)
wincounts$treat <- as.factor(c("ctrl", "treat", NA))
colData(wincounts)

## ----plot-densities-----------------------------------------------------------
library(extraChIPs)
plotAssayDensities(wincounts, colour = "treat", trans = "log1p")

## ----peaks--------------------------------------------------------------------
peaks <- import.bed(
  system.file("extdata", "peaks.bed.gz", package = "extraChIPs")
)
peaks <- granges(peaks)

## ----filtcounts---------------------------------------------------------------
filtcounts <- dualFilter(
  x = wincounts[, !is.na(wincounts$treat)],
  bg = wincounts[, is.na(wincounts$treat)], 
  ref = peaks,
  q = 0.8 # Better to use q = 0.5 on real data
)

## ----plotcpm------------------------------------------------------------------
plotAssayDensities(filtcounts, assay = "logCPM", colour = "treat")
plotAssayPCA(filtcounts, assay = "logCPM", colour = "treat", label = "sample")

## ----dims---------------------------------------------------------------------
dim(wincounts)
dim(filtcounts)

## ----filt-ranges--------------------------------------------------------------
rowRanges(filtcounts)
mean(rowRanges(filtcounts)$overlaps_ref)

## ----voom---------------------------------------------------------------------
v <- voomWeightsFromCPM(
  cpm = assay(filtcounts, "logCPM"), 
  lib.size = filtcounts$totals, 
  isLogCPM = TRUE
)

## ----add-vals-----------------------------------------------------------------
rowRanges(filtcounts)$logCPM <- rowMeans(assay(filtcounts,"logCPM"))
rowRanges(filtcounts)$logFC <- rowDiffs(assay(filtcounts,"logCPM"))[,1]
rowRanges(filtcounts)$PValue <- 1 - pchisq(rowRanges(filtcounts)$logFC^2, 1)

## ----add-ol-------------------------------------------------------------------
res_gr <- mergeByCol(filtcounts, col = "logCPM", pval = "PValue")
res_gr$overlaps_ref <- overlapsAny(res_gr, peaks)

## ----ex-mapping1--------------------------------------------------------------
data("ex_genes")
data("ex_prom")
mapByFeature(
  res_gr, genes = ex_genes, prom = ex_prom, cols = c("gene", "symbol")
)

## ----ex-mapping2--------------------------------------------------------------
data("ex_hic")
res_gr_mapped <- mapByFeature(
  res_gr, 
  genes = ex_genes, 
  prom = ex_prom,
  gi = ex_hic, 
  cols = c("gene", "symbol")
)
res_gr_mapped

## ----plot-pie-----------------------------------------------------------------
res_gr %>% 
  as_tibble() %>% 
  plotPie("overlaps_ref")

## ----plot-dual-pie------------------------------------------------------------
res_gr$Feature <- bestOverlap(
  res_gr, GRangesList(Promoter = ex_prom), missing = "None"
)
res_gr %>% 
  as_tibble() %>% 
  plotPie(x = "Feature", fill = "overlaps_ref") 

## ----bwfl---------------------------------------------------------------------
bwfl <- system.file(
  "extdata", "bigwig", c("ex1.bw", "ex2.bw"), package = "extraChIPs"
) %>% 
  BigWigFileList() %>% 
  setNames(c("ex1", "ex2"))

## ----get-profile--------------------------------------------------------------
pd <- getProfileData(bwfl, res_gr)
pd

## ----profile-heatmap----------------------------------------------------------
plotProfileHeatmap(pd, "profile_data") +
  scale_fill_viridis_c() +
  labs(fill = "CPM") +
  theme_bw()

## ----centred-heatmap----------------------------------------------------------
pd <- getProfileData(bwfl, colToRanges(res_gr, "keyval_range"))
plotProfileHeatmap(pd, "profile_data")  +
  scale_fill_viridis_c() +
  labs(fill = "CPM") +
  theme_bw()

## ----facet-heatmap------------------------------------------------------------
plotProfileHeatmap(
  pd, "profile_data", facetY = "overlaps_ref", linetype = "overlaps_ref"
)  +
  scale_fill_viridis_c() +
  scale_colour_manual(values = c("red", "black")) +
  labs(fill = "CPM") +
  theme_bw()

## ----plot-empty-hfgc----------------------------------------------------------
data("grch37.cytobands")
gr <- range(res_gr)
plotHFGC(gr, cytobands = grch37.cytobands)

## ----add-genes----------------------------------------------------------------
data("ex_trans")
plotHFGC(gr, genes = ex_trans, cytobands = grch37.cytobands)

## ----plot-trans---------------------------------------------------------------
plotHFGC(
  gr, 
  genes = ex_trans, genecol = "wheat",
  collapseTranscripts = FALSE,
  cytobands = grch37.cytobands, zoom = 1.2
)

## ----split-trans--------------------------------------------------------------
status_trans <- splitAsList(ex_trans, ex_trans$status)
plotHFGC(
  gr, 
  genes = status_trans, 
  genecol = list(Up = "forestgreen", Unchanged = "grey", Undetected = "grey80"),
  collapseTranscripts = list(Up = FALSE, Unchanged = FALSE, Undetected = "meta"),
  cytobands = grch37.cytobands, zoom = 1.2
)

## ----plot-hfgc----------------------------------------------------------------
data("ex_prom")
feat_grl <- GRangesList(Promoters = ex_prom)
plotHFGC(
  gr, 
  features = feat_grl, featcol = list(Promoters = "red"),
  genes = status_trans, 
  genecol = list(Up = "forestgreen", Unchanged = "grey", Undetected = "grey80"),
  collapseTranscripts = list(Up = FALSE, Unchanged = FALSE, Undetected = "meta"),
  cytobands = grch37.cytobands, zoom = 1.2
)

## ----plot-with-hic------------------------------------------------------------
plotHFGC(
  gr, 
  hic = ex_hic,
  features = feat_grl, featcol = list(Promoters = "red"),
  genes = status_trans, 
  genecol = list(Up = "forestgreen", Unchanged = "grey", Undetected = "grey80"),
  collapseTranscripts = list(Up = FALSE, Unchanged = FALSE, Undetected = "meta"),
  cytobands = grch37.cytobands, zoom = 1.2
)

## ----plot-with-coverage-------------------------------------------------------
plotHFGC(
  gr, 
  hic = ex_hic,
  features = feat_grl, featcol = list(Promoters = "red"),
  genes = status_trans, 
  coverage = bwfl, linecol = c(ex1 = "#4B0055", ex2 = "#007094"),
  genecol = list(Up = "forestgreen", Unchanged = "grey", Undetected = "grey80"),
  collapseTranscripts = list(Up = FALSE, Unchanged = FALSE, Undetected = "meta"),
  cytobands = grch37.cytobands, zoom = 1.2
)

## ----plot-with-tracks---------------------------------------------------------
cov_list <- list(TF1 = bwfl)
plotHFGC(
  gr, 
  hic = ex_hic,
  features = feat_grl, featcol = list(Promoters = "red"),
  genes = status_trans, 
  coverage = cov_list, 
  linecol = list(TF1 = c(ex1 = "#4B0055", ex2 = "#007094")),
  genecol = list(Up = "forestgreen", Unchanged = "grey", Undetected = "grey80"),
  collapseTranscripts = list(Up = FALSE, Unchanged = FALSE, Undetected = "meta"),
  cytobands = grch37.cytobands, zoom = 1.2
)

## ----cov-annot----------------------------------------------------------------
cov_annot <- splitAsList(res_gr, res_gr$logFC < -1) %>% 
  setNames(c("Unchanged", "Decreased")) %>% 
  endoapply(granges)

## ----plot-annot---------------------------------------------------------------
plotHFGC(
  gr, 
  hic = ex_hic,
  features = feat_grl, featcol = list(Promoters = "red"),
  genes = status_trans, 
  coverage = cov_list, 
  annotation = list(TF1 = cov_annot), 
  annotcol = c(Unchanged = "grey", Decreased = "#3333CC"),
  linecol = list(TF1 = c(ex1 = "#4B0055", ex2 = "#007094")),
  genecol = list(Up = "forestgreen", Unchanged = "grey", Undetected = "grey80"),
  collapseTranscripts = "meta",
  cytobands = grch37.cytobands, zoom = 1.2
)

## -----------------------------------------------------------------------------
sessionInfo()

