A Survey of Regulatory Interactions Among RNA Binding Proteins and MicroRNAs in Cancer

Abstract

Recent advances in genomics and proteomics generated a large amount of trans regulatory data such as those mediated by RNA binding proteins (RBPs) and microRNAs. Since many trans regulators target 3′ UTR of mRNA transcripts, it is likely that there would be interactions, i.e., competitive or cooperative effect, among these trans factors. We compiled the available RBP and microRNA binding sites, mapped them to the mRNA transcripts, and correlated the binding data with mRNA expression data generated by The Cancer Genome Atlas (TCGA). We separated pairs of RBPs and microRNAs into three scenarios: those that have overlapping target sites on the same mRNA transcript (overlapping), those that have target sites on the same mRNA transcript but non-overlapping (neighboring), and those that do not target the same mRNA transcript (independent). Through a regression analysis on expression profiles, we indeed observed interaction effect between RBPs and microRNAs in the majority of the cancer expression data sets. We further discussed implication of such widespread interactions in the context of cancer and diseases.

Introduction

Post-transcriptional regulation is an important yet complex system that modulates expression level of genes inside cells. The commonly studied post-transcriptional regulatory mechanisms include alternative splicing, mRNA modification, regulation by RNA binding proteins (RBP) or by microRNAs. Dysregulation in post-transcriptional regulation has been implicated in many human diseases, especially those mediated by RBPs (Lukong et al., 2008) and by microRNAs (Calin and Croce, 2006). Since many RBPs and microRNAs target the 3′ UTR of the same mRNA transcripts, it is likely that there exist competitive or cooperative interactions among these trans regulators. For example, these trans regulators could compete for the same binding sites on the 3′ UTR; alternatively, binding by one trans regulator could affect the local mRNA conformation and influence binding by other trans regulators. Indeed, a number of earlier studies have investigated interactions between these regulators. Perhaps the most widely studied is the competition between microRNAs, either via competing for the same AGO protein complex or via indirectly modulating the abundance of mRNAs. Califano, Pandolfi, and colleagues proposed the idea that different types of RNAs, either mRNAs or non-coding RNAs, can compete for the same microRNAs, and changes in RNA abundance can modulate the regulatory effect of other microRNAs, hence the concept of ceRNA (competing endogenous RNAs) (Salmena et al., 2011; Sumazin et al., 2011; Li et al., 2014; Tay et al., 2014). Marks and colleagues observed that overexpression of one type of microRNA can overload the AGO proteins thus reduce the possibility of other microRNAs being loaded onto the RISC complex (Khan et al., 2009). In addition to these computational analyses, there are experimental evidence supporting such competitive effects among microRNAs (Hua et al., 2006; Castanotto et al., 2007; Ebert et al., 2007).

Doyle and Tenenbaum (2014) presented a general framework on competition among RBPs; however, there have not been any large-scale experimental survey in this area except for a few cases involving HuR. It was shown that, HuR and AUF1, two RBPs that cause opposing effect on their targets (the binding of HuR stabilizes mRNA while binding of AUF1 promotes degradation), compete for common target sequences and can modulate each other’s effect (Barker et al., 2012). Competition between HuR and wild-type TTP for binding to the HuR transcript has been implicated in HuR regulation and its cytoplasmic localization (Al-Ahmadi et al., 2009). HuR is a widely studied RBP and binds to many target mRNAs, so it is not surprising that many of the reported competition events involved HuR. Zarnack et al. (2013) conducted genome-wide iCLIP experiments and discovered that HNRNP (Heterogeneous Nuclear Ribonucleoprotein C) competes with the splicing factor U2AF65 at many splice sites, which prevents the exonization of Alu elements and protects the integrity of human transcriptome. In another example, it was found that two families of splicing factors, CELF and RBFOX, can antagonize each other in the regulation of global splicing events (Gazzara et al., 2017). Regarding interactions between RBPs and microRNAs, Srikantan et al. (2012) presented an overview of potential competition and cooperation between HuR and microRNAs. For example, it was shown that HuR can recruit miR-19 and the associated RISC complex to the mRNA transcript of small GTPase RhoB and repress protein translation; the interaction between HuR and miR-19 can be decoupled by UV exposure, which prevents the down-regulation of RhoB by miR-19 (Glorian et al., 2011). Such regulatory role of HuR has later been reported for let-7, miR-17-92 and other microRNAs as well (Kim et al., 2009; Kundu et al., 2012; Gunzburg et al., 2015; Mihailovich et al., 2015).

Despite these published studies, to the best of our knowledge, there still lacks a general analysis of global regulatory interactions between RBPs and microRNAs beyond the well-studied proteins such as HuR. In this study, we limited our list of RBPs only to those that target 3′ UTR of mRNAs and excluded those binding events at the 5′ UTR or coding regions (CDS). We also did not include splicing factors that bind to intron and exon junctions, which warrants a separate study. For RBP binding sites, we collected and mapped the binding sites derived from published experimental studies, specifically the POSTAR database (Hu et al., 2017). These RBP binding sites were further refined by considering RNA accessibility to include only those deemed of high confidence. The microRNA target sites were derived from computational predicted databases. Based on the binding sites of RBPs and microRNAs, Figure 1A shows three possible scenarios, which we refer to as “Overlapping,” “Neighboring,” and “Independent,” respectively.

