Many human long noncoding RNAs are localized in the nucleus and can potentially participate in chromatin formation and remodeling ¹ . Recently, technologies such as ChIRP ² , ChRIP ³ , ChOP ⁴ , CHART ⁵ , RAP ⁶ , MARGI ⁷ and GRID ⁸ have been developed to map the genomic interacting sites of various lncRNAs. RNA can interact with chromatin by associating with DNA binding proteins, nascent transcripts, single-stranded or double-stranded DNA, forming R-loops or triple helices, respectively. Growing body of evidence shows that RNA-DNA triplex formation based on the Hoogsteen ⁹ base pairing rules plays important role in RNA-chromatin interactions. Several studies have provided in vitro and in vivo evidence for the existence and biological relevance of triplexes, including pRNA ¹⁰ , Fendrr ¹¹ , Khps1 ¹² , PARTICLE ¹³ , and MEG3 ⁴ .

Computational analyses have revealed that a large population of triplex-forming motifs is present across the human genome with the majority of annotated genes containing at least one triplex target site (TTS), preferentially in regulatory gene regions ¹⁴ , ¹⁵ . Considering the large number of purine-rich sequences in the genome, triplex-mediated targeting of lncRNAs and associated proteins to distinct genomic loci is very likely a commonly used mechanism of gene regulation. Still, there are only a few bioinformatic studies of triplex-based RNA-DNA interactions on the genome-wide scale ¹⁶ .

Here we analyze the genomic regions that are known to interact with the MEG3 lncRNA (the ChOP-seq peaks) or three different short oligos, corresponding to the DNA binding domains (DBDs) of MEG3 and GATA6-AS lncRNAs (the Capture-seq peaks). The current literature usually assumes that the triplex-based interactions have high sequence specificity and each triplex forming oligonucleotide (a transcript region) has a distinct set of genomic binding sites ("triplexome"). We investigate whether all the DNA sites capable of triplex formation are specific enough to be regulated by one particular RNA only or whether different transcripts may have shared TTSs. Our computational analysis revealed a group of genomic regions that may have a very high propensity for triplex formation with a wide range of different RNAs. Therefore, we named such DNA sequences "universal TTSs". We also attempted to reveal the features of these sequences that may be responsible for the observed phenomenon.

The genomic coordinates of the 6837 MEG3 binding sites ⁴ (ChOP-seq peaks) were mapped from the hg19 to the hg38 human genome version using liftOver ¹⁷ (Nov 7, 2017 version). Next, from the 6800 successfully converted peaks we removed two cases corresponding to the genomic regions with ambiguous base-pairs (N) keeping the 6798 ChOP-seq peaks for the analysis. Additionally, to simulate the genomic background 6798 control regions with the lengths matching the selected ChOP-seq peaks were randomly sampled from the human genome using the bedtools ¹⁸ (version 2.27.1, see Supplementary Figure 1).

Triplex-based interactions were predicted by the Triplexator ¹⁴ (version 1.3.2) with the following parameters: -fr off -l 10 -e 10. These values were optimized so that the tool could predict binding between all three RNA-DNA sequence pairs that have been validated in vitro in the original study ⁴ ( Supplementary Table 1). To detect the statistically significant RNA-DNA interactions we developed a method to estimate a p-value from the number of predicted triple helices. Since the MEG3 peaks have different lengths, the expected number of triplexes (i.e. the parameter of the Poisson distribution) is computed based on the lengths of the input RNA and DNA sequences (see below).

The MEG3 binding sites have been identified in the triple negative breast cancer cell line BT-549 ⁴ . To identify all the genes expressed in this cell line we used the RNA-seq data from the control knockdown experiment (ERR652847). The reads were aligned to the human genome (hg38) using HISAT2 ¹⁹ (version 2.1.0) and the number of reads corresponding to each GENCODE ²⁰ (version 28) transcript was calculated by the HTSeq-count tool ²¹ (version 0.10.0). Next, the RPKM values were computed as RPKM = C /( N × L), where C is the number of reads aligned to all the transcript exons, N is the total number of mapped reads (in millions) and L is the transcript length (in kilobases). The most highly expressed (in terms of RPKM) isoform of each gene was considered only. Next, 153 expressed transcripts with the length and GC content similar to the MEG3 lncRNA (NR_002766.2) were selected using the RANN (version 2.6) R package ( Supplementary Figure 2). Additionally, 153 random RNA sequences were obtained by di-nucleotide shuffling the original MEG3 transcript using the uShuffle tool ²² .

All the heatmaps were generated using the Complex-Heatmap ²³ R package. The alignment of the RNA oligo sequences was obtained by the MUSCLE ²⁴ (version 3.8). The locations of the RepeatMasker repeats in the human genome were downloaded from the UCSC Genome Browser ¹⁷ .

