A step-by-step guide to analyzing CAGE data using R/Bioconductor.

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

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 [1– 4].

Cap Analysis of Gene Expression (CAGE) has emerged as one of the dominant high-throughput assays for studying TSSs [5]. 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 [6]. 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 [7, 8]. 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 [9].

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) [7, 10]. Since CAGE captures all capped RNAs, it can also identify long non-coding RNA (lincRNA) [11] 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 [12, 13]. 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 [14]. Only a few packages are dedicated to analyzing 5’-end data in general and CAGE data in particular: [15], [16], [17] and [18], see Table 1.

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	-	-	-	+

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 [19], 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 [10]. 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. [5]). 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 . Both CAGEr, TSRchictet and icetea offers built-in capabilities for differential expression (DE) analysis via the popular or packages [20, 21].

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 ( , , [22, 23]) 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.

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.

Methods

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 [24]. 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:

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:

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

Workflow

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

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

Shows examples of how to perform genomic analyses of CAGE dusters using core Bioconductor packages such as and . 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 CAGEfightR can be used to prepare data for differential expression analysis with popular R packages, including , and [20, 21, 25]. 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 [26]. Once differential expression results have been obtained, they can be combined with other sources of information such as motifs from JASPAR [27] and GO/KEGG terms [28, 29, 30].

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:

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:

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 genome package:

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:

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:

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:

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):

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:

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:

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 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:

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:

This leaves almost 10000 enhancer candidates for analysis.

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:

Then the clusters can be merged: As enhancers are the most complicated type, we keep only enhancers if a TSS and enhancer overlaps:

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

Part 2: Genomic analysis 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 [31]. It takes a bit of code to setup, but the resulting tracks can be reused for later examples:

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:

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:

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.

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:

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

We find that TSSs at annotated promoters are generally highly expressed. Most novel TSSs are expresse d at lower levels, except for some TSSs in 5’-UTRs. Enhancers are expressed at much lower levels than TSSs.

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 [7].

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):

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

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:

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 genome package:

This returns a DNAStringSet-object which we can plot as a sequence logo [32] via the package [33]:

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

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 [13].

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

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):

The output is a GInteractions-object from the package [23]: 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:

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

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.

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” [34].

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):

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

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:

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

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

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 [20] for normalization and EDA since it offers easy to use functions for performing basic analyses. Other popular tools such as [21] and [25] offer similar functionality, as well as more specialized packages for EDA such as .

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:

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

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:

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

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:

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

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):

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:

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 [35] from the 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 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).

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:

Let us redo the PCA-plot, to see the global effect of the batch effect correction:

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:

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

A typical question following identification of differentially expressed TSSs and enhancers, is what TFs might be involved in their regulation. To shed light on this question we can annotate TSSs and enhancers with DNA-binding motifs from the JASPAR database [27].

First we extract the sequences around TSSs and enhancers. Here we simply define it as +/- 500 bp around TSS peak or enhancer midpoint:

Secondly, we use the [36] package to obtain motifs as Position Frequency Matrices (PFMs) from the database:

Thirdly, we use the package [37] to find hits in the sequences:

Finally we can do a simple Fisher’s Exact test to see if a motif co-occurs more with DE TSSs and enhancer than we would expect be chance. Here we will look at the FOS::JUN motif (MA0099.2):

A significant odds ratio above 1 indicate that FOS::JUN is a candidate transcription factor (or, more technically correct, a candidate transcription factor dimer) in regulation of the nanotube response. This is not surprising given that FOS::JUN is part of the TNF-alpha inflammatory pathway (see more below).

Of course, this is a just a very quick and simple analysis of motif enrichment. One could easily have used different regions around TSSs and enhancers and/or split the enrichment analysis between TSSs and enhancers. Other Bioconductor packages like , and implements more advanced statistical methods for calculating enrichment of known motifs. , and can also be used to calculate enrichment of novel motifs, sometimes referred to as motif discovery.

