---
title: "extraChIPs: Differential Signal Using Fixed-Width Windows"
author:
- name: Stevie Pederson
  affiliation: 
  - Black Ochre Data Laboratories, Telethon Kids Institute, Adelaide, Australia
  - Dame Roma Mitchell Cancer Researc Laboratories, University of Adelaide
  - John Curtin School of Medical Research, Australian National University
  email: stephen.pederson.au@gmail.com
package: extraChIPs
bibliography: '`r system.file("references.bib", package = "extraChIPs")`'
output:
  BiocStyle::html_document
vignette: >
  %\VignetteIndexEntry{Differential Signal Analysis (Fixed-Width Windows)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup}
knitr::opts_chunk$set(message = FALSE, crop = NULL)
```


# Introduction

The [GRAVI](https://github.com/smped/GRAVI) workflow, for which this
package is designed, uses sliding windows for differential signal
analysis. However, the use of fixed-width windows, as is common under
DiffBind-style (Ross-Innes et al. 2012) approaches is also possible with
`extraChIPs`. This vignette focusses on using conventional peak calls
and fixed-width approaches to replicate and extend these approaches.

The majority of examples below use heavily reduced datasets to provide
general guidance on using the functions. Some results may appear trivial
as a result, but will hopefully prove far more useful in a true
experimental context. All data, along with this vignette are available
[here](https://github.com/smped/extraChIPs_vignette). Please place all
contents of the data directory in a directory named data in your own
working directory.

# Setup

## Installation

In order to use the package `extraChIPs` and follow this vignette, we
recommend using the package `BiocManager` hosted on CRAN. Once this is
installed, the additional packages required for this vignette
(`tidyverse`, `Rsamtools`, `csaw`, `BiocParallel` and `rtracklayer`) can
also be installed.

``` r 
if (!"BiocManager" %in% rownames(installed.packages()))
  install.packages("BiocManager")
pkg <- c(
  "tidyverse", "Rsamtools", "csaw", "BiocParallel", "rtracklayer", "edgeR", 
  "patchwork", "extraChIPs", "plyranges", "scales", "here", "quantro"
)
BiocManager::install(pkg, update = FALSE)
```

Once these packages are installed, we can load them easily

```{r load-packages}
library(tidyverse)
library(Rsamtools)
library(csaw)
library(BiocParallel)
library(rtracklayer)
library(edgeR)
library(patchwork)
library(extraChIPs)
library(plyranges)
library(scales)
library(glue)
library(ggrepel)
library(here)
library(quantro)
theme_set(theme_bw())
```

## Data

All data for this vignette is expected to be in a sub-directory of the
working directory named “data”, and all paths will be predicated on
this. Please ensure you have all data in this location, obtained from
[here](https://github.com/smped/extraChIPs_vignette).

The data itself is ChIP-Seq data targeting the Estrogen Receptor (ER),
and is taken from the cell-line ZR-75-1 cell-line using data from the
BioProject , Pre-processing was performed using the
[`prepareChIPs`](https://github.com/smped/prepareChIPs) workflow,
written in snakemake (Mölder et al. 2021) and all code is available at
<https://github.com/smped/PRJNA509779>. ER binding was assessed under
Vehicle (E2) and DHT-stimulated (E2DHT) conditions. Using GRCh37 as the
reference genome, a subset of regions found on chromosome 10 are
included in this dataset for simplicity.

First we’ll define our sample data then define our two treatment groups.
Defining a consistent colour palette for all plots is also a good habit
to develop.

```{r samples}
treat_levels <- c("E2", "E2DHT")
treat_colours <- setNames(c("steelblue", "red3"), treat_levels)
samples <- tibble(
  accession = paste0("SRR831518", seq(0, 5)),
  target = "ER",
  treatment = factor(rep(treat_levels, each = 3), levels = treat_levels)
)
samples
```


```{r accessions}
accessions <- samples %>% 
  split(f = .$treatment) %>% 
  lapply(pull, "accession")
```

We’ll eventually be loading counts for differential signal analysis from
a set of BamFiles, so first we’ll create a `BamFileList` with all of
these files. 

```{r bfl, eval = FALSE}
bfl <- here("data", "ER", glue("{samples$accession}.bam")) %>% 
  BamFileList() %>% 
  setNames(str_remove_all(names(.), ".bam"))
