ORIGINAL RESEARCH article

Front. Mol. Biosci., 04 April 2018

Sec. RNA

Volume 5 - 2018 | https://doi.org/10.3389/fmolb.2018.00027

Distinct and Modular Organization of Protein Interacting Sites in Long Non-coding RNAs

  • 1. GN Ramachandran Knowledge Center for Genome Informatics, CSIR Institute of Genomics and Integrative Biology, New Delhi, India

  • 2. CSIR Institute of Genomics and Integrative Biology, Academy of Scientific and Innovative Research, New Delhi, India

Abstract

Background: Long non-coding RNAs (lncRNAs), are being reported to be extensively involved in diverse regulatory roles and have exhibited numerous disease associations. LncRNAs modulate their function through interaction with other biomolecules in the cell including DNA, RNA, and proteins. The availability of genome-scale experimental datasets of RNA binding proteins (RBP) motivated us to understand the role of lncRNAs in terms of its interactions with these proteins. In the current report, we demonstrate a comprehensive study of interactions between RBP and lncRNAs at a transcriptome scale through extensive analysis of the crosslinking and immunoprecipitation (CLIP) experimental datasets available for 70 RNA binding proteins.

Results: Our analysis suggests that density of interaction sites for these proteins was significantly higher for specific sub-classes of lncRNAs when compared to protein-coding transcripts. We also observe a positional preference of these RBPs across lncRNA and protein coding transcripts in addition to a significant co-occurrence of RBPs having similar functions, suggesting a modular organization of these elements across lncRNAs.

Conclusion: The significant enrichment of RBP sites across some lncRNA classes is suggestive that these interactions might be important in understanding the functional role of lncRNA. We observed a significant enrichment of RBPs which are involved in functional roles such as silencing, splicing, mRNA processing, and transport, indicating the potential participation of lncRNAs in such processes.

Background

The recent years have seen the discovery of a large number of novel transcripts which belong to the long non-coding RNA (lncRNA) class in humans and other model organisms (Pauli et al., ). This has been largely contributed by the availability of high-throughput methodologies for transcriptome annotation, including tiling microarrays (Hafner et al., ; Furey, ) and deep sequencing (Roberts et al., ). The recent genome-wide analyzes of lncRNA genes in Humans have annotated over 83,215 transcripts from 32,446 lncRNAs genes (Derrien et al., ; Harrow et al., ). The lncRNA superset presently includes a number of sub-classes which include 3 prime overlapping ncRNA, antisense, bidirectional promoter lncRNA, lincRNA, macro lncRNA, miscRNA, non-coding, processed transcripts, pseudogene, retained intron, sense intronic, sense overlapping, and TEC. By definition lncRNAs encompass all transcripts > 200 nucleotides in length and no ORF coding for more than 30 amino acids (Mercer et al., ). The biogenesis and regulation of lncRNAs have not been studied in great detail, though it is believed that they are transcribed majorly by Polymerase II and are capped and polyadenylated (Goodrich and Kugel, ; Gibb et al., ). One particular class of lncRNAs, the large intergenic non-coding RNA has been primarily discovered through their association with epigenetic marks in the genome (Cabili et al., ; Cao, ). We have recently shown extensive similarities and specific dissimilarities in epigenetic regulation of lncRNAs in comparison to protein-coding genes (Sati et al., ). The precise biological function of many of the lncRNAs are not known, though a handful of the candidates have been recently shown to be mechanistically involved in gene regulation and associated with diseases (Wapinski and Chang, ). Recent reports from our group also suggest processing of a subset of lncRNAs to smaller RNAs (Jalali et al., ), and that a subset of lncRNAs could be potentially targeted by microRNAs (Jalali et al., ), thus constituting an intricate and yet poorly understood network of non-coding RNA mediated regulation.

Mechanistically, the characterization of lncRNA could be generalized as a function of its interactions with other biomolecules in the cell: DNA, RNA, protein, and small-molecules (Bhartiya et al., ). Current studies have showed that molecular and computational biology techniques can act as catalyst in discovering lncRNA-mediated regulation via understanding their interactions with different biomolecules (Jalali et al., ). Recent reports have also suggested the possibility of protein-lncRNA interactions and regulatory interactions mediated through them (Kung et al., ). The present understanding of protein-lncRNA interactions are limited to a handful of candidates associated with proteins involved in epigenetic modifications as in the cases of HOTAIR (Gupta et al., ), Anril (Kogo et al., ), and Xist (Arthold et al., ); splicing as in the case of MALAT1 (or NEAT2) (Tripathi et al., ) conserved nuclear ncRNA; transcriptional regulation through interaction with transcription factors as in the case of Gas5 (Kino et al., ) and few other candidates like Meg3 (Zhao et al., ), DHFR (Blume et al., ), and Gomafu (Sheik Mohamed et al., ). It has been recommended that computational methods for predicting protein-RNA interactions, though less accurate, could be potentially used to guide in experiments (Puton et al., ). Recently experimental methodologies to understand protein-RNA interactions on a genomic-scale, including CLIP-seq (Darnell, ) and variants thereof (Hafner et al., ; Jain et al., ; Konig et al., ) has provided insights into the target-sites of a number of RNA binding proteins with much higher resolution (Popov and Gil, ). The availability of genome-scale maps of RNA binding proteins provide a novel opportunity toward understanding patterns of RNA binding proteins interaction sites in different transcript classes and derive clues on the interaction networks, regulation and functional consequences of these interactions.

Recently, Li and coworkers showed the interaction between protein and lncRNAs, in addition to their association with disease causing SNPs. They have deposited all the interaction data in form of bed files in starBase 2.0 database, the same datasets are also included in our current study (Li et al., ). Tartaglia and coworkers have also employed a novel algorithm catRAPID to evaluate the binding tendency of protein with RNAs (Livi et al., ). A similar study by Park et al. has also attempted to explore the possible functions of lncRNAs by focusing at the RBP-lncRNA interactions. LncRNAtor functionally annotates lncRNA molecules based on their expression profiles and co-expression with mRNAs. It also encompass lncRNA's interaction data with 57 RBPs for 5 organisms (Park et al., ).

