R version: R version 3.6.1 (2019-07-05)

Bioconductor version: 3.10

CAGEfightR version: 1.6.0

1 Background

Transcriptional regulation is one of the most important aspects of gene expression. Transcription Start Sites (TSSs) are focal points of this process: The TSS act as an integration point for a wide range of molecular cues from surrounding genomic areas to determine transcription and ultimately expression levels. These include proximal factors such as chromatin accessibility, chromatin modification, DNA methylation and transcription factor binding, and distal factors such as enhancer activity and chromatin confirmation (Smale and Kadonaga 2003; Kadonaga 2012; Lenhard, Sandelin, and Carninci 2012; Haberle and Stark 2018).

Cap Analysis of Gene Expression (CAGE) has emerged as one of the dominant high-throughput assays for studying TSSs (Adiconis et al. 2018). CAGE is based on “cap trapping”: capturing capped full-length RNAs and sequencing only the first 20-30 nucleotides from the 5’-end, so-called CAGE tags (Takahashi et al. 2012). When mapped to a reference genome, the 5’-ends of CAGE tag identify the actual TSS for respective RNA with basepair-level accuracy. Basepair-accurate TSSs identified this way are referred to as CAGE Transcription Start Sites (CTSSs). RNA polymerase rarely initiates from just a single nucleotide: this is manifested in CAGE data by the fact that CTSSs are mostly found in tightly spaced groups on the same strand. The majority of CAGE studies have merged such CTSSs into genomic blocks typically referred to as Tag Clusters (TCs), using a variety of clustering methods (see below). TCs are often interpreted as TSSs in the more general sense, given that most genes have many CTSSs, but only a few TCs that represent a few major transcripts with highly similar CTSSs (Carninci et al. 2006; Sandelin et al. 2007). Since the number of mapped CAGE tags in a given TC is indicative of the number of RNAs from that region, the number of CAGE tags falling in given TC can be used as a measure of expression (Kawaji et al. 2014).

As CAGE tags can be mapped to a reference genome without the need for transcript annotations, it can detect TSSs of known mRNAs, but also novel alternative TSSs (that might be condition or tissue dependent) (Carninci et al. 2006; Consortium, Pmi, and Dgt 2014). Since CAGE captures all capped RNAs, it can also identify long non-coding RNA (lincRNA) (Hon et al. 2017). It was previously shown that active enhancers are characterized by balanced bidirectional transcription producing enhancer RNAs (eRNAs), making it possible to predict enhancer regions and quantify their expression levels from CAGE data alone (Kim et al. 2010; Andersson et al. 2014). Thus, CAGE data can predict the locations and activity of mRNAs, lincRNAs and enhancers in a single assay, providing a comprehensive view of transcriptional regulation across both inter- and intragenic regions.

Bioconductor contains a vast collection of tools for analyzing transcriptomics datasets, in particular the widely used RNA-Seq and microarray assays(Huber et al. 2015). Only a few packages are dedicated to analyzing 5’-end data in general and CAGE data in particular: TSRchitect (Taylor Raborn, Brendel, and Sridharan, n.d.), icetea (Bhardwaj 2019), CAGEr (Haberle et al. 2015) and CAGEfightR (Thodberg et al. 2019) (Table 1).

CAGEr was the first package solely dedicated to the analysis of CAGE data and was recently updated to more closely adhere to Bioconductor S4-class standards. CAGEr takes as input aligned reads in the form of BAM-files and can identify, quantify, characterize and annotate TSSs. TSSs are found in individual samples using either simple clustering of CTSSs (greedy or distance-based clustering) or the more advanced density-based paraclu clustering method(Frith et al. 2008), and can be aggregated across samples to a set of consensus clusters. Several specialized routines for CAGE data is available, such as G-bias correction of mapped tags, power law normalization of CTSS counts and fine-grained TSS shifts. Finally, CAGEr offers easy interface to the large collection of CAGE data from the FANTOM consortium (Consortium, Pmi, and Dgt 2014). TSRchitect and icetea are two more recent additions to Bioconductor. While being less comprehensive, they aim to be more general and handle more types of 5’-end methods that are conceptually similar to CAGE (RAMPAGE, PEAT, PRO-Cap, etc. (Adiconis et al. 2018)). Both packages can identify, quantify and annotate TSSs, with TSRchitect using an X-means algorithm and icetea using a sliding window approach. icetea offers the unique feature of mapping reads to a reference genome by interfacing with Rsubread. Both CAGEr and icetea offers built-in capabilities for differential expression (DE) analysis via the popular DESeq2 or edgeR packages (Love, Huber, and Anders 2014; Robinson, McCarthy, and Smyth 2010).

CAGEfightR is a recent addition to Bioconductor focused on analyzing CAGE data, but applicable to most 5’-end data. It aims to be general and flexible to allow for easy interfacing with the wealth of other Bioconductor packages. CAGEfightR takes CTSSs stored in BigWig-files as input and uses only standard Bioconductor S4-classes (GenomicRanges, SummarizedExperiment, InteractionSet(Lawrence et al. 2013; Lun, Perry, and Ing-Simmons 2016)) making it easy for users to learn and combine CAGEfightR with functions from other Bioconductor packages (e.g. instead of providing custom wrappers around other packages such as differential expression analysis). In addition to TSS analysis, CAGEfightR is the only package on Bioconductor to also offer functions for enhancer analysis based on CAGE (and similarly scoped) data. This includes enhancer identification and quantification, linking enhancers to TSSs via correlation of expression and finding enhancer clusters, often referred to as stretch- or super enhancers.

Table 1: Comparison of Bioconductor packages for CAGE data analysis.
Analysis icetea TSRchitect CAGEr CAGEfightR
Simplest input FASTQ BAM BAM BigWig
Paired-end input + + - -
CTSS extraction - + + -
TSS calling Sliding window X-means Distance or Paraclu Slice-reduce
TSS shapes - Shape index IQR and TSS shifts IQR, entropy, etc.
Differential Expression + - + -
Enhancer calling - - - +
TSS-enhancer correlation - - - +
Super enhancers - - - +

In this workflow, we illustrate how the CAGEfightR package can be used to orchestrate an end-to-end analysis of CAGE data by making it easy to interface with a wide range of different Bioconductor packages. Highlights include TSS and enhancer candidate identification, differential expression, alternative TSS usage, enrichment of motifs, GO/KEGG terms and calculating TSS-enhancer correlations.

2 Materials and methods

2.1 Dataset

This workflow uses data from “Identification of Gene Transcription Start Sites and Enhancers Responding to Pulmonary Carbon Nanotube Exposure in Vivo” by Bornholdt et al (Bornholdt et al. 2017). This study uses mouse as a model system to investigate how carbon nanotubes affect lung tissue when inhaled. Inhaled nanotubes were previously found to produce a similar response to asbestos, potentially triggering an inflammatory response in the lung tissue leading to drastic changes in gene expression.