FIGURE 1

Materials and Methods

Mapping MicroRNA and RBP Target Sites

The microRNA targeting information was downloaded from TargetScan database (version 7.2) (Agarwal et al., 2015), which as of January 2019 contained 321 human microRNAs and 11,354 target genes. We chose to use TargetScan since it was more frequently updated (last updated in March 2018). The default parameters for target prediction were used. We did not limit our study to only experimentally validated targets such as those curated by TarBase (Karagkouni et al., 2018), since we aimed to finding global competitive effect between trans factors, rather than investigating individual factors. The RBP binding sites were downloaded from POSTAR2¹ (Hu et al., 2017), the coordinates of these binding sites were based on GENCODE Human release 27. The coordinates of microRNA and RBP binding sites were mapped from hg19 to hg38 by using the LiftOver tool. The definition of 3′ UTR from GENCODE was used. Details on the further mapping and processing of these target sites are described in Results.

Randomization of RBP and MicroRNA Target Sites

We conducted randomization on the overlapping RBP (or microRNA) target sites to calculate the baseline probability for two RBPs to have overlapping target sites on an mRNA transcript, provided that they already have target sites on the same mRNA transcript. For a given mRNA, we randomly re-distributed RBP target sites on the 3′ UTR region and compared the number of overlapping RBP target sites with the real data. We conducted randomization 10 times and took the average number of overlapping target sites.

TCGA Datasets

The mRNA and microRNA expression data were downloaded from Broad Institute Firehose website². Genes that have missing value among over 30% of the samples were not considered. We used the knn.impute package to fill in the missing value. The final summary of The Cancer Genome Atlas (TCGA) data is shown in Supplementary Table S1.

Statistical Analysis

We conducted regression analysis as formulated in Eq. 1 to measure the significance of interactions among trans regulators. We used the R package p.adjust to calculate the p-values adjusted by Benjamini–Hochberg method (q-value < 0.05). The p-value for each term tests the null hypothesis that the coefficient is equal to zero (no effect). A low p-value (<0.05) indicates that there is a larger interaction probability between the two trans regulators. We used the R function ks.test (x, y, “greater”) to carry out the Kolmogorov–Smirnov test; x and y represent the p-values calculated from Eq. 1 from Overlapping, Neighboring, and Independent regulator pairs. The “greater” option tests the null hypothesis that x is stochastically smaller than y.

Results

Mapping RBP Targets

The analysis flowchart is shown in Figure 1B and details on how we derived RBP and microRNA binding sites can be found in the section “Materials and Methods.” Supplementary Table S2 summarizes the downloaded RBP binding data. The POSTAR2 database had collected global binding data for 171 RBPs (as of January 2019), mostly were determined by CLIP-Seq or similar approaches and analyzed by software such as PARalyzer or CIMS (Corcoran et al., 2011; Zhang and Darnell, 2011). We noticed that there was often more than one study on a unique RBP; we pooled the binding sites determined from these separate studies. For those genes with multiple alternative mRNA transcripts, we used GENCODE as a reference and selected the transcript with the longest 3′ UTR. Many of these RBP binding sites are likely only bound by RBP in specific cell types due to cell type specific abundance of the RBPs and the target mRNAs, or the competitive or cooperative interaction between RBPs. At this step in our study, we simply mapped all the potential RBP binding sites and their spatial relationship in preparation for the subsequent regression analysis on expression levels.

We noticed that many of the RBP binding sites reported by the PARalyzer or eCLIP were very long and spanned over 50 nucleotides, which was characteristic of the broad peaks called by peak calling software (Van Nostrand et al., 2016) or was the result of multiple overlapping binding sites being merged together (Corcoran et al., 2011). We manually examined the annotation of these RBPs and the RNA binding domains (RBD) within these proteins; the majority of these RBPs or RBDs had evidence to bind to single stranded RNAs (ssRNA) instead of double stranded RNAs (dsRNA). Therefore, we further refined these RBP binding sites by calculating the local RNA structure accessibility and only retained those that were considered to be accessible to RBPs. We used the RNAplfold program in the Vienna Package and used recommended parameters (W = 80 and L = 40) (Lorenz et al., 2011). It was previously shown that inclusion of local RNA accessibly could help improve the identification of RBP binding events (Tafer et al., 2008).

After these filtering steps, we derived a final RBP target set, which included 171 distinct RBPs, 14,520 mRNA transcripts, and 2,865,656 target sites. On average, after considering RNA accessibility, a typical RBP was found to bind 2,077 genes (median 502 genes), and on average each gene was targeted by 24.5 RBPs (median 23 RBPs) (Table 1). We found that a few proteins such as HNRNPC and DDX3X had the most targets (11 RBPs had more than 8000 targets each). These RBPs are among the most studied proteins and have the greatest number of CLIP-Seq experiments done on them. On the other hand, a gene can be targeted by up to 92 RBPs (e.g., HNRNPA2B1 and MAFG, Table 2). The detailed data can be found in Supplementary Tables S3, S4.

