扩增子分析|R分析之微生物生态网络稳定性评估之节点和连接的恒常性、节点持久性以及组成稳定性指数计算

一、引言

周集中老师团队于2021年在Nature climate change发表的文章，阐述了网络稳定性评估的原理算法，并提供了完整的代码。自此对微生物生态网络的评估具有更全面的指标，自此网络稳定性的评估广受大家欢迎。本文将介绍网络稳定性之节点和连接的恒常性、节点持久性以及组成稳定性的原理，计算方法以及代码，如有疑问欢迎讨论交流。

为了理解增温是否以及如何影响构建的微生物生态网络（MEN）的稳定性，我们使用了多种指标来表征网络及其嵌入成员的稳定性，包括鲁棒性、易损性、节点和连接的恒常性、节点持久性以及组成稳定性。我们通过消除模拟和经验观察评估了网络及网络内微生物群落的稳定性。上一节我们讨论了基于模拟的方式评估网络稳定性的稳健性（Robustness）和易损性指标。感兴趣的可以跳转至下述链接阅读：

扩增子分析|微生物生态网络稳定性评估之鲁棒性（Robustness）和易损性（Vulnerability）

本节我们将讨论基于实际观测值评估网络稳定性的节点和连接的恒常性、节点持久性以及组成稳定性指数。

二、基础知识—基于实际观测值的网络稳定性分析

2.1 节点恒常性（Node constancy）

恒常性衡量物种的时间稳定性，定义为μ/σ，其中μ是随时间变化的丰度平均值，σ是标准差。

式中，μi是平均值；

σi是不同时间点网络中物种i丰度的标准差。

仅当物种i在某一时间点处于 MEN 中时，物种i在某一时间点的丰度才为正。否则，物种i在该时间点的丰度被视为零。当σi= 0 时，恒常性不是有限的，因此将从后续分析中移除。最后，所有恒常性值的平均值被报告为不同条件下网络节点的恒常性。使用 Mann-Whitney U检验测试处理之间链接恒定性的差异。

然后，参考Hui等人提出的方法，计算多个网络间的重叠节点数。被包含的时间点称为“顺序”。更高的重叠节点数量表示网络中物种组成的周转速度较慢。

2.2 连接恒常性（Link constancy）

类似于节点恒常性，计算连接恒常性。如果节点i和j在网络中呈正相关，则令= 1；如果节点i和j在网络中呈负相关，则令 = 1。如果i和j 之间没有连接，则令= = 0。连接恒常性（或）的计算方法如下：

式中，

和分别为来自不同时间点所有网络中的均值和标准差，和分别为来自不同时间点所有网络中的均值和标准差。对于= 0或= 0的情况，从后续分析中删除。最后，所有连接恒常性的平均值作为网络连接的恒常性。

2.3 节点持久性（Node persistence）

节点持久性定义为生态系统中共存物种占物种总数的比例。因此，节点的持久性计算为连续多年存在于网络中节点的百分比，计算方法是：

式中，v是在多个连续时间点采集的样本数；

S是网络中的 OTU 总数；

是Dirac delta函数，如果样本i中OTU k的丰度大于 0，则= 1，如果样本i中不存在OUT k ，则 = 0。

2.4 组成稳定性（Compositional stability）

组成稳定性评估群落结构随时间的变化。我们使用样本×OTU矩阵计算了网络中微生物群落的组成稳定性，如下所示

式中，v为在多个连续时间点从同区域采集的样本数。S为网络中的 OTU 总数。为样本i中 OUT k的丰度。若群落结构没有变化，即，稳定性指数等于0。

当v= 2 时，

其中和分别是第i年和第i– 1 年样本中OUT k的丰度。