The dataset consists of CAGE data from mouse lung biopsies: 5 mice whose lungs were instilled with water (Ctrl) and 6 mice whose lungs were instilled with nanotubes (Nano), with CTSSs for each sample stored in BigWig-files, shown in Table 2:

Table 2: Overview of samples in the nanotube exposure experiment.
Group Biological Replicates
Ctrl 5 mice
Nano 6 mice

The data is acquired via the nanotubes data package:


2.2 R-packages

This workflow uses a large number of R-packages: Bioconductor packages are primarily used for data analysis while packages from the tidyverse are used to wrangle and plot results. All these packages are loaded prior to beginning the workflow:

# CRAN packages for data manipulation and plotting

# CAGEfightR and related packages

# Bioconductor packages for differential expression

# Bioconductor packages for enrichment analyses

# Bioconductor data packages

We also set some script-wide settings for later convenience:

# Rename these for easier access
bsg <- BSgenome.Mmusculus.UCSC.mm9
txdb <- TxDb.Mmusculus.UCSC.mm9.knownGene
odb <-

# Script-wide settings
register(MulticoreParam(3)) # Parallel execution when possible
theme_set(theme_light()) # White theme for ggplot2 figures

3 Workflow

The workflow is divided into 3 parts covering different parts of a typical CAGE data analysis:

  1. Shows how to use CAGEfightR to import CTSSs and find and quantify TSS and enhancer candidates.

  2. Shows examples of how to perform genomic analyses of CAGE clusters using core Bioconductor packages such as GenomicRanges and Biostrings. This part covers typical analyses made from CAGE data, from summarizing cluster annotation, TSS shapes and core promoter sequence analysis to more advanced spatial analyses (finding TSS-enhancer correlation links and clusters of enhancers).

  3. Shows how CAGEfightR can be used to prepare data for differential expression analysis with popular R packages, including DESeq2, limma and edgeR (Love, Huber, and Anders 2014; Ritchie et al. 2015; Robinson, McCarthy, and Smyth 2010). Borrowing from RNA-Seq terminology, differential expression can be assessed at multiple different levels: TSS- and enhancer-level, gene-level and differential TSS usage(Soneson, Love, and Robinson 2015). Once differential expression results have been obtained, they can be combined with other sources of information such as motifs from JASPAR (Mathelier et al. 2016) and GO/KEGG terms(Ashburner et al. 2000; The Gene Ontology Consortium 2019; Kanehisa and Goto 2000).

3.1 Part 1: Locating, quantifying and annotating TSSs and enhancers

Before starting the analysis, we recommend gathering all information (Filenames, groups, batches, preparation data, etc.) about the samples to be analyzed in a single data.frame, often called the design matrix. CAGEfightR can keep track of the design matrix throughout the analysis:

      caption = "The initial design matrix for the nanotubes experiment") %>%
    kable_styling(latex_options = "hold_position")
Table 3: The initial design matrix for the nanotubes experiment
Class Name BigWigPlus BigWigMinus
C547 Ctrl C547
C548 Ctrl C548
C549 Ctrl C549
C559 Ctrl C559
C560 Ctrl C560
N13 Nano N13
N14 Nano N14
N15 Nano N15
N16 Nano N16
N17 Nano N17
N18 Nano N18

3.1.1 Obtaining CTSSs

CAGEfightR starts analysis from simple CTSSs, which are the number of CAGE tag 5’-ends mapping to each basepair (bp) in the genome. CTSSs are normally stored as either BED-files, bedGraph-files or BigWig-files. As CTSSs are sparse (only are small fraction of all bps are CTSSs), these files are relatively small and thereby easily shared, many studies will make CTSSs available via online repositories such as GEO.

When preparing CTSSs from new CAGE libraries, tags are normally first barcode split, trimmed and filtered before mapping to a reference genome using standard command line tools. The exact steps will depend on the given CAGE protocol in use. CTSSs can subsequently be extracted from the resulting BAM-files one library at a time, for example using genomecov from bedtools with the -5 setting. TSRchitect and CAGEr also include functions (inputToTSS(), getCTSSs(), respectively) for obtaining CTSSs from BAM-files from within R, with CAGEr having the option of correcting for G-bias when mapping.

CAGEfightR can analyze many types of 5’-end data, as long as they can be represented in a format similar to CTSSs.

3.1.2 Importing CTSSs