Overlap Among RBP Binding Sites

As shown in Figure 1A, we separated pairs of trans factors into three scenarios: Overlapping, Neighboring, and Independent. As described above, in this study we did not consider alternative splicing or alternative 3′ UTRs; instead we choose the longest 3′ UTR for each human gene. We define two trans regulators having overlapping target sites if their target sites, either experimentally determined or predicted, overlap by at least 1 nucleotide. We experimented with alternative definition such as defining overlapping as distance on the mRNA being less than 5 nucleotides, such alterations had no effect on the final conclusions.

Figure 2A shows the distribution of the number of RBP binding sites per gene, after pooling binding sites from different RBPs for the same gene. On average, a human mRNA transcript has 197 RBP binding sites (median 85 binding sites). Figure 2B shows the number of distinct RBPs that target each mRNA. On average, a human mRNA transcript is targeted by 24.5 distinct RBPs (median 23 RBPs). We note that these statistics are based only on the 171 RBPs for which there are in vivo binding data; the number of RBPs that target each mRNA is likely higher.

FIGURE 2

For the 171 RBPs, our analysis show that they have a total of 2,865,656 potential binding sites against the entire human transcriptome. Among these binding sites, 2,579,262 (90%) overlap with another RBP binding site. In contrast, in our random simulations, approximately 2,236,200 (78.1%) of the binding sites overlap with another binding site. We next tested whether such high degree of overlap was caused by the very few RBPs that have many binding sites throughout the transcriptome. To test such possibility, we sequentially removed the RBPs that had the highest number of binding sites and re-calculated the frequency of overlap (Figure 3). The results showed that, even after removing the top 5 RBPs that have the most abundant target sites, on average, 86.3% of the RBPs still have overlapping sites with another RBP (70.1% in random simulations). This indicates that indeed many RBPs share potential overlapping binding sites, at least as determined by pull-down experiments. Similar observation had been made by other authors from the study of in vitro determined binding affinities of RBPs (Van Nostrand et al., 2016). There are two potential explanations for the observed significant overlap between the potential RBP binding sites. A mechanistic explanation is that RBPs preferably recognize RNA regions that have unique nucleotide composition, favorable structural conformation and accessibility. Many of these RBPs share homologous RNA binding domains (RBD), thus it is likely that they bind to similar target sequence motifs. Alternatively, such overlap could be the result of millions of years of adaptive evolution and selection.

FIGURE 3

Table 2 shows some of the genes that are bound by the greatest number of RBPs, also shown are the number of binding sites, and the number of binding sites that overlap with another binding site. The complete list of such data is shown in Supplementary Table S5. Table 3 lists the RBP pairs that have the highest number of overlapping binding sites. HNRNPC, ELAVL1 (HuR), TARDBP (TDP-43), and ATXN2 are among the RBPs that have the highest number of binding sites and target genes (see Table 1); they also have the highest number of overlapping binding sites and overlapping mRNA targets. Figure 4 shows the network of overlapping and neighboring RBPs, represented as a heatmap, detailed data are in Supplementary Table S6.

FIGURE 4

Figure S1 show that, at least as suggested by sequence data, there exists prevalent potential overlap among RBP binding sites. For any RBP-RBP pair, there is 80.1% chance (74% in randomized simulation) that they have overlapping binding sites on at least one mRNA target; and 94.6% chance that they have neighboring binding sites on at least one mRNA. There is 19% chance (15.2% in randomized simulation) that any pair of RBPs have overlapping binding sites on at least 100 mRNAs, and 38% chance that they have neighboring binding sites on at least 100 mRNAs.

Mapping MicroRNA Target Sites

On average, each microRNA is predicted to regulate 557 genes (median 433 genes), and each gene is predicted to be targeted by 16 microRNAs (median 9 microRNAs). Figure 5 shows the distribution of microRNA binding sites among the genes; more detailed information can be found in Supplementary Tables S7, S8. Among the total of 201,235 microRNA binding sites, 157,245 (78%) overlap with another microRNA binding site. There are several previously published studies on overlap and interaction between microRNA regulations including those from our own group (Li et al., 2014), therefore, we focused on interactions between RBPs and microRNAs in this study. Figure S2 shows the overlap between microRNA regulations, represented as a heatmap. Supplementary Table S9 lists the microRNA pairs that have the highest number of overlapping binding sites. Since these top microRNA pairs often belong to the same family and are highly similar in sequence, their binding sites are almost always deemed as overlapping, i.e., no neighboring binding sites between these pairs. The complete list is shown in Supplementary Table S10. It should be noted that human microRNAs consist of many families of homologous microRNAs that share similar seed sequences, e.g., hsa-miR-15a and has-miR-15b; these microRNAs would naturally have overlapping binding sites. We treated these homologous microRNAs as independent regulatory elements since they do not always have similar expression levels across cell types.

FIGURE 5

Overlap Between RBP and MicroRNA Binding Sites

