RNA结构预测竞赛圆满落幕，且看前三名团队飞桨论道(下篇)

点击左上方蓝字关注我们

由百度飞桨举办的螺旋桨RNA结构预测竞赛圆满结束。本次竞赛依托百度研究院在RNA结构预测上的算法优势和百度大脑AI Studio平台优势，汇聚生物计算、人工智能领域从业者和爱好者，探索人工智能技术在生物计算上的技术应用。本次比赛邀请北京大学生命科学学院刘君教授团队提供实验数据，构建测试集，保证比赛的科学性和公正性。各获胜队伍已在比赛页揭晓。今天有幸邀请到冠军队伍分享各自的技术方案，供大家学习参考。快来解锁冠军大佬的冲榜秘籍吧！

冠军-基因疗法团队

成员：邱祥运。邱教授20多年来致力于研究以X射线散射为主要测量手段的科研实验工作, 近年来从事生物大分子力学和结构学基础研究, 近来专注于核酸医学, 借力大数据/深度学习, 探索核酸生物功能、疾病机理和药物开发。

昨天的上篇部分，我们为大家介绍了RNA背景知识，还展示了冠军团队对赛题、数据以及评测指标的理解与分析，内容详见：

RNA结构预测竞赛圆满落幕，且看冠军团队飞桨论道(上篇)

今天我们继续解锁冠军团队的模型搭建及优化过程！

冠军团队方案详解

1. 模型搭建

本方案基于飞桨2.0深度学习平台，过程中参考了大量详细的飞桨文档和开源代码，尤其利用了Paddle.nn中直接可用的TransformerEncoder、LSTM、Conv1D等，损失函数主要采用Paddle.nn提供的函数库。为了更便捷地搭建不同模型并更改模型参数，可以把主要参数存储在一个args结构里，然后包装每个模块(比如LSTM, Conv1D)能接受args为输入。

网络框架构建代码范例：

class Seq2Seq_AttnLSTMConv1DNet(nn.Layer):
    def __init__(self, args):
        super(Seq2Seq_AttnLSTMConv1DNet, self).__init__()
        nn.initializer.set_global_initializer(nn.initializer.KaimingNormal(),
                                              nn.initializer.Constant(0.0))

        self.data_format = args.input_fmt.upper()
        self.feature_dim = int(args.feature_dim)

        self.embed = MyEmbeddingLayer(args, in_features=self.feature_dim)
        self.fc_in = MyLinearTower(args, in_features=self.embed.out_features)
        self.attn = MyAttnTower(args, in_features=self.fc_in.out_features)
        self.lstm = MyLSTMTower(args, in_features=self.attn.out_features)
        self.conv1d = MyConv1DTower(args, in_features=self.lstm.out_features)

模块构建代码范例：

class MyLinearTower(nn.Layer):
    def __init__(self, args, in_features=None, is_return=False):
        """ is_return:True will trun off Act/Norm/Dropout for the last block """
        super(MyLinearTower, self).__init__()
        nn.initializer.set_global_initializer(nn.initializer.KaimingNormal(), nn.initializer.Constant(0.0))

        self.is_return = is_return
        if in_features is None:
            self.in_features = int(args.feature_dim)
        else:
            self.in_features = int(in_features)

        self.data_format = 'NLC'
        self.act_fn = args.act_fn
        self.norm_fn = args.norm_fn
        self.norm_axis = int(args.norm_axis)
        self.dropout = float(args.dropout)

        self.linear_dim = [int(_i) for _i in args.linear_dim] \
                if hasattr(args.linear_dim, '__len__') else [int(args.linear_dim)]
        self.linear_num = int(args.linear_num)
        self.linear_resnet = args.linear_resnet

        in_features = self.in_features
        self.linear_layers = [] # addditional layers if needed
        for i in range(self.linear_num):
            is_return = (i == self.linear_num -1) if self.is_return else False
            self.linear_layers.append(nn.Sequential(*MyLinearBlock(
                    [in_features] + self.linear_dim,
                    dropout = self.dropout,
                    act_fn = self.act_fn,
                    norm_fn = self.norm_fn,
                    norm_axis = self.norm_axis,
                    data_format = self.data_format,
                    is_return = is_return,
            )))

