Zhaolian Lu, Keenan Berry, Zhenbin Hu, Yu Zhan, Tae-Hyuk Ahn, Zhenguo Lin

TSSr was published in NAR Genomics and Bioinformatics in 2021, please cite:

Lu, Z., Berry, K., Hu, Z., Zhan, Y., Ahn, T., & Lin, Z. (2021). TSSr: an R package for comprehensive analyses of TSS sequencing data. NAR Genomics and Bioinformatics, 3.

Workflow functions now return the modified TSSr object and no longer modify the caller's object implicitly. Existing scripts must assign the return value:

# Before 0.99.17 (no longer updates myTSSr)
mergeSamples(myTSSr)

# TSSr 0.99.17 and later
myTSSr <- mergeSamples(myTSSr)

The same assignment pattern applies to getTSS(), normalizeTSS(), filterTSS(), clusterTSS(), consensusCluster(), shapeCluster(), annotateCluster(), deGene(), callEnhancer(), and shiftPromoter(). This change replaces the previous non-standard caller-environment mutation with normal R value semantics.

TSSr is an R package that provides rich functions for mapping transcription start sites (TSSs) and characterizations of structures and activities of core promoters based on various types of sequencing data generated by high-throughput techniques that sequence the very 5' end of RNA transcripts (TSS sequencing), such as cap analysis of gene expression (CAGE), nAnT-iCAGE, NanoCAGE, transcript leader sequencing (TL-seq), transcript isoform sequencing (TIF-seq), TSS-seq, RAMPAGE, single-cell tagged reverse transcription (STRT), global nuclear run-on cap (GRO-cap), and STRIPE-seq.

The main functions of TSSr include identification of TSSs, normalization and filtration of TSS signals, clustering of TSSs to infer core promoters, quantification of core promoter activities and shape, assignment of TSS clusters to genes, gene differential expression, and core promoter shift. TSSr accepts multiple input formats, such as Binary Sequence Alignment Map (BAM) files (single-ended or paired-ended), Browser Extension Data (BED) files, BigWig files, and TSS tables. TSSr generates various types of TSS or core promoter track files that can be visualized in the UCSC Genome Browser or Integrative Genomics Viewer (IGV). TSSr also exports raw and processed result tables and visualization plots. Multi-core processing is supported on Linux and macOS platforms.

• Download R (>= 4.6.0) from CRAN.
  • Windows: https://cran.r-project.org/bin/windows/base/
  • Mac OS X: https://cran.r-project.org/bin/macosx/
  • Linux: https://cran.r-project.org/bin/linux/

• Rsamtools

• GenomicRanges

• GenomicFeature

• Gviz

• rtracklayer

• DESeq2

• BSgenome

  • install the packages by using the following R command:

    install.packages("devtools")
    if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
    BiocManager::install("Rsamtools")
    BiocManager::install("GenomicRanges")
    BiocManager::install("GenomicFeatures")
    BiocManager::install("Gviz")
    BiocManager::install("rtracklayer")
    BiocManager::install("DESeq2")
    BiocManager::install("BSgenome")

• data.table

• stringr

  • install the packages by using the following R command:

  install.packages("data.table")
  install.packages("stringr")

Once TSSr is released through Bioconductor, install with:

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("TSSr")

To install the development version directly from GitHub:

devtools::install_github("Linlab-slu/TSSr", build_vignettes = TRUE)

TSSr accepts two types of input data: read alignment files or TSS tables. The read alignment files in compressed binary alignment map (BAM) format are required if users intend to identify TSSs from mapped sequencing data. BAM files can be derived from mapping of either paired-end or single-end TSS sequencing reads. A TSS table is a matrix of TSS data generated by TSSr or other bioinformatics tools. A TSS table can be a tab-delimited text file, a BigWig (bw) binary type file, or browser extensible data (bed) type file. A tab-delimited TSS table contains chromosome ID, genomic coordinates, strand information, and raw or normalized read counts of each sample. An example TSS table can be downlaoded from http://zlinlab.org/TSSr/ALL.samples.TSS.raw.txt.

