3.1 Getting started

We start the workflow by creating a CAGEexp object, which is a container for storing CAGE datasets and all the results that will be generated by applying specific functions. The CAGEexp objects are an extension of the MultiAssayExperiment class, and therefore can use all their methods. The expression data is stored in CAGEexp using SummarizedExperiment objects, and can also access their methods.

To load the CAGEr package and the other libraries into your R envirnoment type:

library("MultiAssayExperiment")
library("SummarizedExperiment")
library(CAGEr)

3.2 Creating a CAGEexp object

inputFiles = list.files( system.file("extdata", package = "CAGEr")
                       , "ctss$"
                       , full.names = TRUE)

ce <- CAGEexp( genomeName     = "BSgenome.Drerio.UCSC.danRer7"
             , inputFiles     = inputFiles
             , inputFilesType = "ctss"
             , sampleLabels   = sub( ".chr17.ctss", "", basename(inputFiles))
)

ce

## A CAGEexp object of 0 listed
##  experiments with no user-defined names and respective classes. 
##  Containing an ExperimentList class object of length 0:  
## Features: 
##  experiments() - obtain the ExperimentList instance 
##  colData() - the primary/phenotype DataFrame 
##  sampleMap() - the sample availability DataFrame 
##  `$`, `[`, `[[` - extract colData columns, subset, or experiment 
##  *Format() - convert into a long or wide DataFrame 
##  assays() - convert ExperimentList to a SimpleList of matrices

3.3 Reading in the data

In case when the CAGE / TSS data is to be read from input files, an empty CAGEexp object with information about the files is first created as described above in section 3.2. To actually read in the data into the object we use getCTSS() function, that will add an experiment called tagCountMatrix to the CAGEexp object.

getCTSS(ce)
ce

## A CAGEexp object of 1 listed
##  experiment with a user-defined name and respective class. 
##  Containing an ExperimentList class object of length 1: 
##  [1] tagCountMatrix: RangedSummarizedExperiment with 23343 rows and 5 columns 
## Features: 
##  experiments() - obtain the ExperimentList instance 
##  colData() - the primary/phenotype DataFrame 
##  sampleMap() - the sample availability DataFrame 
##  `$`, `[`, `[[` - extract colData columns, subset, or experiment 
##  *Format() - convert into a long or wide DataFrame 
##  assays() - convert ExperimentList to a SimpleList of matrices

This function reads the provided files in the order they were specified in the inputFiles argument. It creates a single set of all TSSs detected across all input datasets (union of TSSs) and a table with counts of CAGE tags supporting each TSS in every dataset. (Note that in case when a CAGEr object is created by coercion from an existing expression table there is no need to call getCTSS()).

Genomic coordinates of all TSSs and numbers of supporting CAGE tags in every input sample can be retrieved using the CTSStagCountSE() function. CTSScoordinatesGR() accesses the CTSS coordinates and CTSStagCountDF() accesses the CTSS expression values.11 Data can also be accessed directly using the native methods of the MultiAssayExperiment and SummarizedExperiment classes, for example ce[["tagCountMatrix"]], rowRanges(ce[["tagCountMatrix"]]) and assay(ce[["tagCountMatrix"]]).

CTSStagCountSE(ce)

## class: RangedSummarizedExperiment 
## dim: 23343 5 
## metadata(0):
## assays(1): counts
## rownames: NULL
## rowData names(0):
## colnames(5): Zf.30p.dome Zf.high Zf.prim6.rep1 Zf.prim6.rep2
##   Zf.unfertilized.egg
## colData names(0):

CTSScoordinatesGR(ce)

## CTSS object with 23343 positions and 0 metadata columns:
##           seqnames       pos strand
##              <Rle> <integer>  <Rle>
##       [1]    chr17  26027430      +
##       [2]    chr17  26050540      +
##       [3]    chr17  26118088      +
##       [4]    chr17  26142853      +
##       [5]    chr17  26166954      +
##       ...      ...       ...    ...
##   [23339]    chr17  45975041      -
##   [23340]    chr17  45975540      -
##   [23341]    chr17  45975544      -
##   [23342]    chr17  45982697      -
##   [23343]    chr17  45999921      -
##   -------
##   seqinfo: 26 sequences (1 circular) from danRer7 genome

CTSStagCountDF(ce)

## DataFrame with 23343 rows and 5 columns
##       Zf.30p.dome Zf.high Zf.prim6.rep1 Zf.prim6.rep2 Zf.unfertilized.egg
##             <Rle>   <Rle>         <Rle>         <Rle>               <Rle>
## 1               0       0             1             0                   0
## 2               0       0             0             0                   1
## 3               0       0             1             0                   0
## 4               0       0             0             1                   0
## 5               0       0             1             0                   0
## ...           ...     ...           ...           ...                 ...
## 23339           1       0             0             0                   0
## 23340           0       2             0             0                   0
## 23341           0       1             0             0                   0
## 23342           0       0             1             0                   0
## 23343           1       0             0             0                   0

For compatiblity with earlier CAGEr works using CAGEset objects, and to provide simpler data formats, the coordinates and expression values can also be accessed as simple data.frames. Note howerver that with large data sets it can cause extreme performance issues.

