Comparison of Methods for Differential Co-expression Analysis for Disease Biomarker Prediction2.286

Comput Biol Med.2019 Oct;113:103380. cdoi: 10.1016/j.compbiomed.2019.103380. Epub 2019 Aug 10.

Abstract

Keywords: Alzheimer's disease; Differential co-expression analysis; Disease biomarkers; Empirical study; Gene expression; Parkinson's disease; miRNA expression.

近年来，芯片和RNA序列技术的广泛使用已导致与基因表达和miRNA表达有关的生物学数据的快速增长，从而增加了对用于更好的转录组分析的计算方法和工具的需求。这些计算方法有助于破译疾病的复杂机制，从而有可能帮助开发有效的药物靶标，以及在早期阶段预防疾病。

在这些计算方法中，基因共表达网络（gene co-expression network，gCEN）广泛用于功能预测，因为共表达网络（co-expression network，CEN）具有挖掘具有相似共表达模式的功能相关基因的能力，擅长在一组样本中鉴定功能相关的基因和miRNA。同样，miRNA共表达网络（miRNA co-expression networks，miCEN）用于研究miRNA与调控基因的相关性，并确定其在表达失调，影响疾病调控等活动中的作用。然而，利用miCEN鉴定疾病相关的mirna的方法尚未得到充分的探索，因此作为构建mi CEN方法的选择变得至关重要。

但是，CEN分析中的一个问题是它无法识别在疾病条件下会发生条件特定变化的基因或miRNA。为了解决相互依赖的问题，出现了一种基于共表达分析的方法，称为差异共表达（differential co-expression，DC）分析，其提供了对控制和疾病状况之间改变机制的见解。

该分析背后的主要思想是，共表达或共调控的基因或miRNA具有相似的表达模式，并且持续在疾病发展过程中，这些共调控的基因或miRNA经历突变，这可能导致它们与其他基因或miRNA的关联发生改变，从而导致跨疾病条件的共调控基因的表达模式有所差异。

下面我们将比较研究四种共表达网络模块检测方法
1. THD-Module Extractor
2. DiffCoEx
3. MODA
4. WGCNA
将这些方法运用到阿尔茨海默病和帕金森病的研究中对基因和miRNA表达数据（3个microarray 和4个RNA-seq数据集）进行如下分析比较
1. 数据预处理，如归一化、补缺失值等
2. 采用上述四种方法机分别进行共表达网络构建和模块挖掘
3. 提取基因或miRNA共表达的模块后，我们用差异共表达分析扩展了共表达分析，具体包括拓扑富集分析、GO/KEGG富集分析、miRNA靶基因富集、验证模块保守性等。
4. 识别hub基因和miRNA
5. 识别疾病的标志物
6. 疾病调控网络分析

方法：
有多种gCEN模块提取技术可以推断基因功能并从基因表达数据中识别疾病与基因的关联。但由于miRNA的长度短，miRNA数据集的变异性高和尺寸小，因此很少有CEN模块提取技术。在我们的研究中，我们使用四种模块提取技术，即THD-Module Extractor，DiffCoEx，MODA和WGCNA，来提取基因和miRNA共表达的模块。
在这些技术中，THD-Module Extractor和WGCNA不假定对照和疾病状况之间存在异常。为了解决这个问题，我们扩展了这些传统的CEN技术，使用DiffCoEx和MODA的差异共表达功能，在控制和疾病条件下识别基因和miRNA之间的变异。

THD-Module Extractor
THD-Module Extractor是一种CEN模块提取方法，在疾病相关分析中考虑边界基因（具有单条边的基因）以及其他共表达的基因。该方法从共表达模块中挖掘出有趣的边界基因，这些边缘基因与核心基因具有较高的语义相关性和较低的表达相似性。该方法可以应用于微阵列AD数据集以鉴定疾病生物标志物。该方法接受两个阈值，即表达相似性阈值（δ）和最小邻域阈值（ρ），从而从基因表达数据构建CEN。并且用户可以调节共表达模块中的最小基因数和生成的模块数。该方法已在MATLAB和R中实现，适用于微阵列数据集，已扩展为通过在疾病状况之间寻找基因和miRNA的协方差来找到差异共表达的基因和miRNA。

