基于线性残差独立性测试的致病基因检测

Causal Gene Identification Based on Linear Residuals Independence Test

摘要

使用因果推断相关的机器学习方法辅助检测致病基因时，作为因果推断的核心工具，条件独立性（CI）测试算法在高维生物数据场景中往往存在时间复杂度高以及准确性低等问题。为此，提出一种融合偏相关测试与线性残差独立性测试算法，压缩CI测试条件集的搜索空间，同时提高准确率。设计一种因果推断策略，在减少冗余CI测试的同时结合V结构与因果函数模型的优点，在应用于真实癌症数据的致病基因检测场景中可以区分Markov等价类，找到真正的因果关系。实验结果表明，提出的算法有较好的致病基因检测性能。

Abstract

In what way to efficiently and accurately look for the causal genes related to the corresponding diseases becomes a hot spot in searching. When no control experiment can be applied, we usually use causal discovery method to detect causal genes. However, the traditional independence tests in high-dimensional data have high time complexity and low accuracy. To alleviate this problem, we propose a residual independence test algorithm that combines partial correlation test and linear residuals independence test to compress the search space of the conditional set of CI(conditional independence) test, and improve the accuracy. Next, we design a causal discovery algorithm based on residual independence test, which can distinguish Markov equivalence classes by V-structure and causal functional model, we apply it to real cancer datasets in the detection of pathogenic genes. The results show that proposed algorithm is significantly better than existing algorithms in many aspects.

Graphical abstract

关键词

因果网络 / 因果推断 / 条件独立性 / 致病基因 / 条件独立性测试

Key words

causal network / causal inference / conditional independence / causal gene / CI（conditional independence） test

0　 引 言

随着生物技术的发展，人类基因组测序和基因组活动监测成本快速下降^[1]，微阵列已经成为检测基因活性的常用工具,若将其应用于疾病诊断和基因治疗，很有可能在未来产生巨大的商业价值^[2]，因而许多生物信息学和生物医学研究人员正着手进行这方面的应用研究。其中，识别与某些病症相关的致病基因是最具有前景的应用之一^[3]，因为它可估计患者潜在患某些病症的风险，并为他们提供进一步的治疗建议。致病基因检测在生物医学领域非常具有挑战性^[4]，是因为在许多情况下都难以进行对照实验，并且仅依靠传统的生物医学发现致病基因代价很高。为了解决这个问题，生物信息学研究者通常借助机器学习和数据挖掘的技术来观察基因表达数据^[5~7]，并设计各种方法以发现与某种特定疾病相关的重要基因。

通常而言，直接发现致病基因是一件十分困难的事情^[8,9]，研究者们经常使用特征选择方法检测与疾病相关的基因。特征选择主要有三种主流的方法：过滤、包装和嵌入。然而，这三种类型的特征选择方法在致病基因检测上都存在不同程度的问题。理论上，这些基于特征选择的方法返回的候选基因集合一般都是由以病症变量为中心的马尔可夫毯^[10]（Markov Blanket，MB）构成的，它是一个包含了致病基因的超集。因此，在假设特征选择结果100%准确的情况下，依然存在一个重要问题——特征集冗余^[11]。

由于特征选择方法存在的缺陷，基于因果推断的致病基因检测方法应运而生^[1,12]。因果网络^[13]是分析变量间因果关系的一种常用工具。因果网络通常是关于一组变量的概率图模型，节点表示变量，有向边表示变量之间蕴含的因果关系。实际上，关于边的信息往往是未知，需要进行因果推断确定节点间是否存在边，以及边的方向。在基于约束的因果推断方法^[14,15]中，通常利用条件独立性（conditional independence，CI）测试^[16,17]学习边的信息，即需要检验变量之间的CI。

通常来说，通过CI测试，因果推断算法如PC算法^[18]可以返回一个与数据实际蕴含的因果网络互为马尔可夫等价类^[12]（Markov equivalence class，MEC）的部分有向无环图^[19]（partial direct acyclic graph，PDAG），缓解了特征集冗余的问题。但通过检测数据中存在的V结构以及通过一致性传播推断因果方向往往无法在一组MEC中找到完整的因果关系。在多种情况下，基于CI测试的方法，包括PC算法和发现MB的一系列算法^[10,20]都不能确定所有被检测出的基因是否为真的致病基因。因此，现有的因果推断方法在检测致病基因过程中耗时长^[21]，查全率低^[22]，准确率低^[23]。

近些年来，研究者们又提出了基于残差的CI检测方法，这些方法能在一组MEC中找到正确的因果结构，其底层模型一般为加性噪声模型(additive noise model，ANM)^[24~29]。ANM可以表示为：y=f\left(X\right)+n，其中n是一个独立于X的加性噪声，X是y的父亲变量集。Grosse-Wentrup等^[27]提出x-f\left(Z\right)\perp \left(y,Z\right)是x\perp y|Z的一个充分非必要条件。Gao等^[15]提出在忠诚性条件和马尔可夫条件下，x-E\left(x|Z\right)\perp y-E\left(y|Z\right)可以推断出x\perp y|Z，不过这个结论仅在线性高斯和非高斯场景下成立，且其中独立性计算往往需要进行高复杂度的基于核方法的密度估算，难以用到包含成千上万个特征变量的基因表达数据。因此，想要通过这些方法来推断因果关系仍然困难重重。为了解决此问题，本文假设这些因果关系是线性的，简化了模型，提出一个全新的方法称为CGD（causal gene discovery），通过基于线性残差独立性测试实现高效的致病基因检测。

