Title: RNA SECONDARY STRUCTURE PREDICTION BY LEARNING UNROLLED ALGORITHMS

Abstract: In this paper, we propose an end-to-end deep learning model, called E2Efold, for RNA secondary structure prediction which can effectively take into account the inherent constraints in the problem. The key idea of E2Efold is to directly predict the RNA base-pairing matrix, and use an unrolled algorithm for constrained programming as the template for deep architectures to enforce constraints. With comprehensive experiments on benchmark datasets, we demonstrate the superior performance of E2Efold: it predicts significantly better structures compared to previous SOTA (especially for pseudoknotted structures), while being as efficient as the fastest algorithms in terms of inference time.

1. INTRODUCTION G G G:

Ribonucleic acid (RNA) is a molecule playing essential roles in numerous cellular processes and regulating expression of genes (Crick, 1970) . It consists of an ordered sequence of nucleotides, with each nucleotide containing one of four bases: Adenine (A), Guanine (G), Cytosine (C) and Uracile (U). This sequence of bases can be represented as

x := (x 1 , . . . , x L ) where x i ∈ {A, G, C, U }, which is known as the primary structure of RNA. The bases can bond with one another to form a set of base-pairs, which defines the secondary structure. A secondary structure can be represented by a binary matrix A * where A * ij = 1 if the i, j-th bases are paired (Fig 1) . Discovering the secondary structure of RNA is important for understanding functions of RNA since the structure essentially affects the interaction and reaction between RNA and other cellular components. Although secondary structure can be determined by experimental assays (e.g. X-ray diffraction), it is slow, expensive and technically challenging. Therefore, computational prediction of RNA secondary structure becomes an important task in RNA research and is useful in many applications such as drug design (Iorns et al., 2007) . Research on computational prediction of RNA secondary structure from knowledge of primary structure has been carried out for decades. Most existing methods assume the secondary structure is a result of energy minimization, i.e., A * = arg min A E x (A). The energy function is either estimated by physics-based thermodynamic experiments (Lorenz et al., 2011; Bellaousov et al., 2013; Markham & Zuker, 2008) or learned from data (Do et al., 2006) . These approaches are faced with a common problem that the search space of all valid secondary structures is exponentially-large with respect to the length L of the sequence. To make the minimization tractable, it is often assumed the base-pairing has a nested structure (Fig 2 left ), and the energy function factorizes pairwisely. With this assumption, dynamic programming (DP) based algorithms can iteratively find the optimal structure for subsequences and thus consider an enormous number of structures in time O(L 3 ).

Although DP-based algorithms have dominated RNA structure prediction, it is notable that they restrict the search space to nested structures, which excludes some valid yet biologically important RNA secondary structures that contain 'pseudoknots', i.e., elements with at least two non-nested base-pairs (Fig 2 right ). Pseudoknots make up roughly 1.4% of base-pairs (Mathews & Turner, 2006) , and are overrepresented in functionally important regions (Hajdin et al., 2013; Staple & Butcher, 2005) . Furthermore, pseudoknots are present in around 40% of the RNAs. They also assist folding into 3D structures (Fechter et al., 2001) and thus should not be ignored. To predict RNA structures with pseudoknots, energy-based methods need to run more computationally intensive algorithms to decode the structures.

In summary, in the presence of more complex structured output (i.e., pseudoknots), it is challenging for energy-based approaches to simultaneously take into account the complex constraints while being efficient. In this paper, we adopt a different viewpoint by assuming that the secondary structure is the output of a feed-forward function, i.e., A * = F θ (x), and propose to learn θ from data in an end-to-end fashion. It avoids the second minimization step needed in energy function based approach, and does not require the output structure to be nested. Furthermore, the feed-forward model can be fitted by directly optimizing the loss that one is interested in.

Despite the above advantages of using a feed-forward model, the architecture design is challenging. To be more concrete, in the RNA case, F θ is difficult to design for the following reasons:

(i) RNA secondary structure needs to obey certain hard constraints (see details in Section 3), which means certain kinds of pairings cannot occur at all (Steeg, 1993) . Ideally, the output of F θ needs to satisfy these constraints. (ii) The number of RNA data points is limited, so we cannot expect that a naive fully connected network can learn the predictive information and constraints directly from data. Thus, inductive biases need to be encoded into the network architecture. (iii) One may take a two-step approach, where a post-processing step can be carried out to enforce the constraints when F θ predicts an invalid structure. However, in this design, the deep network trained in the first stage is unaware of the post-processing stage, making less effective use of the potential prior knowledge encoded in the constraints.

All Binary Structures:

Output Space of E2Efold In this paper, we present an end-to-end deep learning solution which integrates the two stages. The first part of the architecture is a transformer-based deep model called Deep Score Network which represents sequence information useful for structure prediction. The second part is a multilayer network called Post-Processing Network which gradually enforces the constraints and restrict the output space. It is designed based on an unrolled algorithm for solving a constrained optimization. These two networks are coupled together and learned jointly in an end-to-end fashion. Therefore, we call our model E2Efold.