The functional interactions of lncRNAs could be potentially summarized as the sum total of the interactions between other biomolecules independently or in context of one another. The interaction of lncRNAs with genomic DNA and its involvement in chromatin organization (Lee and Bartolomei, ) and with other RNA species (Salmena et al., ; Bhartiya et al., ; Jalali et al., ) including microRNAs (Jalali et al., ) has been explored at length. Though there have been a number of reports characterizing functional roles of lncRNAs through their association with proteins (Wilusz et al., ), no systematic analysis reports has been published on mapping or on characterizing the functional domains of lncRNAs for protein-binding sites. Our study focuses on providing a platform to explore these interactions at a larger scale using computational approaches to functionally indict the lncRNA molecules.

In the present report, we have performed a comprehensive analysis of 70 experimental RNA binding protein datasets available in the public domain. We have derived the peak information (or the most probable site of interaction between protein and RNA) for these RNA binding protein sites at a genome-scale from doRiNA (Blin et al., ), starBase (Yang et al., ; Li et al., ), and CLIPdb (Yang et al., ) and analyzed their binding sites in lncRNAs and protein coding transcripts. Our analysis suggests 6 lncRNA subtypes (viz; antisense, lincRNA, miscRNA, processed transcripts, retained intron, and sense intronic) to be largely enriched for protein-binding sites compared to other subclasses hence potentially contribute to a novel layer of regulatory interactions mediated through protein-RNA interactions in ncRNA transcripts. Our analysis shows the distribution of RBP binding sites on the lncRNA loci as opposed to only protein coding transcripts. In our study, we also reveal an interesting pattern of positional clustering of RBP target sites in lncRNAs suggesting a modular organization of regulatory sites in lncRNAs. We also propose how the functionally similar proteins co-occur in both protein coding and lncRNA transcripts. To our knowledge, this is the most comprehensive study on the comparison of lncRNA-RBP interactions as opposed to protein coding loci.

Methods

Long non-coding RNA datasets

We used the comprehensive compendium of lncRNAs available from GENCODE Version 24 (August 2015 freeze, GRCh38, Ensembl 83, 84) (http://www.gencodegenes.org/) (Harrow et al., ). The lncRNA dataset had a total of 32,446 genes encompassing 83,215 transcripts having 3,14,672 exons comprising of both Ensembl and Havana annotations. LncRNA transcripts were assigned into 13 biotypes, viz, 3prime overlapping ncRNA, antisense, bidirectional promoter lncRNA, lincRNA, macro lncRNA, miscRNA, non-coding, processed transcripts, pseudogene, retained intron, sense intronic, sense overlapping, and TEC. We also extracted the 19,655 protein coding genes with 79,930 transcripts and their 7,11,466 exons.

Genome scale datasets for protein-RNA interactions

We have compiled and analyzed the protein-RNA interaction datasets from public domain for 70 unique proteins derived from 51 publications across 3 databases (detailed in Table 2). The RBP binding sites were downloaded from 3 databases namely: starBase v2.0 (Yang et al., ; Li et al., ), doRiNA 2.0 (Blin et al., ), and ClipDB v1.0 (last updated: April, 2015) (Yang et al., ). The ClipDB database consisted of datasets analyzed using 4 different softwares PARalyzer (Corcoran et al., ), CIMS (Crosslinking induced mutation site) (Moore et al., ), CITS (Weyn-Vanhentenryck et al., ), and Piranha (Uren et al., ).

These datasets comprise of positions of interaction of RNA binding protein and RNA target sites derived after PAR-CLIP (Photoactivatable Ribonucleoside Enhanced Crosslinking and Immunoprecipitation), HITS-CLIP-seq (High Throughput sequencing of RNA isolated by crosslinking immunoprecipitation), RIP-seq (RNA immunoprecipitation), iCLIP (individual nucleotide resolution crosslinking and immunoprecipitation), PAR-iCLIP (Photoactivatable Ribonucleoside Enhanced individual nucleotide resolution crosslinking and immunoprecipitation) and CLASH (cross-linking ligation and sequencing of hybrids) followed by sequencing of the pull-down fraction of RNA. The sequenced RNA is further used to identify exact or probable binding site using various bioinformatic approaches. In case of ClipDB, the peak calling and identification were done using PARalyzer, CIMS, Piranha, and CITS software tools. Hence, we stored each of files derived from all databases in form of peaks as separate files for downstream analysis. Details of all the techniques and methodologies used to process the data used in our analysis is given in Table 1.

Table 1

DatabaseTechnique usedSoftware to predict binding sitesTotal datasetsNo of unique RBP
starBaseCLASH, HITS-CLIP, iCLIP, PAR-CLIPAgo PAR-CLIP raw data were reanalyzed using PARalyzer v1.1, other CLIP-identified binding sites clusters/peaks were used directly from their respective publications8638
CLIPdbHITS-CLIP, iCLIP, PAR-CLIP, PAR-iCLIPPARalyzer, CIMS, Piranha, CITS43059
doRiNAiCLIP, HITS-CLIP,CLIP-Seq, PAR-CLIPBinding sites directly adopted from their respective publications7937

The summary of total datasets examined in our study for analysis.

All these RNA binding sites were liftover to hg38/GRCh38 assembly using the CrossMap-0.2.2 tool (Zhao et al., ). The peak information was available for proteins as shown in Supplementary Tables 1A,B. In total, we considered 7 datasets for our study, namely: (1) starBase; (2) doRiNA; (3) Clipdb-PARalyzer; (4) CLIPdb-CIMS; (5) CLIPdb-CITS; (6) CLIPdb-Piranha-stranded); and (7) CLIPdb-Piranha-non-stranded.

Mapping of RNA binding protein interaction sites

