0 概述

人的精神状况、行为偏好由先天基因组与外界刺激共同决定^[1-2]。基因检测技术的发展使得从全基因组表达水平定量检测基因转录产物mRNA得以实现。“基因→mRNA→蛋白质”为基因表达过程, 可以通过分析mRNA数据来研究发生变化的基因表达。运用机器学习方法对基因表达数据进行研究, 对医学临床诊断具有重要的参考意义。但不同于普通的分类问题, 由于伦理和法律方面的约束与潜在风险, 在医学诊断中对基因表达数据进行研究的目标是得到精确且可解释的分类机制。机器学习领域主流的SVM^[3]、神经网络^[4]等方法对数据具有很强的拟合能力, 但其结果的可解释性较差, 应用于医学诊断领域存在一定局限。

针对上述问题, 有学者将关联规则引入到分类问题中。文献[5]提出CBA算法, 其挖掘基因数据间的关联规则并用于分类。基于关联的分类规则具有直观的可解释性, 但其在高维基因数据上会产生大量的冗余规则, 且其规则挖掘与分类器构建为两阶段算法, 导致即使有巨量规则, 仍会有大量样本未被分类器所覆盖。文献[6]从样本覆盖的角度出发, 提出一种Top-k算法, 其以样本枚举的方式挖掘可以交叉覆盖到全部样本的多个规则组, 以进行分类器构建。以上算法在数据预处理阶段需将连续的基因表达水平值离散化, 即根据阈值的大小将基因表达水平值划分为激活态和抑制态。跨平台基因表达数据的特点是各平台的样本具有相同的特征和标签, 但由于平台的特性, 数据值之间存在尺度差异, 因此不同平台间的阈值也不同, 导致已有分类模式只适用于单一平台。

受k-TSP分类模式中“相对表达反转”概念^[7]的启发, 本文设计一种层级规则树(Hierarchy Rule Tree, HRT)分类模式。传统有监督学习算法推广到基因表达数据领域时, 一般通过基因选择和关联规则发现2种方式^[8-10]。HRT分类模式将上述2种方式进行结合, 基于基因间相对表达反转的特性, 设计相应的数据转换和规则预筛选策略, 进行迭代式分类规则挖掘, 以提高跨平台基因表达数据挖掘的效率。

1 跨平台分类模式 1.1 最高得分对

文献[11]提出一种最高得分对(Top Scoring Pair, TSP)算法, 该算法是纯数据驱动的机器学习方法, 无需使用者设定参数, 避免了精确调整参数所带来的过拟合问题。TSP分类器的出发点是不同基因表达水平之间的相对差异性, 即寻找事件{V_a < V_b}(V_a(V_b)代表基因a(b)的表达值)在不同类样本中发生概率显著不同的那些标志性基因对(a, b)。TSP算法关注事件{V_a < V_b}在每一类样本中发生的概率p_ab(C_m), 其公式定义如下:

$ {P_{ab}}\left( {{C_m}} \right) = \mathit{Prob}\left\{ {{V_a} < {V_b}\left| {Y = {C_m}} \right.} \right\} $

其中, C_m为样本类别, 此处仅考虑二分类问题, m={0, 1}。用Δ_ab代表每一个基因对(a, b)的得分, Δ_ab=|p_ab(C₀)－p_ab(C₁)|, 以得分最高的基因对作为分类器的判定依据。基因数据具有高维度特性, 上万个基因组合出的基因对的数量是亿级, 因此, 有可能出现得分一致的基因对。如果出现这种情况, TSP算法会计算第2个评分指标, 即平均排序差, 以选取唯一的最高得分对。

TSP算法理论上可以作为跨平台分类模式的基础, 但以单个基因对作为全部样本的分类模式过于简单。有研究者基于机器学习集成思想, 以TSP作为基分类器, 提出改进的k-TSP算法, 但仍难以拟合复杂的数据分布。

1.2 层级规则树