1　 基础知识

定义1（因果网络） 因果网络是表示变量间概率依赖关系的有向无环图（direct acyclic graph，DAG），它可表示为一个三元组G=\left(X,E,P\right)。其中，X=\left\{x_{\mathrm{1}},x_{\mathrm{2}},\dots ,x_n\right\}表示DAG中的所有节点的集合；E=e\left(x_i,x_j\right)|x_i,x_j\in X表示DAG中每两个节点间的单向边的集合，e\left(x_i,x_j\right)表示x_i,x_j间存在依赖关系x_i\to x_j，P=\{P(x_i|p{a}_{x_i}),p{a}_{x_i}\in X\}是条件概率的集合，P(x_i|p{a}_{x_i})表示x_i的父节点集p{a}_{x_i}对x_i的概率性影响。

因果网络本质上是联合概率分布P(x_{\mathrm{1}},x_{\mathrm{2}},\cdots ,x_n)中所有CI的一种图形化表示。

定义 2（d-分离） 给定一个DAG G\left(V,E\right)，令P为v_i与v_j之间的一条（先不考虑方向）路径，则称关于v_i与v_j的一条路径P被节点集Z阻塞了，只要满足下面的任意一个条件：

1) P包含一个顺连结构，如v_i\cdots \leftarrow v_k\leftarrow \cdots v_j 或 v_i\cdots \to v_k\to \cdots v_j，且v_k\in Z；

2) P包含一个分连结构，如v_i\cdots \leftarrow v_k\to \cdots v_j，且v_k\in Z；

3) P包含一个V结构，如v_i\cdots \to v_k\leftarrow \cdots v_j，且v_k\notin Z，同时Z也不包含任何v_k的后代节点。

若v_i与v_j之间的所有路径都被Z阻塞，则v_i与v_j被节点集Z给d-分离，Z也被称为v_i与v_j的d-分离集。如果某个Z是所有d-分离集中包含节点数最少的那个，则将称其为最小d-分离集^[8]（minimal d-separator）。

定义3（条件独立性测试） 令X，Y和Z表示三个随机变量的集合，当给定条件Z时，对X和Y进行CI检验，如果X和Y关于Z条件独立，则表示为X\perp Y|Z。

定义3表明，如果X和Y之间是关于Z独立的，即使知道Z的值也无法通过X（或Y）的值得到任何对应Y（或 X）的额外信息。

定义 4（因果Markov假设） 给定一个因果网络G，令P\left(X,Y,Z\right)为三个变量集X，Y和Z的联合概率分布，若X和Y被Z给d-分离 ⇒ X\perp Y|Z，则称P\left(X,Y,Z\right)满足关于G的因果马尔可夫（Markov）假设。

定义 5（因果忠实性假设） 给定一个因果网络G，令P\left(X,Y,Z\right)为三个变量集X，Y和Z的联合概率分布，如果X\perp Y|Z ⇒ X和Y被Z给d-分离，则称P\left(X,Y,Z\right)满足关于G的因果忠实性假设或因果忠实性条件成立。

可以看出，因果Markov假设和忠实性假设将观测数据的组合概率分布与相应的因果网络结构联系起来了。特别地，如果因果忠诚假设得到满足，则可以使用CI测试找到一组与真实因果关系Markov等价的PDAG。一般地，因果推断的时候都会先验地假设因果忠实性条件成立，否则难以进行有效的推断。本文工作也是在这一假设条件成立下进行研究的。

2　 线性残差独立性测试

如前面所说，若对任意两个变量x与y直接进行CI测试，当数据维度较大的时候，运算的时间成本会比较高。在生物信息领域，数据维度通常很高，直接运行CI测试方法很可能在可接受的时间范围内都无法返回结果。本节将详细描述我们设计的高效线性残差独立性测试方法。

2.1　线性加性噪声模型

本文工作是在线性加性噪声模型（linear additive noise model，LANM）基础上进行的。LANM是一个二元组\left(S,P\left(X\right)\right),S=\left\{S_{\mathrm{1}},S_{\mathrm{2}},\dots ,S_n\right\}是n个方程的集合：

S_i:x_i=f_i\left(p{a}_{x_i}\right)+{\epsilon }_i,i=\mathrm{1,2},\dots ,n

其中，每个f_i都是一个线性函数,p{a}_{x_i}对应于DAG G中x_i的父亲节点集,噪声变量{\epsilon }_i服从独立分布,且满足{\epsilon }_i\perp p{a}_{x_i}。

如果一个LANM可以通过残差与自变量的独立关系区分不对称的因果方向，那么称它是可以识别（identifiable）的^[23]。在大部分线性非高斯或非线性高斯场景中，LANM都是可以识别的^[19]。本文假设所有的数据在生成过程中都满足LANM。我们的任务是随机挑选两个节点x和y，然后测试在给定一组变量Z下x和y是否条件独立。下面，介绍不同场景下，CI测试的原理。

2.2　线性场景中的残差独立性与CI的关系

在本节，首先给出关于CI的一个经典理论^[30]，它描述了如何在平方可积函数空间中把相关性与独立性关联起来。

条件独立性特性^[30]（characterization of conditional independence，CCI）：令X,Y,Z为三个随机变量或随机变量集，定义：