By using an unrolled algorithm as the inductive bias to design Post-Processing Network, the output space of E2Efold is constrained (illustrated in Fig 3 ), which makes it easier to learn a good model in the case of limited data and also reduces the overfitting issue. Yet, the constraints encoded in E2Efold are flexible enough such that pseudoknots are not excluded. In summary, E2Efold strikes a nice balance between model biases for learning and expressiveness for valid RNA structures.

We conduct extensive experiments to compare E2Efold with state-of-the-art (SOTA) methods on several RNA benchmark datasets, showing superior performance of E2Efold including:

• being able to predict valid RNA secondary structures including pseudoknots;

• running as efficient as the fastest algorithm in terms of inference time;

• producing structures that are visually close to the true structure;

• better than previous SOTA in terms of F1 score, precision and recall.

Although in this paper we focus on RNA secondary structure prediction, which presents an important and concrete problem where E2Efold leads to significant improvements, our method is generic and can be applied to other problems where constraints need to be enforced or prior knowledge is provided. We imagine that our design idea of learning unrolled algorithm to enforce constraints can also be transferred to problems such as protein folding and natural language understanding problems (e.g., building correspondence structure between different parts in a document).

2. RELATED WORK:

Classical RNA folding methods identify candidate structures for an RNA sequence energy minimization through DP and rely on thousands of experimentally-measured thermodynamic parameters. A few widely used methods such as RNAstructure (Bellaousov et al., 2013) , Vienna RNAfold (Lorenz et al., 2011) and UNAFold (Markham & Zuker, 2008) (Bellaousov & Mathews, 2010) have been made to predict structures with pseudoknots, but the predictive accuracy and efficiency still need to be improved.

Learning-based RNA folding methods such as ContraFold (Do et al., 2006) and ContextFold (Zakov et al., 2011) have been proposed for energy parameters estimation due to the increasing availability of known RNA structures, resulting in higher prediction accuracies, but these methods still rely on the above DP-based algorithms for energy minimization. A recent deep learning model, CDPfold (Zhang et al., 2019) , applied convolutional neural networks to predict base-pairings, but it adopts the dot-bracket representation for RNA secondary structure, which can not represent pseudoknotted structures. Moreover, it requires a DP-based post-processing step whose computational complexity is prohibitive for sequences longer than a few hundreds.

Learning with differentiable algorithms is a useful idea that inspires a series of works (Hershey et al., 2014; Belanger et al., 2017; Ingraham et al., 2018; Chen et al., 2018; Shrivastava et al., 2019) , which shared similar idea of using differentiable unrolled algorithms as a building block in neural architectures. Some models are also applied to structured prediction problems (Hershey et al., 2014; Pillutla et al., 2018; Ingraham et al., 2018) , but they did not consider the challenging RNA secondary structure problem or discuss how to properly incorporating constraints into the architecture.

OptNet (Amos & Kolter, 2017) integrates constraints by differentiating KKT conditions, but it has cubic complexity in the number of variables and constraints, which is prohibitive for the RNA case.

Dependency parsing in NLP is a different but related problem to RNA folding. It predicts the dependency between the words in a sentence. Similar to nested/non-nested structures, the corresponding terms in NLP are projective/non-projective parsing, where most works focus on the former and DP-based inference algorithms are commonly used (McDonald et al., 2005) . Deep learning models (Dozat & Manning, 2016; Kiperwasser & Goldberg, 2016) are proposed to proposed to score the dependency between words, which has a similar flavor to the Deep Score Network in our work.

3. RNA SECONDARY STRUCTURE PREDICTION PROBLEM:

In the RNA secondary structure prediction problem, the input is the ordered sequence of bases x = (x 1 , . . . , x L ) and the output is the RNA secondary structure represented by a matrix A * ∈ {0, 1} L×L . Hard constraints on the forming of an RNA secondary structure dictate that certain kinds of pairings cannot occur at all (Steeg, 1993) . Formally, these constraints are:

(i) Only three types of nucleotides combinations,

B := {AU, U A}∪ {GC, CG} ∪ {GU, U G}, can form base-pairs. ∀i, j, if x i x j / ∈ B, then A ij = 0. (ii) No sharp loops are allowed. ∀|i -j| < 4, A ij = 0.

(iii) There is no overlap of pairs, i.e., it is a matching. ∀i, L j=1 A ij ≤ 1. (i) and (ii) prevent pairing of certain base-pairs based on their types and relative locations. Incorporating these two constraints can help the model exclude lots of illegal pairs. (iii) is a global constraint among the entries of A * .

The space of all valid secondary structures contains all symmetric matrices A ∈ {0, 1} L×L that satisfy the above three constraints. This space is much smaller than the space of all binary matrices {0, 1} L×L . Therefore, if we could incorporate these constraints in our deep model, the reduced output space could help us train a better predictive model with less training data. We do this by using an unrolled algorithm as the inductive bias to design deep architecture.

4. E2EFOLD: DEEP LEARNING MODEL BASED ON UNROLLED ALGORITHM:

