Improved functions for nonlinear sequence comparison using SEEKR

New SEEKR functions including those to estimate significance of k-mer similarity

We developed a series of functions within the SEEKR package that enable end-users to estimate the significance of SEEKR-derived Pearson's r-values. We also created new functions to visualize different aspects of k-mer similarity, and deprecated other functions that were either redundant with new ones or difficult to install using current versions of Python (Fig. 1A). Updates to SEEKR are available for download through the Python Package Index, GitHub, and the Docker Hub (https://pypi.org/project/seekr/; https://github.com/CalabreseLab/seekr; https://hub.docker.com/r/calabreselab/seekr). A separate GitHub page contains the code and console commands used to generate the figures in this manuscript (https://github.com/CalabreseLab/seekr2.0_update_manuscript).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

FIGURE 1.

New SEEKR functions including those to estimate significance of k-mer similarity. (A) Console/command-line functions in SEEKR. Blue text, new or updated functions. Black text, retained functions. Orange text, deprecated functions. “seekr_filter_gencode” replaces “seekr_canonical_gencode.” “seekr_kmer_leiden” replaces “seekr_graph.” The “fasta-shuffle-letters” utility in the MEME Suite (Bailey et al. 2015) can be used in place of “seekr_gen_rand_rnas.” (B) Documentation of the Pearson's r-value threshold at a P-value of P = 0.05 for each probability distribution that best fits the list of r-values derived from 100,000 randomly selected pairwise comparisons from the set of GENCODE v43 Ensembl_canonical lncRNAs at k = 6. “Com.10 best fit”: pairwise comparisons fit to the common10 distributions in scipy.stats. “All best fit”: pairwise comparisons fit to all 124 probability distributions in scipy.stats. The Kolmogorov–Smirnov test was used to assess goodness-of-fit to each distribution. Fitting was repeated 10× for each of the common10 and all distributions. johnsonSU, Johnson's SU distribution. genhyperbolic, generalized hyperbolic distribution. (C) Number of lncRNAs detected as significant at the P < 0.05 threshold when fitting the set of background sequences to the common10 and at the specified subsetting sizes. (D) Number and identity of lncRNAs detected as significantly similar to XIST using k = 4, 5, or 6 when fit to a constant distribution (log-normal). (E) The SEEKR console commands used to download GENCODE lncRNAs and calculate adjusted P-values for a sequence comparison of interest. See also https://github.com/CalabreseLab/seekr2.0_update_manuscript.

As a framework to illustrate the new SEEKR functions, we describe their use in the context of identifying similarities between the lncRNA XIST and other human lncRNAs. We first outline the process of estimating the significance of SEEKR-derived Pearson's r-values before highlighting other functions in SEEKR. To estimate the significance of Pearson's r-values, SEEKR relies on a user-defined set of background sequences to determine what level of similarity would be expected between any two lncRNAs selected by chance. This set of background sequences should be large enough that the total number of pairwise comparisons between sequences in the set will generate a well-populated distribution. For studies in human or mouse, lncRNA sequences curated by GENCODE are a suitable resource (Frankish et al. 2023), and SEEKR provides a function that enables download from current and past versions of the GENCODE database. For the set of background sequences used below, we used SEEKR to download and deduplicate the list of “Ensembl_canonical” GENCODE lncRNAs that are greater than 500 nt in length (henceforth referred to as GENCODE canonical lncRNAs).

After selecting a set of background sequences, the next step in estimating significance using SEEKR is to identify what probability distribution best fits the list of SEEKR-derived Pearson's r-values generated from all possible pairwise comparisons of background sequences. Parameters describing the shape of this probability distribution are then input into the cumulative distribution function to calculate the P-value describing the likelihood that any pair of lncRNAs is more similar to each other than would be expected by chance. The function “seekr_find_dist” is used to find the probability distribution that best fits the set of background sequences, and the function “seekr_find_pval” is used to calculate the P-values from this probability distribution. Lastly, if desired, the function “seekr_adj_pval” can be used to correct P-values for multiple-testing.