The peaks of the RNA binding protein interaction sites were mapped to the lncRNA exons using bespoken perl script and BEDtools (v2.17.0) (Quinlan and Hall, ). The most probable site of interaction (or the peaks) between protein and RNA were derived from datasets taken from doRiNA, starBase, and CLIPdb databases which were processed through standard computational pipelines (as listed in Table 2), offering an easy comparability at the analysis point of view. Further, we tried to analyze the binding sites in each of the individual lncRNA subclasses as defined by GENCODE annotations (i.e., 3 prime overlapping ncRNA, antisense, bidirectional promoter lncRNA, lincRNA, macro lncRNA, miscRNA, non-coding, processed transcripts, pseudogene, retained intron, sense intronic, sense overlapping, and TEC). Similarly, we also plotted the distribution of the binding sites across the protein-coding exons derived from the GENCODE v24 annotation file.

Table 2

S. no.NameFull namedoRiNAstarBaseClipdbPublications
1AGO1Argonaute RISC catalytic component 1✓✓✓Hafner et al., ; Helwak et al., ; Memczak et al.,
2AGO2Argonaute RISC catalytic component 2✓✓✓Hafner et al., ; Gottwein et al., ; Kishore et al., ; Lipchina et al., ; Haecker et al., ; Riley et al., ; Skalsky et al., ; Karginov and Hannon, ; Memczak et al., ; Xue et al.,
3AGO3Argonaute RISC catalytic component 3✓✓✓Hafner et al.,
4AGO4Argonaute RISC catalytic component 4✓✓✓Hafner et al.,
5ALKBH5Alpha-Ketoglutarate-Dependent Dioxygenase AlkB Homolog 5✓✓✓Baltz et al.,
6ATXN2Ataxin 2✓✓Yokoshi et al.,
7C17ORF85/NCBP3Chromosome 17 Open Reading Frame 85/Nuclear Cap Binding Subunit 3✓✓✓Baltz et al.,
8C22ORF28/RTCBChromosome 22 Open Reading Frame 28 /RNA 2,3-Cyclic Phosphate And 5-OH Ligase✓✓Baltz et al.,
9CAPRIN1Cytoplasmic Activation- And Proliferation-Associated Protein 1✓✓✓Baltz et al.,
10CPSF1Cleavage And Polyadenylation Specific Factor 1✓Martin et al.,
11CPSF2Cleavage And Polyadenylation Specific Factor 2✓Martin et al.,
12CPSF3Cleavage And Polyadenylation Specific Factor 3✓Martin et al.,
13CPSF4Cleavage And Polyadenylation Specific Factor 4✓Martin et al.,
14CPSF6Cleavage And Polyadenylation Specific Factor 6✓Martin et al.,
15CPSF7Cleavage And Polyadenylation Specific Factor 7✓Martin et al.,
16CSTF2Cleavage Stimulation Factor Subunit 2✓Martin et al., ; Yao et al.,
17CSTF2TCleavage Stimulation Factor Subunit 2, Tau Variant✓Martin et al.,
18DGCR8DiGeorge Syndrome Critical Region 8✓✓✓Macias et al.,
19EIF4A3Eukaryotic Translation Initiation Factor 4A3✓✓✓Saulière et al.,
20ELAVL1/HUREmbryonic Lethal Abnormal Vision Drosophila)-Like binding protein 1/Human Antigen R✓✓✓Lebedeva et al., ; Mukherjee, ; Friedersdorf and Keene,
21EWSR1Ewing Sarcoma Breakpoint Region 1✓✓✓Hoell et al., ; Paronetto et al.,
22EZH2Enhancer Of Zeste 2 Polycomb Repressive Complex 2 Subunit✓Kaneko et al.,
23FBLFibrillarin✓Kishore et al.,
24FIP1L1Factor Interacting With PAPOLA And CPSF1✓Martin et al.,
25FMR1Fragile X Mental Retardation 1✓✓✓Ascano et al.,
26FOX2RNA binding protein fox-1 homolog 2✓Yeo et al.,
27FUSFused in Sarcoma✓✓✓Hoell et al., ; Lagier-Tourenne et al., ; Nakaya et al., ; Yokoshi et al.,
28FXR1Fragile X Mental Retardation, Autosomal Homolog 1✓✓✓Ascano et al.,
29FXR2Fragile X Mental Retardation Autosomal Homolog 2✓✓✓Ascano et al.,
30HNRNPA1Heterogeneous Nuclear Ribonucleoprotein A1✓Huelga et al.,
31HNRNPA2B1Heterogeneous Nuclear Ribonucleoprotein A2/B1✓Huelga et al.,
32HNRNPCHeterogeneous Nuclear Ribonucleoprotein C (C1/C2)✓✓✓Zarnack et al., ,
33HNRNPDHeterogeneous Nuclear Ribonucleoprotein D✓Yoon et al.,
34HNRNPFHeterogeneous Nuclear Ribonucleoprotein F✓Huelga et al.,
35HNRNPHHeterogeneous Nuclear Ribonucleoprotein H1✓Katz et al.,
36HNRNPLHeterogeneous Nuclear Ribonucleoprotein L✓✓Shankarling et al.,
37HNRNPMHeterogeneous Nuclear Ribonucleoprotein M✓Huelga et al.,
38HNRNPUHeterogeneous Nuclear Ribonucleoprotein U (Scaffold Attachment Factor A)✓Huelga et al., ; Xiao et al.,
39IGF2BP1Insulin Like Growth Factor 2 MRNA Binding Protein 1✓✓✓Hafner et al.,
40IGF2BP2Insulin Like Growth Factor 2 MRNA Binding Protein 2✓✓✓Hafner et al.,
41IGF2BP3Insulin Like Growth Factor 2 MRNA Binding Protein 3✓✓✓Hafner et al.,
42LIN28ALin-28 Homolog A✓✓✓Wilbert et al., ; Hafner et al.,
43LIN28BLin-28 Homolog B✓✓✓Graf et al., ; Hafner et al.,
44METTL3Methyltransferase-Like Protein 3✓✓✓Ping et al.,
45MOV10Moloney Leukemia Virus 10 RISC Complex RNA Helicase✓✓✓Sievers et al.,
46NOP56Nucleolar Protein 5A✓Kishore et al.,
47NOP58Nucleolar Protein 5✓Kishore et al.,
48NUDT21Nudix (Nucleoside Diphosphate Linked Moiety X)-Type Motif 21✓Martin et al.,
49PTBP1Polypyrimidine Tract Binding Protein 1✓✓Xue et al., ; Raj et al.,
50PTBP1 and PTBP2Polypyrimidine Tract Binding Protein 1 and Polypyrimidine Tract Binding Protein 2✓Xue et al.,
51PUM2Pumilio RNA Binding Family Member 2✓✓✓Hafner et al.,
52QKIQuaking homolog KH domain RNA binding✓✓✓Hafner et al.,
53RBM10RNA Binding Motif Protein 10✓Wang et al.,
54RBPMSRNA Binding Protein With Multiple Splicing✓Farazi et al.,
55RTCBRNA 2′,3′-Cyclic Phosphate And 5'-OH Ligase✓Baltz et al.,
56SFRS1Splicing factor arginine/serine-rich 1✓✓Sanford et al.,
57SRRM4Serine/Arginine Repetitive Matrix 4✓Raj et al.,
58TAF15TATA-Box Binding Protein Associated Factor 15✓✓✓Hoell et al., ; Ibrahim et al.,
59TARDBP/TDP-43Transactive response DNA binding protein 43 kDa✓✓✓Tollervey et al., ; Yokoshi et al.,
60TIA1T-Cell-Restricted Intracellular Antigen-1✓✓✓Ule et al.,
61TIAL1TIA-1-Related Protein✓✓Ule et al.,
62TNRC6ATrinucleotide repeat-containing gene 6A✓✓Hafner et al.,
63TNRC6BTrinucleotide repeat-containing gene 6B✓✓Hafner et al.,
64TNRC6CTrinucleotide repeat-containing gene 6C✓✓Hafner et al.,
65U2AF65U2 Small Nuclear Ribonucleoprotein Auxiliary Factor (65kD)✓Zarnack et al.,
66UPF1Up-Frameshift Suppressor 1 Homolog✓Zünd et al.,
67WDR33WD Repeat Domain 33✓Schönemann et al.,
68WTAPWilms Tumor 1 Associated Protein✓Ping et al.,
69YTHDF2YTH N6-Methyladenosine RNA Binding Protein 2✓Wang et al.,
70ZC3H7BZinc Finger CCCH-Type Containing 7B✓✓✓Baltz et al.,