对于一个长度为L的RNA序列，每个碱基字符用一个四维的onehot向量表达，LinearFold的二级结构字符用一个四维的onehot向量表达(三维足够，这里凑了个偶数)，LinearPartition预测可以直接使用。

最终得到的输入矩阵的维度为[N, L, 10]，其中N是batch size，L是序列长度。输出全连接模块的最后一层维度分别为2，所以模型的输出是一个[N, L, 2]维的矩阵。沿最后一维进行softmax操作，得到每个碱基的配对和不配对概率。

损失函数有两种可取的方法，最直接的是对不配对概率计算竞赛采用的均平方差的平方根(RMSD)，同时考虑结果为softmax后的概率，标注为二分类的软标注，我们也可以采用交叉熵做为损失。

两种方法的不同主要在于交叉熵在概率p的梯度上多了一个1/p(1-p)的因子，对接近0或1的概率梯度增强。在实验中我们发现两种损失函数的结果大体相同，猜想是因为概率主要分布在0和1附近，导致这个附加因子的作用不显著。

2. 模型训练

模型训练每个epoch大约5分钟，一般1-2个小时能训练结束，相关代码在fly_paddle.py脚本里。前文提到的args包含了大部分模型和训练的参数, fly_paddle.py设置了缺省参数，同时可以在命令行设置参数取值。

命令行启动模型训练范例：

{
   "data_dir": "data",
   "data_name": null,
   "data_suffix": ".pkl",
   "test_size": 0.1,
   "split_seed": null,
   "batch_size": 4,
   "input_genre": "seq",
   "input_fmt": "NLC",
   "residue_dbn": false,
   "residue_attr": false,
   "residue_extra": false,
   "label_genre": "upp",
   "label_fmt": "NL",
   "label_tone": "soft",
   "label_ntype": 2,
   "label_smooth": false,
   "loss_fn": ["softmax+mse"],
   "loss_sqrt": true,
   "loss_padding": false
}

模型训练中的进度汇报范例：

2021-06-28 14:42:43,270 - INFO - Training, data size: 4500
2021-06-28 14:42:43,272 - INFO -          batch size: 1
2021-06-28 14:42:43,273 - INFO -             shuffle: True
2021-06-28 14:42:43,275 - INFO -        # of batches: 4500
2021-06-28 14:42:43,370 - INFO -      recap interval: 151
2021-06-28 14:42:43,371 - INFO -   validate interval: 450
2021-06-28 14:42:43,371 - INFO -         # of epochs: 21
2021-06-28 14:42:43,372 - INFO -        loss padding: False
2021-06-28 14:42:44,322 - INFO - Epoch/batch: 0/   0, ibatch:    0, loss: 0.8399, std: 0.5991
2021-06-28 14:42:52,789 - INFO - loss: 0.8325, std: 0.5606
2021-06-28 14:42:59,262 - INFO - Epoch/batch: 0/ 151, ibatch:  151, loss: 0.6945, std: 0.6275
2021-06-28 14:43:05,585 - INFO - Epoch/batch: 0/ 302, ibatch:  302, loss: 0.6000, std: 0.6642
2021-06-28 14:43:22,643 - INFO - loss: 0.5584, std: 0.6602
2021-06-28 14:43:22,676 - INFO - Saved model states in: earlystop_0.5584
2021-06-28 14:43:22,678 - INFO - Saved net python code: earlystop_0.5584/paddle_nets.py
2021-06-28 14:43:22,691 - INFO - Saved best model: earlystop_0.5584
2021-06-28 14:43:22,834 - INFO - Epoch/batch: 0/ 453, ibatch:  453, loss: 0.5834, std: 0.6595
2021-06-28 14:43:29,579 - INFO - Epoch/batch: 0/ 604, ibatch:  604, loss: 0.5931, std: 0.6547
2021-06-28 14:43:35,884 - INFO - Epoch/batch: 0/ 755, ibatch:  755, loss: 0.5775, std: 0.6559
2021-06-28 14:43:52,071 - INFO - loss: 0.5481, std: 0.6405
2021-06-28 14:43:52,090 - INFO - Saved model states in: earlystop_0.5481
2021-06-28 14:43:52,092 - INFO - Saved net python code: earlystop_0.5481/paddle_nets.py
2021-06-28 14:43:52,101 - INFO - Saved best model: earlystop_0.5481
2021-06-28 14:43:52,101 - INFO - Removing earlystop model: earlystop_0.5584