We used Triplexator to predict possible triple helices between the full length MEG3 transcript and the 6798 experimentally identified ChOP-seq peaks as well as 6798 control DNA regions (see Methods). The genomic sites with a statistically significant number of triplexes were identified in each DNA set by our custom probabilistic approach (see Methods). As anticipated, more statistically significant (Bonferroni adjusted p-value < 0.01) interactions with the MEG3 lncRNA were predicted for the ChOP-seq peaks than for the control regions. Namely, the interactions with the 3825 (56.3%) ChOP-seq peaks were classified as statistically significant ( Figure 1A, left) while there were only 617 (9.1%, odds ratio test p-value < 2.2 × 10 ⁻¹⁶) such cases among the control DNA regions ( Figure 1B, left). Since the ChOP-seq method has detected RNA contacts with the chromatin (and not the naked DNA) the obtained binding sites can correspond to several different interaction mechanisms including direct RNA-DNA interactions via triple helices or R-loops, RNA-RNA hybridization with nascent transcripts as well as bindings to nuclear proteins. This may be the reason that many MEG3 binding sites did not produce statistically significant predictions with the MEG3 lncRNA. Therefore, these results supported the original conclusion that the MEG3 lncRNA is able to directly interact with the genomic DNA via triple helices.

To check the ability of other RNAs to form triplexes with MEG3 binding sites, we applied Triplexator to a set of 153 expressed transcripts (see Methods). Surprisingly, 65 analyzed RNAs showed results similar to MEG3 lncRNA – they were predicted to form statistically significant interactions with the majority (> 50%) of the ChOP-seq peaks ( Figure 1A, middle). To further investigate possible interactions with the ChOP-seq peaks, 153 artificial sequences were generated by di-nucleotide shuffling of the MEG3 transcript (see Methods). Strikingly, these random "RNAs" produced statistically significant number of triplexes with 39% and 9% of the ChOP-seq peaks and control DNA regions, respectively ( Figure 1A, B, right). These results indicated that the set of the ChOP-seq peaks was different from the randomly sampled genomic sites in that it contained a number of DNA sequences that may be able to interact not only with the MEG3 lncRNA, but with other RNAs as well.

Based on these observations we hypothesized that some of the MEG3-bound genomic sites may be ’universal’, i.e. they may have a potential to form multiple triplexes with a number of different RNAs. Analysis of the Triplexator predictions obtained for the 153 expressed RNAs revealed 1107 (16.3%) ChOP-seq peaks that were predicted to form statistically significant number of triplexes with more than 90% of the analyzed transcripts ( Supplementary Figure 4A). In contrast, the genomic background contained only 41 (0.6%) such sites ( Supplementary Figure 4B). Due to the predicted ability of these genomic regions to form triple helices with various RNAs, we called them "universal triplex target sites (TTSs)". Notably, the identified universal TTSs produced strong p-values for the MEG3 lncRNA as well as for the 153 MEG3 shuffled sequences ( Figure 1C, D). Thus, according to our predictions some of the DNA sequences were more prone to formation of triple helices with different long RNAs and a number of such genomic regions were present among the experimentally identified MEG3 binding sites (ChOP-seq peaks).

To further investigate the predicted ability of the universal TTSs to bind various RNAs we analyzed the results of the recent Capture-seq experiment ²⁵ . This in vitro study has determined genomic binding sites of three different short RNA oligos that corresponded to the DNA-binding domains (DBDs) of the MEG3 and GATA6-AS lncRNAs. Namely, the MEG3_13_41, GATA6_AS_78_118 and MEG3_839_890 oligos were 28, 40 and 48 nt long, respectively. Since the experiment has been performed on the RNA- and protein-free ("naked") genomic DNA, the majority of the identified interactions have been assumed to be direct and mediated by triple helices. Comparison of the genomic coordinates corresponding to the identified target DNA fragments demonstrated that most of the interactions were specific to one oligo only ( Supplementary Figure 5). This can be explained by the fact that the oligo sequences had limited similarity with each other (the identities between the MEG3_13_41- GATA6_AS_78_118, MEG3_839_890- GATA6_AS_78_118 and MEG3_13_41- MEG3_839_890 oligo pairs were 33%, 40% and 25%, respectively – Supplementary Figure 6). Still, 5379 genomic regions were captured by each of the three oligos ( Supplementary Figure 5). Thus, we expected that this set of ’shared Capture-seq peaks’ can be enriched with the potential universal TTSs. To check this we predicted their possible interactions with the 153 RNA sequences that were used in the analysis of the ChOP-seq peaks (see above). Indeed, 2151 (40%) shared Capture-seq peaks were classified as universal TTSs – they had statistically significant number of triplexes with most (> 90%) of the analyzed transcripts ( Figure 1E and Supplementary Figure 4C). Therefore, the fact that the experimentally identified set of shared Capture-seq peaks contained such a high fraction of the universal TTSs indirectly confirmed the predicted property of these special genomic loci.

