Contents

R version: R version 3.6.0 (2019-04-26)

Bioconductor version: 3.10

CAGEfightR version: 1.5.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 seen 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 mRNA from 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) and enhancers RNA (eRNA). It was previously shown that active enhancers are characterized by balanced bidirectional transcription, 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, Thieffry, Vitting-Seerup, et al. 2018), see 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 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, TSRchictet 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
TSS calling sliding window X-means distance or paraclu slice-reduce
TSS shapes - + + +
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 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 wholes 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:

library(nanotubes)

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 the results. All these packages are loaded prior to beginning the workflow:

# CRAN packages for data manipulation and plotting
library(knitr)
library(kableExtra)
library(pheatmap)
library(ggseqlogo)
library(viridis)
library(magrittr)
library(ggforce)
library(ggthemes)
library(tidyverse)

# CAGEfightR and related packages
library(CAGEfightR)
library(GenomicRanges)
library(SummarizedExperiment)
library(GenomicFeatures)
library(BiocParallel)
library(InteractionSet)
library(Gviz)

# Bioconductor packages for differential expression
library(DESeq2)
library(limma)
library(edgeR)
library(sva)

# Bioconductor packages for enrichment analyses
library(TFBSTools)
library(motifmatchr)
library(pathview)

# Bioconductor data packages
library(BSgenome.Mmusculus.UCSC.mm9)
library(TxDb.Mmusculus.UCSC.mm9.knownGene)
library(org.Mm.eg.db)
library(JASPAR2016)

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 <- org.Mm.eg.db

# 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 dusters 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: Tag cluster- 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

CAGEfightR starts analysis from CTSSs, which is the number of CAGE tag 5’-ends mapping to each basepair (bp) in the genome. CTSSs are normally stored as either BED-files or BigWig-files. CAGEfightR works on BigWig-files, since these can be efficiently imported and allow for random access.

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, sometimes called the design matrix. CAGEfightR can keep track of the design matrix throughout the analysis:

data(nanotubes)
kable(nanotubes, 
      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 mm9.CAGE_7J7P_NANO_KON_547.plus.bw mm9.CAGE_7J7P_NANO_KON_547.minus.bw
C548 Ctrl C548 mm9.CAGE_ULAC_NANO_KON_548.plus.bw mm9.CAGE_ULAC_NANO_KON_548.minus.bw
C549 Ctrl C549 mm9.CAGE_YM4F_Nano_KON_549.plus.bw mm9.CAGE_YM4F_Nano_KON_549.minus.bw
C559 Ctrl C559 mm9.CAGE_RSAM_NANO_559.plus.bw mm9.CAGE_RSAM_NANO_559.minus.bw
C560 Ctrl C560 mm9.CAGE_CCLF_NANO_560.plus.bw mm9.CAGE_CCLF_NANO_560.minus.bw
N13 Nano N13 mm9.CAGE_KTRA_Nano_13.plus.bw mm9.CAGE_KTRA_Nano_13.minus.bw
N14 Nano N14 mm9.CAGE_RSAM_NANO_14.plus.bw mm9.CAGE_RSAM_NANO_14.minus.bw
N15 Nano N15 mm9.CAGE_RFQS_Nano_15.plus.bw mm9.CAGE_RFQS_Nano_15.minus.bw
N16 Nano N16 mm9.CAGE_CCLF_NANO_16.plus.bw mm9.CAGE_CCLF_NANO_16.minus.bw
N17 Nano N17 mm9.CAGE_RSAM_NANO_17.plus.bw mm9.CAGE_RSAM_NANO_17.minus.bw
N18 Nano N18 mm9.CAGE_CCLF_NANO_18.plus.bw mm9.CAGE_CCLF_NANO_18.minus.bw

3.1.1 Importing CTSSs

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/Sample1_plus.bw), 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 step in a CAGEfightR analysis, but it can be speed 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: 282 MB

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

# Get a summary
CTSSs
#> 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
rowRanges(CTSSs)
#> GPos 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") %>%
    head
#> 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.2 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 CTSSs count in each sample to Tags-Per-Million (TPM) values, and them sum TPM values across samples:

CTSSs <- CTSSs %>%
    calcTPM() %>%
    calcPooled()
#> 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 only be very lowly expressed. To obtain likely biologically relevant TSSs, we keep only TSSs 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 TSSs candidates for analysis.

Then we turn to bidirectional clustering for identifying bidirectional clusters (BCs), which are candidate for enhancers. Similarly, we can use quickEnhancers to locate and quantify BCs:

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 they are overall lowly expressed, 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.3 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 TSSs 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 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 candidates, 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 enhancer candidates for analysis.

3.1.4 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 TSSs and enhancers do not overlap, so each CAGE tag is only counted 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 are the most complicated type, we keep only enhancers if a TSS and enhancer 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:

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