List of the 70 RNA binding proteins derived from the respective databases.

We further tested the significance of binding frequency for each of the lncRNA biotype when compared to the protein coding transcripts. The normalized frequency of binding was calculated by dividing the unique number of RBP peaks mapped from each dataset by unique number of bases of lncRNA/protein coding/random transcripts per kb. Statistical unpaired t-test was applied using R (version 3.1.3) (R Core Team, ) script, to test if any of the lncRNA biotypes had significantly higher RBP binding frequency as compared to the protein coding transcripts.

Combinatorial patterns for RNA binding protein interaction sites in lncRNAs

We explored the possibility of positional clustering of RNA binding protein interaction sites across the lncRNA and protein coding transcripts. For this, we calculated the co-occurrence binding frequencies for each of the 70 RBPs from the six datasets for each of the lncRNAs and protein coding transcripts in the annotation list. For this analysis we did not consider the CLIPdb-Piranha-non-stranded) dataset due to lack of strand orientation information. Bespoke shell scripts were used to identify RBP sites which co-occurred with each other and were therefore clubbed together.

The coordinates for each RBP peak dataset were intersected separately with both the lncRNA and protein coding exons using BEDtools. These intersecting coordinates were then used to calculate the number of bases which were shared between each of the protein datasets to examine their co-occurrence. The values were further normalized by dividing it with the total number of unique bases of individual RBP datasets which were intersecting with lncRNA and protein coding exons. The mapping percentage in protein coding transcripts provided the baseline for co-occurrence frequency of the binding sites. These co-occurrence frequencies were calculated independently for all the RBP across six datasets.

Positional preference of RNA binding protein interaction sites in lncRNAs

We also examined the positional preference of the RNA binding protein interaction sites across the length of lncRNA transcript. As the length of the transcripts varied considerably in our analysis therefore, we briefly define the length of the transcripts as divided into three equal parts. The length of long non-coding transcripts were normalized to 100 nucleotides and arbitrary divided into three equal parts viz., 5 prime end, the middle region, and 3 prime end for comparisons. The notation 5 prime, middle region, and 3 prime denote the positions of the three equal fragments and have no bearing with 5 prime and 3 prime UTRs. Except for datasets analyzed using Piranha, which did not have strand information of the called RBP peaks, all other datasets were used to check for their positional preference. The unique number of bases intersecting with each of the three lncRNA segments was calculated for each dataset. These were further normalized by dividing these values with the unique number of bases in the respective lncRNA segment. Percentage preference was calculated for each segment and the positional location of RNA protein-binding sites were enumerated and plotted as heatmaps.

Additionally, we also plotted the counts of the RNA binding protein interactions sites in protein coding transcripts derived from GENCODE annotation file and the mappings were divided into 3 regions: 5 prime UTR, coding exons, and 3 prime UTR of the coding genes. The CLIPdb-Piranha-non-stranded dataset were not used for the analysis due to the lack of strand information of the peaks.

Results

Analysis of mapping of RNA binding proteins datasets

