基于遍历基因组合的特征基因选取方法

李 杰

（大理学院数学与计算机学院，云南大理 671003）

随着DNA基因芯片技术的提高，人们研究癌症的方法也得到了快速发展。DNA 微阵列技术使得人们可以在基因水平上了解癌症发生、癌变的机理。通常生物学上认为癌症是由于一个或者某几个基因变异导致的。科学家们正致力于寻找产生癌变的基因或基因组合，近年来各种有效方法如雨后春笋般的涌现出来。

自1999年Golub〔1〕发表了关于急性白血病亚型基因芯片起，从基因芯片的角度研究各类疾病的人越来越多。由于基因芯片数据大部分是高维甚至超高维数据，而样本量却远远少于基因的维数。在统计学上，这样的数据会导致“Curse of Dimensionality”和“Over-Fitting”问题，统计学上解决这一问题的方法通常有两种：数据降维和变量选择。数据降维方法是把高维矩阵在不损失信息或损失少量信息的情况下，把高维矩阵投影到一个低维矩阵的方法，进而可以利用经典的统计模型进行分析，如 PCA〔2〕、SIR〔3〕、MAVE〔4〕等方法；变量选择方法是通过对回归系数的限制，在回归过程中不断把系数变为0 的变量剔除掉，直到达到所需要的变量个数为止，常用的方法有LASSO〔5〕、LAR〔6〕、ALASSO〔7〕、SCAD〔8〕、GLASSO〔9〕等。分类器可以用Logistic 回归、支持向量机、决策树C4.5和随机森林的方法。

从可查的文献显示，选择基因芯片中可能引发癌症的基因，主要思路〔10〕如下：首先剔除无用和冗余基因，其次对剩下的基因进行数据降维或者变量选择，然后把选择出的基因放到分类器中进行分类，最后根据分类器的预测精度进行修正和评价。但是第二步进行变量选择和数据降维时，并没有考虑到基因之间的关系。很可能会出现这样的情况：两个基因单独分析，对癌症并没有影响，但是合在一起对癌症的产生却有很大影响。就如鸡蛋和糖精只能分开吃，不能一起吃的道理一样。本文尝试从这个角度出发，考虑存在此关系的基因对癌症的分类影响。

这方面的文献不是很多。原因大概有3 点：①计算量大，基因个数通常以万计，不同的组合会呈几何级数增长；②计算受到了计算机内存和计算速度的限制；③这种方法比较笨，效率低。但随着计算机多核、多线程、分布式计算和云计算的快速发展，这类的问题似乎可以尝试运行一下。

文章主要对结肠癌数据从考虑成对的基因的相互作用出发，用Logistic分类器对正常人和结肠癌患者进行区分。结肠癌数据共包含62个样本，其中22人为正常人的基因水平，40人为结肠癌的基因水平，共2 000 个基因。考虑成对基因，则有199.9 万个不同的基因组合，遍历所有可能的情况，根据预测精度和AIC准则选择最优的基因组合并与已有的文献进行比较研究，最后根据前1 000个结果，画出基因关系图谱。

1 数据结构及预处理

1.1 数据结构结肠癌数据共包含62 个案例，其中22 个正常人，40 个癌症患者；对每个患者运用DNA 微阵列测得了2 000 个基因的表达水平。用Yi=1 表示第i个人是癌症患者，Yi=0 表示第i个人是正常人，i=1,…,62，则Y是一个二元变量，每一个基因看成一个变量，分别用G1,G2,…,G2000表示。数学模型可以抽象为表1。

表1 结肠癌数据结构

其中Xij，i=1,…,62,j=1,…,2 000 是第i个人第j个基因表达的对数值，该数据可以在BRB-Array-Tools主页上下载。

