CopomuS—Ranking Compensatory Mutations to Guide RNA-RNA Interaction Verification Experiments

In silico RNA-RNA interaction prediction is widely applied to identify putative interaction partners and to assess interaction details in base pair resolution. To verify specific interactions, in vitro evidence can be obtained via compensatory mutation experiments. Unfortunately, the selection of compensatory mutations is non-trivial and typically based on subjective ad hoc decisions. To support the decision process, we introduce our COmPensatOry MUtation Selector CopomuS. CopomuS evaluates the effects of mutations on RNA-RNA interaction formation using a set of objective criteria, and outputs a reliable ranking of compensatory mutation candidates. For RNA-RNA interaction assessment, the state-of-the-art IntaRNA prediction tool is applied. We investigate characteristics of successfully verified RNA-RNA interactions from the literature, which guided the design of CopomuS. Finally, we evaluate its performance based on experimentally validated compensatory mutations of prokaryotic sRNAs and their target mRNAs. CopomuS predictions highly agree with known results, making it a valuable tool to support the design of verification experiments for RNA-RNA interactions. It is part of the IntaRNA package and available as stand-alone webserver for ad hoc application.


Introduction
Many non-coding (nc)RNAs, like bacterial small (s)RNAs or eukayotic micro (mi)RNAs, perform their regulatory functions via direct RNA-RNA interaction (RRI) with their target RNAs [1]. This allows the use of in silico RRI prediction approaches to identify potential targets [2,3], which even provide interaction models in base pair resolution. To verify such predictions in vitro, mutation or deletion experiments are conducted [3,4]. Within deletion experiments, whole (potentially interacting) subsequences are removed from one or both RNAs and the in vitro measured interaction potential is compared to results using wildtype sequences. Since deletion of subsequences can have strong side effects, more sophisticated verification experiments are based on compensatory mutations (CoMs). To this end, both potentially interacting RNAs are mutated, such that wildtype-mutant combinations have a reduced base pairing potential within the interaction site that is regained in mutant-only interactions. Often, only a single base pair is mutated, but multiple concurrent CoMs are also used in literature [5]. A depiction of the setup is given in Figure 1. The in vitro RRI potential can be assessed via various experimental protocols [1], e.g., GFP reporter systems [3]. An RRI is considered verified via CoM if the in vitro RRI signal of the wildtype interaction is lost after mutating one sequence, and recovered when both RNAs are mutated (see Figure 1).

A' B A B
A' B'

A B'
A mutated B mutated A+B mutated A+B wildtype in vitro RRI signal strongly reduced signal recovery base pairing potential with wildtype lost The success of a CoM-based RRI verification strongly depends on the mutation selection. If the base-pair-breaking mutation is not weakening the RRI strong enough, no effect is detectable even if the functional RRI is mutated. On the other hand, if the mutation dramatically changes the intra-molecular structure formation of the respective RNA (e.g. when mutating multiple nucleotides at once), reduced RRI signal might be caused by an inaccessibility of the interaction site rather than the loss of base pairing potential. Thus, the design of potent CoMs is typically done manually based on the personal experience of the experimenter.
Here, we introduce CopomuS, a compensatory mutation selector to support this decision process. CopomuS in silico investigates and compares the effect of CoMs on RRI formation to rank candidate CoMs by their verification potential. To this end, IntaRNA [6,7], a state-of-the-art RRI prediction tool [8], is applied. For each CoM, RRI characteristics of the 4 wildtype-mutant sequence combinations are provided. This annotated list of CoMs is sorted by decreasing verification potential (as assessed by CopomuS) and enables an objective subsequent pruning by the experimenter given his/her expert knowledge. That way, CopomuS very much simplifies CoM selection and reliefs the experimenter of the otherwise necessary manual RRI prediction and comparison. Besides the saving of time, users can easily test the effect of different RRI and CoM constraints on the mutations of interest and the results are simple to reproduce. Finally, CopomuS translates the hypothesis behind compensatory mutation experiments into the respective theoretical model. That is, it implements a meta-strategy that checks for and ensures the intended stability difference pattern of wildtype versus mutant combinations.

