致病氨基酸变异预测的新型融合模型

邵爱斌，杨 洋*

(1. 苏州大学计算机科学与技术学院 江苏 苏州 215006；2. 软件新技术与产业化协同创新中心 南京 210000)

随着测序成本的大幅降低，测序方法在科研和临床医疗中被广泛使用，由此产生了大量测序信息，也包括了越来越多的变异信息。基因变异导致的氨基酸变异可通过多种方式影响蛋白质的结构和功能。当变异发生在蛋白质的某些关键部位，如催化部位或配体相互作用表面，可能导致蛋白质折叠、结构不稳定或蛋白质聚集，进而导致疾病。为实现个性化医疗，追溯疾病发生的机理，预测氨基酸变异的致病性具有很高的研究价值。

与实验方法相比，计算方法具有预测成本低、效率高的突出优势。近年来，研究者们提出并不断改进了多种相关预测模型。其中，PolyPhen-2[1]面对不同预测任务灵敏性的需求，构建了HumDiv和HumVar 两个数据库，运用朴素贝叶斯预测变异致病性。CADD[2]使用支持向量机算法，整合了63 种基因注释，从而得到C 分数来预测致病性。PON-P2[3]利用基因本体(gene ontology, GO)等特征训练，采用自抽样的方式计算置信度。DEOGEN2[4]用可视化的方式提供了每个预测的结果、相关蛋白质背景和起源。MetaSVM[5]开发了基于支持向量机的元分析框架，框架中SVM 的目标函数由铰链损失和稀疏组套索组成。ClinPred[6]首次使用了ClinVar[7]数据库，并训练了两个模型，分别基于随机森林和XGBoost 来获取最高预测结果。PrimateAI[8]结合6 个非人类灵长类动物物种和人类的变异，共收集到38 万条变异数据，训练了一个包含36 层卷积神经网络的深度学习模型。

随着当前可采集的变异数据量的增加，构建一个新的融合模型以提高预测性能变得可行。因此，本文尝试使用深度学习方法从蛋白序列中提取出一些特征，将这些特征与提取并筛选的有效生物特征融合，作为模型的输入，并构建融合模型训练，旨在达到较高的预测性能。

1 数据来源

本文所使用的数据包括人类、动物和植物蛋白序列中的氨基酸变异样本，如表1 所示。其中，人类致病变异取自VariBench[9-10]数据库和ClinVar[7]数据库，共有17631 个。按致病变异与无害变异大约1∶1 的比例，从ExAC[11]数据库获取共18494 个人类无害变异。对于动物变异数据：1)收集OMIA[12]数据库中有“likely causal variants”标记的变异；2)从文献[13]获得其他哺乳动物(狗、鼠和牛)的变异。动物致病氨基酸变异共317 个，中性变异共312 个。进一步从文献[14]取得植物变异数据集，其物种包含拟南芥、水稻和豌豆，数据集由3236 个有害变异和1899 个中性变异构成。

表1 实验数据来源和分布 个

由于某些蛋白序列过长，影响预测性能，本文筛选出长度不超过1500 个的蛋白质序列，共9980条，包含了35179 个氨基酸变异，其中致病变异18521 个，无害变异16658 个。

2 方 法

本文所构建的新型融合模型工作流程如图1所示。

图1 致病氨基酸变异预测的新型融合模型的工作流程

2.1 生物特征提取

综合文献[3, 15]中的特征，本文共提取了1085个可能与氨基酸变异致病性有关的生物特征，包括617 个氨基酸特征、436 个变异类型特征、25 个邻域特征、2 个进化保守特征、3 个蛋白质特征以及GO 频率对数比特征(logical rate, LR)和功能位点特征(functional site, FS)。

617 个氨基酸特征通过AAindex 表[16]提取，对于每个变异，计算氨基酸变异前后得分的差异。

436 个变异类型特征通过两个矩阵获得。其中，400 个特征来自20×20 的矩阵，矩阵的两个维度分别代表变异前和变异后的氨基酸。另外36 个特征来自代表氨基酸物理和化学性质的6×6 矩阵[17]。

