- Software
- Open access
- Published:
RADAR: annotation and prioritization of variants in the post-transcriptional regulome of RNA-binding proteins
Genome Biology volume 21, Article number: 151 (2020)
Abstract
RNA-binding proteins (RBPs) play key roles in post-transcriptional regulation and disease. Their binding sites cover more of the genome than coding exons; nevertheless, most noncoding variant prioritization methods only focus on transcriptional regulation. Here, we integrate the portfolio of ENCODE-RBP experiments to develop RADAR, a variant-scoring framework. RADAR uses conservation, RNA structure, network centrality, and motifs to provide an overall impact score. Then, it further incorporates tissue-specific inputs to highlight disease-specific variants. Our results demonstrate RADAR can successfully pinpoint variants, both somatic and germline, associated with RBP-function dysregulation, which cannot be found by most current prioritization methods, for example, variants affecting splicing.
Background
Dysregulation of gene expression is a hallmark of many diseases, including cancer [1]. In recent years, the accumulation of transcription-level functional characterization data, such as transcriptional factor binding, chromatin accessibility, histone modification, and methylation, has brought great success to annotating and pinpointing deleterious variants. However, beyond transcriptional processing, genes also experience various delicately controlled steps, including the conversion of premature RNA to mature RNA, and then the transportation, translation, and degradation of RNA in the cell. Dysregulation in any one of these steps can alter the final fate of gene products and result in abnormal phenotypes [2,3,4]. Furthermore, the post-transcriptional regulome covers an even larger amount of the genome than coding exons and demonstrates significantly higher cross-population and cross-species conservation. Unfortunately, variant impact in the post-transcriptional regulome has been barely investigated, partially due to the lack of large-scale functional mapping.
RNA-binding proteins (RBPs) have been reported to play essential roles in both co- and post-transcriptional regulation [5,6,7]. RBPs bind to thousands of genes in the cell through multiple processes, including splicing, cleavage and polyadenylation, editing, localization, stability, and translation [8,9,10,11,12]. Recently, scientists have made efforts to complete these post- or co-transcriptional regulome by synthesizing public RBP binding profiles [13,14,15,16], which have greatly expanded our understanding of RBP regulation. Since 2016, the Encyclopedia of DNA Elements (ENCODE) Consortium started to release data from various types of assays on matched cell types to map the functional elements in post-transcriptional regulome. For instance, ENCODE has released large-scale enhanced crosslinking and immunoprecipitation (eCLIP) experiments for hundreds of RBPs [17]. This methodology provides high-quality RBP binding profiles with strict quality control and uniform peak calling to accurately catalog the RBP binding sites at a single-nucleotide resolution. Simultaneously, ENCODE performed expression quantification by RNA-Seq after knocking down various RBPs. Finally, ENCODE has quantitatively assessed the context and structural binding specificity of many RBPs by Bind-n-Seq experiments [18].
In this study, we aimed to construct a comprehensive RBP regulome and a scoring framework to annotate and prioritize variants within it. We collected the full catalog of 318 eCLIP (for 112 RBPs), 76 Bind-n-Seq, and 472 RNA-Seq experiments after RBP knockdown from ENCODE to construct a comprehensive post-transcriptional regulome. By combining polymorphism data from large sequencing cohorts, like the 1000 Genomes Project, we demonstrated that the RBP binding sites showed increased cross-population conservations in both coding and noncoding regions. This strongly indicates the purifying selection on the RBP regulome. Furthermore, we developed a scoring scheme, named RADAR (RNA BinDing Protein regulome Annotation and pRioritization), to investigate variant impact in such regions. RADAR first combines RBP binding, cross-species and cross-population conservation, network, and motif features with polymorphism data to quantify variant impact described by a universal score. Then, it allows tissue- or disease-specific inputs, such as patient expression, somatic mutation profiles, and gene rank list, to further highlight relevant variants (Fig. 1). By applying RADAR to both somatic and germline variants from disease genomes, we demonstrate that it can pinpoint disease-associated variants missed by other methods. In summary, RADAR provides an effective approach to analyze genetic variants in the RBP regulome and can be leveraged to expand our understanding of post-transcriptional regulation. To this end, we have implemented the RADAR annotation and prioritization scheme into a software package for community use (radar.gersteinlab.org).
Results
Defining the RBP regulome using eCLIP data
We used the binding profiles of 112 distinct RBPs from ENCODE to fully explore the human RBP regulome (Additional file 1: Table S1), which has been previously under-investigated. Many of these RBPs are known to play key roles in post-transcriptional regulation, including splicing, RNA localization, transportation, decay, and translation (Additional file 2: Figure S1 and Table S1).
Our definition of the RBP regulome covers 52.6 Mbp of the human genome after duplicate and blacklist removal (Fig. 2a). It is 1.5 and 5.9 times the size of the whole exome and lincRNAs, respectively. In addition, only 53.1% of the RBP regulome has transcription-level annotations, such as transcription factor binding, open chromatin, and active enhancers. 55.1% of the RBP regulome is in the immediate neighborhood of the exome regions, such as coding exons, 3′ or 5′ untranslated regions (UTRs), and nearby introns (Fig. 2c; see the “Methods” section and Additional file 2: Table S2 for more details). Furthermore, we observed a significantly higher cross-species conservation score in the peak regions versus the non-peak regions in almost all annotation categories, providing additional evidence of the regulatory roles of RBPs (Fig. 2c). In summary, the large size of the regulome, the limited overlap with existing annotations, and the elevated conservation level highlight the necessity of computational efforts to annotate and prioritize the RBP regulome.
Using universal features for RADAR score
To annotate and prioritize variants in RBP binding sites, we built a universal score framework for RADAR that includes three components: (1) sequence and structure conservation, (2) network centrality, and (3) nucleotide impact from motif analysis.
Sequence and structure conservation in the RBP regulome
Cross-species sequence comparisons have been widely used to discover regions with biological functions [19, 20]. For example, GERP score maps the human genome to other species to identify nucleotide-level evolutional constraints [21, 22]. Therefore, we used the GERP score in our RADAR universal framework to detect potentially deleterious mutations in the RBP regulome (see the “Methods” section for more details).
Since the enrichment of rare variants indicates a purifying selection in functional regions in the human genome [19, 23, 24], here, we inferred the conservation of RBP binding sites by integrating population-level polymorphism data from large cohorts (i.e., the 1000 Genomes Project) [25, 26]. GC percentage may confound such inference by introducing read coverage variations, which is a sensitive parameter in the downstream variant calling process [27, 28]. Therefore, we calculated the fraction of rare variants, defined as those with derived allele frequencies (DAFs) less than 0.5%, within the binding sites of each RBP. Then, we compared them with those from regions with similar GC content as a background (see the “Methods” section for more details). In total, 88.4% of the RBPs (99 out of 112) showed elevated rare variant fraction in coding regions after GC correction (Fig. 3A). Similarly, in the noncoding part of the binding sites, 93.8% of RBPs (105 out of 112) exhibited an enrichment of rare variants. This observation convincingly demonstrates the accuracy of our RBP regulome definition (Additional file 1: Table S2).
Some well-known disease-causing RBPs demonstrate the largest enrichment of rare variants. For example, the oncogene XRN2, which binds to the 3′ end of transcripts to degrade aberrantly transcribed isoforms, showed significant enrichment of rare variants in its binding sites [29]. Specifically, it demonstrates 12.7% and 10.3% more rare variants in the coding and noncoding regions, respectively (adjusted P values at 1.89 × 10−9 and 2.85 × 10−118 for one-sided binomial tests) [30]. Hence, we used the enrichment of rare variants to infer the selection pressure in RBP binding sites and adjust the universal variant scores in such regulator regions (see the “Methods” section for more details).
RNA secondary structures have been reported to affect almost every step of protein expression and RNA stability [31]. We incorporated structural features predicted by Evofold, which uses a phylogenetic stochastic context-free grammar to identify functional RNAs in the human genome that are deeply conserved across species [32]. We found that the RBP binding sites demonstrated significantly higher conversation after intersecting with conserved structural regions defined by Evofold (Additional file 2: Figure S2). Thus, we used the Evofold regions in our universal scoring system.
Highlighting variants in binding hubs
It has been reported that genes within network hubs demonstrate higher cross-population conservation—a sign of strong purifying selection [23, 24, 33]. We hypothesized that RBP binding hubs could show similar characteristics because once mutated, they might introduce larger regulation alterations. To test this, we separated the regulome based on the number of associated RBPs. Most regulome regions (62%) were associated with only one RBP (Additional file 2: Figure S3). As the number of RBPs increased, we observed a clear trend of larger rare variant enrichment (Fig. 3d). For instance, noncoding regions with at least five or 10 RBPs exhibited 2.2% or 13.4% more rare variants, respectively (top 5% and 1%, Fig. 3d). This observation supports our hypothesis that the RNA regulome hubs are under stronger selection pressure and, therefore, should be given higher priority when evaluating the functional impact of mutations (details in Additional file 2: Figures S4, S5, and S6).
Emphasizing genes differentially expressed after RBP knockdown
RNA-Seq expression profiling before and after shRNA-mediated RBP depletion from ENCODE can help to infer the gene expression changes introduced by RBP knockdown. Variants with disruptive effects on RBP binding may affect or even completely remove the RBP binding and hence affect gene expressions in a similar way. Therefore, we extracted the differentially expressed genes from RNA-Seq before and after shRNA-mediated RBP depletion (Additional file 1: Table S3). Then, we up-weighted all variants that were located near the differentially expressed genes (Additional file 1: Table S4) and simultaneously disrupted the binding of the corresponding RBPs (schematic in Additional file 2: Figure S7).
Using motif analysis to determine nucleotide-level impact
Mutations that change the RBP binding affinity may alter the RBP regulation via motif disruption. We quantified the difference of position weight matrix (PWM) scores of the mutant allele against the reference allele. RADAR consists of two sources of motifs. First, we used the motifs identified from RNA Bind-n-Seq experiments from ENCODE because it has been reported that many RBP binding events in vivo can be captured by binding preferences in vitro. Second, we used the de novo motifs discovered directly from binding peaks using the default settings in DREME (see details in the “Methods” section). For each variant, we quantified the nucleotide effect using the highest motif score from these two sources.
Incorporating user-specific features to re-weight variant impact
Variant prioritization can be improved if informative priors can be appropriately incorporated into the scoring system. Therefore, our RADAR framework allows various types of user inputs to help identify disease-relevant variants. Specifically, we adopted a top-down scheme to incorporate regulator- and element-level information to up-weight factors that are possibly associated with disease of interest.
Highlighting key regulators through expression profiles
Key regulators are often associated with disease progression, so variants that affect such regulation should be prioritized [34]. RADAR finds such key regulators by combining the RBP regulatory network information with expression profiles. Specifically, for cancer, we first constructed the RBP network from eCLIP binding peaks and used the TCGA data to define the gene differential expression status from disease and normal cell types (see the “Methods” section). Then, for each RBP, we quantified its regulation potential by associating its network connectivity with aggregated disease-to-normal differential expressions from many samples using regression. We applied this approach to 19 cancer types from TCGA, and the regulation potentials are given in Fig. 4. The values of the regulation potential (β1, see the “Methods” section) for all cancer types and RBPs are provided in Additional file 1: Table S5. We found that among the RBPs with larger regulation potential, many have been reported as cancer-associated genes (Additional file 1: Table S6). Our regression approach was also performed on a patient level, and survival analysis based on each patient’s regulatory potential was performed (see Fig. 4c). Interestingly, the regulatory potential of two key RBPs PPIL4 and SUB1 was found to be significantly associated with patient survival (Fig. 4c). In our RADAR framework, we further highlight variants that are associated with RBPs with high regulation potential in their corresponding cancer types by adding an extra score to their disease-specific scores (see more details in the “Methods” section). We can easily extend such analysis for other diseases by incorporating differential expression profiles from other cohorts such as GTEx [35, 36].
Up-weighting key elements from either prior knowledge or mutational profiles
RADAR reconsiders the functional impact difference among RBP peaks by their associated genes. For example, genes that undergo significant expression or epigenetic changes are mostly cell type-specific and can be used to highlight more relevant variants. Currently, our RADAR framework can up-weight all the RBP peaks that are close to genes with significant differential expression (DESeq2 [37]).
In addition, RADAR can incorporate somatic variant recurrence, which has been widely used to discover key disease regions, to re-weight different RBP peaks. Peaks with more somatic mutations than expected are often considered to be disease-driving [38,39,40]. Here, we first defined a local background somatic mutation rate from a large cohort of cancer patients to evaluate the mutation burden in each RBP peak (see details in the “Methods” section). Variants that are associated with burdened elements are given higher priority in our scoring scheme.
Prioritizing variants with a RADAR-weighted scoring scheme
By integrating the pre-built and user-specific data context described above, our scoring scheme evaluated the functional impacts of variants that are specific to post-transcriptional regulation (Fig. 1 and Table 1). We used an entropy-based criterion to up-weight rarer annotations. First, RADAR added up the (universal) score of variants for all pre-built features, which include sequence and structure conservation, network binding hub, RBP-gene association, and motif information, using the entropy-based weights. Then, depending on user inputs, RADAR further up-weighted variants’ score based on tissue-specific information from mutations in the key RBP binding sites, nearby genes with differential expression, or the RBP regulatory potential (see Additional file 2: Figure S8 for more details).
Applying RADAR to pathological germline variants
We calculated the universal RADAR scores on all pathological variants from HGMD (version 2015) and compared them with variant scores from 1000 Genomes variants as a background. As expected, the HGMD variants are scored significantly higher (Additional file 2: Figure S9). For example, the mean RADAR score for HGMD variants is 0.589, while it is only 0.025 for 1000 Genomes variants (P value < 2.2 × 10−16 for one-sided Wilcoxon test). We further compared the universal RADAR scores of HGMD and 1000 Genomes variants within the RBP regulome to remove potential bias since HGMD variants may be more likely to be within or nearby to exons. We still observed a significantly higher universal RADAR score in the HGMD ones (1.871 versus 1.337, P value < 2.2 × 10−16 for one-sided Wilcoxon test, Additional file 2: Figure S9).
We further compared RADAR scores of HGMD variants with other methods. In total, 720 HGMD variants were explained by our methods but could not be highlighted by other methods (see details in the “Methods” section, Additional file 1: Table S7). Many of these variants are located nearby the splice junctions. An example is shown in Fig. 5a (more details in Additional file 1: Table S8). This variant is located 4 bp away from the splice junction in BRCA1. eCLIP experiments showed strong binding evidence in 5 RBPs (Fig. 5a). Specifically, the T to C mutation strongly disrupts the binding motif of PRPF8, increasing the possibility of splicing alteration effects. Our finding is not highlighted in previous methods for variant prioritization, such as FunSeq, CADD, and FATHMM-MKL (Additional file 2: Table S3).
Applying RADAR universal score to somatic variants in cancer
We next aimed to leverage our universal RADAR scheme to evaluate the deleteriousness of somatic variants from public datasets. Due to the lack of a gold standard, we evaluated our results from two perspectives. First, we reasoned that since hundreds of cancer-associated genes are known to play essential roles through various pathways [41, 42], variants associated with these genes are likely to have a higher functional impact [23]. To test this hypothesis, we first selected variants within the 1-kb region of the COSMIC [43] genes and compared them with other variants. We tested four cancer types, breast, liver, lung, and prostate cancer, and found in all cases that variants associated with COSMIC genes showed significant enrichment, with a larger RNA-level functional impact (Fig. 6 and Additional file 2: Figure S10). For example, we found a 4.58- and 8.75-fold increase in high-impact variants at a threshold level of 3 and 4, respectively, in breast cancer patients (P < 2.2 × 10−16, one-sided Wilcoxon test).
In our second approach, we hypothesized that variant recurrence could be a sign of functionality and may indicate an association with cancer [19, 23, 24]. Thus, we compared the variants’ scores from RBP binding peaks with or without recurrence. Specifically, we separated the RBP peaks with variants mutated in more than one sample from those that were mutated in only one sample and then compared the universal RADAR scores. We found that in most cancer types, peaks with recurrent variants were associated with a larger fraction of high-impact mutations. For example, in breast cancer, recurrent elements demonstrated a 1.67- and 2.57-fold more high-impact variants with RADAR greater than 3.0 and 4.0, respectively, resulting in a P value of 2.2 × 10−16 in one-sided Wilcoxon test. We observed a similar trend in most of the other cancer types (Additional file 2: Figure S10).
A case study on breast cancer patients using disease-specific features
We applied our method to a set of breast cancer somatic variants from 963 patients released by Alexandrov and Stratton [44]. We used the COSMIC gene list and expression and mutational profiles as additional features. In total, we found that around 3% of the 687,517 variants could alter post-transcriptional regulation to some degree. We incorporated the above disease-specific features and demonstrated how they could help to re-weight the variant scoring process on a coding variant in Fig. 5b. This variant is located within an RBP binding ultra-hot region and showed high sequence conservations (7% more rare variants for its binding RBP). It also demonstrated a strong motif disruption effect (PPIG in Fig. 5b). All such features resulted in a universal RADAR score of 3.67, which is ranked 290 out of all variants. However, we found that it is located in the exon region of the well-known tumor suppressor TP53 (orange track in Fig. 5b), and its binding peaks demonstrated more than expected somatic mutations (purple in Fig. 5). Besides, 3 out of the 6 RBPs binding there showed high regulation potential in breast cancer (green in Fig. 5). Hence, these additional features boost its overall RADAR score to 6.67, which is ranked 47 out of all variants. In comparison, this variant only shows moderate scores for FunSeq2 (3) and CADD (7.46), while it is scored in the top but showed a much lower rank than RADAR.
RADAR aims to prioritize variants relevant to the post-transcriptional regulome, while FunSeq2, FATHMM-MKL, and CADD focus on those that affect the transcriptional regulome. Therefore, we do find many variants that demonstrate a high overall RADAR score but only show moderate FunSeq2, CADD, and FATHMM-MKL scores. For example, 127 coding and 78 noncoding variants that are ranked within the top 1% of overall RADAR scores are not in the top 10% of CADD, FunSeq2, or FATHMM-MKL scores (Additional file 1: Table S9 and T10). Many of such variants are located in RBP binding hubs, undergo strong purifying selection, demonstrate strong motif disruptiveness, and are regulated by key RBPs that are associated with breast cancer from multiple sources of evidence (Additional file 2: Figure S11). We believe the discovery of such events demonstrates the value of RADAR as an important and necessary complement to the existing transcriptional-level function annotation and prioritization tools.
Discussion
In this study, we integrated the full catalog of eCLIP, Bind-n-Seq, and shRNA RNA-Seq experiments from ENCODE to build an RNA regulome for post-transcriptional regulation. Our defined RBP regulome is remarkably larger than one may think. It covers up to 52.6 Mbp of the genome (Fig. 2a), and the majority of it is not covered by previous annotations focusing on transcriptional-level regulation, such as DHS, TFBS, and enhancers. We found that the RBP regulome demonstrated noticeably higher conservation in two aspects: higher cross-species conservation in almost all annotation categories (Fig. 2c) and higher cross-population conservation by showing significant enrichment in rare variants (Fig. 3). These two sources of evidence support the notion that the RBP regulome is under a strong purifying selection and carries out essential biological functions. Also, these results signify the necessity of computational tools to annotate and prioritize variants in the RBP regulome. Furthermore, we showed that the cross-population conservation in eCLIP binding sites demonstrates strong importance when predicting disease versus common germline variants (Additional file 2: Figure S12 and Table S4).
By integrating a variety of regulator-, element-, and nucleotide-level features, we propose an entropy-based scoring frame, RADAR, to investigate the impact of somatic and germline variants. The variant prioritization framework of RADAR contains two parts. First, by incorporating eCLIP, Bind-n-Seq, shRNA RNA-Seq experiments with conservation and structural features, we built a pre-defined data context to quantify the universal variant impact score. This approach is suitable for multiple disease analysis or cases where no other prior information can be used. We applied this RADAR universal score to HGMD pathological variants and highlighted many candidates that cannot be highlighted by other methods. Besides, our RADAR framework provided detailed explanations of the underlying disease-causing mechanism (Fig. 5a). In addition to the universal score, RADAR also allows user-specific inputs such as prior gene knowledge and patient expression and mutation profiles for a re-weighting process to highlight relevant variants in a disease-specific manner. As an example, we applied the RADAR disease-specific scores to variants from several cancer types and showed that RADAR could identify relevant variants in key cancer-associated genes (Fig. 5b). We additionally showed that the RADAR framework trained in a closely matched cell type has better power to pinpoint pathological variants in a particular disease, which is promising as new eCLIP experiments are performed in various other cell types (see Additional file 2: Figures S13 and S14).
It is important to note that as compared to ChIP-Seq experiments which generate peaks with up to kilobase pair resolution, eCLIP experiments provide a higher resolution functional site annotation (even single-nucleotide resolution). Such accurate and compact annotation can greatly improve our variant function interpretations. We also want to mention that most of the current eCLIP peak calling approaches call peaks on the annotated transcribed regions (Additional file 2: Figure S15). With the development of computational approaches for eCLIP peak calling, we hope that the size of our annotated RBP regulome can be further expanded.
Conclusions
In summary, we have shown that RADAR is a useful tool for annotating and prioritizing post-transcriptional regulome for RBPs, which has not been covered by most of the current variant impact interpretation tools. Our method provides additional layers of information to the current gene regulome. Importantly, the RADAR scoring scheme can be used in conjunction with existing transcriptional-level variant impact evaluation tools, such as FunSeq [23, 24], to quantify variant impacts. Given the fast-expanding collection of RBP binding profiles from additional cell types, we envision that our RADAR framework can better tackle the functional consequence of mutations from both somatic and germline genomes.
Methods
eCLIP data processing and quality control
We collected 318 eCLIP experiments of 112 unique RBPs from the ENCODE data portal (encodeprojects.org, released and processed by July 2017). eCLIP data was processed through the ENCODE 3 uniform data processing pipeline, and peaks with score 1000 were used in our analysis. We then removed peaks that overlapped with blacklisted regions. We further separated the peaks into coding regions and the noncoding regions in our analysis to infer the selection pressure. We also provide versions of the eCLIP peaks that are annotated by RBP’s function, such as splicing—which is the most common function aside from RNA binding (see radar.gersteinlab.org).
Universal RADAR score
Cross-population conservation inference
The cross-population conservation score consists of two components. The Shannon entropy considers the length effect of the RBPs while the selection pressure inference aims to determine the conservation of regions. For the Shannon entropy, for each RBP, we define f to be:
where nin represents the number of 1000 Genomes variants falling in that RBP peaks, and ntotal is the total number of 1000 Genomes variants (fixed number). In this way, f considers the binding site coverage of an RBP, since a larger coverage is more likely to have a larger value of nin. The Shannon entropy is therefore equal to:
We then calculate the selection pressure from the enrichment of rare germline variants from the 1000 Genomes Project. Our analysis at each step is separated into coding and noncoding parts. For a given RBP, we suppose its binding peaks contain nr rare variants (DAF ≤ 0.005) and nc common variants. The percentage of rare variants in that RBP’s binding peaks is defined as:
The value of ρ is often confounded by factors such as GC content. To correct for potential GC content bias, we bin the genome into 500 bp bins and group them according to their GC percentage. Then, we compute the background rare variant percentage using the same rare and common variants from the 1000 Genomes Project for each group (see Additional file 2: Figures S16 and S17). For a given RBP with GC percentage g, we select the background group with the closest GC, to obtain a background rare variant percentage \( {\rho}_b^g \). Therefore, after adjusting for GC bias, the enrichment of rare variants is defined to be:
RBPs with a ρadj larger than 1 suggests a higher-than-expected selection pressure. We then adjust the population conservation entropy score as follows:
Given a variant falling in the RBP regulome that intersects a set of RBP eCLIP peaks, set P, the cross-population conservation score of that variant is equal to the maximum Spopulation _ conservation for all RBPs in set P.
Cross-species conservation using GERP
We use the GERP score to measure the cross-species conservation. For each position, a GERP of greater than 2 is often used to define bases that are conserved. The transformation of GERP to a RADAR component score is adapted from Fu et al. Therefore, a sigmoidal transformation is used to fit the GERP scores between 0 and 1, and the parameters are used to force the curve to be sharp at GERP equal to 2 (Additional file 2: Section 13.2).
Structural conservation
We use the output of Evofold as an indicator of cross-species RNA structure. A variant falling in a region given by Evofold as conservative receives a score of 1 while a variant that does not fall in such region receives a score of 0.
RBP binding hubs and networks
We define the number of RBPs binding at a position to as H. We first separated the RBP peaks into the coding and noncoding regions and then grouped the regions on the genome based on H. For each group of regions with hub number H, we calculated the GC-corrected enrichment of rare variants ρadj ∣ H for each group in coding and noncoding regions by Eqs. (3) and (4). We determine the hub numbers, Hnormal, Hhot, and Hultra − hot associated with normal, hot, and ultra-hot regions, respectively. Hhot and Hultra − hot are associated with the top 5% and 1% of binding RBPs (Fig. 3) to represent rare and ultra-rare events, respectively. Our values of ρadj ∣ H are altered in such a way to reflect this phenomenon. The \( {\rho}_{\mathrm{adj}\mid H<{H}_{\mathrm{hot}}} \) associated with hub scores less than that of the hot regions are converted to 0. The \( {\rho}_{\mathrm{adj}\mid {H}_{\mathrm{hot}}\le H<{H}_{\mathrm{ultra}-\mathrm{hot}}} \) associated with hot regions demonstrate a mostly increasing trend and are smoothed using a kernel smoother. The \( {\rho}_{\mathrm{adj}\mid H\ge {H}_{\mathrm{hot}}} \) of the ultra-hot region is kept constant, equal to the max \( {\rho}_{\mathrm{adj}\mid {H}_{\subset \mathrm{hot}}} \) of the hot region.
We also compute f for the set of regions associated with a fixed H (from Eq. (1)). For each value of f, we compute the Shannon entropy from Eq. (2) to be SH. Finally, we multiply the values of ρadj by its respective Shannon entropy, S. A variant falling in the regulome with hub score H would have a score equal to the ρadj ∣ H ∗ SH.
Motif analysis and disruption
We used the changes of PWMs introduced by a variant to quantify the motif disruptiveness effect through motiftools (https://github.com/gersteinlab/MotifTools). Specifically, we defined the disruption score, Dscore, as in Eq. (6) to represent the difference between sequence specificities to an alternative sequence.
where Pref and Palt are the PWM scores from the reference and alternative allele, respectively. Here, the motif scores for reference and alternate sequences are given as:
To quantify a motif-breaking event, we require that the P value for the reference allele is at least 5 × 10−4. There are two motif sources in our analysis. First, we identified RBP motifs using the DREME software (version 4.12.0) directly from RBP peaks [45]. Then, we also incorporated motifs from RNA Bind-n-Seq (RBNS) [18] to characterize sequence and structural specificities of RBPs. For each variant that affected multiple RBP binding profiles, we used the max score (Additional file 2: Table S5). A threshold of Dscore > 3 is used to describe a disruption event that is significant, and a variant having a Dscore less than this threshold receives a score of 0. For variants receiving a Dscore larger than the threshold, we additionally compute the Shannon entropy given for a variant with Dscore = D as:
where f(v, D) represents the number of the 1000 Genomes Project variants, v, that have a Dscore greater than D divided by the total number of 1000 Genomes Project variants. This entropy is then weighted by the population conservation score.
RBP-gene association using shRNA RNA-Seq
To determine if an RBP, R, is associated with a gene, g, we intersect the peaks of R with the transcript annotation of g. If an intersect exists, we form a linkage between the intersected peak of R and g. If some variant falls in that specific peak of R, the variant significantly disrupts the motif of RBP R, and gene g demonstrates at least a 2.5-fold change in its expression after KD of RBP R, we give the variant an additional score of 1.
Tissue-specific score
RBP regulatory potential
RADAR allows inputs in addition to the pre-built context to calculate the disease-specific variant score. In this paper, we used the TCGA expression profiles as an example of the cancer variant prioritization. Specifically, we downloaded expression profiles of 19 cancer patients of 24 types from TCGA. In order to get a robust differential expression analysis, we excluded several cancer types that have less than 10 normal expression profiles and used DESeq2 [37] to find tumor-to-normal differentially expressed genes (corrected P from DESeq2 < 0.05). Let \( {y}_i^k \) represent the differential expression status of gene i of the kth cancer type.
We inferred the regulatory power of each RBP, R, through a regression approach of the above differential expression and RBP network connectivity as:
where \( {\overrightarrow{\boldsymbol{x}}}^R \) is the binary connectivity vector for all genes and R (1 if the gene is a target, otherwise 0). We used the absolute value of β1k, R to indicate the regulation potential of RBP R in cancer type k. If a variant falls in a region with at least one RBP binding, and at least one of the P values associated with β1k, R is significant, then we consider variants falling in that particular RBP to have an additional score of 1.
Recurrence in somatic mutations
We prioritized variants in RBP binding sites are with more-than-expected somatic mutations. To evaluate the somatic mutation burden, we first separated the genome into 1 Mbp bins and calculated a local background mutation rate in each window. Then, for each eCLIP peak, we counted the number of somatic mutations and compared it to the nearest local 1 Mbp context using a one-sided binomial test. If a specific RBP binding site was enriched for somatic mutations, the variant falling in that site was given a higher priority.
Differentially expressed key genes
For each peak of each RBP, we find the associated gene of that peak by intersecting with the GENCODE gene definitions. Using the DESeq2 results, we consider genes with q values that are less than 0.05 to be differentially expressed. If an RBP peak is associated with a gene that is significantly differentially expressed in a tissue type, we increase the score of the variant falling in such peak by 1.
Availability of data and materials
We have made this RNA variant annotation and prioritization tool available as an open-source software at https://github.com/gersteinlab/RADAR [46]. The source code is released under MIT License. The website contains details on usage, examples, resources, and dependencies. We also provided a genome-wide pre-built RADAR baseline score for every base pair on the genome (hg19 and GRCh38 version of the genome). We released all pre-built data context necessary for RADAR and these genome-wide baseline scores at Zenodo which can be found at https://doi.org/10.5281/zenodo.1451838 [47]. Users can also directly query the annotation and functional impact score from radar.gersteinlab.org. See more details in Additional file 2: Section 15.
The eCLIP and shRNA datasets used in RADAR are from the ENCODE Project and available at encodeproject.org [48]. We used germline single-nucleotide polymorphisms from the 1000 Genomes Project, phase 1 [26]. Somatic variants are collected from Alexandrov et al. [49]. The GERP score is from [21, 22]. The RNA-Seq expression data and clinical survival data are from the TCGA Project [50]. Gene annotations are from the GENCODE Project [51, 52].
References
Croce CM. Causes and consequences of microRNA dysregulation in cancer. Nat Rev Genet. 2009;10:704–14. https://doi.org/10.1038/nrg2634.
Romanoski CE, Glass CK, Stunnenberg HG, Wilson L, Almouzni G. Epigenomics: roadmap for regulation. Nature. 2015;518:314–6. https://doi.org/10.1038/518314a.
Yang G, Lu X, Yuan L. LncRNA: a link between RNA and cancer. Biochim Biophys Acta. 2014;1839:1097–109. https://doi.org/10.1016/j.bbagrm.2014.08.012.
Schmitt AM, Chang HY. Gene regulation: long RNAs wire up cancer growth. Nature. 2013;500:536–7. https://doi.org/10.1038/nature12548.
Gerstberger S, Hafner M, Tuschl T. A census of human RNA-binding proteins. Nat Rev Genet. 2014;15:829–45. https://doi.org/10.1038/nrg3813.
van Kouwenhove M, Kedde M, Agami R. MicroRNA regulation by RNA-binding proteins and its implications for cancer. Nat Rev Cancer. 2011;11:644–56. https://doi.org/10.1038/nrc3107.
Swinburne IA, Meyer CA, Liu XS, Silver PA, Brodsky AS. Genomic localization of RNA binding proteins reveals links between pre-mRNA processing and transcription. Genome Res. 2006;16:912–21. https://doi.org/10.1101/gr.5211806.
Dreyfuss G, Kim VN, Kataoka N. Messenger-RNA-binding proteins and the messages they carry. Nat Rev Mol Cell Biol. 2002;3:195–205. https://doi.org/10.1038/nrm760.
Fu XD, Ares M Jr. Context-dependent control of alternative splicing by RNA-binding proteins. Nat Rev Genet. 2014;15:689–701. https://doi.org/10.1038/nrg3778.
Zheng D, Tian B. RNA-binding proteins in regulation of alternative cleavage and polyadenylation. Adv Exp Med Biol. 2014;825:97–127. https://doi.org/10.1007/978-1-4939-1221-6_3.
Fossat N, Tourle K, Radziewic T, Barratt K, Liebhold D, Studdert JB, Power M, Jones V, Loebel DA, Tam PP. C to U RNA editing mediated by APOBEC1 requires RNA-binding protein RBM47. EMBO Rep. 2014;15:903–10. https://doi.org/10.15252/embr.201438450.
Glisovic T, Bachorik JL, Yong J, Dreyfuss G. RNA-binding proteins and post-transcriptional gene regulation. FEBS Lett. 2008;582:1977–86. https://doi.org/10.1016/j.febslet.2008.03.004.
Li JH, Liu S, Zhou H, Qu LH, Yang JH. starBase v2.0: decoding miRNA-ceRNA, miRNA-ncRNA and protein-RNA interaction networks from large-scale CLIP-Seq data. Nucleic Acids Res. 2014;42:D92–7. https://doi.org/10.1093/nar/gkt1248.
Blin K, Dieterich C, Wurmus R, Rajewsky N, Landthaler M, Akalin A. DoRiNA 2.0--upgrading the doRiNA database of RNA interactions in post-transcriptional regulation. Nucleic Acids Res. 2015;43:D160–7. https://doi.org/10.1093/nar/gku1180.
Anders G, Mackowiak SD, Jens M, Maaskola J, Kuntzagk A, Rajewsky N, Landthaler M, Dieterich C. doRiNA: a database of RNA interactions in post-transcriptional regulation. Nucleic Acids Res. 2012;40:D180–6. https://doi.org/10.1093/nar/gkr1007.
Hu B, Yang YT, Huang Y, Zhu Y, Lu ZJ. POSTAR: a platform for exploring post-transcriptional regulation coordinated by RNA-binding proteins. Nucleic Acids Res. 2017;45:D104–14. https://doi.org/10.1093/nar/gkw888.
Van Nostrand EL, Pratt GA, Shishkin AA, Gelboin-Burkhart C, Fang MY, Sundararaman B, Blue SM, Nguyen TB, Surka C, Elkins K, et al. Robust transcriptome-wide discovery of RNA-binding protein binding sites with enhanced CLIP (eCLIP). Nat Methods. 2016;13:508–14. https://doi.org/10.1038/nmeth.3810.
Lambert N, Robertson A, Jangi M, McGeary S, Sharp PA, Burge CB. RNA Bind-n-Seq: quantitative assessment of the sequence and structural binding specificity of RNA binding proteins. Mol Cell. 2014;54:887–900. https://doi.org/10.1016/j.molcel.2014.04.016.
Kircher M, Witten DM, Jain P, O’Roak BJ, Cooper GM, Shendure J. A general framework for estimating the relative pathogenicity of human genetic variants. Nat Genet. 2014;46:310–5. https://doi.org/10.1038/ng.2892.
Ritchie GR, Dunham I, Zeggini E, Flicek P. Functional annotation of noncoding sequence variants. Nat Methods. 2014;11:294–6. https://doi.org/10.1038/nmeth.2832.
Cooper GM, Stone EA, Asimenos G, Program NCS, Green ED, Batzoglou S, Sidow A. Distribution and intensity of constraint in mammalian genomic sequence. Genome Res. 2005;15:901–13 https://doi.org/10.1101/gr.3577405.
Davydov EV, Goode DL, Sirota M, Cooper GM, Sidow A, Batzoglou S. Identifying a high fraction of the human genome to be under selective constraint using GERP++. PLoS Comput Biol. 2010;6:e1001025 https://doi.org/10.1371/journal.pcbi.1001025.
Fu Y, Liu Z, Lou S, Bedford J, Mu XJ, Yip KY, Khurana E, Gerstein M. FunSeq2: a framework for prioritizing noncoding regulatory variants in cancer. Genome Biol. 2014;15:480. https://doi.org/10.1186/s13059-014-0480-5.
Khurana E, Fu Y, Colonna V, Mu XJ, Kang HM, Lappalainen T, Sboner A, Lochovsky L, Chen J, Harmanci A, et al. Integrative annotation of variants from 1092 humans: application to cancer genomics. Science. 2013;342:1235587. https://doi.org/10.1126/science.1235587.
Genomes Project C, Auton A, Brooks LD, Durbin RM, Garrison EP, Kang HM, Korbel JO, Marchini JL, McCarthy S, GA MV, Abecasis GR. A global reference for human genetic variation. Nature. 2015;526:68–74. https://doi.org/10.1038/nature15393.
Genomes Project C, Abecasis GR, Auton A, Brooks LD, MA DP, Durbin RM, Handsaker RE, Kang HM, Marth GT, McVean GA. An integrated map of genetic variation from 1,092 human genomes. Nature. 2012;491:56–65 https://doi.org/10.1038/nature11632.
Garner C. Confounded by sequencing depth in association studies of rare alleles. Genet Epidemiol. 2011;35:261–8. https://doi.org/10.1002/gepi.20574.
Xu C, Nezami Ranjbar MR, Wu Z, DiCarlo J, Wang Y. Detecting very low allele fraction variants using targeted DNA sequencing and a novel molecular barcode-aware variant caller. BMC Genomics. 2017;18:5. https://doi.org/10.1186/s12864-016-3425-4.
Lu Y, Liu P, James M, Vikis HG, Liu H, Wen W, Franklin A, You M. Genetic variants cis-regulating Xrn2 expression contribute to the risk of spontaneous lung tumor. Oncogene. 2010;29:1041–9. https://doi.org/10.1038/onc.2009.396.
Davidson L, Kerr A, West S. Co-transcriptional degradation of aberrant pre-mRNA by Xrn2. EMBO J. 2012;31:2566–78. https://doi.org/10.1038/emboj.2012.101.
Mortimer SA, Kidwell MA, Doudna JA. Insights into RNA structure and function from genome-wide studies. Nat Rev Genet. 2014;15:469–79. https://doi.org/10.1038/nrg3681.
Pedersen JS, Bejerano G, Siepel A, Rosenbloom K, Lindblad-Toh K, Lander ES, Kent J, Miller W, Haussler D. Identification and classification of conserved RNA secondary structures in the human genome. PLoS Comput Biol. 2006;2:e33. https://doi.org/10.1371/journal.pcbi.0020033.
Khurana E, Fu Y, Chen J, Gerstein M. Interpretation of genomic variants using a unified biological network approach. PLoS Comput Biol. 2013;9:e1002886. https://doi.org/10.1371/journal.pcbi.1002886.
Jiang P, Freedman ML, Liu JS, Liu XS. Inference of transcriptional regulation in cancers. Proc Natl Acad Sci U S A. 2015;112:7731–6. https://doi.org/10.1073/pnas.1424272112.
Consortium GT. The Genotype-Tissue Expression (GTEx) project. Nat Genet. 2013;45:580–5. https://doi.org/10.1038/ng.2653.
Consortium GT, Laboratory DA, Coordinating Center -Analysis Working G, Statistical Methods Groups-Analysis Working G, Enhancing Gg, Fund NIHC, Nih/Nci, Nih/Nhgri, Nih/Nimh, Nih/Nida, et al. Genetic effects on gene expression across human tissues. Nature. 2017;550:204–13. https://doi.org/10.1038/nature24277.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550. https://doi.org/10.1186/s13059-014-0550-8.
Lawrence MS, Stojanov P, Polak P, Kryukov GV, Cibulskis K, Sivachenko A, Carter SL, Stewart C, Mermel CH, Roberts SA, et al. Mutational heterogeneity in cancer and the search for new cancer-associated genes. Nature. 2013;499:214–8. https://doi.org/10.1038/nature12213.
Lochovsky L, Zhang J, Fu Y, Khurana E, Gerstein M. LARVA: an integrative framework for large-scale analysis of recurrent variants in noncoding annotations. Nucleic Acids Res. 2015;43:8123–34. https://doi.org/10.1093/nar/gkv803.
Lochovsky L, Zhang J, Gerstein M. MOAT: efficient detection of highly mutated regions with the mutations overburdening annotations tool. Bioinformatics. 2017. https://doi.org/10.1093/bioinformatics/btx700.
Vogelstein B, Kinzler KW. Cancer genes and the pathways they control. Nat Med. 2004;10:789–99. https://doi.org/10.1038/nm1087.
Khurana E, Fu Y, Chakravarty D, Demichelis F, Rubin MA, Gerstein M. Role of non-coding sequence variants in cancer. Nat Rev Genet. 2016;17:93–108. https://doi.org/10.1038/nrg.2015.17.
Forbes SA, Beare D, Bindal N, Bamford S, Ward S, Cole CG, Jia M, Kok C, Boutselakis H, De T, et al: COSMIC: high-resolution cancer genetics using the catalogue of somatic mutations in cancer. Curr Protoc Hum Genet 2016, 91:10 11 11–10 11 37. https://doi.org/10.1002/cphg.21.
Alexandrov LB, Stratton MR. Mutational signatures: the patterns of somatic mutations hidden in cancer genomes. Curr Opin Genet Dev. 2014;24:52–60. https://doi.org/10.1016/j.gde.2013.11.014.
Bailey TL. DREME: motif discovery in transcription factor ChIP-seq data. Bioinformatics. 2011;27:1653–9. https://doi.org/10.1093/bioinformatics/btr261.
Zhang J, Liu J, Lee D, Feng J, Lochovsky L, Lou S, Rutenberg-Schoenberg M, Gerstein M: RADAR GitHub Repository. https://github.com/gersteinlab/RADAR. Accessed 8 Oct 2018.
Zhang J, Liu J, Lee D, Feng J, Lochovsky L, Lou S, Rutenberg-Schoenberg M, Gerstein M: Zenodo repository for RADAR. https://doi.org/10.5281/zenodo.1451838. Accessed 8 Oct 2018.
Consortium EP. An integrated encyclopedia of DNA elements in the human genome. Nature. 2012;489:57–74 https://doi.org/10.1038/nature11247.
Alexandrov LB, Nik-Zainal S, Wedge DC, Aparicio SA, Behjati S, Biankin AV, Bignell GR, Bolli N, Borg A, Borresen-Dale AL, et al. Signatures of mutational processes in human cancer. Nature. 2013;500:415–21 https://doi.org/10.1038/nature12477.
Cancer Genome Atlas Research N, Weinstein JN, Collisson EA, Mills GB, Shaw KR, Ozenberger BA, Ellrott K, Shmulevich I, Sander C, Stuart JM. The Cancer Genome Atlas Pan-Cancer analysis project. Nat Genet. 2013;45:1113–20 https://doi.org/10.1038/ng.2764.
Harrow J, Denoeud F, Frankish A, Reymond A, Chen CK, Chrast J, Lagarde J, Gilbert JG, Storey R, Swarbreck D, et al: GENCODE: producing a reference annotation for ENCODE. Genome Biol 2006, 7 Suppl 1:S4 1–9. https://doi.org/10.1186/gb-2006-7-s1-s4.
Harrow J, Frankish A, Gonzalez JM, Tapanari E, Diekhans M, Kokocinski F, Aken BL, Barrell D, Zadissa A, Searle S, et al. GENCODE: the reference human genome annotation for the ENCODE Project. Genome Res. 2012;22:1760–74 https://doi.org/10.1101/gr.135350.111.
Acknowledgements
We thank Brenton Graveley, Gene Yeo, Peter Freese, and Eric Van Nostrand for the useful discussion on eCLIP and RBNS data.
Peer review information
Anahita Bishop and Kevin Pang were the primary editors of this article and managed its editorial process and peer review in collaboration with the rest of the team.
Funding
This work was supported by the National Institutes of Health (grant number 1U24HG009446-01), AL Williams Professorship, and in part by the facilities and staff of the Yale University Faculty of Arts and Sciences High Performance Computing Center.
Author information
Authors and Affiliations
Contributions
JZ, JL, and MG conceived the study and wrote the manuscript. JZ, JL, DL, LL, and JF wrote the framework and performed the method evaluation. DL and JF developed the website. LL and MR carried out the studies associating RNA structure. JL, DL, and SL participated in the motif analysis. The authors have read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary information
Additional file 1:
Large Supplemental Tables from RADAR providing results and summaries regarding the data used, but are too large to be put in PDF format.
Additional file 2:
Supplemental Figures and Tables for RADAR providing more information on results and methods.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Zhang, J., Liu, J., Lee, D. et al. RADAR: annotation and prioritization of variants in the post-transcriptional regulome of RNA-binding proteins. Genome Biol 21, 151 (2020). https://doi.org/10.1186/s13059-020-01979-4
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s13059-020-01979-4