1.2 预处理由于数据取对数后各个基因的表达水平还是有很大差距，因此需要按基因对数据进行标准化，因此表1中的列需要做标准化变化，以便更好比较。一般在用分类器进行分类之前，需要对冗余基因进行处理，值得注意的是这2 000 个基因的表达水平都是非零，为了能准确反映结果，没有再剔除冗余基因。从2 000 个基因中任意选取2 个基因，不同的组合共有1 999 000 种，任意选择3 个基因，不同的组合共有种。由于选择3 个基因不同的情况达到了13.3 亿，本文暂时不予考虑，而只考虑2个基因共同作用的情形。

2 模型方法

由于Y是一个分类的二元变量，在统计学中，二元离散变量的经典模型是Logistic 模型，假设Y是因变量，x1,x2是自变量，Y取值为1 的概率为p，则Logistic模型为：

其中ε是随机误差项，β0，β1，β2是未知参数，可以通过迭代极大似然估计得到。由于只考虑2个基因联合作用的情形，自变量的个数定位2 个。把上式进行变换可以得到Y取1的概率为：

当β0，β1，β2估计出来后，根据上式就可以计算第i个样本取1的概率pi，当pi＞0.5 时，把第i个样本判定为癌症患者，pi≤0.5 时判为正常人。根据预测值和实际值的比较，可以得到表2。

表2 预测值和实际值

其中n11表示实际为癌症患者并且预测成癌症患者的人数，n10表示实际为正常人但预测成癌症患者的人数，n01表示实际为癌症患者但预测成正常的人数，n00表示实际为正常人并且预测成正常的人数，n11+n10+n01+n00=n=62，预测精度可以表示为：

预测精度越大，效果越好。

3 评价标准

模型评价标准主要有两部分组成：预测精度和统计指标。第一部分主要看预测精度，预测精度越高越好，但是预测精度单个指标并不能反映真实情况，甚至可能会出现错误，例如有100个人，90个正常人，10 个病人，但是预测的结果把这100 个人全部预测成正常人，预测精度也达到了90%，但是它把病人全部预测成正常人，没有一个预测正确，这种方法就没有区分度，需要结合统计指标。统计指标本文采用了AIC准则和β0，β1，β2的显著性水平作为评价标准，在重要性排序上文献中没有统一说法，根据需要拟定重要性大小按照预测精度，AIC信息和β0，β1，β2的显著性水平依次递减。

4 计算流程

计算流程共分为以下几步：①标准化各个基因；②产生2个基因组合的所有可能结果；③遍历所有可能的组合，在每一次遍历过程中计算预测精度、AIC 和未知参数的检验P值，并保存在txt 文档中；④读取最后的数据，按照重要性标准进行排序，找出最合理的模型。流程图见图1。

图1 计算流程图

5 实验结果

本文采用的实验环境配置：Intel（R）Core（TM）i5-3470 CPU@3.20Hz，4GB 内存，Windows 8.1 操作系统，Rstudio 平台，用R 语言进行编程。实验过程中对每个计算过程中的中间结果进行保存，共得到了1 999 000行17列的数据文件，约225.3 MB。用R进行分析得到的主要结果如下。

首先按照预测精确度对所有的数据进行排序，在此基础上再对AIC 从小到大进行排序，同时检验未知参数的t检验P值是否小于0.05，结果显示编号为X55187 和D14812 的基因组合得到最满意的结果，预测精度为0.935 5，预测表格见表3。

表3 实验结果

AIC 为 29.671，在所有 199.9 万中 AIC 中排在第二位，最小的AIC为26.198。实际的模拟结果见表4。

表4 Logistic回归结果

从表4中最后一列可以看出，包括截距项在内，未知参数在置信水平0.01下都通过了检验，参数值都是非零，具有统计学意义。根据表4，回归方程可以写为：

其中x1表示编号为X55187 的基因，x2表示编号为D14812 的基因。零假设下的离差为80.648 4，残差为 23.671 1，广义R2为：（80.6484-23.6711）/80.6484=0.7065。为了进一步说明该回归结果的稳定性，采用交叉留一检验，即依次删除掉62个样本的1个样本，剩下的61 样本作为测试集，统计测试集的预测精度。见图2。