E_{\mathrm{1}}=g\in L_{XZ}^{\mathrm{2}},E\left(g|Z\right)=\mathrm{0}

E_{\mathrm{2}}=h\in L_{YZ}^{\mathrm{2}},E\left(h|Z\right)=\mathrm{0}

E_{\mathrm{3}}=g^{\mathrm{\text{'}}}\in L_X^{\mathrm{2}},E\left(g^{\mathrm{\text{'}}}\right)=\mathrm{0}

E_{\mathrm{4}}=h^{\mathrm{\text{'}}}\in L_Y^{\mathrm{2}},E\left(h^{\mathrm{\text{'}}}\right)=\mathrm{0}

其中，L_{XZ}^{\mathrm{2}}，L_{YZ}^{\mathrm{2}}，L_X^{\mathrm{2}}，L_Y^{\mathrm{2}}分别表示\left(X,Z\right),\left(Y,Z\right),X,Y的平方可积函数空间,则下面四个等式相互等价：

CCI-1：X\perp Y|Z

CCI-2：\forall g\in E_{\mathrm{1}}\mathrm{和}\forall h\in E_{\mathrm{1}},E\left(gh\right)=\mathrm{0}

CCI-3：\forall g\in E_{\mathrm{1}}\mathrm{和}\forall h^{\mathrm{\text{'}}}\in E_{\mathrm{4}},E\left(g{h}^{\mathrm{\text{'}}}\right)=\mathrm{0}

CCI-4：\forall h\in E_{\mathrm{2}}\mathrm{和}\forall g^{\mathrm{\text{'}}}\in E_{\mathrm{3}},E\left(h{g}^{\mathrm{\text{'}}}\right)=\mathrm{0}

特别地，当Z=\mathrm{\varnothing }时,可以得到：

X\perp Y\iff \forall g^{\mathrm{\text{'}}}\in E_{\mathrm{3}}\mathrm{和}\forall h^{\mathrm{\text{'}}}\in E_{\mathrm{4}},E\left(g^{\mathrm{\text{'}}}{h}^{\mathrm{\text{'}}}\right)=\mathrm{0}(1)

CCI描述的是一种在平方可积函数空间中，相关性与独立性的等价关系。下面给出文献[19]的关于条件独立性和残差的独立性的理论结果。

引理1^[19] 给定m+\mathrm{2}个随机变量：x,y以及Z=\left\{z_{\mathrm{1}},z_{\mathrm{2}},\dots ,z_m\right\},它们都是由独立的随机变量s_i\left(i=\mathrm{1,2},\mathrm{3},\dots ,l\right)组成的一个线性组合,如果所有的s_i符合联合高斯分布,当且仅当x-{\rm E}\left(x|Z\right)\perp y-{\rm E}\left(y|Z\right)时，x\perp y|Z成立。

引理2^[19] 给定m+\mathrm{2}个随机变量：x,y以及Z=\left\{z_{\mathrm{1}},z_{\mathrm{2}},\dots ,z_m\right\},它们都是由独立的随机变量s_i\left(i=\mathrm{1,2},\mathrm{3},\dots ,l\right)组成的一个线性组合,如果所有的s_i符合联合高斯分布，那么以下条件相互等价：

1) x\perp y|Z

2) x-E\left(x|Z\right)\perp y

3) y-E\left(y|Z\right)\perp x

4) x-E\left(x|Z\right)\perp y-E\left(y|Z\right)

这个结论意味着在线性高斯分布的条件下，可以将测试x,y和Z之间的条件独立简化为测量x-E\left(x|Z\right)\perp y或者y-E\left(y|Z\right)\perp x。而因为此种高斯场景下，偏相关性^[31]与独立性存在等价关系，即

{\rho }_{xy.Z}=\mathrm{0}\iff x-E\left(x|Z\right)\perp y-E\left(y|Z\right)(2)

因此，通过偏相关系数{\rho }_{xy.Z}的值就可以直接推断CI是否成立。然而，在线性场景中，更加常见的是非高斯情形，我们要作进一步的处理，将根据引理1和引理2，给出线性残差独立性测试方法。

推论1 给定m+\mathrm{2}个随机变量：x,y以及Z=\left\{z_{\mathrm{1}},z_{\mathrm{2}},\dots ,z_m\right\},它们都是由独立的随机变量s_i\left(i=\mathrm{1,2},\mathrm{3},\dots ,l\right)组成的一个线性组合,如果所有的s_i均服从非高斯分布，那么

引理3^[19] 给定m+\mathrm{2}个随机变量：x,y以及Z=\left\{z_{\mathrm{1}},z_{\mathrm{2}},\dots ,z_m\right\},它们都是由独立的随机变量s_i\left(i=\mathrm{1,2},\mathrm{3},\dots ,l\right)组成的一个线性组合,如果所有的s_i均服从非高斯分布，那么

x-E(x|Z)\perp y-E(y|Z)\iff x\perp y|Z(4)

推论2 给定m+\mathrm{2}个随机变量：x,y以及Z=\left\{z_{\mathrm{1}},z_{\mathrm{2}},\dots ,z_m\right\},它们都是由独立的随机变量s_i\left(i=\mathrm{1,2},\mathrm{3},\dots ,l\right)组成的一个线性组合,如果所有的s_i均服从非高斯分布，那么

