RBP-Var: a database of functional variants involved in regulation mediated by RNA-binding proteins.

Transcription factors bind to the genome by forming specific contacts with the primary DNA sequence; however, RNA-binding proteins (RBPs) have greater scope to achieve binding specificity through the RNA secondary structure. It has been revealed that single nucleotide variants (SNVs) that alter RNA structure, also known as RiboSNitches, exhibit 3-fold greater local structure changes than replicates of the same DNA sequence, demonstrated by the fact that depletion of RiboSNitches could result in the alteration of specific RNA shapes at thousands of sites, including 3' UTRs, binding sites of microRNAs and RBPs. However, the network between SNVs and post-transcriptional regulation remains unclear. Here, we developed RBP-Var, a database freely available at http://www.rbp-var.biols.ac.cn/, which provides annotation of functional variants involved in post-transcriptional interaction and regulation. RBP-Var provides an easy-to-use web interface that allows users to rapidly find whether SNVs of interest can transform the secondary structure of RNA and identify RBPs whose binding may be subsequently disrupted. RBP-Var integrates DNA and RNA biology to understand how various genetic variants and post-transcriptional mechanisms cooperate to orchestrate gene expression. In summary, RBP-Var is useful in selecting candidate SNVs for further functional studies and exploring causal SNVs underlying human diseases.

INTRODUCTION

A single nucleotide variant (SNV) can result in complex diseases through various molecular mechanisms, including amino-acid changes at the protein level, alteration of regulatory elements at the transcription level and disruption of RNA–protein interaction at the post-transcription level (1–4). Previous studies mainly concentrated on the effect of SNVs on protein structure (5,6) and the misregulation of transcriptional regulation (7,8), known as cSNVs and rSNVs, respectively. However some SNVs, which are located in non-coding regions and have no relationship with transcriptional regulation (8), can have an observable phenotype or result in a disease. Therefore, some SNVs may play crucial roles in post-transcriptional gene regulation by altering micro (mi)RNA targets, RNA secondary structure and RNA–protein interactions (1,4,9–12).

Intensive attempts have been undertaken to better understand RNA-binding proteins (RBPs) and RNA–protein interactions (13). The advent of next-generation sequencing and relevant molecular techniques have allowed the development of systematic experimental protocols for the identification of protein-binding sites on RNA using the immunoprecipitation of RBPs, with or without in vivo RNA–protein crosslinking, known as crosslinking immunoprecipitation sequencing (CLIP-seq) and RNA immunoprecipitation sequencing (RIP-seq), respectively (14,15). CLIP-seq that combines high-throughput sequencing (HITS)-CLIP (16), photoactivatable ribonucleoside-enhanced (PAR)-CLIP (17) and individual nucleotide resolution crosslinking immunoprecipitation (iCLIP) (18), is a powerful tool for investigating the in vivo binding sites of RBPs on a transcriptome-wide scale, providing higher resolution than previous technologies, such as RIP-seq, for identifying protein-binding sites on RNAs.

Identifying changes in RBP binding induced by SNVs for a messenger RNA (mRNA) in vivo remains challenging because RBPs have much more scope to achieve specificity through the secondary structure of RNA and SNVs could change intramolecular bonds which may twist RNA molecules into hairpins, stem-loops and various other bumps and bulges (13). The impact of SNVs on RNA structure and gene regulation has been widely studied (1–4). Wan et al. (19) studied the initial landscape and variation of RNA secondary structures in a human family trio (mother, father, offspring) and revealed that ≈15% of all transcribed SNVs are RiboSNitches that alter the local RNA structure. These RiboSNitches were also significantly depleted around predicted miRNA target sites and protein-binding sites on RNA. Moreover, Maticzka et al. (20) demonstrated the importance of RNA secondary structure for conferring binding specificity in a subset of RBPs using both RNA structure and sequence features to develop a machine learning-based approach for predicting protein-binding sites on RNA from CLIP data. Fukunaga et al. (21) developed the CapR software that determines the probability of secondary structure throughout an RNA molecule for RBP binding and found that the human RBP Pumilio-2 has a preference for hairpin loop structures. Overall, RNA secondary structures exert a significant effect on RBP binding, and the RNA interactome network can be disrupted once the RNA structures are transformed.