For distribution fitting, “seekr_find_dist” defaults to fit the “common10” distributions defined in Python's Fitter library, and “seekr_find_pval” defaults to selecting the best-fitting distribution from the distributions examined by the user in “seekr_find_dist.” If desired, instead of fitting to the common10, end-users may fit a custom subset or all 124 of the distributions available in scipy.stats (Virtanen et al. 2020); they may also elect to select a distribution other than the best-fitting one for P-value calculations. For background sequences whose comparisons populate a normal-like distribution, such as the GENCODE canonical lncRNAs, we find that fitting to the common10 distributions is efficient and practical. However, if the accuracy of fitting at the threshold of significance is of the highest priority, using all available models in “seekr_find_dist” may be preferable to the common10. Using the set of GENCODE canonical lncRNAs at k = 6, we observed the best-fitting probability distribution from the common10 was most frequently log-normal (Fig. 1B). The best-fitting distribution from the full set of 124 in scipy.stats was most frequently Johnson's SU distribution, but still yielded Pearson's r-values at the P = 0.05 threshold that were on average only ∼3% different from the P-value threshold calculated by the log-normal (Fig. 1B). For all P-value calculations in the work below, we fit our background lncRNAs to the log-normal distribution.

In the case of large background sets, such as the set of GENCODE canonical lncRNAs, which contains 15,500 sequences, fitting a background distribution to all possible pairwise comparisons is computationally intensive and unnecessary for obtaining relatively accurate parameters that describe a distribution's shape (all possible pairwise comparisons from the GENCODE canonical set will yield 15,500², or 240,250,000 data points in a distribution). In these instances, end-users may wish to fit distributions on a subset of pairwise comparisons that are randomly selected from the complete list. Tests on the set of GENCODE canonical lncRNAs indicate that subsetting on 100,000 values provides stable estimates of significance at low computational expense (Fig. 1C). We therefore performed subsetting on 100,000 values for all P-value calculations in our work below.

Lastly, in prior work with SEEKR, we observed that biological signals, including similarities in gene regulatory function and protein-binding, are often detected to similar extents at k-mer lengths of k = 4, 5, or 6 (Kirk et al. 2018; Sprague et al. 2019). We therefore sought to determine the extent to which significance estimates change as a function of k-mer length. We compared the lncRNA XIST to all GENCODE canonical lncRNAs at values of k ranging from k = 4–6. At a threshold of P < 0.05, these comparisons yielded a sum of 3040 distinct lncRNAs across the three values of k; 90% of the 3040 lncRNAs were detected as significantly similar to XIST using at least two different values of k, and 78% were detected as significant at all three values of k (Fig. 1D). Thus, when analyzing the significance of SEEKR-derived Pearson's r-values, different k-mer lengths detect similar sets of lncRNAs.

