CN109033747A - 一种基于pls多扰动集成基因选择及肿瘤特异基因子集的识别方法 - Google Patents

Abstract

本发明涉及一种基于PLS多扰动集成基因选择及肿瘤特异基因子集的识别方法，针对肿瘤微阵列数据的特点，引入不同的扰动机制，给出多扰动集成的基因选择的分析框架；利用PLS多基因度量方法，在该框架下发展出新的基于PLS集成基因选择方法。一方面，本发明的方法是基于子集的整体效应，能够快速识别出具有差异表达的基因，同时还能识别出差异表达信号微弱的基因；另一方面，本发明的方法是基于多重扰动机制，能够识别出一系列不同的，长度小且判别能力强的基因子集。因此，本发明的方法能够识别出的一系列不同的基因子集和弱差异表达基因，通过它们能够更加全面的认识肿瘤基因的特异表达模式。

Description

一种基于PLS多扰动集成基因选择及肿瘤特异基因子集的识 别方法

技术领域

本发明涉及计算科学与生命科学交叉技术领域，特别是一种基于PLS多扰动集成基因选择及肿瘤特异基因子集的识别方法。

背景技术

肿瘤是一种复杂的基因疾病，是由于某些染色体上DNA损伤而导致的细胞内基因异常表达，表现为细胞生长失控，缺乏分化和异常增生的一类复杂疾病。肿瘤基因微阵列(Microarray)能够从分子水平上探索和解释复杂多样的肿瘤疾病的发生与发展形成等过程。针对高通量基因表达谱数据，利用机器学习等技术能够识别出复杂的肿瘤疾病相关的特异基因及其功能，对研究肿瘤的疾病机理、预测肿瘤的疾病类型有着重要的意义。

肿瘤特异基因挖掘是为了寻找肿瘤相关的基因，目前已有方法通常直接或间接使用生物实验结合统计分析方法，来确定一组肿瘤相关的差异表达基因集合。事实上，这种基于统计方法的差异表达基因，对肿瘤的复杂遗传机制的探索是不充分的，或者说，这种所谓差异表达基因子集还不能捕获复杂遗传机制的全貌。

基于微阵列技术的生物学研究中，出现大规模、高通量的实验数据，首先，数据中包含成千上万的基因，表现为数据的维数非常高，相比较于数据的维数，所研究问题的样本容量相对较少，通常只有几十到上千；其次，这类具有高维数少样本特点的数据集中，存在着微弱但有效的基因差异表达信号；最后，由于数据规模通常比较大，很多机器学习所得模型过于复杂，计算所得结果难以从生物意义上解释。这些都对这类数据的探索性分析和信息的提取提出挑战。