While CAGE data is naturally analyzed at the level of clusters (TSSs and enhancers) it is in many cases interesting to also look at gene-level expression estimates. A prime example of this is looking at enrichment of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) terms [28, 29, 30] which are only defined at gene-level. CAGEfightR includes functions for annotating clusters with gene models and summarizing expression to gene-level.

We can annotate clusters with gene IDs in the same manner as Transcript IDs:

And then use CAGEfightR to sum counts of TSSs within genes:

The result is RangedSummarizedExperiment where the ranges are a GRangesList holding the TSSs that were summed within each gene:

The gene IDs in this case is Entrez ID (which is widely used by Bioconductor packages). We can translate these systematic IDs into more human-readable symbols using the annotation package:

Having obtained a gene-level count matrix we can now perform gene-level DE analysis. Here we use limma-voom, since limma makes it easy to perform a subsequent enrichment analysis. Other tools such as DESeq2 (above) or edgeR (see below) could also have been used.

Note: limma is a powerful tool for DE analysis of count-based data. However, since it depends on log transforming counts, it is not always suitable for analyzing datasets where features have very low counts. This is usually not a problem for gene-level analysis, but can be a problem for enhancers, which are generally very lowly expressed.

Similarly to the DESeq2 analysis, we first build the necessary object and then normalize the expression values:

Then we apply the voom-transformation to model the mean-variance trend, for which we also need to specify the design matrix (in this case the design must contain both wanted and unwanted effects!). The same design matrix is then used for fitting the gene-wise models:

We can the both report the overall summary of differential gene expression, and look at the first few top hits:

In addition to looking at individual top genes, we can look at how the differentially expressed genes relate to known databases of gene function to gain insight in what biological processes might be affected in the experiment.

limma makes it easy to perform such an enrichment analysis following a DE analysis. As we have gene indexed by Entrez IDs, we can directly use goana to find enriched GO-terms: goana uses a biased urn-model to estimate enrichment of GO-terms, while taking into account the expression levels of DE genes:

And similarly for KEGG terms:

Both analyses indicate that genes related to the inflammatory response and defense response are upregulated following nanotube exposure. This supports the hypothesis that nanotube induces a response similar to asbestos.

KEGG-terms represents well defined pathways. We can use the package [38] to investigate in more detail the genes in a given enriched pathway. For example, we can look at regulation of gene in the TNF- signalling pathway:

In the above two analyses we looked at whether an individual TSSs or an individual gene was changing expression between experimental groups. However, we might also want to look at whether a gene show differential TSS usage: whether a gene uses different TSSs under different conditions. This problem is similar to differential splicing in RNA-Seq, but looking at TSSs rather than isoforms [26]. Here we will use the edgeR diffSpliceDGE method to find differential TSS usage, although many other packages could have been used, for example diffSplice from limma, , , etc..

Intuitively, diffSpliceDGE tests whether a given TSSs show the same change as other TSSs in the same gene, indicating that TSSs are differentially regulated across the gene. This does however not take into account the relative composition of a given TSSs, e.g. whether a TSS increases from 1%–2% of gene output or 25%–50%. A useful preprocessing step is therefore to filter out TSSs making only a small contribution to total gene expression before analyses.

We use CAGEfightR to remove TSSs that are not expressed as more than 10% of total gene expression in more than 5 samples (We first remove TSSs not assigned to genes):

We can only do differential TSS usage analysis of genes with multiple TSSs. A useful first visualization is therefore to see how many genes use more than one TSS:

While most genes utilize only a single TSSs, many genes use two or more TSSs.

Again, we build the necessary R-objects for running edgeR:

Then we normalize and fit models using the Quasi-likelihood approach, including the diffSpliceDGE step:

Now we can look at differential TSS usage at two-levels: Whether an individual TSS shows differential TSS usage (TSS-level) or whether a gene show differential TSS usage in any way (gene-level). First we can look at individual TSSs (TSS-level differential TSS usage):

The interpretation of log fold changes here is slightly different from before: These log fold changes are relative to the overall log fold change for all TSSs in that gene.

Then we can look at results for each gene (Gene-level differential TSS usage):