基于引理3与推论2，本文提出的线性残差独立性测试（linear residuals independence test，LRIT）方法，见算法1。由于{\rho }_{xy.Z}的计算复杂度远低于x-E\left(x|Z\right)\perp y-E\left(y|Z\right)，而通常推断过程中满足x-E\left(x|Z\right)\perp y-E\left(y|Z\right)的情况要比满足x-E\left(x|Z\right)y-E\left(y|Z\right)的情况少，因此两者可以很好地融合。具体地，给定任意两个变量x和y，以遍历的方式测试条件集Z，先判断是否满足{\rho }_{xy.Z}\ne \mathrm{0}，若不满足，则对x-E\left(x|Z\right)与y-E\left(y|Z\right)是否满足独立性进行测试；若满足，则返回x\perp y|Z。这种测试方法极大地降低了独立性测试的使用次数，提高了CI测试的速度。

要注意的是，在最坏的情况下，该测试方法的算法复杂度等价于x-E\left(x|Z\right)与y-E\left(y|Z\right)的复杂度，即O(n+IT) = O(IT)，其中O(IT)为单次独立性测试的复杂度，其值取决于使用的具体算法。例如，若使用相关性测试，则O(IT) = O(n)。

2.3　致病基因检测算法

基因检测不同于一般的因果推断，因为基因检测只需要确定哪些基因直接影响疾病变量，因此在进行x\perp y|Z的CI测试过程中，可以直接将x看作为病症变量，即单目标因果推断场景。然而，即使在致病基因检测场景下，还需要在上万个基因中遍历搜索其中蕴含的因果关系，计算量依然非常巨大，现有方法在CI测试过程中，不可避免地引入了大量冗余CI测试，若冗余CI测试返回错误的结果，则必然造成真实的致病基因被错误地排除掉，导致查全率低。基于这个观测，我们设计了一种高效的致病基因检测算法CGD（算法2），在每次迭代CI测试中减少了d-分离集的规模，使得该算法能够处理包含数千个特征的基因表达数据。

算法2主要有3个步骤：基于CI测试判断边是否存在，基于V结构原理判断因果方向，以及基于残差独立原理判断因果方向。算法2需要检测两两变量之间的CI，因此复杂度为：O(n²2 ^n-2CI)。

1）　基于CI测试判断边是否存在

在传统的如PC算法中，需要对基因表达数据集V中的每个特征y测试x\perp y|Z。如果CI成立，那么可以由因果忠实性条件判断x与y不相邻，即y必然为非致病基因。一般的基于约束的方法中应用CI测试需要遍历测试V{\backslash }_y中的每个子集Z，本文CGD采取先检查V{\backslash }_{y,N}中的子集Z（算法2，第3行），如果CI不成立，再检测V{\backslash }_y（算法2，第5行）。随着进行CI测试次数增加，集合N的数量也会增长，因此CGD算法可以避免冗余CI测试，算法效率更高。尤其是在大多数生物相关的稀疏的网络场景，大部分非邻接节点会在测试V{\backslash }_{y,N}的时候就被检测出来。

然而，因为Z\subseteq V{\backslash }_{y,N},那么CI是否都能被检测到，即会不会引入冗余特征？为了解答这个问题，我们给出下面这个推论。

推论3 给定两个变量x与y,若Z\left(y\notin Z\right)是x的父子变量集中的任意子集，且满足因果忠诚性条件并且x y|Z，则x和y相邻，或者y是x的后代且两者的最小d-分离集至少包含一个非x父子变量。

推论3易由文献[17]中Theorem 3推导得到。推论3表明，算法CGD在使用CI测试判断边是否存在的过程中（算法2~行），尽管CI测试条件集被缩小为Z\subseteq V{\backslash }_{y,N}，却仅在非常特殊的情况下（y是x的后代节点且两者的最小d-分离集至少包含一个非x父子变量）检测不到CI。因此，很大程度上提高了CI测试效率，若这种特殊情况发生了，也会在检测V{\backslash }_y（算法2，第5行）的时候检测到CI，因此算法是完备的。要注意的是，算法2的第2~10行仅检测出病症变量x的父子变量集，而致病基因仅是父亲集，这是致病基因检测中无法容忍的，因此需要在因果方向学习时，剔除所有子变量。

2）　基于V结构原理判断因果方向

算法2中第11~13行需要对y_{\mathrm{1}}-x-y_{\mathrm{2}}是否构成V结构进行判断，这是非常成熟的算法，本工作遵循PC算法里面的标准流程进行判断即可。如果不存在V结构，则需要作进一步处理。

3）　基于残差独立原理判断因果方向

算法2中的第16行，通过测试x-E\left(x|Z\right)\perp y\left(y\in Z\right),确定y是否为x的直接因节点。事实上，如果x与y相邻，且x-E\left(x|Z\right)\perp y\left(y\in Z\right)成立，那么在线性非高斯场景下，y\to x通常是满足的，这可以由现有的相关理论^[19]保证。

引理 4^[19] 给定变量x以及一个通过可识别LANM生成的变量集Z，对于任意一个邻接于x的y\left(y\in Z\right)，在满足因果忠诚性条件下有以下结论：x-E\left(x|Z\right)\perp y\Rightarrow y\to x\left(\mathrm{或}y\right)。

引理4表明，x的父亲节点可以通过检查x-E\left(x|Z\right)\perp y\left(y\in Z\right)是否成立来判断。基于V结构与基于残差独立检测的因果方向学习策略均不会引入x的后代节点。因此，尽管算法CGD在使用CI测试判断边是否存在的过程中无法检测到所有的CI，即引入了x的后代节点，但在后续因果方向学习中都将这些后代节点一一排除，可见，算法CGD最终返回的是x的完整的父亲节集，即致病基因集。