Recently, great progress has been made in identifying mutations in RBPs associated with disease risk (22) such as mutations in the RBPs FUS and TDP-43 that can cause amyotrophic lateral sclerosis (23). However, it remains unclear what changes occur to RNA substrates in diseases and their effect on RNA structure and RNA–protein interaction induced by SNVs or RNA editing sites (24). Therefore, a comprehensive database is required for the annotation of functional variants involved in post-transcriptional interactions. Here, we developed RBP-Var (Figure 1), an integrated database for the annotation of functional variants involved in RBP-mediated post-transcriptional regulation and interaction. We also determined trait/disease-associations and deciphered gene regulation for SNVs located in potential protein-binding sites on RNA. Furthermore, we developed a heuristic scoring system based on the functional confidence of a variant, including SNV or RNA edit, to assist comparison among annotations.

Figure 1.

Workflow to identify functional rbSNVs by RBP-Var.

DATA COLLECTION AND PROCESSING

Data sources for RBP-Var

RBP-Var provides extensive annotation for the functional variants involved in post-transcriptional interactions. We collected all the available processed data on RBPs from starBase (25) and CLIPdb (26) and extracted sequencing data from Gene Expression Omnibus (GEO) (27). Genomic coordinates of RBPs derived from GEO were converted from human genome assembly hg18 into hg19 using the University of California, Santa Cruz (UCSC) LiftOver tool (http://genome.ucsc.edu). We also retrieved a collection of expression quantitative trait loci (eQTL) from previous studies (28–37), human trait-associated SNVs from the genome-wide association study (GWAS) catalog (38), somatic mutations in cancer from COSMIC (39) and disease-associated mutations from ClinVar (40). In addition, RBP-Var includes collection of RBP binding SNVs (rbSNVs) which may disrupt the binding of RBPs, derived from 112 CLIP-seq data sets (60 RBPs) and 319 motif matching sites (153 RBPs). These functional RBP binding SNVs include deletions, insertions, single nucleotide polymorphism (SNPs), multiple nucleotide polymorphisms (MNPs), microsatellites, insertions/deletions from dbSNP database (version 142) as well as RNA editing (A-to-I editing sites) from RADAR and DARNED (Table 1).

Table 1.

Database content of RBP-Var as of 10 October, 2015

Data type	Sources	Description	Records	Genomic coverage (bp)
CLIP-seq	starBase, CLIPdb, GEO	112 data sets for 60 RBPs	52,528,288	847,717,483
Predicted binding (PWMs)	CISBP-RNA, RBPDB	319 PWMs for 153 RBPs	119,777,778	822,156,225
dbSNP	dbSNP v142	including SNP,InDel and MNP	113,128,211	113,128,211
RNA editing	RADAR,DARNED	A->G editing	2,576,460	2,576,460
miRNA targets	TargetScan, miRanda, miRNASNP	miRNA-SNP pairs	3,361,915	26,958,420
eQTLs	MuTher, SCAN, seeQTL, GTEx, Harvard	SNP-gene pairs	4,546,890	4,546,890
dsQTLs	dsQTL Browser	SNP-gene pairs	214,522	214,522
LD proxies	HapMap, 1000Genome	r2>0.8 by Plink	4,600,670	4,600,670
6mA	MeT-DB	6mA peaks called by MACS	466,037	127,588,574
Splicing sites	H-DBAS	including 7 splicing patterns	1,669,641	158,324,909
trait/disease associations	GWAS, COSMIC, ClinVar, TCGA, DISEASES	SNP-trait or gene-trait pairs	8,955,125	8,955,125

Position weight matrix (PWM) motif matching

Many RBPs interact with mRNAs via a limited set of modular RNA-binding domains, including the RNA recognition motif, heterogeneous nuclear ribonucleoprotein K-homology domain and zinc fingers (11,41). The RNA-binding domains of RBPs initially determine the specificity and preferences of RNA binding with specific sequence motifs (42). Therefore, all PWMs from the catalog of inferred sequence binding preferences of RNA binding proteins (CIS-BP-RNA) (42) and the RNA-binding protein database (RBPDB) (43) deposited in the AURA database (44) were used to call motif matches in the transcriptome. The possible k-mers were aligned to the human transcriptome using MAST in the MEME suite (45) to give a final motif mapping with a match score >0 and P-values < 0.0001.

Impacts of SNVs on RNA structures

Based on human reference transcripts retrieved from UCSC (hg19), we changed the corresponding allele in a given sequence of the reference transcript to an alternative allele as an altered transcript for each SNV. Then, the RNAfold program was employed (46) with default parameters to illustrate the secondary structure and calculate the minimal free energy (MFE, ΔG) for both transcript sequence of wild-type and mutant. Energy changes in RNA structures (ΔΔG) were calculated from MFE differences (ΔΔG = |ΔGalt -ΔGref|), where ΔGref and ΔGalt are the MFE of the reference transcript and altered transcript, respectively. Moreover, for SNP and RNA editing sites, we employed the RNAsnp (2) with default parameters to estimate the mutation effects on local RNA secondary structure in terms of empirical P-values based on global folding in correlation measure by using RNAfold. For insertions, deletions and SNVs of other types, based on changes of MFE,we evaluated the effects on RNA secondary structure in terms of empirical P-values calculated from cumulative probabilities of the Poisson distribution.

Impact of SNPs on miRNA–RNA interactions

Single nucleotide polymorphisms (SNPs) in the target sites of miRNAs may destroy or create miRNA-binding sites on RNAs, which lead to loss-of-function and/or gain-of-function miRNA–RNA interactions. To assess the impact of each SNP on miRNA–RNA interaction, TargetScan (47) and miRanda (48) were used to predict target sites around each SNP with ±25 bp sequences for Ref-transcripts and Alt-transcripts, separately. In addition, we compared miRNA target score versus target score induced by SNV by empirical P-values calculated from cumulative probabilities of the Poisson distribution. We analyzed target alteration of miRNA by SNPs in 3′-UTR of mRNA and miRNA target alteration by SNPs in miRNA seed regions. We found that 888 674 pairs miRNA–RNA interactions was altered significantly (P-value < 0.01) by 253 428 SNPs in 3′-UTR of mRNA in a total of 1 654 644 pairs miRNA–RNA interactions. Meanwhile, we found that 302 723 pairs miRNA–RNA interactions was altered significantly (P-value < 0.01) by 201 SNPs in miRNA seed regions in total 556 008 pairs miRNA–RNA interactions.

Analyzing linkage disequilibrium (LD) proxies

Due to genetic correlation among nearby SNVs, functional analysis on single SNV is unreliable. To extend the post-regulatory feature of a single SNV, LD correlations between SNVs were analyzed, and the set of SNVs that had a strong LD (r2>0.8) with rbSNVs was defined as LD-proxies of rbSNVs. LD data were derived from both merged HapMap phases I+II+III genotype data (49) and integrated 1000-genome phase III release data (50). Data from all populations deposited in these two projects were utilized to perform LD analyses in an LD-window of 500 kb using Plink (51).

Trait/disease-associated label

SNVs occurring at the binding sites of RBPs could contribute to cis-modulation of gene expression by changing the affinity of a regulatory protein to RNA sequence (52). We used trait/disease-associated SNVs (TASs) obtained from GWAS (38), COSMIC (39) and ClinVar (40) to annotate SNVs in potential protein-binding sites identified from 112 CLIP-seq data sets. Expression trait usually differs from other classical complex traits, because the measured mRNA or protein trait is the product of a single gene within a specific chromosomal location. It is known that eQTLs are genomic loci that regulate the expression level of mRNAs and play a crucial role in deciphering gene regulation and spatio-temporal specificity. RBP-Var provides eQTL labels and spatio-temporal labels for the included SNVs. We collected eQTLs from experimentally supported eQTL databases (35–37) and the eQTL browser (http://eqtl.uchicagoedu/cgi-bin/gbrowse/eqtl/) (28–34) to provide association labels for SNVs. Tissue and developmental stage information was labeled according to cell type from which protein-binding sites on RNA were identified.

Classifying functional variants

To better interpret and comprise the catalog of rbSNV information, we developed a heuristic scoring system based on the functional confidence of variants. The scoring system represents with increasing confidence that a variant lies in a functional location and probably results in a functional consequence (i.e. alteration of RBP binding and a gene regulatory effect) (Table 2). For example, we consider variants that are known eQTLs as significant and label them as Category 1. Within Category 1, subcategories indicate additional annotations ranging from the most informational variants (1a, variant may change the motif for RBP binding) to the least informational variants (1e, variant only has a motif for RBP binding).

Table 2.

The variant classification scheme of RBP-Var

Category	Description	# of SNV	# of RNA editing
Likely to affect RBP binding, RNA secondary structure and expression
1a	eQTL+RBP binding+matched RBP motif+riboSNitch+miRNA targets	1	0
1b	eQTL+RBP binding+any RBP motif +riboSNitch+miRNA targets	1,059	0
1c	eQTL+RBP binding+any RBP motif +riboSNitch	2,417	1
1d	eQTL+RBP binding+any RBP motif	334	0
1e	eQTL+RBP binding	916	0
Likely to affect RBP binding, RNA secondary structure
2a	RBP binding+matched RBP motif +riboSNitch+miRNA targets	23	1
2b	RBP binding+any RBP motif +riboSNitch+miRNA targets	29,415	483
2c	RBP binding+any RBP motif +riboSNitch	177,806	5,758
2d	RBP binding+riboSNitch	616,239	46,963
Less likely to affect RBP binding		4,984
3a	RBP binding+any RBP motif+miRNA targets	15	49
3b	RBP binding+matched RBP motif	28,599	0
Minimal possibility to affect RBP binding		11,477
4	RBP binding+any RBP motif	844	534
5	RBP binding+miRNA targets	85	209
6	RBP binding or RBP motif hit	29,415	14

WEB INTERFACE

Data retrieving in RBP-Var could be implemented based on SNP-centric or gene-centric. SNP-centric data retrieving was appropriate to analyze results from genetic studies, especially the results of high-throughput studies, and then provided evidence for further functional studies to identify causal SNPs and shed light on underlying molecular mechanisms. Gene-centric searches were useful in candidate SNP selections that are based on genes of interest. In addition, RBP-Var provided RNA-editing retrieving, which could decipher the potential function of RNA editing involved in post-regulation and gene expression (24). Furthermore, we applied the JBrowse Genome Browser (http://jbrowse.org) to establish the ‘Browser’ page to visualize genome-wide signals of RBP binding for each CLIP-seq data set. Users could browse all binding signals across specific genomic region of their interests. The searching result for SNP rs1802295 in RBP-Var database was used as an example of web-interface (Figure 2).

Figure 2.

Web-interface of RBP-Var. (A–H) The snapshot of searching result for SNP rs1802295 in RBP-Var database.

We found that 874 160 (0.77%) SNVs in the SNP database v.142 were overlapped with RNA binding sites as identified by CLIP-seq. These SNVs included 835 130 SNPs, 26 854 deletions, 11 323 insertions, 166 insertions/deletions, 679 multiple nucleotide polymorphisms (MNPs) and 8 microsatellites. Subsequently, we employed our heuristic scoring system to annotate these SNVs, and we found 828 210 SNVs classified into Category 1/2 (Table 2), which were composed of 802 765 SNPs, 23 570 deletions, 1052 insertions, 152 insertions/deletions, 663 MNPs and 8 microsatellites (Figure 3A).

Figure 3.

The features of RBP-Var. (A) The component of rbSNVs derived from dbSNP. (B) Comparison of SNVs with category score 1/2 deposited in RBP-Var and RegulomeDB with damaging SNVs predicted by SIFT and PolyPhen2.(C) Comparison of RNA editing events with category score 1/2 deposited in RBP-Var with damaging RNA editing events predicted by SIFT and PolyPhen2.

APPLICATIONS

To verify the heuristic system, we compared the functional rbSNVs of RBP-Var with regulatory SNVs from RegulomeDB (8) and deleterious SNVs predicted by SIFT (5) and PolyPhen2 (6). We treated SNVs with category score 1/2 in RBP-Var and as regulatory SNVs in RegulomeDB. A total of 1 712 307 SNVs were identified as deleterious or functional SNVs by at least one of the four classification methods. Among 828 210 functional rbSNPs identified by RBP-Var, 40 904 were classified into functional SNVs by RegulomeDB, 166 647 were presented and predicted as deleterious by PolyPhen2 and 109 364 by SIFT (Figure 3B), whereas, the rest 651 784 rbSNVs were identified uniquely in RBP-Var. Then, we compared the rbSNVs uniquely identified in RBP-Var with GWAS-lead SNVs to determine causal disease variants. As a result, we identified 93 unique rbSNVs that were associated with human diseases (Supplemental Table S1). Similarly, we identified 7108 and 80 928 unique rbSNVs with RBP-Var score 1/2 in ClinVar and COSMIC databases, respectively (Supplemental Tables S2 and S3). These observations indicated that RBP-Var could identify more rbSNVs which could be implicated in human disease and biological functions even though they did not affect protein coding.

As one of the non-coding regulatory regions, untranslated regions (UTRs) are known to have many sequence and structural motifs that can regulate translational and transcriptional efficiency and the stability of mRNAs through interactions with RBPs and other non-coding RNAs like microRNAs (miRNAs) (4). Here, we found that 38 GWAS-lead SNPs located in the UTRs were classified with RBP-Var score 1/2 (Supplemental Table S1). For example, SNP rs1802295 (C/T) which is associated with type 2 diabetes in individuals of South Asian ancestry (53), but its underlying mechanism participated in type 2 diabetes is still unknown (54). Our result revealed SNP rs1802295 is located in the 3′-UTR of VPS26A, a region containing motif of binding by protein PABPC1L and PABPC3 based on PWM matching. This motif is likely disrupted by this SNP as the core component C in the motif would be changed to U (Figure 4A). Meanwhile, the experimental binding intensity around this SNP on mRNA of VPS26A includes binding by protein AGO2, PTBP1 from HeLa cell line, and IGF2BP3, WDR33, ELAVL1 from HEK293 cell line (Figure 4A). All of these RBP bindings may be affected as the secondary structure of mRNA of VPS26A is changed (Figure 4B) and the targeting site of miRNA hsa-miR-510 is potentially lost when the SNP rs1802295 occurs. Thus, we speculated that loss of targeting by miRNA hsa-miR-510 could up-regulate VPS26A and induce a cis-eQTL in blood according to Blood eQTL browser (55) by post-transcriptional regulation mediated by these RBPs. Hence, our results obtained by using the RBP-Var database provide additional insight into mechanisms underlying SNV-induced structural changes in the UTRs as such changes potentially affect the stability of the mRNA or disrupt the interactions with RBPs and targeting of microRNAs present in the UTRs.

Figure 4.

Application of RBP-Var for an SNV related to type-2 diabetes. (A) The RBP binding signal for SNP rs1802295, which is a functional rbSNV associated with type 2 diabetes. The motif of binding by protein PABPC1L and PABPC3 are showed on the top of rs1802295. And the binding intensity around this SNP on mRNA of VPS26A by protein AGO2, PTBP1 from HeLa cell line, and IGF2BP3, WDR33, ELAVL1 from HEK293 cell line is illustrated in different tracks. The targeting site of miRNA hsa-miR-510 is marked by yellow vertical line. (B) The optimal secondary structure of RNA sequence from 200 bp of flanking regions on either side of SNP rs1802295 for wild-type of VPS26A and mutant of VPS26A induced by rs1802295. The green and red nucleotides represent the local sequence in interval of 50 bp for wild-type and mutant of VPS26A, respectively.

RNA editing is a molecular process through which some cells make discrete changes to specific nucleotide sequences within an RNA molecule after it has been generated by RNA polymerase. Adenosine to Inosine (A-to-I) editing is the main form of RNA editing in mammals (56) and occurs in regions of double-stranded RNA. Adenosine deaminases acting on RNA are RNA-editing enzymes involved in the hydrolytic deamination of adenosine to inosine (A-to-I editing). There are many effects of A-to-I editing, arising from the fact that I behaves as if it is guanine in translation and secondary structure formation. These effects include alteration of coding capacity (57), altered miRNA or siRNA target populations, heterochromatin formation, inhibition of miRNA and siRNA processing, and altered splicing. In this study, only a small fraction of RNA editing events deposited in the RADAR/DARNED (58,59) database had deleterious effects on protein structure as predicted by SIFT (5) and PolyPhen2 (6) (Figure 3C). In total, we found that 52,584 unique RNA editing events were classified into category 1/2, which might affect RBP binding as well as RNA secondary structure (Figure 3C).

Taking RNA-editing site chr8:103841636–103841636 for example, it was predicted to be tolerated by SIFT with score 0.55 and be benign by PolyPhen2 with score 0 though this A-to-I RNA editing results in a serine-to-glycine (Ser-to-Gly) amino acid substitution, occurring at residue 367 of AZIN1. However, recent study revealed that this edited form has a stronger affinity to antizyme, and the resultant higher AZIN1 protein stability promotes cell proliferation and predisposes to human hepatocellular carcinoma compared with wild-type AZIN1 protein (60) though underlied mechanism remained unknown. Intriguingly, we found that this A-to-I RNA editing may disrupt the RNA secondary structure of transcript ENST00000518697 of gene AZIN1 (p = 0.0074). And the resultant change of RNA secondary structure may influence the binding of RBPs including RBM47, eIF4AIII, AGO2, ELAVL1, FMR1 and PTBP1. Meanwhile, the searching result of our database showed that this A-to-I RNA editing site was N6-methylated in human hepatoma cell line HepG2. Our observation indicates that this A-to-I RNA editing may deplete the adenine N6-methylation (6mA) as the nucleotide adenine (A) would be transformed into nucleotide inosine (I). Hence, RBP-Var broadens the spectrum for the implication of RNA editing events by evaluating their effects on the alteration of post-regulatory capacity.

DISCUSSION

In this study, we included extensive information on annotated and computed post-regulatory elements in the human genome, and we combined a large array of data sources into a single and integrated database to quickly generate prioritized hypotheses for the function of variants affecting both coding and non-coding regions in a genome. This novel database has a convenient interface that allows the easy query submission and provides an instant classification of significant variants. The SNV summary page allows users to quickly form a hypothesis regarding the true RBP-mediated post-transcriptional consequence of a variant.

Guo et al. (61) constructed the rSNPbase database, which provides a similar annotation of SNVs integrated with RBP-mediated post-transcriptional regulation, using data sets of four RBPs derived from RIP-chip or RIP-seq analysis. Our database provides additional information by prioritizing SNVs within general regulatory regions based on specific RBPs, RNA secondary structure and PWM information.

Our simple heuristic scoring system may be improved over time, as more functional rbSNVs will be validated. Our analysis method and database are centered on a likely disruption of a protein-RNA interaction via alteration of RNA secondary structure and regulation of gene expression by rbSNVs if they are eQTLs. It is likely that there are additional sources of data that reinforce each other in a different manner that need to be explored; however, the scoring system provides significant enrichment along with better category scores, as shown when compared to GWAS-lead, ClinVar and COSMIC rbSNVs. Hence, our database and scoring system provide the best current method for annotating and prioritizing variants involved in RBP-mediated post-transcriptional regulation. Overall, we comprehensively integrated eQTLs, protein-binding sites on RNA, RNA secondary structure, miRNA, and SNV information along with GWAS, COSMIC, ClinVar, 6mA, alternative splicing, and miRNA expression profiles to seek potential post-functional SNVs.

In conclusion, this database provides a user-friendly web interface, allowing users to rapidly find whether SNVs of interest can transform the secondary structure of RNA and identify RBPs whose binding may be subsequently disrupted. Moreover, RBP-Var can assess the impact of each SNV on miRNA–RNA interaction as SNVs may destroy or create miRNA-binding sites, which result in loss-of-function and/or gain-of-function miRNA–RNA interactions. RBP-Var is a useful resource for benchmarking the mutations or RNA-editing events that cause disease by changing post-transcriptional interaction and regulation. In the future, RBP-Var will be updated frequently and extended to other species. We are dedicated to maintaining and improving RBP-Var, since it is a useful resource for the research community.

SUPPLEMENTARY DATA

Supplementary Data are available at NAR Online.