TaRTLEt: Transcriptionally-active Riboswitch Tracer Leveraging Edge deTection

Predicted transcript pools for transcriptionally active riboswitches

Riboswitches are often associated with transcription start sites (TSS), allowing transcription start site analysis to serve as a strategy for identifying novel riboswitches and regulatory RNA (Rosinski-Chupin, Soutourina & Martin-Verstraete, 2014; Rosinski-Chupin et al., 2019; Adams et al., 2017; Yu, Vogel & Förstner, 2018). We reasoned that a transcription-terminating riboswitch in the 5′ UTR of a gene should produce different patterns of RNA-seq read coverage in the ON vs. OFF state (Figs. 1A, 1B). In either case, transcription should begin at the TSS and proceed through the riboswitch itself. In the OFF state, formation of the transcription-terminator structure should truncate transcripts near the riboswitch 3′ end (Fig. 1A top). By contrast, in the ON state, RNA polymerase should read through, producing transcripts that extend well into the downstream gene (Fig. 1A bottom). A riboswitch that produces these two patterns under different conditions exhibits behavior consistent with riboswitch-mediated conditional regulation of transcription termination.

Fragment distributions and terminal pileups

What relation does the coverage map of sequencing reads bear to the pool of transcripts? Each RNA-seq read represents some region of a transcript, but because transcripts are typically fragmented prior to sequencing we know only that the 5′ end of a read must align with a fragment terminus and not whether that fragment is internal or terminal to a transcript. Fragmentation often, as in the validation dataset used here, results in a random biased distribution (Fig. S1) of fragment sizes. Importantly, however, the new 3′ and 5′ termini produced by each fragmentation event come in neighboring pairs. Even though the total number of fragment termini, and thus the total number of read termini, may far exceed the number of original transcript termini, we can still extract a transcript size signal by considering pileups of read termini. Because paired-end sequencing generates reads that represent both ends of each sequenced fragment, the resulting data offers a higher signal:noise ratio than single-end sequencing: there, fragmentation is more extensive and at most one terminus of each fragment is represented in the read set, so that comparable sequencing depth yields greater stochastic fluctuations in the set of termini captured. We therefore designed TaRTLEt v1.0 to restrict analysis to paired-end datasets. Future versions could be extended to analyze appropriately deep single-end datasets.

We reasoned that a convolution approach would allow us to both extract the needed transcript size signal and smooth the overall signal. Convolution is heavily used in applications like edge-detection algorithms (Canny, 1986; Albawi, Mohammed & Al-Zawi, 2017), where it can act as an attenuating filter for neighboring signals. Here, we want neighboring pileups of 3′ and 5′ read termini to drop out, filtering out transcript-internal fragmentation to reveal the transcript-terminal signal. Tracking 3′ ends as +1 and 5′ ends as −1, we can first compute the signed sum of read termini for each position in the reference (Fig. 1C), then calculate a weighted average of surrounding sums for each position. We call the resulting convolution’s features “peaks”, and we seek to identify peaks at which a pileup of 3′ termini is not accompanied by a pileup of 5′ termini immediately downstream.

The choice of kernel function for computation of the weighted average affects the filter’s selectivity for changes in coverage. Rather than a flat uniform distribution, we chose to use kernels constructed by discretizing a Gaussian distribution, such that tuning the spread of this distribution (by the standard deviation σ) tunes the filter’s selectivity for how sharply coverage must change. Because sites of transcription termination, especially intrinsic termination, are constrained (Roberts, 2019; Gusarov & Nudler, 1999), riboswitch-mediated transcription termination should produce a sharp drop in coverage. To match this expectation, we might use a tight convolution (small σ) kernel that attenuates peaks arising from changes over >1–2 nt. However, both biological and technical factors (e.g., thermodynamic fluctuations; processes in library preparation, sequencing, and quality control) are expected to broaden the biological signal to some degree. Through trial and error, we chose a default σ of 1.5 for kernel generation, corresponding to a filter that passes coverage changes localized within a 6-nt region.

Regardless of the number of fragmentation events, the convolution always returns non-zero peaks at positions with asymmetric pileups of 5′ or 3′ termini. These surviving peaks should represent a range of events: transcription termination should be relatively rare and should cause a large, tightly centered coverage drop, giving rise to high, narrow peaks; sporadic processes (e.g., artifacts in library prep or sequencing) should be relatively common and yield lower, broader peaks. While these differences in peak shapes should provide an opportunity to further refine the set of candidate peaks, the probability and details of the processes generating background peaks may be sequence-dependent in as yet poorly understood ways, complicating comparison to a global baseline. Instead, we reasoned that we could collect a per-locus “noise set” of peaks within a defined neighborhood scaled to riboswitch size. Fitting a two-dimensional Gaussian noise distribution to this collection of peaks provides a basis for comparison, allowing us to robustly identify unusually sharp peaks at each locus as candidate sites of riboswitch-mediated transcription termination.