file.exists(path(bfl))
```

    ## [1] TRUE TRUE TRUE TRUE TRUE TRUE

This also enables creation of a `Seqinfo` object based on
the actual reference genome to which the reads were aligned during data
preparation. `Seqinfo` objects are the foundation of working with
GRanges, so defining an object at the start of a workflow is good
practice.
As is common practice for ChIP-Seq analyses, we'll restrict our focus to the 
standard chromosomes.

```{r sq, eval = FALSE}
sq <- seqinfo(bfl)
sq <- keepStandardChromosomes(sq)
isCircular(sq) <- rep(FALSE, length(seqlevels(sq)))
genome(sq) <- "GRCh37"
```

Another key preparatory step for working with peaks is to define a set
of regions as either blacklisted or grey-listed regions. The former are
known problematic regions based on each genome, with data freely
available from
<https://github.com/Boyle-Lab/Blacklist/tree/master/lists>, whilst
grey-listed regions are defined from potentially problematic regions as
detected within the input sample. For our samples code for this is
included in the previously provided repository
(<https://github.com/smped/PRJNA509779>).

```{r omit-ranges, eval = FALSE}
greylist <- import.bed(here("data/chr10_greylist.bed"), seqinfo = sq)
blacklist <- import.bed( here("data/chr10_blacklist.bed"), seqinfo = sq)
omit_ranges <- c(greylist, blacklist)
```

# Working With Peaks

The provided dataset includes six files produced by `macs2 callpeak`
(Zhang et al. 2008) in the `narrowPeak` format, and these are able to be
easily parsed using `extraChIPs`. We’ll immediately pass our black &
grey-listed regions to our parsing function so we can exclude these
regions right from the start. 
By passing the above seqinfo object to this function, we're also telling 
`importPeaks()` to ignore any peaks not on the included chromosomes.

```{r fake-peaks, eval = FALSE}
peaks <- here("data", "ER", glue("{samples$accession}_peaks.narrowPeak")) %>% 
  importPeaks(seqinfo = sq, blacklist = omit_ranges)
```

```{r load-peaks, echo = FALSE}
data("peaks")
```


This will import the peaks from all files as a single `GRangesList`
object, adding the file-name to each element by default. We can easily
modify these names if we so wish.

```{r edit-peaknames}
names(peaks) <- str_remove_all(names(peaks), "_peaks.narrowPeak")
```

Once loaded, we can easily check how similar our replicates are using
`plotOverlaps()`. When three or more sets of peaks are contained in the
`GRangesList`, an UpSet plot will be drawn by default.

```{r plot-overlaps, eval = FALSE, fig.cap = "*UpSet plot showing overlapping peaks across all replicates*"}
plotOverlaps(
  peaks, min_size = 10, .sort_sets = FALSE, 
  set_col = treat_colours[as.character(samples$treatment)]
)
```

Optionally, specifying a column and a suitable function will produce an
additional panel summarising that value. In the following, we’ll show
the maximum score obtained, highlighting that for peaks identified in
only one or two replicates, the overall signal intensity is generally
lower, even in the sample with the strongest signal.

```{r plot-overlaps-score, fig.height = 7, fig.cap = "*UpSet plot showing overlapping peaks across all replicates, with the maximum score across all replicates shown in the upper panel.*"}
plotOverlaps(
  peaks, min_size = 10, .sort_sets = FALSE, var = "score", f = "max",
   set_col = treat_colours[as.character(samples$treatment)]
)
```

A common task at this point may be to define consensus peaks within each
treatment group, by retaining only the peaks found in 2 of the 3
replicates `(p = 2/3)`. The default approach is to take the union of all
ranges, with the returned object containing logical values for each
sample, as well as the number of samples where an overlapping peak was
found.

If we wish to retain any of the original columns, such as the
`macs2 callpeak` score, we can simply pass the column names to
`makeConsensus()`

```{r make-consensus} 
consensus_e2 <- makeConsensus(peaks[accessions$E2], p = 2/3, var = "score")
consensus_e2dht <- makeConsensus(peaks[accessions$E2DHT], p = 2/3, var = "score")
```

Alternatively, we could find the centre of the peaks as part of this
process, by averaging across the estimated peak centres for each sample.
This is a very common step for ChIP-Seq data where the target is a
*transcription factor*, and also forms a key step in the DiffBind
workflow.

In the following code chunk, we first find the centre for each sample
using the information provided by `macs2`, before retaining this column
when calling `makeConsensus()`. This will return each of the individual
centre-position estimates as a list for each merged range, and using
`vapply()` we then take the mean position as our estimate for the
combined peak centre.

```{r make-consensus-centre}
consensus_e2 <- peaks[accessions$E2] %>% 
  endoapply(mutate, centre = start + peak) %>% 
  makeConsensus(p = 2/3, var = "centre") %>% 
  mutate(centre = vapply(centre, mean, numeric(1)))