从规则构成的角度出发, 以基因对(a, b)之间表达值反转为基础的分类模式, 是由2个独立的分类规则$\left\{ { < a, b > \to {C_m}} \right\} 、\left\{ {\neg <a, b > \to {C_m}} \right\}(m = \{ 0, 1\} ) $所构成, 其中, ${ < a, b > } $代表{V_a < V_b}, $ {\neg < a, b > }$代表{V_a≥V_b}。

本文提出的层级规则树以支持度(supp)和置信度(conf)作为规则兴趣度的度量。支持度和置信度是传统关联规则挖掘中的度量标准, 现已扩展到分类规则挖掘中。给定一个数据集R, d为数据集样本, X和Y为样本属性, 对于分类规则X→Y, supp和conf公式定义分别如下:

$ \mathit{supp}\left( {X \to Y} \right) = \frac{{\left| {\left\{ {d \in R;X \subseteq d,Y \subseteq d} \right\}} \right|}}{{\left| R \right|}} $

(1)

$ conf\left( {X \to Y} \right) = \frac{{\left| {\left\{ {d \in R;X \subseteq d,Y \subseteq d} \right\}} \right|}}{{\left| {\left\{ {d \in R;X \subseteq d} \right\}} \right|}} $

(2)

以表 1中的数据分布为例, 有:

$ \mathit{supp}\left( { < a,b > \to {C_0}} \right) = \frac{{{n_1}}}{{{n_1} + {n_2} + {n_3} + {n_4}}} $

$ \mathit{supp}\left( { < a,b > \to {C_1}} \right) = \frac{{{n_3}}}{{{n_1} + {n_2} + {n_3} + {n_4}}} $

$ conf\left( { < a,b > \to {C_0}} \right) = \frac{{{n_1}}}{{{n_1} + {n_2}}} $

$ \mathit{conf}\left( { < a,b > \to {C_1}} \right) = \frac{{{n_3}}}{{{n_3} + {n_4}}} $

可以看出, $ < a, b > 、\neg < a, b > $这2个分类规则基于完全不同的样本子集而得出, 两者之间并不存在关联关系。有可能其中一个规则具有很强的分类能力, 而另外一个规则没有参考价值。将两者相结合进行衡量, 会限制分类器对数据的拟合能力。基于此, 本文以其中的单个规则为基础, 建立一种层级规则树分类模式。

层级规则树首先对置信度最高的规则特征$ < a, b > $所覆盖的样本进行有效分类, 对于剩余样本(特征$\neg < a, b > $所覆盖的样本), 寻找更精确的新规则(如$ < c, d > $)作为分类依据, 迭代挖掘直至样本集中的样本为空或小于设定阈值。后续规则形式可解释为$\left\{ {\neg < a, b > \wedge < c, d > \to {C_m}} \right\}, \{ \neg < a, b > \wedge \neg < c, d > \wedge < e, f > \to {C_m}\} , \cdots $。如图 1所示, 在预测样本类别时, 以挖掘出的顺序作为优先级依次使用多个分类规则, 构成层级规则树分类器CL_HRT。

上述过程是HRT分类模式在二分类问题中的应用, 在解决多分类问题时, 采用如下方式:

给定一个多分类标签$ C=\left\{C_{0}, C_{1}, \cdots, C_{m-1}\right\}$, 将多分类问题分解为m个二分类问题, 即$\left\{ {{C_0}, \sim {C_0}} \right\}, \left\{ {{C_1}, \sim {C_1}} \right\}, \cdots , \left\{ {{C_{m - 1}}, \sim {C_{m - 1}}} \right\} $, 由此得到m个HRT分类器。若m个分类器预测结果无冲突, 则输出此唯一类别; 若存在分类器冲突, 则选择缺省类(样本中的多数类)作为最终结果。

1.3 k-HRT算法

