¹青岛大学，山东 青岛

²北京科技大学，北京

收稿日期：2021年6月7日；录用日期：2021年6月28日；发布日期：2021年7月12日

摘要

为了探究二型糖尿病早期动态变化，本文通过分析小鼠标本的组织特异性表达数据，获得相邻时间段之间的差异表达的证据，以表征基因差异表达的动态变化。该数据集中含有丰富的时空信息，但以往的研究中很难充分利用。我们通过在潜在状态上指定马尔可夫随机场(MRF)合并具有时空结构的复杂数据，用蒙特卡洛期望最大化(MCEM)算法计算模型参数，其关键特征是同时考虑基因表达水平的空间相似性和时间依赖性，仿真研究与实例分析结果都表明该方法具有更高的灵敏度，能识别更多的DE基因，能从数据中提取更多生物学上有意义的结果。

关键词

二型糖尿病早期，马尔可夫随机场，MCEM算法

A Markov Random Field-Based Approach for Early Diagnosis of Type 2 Diabetes

Xuehan Tian¹, Yishu Wang^2*

¹Qingdao University, Qingdao Shandong

²University of Science & Technology Beijing, Beijing

Received: Jun. 7^th, 2021; accepted: Jun. 28^th, 2021; published: Jul. 12^th, 2021

ABSTRACT

In order to explore the early dynamic changes of type 2 diabetes mellitus, this paper analyzed the tissue-specific expression data of mouse specimens to obtain the evidence of differential expression between adjacent time periods, so as to characterize the dynamic changes of differential gene expression. This dataset contains abundant spatiotemporal information, but it is difficult to make full use of it in previous studies. We pass in the potential state specified on Markov random field (MRF) combined with space-time structure of complex data, using the Monte Carlo calculation model for the expectation maximization (MCEM) algorithm parameters, its key characteristic is also considering gene expression level of similarity space and time dependence, simulation research and the example analysis results show that this method has higher sensitivity. More DE genes can be identified and more biologically significant results can be extracted from the data.

Keywords:Type 2 Diabetes, Markov Random Field, MCEM Algorithm

Copyright © 2021 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

糖尿病是一种难以根治的慢性疾病，慢性并发症可遍及全身重要器官，急性并发症严重时可危及生命。其中二型糖尿病(T2D)占糖尿病患者总数的九成以上，糖尿病早期由于因胰岛β细胞的代偿机制代偿机制，血糖变化与正常人群差别较小，做专项检查成本高，容易造成医疗资源浪费。因此，诊断出“糖尿病早期”并积极采取干预治疗手段，将有利于推进治愈二型糖尿病这一难题。

自20世纪90年代中后期，基因微阵列技术问世以来便受到科学界的广泛关注，是基因组学中能有效地测量有机体的基因表达水平所使用最广泛的工具。揭示生物网络是系统生物学中的一个关键目标，现在我们能通过微阵列在各种条件下测量基因组水平的表达情况，这些数据可以从统计上推断基因调控网络。在生物学中，许多不同的生物过程都用图表示，比如：包括蛋白质–蛋白质相互作用网络(PPI)，转录调控网络和基因共表达网络，生物途径等。

我们用由节点 $V=\left\{1,\cdots ,p\right\}$，边 $E\subset V\times V$ 组成的图 $G=\left(V,E\right)$ 表示生物网络。对于条件独立图，当且仅当节点i和节点j在给定所有其他节点都不条件独立的情况下， $\left(i,j\right)\in E$。高斯图形模型(GGMs)则已经被证明是推断条件独立图的最佳方法之一。高斯图模型是一种概率图模型，GGMs起源于20世纪70年代早期，这一概念源于AP Dempster [1] 的经典论文，是研究基因之间关联网络的工具。在GGM中，

假定p维随机变量 $X=\left(X_1,\cdots ,X_p\right)$ 服从多元高斯分布 ${\rm N}\left(\mu ,\Sigma \right)$，$\Theta ={\Sigma }^{-1}$，可以证明 $X_i$ 和 $X_j$ 的条件独立性等效于 ${\Theta }_{ij}\ne 0:X_i\perp X_j|X_{V\backslash \{i,j\}}\iff {\Theta }_{ij}=0$。在GGM中，估计条件独立图等效于估计 $\Theta$ 中的非零条目。