consensus_e2
consensus_e2dht <- peaks[accessions$E2DHT] %>% 
  endoapply(mutate, centre = start + peak) %>% 
  makeConsensus(p = 2/3, var = "centre") %>% 
  mutate(centre = vapply(centre, mean, numeric(1)))
```

We can also inspect these using `plotOverlaps()` provided we use a
`GRangesList` for the input. Now that we only have two elements (one for
each treatment) a VennDiagram will be generated instead of an UpSet
plot.

```{r venn-consensus-overlap, fig.width=6, fig.cap = "*Overlap between consensus peaks identified in a treatment-specific manner*"}
GRangesList(E2 = granges(consensus_e2), E2DHT = granges(consensus_e2dht)) %>% 
  plotOverlaps(set_col = treat_colours[treat_levels])
```


We can now go one step further and define the set of peaks found in
either treatment. Given we’re being inclusive here, we can leave p = 0
so any peak found in *either treatment* is included.

```{r union-peaks}
union_peaks <- GRangesList(
  E2 = select(consensus_e2, centre), 
  E2DHT = select(consensus_e2dht, centre)
) %>% 
  makeConsensus(var = c("centre")) %>% 
  mutate(
    centre = vapply(centre, mean, numeric(1)) %>% round(0)
  ) 
```

Now we have a set of peaks, found in at least 2/3 of samples from either
condition, with estimates of each peak’s centre. The next step would be
to set all peaks as the same width based on the centre position, with a
common width being 500bp.

In the following we’ll perform multiple operations in a single call
mutate, so let’s make sure we know what’s happening.

1.  `glue("{seqnames}:{centre}:{strand}")` uses `glue` syntax to parse
    the seqnames, centre position and strand information as a
    character-like vector with a width of only 1, and using the
    estimated centre as the Range.
2.  We then coerce this to a `GRanges` object, before resizing to the
    desired width.
3.  We also add the original (un-centred) range as an additional column,
    retaining the `GRanges` structure, but discarding anything in the
    `mcols()` element, then
4.  Using `colToRanges()`, we take the centred ranges and place them as
    the core set of GRanges for this object.

This gives a GRanges object with all original information, but with
centred peaks of a fixed width.

```{r centred-peaks, eval = FALSE}
w <- 500
centred_peaks <- union_peaks %>% 
  mutate(
    centre = glue("{seqnames}:{centre}:{strand}") %>% 
      GRanges(seqinfo = sq) %>% 
      resize(width = w),
    union_peak = granges(.)
  ) %>% 
  colToRanges("centre")
```

# Counting Reads

Now we have our centred, fixed-width peaks, we can count reads using
`csaw::regionCounts()` (Lun and Smyth 2016). We know our fragment length
is about 200bp, so we can pass this to the function for a slightly more
sophisticated approach to counting.

```{r se, echo = FALSE}
data("se")
```

```{r fake-se, eval = FALSE}
se <- regionCounts(bfl, centred_peaks, ext = 200)
```

```{r show-se}
se
```


The `colData()` element of the returned object as the columns
*bam.files*, *totals*, *ext* and *rlen*, which are all informative and
can be supplemented with our `samples` data frame. In the following,
we’ll 1) coerce to a `tibble`, 2) `left_join()` the `samples` object, 3)
add the accession as the sample column, 4) set the accession back as the
rownames, then 5) coerce back to the required `DataFrame()` structure.

```{r edit-coldata}
colData(se) <- colData(se) %>% 
  as_tibble(rownames = "accession") %>% 
  left_join(samples) %>% 
  mutate(sample = accession) %>% 
  as.data.frame() %>% 
  column_to_rownames("accession") %>% 
  DataFrame()