集成学习是一种通用的机器学习策略, 与具体的算法模型无关, 其利用多个弱分类器进行集成以提升分类器的最终效果^[12]。常用的集成学习方法有bagging、boosting, 以及由两者衍生出的随机森林、Adaboost、GBDT等。为提升分类器在跨平台基因表达数据集上的泛化能力, 本文借鉴bagging的思想, 以HRT作为基分类器进行集成学习, 得到k-HRT算法。从训练集中随机去除小部分样本, 对多个基分类器CL_i进行k次重复训练, 构成分类器h_k－HRT, 对k个基分类器的预测结果进行无权投票, 将得分最高的类别作为样本的最终预测类别并输出。设x_new是等待预测的样本, 待预测样本类别为$ C = \left\{ {{C_0}, {C_1}, \cdots , {C_{m - 1}}} \right\}$, 则分类器h_k－HRT定义如下:

$ {h_{{\rm{k}} - {\rm{HRT}}}}\left( {{x_{{\rm{new}}}}} \right) = \mathop {\arg \max }\limits_C \sum\limits_{i = 1}^k I \left( {C{L_i}\left( {{x_{{\rm{new}}}}} \right),C} \right) $

(3)

I(CL_i(x_new),C)函数定义如下:

$ I\left( {C{L_i}\left( {{x_{{\rm{new}}}}} \right),C} \right) = \left\{ {\begin{array}{*{20}{l}} {1,C{L_i}\left( {{x_{{\rm{new}}}}} \right) = C}\\ {0,其他} \end{array}} \right. $

(4)

2 快速挖掘算法

在大规模基因表达数据集中, 快速挖掘算法基于相对偏移表(Relative Shifting Table, RST)实现规则预筛选, 以避免对规则空间进行暴力搜索。用N×P的矩阵R表示数据集, 基因表达数据样本数量为N, 每个样本记录的基因(即特征)值的数量为P。算法首先对矩阵R进行数据转换, 然后构建RST进行规则预筛选, 从候选集中选择得分最高的分类规则, 删去该规则所命中的样本, 最后迭代构建层次规则树。完整的层次规则树挖掘算法BuildHRT描述如下。

算法1   BuildHRT算法

输入  数据集R, 最小支持度θ, 候选组数量m

输出  分类器CL

1.DataConvert(R);

2.while |R| ≥N×θ do

3.RST=RelativeShiftTable(R);

4.Can=GetCandidates(T, m);

5.γ=GetTop(RST, Can, θ);

6.if γ!=NULL then

7.insert γ into CL;

8.for each d∈R do

9.if d satisfies γ then

10.delete d from R;

11.else

12.break

13.Set majority class of remain R as default class of CL.

BuildHRT算法的时间复杂度和其在数据集R上的迭代效率有关。在最佳情况下, 迭代一次算法终止, 即HRT算法为一层, 算法时间复杂度为O(mN)。在最差情况下, 迭代$\frac{1}{\theta} $次算法终止, 算法时间复杂度为$O\left(\frac{1+\theta}{2 \theta} m N\right) $, 此时最小支持度θ决定了时间复杂度, θ的理论取值范围为$ \left[\frac{1}{N}, 1\right], \theta$值越小, 算法复杂度越高, 取最小值$ \frac{1}{N}$时算法复杂度为$ O\left(m N^{2}\right)$。

2.1 数据转换

对于以基因对表达反转为基础的分类模式, 基因表达数据上有意义的是各基因之间表达值的大小关系, 而非单个基因的表达值。因此, 在BuildHRT算法的第1行, DataConvert(R)把记录基因表达值的矩阵转换为基因序列矩阵, 以此作为HRT分类模式挖掘算法的基础。具体地, 将矩阵R中每一个样本的基因表达值和基因名组合成一个二元组, 对每个样本内的P个二元组, 以基因表达值为键进行排序, 删去二元组中的基因表达值, 从而将样本转换为由基因名构成的有序序列, 该过程如图 2所示。

2.2 相对偏移表

解决分类问题的关键是发现在不同类样本中表现差异最大的基因表达模式, 并对其进行高效率的挖掘。为避免在高维度、大样本的跨平台基因表达数据集中进行暴力搜索, 减小规则搜索空间, 本文在数据转换的基础上设计一种相对偏移表RST, 以进行规则预筛选。具体过程如下:

1) 计算正、负类样本特征。根据式(5)计算样本的P个基因在矩阵R正类样本C₀中的总得分值$Score_g^{{C_0}} $, 根据总得分值对P个基因由高到低排序, 将排序后得到的基因序列定义为正类样本特征EP。同理可以获取负类样本特征EN。

