群体重测序分析入门：从 FASTQ 到选择清除的完整地图#

做惯了转录组和基因家族挖掘的人，第一次面对群体重测序（population whole-genome resequencing, WGS）项目时容易懵：同样是”测序数据下机”，工具链、统计量和思维方式几乎完全换了一套。这篇是我给自己整理的路线图——把 WGS 群体分析从原始数据到选择清除的完整链条摊开，讲清楚每一步数据长什么样、结果怎么看、坑在哪里，并附上可以直接下载练手的公开数据集和值得逐篇读的论文。我正在用拟南芥 1001 Genomes 的子集实练这条链，跑通之后再写实战篇。

先画一张总图：

FASTQ ──质控──► 比对(BAM) ──► 变异检测(VCF) ──► 过滤(干净 VCF)
                                                    │
        ┌───────────────┬───────────────┬───────────┼───────────────┐
        ▼               ▼               ▼           ▼               ▼
   群体结构          遗传多样性        选择清除     LD 衰减          GWAS
  PCA/ADMIXTURE      π / Fst       XP-CLR/XP-EHH   r² 曲线     MLM + 协变量
   (谁是谁)        (差异多大)       (哪被选中)   (衰减多快)     (位点-表型)

plaintext

一、群体重测序在回答什么问题#

先把三个容易混淆的概念分开：

分析类型	典型问题	数据形态
转录组（RNA-seq）	哪些基因表达变了	读数映射到基因，数表达量
候选基因挖掘 / 比较基因组	某个性状可能由哪些基因控制	少数个体，基因家族/共线性/结构证据链
群体遗传 / WGS	群体的历史与选择压力：驯化来自哪里、哪些区域被人工选择盯上、哪些位点控制表型	大群体（几十到几千个体）的全基因组 SNP

群体分析的本质，是把几百个个体在同一坐标系（参考基因组）上的差异位点收集起来，然后从等位基因频率的空间分布里读出历史：群体分几支（结构）、彼此差多远（多样性）、哪些基因组区域在某个亚群里被”压扁”了多样性（选择清除的信号）、哪些位点的基因型和表型共变（GWAS）。

二、数据长什么样：FASTQ → BAM → VCF 逐层看#

2.1 FASTQ：原始读段#

每条读段四行：序列 ID、碱基序列、”+“、质量字符串。质量字符的 Phred 值减 33（Illumina 1.8+ 编码）就是每个碱基的错误概率，Q20 = 1% 错误率，Q30 = 0.1%。

看数据看什么：fastqc 出报告后重点看三格——Per base sequence quality（箱体何时跌破 Q30，决定修剪长度）、Overrepresented sequences（接头污染）、GC content（异常尖峰提示污染）。多样本用 multiqc 汇总，一眼看出哪个样本是离群点。

2.2 BAM：比对结果#

bwa-mem2 mem 比对 → samtools sort 排序 → samtools markdup 去重（PCR 重复不标记，后面的 GATK 会把扩增伪影当等位基因）。

看数据看什么：

samtools flagstat sample.bam | less -S   # 比对率、properly paired 比例
samtools depth -a sample.bam | awk '{s+=$3} END{print s/NR}'  # 平均深度

bash

经验口径：植物群体重测序样本 10× 以上、比对率 90%+、properly paired 90%+ 属于健康区间；某个样本深度只有 3×，后面变异检测会大量缺失，宁可重测或剔除。

2.3 VCF：变异位点——群体分析真正的原料#

这是最值得花时间学会读的文件。核心列：

CHROM  POS  REF  ALT  QUAL  FILTER  INFO  FORMAT  样本1  样本2 ...

plaintext

CHROM/POS/REF/ALT：参考基因组上哪个位置、参考碱基是什么、变成了什么。
FORMAT 里的关键字段：GT（基因型，0/0 纯合参考、0/1 杂合、1/1 纯合变异、./. 缺失）、DP（该位点覆盖读段数）、GQ（基因型质量）。
判断一个位点靠不靠谱：QUAL 高、FILTER 为 PASS，且多数样本 DP 在 5–100 之间、GQ ≥ 20。一堆 ./. 的位点在群体里没有信息量。

过滤的两条路线：

