6.1.1 Calculating pooled CTSSs

For clustering CTSSs, we first want to calculate the the overall CTSSs signal across all samples, called pooled CTSSs. However, since the different samples have different library sizes, CTSS counts must first be normalized. The simplest way of doing this is to use Tags-Per-Million (TPM, not to be confused with Transcripts-Per-Million or TxPM used in RNA-Seq):

exampleCTSSs <- calcTPM(exampleCTSSs, inputAssay = "counts", outputAssay = "TPM", 
    outputColumn = "subsetTags")
#> Calculating library sizes...
#> Calculating TPM...

calcTPM will calculate the total number of tags for each sample and store them in colData, and then scale counts into TPM and add them as a new assay:

# Library sizes
colData(exampleCTSSs)
#> DataFrame with 3 rows and 5 columns
#>             Name       BigWigPlus       BigWigMinus totalTags subsetTags
#>      <character>      <character>       <character> <numeric>  <numeric>
#> C547        C547 mm9.C547.plus.bw mm9.C547.minus.bw  12436172      63230
#> C548        C548 mm9.C548.plus.bw mm9.C548.minus.bw  12277807      52342
#> C549        C549 mm9.C549.plus.bw mm9.C549.minus.bw  14477276      69728

# TPM values
head(assay(exampleCTSSs, "TPM"))
#> 6 x 3 sparse Matrix of class "dgCMatrix"
#>          C547     C548     C549
#> [1,] 15.81528  .        .      
#> [2,]  .        .       14.34144
#> [3,] 15.81528  .        .      
#> [4,]  .       19.10512  .      
#> [5,] 15.81528  .        .      
#> [6,]  .        .       14.34144

As this is just a subset of the original dataset, we instead use the total number tags from the the complete dataset, by specifying the name of the column in colData (Note that a warning is passed that we are overwritting the previous TPM assay):

exampleCTSSs <- calcTPM(exampleCTSSs, inputAssay = "counts", totalTags = "totalTags", 
    outputAssay = "TPM")
#> Using supplied library sizes...
#> Warning in calcTPM(exampleCTSSs, inputAssay = "counts", totalTags =
#> "totalTags", : object already has an assay named TPM: It will be overwritten!
#> Calculating TPM...
head(assay(exampleCTSSs, "TPM"))
#> 6 x 3 sparse Matrix of class "dgCMatrix"
#>           C547       C548       C549
#> [1,] 0.0804106 .          .         
#> [2,] .         .          0.06907377
#> [3,] 0.0804106 .          .         
#> [4,] .         0.08144777 .         
#> [5,] 0.0804106 .          .         
#> [6,] .         .          0.06907377

calcPooled will calculate pooled CTSSs by summing up TPM across sample for each CTSS and store it in the score-column:

# Library sizes
exampleCTSSs <- calcPooled(exampleCTSSs, inputAssay = "TPM")
#> Warning in calcPooled(exampleCTSSs, inputAssay = "TPM"): object already has a
#> column named score in rowData: It will be overwritten!
rowRanges(exampleCTSSs)
#> GRanges object with 41256 ranges and 1 metadata column:
#>           seqnames    ranges strand |              score
#>              <Rle> <IRanges>  <Rle> |          <numeric>
#>       [1]    chr18  72690884      + | 0.0804105958007014
#>       [2]    chr18  72730496      + | 0.0690737677447056
#>       [3]    chr18  73169663      + | 0.0804105958007014
#>       [4]    chr18  73190583      + | 0.0814477699478417
#>       [5]    chr18  73234436      + | 0.0804105958007014
#>       ...      ...       ...    ... .                ...
#>   [41252]    chr19  61304387      - | 0.0814477699478417
#>   [41253]    chr19  61304427      - | 0.0690737677447056
#>   [41254]    chr19  61332796      - | 0.0690737677447056
#>   [41255]    chr19  61333667      - |  0.162895539895683
#>   [41256]    chr19  61333668      - | 0.0814477699478417
#>   -------
#>   seqinfo: 35 sequences (1 circular) from mm9 genome

6.1.2 Calculating CTSS support to remove excess noise