CAGEfightR uses the convenient BigWigFile and BigWigFileList containers for handling CTSSs stored in BigWig-files (one file for each strand), as these allow for fast random access, inspection and summarization of the genome information stored in the files (e.g. via import(), seqinfo() and summary()). First, we need to tell CAGEfightR where to find the BigWig-files containing CTSSs on the hard drive. Normally, one would supply the paths to each file (e.g. /CAGEdata/BigWigFiles/, but here we will use data from the nanotubes data package:

# Setup paths to file on hard drive
bw_plus <- system.file("extdata", nanotubes$BigWigPlus, 
                        package = "nanotubes", 
                        mustWork = TRUE)
bw_minus <- system.file("extdata", nanotubes$BigWigMinus, 
                        package = "nanotubes", 
                        mustWork = TRUE)

# Save as named BigWigFileList
bw_plus <- BigWigFileList(bw_plus)
bw_minus <- BigWigFileList(bw_minus)
names(bw_plus) <- names(bw_minus) <- nanotubes$Name

The first step is quantifying CTSS usage across all samples. This is often one of the most time-consuming steps in a CAGEfightR analysis, but it can be sped up by using multiple cores (if available, see Materials and Methods). We also need to specify the genome, which we can get from the BSgenome.Mmusculus.UCSC.mm9 genome package:

CTSSs <- quantifyCTSSs(plusStrand = bw_plus,
                       minusStrand = bw_minus,
                       genome = seqinfo(bsg),
                       design = nanotubes)
#> Checking supplied genome compatibility...
#> Iterating over 28 genomic tiles in 11 samples using 3 worker(s)...
#> Importing CTSSs from plus strand...
#> Registered S3 method overwritten by 'pryr':
#>   method      from
#>   print.bytes Rcpp
#> Importing CTSSs from minus strand...
#> Merging strands...
#> ### CTSS summary ###
#> Number of samples: 11
#> Number of CTSSs: 9.339 millions
#> Sparsity: 81.68 %
#> Final object size: 263 MB

The circa 9 million CTSSs are stored as a RangedSummarizedExperiment, which is the standard container of high-throughput experiments in Bioconductor. We can inspect both the ranges and the CTSS counts:

# Get a summary
#> class: RangedSummarizedExperiment 
#> dim: 9338802 11 
#> metadata(0):
#> assays(1): counts
#> rownames: NULL
#> rowData names(0):
#> colnames(11): C547 C548 ... N17 N18
#> colData names(4): Class Name BigWigPlus BigWigMinus

# Extract CTSS positions
#> UnstitchedGPos object with 9338802 positions and 0 metadata columns:
#>                 seqnames       pos strand
#>                    <Rle> <integer>  <Rle>
#>         [1]         chr1   3024556      +
#>         [2]         chr1   3025704      +
#>         [3]         chr1   3025705      +
#>         [4]         chr1   3028283      +
#>         [5]         chr1   3146133      +
#>         ...          ...       ...    ...
#>   [9338798] chrUn_random   5810899      -
#>   [9338799] chrUn_random   5813784      -
#>   [9338800] chrUn_random   5880838      -
#>   [9338801] chrUn_random   5893536      -
#>   [9338802] chrUn_random   5894263      -
#>   -------
#>   seqinfo: 35 sequences (1 circular) from mm9 genome

# Extract CTSS counts
assay(CTSSs, "counts") %>%
#> 6 x 11 sparse Matrix of class "dgCMatrix"
#>    [[ suppressing 11 column names 'C547', 'C548', 'C549' ... ]]
#> [1,] . . 1 . . . . . . . .
#> [2,] . . . 1 . . . . . . .
#> [3,] . . . . 1 . . . . . .
#> [4,] . . . . 1 . . . . . .
#> [5,] . . . . . . 1 . . . .
#> [6,] . 1 . . . . . . . . .

3.1.3 Unidirectional and bidirectional clustering for finding TSS and enhancer candidates:

CAGEfightR finds clusters by calculating the pooled CTSS signal across all samples: We first normalize CTSS counts in each sample to Tags-Per-Million (TPM) values, and then sum TPM values across samples:

CTSSs <- CTSSs %>%
    calcTPM() %>%
#> Calculating library sizes...
#> Calculating TPM...

This will add several new pieces of information to CTSSs: The total number of tags in each library, a new assay called TPM, and the pooled signal for each CTSS.

We can use unidirectional clustering to locate unidirectional clusters, often simply called Tag Clusters (TCs), which are candidates for TSSs. The quickTSSs will both locate and quantify TCs in a single function call:

TCs <- quickTSSs(CTSSs)
#> Using existing score column!
#>  - Running clusterUnidirectionally:
#> Splitting by strand...
#> Slice-reduce to find clusters...
#> Calculating statistics...
#> Preparing output...
#> Tag clustering summary:
#>   Width   Count Percent
#>   Total 3602099 1e+02 %
#>     >=1 2983433 8e+01 %
#>    >=10  577786 2e+01 %
#>   >=100   40842 1e+00 %
#>  >=1000      38 1e-03 %
#>  - Running quantifyClusters:
#> Finding overlaps...
#> Aggregating within clusters...

Note: quickTSSs runs CAGEfightR with default settings. If you have larger or more noisy datasets you most likely want to do a more robust analysis with different settings. See the CAGEfightR vignette for more information.

Many of the identified TCs will be very lowly expressed. To obtain likely biologically relevant TSSs, we keep only TCs expressed at more than 1 TPM in at least 5 samples (5 samples being the size of the smallest experimental group):

TSSs <- TCs %>%
    calcTPM() %>%
    subsetBySupport(inputAssay = "TPM", 
                    unexpressed = 1, 
                    minSamples = 4)
#> Calculating library sizes...
#> Warning in calcTotalTags(object = object, inputAssay = inputAssay,
#> outputColumn = outputColumn): object already has a column named totalTags in
#> colData: It will be overwritten!
#> Calculating TPM...
#> Calculating support...
#> Subsetting...
#> Removed 3573214 out of 3602099 regions (99.2%)

This removed a large number of very lowly expressed TCs, leaving us with almost 30.000 TCs for analysis. For simplicity, we will refer to these TCs as TSS candidates, as each TC can be seen as a measure of the location and activity of the TSS of a transcript or gene. Note that this is a simplification, since a TC technically groups together several bp-accurate CTSSs.

Then we turn to bidirectional clustering for identifying bidirectional clusters (BCs). Similarly, we can use quickEnhancers to locate and quantify BCs (BCs are quantified by summing tags on both strands of the cluster):

BCs <- quickEnhancers(CTSSs)
#> Using existing score column!
#>  - Running clusterBidirectionally:
#> Pre-filtering bidirectional candidate regions...
#> Retaining for analysis: 68.3%
#> Splitting by strand...
#> Calculating windowed coverage on plus strand...
#> Calculating windowed coverage on minus strand...
#> Calculating balance score...
#> Slice-reduce to find bidirectional clusters...
#> Calculating statistics...
#> Preparing output...
#> # Bidirectional clustering summary:
#> Number of bidirectional clusters: 106779
#> Maximum balance score: 1
#> Minimum balance score: 0.950001090872574
#> Maximum width: 1866
#> Minimum width: 401
#>  - Running subsetByBidirectionality:
#> Calculating bidirectionality...
#> Subsetting...
#> Removed 73250 out of 106779 regions (68.6%)
#>  - Running quantifyClusters:
#> Finding overlaps...
#> Aggregating within clusters...

Note: quickEnhancers runs CAGEfightR with default settings. If you have larger or more noisy datasets you most likely want to do a more robust analysis with different settings. See the CAGEfightR vignette for more information.

Again, we are not usually interested in very lowly expressed BCs. As BCs are overall lowly expressed compared to TCs, we will simply filter out BCs without at least 1 count in 5 samples:

BCs <- subsetBySupport(BCs, 
                       inputAssay = "counts", 
                       unexpressed = 0, 
                       minSamples = 4)
#> Calculating support...
#> Subsetting...
#> Removed 20017 out of 33529 regions (59.7%)

3.1.4 Annotating clusters with transcript models

After having located unidirectional and bidirectional clusters, we can annotate them according to known transcript and gene models. These types of annotation are store via TxDb-objects in Bioconductor. Here we will simply use UCSC transcripts included in the TxDb.Mmusculus.UCSC.mm9.knownGene package, but the CAGEfightR vignette includes examples of how to obtain a TxDb object from other sources (GFF/GTF files, AnnotationHub, etc.).

Starting with the TSS candidates, we can not only annotate a TSS with overlapping transcripts, but also in what part of a transcript a TSS lies by using a hierarchical annotation scheme. As some TSS candidates might be very wide, we only use the TSS peak for annotation purposes:

# Annotate with transcript IDs
TSSs <- assignTxID(TSSs, txModels = txdb, swap = "thick")
#> Extracting transcripts...
#> Finding hierachical overlaps...
#> ### Overlap Summary: ###
#> Features overlapping transcripts: 87.65 %
#> Number of unique transcripts: 31898

# Annotate with transcript context
TSSs <- assignTxType(TSSs, txModels = txdb, swap = "thick")
#> Finding hierachical overlaps with swapped ranges...
#> ### Overlap summary: ###
#>       txType count percentage
#> 1   promoter 13395       46.4
#> 2   proximal  2246        7.8
#> 3    fiveUTR  2112        7.3
#> 4   threeUTR  1200        4.2
#> 5        CDS  3356       11.6
#> 6       exon   161        0.6
#> 7     intron  2810        9.7
#> 8  antisense  1294        4.5
#> 9 intergenic  2311        8.0

Almost half of the TSSs were found at annotated promoters, while the other half is novel compared to the UCSC known transcripts.

Transcript annotation is particularly useful for enhancer analysis, as bidirectional clustering might also detect bidirectional promoters. Therefore, a commonly used filtering approached is to only consider BCs in intergenic or intronic regions as enhancer candidates:

# Annotate with transcript context
BCs <- assignTxType(BCs, txModels = txdb, swap = "thick")
#> Finding hierachical overlaps with swapped ranges...
#> ### Overlap summary: ###
#>       txType count percentage
#> 1   promoter   766        5.7
#> 2   proximal  1649       12.2
#> 3    fiveUTR    67        0.5
#> 4   threeUTR   596        4.4
#> 5        CDS   420        3.1
#> 6       exon    71        0.5
#> 7     intron  6815       50.4
#> 8  antisense     0        0.0
#> 9 intergenic  3128       23.1

# Keep only non-exonic BCs as enhancer candidates
Enhancers <- subset(BCs, txType %in% c("intergenic", "intron"))

This leaves almost 10000 BCs for analysis. Again, for simplificity, we will refer to these non-exonic BCs as enhancer candidates for the remainder of the workflow.

3.1.5 Merging into a single dataset

For many downstream analyses, in particular normalization and differential expression, it is useful to combine both TSS and enhancers candidates into a single dataset. This ensures that clusters do not overlap, so that each CAGE tag is counted only once.

We must first ensure that the enhancer and TSS candidates have the same information attached to them, since CAGEfightR will only allow merging of clusters if they have the same sample and cluster information:

# Clean colData
TSSs$totalTags <- NULL
Enhancers$totalTags <- NULL

# Clean rowData
rowData(TSSs)$balance <- NA
rowData(TSSs)$bidirectionality <- NA
rowData(Enhancers)$txID <- NA

# Add labels for making later retrieval easy
rowData(TSSs)$clusterType <- "TSS"
rowData(Enhancers)$clusterType <- "Enhancer"

Then the clusters can be merged: As enhancers could technically be detected as two separate TSSs, we only keep the enhancer if a TSS and enhancer candidate overlaps:

RSE <- combineClusters(object1 = TSSs, 
                       object2 = Enhancers, 
                       removeIfOverlapping = "object1")
#> Removing overlapping features from object1: 374
#> Keeping assays: counts
#> Keeping columns: score, thick, support, txID, txType, balance, bidirectionality, clusterType
#> Merging metadata...
#> Stacking and re-sorting...

We finally calculate the total number of tags and TPM-scaled counts for the final merged dataset:

RSE <- calcTPM(RSE)
#> Calculating library sizes...
#> Calculating TPM...

3.2 Part 2: Genomic analysis of TSSs and enhancers

3.2.1 Genome browser-figures of TSSs and enhancers

First we can simply plot some examples of TSSs and enhancers in a genome browser-style figure using the Gviz package (Hahne and Ivanek 2016). It takes a bit of code to setup, but the resulting tracks can be reused for later examples:

# Genome track
axis_track <- GenomeAxisTrack()

# Annotation track
tx_track <- GeneRegionTrack(txdb, 
                            name = "Gene Models", 
                            col = NA,
                            fill = "bisque4", 
                            shape = "arrow", 
                            showId = TRUE)

A good general strategy for quickly generating genome browser plots is to first define a region of interest, and then only plotting data within that region using subsetByOverlaps. The following code demonstrates this using the first TSS candidate:

# Extract 100 bp around the first TSS
plot_region <- RSE %>% 
    rowRanges() %>% 
    subset(clusterType == "TSS") %>% 
    .[1] %>%
    add(100) %>%

# CTSS track
ctss_track <- CTSSs %>%
    rowRanges() %>%
    subsetByOverlaps(plot_region) %>%
    trackCTSS(name = "CTSSs")
#> Splitting pooled signal by strand...
#> Preparing track...

# Cluster track
cluster_track <- RSE %>%
    subsetByOverlaps(plot_region) %>%
    trackClusters(name = "Clusters", 
                  col = NA, 
                  showId = TRUE)
#> Setting thick and thin features...
#> Merging and sorting...
#> Preparing track...

# Plot tracks together
           from = start(plot_region), 
           to = end(plot_region), 
           chromosome = as.character(seqnames(plot_region)))