colData(se)
```


For QC and visualisation, we can add an additional `logCPM` assay to our
object as well.

```{r add-logcpm}
assay(se, "logCPM") <- cpm(assay(se, "counts"), lib.size = se$totals, log = TRUE)
```

First we might like to check our distribution of counts

```{r plot-assay-densities, fig.cap = "*Count densities for all samples, using the log+1 transformation*"}
plotAssayDensities(se, assay = "counts", colour = "treat", trans = "log1p") +
  scale_colour_manual(values = treat_colours)
```

A PCA plot can also provide insight as to where the variability in the
data lies.

```{r plot-pca, fig.cap = "*PCA plot using logCPM values and showing that replicate variability is larger than variability between treatment groups.*"}
plotAssayPCA(se, assay = "logCPM", colour = "treat", label = "sample") +
  scale_colour_manual(values = treat_colours)
```


# Differential Signal Analysis

## Statistical Testing

In order to perform Differential Signal Analysis, we simply need to
define a model matrix, as for conventional analysis using `edgeR` or
`limma`. We can then pass this, along with our fixed-width counts to
`fitAssayDiff()`. By default normalisation will be *library-size*
normalisation, as is a common default strategy for ChIP-Seq data. In
contrast to sliding window approaches, these results represent our final
results and there is no need for merging windows.

```{r ls-fit}
X <- model.matrix(~treatment, data = colData(se))
ls_res <- fitAssayDiff(se, design = X, asRanges = TRUE)
sum(ls_res$FDR < 0.05)
```


TMM normalisation (Robinson and Oshlack 2010) is another common
strategy, which relies on the data from all treatment groups being drawn
from the same distributions. We can formally test this using the package
`quantro` (Hicks and Irizarry 2015) , which produces p-values for 1)
H<sub>0</sub>: Group medians are drawn from the same distribution, and
2) H<sub>0</sub>: Group-specific distributions are the same.

```{r qtest}
set.seed(100)
qtest <- assay(se, "counts") %>% 
  quantro(groupFactor = se$treatment, B = 1e3)
qtest
```


Here, both p-values are \>0.05, so in conjunction with out visual
inspection earlier, we can confidently apply TMM normalisation. To apply
this, we simply specify the argument `norm = "TMM"` when we call
`fitAssayDiff()`. In the analysis below, we’ve also specified a
fold-change threshold `(fc = 1.2)`, below which, changes in signal are
considered to not be of interest (McCarthy and Smyth 2009). This
threshold is incorporated into the testing so there is no requirement
for *post-hoc* filtering based on a threshold.

```{r tmm-res}
tmm_res <- fitAssayDiff(se, design = X, norm = "TMM", asRanges = TRUE, fc = 1.2)
sum(tmm_res$FDR < 0.05)
```

An MA-plot is a common way of inspecting results and in the following we
use the original ‘union_peak’ in our labelling of points. This serves as
a reminder that the fixed-width windows are in fact *a proxy* for the
entire region for which we have confidently detected ChIP signal, and
that these windows are truly the regions of interest.

```{r plot-ma, fig.cap = "*MA-plot after fitting using TMM normalisation and applying a fold-change threshold during testing. Points are labelled using the original windows obtained when merging replicates and treatment groups.*"}
tmm_res %>% 
  as_tibble() %>% 
  mutate(`FDR < 0.05` = FDR < 0.05) %>% 
  ggplot(aes(logCPM, logFC)) +
  geom_point(aes(colour = `FDR < 0.05`)) +
  geom_smooth(se = FALSE) +
  geom_label_repel(
    aes(label = union_peak), colour = "red",
    data = . %>% dplyr::filter(FDR < 0.05)
  ) +
  scale_colour_manual(values = c("black", "red"))
```


## Mapping to Genes

Whilst knowledge of which regions are showing differential signal, the
fundamental question we are usually asking is about the downstream
regulatory consequences, such as the target gene. Before we can map
peaks to genes, we’ll need to define our genes. In the following, we’ll
use the provided Gencode gene mappings at the gene, transcript and exon
level.

```{r gencode, eval = FALSE}
gencode <- here("data/gencode.v43lift37.chr10.annotation.gtf.gz") %>% 
  import.gff() %>% 
  filter_by_overlaps(GRanges("chr10:42354900-100000000")) %>% 
  split(.$type)