In some cases you might have a very large number of and/or very noisy samples (For example due to poor RNA quality), which can lead to an increase in the number of single tags spread across the genome. To alleviate this issue, CTSSs appearing in only a single or few samples can be discarded. We refer to this as calculating the support: the number of samples expressing a feature above some threshold:

# Count number of samples with MORE ( > ) than 0 counts:
exampleCTSSs <- calcSupport(exampleCTSSs, inputAssay = "counts", outputColumn = "support", 
    unexpressed = 0)
table(rowRanges(exampleCTSSs)$support)
#> 
#>     1     2     3 
#> 33162  4297  3797

The majority of CTSSs are only expressed in a single sample. We can discard these sites using subset, and then recalculate TPM values:

supportedCTSSs <- subset(exampleCTSSs, support > 1)
supportedCTSSs <- calcTPM(supportedCTSSs, totalTags = "totalTags")
#> Using supplied library sizes...
#> Warning in calcTPM(supportedCTSSs, totalTags = "totalTags"): object already
#> has an assay named TPM: It will be overwritten!
#> Calculating TPM...
supportedCTSSs <- calcPooled(supportedCTSSs)
#> Warning in calcPooled(supportedCTSSs): object already has a column named
#> score in rowData: It will be overwritten!

Note, CAGEfightR warned that it was overwritting columns: Another option would be to save the output as new columns. The subsetBySupport function wraps the common task of calling calcSupport and subset.

6.1.3 Unidirectional (tag) clustering: Finding Transcription Start Sites (TSSs)

Once we have obtained pooled CTSSs, we can identify various types of clusters: Transcripts are rarely transcribed from just a single CTSS, but rather from a collection or cluster of nearby CTSSs corresponding to a single TSS.

The simplest way of doing this is to group nearby CTSSs on the same strand into Tag Clusters (TCs), which are likely candidates for real TSSs (although some post-filtering is almost always a good idea, see next section). CAGEfightR does this by a slice-reduce approach: It finds sites above some value (slice) and then merges nearby sites (reduce). In the simplest form this can be run as:

naive_TCs <- clusterUnidirectionally(exampleCTSSs, pooledCutoff = 0, mergeDist = 20)
#> Splitting by strand...
#> Finding tag clusters...
#> Calculating tag cluster statistics...
#> Preparing output...
#> Tag clustering summary:
#>   Width Count Percent
#>   Total 21008 100.0 %
#>     >=1 18387  87.5 %
#>    >=10  2487  11.8 %
#>   >=100   134   0.6 %
#>  >=1000     0   0.0 %

clusterUnidirectionally find TCs and outputs some summary statistics. We see that we have some very wide TCs (> 1000 bp wide). This is often the result of a just few spread-out tags that link several major TCs, which does not represent the TSS structure very well. This can be avoided by tuning the cutoff used for slicing: The smallest value for cutoff value producing the largest number of TCs is chosen:

tuned <- tuneTagClustering(exampleCTSSs, mergeDist = 20)
#> Finding thresholds to be tested...
#> Splitting coverage by strand...
#> Iterating over 11 thresholds using 4 worker(s)...
#> Analyzing plus strand...
#> Analyzing minus strand...
#> Preparing output...
tuned
#>     threshold  nTCs
#> 1  0.00000000 21008
#> 2  0.06907377 15628
#> 3  0.08041060  8861
#> 4  0.08144777  3392
#> 5  0.13814754  3044
#> 6  0.14948436  2694
#> 7  0.15052154  2421
#> 8  0.16082119  2006
#> 9  0.16185837  1743
#> 10 0.16289554  1417
#> 11 0.20722130  1389

Increasing the slice-cutoff by a very small amount greatly increases the number of TCs, by preventing merging of major TCs by weakly expressed CTSSs. We can then use this tuned cutoff value to identify TCs:

optimalCutoff <- tuned[which.max(tuned$nTCs), 1]
TCs <- clusterUnidirectionally(exampleCTSSs, pooledCutoff = optimalCutoff, mergeDist = 20)
#> Splitting by strand...
#> Finding tag clusters...
#> Calculating tag cluster statistics...
#> Preparing output...
#> Tag clustering summary:
#>   Width Count Percent
#>   Total 21008 100.0 %
#>     >=1 18387  87.5 %
#>    >=10  2487  11.8 %
#>   >=100   134   0.6 %
#>  >=1000     0   0.0 %