Compensatory Mutations From Literature
To develop and evaluate CopomuS, we first extracted CoMs of successfully verified RRIs from literature, focusing on sRNA-target RRIs where sRNA binding inhibits translation of the target mRNA. The data was extracted from the benchmark set introduced in [9] and comprises both single-nucleotide CoMs (Table A1; 31 RNA pairs) as well as CoMs involving multiple concurrently mutated nucleotides (Table A2; 28 RNA pairs). To model the workflow based on an sRNA target prediction e.g., following [2,3], we use a genomic context around the start codon for each target. All sequences are provided in Tables A3-A6.

CopomuS Workflow
The central assumptions of CopomuS are that (1) the (main) regulatory RRI between the RNAs is defined by a single interaction site not interrupted by intra-molecular base pairing and (2) IntaRNA's prediction model correctly identifies the important parts of this regulatory RRI.
Only under Assumption-1 we can expect a single CoM to sufficiently alter the RRI potential between the RNAs, since the loss of one site can not be compensated by a concurrently formed second site. A classic example for such a multi-site RRI is the OxyS-fhlA interaction [10]. Furthermore, only under Assumption-2 CopomuS is expected to provide reasonable CoM candidates, since it generates and evaluates them based on the most stable RRIs predicted by IntaRNA.

CoM Generation
Following Assumption-2, CopomuS first computes the most stable RRI that can be formed by the wildtype RNAs using IntaRNA. The (minimal) free energy (MFE) estimate of the RRI computed by IntaRNA provides a stability proxy to assess the RRI potential. That is, the lower the energy the more stable the formed RRI is. Furthermore, the lowest energy RRI is the most likely when assuming a Boltzmann distribution of the energies [11]. The Nearest-Neighbor-model-based energy estimates incorporate both the stability of the inter-molecular base pairing (hybridization energy) and penalty (ED) terms to reflect the accessibility of the interacting subsequences [7]. More precisely, the latter describes the energy needed to break all intra-molecular base pairs that are formed by the RRI's subsequences [12].
CopomuS supports two modes to select inter-molecular base pairs for CoM candidate generation. Per default, all base pairs of the MFE RRI are considered. In addition, one can extend this set with all base pairs of suboptimal RRIs that can be formed by the subsequences covered by the MFE RRI. This softens the dependency on the accuracy of the reported MFE RRI base-pair pattern if alternative patterns within the same site are possible.
Next, identified base pairs are further filtered given user-defined constraints. For instance, specific base pair types, like AU or GU base pairs, can be filtered. Furthermore, lonely base pairs that can not stack with another base pair on either side can be excluded, since they typically provide only low stability contributions. Similarly, base pairs at putative ends of inter-molecular helices can be removed, since they can only form stackings on one side.
Finally, given the set of base pairs to be considered for mutation, the user can define whether all possible CoMs per base pair are to be considered (i.e., 3 non-compatible mutation alternatives for a given base pair like e.g., GU, GC, AU for UA (Note, CG is omitted since wildtype U can form a base pair with mutant G, which needs to be prevented)) or if the respective CoM candidate is only the 'nucleotide flip' of the base pair (as often done in literature), i.e., mutating a GC into a CG base pair.
At the end of the CoM generation, CopomuS outputs a list of CoM candidates to be evaluated and ranked.