# CTSSs 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 at tracks together
plotTracks(list(axis_track, 
                ctss_track,
                cluster_track,
                tx_track),
           from = start(plot_region), 
           to=end(plot_region), 
           chromosome = 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 TC with the thick bar indicating the TSS peak (the overall most used CTSSs within the TC). The bottom track shows the known transcript model at this genomic location. In this case, the CAGE-defined TSS corresponds well to the annotation.

We can also plot the first enhancer:

# Make plotting region
plot_region <- RSE %>% 
    rowRanges %>% 
    subset(clusterType == "Enhancer") %>% 
    .[1] %>%
    add(100) %>%
    unstrand()

# CTSSs 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 at tracks together
plotTracks(list(axis_track, 
                ctss_track,
                cluster_track,
                tx_track),
           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 bidirectional cluster is seen in the middle track, with the midpoint in thick marking the maximally balanced point within the bidirectional cluster.

3.2.2 Location and expression of TSSs and enhancers

In addition to looking at single examples of TSSs and enhancers, 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() %>%
    as.data.frame()

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 TSSs at annotated promoters are generally highly expressed. Most novel TSSs are expressed at lower levels, except for some TSSs in 5’-UTRs. Enhancers are expressed at much lower levels than TSSs.

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 summarizes the shape of a TSS. The Interquartile Range (IQR) can be used to find sharp and broad TSSs. As lowly expressed TSSs cannot show much variation in shape due to their low width and number of tags, we here focused on highly expressed TSSs (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 %>%
    as.data.frame %>%
    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 TSSs are either below or above 10 bp IQR (dashed line), so we use this cutoff to classify TSSs into Sharp and Broad:

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

# Count group sizes
table(rowData(highTSSs)$shape)
#> 
#> Broad Sharp 
#>  9555   812

We can now investigate the core promoters sequences of the two classes of TSSs. We first need to extract the sequences for each TSS: We define this as the TSS peak -40/+10 bp and extract them from using the BSgenome.Mmusculus.UCSC.mm10 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 TSSs tend to have a stronger TATA-box upstream of the TSS compared to Broad TSSs.

3.2.4 Finding candidates for interacting TSSs and enhancers

In addition to simply identifying enhancers, it is also interesting to try infer what genes they might be regulating. CAGE data can itself not provide direct evidence that an enhancer is physically interacting with a TSSs, which would requires 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() %>%
    fct_relevel("TSS")

We can then calculate all pairwise correlations between TSSs and enhancer within a distance of 50 bp. 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):

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
all_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 both the distance and orientation (upstream/downstream relative to TSS) is calculated in addition to the correlation estimate and p-value. For now, we are only interested in 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:

# Make plotting region
plot_region <- cor_links[1] %>% 
    anchors %>% 
    GRangesList() %>% 
    unlist() %>% 
    reduce(ignore.strand=TRUE, 
           min.gapwidth=1e5) %>%
    add(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...


# Cluster 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 at tracks together
plotTracks(list(axis_track, 
                link_track,
                cluster_track,
                tx_track),
           from = start(plot_region), 
           to=end(plot_region), 
           chromosome = as.character(seqnames(plot_region)))
Genome browser example of TSS-enhancer link candidates

Figure 7: Genome browser example of TSS-enhancer link candidates

The top track shows the strength of correlations between 3 TSSs around the Atp1b1 gene. The highest correlation is seen between the upstream TSS and the most distal enhancer.

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 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 into groups 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
stretches <- findStretches(Enhancers, 
                           inputAssay = "TPM",
                           mergeDist = 12500L,
                           minSize = 5,
                           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
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)]

# Inspect
stretches
#> 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 other spanning various parts of genes. Let us plot one of the top intergenic stretches:

# Make plotting region
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...

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

# Plot at tracks together
plotTracks(list(axis_track, 
                stretch_track,
                cluster_track,
                ctss_track),
           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 enhancers, 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 factor we need to include in the final 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.

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:

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 (control vs nano). 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. We do not want to mistake this unwanted variation for biological variation when we test for differential expression. To prevent this, we can include the batch information as a factor in the final model. Let first 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

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(as.data.frame(rowData(RSE))) %>%
    as_tibble

# Show the top hits
res %>%
    top_n(-10, padj) %>%
    dplyr::select(cluster=row, 
                  clusterType, 
                  txType, 
                  baseMean, 
                  log2FoldChange, 
                  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 plot to make sure the DESeq2 analysis was successful. One such example is an MA-plot (another useful plot is 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 differentially expressed TSSs compared to enhancers, which is expected since they are also more highly expressed. Many enhancers 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.

We can tabulate the total number of DE TSSs and enhancers:

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 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, # Unwanted 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

Let us 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 the first principal component (compared to the second principle component before correction).

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

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

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

Finally, we can plot the batch-corrected expression profiles of each individual enhancer:

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`.