In the literature on feed-forward networks for structured prediction, most models are designed using traditional deep learning architectures. However, for RNA secondary structure prediction, directly using these architectures does not work well due to the limited amount of RNA data points and the hard constraints on forming an RNA secondary structure. These challenges motivate the design of our E2Efold deep model, which combines a Deep Score Network with a Post-Processing Network based on an unrolled algorithm for solving a constrained optimization problem.

4.1. DEEP SCORE NETWORK:

The first part of E2Efold is a Deep Score Network U θ (x) whose output is an L × L symmetric matrix. Each entry of this matrix, i.e., U θ (x) ij , indicates the score of nucleotides x i and x j being paired. The x input to the network here is the L × 4 dimensional one-hot embedding. The specific architecture of U θ is shown in 

P i = MLP ψ 1 (i), . . . , ψ (i), ψ +1 (i/L), . . . , ψ n (i/L)

, where {ψ j } is a set of n feature maps such as sin(•), poly(•), sigmoid(•), etc, and MLP(•) denotes multi-layer perceptions. Such position embedding idea has been used in natural language modeling such as BERT (Devlin et al., 2018) , but we adapted for RNA sequence representation;

• a stack of Transformer Encoders (Vaswani et al., 2017) which encode the sequence information and the global dependency between nucleotides; • a 2D Convolution layers (Wang et al., 2017) for outputting the pairwise scores. With the representation power of neural networks, the hope is that we can learn an informative U θ such that higher scoring entries in U θ (x) correspond well to actual paired bases in RNA structure. Once the score matrix U θ (x) is computed, a naive approach to use it is to choose an offset term s ∈ R (e.g., s = 0) and let A ij = 1 if U θ (x) ij > s. However, such entry-wise independent predictions of A ij may result in a matrix A that violates the constraints for a valid RNA secondary structure. Therefore, a postprocessing step is needed to make sure the predicted A is valid. This step could be carried out separately after U θ is learned. But such decoupling of base-pair scoring and post-processing for constraints may lead to sub-optimal results, where the errors in these two stages can not be considered together and tuned together. Instead, we will introduce a Post-Processing Network which can be trained end-to-end together with U θ to enforce the constraints.

4.2. POST-PROCESSING NETWORK:

The second part of E2Efold is a Post-Processing Network PP φ which is an unrolled and parameterized algorithm for solving a constrained optimization problem. We first present how we formulate the post-processing step as a constrained optimization problem and the algorithm for solving it. After that, we show how we use the algorithm as a template to design deep architecture PP φ .

4.2.1. POST-PROCESSING WITH CONSTRAINED OPTIMIZATION:

Formulation of constrained optimization. Given the scores predicted by U θ (x), we define the total score 1 2 i,j (U θ (x) ij -s)A ij as the objective to maximize, where s is an offset term. Clearly, without structure constraints, the optimal solution is to take A ij = 1 when U θ (x) ij > s. Intuitively, the objective measures the covariation between the entries in the scoring matrix and the A matrix. With constraints, the exact maximization becomes intractable. To make it tractable, we consider a convex relaxation of this discrete optimization to a continuous one by allowing A ij ∈ [0, 1]. Consequently, the solution space that we consider to optimize over is A(x) := A ∈ [0, 1] L×L | A is symmetric and satisfies constraints (i)-(iii) in Section 3 .

To further simplify the search space, we define a nonlinear transformation T on R L×L as T ( Â) :=

1 2 Â • Â + ( Â • Â) • M (x)

, where • denotes element-wise multiplication. Matrix M is defined as M (x) ij := 1 if x i x j ∈ B and also |i -j| ≥ 4, and M (x) ij := 0 otherwise. From this definition we can see that M (x) encodes both constraint (i) and (ii). With transformation T , the resulting matrix is non-negative, symmetric, and satisfies constraint (i) and (ii). Hence, by defining A := T ( Â), the solution space is simplified as

A(x) = {A = T ( Â) | Â ∈ R L×L , A1 ≤ 1}.

Finally, we introduce a 1 penalty term Â 1 := i,j | Âij | to make A sparse and formulate the post-processing step as: ( •, • denotes matrix inner product, i.e., sum of entry-wise multiplication)

EQUATION

The advantages of this formulation are that the variables Âij are free variables in R and there are only L inequality constraints A1 ≤ 1. This system of linear inequalities can be replaced by a set of nonlinear equalities relu(A1 -1) = 0 so that the constrained problem can be easily transformed into an unconstrained problem by introducing a Lagrange multiplier λ ∈ R L + :

EQUATION

Algorithm for solving it. We use a primal-dual method for solving Eq. 2 (derived in Appendix B).

In each iteration, Â and λ are updated alternatively by:

(primal) gradient step:

EQUATION

where

EQUATION

EQUATION

(dual) gradient step:

EQUATION