We next investigated overlap between RBP and microRNA binding sites. As described above, our final dataset included a total of 2,865,656 RBP binding sites mediated by 171 RBPs, and 201,235 microRNA binding sites mediated by 321 microRNAs. On average, each mRNA transcript had 197 RBP binding sites (median 85) and 18 microRNA binding sites (median 9). Among the 2,865,656 RBP binding sites, 181,690 (6.3%) (average 243,260, or 8.5% in randomized simulations) overlapped with another microRNA binding site, while 2,361,587 (82.4%) had neighboring microRNA binding sites on the same mRNA transcript. It is surprising that in the randomized simulations more RBP binding sites overlap with another microRNA binding site; one possible explanation is that the real RBP binding sites tend to cluster together in the 3′ UTR thus the randomized RBP binding sites had more chance to overlap with randomized microRNA binding sites. Among the total of 201,235 microRNA binding sites, 111,300 (55.3%) (average 109,503, 54% in randomized simulations) overlapped with an RBP binding site, 11,371 (5.7%) had a neighboring RBP binding site. These results including the simulation results suggests potential prevalent interactions between RBPs and microRNAs.

Figure 6A shows the RBPs and microRNAs that have overlapping binding sites, while Figure 6B shows the RBPs and microRNAs that have neighboring binding sites. Table 4 lists the details of the RBP-microRNA pairs that have the highest number of overlapping or neighboring binding sites. The complete list is in the Supplementary Table S11.

FIGURE 6

For any pair of RBP and microRNA, there is 46.3% chance (average 42.2% in randomized simulation) that they have overlapping binding sites on at least one mRNA transcript, and 93.3% chance that they have neighboring binding sites on at least one mRNA. From Figure 6A and Table 4, we can see that several RBPs such as TARDBP (TDP-43), ATXN2, ELAVL1 (HuR) dominated among the RBP-miRNA pairs that had high number of overlapping binding sites. One potential reason is that these RBPs tend to have high number of binding sites and mRNA targets; indeed, the RBPs that have the highest number of binding sites also tend to have higher number of overlapping microRNA target sites (Pearson correlation 0.97, p-value = 1.36-E116). Also, since many microRNAs belong to the same family, they also have the same number of overlapping or neighboring binding sites with RBPs, e.g., hsa-miR-181a, hsa-miR-181b, hsa-miR-181c, hsa-miR-181d in Table 4.

TARDBP (TDP-43), the top RBP listed in Table 4, is known to bind to both DNA and RNA and can regulate RNA editing as well (Quinones-Valdez et al., 2019); it also interacts with Drosha and Dicer complexes and directly regulates the biogenesis of a subset of microRNAs (Kawahara and Mieda-Sato, 2012; Chen et al., 2018). A recent study has shown that TARDBP can bind to the 3′ UTR of pluripotency factors including Sox 2 and protect these transcripts from degradation mediated by microRNAs such as miR-21 (Modic et al., 2019). Interestingly, both TARDBP and ATXN2 are implicated in neurological disorders such as amyotrophic lateral sclerosis (ALS) (Sreedharan et al., 2008; Elden et al., 2010; Prasad et al., 2019; Yang et al., 2019); in addition, Ataxin-2 has been implicated in microRNA regulation (McCann et al., 2011; Ostrowski et al., 2017).

We next searched the ATtRACT database (Giudice et al., 2016), which has curated RBP binding motifs and positional weight matrices (PWM) determined from pull-down experiments or in vitro methods such as RNAcompete (Ray et al., 2013). For the RBPs shown in Table 4, we only found curated binding motifs for TARDBP (TDP-43) and ELAVL1 (HuR) in this database. TARDBP (TDP-43) has a canonical RNA Recognition Motif (RRM) and it has consensus binding motifs of 5′GAAUGG3′ or 5′GAAUGU3′ as determined by RNAcompete. Hsa-MiR-181 has a mature sequence of 5′AACAUUCAACGCUGUCGGUGAGU3′, which gives rise to seed region of 5′ACAUUCA3′. The complementary target site on the mRNA is 5′UGAAUGU3′, which strongly resembles the curated binding motifs of TARDBP. ELAVL1 (HuR) has consensus U-rich binding motifs. Hsa-miR-340-5p has a mature sequence of 5′UUAUAAAGCAAUGAGACUGAUU3′, which gives rise to seed region of 5′UAUAAAG3′. The complementary target site on the mRNA is 5′CUUUAUA″, which bears resemblance of the U-rich binding motifs of ELAVL1 (HuR). Therefore, the high prevalence of overlap between RBP-microRNA pairs can be explained by the similarity of their binding preferences.

Regression Analysis Between Trans Regulators Based on TCGA Gene Expression Data