3　 实 验

3.1　数据集

本文实验在基因检测领域具有代表性的白血病数据集^[32]上进行，该数据集包含7 129个基因数据，样本来自72个不同的白血病患者，包含47例急性淋巴细胞白血病（ALL）以及25例急性髓细胞白血病（AML）。我们将CGD与另外四个经典的致病基因检测方法CDHD^[12]、SVS^[1]、FRCIT^[26]和CGI^[17]进行比较。CDHD是一种通过PC算法进行骨架学习^[33]，然后利用信息几何模型IGCI^[34]进行因果方向学习的方法；SVS是一种基于V结构的致病基因检测算法；FRCIT是一种基于峰度检测的因果推断算法，它用峰度检测代替了核密度估算，在高维场景中有较好效果；CGI是一种基于非线性回归测试的致病基因识别算法，它采用基于高斯过程回归与核独立性测试的非线性回归测试检测CI。本文实验环境为：11th Gen Intel(R) Core(TM) i7-11700K 3.60 GHz, RAM 32 GB, MATLAB 2020b。

3.2　实验结果及分析

CGD方法从7 129个基因中发现了66个致病基因。我们重点关注在这66个被发现的致病基因中，有多少是被现有研究工作报道过的与急性白血病的发生、诊断以及治疗相关的。表1中列出了通过CGD方法找到的被报道过的相关基因，共有24个。表2为通过 CDHD、SVS、FRCIT和CGI方法得到的基因。

CDHD只对基于返回的Markov等价类中的候选基因进行排序，而不是直接找到确切的致病基因^[14]。该方法返回了前10个基因的排序作为致病基因，而其中5/10的基因被报道与人类急性白血病相关联。SVS只返回了4个致病基因，相对来说，它的召回率非常低。主要原因是SVS仅使用V结构探测致病基因，当非V结构中包含了致病基因时，这些基因无法被检测到。此外，在这个白血病数据集中只有7 129个基因，可能有一些未观察到的基因参与了一些V结构的构成，在缺少未观察到的基因的情况下，SVS就无法恢复相应的V结构，因此无法找到这些致病基因。FRCIT完全依赖于噪声的独立性原理检测基因，返回了16个致病基因，其中有9个被报道过。CGI的表现是这4个对比方法里面最好的，共返回了37个基因，其中有15个被报道过。

对比这些方法，从发现的基因中被报道与急性白血病相关的数量来说，CGD方法在查全率和查准率性能方面明显是最高的。CGD找到的候选基因中，有24个基因是被报道出来的，另外42个基因对于白血病的因果作用尚不明确，但这并不意味着，没有被报道出来的基因就是找错了。如果其他研究学者可以对这些基因给予更多的关注，很多重要且有价值的信息可能就会被发现。这需要运用一些其他方法来验证这个观点。

接下来，我们将从另一个方面评估所提出的方法，检查找到的基因是否可以构成分类特征的子集，即CGD发现的66个基因是否有区分ALL和 AML的能力。从特征选择和因果关系的角度来看，致病基因是分类特征集的一个重要子集，是病症变量的马尔可夫毯一个重要组成部分。因此，致病基因在分类上必然具有较好的能力。

图1展示了CGD发现的基因的表达水平值在AML和ALL的分布情况。每个基因的表达水平在样本中都已进行标准化处理，使其均值为0，标准差为1。对于每个基因而言，基因表达水平高于均值时显示为偏红色（越高越红），低于均值时显示为偏蓝色（越低越蓝）。从图1可以看到，本文发现的66 个基因的表达水平值在ALL和AML间具有明显差异性，其中没有任何基因在两类中均匀地表达。这说明，这66个基因对ALL与AML的分类效果较好。

从统计学的角度来看，致病基因还有一个重要的特征，即通常这些基因比与病症变量没有直接关系的其他基因更接近病症变量。图2给出了CGD所找到的致病基因与病症变量之间的距离关系。距离用\mathrm{1}-\left|{\rho }_{xy}\right|来度量，其中x为基因变量，y为病症变量，{\rho }_{xy}是x和y的Pearson相关系数。首先，将病症变量固定在坐标（1，1）上；然后，随机地将某个基因x放置在以（1，1）为圆心，半径为\mathrm{1}-\left|{\rho }_{xy}\right|的圆周上。

通过观察致病基因与病症变量之间的接近程度，可以发现所找到的致病基因与疾病变量之间的平均距离远小于其他基因与疾病变量之间的平均距离，如表3所示。这一结果验证了CGD方法的有效性。

4　 结 语

本文提出了一种基于线性残差独立性测试的致病基因检测算法。为了克服“维度诅咒”问题，本文的线性残差独立性测试结合了偏相关测试速度快的优点；特别地，所提出的因果方向学习策略融合了V结构识别算法与线性加噪声模型的优点，使得方向学习不受限于Markov等价类与V结构的数量。我们在白血病数据集上，从三个方面对算法进行了测试，实验结果表明提出的方法优于同类其他方法。

由于提出的算法是一种线性方法，一旦数据背后因果关系的非线性程度较高，则算法效果无法得到保障。如何在非线性场景设计高效的致病基因检测算法将是我们后续一个研究重点。

参考文献