图中，d,每两年网络化社区随时间的变化组成稳定性。e ，组成稳定性和节点持久性的皮尔逊相关性。该线表示在变暖（红色）或控制（蓝色）条件下的最小二乘最佳拟合。f ，变暖组（红色框）或控制组（蓝色框）条件下网络复杂性和稳定性指数之间的 Pearson 相关性。（网络复杂性可根据相关代码计算）这里显示显著（P ≤ 0.05 ）相关性，橙色表示正相关，绿色表示负相关。单元格内的数字是相应的相关系数。P > 0.05 的相关性标记为灰色。

三、示例数据和R代码

本文的示例数据和代码均来自于2021年周集中老师团队的Nature climate change文章，感兴趣的同学可以自行去学习。

3.1 组成稳定性指标

🌟准备数据

准备otu.csv表格，该表格为要计算鲁棒性的网络otu数据表。

代码每次计算一个网络的稳健性，因此需要计算几个网络就运行几次代码，每次将输入文件名修改。

🌟完整代码

# 加载自定义的 R 脚本，其中可能包含必要的自定义函数
source("ieggrtools.r")# 定义社区数据和处理数据的文件路径
com.file = "input_file/OTUtable_NetworkedOTUs_AllSamples.txt"
treat.file = "input_file/SampleMap_AllSamples.txt"# 从文件中加载社区数据，并进行转置
comm = t(lazyopen(com.file))
# 加载处理数据
treat = lazyopen(treat.file)# 检查社区数据中是否有缺失值
sum(is.na(comm))
# 将所有的 NA 值替换为 0
comm[is.na(comm)] = 0# 匹配社区数据和处理数据的样本名称
sampc = match.name(rn.list = list(comm = comm, treat = treat))
comm = sampc$comm
treat = sampc$treat# 查看处理数据的前几行，确认哪个列是 Plot_ID，哪个列是 Year
head(treat)# 获取唯一的 Plot_ID 和按顺序排列的 Year
plot.lev = unique(treat$Plot_ID)
year.lev = sort(unique(treat$Year))# 生成一个序列，用于计算多重顺序稳定性（Zeta 多样性水平）
zeta.lev = 2:length(year.lev)# 创建一个包含不同年份窗口的列表，用于计算 Zeta 多样性
year.windows = lapply(1:length(zeta.lev), function(i) {zetai = zeta.lev[i]lapply(1:(length(year.lev) - zetai + 1), function(j) {year.lev[j:(j + zetai - 1)]})
})
names(year.windows) = zeta.lev
year.windows  # 显示生成的年份窗口# 定义一个函数，用于计算一组样本的社区稳定性
comstab <- function(subcom) {((nrow(subcom) * sum(apply(subcom, 2, min))) / sum(subcom))^0.5
}# 计算每个 Plot 在不同年份窗口内的稳定性
stabl = lapply(1:length(year.windows), function(i) {stabi = t(sapply(1:length(plot.lev), function(j) {plotj = plot.lev[j]sapply(1:length(year.windows[[i]]), function(k) {yearwdk = year.windows[[i]][[k]]# 找到满足条件的样本sampijk = rownames(treat)[which((treat$Plot_ID == plotj) & (treat$Year %in% yearwdk))]outijk = NA# 检查样本数量是否匹配年份窗口if (length(sampijk) < length(yearwdk)) {warning("plot ", plotj, " has missing year in year window ", paste(yearwdk, collapse = ","))} else if (length(sampijk) > length(yearwdk)) {warning("plot ", plotj, " has duplicate samples in at least one year of window ", paste(yearwdk, collapse = ","))} else {# 计算样本的社区稳定性comijk = comm[which(rownames(comm) %in% sampijk), , drop = FALSE]outijk = comstab(comijk)}outijk})}))# 确保行数和 Plot_ID 的数量匹配if (nrow(stabi) != length(plot.lev) & nrow(stabi) == 1) {stabi = t(stabi)}rownames(stabi) = plot.levcolnames(stabi) = sapply(year.windows[[i]], function(v) {paste0("Zeta", zeta.lev[i], paste0(v, collapse = ""))})stabi
})
# 将所有的稳定性数据合并成一个矩阵
stabm = Reduce(cbind, stabl)# 移除特定年份 (Y09) 并注释每个 Plot 的处理信息
treat.rm09 = treat[which(treat$Year != "Y09"), , drop = FALSE]
output = data.frame(treat.rm09[match(rownames(stabm), treat.rm09$Plot_ID), 1:5, drop = FALSE], stabm, stringsAsFactors = FALSE)
output = output[order(output$Warming, output$Clipping, output$Precip, output$Plot_ID), , drop = FALSE]
head(output)  # 显示输出的前几行# 加载 tidyr 包，用于数据转换
library(tidyr)# 将稳定性数据转换为长格式，方便绘图
cs.mo = output
cs.long = gather(cs.mo, year, cs, Zeta2Y09Y10:Zeta6Y09Y10Y11Y12Y13Y14, factor_key = TRUE)
cs.long = cs.long[-which(is.na(cs.long$cs)), ]
cs.long$EndYear = as.numeric(gsub(".*Y", "", cs.long$year))# 定义一个函数，用于生成线性模型文本
get_text <- function(y, x) {lm = lm(y ~ x)lm_p = round(summary(lm)$coefficients[2, 4], 3)lm_r = round(summary(lm)$adj.r.squared, 3)lm_s = round(summary(lm)$coefficients[2, 1], 3)txt = paste(lm_s, " (", lm_r, ", ", lm_p, ")", sep = "")return(txt)
}# 定义一个函数，用于绘制完整的社区稳定性图
plot_cs_full <- function(xc, yc, xw, yw) {plot(yw ~ xw, ylim = c(0.2, 1), xlim = c(10, 14), ylab = "Compositional stability", xlab = "Year", pch = 19, cex.axis = 0.7, cex.lab = 0.7, tck = -0.03, col = "#e7211f")abline(lm(yw ~ xw), col = "#e7211f")points(yc ~ xc, col = "#214da0")abline(lm(yc ~ xc), col = "#214da0")txt1 = get_text(yw, xw)txt2 = get_text(yc, xc)mtext(txt1, side = 3, line = 0.8, cex = 0.5, col = "#e7211f")mtext(txt2, side = 3, line = 0.1, cex = 0.5, col = "#214da0")
}# 定义一个函数，用于绘制 Zeta6 稳定性图，并进行 Wilcoxon 检验
plot_cs_zeta6 <- function(xc, yc, xw, yw) {plot(yw ~ xw, ylim = c(0.2, 1), xlim = c(10, 14), ylab = "Compositional stability", xlab = "Year", pch = 19, cex.axis = 0.7, cex.lab = 0.7, tck = -0.03, col = "#e7211f")points(yc ~ xc, col = "#214da0")wt = wilcox.test(yw, yc)mtext(paste(wt$statistic, wt$p.value, sep = ","), cex = 0.5)
}# 绘制 Figure 3d：相邻两年的组分稳定性
# 设置颜色和符号
colors = c(rep("#214da0", 5), rep("#e7211f", 5))
syms = c(1, 1, 1, 1, 1, 19, 19, 19, 19, 19)
# 定义一个函数 plot_cs，用于绘制稳定性图，带误差条
plot_cs <- function(yw, yc, ebw, ebc) {x = c(1, 2, 3, 4, 5)yl = c(min(c(yw - ebw, yc - ebc)), max(c(yw + ebw, yc + ebc)))# 绘制加温处理 (W) 的数据及回归线plot(yw ~ x, ylim = yl, ylab = "Compositional stability", xlab = "Year of warming", pch = 19, cex.axis = 0.7, cex.lab = 0.7, tck = -0.03, col = "red")abline(lm(yw ~ x), col = "red")arrows(x, yw - ebw, x, yw + ebw, length = 0.05, angle = 90, code = 3, col = "red")# 绘制无加温处理 (N) 的数据及回归线points(yc ~ x, col = "blue")abline(lm(yc ~ x), lty = 2, col = "blue")arrows(x, yc - ebc, x, yc + ebc, length = 0.05, angle = 90, code = 3, col = "blue")# 获取线性模型文本txt1 = get_text(yw, x)txt2 = get_text(yc, x)# 添加文本到图表中mtext(txt1, side = 3, line = 0.8, cex = 0.5, col = "red")mtext(txt2, side = 3, line = 0.1, cex = 0.5, col = "blue")
}# 筛选 Zeta2 数据并计算均值和标准差
zeta2 = cs.long[grep("Zeta2", cs.long$year), ]
zeta2_means = aggregate(zeta2$cs, list(zeta2$Warming, zeta2$EndYear), FUN = mean)
zeta2_sds = aggregate(zeta2$cs, list(zeta2$Warming, zeta2$EndYear), FUN = sd)# 绘制相邻年份的组分稳定性图
plot_cs(yw = zeta2_means$x[which(zeta2_means$Group.1 == "W")],yc = zeta2_means$x[which(zeta2_means$Group.1 == "N")],ebw = zeta2_sds$x[which(zeta2_sds$Group.1 == "W")],ebc = zeta2_sds$x[which(zeta2_sds$Group.1 == "N")])#################即可完成组合稳定性随时间变化图#####################

🌟输出结果

图片解读，红色是变温，蓝色是对照。数字0.071是线性拟合斜率，0.862和0.015是调整后的r2和P值。

3.2 节点恒常性指标

🌟准备数据

准备OTUtable_NetworkedOTUs_AllSamples.csv表格，该表格为要节点恒常性的otu数据表

准备SampleMap_AllSamples.csv表格，样本基本映射文件

🌟完整代码

# 清理工作环境中的所有对象
rm(list = ls())# 加载所需的 R 包
library(ggplot2)  # 用于数据可视化
library(gridExtra)  # 用于排列多个图表# 读取 OTU 表和样本映射文件
otu <- read.csv("OTUtable_NetworkedOTUs_AllSamples.csv",  header=T, row.names=1)
map <- read.csv("SampleMap_AllSamples.csv",  header=T)# 定义每年样本的索引
id_09 = rep(which(map$Year == "Y09"), each=2)  # 重复每个索引两次
id_10 = which(map$Year == "Y10")
id_11 = which(map$Year == "Y11")
id_12 = which(map$Year == "Y12")
id_13 = which(map$Year == "Y13")
id_14 = which(map$Year == "Y14")# 检查样本的 plot 顺序
data.frame(map$Plot_full_name[id_09],map$Plot_full_name[id_10],map$Plot_full_name[id_11],map$Plot_full_name[id_12],map$Plot_full_name[id_13],map$Plot_full_name[id_14])# 获取第 14 年的温度处理（是否加温）
warming_trt = map$Warming[id_14]# 将 OTU 表按年份拆分为不同数据框
otu_09 = as.data.frame(otu[, id_09])
otu_10 = as.data.frame(otu[, id_10])
otu_11 = as.data.frame(otu[, id_11])
otu_12 = as.data.frame(otu[, id_12])
otu_13 = as.data.frame(otu[, id_13])
otu_14 = as.data.frame(otu[, id_14])# 检查 OTU 表的样本顺序是否与映射文件一致
sum(names(otu_09) != map$Sample[id_09])  # 检查列名是否与样本名称匹配
sum(names(otu_10) != map$Sample[id_10])
sum(names(otu_11) != map$Sample[id_11])
sum(names(otu_12) != map$Sample[id_12])
sum(names(otu_13) != map$Sample[id_13])
sum(names(otu_14) != map$Sample[id_14])# 计算节点稳定性（constancy）
# 计算每个 OTU 的平均值
otu_mean = (otu_09 + otu_10 + otu_11 + otu_12 + otu_13 + otu_14) / 6
# 计算标准差
otu_sd = sqrt(((otu_09 - otu_mean)^2 + (otu_10 - otu_mean)^2 + (otu_11 - otu_mean)^2 +(otu_12 - otu_mean)^2 + (otu_13 - otu_mean)^2 + (otu_14 - otu_mean)^2) / 5)
# 计算节点稳定性为平均值除以标准差
otu_constancy = otu_mean / otu_sd# 将节点稳定性结果保存为 CSV 文件
write.table(otu_constancy, "observed_constancy_of_each_node.csv", sep=",")# 根据温度处理分割节点稳定性数据
otu_constancy_w = otu_constancy[, which(warming_trt == "W")]  # 加温处理
otu_constancy_c = otu_constancy[, which(warming_trt == "N")]  # 控制处理# 计算每个 OTU 在加温和控制处理下的平均节点稳定性
otu_constancy_w_avg = rowMeans(otu_constancy_w)
otu_constancy_c_avg = rowMeans(otu_constancy_c)# 提取有限值的节点稳定性
con_w = otu_constancy_w_avg[is.finite(otu_constancy_w_avg)]
con_c = otu_constancy_c_avg[is.finite(otu_constancy_c_avg)]# 创建一个数据框，用于绘图
nc_df = data.frame(warming = c(rep("control", length(con_c)), rep("warming", length(con_w))),nc = c(con_c, con_w))# 使用 ggplot2 绘制箱线图和散点图
ggplot(nc_df, aes(x = warming, y = nc, fill = warming)) +geom_boxplot(alpha = 1, width = 0.4, outlier.shape = NA) +  # 绘制箱线图，不显示离群点geom_jitter(shape = 16, size = 0.5, position = position_jitterdodge(jitter.width = 0.01, dodge.width = 0.8)) +  # 绘制抖动的散点图scale_fill_manual(values = c("#214da0", "#e7211f")) +  # 设置填充颜色labs(y = "Node constancy") +  # 设置 y 轴标签theme(axis.line = element_line(colour = "black"),  # 轴线颜色为黑色panel.grid.major = element_blank(),  # 移除主要网格线panel.grid.minor = element_blank(),  # 移除次要网格线panel.background = element_blank(),  # 移除背景panel.border = element_rect(colour = "black", fill = NA, size = 0.6),  # 设置边框legend.position = "none"  # 不显示图例)

🌟输出结果

3.3 连接恒常性指标

🌟准备数据

准备NetworkEdgeTable_AllNetworks.csv表格，即边属性表，该表可由网络构建上游分析获得。

🌟完整代码

# 清理工作环境中的所有对象
rm(list = ls())
# 读取网络连接表
edge_prop = read.csv("NetworkEdgeTable_AllNetworks.csv", header=T)
edge_prop$edge_sign = paste(edge_prop$edge, edge_prop$V5, sep="_")  # 创建带符号的边缘标识
# 计算唯一连接的总数和带符号的总数
total_e = length(unique(edge_prop$edge))
total_e_sign = length(unique(edge_prop$edge_sign))  # 注意：符号是一致的
# 创建连接表
edge_table = data.frame(matrix(NA, nrow=total_e, ncol=11))
row.names(edge_table) = unique(edge_prop$edge)
names(edge_table) = unique(edge_prop$year_trt)
# 填充连接表
for (i in 1:ncol(edge_table)) {edge_prop_year = edge_prop[which(edge_prop$year_trt == names(edge_table[i])), ]edge_table[which(rownames(edge_table) %in% edge_prop_year$edge), i] = 1
}
edge_table[is.na(edge_table)] <- 0  # 将 NA 替换为 0# 查看连接表的前几行
head(edge_table)# 计算平均值和标准差
edge_table$avg_c = rowMeans(edge_table[, c(1, 2, 4, 6, 8, 10)])
edge_table$sd_c = apply(edge_table[, c(1, 2, 4, 6, 8, 10)], 1, FUN=sd)
edge_table$avg_w = rowMeans(edge_table[, c(1, 3, 5, 7, 9, 11)])
edge_table$sd_w = apply(edge_table[, c(1, 3, 5, 7, 9, 11)], 1, FUN=sd)# 计算恒常性
edge_table$constancy_c = edge_table$avg_c / edge_table$sd_c
edge_table$constancy_w = edge_table$avg_w / edge_table$sd_w# 处理无穷大值
edge_table$constancy_c[is.infinite(edge_table$constancy_c)] = NA
edge_table$constancy_w[is.infinite(edge_table$constancy_w)] = NA# edge_table[1:100,]# 计算链接恒常性的相关统计
v_constancy_c = edge_table$constancy_c[which(!is.na(edge_table$constancy_c))]
n_constancy_c = sum(!is.na(edge_table$constancy_c))
avg_constancy_c = mean(edge_table$constancy_c, na.rm=T)
sd_constancy_c = sd(edge_table$constancy_c, na.rm=T)v_constancy_w = edge_table$constancy_w[which(!is.na(edge_table$constancy_w))]
n_constancy_w = sum(!is.na(edge_table$constancy_w))
avg_constancy_w = mean(edge_table$constancy_w, na.rm=T)
sd_constancy_w = sd(edge_table$constancy_w, na.rm=T)# 进行 t 检验
t.test(v_constancy_w, v_constancy_c)library(ggplot2)# 绘制扩展数据图 5d - 连接恒常性
df_lc_uw = data.frame(Treatment = c("control", "warming"), Link_constancy = c(avg_constancy_c, avg_constancy_w), se = c(sd_constancy_c / sqrt(n_constancy_c), sd_constancy_w / sqrt(n_constancy_w))
)# 创建条形图
ggplot(df_lc_uw, aes(x=Treatment, y= Link_constancy)) +geom_bar(stat="identity", fill=c("#214da0", "#e7211f"), width=0.5) +geom_errorbar(aes(ymin=Link_constancy - se, ymax=Link_constancy + se), width=0.2) +labs(x="Treatment", y="Unweighted link constancy") +  # 添加中文标签theme(axis.line = element_line(colour="black"),panel.grid.major = element_blank(),panel.grid.minor = element_blank(),panel.background = element_blank(),panel.border = element_rect(colour="black", fill=NA, size=0.6),legend.position="none")  # 隐藏图例

🌟输出结果

3.4 节点持久性指标

🌟准备数据

准备OTUtable_NetworkedOTUs_AllSamples.csv表格，该表格为要节点恒常性的otu数据表。

准备SampleMap_AllSamples.csv表格，样本基本映射文件

🌟完整代码

# 清理工作环境中的所有对象
rm(list = ls())# 读取 OTU 表和样本映射文件
otutab <- read.csv("OTUtable_NetworkedOTUs_AllSamples.csv",  header=T, row.names=1)
treatused <- read.csv("SampleMap_AllSamples.csv",  header=T)### 计算每个样地的持久性
## 这里持久性 = 在特定年份中物种的比例
comm <- t(otutab)  # 转置OTU表
sum(row.names(comm) == row.names(treatused))  # 确保样本名称匹配fake09.N <- comm[1:24, ]
row.names(fake09.N) <- gsub("S", "N", row.names(fake09.N))  # 替换样本名称
comm2 <- rbind(fake09.N, comm)  # 合并数据fake09.treat <- treatused[1:24, ]
row.names(fake09.treat) <- gsub("S", "N", row.names(fake09.treat))  # 替换样本名称
fake09.treat$Plot_full_name <- gsub("S", "N", fake09.treat$Plot_full_name)  # 替换样地名称
fake09.treat$Plot_ID <- gsub("S", "N", fake09.treat$Plot_ID)
treat2 <- rbind(fake09.treat, treatused)  # 合并处理数据
treat2$Plot_ID = factor(treat2$Plot_ID)  # 转换为因子# 组织09年的升温处理
treat2$Warming[1:48] <- treat2$Warming[49:96][match(treat2$Plot_ID[1:48], treat2$Plot_ID[49:96])]# 按照样地ID分割数据
test1 <- split(data.frame(comm2), treat2$Plot_ID)# 定义年份
years = c("Y09", "Y10", "Y11", "Y12", "Y13", "Y14")
l <- rep(list(0:1), 6)
allcombines <- expand.grid(l)  # 生成所有组合
allcombines <- allcombines[c(4, 7, 13, 25, 49, 8, 15, 29, 57, 16, 31, 61, 32, 63, 64), ]  # 仅选择相邻年份
allcombines <- t(allcombines)# 计算持久性
test <- sapply(test1, function(x) {test3 <- sapply(1:ncol(allcombines), function(i) {coms <- allcombines[, i] * xtotal.sp <- sum(colSums(coms > 0) > 0)  # 计算总物种数persist.sp <- sum(colSums(coms > 0) == sum(allcombines[, i]))  # 计算持久物种数prop = persist.sp / total.sp  # 计算比例names(prop) = paste0(years[allcombines[, i] > 0], collapse = "")prop})test3
})# 整理结果
persist.result <- t(test)
output = data.frame(treat2[match(rownames(persist.result), treat2$Plot_ID), 1:5, drop=FALSE], persist.result, stringsAsFactors = FALSE)
names(output) <- paste(c(rep("", 5),rep("Zeta2", 5),rep("Zeta3", 4),rep("Zeta4", 3),rep("Zeta5", 2),rep("Zeta6", 1)), names(output), sep="")
head(output)
#输出结果文件
write.csv(output, "MultiOrderPersistence.csv")library(tidyr)
## 绘制图 S7fghij - 多级持久性
op.mo = outputop.long <- gather(op.mo, year, op, Zeta2Y09Y10:Zeta6Y09Y10Y11Y12Y13Y14, factor_key=TRUE)  # 转换为长格式
op.long$EndYear = as.numeric(gsub(".*Y", "", op.long$year))  # 提取结束年份# 定义文本提取函数
get_text <- function(y, x) {lm = lm(y ~ x) lm_p = round(summary(lm)$coefficients[2, 4], 3)lm_r = round(summary(lm)$adj.r.squared, 3)lm_s = round(summary(lm)$coefficients[2, 1], 3)	txt = paste(lm_s, " (", lm_r, ", ", lm_p, ")", sep="")	return(txt)
}# 绘制持久性图
plot_op_full <- function(xc, yc, xw, yw) {plot(yw ~ xw, ylim=c(0, 0.3), xlim=c(10, 14), ylab="Node persistence", xlab="Year", pch=19, cex.axis=0.7, cex.lab=0.7, tck=-0.03, col="#e7211f")abline(lm(yw ~ xw), col="#e7211f")  # 添加回归线points(yc ~ xc, col="#214da0")  # 绘制控制组abline(lm(yc ~ xc), col="#214da0")  # 添加回归线txt1 = get_text(yw, xw)txt2 = get_text(yc, xc)mtext(txt1, side=3, line=0.8, cex=0.5, col="#e7211f")  # 显示文本mtext(txt2, side=3, line=0.1, cex=0.5, col="#214da0")
}# 绘制每个Zeta的图
par(mfrow=c(2, 3))
plot_op_full(xc=op.long$EndYear[intersect(which(op.long$Warming=="N"), grep("Zeta2", op.long$year))], yc=op.long$op[intersect(which(op.long$Warming=="N"), grep("Zeta2", op.long$year))], xw=op.long$EndYear[intersect(which(op.long$Warming=="W"), grep("Zeta2", op.long$year))], yw=op.long$op[intersect(which(op.long$Warming=="W"), grep("Zeta2", op.long$year))])plot_op_full(xc=op.long$EndYear[intersect(which(op.long$Warming=="N"), grep("Zeta3", op.long$year))], yc=op.long$op[intersect(which(op.long$Warming=="N"), grep("Zeta3", op.long$year))], xw=op.long$EndYear[intersect(which(op.long$Warming=="W"), grep("Zeta3", op.long$year))], yw=op.long$op[intersect(which(op.long$Warming=="W"), grep("Zeta3", op.long$year))])plot_op_full(xc=op.long$EndYear[intersect(which(op.long$Warming=="N"), grep("Zeta4", op.long$year))], yc=op.long$op[intersect(which(op.long$Warming=="N"), grep("Zeta4", op.long$year))], xw=op.long$EndYear[intersect(which(op.long$Warming=="W"), grep("Zeta4", op.long$year))], yw=op.long$op[intersect(which(op.long$Warming=="W"), grep("Zeta4", op.long$year))])plot_op_full(xc=op.long$EndYear[intersect(which(op.long$Warming=="N"), grep("Zeta5", op.long$year))], yc=op.long$op[intersect(which(op.long$Warming=="N"), grep("Zeta5", op.long$year))], xw=op.long$EndYear[intersect(which(op.long$Warming=="W"), grep("Zeta5", op.long$year))], yw=op.long$op[intersect(which(op.long$Warming=="W"), grep("Zeta5", op.long$year))])plot_op_zeta6(xc=op.long$EndYear[intersect(which(op.long$Warming=="N"), grep("Zeta6", op.long$year))], yc=op.long$op[intersect(which(op.long$Warming=="N"), grep("Zeta6", op.long$year))], xw=op.long$EndYear[intersect(which(op.long$Warming=="W"), grep("Zeta6", op.long$year))], yw=op.long$op[intersect(which(op.long$Warming=="W"), grep("Zeta6", op.long$year))])# 绘制相邻两年的观察到的持久性
colors = c(rep("#214da0", 5), rep("#e7211f", 5))
syms = c(1, 1, 1, 1, 1, 19, 19, 19, 19, 19)plot_op <- function(yw, yc, ebw, ebc) {x = c(1, 2, 3, 4, 5)yl = c(min(c(yw - ebw, yc - ebc)), max(c(yw + ebw, yc + ebc)))plot(yw ~ x, ylim=yl, ylab="Observed Node Persistence", xlab="Year of warming", pch=19, cex.axis=0.7, cex.lab=0.7, tck=-0.03, col="red")abline(lm(yw ~ x), col="red")  # 添加回归线arrows(x, yw - ebw, x, yw + ebw, length=0.05, angle=90, code=3, col="red")  # 添加误差条points(yc ~ x, col="blue")  # 绘制控制组abline(lm(yc ~ x), lty=2, col="blue")  # 添加回归线arrows(x, yc - ebc, x, yc + ebc, length=0.05, angle=90, code=3, col="blue")  # 添加误差条txt1 = get_text(yw, x)txt2 = get_text(yc, x)mtext(txt1, side=3, line=0.8, cex=0.5, col="red")  # 显示文本mtext(txt2, side=3, line=0.1, cex=0.5, col="blue")
}# 提取Zeta2的数据
zeta2 = op.long[grep("Zeta2", op.long$year), ]
zeta2_means = aggregate(zeta2$op, list(zeta2$Warming, zeta2$EndYear), FUN=mean)  # 计算均值
zeta2_sds = aggregate(zeta2$op, list(zeta2$Warming, zeta2$EndYear), FUN=sd)  # 计算标准差# 绘制Zeta2的观察到的持久性
plot_op(yw=zeta2_means$x[which(zeta2_means$Group.1 == "W")], yc=zeta2_means$x[which(zeta2_means$Group.1 == "N")], ebw=zeta2_sds$x[which(zeta2_sds$Group.1 == "W")], ebc=zeta2_sds$x[which(zeta2_sds$Group.1 == "N")])

🌟输出结果

图片解读，红色是变温，蓝色是对照。数字0.048是线性拟合斜率，0.866和0.014是调整后的r2和P值。