Probing fragment-end distributions across treatments

For any of these sites to reveal conditional regulation by a transcription-terminating riboswitch, the efficiency of termination (that is, the proportion of transcripts that terminate) near the riboswitch 3′ end must vary across treatments that evoke ON and OFF states. For a riboswitch that acts by transcription termination, termination efficiency should be low (low or no coverage drop) in the ON state and high (large coverage drop) in the OFF state. We bound the region of interest for this search to a range of 0.5 ×riboswitch size both upstream and downstream of the riboswitch 3′ end and examine the fractional coverage change (that is, (change in coverage across peak)/(average coverage depth across entire riboswitch)) at each convolution peak in the region under each treatment condition. To compare results at coincident peaks across treatment conditions, we account for the possibility that peak location might shift slightly due to small variations in process noise by clustering peak locations across treatment conditions (Fig. 2). Any cluster that captures both ON and OFF states of a transcription-terminating riboswitch should exhibit fractional coverage change with a mean that is significantly more negative and a variance that is significantly larger than baseline (Figs. 2 and 3).

Identification of riboswitches and open reading frames

Input sequences in FASTA format are processed through Infernal v1.1.5 (http://eddylab.org/infernal) (Nawrocki & Eddy, 2013b) and Prodigal v2.6.3 (https://github.com/hyattpd/Prodigal) (Hyatt et al., 2010). For Infernal, each input sequence is passed to cmscan using the - - cut_ga - - rfam - - nohmmonly - - noali - - tblout options against a comprehensive riboswitch covariance model of families compiled from the latest Rfam collection (version 14.9 as of writing) (https://rfam.org/) (Griffiths-Jones et al., 2003). cmscan output is parsed to identify riboswitch loci for further processing. Prodigal is called with the default output and translation options. For each identified riboswitch, the downstream gene is identified from the Prodigal output and recorded for downstream plot annotations. The downstream gene is defined as the first open reading frame past the riboswitch 5′ end and on the same strand. This definition captures open reading frames (ORFs) that overlap with part of the identified riboswitch sequence, in order to maximize the likelihood of capturing ORFs that are biologically relevant to the riboswitch.

Riboswitch reference generation

To reduce computational requirements for transcript mapping, TaRTLEt trims input genomic sequences back to riboswitch neighborhoods, producing reference sequences that stretch from 1,000 nt upstream of each identified riboswitch locus to 1,000 nt downstream. (For metagenomic inputs, where riboswitch loci may lie less than 1,000 nt from contig ends, neighborhoods are permitted to be smaller or asymmetric.) We chose this neighborhood size on the basis of common practice in paired-end sequencing library preparation, where the target fragmentation size is often less than 1,000 nt to preserve base quality (Tan et al., 2019). Consequently, the generated reference neighborhoods are typically (96.61% of riboswitch-aligned reads in our validation dataset) large enough to capture both a riboswitch-aligned read and its paired-end mate. All references generated by TaRTLEt are oriented 5′→ 3′; references derived from (-) strand riboswitches are reverse complemented before saving to file.

Reference alignment

RNA-seq reads are aligned to the generated riboswitch reference using HISAT2 (http://daehwankimlab.github.io/hisat2/) (Kim et al., 2019). HISAT2 was chosen for its fast per-read alignment time, good parallel scaling, and acceptable alignment rates (Musich, Cadle-Davidson & Osier, 2021). HISAT2 calls are made using the default options, but any additional valid alignment options may be passed through the TaRTLEt interface. The current TaRTLEt implementation requires that the input RNA-seq datasets be paired-end. The resulting SAM files are converted to sorted BAMs for downstream processing.

We recommend running HISAT2 using the options - - no-unal - - score-min L,0,−0.4. The use of - - no-unal dramatically reduces output file sizes by omitting unaligned reads. The use of - - score-min L,0,−0.4 decreases stringency to allow reads with alignment mismatches/soft-clipping to be included in the output; any alignment errors are handled during downstream processing. The validation dataset presented here was run using these parameters in addition to -p 40 (performance tuning with multiple threads) and -t (wall-time logging).

Processing mapped reads

TaRTLEt heavily utilizes the pysam (https://github.com/pysam-developers/pysam) library (a python interface for samtools (Li et al., 2009) and htslib (Bonfield et al., 2021)) to parse BAM files and generate alignment data objects for downstream data transformation. Once mapped reads/read pairs have been collected from the sorted BAMs, they are filtered to discard pairs in invalid orientations (RF, TANDEM, FFR, RRF), which may reflect sequencing errors. Only pairs in the EQ and FR orientations, and pairs with only one mapped read, are retained as valid. The EQ (“equal”) orientation denotes reads aligned in opposite directions that overlap across their entire length, as is often the case when the fragment sequenced is ∼150 bp, comparable in scale to individual short reads. In the FR orientation, the left-most aligned base (looking 5′→ 3′ in the reference) is either shared between the F and R oriented reads or belongs only to the F oriented read, and the right-most aligned base belongs exclusively to the R read.

Next, TaRTLEt generates fragment coverage pileups using the filtered set of reads for every position in the reference sequence. Single reads (with unmapped mates) are discarded during counting unless the - - allow-single-reads flag is passed to TaRTLEt. TaRTLEt tracks five types of coverage. Two can incorporate data from single reads if included: “read coverage” counts bases directly mapped by a read, while “clipped fragment coverage” counts bases mapped by HISAT2’s soft-clipped read regions. The remaining three types require both mates to be mapped. “Overlapped fragment coverage” counts bases overlapped by both mates; “inferred fragment coverage” counts unmapped bases that lie between the mapped regions of the two mates; and “fragment termini coverage” counts only the mapped base that is the origin of read synthesis, i.e., looking 5′→ 3′ in the reference, the left-most mapped base for an F read and right-most mapped base for an R read. This last category captures the ends of sequenced fragments to feed into the signed sum of termini for convolution analysis; we count the fragment terminus captured by the F read as −1 and the fragment terminus captured by the R read as +1. Reads that have soft-clipped regions are discarded from this step unless the - - allow-soft-clips flag is specified. Fragment termini counts are removed from the read coverage and/or overlapped coverage counts, and read coverage counts are removed from the overlapped coverage counts, such that no mapped base is counted more than once. This ensures that counts across all coverage types at a given position sum to the read depth at that position.

Finally, to ensure sufficient coverage for further analysis, riboswitch loci are filtered based on the average read coverage across the length of the riboswitch. The coverage threshold can be adjusted with the - - min-cov-depth option; the validation presented here uses the default value of 15.

Transformation and feature identification

Fragment termini coverage pileups are transformed using a convolution to extract information about transcription starts and stops. TaRTLEt generates a 51-element one-dimensional weighted kernel by discretizing a Gaussian probability density function with σ = 1.5. The weighted kernel is then applied to the fragment termini coverage array through the convolution to generate an array of peaks corresponding to transcription events. The convolution transforms coverage at each position into convolution amplitude (arbitrary units).

TaRTLEt conducts simple peak detection as follows for a positive peak: a local convolution amplitude maximum is defined as the peak “summit”. The peak region covers all the bases on either side of the summit at which convolution amplitude (a) is monotonically decreasing and (b) does not enter or cross a zero-region, defined as the interval $\left(-0.1,0.1\right)$. To delineate ambiguous shoulder regions between peaks, peak bounds follow pythonic inclusive-exclusive convention where the left peak bound represents the left-most valid base according to the monotonicity and zero-crossing criteria and the right peak bound represents the first base that fails the criteria. Negative- and positive-amplitude peaks are considered candidate transcription start and termination events, respectively.

Next, TaRTLEt assesses candidate transcription-termination peaks individually to determine which rise above background. Working with per-locus, per-condition convolution plots, TaRTLEt compares each positive peak within a user-defined region to a two-dimensional multivariate normal (MVN) distribution fitted to the widths and amplitudes of each other peak in the region (the “noise set”). The region that feeds into the noise set is defined relative to riboswitch length, with parameters defining the proportion of the riboswitch that should be excluded at the 5′ end and how far past the riboswitch 3′ end the region should stretch. (Our validation analysis used - - ext-prop −0.3 1.0 to define a noise set region that excludes the first 30% of the riboswitch, includes the remainder, and continues beyond the riboswitch 3′ end to include a region equal in length to the complete riboswitch.) The likelihood of drawing a given peak from the MVN is used as a pseudo p-value, with a significance threshold of <0.05. Peaks that pass this test and exceed a user-configurable coverage change threshold (by default, a drop of ≥20% of the mean coverage across the locus) are called as transcription termination events.

To facilitate biological interpretation of riboswitch activity, TaRTLEt then asks, for each condition and each locus, whether any passing transcription-termination peaks were observed. This within-condition check does not in itself offer evidence for riboswitch regulation by transcription termination. Instead, it (1) feeds into the cluster-level calls described below; (2) supports manual curation of results; and (3) allows the user to ask, among the transcription-terminating riboswitches identified below, what set of experimental conditions evoke the ON state vs. the OFF state.

Peak clustering and statistics

To detect differential transcription termination, TaRTLEt compares data from all treatment conditions at each riboswitch locus. All peaks that fall within the riboswitch region of interest (user-configurable; by default, within ± 0.5 ×  riboswitch size of the riboswitch 3′ end) under any condition are clustered by position to enable cross-condition comparisons. The riboswitch 3′ end defines position 0, and peak positions are given relative to this site as a fraction of the riboswitch size. For example, a position of −0.15 would indicate that the peak is located 15% of the riboswitch size upstream of the riboswitch 3′ end. Clustering is implemented via fclusterdata from the SciPy package and follows a cophenetic distancing algorithm with complete linkage. The default (user-configurable) cophenetic distance threshold of 0.04 was chosen through trial and error (Fig. S3).

Once peaks are assigned to clusters, cross-cluster comparisons at each locus are used to test for condition-dependent transcription termination. Each peak is characterized by the change in total fragment coverage across its bounds. The fractional coverage change is defined as the difference in total coverage between the left and right bounds of the peak, expressed as a fraction of the average total read coverage between the 5′ and 3′ ends of the riboswitch. TaRTLEt performs a Mann–Whitney one-tailed U-test (scipy.stats.mannwhitneyu) to compare each cluster’s fractional coverage change to the mean fractional coverage change of all other peaks, to test whether the cluster mean is more negative than the mean of the set of all other peaks.

Testing for differences in cluster variance must account for bias from differences in cluster sizes. TaRTLEt uses Levene’s test (scipy.stats.levene) to run pairwise comparisons between the variance of a given cluster enclosing n peaks and that of random samples of n peaks chosen without replacement from all other peaks at the same riboswitch locus until the superset of peaks is exhausted or at least 60 comparisons have been made. If the superset is exhausted before 60 comparisons are made, all peaks are replaced and the clusters are randomly resampled to reach 60 comparisons. TaRTLEt returns the median of the resulting 60 p-values for each peak cluster.

Finally, TaRTLEt combines peak- and cluster-level tests to assess significance. To be called as showing evidence of condition-dependent transcription termination, a riboswitch locus must have at least one cluster that (i) includes at least one significant transcription-termination peak (i.e., distinct from the noise-set MVN at p <0.05 and coverage drop greater than or equal to the user-defined threshold; see previous section); (ii) passes the Mann–Whitney U-test for mean fractional coverage change; and (iii) passes the Levene test for increased variance.

Description of outputs

TaRTLEt outputs data on two layers of abstraction: single-condition convolution peaks, and cross-condition peak clusters. At each level, TaRTLEt generates both tabular output on each feature and plots of each riboswitch locus. For each candidate peak in the convolution, the single-condition tabular output (peak_log.csv) captures the riboswitch locus, the RNA-seq dataset, the location and shape of the peak, the “noise set” of comparator peaks, TaRTLEt’s pass/fail call for that peak as a candidate transcription termination event, and the reason for that call. Single-condition locus plots are sorted into pass/ and fail/ subdirectories, where pass/ contains locus plots for all conditions in which any peak, regardless of location within the locus, was called as a transcription termination event. A sample locus plot of the molybdenum cofactor (MOCO) riboswitch from E. coli is shown in Fig. 4; for a sample peak log, see Table S1.

At the level of cross-condition cluster comparisons, TaRTLEt records each cluster’s mean and variance test statistics in cluster_stats.csv, then uses this file in conjunction with peak_log.csv to generate a summary peak plot that allows easy inspection of condition-dependent modulation of transcription termination efficiency near the 3′ ends of all riboswitches in the dataset. Figure 5 shows a sample summary peak plot. The user may simply scan the summary peak plot to identify loci in which a cluster was marked significant (*), optionally using the locus plots to identify particular experimental conditions that evoked a change in riboswitch transcription termination efficiency (at loci identified as modulating transcription termination) or to spot-check peak call quality.