Genome browser example of TSS candidate

Figure 1: Genome browser example of TSS candidate

The top track shows the pooled CTSS signal and the middle track shows the identified TSS candidate. The thick bar indicates the TSS candidate peak (the overall most used CTSSs within the TSS candidate). The bottom track shows the known transcript model at this genomic location. In this case, the CAGE-defined TSS candidate corresponds well to the annotation.

We can also plot the first enhancer candidate:

# Extract 100 bp around the first enhancer
plot_region <- RSE %>% 
    rowRanges() %>% 
    subset(clusterType == "Enhancer") %>% 
    .[1] %>%
    add(100) %>%

# CTSS track
ctss_track <- CTSSs %>%
    rowRanges() %>%
    subsetByOverlaps(plot_region) %>%
    trackCTSS(name = "CTSSs")
#> Splitting pooled signal by strand...
#> Preparing track...

# Cluster track
cluster_track <- RSE %>%
    rowRanges %>%
    subsetByOverlaps(plot_region) %>%
    trackClusters(name = "Clusters", 
                  col = NA, 
                  showId = TRUE)
#> Setting thick and thin features...
#> Merging and sorting...
#> Preparing track...

# Plot tracks together
           from = start(plot_region), 
           to = end(plot_region), 
           chromosome = as.character(seqnames(plot_region)))
Genome browser example of enhancer candidate

Figure 2: Genome browser example of enhancer candidate

Here we see the bidirectional pattern characteristic of active enhancers. The enhancer candidate is seen in the middle track. The midpoint in thick marks the maximally balanced point within the enhancer candidate.

3.2.2 Location and expression of TSSs and enhancers