Note, we no longer have any very long TCs (For a real dataset it’s not uncommon to have a few very long TCs left). Let’s take a closer look at what’s returned:

TCs
#> GRanges object with 21008 ranges and 2 metadata columns:
#>                             seqnames            ranges strand |
#>                                <Rle>         <IRanges>  <Rle> |
#>   chr18:72690884-72690884;+    chr18          72690884      + |
#>   chr18:72730496-72730496;+    chr18          72730496      + |
#>   chr18:73169663-73169663;+    chr18          73169663      + |
#>   chr18:73190583-73190583;+    chr18          73190583      + |
#>   chr18:73234436-73234436;+    chr18          73234436      + |
#>                         ...      ...               ...    ... .
#>   chr19:61304092-61304092;-    chr19          61304092      - |
#>   chr19:61304278-61304387;-    chr19 61304278-61304387      - |
#>   chr19:61304427-61304427;-    chr19          61304427      - |
#>   chr19:61332796-61332796;-    chr19          61332796      - |
#>   chr19:61333667-61333668;-    chr19 61333667-61333668      - |
#>                                          score     thick
#>                                      <numeric> <IRanges>
#>   chr18:72690884-72690884;+ 0.0804105958007014  72690884
#>   chr18:72730496-72730496;+ 0.0690737677447056  72730496
#>   chr18:73169663-73169663;+ 0.0804105958007014  73169663
#>   chr18:73190583-73190583;+ 0.0814477699478417  73190583
#>   chr18:73234436-73234436;+ 0.0804105958007014  73234436
#>                         ...                ...       ...
#>   chr19:61304092-61304092;- 0.0690737677447056  61304092
#>   chr19:61304278-61304387;-   86.4736312947607  61304320
#>   chr19:61304427-61304427;- 0.0690737677447056  61304427
#>   chr19:61332796-61332796;- 0.0690737677447056  61332796
#>   chr19:61333667-61333668;-  0.244343309843525  61333667
#>   -------
#>   seqinfo: 35 sequences (1 circular) from mm9 genome

The GRanges includes the chromosome, start and end positions of each TCs, as well as two other key values: The TC peak in the thick column: This the single position within the TC with the highest pooled TPM, and the sum of pooled TPM within the TC in the score column. In case you want to save TCs to a bed file, you can use the export function from the rtracklayer package.

6.1.4 Bidirectional clustering: Finding Enhancers

CAGE is a unique technology in that it allows for the robust identification of enhancers by transcription of enhancer RNAs (eRNAs). Active enhancers produce weak but consistent bidirectional transcription of capped eRNAs, resulting in a characteristic CTSS pattern of two diverging peaks approximally 400 basepairs apart.

CAGEfightR can identify this pattern by calculating a balance score ranging from 0 to 1, where 1 corresponds to perfectly divergent site (Technically, the balance score is the Bhattacharyya coefficient measuring agreement with the “ideal” divergent site, see the CAGEfightR paper). Again, a slice-reduce approach can be use to identify invidividual enhancers based on the balance scores. We refer to this approach as bidirectional clustering as opposed to the conventional unidirectonal cluster used to identify TSSs: and hence the name of the cluster functions: clusterBidirectionally and clusterUnidirectionally:

enhancers <- clusterBidirectionally(exampleCTSSs, balanceThreshold = 0.95)
#> 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...

Let’s look closer at the returned GRanges:

enhancers
#> GRanges object with 674 ranges and 3 metadata columns:
#>                           seqnames            ranges strand |
#>                              <Rle>         <IRanges>  <Rle> |
#>   chr18:73731202-73731784    chr18 73731202-73731784      * |
#>   chr18:73731956-73732359    chr18 73731956-73732359      * |
#>   chr18:73750120-73750564    chr18 73750120-73750564      * |
#>   chr18:73752492-73752953    chr18 73752492-73752953      * |
#>   chr18:73802190-73802733    chr18 73802190-73802733      * |
#>                       ...      ...               ...    ... .
#>   chr19:61182777-61183299    chr19 61182777-61183299      * |
#>   chr19:61193381-61193912    chr19 61193381-61193912      * |
#>   chr19:61194241-61194692    chr19 61194241-61194692      * |
#>   chr19:61202166-61202703    chr19 61202166-61202703      * |
#>   chr19:61260539-61260971    chr19 61260539-61260971      * |
#>                                       score     thick           balance
#>                                   <numeric> <IRanges>         <numeric>
#>   chr18:73731202-73731784 0.160821191601403  73731493                 1
#>   chr18:73731956-73732359  1.22166008691667  73732159 0.997127396724092
#>   chr18:73750120-73750564 0.242268961549244  73750342 0.985341340592472
#>   chr18:73752492-73752953 0.149484363545407  73752722 0.999279749024173
#>   chr18:73802190-73802733 0.150521537692547  73802461 0.999153450069061
#>                       ...               ...       ...               ...
#>   chr19:61182777-61183299  0.31134272929395  61183042 0.999802454365801
#>   chr19:61193381-61193912 0.762907342371592  61193599 0.999770184445126
#>   chr19:61194241-61194692 0.218558131290113  61194460 0.999279749024173
#>   chr19:61202166-61202703 0.149484363545407  61202434 0.999279749024173
#>   chr19:61260539-61260971 0.138147535489411  61260755                 1
#>   -------
#>   seqinfo: 35 sequences (1 circular) from mm9 genome

The GRanges includes the chromosome, start and end positions of each enhancer (enhancers are unstranded), as well as two other key values: the score column holds the sum of pooled CTSSs of the enhancer (on both strands) and the enhancer midpoint in the thick column: This is the maximally balanced site in the enhancers. Again, the export function from the rtracklayer package can be used to export enhancers to a BED-file.

As balance here is solely defined on the pooled CTSSs, a useful filter is to make sure the bidirectional site is also bidirectional in at least a single sample. Similarly to how we calculated support, we can calculate the sample-wise bidirectionality of enhancers:

# Calculate number of bidirectional samples
enhancers <- calcBidirectionality(enhancers, samples = exampleCTSSs)

# Summarize
table(enhancers$bidirectionality)
#> 
#>   0   1   2   3 
#> 284 319  44  27

Many enhancers are not observed to be bidirectional in one or more samples. We can remove these using subset (The subsetBySupport function wraps the common task of calling calcBidirectionality and subset):

enhancers <- subset(enhancers, bidirectionality > 0)

This bidirectional clustering approach identifies any bidirectional site in the genome. This means that for example bidirectional promoters will also detected (given that they are sufficiently balanced). To remove these cases, annotation with transcript models can be used to remove bidirectional clusters overlapping known promoters and exons. This is described in the next section.

6.2.1 Quantifying expression at cluster level

To look at differential expression, we must quantify the expression of each TC in each sample (in almost all cases you want to quantify the raw CTSS counts). This can be done using the quantifyClusters function (The example datasets have already been quantified, so here we will re-quantify the clusters):

requantified_TSSs <- quantifyClusters(exampleCTSSs, clusters = rowRanges(exampleUnidirectional), 
    inputAssay = "counts")
#> Finding overlaps...
#> Aggregating within clusters...
requantified_enhancers <- quantifyClusters(exampleCTSSs, clusters = rowRanges(exampleBidirectional), 
    inputAssay = "counts")
#> Finding overlaps...
#> Aggregating within clusters...

This returns a RangedSummarizedExperiment with clusters in rowRanges and count matrix in assays.

6.2.2 Removing weakly expressed clusters

We often wish to remove weakly expressed clusters. Similarly to CTSSs, a simple approach is to remove cluster based on counts using the subsetBySupport function:

# Only keep enhancers expressed in more than one sample
supported_enhancers <- subsetBySupport(exampleBidirectional, inputAssay = "counts", 
    unexpressed = 0, minSamples = 1)
#> Calculating support...
#> Subsetting...
#> Removed 105 out of 390 regions (26.9%)

Similarly, we can also first calculate TPM values for clusters, and then filter on TPM values:

# Calculate TPM using pre-calculated total tags:
exampleUnidirectional <- calcTPM(exampleUnidirectional, totalTags = "totalTags")
#> Using supplied library sizes...
#> Calculating TPM...

# Only TSSs expressed at more than 1 TPM in more than 2 samples
exampleUnidirectional <- subsetBySupport(exampleUnidirectional, inputAssay = "TPM", 
    unexpressed = 1, minSamples = 2)
