--- title: "TSSr Vignette" author: "Zhaolian Lu, Keenan Berry, Zhenbin Hu, Yu Zhan, Tae-Hyuk Ahn, Zhenguo Lin" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteEngine{knitr::rmarkdown} %\VignetteIndexEntry{TSSr Vignette} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` # Introduction TSSr is a comprehensive R/Bioconductor package for analyzing transcription start site (TSS) sequencing data. It supports multiple input formats including BAM, BED, BigWig, and TSS tables, and provides a complete workflow from TSS identification through downstream analyses such as core promoter shape analysis, gene annotation, differential expression, and promoter shifting. ## Installation Install TSSr through Bioconductor with `BiocManager`: ```r if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("TSSr") ``` And load `TSSr`: ```{r library, results='hide', message=FALSE} library(TSSr) ``` # Citation If you use TSSr, please cite the following article: ```{r citation, eval=TRUE} citation("TSSr") ``` # Getting help For general questions about the usage of TSSr, use the [official Bioconductor support forum](https://support.bioconductor.org) and tag your question "TSSr". We strive to answer questions as quickly as possible. For technical questions, bug reports and suggestions for new features, we refer to the [TSSr github page](https://github.com/Linlab-slu/TSSr/issues). # Quick start This vignette runs an analysis from input files rather than beginning with a precomputed `TSSr` object. The same object, named `myTSSr`, is returned and reassigned at every stage of the workflow. ## Bundled input files The example TSS table contains four samples. The accompanying GFF file provides gene coordinates for annotation. `system.file()` locates these files inside the installed package, and `mustWork = TRUE` ensures that a missing file is reported immediately. ```{r bundled-input-files} exampleInput <- system.file( "extdata", "example-tss-table.tsv", package = "TSSr", mustWork = TRUE ) exampleAnnotation <- system.file( "extdata", "example-annotation.gff3", package = "TSSr", mustWork = TRUE ) ``` We create one `TSSr` object named `myTSSr` and retain this name throughout the vignette. The four input samples are assigned to two merged groups: `SL01` and `SL02` form the control group, while `SL03` and `SL04` form the treat group. ```{r construct-object} myTSSr <- TSSr( genomeName = "BSgenome.Scerevisiae.UCSC.sacCer3", inputFiles = exampleInput, inputFilesType = "TSStable", sampleLabels = c("SL01", "SL02", "SL03", "SL04"), sampleLabelsMerged = c("control", "treat"), mergeIndex = c(1, 1, 2, 2), refSource = exampleAnnotation ) myTSSr ``` ## Using your own data To analyse your own data, replace the values in the `TSSr()` call above while keeping the object name `myTSSr`. The supplied values must satisfy the requirements below. ### Constructor arguments | Argument | Required | Accepted values and restrictions | |:--|:--:|:--| | `genomeName` | Yes | A single character string naming an installed `BSgenome` package that matches the genome assembly used to generate the input files, for example `"BSgenome.Scerevisiae.UCSC.sacCer3"`. | | `inputFiles` | Yes | A character vector containing paths to existing regular files. For `"bam"`, `"bamPairedEnd"`, `"bed"`, `"tss"`, and `"BigWig"`, supply one file per original sample. For `"TSStable"`, supply exactly one table containing all samples. | | `inputFilesType` | Yes | One case-sensitive value chosen from `"bam"`, `"bamPairedEnd"`, `"bed"`, `"tss"`, `"BigWig"`, or `"TSStable"`. See the table below. | | `sampleLabels` | Yes | A non-empty character vector naming the original samples. Except for `"TSStable"`, its length must equal `length(inputFiles)`. For `"TSStable"`, the labels identify sample columns in the single input table. | | `sampleLabelsMerged` | No | A character vector naming the groups produced by replicate merging, in the order indicated by `mergeIndex`. It must be supplied together with `mergeIndex`. If both arguments are omitted, each sample is treated as a separate group. | | `mergeIndex` | No | A numeric vector assigning every original sample to a merged group. Its length must equal `length(sampleLabels)`. Use consecutive positive integers beginning with `1`, such as `c(1, 1, 2, 2)`. The number of unique indices must equal `length(sampleLabelsMerged)`. | | `refSource` | No | Zero or one path to an existing genome-annotation file used by `annotateCluster()`. Supply it when cluster annotation will be performed and no annotation table has otherwise been added to the object. | | `organismName` | No | Zero or one character string describing the organism. This optional metadata is not required for TSS clustering. | ### Accepted input types The value of `inputFilesType` is case-sensitive and must be one of the following: | `inputFilesType` | Input represented | Required file arrangement | |:--|:--|:--| | `"bam"` | Single-end BAM alignments | One BAM file per sample | | `"bamPairedEnd"` | Paired-end BAM alignments | One BAM file per sample | | `"bed"` | BED-formatted TSS data | One BED file per sample | | `"tss"` | CTSS/TSS position-and-signal data | One file per sample | | `"BigWig"` | BigWig signal data | One BigWig file per sample | | `"TSStable"` | A tabular TSS matrix containing multiple samples | Exactly one table containing all sample columns | For example, four original samples can be merged into two experimental groups using `sampleLabelsMerged = c("control", "treat")` and `mergeIndex = c(1, 1, 2, 2)`. Index `1` corresponds to `"control"`, and index `2` corresponds to `"treat"`; the ordering of `sampleLabelsMerged` must therefore match the numeric group indices. If replicate merging is not required, omit both `sampleLabelsMerged` and `mergeIndex`. `TSSr()` will automatically retain each original sample as a separate group. Do not supply only one of these two arguments. ## Importing TSS data `getTSS()` reads the input file and stores the imported signal in the raw TSS matrix. Displaying its dimensions provides a quick confirmation that the data were imported before later analysis steps are run. ```{r import-tss, message=FALSE} myTSSr <- getTSS(myTSSr) dim(TSSmatrix(myTSSr, data = "raw")) ``` For formats other than `"TSStable"`, the import step may require additional format-specific settings. See `?getTSS` before adapting the example to BAM or other position-based files. # TSS data processing The `myTSSr` object now contains the chromosome I signals imported above. The next chunk merges the original samples into their experimental groups, normalizes the merged signals, and applies the TPM filter in sequence. ```{r TSSprocessing, tidy=FALSE} # Merge replicates myTSSr <- mergeSamples(myTSSr) # Normalization myTSSr <- normalizeTSS(myTSSr) # TSS filtering rowsBeforeFiltering <- nrow(TSSmatrix(myTSSr, data = "processed")) myTSSr <- filterTSS(myTSSr, method = "TPM", tpmLow = 2) c( before = rowsBeforeFiltering, after = nrow(TSSmatrix(myTSSr, data = "processed")) ) ``` These functions populate the processed TSS matrix in the returned object. The reported row counts show how many genomic TSS positions remain after the TPM filter. Both values are 482 for the bundled fixture because no positions are removed at the selected threshold; user datasets may show a decrease at this step. Use `TSSmatrix()` to inspect this matrix through the public accessor: ```{r TSSprocessing-show} # Access the processed TSS matrix without exposing internal storage head(TSSmatrix(myTSSr, data = "processed")) ``` The first three columns identify the chromosome, genomic position, and strand; the remaining columns contain processed signals for the merged sample groups. This matrix is the input for TSS clustering. # TSS clustering Nearby TSSs can represent transcription initiation from the same core promoter. `clusterTSS()` groups these positions into sample-specific tag clusters using the peak-based `peakclu` method. Here, candidate peaks must be at least 100 bp apart, a cluster can extend across gaps of at most 30 bp, TSSs below 2% of the local peak are removed, and clusters with signal below 1 TPM are discarded. `consensusCluster()` then relates tag clusters across sample groups whose dominant TSSs are within 50 bp, providing common regions for comparative analyses. ```{r TSSclustering, tidy=FALSE} # TSS clustering myTSSr <- clusterTSS(myTSSr, method = "peakclu", peakDistance = 100, extensionDistance = 30, localThreshold = 0.02, clusterThreshold = 1, useMultiCore = FALSE, numCores = NULL ) # Aggregating consensus clusters myTSSr <- consensusCluster(myTSSr, dis = 50, useMultiCore = FALSE) ``` The returned object now contains both `tagClusters` and `consensusClusters`. Their accessors expose the sample-specific calls and the cross-sample regions, respectively: ```{r TSSclustering-show} # Tag clusters per sample head(tagClusters(myTSSr, sample = "control")) # Consensus clusters across samples head(consensusClusters(myTSSr, sample = "control")) ``` Each table reports cluster coordinates, the dominant TSS, TSS signal, and inter-quantile boundaries. Consensus clusters are used below so that promoter shape, annotation, and between-group comparisons refer to the same genomic regions. See `?clusterTSS` and `?consensusCluster` for parameter details. # Core promoter shape Core promoter shape describes how initiation signals are distributed within a cluster. `shapeCluster()` can calculate either the Shape Index (SI) or Promoter Shape Score (PSS); this example calculates PSS for the consensus clusters. PSS combines the inter-quantile width with the distribution of TSS signals. Smaller PSS values indicate sharper promoters, and a singleton has a score of zero. ```{r shapeCluster, tidy=FALSE} # Calculating core promoter shape score myTSSr <- shapeCluster(myTSSr, clusters = "consensusClusters", method = "PSS", useMultiCore = FALSE, numCores = NULL ) ``` The scores are stored in `clusterShape` for each merged sample group and can be retrieved without accessing the object's internal slots: ```{r shapeCluster-show} # Shape scores per cluster head(clusterShape(myTSSr, sample = "control")) ``` These values can be used to compare promoter architectures across samples or exported with `exportShapeTable()`. Only the most recently calculated shape method is retained; see `?shapeCluster` for the SI alternative. # Annotation of core promoters Annotation connects consensus clusters to downstream genes so that subsequent analyses can operate at gene level. The `refSource` supplied when `myTSSr` was constructed points to the bundled GFF annotation. With the settings below, a cluster is considered up to 1000 bp upstream of a gene start; overlapping upstream features use the more restrictive 500-bp rule. Weak downstream clusters below 2% of the strongest cluster in the promoter region are filtered to reduce recapping and other low-signal events. ```{r annotateCluster, tidy=FALSE} # Assign clusters to the annotated features myTSSr <- annotateCluster(myTSSr, clusters = "consensusClusters", filterCluster = TRUE, filterClusterThreshold = 0.02, annotationType = "genes", upstream = 1000, upstreamOverlap = 500, downstream = 0) ``` The returned object separates clusters into assigned, unassigned, and, when filtering is enabled, filtered tables. The following accessor displays clusters that were assigned to genes in the control group: ```{r annotateCluster-show} # Clusters assigned to genes head(assignedClusters(myTSSr, sample = "control")) ``` The `gene` column records the assigned feature, while the cluster columns retain the genomic position and signal information. Assigned clusters support the gene-level differential expression and promoter-shift analyses below; unassigned clusters are considered during enhancer identification. See `?annotateCluster` when adapting promoter windows to another organism. # Enhancer identification Active enhancers can produce nearby, divergent transcription signals. `callEnhancer()` searches the unassigned clusters for opposite-strand pairs no more than 400 bp apart, applies a directionality filter, and excludes candidates within 2000 bp of a gene's main annotated promoter. The result should therefore be interpreted as a set of putative enhancers rather than independently validated regulatory elements. ```{r callEnhancer, tidy=FALSE} myTSSr <- callEnhancer(myTSSr, flanking = 400, dis2gene = 2000) ``` Candidate pairs and their strand-specific signals are stored in the `enhancers` slot. Use the public accessor to inspect one sample group: ```{r callEnhancer-show} # Putative enhancers per sample head(enhancers(myTSSr, sample = "control")) ``` The output identifies the paired clusters, their dominant TSS positions and signals, and the directionality score `D`. See `?callEnhancer` for the candidate selection parameters and `exportEnhancerTable()` for saving the complete table. # Differential expression analysis After annotation, TSSr can aggregate raw TSS counts from assigned clusters by gene and use DESeq2 to compare sample groups. `comparePairs` sets the order of the comparison, here `control` versus `treat`, and `pval = 0.01` retains genes with an adjusted p value below 0.01 in the significant-results table. DESeq2 is a suggested dependency and must be installed to run this step. The bundled example contains only a few annotated genes, so it uses `fitType = "mean"` to avoid an unstable local dispersion-trend fit. The default remains `fitType = "parametric"` for ordinary analyses. ```{r deGene, tidy=FALSE} # Gene-level differential expression using DESeq2 myTSSr <- deGene(myTSSr, comparePairs = list(c("control", "treat")), pval = 0.01, useMultiCore = FALSE, numCores = NULL, fitType = "mean") ``` The returned object stores the gene-level count tables and both complete and significant DESeq2 results. `DEtables()` selects a comparison and result set: ```{r deGene-show} # All differential expression results head(DEtables( myTSSr, comparison = "control_VS_treat", result = "all" )) ``` Each row in the displayed table corresponds to one annotated gene and includes the DESeq2 effect-size and significance statistics. `result = "significant"` returns the subset with an adjusted p value below the `pval` threshold. At `pval = 0.01`, none of the six genes in this small bundled dataset meet that criterion, so its significant-results table is empty. Use `plotDE()` to visualize the results or `exportDETable()` to save them. # Core promoter shifts Genes with multiple core promoters may use them at different proportions between conditions. `shiftPromoter()` tests the annotated promoter profiles for the requested `control`-versus-`treat` comparison. Sparse gene-level 2-by-2 tables can have small expected counts; in that case, the function emits one summary warning for all affected tests and their approximate p values may be unreliable. See `?shiftPromoter` for details. ```{r shiftPromoter, tidy=FALSE} # Calculate core promoter shifts myTSSr <- shiftPromoter(myTSSr, comparePairs = list(c("control", "treat")), pval = 0.01) ``` The test results are stored by comparison in the `PromoterShift` slot. The accessor below retrieves the table for this pair: ```{r shiftPromoter-show} # Promoter shift results head(PromoterShift(myTSSr, comparison = "control_VS_treat")) ``` The table reports gene-level promoter-shift statistics and associated p values. Because small expected counts weaken the chi-squared approximation, inspect the underlying promoter signals before interpreting affected genes. # Exporting results TSSr can save processed TSS values and cluster tables for downstream analysis, or write genome-browser tracks for visualization. The functions below write to the current working directory. The vignette temporarily switches to `tempdir()` so that rendering does not add result files to the package source. ```{r export} oldwd <- setwd(tempdir()) # Export processed TSS values to ALL.samples.TSS.processed.txt exportTSStable(myTSSr, data = "processed") # Export cluster tables exportClustersTable(myTSSr, data = "assigned") # Export to BED format for genome browsers exportClustersToBed(myTSSr, data = "consensusClusters") # Export TSS to bedGraph exportTSStoBedgraph(myTSSr, data = "processed") setwd(oldwd) ``` This creates a tab-delimited processed TSS table, per-sample assigned-cluster tables, BED files for assigned consensus clusters, and strand-specific bedGraph tracks for the processed TSS signal. BED and bedGraph files can be loaded into genome browsers such as UCSC Genome Browser or IGV. In an interactive analysis, set the working directory to the desired output location before calling these functions; see the individual export help pages for file-naming options. # Precomputed example object The package also includes `exampleTSSr`, a precomputed object for users who want to inspect a populated `TSSr` object without rerunning the workflow. It is not used as the input to the analysis above and does not replace `myTSSr`. ```{r precomputed-example} data("exampleTSSr") exampleTSSr ``` # Complete workflow in one pipeline The step-by-step assignments above make intermediate results easy to inspect. Once the settings have been established, the same analysis can be expressed as one base R pipeline. The object is still named `myTSSr`; only the style of the assignment changes. Replace the constructor values with the paths, labels, genome, and annotation for your own experiment. ```{r complete-pipeline, eval=FALSE} myTSSr <- TSSr( genomeName = "BSgenome.Scerevisiae.UCSC.sacCer3", inputFiles = exampleInput, inputFilesType = "TSStable", sampleLabels = c("SL01", "SL02", "SL03", "SL04"), sampleLabelsMerged = c("control", "treat"), mergeIndex = c(1, 1, 2, 2), refSource = exampleAnnotation ) |> getTSS() |> mergeSamples() |> normalizeTSS() |> filterTSS(method = "TPM", tpmLow = 2) |> clusterTSS( method = "peakclu", peakDistance = 100, extensionDistance = 30, localThreshold = 0.02, clusterThreshold = 1, useMultiCore = FALSE, numCores = NULL ) |> consensusCluster( dis = 50, useMultiCore = FALSE ) |> shapeCluster( clusters = "consensusClusters", method = "PSS", useMultiCore = FALSE, numCores = NULL ) |> annotateCluster( clusters = "consensusClusters", filterCluster = TRUE, filterClusterThreshold = 0.02, annotationType = "genes", upstream = 1000, upstreamOverlap = 500, downstream = 0 ) |> callEnhancer( flanking = 400, dis2gene = 2000 ) |> deGene( comparePairs = list(c("control", "treat")), pval = 0.01, useMultiCore = FALSE, numCores = NULL, fitType = "mean" ) |> shiftPromoter( comparePairs = list(c("control", "treat")), pval = 0.01 ) ``` The export functions are intentionally kept outside this pipeline because they write files rather than add a new analysis stage to the returned object. # Session info Recording package and R versions makes the rendered analysis easier to reproduce and diagnose: ```{r sessionInfo, tidy=FALSE} sessionInfo() ```