We see that the two lists agree, which is not surprising given that the gene-level results are obtained by aggregating TSS-level p-values across genes.

We can look at closer at the TSS usage in on of the top hits: We can visualize the batch-corrected expression (See above) of each TSS in the Fblim1 gene via a heatmap:

Fblim1 has 3 TSSs, with 2 of them being used in the Ctrl samples, while the Nano samples also uses the chr4:141154044-141154185;- TSS, as also seen in the TSS-level table above. While a heatmap is useful for seeing expression changes, a genome browser view is better to inspect the genomic context of each TSSs:

The Fblim1 gene uses two annotated TSSs, but the Nano samples also uses a novel intronic TSS.

Discussion

This workflow is intended as providing an outline of the basic building blocks of CAGE data analysis, going from clustering, to spatial analyses to differential expression. More advanced analyses can be strung together from these basic elements: Finding enhancers linked to DE TSSs, enhancer stretches composed of DE enhancer, comparing DNA binding motif enrichments between DE enhancers and TSSs, etc.

One aspect not covered in this workflow is the utility of CAGE data (and 5’-end data in general) in providing accurate TSSs for studying other types of data. For example, having accurate TSSs is highly beneficial in chromatin research, since the location and nucleosome and TSSs are closely related [13, 39, 40]. CAGE can be combined with chromatin confirmation assays such as HiC to find new enhancers that are both co-expressed and physically interacting with TSSs. Many genome-wide association studies are finding that disease-related genetic variants are found in intergenic regions, that are often poorly annotated. The accurate enhancer locations provided by CAGE can greatly aid interpretation of such variants [41]. The adherence of CAGEfightR to standard Bioconductor classes facilitates these inter-assay analyses by making it easy to mix-and-match multiple packages developed for different experimental assays.

Software and data availability

The following software versions were used in this article:

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

Bioconductor version: 3.9

CAGEfightR version: 1.4.0

CAGEWorkflow: https://doi.org/10.18129/B9.bioc.CAGEWorkflow [42]

License: GPL-3

Mouse nanotube CAGE data: https://doi.org/10.18129/B9.bioc.nanotubes [24]

License: GPL-3

Thodberg and Sandelin present CAGEfightR, a new Bioconductor package suited for analysis of 5'end oriented datasets, derived from CAGE-seq and similar techniques.

Although similar Bioconductor packages exist (icetea, CAGEr, TSRchitect), the greatest strength of CAGEfightR is unique in its ability to call or predict putative enhancers based on the bidirectional transcription initiation signature, thus filling the gap in the current R-based toolkit for CAGE-like data analysis.

In addition, using CAGEfightR, hypotheses of enhancer-promoter interactions based on co-expression levels can be easily set and visualised.

The package is well-documented, and the step-by-step protocol well written and easy to follow. I only have a few minor comments which could benefit the general reader:

The introduction and abstract could state more clearly that CAGE allows TSS mapping of only RNA polymerase II transcripts. Though this is implied through usage of cap-trapping, I would keep in mind that this workflow might be used by general readers not so familiar with TSS-mapping techniques. Considering the advent of technologies that capture RNA polI-RNApolIII transcripts, expected to be much noisier, it would be better to make it as clear as possible. If the authors believe CAGEfightR could be of use on noisier data, from experimental methods based on negative selection (such as TSS-seq), it would be worth testing this and including a few sentences, as this would promote CAGEfightR usage on any TSS mapping technique.

I support comparison of existing packages in the form of Table 1; however, I would expand this comparison to include unique features that perhaps CAGEfightR does not have - e.g. icetea and TSRchitect support paired-end data, CAGEr has TSS-shifting discovery function and implemented G-correction function to remove mismatching G’s from 5’ends of reads.

Common problem with CAGE and CAGE-like data which is obtained through reverse transcription, is the addition of a G, or so called G-bias upstream of the true transcription start site.

It would be beneficial for general readers as this is a step-by-step protocol to discuss how to generate BigWig files from fastq files, and how to address the G-bias problem/i.e. remove the mismatching G’s at the 5’end of the reads.