VQSR（GATK 变异质量重校准）：需要大量已知可靠位点做训练集，适合有人类、拟南芥、水稻这类注释成熟物种的大项目。
硬过滤：无训练集时的标准做法（多数非模式物种走这条）。GATK 官方 SNP 推荐阈值：QD < 2.0、FS > 60.0、MQ < 40.0、SOR > 3.0、ReadPosRankSum < -8.0。之后还要做群体级过滤——缺失率（--max-missing 0.8）和最小等位基因频率（--maf 0.05）：MAF 过低的位点在群体统计里全是噪声。

工具选择上的一句话答案：有参考基因组、样本量不大时 GATK HaplotypeCaller（gVCF 流程）是金标准；样本多、赶进度或做初步筛选时 bcftools mpileup | bcftools call 更快，两者结果主体一致、边界位点有出入——重要的是整个项目用同一条管线，别混。

三、全流程地图：每步做什么、结果怎么看#

3.1 群体结构：PCA、ADMIXTURE、系统发育树#

PCA（plink --pca 或 EIGENSOFT smartpca）：把几千个样本 × 百万 SNP 压到二维。怎么看：PC1 通常是最大的分层（物种/生态型/地理），如果 PC1 把样品分成了栽培种和野生种两团，后续所有统计都必须考虑这个分层，否则全是混杂效应。
ADMIXTURE：每个个体按 K 个祖先成分画堆积条形图。K 怎么定：跑 K=2 到 K=10，取交叉验证误差（CV error）最低点；同时看生物学合理性——CV 最低是 K=7 但 K=3 就对应了三个地理亚群，通常 K=3 是更可解释的答案，两者都要报告。
系统发育树（IQ-TREE，用 SNP 或全基因组联配）：看拓扑是否与 PCA/ADMIXTURE 一致。三件事互相印证，群体结构才算立住。

3.2 遗传多样性：π 与 Fst#

** nucleotide diversity (π)**：群体内平均每位点杂合度。驯化群体的 π 显著低于野生群体——瓶颈效应的直接证据。
Fst（Weir & Cockerham 1984）：群体间分化度，0 到 1。看的是滑窗值（比如 50 kb 窗口、10 kb 步长），全基因组 Fst 是个没有定位能力的平均数，有意义的永远是某段窗口显著高于背景。
工具：vcftools --window-pi / --weir-fst-pop，或 VCFtools 的现代替代 pixy（处理缺失数据的偏差更小）。

3.3 选择清除：驯化和改良的指纹#

人工选择会让受选择区域多样性骤降、等位基因频率被推向固定。三类互补信号：

Fst / πratio（组间比较）：栽培 vs 野生的高 Fst 窗口。
XP-CLR（Chen et al. 2010）：多等位基因连锁 haplotype 模式，对软选择也敏感。
XP-EHH（Sabeti et al. 2007）：长单倍型纯合的扩展，对接近固定的硬选择敏感。

实操要点：窗口与步长的选择决定分辨率和噪声的平衡（常用 10–50 kb 窗、1/10 窗长为步长）；单指标峰值不算证据，至少两个指标重叠 + 基因注释落到已知驯化基因附近（比如水稻 sh4、progs1 级别的验证逻辑），才是能写进论文的选择区间。

3.4 LD 衰减#

位点间连锁不平衡（r²）随距离衰减，衰减距离决定 GWAS 需要的标记密度和关联信号的定位精度。工具 PopLDdecay 或 plink --r2。怎么看：r² 降到背景值一半时的物理距离就是衰减距离；自交物种（拟南芥、水稻粳稻）衰减慢，异交物种衰减快——荞麦恰好处在异交端，意味着 GWAS 需要更密的标记，这也是荞麦 GWAS 相对少、BSA/QTL-seq 相对多的原因之一。

3.5 GWAS：基因型和表型的全基因组关联#

模型：GLM（不考虑结构）基本会被混杂效应淹没，标准做法是 MLM——群体结构（PCA 前几维或 ADMIXTURE 的 Q 矩阵）作固定效应协变量 + 亲缘关系矩阵（K）控制家系（Yu et al. 2006 统一模型）。工具：GEMMA / EMAXX / GAPIT(R)。
结果怎么看：曼哈顿图上显著峰（Bonferroni 或 FDR 校正线以上）+ QQ 图看膨胀因子——观测 P 值整体偏离期望线说明结构没控干净，先回头修协变量，别急着报位点。
表型要多年多点、至少有重复，单点单年的 GWAS 关联出来一半是环境噪声。

四、示例数据集：从哪里开始练手#