where α, β are step sizes and γ α , γ β are decaying coefficients. When it converges at T , an approximate solution Round A T = T ( ÂT ) is obtained. With this algorithm operated on the learned U θ (x), even if this step is disconnected to the training phase of U θ (x), the final prediction works much better than many other existing methods (as reported in Section 6). Next, we introduce how to couple this post-processing step with the training of U θ (x) to further improve the performance.

4.2.2. POST-PROCESSING NETWORK VIA AN UNROLLED ALGORITHM:

We design a Post-Processing Network, denoted by PP φ , based on the above algorithm. After it is defined, we can connect it with the deep score network U θ and train them jointly in an end-to-end fashion, so that the training phase of U θ (x) is aware of the post-processing step.

Algorithm 1: Post-Processing Network PP φ (U, M )

Parameters φ := {w, s, α, β, γ α , γ β , ρ} U ← softsign(U -s) • U Â0 ← softsign(U -s) • sigmoid(U ) A 0 ← T ( Â0 ); λ 0 ← w • relu(A 0 1 -1) For t = 0, . . . , T -1 do λ t+1 , A t+1 , Ât+1 = PPcell φ (U, M, λ t , A t , Ât , t) return {A t } T t=1 Algorithm 2: Neural Cell PPcell φ Function PPcell φ (U, M, λ, A, Â, t): G ← 1 2 U -(λ • softsign(A1 -1)) 1 Ȧ ← Â + α • γ α t • Â • M • (G + G ) Â ← relu(| Ȧ| -ρ • α • γ α t ) Â ← 1 -relu(1 -Â) [i.e.,min( Â, 1)] A ← T ( Â); λ ← λ+β•γ β t •relu(A1-1) return λ, A, Â

The specific computation graph of PP φ is given in Algorithm 1, whose main component is a recurrent cell which we call PPcell φ . The computation graph is almost the same as the iterative update from Eq. 3 to Eq. 6, except for several modifications:

• (learnable hyperparameters) The hyperparameters including step sizes α, β, decaying rate γ α , γ β , sparsity coefficient ρ and the offset term s are treated as learnable parameters in φ, so that there is no need to tune the hyperparameters by hand but automatically learn them from data instead. • (fixed # iterations) Instead of running the iterative updates until convergence, PPcell φ is applied recursively for T iterations where T is a manually fixed number. This is why in Fig 3 the output space of E2Efold is slightly larger than the true solution space. • (smoothed sign function) Resulted from the gradient of relu(•), the update step in Eq. 4 contains a sign(•) function. However, to push gradient through PP φ , we require a differentiable update step. Therefore, we use a smoothed sign function defined as softsign(c) := 1/(1 + exp(-kc)), where k is a temperature. • (clip Â) An additional step, Â ← min( Â, 1), is included to make the output A t at each iteration stay in the range [0, 1] L×L . This is useful for computing the loss over intermediate results {A t } T t=1 , for which we will explain more in Section 5. With these modifications, the Post-Processing Network PP φ is a tuning-free and differentiable unrolled algorithm with meaningful intermediate outputs. Combining it with the deep score network, the final deep model is

E2Efold : {A t } T t=1 = Post-Process Network PP φ ( U θ (x)

Deep Score Network , M (x)) .

(7)

5. END-TO-END TRAINING ALGORITHM:

Given a dataset D containing examples of input-output pairs (x, A * ), the training procedure of E2Efold is similar to standard gradient-based supervised learning. However, for RNA secondary structure prediction problems, commonly used metrics for evaluating predictive performances are F1 score, precision and recall, which are non-differentiable.

Differentiable F1 Loss. To directly optimize these metrics, we mimic true positive (TP), false positive (FP), true negative (TN) and false negative (FN) by defining continuous functions on [0, 1] L×L :

TP = A, A * , FP = A, 1 -A * , FN = 1 -A, A * , TN = 1 -A, 1 -A * .

Since F1 = 2TP/(2TP + FP + FN), we define a loss function to mimic the negative of F1 score as:

EQUATION

Assuming that ij A * ij = 0, this loss is well-defined and differentiable on [0, 1] L×L . Precision and recall losses can be defined in a similar way, but we optimize F1 score in this paper.

It is notable that this F1 loss takes advantages over other differentiable losses including 2 and cross-entropy losses, because there are much more negative samples (i.e. A ij = 0) than positive samples (i.e. A ij = 1). A hand-tuned weight is needed to balance them while using 2 or crossentropy losses, but F1 loss handles this issue automatically, which can be useful for a number of problems (Wang et al., 2016; Li et al., 2017) .

Overall Loss Function. As noted earlier, E2Efold outputs a matrix A t ∈ [0, 1] L×L in each iteration. This allows us to add auxiliary losses to regularize the intermediate results, guiding it to learn parameters which can generate a smooth solution trajectory. More specifically, we use an objective that depends on the entire trajectory of optimization:

EQUATION

where {A t } T t=1 = PP φ (U θ (x), M (x)) and γ ≤ 1 is a discounting factor. Empirically, we find it very useful to pre-train U θ using logistic regression loss. Also, it is helpful to add this additional loss to Eq. 9 as a regularization.

6. EXPERIMENTS:

We compare E2Efold with the SOTA and also the most commonly used methods in the RNA secondary structure prediction field on two benchmark datasets. It is revealed from the experimental results that E2Efold achieves 29.7% improvement in terms of F1 score on RNAstralign dataset and it infers the RNA secondary structure as fast as the most efficient algorithm (LinearFold) among existing ones. An ablation study is also conducted to show the necessity of pushing gradient through the post-processing step. The codes for reproducing the experimental results are released.1 2019), for a base pair (i, j), the following predictions are also considered as correct: (i + 1, j), (i -1, j), (i, j + 1), (i, j -1), so we also reported the metrics when one-position shift is allowed. As shown in Table 2 , traditional methods can achieve a F1 score ranging from 0.433 to 0.624, which is consistent with the performance reported with their original papers. The two learning-based methods, CONTRAfold and CDPfold, can outperform classical methods with reasonable margin on some criteria. E2Efold, on the other hand, significantly outperforms all previous methods across all criteria, with at least 20% improvement. Notice that, for almost all the other methods, the recall is usually higher than precision, while for E2Efold, the precision is higher than recall. That can be the result of incorporating constraints during neural network training. Fig 5 shows the distributions of F1 scores for each method. It suggests that E2Efold has consistently good performance.