#> Calculating support...
#> Subsetting...
#> Removed 17141 out of 17313 regions (99%)

6.2.3 Annotation with transcript models

While CAGE can identify TSSs and enhancers completely independent of annotation, it is often useful to compare these to existing annotations. Bioconductor stores transcript models using the TxDb class from the GenomicFeatures package. There are multiple ways to obtain a TxDb object:

• Use the TxDb-objects included as Bioconductor packages: For this vignette we use the TxDb.Mmusculus.UCSC.mm9.knownGene package.
• Import a GTF or GFF3 file using the makeTxDbFromGFF function from GenomicFeatures.
• Use the AnnotationHub package to directly download TxDb objects for a wide range of organisms.
• Use makeTxDbFromBiomart, makeTxDbFromUCSC or makeTxDbFromEnsembl from from GenomicFeatures to create TxDb objects from various online sources.

In case you really do not want to use a TxDb object, all annotation functions in CAGEfightR also accepts GRanges or GRangesList as input.

library(TxDb.Mmusculus.UCSC.mm9.knownGene)
txdb <- TxDb.Mmusculus.UCSC.mm9.knownGene
txdb
#> TxDb object:
#> # Db type: TxDb
#> # Supporting package: GenomicFeatures
#> # Data source: UCSC
#> # Genome: mm9
#> # Organism: Mus musculus
#> # Taxonomy ID: 10090
#> # UCSC Table: knownGene
#> # Resource URL: http://genome.ucsc.edu/
#> # Type of Gene ID: Entrez Gene ID
#> # Full dataset: yes
#> # miRBase build ID: NA
#> # transcript_nrow: 55419
#> # exon_nrow: 246570
#> # cds_nrow: 213117
#> # Db created by: GenomicFeatures package from Bioconductor
#> # Creation time: 2015-10-07 18:13:02 +0000 (Wed, 07 Oct 2015)
#> # GenomicFeatures version at creation time: 1.21.30
#> # RSQLite version at creation time: 1.0.0
#> # DBSCHEMAVERSION: 1.1

First, we might simply inquire as to what transcripts overlaps our clusters, which is done via the assignTxID function:

exampleUnidirectional <- assignTxID(exampleUnidirectional, txModels = txdb, outputColumn = "txID")
#> Extracting transcripts...
#> Finding hierachical overlaps...
#> ### Overlap Summary: ###
#> Features overlapping transcripts: 87.79 %
#> Number of unique transcripts: 232

Note, that a single TC can overlap multiple transcripts (names are separated with a ;):

rowRanges(exampleUnidirectional)[5:6]
#> GRanges object with 2 ranges and 4 metadata columns:
#>                             seqnames            ranges strand |
#>                                <Rle>         <IRanges>  <Rle> |
#>   chr18:74427830-74427996;+    chr18 74427830-74427996      + |
#>   chr18:74442758-74442837;+    chr18 74442758-74442837      + |
#>                                        score     thick   support
#>                                    <numeric> <IRanges> <integer>
#>   chr18:74427830-74427996;+ 63.4299739643187  74427933         3
#>   chr18:74442758-74442837;+ 11.2621271351077  74442817         3
#>                                              txID
#>                                       <character>
#>   chr18:74427830-74427996;+ uc008fph.2;uc008fpj.2
#>   chr18:74442758-74442837;+ uc008fpl.2;uc008fpm.2
#>   -------
#>   seqinfo: 35 sequences (1 circular) from mm9 genome

In many cases, we are not so much interested in what specific transcripts a cluster overlaps, but rather where in a transcript a cluster lies - for example whether it is located at the annotated TSSs, in the 5’-UTR or in the CDS region. However, due to overlapping transcripts, these categories might overlap, for example a region might be an exon in one transcripts, but skipped in another due to alternative splicing. To get around this, CAGEfightR implements a hierachical annotation scheme (See the CAGEfightR paper for a detailed description), implemented in the assignTxType function:

exampleUnidirectional <- assignTxType(exampleUnidirectional, txModels = txdb, 
    outputColumn = "txTyp")