seqlevels(gencode) <- seqlevels(sq)
seqinfo(gencode) <- sq
```

Mapping to genes using `mapByFeature()` uses additional annotations,
such as whether the peak overlaps a promoter, enhancer or long-range
interaction. Here we’ll just use promoters, so let’s create a set of
promoters from our transcript-level information, ensuring we incorporate
all possible promoters within a gene, and merging any overlapping ranges
using `reduceMC()`

```{r promoters, eval = FALSE}
promoters <- gencode$transcript %>% 
    select(gene_id, ends_with("name")) %>% 
    promoters(upstream = 2500, downstream = 500) %>% 
    reduceMC(simplify = FALSE)
promoters
```

    ## GRanges object with 1678 ranges and 3 metadata columns:
    ##          seqnames            ranges strand |
    ##             <Rle>         <IRanges>  <Rle> |
    ##      [1]    chr10 42678287-42681286      + |
    ##      [2]    chr10 42702938-42705937      + |
    ##      [3]    chr10 42735669-42738668      + |
    ##      [4]    chr10 42743933-42746932      + |
    ##      [5]    chr10 42968428-42973155      + |
    ##      ...      ...               ...    ... .
    ##   [1674]    chr10 99635155-99638154      - |
    ##   [1675]    chr10 99643500-99646805      - |
    ##   [1676]    chr10 99695536-99698535      - |
    ##   [1677]    chr10 99770595-99773594      - |
    ##   [1678]    chr10 99789879-99793085      - |
    ##                                                                     gene_id
    ##                                                             <CharacterList>
    ##      [1]                                                ENSG00000237592.2_5
    ##      [2]                                                ENSG00000271650.1_7
    ##      [3]                                                ENSG00000290458.1_2
    ##      [4]                                                ENSG00000274167.5_8
    ##      [5] ENSG00000185904.12_9,ENSG00000185904.12_9,ENSG00000185904.12_9,...
    ##      ...                                                                ...
    ##   [1674]                                                  ENSG00000265398.1
    ##   [1675]                        ENSG00000095713.14_13,ENSG00000095713.14_13
    ##   [1676]                                              ENSG00000095713.14_13
    ##   [1677]                                              ENSG00000095713.14_13
    ##   [1678]                        ENSG00000095713.14_13,ENSG00000095713.14_13
    ##                                  gene_name
    ##                            <CharacterList>
    ##      [1]                       IGKV1OR10-1
    ##      [2]                   ENSG00000271650
    ##      [3]                   ENSG00000290458
    ##      [4]                   ENSG00000274167
    ##      [5] LINC00839,LINC00839,LINC00839,...
    ##      ...                               ...
    ##   [1674]                        AL139239.1
    ##   [1675]                     CRTAC1,CRTAC1
    ##   [1676]                            CRTAC1
    ##   [1677]                            CRTAC1
    ##   [1678]                     CRTAC1,CRTAC1
    ##                                        transcript_name
    ##                                        <CharacterList>
    ##      [1]                               IGKV1OR10-1-201
    ##      [2]                               ENST00000605702
    ##      [3]                               ENST00000622823
    ##      [4]                               ENST00000622650
    ##      [5] LINC00839-204,LINC00839-203,LINC00839-202,...
    ##      ...                                           ...
    ##   [1674]                                AL139239.1-201
    ##   [1675]                         CRTAC1-205,CRTAC1-206
    ##   [1676]                                    CRTAC1-204
    ##   [1677]                                    CRTAC1-201
    ##   [1678]                         CRTAC1-203,CRTAC1-202
    ##   -------
    ##   seqinfo: 25 sequences from GRCh37 genome

Now we’ll pass these to `mapByFeature()`, but first, we’ll place the
original ‘union_peak’ back as the core of the GRanges object. This will
retain all the results from testing, but ensures the correct region is
mapped to genes.

```{r tmm-mapped-res, eval = FALSE}
tmm_mapped_res <- tmm_res %>% 
  colToRanges("union_peak") %>% 
  mapByFeature(genes = gencode$gene, prom = promoters) %>% 
  mutate(
    status = case_when(
      FDR >= .05 ~ "Unchanged",
      logFC > 0 ~ "Increased",
      logFC < 0 ~ "Decreased"
    )
  )
