如何用R的hierfstat计算各种群的种群特异性FST及基础统计量

一、核心结论 #

你可以通过单个种群vs其余所有种群合并组的对比思路，用hierfstat计算种群特异性FST；你之前尝试拆分种群的方法是合理的，但需要明确统计逻辑并注意计算细节。

二、具体实现步骤 #

假设你的hierfstat对象为hs_obj，种群分组信息存储在hs_obj$pop中：

1. 提取所有唯一种群名称

pop_names <- unique(hs_obj$pop)

2. 循环计算每个种群与其余种群的特异性FST及配套统计量

# 初始化结果容器
pop_specific_stats <- list()

for (pop in pop_names) {
  # 构建临时分组：目标种群为一组，其余种群合并为"Rest"组
  temp_pop_vector <- ifelse(hs_obj$pop == pop, pop, "Rest")
  # 更新临时hierfstat对象的分组信息
  temp_hs <- hs_obj
  temp_hs$pop <- temp_pop_vector
  
  # 计算两组间的统计量
  bs_result <- basic.stats(temp_hs)
  # 提取目标种群的FIS、Ho、He，以及两组间的FST
  pop_specific_stats[[pop]] <- data.frame(
    FST = bs_result$overall$FST,
    FIS = bs_result$Fis[names(bs_result$Fis) == pop],
    Ho = bs_result$Ho[names(bs_result$Ho) == pop],
    He = bs_result$Hs[names(bs_result$Hs) == pop]
  )
}

# 转换为数据框便于查看和后续分析
pop_stats_summary <- do.call(rbind, pop_specific_stats)

3. 可选：计算种群间两两FST
   如果需要单个种群与其他每个种群的两两分化程度，直接使用pairwise.fst()：

pairwise_fst_matrix <- pairwise.fst(hs_obj)
# 提取每个种群与其他种群的平均FST作为另一种特异性指标
pop_avg_pairwise_fst <- apply(pairwise_fst_matrix, 1, mean, na.rm = TRUE)

三、拆分种群方法的合理性与注意事项 #

• 合理性：拆分种群后将单个种群与其余合并组对比计算FST，本质是衡量该种群相对于整体的分化程度，是种群特异性分化的常用量化方式，逻辑上成立。
• 样本量影响：若不同种群样本量差异极大，合并其余种群可能稀释小种群的信号，建议单独标注样本量极小的种群结果，谨慎解读。
• 大数据集优化：你的数据包含38000个位点，循环计算可能较慢，可改用lapply替代循环，或用parallel包实现并行计算提升效率。
• FIS的计算：FIS是种群内近交系数，与其他种群无关，直接用basic.stats(hs_obj)$Fis即可得到每个种群的FIS，无需拆分种群重复计算。

内容的提问来源于stack exchange，提问作者Pei-Wei Sun