25 个邻域特征表示变异位点的序列位置特征。20 维向量在23 个位置，即变异位点前后11个位置的窗口中统计每种氨基酸类型的频数[18]。此外，在23 个位置的邻域窗口中统计了5 种氨基酸类别，包括非极性、极性、带电、带正电和带负电的频数。

使用DIAMOND[19]将每个蛋白质序列在SwissProt[20]中进行同源序列比对，找到同源序列并统计命中数。比对后，使用SIFT 4G[21]计算每个变体的SIFT 分数。

3 个蛋白质特征包括蛋白质序列的长度，变异是否发生在序列的第一位氨基酸，以及变异在序列中的位置。

利用朴素贝叶斯算法计算特征LR 和FS 如下：

1)使用R Bioconductor 工具GO.db(https://bioconductor.org/packages/GO.db/)，从AmiGO[22]与QuickGO[23]中收集蛋白质的所有GO 术语以及这些术语对应的全部上层GO 术语，并过滤以避免重复。然后，在数据集中对收集到的GO 术语分别统计两个类别，即致病和无害的频数，计算该蛋白GO 的LR 如下：

式中，f(Pi)是第i个GO 术语在致病集的频数；f(Ni)是第i个GO 术语在无害集中的频数。为了避免分数无意义，频数都加1。如果一个蛋白质没有GO 术语，则LR 为0。

2)根据UniProtKB/SwissProt[20]查找变异位点的Site 术语，以类似LR 的计算方法得到该蛋白Site 的功能位点特征。

2.2 生物特征选择

特征选择也称属性选择，通过从当前全部特征中筛选出部分特征的方式，达到模型某项指标最优化的目的。较多的特征中可能存在无关特征和冗余特征，这也许会导致维度灾难，使得模型的性能下降，训练的时间开销增加。当特征较多而样本数量较少时，模型往往容易出现过拟合的问题。为了避免特征冗余和模型过拟合的风险，减少算力资源浪费，增加模型可解释性，需要筛选出与分类性能密切相关的生物特征。

参照文献[1, 3]的特征选择方法，本文使用LightGBM 算法结合递归特征消除(recursive feature elimination, RFE)方法[24]进行特征选择。特征选择的过程是以全部特征构建模型，为每个特征分配权重，并去除权重最小的特征，然后在剩余的特征上重复上述操作，直到达到预设的特征数量。

2.3 通过神经网络提取特征

随着数据量增加，有望通过深度学习方法提取出一些没有明确生物意义但对分类有帮助的特征。本文构建的深度学习网络主要包括两个部分：卷积神经网络(convolutional neural networks, CNN)和双向长短期记忆神经网络(bi-directional long-short term memory, Bi-LSTM)。

CNN 提取特征的基本思路是通过不断平移卷积核，对原始矩阵的每个元素进行卷积操作，从而提取到某些细节特征。卷积核的感受野与卷积核大小有关，感受野不同，卷积网络能提取到的特征细节也不一样。本模型设置了多条卷积线路，不同线路对应不同卷积核。通过融合多种卷积核的卷积结果，可以提取到较全面的特征。本模块使用了多层的小卷积核，代替少层的大卷积核。在相同大小感受野的前提下，能够缩小计算量，提高训练效率。对于每层卷积提取的特征，再利用激活函数进行非线性变换，以提取更具有泛化性的特征。

特征提取中，本文进一步使用了Bi-LSTM，其构造中有记忆门和遗忘门，经过训练，能够选择记忆和遗忘某些信息，探索到CNN 输出的特征向量中距离较远的关系。Bi-LSTM 由正向LSTM 和逆向LSTM 组成，氨基酸变异位点特征同时与该位点前后序列都有关。经过特征累计，Bi-LSTM能够提取出位点前后向的双向特征。

2.4 融合预测算法