The authors prefer to use the term TSS regions or TSSs in place of tag clusters as aggregated CTSSs, and even TSSs as a broader term for CTSS while CTSS is just a CAGE-supported TSS. This becomes very confusing for a general reader, especially page 36:

“Now we look at differential TSS usage at two-levels: Whether an individual TSS shows differential TSS usage (TSS-level) or whether a gene show differential TSS usage in any way (gene-level). First we can look at individual TSSs (TSS-level differential TSS usage).”

I would suggest usage of CTSS for individual CAGE-supported TSSs, tag cluster - for an aggregate of individual CTSSs based on distance based clustering or whatever methodology, and tag cluster can be interchangeably used with promoter where needed.

Interquartile range should probably be interquantile range as it spans middle 10-90th percentile of the signal. It would also be beneficial to explain why is it used instead of all signal (more robust measure that excludes outliers and is less sensitive to sequencing depth etc).

I am a bit surprised the authors use such harsh filtering step before plotting distributions of IQR (Figure 5, TPM >=10), I would assume that the problem is in tag clusters which seem sharp - single bp, and therefore it would perhaps be beneficial to add a more stringent filtering step only to single bp tag clusters to be retained only if highly expressed (>= 5 TPM), while a lower threshold can be applied on broad tag clusters as multiple CTSS within a tag cluster give more certainty that it is not just noise.

I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.

In this paper, the authors present a cookbook for analyzing CAGE data mainly through their own Bioconductor package, CAGEfightR, applied to sample data, which were analyzed and published by their group previously. By following the given source codes, readers will be able to learn how they can obtain various kinds of information rather easily. Thus, I feel that this publication will be useful for those who want to learn how to analyze CAGE data quickly. On the other hand, it will not tell us how to solve our deeper problems in research. I know that this is out of the scope of this tutorial but I can’t help but feeling, for example, how the given procedure for reducing batch effects can be justified (see below). Thus, here are some of my suggestions for its further improvement:

In our realistic situations, the EDA approach is quite important. In this sense, I appreciate their demonstration on how to remove batch effects from the expression data between various samples, using ComBat (Figs. 9 and 11). However, most researchers will not be satisfied with just seeing that the PC1 has become to separate positive and negative groups; it is natural that they would like to confirm if the correction was enough or not; they would also wish to see what the new PC2 as well as the old PC1 represent. Therefore, I recommend the authors to extend Table 2 for characterizing each sample from various features (e.g., experimental conditions and data size) and to use such features for the interpretation. More discussion and/or additional attempts to clarify the most probable main reason for the initial batch effects would be desirable.

Similarly, since the first author does not seem to have been a member of the previous analysis, it is interesting to see the consistency between the two studies. For example, the observation that inflammation-related genes were activated seems to be the same in both analyses. Then, are the genes with differential TSSs likely to explain the phenomenon? How much are the detected enhancers contributing to the differential expression? Do these enhancers (or newly activated TSSs) share any over-represented motifs? From the same reason, I recommend the authors to avoid using (ugly) chromosomal coordinates to represent genes/promoters/enhancers, wherever possible. It would be great if the authors can show that they could perform deeper analyses this time.

For the convenience of wider readers, it might be useful to show the way how to obtain BigWig files from rawer data (BAM or fastq, if possible).

Similarly, a summary table of used tools (except CAGEfightR), containing their input file information as well as their main purposes might be useful.

Also, it might be useful if there is a summary on what CAGEfightR can do/cannot do. For example, is it possible to combine different sources of CAGE data with this workflow?

One of my students tried to follow the workflow. At first, she failed to install some of the packages. It was because the version of R she used was R3.5. Thus, this point should be clearly noted. In addition, she reports that to run the code “Fit and shrink DE model”, the installation of “statmod” and “BiasedUrn” was necessary. Perhaps, some additional information on how to setup initial environment would be useful for beginners.

