从Fst到Tajima's D：WGS重测序数据群体遗传分析实战指南
当面对海量的全基因组重测序数据时，如何从中提取出有价值的群体遗传信息？本文将带你从VCF文件开始，逐步完成Fst、Pi和Tajima's D等关键参数的计算与分析。不同于理论教科书，我们聚焦于实际代码操作与结果生物学解读，适合刚接触群体遗传分析的生物信息学研究者。
1. 分析前的准备工作
在开始计算前，我们需要确保环境配置正确并理解数据的基本特征。群体遗传分析通常从VCF文件开始，这是存储样本基因型变异信息的标准格式。
首先检查vcftools是否安装：
vcftools --version若未安装，可通过conda快速配置环境：
conda create -n popgen -c bioconda vcftools plink r-tidyverse conda activate popgen对原始VCF文件进行质量检查是必不可少的步骤：
vcftools --vcf input.vcf --check-indels --out quality_check注意：在实际分析中，我们常遇到以下VCF文件问题：
- 样本ID重复或格式不规范
- 染色体命名不一致（如chr1 vs 1）
- 缺失基因型比例过高
提示：使用
bcftools stats input.vcf > vcf_stats.txt可生成详细的变异统计报告，帮助评估数据质量
2. 群体遗传分化分析：Fst计算与解读
Fst是衡量群体间遗传分化的核心指标，其值域为0-1，越接近1表示分化程度越高。我们使用vcftools计算滑动窗口Fst：
vcftools --vcf input.vcf --weir-fst-pop pop1.txt --weir-fst-pop pop2.txt --fst-window-size 50000 --fst-window-step 10000 --out pop1_vs_pop2关键参数说明：
--fst-window-size：窗口大小（建议参考连锁不平衡衰减距离）--fst-window-step：滑动步长（通常为窗口大小的1/5-1/2）
计算结果包含三列：
CHROM BIN_START BIN_END N_VARIANTS WEIGHTED_FST MEAN_FSTFst值生物学解释指南：
| Fst范围 | 分化程度 | 可能生物学意义 |
|---|---|---|
| 0-0.05 | 低 | 频繁基因流或近期分化 |
| 0.05-0.15 | 中等 | 适度隔离或生态分化 |
| >0.15 | 高 | 长期隔离或强选择压力 |
在R中可视化Fst分布：
library(ggplot2) fst_data <- read.table("pop1_vs_pop2.windowed.weir.fst", header=T) ggplot(fst_data, aes(x=MEAN_FST)) + geom_histogram(bins=30, fill="steelblue") + labs(x="Fst值", y="窗口数量", title="群体间Fst值分布")3. 遗传多样性分析：π与θ的计算
核苷酸多样性π反映群体内遗传变异水平，计算公式为：
π = Σ(2n/(n-1)) * p(1-p)其中n为样本数，p为等位基因频率。使用vcftools计算：
vcftools --vcf input.vcf --window-pi 50000 --window-pi-step 10000 --out pop_piθ是另一个多样性指标，基于分离位点数目估算：
vcftools --vcf input.vcf --site-pi --out pop_theta常见问题解答：
- 为何π和θ值有时差异很大？
- π对低频变异更敏感
- θ受测序深度影响更大
- 水稻群体的典型π值范围：0.001-0.005
- 果蝇群体的典型π值范围：0.01-0.02
4. 中性检验：Tajima's D的计算与解释
Tajima's D检验群体是否符合中性进化模型，计算公式为：
D = (π - θ)/√(Var(π - θ))计算命令：
vcftools --vcf input.vcf --TajimaD 50000 --out pop_tajima结果解读框架：
| Tajima's D值 | 可能进化场景 | 选择类型暗示 |
|---|---|---|
| 显著>0 | 群体收缩或平衡选择 | 维持多态性 |
| ≈0 | 符合中性模型 | 无显著选择 |
| 显著<0 | 群体扩张或正选择 | 定向选择 |
在R中绘制曼哈顿图：
tajima <- read.table("pop_tajima.Tajima.D", header=T) ggplot(tajima, aes(x=BIN_START/1e6, y=TajimaD, color=CHROM)) + geom_point(alpha=0.6) + facet_grid(.~CHROM, scales="free_x") + labs(x="位置(Mb)", y="Tajima's D") + geom_hline(yintercept=0, linetype="dashed")5. 进阶分析：选择信号检测与结果整合
结合多种统计量可提高选择信号检测的可靠性。以下是常见分析流程：
- 计算XP-CLR（跨群体复合似然比）：
xpclr --input input.vcf --out pop_xpclr --samplesA pop1.txt --samplesB pop2.txt- 整合分析结果：
library(qvalue) fst <- fst_data[,c("CHROM","BIN_START","BIN_END","MEAN_FST")] pi <- read.table("pop_pi.windowed.pi", header=T) tajima <- read.table("pop_tajima.Tajima.D", header=T) merged <- merge(fst, pi, by=c("CHROM","BIN_START","BIN_END")) merged <- merge(merged, tajima, by=c("CHROM","BIN_START","BIN_END")) # 计算综合得分 merged$composite_score <- scale(merged$MEAN_FST) + scale(-log10(merged$PI)) + scale(abs(merged$TajimaD))- 候选基因注释示例：
bedtools intersect -a candidate_regions.bed -b gene_annotation.gff -wo > candidate_genes.txt实战建议：
- 不同物种需调整窗口参数
- 建议结合至少3种方法的结果进行判断
- 关注多个统计量同时异常的基因组区域
