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).

猜你喜欢

系列教程：微生物组入门 Biostar 微生物组  宏基因组

专业技能：学术图表 高分文章 生信宝典 不可或缺的人

一文读懂：宏基因组 寄生虫益处 进化树

必备技能：提问 搜索  Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

扩增子分析：图表解读 分析流程 统计绘图

在线工具：16S预测培养基 生信绘图

科研经验：云笔记  云协作 公众号

编程模板: Shell  R Perl

生物科普:  肠道细菌 人体上的生命 生命大跃进  细胞暗战 人体奥秘  

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。PI请明示身份，另有海内外微生物相关PI群供大佬合作交流。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍未解决群内讨论，问题不私聊，帮助同行。

学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

点击阅读原文，跳转最新文章目录阅读