In addition to looking at single examples of TSS and enhancer candidates, we also want to get an overview of the number and expression of clusters in relation to transcript annotation. First we extract all of the necessary data from the RangedSummarizedExperiment into an ordinary data.frame:

cluster_info <- RSE %>%
    rowData() %>%

Then we use ggplot2 to plot the number and expression levels of clusters in each annotation category:

# Number of clusters
ggplot(cluster_info, aes(x = txType, fill = clusterType)) +
    geom_bar(alpha = 0.75, position = "dodge", color = "black") +
    scale_fill_colorblind("Cluster type") +
    labs(x = "Cluster annotation", y = "Frequency") +
    theme(axis.text.x = element_text(angle = 90, hjust = 1))
Number of clusters within each annotation category

Figure 3: Number of clusters within each annotation category

# Expression of clusters
ggplot(cluster_info, aes(x = txType, 
                         y = log2(score / ncol(RSE)), 
                         fill = clusterType)) +
    geom_violin(alpha = 0.75, draw_quantiles = c(0.25, 0.50, 0.75)) +
    scale_fill_colorblind("Cluster type") +
    labs(x = "Cluster annotation", y = "log2(TPM)") +
    theme(axis.text.x = element_text(angle = 90, hjust = 1))
#> Warning in regularize.values(x, y, ties, missing(ties)): collapsing to unique
#> 'x' values
Expression of clusters within each annotation category

Figure 4: Expression of clusters within each annotation category

We find that TSS candidates at annotated promoters are generally highly expressed. Most novel TSS candidates are expressed at lower levels, except for some TSS candidates in 5’-UTRs. Enhancer candidates are expressed at much lower levels than TSS candidates.

3.2.3 Analysing TSS shapes and sequences

A classic analysis of CAGE data is to divide TSSs into Sharp and Broad classes, which show different core promoter regions and different expression patterns across tissues(Carninci et al. 2006).

CAGEfightR can calculate several shape statistics that summarize the shape of a TSS. The Interquantile Range (IQR) can be used to find sharp and broad TSSs. The IQR measures the bp-distance holding e.g. 10-90% of the pooled expression in a TSS candidate, which dampens the effect of possible straggler tags that can greatly extend the width of a TSS candidate without contributing much to expression. As lowly expressed TSS candidates cannot show much variation in shape due to their low width and number of tags, we here focus on highly expressed TSS candidates (average TPM >= 10):

# Select highly expressed TSSs
highTSSs <- subset(RSE, clusterType == 'TSS' & score / ncol(RSE) >= 10)

# Calculate IQR as 10%-90% interval 
highTSSs <- calcShape(highTSSs, 
                      pooled = CTSSs, 
                      shapeFunction = shapeIQR, 
                      lower = 0.10, 
                      upper = 0.90)
#> Splitting by strand...
#> Applying function to each cluster...
#> Preparing output output...

We can then plot the bimodal distribution of IQRs. We use a zoom-in panel to highlight the distinction between the two classes:

highTSSs %>%
    rowData() %>% %>%
    ggplot(aes(x = IQR)) +
    geom_histogram(binwidth = 1, 
                   fill = "hotpink", 
                   alpha = 0.75) +
    geom_vline(xintercept = 10, 
               linetype = "dashed", 
               alpha = 0.75, 
               color = "black") +
    facet_zoom(xlim = c(0,100)) +
    labs(x = "10-90% IQR", 
         y = "Frequency")
Bimodal distribution of Interquartile Ranges (IQRs) of highly expressed TSSs

Figure 5: Bimodal distribution of Interquartile Ranges (IQRs) of highly expressed TSSs

We see most TSS candidates are either below or above 10 bp IQR (dashed line), so we use this cutoff to classify TSS candidates into Sharp and Broad:

# Divide into groups
rowData(highTSSs)$shape <- ifelse(rowData(highTSSs)$IQR < 10, "Sharp", "Broad")

# Count group sizes
#> Broad Sharp 
#>  9555   812

We can now investigate the core promoter sequences of the two classes of TSS candidates. We first need to extract the surrounding promoter sequence for each TSS candidate: We define this as the TSS candidate peak -40/+10 bp and extract them using the BSgenome.Mmusculus.UCSC.mm9 genome package:

promoter_seqs <- highTSSs %>%
    rowRanges() %>%
    swapRanges() %>%
    promoters(upstream = 40, downstream = 10) %>%
    getSeq(bsg, .)

This returns a DNAStringSet-object which we can plot as a sequence logo (Schneider and Stephens 1990) via the ggseqlogo package(Wagih 2017):

promoter_seqs %>%
    as.character %>%
    split(rowData(highTSSs)$shape) %>%
    ggseqlogo(data = ., ncol = 2, nrow = 1) +
    theme_logo() +
    theme(axis.title.x = element_blank(),
          axis.text.x = element_blank(),
          axis.ticks.x = element_blank())
Sequence logos of core promoter regions of Sharp and Broad classes of TSSs

Figure 6: Sequence logos of core promoter regions of Sharp and Broad classes of TSSs

As expected, we observe that Sharp TSS candidates tend to have a stronger TATA-box upstream of the TSS peak compared to Broad TSS candidates.

3.2.4 Finding candidates for interacting TSSs and enhancers

In addition to simply identifying enhancers, it is also interesting to try and infer what genes they might be regulating. CAGE data can itself not provide direct evidence that an enhancer is physically interacting with a TSS. This would require specialized chromatin confirmation capture assays such as HiC, 4C, 5C, etc. However, previous studies have shown that TSSs and enhancers that are close to each other and have highly correlated expression are more likely to be interacting. We can therefore use distance and correlation of expression between TSSs and enhancers to identify TSSs-enhancer links as candidates for physical interactions(Andersson et al. 2014).

To do this with CAGEfightR, we first need to indicate the two types of clusters as a factor with two levels:

rowData(RSE)$clusterType <- RSE %>%
    rowData() %>%
    use_series("clusterType") %>%
    as_factor() %>%

We can then calculate all pairwise correlations between TSSs and enhancer within a distance of 50 kbp. Here we use the non-parametric Kendall’s tau as a measure of correlation, but other functions for calculating correlation can be supplied (e.g. one could calculate Pearson’s r on log-transformed TPM values to only capture linear relationships):

# Find links and calculate correlations
all_links <- RSE %>%
    swapRanges() %>%
    findLinks(maxDist = 5e4L,
              directional = "clusterType",
              inputAssay = "TPM",
              method = "kendall")
#> Finding directional links from TSS to Enhancer...
#> Calculating 41311 pairwise correlations...
#> Preparing output...
#> # Link summary:
#> Number of links: 41311
#> Summary of pairwise distance:
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
#>     205    8832   21307   22341   35060   50000

