Long-range RNA structures in the human transcriptome beyond evolutionarily conserved regions

RIC-seq and RNA-seq experiments

The results of RIC-seq experiments conducted in seven human cell lines including GM12878, H1, HeLa, HepG2, IMR90, K562, and hNPC (two bioreplicates each) were downloaded from the Gene Expression Omnibus under the accession numbers GSE127188 and GSE190214 in FASTQ format. The matched control RNA-seq experiments were downloaded from the ENCODE consortium under the accession numbers listed in Table S1. RIC-seq and the matched RNA-seq data were processed by RNAcontacts pipeline (Margasyuk et al., 2023d) using February 2009 (hg19) assembly of the human genome and GENCODE transcript annotation v41lift37, which were downloaded from Genome Reference Consortium (Church et al., 2011) and GENCODE website (Harrow et al., 2012), respectively.

Pairs of nested contact clusters

PHRIC pipeline starts with the output files of the RNAcontacts, which contain the coordinates of ligation junctions extracted from RIC-seq read alignments. These input files are combined into one table with an additional field listing the experiment identifier. Next, the intrachromosomal junctions below 200 nts and junctions between different chromosomes are removed, and the remaining junctions are aggregated into contacts by computing the number of supporting reads in each experiment. To cluster closely located contacts, we first cluster their left and right split points using bedtools cluster program (Quinlan, 2014), which merges split points into one cluster if the distance between them is not greater than 10 nts, as in Margasyuk et al. (2023d). The contacts are combined into a cluster if their left split points belong to the same cluster, and the right split points belong to the same cluster. Each contact cluster is characterized by the number of supporting experiments and by the total number of supporting reads. Contact clusters are assigned to the positive or negative strand based on the orientation of transcripts, to which they belong.

To identify pairs of nested contact clusters (PNCC) we devised an efficient procedure based on bedtools intersect program (Quinlan, 2014). First, we identify two lists of pairs of contact clusters A and B such that (1) the left segment of B belongs to the window $[10;100]$ nts from the left segment of A, and (2) the right segment of B belongs to the window $[-100;-10]$ from right segment of A. These lists are intersected using bedtools intersect as follows. We create two files in BED format that contain the coordinates of the left and the right segments of each contact cluster, respectively. Similarly, we create two additional files in BED format that store the coordinates of the left and the right segments with indents $[10;100]$ and $[-100;-10]$, respectively. Pairwise intersection of these files and additional matching the identifiers of the constituent contacts yield the desired list of PNCC. This procedure is illustrated in Fig. S1.

The list of PNCC is filtered according to the following conditions: (1) read support of both clusters is $\ge$3; (2) distance between clusters is between 10 and 50 nts; and (3) none of the folding intervals (intervals 1–2 and 3–4 in Fig. 1) intersects with genomic repeats annotated within RepeatMasker track of UCSC Genome Browser including SINE (short interspersed nuclear elements) and LINE (long interspersed nuclear elements), and other repeat types (Jurka, 2000). The condition (1) roughly corresponds to the significance cutoff for PCCRs by the frequency of forked eCLIP peaks (Margasyuk et al., 2023a); the condition (2) represents the range of lengths of PCCRs (Kalmykova et al., 2021); the condition (3) was introduced to filter out abundant RNA structures formed by low complexity regions as in previous works (Kalmykova et al., 2021; Pervouchine, 2014). The nucleotide sequences of the folding intervals are extracted from the genome using bedtools getfasta program (Quinlan, 2014).

RNA structure prediction and classification

The secondary structure between folding intervals were predicted for each PNCC by calling the PREPH program as fold.py -k 3 -a 3 -e -1 -u False -d 2. That is, RNA structure was required to have a minimal helix length of 3 nts and at most two GT base pairs in each such helix, similar settings that were used earlier to predict PCCRs (Kalmykova et al., 2021). There was no default threshold on the free energy and no prediction of suboptimal structures. The results of PREPH were parsed by custom scripts to extract the coordinates of the complementary regions (called left and right “handles”), their base-pairing scheme in dot-bracket notation, and the free energy of hybridization ( $\mathrm{\Delta }G$). The predicted RNA structure was categorized as conserved if both complementary parts were located within the set of conserved RNA elements (phastConsElements track for the alignment of 100 vertebrates genomes to the human genome (Siepel et al., 2005)), or non-conserved otherwise.

Compensatory substitutions