To estimate the performance of E2Efold on long sequences, we also compute the F1 scores weighted by the length of sequences, such that the results are more dominated by longer sequences. Detailed results are given in Appendix D.3. 3 . It is understandable that the performances of classical methods which are not learningbased are consistent with that on RNAStralign. The performance of E2Efold, though is not as good as that on RNAStralign, is still better than all the other methods across different evaluation criteria. In addition, since the original ArchiveII dataset contains domain sequences (subsequences), we remove the domains and report the results in Appendix D.4, which are similar to results in Table 3 .

Inference Time Comparison. We record the running time of all algorithms for predicting RNA secondary structures on the RNAStralign test set, which is summarized in Table 4 . LinearFold is the most efficient among baselines because it uses beam pruning heuristic to accelerate DP. CDPfold, which achieves higher F1 score than other baselines, however, is extremely slow due to its DP post-processing step. Since we use a gradient-based algorithm which is simple to design the Post-Processing Network, E2Efold is fast. On GPU, E2Efold has similar inference time as LinearFold. Pseudoknot Prediction. Even though E2Efold does not exclude pseudoknots, it is not sure whether it actually generates pseudoknotted structures. Therefore, we pick all sequences containing pseudoknots and compute the averaged F1 score only on this set. Besides, we count the number of pseudoknotted sequences that are predicted as pseudoknotted and report this count as true positive (TP). Similarly we report TN, FP and FN in Table 5 along with the F1 score. Most tools exclude pseudoknots while RNAstructure is the most famous one that can predict pseudoknots, so we choose it for comparison. Ablation Study. To exam whether integrating the two stages by pushing gradient through the postprocess is necessary for performance of E2Efold, we conduct an ablation study (Table 6 ). We test the performance when the post-processing step is disconnected with the training of Deep Score Network U θ . We apply the post-processing step (i.e., for solving augmented Lagrangian) after U θ is learned (thus the notation "U θ + PP" in Table 6 ). Although "U θ + PP" performs decently well, with constraints incorporated into training, E2Efold still has significant advantages over it.

Discussion. To better estimate the performance of E2Efold on different RNA types, we include the per-family F1 scores in Appendix D.5. E2Efold performs significantly better than other methods in 16S rRNA, tRNA, 5S RNA, tmRNA, and telomerase. These results are from a single model. In the future, we can view it as multi-task learning and further improve the performance by learning multiple models for different RNA families and learning an additional classifier to predict which model to use for the input sequence.

A MORE DISCUSSION ON RELATED WORKS

Here we explain the difference between our approach and other works on unrolling optimization problems.

First, our view of incorporating constraints to reduce output space and to reduce sample complexity is novel. Previous works (Hershey et al., 2014; Belanger et al., 2017; Ingraham et al., 2018) did not discuss these aspects. The most related work which also integrates constraints is OptNet (Amos & Kolter, 2017) , but its very expensive and can not scale to the RNA problem. Therefore, our proposed approach is a simple and effective one.

Second, compared to (Chen et al., 2018; Shrivastava et al., 2019) , our approach has a different purpose of using the algorithm. Their goal is to learn a better algorithm, so they commonly make their architecture more flexible than the original algorithm for the room of improvement. However, we aim at enforcing constraints. To ensure that constraints are nicely incorporated, we keep the original structure of the algorithm and only make the hyperparameters learnable.

Finally, although all works consider end-to-end training, none of them can directly optimize the F1 score. We proposed a differentiable loss function to mimic the F1 score/precision/recall, which is effective and also very useful when negative samples are much fewer than positive samples (or the inverse).

B DERIVATION OF THE PROXIMAL GRADIENT STEP:

The maximization step in Eq. 1 can be written as the following minimization:

EQUATION

Consider the quadratic approximation of -f ( Â) centered at Ât :

