基于 Mantel 检验的网络图与相关性热力图结合的顶刊可视化图表实现

文摘   2024-10-12 18:00   中国香港  

背景

在生态学、环境科学和微生物组学等领域,复杂的环境因子与生物群落之间的相互关系是研究的重点,为了可视化这种多维度的相互作用,文献中经常使用诸如Mantel检验、相关性矩阵和网络图来揭示变量之间的联系,在R语言中,丰富的可视化库提供了强大的工具来构建这些图表,例如通过ggplot2、corrplot等包,然而,在Python中虽然缺少类似的高集成度工具,但可以通过组合多个Python库来灵活实现同样的可视化效果

文献展示的B图通过结合Mantel检验和相关性分析,构建了一个环境变量和生物群落之间的关系网络,采用不同的线条颜色和粗细来表示p值显著性,这一图形能够有效地展示环境因子与微生物组分类及功能组成之间的复杂关系,在本文中,将介绍如何利用Python中的matplotlib、seaborn和networkx等工具包,实现类似于文献中B图的复现,展示变量之间的关系

代码实现

网络图和相关矩阵的初步可视化

import randomimport numpy as npimport pandas as pdimport stringimport matplotlib.pyplot as pltimport networkx as nximport seaborn as snsfrom matplotlib.lines import Line2Dimport warningswarnings.filterwarnings("ignore")
# 1. 初始化随机数种子seed = random.randint(1, 100)random.seed(seed)
# 2. 设置矩阵大小num = 15
# 3. 创建随机矩阵和DataFramematrix = np.random.rand(num, num)df = pd.DataFrame(matrix, columns=list(string.ascii_uppercase[:num]), # 列名称为前15个大写字母 index=list(string.ascii_uppercase[:num])) # 行名称为前15个大写字母
# 4. 定义p值列表p = [0.001, 0.005, 0.02, 0.05, 0.1, 0.2]p_values = pd.DataFrame({ 'spec01': [random.choice(p) for _ in range(num)], 'spec02': [random.choice(p) for _ in range(num)], 'spec03': [random.choice(p) for _ in range(num)], 'spec04': [random.choice(p) for _ in range(num)],}, index=df.columns)
# 5. 生成r值的随机矩阵r_values = pd.DataFrame({ 'spec01': [random.random() * 2 - 1 for _ in range(num)], # 生成 [-1, 1] 之间的随机数 'spec02': [random.random() * 2 - 1 for _ in range(num)], 'spec03': [random.random() * 2 - 1 for _ in range(num)], 'spec04': [random.random() * 2 - 1 for _ in range(num)],}, index=df.columns)
# 6. 计算相关矩阵,并仅保留下三角矩阵,掩盖对角线corr_matrix = df.corr()
# 将对角线部分(相关系数为1)设置为NaN,使其不显示np.fill_diagonal(corr_matrix.values, np.nan)
# 掩盖上三角部分mask = np.triu(np.ones_like(corr_matrix, dtype=bool)) # 仅掩盖上三角部分
# 7. 使用Seaborn绘制仅显示下三角的相关性矩阵热力图,保留对角线但为空白fig, ax = plt.subplots(figsize=(8, 8))sns.heatmap(corr_matrix, annot=True, cmap='RdYlGn_r', ax=ax, square=True, cbar=True, fmt=".2f", mask=mask, linewidths=1, linecolor='white', annot_kws={"size": 8}, cbar_kws={"shrink": 0.8, "label": "Correlation"})
# 8. 构建网络图G = nx.Graph()G.add_nodes_from(p_values.columns) # 添加spec列作为节点G.add_nodes_from(p_values.index) # 添加index作为节点
# 9. 根据r值和p值构建边,并设置规则,包括透明度for spec in p_values.columns: for element in p_values.index: # 根据p值设置颜色和透明度 p_val = p_values.loc[element, spec] r_val = abs(r_values.loc[element, spec]) # 获取对应的r值并取绝对值 if p_val < 0.001: color = '#b2d8b2' # 淡绿色 alpha = 0.9 # 接近不透明 elif 0.001 <= p_val < 0.01: color = '#d3d3d3' # 淡灰色 alpha = 0.7 # 半透明 elif 0.01 <= p_val < 0.05: color = '#add8e6' # 淡蓝色 alpha = 0.5 # 更透明 else: color = '#dda0dd' # 淡紫色 alpha = 0.3 # 非常透明
# 根据r值的绝对值调整边的权重(粗细),例如r值范围在[0, 1]之间,可以映射为[2, 6]的边宽度 if r_val > 0.5: weight = 6 # 粗线条 elif 0.25 <= r_val <= 0.5: weight = 4 # 中等线条 else: weight = 2 # 细线条
G.add_edge( spec, element, weight=weight, # 根据r值设置边的权重 color=color, # 根据p值设置边的颜色 alpha=alpha # 根据p值设置透明度 )
# 10. 提取边的权重、颜色和透明度weights = [G[u][v]['weight'] for u, v in G.edges()]colors = [G[u][v]['color'] for u, v in G.edges()]alphas = [G[u][v]['alpha'] for u, v in G.edges()]
# 11. 定义节点位置# 将矩阵内部的节点放置在对角线上pos = {}for i, node in enumerate(p_values.index): # 将节点放置在对角线上,即(i, i)处,对应热力图的网格坐标 pos[node] = (i + 0.5, i + 0.5)
# 定义spec01, spec02, spec03, spec04的位置,放在指定的坐标点spec_positions = { 'spec01': (9, 1), # (9, 1) 'spec02': (11, 3), # (11, 3) 'spec03': (13, 5), # (13, 5) 'spec04': (14.5, 7), # (15.5, 7)}
# 更新spec节点的坐标pos.update(spec_positions)
# 12. 绘制节点nx.draw_networkx_nodes(G, pos, node_color='pink', # 节点颜色 node_size=300, # 节点大小 edgecolors='black', # 节点边框颜色 ax=ax)
# 13. 绘制边,包含透明度for i, (u, v) in enumerate(G.edges()): nx.draw_networkx_edges(G, pos, edgelist=[(u, v)], # 单独绘制每条边 width=weights[i], # 根据r值设置边的宽度 edge_color=colors[i], # 边的颜色 alpha=alphas[i], # 边的透明度 ax=ax, connectionstyle='arc3,rad=0.2') # 弧形边的连接样式
# 14. 添加文本标注for node, (x, y) in pos.items(): ax.text(x, y, node, fontsize=10, ha='center', va='center')
# 15. 添加线条颜色和宽度图例# P-value 图例元素p_value_legend_elements = [ Line2D([0], [0], color='#b2d8b2', lw=4, label='p < 0.001'), # 绿色 Line2D([0], [0], color='#d3d3d3', lw=4, label='0.001 <= p < 0.01'), # 灰色 Line2D([0], [0], color='#add8e6', lw=4, label='0.01 <= p < 0.05'), # 蓝色 Line2D([0], [0], color='#dda0dd', lw=4, label='p >= 0.05') # 紫色]
# Mantel's r 图例元素r_value_legend_elements = [ Line2D([0], [0], color='grey', lw=6, label='|r| > 0.5'), # 线条宽度 6 Line2D([0], [0], color='grey', lw=4, label='0.25 <= |r| <= 0.5'), # 线条宽度 4 Line2D([0], [0], color='lightgrey', lw=2, label='|r| < 0.25') # 线条宽度 2]
# 合并两个图例元素combined_legend_elements = p_value_legend_elements + r_value_legend_elements
# 创建一个综合的图例ax.legend(handles=combined_legend_elements, loc='upper left', bbox_to_anchor=(1.2, 1), title="P-value and |Mantel's r|")plt.savefig("3.pdf", format='pdf', bbox_inches='tight', dpi=1200)plt.show()

通过热力图和网络图直观展示了随机生成矩阵的相关性与节点间的关联。首先,代码直接生成相关矩阵并绘制热力图,显示矩阵中各元素的相关性。随后,基于p值和r值构建网络图,节点间的边根据统计显著性(p值)和相关强度(r值)来设置颜色、透明度和粗细,展示了不同节点之间的关联强度和统计显著性,最终,通过图例帮助解释不同边的特性,使得数据的相关性和统计意义更加清晰易懂,图中的线条其实应该根据Mantel检验的p值和r值绘制,p值决定了线条的颜色梯度,而r值则控制线条的粗细,接下来生成模拟数据进行Mantel检验并进行绘图

真实场景下绘图

模拟数据生成

import numpy as npimport pandas as pdimport seaborn as snsimport matplotlib.pyplot as pltimport networkx as nx
# 设置随机种子,保证结果可重复np.random.seed(42)
# 生成模拟的环境变量数据env_variables = ['N', 'P', 'K', 'Ca', 'Mg', 'S', 'Al', 'Fe', 'Mn', 'Zn', 'Mo', 'BS', 'HD', 'pH']n_samples = 100
# 创建一个包含环境变量的随机数据集env_data = pd.DataFrame(np.random.rand(n_samples, len(env_variables)), columns=env_variables)
# 生成模拟的物种数据species = ['Spec01', 'Spec02', 'Spec03']species_data = pd.DataFrame({ 'Spec01': np.random.rand(n_samples), 'Spec02': np.random.rand(n_samples), 'Spec03': np.random.rand(n_samples)})
# 将环境变量和物种数据合并combined_data = pd.concat([env_data, species_data], axis=1)combined_data

通过生成包含14种环境变量和3种物种丰度的随机数据集,并将它们合并为一个完整的数据框,为后续的生态学分析、物种与环境变量之间关系的研究以及可视化提供了初步的基础数据支持

Mantel检验

import mantelfrom scipy.spatial.distance import pdist, squareform
# 初始化存储p值和相关系数的DataFramep_values_df = pd.DataFrame(index=env_variables, columns=species)r_values_df = pd.DataFrame(index=env_variables, columns=species)
# 逐个环境变量计算与每个物种之间的Mantel检验p值和相关系数for env_var in env_variables: # 环境变量的距离矩阵 env_distance_matrix = squareform(pdist(combined_data[[env_var]], metric='euclidean')) for spec in species: # 物种的距离矩阵 species_distance_matrix = squareform(pdist(combined_data[[spec]], metric='euclidean')) # 使用mantel库计算Mantel检验 mantel_stat, p_value, _ = mantel.test(env_distance_matrix, species_distance_matrix, perms=999, method='pearson')
# 将p值和Mantel's r值存储在对应的DataFrame中 p_values_df.loc[env_var, spec] = p_value r_values_df.loc[env_var, spec] = mantel_stat
p_values_df

r_values_df

计算每个环境变量与每个物种之间的Mantel检验,并将计算得到的p值和Mantel的相关系数(r值)存储在两个DataFrame中,以便后续绘图

图形绘制

# 6. 计算相关矩阵,并仅保留下三角矩阵,掩盖对角线corr_matrix = env_data.corr()
# 将对角线部分(相关系数为1)设置为NaN,使其不显示np.fill_diagonal(corr_matrix.values, np.nan)
# 掩盖上三角部分mask = np.triu(np.ones_like(corr_matrix, dtype=bool)) # 仅掩盖上三角部分
# 7. 使用Seaborn绘制仅显示下三角的相关性矩阵热力图,保留对角线但为空白fig, ax = plt.subplots(figsize=(8, 8))sns.heatmap(corr_matrix, annot=True, cmap='RdYlGn_r', ax=ax, square=True, cbar=True, fmt=".2f", mask=mask, linewidths=1, linecolor='white', annot_kws={"size": 8}, cbar_kws={"shrink": 0.8, "label": "Correlation"})
# 8. 构建网络图G = nx.Graph()G.add_nodes_from(p_values_df.columns) # 添加spec列作为节点G.add_nodes_from(p_values_df.index) # 添加index作为节点
# 9. 根据r值和p值构建边,并设置规则,包括透明度for spec in p_values_df.columns: for element in p_values_df.index: # 根据p值设置颜色和透明度 p_val = p_values_df.loc[element, spec] r_val = abs(r_values_df.loc[element, spec]) # 获取对应的r值并取绝对值 if p_val < 0.001: color = '#b2d8b2' # 淡绿色 alpha = 0.9 # 接近不透明 elif 0.001 <= p_val < 0.01: color = '#d3d3d3' # 淡灰色 alpha = 0.7 # 半透明 elif 0.01 <= p_val < 0.05: color = '#add8e6' # 淡蓝色 alpha = 0.5 # 更透明 else: color = '#dda0dd' # 淡紫色 alpha = 0.3 # 非常透明
# 根据r值的绝对值调整边的权重(粗细),例如r值范围在[0, 1]之间,可以映射为[2, 6]的边宽度 if r_val > 0.5: weight = 6 # 粗线条 elif 0.25 <= r_val <= 0.5: weight = 4 # 中等线条 else: weight = 2 # 细线条
G.add_edge( spec, element, weight=weight, # 根据r值设置边的权重 color=color, # 根据p值设置边的颜色 alpha=alpha # 根据p值设置透明度 )
# 10. 提取边的权重、颜色和透明度weights = [G[u][v]['weight'] for u, v in G.edges()]colors = [G[u][v]['color'] for u, v in G.edges()]alphas = [G[u][v]['alpha'] for u, v in G.edges()]
# 11. 定义节点位置# 将矩阵内部的节点放置在对角线上pos = {}for i, node in enumerate(p_values_df.index): # 将节点放置在对角线上,即(i, i)处,对应热力图的网格坐标 pos[node] = (i + 0.5, i + 0.5)
# 定义spec01, spec02, spec03, spec04的位置,放在指定的坐标点spec_positions = { 'Spec01': (9, 1), # (9, 1) 'Spec02': (11, 3), # (11, 3) 'Spec03': (13, 5), # (13, 5)}
# 更新spec节点的坐标pos.update(spec_positions)
# 12. 绘制节点nx.draw_networkx_nodes(G, pos, node_color='pink', # 节点颜色 node_size=300, # 节点大小 edgecolors='black', # 节点边框颜色 ax=ax)
# 13. 绘制边,包含透明度for i, (u, v) in enumerate(G.edges()): nx.draw_networkx_edges(G, pos, edgelist=[(u, v)], # 单独绘制每条边 width=weights[i], # 根据r值设置边的宽度 edge_color=colors[i], # 边的颜色 alpha=alphas[i], # 边的透明度 ax=ax, connectionstyle='arc3,rad=0.2') # 弧形边的连接样式
# 14. 添加文本标注for node, (x, y) in pos.items(): ax.text(x, y, node, fontsize=10, ha='center', va='center')
# 15. 添加线条颜色和宽度图例# P-value 图例元素p_value_legend_elements = [ Line2D([0], [0], color='#b2d8b2', lw=4, label='p < 0.001'), # 绿色 Line2D([0], [0], color='#d3d3d3', lw=4, label='0.001 <= p < 0.01'), # 灰色 Line2D([0], [0], color='#add8e6', lw=4, label='0.01 <= p < 0.05'), # 蓝色 Line2D([0], [0], color='#dda0dd', lw=4, label='p >= 0.05') # 紫色]
# Mantel's r 图例元素r_value_legend_elements = [ Line2D([0], [0], color='grey', lw=6, label='|r| > 0.5'), # 线条宽度 6 Line2D([0], [0], color='grey', lw=4, label='0.25 <= |r| <= 0.5'), # 线条宽度 4 Line2D([0], [0], color='lightgrey', lw=2, label='|r| < 0.25') # 线条宽度 2]
# 合并两个图例元素combined_legend_elements = p_value_legend_elements + r_value_legend_elements
# 创建一个综合的图例ax.legend(handles=combined_legend_elements, loc='upper left', bbox_to_anchor=(1.2, 1), title="P-value and |Mantel's r|")plt.savefig("4.pdf", format='pdf', bbox_inches='tight', dpi=1200)plt.show()

至此,该图形已成功绘制。如果读者使用自己的数据,可能由于数据维度不同,需要调整节点位置,即可生成相应的图表。然而,为了更简便地生成类似文献中图B的复杂网络图,建议使用R语言。R具备针对这类可视化的专用绘图包,如ggplot2、vegan和corrplot,能够直接处理Mantel检验和环境变量的可视化,操作简便。而Python虽然可以通过networkx和matplotlib生成类似图表,但过程较为复杂,且代码量较大,效率不及R语言

往期推荐

SCI图表复现:整合数据分布与相关系数的高级可视化策略
复现顶刊Streamlit部署预测模型APP
树模型系列:如何通过XGBoost提取特征贡献度
SHAP进阶解析:机器学习、深度学习模型解释保姆级教程
特征选择:Lasso和Boruta算法的结合应用
从基础到进阶:优化SHAP力图,让样本解读更直观
SCI图表复现:优化SHAP特征贡献图展示更多模型细节
多模型中的特征贡献度比较与可视化图解

基于相关性与标准差的多模型评价指标可视化比较 —— 泰勒图应用解析

微信号|deep_ML

欢迎添加作者微信进入Python、ChatGPT群

进群请备注Python或AI进入相关群

无需科学上网、同步官网所有功能、使用无限

如果你对类似于这样的文章感兴趣。

欢迎关注、点赞、转发~

个人观点,仅供参考

RPython
人生苦短,R和Python。
 最新文章