arrange(tmm_mapped_res, PValue)
```

    ## GRanges object with 188 ranges and 10 metadata columns:
    ##         seqnames            ranges strand |        E2     E2DHT         n
    ##            <Rle>         <IRanges>  <Rle> | <logical> <logical> <numeric>
    ##     [1]    chr10 81101906-81102928      * |      TRUE      TRUE         2
    ##     [2]    chr10 79629641-79630271      * |      TRUE      TRUE         2
    ##     [3]    chr10 89407752-89408138      * |     FALSE      TRUE         1
    ##     [4]    chr10 52233596-52233998      * |      TRUE      TRUE         2
    ##     [5]    chr10 91651138-91651433      * |     FALSE      TRUE         1
    ##     ...      ...               ...    ... .       ...       ...       ...
    ##   [184]    chr10 89395292-89395626      * |      TRUE      TRUE         2
    ##   [185]    chr10 82723984-82724328      * |      TRUE      TRUE         2
    ##   [186]    chr10 93120411-93121224      * |      TRUE      TRUE         2
    ##   [187]    chr10 95755308-95755721      * |      TRUE      TRUE         2
    ##   [188]    chr10 79190987-79191351      * |     FALSE      TRUE         1
    ##               logFC    logCPM      PValue         FDR
    ##           <numeric> <numeric>   <numeric>   <numeric>
    ##     [1]    1.880493   7.91492 2.50406e-25 4.70764e-23
    ##     [2]    0.940221   8.07204 8.07816e-08 7.59347e-06
    ##     [3]    1.578056   6.16574 4.18078e-07 2.61995e-05
    ##     [4]    1.057814   6.67466 5.73897e-05 2.69732e-03
    ##     [5]    1.546530   5.15431 1.75562e-04 6.60112e-03
    ##     ...         ...       ...         ...         ...
    ##   [184] -0.01346141   6.36140    0.961749    0.982656
    ##   [185] -0.00855862   6.91968    0.972422    0.988154
    ##   [186] -0.01173559   9.32397    0.979561    0.988154
    ##   [187]  0.00721240   5.88042    0.982898    0.988154
    ##   [188] -0.00157801   6.41182    0.994914    0.994914
    ##                                            gene_id
    ##                                    <CharacterList>
    ##     [1]                       ENSG00000108179.14_6
    ##     [2]                      ENSG00000151208.17_11
    ##     [3]   ENSG00000225913.2_9,ENSG00000196566.2_10
    ##     [4] ENSG00000198964.14_10,ENSG00000225303.2_10
    ##     [5]                        ENSG00000280560.3_9
    ##     ...                                        ...
    ##   [184]   ENSG00000225913.2_9,ENSG00000196566.2_10
    ##   [185]                        ENSG00000227209.1_5
    ##   [186]                        ENSG00000289228.2_2
    ##   [187]                      ENSG00000138193.17_12
    ##   [188]                      ENSG00000156113.25_17
    ##                               gene_name      status
    ##                         <CharacterList> <character>
    ##     [1]                            PPIF   Increased
    ##     [2]                            DLG5   Increased
    ##     [3] ENSG00000225913,ENSG00000196566   Increased
    ##     [4]           SGMS1,ENSG00000225303   Increased
    ##     [5]                       LINC01374   Increased
    ##     ...                             ...         ...
    ##   [184] ENSG00000225913,ENSG00000196566   Unchanged
    ##   [185]                         WARS2P1   Unchanged
    ##   [186]                 ENSG00000289228   Unchanged
    ##   [187]                           PLCE1   Unchanged
    ##   [188]                          KCNMA1   Unchanged
    ##   -------
    ##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

# Inspection of Results

## Profile Heatmaps

When analysing a transcription factor, checking the binding profile
across our treatment groups can be informative, and is often performed
using ‘Profile Heatmaps’ where coverage is smoothed within bins
surrounding our peak centre.

The function `getProfileData()` takes a set of ranges and a
BigWigFileList, and performs the smoothing, which is then passed to the
function `plotProfileHeatmap()`.

The following shows the three steps of 1) defining the ranges, 2)
obtaining the smoothed binding profiles, and 3) drawing the heatmap.
Note that we can facet the heatmaps by selecting the ‘status’ column to
separate any Increased or Decreased regions. By default, this will also
draw the smoothed lines in the top panel using different colours.

These plots can be used to show coverage-like values (SPMR or CPM) or we
can use fold-enrichment over the input sample(s), as is also produced by
`macs2 bdgcmp`. This data isn’t generally visualised using
log-transformation so we’ll set `log = FALSE` in our call to
`getProfileData()`

```{r profile-heatmap, eval = FALSE}
fe_bw <- here("data", "ER", glue("{treat_levels}_FE_chr10.bw")) %>% 
  BigWigFileList() %>% 
  setNames(treat_levels)