EQUATION

and rewrite the optimization in Eq. 10 as

EQUATION

where

EQUATION

Next, we define proximal mapping as a function depending on α as follows:

EQUATION

Since we always use 

EQUATION

The last equation holds since Ât will remain symmetric in our algorithm if the initial Â0 is symmetric. Moreover, in the main text, α is replaced by α • γ t α .

C IMPLEMENTATION AND TRAINING DETAILS:

We used Pytorch to implement the whole package of E2Efold.

Deep Score Network. In the deep score network, we used a hyper-parameter, d, which was set as 10 in the final model, to control the model capacity. In the transformer encoder layers, we set the number of heads as 2, the dimension of the feed-forward network as 2048, the dropout rate as 0.1. As for the position encoding, we used 58 base functions to form the position feature map, which goes through a 3-layer fully-connected neural network (the number of hidden neurons is 5 * d) to generate the final position embedding, whose dimension is L by d. In the final output layer, the pairwise concatenation is carried out in the following way: Let X ∈ R L×3d be the input to the final output layers in Figure 4 (which is the concatenation of the sequence embedding and position embedding). The pairwise concatenation results in a tensor Y ∈ R L×L×6d defined as

EQUATION

where Y (i, j, :) ∈ R 6d , X(i, :) ∈ R 3d , and X(j, :) ∈ R 3d .

In the 2D convolution layers, the the channel of the feature map gradually change from 6 * d to d , and finally to 1. We set the kernel size as 1 to translate the feature map into the final score matrix. Each 2D convolution layer is followed by a batch normalization layer. We used ReLU as the activation function within the whole score network.

Post-Processing Network. In the PP network, we initialized w as 1, s as log(9), α as 0.01, β as 0.1, γ α as 0.99, γ β as 0.99, and ρ as 1. We set T as 20.

Training details. During training, we first pre-trained a deep score network and then fine-tuned the score network and the PP network together. To pre-train the score network, we used binary crossentropy loss and Adam optimizer. Since, in the contact map, most entries are 0, we used weighted loss and set the positive sample weight as 300. The batch size was set to fully use the GPU memory, which was 20 for the Titan Xp card. We pre-train the score network for 100 epochs. As for the fine-tuning, we used binary cross-entropy loss for the score network and F1 loss for the PP network and summed up these two losses as the final loss. The user can also choose to only use the F1 loss or use another coefficient to weight the loss estimated on the score network U θ . Due to the limitation of the GPU memory, we set the batch size as 8. However, we updated the model's parameters every 30 steps to stabilize the training process. We fine-tuned the whole model for 20 epochs. Also, since the data for different RNA families are imbalanced, we up-sampled the data in the small RNA families based on their size. For the training of the score network U θ in the ablation study, it is exactly the same as the training of the above mentioned process. Except that during the fine-tune process, there is the unrolled number of iterations is set to be 0. (iii) Samples in ArchiveII dataset are i.i.d. from the distribution P(ArcII).

To compare the differences among these data distributions, we can test the following hypothesis: 

:

The codes for reproducing the experimental results are released at https://github.com/ml4bio/e2efold.

CONCLUSIONWe propose a novel DL model, E2Efold, for RNA secondary structure prediction, which incorporates hard constraints in its architecture design. Comprehensive experiments are conducted to show the superior performance of E2Efold, no matter on quantitative criteria, running time, or visualization. Further studies need to be conducted to deal with the RNA types with less samples. Finally, we believe the idea of unrolling constrained programming and pushing gradient through post-processing can be generic and useful for other constrained structured prediction problems.

Figures and Tables:

Figure 1: Graph and matrix representations of RNA secondary structure.

Figure 2: Nested and non-nested structures.

Figure 3: Output space of E2Efold.

Fig 4. It mainly consists of • a position embedding matrix P which distinguishes {x i } L i=1 by their exact and relative positions:

Figure 4: Architecture of Deep Score Network.

Figure 5: Distribution of F1 score.

Figure 6: The RNAStralign length distribution.

(a) P(RNAStr train ) = P(RNAStr test ) (b) P(RNAStr train ) = P(ArchiveII)

Dataset Statistics

Table content: <table><tr><td>Type</td><td colspan="2">ArchiveII</td><td colspan="2">RNAStralign</td></tr><tr><td/><td colspan="4">length #samples length #samples</td></tr><tr><td>All</td><td colspan="4">28∼2968 3975 30∼1851 30451</td></tr><tr><td colspan="3">16SrRNA 73∼1995 110</td><td colspan="2">54∼1851 11620</td></tr><tr><td colspan="5">5SrRNA 102∼135 1283 104∼132 9385</td></tr><tr><td>tRNA</td><td>54∼93</td><td>557</td><td colspan="2">59∼95 6443</td></tr><tr><td>grp1</td><td>210∼736</td><td>98</td><td colspan="2">163∼615 1502</td></tr><tr><td>SRP</td><td>28∼533</td><td>928</td><td colspan="2">30∼553 468</td></tr><tr><td colspan="3">tmRNA 102∼437 462</td><td colspan="2">102∼437 572</td></tr><tr><td colspan="3">RNaseP 120∼486 454</td><td colspan="2">189∼486 434</td></tr><tr><td colspan="2">telomerase 382∼559</td><td>37</td><td colspan="2">382∼559 37</td></tr><tr><td colspan="3">23SrRNA 242∼2968 35</td><td>-</td><td>-</td></tr><tr><td>grp2</td><td>619∼780</td><td>11</td><td>-</td><td>-</td></tr></table>