We analyzed publicly available datasets for 70 RNA binding proteins derived from seven datasets encompassing five technologies viz. PAR-CLIP, HITS-CLIP, iCLIP, RIP-seq, and CLASH. The experimental datasets were downloaded for RNA binding proteins from three databases (details in Table 1). The experiments briefly included high-throughput genome-scale analysis of RNA protein interactions through pull down and sequencing. The derived data in form of interaction sites (or peaks) which were pre-processed using different computational pipelines including PARalyzer, CIMS, Piranha, and CITS for each of the proteins and were mapped onto the hg38 build of the Human reference genome. The total number of peaks mapping to the genome for respective datasets corresponding to each RNA binding protein has been detailed in Supplementary Tables 1A,B. Each of the dataset was kept as a separate file even if the name of the RNA binding protein was same. This was followed to maintain the identity of each dataset as there were differences in number of peaks for same proteins across different databases which could be attributed to the different experimental protocols used for processing including difference in cell lines, conditions or end points, or downstream computational processing. As same protein was present in more than one dataset, we did not group them as one because different databases had differences in the number and position of peaks owing to the differences in the peak calling softwares and computational pipelines adopted by the users. Nevertheless, the differences in the global frequencies have not been influenced by these.

Comparison of RNA binding protein interaction sites within lncRNAs and protein coding genes

We compared the interaction sites for each of the RNA binding proteins in lncRNAs as well as protein-coding transcripts. Toward this end, we used the transcript annotations as provided by GENCODE V24 (Harrow et al., ) for protein-coding transcripts and lncRNAs. In total the dataset comprised of 79,930 protein-coding transcripts from 19,655 genes and 83,215 lncRNA transcripts arising out of 32,446 genic loci. We analyzed the distribution of RNA binding protein interaction sites across lncRNAs and protein coding transcripts.

All proteins showed distinct frequency distribution across both protein-coding and long non-coding transcripts. In general, RBP binding was higher in protein coding transcripts when compared to long non-coding transcripts. But when we looked closely, few of RBPs showed higher enrichment for lncRNA subclass when compared to protein coding transcripts. We tested the significance of the enrichment of RBP sites across lncRNA subtypes as opposed to protein coding transcripts using paired t-test. We observed that six of the biotypes including antisense, lincRNA, miscRNA, processed transcripts, retained intron, and sense intronic were more enriched (p-value ≤ 0.05) for RBP sites as opposed to protein coding transcripts in some or the other RBP dataset.

We plotted the binding frequencies of RBPs in lncRNAs and protein coding transcripts for each of the seven datasets as separate graphs. Those datasets and biotypes which had a significantly higher binding for RBPs have been plotted (Figure 1, Supplementary Figures 1, 2). The RBP binding frequency for CLIPdb-CIMS dataset was significantly higher in lincRNA class when compared to protein coding transcripts for all proteins, while HNRNP (F, H, and U) protein had consistent enrichment for miscRNA class (Figure 1). HNRNP complexes help in processing of pre-mRNAs into functional, translatable mRNAs in the cytoplasm. AGO group from CLIPdb-Piranha-non-stranded dataset were mostly enriched for miscRNA, sense intronic, and lincRNA class compared to protein coding transcript while most of proteins showed enrichment for miscRNA and lincRNA classes (Supplementary Figure 1). In Supplementary Figure 1B, we observed miscRNA and lincRNA class to be mostly enriched for most of proteins including AGO proteins, CSTF2 in sense intronic and DGCR8 in retained intron class. AGO2 protein is an important part of RNA-induced silencing complex (RISC) and is required for RNA-mediated gene silencing (RNAi). CSTF2 plays role in polyadenylation and 3'-end cleavage of mammalian pre-mRNAs. DGCR8 is a component of the microprocessor complex that acts as a RNA- and heme-binding protein that is involved in the initial step of microRNA (miRNA) biogenesis. For the starBase, CLIPdb-CITS, doRiNA, Clipdb-PARalyzer datasets RBPs showed higher frequency distribution for lncRNAs (miscRNA, retained intron processed transcript) compared to protein coding transcripts (Supplementary Figures 2A–D), ATXN2 protein from Supplementary Figure 2D had a comparable binding frequency in miscRNA class to protein coding transcripts. This protein is involved in EGFR trafficking, acting as negative regulator of endocytic EGFR internalization at the plasma membrane. Proteins from CLIPdb-Piranha-stranded had enrichment for miscRNA class when compared to protein coding transcripts (Supplementary Figure 2E).

Figure 1

We additionally chose a random set of 1,000,000 (1 million) genomic loci as a control set with an average length of 240 bases and mapped the RBP sites across this control set. The frequencies of protein binding sites across these random genomic loci, lncRNA, and protein coding transcripts of randomly chosen RBPs from each of the six datasets have been depicted in the Supplementary Figure 3, to illustrate that the frequency of protein binding sites in lncRNAs is not an arbitrary event. The observed RBP frequency was significantly lower for these random positions when compared to protein coding transcripts and lncRNAs. This clearly substantiates the fact that the observed RBP distribution frequencies are not just due to randomness but are inherently due to the class of RNA they bind.

Combinatorial patterns for protein-binding sites in lncRNAs show similar proteins have overlapping binding sites

The seven datasets considered in this study were observed to map onto lncRNA transcripts as well as protein-coding transcripts. To understand whether they map to common subset of loci in the respective transcripts, we evaluated the positional overlaps of the binding sites for each protein from these seven datasets individually. The counts of overlaps were measured as proportion of the total number of independent occurrences of binding sites for each protein. The overlaps were counted separately for all positions in the protein coding transcripts and in lncRNAs. The mapping in protein coding transcripts served as the control set which provided a fair idea of the general overlap in the genomic scale.

Four proteins from the CLIPdb-CITS dataset CSTF2, HNRNPC, TARDB, and TIA1 showed maximum co-occurrence with their respective set of proteins both in protein coding and lncRNAs transcripts while CSTF2, HNRNPC, and TIAL1 co-occurred with each other as well. Our analysis revealed that similar functioning proteins have significantly higher overlapping binding sites with each other, as expected, while EZH2 was an exception in this dataset (Figure 2).

Figure 2

Similarly, RBPs from other five datasets also showed same behavior of co-occurrence between the same set of proteins as shown in Supplementary Figures 4–8 as heatmap. ELAVL1 co-occurred with HUR proteins from doRiNA dataset with high co-occurrence binding frequency as both being the alternate name of same protein. HNRNPF co-occurred with HNRNPU; both are part of the same HNRNP complex, infact all the HNRNP proteins are related to each other.

