火山图是代谢组学数据分析中用于可视化差异代谢物的核心工具，它通过散点图形式，将每个代谢物的变化倍数（通常为log2(Fold Change)）与统计显著性（-log10(p-value)）相结合，从而直观地区分上调、下调及不显著的代谢物 [ref_1][ref_4]。以下是如何使用Python处理代谢组学数据并绘制高质量火山图的详细教程，涵盖数据预处理、差异分析、可视化及高级定制。 ### 一、 代谢组学数据准备与预处理 典型的代谢组学差异分析结果数据应至少包含代谢物名称、变化倍数（FC或log2FC）和显著性p值（或校正后的p值如FDR）[ref_5]。 ```python import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns from scipy import stats import warnings warnings.filterwarnings('ignore') # 1. 载入代谢组学差异分析结果 # 假设数据格式为CSV，包含以下列：'Metabolite', 'log2FC', 'pvalue', 'VIP' (可选) df = pd.read_csv('metabolomics_diff_results.csv') print("数据预览：") print(df.head()) print(f"\n数据形状：{df.shape}") # 2. 数据预处理与检查 # 确保必要的列存在，VIP值常见于基于PLS-DA模型的代谢组学分析 [ref_4] required_cols = ['Metabolite', 'log2FC', 'pvalue'] if not all(col in df.columns for col in required_cols): # 尝试常见列名变体 col_map = {'log2FoldChange': 'log2FC', 'logFC': 'log2FC', 'P.Value': 'pvalue', 'p.value': 'pvalue', 'FDR': 'pvalue'} df = df.rename(columns=col_map) # 如果还没有log2FC，但从FC计算 if 'FC' in df.columns and 'log2FC' not in df.columns: df['log2FC'] = np.log2(df['FC']) # 处理缺失值或无限值 df = df.replace([np.inf, -np.inf], np.nan).dropna(subset=['log2FC', 'pvalue']) df['pvalue'] = pd.to_numeric(df['pvalue'], errors='coerce') df = df.dropna(subset=['pvalue']) # 计算 -log10(pvalue) 用于绘图 df['neg_log10_pvalue'] = -np.log10(df['pvalue']) print(f"处理后数据量：{len(df)}") ``` ### 二、 定义差异代谢物筛选阈值与分类 在代谢组学中，阈值设定需结合统计学显著性和生物学变化幅度，VIP值（Variable Importance in Projection）常作为多元统计分析的重要补充指标 [ref_4]。 ```python # 3. 设定差异筛选阈值 log2fc_thresh = 1.0 # 对应原始FC=2倍变化 pval_thresh = 0.05 vip_thresh = 1.0 if 'VIP' in df.columns else None # VIP阈值通常为1.0或1.5 [ref_4] # 4. 根据阈值对代谢物进行分类 def classify_metabolite(row): sig = row['pvalue'] < pval_thresh fc_up = row['log2FC'] >= log2fc_thresh fc_down = row['log2FC'] <= -log2fc_thresh if vip_thresh and 'VIP' in row: vip_sig = row['VIP'] >= vip_thresh # 结合VIP值、p值和FC进行综合筛选 [ref_4] if sig and vip_sig and fc_up: return 'Up-regulated (VIP)' elif sig and vip_sig and fc_down: return 'Down-regulated (VIP)' elif sig and fc_up: return 'Up-regulated' elif sig and fc_down: return 'Down-regulated' else: return 'Not Significant' else: # 仅使用p值和FC if sig and fc_up: return 'Up-regulated' elif sig and fc_down: return 'Down-regulated' else: return 'Not Significant' df['Regulation'] = df.apply(classify_metabolite, axis=1) regulation_counts = df['Regulation'].value_counts() print("差异代谢物分类统计：") print(regulation_counts) ``` ### 三、 绘制基础火山图 使用`matplotlib`和`seaborn`绘制清晰直观的火山图。 ```python # 5. 设置绘图风格与颜色 plt.style.use('seaborn-v0_8-whitegrid') # 定义颜色方案：上调-红，下调-蓝，VIP显著-深色，不显著-灰 if vip_thresh and 'VIP' in df.columns: palette = { 'Up-regulated (VIP)': '#CB181D', # 深红 'Down-regulated (VIP)': '#2171B5', # 深蓝 'Up-regulated': '#FB9A99', # 浅红 'Down-regulated': '#9ECAE1', # 浅蓝 'Not Significant': '#D9D9D9' # 灰 } else: palette = { 'Up-regulated': '#E64B35', 'Down-regulated': '#3182BD', 'Not Significant': '#BDBDBD' } # 6. 创建火山图 fig, ax = plt.subplots(figsize=(10, 8)) # 按分类绘制散点，便于图例顺序控制 categories = sorted(df['Regulation'].unique(), key=lambda x: ('Not Significant' in x, x)) for category in categories: if category not in palette: continue subset = df[df['Regulation'] == category] ax.scatter( subset['log2FC'], subset['neg_log10_pvalue'], c=palette[category], label=category, alpha=0.7, s=40, # 点大小 edgecolor='black', linewidth=0.5 ) # 7. 添加阈值线 # 垂直阈值线 (log2FC) ax.axvline(x=log2fc_thresh, color='gray', linestyle='--', linewidth=1.2, alpha=0.8) ax.axvline(x=-log2fc_thresh, color='gray', linestyle='--', linewidth=1.2, alpha=0.8) # 水平阈值线 (-log10(p-value)) pval_thresh_log10 = -np.log10(pval_thresh) ax.axhline(y=pval_thresh_log10, color='gray', linestyle='--', linewidth=1.2, alpha=0.8) # 可选：添加VIP阈值区域（如果使用VIP值） if vip_thresh and 'VIP' in df.columns: # 可以用不同形状或添加文本标注来提示VIP筛选 ax.text(ax.get_xlim()[1]*0.7, ax.get_ylim()[1]*0.9, f'VIP ≥ {vip_thresh}', fontsize=10, style='italic', bbox=dict(boxstyle="round,pad=0.3", facecolor='yellow', alpha=0.3)) # 8. 标记Top N显著代谢物（避免标签重叠） top_n = 15 # 选择策略：p值最小且|log2FC|较大的代谢物 df['abs_log2fc'] = df['log2FC'].abs() top_metabolites = df[df['Regulation'].str.contains('regulated')].nsmallest(top_n, 'pvalue') # 简单文本标注（生产环境建议使用adjustText库优化位置） for idx, row in top_metabolites.iterrows(): ax.annotate( row['Metabolite'], xy=(row['log2FC'], row['neg_log10_pvalue']), xytext=(5, 5), # 偏移量 textcoords='offset points', fontsize=8, alpha=0.85, arrowprops=dict(arrowstyle='->', color='black', lw=0.5, alpha=0.5) ) ``` ### 四、 图表美化与信息增强 优化图表元素，使其符合出版要求并包含完整信息 [ref_2]。 ```python # 9. 设置坐标轴与标签 ax.set_xlabel('log$_2$(Fold Change)', fontsize=14, fontweight='bold') ax.set_ylabel('-log$_{10}$(p-value)', fontsize=14, fontweight='bold') ax.set_title('Volcano Plot of Metabolomics Data', fontsize=16, fontweight='bold', pad=20) # 设置坐标轴范围，留出一些边距 x_min, x_max = df['log2FC'].min(), df['log2FC'].max() y_max = df['neg_log10_pvalue'].max() ax.set_xlim(x_min * 1.1, x_max * 1.1) ax.set_ylim(0, y_max * 1.05) # 10. 添加网格与图例 ax.grid(True, which='both', linestyle=':', linewidth=0.5, alpha=0.4) # 将图例放在图表外侧，避免遮挡数据 ax.legend( title='Metabolite Regulation', title_fontsize=12, fontsize=10, loc='center left', bbox_to_anchor=(1.02, 0.5), # 将图例放在图表右侧外部 frameon=True, fancybox=True, shadow=False, borderaxespad=0. ) # 11. 在图表内添加统计信息摘要框 summary_text = (f'Total Metabolites: {len(df):,}\n' f'Up-regulated: {regulation_counts.get("Up-regulated", 0) + regulation_counts.get("Up-regulated (VIP)", 0):,}\n' f'Down-regulated: {regulation_counts.get("Down-regulated", 0) + regulation_counts.get("Down-regulated (VIP)", 0):,}\n' f'Thresholds: |log2FC|>{log2fc_thresh}, p<{pval_thresh}') if vip_thresh and 'VIP' in df.columns: summary_text += f'\nVIP≥{vip_thresh}' props = dict(boxstyle='round', facecolor='lightgrey', alpha=0.8) ax.text(0.02, 0.98, summary_text, transform=ax.transAxes, fontsize=9, verticalalignment='top', bbox=props, family='monospace') plt.tight_layout(rect=[0, 0, 0.85, 1]) # 为右侧图例留出空间 plt.show() ``` ### 五、 高级定制与交互式可视化 对于更复杂的分析或探索性数据挖掘，可以使用`plotly`库创建交互式火山图。 ```python # 12. 使用Plotly创建交互式火山图（可选） try: import plotly.express as px import plotly.graph_objects as go # 为交互式图表准备数据 plot_df = df.copy() # 添加悬停文本信息 plot_df['hover_text'] = plot_df.apply( lambda row: f"<b>{row['Metabolite']}</b><br>" f"log2FC: {row['log2FC']:.3f}<br>" f"p-value: {row['pvalue']:.2e}<br>" f"-log10(p): {row['neg_log10_pvalue']:.2f}<br>" f"Regulation: {row['Regulation']}" + (f"<br>VIP: {row['VIP']:.2f}" if 'VIP' in row else ""), axis=1 ) fig_interactive = px.scatter( plot_df, x='log2FC', y='neg_log10_pvalue', color='Regulation', color_discrete_map=palette, hover_name='Metabolite', hover_data={'log2FC': ':.3f', 'pvalue': ':.2e', 'neg_log10_pvalue': ':.2f', 'Regulation': True}, title='Interactive Volcano Plot - Metabolomics Data', labels={'log2FC': 'log2(Fold Change)', 'neg_log10_pvalue': '-log10(p-value)'}, height=700 ) # 添加阈值线 fig_interactive.add_hline(y=pval_thresh_log10, line_dash="dash", line_color="gray", opacity=0.7, annotation_text=f"p={pval_thresh}", annotation_position="top left") fig_interactive.add_vline(x=log2fc_thresh, line_dash="dash", line_color="gray", opacity=0.7, annotation_text=f"FC={2**log2fc_thresh:.1f}", annotation_position="top right") fig_interactive.add_vline(x=-log2fc_thresh, line_dash="dash", line_color="gray", opacity=0.7, annotation_text=f"FC={2**(-log2fc_thresh):.2f}", annotation_position="top left") # 更新布局 fig_interactive.update_layout( legend_title_text='Regulation', plot_bgcolor='white', hoverlabel=dict(bgcolor="white", font_size=12), xaxis=dict(showgrid=True, gridcolor='lightgray', zerolinecolor='black'), yaxis=dict(showgrid=True, gridcolor='lightgray', zerolinecolor='black') ) fig_interactive.show() # 保存交互式HTML文件 # fig_interactive.write_html("interactive_volcano_plot.html") except ImportError: print("Plotly未安装，如需交互式图表请运行: pip install plotly") ``` ### 六、 结果解读与下游分析 绘制火山图后，关键的差异代谢物列表可用于后续通路富集分析或生物学解释 [ref_5]。 ```python # 13. 提取并保存显著差异代谢物列表 significant_up = df[df['Regulation'].str.contains('Up-regulated')] significant_down = df[df['Regulation'].str.contains('Down-regulated')] print(f"\n显著上调代谢物数量：{len(significant_up)}") print(f"显著下调代谢物数量：{len(significant_down)}") # 保存结果 output_cols = ['Metabolite', 'log2FC', 'pvalue', 'neg_log10_pvalue', 'Regulation'] if 'VIP' in df.columns: output_cols.append('VIP') significant_up[output_cols].sort_values('pvalue').to_csv('significant_up_metabolites.csv', index=False) significant_down[output_cols].sort_values('pvalue').to_csv('significant_down_metabolites.csv', index=False) print("差异代谢物列表已保存为CSV文件。") # 14. 可选：生成简单的统计摘要 summary_stats = pd.DataFrame({ 'Metric': ['Total Metabolites', 'Significant (p<0.05)', 'Up-regulated', 'Down-regulated', 'Max -log10(p)', 'Max |log2FC|'], 'Value': [len(df), len(significant_up) + len(significant_down), len(significant_up), len(significant_down), df['neg_log10_pvalue'].max(), df['log2FC'].abs().max()] }) print("\n分析摘要：") print(summary_stats.to_string(index=False)) ``` ### 七、 关键参数与注意事项总结 下表总结了代谢组学火山图绘制中的关键参数及其意义： | 参数/步骤 | 说明与建议 | 代谢组学特定考量 | | :--- | :--- | :--- | | **数据输入** | 需包含代谢物ID、log2FC、p值。VIP值作为重要补充 [ref_4]。 | 数据应经过适当的预处理和标准化 [ref_5]。 | | **log2FC阈值** | 通常设为0.5 (FC~1.4) 或 1.0 (FC=2)。 | 可根据实验体系和代谢物特性调整，脂质组学可能使用不同阈值。 | | **p值阈值** | 常设为0.05。建议同时考虑FDR校正 [ref_5]。 | 对于高通量数据，使用校正后的p值 (如FDR, q-value) 可控制假阳性。 | | **VIP阈值** | 常用于PLS-DA模型，阈值常设为1.0或1.5 [ref_4]。 | VIP值高不一定代表生物学意义重大，需结合p值和FC综合判断 [ref_4]。 | | **可视化优化** | 使用区分度高的颜色，添加阈值线和统计摘要。 | 可考虑用点大小映射VIP值，用形状映射代谢物类别（如脂质、氨基酸）。 | | **结果解读** | 关注落在阈值线右上/左上象限的代谢物。 | 结合已知通路，对显著差异代谢物进行通路富集分析，挖掘生物学意义 [ref_5]。 | 通过以上步骤，您可以系统性地使用Python对代谢组学数据进行差异分析可视化。核心在于根据生物学问题合理设置阈值，并利用火山图直观地筛选出在变化幅度和统计学意义上均显著的潜在标志物 [ref_3][ref_4]。