Putting these recommendations to use, selecting a k-mer length k = 6 and using the set of GENCODE canonical lncRNAs as both the background and the set of lncRNAs to search, we identified 2768 lncRNAs whose overall 6-mer contents are more similar to XIST than would be expected by chance; after correcting for multiple-testing using the Benjamini–Hochberg method, this list was reduced in size to 1265 lncRNAs (Fig. 1E; Table 1; Supplemental Table S1; https://github.com/CalabreseLab/seekr2.0_update_manuscript).

SEEKR functions to visualize features of k-mer similarity

To demonstrate graphing functions within SEEKR, we next investigated LINC00632, DLX6-AS1, and PCDH10-DT, which are highly similar to XIST (Table 1) and whose Ensembl_canonical transcript isoforms are strongly supported by short-read RNA-seq data (de Goede et al. 2021). XIST is thought to encode its repressive functions through the cumulative action of separate domains that interact with different subsets of proteins (Trotman et al. 2021). While full repression by XIST requires its entire sequence, some of its regions most critical for repression are comprised of tandemly repeated sequences, referred to as Repeats A, B, D, E, and F (Dixon-McDougall and Brown 2021, 2022). Therefore, to determine whether our select XIST-like lncRNAs harbored regional similarity to XIST, we separated the lncRNAs into ∼500 nt fragments and compared the fragments in each lncRNA to each fragment within XIST. Because data suggest Repeats A, B, D, E, and F function as discrete domains, we evaluated each XIST Repeat as a single intact fragment and separated the intervening XIST intervals into ∼500 fragments (Fig. 2A; Supplemental Files S1 and S2 contain fragments of XIST and all other GENCODE canonical lncRNAs used in this study, respectively). From these analyses, we found that LINC00632, DLX6-AS1, and PCDH10-DT each harbor significant similarity to Repeat E but lack similarity to the other Repeats in XIST (Fig. 2B–D). Moreover, these three lncRNAs also harbor significant similarity to fragments distributed across the final exon of XIST (Figs. 2B–D). A closer examination revealed that the similarities could be attributed to a uniform enrichment of k-mers rich in A and T nucleotides, whereas k-mers rich in G and C nucleotides were among the most variably enriched (Figs. 2E–G). We also observed that many of the 500 nt fragments within XIST’s final exon (i.e., Int.6 in Fig. 2A) were significantly more similar to each other than would be expected by chance (Fig 2H,I; Supplemental Fig. S1). We selected a k-mer length of k = 4 for these analyses, so as to reduce the number of k-mers that have zero count values in the analyzed fragments, following prior guidance described in Kirk et al. (2018) and Sprague et al. (2019). However, similar trends were detected at k-mer lengths k = 5 and 6 (Supplemental Fig. S2).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

FIGURE 2.

SEEKR functions to visualize features of k-mer similarity. (A) Intervals/chunks of XIST used to identify regional similarities. Repeats A–F in XIST are colored in non-gray shades and were each evaluated as single intervals, regardless of length. “Repeat F” is annotated here as a region that spans nt 743–1807 in XIST (ENST00000429829.6). (B–D) Regional similarities between XIST and LINC00632, DLX6-AS1, and PCDH10-DT. The latter three lncRNAs were separated into ∼500 nt consecutive intervals and compared to the XIST intervals in (A). Heatmaps show P-values for each comparison. Comparisons were made using k = 4. (E) Heatmap displaying log2-transformed z-scores at k = 4, for the 3′ portion of XIST (Interval 6.2 to its transcript end), Repeat E, and LINC00632, DLX6-AS1, and PCDH10-DT relative to the set of nonidentical “Ensembl_canonical” GENCODE v43 lncRNA transcripts ≥500 nt in length. Log2-transformed z-scores were derived using seekr_kmer_counts. Only those k-mers beginning with “T” nucleotides are displayed for clarity. (F) Top 10 4-mers that are the most enriched in the 3′ portion of XIST and the full-length lncRNAs from (E) relative to the set of “Ensembl_canonical” also from (E). (G) Top 10 4-mers that exhibit the most variable enrichment in the 3′ portion of XIST and the full-length lncRNAs from (E) (i.e., highest standard deviation across the RNAs). (H) Regional similarities within XIST. Heatmap shows P-values for each pairwise comparison of XIST intervals shown in (A). Comparisons were made using k = 4. (I) Dendrogram showing clusters of XIST intervals from (A) that harbor related k-mer contents. Clusters (different colors) were defined as the nodes whose correlation distances are <0.7*[max_distance_in_dendrogram]. Comparisons were made using k = 4. The SEEKR console/command line functions that were used to generate graphs are displayed above each figure panel.

Domain-based search identifies chromatin-associated lncRNAs that harbor XIST-like fragments

Given that the tandem repeats within XIST are some of the regions most essential for its ability to induce and maintain gene silencing, we next used SEEKR to perform a parallel search for XIST-like lncRNAs. In this search, we separated all 15,550 GENCODE canonical lncRNAs into ∼500 nt fragments (including XIST, as a positive control), and then used SEEKR to identify lncRNAs that contain fragments with significant k-mer similarity to XIST Repeats A, B, D, E, and F (at an unadjusted P-value of <0.05). We then summed the number of XIST-similar fragments in each lncRNA and used these sums to rank lncRNAs by their overall XIST-likeness. This analysis identified several intriguing lncRNAs (Table 2; Supplemental Table S2). As might have been expected, XIST ranked high (third). However, the known repressive lncRNA KCNQ1OT1 ranked eighth (in a tie with five other transcripts) (Schertzer et al. 2019; Quinodoz et al. 2021). Moreover, 18 of the top 22 lncRNAs were expressed at detectable levels in K562 cells, and of those 18 lncRNAs, 17 were chromatin-enriched, including the architectural lncRNA NEAT1 (Table 2; Dunham et al. 2012; Obuse and Hirose 2023). Thus, with this simple fragment-based search, we identified several chromatin-enriched lncRNAs that harbor domains that resemble those required for repression by XIST, including KCNQ1OT1 and NEAT1. Additional XIST-like transcripts of interest include ENST00000605862.6, ENST00000506640.3, and ENST00000622550.2, functionally uncharacterized lncRNAs which are both spliced and expressed in multiple cell types (Dunham et al. 2012; de Goede et al. 2021). We also identified RENO1, a conserved lncRNA whose depletion in mouse embryonic stem cells causes changes in gene expression and a failure to differentiate properly into neurons (Hezroni et al. 2020).

The tandem repeats in XIST appear in a specific order and are thought to serve as recruitment centers for specific RNA-binding proteins (RBPs) (Trotman et al. 2021). With the possible exception of KCNQ1OT1, XIST-like fragments within the lncRNAs of Table 2 appeared in an order that was different from that found in XIST (Fig. 3). However, we hypothesized that at least some of the fragments might still associate with the same RBPs as their cognate domains in XIST. Therefore, in three lncRNAs of interest, we examined the eCLIP read density profiles of six RBPs that have been subject to eCLIP in K562 cells and exhibit characteristic enrichment in the tandem repeats in XIST (Van Nostrand et al. 2020). These RBPs included RBM15 (enriched over Repeat A); HNRNPM (enriched over Repeat F region); HNRNPK (enriched over Repeats B and D); and MATR3, PTBP1, and TIA1 (enriched over Repeat E; Fig. 3A). In all three lncRNAs, we observed XIST-like fragments that colocalized with the eCLIP signal from the expected RBPs. In KCNQ1OT1, Repeat A-, F-, and B-like fragments colocalized with RBM15, HNRNPM, and HNRNPK, respectively (Fig. 3B). In NEAT1, Repeat F-, B/D-, and E-like fragments colocalized with HNRNPM, HNRNPK, and MATR3/PTBP1/TIA1, respectively (Fig. 3C). In ENST00000605862.6, A- and B-like fragments colocalized with RBM15 and HNRNPK, respectively (Fig. 3D). Thus, while the arrangement and number of XIST-like fragments in each of those lncRNAs differs from XIST, several of their XIST-like fragments exhibited enriched association for the expected RBPs.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

FIGURE 3.

Domain-based search identifies chromatin-associated lncRNAs that harbor XIST-like fragments. (A–D) Screen images from the UCSC Genome Browser displaying gene annotations, location and identity of XIST-like fragments, read density from chromatin-associated RNA-seq and from Bru-seq, and background-corrected eCLIP signal for a subset of RBPs enriched over the tandem repeats of XIST.

XIST-like fragments in K562-expressed lncRNAs associate with XIST-binding proteins

We next examined the extent to which SEEKR-defined XIST-like fragments colocalized with XIST-binding proteins in all lncRNAs expressed in K562 cells (n = 3068; expression defined as >0.0625 TPM from total RNA-seq data [Dunham et al. 2012]). We used the Wilcoxon signed-rank test to compare eCLIP read counts in each XIST-like fragment to eCLIP read counts in a paired set of fragments that were randomly shuffled among the 3068 K562-expressed lncRNAs (excluding the regions defined by SEEKR to be XIST-like). Strikingly, we found that the cognate binding proteins of each XIST Repeat were significantly enriched over XIST-like fragments when compared to randomly shuffled controls (Fig. 4). Moreover, for all Repeats except Repeat A, the expected cognate protein(s) was the one that was the most significantly enriched over the set of corresponding XIST-like fragments. Specifically, HNRNPK was the most significantly enriched protein over Repeat B- and D-like fragments, MATR3 and PTBP1 were the most significantly enriched over Repeat E-like fragments, and HNRNPM was the most significantly enriched over Repeat F-like fragments (Fig. 4). For Repeat A-like fragments, we observed that all five proteins were significantly enriched over shuffled controls (Fig. 4); these findings can be rationalized because Repeat A is arguably the most complex of all XIST repeats, and it contains consensus binding motifs for each of the proteins we analyzed. Thus, fragment-based analyses with SEEKR can be used to identify domains in other lncRNAs that associate with proteins that exhibit enriched binding to a domain of interest.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

FIGURE 4.

XIST-like fragments in K562-expressed lncRNAs associate with XIST-binding proteins. (A) Summary table displaying the number of K562 lncRNA fragments detected as significantly similar (P < 0.05) to each XIST Repeat (no. of Frags.), along with Bonferroni-corrected adjusted P-values from Wilcoxon signed-rank tests describing the comparison between eCLIP read density in the set of XIST-like fragments versus shuffled control fragments, for each of five XIST-binding proteins. (B) Box plots displaying eCLIP read density for select XIST-like proteins under XL, the set of detected XIST-like fragments, and C, the set of shuffled controls, for each of five XIST-binding proteins. Yellow highlights, the protein(s) that associate most robustly with each XIST Repeat (Van Nostrand et al. 2020).

Not all expressed, XIST-like lncRNA annotations are supported by RNA-seq data in K562 cells

In carrying out these analyses, we also noted that transcript annotations for certain highly ranked, expressed lncRNAs in Table 2 were not strongly supported by short-read RNA-seq data from K562 or HepG2 cells (Dunham et al. 2012). The most highly ranked XIST-similar lncRNA, HELLPAR, is annotated as an unspliced, 200 kb lncRNA that begins at the 3′ end of the protein-coding gene PARPBP and terminates near the 5′ end of the protein-coding gene IGF1 (van Dijk et al. 2012; Frankish et al. 2021). Using standard approaches to quantify transcript abundance from RNA-seq data, HELLPAR registers as being expressed and chromatin-associated in K562 cells (Table 2; Supplemental Table S2). However, RNA-seq read density suggests that transcription across the HELLPAR locus is not due to the expression of a single lncRNA. Rather, it appears to result from a combination of imprecise termination of the upstream PARPBP gene, leading to read-through transcription in the 5′ half of the HELLPAR locus, and the transcription of a separate lncRNA, LINC02456, whose promoter is found in the center of the locus and whose transcribed product runs antisense through the IGF1 protein-coding gene (Fig. 5A; Dunham et al. 2012; de Goede et al. 2021; Zhu et al. 2023).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

FIGURE 5.

Not all expressed, XIST-like lncRNA annotations are supported by RNA-seq data in K562 cells. (A–F) Screen images from the UCSC Genome Browser displaying gene annotations, location and identity of XIST-like fragments, and read density from ribosome-depletion total RNA-seq (total), chromatin-associated RNA-seq (chromatin), and Bru-seq experiments (bruseq) performed in K562 cells (Dunham et al. 2012; Luo et al. 2020).

Likewise, again visualizing short-read RNA-seq data from Dunham et al. (2012), we found examples of other XIST-like lncRNAs that register as being expressed and chromatin-associated in K562 cells but whose transcript annotations were not strongly supported. These included TSIX, whose only RNA-seq read density coincided with the gene structure of XIST on the opposite strand, and may have arisen due to imperfect strand-specificity of the dUTP RNA-seq protocol (Fig. 5B; Levin et al. 2010); LINC02604, which is annotated as a monoexonic lncRNA but appears to be part of a much larger transcribed and spliced region (Fig. 5C); and monoexonic lncRNAs that overlap in the sense direction relative to longer protein-coding or lncRNA genes: ENST00000623075.1, PRNCR1, and ENST00000624919.1 (Fig. 5D–F). In each of these examples, RNA-seq abundance estimators register the lncRNAs as being expressed, underscoring the importance of examining RNA-seq read density relative to lncRNA transcript annotations before investigating SEEKR-derived sequence similarities.

Using the SEEKR P-value function to augment evolutionary studies of lncRNAs

We previously demonstrated that SEEKR can be applied to detect k-mer similarity in lncRNAs between different species (Kirk et al. 2018; Sprague et al. 2019). We reasoned that the newly developed P-value function in SEEKR could help to prioritize results from between-species lncRNA comparisons, just as above we demonstrated that it can help prioritize results from within-species lncRNA comparisons. Therefore, we investigated the SEEKR similarity scores of pairs of human and mouse lncRNAs that a recent study identified as potentially homologous owing to their conserved genomic locations and shared patterns of protein-binding motifs (Huang et al. 2024). Specifically, we examined two sets of lncRNAs: a set of 5564 human lncRNAs that have predicted homologs in mouse (the “HM” set); and a set of 570 human lncRNAs that have predicted homologs in mouse and zebrafish (the “HMZ” set) (Huang et al. 2024). Comparing human and mouse lncRNA pairs against the distribution of pairwise comparisons between all human lncRNAs at k-mer length k = 6, we observed that the distributions of SEEKR similarity scores in both the HM and HMZ sets exhibited significant rightward shifts, as would be expected for conserved sequences (P < 2×10⁻¹⁶ for both comparisons; Wilcoxon rank sum test [Fig. 6; Kirk et al. 2018]). About 20% of lncRNA pairs from the HM and HMZ sets harbored SEEKR-derived 6-mer contents that were more similar to each other than would be expected by chance, using the set of all human lncRNAs as background (P < 0.05; 1222 of 7355 pairs in the HM set and 165 of 570 lncRNA pairs in the HMZ set; Fig. 6; orange bars). At the level of 6-mer content, most of the significantly similar pairs were less similar to each other than lncRNAs known to exhibit strong conservation, such as H19, NEAT1, XIST, MALAT1, and KCNQ1OT1 (Fig. 6; orange bars vs. green lines). Examining further the set of 165 HMZ lncRNA pairs that harbored significant k-mer similarity to each other, we identified two pairs of human/mouse lncRNAs that each harbored a high number of fragments similar to XIST’s tandem repeat regions: HRAT92 (human) and 6330403L08Rik (mouse); and SNHG14 (human) and 9330162G02Rik (mouse). We defined a “high” number of fragments in this analysis as any lncRNA that had greater than seven XIST-like fragments, which corresponds to a fragment count that ranks in the 99.5% or greater when considering all GENCODE canonical human lncRNAs (Supplemental Tables S2 and S3). These analyses demonstrate ways in which SEEKR can be used to prioritize similarity between lncRNAs across species to augment other approaches to study evolutionary relationships.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

FIGURE 6.

Using the SEEKR P-value function to augment evolutionary studies of lncRNAs. (A,B) Distributions of SEEKR-derived r-values calculated using k = 6 for all human–mouse lncRNA pairs defined as coPARSE in human and mouse (A; HM set) and human, mouse, and zebrafish (B; HMZ) set (Huang et al. 2024). HM and HMZ r-values are binned in discrete bars, with corresponding counts shown on the right-hand y-axes. Orange bars, binned counts of lncRNA pairs whose P-values of similarity are <0.05; blue bars, all other bins. Green lines mark SEEKR-derived r-values between exemplar pairs of conserved human and mouse lncRNAs. SEEKR-derived r-values for all pairwise comparisons of human lncRNAs are shown as continuous gray distributions whose densities are displayed on the left-hand y-axes.