While protein having similar function such as AGO and DGCR8 proteins were co-occurring in both the doRiNA and CLIPdb-CIMS datasets. Similarly, TNRC6 (A-C) proteins co-occurred with AGO proteins from CLIPdb-Piranha-stranded, Clipdb-PARalyzer, and starBase datasets, from previous observations it is has been seen that functionally related proteins co-occur as in case of TNRC6 with Argonautes, as they have shown to be to play important roles in microRNA mediated regulation of transcripts (Baillat and Shiekhattar, ; Chen et al., ). ATXN2 and TARDB from Clipdb-PARalyzer are known to associate in one complex depending on RNA where they bind, we observed them to co-occur in our analysis (Elden et al., ). From Clipdb-PARalyzer dataset CSTF2 co-occurred with CPSF proteins. Argonaute protein was observed to co-occur with FUS, HNRNP, PTBP1, and PTBP2 from CLIPdb-CIMS datasets and from literature it has been reported that all these proteins interact with each except AGO, hence we believe if other proteins co-occur then AGO should also functionally correlate with these proteins. From starBase dataset, we also observed TAF15 and FUS co-occurred. In addition, we also observed that FUS and TARDB proteins co-occurred from Clipdb-PARalyzer dataset and AGO group of proteins from CLIPdb-CIMS dataset co-occured with HNRNP2B1, HNRNPF, HNRNPM, and HNRNPU proteins. There were other proteins also which co-occurred but with low co-occurrence binding frequency. There was no stark difference in the overlaps of the binding sites between protein coding transcripts and lncRNA sites for each of the proteins considered in our analysis.

Positional clustering of the protein-binding sites

Positional preferences of the RNA binding protein interaction sites were examined across the entire length lncRNAs. The entire length of transcript was calculated by summing up the lengths of individual exons falling in a transcript and then calculating the position of the mapped RNA binding protein interaction site across this calculated length. As the length of the transcript varied therefore, the entire length was arbitrarily divided into three equal parts viz. 5 prime end, middle region, and 3 prime end. Our analysis revealed that the number of RNA binding protein interaction sites for most of the proteins were in majorly mapping to the 3 prime end and the mid segment of the transcripts as shown in Figure 3 and Supplementary Figure 9. To observe the frequencies of binding sites in protein coding transcripts, we mapped and analyzed the RNA binding protein interaction sites in the protein coding transcripts. The binding frequencies for RBPs were evaluated in protein coding transcripts which were divided as 5 prime UTR, CDS, and 3 prime UTR. The data for the same was derived from GENCODE annotation file in form of bed files. We observed that RNA binding protein interaction sites were distributed in 3 prime UTR, 5 prime UTR, and coding exons and frequencies varied for each protein. The HUR/ELAV1 protein showed a positional preference toward the 3 prime end across the lncRNA transcript and the same has been reported recently by Wang and group (Wang et al., ) (Figure 4, Supplementary Figures 10, 11).

Figure 3

Figure 4

We further observed that AGO proteins across the three datasets, namely; Clipdb-PARalyzer, starBase, and doRiNA showed to have a positional preference in protein coding and lncRNA transcripts (Figures 3, 4). When we examined the mapping for the three datasets in protein coding transcripts, we observed that AGO protein showed preference toward the 3 prime UTR. Previous reports have shown AGO proteins bound to miRNAs to target toward 3 prime end of mRNA thereby affecting its translation (Pillai et al., ). Such positional preference for AGO proteins is an established fact when targeting the 3′ end of mRNAs leading to post-transcriptional silencing. We observed similar positional preference for AGO protein in lncRNAs, thereby suggesting certain regulatory roles.

High frequencies of RNA binding protein interaction sites in a subset of transcripts

We also observed that many well-known lncRNAs including XIST, NEAT1, OIP5-AS1, and MALAT1 had large number of RNA binding protein sites across their length. A subset of well-annotated lncRNA genes had consistently large number of binding sites for majority of the proteins considered. MALAT1 (metastasis associated lung adenocarcinoma transcript 1), a well-studied lncRNA with intricate roles in the pathophysiology of cancer Metastases is one of such candidate (Gutschner et al., ). MALAT1 is highly conserved amongst mammals and is known to be localized in nucleus. We plotted the binding sites for all RBPs to the full-length of MALAT1 transcripts and the same is shown in Figure 5 for ClipDB-CIMS, CLIPdb-CITS, and CLIPdb-Piranha-stranded datasets. We combined all the datasets for each protein within a database and divided them into three classes (Cytoplasmic, Nuclear, or Both) based on their cellular localization. The distribution profiles for all the RBPs across the MALAT1 gene was derived using UCSC Genome Browser (Meyer et al., ).

Figure 5

We observed that the RBPs known to be localizing in nucleus were shown to have higher binding sites across MALAT1 when compared to other RBPs. The functional interaction of MALAT1 with a number of RNA binding proteins have been previously studied (Tripathi et al., ), suggesting extensive functional link to the interactions and thereby providing interesting insights for lncRNA functions and biological regulatory networks they take part in. The mapping for all other datasets across the MALAT1 lncRNA is shown in Supplementary Figures 12, 13.

Discussion

LncRNAs have lately emerged as one of the major transcript forms encoded by the human genome, the numbers growing as much as the number of protein-coding transcripts over the years. GENCODE v24 has 83,215 lncRNA loci compared to 79,930 protein-coding gene loci. The functional role of many candidate lncRNAs have been extensively studied in the recent past, nevertheless the general lack of conservation of lncRNAs, even between closely related organisms, barring a handful of candidate lncRNAs has restricted the possibility to model functionalities of lncRNAs in model systems.