在统计分析方法中，基因选择大多是基于单个基因进行疾病表型相关的分析和研究，常见的差异基因表达分析(difference expression,DE)有：倍数分析法(Ratio值)，参数分析法(调节后t检验)与方差分析(ANOVA)，非参数分析法(芯片显著性分析(SAM)与经验贝叶斯法(eBayes)、排序制造(RP)等。总体而言，这类方法都默认为每个基因对疾病表型的作用是相互独立的。在机器学习中，基因选择方法主要有基于Filter基因选择和基于Wrapper基因选择。Filter方法评估每个基因仍是依赖于该基因自身的固有信息，与其它基因的信息无关，其缺点是忽略了基因-基因的相关性。Wrapper方法以最大化分类器识别率为目标进行基因子集选择，其缺点是直接或间接的删除在生物学上可能极为重要的“冗余”基因。

利用基因表达谱识别不同肿瘤疾病状态之间发生差异表达的基因及其功能是研究肿瘤疾病发病机理的一项重要任务。肿瘤等复杂疾病是一种系统性疾病，其基因表达改变可能涉及广泛的功能。在实际的生物遗传过程中，通常普遍地存在着基因与表型、基因与基因复杂相关作用和交互作用。因此，更为合理和准确的分析方法是要同时考虑基因与表型、基因与基因之间的相关性以及交互作用。有效的肿瘤特异基因选择方法应该是能够考虑基因间的联合分布，即同时考虑所有的基因，允许检测那些具有较小主效应，但存在有较强交互效应的基因。实现对强主效应基因捕获的同时，也能捕获“较小主效应但较强交互效应”的基因。然而，传统的差异基因表达分析法与Filter方法，它们都可能会直接删除这些在生物学上极为重要的“无关”基因。另一方面，肿瘤基因芯片中有许多信息基因的表达谱会呈现高度相关，表现出具有共同表达某功能的基因组，而Wrapper方法通常直接或间接的删除这些在生物学上极为重要的“冗余”基因。事实上，这些冗余基因能够用来识别肿瘤组织中异常表达的基因，它们可能是生物体生化通路上的共调节基因，或者是在不同平行通道上的基因，对这些基因表达行为的分析能够更加全面认识肿瘤特异基因表达的特点。

传统基因选择方法的另一问题是选择结果的不稳定。由于数据的轻微扰动，或对同一个数据集采用不同的方法，基因选择的结果通常也不一致。于是，为了提高基因选择结果的稳定性和鲁棒性，集成学习的思想已逐渐被大家所共识。近年来，以多学习器集成的信息融合技术，在模式识别和机器学习领域中受到广泛关注。集成学习是一种机器学习范式，训练多个(同质或异质)学习器求解同一问题。集成学习算法对所训练的多个基学习器结果进行统计整合(如多数投票)。为了提升集成学习效果，增强成员学习器之间的差异性，不同的扰动机制被提出，包括有样本空间扰动，特征空间扰动以及模型参数的扰动等。通过多种扰动机制以加大学习器成员之间的差异性。

排序聚合(Ranking aggregation,RA)方法采用集成学习的思想，通过执行多次基因排序，把各个排序结果按照某一方式融合后选择基因子集，能够有效提高基因选择的稳定性和鲁棒性。排序聚合方法采用平均聚合等方式把基因的排名(Rank)或得分值(Score)融合在一起，并按照大小进行排序，得到稳定的基因排名，最后从中选择出若干个基因构成一组基因子集。使用合理的多基因度量方法时，从扰动所得的数据子集中，算法能够获得表达局部数据分布的相关基因，在整体数据分布中可能表现为弱差异表达的基因(或无关基因)。然而，这些基因通常能够增强最终子集整体的类判别能力。这类排序聚合方法主要针对基因子集的稳定性，算法最终目标是要获得一个稳定的基因子集。显然，所得的单个基因子集通常还不足以洞察复杂遗传机制的全貌。

发明内容

有鉴于此，本发明的目的是提出一种基于PLS多扰动集成基因选择及肿瘤特异基因子集的识别方法，利用该方法能够识别出一系列不同的基因子集和弱差异表达基因，通过这些不同基因子集以及弱差异表达基因能够更加全面的认识肿瘤基因的特异表达模式。

本发明采用以下方案实现：一种基于PLS多扰动集成基因选择及肿瘤特异基因子集的识别方法，包括以下步骤：

步骤S1：建立多响应变量PLS模型，使用SIMPLS算法求解多响应变量PLS模型，实现基于PLS多基因度量；

步骤S2：利用基于PLS多基因度量方法，在多扰动集成的基因选择的框架下，对样本数据进行基于PLS集成基因选择，得到样本数据的基因列表；

步骤S3：利用基分类器，从上述排序好的基因列表中识别出识别率最高的前k个基因，形成肿瘤特异基因子集。

本发明针对肿瘤微阵列数据的特点，引入不同的扰动机制，给出多扰动集成的基因选择的分析框架；利用PLS多基因度量方法，在该框架下发展出新的基于PLS集成基因选择方法。一方面，本发明的方法是基于子集的整体效应，能够快速识别出具有差异表达的基因，同时还能识别出差异表达信号微弱的基因；另一方面，本发明的方法是基于多重扰动机制，能够识别出一系列不同的，长度小且判别能力强的基因子集。在传统的差异表达基因识别过程中，通常没有考虑基因与表型、基因与基因之间的相关性及其交互作用，忽略了那些具有较小主效应，但存在有较强交互效应的基因。同时，针对有限样例的数据集，传统方法通常只能得到一组的差异表达基因。而通过本发明的方法能够识别出的一系列不同的基因子集和弱差异表达基因，这些不同的基因子集和弱差异表达基因能够更加全面的认识肿瘤基因的特异表达模式，对辅助理解肿瘤发生发展的机理、寻找肿瘤药物治疗的分子靶标，以及可靠的分子诊断与治疗，提供了新的视角。

进一步地，所述步骤S1中，建立多响应变量PLS模型具体为：

给定一个已知类别的样本数据集{(X_i,y_i)|X_i∈R^p,y_i∈Y_C,i＝1,…,n}，其中Y_C＝{c₁,c₂,…,c_g}是类标签集合，g表示类别的个数，每个样本的y_i对应样本可能所属的某个类别Y_C；

对原始的样本类标签(y)_n×1进行编码，定义模式分类PLS模型中的因变量为Y＝(y_ij)_n×g∈{0,1}^n×g，其中包括n个观测样本，g个类别，则

也即经编码后的因变量为矩阵Y＝(y_ij)_n×g；

通过利用类编码扩展方法，将原始的样本类标签(y)_n×1扩展成多响应变量矩阵Y＝(y_ij)_n×g，建立多响应变量PLS模型。

进一步地，所述步骤S1中，使用SIMPLS算法求解多响应变量PLS模型时，设置PLS模型中成分个数等于类别数g，以实现基于PLS多基因度量，具体实现步骤为：

步骤S11：输入n行p列的样本数据集X_n×p，以及对应的类标签列向量Y_n×1；

步骤S12：编码类标签列向量Y_n×1生成矩阵ClsY_n×g，其中g为类别数目；

步骤S13：调用算法SIMPLS(X,ClsY,g)，得到成分矩阵T和权重矩阵W，并计算方差解释量Rd；

步骤S14：设定变量i＝1,2,...p，由i＝1开始进行计算；

步骤S15：计算得到第i个基因的vip值，vip值即为该基因的权值，再令i＝i+1后再次进行计算；

步骤S16：所有计算完成后输出每个列的基因的权值。

进一步地，所述步骤S2中，对样本数据进行基于PLS集成基因选择具体包括以下步骤：

步骤S21：生成不同的数据子集：通过对样本数据集的随机扰动，生成B个不同的数据子集；

步骤S22：计算局部的基因权值：使用多基因度量方法，在每个数据子集上计算相应的局部基因权值；

步骤S23：集成所有的基因权值：对全部nB组局部基因权值依概率加权平均，返回达成共识的基因重要性列表。

进一步地，所述步骤S21具体为：对样本数据集X_n×p做bootstrap重抽样，生成B个样本子集，并且对每一样本子集的p个基因进行不放回随机抽样，抽样数为(取整)，最终生成B个的n行列的数据子集_i(1≤i≤B)，数据子集_i记X(n⁽ⁱ⁾,p⁽ⁱ⁾)；

其中，n⁽ⁱ⁾＝{n₁ ⁽ⁱ⁾,n₂ ⁽ⁱ⁾,...,n_n ⁽ⁱ⁾}表示子集_i的行标，也即n₁ ⁽ⁱ⁾,n₂ ⁽ⁱ⁾,...,n_n ⁽ⁱ⁾是从原始样本序号1,2,...,n中第i次bootstrap抽样的样本序号，并且1≤n₁ ⁽ⁱ⁾≤n₂ ⁽ⁱ⁾≤...≤n_n ⁽ⁱ⁾≤n；p⁽ⁱ⁾＝{p₁ ⁽ⁱ⁾,p₂ ⁽ⁱ⁾,...,p_M ⁽ⁱ⁾}表示子集_i的列标，也即p₁ ⁽ⁱ⁾,p₂ ⁽ⁱ⁾,...,p_M ⁽ⁱ⁾是从原始基因序号1,2,...,p中第i次的不放回随机抽样，并且1≤p₁ ⁽ⁱ⁾<p₂ ⁽ⁱ⁾<...<p_M ⁽ⁱ⁾≤p，其中抽样数M等于取整数。

进一步地，所述步骤S22具体为：对每个数据子集子集_i(1≤i≤B)计算基因权值，使用基于PLS多基因度量，求解在该子集中每个基因的权重值w_ij(1≤i≤B,1≤j≤p)，其中若则w_ij＝0，得到基因权重向量w_i＝{w_i1,w_i2,...,w_ip}。

进一步地，所述步骤S23具体为：计算所有基因的最终权值向量，即对B个数据子集求解得到B组权值向量后，按每个基因随机出现的次数求平均，对所有基因依平均权值降序排列，返回最后的集成的基因重要性列表。

进一步地，所述步骤S3中，识别肿瘤特异基因子集具体包括以下步骤：

步骤S31：根据步骤S2中对样本数据集中全部基因的排名，筛选出前k个基因；

步骤S32：在所选前k个基因构成的子集上，使用支持向量机为基分类器，对样本数据集进行模型训练；

步骤：利用模型训练得到的SVM分类模型，对测试样例在相应的k个基因上进行测试，计算识别率，记录测试识别率达到最大时的基因子集。

与现有技术相比，本发明有以下有益效果：本发明针对肿瘤微阵列数据的特点，引入不同的扰动机制，给出多扰动集成的基因选择的分析框架；利用PLS多基因度量方法，在该框架下发展出新的基于PLS集成基因选择方法。一方面，本发明的方法是基于子集的整体效应，能够快速识别出具有差异表达的基因，同时还能识别出差异表达信号微弱的基因；另一方面，本发明的方法是基于多重扰动机制，能够识别出一系列不同的，长度小且判别能力强的基因子集。因此，通过本发明的方法能够识别出的一系列不同的基因子集和弱差异表达基因，也能够更全面的认识肿瘤基因的特异表达模式，对辅助理解肿瘤发生发展的机理、寻找肿瘤药物治疗的分子靶标，以及可靠的分子诊断与治疗，提供了新的视角。

附图说明

图1是本发明实施例中基于多扰动集成的基因选择方法的原理示意图。

图2是本发明实施例中不同基因在两类样例(relapse和non-relapse)的表达情况示意图。

图3是本发明实施例中基因子集对不同样例的分离能力，也即在所选基因子集所张成的特征空间上，两类样例relapse和non-relapse的分布情况示意图。

图4是本发明实施例中弱差异表达基因的表达丰度图在癌组织和正常组织中的示意图。

图5是本发明实施例中弱差异表达基因的Kaplan-Meier生存曲线图。

具体实施方式

下面结合实施方式与实施例对本发明做进一步说明。

如图1所示，本实施例提供了一种基于PLS多扰动集成基因选择及肿瘤特异基因子集的识别方法，包括以下步骤：

步骤S1：建立多响应变量PLS模型，使用SIMPLS算法求解多响应变量PLS模型，实现基于PLS多基因度量；

步骤S2：利用基于PLS多基因度量方法，在多扰动集成的基因选择的框架下，对样本数据进行基于PLS集成基因选择，得到样本数据的基因列表；

步骤S3：利用基分类器，从上述排序好的基因列表中识别出识别率最高的前k个基因，形成肿瘤特异基因子集。

在本实施例中，基因选择是从原始的高维基因集合X＝{X¹,X²,L,X^p}中选择基因子集满足k＝p，其中p为原始基因的个数，k为所选基因的个数。基于机器学习的基因选择方法通常是对子集X^(k)进行整体评估，而传统的差异基因选择方法是对单个基因以统计显著性进行评估，观测p值是否显著。在肿瘤基因微阵列数据分析中，有效识别出肿瘤特异表达基因至关重要。在对基因的重要性进行度量时，本发明引入多基因度量的定义，并提出单基因度量和多基因度量。

单基因度量，是指第i个基因Xⁱ的权重值(或重要性得分)仅由单个基因Xⁱ自身决定，而与其它基因无关，也即假设基因之间存在独立性。多基因度量，是指第i个基因Xⁱ的权重值是由基因子集X^(k)决定，其中X^(k)是包含Xⁱ且长度为k(k>1)的基因子集。也即基因Xⁱ的权值是由k个(甚至k＝p)基因共同决定，而非仅仅单个基因Xⁱ决定其自身的权值，我们称这样的基因度量方法为多基因度量。

因此，把变量投影重要性指标(VIP)指标应用于基因权重值的计算中，给出类编码扩展方法，则所述步骤S1中，建立多响应变量PLS模型具体为：

给定一个已知类别的样本数据集{(X_i,y_i)|X_i∈R^p,y_i∈Y_C,i＝1,…,n}，其中Y_C＝{c₁,c₂,…,c_g}是类标签集合，g表示类别的个数，每个样本的y_i对应样本可能所属的某个类别Y_C；

对原始的样本类标签(y)_n×1进行编码，定义模式分类PLS模型中的因变量为Y＝(y_ij)_n×g∈{0,1}^n×g，其中包括n个观测样本，g个类别，则

也即经编码后的因变量为矩阵Y＝(y_ij)_n×g；

通过利用类编码扩展方法，将原始的样本类标签(y)_n×1扩展成多响应变量矩阵Y＝(y_ij)_n×g，建立多响应变量PLS模型。

在本实施例中，所述步骤S1中，使用SIMPLS算法求解多响应变量PLS模型时，设置PLS模型中成分个数等于类别数g，以实现基于PLS多基因度量，具体实现步骤为：

步骤S11：输入n行p列的样本数据集X_n×p，以及对应的类标签列向量Y_n×1；

步骤S12：编码类标签列向量Y_n×1生成矩阵ClsY_n×g，其中g为类别数目；

步骤S13：调用算法SIMPLS(X,ClsY,g)，得到成分矩阵T和权重矩阵W，并计算方差解释量Rd；

步骤S14：设定变量i＝1,2,...p，由i＝1开始进行计算；

步骤S15：计算得到第i个基因的vip值，vip值即为该基因的权值，再令i＝i+1后再次进行计算；

步骤S16：所有计算完成后输出每个基因的权值。

根据上述步骤，以下为基于PLS多基因度量方法(PLSMGM)在程序中的具体实现步骤：

输入：n行p列数据集X_n×p，及对应的类标签列向量Y_n×1

输出：每个基因的权值得分Score

(1):初始化：编码类标签列向量Y_n×1生成矩阵ClsY_n×g，其中g为类别数目(2):调用函数SIMPLS(X,ClsY,g)得到成分矩阵T和权重矩阵W，计算方差解释量RD

(3):For i＝1,2,...p do

(4):计算得到第i个基因的vip值

(5):Score(i)＝vip

(6):EndFor

(7):返回Score

在本实施例中，利用PLS多基因度量方法(PLSMGM)，在多扰动集成的基因选择的框架下，给出基于PLS集成基因选择方法(PLSEGS)，则所述步骤S2中，对样本数据进行基于PLS集成基因选择具体包括以下步骤：

步骤S21：生成不同的数据子集：通过对样本数据集的随机扰动，生成B个不同的数据子集；

步骤S22：计算局部的基因权值：使用多基因度量方法，在每个数据子集上计算相应的局部基因权值；

步骤S23：集成所有的基因权值：对全部B组局部基因权值依概率加权平均，返回达成共识的基因重要性列表。

在本实施例中，所述步骤S21具体为：对样本数据集X_n×p做bootstrap重抽样，生成B个样本子集，并且对每一样本子集的p个基因进行不放回随机抽样，抽样数为(取整)，最终生成B个的n行列的数据子集_i(1≤i≤B)，数据子集_i记为X(n⁽ⁱ⁾,p⁽ⁱ⁾)；

其中，n⁽ⁱ⁾＝{n₁ ⁽ⁱ⁾,n₂ ⁽ⁱ⁾,...,n_n ⁽ⁱ⁾}表示子集_i的行标，也即n₁ ⁽ⁱ⁾,n₂ ⁽ⁱ⁾,...,n_n ⁽ⁱ⁾是从原始样本序号1,2,...,n中第i次bootstrap抽样的样本序号，并且1≤n₁ ⁽ⁱ⁾≤n₂ ⁽ⁱ⁾≤...≤n_n ⁽ⁱ⁾≤n；p⁽ⁱ⁾＝{p₁ ⁽ⁱ⁾,p₂ ⁽ⁱ⁾,...,p_M ⁽ⁱ⁾}表示子集_i的列标，也即p₁ ⁽ⁱ⁾,p₂ ⁽ⁱ⁾,...,p_M ⁽ⁱ⁾是从原始基因序号1,2,...,p中第i次的不放回随机抽样，并且1≤p₁ ⁽ⁱ⁾<p₂ ⁽ⁱ⁾<...<p_M ⁽ⁱ⁾≤p，其中抽样数M等于取整数。

在本实施例中，所述步骤S22具体为：对每个数据子集子集_i(1≤i≤B)计算基因权值，使用基于PLS多基因度量方法(PLSMGM)，求解在该子集中每个基因的权重值w_ij(1≤i≤B,1≤j≤p)，其中若则w_ij＝0，得到基因权重向量w_i＝{w_i1,w_i2,...,w_ip}。

在本实施例中，所述步骤S23具体为：计算所有基因的最终权值向量，即对B个数据子集求解得到B组权值向量后，按每个基因随机出现的次数求平均，对所有基因依平均权值降序排列，返回最后的集成的基因重要性列表。

根据上述步骤，以下给出基于PLS集成基因选择(PLSEGS)在程序中的具体实现步骤：

输入：n行p列数据集X_n×p，及对应的类标签列向量Y_n×1，抽样次数B

输出：基因重要性排名Ranked

(01)：初始化：[n,p]＝size(X)，基因的权重向量W＝[0,0,...,0]选中的基因个数向量nsel＝[0,0,...,0]

(02)：For i＝1,2,...B do

(03)：1到p中采样个基因，序号为g_sel

(04)：s_sel＝bootstrap(1:n)，1到n中bootstrap采样n个样本，序号为g_sel

(05)：sub_X＝X(s_sel,g_sel)，生成X的数据子集sub_X

(06)：sub_Y＝Y(s_sel,:)，生成对应于sub_X的类标签向量sub_Y

(07)：loc_w＝PLSMGM(sub_X,sub_Y)，调用PLSMGM生成基因的权重值

(08)：W(:,g_sel)＝W(:,g_sel)+loc_w，更新权重值W在基因位置g_sel上

(09)：nsel(:,g_sel)＝nsel(:,g_sel)+1，更新选中的基因个数nsel

(10)：EndFor

(11)：ens_w＝W/nsel

(12)：从大到小排序ens_w，记录相应的排名位置Ranked

(13)：返回Ranked

在本实施例中，根据最新研究可知，导致肿瘤发生的基因变异数量通常介于1到10之间，研究人员还发现，一些关键性致癌突变几乎有一半发生在尚未被确定与肿瘤相关的基因中，这些也说明更多致癌基因仍有待被发现；

因此，针对有限样例的肿瘤基因微阵列数据集，利用机器学习方法，挖掘出一系列的不同的特异基因子集，这些子集仅含有少数基因(设置不超过10)，且这些子集能实现对测试样例集的较高识别(甚至100％识别)，这些子集中除了关键性的差异表达基因外，还存在一些重要的弱差异表达基因，则所述步骤S3中，识别肿瘤特异基因子集具体包括以下步骤：

步骤S31：根据步骤S2中对样本数据集中全部基因的排名，筛选出前k个基因；

步骤S32：在所选前k个基因构成的子集上，使用支持向量机为基分类器，对样本数据集进行模型训练；

步骤S33：利用模型训练得到的SVM分类模型，对测试样例在相应的k个基因上进行测试，计算识别率，记录测试识别率达到最大时的基因子集。

该步骤在程序中的详细的计算过程如下：

步骤1:初始化：nRun＝0,k＝0,max_nRun＝5000,max_k＝10,B＝2000

步骤2:调用本发明的方法PLSEGS，在训练集上计算每个基因的排名

步骤3:选取前k个基因，在训练集上训练SVM分类模型

步骤4:使用分类模型在测试集的对应k个基因上进行分类，计算识别率

步骤5:k＝k+1，如果k<max_k，跳转到(步骤3)

步骤6:nRun＝nRun+1；如果nRun≤max_nRun，重复(步骤2)到(步骤5)

步骤7:返回nRun个子集，对应于每次内部循环中识别率达到最高的前k个基因。

综上所述，本实施例针对肿瘤微阵列数据的特点，引入不同的扰动机制，给出多扰动集成的基因选择的分析框架；利用PLS多基因度量方法，在该框架下发展出新的基于PLS集成基因选择方法。一方面，该方法是基于子集的整体效应，能够快速识别出具有差异表达的基因，同时还能识别出差异表达信号微弱的基因；另一方面，该方法是基于多重扰动机制，能够识别出一系列不同的，长度小且判别能力强的基因子集。通过该方法能够识别出的一系列不同的基因子集和弱差异表达基因，也能够更全面的认识肿瘤基因的特异表达模式，对辅助理解肿瘤发生发展的机理、寻找肿瘤药物治疗的分子靶标，以及可靠的分子诊断与治疗，提供了新的视角。

在本实施例中，为了验证该识别方法的有益效果，选取两个人类肿瘤微阵列数据集。数据集包含有限样例，但基因数量较多，属于典型的高维小样本数据。在生物信息学、模式识别和机器学习等领域中它们被大量研究，经常被用来验证基因选择方法和机器学习的性能，其中所选的乳腺癌数据集是公认的较难分类判别的数据。相关数据来源及其背景描述如下：

乳腺癌是影响女性健康的最主要的恶性肿瘤之一，最新研究报告显示，全球肿瘤5年生存率在稳步提升，乳腺癌的5年生存率超过80％，但不同地区患者的生存情况仍然不容乐观。这里所使用的乳腺癌数据集共包含24481个基因和97个病例样本的信息，其中46个样本在初次诊断治疗后5年内癌症发生转移(数据集中被标记为“relapse”)，而剩余的51个病例(标记为“non-relapse”的样本)初次治疗后至少在5年时间间隔内仍然保持健康的状态。该数据集分成训练样本集和测试样本集，训练样本集包括78例训练样本(38例relapse+44例non-relapse)，测试样本集包括19例测试样本(12例relapse+7例non-relapse)。

前列腺癌是世界范围内高发和主要的致死癌症之一，前列腺癌具有遗传易感性。这里所使用的前列腺癌数据集主要用于预测前列腺癌的临床诊断结果。它总共收集136例前列腺病例样本信息，每个样本的表达基因数目为12600。其中75例为前列腺癌肿瘤样本(PTS:prostate tumor sample)，59例为正常前列腺组织(NPS:normal prostate sample)。该数据集分成训练样本集和测试样本集，训练样本集包括102例训练样本(52例PTS+50例NPS)，测试样本集包括34例测试样本(25例PTS+9例NPS)。

采用以上所述样本进行后续测试，具体分为以下几点：

1)基于不同分类器的结果比较

