群体分化分析 — Fst 计算与解读¶
一句话概述：Fst 衡量群体间遗传分化程度（0=无分化，1=完全分化），是群体遗传学中最基本、最常用的统计量，可用 VCFtools、PLINK 和 PopGenome 计算。
核心知识点速查表¶
| 概念 | 白话解释 |
|---|---|
| Fst（固定指数） | 群体间遗传变异占总变异的比例 |
| Fst = 0 | 两个群体无分化（像一个群体） |
| Fst = 1 | 两个群体完全分化（无共享等位基因） |
| Weir & Cockerham Fst | 校正了样本量偏差的 Fst 估计方法 |
| Hudson Fst | 另一种 Fst 估计，适合不等样本量 |
| 滑动窗口 Fst | 沿基因组滑动计算局部 Fst，检测分化区域 |
| Fst outlier | Fst 异常高的基因组区域，可能受歧化选择 |
| PBS | 群体分支统计量，基于 Fst 检测特定谱系的选择 |
一、Fst 的概念（白话版）¶
比喻：假设有两个学校（群体），每个学生有不同的手环颜色（等位基因）。如果两个学校的手环颜色分布一样（红色30%、蓝色70%），那 Fst=0。如果学校A全是红色、学校B全是蓝色，那 Fst=1。
Fst 的经验解读： - Fst < 0.05 → 分化很小（同一物种不同群体） - 0.05~0.15 → 中等分化 - 0.15~0.25 → 大分化 - Fst > 0.25 → 非常大的分化（可能不同亚种）
二、方法 1：VCFtools 计算 Fst¶
# VCFtools 是计算 Fst 最简单的方法
# 安装
conda install -c bioconda vcftools # conda安装
# 准备群体标签文件（每行一个样本名）
echo -e "sample1\nsample2\nsample3" > pop1.txt # 群体1样本列表
echo -e "sample4\nsample5\nsample6" > pop2.txt # 群体2样本列表
# 计算全基因组平均 Fst（Weir & Cockerham方法）
vcftools --gzvcf input.vcf.gz \ # 输入VCF
--weir-fst-pop pop1.txt \ # 群体1
--weir-fst-pop pop2.txt \ # 群体2
--out pop1_vs_pop2 # 输出前缀
# 输出文件：pop1_vs_pop2.weir.fst
# 格式：CHROM POS WEIR_AND_COCKERHAM_FST
# 最后一行显示全基因组加权平均Fst
# 滑动窗口 Fst（检测高分化区域）
vcftools --gzvcf input.vcf.gz \
--weir-fst-pop pop1.txt \
--weir-fst-pop pop2.txt \
--fst-window-size 50000 \ # 窗口大小50kb
--fst-window-step 25000 \ # 步长25kb
--out pop1_vs_pop2_window # 输出前缀
# 输出文件：pop1_vs_pop2_window.windowed.weir.fst
# 格式：CHROM BIN_START BIN_END N_VARIANTS WEIGHTED_FST MEAN_FST
三、方法 2：PLINK 计算 Fst¶
# PLINK 2.0 支持 Fst 计算
# 准备群体信息（在 .fam 文件的 FID 列标注群体）
# 或使用 --within 参数指定群体分组
# 生成群体分配文件
cat << 'EOF' > clusters.txt
sample1 sample1 pop1
sample2 sample2 pop1
sample3 sample3 pop2
sample4 sample4 pop2
EOF
# 计算 Fst
plink2 --vcf input.vcf.gz \
--fst CATPHENO \ # 分类表型（群体分组）
--within clusters.txt \ # 群体分配文件
--out fst_result # 输出前缀
# PLINK 1.9 版本
plink --bfile mydata \
--fst \ # 计算Fst
--within clusters.txt \
--out fst_result
四、Fst Manhattan 图可视化¶
# ========== R脚本：Fst Manhattan plot ==========
# 读取滑动窗口Fst结果
fst <- read.table("pop1_vs_pop2_window.windowed.weir.fst",
header=TRUE) # 读取VCFtools输出
# 处理负值（Weir & Cockerham Fst 可以为负，设为0）
fst$WEIGHTED_FST[fst$WEIGHTED_FST < 0] <- 0 # 负值设为0
# 为每条染色体分配颜色
fst$CHR <- as.numeric(gsub("chr", "", fst$CHROM)) # 染色体编号
# 计算全基因组位置
fst$cum_pos <- 0 # 初始化
chr_sizes <- tapply(fst$BIN_END, fst$CHR, max) # 每条染色体的大小
offset <- cumsum(c(0, chr_sizes[-length(chr_sizes)])) # 偏移量
for(i in seq_along(unique(fst$CHR))) {
chr <- unique(fst$CHR)[i]
fst$cum_pos[fst$CHR == chr] <- fst$BIN_START[fst$CHR == chr] + offset[i]
}
# 画 Manhattan plot
pdf("Fst_manhattan.pdf", width=14, height=5)
plot(fst$cum_pos / 1e6, # X轴：位置(Mb)
fst$WEIGHTED_FST, # Y轴：Fst
pch=20, cex=0.5, # 小点
col=ifelse(fst$CHR %% 2 == 0, # 奇偶染色体不同颜色
"steelblue", "gray50"),
xlab="Genomic position (Mb)", # X轴标签
ylab=expression(F[ST]), # Y轴标签（下标ST）
main="Genome-wide Fst: Pop1 vs Pop2") # 标题
# 添加 top 1% 阈值线
threshold <- quantile(fst$WEIGHTED_FST, 0.99, na.rm=TRUE) # top 1%
abline(h=threshold, col="red", lty=2) # 红色虚线
dev.off()
五、成对 Fst 矩阵与热图¶
# 当有多个群体时，计算成对 Fst 矩阵
# 假设已经用 VCFtools 计算了所有群体对的 Fst
# 构建 Fst 矩阵
pops <- c("PopA", "PopB", "PopC", "PopD", "PopE") # 群体名列表
n <- length(pops)
fst_matrix <- matrix(0, n, n, dimnames=list(pops, pops)) # 初始化矩阵
# 填入 Fst 值（从 VCFtools 输出中提取加权平均 Fst）
fst_matrix["PopA", "PopB"] <- 0.05
fst_matrix["PopA", "PopC"] <- 0.12
fst_matrix["PopA", "PopD"] <- 0.08
fst_matrix["PopB", "PopC"] <- 0.15
fst_matrix["PopB", "PopD"] <- 0.10
fst_matrix["PopC", "PopD"] <- 0.03
# ... 填入所有成对 Fst
fst_matrix <- fst_matrix + t(fst_matrix) # 对称化
# 画热图
library(pheatmap) # 加载pheatmap包
pdf("Fst_heatmap.pdf", width=8, height=7)
pheatmap(fst_matrix, # Fst矩阵
display_numbers=TRUE, # 显示数值
number_format="%.3f", # 3位小数
color=colorRampPalette(c("white","orange","red"))(50), # 颜色
main="Pairwise Fst between populations") # 标题
dev.off()
六、PBS（群体分支统计量）¶
# PBS 基于三个群体的 Fst 计算特定谱系的分化
# PBS 高值 → 该谱系经历了强选择
# PBS 计算公式：
# T_AB = -log(1 - Fst_AB) # 分化时间的代理
# PBS_A = (T_AB + T_AC - T_BC) / 2
# 用R计算PBS
cat << 'EOF' > calc_pbs.R
# 读取三个群体对的窗口 Fst
fst_AB <- read.table("popA_vs_popB.windowed.weir.fst", header=TRUE)
fst_AC <- read.table("popA_vs_popC.windowed.weir.fst", header=TRUE)
fst_BC <- read.table("popB_vs_popC.windowed.weir.fst", header=TRUE)
# Fst 负值设为 0
fst_AB$WEIGHTED_FST[fst_AB$WEIGHTED_FST < 0] <- 0
fst_AC$WEIGHTED_FST[fst_AC$WEIGHTED_FST < 0] <- 0
fst_BC$WEIGHTED_FST[fst_BC$WEIGHTED_FST < 0] <- 0
# 确保 Fst < 1（避免 log(0)）
cap_fst <- function(x) pmin(x, 0.9999) # 限制最大值
T_AB <- -log(1 - cap_fst(fst_AB$WEIGHTED_FST)) # 转换为T
T_AC <- -log(1 - cap_fst(fst_AC$WEIGHTED_FST))
T_BC <- -log(1 - cap_fst(fst_BC$WEIGHTED_FST))
# 计算 PBS（以群体A为焦点）
PBS_A <- (T_AB + T_AC - T_BC) / 2 # PBS公式
# 输出结果
result <- data.frame(CHROM=fst_AB$CHROM,
BIN_START=fst_AB$BIN_START,
BIN_END=fst_AB$BIN_END,
PBS=PBS_A)
write.table(result, "PBS_popA.txt", row.names=FALSE, quote=FALSE, sep="\t")
EOF
Rscript calc_pbs.R
七、常见报错与解决¶
| 报错信息 | 原因 | 解决方法 |
|---|---|---|
| Fst 为负值 | Weir & Cockerham 方法在分化极小时可以为负 | 正常现象，分析时设为 0 |
No individuals in population | 群体文件中的样本名与VCF不匹配 | 检查样本名拼写 |
| 全基因组 Fst 接近 0 | 两个群体实际没什么分化 | 正常结果，不是错误 |
| 滑动窗口中有大量 NA | 某些窗口中 SNP 太少 | 减小窗口大小或用更密的SNP |
| Fst > 0.5 但物种一样 | 群体确实分化很大或数据有问题 | 检查是否有批次效应 |
八、面试高频问题¶
Q1：Fst 和遗传距离的区别？
Fst 衡量群体间分化比例，遗传距离（如 Nei's D）衡量两个群体间的绝对差异。Fst 受群体大小和迁移率影响，遗传距离更多反映分化时间。Fst 适合比较同一物种的群体，遗传距离适合物种间比较。
Q2：怎么用 Fst 检测自然选择？
Fst outlier 方法：在全基因组计算 Fst，找出 Fst 异常高的区域。这些区域可能受歧化选择（不同群体适应不同环境）。常用工具：BayeScan、OutFLANK、pcadapt。
Q3：Weir & Cockerham Fst 为什么可以是负数？
这是统计估计量，当群体内变异大于群体间变异时（即实际无分化），估计值可以为负。负值应解读为 Fst=0（无分化）。
九、速查表¶
# === Fst 计算速查 ===
# VCFtools 单SNP Fst
vcftools --gzvcf in.vcf.gz --weir-fst-pop pop1.txt --weir-fst-pop pop2.txt --out fst
# VCFtools 滑动窗口 Fst
vcftools --gzvcf in.vcf.gz --weir-fst-pop pop1.txt --weir-fst-pop pop2.txt \
--fst-window-size 50000 --fst-window-step 25000 --out fst_window
# 多群体成对 Fst（批量计算）
for p1 in popA popB popC; do
for p2 in popA popB popC; do
if [ "$p1" \< "$p2" ]; then
vcftools --gzvcf in.vcf.gz --weir-fst-pop ${p1}.txt --weir-fst-pop ${p2}.txt --out ${p1}_${p2}
fi
done
done
# Fst 解读标准：
# < 0.05 → 很小分化
# 0.05~0.15 → 中等分化
# 0.15~0.25 → 大分化
# > 0.25 → 非常大分化
参考资料：Weir & Cockerham 1984、VCFtools文档、PLINK 2.0文档、Yi et al. 2010 PBS方法