Table content: <table/>

Performance comparison on ArchiveII

Table content: <table><tr><td>Method</td><td>Prec Rec</td><td>F1 Prec(S) Rec(S) F1(S)</td></tr><tr><td>E2Efold</td><td colspan="2">0.734 0.66 0.686 0.758 0.676 0.704</td></tr><tr><td>CDPfold</td><td colspan="2">0.557 0.535 0.545 0.612 0.585 0.597</td></tr><tr><td colspan="3">LinearFold 0.641 0.617 0.621 0.668 0.644 0.647</td></tr><tr><td>Mfold</td><td colspan="2">0.428 0.383 0.401 0.450 0.403 0.421</td></tr><tr><td colspan="3">RNAstructure 0.563 0.615 0.585 0.590 0.645 0.613</td></tr><tr><td>RNAfold</td><td colspan="2">0.565 0.627 0.592 0.586 0.652 0.615</td></tr><tr><td colspan="3">CONTRAfold 0.607 0.679 0.638 0.629 0.705 0.662</td></tr></table>

Inference time on RNAStralignTest On ArchiveII Without Re-training. To mimic the real world scenario where the users want to predict newly discovered RNA's structures which may have a distribution different from the training dataset, we directly test the model learned from RNAStralign training set on the ArchiveII dataset, without re-training the model. To make the comparison fair, we exclude sequences that are overlapped with the RNAStralign dataset. We then test the model on sequences in ArchiveII that have overlapping RNA types (5SrRNA, 16SrRNA, etc) with the RNAStralign dataset. Results are shown in Table

Table content: <table><tr><td>Method</td><td/><td colspan="2">total run time time per seq</td></tr><tr><td colspan="3">E2Efold (Pytorch) 19m (GPU)</td><td>0.40s</td></tr><tr><td colspan="4">CDPfold (Pytorch) 440m*32 threads 300.107s</td></tr><tr><td>LinearFold</td><td>(C)</td><td>20m</td><td>0.43s</td></tr><tr><td>Mfold</td><td>(C)</td><td>360m</td><td>7.65s</td></tr><tr><td colspan="2">RNAstructure (C)</td><td>3 days</td><td>142.02s</td></tr><tr><td>RNAfold</td><td>(C)</td><td>26m</td><td>0.55s</td></tr><tr><td colspan="2">CONTRAfold (C)</td><td>1 day</td><td>30.58s</td></tr></table>

Evaluation of pseudoknot prediction

Table content: <table><tr><td>Method Set F1 TP FP TN FN</td></tr><tr><td>E2Efold 0.710 1312 242 1271 0</td></tr><tr><td>RNAstructure 0.472 1248 307 983 286</td></tr></table>

RNAstructure and CONTRAfold produce comparatively more reasonable visualizations among all baselines, so we compare with them. These two methods can capture a rough sketch of the structure, but not good enough. For most cases, E2Efold produces structures most similar to the ground-truths. Moreover, it works surprisingly well for some RNA sequences that are long and very difficult to predict. Ablation study (RNAStralign test set)

Table content: <table><tr><td>true structure</td><td>E2Efold</td><td>RNAstructure</td><td>CONTRAfold</td><td>true structure</td><td>E2Efold</td><td>RNAstructure</td><td>CONTRAfold</td></tr><tr><td colspan="5">Visualization. We visualize predicted structures of three</td><td/><td/><td/></tr><tr><td colspan="5">RNA sequences in the main text. More examples are</td><td/><td/><td/></tr><tr><td colspan="5">provided in appendix (Fig 8 to 14). In these figures,</td><td/><td/><td/></tr><tr><td colspan="5">purple lines indicate edges of pseudoknotted elements.</td><td/><td/><td/></tr><tr><td colspan="5">Although CDPfold has higher F1 score than other base-</td><td/><td/><td/></tr><tr><td colspan="5">true structure lines, its predictions are visually far from the ground-truth. Instead,</td><td>true structure E2Efold</td><td/><td>E2Efold</td></tr></table>

RNAStralign dataset splits statistics Samples in RNAStralign training set are i.i.d. from the distribution P(RNAStr train ); (ii) Samples in RNAStralign testing set are i.i.d. from the distribution P(RNAStr test );

