基于集成学习方法识别导致人类疾病的遗传变异
1. 引言人类疾病大多与遗传因素有关，这些疾病依据其致病机理的不同通常划分为单基因病与多基因病两类。前者也称为孟德尔疾病，由单个基因上的单个遗传变异所导致并符合孟德尔遗传方式，致病机理简单但比较少见；后者又称为复杂疾病，由多个基因上的多个遗传变异、它们之间的相互作用、以及它们与环境的相互作用所导致，致病机理复杂且广泛存在。随着医学遗传学的发展，人们认识到虽然环境和行为等因素对复杂疾病有较大影响，但对肝癌、糖尿病、心脏病等几乎每一种常见的复杂疾病都存在一些关键的致病遗传变异，它们的发生将使个体患病的风险大大增加。因而，识别与复杂疾病相关的遗传变异成为理解这些疾病的遗传学基础，进而预防、诊断和治疗这些疾病的首要步骤[1-3]。
单核苷酸多态性(Single Nucleotide Polymorphisms，SNPs)是指在DNA 水平上由于单个核苷酸位置上的转换、颠换、插入、缺失等变异事件引起的序列多态性，是人类基因组最丰富的一类遗传变异[4]。根据对生物遗传性状的影响，发生在蛋白质编码区的单核苷酸多态性又可分为同义单核苷酸多态性(Synonymous SNPs)和非同义单核苷酸多态性(Nonsynonymous SNPs, nsSNPs)两类。nsSNPs可使其编码的蛋白质序列发生改变，在蛋白质层次上造成氨基酸置换(Amino acid substitution)，从而潜在地影响蛋白质的结构和功能，并可能导致生物性状的改变[5]。在人类遗传疾病的研究中，通常将nsSNPs分为致病和非致病两类。通过对复杂疾病相关联的nsSNPs进行分析，可以帮助人们理解这些疾病的遗传学机理，从而有效地指导这些疾病的治疗和预防。另一方面，对疾病相关nsSNPs的分析也可用于指导个性化临床用药[5]。然而，传统的通过体外定点突变实验和基因敲除／插入实验从庞大的nsSNPs数据库中挖掘出致病nsSNPs往往需要耗费巨大的人力物力[5]。因此，利用机器学习的方法准确快速的识别致病nsSNPs越来越受到人们的重视。
致病nsSNPs的识别通常被建模为一个二值分类问题。人们根据nsSNPs及其相应氨基酸置换的物理化学特性计算出数值特征，然后采用一定的分类方法预测出表征nsSNPs致病可能性的数值[6-15]，进而根据这一数值特征对候选nsSNPs进行分类。虽然这些方法获得了显着的成功，但是它们几乎都是基于通过蛋白质序列、结构、注释和进化信息计算得出的特征。由于已知的人类蛋白质结构信息非常有限，利用结构信息计算特征的方法在应用上受到很大的限制，因此，仅基于蛋白质序列提取的特征具有更广泛的应用范围。另一方面，已有方法大多采用较为传统的基于单一分类器的识别模型，包括Logistic回归[10]，贝叶斯模型[11]，神经网络[16]，决策树[17,18]，支持向量机[19,20]等等，而这些方法在对nsSNPs进行分类时的效果都不是很理想。
本文中，我们基于一组仅从蛋白质序列中提取的26维特征[21]对致病nsSNPs进行识别。我们使用5种常见的集成学习方法（AdaBoost[22]、LogitBoost[23]、随机森林[24]、L2boosting算法[25]和随机梯度回归算法[26]）和两种传统分类方法（决策树、支持向量机）对nsSNPs进行分类，并且对比了这些方法的性能。结果表明，基于集成学习方法的机器学习方法可以更有效地识别致病nsSNPs，其中又以LogitBoost的性能最佳。
和置换后的）。
2. 数据和方法
2.1 数据来源
Swiss-Prot数据库（57.13版，发布于2010年1月19日）收集了发生在12,189个人类蛋白质中的61,415个氨基酸置换。每个置换项都被标注为—疾病‖，—多态‖或—未分类‖。在本文中，我们将标注为—疾病‖的氨基酸置换称为致病变异，将标注为—多态‖的置换称为中性多态。Pfam 数据库[27]收集了蛋白质家族和多序列比对信息。该数据库第24.0版（2009年10月发布）包含了11,912个蛋白质家族的序列比对模型。在本文中，我们研究在Pfam数据库中至少有20个同源蛋白质的人类蛋白质，从而收集到发生在6,578个人类蛋白质上的15,150个疾病突变，13,558个中性多态。
2.2 特征集
我们采用了由Jiang et al. [21]提出的仅基于蛋白质序列计算的26维特征：包括基于氨基酸三个物理化学性质（分子量，pI值，疏水性）的特征，基于3个蛋白质二级结构（α-螺旋、β-折叠、β-转角）中氨基酸发生频率与已知的蛋白质二级结构信息计算所得的相对频率特征和两个进化保守性评分。前6个特征可以从文献中[28,29]或者只用蛋白质序列信息计算[30]得到。进化保守性评分可以定义为在Pfam多序列比对相应位置中氨基酸出现的频率（原始
对于查询蛋白质中给定的氨基酸置换对（原始->置换后），上述理化性质和相对频率可以由原始氨基酸（Org）和置换后的氨基酸（Sub）计算得到，也可由一个窗口内（W 是窗口的大小）查询蛋白质包括的所有氨基酸序列或查询蛋白质与其同源蛋白质序列比对后的列加权情况计算得到。在本文中，我们采用了与Jiang et al [21]相同设置的窗口参数，即窗口的大小为W=9。
综上所述，我们定义了一个包括24 个理化和相对频率特性（每六个氨基酸特征在四种不同情况下分别计算）和2 个保守性评分（原始氨基酸和置换后的氨基酸）的特征集。
2.3 分类方法
理想的分类方法应具有最好的预测能力和最少的计算复杂度。最近有关统计学习的研究表明，集成学习算法有很好的分类效果，并在数据挖掘、模式识别、文本分类、预测等方面得到了重要应用。与传统的机器学习方法相比，集成学习的方法不追求在一次分类中取得最优的结果，而是通过构建一系列分类效果仅好于随机猜测的—弱‖分类器，并将它们通过某种方式组合成一个—强‖分类器，使组合而成的—强‖分类器具有更好的分类效果。本文中我们选择了近几年很流行的5 种集成算法——AdaBoost 算法[22]，LogitBoost 算法[23]，随机森林算法[24]，L2boosting 算法[25]和随机梯度回归算法[26]——作为训练与预测的分类方法。
2.3.1 AdaBoost 算法
AdaBoost算法[23]，即自适应boosting算法，是由Freund和Schapire在1995年提出的现阶段最流行的boosting算法之一。这种算法的特点是需要迭代很多次，在每一次的迭代中，基于现有样本数据集构建一个基础的—弱‖分类器，反复调用该—弱‖分类器，通过对每轮中错判的样本赋予更大的权重，使其更关注那些难判的样本，经过多轮循环，最后采用加权的方法将各轮的—弱‖分类器合成—强‖分类器，从而得到较高精度的预测模型。
具体而言，我们定义1 1 ( , ),..., ( , ) N N x y x y 为训练集，其中向量i x 为特征值， i y 为第i个样本的分类，取值于{0，1}。我们定义( ) ( ) ( )1( ) ( )m M m mi m iF x c f x== ? ，其中( )( ) mi F x 是M 个弱分类器( )( ) mi f x 的加权集成， (m) c 为常系数。首先，我们对所有样本设定相同的初始权重(0) 1/ i w = N ，在每一轮迭代中，我们计算集成系数(m) c 并利用上一轮分类器的权重误差更新样本权重(m)i w 。随后，我们通过以下公式估计此轮中y =1 的概率并且在每一轮中更新它直到输出最后分类概率为止。最后，我们输出分类器( ( )) M sign F x 和分类概率( ) M p x 。
由于boosting是一种元算法，它可以与许多其他学习方法结合使用，从而提高其性能。
人们普遍认为，结合决策树和AdaBoost的算法（Breiman, 1998）是现在世界上最好的AdaBoost分类器。所以，本文中的AdaBoost算法采用了R软件中的决策树函数（tree）作为基础—弱‖分类器。
2.3.2 LogitBoost 算法
LogitBoost 算法是由Frieman、Hastie、Tibshirani 于2000 年在对AdaBoost 算法进行改进的基础上提出的。与Adaboost 算法相似，我们定义1 1 ( , ),..., ( , ) N N x y x y 为训练集，其中向量i x 为特征值， i y 为第i 个样本的分类，取值于{0，1}。通过公式(1) 我们估计每轮中y =1的概率( ) i p x 。我们同样定义( ) ( ) ( )1( ) ( )m M m mi m iF x c f x== ? ，其中( )( ) mi F x 是M 个弱分类器( )( ) mi f x 的加权集成， (m) c 为常系数，取值为1/2。LogitBoost 算法是利用自适应Newton步骤并根据最大化Bernoulli 似然函数来寻找合适的对数对称模型，并且在每次迭代中优化了关于f的函数E(-y(F(x) + f (x))) (Friedman et al., 1998)。具体来说，首先，我们对所有样本设定相同的初始权重(0) 1/ i w = N ，在每一轮迭代中，我们计算新的权重和工作响应。
随后，我们通过最小化加权的最小二乘寻找一个合适的回归函数(m) f 。最后，我们对集成函数和每次迭代的分类概率进行更新，并最终输出分类器( ( )) M sign F x 和分类概率( ) M p x 。
AdaBoost 算法与LogitBoost 算法的不同在于，AdaBoost 算法采用指数型损失函数，而LogitBoost 算法采用负对数似然损失函数[23]。相对于优化指数损失的Adaboost 算法，LogitBoost 算法更具有鲁棒性，可以显着地降低训练误差的同时提高泛化能力，并改善AdaBoost 算法存在的过拟合以及对噪声敏感等问题。LogitBoost 算法可以采用Logit 模型、Decision Stump 或者决策树作为基础分类器。
2.3.3 随机森林算法
随机森林算法[24]是Leo Breiman于2001年提出的，通过生成许多决策树并遵循少数服从多数的原则由他们—投票‖进行最终分类的一类集成学习算法。它具有需要调整的参数较少、不必担心过度拟合、分类速度很快，能高效处理大样本数据以及较强的抗噪音能力等特点，因此，在基因芯片数据挖掘、代谢途径分析及药物筛选等生物学领域得到应用并取得了较好的效果。此外，随机森林还有一个很大的优点，即在对数据进行分类的同时，可以给出各个特征在分类过程中的重要性评分，根据该评分能够筛选出相对重要的变量。
我们定义1 1 ( , ),..., ( , ) N N x y x y 为训练集，其中向量i x 为特征值， i y 为第i 个样本的分类，取值于{0，1}。随机森林算法的目标为产生M 棵决策树。在此算法中，每棵树的生长须符合以下规则：
1） 如果训练集中样本数为N，从原来数据中有替换的随机抽取N 个样本。这些样本是用来生成树的样本集。
2） 对于 M 个输入变量，在每个节点随机地选定m 个变量，m<<M，用能将这m 个变量进行最好分类的分类点作为分类的节点。在生成森林的过程中，m 的值保持不变。
3） 每一棵树都尽可能的延伸。没有剪枝。
而未被选定的样本称为袋外样本（OOB）。根据决策树算法，我们使用这种新的训练集最大限度的生成一棵没有任何剪枝的决策树，在每个节点只随机选取少量的变量（m 个）作为特征向量进行分类，而不是考虑所有的预测变量（随机特征选择）。重复该过程M 次，由此建立M 棵独立的决策树，从而形成一个随机生成的—森林‖。
2.3.4 L2boosting 算法
L2boosting算法是由Peter Buhlmann和Bin Yu于2003年首先提出的一种优化的boosting算法，其主要思想是：进行最小二乘拟合残差的迭代梯度boosting算法[25]。
我们定义1 1 ( , ),..., ( , ) N N x y x y 为训练集，其中向量i x 为特征值， i y 为第i个样本的分类，取值于{0，1}。定义F(x) = E(Y | X = x)为决策函数，损失函数为2 C(y, f ) = (y - f ) / 2。
L2boosting算法流程为：
1） （初始化）根据训练集合，我们利用最小二乘拟合出一个普通的梯度下降函数F0(x)作为初始的实值学习器。
2） 计算残差( ) i i m i U = Y - F X ，其中i=1,2…,n，m为boosting的次数。利用当前残差拟合实值学习器，并进行更新： 1 1 ( ) ( ) ( ) m m m F F f + + × = × + × ，其中1( ) m f + × 为梯度向量。
3） （迭代）增加一次迭代次数，使m=m+1，并重复第2步。
本文中，L2boosting算法是用R软件中mboost软件包的glmboost函数实现的。
2.3.5 随机梯度回归算法
随机梯度回归算法[26]是斯坦福大学教授Friedman于2001年提出的一种回归预测方法，也是一种Boosting算法，基础学习器一般选取回归树。借助梯度下降优化思想，利用损失函数负梯度产生的伪残差建立迭代回归树。按照Bagging的思想，在每次利用伪残差建立回归树的过程中不是使用所有的伪残差，而是随机选取一部分残差用来得到回归树。最后得到的模型就是一些回归树的线性组合。即其中 ( ) m T x 是回归树， m a 是线性组合的系数。这种方法具有回归树的统计特性，而且通常比单棵回归树预测的精度高[31]。R软件中mboost软件包的blackboost函数就是基于该思想编写的R代码。
2.4 验证方法和评估
标准基于从蛋白质序列中提取的26维特征，我们对这5种常见的集成学习方法（AdaBoost、LogitBoost、随机森林、L2boosting算法和随机梯度回归算法）和两种传统分类方法（决策树、支持向量机）在这一分类问题中的性能进行比较。验证实验时，我们分别用到了R软件中的Tree、e1071、randomForest、ada、caTools、mboost软件包，并参照软件包中使用手册的建议设置了具体参数，通过多次实验选择了能给予我们最好预测结果的参数值。
为了验证分类方法的稳定性及泛化能力，我们采用了10倍交叉验证方法(10-fold crossvalidation，10-cv)，具体做法是：将训练的28708个nsSNPs样本随机分为10组(每组约包含2871个样本)，然后每次留出1组作为测试数据，另9组作为训练数据，这样轮流进行10次实验，使得每组数据都恰能作为一次测试数据进行预测。最后，计算10次预测结果的平均值（方差）作为最终结果。
结合10倍交叉验证结果和确定的决策阈值，我们计算在该决策阈值下的真阳性（TP），真阴性（TN），假阳性（FP）和假阴性（FN），进而计算以下4个标准对分类方法的性能进行评价：预测的准确率（ACC），ROC曲线的面积（AUC），平衡错误率（BER），和马修斯相关系数（MCC）。其中，在给定决策阈值下的平衡误码率（BER）和马修斯相关系数（MCC）定义为：
这四个指标中，前两个指标为我们提供了对给定预测方法能力的全面评估，后两个指标考虑了两类样本的不平衡性，给我们详细地提供了该方法在给定的决策阈值下的预测能力。一般来说，BER越小，ACC、AUC和MCC越大，分类方法越好。
3. 结果
通过10 倍交叉验证，我们对5 种常见的集成学习方法（AdaBoost、LogitBoost、随机森林、L2boosting 算法和随机梯度回归算法）和两种传统方法（决策树、支持向量机）的分类性能进行了比较。结果表明，以上7 种方法都可以对nsSNPs 数据进行分类（1 为致病nsSNPs，0 为中性nsSNPs）。根据预测结果我们也可以判断该nsSNPs 为致病nsSNPs 的概率。
由于AdaBoost 算法、随机梯度回归算法、Logitboost 算法和随机森林（RF）算法均采用决策树算法作为—弱‖分类器或基础分类器，所以对于最终得到的集成—强‖分类器，其预测结果的四项指标均优于单独的决策树算法（如图1 所示）。对于分类效果较好的SVM 算法，它的ACC 可以达到0.7322、AUC 可以达到0.7927，而随机森林（RF）算法的ACC 可以达到0.7498、AUC 可以达到0.8057，Logitboost 算法的ACC 可以达到0.7565、AUC 可以达到0.8203（见表1）。由此可知，相比于SVM 算法，Logitboost 算法和随机森林（RF）算法的稳定性更好、准确率更高。因此，我们可以认为基于集成学习方法的机器学习方法可以更有效地识别致病nsSNPs。
同样的排序结果，即LogitBoost > 随机森林> SVM > AdaBoost >随机梯度回归> 决策树> L2boosting，其中—>‖代表性能优于。可以看出，Logitboost算法的ACC，AUC，MCC指标均明显高于其余方法，而Logitboost算法的BER指标明显低于其余方法。因此Logitboost算法具有7种分类方法中最好的准确率和平衡性。此外，由于ROC曲线可以全面的评价在不同阈值下各个分类方法的准确性，从图1的ROC曲线中我们可以进一步确认集成学习算法中的Logitboost算法性能最优，预测结果最准确。
4. 结论
本文中，我们将致病nsSNPs的识别问题建模为二值分类问题，基于从蛋白质序列中提取的26维特征，利用多种集成学习方法对nsSNPs进行分类，从而识别致病nsSNPs。我们的方法有以下几个优点。首先，识别致病nsSNPs将有助于人们理解人类疾病的遗传学基础，从而对这些疾病的预防、诊断、治疗起到一定的促进作用。第二，我们的方法不依赖于蛋白质结构的信息，它是从蛋白质序列和多序列比对中计算得出的简单的数值特征。第三，我们的方法采用了简单而强大的集成学习方法，较传统的基于单一机器学习算法的识别模型，它具有更高的预测准确性。最后，我们的方法可以结合大规模疾病项目中的病例对照数据。对于大规模数据处理和多基因疾病的问题，传统的统计方法需要做多假设检验。相比之下，我们的方法是基于纯粹的疾病突变的生物物理和生物化学信息以及相应蛋白质序列的进化信息，因此，它弥补了传统方法的一些不足之处。当然，我们的方法也可以按以下方式进行改善。首先，我们目前使用Pfam多序列比对方法提取查询蛋白质序列的保守蛋白质区域。因此，我们研究的突变受限于Pfam已知的蛋白质区域。此限制可以使用一些其他的多序列比对方法得以克服，如BLAST、PSI-BLAST和PANTHER。其次，现在大部分使用分类模型的研究（包括我们的）给出的预测结果是分级化的突变影响，这是由于大多数突变数据库中突变影响是按分类值给出的。作为一个长期的目标，我们应在预测基因突变对常见疾病造成的微小影响这一领域迈出重要的一步。其主要的困难是缺乏足够的样品来建立统计学习模型。当数据库资源充足时，目前的许多分类模型（包括我们的）可以延伸为回归模型，从而建立回归模型与微小突变影响之间的联系。最后，我们的方法仅限于识别基因编码区的致病nsSNPs。但是，在基因组其他区域的突变也可能会导致疾病，例如基因的转录因子结合位点和启动子区域，因此，我们需要进一步研究位于其它区域的基因突变。