To assess the abundance and statistical significance of compensatory substitutions for PHRIC predictions, including those outside of conserved 100-vertebrate alignment blocks, we employed the procedure that was outlined previously (Kalmykova et al., 2021). Multiple sequence alignments (MSA) of 46 mammalian genomes were downloaded from the UCSC Genome Browser website in MAF format (Kent et al., 2002). We analyzed pairs of alignment blocks that were cut out from the MSA by the predicted complementary regions. MSA columns containing more than 80% gaps and rows containing more than 10% gaps were removed. The resulting MSA pairs were passed to the R-scape v1.2.340 software as explained earlier (Kalmykova et al., 2021) after merging through a spacer containing 10 adenine nucleotides along with the restricted phylogenetic tree. The E-value of the structure was calculated as the product of E-values of its constituent base pairs that were reported by R-scape.

PHRIC pipeline

PHRIC is implemented in a reproducible and scalable workflow management system Snakemake. It is available publicly through the GitHub repository (Margasyuk et al., 2023b).

Pairs of nested contact clusters

RIC-seq experiments in seven human cell lines (see Methods) containing 55–170 million raw reads per replicate were analyzed to extract split reads encoding RNA contacts (see Methods). In total, 2–10 million split positions supported by 15–40 million individual split reads were obtained per cell line (Fig. S2). Further analysis was confined to intragenic split reads with distance between split positions of at least 200 nts. Split points from all RIC-seq experiments were pooled and clustered using single-linkage clustering with the distance threshold of 10 nts, resulting in ~35 millions of RNA contact clusters. Due to inherent sparsity of RIC-seq data, we chose to merge split points from all RIC-seq experiments rather than to analyze them separately (Cao et al., 2021; Cai et al., 2020). Consequently, the variation in the number of split reads across cell lines did not affect the downstream analysis. Each cluster of RNA contacts was characterized by the set of RIC-seq experiments, in which it was supported, and by the total number of supporting split reads from all RIC-seq experiments. RNA contact clusters were observed in on average $1.34$ RIC-seq experiments and were supported by on average $2.05$ split reads. Then, we found PNCC with distances between contacts ranging from 10 to 100 nts (Fig. 1). The resulting set of ~2.6 million PNCC was further filtered to constrain the distance between the outer and the inner contacts to be between 10 and 50 nts while requiring them to be supported by at least three split reads and excluding clusters that intersect annotated genomic repeats (see Methods). This resulted in ~29,000 PNCC, which were next passed to PREPH (Kalmykova et al., 2021) to predict RNA structure.

RNA structure properties

PREPH predicts nearly-perfect stretches of complementary nucleotides in a pair of input sequences using the dynamic programming matrix based on precomputed helix energies for all $k$-mers and energies of short internal loops and bulges (Kalmykova et al., 2021). In application to the set of ~29,000 PNCC supported by at least three split reads, it yielded 11,998 predicted RNA structures at the equilibrium free energy cutoff $\mathrm{\Delta }G<-15$ kcal/mol having a decaying $\mathrm{\Delta }G$ distribution with the median of $-23.1$ kcal/mol (Fig. 2A). We subdivided the predicted RNA structures into three categories: intronic, in which both sequences were located entirely in introns; exonic, in which both sequences were located entirely in exons; and mixed, which was comprised of cases when one of the sequences overlapped a splice site or an exonic sequence contacted an intronic sequence. In all three groups, approximately 40% of the predictions were supported by 5–10 reads, roughly 40% of the predictions were supported by 10–20 reads, and 20% of the predictions were supported by more than 20 reads (Table 1). The boundaries defining these groups represent natural cutoffs subdividing the set of predictions into subsets of roughly the same cardinality with increasing read support level. Remarkably, more than 70% of the predicted structures were located outside of genomic blocks conserved in 100 vertebrates.

Our expectation was that PNCC with higher levels of RIC-seq support yield more stable RNA structures. Indeed, the median absolute value of $\mathrm{\Delta }G$ increased with increasing $r$, the total number of supporting reads (Fig. 2B), but the magnitude of this increase was very small (on average, 0.03 kcal/mol with each additional supporting read). We next asked if the observed free energy distribution is non-random with respect to the rewired control (Pervouchine et al., 2012), in which the interacting sequences were randomly exchanged (see Methods). To control for the dinucleotide frequencies, which affect $\mathrm{\Delta }G$ values, we separately estimated the free energy of interaction in PNCC, in which one of the nucleotide sequences was reverse complemented. Both the rewired pairs and the reverse complemented control resulted in significantly lower median $\mathrm{\Delta }G$ values (Fig. 2C). This findings indicate that RNA structures predicted by the PHRIC pipeline are more stable than randomly occurring structures, and that their stability correlates with RIC-seq read support.

