### Consensus folding energy and covariance score in RNAalifold

The basic approach of RNAalifold [25] is to integrate covarying mutation into the thermodynamic model to predict consensus structures. First, covariance scores are computed for all pairs of columns to determine possible pairing positions in the consensus structure. Then, based on the MFE model, the minimum average folding energy is computed with dynamic programming. Assume the given alignment is denoted by \mathbb{A}, which contains *N* sequences \mathbb{A}=\{{s}^{1},{s}^{2},\dots ,{s}^{N}\}. Each sequence contains *L* symbols, including nucleotides and gaps, and {s}_{i}^{k} represents the *i*^{th} symbol (1≤*i*≤*L*) at the *k*^{th} (1≤*k*≤*N*) RNA sequence. The minimization of free energy is computed by using the following recursive functions:

\begin{array}{l}\phantom{\rule{1.9em}{0ex}}{F}_{i,j}=\text{min}({F}_{i+1,j},\underset{i<k\le j}{\text{min}}({C}_{i,k}+{F}_{k+1,j}))\\ \phantom{\rule{2.1em}{0ex}}{C}_{i,j}={\varphi}_{2}{\gamma}_{i,j}+\text{min}\phantom{\rule{3pt}{0ex}}\left\{\begin{array}{l}\sum _{{s}^{k}\in \mathbb{A}}H(i,j,{s}^{k})\\ \underset{i<p<q<j}{\text{min}}\left(\sum _{{s}^{k}\in \mathbb{A}}I(i,j,p,q,{s}^{k})+{C}_{p,q}\right)\\ \underset{i<p<j}{\text{min}}\left(F{M}_{i,p}+\mathit{\text{FM}}{1}_{p+1,j}+{M}_{a}\right)\end{array}\right.\phantom{\rule{3pt}{0ex}}\\ \phantom{\rule{1em}{0ex}}F{M}_{i,j}=\text{min}\phantom{\rule{3pt}{0ex}}\left\{\begin{array}{l}F{M}_{i+1,j}+{M}_{c}\\ \underset{i<p<j}{\text{min}}{C}_{i,p}+F{M}_{p+1,j}+{M}_{b}\\ \mathit{\text{FM}}{1}_{i,j}\end{array}\phantom{\rule{3pt}{0ex}}\right.\\ \phantom{\rule{3pt}{0ex}}\mathit{\text{FM}}{1}_{i,j}=\text{min}\phantom{\rule{3pt}{0ex}}\left(\mathit{\text{FM}}{1}_{i,j-1}+{M}_{c},{C}_{i,k}\phantom{\rule{3pt}{0ex}}\right)\end{array}

(1)

where *F*_{i,j},*C*_{i,j},*F* *M*_{i,j},*F* *M*1_{i,j} denote the minimum free energies for the region between *i*^{th} column and *j*^{th} column with unconstrained structure, with enclosed structure, with a multi-loop, and with a multi-loop containing a single branch, respectively. *H*(*i*,*j*,*s*) is the free energy for a hairpin loop enclosed by *s*_{
i
} and *s*_{
j
}, and *I*(*i*,*j*,*p*,*q*,*s*) is the free energy for an internal loop containing two base-pairs, one is between *s*_{
i
} and *s*_{
j
} and the other is between *s*_{
p
} and *s*_{
q
}. *M*_{
a
},*M*_{
c
} are penalties for closing bases and non-pairing bases in multi-loops. *M*_{
b
} is the bonus for branch bases in multi-loops.

The recursive functions were derived from the Turner’s model [35]. One major change made by RNAalifold for consensus folding is the usage of covariance score *γ*. It is not only a factor in the computing of free energy, but also determines the possible pairing columns in the alignment. Two parts, one is bonus and the other is penalty, are in this score. The first part of the covariance score is called the conservation score. For ({s}_{i}^{k},{s}_{j}^{k}) and ({s}_{i}^{l},{s}_{j}^{l}), three levels of confidence for pairing are assessed: base-pairs without mutation, base-pairs with one mutation, and base-pairs with two mutations. In the latest version of Vienna RNA package (v2.0) [36], the recursive function for computing conservation score is:

\begin{array}{l}{V}_{i,j}\phantom{\rule{0.3em}{0ex}}=\phantom{\rule{0.3em}{0ex}}\frac{1}{N}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\sum _{1\le k<l\le N}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\left\{\begin{array}{ll}\phantom{\rule{0.3em}{0ex}}h({s}_{i}^{k},{s}_{i}^{l})\phantom{\rule{0.3em}{0ex}}+h({s}_{j}^{k},{s}_{j}^{l})& \text{if}\phantom{\rule{1em}{0ex}}({s}_{i}^{k},{s}_{j}^{k})\phantom{\rule{0.3em}{0ex}}\in \phantom{\rule{0.3em}{0ex}}\mathbb{B}\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}({s}_{i}^{l},{s}_{j}^{l})\phantom{\rule{0.3em}{0ex}}\in \phantom{\rule{0.3em}{0ex}}\mathbb{B}\\ 0& \text{otherwise}\end{array}\right.\end{array}

(2)

where *h*(*x*,*y*) is the Hamming distance between base *x* and base *y*, and \mathbb{B} ={‘AU’, ‘UA’, ‘CG’, ‘GC’, ‘GU’, ‘UG’} is the set of all possible base-pairs. The second part is the penalty score *Q*_{i,j}, which deals with a pair of symbols that cannot form a base-pair:

{Q}_{i,j}=\sum _{1\le k\le N}\left\{\begin{array}{ll}0& \phantom{\rule{1em}{0ex}}\text{if}\phantom{\rule{3pt}{0ex}}({s}_{i}^{k},{s}_{j}^{k})\in \mathbb{B}\\ 0.25& \phantom{\rule{1em}{0ex}}\text{if}\phantom{\rule{3pt}{0ex}}{s}_{i}^{k}\phantom{\rule{3pt}{0ex}}\text{and}\phantom{\rule{3pt}{0ex}}{s}_{j}^{k}\text{are gaps}\\ 1& \phantom{\rule{6pt}{0ex}}\text{otherwise}\end{array}\right.

(3)

Overall, the covariance score is:

{\gamma}_{i,j}={V}_{i,j}-{\varphi}_{1}\times {Q}_{i,j}

(4)

where *ϕ*_{1}=*ϕ*_{2}=1. A threshold value *γ*_{
t
}=−2 is defined for *γ*_{i,j}. If *γ*_{i,j}>*γ*_{
t
}, *i*^{th} column and *j*^{th} column are considered to be pairing columns. In the final output, the minimum average folding energy, including the covariance score, is normalized by dividing *N*.

### Phylogenetic-based covarying mutation

RNAalifold incorporates covarying mutations into consensus folding to improve the detection of pairing columns. From Equation (2), it can be seen that RNAalifold counts the level of covariance by treating all sequences equally and try all possible combinations of base-pairs. In short, RNAalifold models the relationship of sequences as a complete graph. As a result, the specific evolutionary relationship among sequences in the phylogenetic history is ignored. Take the RNA structural alignment in Figure 1 as an example. The red and green columns achieve the same covariance score (2) in RNAalifold. However, as described in [37], the conservation evidence in Figure 1(c) is stronger than that in Figure 1(b) because at least two mutations occur at the green columns while only one is required to form the red ones.

PhyloRNAalifold models the relationship of aligned sequences as a tree by introducing the phylogenetic history of the alignment into the computation of covariance scores. The level of structural conservation is measured by the number of covarying mutations on the tree. Our assumption is that more covarying mutations on the tree mean stronger evidence of conservation. In addition, PhyloRNAalifold does not discard the original scoring scheme of RNAalifold, because experimental results showed this scheme can infer significant RNA structural aspects with high sensitivity and selectivity [38]. Assume *m*_{i,j} covarying mutations occur between *i*^{th} and *j*^{th} columns on the alignment \mathbb{A}’s phylogenetic tree and the number of base-pairs on those columns is *b*_{i,j}. The value of *m*_{i,j} depends on the size of the alignment. Since our approach focuses on improving the bonus part of the covariance scores, the number of covarying mutations is normalized with its upper bound: \frac{{m}_{i,j}}{{b}_{i,j}-1}. A new factor for the conservation score is proposed:

{\u220a}_{i,j}=1+\beta \times \frac{{m}_{i,j}}{{b}_{i,j}-1}

(5)

where *β* is the scale parameter for the normalized covarying mutation numbers. PhyloRNAalifold computes covariance scores with the following formula:

{\gamma}_{i,j}^{p}={\u220a}_{i,j}\times {V}_{i,j}-{\varphi}_{1}\times {Q}_{i,j}

(6)

All the other parameters and their default values in RNAalifold are retained. Due to the fact that {\gamma}_{i,j}^{p}\ge {\gamma}_{i,j} (*∊*_{i,j}≥1), two columns would be marked as pairing in PhyloRNAalifold if their covariance score in RNAalifold is greater than the threshold *γ*_{
t
} (the default value of *γ*_{
t
} is -2). Thus the advantage of PhyloRNAalifold is to import more potentially pairing positions with high mutation numbers.

### Computing the number of covarying mutations

Given a phylogenetic tree and labels at its leaves, the Fitch algorithm can optimize nucleotide assignment of the internal nodes to minimize the number of mutations [39]. If we model solving phylogeny as a maximum parsimony problem, this number can be taken as the actual number of mutations. The Fitch algorithm consists of a forward phase and a backward phase. In the forward phase, all possible labels at each internal node are inferred. In addition, the number of mutations is estimated during a bottom-up traversal. In the backward phase, a top-down pass is performed to find the optimal label at each internal node. Only the forward algorithm is applied to PhyloRNAalifold, since we do not need the exact labels at the internal nodes, but only the number of mutations on the tree. Without loss of generality, we require \mathbb{T} to be a rooted binary tree. *r* denotes the root of \mathbb{T} and *v*, *v*_{
l
}, *v*_{
r
} denote a node, left child of *v*, and right child of *v* respectively. *F*(*v*) is the set of possible labels at node *v*, and *c* *o* *s* *t*(*v*) is the number of mutations on the sub-tree which is rooted at *v*. Then the forward phase can be described with the following recursive functions:

\begin{array}{l}\begin{array}{l}\phantom{\rule{1.2em}{0ex}}F(v)=\left\{\begin{array}{ll}F({v}_{l})\cap F({v}_{r})& \phantom{\rule{1em}{0ex}}\text{if}\phantom{\rule{2.77626pt}{0ex}}F({v}_{l})\cap F({v}_{r})\ne \varnothing \\ F({v}_{l})\cup F({v}_{r})& \phantom{\rule{1em}{0ex}}\text{otherwise}\end{array}\right.\\ \phantom{\rule{-6em}{0ex}}\mathit{\text{cost}}(v)=\left\{\begin{array}{ll}\mathit{\text{cost}}({v}_{l})+\mathit{\text{cost}}({v}_{r})& \phantom{\rule{2.77626pt}{0ex}}\phantom{\rule{2.77626pt}{0ex}}\text{if}\phantom{\rule{2.77626pt}{0ex}}F({v}_{l})\cap F({v}_{r})\ne \varnothing \\ \mathit{\text{cost}}({v}_{l})+\mathit{\text{cost}}({v}_{r})+1& \phantom{\rule{2.77626pt}{0ex}}\phantom{\rule{2.77626pt}{0ex}}\text{otherwise}\end{array}\right.\end{array}\end{array}

(7)

For each leaf, *F*(*v*) is a base at the corresponding sequence. After the computation is finished, *c* *o* *s* *t*(*r*) shows the minimum number of mutations on the phylogenetic tree. The optimization of this algorithm was proved in [40].

In Equation 5, the computation of *∊*_{i,j} does not depend on non-pairing bases. Therefore, in the revised Fitch algorithm non-pairing bases need not to be considered when the number of covarying mutations is computed. We changed the original Fitch algorithm in two ways: (1) at any leaf node, if ({s}_{i}^{k},{s}_{j}^{k})\notin \mathbb{B}, set ({s}_{i}^{k},{s}_{j}^{k})= (‘-’, ‘-’); (2) for one internal node *v*, if the bases at *v*_{
l
}(*v*_{
r
}) is (‘-’, ‘-’), *v* will obtain *F*(*v*_{
r
})(*F*(*v*_{
l
})) as its label. One example of this algorithm is shown in Figure 1(d). The revised Fitch algorithm can be described by using the following functions:

\begin{array}{c}\begin{array}{l}\phantom{\rule{1.2em}{0ex}}F(v)=\left\{\begin{array}{ll}F({v}_{l})\cap F({v}_{r})& \phantom{\rule{1em}{0ex}}\text{if}\phantom{\rule{2.77626pt}{0ex}}F({v}_{l})\cap F({v}_{r})\ne \varnothing \phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}F({v}_{l})\ne \text{(\u2018-\u2019, \u2018-\u2019)}\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}F({v}_{r})\ne \text{(\u2018-\u2019, \u2018-\u2019)}\\ F({v}_{l})& \phantom{\rule{1em}{0ex}}\text{if}F\phantom{\rule{2.77626pt}{0ex}}({v}_{r})=\text{(\u2018-\u2019, \u2018-\u2019)}\\ F({v}_{r})& \phantom{\rule{1em}{0ex}}\text{if}\phantom{\rule{2.77626pt}{0ex}}F({v}_{l})=\text{(\u2018-\u2019, \u2018-\u2019)}\\ F({v}_{l})\cup F({v}_{r})& \phantom{\rule{1em}{0ex}}\text{otherwise}\end{array}\right.\\ \mathit{\text{cost}}(v)=\left\{\begin{array}{ll}\mathit{\text{cost}}({v}_{l})+\mathit{\text{cost}}({v}_{r})& \text{if}F({v}_{l})\cap \phantom{\rule{2.77626pt}{0ex}}F({v}_{r})\ne \varnothing \phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}F({v}_{l})\ne \text{(\u2018-\u2019, \u2018-\u2019)}\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}F({v}_{r})\ne \text{(\u2018-\u2019, \u2018-\u2019)}\\ \mathit{\text{cost}}({v}_{l})& \text{if}F({v}_{r})=\text{(\u2018-\u2019, \u2018-\u2019)}\\ \mathit{\text{cost}}({v}_{r})& \text{if}\phantom{\rule{2.77626pt}{0ex}}F({v}_{l})=\text{(\u2018-\u2019, \u2018-\u2019)}\\ \mathit{\text{cost}}({v}_{l})+\mathit{\text{cost}}({v}_{r})+1& \text{otherwise}\end{array}\right.\end{array}\end{array}

(8)

It is easy to see that our algorithm is optimal, because it only excludes non-pairing bases from the computation of the original Fitch algorithm.

In PhyloRNAalifold, the tree structure is an input variable and the clients can use any phylogenetic tree construction algorithm to build it. The time complexity of the original RNAalifold algorithm is *O*(*m*×*n*^{2}+*n*^{3}) [41], where *n* is the length of the alignment and *m* is the number of sequences in the alignment. The extra computation in PhyloRNAalifold is caused by the revised Fitch algorithm, whose time complexity ranges from *O*(log*m*) to *O*(*m*). In addition, PhyloRNAalifold needs to compute *∊*_{i,j} for each pair of columns in the alignment. Thus the overall time consumption of the revised Fitch algorithm is *O*(log*m*×*n*^{2}) or *O*(*m*×*n*^{2}). Neither of them increases the time complexity of RNAalifold.