$ Score_g^{{C_0}} = \sum\limits_{d \in {C_0}} I (d,g) $

(5)

2) 根据式(6)计算每个基因g在正类样本特征EP上相对于负类样本特征EN中的相对偏移量Dis_g。

$ Di{s_g} = I\left( {EP,g} \right) - I\left( {EN,g} \right) $

(6)

3) 对P个基因按照偏移量得分值进行排序, 得到的基因序列即为相对偏移表RST。假设在数据矩阵R上得到正类样本特征EP=daceb、负类样本特征EN=dbeca, 则计算的相对偏移量Dis_g如表 2所示, 最终求得相对偏移表RST=acdeb。

2.3 层级规则树挖掘

利用相对偏移表RST进行规则预筛选, 目标是找到在不同类中相对表达反转的基因对(a, b)。令x为负类样本特征EN相对于正类样本特征EP向前偏移最多的基因, y为EN相对于EP向后偏移最多的基因, 则基因对(x, y)有极大概率在不同类样本上完成相对表达反转。向前、向后偏移最多的基因分别位于相对偏移表T的首、尾位置。因此, 为了提高算法的稳定性, 本文从相对偏移表首、尾各取出m个基因组合为候选基因对Can(算法1中的第4行), 并从中选取置信度最高的基因对作为目标分类规则(算法1中的第5行)。

在每一层寻找分类规则γ的具体过程如下:

1) 从相对偏移表RST的首、尾各随机取出m个基因, 按照位置对应组合为m个候选基因对。

2) 每一个候选基因对(a, b)可以产生4个分类规则: $ \left\{ { < a, b > \to {C_0}} \right\}, \left\{ { < a, b > \to {C_1}} \right\}, \left\{ {\neg <a, b > \to {C_0}} \right\}, \left\{ {\neg < a, b > \to {C_1}} \right\}$。从中选取置信度得分最高的作为基因对(a, b)的唯一分类规则, 共生成m个候选分类规则。

3) 从m个候选分类规则中选取置信度得分最高且满足最小支持度θ的规则, 作为目标分类规则γ。

将分类规则γ加入层次规则树分类器中, 从训练数据集上删去被分类规则γ前项所命中的样本, 重复算法1的第3行~第12行, 直至剩余样本数量小于最小支持度θ, 将剩余样本中的多数类作为层次规则树分类器中的缺省类。

3 实验结果与分析

为验证本文分类模式的性能, 从NCBI上获取多个平台的基因数据集作为原始实验数据, 进行仿真对比与分析, 多个数据集所共有的标签是性别与年龄。

3.1 数据预处理

原始基因表达数据样本的维度为20 660, 在这些属性上存在很多缺失值。由于这些样本来自不同的GPL平台, 各样本上缺失值的数量与位置不同, 因此需要找到缺失值最少的样本与维度。本文受数据库搜索方法ripple join^[13]的启发, 使用类似的方法进行缺失值数据的消除。首先将样本及其维度根据NA值的数量从高到低分别进行排序, 然后记录维度映射表。初始化一个空矩阵, 迭代扩展行和列, 找到最大的无缺失样本矩阵。经过数据清洗后, 最终可以得到无缺失值的2个样本集, 数据分布如表 3所示。其中, 性别数据集来自7个平台:GPL10558, GPL6102, GPL6884, GPL4133, GPL6947, GPL6480, GPL570;年龄数据集来自6个平台:GPL10558, GPL6102, GPL6884, GPL6947, GPL6480, GPL570。