# Show links
#> GInteractions object with 41311 interactions and 4 metadata columns:
#>              seqnames1   ranges1        seqnames2   ranges2 | orientation
#>                  <Rle> <IRanges>            <Rle> <IRanges> | <character>
#>       [1]         chr1   6204746 ---         chr1   6226837 |  downstream
#>       [2]         chr1   7079251 ---         chr1   7083527 |  downstream
#>       [3]         chr1   9535519 ---         chr1   9554735 |  downstream
#>       [4]         chr1   9538162 ---         chr1   9554735 |  downstream
#>       [5]         chr1  20941781 ---         chr1  20990601 |  downstream
#>       ...          ...       ... ...          ...       ... .         ...
#>   [41307]  chr9_random    193165 ---  chr9_random    217926 |    upstream
#>   [41308]  chr9_random    193165 ---  chr9_random    242951 |    upstream
#>   [41309]  chr9_random    223641 ---  chr9_random    217926 |  downstream
#>   [41310]  chr9_random    223641 ---  chr9_random    242951 |    upstream
#>   [41311] chrUn_random   3714359 --- chrUn_random   3718258 |    upstream
#>            distance            estimate            p.value
#>           <integer>           <numeric>          <numeric>
#>       [1]     22090 -0.0603022689155527  0.805433562909099
#>       [2]      4275   0.365994211051474  0.128612838399956
#>       [3]     19215   -0.21320071635561  0.392330339776564
#>       [4]     16572   0.341121146168977   0.17111237306132
#>       [5]     48819    0.14070529413629  0.565460671338501
#>       ...       ...                 ...                ...
#>   [41307]     24760   0.477084298221423 0.0423302291213607
#>   [41308]     49785   0.180906806746658  0.459929012970529
#>   [41309]      5714 -0.0366987921708787  0.875896057922941
#>   [41310]     19309  -0.261309831967395   0.28579482541369
#>   [41311]      3898  -0.170560573084488  0.493773664508106
#>   -------
#>   regions: 38454 ranges and 8 metadata columns
#>   seqinfo: 35 sequences (1 circular) from mm9 genome

The output is a GInteractions-object from the InteractionSet package(Lun, Perry, and Ing-Simmons 2016): For each TSS-enhancer link, both the distance and orientation (upstream/downstream relative to TSS) is calculated in addition to the correlation estimate and p-value. If one were to extract a set of highly correlated links for further analysis, it would be appropriate to correct the p-values for multiple testing, e.g. with the p.adjust(). For now, we are only interested in the top positive correlations, so we subset and sort the links:

# Subset to only positive correlation
cor_links <- subset(all_links, estimate > 0)

# Sort based on correlation
cor_links <- cor_links[order(cor_links$estimate, decreasing = TRUE)]

We can then visualize the correlation patterns across a genomic region, here using the most correlated TSS-enhancer link:

# Extract region around the link of interest
plot_region <- cor_links[1] %>% 
    boundingBox() %>% 
    linearize(1:2) %>%

# Cluster track
cluster_track <- RSE %>%
    subsetByOverlaps(plot_region) %>%
    trackClusters(name = "Clusters", 
                  col = NA, 
                  showId = TRUE)
#> Setting thick and thin features...
#> Merging and sorting...
#> Preparing track...

# Link track
link_track <- cor_links %>%
    subsetByOverlaps(plot_region) %>%
    trackLinks(name = "Links",
               interaction.measure = "p.value",
               interaction.dimension.transform = "log",
               col.outside = "grey",
               plot.anchors = FALSE,
               col.interactions = "black")

# Plot tracks together
           from = start(plot_region), 
           to = end(plot_region), 
           chromosome = as.character(seqnames(plot_region)))
Genome browser example of TSS-enhancer link

Figure 7: Genome browser example of TSS-enhancer link

The top track shows the correlation between 3 TSS candidates around the Atp1b1 gene. The most significant correlation is seen between the upstream TSS candidate and the most distal enhancer candidate.

A word of caution on calculating correlations between TSSs and enhancers in this manner: As we are calculating the correlation of expression across biological replicates from two conditions (Ctrl and Nano), high correlations could also arise from a TSS and enhancer candidate responding in the same direction in response to the treatment. This means that the correlation observed when combining all samples across conditions can be different from the correlation calculated within each condition (This unintuitive phenomenon is known as Simpson’s Paradox). To avoid including such cases, one could analyze each condition separately to find TSS-enhancer links within each state. As an extension of this approach, one could also look at TSS-enhancer links that show different strengths of correlation in different states. Analyses of this type are referred to as differential coexpression analysis.

3.2.5 Finding stretches of enhancers

Several studies have found that groups or stretches of closely spaced enhancers tend to show different chromatin characteristics and functions compared to singleton enhancers. Such groups of enhancers are often referred to as “super enhancers” or “stretch enhancers”(Pott and Lieb 2015).

CAGEfightR can detect such enhancer stretches based on CAGE data. CAGEfightR groups nearby enhancers and calculates the average pairwise correlation between them, shown below (again using Kendall’s tau):

# Subset to only enhancers
Enhancers <- subset(RSE, clusterType == "Enhancer")

# Find stretches within 12.5 kbp
stretches <- findStretches(Enhancers, 
                           inputAssay = "TPM",
                           mergeDist = 12500L,
                           minSize = 5L,
                           method = "kendall")
#> Finding stretches...
#> Calculating correlations...
#> # Stretch summary:
#> Number of stretches: 95
#> Total number of clusters inside stretches: 587 / 9943
#> Minimum clusters: 5
#> Maximum clusters: 15
#> Minimum width: 7147
#> Maximum width: 92531
#> Summary of average pairwise correlations:
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
#> -0.10038  0.01351  0.08107  0.09097  0.16171  0.37105

Similarly to TSSs and enhancers, we can also annotate stretches based on their relation with known transcripts:

# Annotate transcript models
stretches <- assignTxType(stretches, txModels = txdb)
#> Finding hierachical overlaps...
#> ### Overlap summary: ###
#>       txType count percentage
#> 1   promoter    50       52.6
#> 2   proximal     0        0.0
#> 3    fiveUTR     6        6.3
#> 4   threeUTR     5        5.3
#> 5        CDS     3        3.2
#> 6       exon     2        2.1
#> 7     intron    15       15.8
#> 8  antisense     0        0.0
#> 9 intergenic    14       14.7

# Sort by correlation
stretches <- stretches[order(stretches$aveCor, decreasing = TRUE)]