为了估计GGM，Meinshausen & Buhlmann [2] 证明了使用套索的邻域选择是一个计算吸引力的标准协方差稀疏高维图选择模型。邻域选择分别估计了图中的每个节点的条件独立性限制，因此等价于高斯线性模型的变量选择。同时证明了稀疏的邻域选择方案对于高维图是一致，其一致性取决于惩罚参数的选择。Yuan & Lin [3] 提出了一种有效的算法来扩展Lasso方法的特征选择，并表明这些扩展具有优越的性能，在一定条件下，Lasso变量的选择具有一致性。Friedman et al. [4] 考虑用逆协方差矩阵上的Lasso惩罚来估计稀疏图的问题.，利用Lasso的坐标下降过程，开发了一种简单的图形Lasso算法，命名为图形套索(glasso)。Cai et al. [5]，Dobra et al. [6]，Wang et al. [7] 等三人同样对估计GGM做出了许多贡献。

本文提出的工作是基于对小鼠标本的组织特异性表达数据的分析所激发的，我们希望通过分析小鼠标本的组织特异性表达数据，在小鼠发育的不同时期，不同的器官组织基因表达动态发育过程不同，为了有效利用该数据集中丰富的时空信息，我们的方法基于马尔可夫随机场(MRF)模型识别相邻时间点的差异表达基因(DE)，用蒙特卡洛期望最大化(MCEM)算法计算模型参数，考虑不同时间不同器官组织的病理性变化。我们的方法的关键特征是同时考虑基因表达水平的空间相似性和时间依赖性，获得相邻时间段之间的差异表达的证据，以表征基因差异表达的动态变化，从而探究二型糖尿病早期动态变化。以更好地从数据中提取生物学上有意义的结果。

2. 模型构建

2.1. 差异基因的潜状态模型

对于二型糖尿的早期患者，因胰岛β细胞的代偿机制，血糖波动于正常人群差异很小，通过现有检测方式诊断时患者已经不可能自愈。为了探究二型糖尿病早期动态变化，我们通过分析小鼠标本的组织特异性表达数据，获得相邻时间段之间的差异表达的证据，以表征基因差异表达的动态变化。

首先我们假设 $B=3$ 为脂肪，肝脏及胰腺这三处采样组织， $T=3$ 为小鼠发育的第1周、9周、18周三个年龄阶段， $G=45101$ 为基因数。令 $n_{bt}$ 表示在区域b处，年龄为t的小鼠的重复次数， $N_b={\left(n_{b1},\cdots ,n_{bt},\cdots ,n_{bT}\right)}^{\prime }$ 为区域b中重复次数的列向量， $N=\left(N_1,\cdots ,N_b,\cdots ,N_B\right)$ 为重复次数

我们假设 $y_{bgtk},k=1,\cdots ,n_{bt}$ 服从均值为 ${\mu }_{bgt}$ 和方差为 ${\sigma }_0^2$ 的正态分布 $y_{bgtk}~{\rm N}\left({\mu }_{bgt},{\sigma }_0^2\right)$。 $x_{bgt}$ 为二进制潜状态，表示基因g是否在区域b时刻t中表达，以 $x_{bgt}$ 为条件，我们假设 ${\mu }_{bgt}$ 遵循高斯分布：

$\begin{array}{l}{\mu }_{bgt}|x_{bgt}=0~{\rm N}\left({\mu }_{1b},{\sigma }_{1b}^2\right)\\ {\mu }_{bgt}|x_{bgt}=1~{\rm N}\left({\mu }_{2b},{\sigma }_{2b}^2\right)\end{array}$

我们假设在特定的区域，混合成分的均值和方差不同，用 ${\mu }_1,{\mu }_2,{\sigma }_1,{\sigma }_2$ 表示所有区域的参数向量。 $y_{bgtk}$ 在 $x_{bgt}$ 条件下的分布具有以下形式：