CoM Characteristics and Ranking
Each CoM candidate is evaluated, using IntaRNA MFE-RRI predictions for all 4 sequence combinations for the current CoM, i.e., wildtype-only, wildtype-mutant as well as mutant-only combinations (see Figure 2). Depending on which ranking the user has specified, various RRI characteristics are aggregated for each combination, e.g., RRI stability in terms of MFE, RRI position, base pairs, or the accessibility of the mutated positions.
CopomuS implements a hierarchical CoM ranking. To this end, different classification and sorting functions are provided that can be sequentially combined to select for specific CoM characteristics. The most important classifier mfeCover checks whether the CoM is also covered by the mutant-only RRI. If not, the mutant-only RRI was found in a different location or no stable RRI was found at all. Next is the E classifier, which evaluates the desired reduction in RRI potential of wildtype-mutant combinations compared to wildtype-only or mutant-only RRIs. As discussed, the RRI MFE provides a proxy to assess interaction stability. Thus, CopomuS demands that both wildtype-only (ww) and mutant-only (mm) MFE are below zero and both wildtype-mutant MFEs are worse (higher), using two thresholds α and β, respectively. More formally, MFE(ww) + α < min(MFE(mw),MFE(wm)) and MFE(mm) + β < min(MFE(mw),MFE(wm)) have to be satisfied. This way, CopomuS can identify CoM candidates that have a high chance to show the expected in vitro RRI signal pattern from Figure 1 needed for RRI verification. The combination and hierarchy of classifier functions already defines a ranking of CoM subsets with decreasing level of constraint satisfaction.
Finally, each subset can be sorted via the minDeltaE function to get a final ranking of CoM candidates. minDeltaE favors CoMs with higher minimal MFE difference between wildtype-mutant combinations to wildtype-only or mutant-only combinations, i.e., it calculates and compares (min(MFE(mw),MFE(wm)) − max(MFE(ww),MFE(mm))). Therefore, the top-ranked CoM candidate will show the strongest RRI stability reduction for both wildtype-mutant combinations.