最新研究表明，导致肿瘤发生的基因变异数量通常介于1到10之间，研究人员还发现，一些关键性致癌突变几乎有一半发生在尚未被确定与肿瘤相关的基因中，这些也说明更多致癌基因仍有待被发现。以下实验主要考查本发明的方法与不同基分类器的结合使用，并以基分类器对测试集的最大识别率(这里识别率被设置为100％)为目标，考查该方法所选出包含基因个数较少的基因子集。

在对本发明的集成基因选择方法进行评价过程中，仅仅关注本发明的方法对不同基因子集的总体判别能力，使用5个不同类型的分类器，支持向量机SVM(线性核)，Fisher线性判别法(FDA)，K近邻分类器(K＝10)以及最近邻分类器(1NN)，朴素Bayesian分类器(NBC)对测试集数据进行分类识别。以下是利用这些不同的分类器，本发明的方法在2个肿瘤数据上挖掘出长度最小的基因子集，表1和表2给出具有最小长度的基因子集，表中阴影标注的基因是弱差异表达基因。

表1.本发明方法所筛选具有最小长度的基因子集，使用不同基分类器，乳腺癌的测试数据集在这些子集上均能够被正确识别(识别率100％)

表2.本发明方法所筛选具有最小长度的基因子集，使用不同基分类器，前列腺癌的测试数据集在这些子集上均能够被正确识别(识别率100％)