The above sections only investigated the potential of interactions between trans factors based on the location of their binding sites on the 3′ UTR of mRNAs. Next, we set out to analyze gene expression data to determine whether there are indeed interactions between these overlapping or neighboring trans regulators. The most definitive evidence to validate interactions between two trans factors is to experimentally remove one trans factor or its target site, then ascertain whether the regulatory effect of the other trans factor is elevated (competition) or diminished (cooperation). Such approaches were previously taken in the study of competitive effect between microRNAs (van Dongen et al., 2008; Khan et al., 2009), and between selected microRNA and RBPs (Kim et al., 2009; Glorian et al., 2011). Since the goal of our study is to conduct a global analysis and identify potential interactions between trans factors, we need an effective statistical method for this purpose. Toward this goal, we adopted a simplified statistical approach (Eq. 1), in which x₁ and x₂ represent the expression level of two trans regulators (RBPs or miRNAs), a₁ represents the “regulatory strength” that y is regulated by x₁, a₂ represents the regulatory strength that y is regulated by x₂, a₃ is the interaction effect between x₁ and x₂, and a₄ is the random error term, representing regulatory effect by other factors. If indeed there are interactions between two trans factors, either competitive or cooperative in a specific cancer type, we are likely to see statistically significant interactions between these two trans factors.

Equation 1:

We downloaded the TCGA gene expression data from the Broad Institute FIREHOSE website (version 2016_01_28), and tested Eq. 1 on the following cancer types that have the greatest number of expression datasets: BRCA, THCA, LIHC, LUSC, LUAD, PRAD, KIRC. The number of gene expression datasets for these cancer types is listed in Supplementary Table S1.

For each cancer type, TCGA contains matched data for both cancer patients and normal controls. We conducted separate studies on three different gene expression datasets: tumor samples only, normal control samples only, and a union of tumor and normal control samples (results are summarized in Tables 5A–C, respectively). We first conducted regression analysis on RBP-microRNA pairs, i.e., x₁ represents expression of an RBP, x₂ represents expression level of a microRNA in Eq. 1. For each dataset and for each cancer type, we enumerated all the Overlapping, Neighboring, and Independent RBP-microRNA pairs, and conducted regression analysis on their expression profiles to identify interacting pairs. We used the R package p.adjust to calculate the p-values adjusted by Benjamini–Hochberg method (q-value < 0.05). If our hypothesis is correct, we expect to see more RBP-microRNA pairs that have statistically significant p-value on a₃ term; in other words, more interacting pairs among Overlapping factors, fewer among Neighboring factors and the fewest among Independent factors. To test such hypothesis, we used Kolmogorov–Smirnov (K–S) test to compare between these three groups, i.e., the Overlapping pairs, the Neighboring pairs and the Independent pairs, to see which group has collectively lower p-values, i.e., which group has more pairs of interacting trans factors.

Table 5A lists the analysis results from the tumor samples of the TCGA dataset. We found that, in most of the tumor types except for BRCA, collectively, the “Overlapping” pairs of RBP and microRNA had significantly lower p-values than “Independent” pairs (Column 3). It is likely that the breast cancer samples consist of subtypes and larger differences between tumor samples, which makes it less efficient in the regression analysis as outlined in Eq. 1. When comparing between “Neighboring” and “Independent” pairs (Column 4), in four cancer types (LUAD, LUSC, PRAD, THCA), Neighboring pairs have collectively lower p-values than Independent pairs, which is consistent with our expectations. In contrast, Column 2 shows that the Overlapping RBP-microRNA pairs have collectively lower p-values than Neighboring pairs only in three cancer types (KIRC, LIHC, LUAD). In our model (Figure 1A), both the Overlapping and Neighboring pairs can interact with each other, although in different mode. The Overlapping factors directly compete for the same target sequence in the mRNA transcript, while the Neighboring factors interact with each other in an indirect manner. Taken this altogether, the results in Table 5A is consistent with our hypothesis that Overlapping and Neighboring pairs of RBP and microRNA can interact with each other on regulating common gene targets.

Table 5B lists a similar comparison on the normal control samples from TCGA, while Table 5C lists the results from the pooled tumor and control samples, on pairs of RBPs and microRNAs. The results in Table 5C is mostly consistent with the results in Table 5A, while the results from normal control samples (Table 5B) reveals fewer interacting RBP-microRNA pairs as indicated by the modest p-values between Overlapping and Independent pairs, and between Neighboring and Independent pairs. It is not clear why the normal control samples are not as effective as cancer samples in revealing interacting RBP-microRNA pairs. One potential reason could be the relative homogeneity in gene expression among the normal control samples, which reduces the statistical power in the regression analysis (Eq. 1). We also note that such widespread interactions between RBPs and microRNAs are observed in some cancer types but not in others. It could be because the gene expression data from certain cancer types have limited statistical power, due to study design, or the homogeneity of the cohort. Regardless, as shown in Table 5A, we indeed see the definitive evidence of interacting trans acting pairs in the tumor samples.

Table 5 shows the results of regression between RBPs and microRNAs, we next conducted regression analysis between RBP-RBP pairs (Table 6) and between microRNA-microRNA pairs (Table 7). Consistent with our hypothesis, Table 6 shows that, collectively, there are significantly more regulatory interactions among Overlapping and Neighboring pairs of RBPs than Independent pairs (Column 3 and 4); the prevalence of interactions between Overlapping pairs and Neighboring pairs of RBPs are statistically indistinguishable (Column 2). Regression analysis of microRNA pairs from TCGA tumor samples are shown in Table 7, which shows Overlapping pairs collectively have more interactions than Neighboring and Independent pairs, and Neighboring pairs have more prevalent interactions than Independent pair.