As a cookbook, it is desirable that users can find their necessary items more easily (via a table of contents with clearer headers, perhaps?).

There seems to be some confusion on the versions used:

As for the mouse genome sequence, (1) both mm9 and mm10 are used. (2) Why didn’t they use the latest version of JASPAR (Jaspar2018)?

The manuscript seems to contain many typos. Here are some that we found (there are likely to be more): (1)This is workflow is a case study on (2) CAGE dusters (3) can be speed up (4) novel TSSs are expresse d (5) to try infer (6 )this is a just a very quick (7) The returned GRanges contains the the location (8) going from clustering, to spatial analyses to differential expression.

We confirm that we have read this submission and believe that we have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.

Thodberg and Sandelin describe a comprehensive workflow for the basic analysis of CAGE data. The workflow is simple, well-presented and the code mostly runs without problems. I only have a few minor comments to improve its usefulness for the general reader:

There should be a general overview on how the CTSS BigWig files are generated from raw FASTQ files.

Consider using ExperimentHub for the nanotubes package. This provides greater control over which files are downloaded, rather than forcing the user to obtain all files at installation. This is especially useful if you have multiple data sets - see, for example, the chipseqDBData package. [1]

Explain what a BigWigFileList is, and why it needs to be used instead of having a simple character vector.

Consider not having string'ified intervals as the row names in the output of quickTSSs(). In some situations, generation of the strings can use more RAM than the actual data itself. It definitely slows down any attempt to 'show' the output. I would suggest that this be deferred as late as possible, e.g., until someone needs the strings as row names of a data frame to save to file.

More details are required on how quantification is performed. For example, I assume counts are summed directly from single strands for TSSs. For enhancers, are counts summed from both strands?

"As enhancers are the most complicated type, we keep only enhancers if a TSS and enhancer overlaps:" The complexity of the enhancers doesn't really provide a motivation for only keeping enhancers in cases of overlaps. The better reason is that all enhancers would be detected as two TSSs if the strandedness was ignored; if they do overlap, it would be more appropriate to interpret them as a single enhancer rather than as two distinct TSS events.

The Interquartile Range (IQR)... of what? I assume that the range refers to the length of the interval that contains 10 to 90% of a TSS's counts. Incidentally, the IQR is no longer an IQR if it's redefined from 10-90%.

The single quotes in the highTSSs call are malformed, which prevents copy-pasting.

One could consider using a 2-component mixture model to classify elements into sharp/broad in a more automated manner.

I presume that the pairwise correlations for the TSS-enhancer interactions are computed across samples for each TSS/enhancer pair. If so, is this done after blocking on the design? Otherwise it is possible to obtain strong positive correlations simply because a TSS and the enhancer happen to respond in the same direction to a given treatment. If there is a genuine physical interaction, it should manifest as correlations within each treatment condition, where stochastic differences in RNA polymerase activity affect both the TSS and its interacting enhancer.

There is no correction for multiple testing in the p-values from the links. While I recognise that this is a limitation of the small sample size, it should still be pointed out as a caveat of the analysis.

All normalization steps in the DE analyses assume that most of the input features are not differentially expressed between conditions. This is usually not a concern, but if aggressive feature selection is performed, it may become a problem. For example, if one were to perform the DE analysis on superenhancers only, it would give incorrect results in situations where the superenhancer activity increases globally upon treatment.

I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.

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

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

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.345412	0
chr19:56668148-56668332;+	TSS	CDS	103.8795	-2.254371	0

Table 6.

Global summary of differentially expressed genes.

	(Intercept)	BatchB	ClassNano
Down	51	2572	1505
NotSig	463	8278	10373
Up	13053	2717	1689

Table 7.

Top differentially expressed genes.

	symbol	nClusters	AveExpr	logFC	adj.P.Val