图2 交叉留一检验预测精度趋势图

图2显示，预测精度大部分是0.934，其中有4次试验的预测精度达到了0.951，这说明选择出的基因组合具有较好的稳健性。不会因为样本的变化，而对预测精度产生较大波动。综上所述，编号为X55187和D14812的基因组合是比较适合区分结肠癌患者和正常人的。

因为遍历了所有可能情况，得到的结果丰富，为了进一步展示其它结果，我们截取了199.9万中预测精度大于0.9的所有基因组合，共有640对基因组合，对这些基因组合进行频数分析，得到结果见表5。

表5 基因频率分析结果

张靖等〔11〕研究结果显示，基因R87126，H08393，M63391，X12671是引起结肠癌的可疑基因，跟本文中频率分析结果相同的有两个基因R87126，M63391，注意这里考察的是成对基因的频率分析，也就是说本文只是统计了各个基因出现的频率，并没有考虑成对基因之间内部的关系。虽然R87126，M63391 在统计结果中出现的频率最高，但是这两者组合成的基因对并没有出现在640 对当中，重新考察R87126，M63391组成的基因对，结果显示预测精度仅为87.10%，AIC为43.234，广义R2为0.538 3，并且M63391基因的系数在0.05的置信水平下没有通过检验。因此并不是基因出现频率越大，分类的预测精度就越高。

6 结论

本文在癌症是由某个或者某几个基因共同变异导致的假设下，采用遍历所有基因组的方法，运用Logistic 分类器，对结肠癌的基因进行了筛选，按照预测精度和AIC 准则得到了最优的基因组合（X55187，D14812），并与已有文献结果进行比较，得出频率高的基因组合预测精度并不一定高的结论。由于实验条件限制，本文中并没有考察3 个以及3个以上基因组合对分类的影响。

〔1〕Alon U，Barkai N，Notterman D A，et al. Broad patterns of gene expression revealed by clustering analysis of tumor and normal colon tissues probed by oligonucleotide arrays〔J〕.Proc Natl Acad Sci，1999，96：6745-6750.

〔2〕Pearson K. On lines and planes of closest fit to systems of points in space〔J〕. Philosophical Magazine，1901，2（6）：559-572.

〔3〕Li K C. Sliced inverse regression for dimension reduction〔J〕. Journal of the American Statistical Association，1991，86：316-327.

〔4〕Xia Yingcun，Tong Howell，Li W K.An adaptive estimation of dimension reduction space〔J〕. Journal of the Royal Statistical Society：Series B（Statistical Methodology），2002，64（3）：363-410.

〔5〕Tibshirani R. Regression shrinkage and selection via the lasso〔J〕.J Royal Statist Soc B，1996，58（1）：267-288.

〔6〕Zhou Hui. The Adaptive Lasso and Its Oracle Properties〔J〕. Journal of the American Statistical Association，2006，101（476）：1418-1429.

〔7〕Fan J，Li R. Variable selection via nonconcave penalized likelihood and its oracle properties〔J〕.J Amer Statist Assoc，2001，96：1348-1360.

〔8〕Noah S，Jerome F，Trevor H，et al.A sparse-group lass〔J〕.Journal of Computational and Graphical Statistics，2013，22（2）：231-245.

〔9〕Efron B，Hastie T J，Johnstone I，et al. Least Angle Regression〔J〕.Annals of Statistics，2004，32（2）：407-499.

〔10〕李杰，邓丽君，唐胜男.基于BB—SIR方法的结肠癌特征基因提取〔J〕.世界科技研究与发展，2011，33（4）：588-591.

〔11〕张靖，胡学钢，李培培，等.基于迭代Lasso 的肿瘤分类信息基因选择方法研究〔J〕. 模式识别与人工智能，2014，27（1）：49-59.