#> Finding hierachical overlaps...
#> ### Overlap summary: ###
#>       txType count percentage
#> 1   promoter   101       58.7
#> 2   proximal    11        6.4
#> 3    fiveUTR    15        8.7
#> 4   threeUTR     3        1.7
#> 5        CDS     8        4.7
#> 6       exon     0        0.0
#> 7     intron    13        7.6
#> 8  antisense     6        3.5
#> 9 intergenic    15        8.7

The function returns a small summary of the annotation: At the top of the hierachy we have overlap with annotated promoter (+/- 100 basepairs from the annotated TSSs, with this distance being changeable when running the function), followed by proximal (- 1000 upstrema of the annotated TSS). Coding transcripts can be annotated as 5’-UTR, 3-’UTR, CDS and intron, while non-conding transcripts are simply annotated as exon or intron. At the bottom of the hierachy, clusters overlapping genes on the opposite strand are annotated as antisense. Finally, in case a cluster does not overlap any transcript, it is annotated as intergenic.

Highly expressed TSSs can be quite wide, and overlap many different categories. It can be useful to only annotate the TSSs peak rather than the the entire TSSs to obtain more representative annotations:

exampleUnidirectional <- assignTxType(exampleUnidirectional, txModels = txdb, 
    outputColumn = "peakTxType", swap = "thick")
#> Finding hierachical overlaps with swapped ranges...
#> ### Overlap summary: ###
#>       txType count percentage
#> 1   promoter    94       54.7
#> 2   proximal    15        8.7
#> 3    fiveUTR    17        9.9
#> 4   threeUTR     3        1.7
#> 5        CDS     6        3.5
#> 6       exon     0        0.0
#> 7     intron    16        9.3
#> 8  antisense     6        3.5
#> 9 intergenic    15        8.7

Annotation with transcript models are particular useful for enhancer detecting, since bidirectional cluster might in some cases be very balanced bidirectional promoters. We can remove enhancers overlapping known transcripts:

# Annotate with TxTypes
exampleBidirectional <- assignTxType(exampleBidirectional, txModels = txdb, outputColumn = "txType")
#> Finding hierachical overlaps...
#> ### Overlap summary: ###
#>       txType count percentage
#> 1   promoter    23        5.9
#> 2   proximal    15        3.8
#> 3    fiveUTR     4        1.0
#> 4   threeUTR    29        7.4
#> 5        CDS    45       11.5
#> 6       exon     5        1.3
#> 7     intron   231       59.2
#> 8  antisense     0        0.0
#> 9 intergenic    38        9.7

# Only keep intronic and intergenic enhancers
exampleBidirectional <- subset(exampleBidirectional, txType %in% c("intron", "intergenic"))

Of course, this way of filtering relies on having accurate transcript models!

6.2.4 Quantifying shape of Tag Clusters

We have already looked at two characteristics of Tag Clusters (TCs) in addition their genomic location: the expression level (score column) and TSSs peak (thick column). One can further describe TCs beyond these two measures. For example, it was found that mammalian TSSs can be divided into “sharp” or “broad” TSSs, by comparing the Interquartile Range (IQR). We refer to such a measure as a shape statistic and CAGEfightR includes a few built-in functions for shape statistics (IQR, entropy and multimodality) as well as an easy framework for implementing your own.

Let’s see how this works by calculating the IQR for TSSs. First, we need to calculate pooled CTSS:

# Recalculate pooled signal
exampleCTSSs <- calcTPM(exampleCTSSs, totalTags = "totalTags")
#> Using supplied library sizes...
#> Calculating TPM...
exampleCTSSs <- calcPooled(exampleCTSSs)
#> Warning in calcPooled(exampleCTSSs): object already has a column named score
#> in rowData: It will be overwritten!

# Calculate shape
exampleUnidirectional <- calcShape(exampleUnidirectional, pooled = exampleCTSSs, 
    outputColumn = "IQR", shapeFunction = shapeIQR, lower = 0.25, upper = 0.75)
#> Splitting coverage by strand...
#> Applying function to each TC...
#> Assembling output...

A plot of the distribution of IQR values reveals the distinction between broad and sharp TSSs:

hist(rowRanges(exampleUnidirectional)$IQR, breaks = max(rowRanges(exampleUnidirectional)$IQR), 
    xlim = c(0, 100), xlab = "IQR", col = "red")

shapeEntropy works in the same way. Advanced users may implement their own shape statistic by writing a new function:

# Write a function that quantifies the lean of a TSS
shapeLean <- function(x, direction) {
    # Coerce to normal vector
    x <- as.vector(x)
    
    # Find highest position:
    i <- which.max(x)
    
    # Calculate sum
    i_total <- sum(x)
    
    # Calculate lean fraction
    if (direction == "right") {
        i_lean <- sum(x[i:length(x)])
    } else if (direction == "left") {
        i_lean <- sum(x[1:i])
    } else {
        stop("direction must be left or right!")
    }
    
    # Calculate lean
    o <- i_lean/i_total
    
    # Return
    o
}

# Calculate lean statistics, additional arguments can be passed to calcShape
# via '...'
exampleUnidirectional <- calcShape(exampleUnidirectional, exampleCTSSs, outputColumn = "leanRight", 
    shapeFunction = shapeLean, direction = "right")
#> Splitting coverage by strand...
#> Applying function to each TC...
#> Assembling output...

exampleUnidirectional <- calcShape(exampleUnidirectional, exampleCTSSs, outputColumn = "leanLeft", 
    shapeFunction = shapeLean, direction = "left")
#> Splitting coverage by strand...
#> Applying function to each TC...
#> Assembling output...

6.3.1 Annotation with gene models

In a similar fashion to transcript annotation described above Annotation with transcript models, gene models are obtained via TxDb objects:

# Load example TSS
data(exampleUnidirectional)

# Keep only TCs expressed at more than 1 TPM in more than 2 samples:
exampleUnidirectional <- calcTPM(exampleUnidirectional, totalTags = "totalTags")
#> Using supplied library sizes...
#> Calculating TPM...
exampleUnidirectional <- subsetBySupport(exampleUnidirectional, inputAssay = "TPM", 
    unexpressed = 1, minSamples = 2)
#> Calculating support...
#> Subsetting...
#> Removed 17141 out of 17313 regions (99%)

# Use the Bioconductor mm9 UCSC TxXb
library(TxDb.Mmusculus.UCSC.mm9.knownGene)
txdb <- TxDb.Mmusculus.UCSC.mm9.knownGene

The assignGeneID function annotate ranges with gene IDs:

exampleUnidirectional <- assignGeneID(exampleUnidirectional, geneModels = txdb, 
    outputColumn = "geneID")
#> Extracting genes...
#> Overlapping while taking distance to nearest TSS into account...
#> Finding hierachical overlaps...
#> ### Overlap Summary: ###
#> Features overlapping genes: 88.37 %
#> Number of unique genes: 89

TSSs not overlappping any genes receive NA, and in the (usually rare) case of a TC overlapping multiple genes, the gene with the closest annotated TSS is chosen:

What type of gene ID is returned of course depends on the source of the TxDb object supplied. In this case, it is Entrez IDs, the most commonly used ID system in Bioconductor. Entrez IDs can be used to find the corresponding gene symbols using an OrgDb object:

# Use Bioconductor OrgDb package
library(org.Mm.eg.db)
#> 
odb <- org.Mm.eg.db

# Match IDs to symbols
symbols <- mapIds(odb, keys = rowRanges(exampleUnidirectional)$geneID, keytype = "ENTREZID", 
    column = "SYMBOL")
#> 'select()' returned 1:1 mapping between keys and columns

# Add to object
rowRanges(exampleUnidirectional)$symbol <- as.character(symbols)

For most gene-level operations in CAGEfightR NAs (TCs not overlapping genes) are simply removed. If you wish to keep these TCs, you can relabel them as “Novel” genes via the assignMissingID function:

exampleUnidirectional <- assignMissingID(exampleUnidirectional, outputColumn = "symbol")
#> Assigned 2 missing IDs

6.3.2 Quantify expression at gene-level

Once TCs have been assigned to genes, TC expression can be summed up within genes to obtain a gene-level expression matrix via the quantifyGenes function:

genelevel <- quantifyGenes(exampleUnidirectional, genes = "geneID", inputAssay = "counts")

This returns a RangedSummarizedExperiment with the same number of rows as genes. The ranges are returned as a GRangesList of the individual TSSs making up the expression of the genes:

rowRanges(genelevel)
#> GRangesList object of length 88:
#> $100217439 
#> GRanges object with 5 ranges and 5 metadata columns:
#>                             seqnames            ranges strand |
#>                                <Rle>         <IRanges>  <Rle> |
#>   chr19:53899602-53899622;-    chr19 53899602-53899622      - |
#>   chr19:55817615-55817733;-    chr19 55817615-55817733      - |
#>   chr19:56738135-56738213;-    chr19 56738135-56738213      - |
#>   chr19:57432920-57432984;-    chr19 57432920-57432984      - |
#>   chr19:60842388-60842675;-    chr19 60842388-60842675      - |
#>                                        score     thick   support      geneID
#>                                    <numeric> <IRanges> <integer> <character>
#>   chr19:53899602-53899622;- 5.20115106166876  53899605         3   100217439
#>   chr19:55817615-55817733;- 7.01452555772393  55817695         3   100217439
#>   chr19:56738135-56738213;- 11.9559511230901  56738161         3   100217439
#>   chr19:57432920-57432984;- 7.02786464660017  57432963         3   100217439
#>   chr19:60842388-60842675;- 14.6610146789778  60842645         3   100217439
#>                                  symbol
#>                             <character>
#>   chr19:53899602-53899622;-     Snora19
#>   chr19:55817615-55817733;-     Snora19
#>   chr19:56738135-56738213;-     Snora19
#>   chr19:57432920-57432984;-     Snora19
#>   chr19:60842388-60842675;-     Snora19
#> 
#> $100217466 
#> GRanges object with 1 range and 5 metadata columns:
#>                             seqnames            ranges strand |
#>   chr18:74938323-74938515;-    chr18 74938323-74938515      - |
#>                                        score    thick support    geneID
#>   chr18:74938323-74938515;- 243.886676380891 74938434       3 100217466
#>                               symbol
#>   chr18:74938323-74938515;- Scarna17
#> 
#> $100502841 
#> GRanges object with 1 range and 5 metadata columns:
#>                             seqnames            ranges strand |
#>   chr18:78135139-78135265;+    chr18 78135139-78135265      + |
#>                                        score    thick support    geneID
#>   chr18:78135139-78135265;+ 101.861136805853 78135206       3 100502841
#>                             symbol
#>   chr18:78135139-78135265;+   Epg5
#> 
#> ...
#> <85 more elements>
#> -------
#> seqinfo: 35 sequences (1 circular) from mm9 genome

Gene-level values are accesible via the rowData function:

rowData(genelevel)
#> DataFrame with 88 rows and 1 column
#>     nClusters
#>     <integer>
#> 1           5
#> 2           1
#> 3           1
#> 4           1
#> 5           2
#> ...       ...
#> 84          1
#> 85          2
#> 86          1
#> 87          1
#> 88          1

The gene-level expression matrix can now be supplied to many other Bioconductor packages, including DESeq2, edgeR, limma, etc.

6.3.3 Filtering clusters based on gene composition

When looking at Differential TSS Usage (DTU) we are interested in changes in the contribution of individual TCs to overall gene expression. Since CAGE can detect even very lowly expressed TSSs, it is common to see TSSs making up only a very small fraction of the total gene expression. Removing these TSSs can improve power and clarity of downstream analyses.

CAGEfightR can calculate the composition value of intragenic TCs, defined as the number of samples expressing a TC at more than a certain fraction of total gene expression. This is implemented in the calcComposition function:

# Remove TSSs not belonging to any genes
intragenicTSSs <- subset(exampleUnidirectional, !is.na(geneID))

# Calculate composition: The number of samples expressing TSSs above 10% of
# the total gene expression.
intragenicTSSs <- calcComposition(intragenicTSSs, outputColumn = "composition", 
    unexpressed = 0.1, genes = "geneID")

# Overview of composition values:
table(rowRanges(intragenicTSSs)$composition)
#> 
#>   0   1   2   3 
#>  22   8   6 116

We can see that many TSSs makes up less than 10% in all samples. These can be removed with subset:

# Remove TSSs with a composition score less than 3
intragenicTSSs <- subset(intragenicTSSs, composition > 2)

The subsetByComposition function wraps the common task of calling calcComposition and subset. After filtering, the resulting expression matrix can be used with popular Bioconductor packages for DTU calling, such as DEXSeq, DRIMSeq, edgeR, limma, etc.