66938	Sh3d21	3	5.871004	3.075745	0.0e+00
245049	Myrip	2	4.371325	2.414055	7.0e-07
12722	Clca3a1	1	3.020528	3.692198	7.0e-07
382864	Colq	3	2.770158	-3.426911	1.1e-06
20716	Serpina3n	5	6.384175	1.872782	3.0e-06
72275	2200002D01Rik	2	7.208031	1.693257	5.5e-06
381813	Prmt8	4	4.553612	1.409006	5.8e-06
170706	Tmem37	2	5.503908	1.679690	5.8e-06
18654	Pgf	1	4.862055	2.337045	5.8e-06
20361	Sema7a	1	7.612236	1.473680	5.9e-06

Table 8.

Top enriched or depleted GO-terms.

	Term	Ont	N	Up	Down	P.Up	P.Down
GO:0006954	inflammatory response	BP	556	142	51	0	0.9562685
GO:0006952	defense response	BP	1072	224	99	0	0.9878373
GO:0097529	myeloid leukocyte migration	BP	170	61	14	0	0.9359984
GO:0010033	response to organic substance	BP	2074	370	196	0	0.9987104
GO:0006950	response to stress	BP	2755	464	246	0	0.9999946
GO:0006955	immune response	BP	1034	210	96	0	0.9833226
GO:0042221	response to chemical	BP	2762	467	292	0	0.9178712
GO:0050900	leukocyte migration	BP	288	83	23	0	0.9792828
GO:0001816	cytokine production	BP	634	143	45	0	0.9998658
GO:0001817	regulation of cytokine production	BP	570	132	39	0	0.9998856

Table 9.

Top enriched of depleted KEGG-terms.

	Pathway	N	Up	Down	P.Up	P.Down
path:mmu04060	Cytokine-cytokine receptor interaction	173	56	13	0.0000000	0.9579351
path:mmu04668	TNF signaling pathway	105	31	8	0.0000037	0.9186628
path:mmu00600	Sphingolipid metabolism	41	17	2	0.0000051	0.9583011
path:mmu00980	Metabolism of xenobiotics by cytochrome P450	48	4	17	0.8857194	0.0000137
path:mmu03010	Ribosome	122	32	2	0.0000226	0.9999900
path:mmu04064	NF-kappa B signaling pathway	85	24	5	0.0000704	0.9655534
path:mmu04657	IL-17 signaling pathway	74	22	2	0.0000806	0.9985563
path:mmu00982	Drug metabolism - cytochrome P450	46	5	15	0.7266916	0.0001238
path:mmu04630	JAK-STAT signaling pathway	112	29	7	0.0001453	0.9785951
path:mmu04512	ECM-receptor interaction	69	21	13	0.0001488	0.0577601

Table 10.

Top differentially used TSSs.

	txType	geneID	symbol	logFC	FDR
chr17:13840650-13840851;-	intron	21646	Tcte2	1.7889344	0e+00
chr10:57857044-57857314;+	promoter	110829	Lims1	-1.0651946	0e+00
chr14:70215678-70215876;-	intron	246710	Rhobtb2	2.4933979	0e+00
chr4:141154044-141154185;-	intron	74202	Fblim1	1.7018062	0e+00
chr17:33966135-33966308;+	intron	66416	Ndufa7	2.1612127	0e+00
chr15:76428030-76428201;-	intron	94230	Cpsf1	1.4598815	0e+00
chr19:57271818-57272125;-	promoter	226251	Ablim1	1.1456163	0e+00
chr9:77788968-77789200;+	intron	68801	Elovl5	0.9810692	1e-07
chr11:116395161-116395462;+	proximal	20698	Sphk1	1.7471930	1e-07
chr2:91496305-91496449;+	intron	228359	Arhgap1	0.9809491	3e-07

Table 11.

Top genes showing any differential TSS usage.

geneID	symbol	NExons	FDR
21646	Tcte2	4	0e+00
110829	Lims1	3	0e+00
246710	Rhobtb2	3	0e+00
74202	Fblim1	3	0e+00
66416	Ndufa7	3	0e+00
94230	Cpsf1	2	0e+00
226251	Ablim1	3	0e+00
68801	Elovl5	2	1e-07
20698	Sphk1	3	1e-07
228359	Arhgap1	2	2e-07