从表1和表2容易看出，在乳腺癌数据集上，本发明的方法能够识别出不同基因子集，分类器SVC、FDA和KNN(K＝10)只需要三个基因就能够完全识别测试集，同样，在分类器NBC和1NN上只需要四个基因也能完全识别测试集。类似地，在前列腺癌数据集上，算法能够识别出不同基因子集，所有基分类器都只需要二个基因就能够完全识别测试集。这些均表明，本发明方法能够挖掘出具有强分离能力的长度较小的不同的基因子集。同时，在这些具有较高识别率的基因子集中，还分布着不同的弱差异表达基因。这是由于本发明方法是基于子集的整体效应，所挖掘出的基因间存在有交互效应。从而使得该方法能够识别出这些弱差异表达基因，这些弱差异表达基因将有助于从整体上去理解基因之间的特异表达模式。

2)和非集成方法的比较分析

这里的实验主要考查传统非集成基因选择方法与本发明的集成基因选择方法之间的性能差异。通过分类器的不同性能指标，间接刻画集成与非集成基因选择方法之间的差异。对于非集成基因选择方法，选择6个不同的经典方法，包括单变量方法和多变量方法，其中单变量方法有：基于参数统计检验T-test和非参数统计检验Ranksum；多变量方法有：基于集成的随机森林、基于近邻信息的ReliefF和基于递归特征消除的PLSRFE和SVMRFE。利用这些不同的基因选择方法，在2个肿瘤数据上筛选长度不超过100的基因子集，表3列出所有方法在最大识别率时相关的指标信息。

