采用网络分析方法评估不同生态微生物群落（如土壤、水体和根际）之间的关联。本文介绍了一种利用WGCNA算法分析因不同生态环境而产生的微生物群落中可能存在的各种共现网络的实验方案。
需要JoVE订阅才能观看此内容。 请登录或开始免费试用
方法文章
基于加权相关网络的不同生境下根系微生物群落的分化
* These authors contributed equally
本文内容
摘要
摘要
根部微生物组在植物生长和环境适应中发挥着重要作用。网络分析是研究群落的重要工具，能够有效探究不同环境中不同微生物物种之间的相互作用关系或共现模式。本文旨在详细介绍如何利用加权相关网络算法，分析因不同生态环境而在微生物群落中可能出现的各类共现网络。本实验的所有分析均在WGCNA软件包中进行。WGCNA是一个用于加权相关网络分析的R语言软件包。用于演示这些方法的实验数据来自NCBI（国家生物技术信息中心）数据库中的水稻（Oryza sativa）根系三个生态位的微生物群落数据。我们采用加权相关网络算法，分别构建了这三个生态位中微生物群落的共丰度网络，并鉴定了内生圈、根表层和根际土壤之间的差异共丰度网络。此外，通过“WGCNA”软件包获得了网络中的核心属，这些核心属在网络功能中发挥重要的调控作用。这些方法使研究人员能够分析微生物网络对环境扰动的响应，并验证不同的微生物生态响应理论。结果表明，在水稻的内生圈、根表层和根际土壤中均鉴定出显著差异的微生物网络。
引言
微生物组研究对于理解与调控生态系统过程具有重要意义1,4。微生物种群通过相互作用的生态网络彼此关联，这些网络的特征可影响微生物对环境变化的响应3,4。此外，这些网络的特性影响着微生物群落的稳定性，并与土壤功能密切相关5。加权基因共表达网络分析（WGCNA）目前已广泛应用于基因与微生物群落关系的研究6。以往的研究主要关注不同基因或种群网络与外部环境之间的关联7。然而，在不同环境条件下由微生物种群形成的共现网络差异却鲜有研究。本文所呈现的研究旨在提供深入见解和具体细节，介绍如何快速实施WGCNA算法，以构建在不同环境条件下采集的微生物组样本的共现网络。基于分析结果，我们评估了种群的组成及差异，并进一步探讨了不同微生物种群之间的关系。本研究应用了如下加权相关性网络算法的基本流程8：首先，通过计算操作分类单元（OTU）表达谱之间的皮尔逊相关系数构建相似性矩阵；随后，依据无标度拓扑标准选择邻接函数参数（幂函数或S型函数），将相似性矩阵转化为邻接矩阵，每个共现网络对应一个邻接矩阵；接着，采用基于TOM的不相似性指标结合平均连接层次聚类方法，将表达模式一致的OTU聚类为模块；进一步计算保守性统计量与相关参数分析模块之间的关系，最终识别出模块中的核心OTU。这些方法特别适用于分析不同环境条件下多种微生物种群间网络结构的差异。在本文中，我们详细描述了共表达网络构建的方法、模块间差异的分析过程，并简要概述了用于识别不同模块网络中核心物种的操作步骤。
访问受限。请登录或开始试用以查看此内容。
方案
1. 数据下载
- 从NCBI数据库下载登录号PRJNA386367的数据。从该登录号的数据中，选取2014年在美国加利福尼亚州阿布克尔（Arbuckle）的淹水水稻田中种植14周的水稻植株的根际、根表和内生微生物组数据。
注：根际、根表和内生微生物组数据在登录号PRJNA386367中以OTU表形式提供。
2. 最佳功率值的确定
注意：WGCNA 软件包包含以下所有功能参数。WGCNA 是一个用于加权相关网络分析的 R 软件包。主要命令行参见补充材料 S1。
- 在 R 语言环境中，打开 Rstudio 软件并安装 WGCNA 软件包。
- 加载数据，并使用 goodSamplesGenes 函数检查数据的正确性。执行以下命令行：
"gsg = goodSamplesGenes(datExpr0, verbose = 3)
gsg$allOK "
点击 Run。 - 检查是否存在离群值，并保存符合要求的样本。当检查结果为 TRUE 时，继续下一步操作。保存结果。
- 使用 PickSoftThreshold 函数计算两组数据在不同幂值下的无标度指数 R2。执行以下命令行：
"sft = pickSoftThreshold(datExpr0, powerVector = powers, verbose = 5)"
点击 Run。 - 可视化结果（图 1）。执行以下命令行：
"plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
xlab="软阈值（power）",ylab="无标度拓扑模型拟合度，有符号 R^2",type="n",
main = paste("ES_无标度性"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
labels=powers,cex=cex1,col="red");
abline(h=0.9,col="red")
plot(sft$fitIndices[,1], sft$fitIndices[,5],
xlab="软阈值（power）",ylab="平均连接度", type="n",
main = paste("ES_平均连接度"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")"
点击 Run。
注：加权相关网络算法的前提是所构建的共表达网络结构符合无标度拓扑标准，从而提高其鲁棒性。无标度指数越接近 1，表明网络结构越接近无标度网络。 - 选择无标度指数 R2 大于 0.9 时对应的幂值，并进入下一步分析。
注：当无标度指数接近 1 时，网络结构更接近无标度网络。在分析两个或多个网络时，需选择使每个网络均接近无标度网络的幂值，以满足共表达网络之间的可比性。
3. 共表达网络的构建与模块识别
注意：根据上述计算得到的功率值构建共现网络。关键命令行参见补充材料 S2。
- 使用 WGCNA 软件包中的 adjacency 函数为构建符号共现网络添加有符号参数。执行命令行：
"adjacency = adjacency(datExpr0, power = softPower)"
点击 Run。 - 应用 TOM-similarity 函数构建拓扑重叠网络，并计算不相似性网络。执行命令行：
"TOM = TOMsimilarity(adjacency);
dissTOM = 1-TOM"
点击 Run。
注意：添加有符号参数用于设定拓扑重叠网络的类型。 - 使用 hclust 函数选择平均连接层次聚类方法进行层次聚类。执行命令行：
"geneTree = hclust(as.dist(dissTOM), method = "average");"
点击 Run。 - 使用 cutreeDynamic 函数执行动态分支剪切，并将 minClusterSize 参数设置为 30，以获得模块识别结果。执行命令行：
"dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM, deepSplit = 2, pamRespectsDendro = FALSE, minClusterSize = minModuleSize);"
点击 Run。
注意：最小模块大小不得低于 30。 - 通过 moduleEigengenes 函数计算每个 OTU 模块的模块特征值。执行命令行：
"MEList = moduleEigengenes(datExpr0, colors = dynamicColors)
MEs = MEList$eigengenes"
点击 Run。
注意：模块特征值代表该模块中 OTU 的整体表达水平，它并非某个特定 OTU，而是通过奇异网络值分解获得的每个聚类的第一主成分。 - 基于模块特征值的相关系数执行聚类函数。使用 mergeCloseModules 函数合并相关性低于 0.25 的模块。执行命令行：
"merge = mergeCloseModules(datExpr0, dynamicColors, cutHeight = MEDissThres, verbose = 3)"
点击 Run。 - 最后，使用 plotDendroAndColors 函数进行可视化，以获得每个共表达网络的模块分配示意图（图 2）。使用 table 函数从模块分配表中提取每个 OTU 对应的模块归属信息。执行命令行：
"plotDendroAndColors(geneTree, mergedColors, "Merged dynamic",dendroLabels = FALSE,
hang = 0.03,addGuide = TRUE, guideHang = 0.05,
main = "ES_Gene dendrogram and module colors")"
点击 Run。
注意：在共表达网络的模块分配图中，不同颜色代表不同的模块，灰色代表无法归入任何模块的 OTU。灰色模块中 OTU 数量越多，表明表达矩阵前期预处理质量较差。
4. 模块比较
注意：该方法可用于比较两个生态微生物群落的网络模块。本文中，比较内生环境与根面、内生环境与根际、根际与根面之间的微生物网络模块差异。
- 保存性检验
- 加载前几步中保存的两组数据集的参数和结果。
- 将一组微生物数据的网络模块划分结果设为参考组，另一组设为测试组。
- 使用 modulePreservation 函数计算保守性统计参数 Z_summary 和 medianRank 的值。执行命令行：
"system.time({mp=modulePreservation(multiExpr,
multiColor,referenceNetworks=1,
nPermutation=100, randomSeed=1,quickCor=0,verbose=3)})"
点击 Run。
注意：该结果可用于量化模块之间的保守性。Z_summary>10 表示两个模块高度保守，而 Z_summary<2 表示模块未被保守。medianRank 通过排序反映模块的相对保守性，medianRank 值越高表示模块越未被保守。（关键命令行参见 Supplement S3。） - 使用 plot 函数可视化结果（图 3），获取 Z_summary 和 medianRank 参数值（表 1）。
注意：同时满足 Z_summary 值小于 2 且 medianRank 值处于最高水平的网络模块，是两个生态微生物群落中最未被保守的模块。 - 根据上述两个统计参数的结果，确定两个网络中最为未被保守的模块。
- 模块成员关系的相关性分析
- 将两个网络的模块划分结果分别设为参考组和测试组。
注意：设置需与保存性检验一致。 - 使用 corPvalueStudent 函数提取若干候选模块中每个 OTU 的 kME（模块成员度）值。
执行命令行：
"Pvalue = as.data.frame(corPvalueStudent(as.matrix
(ModuleMembership), Samples))"
点击 Run。
注意：kME 表示模块成员度，ME 表示模块特征向量（module eigen），代表模块中 OTU 表达的整体水平。kME 是每个 OTU 与 ME 之间的相关系数，通过 OTU 的 kME 值可量化其在网络中的重要性。（关键命令行参见 Supplement S4。） - 随后，使用 verboseScatterplot 函数计算两个网络中对应 OTU 的 kME 值之间的相关系数，并绘制相关性分析图（图 4）。
执行命令行：
"verboseScatterplot(abs(TModuleMembership
[TmoduleGenes, Tcolumn]),
abs(NModuleMembership[NmoduleGenes, Ncolumn]),
xlab = paste("kME in", "ES"),
ylab = paste("kME in", "RP"),
main = paste("lightyellow"),
cex.main = 1.7, cex.lab = 1.6, cex.axis = 1.6, col = modulecolor)"
点击 Run。 - 选择两个网络中 OTU 的 kME 值相关系数最小的模块，认为该模块在两个网络间的差异最大。
- 将两个网络的模块划分结果分别设为参考组和测试组。
5. 微生物差异网络模块的分析
- 通过统计分析差异最大模块的OTU序列集，获得优势细菌门的数据。
注意：差异最大模块的OTU序列集按门水平分类进行汇总，优势细菌门所占比例超过10%。 - 然后，使用exportNetworkToCytoscape函数获取最大差异模块中OTU相互作用关系信息的文件。
执行命令行：
"cyt = exportNetworkToCytoscape(modTOM,
edgeFile = paste("NEW-ES_CytoscapeInput-edges-", modules , ".txt", sep=""),
nodeFile = paste("NEW-ES_CytoscapeInput-nodes-", modules, ".txt", sep=""),
weighted = TRUE,threshold = 0.5, nodeNames = modProbes,
altNodeNames = modGenes, nodeAttr = moduleColors[inModule])"
点击 运行。 - 将文件导入Cytoscape，设置阈值为0.5，并根据需要调整其他参数。
- 构建差异微生物的共现网络（图5）。
- 获得在网络中具有最重要调控作用的核心属的信息。
注意：根据OTU的kME值，可定义核心属。 - 最后，评估核心属的功能，并分析其对整个差异网络的影响。
访问受限。请登录或开始试用以查看此内容。
结果
本文中的代表性结果下载自NCBI数据库中2014年加利福尼亚Abaker水稻根系微生物组数据（PRJNA386367）9。该数据包括在淹水水稻田中生长14周的水稻植株的根际、根表和内生圈微生物组样本。我们采用WGCNA算法选择满足接近无尺度网络的三个网络的幂值（图1），并构建了三个共表达网络（图2）。在内生圈、根表和根际土壤微生物共表达网络中，分别鉴定了23、22和21个模块。这些结果表明，这三个生态位中的微生物相互作用网络数量基本相等。
我们进一步比较了内生圈、根表层和根际土壤中微生物网络模块的差异。获得了三个生态位组模块的以下保真性检验结果。在根际土壤与根表层之间存在三个极不保守的模块（图3a，表1）；在根际土壤与内生圈之间存在九个极不保守的模块（图3b
访问受限。请登录或开始试用以查看此内容。
讨论
相关性网络在生物信息学应用中正被越来越多地使用。WGCNA 是一种系统生物学方法，用于描述生物系统中各个元素之间关系的分析12。早期关于 WGCNA 的研究使用了 R 软件包13,14,15。该软件包包含用于网络构建、模块检测、拓扑性质计算、数据模拟、可视化以及与外部软件交互的功能。WGCNA 已被广泛应用于分析脑癌16、酵母细胞周期17、小鼠遗传学18,19、灵长类动物脑组织20,21、糖尿病22以及植物
访问受限。请登录或开始试用以查看此内容。
披露
作者无任何利益冲突需要披露。
致谢
本论文的撰写得到了中国国家自然科学基金-贵州省人民政府喀斯特科学研究中心项目（U1812401）、贵州师范大学博士研究项目（GZNUD[2017]1）、贵州省科技支撑计划项目（QKHZC[2021]YB459）以及贵阳市科技计划项目（[2019]2-8）的资助。
作者感谢 Edwards J.A 等人提供了水稻微生物组数据至公共数据库，并感谢 TopEdit（www.topeditsci.com）在本手稿撰写过程中提供的语言学协助。
访问受限。请登录或开始试用以查看此内容。
材料
| 姓名 | 公司 | 目录编号 | 评论 |
|---|---|---|---|
| R | The University of Auckland | version 4.0.2 | R 是一个用于统计计算和图形的自由软件环境，可在多种 UNIX 平台、Windows 和 MacOS 系统上编译和运行。 |
| RStdio | JJ Allaire | version 1.4.1103 | RStudio 集成开发环境是一组集成工具，旨在帮助用户更高效地使用 R 和 Python。 |
| Cytoscape | version 3.7.1 | Cytoscape 是一个开源软件平台，用于可视化复杂网络，并将其与任意类型的属性数据进行整合。 | |
| NCBI database | 美国国家生物技术信息中心通过提供生物医学和基因组信息，推动科学与健康事业的发展。 |
参考文献
- Philippot, L., Raaijmakers, J. M., Lemanceau, P., vander Putten, W. H. Going back to the roots: the microbial ecology of the rhizosphere. Nature Reviews Microbiology. 11, 789-799 (2013).
- Fierer, N. Embracing the unknown: disentangling the complexities of the soil microbiome. Nature Review Microbiology. 15 (10), 579-590 (2017).
- Jin, J., Wang, G. H., Liu, X. B., Liu, J. D., Chen, X. L., Herbert, S. J. Temporal and spatial dynamics of bacterial community in the rhizosphere of soybean genotypes grown in a black soil. Pedosphere. 19 (6), 808-816 (2009).
- Ma, B., et al. Genetic correlation network prediction of forest soil microbial functional organization. ISME J. 12 (10), 2492-2505 (2018).
- de Vries, F. T., et al. Soil bacterial networks are less stable under drought than fungal networks. Nature Communications. 9 (1), 3033(2018).
- Colin, C., et al. Correlating transcriptional networks to breast cancer survival: a large-scale coexpression analysis. Carcinogenesis. (10), 2300-2308 (2013).
- Ma, B., Zhao, K., Lv, X., et al. Genetic correlation network prediction of forest soil microbial functional organization. ISME J. 12, 2492-2505 (2018).
- Zhang, B., Horvath, S. A general framework for weighted gene co-expression network analysis. Statistical applications in genetics and molecular biology. 4 (1), (2005).
- Edwards, J. A., et al. Compositional shifts in root-associated bacterial and archaeal microbiota track the plant life cycle in field-grown rice. PLoS Biology. 16 (2), 2003862(2018).
- Bashan, Y., De-Bashan, L. E. How the Plant Growth-Promoting Bacterium Azospirillum Promotes Plant Growth-A Critical Assessment. Advances in Agronomy. 108, 77-136 (2010).
- Lovley, D. R., et al. Geobacter: The Microbe Electric's Physiology, Ecology, and Practical Applications. Advances in Microbial Physiology. 59, 1(2011).
- Langfelder, P., Horvath, S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 9, 559(2008).
- Zhang, B., Horvath, S. A General Framework for Weighted Gene Co-expression Network Analysis. Statistical Applications in Genetics and Molecular Biology. 4 (1), 17(2005).
- Horvath, S., Dong, J. Geometric interpretation of Gene Co-expression Network Analysis. PLoS Computational Biology. 4 (8), 1000117(2008).
- Langfelder, P., Horvath, S. Eigengene networks for studying the relationships between co-expression modules. BMC Systems Biology. 1, 54(2007).
- Horvath, S., et al. Analysis of Oncogenic Signaling Networks in Glioblastoma Identifies ASPM as a Novel Molecular Target. Proceedings of the National Academy of Sciences of the United States of America. 103 (46), 17402-17407 (2006).
- Carlson, M. R., et al. and Sequence Conservation: Predictions from Modular Yeast Co-expression Networks. BMC Genomics. 7 (1), 40(2006).
- Fuller, T., et al. Weighted Gene Co-expression Network Analysis Strategies Applied to Mouse Weight. Mammalian Genome. 6 (18), 463-472 (2007).
- Yin, L., Wang, Y., Lin, Y., et al. Explorative analysis of the gene expression profile during liver regeneration of mouse: a microarray-based study[J]. Artificial Cells Nanomedicine & Biotechnology. 47 (1), 1113-1121 (2019).
- Oldham, M., Horvath, S., Geschwind, D. Conservation and Evolution of Gene Co-expression Networks in Human and Chimpanzee Brains. Proceedings of the National Academy of Sciences of the United States of America. 103 (47), 17973-17978 (2006).
- Oldham, M. C., et al. Functional organization of the transcriptome in human brain. Nature Neuroscience. 11 (11), 1271-1282 (2008).
- Keller, M. P., et al. A gene expression network model of type 2 diabetes links cell cycle regulation in islets with diabetes susceptibility. Genome Research. 18 (5), 706-716 (2008).
- Weston, D., Gunter, L., Rogers, A., Wullschleger, S. Connecting genes, coexpression modules, and molecular signatures to environmental stress phenotypes in plants. BMC Systems Biology. 2 (1), 16(2008).
- Jorda´n, F. Keystone species and food webs. Biological Sciences. 364, 1733-1741 (2009).
- Backhed, F., Ley, R. E., Sonnenburg, J. L., Peterson, D. A., Gordon, J. I. Host-Bacterial Mutualism in the Human Intestine. Science. 307, 1915-1920 (2009).
访问受限。请登录或开始试用以查看此内容。
重印与许可