# Show the results
#> GRanges object with 95 ranges and 4 metadata columns:
#>                             seqnames              ranges strand |
#>                                <Rle>           <IRanges>  <Rle> |
#>     chr11:98628005-98647506    chr11   98628005-98647506      * |
#>    chr7:139979437-140003112     chr7 139979437-140003112      * |
#>     chr15:31261340-31279984    chr15   31261340-31279984      * |
#>   chr11:117733009-117752208    chr11 117733009-117752208      * |
#>      chr7:97167988-97188451     chr7   97167988-97188451      * |
#>                         ...      ...                 ...    ... .
#>   chr15:101076561-101093429    chr15 101076561-101093429      * |
#>     chr16:91373912-91399202    chr16   91373912-91399202      * |
#>    chr7:132619265-132644381     chr7 132619265-132644381      * |
#>     chr15:79181690-79208915    chr15   79181690-79208915      * |
#>     chr10:94708643-94729408    chr10   94708643-94729408      * |
#>                                         revmap nClusters              aveCor
#>                                  <IntegerList> <integer>           <numeric>
#>     chr11:98628005-98647506 6600,6601,6602,...         6   0.371052840516797
#>    chr7:139979437-140003112 4220,4221,4222,...         5   0.328630841442886
#>     chr15:31261340-31279984 7962,7963,7964,...         5   0.301603791540209
#>   chr11:117733009-117752208 6785,6786,6787,...         6   0.284399425439616
#>      chr7:97167988-97188451 4022,4023,4024,...         6   0.262199740521045
#>                         ...                ...       ...                 ...
#>   chr15:101076561-101093429 8320,8321,8322,...         5 -0.0549688493223916
#>     chr16:91373912-91399202 8643,8644,8645,...         7 -0.0598361076236999
#>    chr7:132619265-132644381 4160,4161,4162,...         5 -0.0626248504104628
#>     chr15:79181690-79208915 8144,8145,8146,...         5 -0.0981772309926707
#>     chr10:94708643-94729408 5823,5824,5825,...         5  -0.100380656957041
#>                                 txType
#>                               <factor>
#>     chr11:98628005-98647506   promoter
#>    chr7:139979437-140003112   promoter
#>     chr15:31261340-31279984     intron
#>   chr11:117733009-117752208   promoter
#>      chr7:97167988-97188451   promoter
#>                         ...        ...
#>   chr15:101076561-101093429 intergenic
#>     chr16:91373912-91399202    fiveUTR
#>    chr7:132619265-132644381   promoter
#>     chr15:79181690-79208915   promoter
#>     chr10:94708643-94729408     intron
#>   -------
#>   seqinfo: 35 sequences (1 circular) from mm9 genome

The returned GRanges contains the the location, number of enhancers and average correlation for each stretch. Stretches are found in a variety of context, some being intergenic and others spanning various parts of genes. Let us plot one of the top intergenic stretches:

# Extract region around a stretch of enhancers
plot_region <- stretches["chr17:26666593-26675486"] + 1000

# Cluster track
cluster_track <- RSE %>%
    subsetByOverlaps(plot_region) %>%
    trackClusters(name = "Clusters", 
                  col = NA, 
                  showId = TRUE)
#> Setting thick and thin features...
#> Merging and sorting...
#> Preparing track...

# CTSS track
ctss_track <- CTSSs %>%
    subsetByOverlaps(plot_region) %>%
    trackCTSS(name = "CTSSs")
#> Splitting pooled signal by strand...
#> Preparing track...

# Stretch enhancer track
stretch_track <- stretches %>%
    subsetByOverlaps(plot_region) %>%
    AnnotationTrack(name = "Stretches", fill = "hotpink", col = NULL)

# Plot tracks together
           from = start(plot_region), 
           to = end(plot_region), 
           chromosome = as.character(seqnames(plot_region)))
Genome browser example of enhancer stretch

Figure 8: Genome browser example of enhancer stretch

This stretch is composed of at least 5 enhancer candidates, each of which shows bidirectional transcription.

3.3 Part 3: Differential Expression analysis of TSSs, enhancers and genes

3.3.1 Normalization of expression and EDA

Before performing statistical tests for various measures of Differential Expression (DE), it is important to first conduct a thorough Exploratory Data Analysis (EDA) to identify what factors we need to include in the final DE model.

Here we will use DESeq2 (Love, Huber, and Anders 2014) for normalization and EDA since it offers easy to use functions for performing basic analyses. Other popular tools such as edgeR (Robinson, McCarthy, and Smyth 2010) and limma (Ritchie et al. 2015) offer similar functionality, as well as more specialized packages for EDA such as EDASeq.

DESeq2 offers sophisticated normalization and transformation of count data in the form of the variance stabilized transformation: this adds a dynamic pseudo-count to normalized expression values before log transforming to dampen the inherent mean-variance relationship of count data. This is particularly useful for CAGE data, as CAGE can detect even very lowly expressed TSSs and enhancers.

Note: Due to their overall lower expression, enhancer candidate tags make up only a small proportion of the total number of tags. As proper estimation of normalization factors assumes a large number non-DE features, both TSS and enhancer candidates should normally be included in a DE analysis.

First, we fit a “blind” version of the variance stabilizing transformation, since we do not yet know what design is appropriate for this particular study:

# Create DESeq2 object with blank design
dds_blind <- DESeqDataSet(RSE, design = ~ 1)

# Normalize and log transform
vst_blind <- vst(dds_blind, blind = TRUE)

A very useful first representation is a Principal Components Analysis (PCA) plot summarizing variance across the entire experiment as Principle Components (PCs):

plotPCA(vst_blind, "Class")
PCA-plot of variance stabilized expression.

Figure 9: PCA-plot of variance stabilized expression

We observe that PC2 separates the samples according to the experimental group (Nano vs Ctrl). However, PC1 also separates samples into two groups. This is suggestive of an unwanted yet systematic effect on expression, often referred as a batch effect. Batch effect can arise for a multitude of reasons, e.g. libraries being prepared by different labs or people or using slightly different reagent pools. Often, batch effects co-ocur with the date libraries are prepared, and indeed Bornholdt et al suggests this as the cause of the batch effect in the original study.

We do not want to mistake this unwanted variation for biological variation when we test for DE. To prevent this, we can include the batch effect as a factor in the final DE model. First, we define the batch variable:

# Extract PCA results
pca_res <- plotPCA(vst_blind, "Class", returnData = TRUE)

# Define a new variable using PC1
batch_var <- ifelse(pca_res$PC1 > 0, "A", "B")

# Attach the batch variable as a factor to the experiment
RSE$Batch <- factor(batch_var)

# Show the new design
RSE %>%
    colData() %>%
    subset(select = c(Class, Batch)) %>%
    kable(caption = "Design matrix after adding new batch covariate.") %>%
    kable_styling(latex_options = "hold_position")
