SparCC的微生物网络构建示例
续前文“微生物共发生网络”，本篇继续简介SparCC的网络构建方法。
基于高通量测序的技术，例如16S rRNA分析，为阐明天然微生物群落的复杂结构提供了便利。对于识别群落中微生物相互作用，基于物种丰度组成数据的相关性分析是一种常见方法，但是这种分析可能会产生与真实情况相违背的虚假关系，归因于测序获得的物种丰度或基因丰度等很难绝对定量。群落多样性是调节这种组合效应的剧烈程度的关键因素，为此SparCC方法被开发出来，它能够根据成分数据估算相关值（Friedman and Alm, 2012）。已经表明，SparCC是一种适合测序数据特征的新颖方法，可以推断物种或基因之间的相关性，以及构建微生物物种或基因功能相互作用网络。
接下来展示SparCC工具计算相关性的方法，并通过相关矩阵构建关联网络。
SparCC安装
在该链接中访问SparCC资源，下载程序及查看操作文档：
https://bitbucket.org/yonatanf/sparcc/src/default/
在Linux平台下，命令行操作如下。
#克隆库
hg clone https://bitbucket.org/yonatanf/sparcc
#添加环境变量（例如我将 sparcc 存放至 /home/ly/software/sparcc/）
export PATH="/home/ly/software/sparcc/:$PATH"
#添加可执行权限
chmod -R 755 /home/ly/software/sparcc/
#sparcc 需要 python2 环境支持（python3 不行）
#如果出来帮助选项就代表成功了
#若提示缺少 python 模块（numpy、pandas），额外安装下即可
SparCC.py -h
SparCC构建微生物共发生网络示例
接下来展示SparCC的使用，以“微生物共发生网络”为例。
备注：SparCC安装路径中的“example”，提供了示例数据和示例结果文件，可帮助了解。这些文件的具体操作细节来自“readme.rst”中的示例流程。
我们换个数据走一下。
首先准备一个物种丰度表，以OTU丰度表为例，其中每一列代表一个样本，每一行代表一种OTU，交叉区域为各OTU在各样本中的丰度。
备注：该OTU表可以提前作些预处理，例如剔除低丰度、低频物种，或者执行有效的丰度预转化等。
1、计算相关矩阵
如下示例，默认使用SparCC相关性计算OTU间的关联程度，SparCC将根据观测值的Dirichlet分布对真实得分进行估计，并对5次估计取平均获得观测得分。若期望使用其它相关系数，可通过-a或--algo参数指定。
#第 1 步，计算观测值的相关矩阵
#SparCC.py -h
SparCC.py otu_table.txt -i 5 --cor_file=cor_sparcc.out.txt > sparcc.log
对于输出的矩阵“cor_sparcc.out.txt”中的数值，可将它理解为数据集中各OTU两两之间的关联程度，正、负值分别表示了正、负关联（丰度的相同或相反趋势改变），值越大代表关联强度越高。
接下来，就需要评估这种关联程度是否是显著（有效）的。
2、获得自举分布
基于重抽样的随机替换为评估显著性提供了方法。SparCC通过bootstrap方法对原数据集重抽样，获得随机（替换）的数据集。在后续，将通过这些随机数据集计算伪p值，以用于评估初始观测得分的显著性。
如下示例，将生成100个重抽样后的随机数据集。
#第 2 步，通过自举抽样获得随机数据集
#MakeBootstraps.py -h
MakeBootstraps.py otu_table.txt -n 100 -t bootstrap_#.txt -p pvals/ >> sparcc.log
注：SparCC的文档中建议不少于100次抽样，以在后续获得足够可信的伪p值。
3、计算伪p值
在获得了自举分布（多次重抽样的随机数据集）后，计算这些随机数据集中OTU的相关程度，即获得随机值的相关矩阵。随后，通过比较观测值的相关矩阵中数值在随机值的相关矩阵中的分布，从中生成伪p值，由于相关性有正也有负，需计算双侧区间。
#第 3 步，计算伪 p 值，作为评估相关性显著的依据
#首先通过循环语句批处理，获得各随机数据集中变量的相关矩阵（随机值的相关矩阵）
for n in {0..100}; do SparCC.py pvals/bootstrap_${n}.txt -i 5 --cor_file=pvals/bootstrap_cor_${n}.txt >> sparcc.log; done
#通过观测值的相关矩阵中系数（cor0），以及随机值的相关矩阵中系数（corN），考虑 |cor0|>|corN| 的频率，获得伪 p 值（我猜的应该是这样......）
#PseudoPvals.py -h
PseudoPvals.py cor_sparcc.out.txt pvals/bootstrap_cor_#.txt 100 -o pvals/pvals.two_sided.txt -t two_sided >> sparcc.log
最后获得的矩阵“pvals.two_sided.txt”，为伪p值矩阵，代表了两两OTU之间相关程度的显著性，值越低越显著。该矩阵与初始获得的观测值的相关矩阵“cor_sparcc.out.txt”对应。
随后，同时考虑相关性的强度以及显著水平，对相关矩阵中的数值作下筛选。
4、合并相关矩阵和p值矩阵，获得最终网络
不妨考虑使用R进行矩阵操作，根据相关性的强度以及显著水平自定义筛选，只保留具有显著的强相关关系，如下示例，最终获得邻接矩阵类型的网络文件。
#观测值的相关矩阵
cor_sparcc <- read.delim('cor_sparcc.out.txt', row.names = 1, sep = '\t', check.names = FALSE)
#伪 p 值矩阵
pvals <- read.delim('pvals.two_sided.txt', row.names = 1, sep = '\t', check.names = FALSE)
#保留 |相关性|≥0.8 且 p<0.01的值
cor_sparcc[abs(cor_sparcc) < 0.8] <- 0
pvals[pvals>=0.01] <- -1
pvals[pvals<0.01 & pvals>=0] <- 1
pvals[pvals==-1] <- 0
#筛选后的邻接矩阵
adj <- as.matrix(cor_sparcc) * as.matrix(pvals)
diag(adj) <- 0 #将相关矩阵中对角线中的值（代表了自相关）转为 0
write.table(data.frame(adj, check.names = FALSE), 'neetwork.adj.txt', col.names = NA, sep = '\t', quote = FALSE)
图示邻接矩阵，不存在相关就是0，存在相关就是非0的数值，正值表示正相关，负值表示负相关，数值的绝对值大小代表相关强度。
R包igraph的网络操作
随后，不妨继续使用R，通过邻接矩阵构建网络，并对网络格式进行转换，以便能够使用更多工具（如Cytoscape、Gephi等）进行统计分析、可视化操作等。
igraph包提供了灵活的网络操作方法，首先通过它转换网络格式。
##网络格式转换
library(igraph)
#输入数据，邻接矩阵
neetwork_adj <- read.delim('neetwork.adj.txt', row.names = 1, sep = '\t', check.names = FALSE)
head(neetwork_adj)[1:6] #邻接矩阵类型的网络文件
#邻接矩阵 -> igraph 的邻接列表，获得含权的无向网络
g <- graph_from_adjacency_matrix(as.matrix(neetwork_adj), mode = 'undirected', weighted = TRUE, diag = FALSE)
g #igraph 的邻接列表
#这种转换模式下，默认的边权重代表了 sparcc 计算的相关性（存在负值）
#由于边权重通常为正值，因此最好取个绝对值，相关性重新复制一列作为记录
E(g)$sparcc <- E(g)$weight
E(g)$weight <- abs(E(g)$weight)
#再转为其它类型的网络文件，例如
#再由 igraph 的邻接列表转换回邻接矩阵
adj_matrix <- as.matrix(get.adjacency(g, attr = 'sparcc'))
write.table(data.frame(adj_matrix, check.names = FALSE), 'network.adj_matrix.txt', col.names = NA, sep = '\t', quote = FALSE)
#graphml 格式，可使用 gephi 软件打开并进行可视化编辑
write.graph(g, 'network.graphml', format = 'graphml')
#gml 格式，可使用 cytoscape 软件打开并进行可视化编辑
write.graph(g, 'network.gml', format = 'gml')
#边列表，也可以直接导入至 gephi 或 cytoscape 等网络可视化软件中进行编辑
edge <- data.frame(as_edgelist(g))
edge_list <- data.frame(
source = edge[[1]],
target = edge[[2]],
weight = E(g)$weight,
sparcc = E(g)$sparcc
)
head(edge_list)
write.table(edge_list, 'network.edge_list.txt', sep = '\t', row.names = FALSE, quote = FALSE)
#节点属性列表，对应边列表，记录节点属性，例如
node_list <- data.frame(
nodes_id = V(g)$name, #节点名称
degree = degree(g) #节点度
)
head(node_list)
write.table(node_list, 'network.node_list.txt', sep = '\t', row.names = FALSE, quote = FALSE)
随后，可继续在R中编辑网络。
例如，先前的博文简介过在R中进行网络拓扑特征计算的方法，如有需要可点击查看：
或者使用其它图形界面的工具处理更方便。
例如后续使用Gephi打开上述转化后的文件“network.graphml”，计算拓扑属性，可视化等。
对于Gephi软件，或者Cytoscape等，还请自行了解这些软件的使用了。
参考文献
Friedman J, Alm E J. Inferring Correlation Networks from Genomic Survey Data. PLOS Computational Biology, 2012, 8(9).
猜你喜欢
10000+：菌群分析 宝宝与猫狗 梅毒狂想曲 提DNA发Nature Cell专刊 肠道指挥大脑
文献阅读 热心肠 SemanticScholar Geenmedical
16S功能预测 PICRUSt FAPROTAX Bugbase Tax4Fun
生物科普: 肠道细菌 人体上的生命 生命大跃进 细胞暗战 人体奥秘
写在后面
为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。PI请明示身份，另有海内外微生物相关PI群供大佬合作交流。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍未解决群内讨论，问题不私聊，帮助同行。
学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”
点击阅读原文，跳转最新文章目录阅读
扫一扫
1491
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