The availability of genome-scale assays for evaluating protein-binding sites in RNA (Kishore et al., ), has offered new opportunities to address this issue at much higher confidence and resolution than which were provided by computational approaches (Bellucci et al., ; Puton et al., ). To date, seven datasets for genome-scale protein-RNA interactions are available in public domain (i.e., doRiNA, Clipdb, starBase) and the present analysis makes use of all these available datasets. We show such approaches involving repurposing of datasets could provide immense insights into the biological functions with potential regulation of lncRNAs.

In the present study, we have used the peak information (or the most probable site of interaction between protein and RNA) from seven datasets processed through standardized computational pipeline for accurate assessment of protein-RNA interaction sites (doRiNA, Clipdb, starBase). This allowed us to compare the frequencies of the protein binding sites in systematic fashion. It has not escaped our attention that the datasets encompass a diverse set of experiments; cell line, and experimental protocols, nevertheless; our findings hold true despite these differences available in public domain as part of this analysis encompassing six experimental databases of RNA binding proteins. For instance, one of the most studied RBP, the Argonaute datasets showed similar trends regardless of the diverse experimental protocols (HITS-CLIP, iCLIP, PAR-CLIP) and analysis methodologies employed.

The RBPs considered in our study are known to be involved in varied types of functional roles including silencing, splicing, stability, mRNA processing, and transport. In the current study, we observed RBPs enriched for specific lncRNA biotypes are involved in diverse functions, suggesting their probable functional mechanism of action. RBPs such as AGO, DGCR8, EWSR1, TNRC6A/B/C, and FUS, involved in maintenance of the stability of RNA, were having significant enrichment for the lincRNA, miscRNA, retained intron subclasses suggesting they might be acting as either transporters or as sponges for these RBPs. Another set of RBPs such as CPSF complex, FBL, TAF15, and HNRNP family, playing a role in mRNA processing were shown to be enriched in lncRNA subclasses, signifying that lncRNAs inturn might be acting as guides. These proteins might be also involved in mechanism of lncRNA biogenesis. Enrichment was also observed for proteins such ATXN2, C17ORF85, and HNRNPs which predominantly are involved in the export and transporting of RNA moieties, in addition to proteins such as EIF4A3, FOX2, PTBP1, QKI, SFRS1, SRRM4 among others which are predominantly involved in splicing. Hence our analysis suggests that interaction of lncRNA with such types of RBPs surely provide hints about the possible functional roles lncRNAs might be playing which can be validated by experimental approaches.

We also highlight the localization of lncRNAs and RBPs within a cell. We classified the RBPs based on their known localization within the cells and overlapped it with MALAT1, which is an established nuclear enriched lncRNA. The results indicated that the intensity of nuclear localized RBPs were higher for MALAT1 across all the seven datasets. This further strengthened the fact that these bindings were not an arbitrary event and are indeed interacting with the co-localized lncRNAs.

The present analysis reveals a set of interesting characteristics of protein-RNA interaction in the context of lncRNAs: (1) high frequency of RNA-protein interaction sites in lncRNAs subclasses; (2) co-occurrence of RNA binding protein interaction sites; and (3) positional preference for the binding sites across the transcript length. This analysis, to our best of knowledge is the most comprehensive analysis of RNA binding protein interaction sites in lncRNAs, and provides the basis for further analysis on the functional consequences of these patterns. It has also not escaped our attention that targeting protein-interaction sites and thus the functionalities could be in the future therapeutically explored. Recent reports from other laboratories have explored the possibility of targeting RNA structures using small molecules (Jamal et al., ; Bose et al., ). Further availability of genome-scale protein-RNA interaction datasets and availability of tools to query RNA secondary structures at genome scale (Hofacker, ) would provide us with immense opportunities toward understanding the entire repertoire of functional RNA interactions and phenotypic correlates at a genome-scale level. This would also form the much-needed resource of knowledge to potentially query and understand consequences of genomic variations at these loci.

Conclusion

The interactions between proteins and RNA molecules can provide the essential insights into the functioning of the lncRNAs. In this study, we highlight the enrichment of RBP sites across some of the lncRNA transcript classes in comparison with protein coding transcripts. We have systematically demonstrated that proteins having similar functional roles showed a higher co-occurrence across both lncRNA and protein coding transcripts. Also, the positional preference of most of RBPs agreed with their possible functional roles. Our study gives a compendium of lncRNA and RBP interactions suggesting a large number of functional roles which they can play including silencing, splicing, mRNA processing, export or transport.

Statements

Author contributions

VS conceptualized the analysis. Data analysis was performed by SJ and SG. SJ prepared the data summaries and visualization. SJ and SG wrote the manuscript. All authors reviewed the manuscript.

Funding

This work was funded by the Council of Scientific and Industrial Research (CSIR), India through Grant GENCODE-C (BSC0123).

Acknowledgments

Authors would like to acknowledge Dr. S. Ramachandran and Dr. Sheetal Gandotra for their valuable discussions which helped in compiling the analysis of this study and writing of the manuscript.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer NA and handling Editor declared their shared affiliation.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2018.00027/full#supplementary-material

Supplementary Figure 1

(A,B) Distribution of RNA binding proteins from CLIPdb-Piranha-stranded across six biotypes of lncRNA genes and protein-coding genes. X-axis of the graph shows the distribution of RNA binding protein interaction sites in subclasses of lncRNAs and protein coding genes frequency of binding sites.

Supplementary Figure 2

Distribution of RNA binding proteins from (A) starBase, (B) CLIPdb-CITS, (C) doRiNA, (D) Clipdb-PARalyzer, and (E) CLIPdb-Piranha-non-stranded across 6 biotypes of lncRNA genes and protein-coding genes. X-axis of the graph shows the distribution of RNA binding protein interaction sites in subclasses of lncRNAs and protein coding genes frequency of binding sites.

Supplementary Figure 3