$\begin{array}{l}{y}_{bgtk}|x_{bgt}=0~N\left({\mu }_{1b},{\sigma }_{1b}^2+{\sigma }_0^2\right)\\ y_{bgtk}|x_{bgt}=1~N\left({\mu }_{2b},{\sigma }_{2b}^2+{\sigma }_0^2\right)\end{array}$

给定潜变量数组 $X$，假定 $f\left(Y|X\right)={\displaystyle \underset{b=1}{\overset{B}{\prod }}{\displaystyle \underset{g=1}{\overset{G}{\prod }}{\displaystyle \underset{t=1}{\overset{T}{\prod }}f\left(y_{bgt}|x_{bgt}\right)}}}$ 条件独立，则

$f\left(y_{bgt}|x_{bgt}\right)=\underset{k=1}{\overset{n_{bt}}{\prod }}f\left(y_{bgt}|x_{bgt}\right)$.

下面我们对于差异表达基因(DE)的分析，我们首先将观察到的数据转换为数组 $B\times G\times \left(T-1\right)$，在相邻时段之间执行t检验并将t统计量转换为 $z~\text{scores}$。设 $y_{bg(t-1)}$ 和 $y_{bgt}$ 分别为t-1与t时刻区域b中基因g的表达值的向量，双样本t统计量为：

$t_{bg(t-1)}=\frac{{\bar{y}}_{bgt}-{\bar{y}}_{bg(t-1)}}{s}$

其中s为 ${\bar{y}}_{bgt}-{\bar{y}}_{bg(t-1)}$ 的标准误差的估计值。然后将检验统计量 $t_{bg(t-1)}$ 转换为 $z~\text{scores}$ ：

$z_{bg(t-1)}={\Phi }^{-1}\left(F_{n_{bt}+n_{b(t-1)-2}}\left(t_{bg(t-1)}\right)\right)$.

其中 $n_{bg(t-1)}$ 和 $n_{bgt}$ 是 $y_{bg(t-1)}$ 和 $y_{bgt}$ 中重复的数目。

基因表达数据由 $B\times G\times \left(T-1\right)$ 数组 $Z$ 表示。 $s_{bgt}$ 表示表示基因g在区域b是否在时段t和t+1之间差异表达的二进制潜状态，这是我们推断的目的。设 $S$ 为维数 $B\times G\times \left(T-1\right)$ 的潜状态数组。在 $s_{bgt}$ 条件下，我们假设 $z_{bgt}$ 遵循混合分布：

$f\left(z_{bgt}|s_{bgt}\right)=\left(1-s_{bgt}\right)f_0\left(z_{bgt}\right)+s_{bgt}{f}_1\left(z_{bgt}\right)$.

其中 $f_0\left(z\right)$ 是零密度， $f_1\left(z\right)$ 是非零密度。我们假设零密度遵循标准的 $N\left(0,\text{1}\right)$ 分布。我们采用非参数经验贝叶斯框架(Efron [8] )，通过使用R包locfdr将非空密度与自然样条拟合。给定 $S$，假定条件独立：

$f\left(Z|S\right)={\displaystyle \underset{b=1}{\overset{B}{\prod }}{\displaystyle \underset{g=1}{\overset{G}{\prod }}{\displaystyle \underset{t=1}{\overset{T-1}{\prod }}f\left(z_{bgt}|s_{bgt}\right)}}}$.

2.2. MRF模型的先验概率

上述模型和推理目标中的一个关键组成部分是我们未知的潜在状态数组 $X$,下面我们用MRF模型生成 $X$ 的先验值 $p\left(X\right)$。对于每个基因 $g$，我们构建一个无向图 $G_g=\left(V_g,E_g\right)$，其中

$V_g=\left\{x_{bgt}:b=1,\cdots ,B,t=1,\cdots ,T\right\}$ 是节点的集合， $E_g$ 是边的集合， $E_g$ 由 $E_{g\text{1}}$ 与 $E_{g\text{2}}$ 组成，