RNA structure in exons and introns

Previous studies of RNA secondary structure and RBP interaction landscapes in eukaryotic nuclei demonstrated that introns tend to be more structured than exons (Gosai et al., 2015; Sun et al., 2019; Zafrir & Tuller, 2015). However, these studies assessed the propensity of RNA bases to be involved in local RNA structure and did not consider long-range RNA structure organization. We revisited this question by subdividing the exonic, intronic, and mixed into four free energy groups, 15–20, 20–25, 25–30, and >30 kcal/mol (by absolute value) as in Kalmykova et al. (2021) and, on one hand, into groups with high ( $r\ge 12$) and low ( $r<12$) read support according to PNCC, to which they belonged. The boundary of $r=12$ is equal to the median of the read support distribution.

The intronic RNA structures were characterized by a higher proportion of predictions with free energies exceeding $25$ kcal/mol by absolute value (Table 2). Furthermore, the intronic RNA structures with high read support had a significantly larger $\mathrm{\Delta }G$ values than those with low read support (P-value < 0.1%, Mann-Whitney test), while in exonic and mixed groups the difference between the high and the low read support groups was not significant (Fig. 3A). Next, we compared the $\mathrm{\Delta }G$ values between the observed and the rewired sets of PNCC in the intronic, exonic, and other groups. The median free energies in all three groups were significantly larger as compared to the rewired control set (P-value < 0.1%, Mann-Whitney test), with a larger magnitude of difference for the intronic group than for the exonic group (Fig. 3B). Finally, we considered the distribution of $\mathrm{\Delta }\mathrm{\Delta }G=\mathrm{\Delta }{G}_{obs}-\mathrm{\Delta }{G}_{RC}$ values in the matched set of the RNA structures that were actually observed (with equilibrium free energy $\mathrm{\Delta }{G}_{obs}$) and the control set, in which one of the sequences was reverse complemented (with equilibrium free energy $\mathrm{\Delta }{G}_{RC}$). Again, intronic RNA structures were more stable with respect to the reverse complemented control compared to exonic and mixed groups (Fig. 3C), thus confirming that introns have a higher propensity to form stable long-range RNA structures between each other.

Next, we focused on a subset of intronic RNA structures that loop out at least one annotated exon and computed the inclusion rates $\mathrm{\Psi }$ of these exons across all cell lines (Fig. 4A). The median $\mathrm{\Psi }$ value decreased significantly for exons that are looped out by intronic RNA structures with increasing the read support (P-value < 1%, Mann-Whitney test), in full agreement with previous findings that exon inclusion drops with increasing stability of the surrounding RNA structure (Kalmykova et al., 2021; Margasyuk et al., 2023a), which positively correlates with the read support as we showed before (Fig. 2B).

A comparison of $\mathrm{\Delta }G$ values between RNA structures located in conserved and non-conserved parts of exonic and intronic regions revealed that non-conserved intronic RNA structures were as stable as conserved intronic RNA structures, while non-conserved exonic RNA structures were on average even more stable than conserved exonic RNA structures (Fig. 4B). This indicates that exonic sequences lacking constraints on maintaining the encoded aminoacid sequence can evolve more stable RNA structures.

Finally, we subdivided intronic and exonic RNA structures into three classes corresponding to interactions within the same intron or exon, adjacent (i.e., consecutive) introns or exons, and distant (i.e., non-consecutive) introns or exons. Exonic RNA structures distributed almost equally likely between these groups, while intronic RNA structures have a strong preference to concentrate in the same intron (Fig. 4C). The distance between complementary sequences did not confound this result as the distributions of distances were almost the same (Fig. S3). This remarkable observation now provides an experiment-derived evidence for the tendency of RNA structures to concentrate within the same intron, which was suggested earlier by several bioinformatic studies (Kalmykova et al., 2021; Pervouchine et al., 2012).

Case studies

In this section, we discuss a few examples of long-range RNA structures from the list of 11,998 PHRIC predictions. The entire list is available as File S1, which can be visualized through UCSC genome browser along with RIC-seq contacts in File S2 (Margasyuk et al., 2023c).