Table content: <table><tr><td>RNA type</td><td>All</td><td colspan="3">Training Validation Testing</td></tr><tr><td>16SrRNA</td><td>11620</td><td>9325</td><td>1145</td><td>1150</td></tr><tr><td>5SrRNA</td><td>9385</td><td>7687</td><td>819</td><td>879</td></tr><tr><td>tRNA</td><td>6443</td><td>5412</td><td>527</td><td>504</td></tr><tr><td>grp1</td><td>1502</td><td>1243</td><td>123</td><td>136</td></tr><tr><td>SRP</td><td>468</td><td>379</td><td>36</td><td>53</td></tr><tr><td>tmRNA</td><td>572</td><td>461</td><td>50</td><td>61</td></tr><tr><td>RNaseP</td><td>434</td><td>360</td><td>37</td><td>37</td></tr><tr><td>telomerase</td><td>37</td><td>28</td><td>4</td><td>5</td></tr><tr><td colspan="2">RNAStralign 30451</td><td>24895</td><td>2702</td><td>2854</td></tr><tr><td colspan="2">D.2 TWO-SAMPLE HYPOTHESIS TESTING</td><td/><td/><td/></tr><tr><td colspan="5">To better understand the data distribution in different datasets, we provide statistical hypothesis test</td></tr><tr><td>results in this section.</td><td/><td/><td/><td/></tr><tr><td>We can assume that</td><td/><td/><td/><td/></tr><tr><td>(i)</td><td/><td/><td/><td/></tr></table>

Figure 10: Visualization of Group I intron, IC3, Kaf.c.trnL.Figure 13: Visualization of tmRNA, uncu.bact. AF389956.

Table content: <table><tr><td>true structure</td><td>E2Efold</td><td>RNAstructure</td><td>CDPfold</td></tr><tr><td>true structure</td><td>E2Efold</td><td>RNAstructure</td><td>CDPfold</td></tr><tr><td>CONTRAfold</td><td>LinearFold</td><td>Mfold</td><td>RNAfold</td></tr><tr><td>CONTRAfold</td><td>LinearFold</td><td>Mfold</td><td>RNAfold</td></tr><tr><td>true structure</td><td>E2Efold</td><td>RNAstructure</td><td>CDPfold</td></tr><tr><td>CONTRAfold</td><td>LinearFold</td><td>Mfold</td><td>RNAfold</td></tr><tr><td>true structure</td><td>E2Efold</td><td>RNAstructure</td><td>CDPfold</td></tr><tr><td/><td colspan="2">Figure 11: Visualization of RNaseP, A.salinestris-184.</td><td/></tr><tr><td>true structure</td><td>E2Efold</td><td>RNAstructure</td><td>CDPfold</td></tr><tr><td>CONTRAfold</td><td>LinearFold</td><td>Mfold</td><td>RNAfold</td></tr><tr><td>CONTRAfold</td><td>LinearFold</td><td>Mfold</td><td>RNAfold</td></tr></table>

Peer review (iclr2020):
This paper introduces an end-to-end method to predict the secondary structure of RNA, by mapping the nucleotide sequence to a binary affinity matrix. The authors decompose this problem into two part: (i) predicting an affinity score between each base pair in the input sequence, using a combination of a transformer sequence encoder network and a convolutional decoder, and (ii) a post-processing step that ensures that structural local and global constraints are enforced. An innovation is to express this post-processing as an unrolled sequence of proximal gradient descent steps, which are fully differentiable, and allow the full combination of (i)+(ii) to be trained end-to-end. A thorough set of experiments validate the approach. Overall, the paper is well written and easy to follow. The approach of unrolling structural constraints as shown in the paper is interesting and applicable to much wider domains than secondary structure prediction. The proposed approach appears to provide a novel, convincing and non-obvious solution to RNA secondary structure prediction, and subject to suggestions below, would represent a valuable contribution to ICLR. The principal area for improvement would be to include additional detail (perhaps in appendix) on the model hyperparameter configurations that were used in the experiments. Moreover, more details on the set of psi functions, and the MLP details for P_i (e.g. number of hidden units, activation function, the use of dropout, batch normalization, etc) should be given, as well as more information on the specifics how how the “pairwise concatenation” is carried out in the output layer. What unrolling constant T is used? Finally, in the ablation study (p. 8) details on how U_	heta is trained by itself (without the post-processing step) should be given. Detailed comments: * Overall, the whole paper should be thoroughly reviewed for English grammar and writing style; a subset of suggested changes follow. * p. 1: structure a result ==> structure is a result * p. 2: energy based methods ==> energy-based methods * p. 2: energy function based approaches ==> energy function-based approaches * p. 2: view point ==> viewpoint * p. 2: E2Efold is flexible ==> E2Efold are flexible * p. 2: nearly efficient ==> nearly efficiently * p. 3: typically scale ==> typically scale as * p. 3: few hundred. ==> few hundreds. * p. 4: all binary matrix ==> all binary matrices * p. 4: output space can help ==> output space could help * p. 5: formulation are the ==> formulation are that the * p. 6: eq. (7) should contain quantities indexed by in the RHS * p. 8: pesudoknotted ==> pseudoknotted * p. 9 ff: in the bibliography, all lowercase rna should be uppercase RNA. Use {RNA} in bibtex entries.