[1] CAI R C， ZHANG Z J， HAO Z F. Causal gene identification using combinatorial V-structure search［J］. Neural Networks， 2013， 43： 63-71. DOI： 10.1016/j.neunet.2013.01.025 .

[2] MOLLA M， WADDELL M， PAGE D， et al. Using machine learning to design and interpret gene-expression microarrays［J］. AI Magazine， 2004， 25（1）： 23-44.

[3] NOBLE D. Genes and causation［J］. Philosophical Transactions of the Royal Society A： Mathematical， Physical and Engineering Sciences， 2008， 366（1878）： 3001-3015. DOI： 10.1098/rsta.2008.0086 .

[4] COOKSON W， LIANG L M， ABECASIS G， et al. Mapping complex disease traits with global gene expression［J］. Nature Reviews Genetics， 2009， 10（3）： 184-194. DOI： 10.1038/nrg2537 .

[5] DING C， PENG H. Minimum redundancy feature selection from microarray gene expression data［J］. Journal of Bioinformatics & Computational Biology， 2008，3（2）： 185-205.

[6] BLEULER S， PRELIC A， ZITZLER E. An EA framework for biclustering of gene expression data［C］//Proceedings of the 2004 Congress on Evolutionary Computation. New York： IEEE Press， 2004： 166-173. DOI： 10.1109/CEC.2004.1330853 .

[7] CHAI L E， MOHAMAD M S， DERIS S， et al. Inferring gene regulatory networks from gene expression data by a dynamic bayesian network-based model［C］//Distributed Computing and Artificial Intelligence： 9th International Conference. Heidelberg：Springer， 2012： 379-386. DOI： 10.1007/978-3-642-28765-7_45 .

[8] GUYON I， ELISSEEFF A. An introduction to variable and feature selection［J］. Journal of Machine Learning Research， 2003， 3： 1157-1182. DOI：10.1162/153244303322753616 .

[9] CAI J， LUO J W， WANG S L， et al. Feature selection in machine learning： A new perspective［J］. Neurocomputing， 2018， 300： 70-79. DOI： 10.1016/j.neucom.2017.11.077 .

[10] WANG H， LING Z， YU K， et al. Towards efficient and effective discovery of Markov blankets for feature selection［J］. Information Sciences， 2020， 509： 227-242. DOI： 10.1016/j.ins.2019.09.010 .

[11] 龚俭， 梅海彬， 丁勇， 等. 多特征关联的入侵事件冗余消除［J］. 东南大学学报（自然科学版）， 2005， 35（3）： 366-371. DOI： 10.3321/j.issn： 1001-0505.2005.03.010 .

[12] GONG J， MEI H B， DING Y， et al. Multi-feature correlation redundance elimination of intrusion event［J］. Journal of Southeast University （Natural Science Edition）， 2005， 35（3）： 366-371. DOI： 10.3321/j.issn： 1001-0505.2005.03.010（Ch ）.

[13] HAO Z F， ZHANG H， CAI R C， et al. Causal discovery on high dimensional data［J］. Applied Intelligence， 2015， 42（3）： 594-607. DOI： 10.1007/s10489-014-0607-0 .

[14] PEARL J. Causality［M］. Cambridge：Cambridge University Press， 2009. DOI: 10.1017/cbo9780511803161 .

[15] EDWARDS D. Introduction to Graphical Modelling［M］. Berlin：Springer Science & Business Media， 2012.

[16] GAO T， JI Q. Local causal discovery of direct causes and effects［J］. Advances in Neural Information Processing Systems， 2015， 2： 2512–2520. DOI： 10.5555/2969442.2969520 .

[17] TSAMARDINOS I， BROWN L E， ALIFERIS C F. The max-min hill-climbing Bayesian network structure learning algorithm［J］. Machine Language， 2006， 65（1）： 31-78. DOI： 10.1007/s10994-006-6889-7 .

[18] ZHANG H， YAN C X， XIA Y W， et al. Causal gene identification using non-linear regression-based independence tests［J］. IEEE/ACM Transactions on Computational Biology and Bioinformatics， 2023， 20（1）： 185-195. DOI： 10.1109/TCBB.2022.3149864 .

[19] SPIRTES P， GLYMOUR C N， Causation SCHEINES R.， Prediction， and Search ［M］. Cambridge：MIT Press， 2000.

[20] ZHANG H， ZHOU S G， GUAN J H， et al. Measuring conditional independence by independent residuals for causal discovery［J］. ACM Transactions on Intelligent Systems and Technology， 2019， 10（5）： 1-19. DOI： 10.1145/3325708 .

[21] PELLET J P， ELISSEEFF A. Using Markov Blankets for causal structure learning［J］. Journal of Machine Learning Research， 2008， 9（7）： 1295–1342. DOI： 10.1145/1390681.1442776 .

[22] ZHANG H， ZHOU S G， YAN C X， et al. Recursively learning causal structures using regression-based conditional independence test［J］. Proceedings of the AAAI Conference on Artificial Intelligence， 2019， 33（1）： 3108-3115. DOI： 10.1609/aaai.v33i01.33013108 .

[23] ZHANG H， ZHOU S G， GUAN J H. Measuring conditional independence by independent residuals： Theoretical results and application in causal discovery［J］. Proceedings of the AAAI Conference on Artificial Intelligence， 2018， 32（1）： 2029 -2036. DOI： 10.1609/aaai.v32i1.11555 .