数据集	物种	规模	获取	适合练什么
1001 Genomes	拟南芥	1135 个体	1001genomes.org，可直接下载现成 VCF	最佳入门：注释成熟，可取 50–100 个体子集练 PCA/ADMIXTURE/π，下载量可控
水稻 3K	栽培稻	3024 份	SNP-Seek 数据库（irri.org）	大群体结构分析、籼粳分化、驯化信号复现
1000 Genomes	人	2504 个体	官网/ENA	群体分析方法的原产地，教程最丰富
DPGP	果蝇	非洲群体	DPGP 官网	野生群体、高重组率的对照练习
荞麦重测序	苦荞/甜荞	510–572 份	见 BuckwheatGPDB（buckwheat-gpdb.cn），收录荞麦基因组、基因型、表型与 GWAS 数据集	我的课题方向：荞麦群体结构、驯化与性状定位的第一手参考

建议路径：先用 1001 Genomes 的 50 个体子集把 PCA → ADMIXTURE → π → Fst 跑通（一台普通工作站一两天），再换荞麦数据复现一篇论文的图，最后才碰 GWAS。

五、论文清单：按阅读顺序#

第一批：经典应用（先看故事怎么讲圆）

Huang et al., 2010, Nature Genetics —— 水稻 517 份地方种 WGS，群体结构与驯化的模板论文。
Huang et al., 2012, Nature —— 950 份水稻重测序 + GWAS 定位驯化基因，“群体分析讲故事”的天花板。
1001 Genomes Consortium, 2016, Cell —— 1135 份拟南芥全球多态性图谱，看图学作图。
The 3,000 Rice Genomes Project, 2014, GigaScience —— 大群体资源型论文怎么写。

第二批：荞麦专题（我的方向）

Zhang et al., 2021, Plant Biotechnology Journal —— 苦荞 510 份全球种质重测序，GWAS 定位农艺性状，含多重起源证据。
Zhang et al., 2023, Molecular Plant —— 甜荞 572 份泛变异图谱，比较基因组 + 群体遗传结合。
Li et al., 2025, Cell Reports —— 苦荞完整参考基因组 + 突变群体重测序。
Advanced Science（2025）—— 苦荞离子组 GWAS，DOI: 10.1002/advs.202412291。

第三批：方法溯源（按工具挂）

工具/统计量	原文
PCA（EIGENSOFT）	Patterson, Price & Reich, 2006, PLoS Genetics
ADMIXTURE	Alexander, Novembre & Lange, 2009, Genome Biology
Fst	Weir & Cockerham, 1984, Evolution
XP-CLR	Chen, Patterson & Reich, 2010, PLoS Genetics
XP-EHH / iHS	Sabeti et al., 2007, Nature
GWAS 统一 MLM	Yu et al., 2006, Nature Genetics
GEMMA	Zhou & Stephens, 2012, Nature Genetics

读法：方法原文不需要精读推导，把”统计量在度量什么、什么条件下会失效”两件事读明白就够用了。

六、自检清单：能不能独立扛一个 WGS 群体项目#

以下每条都要能用一句话回答，答不上来就回到对应章节：

变异检测选 GATK 还是 bcftools，依据是什么？
硬过滤的参数为什么是那几个？VQSR 什么时候不可用？
ADMIXTURE 的 K 值怎么定？CV error 和生物学解释冲突时怎么办？
Fst 为什么要做滑窗？窗口和步长怎么权衡？
XP-CLR、XP-EHH、Fst/πratio 各自对什么类型的选择敏感？
LD 衰减距离怎么读？它如何决定 GWAS 的标记密度和定位精度？
GWAS 的 MLM 里，群体结构和亲缘关系分别控什么？QQ 图膨胀说明什么？

七条全能答上来，入职第一个月接手一个 WGS 群体项目就不会卡住；能答五条，说明还处在”候选基因挖掘”层，把第一、三批论文读完再自测一遍。

写在最后#

这条工具链和转录组是两套完全不同的肌肉：转录组的核心是”条件间比较”，群体遗传的核心是”频率与时间”。我给自己的实练计划是从 1001 Genomes 五十个体子集开始，两周内把 PCA、ADMIXTURE、π、Fst 跑通并画出可发表级别的图，然后切到荞麦数据复现 Zhang 2021 的群体结构图。跑通之后，实战篇会把每一步的命令、参数和踩过的坑完整记录下来——和这篇地图放在一起，就是一条从入门到能扛项目的完整路径。