表3.不同基因选择方法的性能比较(非集成方法vs.集成方法)

注：基分类器为线性核支持向量分类器(linear SVC)

从表3容易看出，在乳腺癌数据集上，本发明的集成基因选择方法能够选择出识别率100％的仅包含三个基因的子集，并且集成方法所生成的不同子集，它们的性能指标的中数也优于所有非集成方法的指标值。同样的，在前列腺癌数据集上，本发明的集成基因选择方法能够选择出识别率100％的仅包含二个基因的子集，并且集成方法所生成的不同子集，它们的性能指标的中数也达到非集成方法中的最大指标值，但集成方法却有更小的子集长度。这些也都表明本发明的集成基因选择方法能够生成具有更强判别能力的基因子集。

3)肿瘤特异基因子集的识别与分析

评价方法优劣还要视其结果是否具有生物学含义。以下仅针对乳腺癌数据集，应用生物信息学并结合文献检索，对本发明的集成基因选择方法所识别不同基因子集做进一步的验证分析。

(一)特异基因及基因子集的类分离能力

以下给出本发明方法所识别出的长度小于10的不同基因子集，所有这些子集在测试集上均实现对两类样例的正确分类(识别率100％)，也即这些基因子集对relapses和non-relapsed两类样例具有极强的判别能力。从表4容易看出，本发明方法能够挖掘出众多不同的基因子集，并且这些子集中存在有冗余基因，也即每个基因子集中几乎都包含有强差异表达基因Feat#10889(TSPYL5)和Feat#13800(ATP5E)。同时，表4中不同基因子集还存在弱差异表达基因(阴影标注)。