Finally, we attempted to reveal the features of the universal TTSs that may allow them to interact with several different RNAs. For this purpose we compared sequence composition of the universal and all the other (i.e. non-universal) genomic regions from each set. While the GC content of the universal and non-universal DNA sequences were similar, the universal TTSs had higher purine (G or A) and, especially, poly-purine content (see Supplementary Figure 7 and Methods for the definitions). To find out the origin of these poly-purine elements we analyzed the classes of the overlapping genomic repeats. All three sets of the universal TTSs were enriched with the purine-rich low complexity regions, LCRs ( Figure 1F and Supplementary Figure 8). Additional analysis of several universal TTSs confirmed that these LCRs were predicted to form multiple triple helices with the majority of the analyzed transcripts (see Supplementary Figure 9 for representative cases). Therefore, the presence of the purine-rich low complexity elements was the characteristic property of the universal TTSs that potentially allowed them to interact with a wide range of different RNAs. All together our results suggested the existence of a special type of genomic loci that may function as RNA-binding hubs.

The importance of triplex-dependent gene regulation in the genomes of higher organisms is becoming a generally accepted concept. Here we performed a large-scale bioinformatic analysis of the genomic regions (ChOP-seq and Capture-seq peaks) that have been shown experimentally to interact with particular RNAs (MEG3 lncRNA or short oligos). To filter out not significant Triplexator predictions, the statistical significance of every RNA-DNA interaction was estimated from the Poisson distribution. To our surprise for some genomic regions (that we called "universal triplex target sites") Triplexator predicted statistically significant (adjusted p-value < 0.01) number of triplexes with the majority (> 90%) of the analyzed transcripts. According to our analysis universal TTSs are quite rare in the human genome – there were only 0.6% of them among the 6798 randomly sampled regions. On the other hand, 16.3% of the experimentally identified MEG3 binding sites (ChOP-seq peaks) were classified as universal TTSs. Additionally, genomic regions that have been shown to form triplexes with three different oligos (shared Capture-seq peaks) contained 40% of the universal TTSs. All three sets of the identified universal TTSs were enriched with the purine rich low complexity regions.

The theoretical possibility of the universal TTS existence comes from the degeneracy of the Hoogsteen rules ⁹ . In fact, the triplex-based interaction can be formed in both orientations (parallel and anti-parallel) and it involves only purines (G or A) in the DNA. Additionally, the DNA guanine and adenine can bind to RNA guanine and uracil, respectively, in both orientations while the A::A pairing occurs in the anti-parallel orientation only. This makes the long poly-purine elements a possible targets for a number of different RNA oligos.

One of the possible and actively discussed roles of the chromatin bound RNAs (including lncRNAs) is to bring different chromosomal parts together to enable the remote DNA-DNA contacts ⁸ . Moreover, it has recently been shown that RNAs originating from super-enhancers form triplexes at distant regions ²⁶ . Therefore, it is possible that universal TTSs may facilitate distal enhancer-promoter interactions via engagement with the same enhancer RNA. In line with this hypothesis, we observed the statistical significant enrichment of the universal Capture-seq peaks near (< 1 kb) the transcription start sites (TSSs) of the annotated genes ( Supplementary Figure 10A). However, the computationally predicted universal ChOP-seq and background TTSs did not have such trend ( Supplementary Figure 10B,C). Thus, the experimentally identified shared Capture-seq peaks may be more suitable for subsequent functional validation of the universal TTSs in living cells.

Importantly, the current computational analysis has a number of limitations. Namely, the triplex-based interactions of the full length transcripts were predicted without taking their secondary structure into account. We are not aware of any bioinformatics tools that would be able to produce such predictions. Moreover, cellular localization of the 153 selected expressed transcripts as well as DNA binding proteins and chromatin compaction were not considered. Therefore, our simulations are more similar to the in vitro Capture-seq experiments with short oligos than to the interactions of long transcripts with the chromatin inside the nucleus. Comprehensive identification of all the RNA-DNA interactions obtained by high throughput experimental methods may clarify the predicted functionality of the universal TTSs in the cell. Although a few methods for this task have recently been developed, the length of the sequencing reads (e.g. about 40 bp of DNA in case of GRID-seq) does not allow to reliably determine interactions with the long low complexity regions (including universal TTSs). We are looking forward to the new high quality experimental data to gain further insight into the triplex-based RNA-chromatin interactions in vivo.