DiffCoEx
DiffCoEx基于WGCNA构建CEN，并遵循五个步骤来鉴定差异共表达的基因。首先，对于每个条件，利用Pearson相关构造一个邻接矩阵，然后计算条件对的邻接差矩阵。从邻接差矩阵出发，利用拓扑重叠得分计算不同矩阵，其中高得分表示具有相同邻居的基因，最后，从不同矩阵中找到共表达的模块。该方法取决于对用于构造邻接矩阵和tree cutting阈值的参数的调整。该方法在R中实现，仅用于微阵列数据集。

MODA
MODA（MOdule Differential Analysis for weighted gene co-expression network）是基于WGCNA的CEN模块提取技术，用于识别与疾病相关的基因差异表达模块。MODA使用一种节省样本的方法从单个或多个样本构建CEN，并通过比较跨条件的网络来识别跨条件的差异共表达的基因模块。为了检测来自每个网络的模块，通过调整tree cutting阈值来应用分层聚类。最佳切割阈值是基于加权网络的平均模块化确定的。该方法在R中实现，结果受参数值的影响很大。

WGCNA
加权基因共表达网络分析（WGCNA）使用Pearson相关来计算所有样本或条件下基因之间的相关模式。WGCNA通过使用模块特征基因对每个模块进行汇总，使用层次聚类从CEN使用不同的tree cutting阈值来识别模块。特征基因通过计算模块的成员关系来帮助关联模块。WGCNA是使用最广泛的CEN技术，并用于许多应用程序，例如癌症的遗传分析，小鼠和酵母的基因组分析以及脑MRI数据分析。该方法在R中实现，可用于微阵列基因表达数据和RNA-Seq数据。要找到差异共表达的基因或miRNA，我们使用DGCL软件包中的WGCNA()函数。
表1突出显示了四种CEN模块提取技术的主要功能

在我们的研究中，我们使用具有基因和miRNA表达谱的七个数据集（四个RNA-seq数据集和三个microarray数据集）来进行模块的（差异）共表达分析。我们扩展了差异共表达基因的分析，包括GO/KEGG富集分析，拓扑富集（针对差异共表达的基因），miRNA靶标富集（针对差异共表达的miRNA）和中枢基因的识别。从差异共表达的基因或miRNA模块中提取的中枢基因在整个疾病条件下都会发生调节变化。根据网络统计数据，如中心基因的程度，我们确定了中心基因和中心mirna，它们在疾病网络中发生定向障碍。我们进一步评估条件特定的枢纽基因和枢纽miRNA，并确定它们在AD和PD进展中的潜在作用。

表2概述了我们实验研究中使用的微阵列和RNA-seq数据集，包括microarray数据（两个基因表达AD数据集（GSE4757，GSE5281）和一个miRNA表达AD数据集（GSE16759））和RNA-seq数据（两个基因表达PD数据集（GSE62642，GSE68719）和两个miRNA表达PD数据集（GSE72962，GSE77668））。并对microarray数据和RNA-seq数据进行预处理。

结果：
我们将对GO富集分析，拓扑富集分析，miRNA靶标富集和KEGG通路富集分析等方面的四种模块提取方法（THD-Module Extractor，DiffCoEx，MODA和WGCNA）进行评估。此外，我们分别通过识别GSE68719和GSE72962 PD数据集的对照和疾病网络中的中枢基因和中枢miRNAs来扩展分析。我们评估了从GSE4757 AD数据集中获得的差异共表达基因中鉴定的hub基因的基因调控网络，以及在GSE4757 AD数据集的对照和疾病网络中差异共表达基因模块的模块保守性。