表4.长度不超过10的不同的基因子集，它们在测试数据集上实现100％识别

注：阴影标记的基因为弱差异表达基因，也即没有通过统计显著性检验。(基分类器：LinearSVC)

以下观测基因子集中每个基因在两类样例上的表达情况，以及所选基因子集整体的类分离能力。这里仅对表4中的第二组基因子集做可视化分析。如图2给出两类样例(relapse和non-relapse)在3个基因上的表达值分布，左边为训练集，右边为测试集，从图中容易看出，在训练集和测试集上，前两个基因TSPYL5和ATP5E在两类样例上存在显著的差异表达，也即它们是所谓的强差异表达基因。而基因STK3的散点图(尤其在测试集上)显示，其在两类样例上并不存在差异表达，也即该基因为弱差异表达基因。

更进一步观测这三个基因所组成的子集的整体效应，其对两类样例的分离能力。如图3所示，首先观测两个强差异表达基因子集对两类样例的分离能力，图3的上半部分为两类样例在这两个基因(Feat#10889和Feat#13800)上的分布情况，容易看出，对于训练集和测试集，它们均未能实现线性可分离。图3的下半部分显示，在引入弱差异表达基因Feat#9739(STK3)后，该子集在训练集上两类样例的分离性有所提升，特别在测试集上，可以实现对两类样例的线性分离。这间接表明，这个弱差异表达基因和其它两个基因存在有交互效应，从而使得整个基因子集对类别的分离能力表现出更好的整体效应。

(二)弱差异表达基因及其功能分析

针对本发明所识别出的弱差异表达基因进行分析，表5给出这些基因的相关信息，具体包括基因的序号、存取号、基因名称及其功能描述、基因在两类样例中的表达差异的显著性p值(这些基因在两类样例上均未能通过统计显著性检验，显著水平为0.01)。

表5.本发明所识别出的部分的弱差异表达基因的相关信息的描述

在对基因微阵列数据分析过程中，不同的数据分析方法对结果有显著影响。强差异表达基因通常易于实验验证，然而基于不同统计学方法以及不同因素所得到的弱差异表达基因)通常难以验证。因此，为了说明本发明方法所得弱差异表达基因的有效性，利用TCGA数据库并结合PubMed数据库做进一步的验证，首先在TCGA数据库上对这些弱差异基因进行分析，图4给出它们在癌组织和正常组织上的表达丰度图。容易看出，除了基因APOLD1外，其它基因在两类组织中的表达丰度都没有显著差异，如基因ERI1在癌组织和正常组织上的表达丰度都很小，且丰度值几乎一致。事实上，这类弱差异表达基因很少被关注，传统实验方法通常是直接忽略它们，而更加专注于那些所谓的强差异表达基因的研究。