$E_{g1}=\left\{\left(x_{bgt,}{x}_{b^{\prime }g{t}^{\prime }}\right):b\ne b^{\prime },t=t^{\prime }\right\}$ 为包含捕获的区域之间空间依赖性的边缘 $E_{g\text{2}}=\left\{\left(x_{bgt,}{x}_{b^{\prime }g{t}^{\prime }}\right):b=b^{\prime },|t-t^{\prime }|=1\right\}$

为包含捕获的相邻周期之间时间依赖性的边缘。

对于 $p\left(X\right)$ 的联合分布概率，我们建立下列成对交互的MRF模型(Besag [9], Lin et al. [10] )：

$\begin{array}{l}p\left(X|\Phi \right)\propto {\displaystyle \underset{g=1}{\overset{G}{\prod }}\exp \{{\gamma }_0{\displaystyle \underset{V_g}{\sum }{I}_0}\left(x_{bgt}\right)+{\gamma }_1{\displaystyle \underset{V_g}{\sum }{I}_1}\left(x_{bgt}\right)}\\ +{\beta }_1{\displaystyle \underset{E_{g1}}{\sum }\left[I_0\left(x_{bgt}\right)I_0\left(x_{b^{\prime }g{t}^{\prime }}\right)+I_1\left(x_{bgt}\right)I_1\left(x_{b^{\prime }g{t}^{\prime }}\right)\right]}\\ +{\beta }_2{\displaystyle \underset{E_{g2}}{\sum }\left[I_0\left(x_{bgt}\right)I_0\left(x_{b^{\prime }g{t}^{\prime }}\right)+I_1\left(x_{bgt}\right)I_1\left(x_{b^{\prime }g{t}^{\prime }}\right)\right]}\}\end{array}$ 2-(1)

其中 $I_0()$ 与 $I_{\text{1}}()$ 为指标函数，令 $\gamma ={\gamma }_{\text{1}}-{\gamma }_{\text{0}}$，可得条件概率：

$\begin{array}{l}p\left(x_{bgt}|X;\Phi \right)=\frac{\exp \left\{x_{bgt}F\left(x_{bgt},\Phi \right)\right\}}{1+\exp \left\{F\left(x_{bgt},\Phi \right)\right\}}\\ F\left(x_{bgt},\Phi \right)=\gamma +{\beta }_1{\displaystyle \underset{b^{\prime }\ne b}{\sum }\left(2x_{b^{\prime }gt}-1\right)}+{\beta }_2\left\{I_{t\ne 1}\left[2x_{bg(t-1)}-1\right]+I_{t\ne T}\left[2x_{bg(t+1)}-1\right]\right\}\end{array}$ 2-(2)

其中 ${\beta }_1$ 为捕获空间依赖性的参数， ${\beta }_2$ 为捕获时间间依赖性的参数。

接下来，我们考虑先验分布 $p\left(S\right)$ 的MRF模型，同时考虑了时间依赖性和空间相似性。我们将3个区域分成两组：用 $B_p$ 表示胰腺区域， $B_n$ 非胰腺区域。联合概率类似于2-(1)，可以计算条件概率并具有以下形式：

$p\left(s_{bgt}|S/s_{bgt};{\Phi }_{DE}\right)=\frac{\exp \left\{s_{bgt}{F}_{DE}\left(s_{bgt},{\Phi }_{DE}\right)\right\}}{1+\exp \left\{F_{DE}\left(s_{bgt},{\Phi }_{DE}\right)\right\}}$ 2-(3)

如果 $b\in B_p$ ：

$\begin{array}{c}{F}_{DE}\left(s_{bgt},{\Phi }_{DE}\right)={\gamma }_{DE}+{\beta }_{pp}{\displaystyle \underset{b^{\prime }\in B_{n/b}}{\sum }\left(2s_{b^{\prime }gt}-1\right)}+{\beta }_{pn}{\displaystyle \underset{b^{\prime }\in B_n}{\sum }\left(2s_{b^{\prime }gt}-1\right)}\\ +{\beta }_t\left\{I_{t\ne 1}\left[2s_{bg(t-1)}-1\right]+I_{t\ne T}\left[2s_{bg(t+1)}-1\right]\right\}\end{array}$

如果 $b\in B_n$ ：