3.2 实验设置

性别数据集的2个标签为Male、Female, 在数据预处理阶段, 去除两性分别独有的基因特征, 以男性、女性所共有的基因作为样本特征。年龄数据集的正、负类标签为Young、Old, 由于年龄是连续的数值, 年龄数据集对Old标签的定义是年龄大于60岁的人群。

实验中的对比算法为本文k-HRT算法、k-TSP算法^[7]、SVM-RFE算法^[3]。实验分为2组:

实验1   2个数据集上的无平台信息实验。

实验2   2个数据集上的平台迁移实验。

本文自行从NCBI获取并整合实验数据, 无训练集与测试集的分别。对于第1组无平台信息实验, 训练集与测试集按8:2的比例进行无偏划分; 对于第2组平台迁移实验, 将源平台作为训练集、目标平台作为测试集。对于BuildHRT算法中的2个自定义参数, 最小支持度θ设定为20/N, 候选基因对数量m为50。用于集成的层级规则树数目k设定为3, 对比算法k-TSP中的k值也为3。

3.3 结果分析

在无平台信息实验中, 图 3、图 4分别为性别、年龄2个数据集上各算法的准确度随样本量的变化情况。由图 3、图 4可以看出, 在类均衡的性别数据集上, k-HRT算法与k-TSP算法均优于无法消除平台差异性的SVM-RFE算法; 在类不均衡的年龄数据集上, k-HRT算法的稳定性与准确性相比其他算法均具有较大优势。图 5、图 6分别为性别、年龄2个数据集上各算法的CPU运行时间随样本量的变化情况。由图 5、图 6可以看出, 相比其他2种算法, k-TSP算法效率较差, 难以应对大规模数据集。表 4所示为2个数据集上最大样本量时各算法的准确度对比。

对于平台迁移实验, 表 5所示为性别、年龄2个数据集上各算法的准确度对比, 其中, 每个数据集上准确度最高的值用加粗标出。由表 5可以看出, 该实验的结论与无平台信息实验基本一致。

3.4 生物学解释

生物的衰老过程和某些基因表达水平的变化有关^[14]。本文对年龄组实验中k-HRT分类器上的如下分类规则进行生物学分析:

If IFNA17 < ALKBH1, then Young.

If BCORL1 < COX17, then Young.

If BCORL1≥COX17 and MRAS < ZNF75D, then Old.

ALKBH1与DNA烷基化损伤修复机制有关^[15]。文献[16]研究表明, ALKBH1在胚胎干细胞转录网络中扮演着核心角色。IFNA17由巨噬细胞产生, 具有抗病毒活性, 文献[17]发现其和细胞凋亡相关。文献[18]在CEPH Uath家族数据集上的研究结果也表明, IFNA17的表达水平和年龄有关。COX17则与生物的新陈代谢功能相关^[19]。BCORL1编码的蛋白质是转录辅阻遏物, 目前尚未发现BCORL1在生物衰老方面的相关特性。

4 结束语

传统分类模式在跨平台基因表达数据平台间难以有效迁移。为此, 本文提出一种基于层级规则树的k-HRT算法。针对由跨平台特性所带来的大规模数据问题, 设计与分类模式相切合的数据转换和规则预筛选策略。同步进行规则挖掘与分类器构建, 以解决规则冗余与样本覆盖等问题。在真实基因表达数据集上的实验结果验证了k-HRT算法的可行性与高效性。本文所提出的数据转换与规则预筛选策略能够为相关的数据集挖掘工作提供一种新思路。下一步考虑将k-HRT算法应用到连续数值型跨平台数据集中, 以进行数据分类。