##欢迎批评指正！！

基于高通量测序的技术，例如16S rRNA分析，为阐明天然微生物群落的复杂结构提供了便利。对于识别群落中微生物相互作用，基于物种丰度组成数据的相关性分析是一种常见方法，但是这种分析可能会产生与真实情况相违背的虚假关系，归因于测序获得的物种丰度或基因丰度等很难绝对定量。群落多样性是调节这种组合效应的剧烈程度的关键因素，为此SparCC方法被开发出来，它能够根据成分数据估算相关值（Friedman and Alm, 2012）。已经表明，SparCC是一种适合测序数据特征的新颖方法，可以推断物种或基因之间的相关性，以及构建微生物物种或基因功能相互作用网络。

首先准备一个物种丰度表，以OTU丰度表为例，其中每一列代表一个样本，每一行代表一种OTU，交叉区域为各OTU在各样本中的丰度。

备注：该OTU表可以提前作些预处理，例如剔除低丰度、低频物种，或者执行有效的丰度预转化等。该文件无最后一列注释信息，使用时需要将物种注释结果剔除。

1、计算相关矩阵 如下示例，默认使用SparCC相关性计算OTU间的关联程度，SparCC将根据观测值的Dirichlet分布对真实得分进行估计，并对5次估计取平均获得观测得分。若期望使用其它相关系数，可通过-a或–algo参数指定。

/TJPROJ5/META_ASS/script/yuxi/work/miniconda3/bin/python\
/TJPROJ5/META_ASS/script/yuxi/Share/sparCC/SparCC3/SparCC.py example/fake_data.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 /TJPROJ5/META_ASS/script/yuxi/work/miniconda3/bin/python /TJPROJ5/META_ASS/script/yuxi/Share/sparCC/SparCC3/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
/TJPROJ5/META_ASS/script/yuxi/work/miniconda3/bin/python /TJPROJ5/META_ASS/script/yuxi/Share/sparCC/SparCC3/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,
    interaction_type = ifelse(E(g)$sparcc>0, 'positive', 'negative')
    )
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中编辑网络。

1. sh step1.sh
2. sh step2.sh
3. sh step3.sh
4. R -f 1_deal.R
5. R -f 2_deal.R