$\begin{array}{c}{F}_{DE}\left(s_{bgt},{\Phi }_{DE}\right)={\gamma }_{DE}+{\beta }_{nn}{\displaystyle \underset{b^{\prime }\in B_{n/b}}{\sum }\left(2s_{b^{\prime }gt}-1\right)}+{\beta }_{np}{\displaystyle \underset{b^{\prime }\in B_p}{\sum }\left(2s_{b^{\prime }gt}-1\right)}\\ +{\beta }_t\left\{I_{t\ne 1}\left[2s_{bg(t-1)}-1\right]+I_{t\ne T}\left[2s_{bg(t+1)}-1\right]\right\}\end{array}$.

其中 ${\Phi }_{DE}=\left({\beta }_{pp},{\beta }_{nn},{\beta }_{pn},{\beta }_{np}\right)$，${\beta }_{pp}$ 为胰腺区之间的空间系数， ${\beta }_{nn}$ 为非胰腺区之间的空间系数， ${\beta }_{pn}$ 为胰腺区到非胰腺区之间的空间系数， ${\beta }_{np}$ 为非胰腺区到胰腺区之间的空间系数，由对称性可得， ${\beta }_{pn}={\beta }_{np}$。

3. 参数估计与后验概率估计

3.1. 参数估计

3.1.1. p(X)的参数估计算法

估计MRF参数 $\Phi =\left({\gamma }_1,{\beta }_1,{\beta }_2\right)$ 和高斯混合模型参数 $\Theta =\left({\mu }_1,{\sigma }_1;{\mu }_2,{\sigma }_2\right)$ ：

如果给定数组X， $\Phi$，$\Theta$,可以通过最大似然估计(MLE)进行估计。但是，潜伏状态X是未观察到的，因此也需要进行估计。通常针对缺失数据估计实施期望最大化(EM)算法，但由于期望项难以控制，因此不适用于我们的模型。因此，我们提出以下MCEM算法 [11] 来估计：

1) 通过无偏估计量估计 ${\sigma }_0$ ：

${\hat{\sigma }}_0{}^2\frac{1}{G\times {\displaystyle {\sum }_{b=1}^B{\displaystyle {\sum }_{t=1}^T\left(n_{bt}-1\right)}}}{\displaystyle \underset{g=1}{\overset{G}{\sum }}{\displaystyle \underset{b=1}{\overset{B}{\sum }}{\displaystyle \underset{t=1}{\overset{T}{\sum }}{\displaystyle \underset{k=1}{\overset{n_{bt}}{\sum }}{\left(y_{bgtk}-{\bar{y}}_{bgt}\right)}^2}}}}$

2) 在不考虑空间和时间依赖性的情况下，通过简单的高斯混合模型获得初始估计值 $\hat{X}$ 和。

3) 因为没有针对的显式MLE，所以选择初始估计，将该估计将最大化得到以下伪似然函数(Besag [12] )：

4) 令，EM算法中预期的完整数据对数似然可通过蒙特卡罗方法(Wei and Tanner [11] )进行近似：

3-(1)

其中通过吉布斯采样获得。从到，中的所有条目都被下列公式按顺序更新，

3-(2)

5) 用更新，使3-(1)最大化：

与步骤3相同，我们用中的伪似然函数替换似然函数。包含，的项是可分离的，因此可以分别进行优化。

6) 重复步骤4，步骤5，直至收敛。

3.1.2. p(S)的参数估计算法

为了随着时间的推移识别DE基因。在该模型中，只需要迭代地更新MRF先验中的参数。该算法与上面的算法有一些相似之处：

1) 汇总Z中的并通过locfdr程序估算。

2) 不考虑空间和时间依赖性的情况下，通过简单的混合模型获得初始估。

3) 获得初始估计，其最大化的伪似然函数：

其中如2-(3)所定义。

4) 通过蒙特卡罗近似估计预期的完整数据对数似然：

3-(3)

其中由Gibbs采样获得，从到中的所有条目都下列公式被按顺序更新。

3-(4)

5) 通过更新，直至得到最大化3-(3)。

6) 重复步骤4，步骤5，直至收敛。

3.2. 后验概率估计

