技术问询：如何在PCoA分析图中添加性状向量/箭头？

嘿，看你已经把PCoA性状分析的前期流程走得很顺畅了——从构建物种×性状矩阵、考虑性状权重拆分，到用lingoes法校正负特征值，这几步都踩得很准！接下来要绘制带性状向量的PCoA图，我给你分享两个常用工具（R和Python）的实现方法，都是生态领域里常用的套路：

R里的vegan包是生态分析的利器，能轻松搞定PCoA计算和性状向量拟合：

• 首先，用cmdscale完成带lingoes校正的PCoA计算：

  # 假设你的距离矩阵已经构建好，命名为dist_mat
  pcoa_result <- cmdscale(dist_mat, k = 2, eig = TRUE, add = TRUE, x.ret = TRUE)
  # 参数说明：k=2指定保留前2个轴，add=TRUE就是启用lingoes校正
  

• 接着用envfit函数拟合性状与PCoA轴的相关性：

  library(vegan)
  # trait_mat是你的物种×性状矩阵（行=物种，列=性状）
  trait_fit <- envfit(pcoa_result$points, trait_mat, permutations = 999)
  # permutations设置置换次数，用于计算显著性
  

• 最后绘制PCoA图并添加性状向量：

  # 先绘制物种点
  plot(pcoa_result$points, pch = 16, col = "darkblue", 
       xlab = paste0("PCoA 1 (", round(pcoa_result$eig[1]/sum(pcoa_result$eig)*100, 1), "%)"),
       ylab = paste0("PCoA 2 (", round(pcoa_result$eig[2]/sum(pcoa_result$eig)*100, 1), "%)"),
       main = "PCoA with Trait Vectors")
  # 添加显著的性状向量（p<0.05）
  plot(trait_fit, p.max = 0.05, col = "red", lwd = 2)
  # 加个图例更清晰
  legend("topright", legend = c("Species", "Significant Traits"), 
         pch = c(16, NA), lty = c(NA, 1), col = c("darkblue", "red"))
  

注意：性状矩阵的物种顺序必须和PCoA分析中的物种顺序完全一致，否则相关性计算会出错！如果是分类性状，记得先转换成哑变量再输入envfit。

Python可以用scikit-bio做PCoA，结合matplotlib手动绘制性状向量：

• 第一步，用scikit-bio完成lingoes校正的PCoA：

  import skbio
  # dist_mat是你的距离矩阵（可以是skbio的DistanceMatrix对象，或者普通二维数组）
  pcoa = skbio.stats.ordination.pcoa(dist_mat, correction='lingoes')
  

• 第二步，计算性状与PCoA轴的相关性（这里用皮尔逊相关）：

  import pandas as pd
  import scipy.stats as stats
  import matplotlib.pyplot as plt
  
  # trait_mat是DataFrame格式，行=物种，列=性状；pcoa.samples是PCoA坐标表
  trait_correlations = {}
  for trait in trait_mat.columns:
      # 计算性状与PC1、PC2的相关性和p值
      corr_pc1, p_pc1 = stats.pearsonr(pcoa.samples['PC1'], trait_mat[trait])
      corr_pc2, p_pc2 = stats.pearsonr(pcoa.samples['PC2'], trait_mat[trait])
      trait_correlations[trait] = (corr_pc1, corr_pc2, p_pc1, p_pc2)
  
  # 转换成DataFrame方便筛选显著性状
  corr_df = pd.DataFrame.from_dict(trait_correlations, orient='index', 
                                   columns=['PC1_corr', 'PC2_corr', 'PC1_p', 'PC2_p'])
  # 筛选p<0.05的显著性状
  sig_traits = corr_df[(corr_df['PC1_p'] < 0.05) | (corr_df['PC2_p'] < 0.05)]
  

• 第三步，绘制PCoA图并添加性状箭头：

  fig, ax = plt.subplots(figsize=(8, 6))
  # 绘制物种点
  ax.scatter(pcoa.samples['PC1'], pcoa.samples['PC2'], c='darkblue', alpha=0.7)
  # 计算轴的范围，用于缩放箭头长度
  x_range = pcoa.samples['PC1'].max() - pcoa.samples['PC1'].min()
  y_range = pcoa.samples['PC2'].max() - pcoa.samples['PC2'].min()
  
  # 添加性状箭头
  for trait, row in sig_traits.iterrows():
      # 箭头长度用相关性乘以轴范围的0.8倍，避免超出图范围
      ax.arrow(0, 0, row['PC1_corr'] * x_range * 0.8, row['PC2_corr'] * y_range * 0.8,
               head_width=0.05*x_range, head_length=0.05*y_range, fc='red', ec='red', lw=2)
      # 标注性状名称
      ax.text(row['PC1_corr'] * x_range * 0.9, row['PC2_corr'] * y_range * 0.9, 
              trait, color='red', fontweight='bold')
  
  # 设置轴标签和标题
  ax.set_xlabel(f"PCoA 1 ({pcoa.proportion_explained['PC1']*100:.1f}%)")
  ax.set_ylabel(f"PCoA 2 ({pcoa.proportion_explained['PC2']*100:.1f}%)")
  ax.set_title("PCoA with Trait Vectors")
  plt.grid(alpha=0.3)
  plt.show()
  

小提示：如果是分类性状，需要先做one-hot编码（比如用pandas.get_dummies），再计算相关性，不然结果会不准确。

内容的提问来源于stack exchange，提问作者Vincent H.