head(CTSScoordinates(ce))

##     chr      pos strand
## 1 chr17 26027430      +
## 2 chr17 26050540      +
## 3 chr17 26118088      +
## 4 chr17 26142853      +
## 5 chr17 26166954      +
## 6 chr17 26222417      +

head(CTSStagCountDf(ce))

##   Zf.30p.dome Zf.high Zf.prim6.rep1 Zf.prim6.rep2 Zf.unfertilized.egg
## 1           0       0             1             0                   0
## 2           0       0             0             0                   1
## 3           0       0             1             0                   0
## 4           0       0             0             1                   0
## 5           0       0             1             0                   0
## 6           1       1             0             0                   0

head(CTSStagCount(ce))

##     chr      pos strand Zf.30p.dome Zf.high Zf.prim6.rep1 Zf.prim6.rep2
## 1 chr17 26027430      +           0       0             1             0
## 2 chr17 26050540      +           0       0             0             0
## 3 chr17 26118088      +           0       0             1             0
## 4 chr17 26142853      +           0       0             0             1
## 5 chr17 26166954      +           0       0             1             0
## 6 chr17 26222417      +           1       1             0             0
##   Zf.unfertilized.egg
## 1                   0
## 2                   1
## 3                   0
## 4                   0
## 5                   0
## 6                   0

Note that the samples are ordered in the way they were supplied when creating the CAGEexp object and will be presented in that order in all the results and plots. To check sample labels and their ordering type:

sampleLabels(ce)

##             #FF0000FF             #CCFF00FF             #00FF66FF 
##         "Zf.30p.dome"             "Zf.high"       "Zf.prim6.rep1" 
##             #0066FFFF             #CC00FFFF 
##       "Zf.prim6.rep2" "Zf.unfertilized.egg"

In addition, a colour is assigned to each sample, which is consistently used to depict that sample in all the plots. By default a rainbow palette of colours is used and the hexadecimal format of the assigned colours can be seen as names attribute of sample labels shown above. The colours can be changed to taste at any point in the workflow using the setColors() function.

3.4 Quality controls and preliminary analyses

annotateCTSS(ce, exampleZv9_annot)

colData(ce)[,c("librarySizes", "promoter", "exon", "intron", "unknown")]

## DataFrame with 5 rows and 5 columns
##                     librarySizes  promoter      exon    intron   unknown
##                        <integer> <integer> <integer> <integer> <integer>
## Zf.30p.dome                41814     37843      2352       594      1025
## Zf.high                    45910     41671      2848       419       972
## Zf.prim6.rep1              34053     29531      2714       937       871
## Zf.prim6.rep2              34947     30799      2320       834       994
## Zf.unfertilized.egg        56140     51114      2860       400      1766

CTSScoordinatesGR(ce)

## CTSS object with 23343 positions and 2 metadata columns:
##           seqnames       pos strand |  genes annotation
##              <Rle> <integer>  <Rle> |  <Rle>      <Rle>
##       [1]    chr17  26027430      + |           unknown
##       [2]    chr17  26050540      + | grid1a   promoter
##       [3]    chr17  26118088      + | grid1a       exon
##       [4]    chr17  26142853      + | grid1a     intron
##       [5]    chr17  26166954      + | grid1a       exon
##       ...      ...       ...    ... .    ...        ...
##   [23339]    chr17  45975041      - |           unknown
##   [23340]    chr17  45975540      - |           unknown
##   [23341]    chr17  45975544      - |           unknown
##   [23342]    chr17  45982697      - |           unknown
##   [23343]    chr17  45999921      - |           unknown
##   -------
##   seqinfo: 26 sequences (1 circular) from danRer7 genome

plotAnnot(ce, "counts")

## Warning: Removed 20 rows containing missing values (geom_segment).

## Warning: Removed 20 rows containing missing values (geom_point).

3.4.2 Correlation between samples

As part of the basic sanity checks, we can explore the data by looking at the correlation between the samples. The plotCorrelation2() function will plot pairwise scatter plots of expression scores per TSS or consensus cluster and calculate correlation coefficients between all possible pairs of samples22 Alternatively, the plotCorrelation() function does the same and colors the scatterplots according to point density, but is much slower.. A threshold can be set, so that only regions with an expression score (raw or normalized) above the threshold (either in one or both samples) are considered when calculating correlation. Three different correlation measures are supported: Pearson’s, Spearman’s and Kendall’s correlation coefficients. Note that while the scatterplots are on a logarithmic scale with pseudocount added to the zero values, the correlation coefficients are calculated on untransformed (but thresholded) data.

corr.m <- plotCorrelation2( ce, samples = "all"
                          , tagCountThreshold = 1, applyThresholdBoth = FALSE
                          , method = "pearson")

Figure 2: Correlation of raw CAGE tag counts per TSS

3.5 Merging of replicates

Based on calculated correlation we might want to merge and/or rearrange some of the datasets. To rearrange the samples in the temporal order of the zebrafish development (unfertilized egg -> high -> 30 percent dome -> prim6) and to merge the two replicas for the prim6 developmental stage we use the mergeSamples() function:

mergeSamples(ce, mergeIndex = c(3,2,4,4,1), 
            mergedSampleLabels = c("zf_unfertilized_egg", "zf_high", "zf_30p_dome", "zf_prim6"))
annotateCTSS(ce, exampleZv9_annot)

The mergeIndex argument controls which samples will be merged and how the final dataset will be ordered. Samples labeled by the same number (in our case samples three and four) will be merged together by summing number of CAGE tags per TSS. The final set of samples will be ordered in the ascending order of values provided in mergeIndex and will be labeled by the labels provided in the mergedSampleLabels argument. Note that mergeSamples function resets all slots with results of downstream analyses, so in case there were any results in the CAGEexp object prior to merging, they will be removed. Thus, annotation has to be redone.

3.6 Normalization

Library sizes (number of total sequenced tags) of individual experiments differ, thus normalization is required to make them comparable. The librarySizes function returns the total number of CAGE tags in each sample:

librarySizes(ce)

## [1] 56140 45910 41814 69000

The CAGEr package supports both simple tags per million normalization and power-law based normalization. It has been shown that many CAGE datasets follow a power-law distribution (Balwierz et al. 2009). Plotting the number of CAGE tags (X-axis) against the number of TSSs that are supported by <= of that number of tags (Y-axis) results in a distribution that can be approximated by a power-law. On a log-log scale this reverse cumulative distribution will manifest as a monotonically decreasing linear function, which can be defined as

\[y = -1 * \alpha * x + \beta\]

and is fully determined by the slope \(\alpha\) and total number of tags T (which together with \(\alpha\) determines the value of \(\beta\)).

To check whether our CAGE datasets follow power-law distribution and in which range of values, we can use the plotReverseCumulatives function:

plotReverseCumulatives(ce, fitInRange = c(5, 1000), onePlot = TRUE)

Figure 3: Reverse cumulative distribution of CAGE tags

In addition to the reverse cumulative plots (Figure 3), a power-law distribution will be fitted to each reverse cumulative using values in the specified range (denoted with dashed lines in Figure 3) and the value of \(\alpha\) will be reported for each sample (shown in the brackets in the Figure 3 legend). The plots can help in choosing the optimal parameters for power-law based normalization. We can see that the reverse cumulative distributions look similar and follow the power-law in the central part of the CAGE tag counts values with a slope between -1.1 and -1.3. Thus, we choose a range from 5 to 1000 tags to fit a power-law, and we normalize all samples to a referent power-law distribution with a total of 50,000 tags and slope of -1.2 (\(\alpha = 1.2\)).33 Note that since this example dataset contains only data from one part of chromosome 17 and the total number of tags is very small, we normalize to a referent distribution with a similarly small number of tags. When analyzing full datasets it is reasonable to set total number of tags for referent distribution to one million to get normalized tags per million values.

To perform normalization we pass these parameters to the normalizeTagCount function.

normalizeTagCount(ce, method = "powerLaw", fitInRange = c(5, 1000), alpha = 1.2, T = 5*10^4)
ce[["tagCountMatrix"]]

## class: RangedSummarizedExperiment 
## dim: 23343 4 
## metadata(0):
## assays(2): counts normalizedTpmMatrix
## rownames: NULL
## rowData names(2): genes annotation
## colnames(4): zf_unfertilized_egg zf_high zf_30p_dome zf_prim6
## colData names(0):

The normalization is performed as described in (Balwierz et al. 2009):

• Power-law is fitted to the reverse cumulative distribution in the specified range of CAGE tags values to each sample separately.
• A referent power-law distribution is defined based on the provided alpha (slope in the log-log representation) and T (total number of tags) parameters. Setting T to 1 million results in normalized tags per million (tpm) values.
• Every sample is normalized to the defined referent distribution, i.e. given the parameters that approximate its own power-law distribution it is calculated how many tags would each TSS have in the referent power-law distribution.

In addition to the two provided normalization methods, a pass-through option none can be set as method parameter to keep using raw tag counts in all downstream steps. Note that normalizeTagCount() has to be applied to CAGEr object before moving to next steps. Thus, in order to keep using raw tag counts run the function with method="none". In that case, all results and parameters in the further steps that would normally refer to normalized CAGE signal (denoted as tpm), will actually be raw tag counts.

3.7 CTSS clustering

Transcription start sites are found in the promoter region of a gene and reflect the transcriptional activity of that promoter (Figure 5). TSSs in the close proximity of each other give rise to a functionally equivalent set of transcripts and are likely regulated by the same promoter elements. Thus, TSSs can be spatially clustered into larger transcriptional units, called tag clusters (TCs) that correspond to individual promoters. CAGEr supports three methods for spatial clustering of TSSs along the genome, two ab initio methods driven by the data itself, as well as assigning TSSs to predefined genomic regions:

• Simple distance-based clustering in which two neighbouring TSSs are joined together if they are closer than some specified distance (greedy algorithm);

• Parametric clustering of data attached to sequences based on the density of the signal (Frith et al. 2007), http://www.cbrc.jp/paraclu/;