模型搭建和优化共分为三个步骤:

第一步，培养对深度学习任务的直觉。先用最简单的全连接模型看一下不考虑碱基间相互作用时的RMSD。如果输入的碱基信息只有其字符编码的四维onehot向量，RMSD为～0.42，和随机猜测概率相差无几。加入LinearFold和LinearPartition后，预测信息能降低RMSD到0.25左右。可见LinearFold和LinearPartition给出了很好的初始点，如果要进一步降低RMSD, 就需要考虑碱基间的相互作用。

第二步，分步加入TransformerEncoder、LSTM和Conv1D。发现单一的TransformerEncoder或者LSTM足够降低RMSD到[0.21, 0.22]区间，Conv1D模块对RMSD影响较小, 与期望大致相符，因为一层的Conv1D只是考虑到了紧邻碱基间的作用。最后三个模块一起能降低RMSD到[0.20, 0.21]区间。

第三步，摸索超参数并改进模型来进一步降低RMSD。首先对learning_rate、dropout和weight decay做了一系列摸索，发现在常用的超参数区间里模型表现相当，可能原因有2个：(1)参数空间相对较小；(2)起始的数据和标注相距不远。因为训练集和验证集的RMSD在所有的训练中RMSD相当, 故没有采用参数正则化。

3. 方案总结

图1展示了在训练集上模型预测的碱基不配对概率分布，第一个突出的问题是有39%碱基的不配对概率在[0.1, 0.9]区间，而标注集中只有22%碱基的不配对概率在此区间。一个简单直接的方法是增加[0.1, 0.9]区间对损失的贡献，比如从损失中减去(概率-0.5)的平方, 然而这方面的尝试未能减小RMSD。

另一个明显的问题是RMSD过大，常见的解决方法是增加模型的深度或宽度。更大的模型能更好地拟合训练数据(比如模块维度为128、层数为3时可达到RMSD< 0.1)，验证数据的RMSD停留在0.2左右，然而这是一个典型的过拟合现象，故没有采用更大的模型。

图1 模型预测碱基不配对概率分布

4. 改进方向

最后递交的预测结果采用了模型平均。在每次训练中从训练集中随机选取10%的序列做为验证集并存取验证集最小RMSD的模型参数，并通过多次训练得到同一模型的不同参数集，最后平均三个最好参数集的预测结果。

训练集得到的RMSD~0.20, 测试集得到的RMSD差一些, ~0.24，一个可能的原因是测试集有31%的序列长度超过训练和验证集中的最长序列长度, 增加了预测难度。无论是0.20还是0.24, 和实验结果仍有差距, 于是尝试使用以下方法提高预测精度：

(1) 增强输入信息。比如在LinearFold和LinearPartition基础上加入更多计算方法的预测数据。在RNA结构预测领域有很多可用的软件, 比如基于物理的MFOLD[1]、VIENNA RNAFOLD[2]、RNASTRUCTURE[3]，基于机器学习的CONTRAFOLD[4]、CDPFOLD[5]，和近年来快速发展的深度学习方法SPOT-RNA[6]、E2EFOLD[7]、MXFOLD2[8]、DIRECT[9]、RNACONCAT[10]、DMFOLD[11]等。

(2)训练更大模型。较小的模型难于捕捉RNA碱基配对中受支配的碱基间的多体相互吸引与竞争。尤其是当RNA序列长度更长(比如新冠病毒的基因有~30,000个碱基, mRNA疫苗也有几千个碱基), 并且其碱基经过大量的化学修饰后, 预测其结构更是难上加难。要能够预测更长更复杂的RNA, 更大的模型几乎是必须的。图2展示了模型和标注之间的对比, 上下两组的长度分别为110和480, 左列两图展示了标注(蓝)、预测(红)和它们的差(绿)对序列位置(X), 右列两图展示了预测(Y)对标注(X)。可以看到较短的序列预测准确度更高, 而长序列的预测差距较大。更大的模型需要更多更全面的训练数据, 尤其对于有化学修饰的碱基, 实验测量很有可能是一个瓶颈，需要另辟蹊径, 比如用模拟的方法得到更多的数据。