In the following we focus our analysis to the interactions between RBP-microRNA pairs. We ask whether the statistically significant RBP-microRNA pairs and significant RBP-microRNA-mRNA trios that we observe in one cancer type are also found in other cancer types. If indeed this is the case, then it would be considered as additional evidence for the validity of such interactions. Indeed, many of such interacting pairs or trios are found in more than one cancer type. However, we need to remove one confounding factor, as many of these RBP-microRNA pairs are present more frequently than other pairs in the input data (Eq. 1). We conducted a hypergeometric Fisher’s exact test and calculated the p-value of these frequently occurring RBP-microRNA pairs in each cancer types. Supplementary Table S12 lists the top occurring RBP-microRNA pairs and their p-values (Fisher’s exact test), and Supplementary Table S13 lists the top significant RBP-microRNA-mRNA trios and their p-values and q-values after Benjamini–Hochberg correction. The complete list of significant RBP-microRNA pairs and RBP-microRNA-mRNA trios from each cancer type can be found in the Supplementary Tables S14, S15 (also see the Venn diagrams in Figure S3). The detail of RBP-RBP pairs and RBP-RBP-mRNA trios are in Supplementary Tables S16–S19 (also see the Venn diagram in Figure S4), and the detail of microRNA-microRNA pairs and microRNA-microRNA-mRNA trios are in Supplementary Tables S20–S23 (also see the Venn diagram in Supplementary Figure S5).

Discussion