• Counting TSSs and their signal in a set of user supplied genomic regions (e.g. annotation derived promoter regions or other regions of interest).

These functionalities are provided in the function clusterCTSS(), which accepts additional arguments for controlling which CTSSs will be included in the clustering as well as for refining the final set of tag clusters.

We will perform a simple distance-based clustering using 20 bp as a maximal allowed distance between two neighbouring TSSs. Prior to clustering we will filter out low-fidelity TSSs - the ones supported by less than 1 normalized tag counts in all of the samples.

clusterCTSS( object = ce
           , threshold = 1
           , thresholdIsTpm = TRUE
           , nrPassThreshold = 1
           , method = "distclu"
           , maxDist = 20
           , removeSingletons = TRUE
           , keepSingletonsAbove = 5)

Our final set of tag clusters will not include singletons (clusters with only one TSS), unless the normalized signal is above 5, it is a reasonably supported TSS. The clusterCTSS function creates a set of clusters for each sample separately; for each cluster it returns the genomic coordinates, counts the number of TSSs within the cluster, determines the position of the most frequently used (dominant) TSS, calculates the total CAGE signal within the cluster and CAGE signal supporting the dominant TSS only. We can extract tag clusters for a desired sample from CAGEexp object by calling the tagClusters function:

head(tagClusters(ce, sample = "zf_unfertilized_egg"))

##   cluster   chr    start      end strand   tpm nr_ctss dominant_ctss
## 1       1 chr17 26453631 26453708      +  27.0      12      26453667
## 2       2 chr17 26564507 26564610      + 128.6      24      26564585
## 3       3 chr17 26595636 26595793      + 217.0      35      26595750
## 4       4 chr17 26596032 26596091      +  10.4       9      26596070
## 5       5 chr17 26596117 26596127      +  12.2       4      26596118
## 6       6 chr17 26596149 26596175      +  11.0       5      26596153
##   tpm.dominant_ctss
## 1              8.25
## 2             29.28
## 3            100.97
## 4              3.22
## 5              5.74
## 6              3.85

3.8 Promoter width

Genome-wide mapping of TSSs using CAGE has initially revealed two major classes of promoters in mammals (Carninci et al. 2006), with respect to the number and distribution of TSSs within the promoter. They have been further correlated with differences in the underlying sequence and the functional classes of the genes they regulate, as well as the organization of the chromatin around them.

• “broad” promoters with multiple TSSs characterized by a high GC content and overlap with a CpG island, which are associated with widely expressed or developmentally regulated genes;
• “sharp” promoters with one dominant TSS often associated with a TATA-box at a fixed upstream distance, which often regulate tissue-specific transcription.

Thus, the width of the promoter is an important characteristic that distinguishes different functional classes of promoters. CAGEr analyzes promoter width across all samples present in the CAGEexp object. It defines promoter width by taking into account both the positions and the CAGE signal at TSSs along the tag cluster, thus making it more robust with respect to total expression and local level of noise at the promoter. Width of every tag cluster is calculated as following:

1. Cumulative distribution of CAGE signal along the cluster is calculated.
2. A “lower” (qLow) and an “upper” (qUp) quantile are selected by the user.
3. From the 5′ end the position, the position of a quantile \(q\) is detemined as the first base in which of the cumulative expression is higher or equal to \(q\%\) of the total CAGE signal of that cluster.
4. Promoter interquantile width is defined as the distance (in base pairs) between the two quantile positions.

The procedure is schematically shown in Figure 4.

Figure 4: Cumulative distribution of CAGE signal and definition of interquantile width

Required computations are done using cumulativeCTSSdistribution and quantilePositions functions, which calculate cumulative distribution for every tag cluster in each of the samples and determine the positions of selected quantiles, respectively:

cumulativeCTSSdistribution(ce, clusters = "tagClusters", useMulticore = T)
quantilePositions(ce, clusters = "tagClusters", qLow = 0.1, qUp = 0.9)

Tag clusters and their interquantile width can be retrieved by calling tagClusters function:

tagClustersGR( ce, "zf_unfertilized_egg"
             , returnInterquantileWidth = TRUE,  qLow = 0.1, qUp = 0.9)

## TagClusters object with 517 ranges and 7 metadata columns:
##       seqnames            ranges strand |  score   nr_ctss dominant_ctss
##          <Rle>         <IRanges>  <Rle> |  <Rle> <integer>     <integer>
##     1    chr17 26453632-26453708      + |   27.0        12      26453667
##     2    chr17 26564508-26564610      + |  128.6        24      26564585
##     3    chr17 26595637-26595793      + |  217.0        35      26595750
##     4    chr17 26596033-26596091      + |   10.4         9      26596070
##     5    chr17 26596118-26596127      + |   12.2         4      26596118
##   ...      ...               ...    ... .    ...       ...           ...
##   513    chr17 45700182-45700187      - |   9.62         3      45700182
##   514    chr17 45901329-45901330      - |   1.96         2      45901329
##   515    chr17 45901698-45901710      - |  27.65         4      45901698
##   516    chr17 45901732-45901784      - | 119.97        15      45901749
##   517    chr17 45901814-45901816      - |   3.25         2      45901816
##       tpm.dominant_ctss q_0.1 q_0.9 interquantile_width
##                   <Rle> <Rle> <Rle>               <Rle>
##     1              8.25    36    72                  37
##     2             29.28    17    81                  65
##     3            100.97    37   114                  78
##     4              3.22     1    50                  50
##     5              5.74     1    10                  10
##   ...               ...   ...   ...                 ...
##   513              6.37     1     6                   6
##   514              1.30     1     2                   2
##   515             23.75     1     2                   2
##   516             83.45     2    21                  20
##   517              1.94     1     3                   3
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