图2 模型预测碱基不配对概率和标注的对比

GitHub源代码: 

https://github.com/genetic-medicine/PaddleHelix_RNA_UPP

5. 飞桨使用体验

从纯科学角度来看, 一个理想的模型只需要RNA序列信息就能准确预测碱基不配对概率和其它的结构信息。这并不是可望而不可及的, 在现有的方法里, 百度螺旋桨的LinearFold /LinearPartition和前面提到的第一个端到端深度学习方法SPOT-RNA，都只需要序列就可以非常好地预测RNA的二级结构。

从实用角度来看, 即使一个方法需要综合很多方面的信息(比如需要聚合很多现有软件的预测结果)，能更快更准确的预测RNA结构也已足够。飞桨的API风格与tensorflow和pytorch非常接近，所以对有深度学习经验的同学来说用飞桨能很快上手，速度与灵活性也不比其他深度学习框架差而且还可免费领取GPU算力，对入门AI的朋友来说比较友好。

更多PaddleHelix介绍，扫码即可查看⬇️

参考文献

[1] Zuker, M., 2003. Mfold web server for nucleic acid folding and hybridization prediction. Nucleic Acids Research, 31(13), pp. 3406-3415.

[2] Lorenz, R., Bernhart, S.H., Höner zu Siederdissen C, Tafer, H., Flamm, C., Stadler, P.F. and Hofacker, I.L.,2011. ViennaRNA Package 2.0. Algorithms for Molecular Biology, 6(1),pp.26.

[3] Reuter, J.S. and Mathews, D.H.,2010. RNAstructure: software for RNA secondary structure prediction and analysis. BMC Bioinformatics,11(1), pp.129.

[4] Do, C.B., Woods, D.A. and Batzoglou, S.,2006. CONTRAfold: RNA secondary structure prediction without physics-based models. Bioinformatics,22(14), pp. e90-e98.

[5] Zhang, H., Zhang, C., Li, Z., Li, C., Wei, X., Zhang, B. and Liu, Y., 2019. A New Method of RNA Secondary Structure Prediction Based on Convolutional Neural Network and Dynamic Programming. Frontiers in Genetics,10(467).

[6] Singh, J., Hanson, J., Paliwal, K. and Zhou, Y., 2019. RNA secondary structure prediction using an ensemble of two-dimensional deep neural networks and transfer learning. Nature Communications,10(1), pp.5407.

[7] Chen, X., Li, Y., Umarov, R., Gao, X. and Song, L., 2020. RNA Secondary Structure Prediction by Learning Unrolled Algorithms. [arXiv:2002.05810 p.]. Available from: https://ui.adsabs.harvard.edu/abs/2020arXiv200205810C.

[8] Sato, K., Akiyama, M., and Sakakibara, Y., 2021. RNA secondary structure prediction using deep learning with thermodynamic integration. Nature Communications,12(1), pp. 941.

[9] Jian, Y., Wang, X., Qiu, J., Wang, H., Liu, Z., Zhao, Y. and Zeng, C., 2019. DIRECT: RNA contact predictions by integrating structural patterns. BMC Bioinformatics,20(1), pp. 497.

[10] Sun, S., Wang, W., Peng, Z. and Yang, J., 2020. RNA inter-nucleotide 3D closeness prediction by deep residual neural networks. Bioinformatics,37(8), pp.1093-1098.

[11] Wang, L., Liu, Y., Zhong, X., Liu, H., Lu, C., Li, C. and Zhang, H., 2019. DMfold: A Novel Method to Predict RNA Secondary Structure with Pseudoknots Based on Deep Learning and Improved Base Pair Maximization Principle. Frontiers in Genetics,10(143).

飞桨(PaddlePaddle)以百度多年的深度学习技术研究和业务应用为基础，集深度学习核心训练和推理框架、基础模型库、端到端开发套件和丰富的工具组件于一体，是中国首个自主研发、功能丰富、开源开放的产业级深度学习平台。飞桨企业版针对企业级需求增强了相应特性，包含零门槛AI开发平台EasyDL和全功能AI开发平台BML。EasyDL主要面向中小企业，提供零门槛、预置丰富网络和模型、便捷高效的开发平台；BML是为大型企业提供的功能全面、可灵活定制和被深度集成的开发平台。

END