The human GANAB gene encodes the glucosidase II $\alpha$ subunit and is involved in autosomal-dominant polycystic kidney and liver disease (Porath et al., 2016; Besse et al., 2018). One of its internal exons, exon 6, is spliced alternatively. We detected two PNCC surrounding exon 6 that are supported by a large number of RIC-seq split reads (Fig. 5A). PHRIC pipeline predicts two pairs of complementary regions that constitute long-range RNA structures with free energies of $-26.3$ and $-22.1$ kcal/mol, respectively. We hypothesize that complementary base pairing mediated by these RNA structures are responsible for alternative splicing of exon 6.

Another example of PHRIC predictions is the RNA structure in the human NDUFB5 gene, which encodes a subunit of the multisubunit NADH-ubiquinone oxidoreductase (complex I). Three transcript variants encoding different splice isoforms have been found for this gene, and two of them differ by alternative inclusion of exon 7 (Fig. 5B). We detected two pairs of complementary intronic sequences that are supported by clusters of RIC-seq contacts. Mutually exclusive splicing of exon 7 could be modulated by these RNA structures with free energies $-28.8$ and $-22.1$ kcal/mol.

Finally, we discuss the ZNF655 gene, which encodes a zinc finger protein that is possibly involved in transcriptional regulation. It accelerates the progression of pancreatic cancer by promoting the binding of E2F1 and CDK1 (Shao et al., 2022). It contains a cassette exon 4, which is surrounded by a pair of complementary sequences capable forming a duplex with free energy $-25.9$ kcal/mol. Alternative splicing of this exon could also be modulated by RNA structure.

Evolutionary signatures beyond sequence conservation in vertebrates

Our earlier study has identified pairs of conserved complementary regions (Kalmykova et al., 2021) within the so-called conserved RNA elements, which are derived from Multiz sequence alignments of 100 vertebrate genomes using phylo-HMM (Felsenstein & Churchill, 1996; Blanchette et al., 2004). The generative model of phylo-HMM contains a state for conserved sites and a state for non-conserved sites, transitions between which determine the borders of conserved RNA elements. These borders and the set of conserved RNA elements itself vary depending on the set of species that were passed to the model as an input.

In this section, we asked whether the evolutionary signatures of RNA structure can be extracted directly from genomic multiple sequence alignments. Towards this goal, we restricted our analysis to 46 mammalian genomes including the human genome (Blanchette et al., 2004; Kent et al., 2002), and applied the R-scape program (Rivas, Clements & Eddy, 2017), which scores independent occurrence of compensatory substitutions on different branches of the phylogenetic tree, to the alignment blocks cut out by PHRIC predictions. The statistical significance of pairwise covariations in each predicted structure was estimated as a product of E-values reported by R-scape for all its base pairs. Out of 11,998 pairs of complementary regions originally predicted by PHRIC, 11,224 had at least one base pair with E-value < 1, and only 308 pairs had a significant E-value (below 5%) after Benjamini-Hochberg adjustment for multiple testing.

As in Kalmykova et al. (2021), we found that RNA structures with significant compensatory substitutions had significantly larger equilibrium free energies (by absolute value) than RNA structures with non-significant covariations (P-value < 0.001, Mann-Whitney test) (Fig. 6A). The median lengths of the base-paired regions did not differ significantly between these two sets (P-value = 0.15, Mann-Whitney test), thus excluding the possibility that longer RNA structures contribute simultaneously to both $\mathrm{\Delta }G$ and the E-value.

An RNA structure with many independent compensatory substitutions ( $E-\mathrm{v}\mathrm{a}\mathrm{l}\mathrm{u}\mathrm{e}=1.6\cdot {10}^{-14}$) was detected in NFAT5 gene, a member of the nuclear factors of activated T cells family of transcription factors, which plays a central role in inducible gene transcription during the immune response (Neuhofer, 2010). The intron spanning between exons 7 and 8 of NFAT5 contains a pair of complementary sequences, R1 and R2, which are strongly supported by RIC-seq as PNCC, but fall outside of conserved RNA elements in vertebrates (Fig. 6B). Manual inspection of the multiple sequence alignment of mammalian homologs revealed that intronic sequences corresponding to R1 and R2 are missing in rodents, but present in primates, canids, and African elephant, with multiple independent compensatory substitutions. Remarkably, the structure formed by R1 and R2 is nested within a larger conserved structure (id355386, $\mathrm{\Delta }G=-40.1$ kcal/mol) that was predicted for this gene earlier (Kalmykova et al., 2021).