Once the cumulative distributions and the positions of quantiles have been calculated, the histograms of interquantile width can be plotted to globally compare the promoter width across different samples (Figure ??:

plotInterquantileWidth(ce, clusters = "tagClusters", tpmThreshold = 3, qLow = 0.1, qUp = 0.9)

Significant difference in the promoter width might indicate global differences in the modes of gene regulation between the two samples. The histograms can also help in choosing an appropriate width threshold for separating sharp and broad promoters.

3.9 Creating consensus promoters across samples

Tag clusters are created for each sample individually and they are often sample-specific, thus can be present in one sample but absent in another. In addition, in many cases tag clusters do not coincide perfectly within the same promoter region, or there might be two clusters in one sample and only one larger in the other. To be able to compare genome-wide transcriptional activity across samples and to perform expression profiling, a single set of consensus clusters needs to be created. This is done using the aggregateTagClusters function, which aggregates tag clusters from all samples into a single set of non-overlapping consensus clusters:

aggregateTagClusters(ce, tpmThreshold = 5, qLow = 0.1, qUp = 0.9, maxDist = 100)
ce$outOfClusters / ce$librarySizes *100

## zf_unfertilized_egg             zf_high         zf_30p_dome 
##                33.3                32.9                32.1 
##            zf_prim6 
##                34.4

Tag clusters can be aggregated using their full span (from start to end) or using positions of previously calculated quantiles as their boundaries. Only tag clusters above given tag count threshold will be considered and two clusters will be aggregated together if their boundaries (i.e. either starts and ends or positions of quantiles) are <= maxDist apart. Final set of consensus clusters can be retrieved by:

consensusClustersGR(ce)

## ConsensusClusters object with 253 ranges and 3 metadata columns:
##                             seqnames            ranges strand |
##                                <Rle>         <IRanges>  <Rle> |
##   chr17:26379388-26379679:-    chr17 26379388-26379679      - |
##   chr17:26446559-26446651:-    chr17 26446559-26446651      - |
##   chr17:26452451-26452571:-    chr17 26452451-26452571      - |
##   chr17:26453648-26453757:+    chr17 26453648-26453757      + |
##   chr17:26555294-26555359:-    chr17 26555294-26555359      - |
##                         ...      ...               ...    ... .
##   chr17:45700093-45700189:-    chr17 45700093-45700189      - |
##   chr17:45806843-45806894:+    chr17 45806843-45806894      + |
##   chr17:45901696-45901758:-    chr17 45901696-45901758      - |
##   chr17:45905581-45905664:+    chr17 45905581-45905664      + |
##   chr17:45975253-45975297:+    chr17 45975253-45975297      + |
##                                        score consensus.cluster
##                                    <numeric>         <integer>
##   chr17:26379388-26379679:- 179.803895280803               196
##   chr17:26446559-26446651:- 48.1654641622959               197
##   chr17:26452451-26452571:- 184.204291454949               198
##   chr17:26453648-26453757:+ 129.733525462687                 1
##   chr17:26555294-26555359:- 19.2388260848655               199
##                         ...              ...               ...
##   chr17:45700093-45700189:- 9.21466887879838               368
##   chr17:45806843-45806894:+  108.34520829738               190
##   chr17:45901696-45901758:-  789.81397211728               372
##   chr17:45905581-45905664:+ 321.125346126617               192
##   chr17:45975253-45975297:+ 35.8163066635267               193
##                                          tpm
##                                    <numeric>
##   chr17:26379388-26379679:- 287.247912339468
##   chr17:26446559-26446651:- 85.7757116833182
##   chr17:26452451-26452571:- 269.689488829474
##   chr17:26453648-26453757:+ 185.967802102194
##   chr17:26555294-26555359:- 47.6676600011559
##                         ...              ...
##   chr17:45700093-45700189:- 23.2038219025773
##   chr17:45806843-45806894:+ 173.614212890553
##   chr17:45901696-45901758:- 875.831472947622
##   chr17:45905581-45905664:+ 392.236067622669
##   chr17:45975253-45975297:+ 55.6241013192924
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

which will return genomic coordinates and sum of CAGE signal across all samples for each consensus cluster (the tpm column).

Analogously to tag clusters, analysis of promoter width can be performed for consensus clusters as well, using the same cumulativeCTSSdistribution, quantilePositions and plotInterquantileWidth functions as described above, but by setting the clusters parameter to "consensusClusters". Like the CTSS, the consensus clusters can also be annotated:

annotateConsensusClusters(ce, exampleZv9_annot)
cumulativeCTSSdistribution(ce, clusters = "consensusClusters", useMulticore = TRUE)
quantilePositions(ce, clusters = "consensusClusters", qLow = 0.1, qUp = 0.9, useMulticore = TRUE)

Although consensus clusters are created to represent consensus across all samples, they obviously have different CAGE signal and can have different width or position of the dominant TSS in the different samples. Sample-specific information on consensus clusters can be retrieved with the function, by specifying desired sample name (analogous to retrieving tag clusters):

consensusClustersGR( ce, sample = "zf_unfertilized_egg"
                       , returnInterquantileWidth = TRUE,  qLow = 0.1, qUp = 0.9)

## ConsensusClusters object with 253 ranges and 8 metadata columns:
##                             seqnames            ranges strand |
##                                <Rle>         <IRanges>  <Rle> |
##   chr17:26379388-26379679:-    chr17 26379388-26379679      - |
##   chr17:26446559-26446651:-    chr17 26446559-26446651      - |
##   chr17:26452451-26452571:-    chr17 26452451-26452571      - |
##   chr17:26453648-26453757:+    chr17 26453648-26453757      + |
##   chr17:26555294-26555359:-    chr17 26555294-26555359      - |
##                         ...      ...               ...    ... .
##   chr17:45700093-45700189:-    chr17 45700093-45700189      - |
##   chr17:45806843-45806894:+    chr17 45806843-45806894      + |
##   chr17:45901696-45901758:-    chr17 45901696-45901758      - |
##   chr17:45905581-45905664:+    chr17 45905581-45905664      + |
##   chr17:45975253-45975297:+    chr17 45975253-45975297      + |
##                                        score consensus.cluster
##                                    <numeric>         <integer>
##   chr17:26379388-26379679:-  49.514684740936               196
##   chr17:26446559-26446651:- 10.9138520447868               197
##   chr17:26452451-26452571:- 61.8658559845496               198
##   chr17:26453648-26453757:+ 19.1354399087033                 1
##   chr17:26555294-26555359:- 6.37015756741174               199
##                         ...              ...               ...
##   chr17:45700093-45700189:- 6.37015756741174               368
##   chr17:45806843-45806894:+ 23.0404872540061               190
##   chr17:45901696-45901758:- 132.654636120977               372
##   chr17:45905581-45905664:+   106.7961601026               192
##   chr17:45975253-45975297:+                0               193
##                                          tpm annotation           genes
##                                    <numeric>      <Rle>           <Rle>
##   chr17:26379388-26379679:-  49.514684740936   promoter CCSER2 (1 of 2)
##   chr17:26446559-26446651:- 10.9138520447868   promoter        fam149b1
##   chr17:26452451-26452571:- 61.8658559845496   promoter           mrp63
##   chr17:26453648-26453757:+ 19.1354399087033   promoter           ttc7b
##   chr17:26555294-26555359:- 6.37015756741174       exon          calm1a
##                         ...              ...        ...             ...
##   chr17:45700093-45700189:- 6.37015756741174       exon           susd4
##   chr17:45806843-45806894:+ 23.0404872540061   promoter          vcpkmt
##   chr17:45901696-45901758:- 132.654636120977   promoter           arf6b
##   chr17:45905581-45905664:+   106.7961601026   promoter            aspg
##   chr17:45975253-45975297:+                0    unknown                
##                             q_0.1 q_0.9 interquantile_width
##                             <Rle> <Rle>               <Rle>
##   chr17:26379388-26379679:-    20    56                  37
##   chr17:26446559-26446651:-    31    65                  35
##   chr17:26452451-26452571:-    48    77                  30
##   chr17:26453648-26453757:+    37   271                 235
##   chr17:26555294-26555359:-     3   109                 107
##                         ...   ...   ...                 ...
##   chr17:45700093-45700189:-     1     1                   1
##   chr17:45806843-45806894:+    26    74                  49
##   chr17:45901696-45901758:-    25    31                   7
##   chr17:45905581-45905664:+     3    92                  90
##   chr17:45975253-45975297:+     3    54                  52
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

This will, in addition to genomic coordinates of the consensus clusters, which are constant across all samples, also return the position of the dominant TSS, the CAGE signal (tpm) and the interquantile width specific for a given sample. Note that when specifying individual sample, only the consensus clusters that have some CAGE signal in that sample will be returned (which will be a subset of all consensus clusters).

3.10 Data export for genome browsers

CAGE data can be visualized in the genomic context by exporting raw or normalized CAGE tag counts to a bedGraph (or BigWig) file and uploading (or linking) the file to a genome browser. Positions of TSSs and tag counts supporting them are exported using exportCTSStoBedGraph()44 Note that the ZENBU genome browser can display natively data from BAM or BED files as coverage tracks.:

exportCTSStoBedGraph(ce, values = "normalized", format = "bedGraph", oneFile = TRUE)

This will produce a single bedGraph file with multiple annotated tracks that can be directly visualized as custom tracks in the genome browser (Figure 5).

There are two tracks per sample; one for TSSs on the plus strand and the other for the minus strand. Values for TSSs on minus strand are shown as negative and are pointing downwards in the browser.

Alternatively, the tracks can be exported to a binary BigWig format:

exportCTSStoBedGraph(ce, values = "normalized", format = "BigWig")

which will produce two BigWig files per sample (one for TSSs on the plus strand and the other for the minus strand) and an accompanying text file with track headers.

Figure 5: CAGE data bedGraph track visualized in the UCSC Genome Browser

Interquantile width can also be visualized in a gene-like representation in the UCSC genome browser by exporting the data into a BED file:

exportToBed(object = ce, what = "tagClusters", qLow = 0.1, qUp = 0.9, oneFile = TRUE)

In this gene-like representation (Figure 6), the oriented line shows the full span of the cluster, filled block marks the interquantile width and a single base-pair thick block denotes the position of the dominant TSS.

Figure 6: Tag clusters visualization in the genome browser

3.11 Expression profiling

Since CAGE signal reflects level of transcription from a given TSS or promoter it can be used for 5′ centered expression profiling. Expression clustering can be done at level of individual CTSSs or at level of entire promoters (consensus clusters). In the former case, feature vector containing log transformed and scaled normalized CAGE signal at individual TSS across multiple samples is used as input for clustering algorithm, whereas in the latter case CAGE signal within the entire consensus cluster is used. CAGEr supports two unsupervised clustering algorithms: kmeans and self-organizing maps (SOM). Both algorithms require to specify a number of clusters in advance.

We will perform expression clustering at the level of entire promoter using SOM algorithm and applying it only to promoters with normalized CAGE signal \(>= 15\) in at least one sample.

# Not supported yet for CAGEexp objects, sorry.
# getExpressionProfiles(ce, what = "consensusClusters", tpmThreshold = 10, 
#       nrPassThreshold = 1, method = "som", xDim = 4, yDim = 2)

Distribution of expression across samples for 8 clusters returned by SOM (4 \(\times\) 2 map) can be visualized using the plotExpressionProfiles function as shown in Figure 7:

# Not supported yet for CAGEexp objects, sorry.
# plotExpressionProfiles(ce, what = "consensusClusters")

Figure 7: Expression clusters

Each cluster is shown in different color and is marked by its label and the number of elements (promoters) in the cluster. We can extract promoters belonging to a specific cluster by typing:

# Not supported yet for CAGEexp objects, sorry.
# class3_1 <- extractExpressionClass(ce, 
#       what = "consensusClusters", which = "3_1")
# head(class3_1)

Consensus clusters and information on their expression profile can be exported to a BED file, which allows visualization of the promoters in the genome browser colored in the color of the expression cluster they belong to (Figure 8:

# Not supported yet for CAGEexp objects, sorry.
# exportToBed(ce, what = "consensusClusters", 
#       colorByExpressionProfile = TRUE)

Figure 8: Consensus clusters colored by expression profile in the genome browser

Expression profiling of individual TSSs is done using the same procedure as described above for consensus clusters, only by setting wha = "CTSS" in all of the functions.

3.12 Differential expression analysis

The raw expression table for the consensus clusters can be exported to the DESeq2 package for differential expression analysis. For this, the column data needs to contain factors that can group the samples. They can have any name.

ce$group <- factor(c("a", "a", "b", "b"))
dds <- consensusClustersDESeq2(ce, ~group)

3.13 Shifting promoters

As shown in Figure 6, TSSs within the same promoter region can be used differently in different samples. Thus, although the overall transcription level from a promoter does not change between the samples, the differential usage of TSSs or promoter shifting may indicate changes in the regulation of transcription from that promoter, which cannot be detected by expression profiling. To detect this promoter shifting, a method described in @[Haberle:2014] has been implemented in CAGEr. Shifting can be detected between two individual samples or between two groups of samples. In the latter case, samples are first merged into groups and then compared in the same way as two individual samples. For all promoters a shifting score is calculated based on the difference in the cumulative distribution of CAGE signal along that promoter in the two samples. In addition, a more general assessment of differential TSS usage is obtained by performing Kolmogorov-Smirnov test on the cumulative distributions of CAGE signal, as described below. Thus, prior to shifting score calculation and statistical testing, we have to calculate cumulative distribution along all consensus clusters:

cumulativeCTSSdistribution(ce, clusters = "consensusClusters")

Next, we calculate a shifting score and P-value of Kolmogorov-Smirnov test for all promoters comparing two specified samples:

# Not supported yet for CAGEexp objects, sorry.
# scoreShift(ce, groupX = "zf_unfertilized_egg", groupY = "zf_prim6",
#       testKS = TRUE, useTpmKS = FALSE)

This function will calculate shifting score as illustrated in Figure 9. Values of shifting score are in range between -Inf and 1. Positive values can be interpreted as the proportion of transcription initiation in the sample with lower expression that is happening “outside” (either upstream or downstream) of the region used for transcription initiation in the other sample. In contrast, negative values indicate no physical separation, i.e. the region used for transcription initiation in the sample with lower expression is completely contained within the region used for transcription initiation in the other sample. Thus, shifting score detects only the degree of upstream or downstream shifting, but does not detect more general changes in TSS rearrangement in the region, e.g. narrowing or broadening of the region used for transcription.

To assess any general change in the TSS usage within the promoter region, a two-sample Kolmogorov-Smirnov (K-S) test on cumulative sums of CAGE signal along the consensus cluster is performed. Cumulative sums in both samples are scaled to range between 0 and 1 and are considered to be empirical cumulative distribution functions (ECDF) reflecting sampling of TSS positions during transcription initiation. K-S test is performed to assess whether the two underlying probability distributions differ. To obtain a P-value i.e. the level at which the null-hypothesis can be rejected), sample sizes that generated the ECDFs are required, in addition to actual K-S statistics calculated from ECDFs. These are derived either from raw tag counts, i.e. exact number of times each TSS in the cluster was sampled during sequencing (when useTpmKS = FALSE), or from normalized tpm values (when useTpmKS = TRUE). P-values obtained from K-S tests are further corrected for multiple testing using Benjamini and Hochenberg (BH) method and for each P-value a corresponding false-discovery rate (FDR) is also reported.

Figure 9: Calculation of shifting score

We can select a subset of promoters with shifting score and/or FDR above specified threshold:

# Not supported yet for CAGEexp objects, sorry.
# shifting.promoters <- getShiftingPromoters(ce, 
#       tpmThreshold = 5, scoreThreshold = 0.6,
#       fdrThreshold = 0.01)
# head(shifting.promoters)

The getShiftingPromoters function returns genomic coordinates, shifting score and P-value (FDR) of the promoters, as well as the value of CAGE signal and position of the dominant TSS in the two compared (groups of) samples. Figure 10 shows the difference in the CAGE signal between the two compared samples for one of the selected high-scoring shifting promoters.

Figure 10: Example of shifting promoter

4.1 Creating a CAGEexp object by coercing a data frame

A CAGEexp object can also be created directly by coercing a data frame containing single base-pair TSS information. To be able to do the coercion into a CAGEexp, the data frame must conform with the following:

• The data frame must have at least 4 columns;

• the first three columns must be named chr, pos and strand, and contain chromosome name, 1-based genomic coordinate of the TSS (positive integer) and TSS strand information (+ or -), respectively;

• these first three columns must be of the class character, integer and character, respectively;

• all additional columns must be of the class integer and should contain raw CAGE tag counts (non-negative integer) supporting each TSS in different samples (columns). At least one such column with tag counts must be present;

• the names of the columns containing tag counts must begin with a letter, and these column names are used as sample labels in the resulting CAGEexp object.

An example of such data frame is shown below:

TSS.df <- read.table(system.file( "extdata/Zf.unfertilized.egg.chr17.ctss"
                                , package = "CAGEr"))
# make sure the column names are as required
colnames(TSS.df) <- c("chr", "pos", "strand", "zf_unfertilized_egg")
# make sure the column classes are as required
TSS.df$chr <- as.character(TSS.df$chr)
TSS.df$pos <- as.integer(TSS.df$pos)
TSS.df$strand <- as.character(TSS.df$strand)
TSS.df$zf_unfertilized_egg <- as.integer(TSS.df$zf_unfertilized_egg)
head(TSS.df)

##     chr      pos strand zf_unfertilized_egg
## 1 chr17 26050540      +                   1
## 2 chr17 26074127      -                   2
## 3 chr17 26074129      -                   3
## 4 chr17 26222545      -                   1
## 5 chr17 26322780      -                   1
## 6 chr17 26322832      -                   2

This data.frame can now be coerced to a CAGEexp object, which will fill the corresponding slots of the object with provided TSS information:

ce.coerced <- as(TSS.df, "CAGEexp")
ce.coerced

## A CAGEexp object of 1 listed
##  experiment with a user-defined name and respective class. 
##  Containing an ExperimentList class object of length 1: 
##  [1] tagCountMatrix: RangedSummarizedExperiment with 8395 rows and 1 columns 
## Features: 
##  experiments() - obtain the ExperimentList instance 
##  colData() - the primary/phenotype DataFrame 
##  sampleMap() - the sample availability DataFrame 
##  `$`, `[`, `[[` - extract colData columns, subset, or experiment 
##  *Format() - convert into a long or wide DataFrame 
##  assays() - convert ExperimentList to a SimpleList of matrices

4.2 Summary of the CAGEr accessor functions

Originally there was one accessor per slot in CAGEset objects (the original CAGEr format), but now that I added the CAGEexp class, that uses different underlying formats, the number of accessors increased because a) I provide the old ones for backwards compatibility and b) there multiple possible output formats.

Before releasing this CAGEr update in Bioconductor, I would like to be sure that the number of accessors and the name scheme are good enough.

Please let me know your opinion about the names and scope of the accessors below:

4.3 Summary of the CAGEexp experiment slots and assays

A CAGEexp object can contain the following experiments.

Please let me know your opinion about the names

4.4 Summary of the CAGEr classes

The CTSS, CTSS.chr, TagCluster and ConsensusClsuters are mostly used internally or type safety and preventing me (Charles) from mixing up inputs. They are visible from the outside. Should they be used more extensively ? Can they be replaced by more “core” Bioconductor classes ?