In this study, we analyzed 171 RBPs for which there are high quality CLIP-Seq experimental data available. This represents only a small fraction of the human RBP repertoire (424 known and predicted RBPs (Ray et al., 2013) or 1,542 manually curated RBPs (Gerstberger et al., 2014). We also used the predicted microRNA target sites instead of experimentally validated microRNA target sites, which is still scarce at the moment. Despite such limitations, we think our combined two-prong analysis on both gene expression and sequence data is still robust to discover regulatory interactions between these trans regulators. There are two lines of evidence to support this. Many of the trans interactions we discovered based on regression analysis were previously reported in the literature, such as ELAVL1 (HuR) and miR-19 (Kim et al., 2009; Glorian et al., 2011), F6 and YTHDC1 (Kasowitz et al., 2018). In addition, many of these interactions between trans factors are discovered in multiple cancer types, after controlling for false discovery rate.

Many of the identified trans regulators such as HuR have a large number of targets in the cell, so there is a possibility that some of these predicted or observed interactions between trans regulators are indirect effects instead of direct effects manifested on the same mRNA. The only definitive way to validate these trans regulatory effects is to quantify the individual regulatory effect of these trans factors and ascertain whether the combined regulatory effect is additive or not, similar to the study of epistasis or genetic interaction between regulatory elements. With the accumulation of genetic variation and individual gene expression data from projects such as GTEx, it may be feasible in the near future to achieve such a goal.

We showed evidence that some of these RBPs and microRNAs have cooperative or competitive effect in gene regulation; some of these trans regulators or target genes are involved in cancer or neurological disorders. A potentially intriguing and useful application of such trans co-regulation effect is to manipulate one trans regulator to enhance or attenuate the effect of another trans regulator. Several intriguing examples had been reported in the literature, for example knocking down a microRNA to sensitize cancer cell’s response to drugs that target RBPs (Gabra and Salmena, 2017). We hope our systematic analysis will provide some insight on understanding the interaction and mechanism of trans factors in human diseases and on designing on effective therapeutic approaches.

References

• 1

  AgarwalV.BellG. W.NamJ. W.BartelD. P. (2015). Predicting effective microRNA target sites in mammalian mRNAs.eLife4:e05005.

  • Google Scholar
  • View reference in article

• 2

  Al-AhmadiW.Al-GhamdiM.Al-HajL.Al-SaifM.KhabarK. S. (2009). Alternative polyadenylation variants of the RNA binding protein, HuR: abundance, role of AU-rich elements and auto-regulation.Nucleic Acids Res.373612–3624. 10.1093/nar/gkp223

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  BarkerA.EpisM. R.PorterC. J.HopkinsB. R.WilceM. C.WilceJ. A.et al (2012). Sequence requirements for RNA binding by HuR and AUF1.J. Biochem.151423–437. 10.1093/jb/mvs010

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  CalinG. A.CroceC. M. (2006). MicroRNA signatures in human cancers.Nat. Rev. Cancer6857–866. 10.1038/nrc1997

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  CastanottoD.SakuraiK.LingemanR.LiH.ShivelyL.AagaardL.et al (2007). Combinatorial delivery of small interfering RNAs reduces RNAi efficacy by selective incorporation into RISC.Nucleic Acids Res.355154–5164. 10.1093/nar/gkm543

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  ChenX.FanZ.McGeeW.ChenM.KongR.WenP.et al (2018). TDP-43 regulates cancer-associated microRNAs.Protein Cell9848–866. 10.1007/s13238-017-0480-9

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  CorcoranD. L.GeorgievS.MukherjeeN.GottweinE.SkalskyR. L.KeeneJ. D.et al (2011). PARalyzer: definition of RNA binding sites from PAR-CLIP short-read sequence data.Genome Biol.12R79–R140.

  • Google Scholar
  • View reference in article

• 8

  DoyleF.TenenbaumS. A. (2014). Trans-regulation of RNA-binding protein motifs by microRNA.Front. Genet.5:79. 10.3389/fgene.2014.00079

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 9

  EbertM. S.NeilsonJ. R.SharpP. A. (2007). MicroRNA sponges: competitive inhibitors of small RNAs in mammalian cells.Nat. Methods4721–726. 10.1038/nmeth1079

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  EldenA. C.KimH. J.HartM. P.Chen-PlotkinA. S.JohnsonB. S.FangX.et al (2010). Ataxin-2 intermediate-length polyglutamine expansions are associated with increased risk for ALS.Nature4661069–1075. 10.1038/nature09320

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  GabraM. M.SalmenaL. (2017). microRNAs and acute myeloid leukemia chemoresistance: a mechanistic overview.Front. Oncol.7:255. 10.3389/fonc.2017.00255

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  GazzaraM. R.MalloryM. J.RoytenbergR.LindbergJ. P.JhaA.LynchK. W.et al (2017). Ancient antagonism between CELF and RBFOX families tunes mRNA splicing outcomes.Genome Res.271360–1370. 10.1101/gr.220517.117

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  GerstbergerS.HafnerM.TuschlT. (2014). A census of human RNA-binding proteins.Nat. Rev. Genet.15829–845. 10.1038/nrg3813

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  GiudiceG.Sanchez-CaboF.TorrojaC.Lara-PezziE. (2016). ATtRACT-a database of RNA-binding proteins and associated motifs.Database2016:baw035. 10.1093/database/baw035

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  GlorianV.MaillotG.PolesS.IacovoniJ. S.FavreG.VagnerS. (2011). HuR-dependent loading of miRNA RISC to the mRNA encoding the Ras-related small GTPase RhoB controls its translation during UV-induced apoptosis.Cell Death Diff.181692–1701. 10.1038/cdd.2011.35

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  GunzburgM. J.SivakumaranA.PendiniN. R.YoonJ. H.GorospeM.WilceM. C.et al (2015). Cooperative interplay of let-7 mimic and HuR with MYC RNA.Cell Cycle142729–2733. 10.1080/15384101.2015.1069930

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  HuB.YangY. T.HuangY.ZhuY.LuZ. J. (2017). POSTAR: a platform for exploring post-transcriptional regulation coordinated by RNA-binding proteins.Nucleic Acids Res.45D104–D114.

  • Google Scholar
  • View reference in article

• 18

  HuaZ.LvQ.YeW.WongC. K.CaiG.GuD.et al (2006). MiRNA-directed regulation of VEGF and other angiogenic factors under hypoxia.PLoS One1:e116. 10.1371/journal.pone.0000116

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  KaragkouniD.ParaskevopoulouM. D.ChatzopoulosS.VlachosI. S.TastsoglouS.KanellosI.et al (2018). DIANA-TarBase v8: a decade-long collection of experimentally supported miRNA-gene interactions.Nucleic Acids Res.46D239–D245.

  • Google Scholar
  • View reference in article

• 20

  KasowitzS. D.MaJ.AndersonS. J.LeuN. A.XuY.GregoryB. D.et al (2018). Nuclear m6A reader YTHDC1 regulates alternative polyadenylation and splicing during mouse oocyte development.PLoS Genet.14:e1007412. 10.1371/journal.pgen.1007412

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  KawaharaY.Mieda-SatoA. (2012). TDP-43 promotes microRNA biogenesis as a component of the Drosha and Dicer complexes.Proc. Natl. Acad. Sci. U.S.A.1093347–3352. 10.1073/pnas.1112427109

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  KhanA. A.BetelD.MillerM. L.SanderC.LeslieC. S.MarksD. S. (2009). Transfection of small RNAs globally perturbs gene regulation by endogenous microRNAs.Nat. Biotechnol.27549–555. 10.1038/nbt.1543

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 23

  KimH. H.KuwanoY.SrikantanS.LeeE. K.MartindaleJ. L.GorospeM. (2009). HuR recruits let-7/RISC to repress c-Myc expression.Genes Dev.231743–1748. 10.1101/gad.1812509

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  KunduP.FabianM. R.SonenbergN.BhattacharyyaS. N.FilipowiczW. (2012). HuR protein attenuates miRNA-mediated repression by promoting miRISC dissociation from the target RNA.Nucleic Acids Res.405088–5100. 10.1093/nar/gks148

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  LiY.LiangC.WongK. C.JinK.ZhangZ. (2014). Inferring probabilistic miRNA-mRNA interaction signatures in cancers: a role-switch approach.Nucleic Acids Res.42:e76. 10.1093/nar/gku182

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  LorenzR.BernhartS. H.Honer Zu SiederdissenC.TaferH.FlammC.StadlerP. F.et al (2011). ViennaRNA Package 2.0.Algorithms Mol. Biol.6:26.

  • Google Scholar
  • View reference in article

• 27

  LukongK. E.ChangK.-W.KhandjianE. W.RichardS. (2008). RNA-binding proteins in human genetic disease.Trends Genet.24416–425. 10.1016/j.tig.2008.05.004

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  McCannC.HolohanE. E.DasS.DervanA.LarkinA.LeeJ. A.et al (2011). The Ataxin-2 protein is required for microRNA function and synapse-specific long-term olfactory habituation.Proc. Natl. Acad. Sci. U.S.A.108E655–E662.

  • Google Scholar

• 29

  MihailovichM.BremangM.SpadottoV.MusianiD.VitaleE.VaranoG.et al (2015). miR-17-92 fine-tunes MYC expression and function to ensure optimal B cell lymphoma growth.Nat. Commun.6:8725.

  • Google Scholar

• 30

  ModicM.GroschM.RotG.SchirgeS.LepkoT.YamazakiT.et al (2019). Cross-regulation between TDP-43 and paraspeckles promotes pluripotency-differentiation transition.Mol. Cell.74951.e13–965.e13.

  • Google Scholar

• 31

  OstrowskiL. A.HallA. C.MekhailK. (2017). Ataxin-2: from RNA control to human health and disease.Genes8:157. 10.3390/genes8060157

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 32

  PrasadA.BharathiV.SivalingamV.GirdharA.PatelB. K. (2019). Molecular mechanisms of TDP-43 misfolding and pathology in amyotrophic lateral sclerosis.Front. Mol. Neurosci.12:25. 10.3389/fnmol.2019.00025

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 33

  Quinones-ValdezG.TranS. S.JunH. I.BahnJ. H.YangE. W.ZhanL.et al (2019). Regulation of RNA editing by RNA-binding proteins in human cells.Commun. Biol.2:19.

  • Google Scholar

• 34

  RayD.KazanH.CookK. B.WeirauchM. T.NajafabadiH. S.LiX.et al (2013). A compendium of RNA-binding motifs for decoding gene regulation.Nature499172–177.

  • Google Scholar

• 35

  SalmenaL.PolisenoL.TayY.KatsL.PandolfiP. P. (2011). A ceRNA hypothesis: the rosetta stone of a hidden RNA language?Cell146353–358. 10.1016/j.cell.2011.07.014

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 36

  SreedharanJ.BlairI. P.TripathiV. B.HuX.VanceC.RogeljB.et al (2008). TDP-43 mutations in familial and sporadic amyotrophic lateral sclerosis.Science3191668–1672.

  • Google Scholar

• 37

  SrikantanS.TominagaK.GorospeM. (2012). Functional interplay between RNA-binding protein HuR and microRNAs.Curr. Protein Pept. Sci.13372–379. 10.2174/138920312801619394

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 38

  SumazinP.YangX.ChiuH. S.ChungW. J.IyerA.Llobet-NavasD.et al (2011). An extensive microRNA-mediated network of RNA-RNA interactions regulates established oncogenic pathways in glioblastoma.Cell147370–381. 10.1016/j.cell.2011.09.041

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 39

  TaferH.AmeresS. L.ObernostererG.GebeshuberC. A.SchroederR.MartinezJ.et al (2008). The impact of target site accessibility on the design of effective siRNAs.Nat. Biotechnol.26578–583. 10.1038/nbt1404

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 40

  TayY.RinnJ.PandolfiP. P. (2014). The multilayered complexity of ceRNA crosstalk and competition.Nature505344–352. 10.1038/nature12986

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 41

  van DongenS.Abreu-GoodgerC.EnrightA. J. (2008). Detecting microRNA binding and siRNA off-target effects from expression data.Nat. Methods51023–1025. 10.1038/nmeth.1267

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 42

  Van NostrandE. L.PrattG. A.ShishkinA. A.Gelboin-BurkhartC.FangM. Y.SundararamanB.et al (2016). Robust transcriptome-wide discovery of RNA-binding protein binding sites with enhanced CLIP (eCLIP).Nat. Methods13508–514. 10.1038/nmeth.3810

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 43

  YangY.HallidayG. M.KiernanM. C.TanR. H. (2019). TDP-43 levels in the brain tissue of ALS cases with and without C9ORF72 or ATXN2 gene expansions.Neurology93e1748–e1755.

  • Google Scholar

• 44

  ZarnackK.KonigJ.TajnikM.MartincorenaI.EustermannS.StevantI.et al (2013). Direct competition between hnRNP C and U2AF65 protects the transcriptome from the exonization of Alu elements.Cell152453–466. 10.1016/j.cell.2012.12.023

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 45

  ZhangC.DarnellR. B. (2011). Mapping in vivo protein-RNA interactions at single-nucleotide resolution from HITS-CLIP data.Nat. Biotechnol.29607–614. 10.1038/nbt.1873

  • Pubmed Abstract
  • CrossRef
  • Google Scholar