提取并筛选出的生物特征与深度学习提取到的特征以拼接的形式相融合，作为融合预测算法的输入。融合特征与预测结果之间的关系比较复杂，多个模型的融合或许能取得较好的性能。不同的模型从不同的角度去观测数据。堆叠法(stacking)可以整合多个模型观测数据的角度，相互取长补短，能够得到一个更加全面和优秀的结果。LightGBM[25]采用了基于梯度的单边采样，留下梯度大的样本，随机抽样梯度小的样本，来分割数据样本。XGBoost[26]通过直方图算法，将全部样本按照某个特征分成离散区域，并根据这些区域确定最佳分割。CatBoost 采用贪婪策略，整合不同类别型特征得到新的特征，能够直接处理类别型特征。所以，XGBoost、CatBoost 和LightGBM 这3 种提升boosting 算法各有特点。

本融合模型选择这3 个算法作为第一层模型，训练并保留预测结果，得到3 列新的特征表达。第二层随机森林用得到的特征表达进行训练，得到最终的分类结果。以此达到平衡各算法的优缺点，提高分类性能的目的。

2.5 模型评价

本文引入了一些评价指标对模型进行评估，来直观辨别分类性能的优劣，其中包括阳性预测值(PPV)、阴性预测值(NPV)、敏感性(TPR)、特异性(TNR)、准确性(ACC)、马修斯相关系数(MCC)和总体绩效指标(overall performance measure, OPM)。这些评价指标的数学定义如下：

MCC 在TP、TN、FP 和FN 基础上评价二元分类性能，其值在−1～+1 之间，值越大分类性能越好。OPM 由PON-P2[3]定义，将MCC 归一化为0～1 值(nMCC)，结合PPV、NPV、TPR、TNR和ACC，如式(9)，两两之和为棱长作长方体，长方体体积的1/8 作为OPM 的值。OPM 的取值范围为0～1，值越接近1，分类性能越好。

3 实验结果

