A step-by-step guide to analyzing CAGE data using R/Bioconductor
25 October 2023
Abstract
Cap Analysis of Gene Expression (CAGE) is one of the most popular 5’-end sequencing methods. In a single experiment, CAGE can be used to locate and quantify the expression of both Transcription Start Sites (TSSs) and enhancers. This workflow is a case study on how to use the CAGEfightR package to orchestrate analysis of CAGE data within the Bioconductor project. This workflow starts from BigWig-files and covers both basic CAGE analyses such as identifying, quantifying and annotating TSSs and enhancers, advanced analysis such as finding interacting TSS-enhancer pairs and enhancer clusters, to differential expression analysis and alternative TSS usage. R-code, discussion and references are intertwined to help provide guidelines for future CAGE studies of the same kind.Contents
- 1 Background
- 2 Materials and methods
- 3 Workflow
- 4 Discussion
- 5 Author information
- 6 Competing interests
- 7 Grant information
- 8 Acknowledgments
R version: R version 4.3.1 (2023-06-16)
Bioconductor version: 3.18
CAGEfightR version: 1.22.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.
| Analysis | icetea | TSRchitect | CAGEr | CAGEfightR |
|---|---|---|---|---|
| Simplest input | FASTQ | BAM | BAM | BigWig/bedGraph |
| 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:
| Group | Biological Replicates | |
|---|---|---|
| Ctrl | 5 mice | |
| Nano | 6 mice |
The data is acquired via the nanotubes data package:
library(nanotubes)
data(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 results. All these packages are loaded prior to beginning the workflow:
# CRAN packages for data manipulation and plotting
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(statmod)
library(BiasedUrn)
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 figures3 Workflow
The workflow is divided into 3 parts covering different parts of a typical CAGE data analysis:
Shows how to use
CAGEfightRto import CTSSs and find and quantify TSS and enhancer candidates.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).
Shows how
CAGEfightRcan 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:
knitr::kable(nanotubes,
caption = "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 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/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$NameThe 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 design...
#> Checking supplied genome compatibility...
#> Iterating over 1 genomic tiles in 11 samples using 3 worker(s)...
#> Importing CTSSs from plus strand...
#> Importing CTSSs from minus strand...
#> Merging strands...
#> Formatting output...
#> ### CTSS summary ###
#> Number of samples: 11
#> Number of CTSSs: 9.339 millions
#> Sparsity: 81.68 %
#> Type of rowRanges: StitchedGPos
#> Final object size: 280.88 MBThe 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
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)
#> StitchedGPos 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.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() %>%
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 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 support...
#> 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...
#> 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%
#> 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.0Almost 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) %>%
unstrand()
# 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
plotTracks(list(axis_track,
ctss_track,
cluster_track,
tx_track),
from = start(plot_region),
to = end(plot_region),
chromosome = as.character(seqnames(plot_region)))
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) %>%
unstrand()
# 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
plotTracks(list(axis_track,
ctss_track,
cluster_track,
tx_track),
from = start(plot_region),
to = end(plot_region),
chromosome = as.character(seqnames(plot_region)))
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() %>%
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))
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))
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() %>%
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")
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
table(rowData(highTSSs)$shape)
#>
#> Broad Sharp
#> 9555 812We 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())
#> Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
#> of ggplot2 3.3.4.
#> ℹ The deprecated feature was likely used in the ggseqlogo package.
#> Please report the issue at <https://github.com/omarwagih/ggseqlogo/issues>.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.
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() %>%
fct_relevel("TSS")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
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.0603023 0.805434
#> [2] 4275 0.3659942 0.128613
#> [3] 19215 -0.2132007 0.392330
#> [4] 16572 0.3411211 0.171112
#> [5] 48819 0.1407053 0.565461
#> ... ... ... ...
#> [41307] 24760 0.4770843 0.0423302
#> [41308] 49785 0.1809068 0.4599290
#> [41309] 5714 -0.0366988 0.8758961
#> [41310] 19309 -0.2613098 0.2857948
#> [41311] 3898 -0.1705606 0.4937737
#> -------
#> regions: 38454 ranges and 8 metadata columns
#> seqinfo: 35 sequences (1 circular) from mm9 genomeThe 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) %>%
add(1000L)
# 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
plotTracks(list(axis_track,
link_track,
cluster_track,
tx_track),
from = start(plot_region),
to = end(plot_region),
chromosome = as.character(seqnames(plot_region)))
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.37105Similarly 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
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 txType
#> <IntegerList> <integer> <numeric> <factor>
#> chr11:98628005-98647506 6600,6601,6602,... 6 0.371053 promoter
#> chr7:139979437-140003112 4220,4221,4222,... 5 0.328631 promoter
#> chr15:31261340-31279984 7962,7963,7964,... 5 0.301604 intron
#> chr11:117733009-117752208 6785,6786,6787,... 6 0.284399 promoter
#> chr7:97167988-97188451 4022,4023,4024,... 6 0.262200 promoter
#> ... ... ... ... ...
#> chr15:101076561-101093429 8320,8321,8322,... 5 -0.0549688 intergenic
#> chr16:91373912-91399202 8643,8644,8645,... 7 -0.0598361 fiveUTR
#> chr7:132619265-132644381 4160,4161,4162,... 5 -0.0626249 promoter
#> chr15:79181690-79208915 8144,8145,8146,... 5 -0.0981772 promoter
#> chr10:94708643-94729408 5823,5824,5825,... 5 -0.1003807 intron
#> -------
#> seqinfo: 35 sequences (1 circular) from mm9 genomeThe 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
plotTracks(list(axis_track,
stretch_track,
cluster_track,
ctss_track),
from = start(plot_region),
to = end(plot_region),
chromosome = as.character(seqnames(plot_region)))
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")
#> using ntop=500 top features by variance
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)
#> using ntop=500 top features by variance
# 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)) %>%
knitr::kable(caption = "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 testingWe 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(n = -10, wt = padj) %>%
dplyr::select(cluster = row,
clusterType,
txType,
baseMean,
log2FoldChange,
padj) %>%
knitr::kable(caption = "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 ~ .)
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 1993.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)
#> Found 253 genes with uniform expression within a single batch (all zeros); these will not be adjusted for batch.
#> Found2batches
#> Adjusting for1covariate(s) or covariate level(s)
#> Standardizing Data across genes
#> Fitting L/S model and finding priors
#> Finding parametric adjustments
#> Adjusting the DataWe 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")
#> using ntop=500 top features by variance
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)) %>%
pull(row)
# 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")
#> Bin width defaults to 1/30 of the range of the data. Pick better value with
#> `binwidth`.