[24] ZHANG H， ZHOU S G， ZHANG K， et al. Causal discovery using regression-based conditional independence tests［J］. Proceedings of the AAAI Conference on Artificial Intelligence， 2017， 31（1）： 1250-125. DOI： 10.1609/aaai.v31i1.10698 .

[25] ZHANG K， HYVÄRINEN A. On the identifiability of the post-nonlinear causal model［C］//Proceedings of the 25th Conference on Uncertainty in Artificial Intelligence. New York： ACM， 2009： 647-655. DOI： 10.5555/1795114.1795190 .

[26] PETERS J， JANZING D， SCHÖLKOPF B. Causal inference on discrete data using additive noise models［J］. IEEE Transactions on Pattern Analysis and Machine Intelligence， 2011， 33（12）： 2436-2450. DOI： 10.1109/TPAMI.2011.71 .

[27] ZHANG H， ZHANG K， ZHOU S G， et al. Testing independence between linear combinations for causal discovery［J］. Proceedings of the AAAI Conference on Artificial Intelligence， 2021， 35（7）： 6538-6546. DOI： 10.1609/aaai.v35i7.16810 .

[28] GROSSE-WENTRUP M， JANZING D， SIEGEL M， et al. Identification of causal relations in neuroimaging data with latent confounders： An instrumental variable approach［J］. NeuroImage， 2016， 125： 825-833. DOI： 10.1016/j.neuroimage.2015.10.062 .

[29] RAMSEY J D. A scalable conditional independence test for nonlinear， non-Gaussian data［EB/OL］. 2014： arXiv： 1401.5031.

[30] FLAXMAN S R， NEILL D B， SMOLA A J. Gaussian processes for independence tests with non-iid data in causal inference［J］. ACM Transactions on Intelligent Systems and Technology， 2016， 7（2）： 1-23. DOI： 10.1145/2806892 .

[31] PDAUDIN J J. Partial association measures and an application to qualitative regression［J］. Biometrika， 1980， 67（3）： 581-590. DOI： 10.1093/biomet/67.3.581 .

[32] PENG J， WANG P， ZHOU N F， et al. Partial correlation estimation by joint sparse regression models［J］. Journal of the American Statistical Association， 2009， 104（486）： 735-746. DOI： 10.1198/jasa.2009.0126 .

[33] GOLUB T R， SLONIM D K， TAMAYO P， et al. Molecular classification of cancer： Class discovery and class prediction by gene expression monitoring［J］. Science， 1999， 286（5439）： 531-537. DOI： 10.1126/science.286.5439.531 .

[34] CHEN A H， TSAU Y W， LIN C H. Novel methods to identify biologically relevant genes for leukemia and prostate cancer from gene expression profiles［J］. BMC Genomics， 2010， 11： 274. DOI： 10.1186/1471-2164-11-274 .

[35] JANZING D， STEUDEL B， SHAJARISALES N， et al. Justifying information-geometric causal inference［M］//Measures of Complexity. Cham： Springer International Publishing， 2015： 253-265. DOI： 10.1007/978-3-319-21852-6_18 .

[36] KAGER L， LION T， ATTARBASCHI A， et al. Incidence and outcome of TCF3-PBX1-positive acute lymphoblastic leukemia in Austrian children［J］. Haematologica， 2007， 92（11）： 1561-1564. DOI： 10.3324/haematol.11239 .

[37] KEIM D， HAILAT N， MELHEM R， et al. Proliferation-related expression of p19/nm23 nucleoside diphosphate kinase［J］. The Journal of Clinical Investigation， 1992， 89（3）： 919-924. DOI： 10.1172/JCI115672 .

[38] YANG P， ZHOU B B， ZHANG Z， et al. A multi-filter enhanced genetic ensemble system for gene selection and sample classification of microarray data［J］. BMC bioinformatics， 2010， 11（1）： 1-12. DOI： 10.1186/1471-2105-11-S1-S5 .

[39] GRISOLANO J L， WESSELSCHMIDT R L， PELICCI P G， et al. Altered myeloid development and acute leukemia in transgenic mice expressing PML-RARα under control of cathepsin G regulatory sequences［J］. Blood， 1997， 89（2）： 376-387. DOI： 10.1006/bcmd.1997.0116 .

[40] BERTOLI S， PAUBELLE E， BÉRARD E， et al. Ferritin heavy/light chain （FTH1/FTL） expression， serum ferritin levels， and their functional as well as prognostic roles in acute myeloid leukemia［J］. European journal of haematology， 2019， 102（2）： 131-142. DOI： 10.1111/ejh.13183 .

[41] STORR J， DOLAN G， COUSTAN-SMITH E， et al. Value of monoclonal anti-myeloperoxidase （MPO7） for diagnosing acute leukaemia［J］. Journal of clinical pathology， 1990， 43（10）： 847-849. DOI： 10.1136/jcp.43.10.847 .

[42] LINENBERGER M L. CD33-directed therapy with gemtuzumab ozogamicin in acute myeloid leukemia： Progress in understanding cytotoxicity and potential mechanisms of drug resistance［J］. Leukemia， 2005， 19（2）： 176-182. DOI： 10.1038/sj.leu.2403598 .

[43] HATFIELD K J， OLSNES A M， GJERTSEN B T， et al. Antiangiogenic therapy in acute myelogenous leukemia： Targeting of vascular endothelial growth factor and interleukin 8 as possible antileukemic strategies［J］. Current Cancer Drug Targets， 2005， 5（4）： 229-248. DOI： 10.2174/1568009054064651 .

