Keywords
MEG3 lncRNA, triple helix, triplex target sites (TTS)
MEG3 lncRNA, triple helix, triplex target sites (TTS)
As it was suggested by the reviewers, we significantly increased the number of the RNA sequences (from 6 to 306) that were used for the analysis and used the exact locations for all the peaks (instead of 3kb bins that were originally used). Additionally, we developed a probabilistic approach to estimate the statistical significance of each RNA-DNA interaction. Finally, we added the analysis of the recently published data obtained using the Capture-seq experiment.
The data has been expanded and uploaded to Zenodo and GitHub.
We have updated the title of the article to: ‘Purine-rich low complexity regions are potential RNA binding hubs in the human genome’.
Nevertheless, with all these changes our main conclusion remained the same. Namely, we proposed that the purine rich low complexity genomic regions may have an ability to interact with various different RNAs.
We also speculate about possible biological role of these special genomic loci.
See the authors' detailed response to the review by Andrey A. Mironov
See the authors' detailed response to the review by Ingrid Grummt
See the authors' detailed response to the review by Hao Zhu
Many human long noncoding RNAs are localized in the nucleus and can potentially participate in chromatin formation and remodeling1. Recently, technologies such as ChIRP2, ChRIP3, ChOP4, CHART5, RAP6, MARGI7 and GRID8 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 Hoogsteen9 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 pRNA10, Fendrr11, Khps112, PARTICLE13, and MEG34.
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 regions14,15. 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 scale16.
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 sites4 (ChOP-seq peaks) were mapped from the hg19 to the hg38 human genome version using liftOver17 (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 bedtools18 (version 2.27.1, see Supplementary Figure 1).
Triplex-based interactions were predicted by the Triplexator14 (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 study4 (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-5494. 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 HISAT219 (version 2.1.0) and the number of reads corresponding to each GENCODE20 (version 28) transcript was calculated by the HTSeq-count tool21 (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 tool22.
All the heatmaps were generated using the Complex-Heatmap23 R package. The alignment of the RNA oligo sequences was obtained by the MUSCLE24 (version 3.8). The locations of the RepeatMasker repeats in the human genome were downloaded from the UCSC Genome Browser17.
For a given pair of RNA and DNA sequences, Triplexator outputs all the possible triple helices that satisfy the user-defined thresholds. Notably, the number of the predicted triplexes increases with the lengths of input sequences (Supplementary Figure 3A). To account for this dependence the normalized number of triplexes (i.e. the "triplex potential" or tpot) is also computed by the Triplexator. Although this allows to compare triplexes predicted for RNA-DNA pairs with different lengths, it does not provide information about significance of these interactions.
To estimate the probability to observe a particular number of triplexes by chance (e.g. from the random sequences with the same lengths) we analyzed the average number of predicted triplexes between random sequences of various lengths. Namely, we considered four different RNA lengths (LRNA = {500, 1000, 1500, 2000}) and ten different DNA lengths (LDNA = {200, 400, ..., 1800, 2000}). For each of the 40 combinations of (LRNA, LDNA), 100 random sequences with the length of LRNA and 100 random sequences with the length of LDNA were generated (with the equal frequencies for all the four nucleotides).
For each RNA-DNA pair triple helices were predicted by the Triplexator with the parameters optimized for MEG3 lncRNA (-fr off -l 10 -e 10). Next, for every combination of (LRNA, LDNA) the average number of predicted triplexes (λ) was computed from all the 10000 predictions (Supplementary Table 2). Finally, a linear regression model for λ was fitted to all the obtained values (adjusted R2 = 87%, see Supplementary Figure 3B):
where θ0 = −0.688, θ1 = 5.37 × 10−4 and θ2 = 6.03 × 10−4.Thus, the statistical significance of the number of triple helices Ntpx, predicted between RNA of length LRNA and DNA of length LDNA, can be estimated as follows. First, the expected average number of predicted triplexes (λ) is computed from the equation (1). Next, the expected distribution of the number of predicted triplexes (H0) is simulated by the Poisson distribution with the obtained value of λ (Supplementary Figure 3C). Finally, the p-value of the observed number of triple helices (Ntpx) can be estimated as follows:
Importantly, the Ntpx value is taken from the "Total (abs)" column of the triplex_search.summary file. The same value is used by the Triplexator to compute its "triplex potential" (tpot). The ’Total (abs)’ is the total number of all possible triplexes that satisfy the user-defined thresholds (overlaps are allowed). Thus, for a single triplex longer than the minimal length (10 in our settings), the ’Total (abs)’ value may be greater than 1. For example, 11 bp DNA fragment 5’-GAGAGAGAGAG-3’ and 11 nt RNA oligo 5’-GAGAGAGAGAG-3’ can interact with each other forming one long anti-parallel triplex without mismatches. However, with the minimal allowed triplex length set to 10, the Ntpx is equal to 3. This includes the long triplex of length 11 as well as the two triplexes of length 10 without the first or the last position of the long triplex. Therefore, a single long triplex is likely to produce a large Ntpx value and, consequently, a statistically significant p-value.
Due to the properties of the Watson-Crick base pairing model, the GC content of a sequence corresponding to the forward (+) DNA strand is equal to the GC content of the reverse (-) DNA strand. However, the GA content is more important for triplex based interactions because RNA can only form triple helices with the purines in the DNA. In contrast to the GC content, the GA content can be different between the DNA strands. Moreover, if one strand is purine rich, the other strand is automatically purine poor. For example, for the DNA sequence 5’-GGGGGAGA-3’ the purine content of the direct strand is 100%, while the other strand (3’-CCCCCTCT-5’) has no purines at all (i.e. its purine content is 0%).
Thus, we define the purine content of a given DNA fragment as the maximum value between the two strands, i.e.:
It should be noted that the purine content computed by the formula (3) is always ≥ 50%. To work with a measure that is defined from 0% to 100%, we introduce the poly-purine content. We define a poly-purine element Ps as a continuous stretch of 10 or more purines located on the DNA strand s (where s = {+, −}) . For a given DNA sequence that has N+ poly-purine elements on the forward strand and N− poly-purine elements on the reverse strand, the poly-purine content is computed as follows:
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−16) 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 experiment25. 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 rules9. 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 contacts8. Moreover, it has recently been shown that RNAs originating from super-enhancers form triplexes at distant regions26. 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.
Zenodo: vanya-antonov/universal_tts: The initial release of the code, data files and images related to universal TTSs. http://doi.org/10.5281/zenodo.265480027
This project contains the following underlying data:
universal_tts-v1.0.0.zip?download=1.zip
– data (folder containing underlying data, description of individual files can be found in Supplementary Table 3)
Zenodo: vanya-antonov/universal_tts: The initial release of the code, data files and images related to universal TTSs. http://doi.org/10.5281/zenodo.265480027
This project contains the following extended data:
Data and code are available under the terms of GNU General Public License version 3 (GPL-3.0).
This work was partially supported by Russian Science Foundation (grant 14-15-30002 to YAM).
The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Supplementary File 1. File containing Supplementary figure 1–Supplementary figure 10, and Supplementary table 1–Supplementary table 3.
We are thankful to Dr. Chandrasekhar Kanduri (University of Gothenburg, Sweden) for providing the original coordinates of the ChOP-seq peaks.
Views | Downloads | |
---|---|---|
F1000Research | - | - |
PubMed Central
Data from PMC are received and updated monthly.
|
- | - |
Competing Interests: No competing interests were disclosed.
Competing Interests: No competing interests were disclosed.
Reviewer Expertise: Regulation of gene expression by noncoding RNA
Is the work clearly and accurately presented and does it cite the current literature?
Yes
Is the study design appropriate and is the work technically sound?
No
Are sufficient details of methods and analysis provided to allow replication by others?
Partly
If applicable, is the statistical analysis and its interpretation appropriate?
No
Are all the source data underlying the results available to ensure full reproducibility?
Yes
Are the conclusions drawn adequately supported by the results?
No
Competing Interests: No competing interests were disclosed.
Reviewer Expertise: Genome analysis, lncRNA analysis, molecular evolution
Is the work clearly and accurately presented and does it cite the current literature?
Yes
Is the study design appropriate and is the work technically sound?
Yes
Are sufficient details of methods and analysis provided to allow replication by others?
Yes
If applicable, is the statistical analysis and its interpretation appropriate?
Yes
Are all the source data underlying the results available to ensure full reproducibility?
Yes
Are the conclusions drawn adequately supported by the results?
Yes
References
1. Mondal T, Subhash S, Vaid R, Enroth S, et al.: MEG3 long noncoding RNA regulates the TGF-β pathway genes through formation of RNA-DNA triplex structures.Nat Commun. 2015; 6: 7743 PubMed Abstract | Publisher Full TextCompeting Interests: No competing interests were disclosed.
Is the work clearly and accurately presented and does it cite the current literature?
No
Is the study design appropriate and is the work technically sound?
No
Are sufficient details of methods and analysis provided to allow replication by others?
Partly
If applicable, is the statistical analysis and its interpretation appropriate?
No
Are all the source data underlying the results available to ensure full reproducibility?
Partly
Are the conclusions drawn adequately supported by the results?
No
Competing Interests: No competing interests were disclosed.
Reviewer Expertise: Regulation of gene expression by noncoding RNA
Alongside their report, reviewers assign a status to the article:
Invited Reviewers | |||
---|---|---|---|
1 | 2 | 3 | |
Version 2 (revision) 09 May 19 |
read | read | |
Version 1 17 Jan 18 |
read | read | read |
Provide sufficient details of any financial or non-financial competing interests to enable users to assess whether your comments might lead a reasonable person to question your impartiality. Consider the following examples, but note that this is not an exhaustive list:
Sign up for content alerts and receive a weekly or monthly email with all newly published articles
Already registered? Sign in
The email address should be the one you originally registered with F1000.
You registered with F1000 via Google, so we cannot reset your password.
To sign in, please click here.
If you still need help with your Google account password, please click here.
You registered with F1000 via Facebook, so we cannot reset your password.
To sign in, please click here.
If you still need help with your Facebook account password, please click here.
If your email address is registered with us, we will email you instructions to reset your password.
If you think you should have received this email but it has not arrived, please check your spam filters and/or contact for further assistance.
Comments on this article Comments (0)