TSSr: an R/Bioconductor package for comprehensive analyses of transcription start site (TSS) data

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

How to cite TSSr?

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.

Important migration note for TSSr 0.99.17 and later

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.

1. Introduction

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.

TSSr flowchart

2. Pre-requisites:

R version

R packages

3. Installing TSSr Package:

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)

4. Input data for TSSr

Input data files

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.

Reference genome

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.

S4 object

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

5. Usage:

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

01_TSS_correlation_plot_of_all_samples

02_PCA_plot

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

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"))
exportClustersTable(myTSSr, data = "tagClusters")
exportClustersToBed(myTSSr, data = "tagClusters") 
exportClustersTable(myTSSr, data = "consensusClusters")
exportClustersToBed(myTSSr, data = "consensusClusters")	

03_Interquantile_plot_of_ALL_samples

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)

04_Shape_plot_of_ALL_samples

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

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

05_Volcano_plot

06_TSS_graphs

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

    exportShiftTable(myTSSr)

7. Contact:

Zhenguo Lin, PhD

Department of Biology, Saint Louis University, USA.

Email: zhenguo.lin@slu.edu

8. Reference:

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.