[44] ZAGHLOUL A， KAMAL H， ISMAIL M M， et al. Soluble carcinoembryonic antigen cell adhesion molecule 1， 6 and 8 in acute myeloid leukemia： Their relation to survival and prognosis［J］. Life Science Journal， 2018， 15（4）：1-11. DOI： 10.7537/marslsj150418.01 .

[45] SELLIN M E， HOLMFELDT P， STENMARK S， et al. Op18/Stathmin counteracts the activity of overexpressed tubulin-disrupting proteins in a human leukemia cell line［J］. Experimental Cell Research， 2008， 314（6）： 1367-1377. DOI： 10.1016/j.yexcr.2007.12.018 .

[46] JANOWSKA-WIECZOREK A， BELCH A R， JACOBS A， et al. Increased circulating colony-stimulating factor-1 in patients with preleukemia， leukemia， and lymphoid malignancies［J］. Blood， 1991， 77（8）： 1796-1803. DOI: 10.1182/blood.v77.8.1796.bloodjournal7781796 .

[47] CHENG C L， HOU H A， LEE M C， et al. Higher bone marrow LGALS3 expression is an independent unfavorable prognostic factor for overall survival in patients with acute myeloid leukemia［J］. Blood， 2013， 121（16）： 3172-3180. DOI： 10.1182/blood-2012-07-443762 .

[48] NAGAI T， HARIGAE H， ISHIHARA H， et al. Transcription factor GATA-2 is expressed in erythroid， early myeloid， and CD34+ human leukemia-derived cell lines［J］. Blood， 1994， 84（4）：1074-84. DOI： 10.1016/0955-3886（94）90139-2 .

[49] HOYOS V， SAVOLDO B， QUINTARELLI C， et al. Engineering CD19-specific T lymphocytes with interleukin-15 and a suicide gene to enhance their anti-lymphoma/leukemia effects and safety［J］. Leukemia， 2010， 24（6）： 1160-1170. DOI： 10.1038/leu.2010.60 .

[50] RUSINIAK M E， YU M， ROSS D T， et al. Identification of B94 （TNFAIP2） as a potential retinoic acid target gene in acute promyelocytic leukemia［J］. Cancer Research， 2000， 60（7）： 1824-1829. DOI： 10.1097/00002820-200004000-00012 .

[51] RUTELLA S， RUMI C， PUGGIONI P， et al. Expression of thrombospondin receptor （CD36） in B-cell chronic lymphocytic leukemia as an indicator of tumor cell dissemination［J］. Haematologica， 1999， 84（5）： 419-424. DOI： 10.1016/S0268-9499（99）90103-0 .

[52] SEIKI M， INOUE J， TAKEDA T， et al. Direct evidence that p40x of human T-cell leukemia virus type I is a trans-acting transcriptional activator［J］. The EMBO Journal， 1986， 5（3）： 561-565. DOI： 10.1002/j.1460-2075.1986.tb04247.x .

[53] TANAKA H， ITO T， KYO T， et al. Treatment with IFN？ in vivo up-regulates serum-soluble TNF-related apoptosis inducing ligand （sTRAIL） levels and TRAIL mRNA expressions in neutrophils in chronic myelogenous leukemia patients［J］. European Journal of Haematology， 2007， 78（5）： 389-398. DOI： 10.1111/j.1600-0609.2007.00834.x .

[54] HU X L， YANG D F， ZIMMERMAN M， et al. IRF8 regulates acid ceramidase expression to mediate apoptosis and suppresses myelogeneous leukemia［J］. Cancer Research， 2011， 71（8）： 2882-2891. DOI： 10.1158/0008-5472.CAN-10-2493 .

[55] KIDSON C. Two leukocyte enzymes with catalase activity in acute leukemia［J］. Biochemical and Biophysical Research Communications， 1962， 9（1/2）： 138-141. DOI： 10.1016/0006-291X（62）90102-X .

[56] 陈玲. 急性白血病患者骨髓中细胞间粘附分子-3的表达与临床关系的研究［D］. 济南： 山东大学， 2008.

[57] CHEN L. Study on the Expression of Intercellular Adhesion Molecule⁃3 in Bone Marrow of Patients with Acute Leukemia and its Clinical Relationship ［D］. Jinan： Shandong University， 2008 （Ch）.

[58] 刘继明，曾慧兰，张毅，等. 白细胞介素8及其受体在急性白血病的表达及其意义［J］. 中华血液学杂志， 1999， 20（1）： 24-26. DOI: 10.3760/j:issn:0253-2727.1999.01.007 .

[59] LIU J M， ZENG H L， ZHANG Y，et al. Study on the expression of interleukin-8 and its receptors in acute leukemia［J］. Chinese Journal of Hematology， 1999， 20（1）： 24-26 （Ch）. DOI: 10.3760/j:issn:0253-2727.1999.01.007 .

[60] KONOPLEVA M， MIKHAIL A， ESTROV Z， et al. Expression and function of leptin receptor isoforms in myeloid leukemia and myelodysplastic syndromes： Proliferative and anti-apoptotic activities［J］. Blood， 1999， 93（5）： 1668-1676. DOI: 10.1182/blood.v93.5.1668.405a15_1668_1676 .

基金资助

国家自然科学基金(U1936205)

国家自然科学基金(61972100)

国家自然科学基金(62006051)