Availability
CopomuS is implemented in Python and part of the IntaRNA package version 3.2.0 or higher (freely available at https://github.com/BackofenLab/IntaRNA/). Due to its modular implementation, it can be easily expanded by further classification or sorting functions. Given its flexible command line interface and the provided CSV-based output, CopomuS can be easily integrated in automated pipelines and workflows. Its webserver (freely available at http://rna.informatik.uni-freiburg.de/ CopomuS/) for ad hoc usage is part of the Freiburg RNA tools framework [13] version 4.8.0 or higher.

Results
To evaluate the rationals behind CopomuS and to better understand CoM characteristics, we studied RRI characteristics of CoMs from literature. Subsequently, these CoMs were also used to benchmark CopomuS.

Statistics of CoMs from Literature
First, the distribution of mutation types was assessed, which is shown in Figure 3A. Within the single-nt CoM data set, mainly GC base pairs (29/31) have been mutated and all base pairs were 'flipped', i.e., mutated into their wildtype pairing partner. Multi-nt CoMs are also dominated by flipped mutations (52/86) and mutations of GC wildtype base pairs (42/86).
Next, we studied whether or not IntaRNA is able to correctly identify the CoM. We therefore identified the rank of the RRI predicted by IntaRNA (including suboptimals; sorted by energy) that includes at least one of the base pairs of the CoM. Respective statistics are provided in Figure 3B.  All base pairs of the single-nt CoMs can form base pair stackings (to at least one side), i.e., no lonely base pairs have been mutated. Most (24/31) single-nt CoMs are within helices or can stack to both sides.

Energy Profiling of CoMs from Literature
Next, we evaluated if IntaRNA MFE predictions are a useful proxy for RRI stability in order to identify highly potent CoM candidates based on the RRI stability differences between wildtype-/mutant-only (ww, mm) and wildtype-mutant combinations (wm, mw). We restricted the investigation to the 22 single-nt GC-mutating CoMs found in MFE RRIs. As a background model, we identified all additional GC base pairs (in total 207) from their respective MFE RRIs and treated them as flipped CGCG mutations. For each CoM and wildtype-mutant combination, we computed the respective MFE. Figure 4 visualizes the MFE distributions as well as MFE differences for both the CoMs known from literature as well as the background CoMs.
First, we compared the energy distributions of known CoMs and the background shown in Figure 4A. That is, for each wildtype-mutant combination (ww, wm, mw, mm), we compare the set of known CoM MFE values with the set of background CoM MFE values. This is done based on a sample t-test (known CoMs vs. background CoMs, two-tailed, unequal variance, unpaired). Respective p-values are reported in Figure 4A. Here, we only find the mutant-only combinations (mm) to be significantly different.
We also compared for each CoM the MFE of a given RNA combination (e.g., wm) with the respective value of the wildtype-only (ww) or mutant-only (mm) combination from the same data set. Beside visualizing the differences, Figure 4B also shows the p-values of respective paired sample t-tests (combination vs. ww|mm, two-tailed, unequal variance). With regard to wildtype-only (ww) MFE, both the known CoMs as well as the background CoMs show on average higher wildtype-mutant energies (i.e., wm-ww, mw-ww are positive). For the known CoMs, this also holds for mutant-only combinations (mm-ww), while the background CoMs do not show a significant difference between wildtype-only and mutant-only energies (mm-ww). That is, known CoMs show on average a stronger reduction of RRI potential of compensatory mutations (mm-ww) compared to the background. This difference also manifests in the energy differences when relating to the MFE of mutant-only (mm) combinations, i.e., the wildtype-mutant combinations of known CoMs (wm-mm, mw-mm) show in contrast to the background model on average no significant energy change. This is also reflected in the minDeltaE distributions, i.e., the mean minDeltaE difference is about zero for known CoMs while a slightly positive average value is observed for the background.

CopomuS Parameter Optimization
To find suitable values for the α/β parameters of the MFE-based CoM classification, we performed a parameter sweep in the range [1,5] given the observations from Figure 4. Here, we used the whole single-nt CoM data set and restricted the CoM candidate generation to the flipping of the MFEs' GC base pairs. We assessed whether or not the known single-nt CoM was among the valid candidates which fulfilled the MFE-difference constraints for the given α/β values and how many such candidates are left. A low number implies a clear separation of the known CoM from other candidates.
Results are depicted in Figure 5A. An α/β threshold of 1 prunes the candidate list already quite strongly. Increasing the thresholds results in a nearly linear loss of known CoMs among the remaining CoM candidates such that for value 2.5 already 10/31 known CoMs are considered invalid. On the other hand, no significant reduction of the average number of valid candidates can be observed for values between 1 and 2.5.
We thus performed another parameter screen to investigate non-uniform combinations of α and β, also summarized in Figure 5B. We observe that the effect of the energy-based candidate pruning is mainly governed by β, i.e., the energy difference to the mutant-only MFE. This is in accordance with the MFE difference statistics from Figure 4. As a result, we fix the default values of CopomuS's E function to α = 2 and β = 1.

CopomuS Benchmark
Finally, we evaluated CopomuS's ranking of known single-nt CoMs among its generated CoM candidates for different combinations of classification and sorting functions. We are interested in combinations for which many known CoMs fulfill all constraints with low (mean) ranking. Given our preliminary studies, we tested the effect of classification based on mfeCover and E (using α = 2 and β = 1) as well as sorting by minDeltaE. We restrict candidate generation to the flipping of GC and AU base pairs of the MFE RRI that are part of helices (no lonely base pairs or helix ends). Figure 6 summarizes the results. Classification based on mfeCover has only minor effects, while E-based pruning shows much stronger candidate set reductions without altering the number of known CoMs that are not fulfilling the constraints. minDeltaE sorting alone provides already very good results, which can be slightly enhanced when combined with E classification.

Statistics of CoMs from Literature
The high abundance of GC mutations is related to the higher base pairing strength of GC base pairs compared to AU or GU base pairs. Thus, preventing the formation of a GC base pair via mutation will have on average a higher impact on RRI stability and thus interaction potential.
Within multi-nt CoMs, GC mutations are not as dominant as in single-nt CoMs. Here, the set of mutations of the multi-nt CoM often aims at the core of a long inter-molecular helix to prevent their formation in wildtype-mutant combinations. Thus, less focus on GC base pairs is possible and even GU base pairs are part of the CoM. If possible, no GU base pair is introduced in the mutant (only 1/86), since GU base pairs have the lowest stability contributions.
Flipping wildtype base pairs in mutants is an easy way to ensure base pair incompatibility in wildtype-mutant combinations and enables (on average) a similar stability of the wildtype-only and mutant-only RRI. Thus, we chose flipping as the standard mode of CoM candidate generation.
The high rate of CoMs found within the IntaRNA MFE predictions supports our choice to build CopomuS's CoM selection based on IntaRNA RRI predictions. The prediction can fail e.g., if additional interaction partners, like RNA-binding proteins such as Hfq [14], are needed to guide or stabilize the RRI formation. In that case, IntaRNA's prediction model has to be guided with respective constraints, which can be incorporated in the CopomuS workflow using its optional IntaRNA parameter file. Examples for such constraints are location constraints where RRIs are assumed, structure probing data from e.g., SHAPE experiments [15] or explicit seed interaction information. Another reason for missing a known CoM is that only one base pair pattern of an RRI is investigated, but sometimes slightly different patterns with equal energy within the same boundaries are possible. Since IntaRNA reports only one pattern per RRI boundaries, a known CoM present in an alternative pattern is missed.
The high abundance of stacked base pairs within the data set motivated the introduced CoM candidate filter capabilities of CopomuS, to exclude lonely and helix-end base pairs from the candidate generation.

Energy Profiling of CoMs from Literature
The pattern of average MFE differences of known CoMs (inversely) follows the desired relation of in vitro RRI signals depicted in Figure 1. That is, wildtype-mutant RRIs show on average lower stability (higher MFEs) compared to wildtype-only or mutant-only interactions. This, in concert with the less prominent pattern within the background model, strongly supports our hypothesis that we can indeed use IntaRNA MFE predictions to select highly potent CoM candidates. Furthermore, the higher mean minDeltaE values of CoMs motivate our final ranking of CoM candidates in CopomuS to identify mutations that show the desired pattern most prominently.

CopomuS Benchmark
Nine of all known CoMs from literature are not among the remaining valid candidates for all measures and combinations. Seven of these are considered invalid by mfeCover due to a lack of MFE RRI coverage (compare Figure 3B) and two are at helix ends and thus filtered during candidate generation.
As expected, the combination of mfeCover with the E classifier impacts the ranking less than combining mfeCover with minDeltaE; since the latter provides a high-resolution sorting instead. Nevertheless, a combination of all three functions provides the best average rank of about 4 and furthermore instantiates all rationales underlying CopomuS. Thus, we make the combination of mfeCover and E with a final sorting by minDeltaE the default ranking of CopomuS.

Conclusions
CopomuS implements an automated, objective evaluation strategy to identify compensatory mutations (CoMs) based on IntaRNA-based RRI stability analyses. That is, top-ranked CoMs show an in silico RRI stability pattern that follows the desired pattern of in vitro RRI potentials. The required scoring functions were derived from characteristics of CoMs used in successful RRI verification experiments known from literature. We could show that the introduced measures efficiently reduce the set of CoM candidates, while the known CoM was found on average within the top-5 candidates. That way, experimenters can easily and reproducibly pick promising candidates from the provided list without the need of time consuming manual RRI prediction and comparison. Therefore, we consider CopomuS a valuable tool to guide the difficult design of highly promising CoM-based RRI verification experiments.
Currently, CopomuS supports only single-nt CoM generation and evaluation. While potent multi-nt CoMs could be derived from top-ranked single-nt CoMs that can stack, no correct ranking of their combination can be done based on the current output. We therefore work on the extension of CopomuS to multi-nt CoM generation and testing, to also provide reliable selection criteria for this setup. Furthermore, we would like to evaluate CopomuS's ranking order, which is currently not possible due to a lack of experimental data. That is, published CoM data is biased towards "valid" CoMs that were successfully used to verify an RRI. Non-successful CoMs or CoMs with low experimental effect are typically not published. Eventually, an extensive CoM feature and ranking evaluation would require in vitro measurements of multiple CoMs for the same RRI within the same experimental setup to be comparable.

Abbreviations
The following abbreviations are used in this manuscript:

RRI
RNA-RNA interaction CoM compensatory mutation MFE minimum free energy nt nucleotide ww, wm, mw, mm sRNA-mRNA combinations of (w)ildtype and (m)utated variants AUAU, AUCG, ... CoM classes where the first two letters encode the wildtype and the last two letters the mutated nts, e.g., AUCG encodes the mutation of either an AU or UA RRI base pair into either CG or GC Appendix A. Data Table A1. Single-nucleotide CoMs from literature used for this study.