其次，在PubMed数据库结合Mesh，KEGG等数据库对这些弱差异表达基因进行分析，表6给出这些基因在以上数据库中的相关的统计汇总等信息，具体包括：乳腺癌相关文献中出现该基因的文献数、对应乳腺癌疾病相关化合物中该基因的文献数以及该基因本身所包含的致病突变数。表6容易看出，它们中间有6个基因在PubMed等数据库中有相应的研究及文献记载，而另外6个基因(阴影标注)在PubMed等数据库中还没有其关于乳腺癌的直接报告。

表6.弱差异表达基因在PubMed数据库中的相关信息

进一步通过PubMed文献检索以及Mesh，KEGG等数据库，我们发现这些无正面报道的弱差异表达基因，它们存在与其它相关联的基因，如表7所示，这些与其关联的其它基因，在乳腺癌研究文献中却存在有相对应的研究报道，这也间接表明这些弱差异表达基因和乳腺癌并非毫无关系，它们可能存在有对乳腺癌的间接表达，可能是乳腺癌的潜在基因。

表7.弱差异表达基因与其它基因间的相互作用关系及其证据来源.

最后，使用在线分析工具：www.kmplot.com/mirpower，分析了这6个无正面报道的弱差异表达基因，图5给出生存分析中的Kaplan-Meier图，它描述了不同基因的表达对乳腺癌患者生存概率随时间的变化，它能够很好地描述生存过程。容易看出这6个基因中的5个基因(除了基因DCAF4)的对数秩(Log-rank)的p值小于0.05，表明这些基因对乳腺癌生存有显著差异。

因此，本发明的方法所识别这些基因子集可能与乳腺癌预后相关，特别的，通过查找相关文献，发现12个弱差异表达基因中有6个基因已被证实与乳腺癌相关，它们在PubMed等数据库中存在有相关研究报道，而另外6个基因尽管在PubMed等数据库中没有直接的研究报道，但它们和其它基因存在关联作用，对乳腺癌存在间接表达。并且生存分析的Kaplan-Meier图也表明这些弱差异表达基因对乳腺癌的生存率有关，它们也可能是乳腺癌病人的预后标志基因。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种基于PLS多扰动集成基因选择及肿瘤特异基因子集的识别方法，其特征在于：包括以下步骤：

步骤S1：建立多响应变量PLS模型，使用SIMPLS算法求解多响应变量PLS模型，实现基于PLS多基因度量；

步骤S2：利用基于PLS多基因度量方法，在多扰动集成的基因选择的框架下，对样本数据进行基于PLS集成基因选择，得到样本数据的基因列表；

步骤S3：利用基分类器，从上述排序好的基因列表中识别出识别率最高的前k个基因，形成肿瘤特异基因子集。

2.根据权利要求1所述的一种基于PLS多扰动集成基因选择及肿瘤特异基因子集的识别方法，其特征在于：所述步骤S1中，建立多响应变量PLS模型具体为：