Table 4: Design matrix after adding new batch covariate.
Class Batch
C547 Ctrl B
C548 Ctrl B
C549 Ctrl B
C559 Ctrl A
C560 Ctrl A
N13 Nano B
N14 Nano A
N15 Nano B
N16 Nano A
N17 Nano A
N18 Nano A

As an alternative to manually defining the batch variable, tools such as sva and RUVSeq can be used to estimate unknown batch effects from the data.

3.3.2 Cluster-level differential expression

Following our short EDA above, we are ready to specify the final design for the experiment: We want to take into account both the Class and Batch of samples:

# Specify design
dds <- DESeqDataSet(RSE, design = ~ Batch + Class)

# Fit DESeq2 model
dds <- DESeq(dds)
#> estimating size factors
#> estimating dispersions
#> gene-wise dispersion estimates
#> mean-dispersion relationship
#> final dispersion estimates
#> fitting model and testing

We can now extract estimated effects (log fold changes) and statistical significance (p-values) for the Nano-vs-Ctrl comparison, implicitly correcting for the batch effect:

# Extract results
res <- results(dds,
               contrast = c("Class", "Nano", "Ctrl"),
               alpha = 0.05, 
               independentFiltering = TRUE, 
               tidy = TRUE) %>%
    bind_cols( %>%

# Show the top hits
res %>%
    top_n(n = -10, wt = padj) %>%
    dplyr::select(cluster = row, 
                  padj) %>%
    kable(caption = "Top differentially expressed TSS and enhancer candidates") %>%
    kable_styling(latex_options = "hold_position")
Table 5: Top differentially expressed TSS and enhancer candidates
cluster clusterType txType baseMean log2FoldChange padj
chr1:73977049-73977548;- TSS intron 1183.3740 2.838367 0
chr2:32243097-32243468;- TSS promoter 30799.5953 3.741789 0
chr3:144423689-144423778;- TSS promoter 191.0431 3.709530 0
chr4:125840648-125840820;- TSS proximal 1063.4328 3.867574 0
chr4:137325466-137325712;- TSS intron 176.7636 3.912592 0
chr7:53971039-53971170;- TSS promoter 8720.5204 6.696838 0
chr9:120212846-120213294;+ TSS promoter 316.0582 2.404706 0
chr11:83222553-83222887;+ TSS proximal 228.5560 6.098838 0
chr12:105649334-105649472;+ TSS CDS 175.1364 3.345411 0
chr19:56668148-56668332;+ TSS CDS 103.8795 -2.254371 0

It always a good idea to inspect a few diagnostics plots to make sure the DESeq2 analysis was successful. One such example is an MA-plot (another useful plot is the p-value histogram):

ggplot(res, aes(x = log2(baseMean), 
                y = log2FoldChange, 
                color = padj < 0.05)) +
    geom_point(alpha = 0.25) +
    geom_hline(yintercept = 0, 
               linetype = "dashed", 
               alpha = 0.75) +
    facet_grid(clusterType ~ .)
Diagnostic MA plot of the differential expression analysis

Figure 10: Diagnostic MA plot of the differential expression analysis

We can see that we overall find more DE TSS compared to enhancer candidates, which is expected since TSS candidates are more highly expressed. Many enhancer candidates are filtered away for the final DESeq2 analysis (The “Independent Filtering” Step), as their expression level is too low to detect any DE: This increases power for detecting DE at higher expression levels for the remaining TSS and enhancer candidates.

We can tabulate the total number of DE TSS and enhancer candidates:

table(clusterType = rowRanges(RSE)$clusterType, 
      DE = res$padj < 0.05)
#>            DE
#> clusterType FALSE  TRUE
#>    TSS      22071  6385
#>    Enhancer  3034   199

3.3.3 Correcting expression estimates for batch effects

In addition to looking at estimates and significance for each cluster, we might also want to look at individual expression values for some top hits. However, we then need to also correct the expression estimates themselves for batch effects, just like we did for log fold changes and p-values (using the same model of course).

Here we use ComBat(Johnson, Li, and Rabinovic 2007) from the sva package which is suitable for removing simple batch effects from small experiments. For more advanced setups, removeBatchEffect from limma can remove arbitrarily complex batch effects. The RUVSeq package and fsva from sva can be used to remove unknown batch effects.

We again use the variance stabilizing transformation to prepare the data for ComBat (this makes count data resemble expression estimates obtained from microarrays, as ComBat was originally developed for microarrays):

# Guided / non-blind  variance stabilizing transformation
vst_guided <- varianceStabilizingTransformation(dds, blind = FALSE)

To run ComBat we need two additional pieces of information: i) A design matrix describing the biological or wanted effects and ii) the known but unwanted batch effect. We first specify the design matrix, and then run ComBat:

# Design matrix of wanted effects
bio_effects <- model.matrix(~ Class, data = colData(RSE))

# Run ComBat
assay(RSE, "ComBat") <- ComBat(dat = assay(vst_guided),
                                    batch = RSE$Batch,
                                    mod = bio_effects)
#> Found2batches
#> Adjusting for1covariate(s) or covariate level(s)
#> Standardizing Data across genes
#> Fitting L/S model and finding priors
#> Finding parametric adjustments
#> Adjusting the Data

We can redo the PCA-plot, to see the global effect of the batch effect correction:

# Overwrite assay 
assay(vst_guided) <- assay(RSE, "ComBat")

# Plot as before
plotPCA(vst_guided, "Class")
PCA-plot of batch corrected expression.

Figure 11: PCA-plot of batch corrected expression

Now Nano and Ctrl are separated along PC1 (compared to PC2 before correction). As PC1 captures the most variance, this indicates that the batch effect has been removed and that the experimental group is now the main contributor to variance of expression.

Then we extract the top 10 DE enhancer candidates using the following tidyverse code:

# Find top 10 DE enhancers
top10 <- res %>%
    filter(clusterType == "Enhancer", padj < 0.05) %>%
    group_by(log2FoldChange >= 0) %>%
    top_n(n = 5, wt = abs(log2FoldChange)) %>%

# Extract expression values in tidy format
tidyEnhancers <- assay(RSE, "ComBat")[top10, ] %>%
    t() %>%
    as_tibble(rownames = "Sample") %>%
    mutate(Class = RSE$Class) %>%
    gather(key = "Enhancer", 
           value = "Expression", 
           -Sample, -Class, 
           factor_key = TRUE)

Finally, we can plot the batch-corrected expression of each top enhancer candidate:

ggplot(tidyEnhancers, aes(x = Class, 
                          y = Expression, 
                          fill = Class)) +
    geom_dotplot(stackdir = "center", 
                 binaxis = "y", 
                 dotsize = 3) +
    facet_wrap(~ Enhancer, 
               ncol = 2, 
               scales = "free_y")
#> `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.