为了获得后验概率的估计，我们实现了一个单独的Gibbs采样，并通过MCEM算法估计模型参数。潜在状态根据3-(2)和3-(4)依次更新，关于对DE基因的推断，我们采用基于后验概率的FDR定义(Newton

et al. [13]；Li，Wei and Maris [14] )。局部后验f.d.r由Gibbs采样估计得来。为的升序排序值。让，然后我们拒绝所有零假设，。

4. 仿真模拟与实例分析

4.1. 仿真模拟

为了评估我们提出的MRF模型的性能本节进行了仿真模拟。在模拟中通过3轮重复次数生成了在3个区域，3个时间点的500个基因的数据。我们考虑了两种模拟设置模拟的：

模拟设置1：用Gibbs采样器模拟潜状态阵列。采样器从一个可能是EE或DE随机数组开始，在每一轮吉布斯抽样中，每个潜在状态都根据3-(7)顺序更新一次。我们进行了三轮吉布斯抽样，以获得潜在的状态数组S。参数为：，，，，。

模拟设置2：在时间点1中，所有基因的死亡概率均为0.1。隐马尔可夫状态随时间模型改变：如果一个基因在t−1时刻为DE，t时刻变成EE概率为0.5；为了保持DE基因的数量不变，我们随机选择相同数量的EE基因在t时刻切换到DE。三层的潜在状态最初设置是相同的。然后随机选择DE状态的不同比例(0.1、0.2、0.5)切换为EE；相同数量的EE状态也切换为DE，以保持DE基因总数不变。

我们在敏感性、特异性和FDR这三个方面比较了这两个模型(见表1)。对于MRF和EB模型，我们选择了，满足，控制FDR小于0.05，与EB相比，我们的方法(MRF)了更高的灵敏度。而当相似度降低时，即，从HMM(0.1)到HMM(0.5)，FDR在我们的模型中略有增加，但在EB模型中并没有完全增加。

表1. MRF与EB的比较

4.2. 实例分析

本文数据来自小鼠标本的组织特异性表达数据数据集，编号为GSE77943。

我们应用MRF模型来识别相邻时期之间的DE基因。运行20次MCEM算法的迭代，Gibbs采样器的设置为500/1500，估计的MRF参数为，，，，。胰腺与非胰腺区域系数远小于胰腺与胰腺和非胰腺与非胰腺系数，表明胰腺与非胰腺区域之间的组间差异。当假设没有空间和时间依赖性时，模型简化为简单的经验贝叶斯(EB)模型，基于局部FDR控制程序，MRF和EB模型中的阈值分别为0.46和0.32。在两个模型中被鉴定为DE的基因数量是30,605 (MRF)和13,273 (EB)，其中11,149 (84%)重叠。较高的阈值导致在MRF模型中识别为DE的更多基因。

后验概率对超参数的选择敏感，我们通过改变后局部FDR的阈值来计算灵敏度和特异性，我们将本文提出的MRF模型与EB模型进行了比较，该模型假设没有时间和空间依赖性。对于两组分别绘制ROC曲线如图1所示，与EB模型相比，MRF模型AUC面积更大，MRF模型的性能更好，MRF模型从空间相似性中获益更多。

图1. EB与MRF模型的ROC曲线

5. 结论

图形模型在基因表达数据的分析中应用广泛，基于基因表达数据的图形模型可以为可视化基因之间的关系和生成生物假设提供一个有用的工具，我们对图形模型之间的差异很感兴趣。在本文中，我们使用了一种贝叶斯邻域选择程序来估计高斯图形模型，结合了马尔可夫随机场来合并具有时空结构的数据。仿真研究和实例分析表明，将复杂的数据结构合并到联合建模框架中有助于估计，与常用的高斯混合模型相比，该模型具有灵敏度更高，能识别更多的DE基因。

文章引用

田雪寒,王艺舒. 基于马尔科夫随机场对二型糖尿病早期诊断
A Markov Random Field-Based Approach for Early Diagnosis of Type 2 Diabetes[J]. 应用数学进展, 2021, 10(07): 2339-2347. https://doi.org/10.12677/AAM.2021.107244

参考文献