1.差异共表达基因模块
1.1.GO和拓扑富集分析
我们研究了从AD和PD数据集中提取的GO term和共表达模块的拓扑意义。对于GO富集分析，p值定义为在与特定GO项相关的大小为n的簇中获得k个或更多基因的可能性，而q值它提供了最小的错误发现率（FDR）。p和q的值越小，表示模块的统计意义越高，即该模块不是偶然发现的。对于拓扑属性，如最短路径长度（SPL），节点的介数（NB）和连通度（Degree）等，确定模块的拓扑富集情况。
对于AD数据集：
THD-Module Extractor的性能：对于GSE4757 AD数据集，提取的最佳差异共表达基因模块的p值和q值分别为1.96E-06和1.38E-03。该模块还拓扑属性SPL，NB，Degree值分别为3.97，20666和11.19。对于GSE5281 AD数据集，提取出差异共表达基因模块的p、q值分别为7.72E-07、9.50E-03，拓扑属性SPL，NB，Degree度值分别为4.21、11895、8.03。
DiffCoEx的性能：对于GSE4757 AD数据集，提取的最佳差异共表达基因模块的p值和q值分别为2.79E-05和2.75E-02。相对而言，该模块的拓扑属性SPL，NB，Degree值分别为4.03、20585和10.79。同样，对于GSE5281 AD数据集，提取的最佳差分共表达模块的p、q值分别为5.10E-06、3.06E-02，模块的拓扑属性SPL，NB，Degree值分别为4.06、21094和10.58。
MODA的性能：从GSE4757 AD数据集中获得的最佳差异共表达基因模块的p、q值分别为6.85E-06、2.45E-03，其拓扑属性SPL，NB，Degree值分别为、3.99、44151和21。对于GSE5281 AD数据集，提取最佳差分共表达模块的p、q值分别为1.03E-04、2.85E-02，而拓扑属性SPL，NB，Degree值分别为3.92、10897和9.03。
WGCNA的性能：从GSE5281 AD数据集中提取的最佳差异共表达基因模块，其p、q值分别为6.71E-05、1.36E-02，其拓扑属性SPL，NB，Degree值分别为4.19、2425、4.94。同样，对于GSE5281 AD数据集，提取的差异共表达基因模块的p、q值分别为7.10E-06、2.24E-03，而拓扑属性SPL，NB，Degree值分别为4.03、9670、8.41。
比较:我们分别从GSE4757和GSE5281 AD数据集中提取的差异共表达基因模块，对其p值、q值、SPL、NB、Degree进行了比较。在我们对GSE4757的实验研究中，我们发现THD-Module Extractor能够获得最具统计意义和拓扑属性的差异共表达模块。来自MODA的差异共表达模的拓扑值也很重要。同样，对于GSE5281数据集，我们观察到使用THD-Module Extractor提取出统计意义上最具差异的共表达基因模块。

对于PD数据集进行同样的分析比较：
比较: 我们观察到，对于GSE62642， DiffCoEx获得的差异共表达的基因模块的p值和q值分别为2.24E-07和2.93E-03，具有较高的统计学意义，优于其他三种方法。MODA获得的差异共表达的基因模块在拓扑属性上有更好的表现，SPL值较低为4.26，NB值较高为50394，度较高为18.33。对于WGCNA，我们没有发现统计学意义和拓扑丰富的差异共表达基因模块。另一方面，利用WGCNA提取GSE68719数据集在统计学意义上的最佳差异共表达基因模块，其p和q值分别为4.58E-08和5.79E-05。与其他三种方法提取的差异共表达模组相比，利用MODA从GSE68719中提取的差异共表达模组在拓扑统计(即SPL、NB和Degree)方面具有最佳的价值。

在表3中，我们总结了根据上述GO和拓扑富集分析，使用THD-Module Extractor，DiffCoEx，MODA和WGCNA提取的最佳差异共表达基因模块。

1.2.通路富集分析
通路富集分析是对在疾病发生过程中发生差异表达的基因或基因组的重要性的测量，从p值和q值上展示了与差异共表达基因模块相关的KEGG通路的重要性。