The reference genome stored as a BSgenome data package must be provided for identification of TSSs from bam files. For available BSgenome genome objects, please check (https://bioconductor.org/packages/BSgenome). In case the BSgenome object of a species is unavailable in the BSgenome package (not in the list returned by BSgenome::available.genomes() function), a custom BSgenome object has to be built and installed before running this package. (See the vignette “How to forge a BSgenome data package” within the BSgenome package for instructions). A genome annotation file (GTF or GFF file) is required for assigning TSS clusters to genes or other genomic features by using the annotateCluster function of TSSr. The chromosome names in BSgenome, genome annotation file as well as the reference genome used for mapping sequencing reads must be consistent.

TSSr uses S4 object to store all input files and arguments and generate downstream analysis results.

• Launch TSSr

    library(TSSr)
  

• Example TSSr object. Users may explore the structure of a TSSr object and TSSr functions using provided example TSSr object "exampleTSSr". Skip this step if the user plan to generate a new TSSr object from bam or TSS table files.

    data(exampleTSSr)
    myTSSr <- exampleTSSr 
  

• Generating a new TSSr object from BAM or TSS table files. To generate a new TSSr object, several arguments of the constructer "new" function must be specified. The first step is to provide the path and file names of input files for argument "inputFiles". Four example BAM files (S01.sorted.bam, S02.sorted.bam, S03.sorted.bam, S04.sorted.bam) can be downloaded from http://www.zlinlab.org/TSSr.html. These BAM files were generated using nAnT-iCAGE from a budding yeast Saccharomyces cerevisiae (Lu and Lin 2019). The nAnT-iCAGE reads were mapped to the reference genome of S. cerevisiae (R64-2-1) using HISAT2 (Kim, Langmead et al. 2015). To reduced sample file size, each BAM file only includes sequencing reads mapped to two chromosomes (Chr I and Chr II). A TSS table genereated from the four BAM files can be downloaded from http://zlinlab.org/TSSr/ALL.samples.TSS.raw.txt, which can also be used as a input file. Assuming the bam or TSS files are saved to the current working directory:

    inputFiles <- c("S01.sorted.bam", "S02.sorted.bam", "S03.sorted.bam", "S04.sorted.bam") # if use a TSS table: inputFiles <- "ALL.samples.TSS.raw.txt"
    inputFilesType <- "bam" # set “inputFilesType” as “bamPairedEnd” for paired-end BAM files, and as "TSStable" if the input file is a TSS table 
  

• Install and launch required BSgenome object. The corresponding BSgenome object of S. cerevisiae “BSgenome.Scerevisiae.UCSC.sacCer3” is available from BSgenome Bioconductor package (https://bioconductor.org/packages/BSgenome).

    if (!requireNamespace("BSgenome.Scerevisiae.UCSC.sacCer3", quietly = TRUE))
    BiocManager::install("BSgenome.Scerevisiae.UCSC.sacCer3")
    library(BSgenome.Scerevisiae.UCSC.sacCer3) 
  

• Provide path and file name of genome annotation file. The example genome annotation of S. cerevisiae (R64-2-1) was originally obtained from the Saccharomyces Genome Database. It ban be downlaoded from http://zlinlab.org/TSSr/saccharomyces_cerevisiae_R64-2-1.gff. Assuming the gff file is saved to the current working directory:

    refSource <- "saccharomyces_cerevisiae_R64-2-1.gff" 
  

• Creating a new TSSr object using the constructor function. In addition to inputFiles, inputFilesType, genomeName, and the optional annotation file refSource, users provide sampleLabels, sampleLabelsMerged, and mergeIndex. In this example, SL01 and SL02 are biological replicates from the control condition, while SL03 and SL04 are replicates from the treatment condition. Therefore, sampleLabelsMerged = c("control", "treat") and mergeIndex = c(1, 1, 2, 2) merge the four samples into two groups.

  (organismName is retained as optional descriptive metadata for backward compatibility. It is not required by annotateCluster() and does not affect cluster-to-gene assignment, which is determined by the genomic coordinates and strands in refSource or refTable.)

    myTSSr <- new("TSSr", genomeName = "BSgenome.Scerevisiae.UCSC.sacCer3"
              ,inputFiles = inputFiles
              ,inputFilesType = inputFilesType
              ,sampleLabels = c("SL01","SL02","SL03","SL04")
              ,sampleLabelsMerged = c("control","treat")
              ,mergeIndex = c(1,1,2,2)
              ,refSource = refSource)
  

  Typing the object name displays a compact summary. Analysis results are initially marked as not run:

     myTSSr
      
    # TSSr object
    #   Genome: BSgenome.Scerevisiae.UCSC.sacCer3
    #   Samples: 4 (SL01, SL02, SL03, SL04)
    #   Merged samples: 2 (control, treat)
    #   TSSs: raw 0; processed 0
    #   Analyses:
    #     Tag clusters: <not run>
    #     Consensus clusters: <not run>
    #     Cluster shapes: <not run>
    #     Assigned clusters: <not run>
    #     Unassigned clusters: <not run>
    #     Filtered clusters: <not run>
    #     Enhancers: <not run>
    #     DE comparisons: <not run>
    #     TAG tables: <not run>
    #     Promoter shifts: <not run>
  

• TSS calling from bam files or retrieving TSS data from TSS table using "getTSS" function. The "getTSS" function identifies genomic coordinates of all TSSs and calculates read counts supporting each TSS from bam file. The results are available through TSSmatrix(myTSSr, data = "raw"). Before TSS calling, TSSr removes reads that are below certain sequencing quality and mapping quality. The default threshold for Phred quality score is 10, and mapping quality (MAPQ score) is 20. Users may change these arguments by setting different values for “sequencingQualityThreshold” and “mappingQualityThreshold” when running the “getTSS" function. If a mapped TSS sequencing read starts with a G that is a mismatch to the reference genome, the uncoded 5’ end G is likely the m7G cap, and thus the uncoded G will be removed from TSS calling by "getTSS". If a matched G at the 5’end of a tag is considered as an added cap, TSSr treats the 5’end of reads with matched G as genome-coded G, and the first G is not removed when calling TSS. This strategy is based on a stronge preference of PyPu dinucleotide at the [-1, +1] sites. This strategy also makes TSSr suitable for calling TSSs from 5’end sequencing reads that are not based on cap capture techniques.

Update in TSSr 0.99.16 (2026-06-16): The BAM implementation keeps the original G-specific biological rule for uncoded cap bases. For BAM input with softclippingAllowed = FALSE, TSSr removes only consecutive terminal uncoded G bases that mismatch the reference genome. On the plus strand, this means leading read G bases at the 5' end. On the minus strand, this means trailing read C bases that correspond to transcript-sense 5' G bases. Matched terminal G bases are treated as genome-coded and are not removed, and non-G mismatches are not trimmed. Consecutive mismatched terminal G bases are removed until the first matched G or non-G read base.

Considering that soft-clipping during mapping of sequencing reads to reference genome could also exclude non-matching 5’end bases other than unmatched Gs, which introduces false positive TSSs. Therefore, we’d recommend not to use soft-clipping for sequencing data. In the case that soft-clipped bam files are used, the users should set "softclippingAllowed = TRUE" for getTSS function.

Update in TSSr 0.99.16 (2026-06-16): When softclippingAllowed = TRUE, TSSr uses the aligner's aligned 5' boundary directly and does not apply uncoded G correction. When softclippingAllowed = FALSE, TSSr applies the corrected G-only uncoded G correction described above. BAM CIGAR widths are now calculated in reference coordinates: soft clips, hard clips, and insertions do not extend the reference interval, whereas reference-consuming operations determine the aligned interval and the minus-strand end coordinate.

    myTSSr <- getTSS(myTSSr)

Note: If you ran STAR with the default alignment parameters (without using --alignEndsType Extend5pOfRead1), please ensure that you set softclippingAllowed = TRUE when running TSSr.

TSS calling from bam files or retrieving TSS data from TSS table

    TSSmatrix(myTSSr, data = "raw")

    #          chr    pos strand SL01 SL02 SL03 SL04
    #     1:  chrI   1561      +    0    0    0    1
    #     2:  chrI   5759      +    1    0    0    0
    #     3:  chrI   5765      +    1    0    0    0
    #     4:  chrI   5773      +    1    0    0    0
    #     5:  chrI   5925      +    0    1    0    0
    #    ---
    # 29452:  chrI 227276      -    0    2    2    4
    # 29453:  chrI 227279      -    0    0    0    1
    # 29454:  chrI 227289      -    0    0    0    1
    # 29455:  chrI 227300      -    0    0    0    1
    # 29456:  chrI 227310      -    0    0    0    2

• TSSrawMatrix contains genomic coordinates and read counts of each called TSS. exportTSStable(..., data = "raw", merged = FALSE) writes these data to "ALL.samples.TSS.raw.txt" for reuse in TSSr or other tools. Setting merged = TRUE instead combines raw counts according to mergeIndex and writes "ALL.samples.TSS.raw.merged.txt".

    exportTSStable(myTSSr, data = "raw", merged = FALSE)
  

• The "plotCorrelation" function is used to calculate the pairwise correlation coefficients and plot pairwise scatter plots of TSS counts. A subset of samples can also be specified to display the pairwise correlations. Three correlation methods are supported: “pearson”, “kendall”, or “spearman”. The default setting generates scatter plots for all samples (samples = "all"). For plotting a subset of samples, users may provide a list of sample labels for "sample", e.g., samples = c("SL01","SL02","SL03").

    plotCorrelation(myTSSr, samples = "all")
  

• TSS includes "plotPCA" function to plot results principle component analysis (PCA) of all samples to obtain an overview of the TSS data in a low-dimensional subspace and for quality assessment. plotPCA will make a biplot which visualizes both how samples relate to each other in terms of the first principal component (PC1) and the second principal component (PC2) and simultaneously show how each variable contributes to each principal component.

   plotTssPCA(myTSSr, TSS.threshold=10)
  

• Merging samples (biological replicates). Users can merge multiple samples (e.g., biological replicates) into previously defined groups with mergeSamples function. The "mergeIndex" argument directs which samples will be merged and how the final dataset will be ordered accordingly. Retrieve the merged read counts and genomic coordinates with TSSmatrix(myTSSr, data = "processed").

    myTSSr <- mergeSamples(myTSSr)
    
    TSSmatrix(myTSSr, data = "processed")
    
    #          chr    pos strand control treat
    #     1:  chrI   1561      +       0     1
    #     2:  chrI   5759      +       1     0
    #     3:  chrI   5765      +       1     0
    #     4:  chrI   5773      +       1     0
    #     5:  chrI   5925      +       1     0
    #    ---                                  
    # 29452:  chrI 227276      -       2     6
    # 29453:  chrI 227279      -       0     1
    # 29454:  chrI 227289      -       0     1
    # 29455:  chrI 227300      -       0     1
    # 29456:  chrI 227310      -       0     2
  

  To return library sizes (number of total read counts) of merged samples in TSSr object in the specified order (note: order is specified in mergeSample function):

    librarySizes(myTSSr)
    # control   treat 
    #  596280  982703
  

  Library sizes are usually different among samples. To provide between-sample comparability, the raw read counts of each TSS need to be scaled as tags per million mapped reads (TPM) with normalizeTSS function. TSSr also provides two options to exclude TSSs with a low support: "poisson" and "tpm". If you intend to remove TSSs with low support using the “poisson” option, you should skip the "normalizeTSS" because the "poisson" method (method = "poisson") use the raw read count data to calculate the probability of observing k numbers of reads supporting each TSS based on the sequencing depth of the sample per the Poisson distribution. Only TSSs with a significantly larger number of supporting reads than expected (default threshold p < 0.01) are considered as qualified TSSs. Non-significant TSSs are thus filtered by TSSr. TSSr will normalize the raw read counts as TPM after filtering by the "poisson" method.

    poissonTSSr <- filterTSS(myTSSr, method = "poisson")
  

  If you use the "normalizeTSS" function, you may use the “filterTSS” option to remove TSSs with low support from mapped reads based on TPM value (method = "TPM"), and any TSS that has a lower TPM value than the user-defined threshold "tpmLow" will be removed (default TPM threshold = 0.1). For a threshold to remove a nonzero count, it must be greater than the smallest possible nonzero TPM for that sample, approximately 1e6 / librarySize. The reduced chromosome I example therefore uses tpmLow = 2 so that filtering has a measurable effect.

    myTSSr <- normalizeTSS(myTSSr)
    myTSSr <- filterTSS(myTSSr, method = "TPM", tpmLow = 2)
    TSSmatrix(myTSSr, data = "processed")
  
    
    #          chr    pos strand   control     treat
    #     1:  chrI   6163      +  0.000000  2.035203
    #     2:  chrI   6530      +  5.031193  0.000000
    #     3:  chrI   6534      +  0.000000  2.035203
    #     4:  chrI   6545      +  0.000000  2.035203
    #     5:  chrI   6548      + 20.124774  4.070406
    #    ---
    # 14891:  chrI 227270      -  3.354129 11.193616
    # 14892:  chrI 227271      -  0.000000  2.035203
    # 14893:  chrI 227274      -  5.031193  5.088007
    # 14894:  chrI 227276      -  3.354129  6.105609
    # 14895:  chrI 227310      -  0.000000  2.035203
  

  Similar to the raw read counts, the normalized and filtered (processed) TSS matrix can be exported to the delimited text file "ALL.samples.TSS.processed.txt" with "exportTSStable" or to bedGraph/BigWig files with "exportTSStoBedgraph". bedGraph/BigWig files can be visualized in the UCSC Genome Browser or Integrative Genomics Viewer (IGV) or Genome Browser at YeasTSS.org for selected yeast species.

exportTSStable(myTSSr, data = "processed") 
exportTSStoBedgraph(myTSSr, data = "processed", format = "bedGraph") 
exportTSStoBedgraph(myTSSr, data = "processed", format = "BigWig")

• Clustering TSSs to infer core promoters with “clusterTSS” function.

  The “clusterTSS” function was designed to group neighboring TSSs into distinct TSS clusters (TCs), representing putative core promoters. It implements a TSS clustering algorithm based on peaking identification, namely “peakclu” (peak clustering) (Lu and Lin, 2021). Briefly, peakclu applies a sliding window approach (default window size = 100 bp with step size = 1) to scan TSS signals from the 5′ end of both strands of each chromosome. In each window, the TSS with the highest TPM value was identified as the peak. The surrounding TSSs are grouped with the peak into a TC. The clustering process of a TC terminates if a TSS is ≥ n bp (default n = 25) away from the nearest upstream TSS. In addition to setting a minimal allowed distance between peaks, TSSr offers another option to set maximal allowed extension distance between neighboring TSSs around peaks, which enables users to define the boundaries between neighboring core promoters. clusterTSS calculates inter-quantile width of a core promoter based on the cumulative distribution of TSS signals within the promoter. The positions of the 10th to 90th quantiles of TSS signals, which include at least 80% transcription initiation signals within a cluster, were defined as the 5’ and 3’ boundaries of the core promoter. The clustering step might be slow especially when the number of TSSs is in millions. Using multicores is highly recommended.

Note: The variable localThreshold = 0.02 means that only signals with a minimum TSS intensity greater than 2% of the dominant TSS will be included in the cluster; otherwise, they will be filtered out. Meanwhile, clusterThreshold = 1 means that a cluster will only be retained if the total TPM (Transcripts Per Million) value exceeds 1.

        myTSSr <- clusterTSS(myTSSr, method = "peakclu",peakDistance=100,extensionDistance=30
	           ,localThreshold = 0.02, clusterThreshold = 1
	           ,useMultiCore=FALSE, numCores=NULL)

        head(tagClusters(myTSSr, sample = "control"))
        head(tagClusters(myTSSr, sample = "treat"))

• The results of TSS clustering can be exported to delimited text file with "exportClustersTable" function or bedGraph files with "exportClustersToBed" function.

exportClustersTable(myTSSr, data = "tagClusters")
exportClustersToBed(myTSSr, data = "tagClusters") 

• Aggregating consensus TSS clusters

  TSSr infers a set of consensus core promoters using the “consensusCluster” function to assign the same ID for TCs belong to the same core promoter, which allows subsequent comparative studies across samples. TCs from different samples are considered to belong to the same consensus core promoter if the distance of their dominant TSSs is smaller than a user-defined distance (default = 50 bp). Similarly to clusterTSS function, consensusCluster function also returns genomic coordinates, sum of TSS tags, dominate TSS coordinate, a lower (q0.1) and an upper (q0.9) quantile coordinates, and interquantile widths for each consensus cluster in each sample.

    myTSSr <- consensusCluster(myTSSr, dis = 50, useMultiCore = FALSE, numCores = NULL)
  

  Similarly, the detailed information of consensus clusters can be exported to delimited text files with "exportClustersTable" function or bedGraph files with "exportClustersToBed" function.

exportClustersTable(myTSSr, data = "consensusClusters")
exportClustersToBed(myTSSr, data = "consensusClusters")	

• Quantification of core promoter shape

  Core promoter shape reflects the distribution of TSS signals within a core promoter. TSSr provides three different options, inter-quantile width, promoter shape score (PSS) and shape index (SI), to quantify core promoter shape. Inter-quantile width refers to the distance between the locations of the 10th percentile to the 90th percentile TSS signals within a TSS cluster. Thus, it measures the width of a core promoter, but lacks the information of distribution patterns of TSS signals within a core promoter. Inter-quantile width could be significantly affected by different clustering methods. plotInterQuantile function plots interquantile width of each sample.

    plotInterQuantile(myTSSr,samples = "all",tagsThreshold = 1)
  

Promoter shape score (PSS) integrates both inter quantile width and the observed probabilities of tags at every TSSs within a cluster (Lu and Lin 2019). PSS can be calculated using using shapeCluster function with method set as “PSS”. The smaller value represents the sharper core promoter. PSS is 0 representing singletons. SI is determined by the probabilities of tags at every TSSs within one cluster (Hoskins, Landolin et al. 2011). SI is also calculated using shapeCluster function with method set as “SI”. The greater value represents the sharper core promoter. The SI = 2 represents singletons.

    myTSSr <- shapeCluster(myTSSr,clusters = "consensusClusters", method = "PSS",useMultiCore= FALSE, numCores = NULL)

    head(clusterShape(myTSSr, sample = "control"))
    head(clusterShape(myTSSr, sample = "treat"))

The distributions of core promoter shape (PSS or SI) can be illustrated by a histograms with plotShape function. Only the most recently calculated shape method is retained and can be retrieved with clusterShape().

    plotShape(myTSSr)

The complete list of PSS can be exported to delimited text file with "exportShapeTable" function.

exportShapeTable(myTSSr)

• Annotation (Assigning TCs to genes)

  Assigning TCs to downstream genes as their core promoters is required for annotation of the 5’ boundaries of genomic features. This process is also a prerequisite for further interrogations of regulated transcription initiation at the gene level. TSSr offers the “annotateCluster” function to assign TCs to their downstream genes. The assignment of a TC to a gene is based on the distance between the position of the dominant TSS of a TC and the annotated 5’ends of coding sequences (CDS) or transcripts. The default maximum distance between the dominant TSS and CDS is 1000 bp (“upstream = 1000”). If a TC overlaps with the CDS of an upstream gene, the dominant TSS of the TC must be within 500 bp to the 3’end of the overlapping CDS by default (“upstreamOverlap = 500”) to be eligible to assign to its downstream gene. If the 5’ends of annotated transcripts, instead of CDS, are used for TC assignment, users should set “annotationType” as “transcript,” and the default distance parameter is 500 bp. Because the genomes size and the number of introns vary substantially among organisms, it is necessary to apply customized criteria for TC assignment for different organisms. Users are advised to adjust the assignment criteria for core promoter assignment in TSSr. By default, only TCs with ≥ 0.02 TPM are used for the annotation process. To reduce transcriptional or technical noise of small clusters downstream a strong cluster, the filterCluster argument was set as "filterClusterThreshold = 0.02", indicating that any TC with TPM value < 0.02 TPM will be excluded from assigning to genes.

    myTSSr <- annotateCluster(myTSSr,clusters = "consensusClusters",filterCluster = TRUE,
              filterClusterThreshold = 0.02, annotationType = "genes"
              ,upstream=1000, upstreamOverlap = 500, downstream = 0)
              
    head(assignedClusters(myTSSr, sample = "control"))
    head(assignedClusters(myTSSr, sample = "treat"))
  

  Instead of visualizing TSSs and core promoters in the UCSC Genome Browser or IGV, plotTSS function is able to generate publish ready figures when list of interested genes are provided and plotting region is specified.

  The results of TC annotation can be exported to delimited text files with "exportClustersTable" function.

exportClustersTable(myTSSr, data = "assigned")
exportClustersTable(myTSSr, data = "unassigned")

• Analysis of enhancers

  TSS data allow for the robust identification of enhancers by transcription of enhancer RNAs (eRNAs). Active enhancers produce bidirectional transcription of capped eRNAs, resulting in two diverging tag clusters by at most 400 bp. TSSr can identify this bidirectional cluster pairs and calculate a sample-set wide directionality score D for each locus (Andersson et al., 2014). D = (F-R)/(F+R), where F is the aggregated normalized tag counts in forward strandm and R is the aggregated normalized tag counts in reverse strand. Putative enhancers were then filtered with |D| < 0.8.

    myTSSr <- callEnhancer(myTSSr, flanking = 400)
  

  The results of putative enhancers can be exported to delimited text files with "exportEnhancerTable" function.

    exportEnhancerTable(myTSSr)
  

• Differential expression analysis

  The number of tags at each TSS reflects the number of transcripts initiated at the TSS. Thus, TSS data can be used for expression profiling. With specified sample pairs for comparison, deGene function counts raw tags of each consensus clusters and utilizes the DESeq2 package (Love, Huber et al. 2014) for differential expression analysis.

    myTSSr <- deGene(myTSSr,comparePairs=list(c("control","treat")), pval = 0.01,useMultiCore=FALSE, numCores=NULL)
  

  Differential expression analysis results can be visualized by plotDE function which generates a volcano plots. Names of genes differential expressed between the compared pairs are displayed on the dots when the withGeneName argument is set as TRUE.

    plotDE(myTSSr, withGeneName = "TRUE",xlim=c(-2.5, 2.5),ylim=c(0,10))
  

  Differential expression analysis results can also be exported to a text file with exportDETable function.

    exportDETable(myTSSr, data = "sig")
  

• Core promoter shifts

  One gene might have multiple core promoters which can be used differently in different samples. TSSr implements degree of shift (Ds) algorithm (Lu and Lin 2019) to quantify the degree of promoter shift across different samples.

    myTSSr <- shiftPromoter(myTSSr,comparePairs=list(c("control","treat")), pval = 0.01)
            
    head(PromoterShift(myTSSr, comparison = "control_VS_treat"))
  

  Below is an example of core promoter shift in genes YBR168W and YBL067C. The two major core promoters are differently used in control and treat samples.

    plotTSS(myTSSr,samples=c("control","treat"),tssData = "processed"
    ,clusters = "assigned"
    ,clusterThreshold = 0.02 
    ,genelist=c("YBR168W","YBL067C")
    ,up.dis =500,down.dis = 100,yFixed=TRUE)
  

Results of core promoter shift analysis can also be exported to a text file with exportShiftTable function.

    exportShiftTable(myTSSr)

Zhenguo Lin, PhD

Department of Biology, Saint Louis University, USA.

Email: [email protected]

Andersson, R., Gebhard, C., Miguel-Escalada, I., Hoof, I., Bornholdt, J., Boyd, M., Chen, Y., Zhao, X., Schmidl, C., Suzuki, T., Ntini, E., Arner, E., Valen, E., Li, K., Schwarzfischer, L., Glatz, D., Raithel, J., Lilje, B., Rapin, N., Bagger, F. O., … Sandelin, A. (2014). "An atlas of active enhancers across human cell types and tissues". Nature 507(7493), 455–461.

Arribere, J. A. and W. V. Gilbert (2013). "Roles for transcript leaders in translation and mRNA decay revealed by transcript leader sequencing." Genome Res 23(6): 977-987.

Cumbie, J. S., M. G. Ivanchenko and M. Megraw (2015). "NanoCAGE-XL and CapFilter: an approach to genome wide identification of high confidence transcription start sites." BMC Genomics 16: 597.

Cvetesic, N., H. G. Leitch, M. Borkowska, F. Muller, P. Carninci, P. Hajkova and B. Lenhard (2018). "SLIC-CAGE: high-resolution transcription start site mapping using nanogram-levels of total RNA." Genome Res 28(12): 1943-1956.

Kruesi, W. S., L. J. Core, C. T. Waters, J. T. Lis and B. J. Meyer (2013). "Condensin controls recruitment of RNA polymerase II to achieve nematode X-chromosome dosage compensation." Elife 2: e00808.

Mahat, D. B., H. Kwak, G. T. Booth, I. H. Jonkers, C. G. Danko, R. K. Patel, C. T. Waters, K. Munson, L. J. Core and J. T. Lis (2016). "Base-pair-resolution genome-wide mapping of active RNA polymerases using precision nuclear run-on (PRO-seq)." Nat Protoc 11(8): 1455-1476.

Mahat, D. B., H. Kwak, G. T. Booth, I. H. Jonkers, C. G. Danko, R. K. Patel, C. T. Waters, K. Munson, L. J. Core and J. T. Lis (2016). "Base-pair-resolution genome-wide mapping of active RNA polymerases using precision nuclear run-on (PRO-seq)." Nat Protoc 11(8): 1455-1476.

Malabat, C., F. Feuerbach, L. Ma, C. Saveanu and A. Jacquier (2015). "Quality control of transcription start site selection by nonsense-mediated-mRNA decay." Elife 4.

Murata, M., H. Nishiyori-Sueki, M. Kojima-Ishiyama, P. Carninci, Y. Hayashizaki and M. Itoh (2014). "Detecting expressed genes using CAGE." Methods Mol Biol 1164: 67-85.

Pelechano, V., W. Wei and L. M. Steinmetz (2013). "Extensive transcriptional heterogeneity revealed by isoform profiling." Nature 497(7447): 127-131.

Takahashi, H., T. Lassmann, M. Murata and P. Carninci (2012). "5' end-centered expression profiling using cap-analysis gene expression and next-generation sequencing." Nat Protoc 7(3): 542-561.