Distribution of RNA binding proteins CLIPdb-CIMS (CSTF2), starBase (AGO1), CLIPdb-CITS (HNRNPC), doRiNA (AGO2), Clipdb-PARalyzer (AGO2), and CLIPdb-Piranha-non-stranded (AGO2) across lncRNA, Protein Coding Transcript and Random Genomic Loci. X-axis of the graph represents random RBPs selected from each dataset and Y-axis depicts the normalized frequency of RNA binding protein interaction sites. The frequency is calculated as the number of unique RBP peaks per unique number of exonic bases per kilobase mapped.

Supplementary Figure 4

The Heatmap depicts the combinatorial patterns of clustered protein-binding sites across lncRNAs (blue in color) and protein coding transcripts (red in color) for doRiNA dataset RBPs. The scale here signifies the number of overlapping binding sites per total number of occurrences for the independent proteins.

Supplementary Figure 5

The Heatmap depicts the combinatorial patterns of clustered protein-binding sites across lncRNAs (blue in color) and protein coding transcripts (red in color) for CLIPdb-CIMS dataset RBPs. The scale here signifies the number of overlapping binding sites per total number of occurrences for the independent proteins.

Supplementary Figure 6

The Heatmap depicts the combinatorial patterns of clustered protein-binding sites across lncRNAs (blue in color) and protein coding transcripts (red in color) for CLIPdb-Piranha-non-stranded dataset) RBP. The scale here signifies the number of overlapping binding sites per total number of occurrences for the independent proteins.

Supplementary Figure 7

The Heatmap depicts the combinatorial patterns of clustered protein-binding sites across lncRNAs (blue in color) and protein coding transcripts (red in color) for starBase dataset RBPs. The scale here signifies the number of overlapping binding sites per total number of occurrences for the independent proteins.

Supplementary Figure 8

The Heatmap depicts the combinatorial patterns of clustered protein-binding sites across lncRNAs (blue in color) and protein coding transcripts (red in color) for Clipdb-PARalyzer dataset) RBPs. The scale here signifies the number of overlapping binding sites per total number of occurrences for the independent proteins.

Supplementary Figure 9

Positional preference of protein-binding sites in lncRNAs transcripts for (A) Clipdb-PARalyzer, (B) CLIPdb-CIMS, (C) starBase, (D) doRiNA, (E) CLIPdb-CITS, and (F) CLIPdb-Piranha-stranded.

Supplementary Figure 10

(A) Distribution of RNA binding proteins sites from (A) Clipdb-PARalyzer, (B) CLIPdb-CIMS, and (C) starBase datasets across Refseq genes. X-axis of the graph depicts the distribution of RNA binding protein interaction sites in refseq genes and Y-axis is the frequency of binding sites.

Supplementary Figure 11

(B) Distribution of RNA binding proteins sites from (A) doRiNA, (B) CLIPdb-CITS and (C) CLIPdb-Piranha-stranded datasets across Refseq genes. X-axis of the graph depicts the distribution of RNA binding protein interaction sites in refseq genes and Y-axis is the frequency of binding sites.

Supplementary Figure 12

Depiction of the mapping of RNA binding protein interaction sites from Clipdb-PARalyzer datasets across the length of MALAT1 lncRNA. The RBP highlighted in gray box are the ones generally localized to cytoplasm (C). The RBP generally localized to nucleus (N) are marked as yellow box. C/N labeled RBPs is the ones which are present in both Nucleus and Cytoplasm.

Supplementary Figure 13

Depiction of the mapping of RNA binding protein interaction sites from ClipDB (doRiNA and starBase datasets) across the length of MALAT1 lncRNA. The RBP highlighted in gray box are the ones generally localized to cytoplasm (C). The RBP generally localized to nucleus (N) are marked as yellow box. C/N labeled RBPs is the ones which are present in both Nucleus and Cytoplasm.

Supplementary Tables 1

(A) Detailed list of publically available datasets derieved from Starbase and Dorina databases. (B) Detailed list of publically available datasets derieved from CLIPdb database.

    Abbreviations

  • lncRNA

    long non-coding RNA

  • DNA

    Deoxyribonucleic acid

  • CLIP

    cross-linking immunoprecipitation

  • HITS-CLIP

    UV cross-linking and immunoprecipitation with high-throughput sequencing

  • PAR-CLIP

    Photoactivatable Ribonucleoside-Enhanced Crosslinking and Immunoprecipitation

  • iCLIP

    individual-nucleotide resolution Cross-Linking and Immuno Precipitation

  • ncRNA

    Non-coding RNA

  • lincRNA

    long intergenic RNA

  • TEC

    To be Experimentally Confirmed

  • RIP-seq

    RNA Immunoprecipitation sequencing

  • CLASH

    cross-linking ligation and sequencing of hybrids

  • CIMS

    Crosslinking induced mutation site

  • CITS

    crosslinking induced truncation analysis

  • RISC

    RNA-induced silencing complex

  • RNAi

    RNA interference

  • miRNA

    microRNA

  • RBP

    RNA binding protein

  • UTR

    untranslated region

  • CDS

    coding sequence

  • UCSC

    University of California, Santa Cruz.

References

Summary

Keywords

long non-coding RNAs, RNA binding proteins, protein-lncRNA interactions, Argonaute (ago), MALAT1

Citation

Jalali S, Gandhi S and Scaria V (2018) Distinct and Modular Organization of Protein Interacting Sites in Long Non-coding RNAs. Front. Mol. Biosci. 5:27. doi: 10.3389/fmolb.2018.00027

Received

27 September 2017

Accepted

14 March 2018

Published

04 April 2018

Volume

5 - 2018

Edited by

Naoyuki Kataoka, The University of Tokyo, Japan

Reviewed by

Keith W. Vance, University of Bath, United Kingdom; Nobuyoshi Akimitsu, The University of Tokyo, Japan

Updates

Copyright

© 2018 Jalali, Gandhi and Scaria.

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Saakshi Jalali Vinod Scaria

This article was submitted to RNA, a section of the journal Frontiers in Molecular Biosciences

Disclaimer

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.

Outline

Figures

Cite article

Copy to clipboard


Export citation file


Share article

Share on WeChat

Scan with WeChat to share this article

Article metrics