sig_ranges <- filter(tmm_mapped_res, FDR < 0.05)
pd_fe <- getProfileData(fe_bw, sig_ranges, log = FALSE) 
pd_fe %>% 
  plotProfileHeatmap("profile_data", facetY = "status") +
  scale_fill_gradient(low = "white", high = "red") +
  labs(fill = "Fold\nEnrichment")
```


## Session Info

```{r session-info}
sessionInfo()
```

 
## References

<div id="refs" class="references csl-bib-body hanging-indent">

<div id="ref-gviz" class="csl-entry">

Hahne, Florian, and Robert Ivanek. 2016. “Statistical Genomics: Methods
and Protocols.” In, edited by Ewy Mathé and Sean Davis, 335–51. New
York, NY: Springer New York.
<https://doi.org/10.1007/978-1-4939-3578-9_16>.

</div>

<div id="ref-Hicks2015-ee" class="csl-entry">

Hicks, Stephanie C, and Rafael A Irizarry. 2015. “Quantro: A Data-Driven
Approach to Guide the Choice of an Appropriate Normalization Method.”
*Genome Biol.* 16 (1): 117.

</div>

<div id="ref-csaw2016" class="csl-entry">

Lun, Aaron T L, and Gordon K Smyth. 2016. “Csaw: A Bioconductor Package
for Differential Binding Analysis of ChIP-Seq Data Using Sliding
Windows.” *Nucleic Acids Res.* 44 (5): e45.

</div>

<div id="ref-McCarthy2009-qf" class="csl-entry">

McCarthy, Davis J, and Gordon K Smyth. 2009. “Testing Significance
Relative to a Fold-Change Threshold Is a TREAT.” *Bioinformatics* 25
(6): 765–71.

</div>

<div id="ref-Molder2021-mo" class="csl-entry">

Mölder, Felix, Kim Philipp Jablonski, Brice Letcher, Michael B Hall,
Christopher H Tomkins-Tinch, Vanessa Sochat, Jan Forster, et al. 2021.
“Sustainable Data Analysis with Snakemake.” *F1000Res.* 10 (January):
33.

</div>

<div id="ref-Robinson2010-qp" class="csl-entry">

Robinson, Mark D, and Alicia Oshlack. 2010. “A Scaling Normalization
Method for Differential Expression Analysis of
<span class="nocase">RNA-seq</span> Data.” *Genome Biol.* 11 (3): R25.

</div>

<div id="ref-DiffBind2012" class="csl-entry">

Ross-Innes, Caryn S., Rory Stark, Andrew E. Teschendorff, Kelly A.
Holmes, H. Raza Ali, Mark J. Dunning, Gordon D. Brown, et al. 2012.
“Differential Oestrogen Receptor Binding Is Associated with Clinical
Outcome in Breast Cancer.” *Nature* 481: –4.
<http://www.nature.com/nature/journal/v481/n7381/full/nature10730.html>.

</div>

<div id="ref-Zhang2008-ms" class="csl-entry">

Zhang, Yong, Tao Liu, Clifford A Meyer, Jérôme Eeckhoute, David S
Johnson, Bradley E Bernstein, Chad Nusbaum, et al. 2008. “Model-Based
Analysis of ChIP-Seq (MACS).” *Genome Biol.* 9 (9): R137.

</div>

</div>

<br>