2.差异共表达miRNA模块
2.1.GO富集和miRNA靶富集分析
对于AD数据集：
比较:我们比较了从GSE16759 AD、GSE72962 PD和GSE77668 PD中提取的差异共表达miRNA模块，从p值和O/E比值上比较了它们的GO富集和miRNA靶富集。
从这四种方法中都可以看出，从GSE16759 AD数据集中提取的这些差异共表达miRNA模块都富含与AD相关的低p值靶基因。然而，在靶基因富集方面， DiffCoEx、MODA和WGCNA的表现较差（GSE16759 AD的样本量较小，不适合用于具有生物学意义的模块提取）。对于THD-Module Extractor，从GSE16759 AD miRNA数据集中映射得到的最佳差异共表达模块的p值和靶基因富集方面O/E比值分别为9.55E-205和54.6154(最高)，而其他差异共表达的模块没有显著的统计和O/E比值。
对于PD数据集进行同样的分析比较：
比较:对于GSE77668 PD miRNA数据集，THD-Module Extractor和WGCNA都有统计学意义的模块。然而，使用WGCNA提取的这个具有统计意义的模块的O/E比率低于THD-Module Extractor。因此，从WGCNA、DiffCoEx、MODA中提取的差异共表达模块的生物学意义低于THD- Module Extractor。同样，对于GSE72962 PD miRNA数据集，DiffCoEx产生了一个最低p值为2.94E-111的模块。然而，这种差异共表达模块的O/E比低于THD-Module Extractor。

2.2.通路富集分析
比较：从GSE72962 PD中提取microrna的模块分析，THD-Module Extractor和DiffCoEx提取的microrna的模块的较好（p值较低）。并根据相关富集的通路信息推断提取的miRNA模块也与生物学意义上的通路或生物学过程相关。

3.识别hub基因/miRNA
模块中连接基因/miRNA最高程度的基因/miRNA被称为该模块的hub基因/miRNA。hub基因和hub miRNA在gCEN和miCEN中的表达度最高，具有较高的生物学意义。对这些hub基因的进一步分析已经证明，这些hub基因在不同的疾病条件下变得不活跃或发生变异，从而导致新神经元细胞死亡和神经退行性疾病，如PD、AD和多发性硬化(MS)。例如，RAD52是DNA修复和保护端粒结构的主要候选基因。在衰老过程中，RAD52的失活或微扰会导致细胞损伤和神经元死亡。

4.hub基因的基因调控网络分析
我们评估了由转运因子介导的hub基因的调控网络。这个调控网络被称为基因调控网络(GRN)，它提供了关于基因-基因相互作用的见解，并帮助识别和确定疾病中的治疗靶点。这些基因调控mRNA和其他蛋白的基因表达，在疾病的形态发生和发病机理中起着至关重要的作用。我们使用GENIE3算法推断调控网络。使用THD-Module Extractor、DiffCoEx、MODA和WGCNA从GSE4757 AD数据集中提取的hub基因的grn的邻接矩阵。我们发现，DDR1、ABCA1、INPP5D等hub基因是基因调控网络中的调控基因，并已被证实与AD相关。

5.评估模块preservation统计信息
差异共表达基因或miRNA模块是条件特异性的，在对照组和疾病组条件之间发生调控变化和基因连接的差异。我们使用模块preservation统计信息来探索AD进展过程中基因共表达网络(gCEN)拓扑结构的变化。通常，模块preservation统计信息用于查找参考和测试网络之间是否保存了基因模块。我们使用modulePreservation()将疾病网络的扰动和失调与控制网络相关联。
下表给出了使用WGCNA从GSE4757 AD数据集中提取的差异共表达基因模块preservation统计信息。黑色，黄色，棕色和红色模块中的差异共表达基因表现出较低preservation Zsummary，即共表达结构的缺失。

类似地，在图B中，我们确定以黑色，黄色，棕色和红色显示的这四个基因模块的preservation Zsummary在2和19之间，这表明该模块在疾病（AD）网络中被适当保存。

结论：

在详细分析的基础上，我们得出结论，所有四种方法均产生了生物学上丰富且具有统计学意义的基因和miRNA差异共表达模块。但是，与其他三种方法在GO富集分析，拓扑意义，通路富集分析，hub基因鉴定，基因调控网络分析和疾病生物标志物鉴定等方面相比，使用THD-Module Extractor提取的差异共表达基因和miRNA模块的性能要好。THD-Module Extractor的工作原理是，它既包括基因表达相似性又包括语义相似性来识别共表达的模块，这导致了挖掘到的差异共表达的基因或miRNA模块可以包括功能上与AD或PD相关的通路相关的重要基因或miRNA。

转自生信人