在数据集中随机选取90%的变异样本作为训练集进行10 重交叉验证。实验用的部分数据和代码已经上传GitHub (https://github.com/s1mpleCN/PONDML)。

3.1 生物特征筛选结果

参照文献[1, 3]的特征筛选经验，本文用RFE 算法分别筛选了100、50、20、10 个特征进行比较。其交叉验证结果如表2 所示，其中仅基于20 个生物特征，预测准确性ACC 即可达到89.7%，MCC 为79.4%，OPM 为72.2%。综合考虑特征数量和模型性能的平衡，本文最终选择了筛选出的20 个生物特征来进行特征融合。

表2 筛选不同数量特征的CV 结果

3.2 通过神经网络提取特征的结果

考虑到原始蛋白序列是由英文字母组成的，无法直接用深度学习方法训练。借鉴自然语言处理的做法，将原始蛋白质序列按照氨基酸单字母缩写的顺序从小到大进行编码[27]。神经网络中的Bi-LSTM要求输入的序列长度限制为固定长度，在长度短的序列后面填充0，使序列最终长度均为1500 个氨基酸。通过词嵌入方法将序列中每个数字编码，转化为一个实数向量。本文设置嵌入参数为128，输出为1500×128 的矩阵。

本文构建的CNN 使用了一维卷积函数，只需要定义核大小的一个参数，另一参数默认取词嵌入部分设置的长度128。CNN 共构造了3 条卷积路线，分别对应3 种大小的卷积核(7×128、5×128和3×128)。每条路线中都包含5 个卷积层，每个卷积层之后都进行最大池化操作，从而提取出更显著的特征。将3 条路线提取的特征以拼接的方式融合，作为Bi-LSTM 的输入。经过实验比较，Bi-LSTM 层的输出参数设为70 时，预测性能最好。此时，正向LSTM 和逆向LSTM 都可以提取到70 个特征。

在构建的神经网络中加入丢弃层(dropout)，每次随机选择训练过程中的部分神经元进行计算，避免模型过拟合。网络中加入批标准化层(batch normalization)，可以避免梯度爆炸，加快迭代速度。综合考虑模型的准确性和泛化性，利用早停机制(early stopping)，共进行30 次迭代。深度学习网络参数详见表3。

表3 深度学习网络中的参数

神经网络训练的模型性能如表4 所示，其准确性ACC 为84.9%，MCC 为69.7%，OPM 达到61.2%。

表4 神经网络训练的模型性能结果

3.3 融合模型交叉验证性能结果

如图2 所示，本融合模型将每次交叉验证的训练集随机分成3 份，分两层进行训练。首先选择XGBoost，每一次用其中两份数据预测剩下的1 份数据以及验证集的数据，保留得到的预测结果。每份数据都预测一次，共可以保留3 份训练集的预测结果和3 份验证集的预测结果。整合训练集的预测结果作为第二层的一列新的特征表达。将3 份验证集的预测结果取平均值用作第二层验证集的一列特征表达。然后用CatBoost 和LightGBM 重复上面步骤。最后，在第二层中用随机森林进行训练预测，相当于对第一层的模型做了融合。

图2 融合预测算法工作流程

本文融合模型交叉验证结果如表5 所示，其中准确性ACC 为92.8%，MCC 达到85.5%，OPM 值为79.8%。

表5 融合模型10 重交叉验证结果

3.4 盲测结果及与其他工具比较

以剩余的10%数据作为盲测集，使用构建好的融合模型进行预测。其结果如表6 第二列所示，本文融合模型准确性ACC 达到93.1%，MCC 为86.1%，OPM 为80.6%。

本文融合模型在盲测集上的ROC 曲线如图3所示，曲线下面积(AUC)为0.921。

图3 盲测集ROC 曲线

本文还在该盲测集上比较了几种常见的预测工具性能，包括Polyphen2[1]、CADD[2]、PON-P2[3]、DEOGEN2[4]、MetaSVM[5]、ClinPred[6]和PrimateAI[8]。其中CADD 阈值设为20 时，其性能取得最佳。

由表6 可见，Polyphen2 准确性较低，为75.0%。考虑到Polyphen2 训练的数据量仅为2 万条，限制了它的性能。CADD 准确性达到77.5%。该模型以36 个其他预测模型的得分作为特征，或许某些预测模型并不能有效提高氨基酸变异致病性的预测性能。定义覆盖率为模型能够预测的变异数/测试集总样本数量。PON-P2 准确性最高，达到93.7%，但覆盖率仅为54.2%，这是因为PON-P2 舍弃了置信度不高的变异预测。而本融合模型覆盖率达到了100.0%。由于综合了两个预测模型，ClinPred 准确性较高，为91.1%，但覆盖率也仅为76.5%。PrimateAI 是一个基于深度学习的预测模型，直接从序列中提取特征，准确性较低，仅为56.0%。可见生物特征在变异致病性预测中扮演着更为重要的角色。

表6 与常用的预测模型在盲测集上的性能结果比较

4 结束语

氨基酸变异常常会对蛋白质的结构和功能造成影响，进而导致疾病。基于计算的方法作为一种预测氨基酸变异致病性的有效途径，被研究者广泛使用。为了提高预测准确性和泛化性，本文构建了一个新型融合模型。首先，收集到包含人类、动物和植物种群的氨基酸变异作为数据集。接着，提取有效生物特征，并用RFE 筛选出最优特征子集。然后，使用深度学习网络提取特征，深度学习网络由CNN 和Bi-LSTM 组成。将筛选完的生物特征和深度学习网络提取的特征以拼接的方式融合，作为预测输入。最后，构建一个基于XGBoost、CatBoost、LightGBM 和随机森林的融合模型，得到最终的预测结果。融合模型的交叉验证准确性ACC 达到92.8%，MCC 达到85.5%，OPM 值为79.8%。与其他工具相比，本文模型具有更高的准确性和泛化性。该模型工具可用于辅助临床诊断和药物设计，降低研发成本。