给定一个已知类别的样本数据集{(X_i,y_i)|X_i∈R^p,y_i∈Y_C,i＝1,…,n}，其中Y_C＝{c₁,c₂,…,c_g}是类标签集合，g表示类别的个数，每个样本的y_i对应样本可能所属的某个类别Y_C；

对原始的样本类标签(y)_n×1进行编码，定义模式分类PLS模型中的因变量为Y＝(y_ij)_n×g∈{0,1}^n×g，其中包括n个观测样本，g个类别，则

也即经编码后的因变量为矩阵Y＝(y_ij)_n×g；

通过利用类编码扩展方法，将原始的样本类标签(y)_n×1扩展成多响应变量矩阵Y＝(y_ij)_n×g，建立多响应变量PLS模型。

3.根据权利要求1所述的一种基于PLS多扰动集成基因选择及肿瘤特异基因子集的识别方法，其特征在于：所述步骤S1中，使用SIMPLS算法求解多响应变量PLS模型时，设置PLS模型中成分个数等于类别数g，以实现基于PLS多基因度量，具体实现步骤为：

步骤S11：输入n行p列的样本数据集X_n×p，以及对应的类标签列向量Y_n×1；

步骤S12：编码类标签列向量Y_n×1生成矩阵ClsY_n×g，其中g为类别数目；

步骤S13：调用算法SIMPLS(X,ClsY,g)，得到成分矩阵T和权重矩阵W，并计算方差解释量Rd；

步骤S14：设定变量i＝1,2,...p，由i＝1开始进行计算；

步骤S15：计算得到第i个基因的vip值，vip值即为该基因的权值，再令i＝i+1后再次进行计算；

步骤S16：所有计算完成后输出每个基因的权值。

4.根据权利要求1所述的一种基于PLS多扰动集成基因选择及肿瘤特异基因子集的识别方法，其特征在于：所述步骤S2中，对样本数据进行基于PLS集成基因选择具体包括以下步骤：

步骤S21：生成不同的数据子集：通过对样本数据集的随机扰动，生成B个不同的数据子集；

步骤S22：计算局部的基因权值：使用多基因度量方法，在每个数据子集上计算相应的局部基因权值；

步骤S23：集成所有的基因权值：对全部B组局部基因权值依概率加权平均，返回达成共识的基因重要性列表。

5.根据权利要求4所述的一种基于PLS多扰动集成基因选择及肿瘤特异基因子集的识别方法，其特征在于：所述步骤S21具体为：对样本数据集X_n×p做bootstrap重抽样，生成B个样本子集，并且对每一样本子集的p个基因进行不放回随机抽样，抽样数为(取整)，最终生成B个的n行列的数据子集_i(1≤i≤B)，数据子集_i记为X(n⁽ⁱ⁾,p⁽ⁱ⁾)；

其中，n⁽ⁱ⁾＝{n₁ ⁽ⁱ⁾,n₂ ⁽ⁱ⁾,...,n_n ⁽ⁱ⁾}表示子集_i的行标，也即n₁ ⁽ⁱ⁾,n₂ ⁽ⁱ⁾,...,n_n ⁽ⁱ⁾是从原始样本序号1,2,...,n中第i次bootstrap抽样的样本序号，并且1≤n₁ ⁽ⁱ⁾≤n₂ ⁽ⁱ⁾≤...≤n_n ⁽ⁱ⁾≤n；p⁽ⁱ⁾＝{p₁ ⁽ⁱ⁾,p₂ ⁽ⁱ⁾,...,p_M ⁽ⁱ⁾}表示子集_i的列标，也即p₁ ⁽ⁱ⁾,p₂ ⁽ⁱ⁾,...,p_M ⁽ⁱ⁾是从原始基因序号1,2,...,p中第i次的不放回随机抽样，并且1≤p₁ ⁽ⁱ⁾<p₂ ⁽ⁱ⁾<...<p_M ⁽ⁱ⁾≤p，其中抽样数M等于取整数。

6.根据权利要求4所述的一种基于PLS多扰动集成基因选择及肿瘤特异基因子集的识别方法，其特征在于：所述步骤S22具体为：对每个数据子集子集_i(1≤i≤B)计算基因权值，使用基于PLS多基因度量，求解在该子集中每个基因的权重值w_ij(1≤i≤B,1≤j≤p)，其中若则w_ij＝0，得到基因权重向量w_i＝{w_i1,w_i2,...,w_ip}。

7.根据权利要求4所述的一种基于PLS多扰动集成基因选择及肿瘤特异基因子集的识别方法，其特征在于：所述步骤S23具体为：计算所有基因的最终权值向量，即对B个数据子集求解得到B组权值向量后，按每个基因随机出现的次数求平均，对所有基因依平均权值降序排列，返回最后的集成的基因重要性列表。

8.根据权利要求1所述的一种基于PLS多扰动集成基因选择及肿瘤特异基因子集的识别方法，其特征在于：所述步骤S3中，识别肿瘤特异基因子集具体包括以下步骤：

步骤S31：根据步骤S2中对样本数据集中全部基因的排名，筛选出前k个基因；

步骤S32：在所选前k个基因构成的子集上，使用支持向量机为基分类器，对样本数据集进行模型训练；

步骤S33：利用模型训练得到的SVM分类模型，对测试样例在相应的k个基因上进行测试，计算识别率，记录测试识别率达到最大时的基因子集。