企业官网建设流程全解析

总有刚入门的同学拿着一份VCF来找我，开口就是：“帮我跑个Fst和Tajima's D，两个群体的SNP我都提出来了。”听到这话我就知道，后面至少埋着五六个坑。不是说这些统计量难算——现成工具一大堆，半小时就能出一堆数字——难的是你根本不知道这些数字在回答什么问题，更不知道哪些数字是从有问题的输入数据里硬算出来的。

群体基因组学入门，绕不开三样东西：等位基因频谱、群体分化、选择信号检测。这三者其实是一套组合拳：等位基因频谱是群体内部变异的基础人口普查，Fst这类指标量化群体之间的分歧程度，选择信号检测则是在多样性异常和分化异常的地方找生物学故事。这篇文章不堆教科书定义，我直接从实际分析的角度，把三个概念怎么算、怎么读、怎么配着用，以及我踩过的坑一次说清楚。适合刚跑通比对和变异检测、准备正式进入群体遗传分析的读者。

1. 等位基因频谱、群体分化、选择信号：先理清这一套逻辑

1.1 SFS是群体内部的“变异人口普查”

等位基因频谱（Site Frequency Spectrum，SFS）做的事情很简单：从某个群体里取样n条同源染色体，对每个双等位基因位点统计“衍生等位基因出现了几次”，然后按出现次数从1到n-1做一个频数计数。听起来像数人头，但它其实是群体内部变异结构最浓缩的一张快照。

这张快照之所以重要，是因为群体历史的几乎每个事件都会在它的形状上留下痕迹。群体扩张时突变累积加快但还没有被漂变淘汰，频谱左边（稀有变异）会异常高；群体经历瓶颈时低频变异大量丢失，频谱中间和右边会相对上翘。选择也一样：一个刚完成扫荡的位点附近，稀有变异往往还没恢复，频谱整体左偏。所以读懂SFS，是后面判断选择信号的基础功。

1.2 Fst和dxy量化的是“群体之间”的分歧

如果说SFS看的是“一个群体内部有多少差异”，那群体分化看的就是“两个群体的组成到底隔了多远”。最常用的相对指标是Fst，它的本质是比较群体间的等位基因频率差异占总的遗传变异多少比例。0.01说明两个群体几乎不分家，0.2以上在自然群体里已经是相当可观的隔离程度。

但这里有个很多人一开始没意识到的点：Fst是相对差异，不是绝对差异。一个多样性极低的群体和一个多样性极高的群体，Fst的数值会受各自内部多样性牵制。所以现在分析里普遍把dxy（绝对分歧）也带上。dxy是直接数两个群体之间任意两条序列的平均差异，不关心“相对”二字。两个指标配合，才能判断一块区域的高Fst到底是真分化，还是仅仅因为多样性被压低了。

1.3 选择信号检测是在找异常，但异常不等于证据

选择信号检测的思路本质上是一个“异常排查”的思路：在滑窗里计算各种统计量，然后找出分布尾部的极端窗口，再去看这些窗口对应的基因有没有故事。常用的统计量包括核苷酸多样性π、Tajima's D、Fst、群体间分枝统计量PBS，以及需要单倍体相位信息的iHS、XP-EHH这类LD相关统计量。

需要提醒的是，选择信号检测输出的永远是“候选”，不是“结论”。重测序数据本身带不进表型，一个窗口出现极端Fst，可能是局部适应，也可能是参考基因组拼接错误、拷贝数变异、背景选择，甚至只是样本量不均衡造成的假象。所以这一篇里我会反复强调一件事：多个独立统计量交叉验证，再加功能注释，才有可能把候选缩小到值得花钱做实验验证的位点。

2. 等位基因频谱：从读懂一张图开始

2.1 展开频谱和折叠频谱，差的是一个祖先状态

SFS严格定义里统计的是“衍生等位基因”的频率，也就是突变产生以后的那个新等位基因。可问题是，序列数据本身不会告诉你哪个等位基因是祖先状态，哪个是突变来的。这个信息通常要依赖一个亲缘关系足够远的外群物种：如果外群和群体里的一条序列一致，那这个等位基因大概率是祖先型；不一致，就是衍生型。

有外群、能确定祖先状态的时候，得到的是“展开频谱”（unfolded SFS），从频率1到n-1每个bin都独立有意义。没有可靠外群，或者外群分歧太远导致测序比对偏差很大，那就只能把频率为i和n-i的bin折叠起来，得到“折叠频谱”（folded SFS）。折叠后信息量会损失，例如Tajima's D这类依赖衍生等位基因方向的统计量就没法直接从折叠频谱里推，只能另算。

我个人的建议是：能做展开就不要用折叠。因为折叠相当于把两个bin合并，等于把高频和低频的信息强行混在一起，很多细微的历史信号会被抹掉。但前提是你手上确实有一个高质量外群。如果外群比对率低、覆盖度参差不齐，那折叠频谱反而更稳。

2.2 低深度数据算SFS的正确姿势

很多人习惯先做完变异检测拿到VCF，再从VCF里数等位基因频率。这个流程在测序深度足够高（比如30x以上）时问题不大，但重测序项目为了控制成本，动辄只有5x到10x的深度。这种数据里，单个位点的基因型本身就不确定，硬call一个基因型Ignores不确定性，低频变异——恰恰是SFS最看重的部分——会被成批漏掉。

正确的做法是用基因型似然（genotype likelihood）框架，在每一个位点上保留四种可能的似然值，而不是先做一个稳健的基因型判断。这个思路最有代表性的实现就是ANGSD。流程大致是：先对每个群体分别生成位点等位基因频率的似然文件（.saf），再用realSFS做最大似然估计得到SFS。整个过程不产出VCF，却能直接用上所有reads的信息，这是低深度数据算SFS最稳的方案。

2.3 一张SFS图该怎么解读

拿到SFS之后，别急着套软件推断历史，先自己看图。横轴是衍生等位基因频率类别（或者折叠后的最小频率），纵轴是位点数量或者占比。我一般会先看三个位置：

• 左边第一个bin（singleton）占比高不高。如果单拷贝变异占了相当大比例，往往是群体扩张、近期选择扫荡或者较强的纯化选择。当然，测序错误也会堆在左边，所以低深度数据要先确认错误率控制得住。
• 中间频率的bin有没有异常上翘。中频变异过多，常见的解释是平衡选择（比如与病原抗性相关的多态长期维持），也可能是分析时混入了两个有结构、但没有被正确区分的亚群——这个混淆经常被忽略。
• 最右端高频bin是不是接近零。折叠频谱里最右端对应的是频率接近0.5的位点，这部分如果异常高，就要怀疑参考基因组偏向或者样本间亲缘关系过近。

看完形状，再用SFS做正式推断才有底。比如用∂a∂i或者fastsimcoal这类工具去拟合并估计种群历史参数时，SFS的bin值会被直接当作多项式分布的观测值，任何方向性偏差最终都会变成错误的历史参数。

3. 群体分化：Fst只是半个故事

3.1 不同软件算出来的Fst其实不是同一个东西

Fst这个名字背后有一堆估算公式，实践里最常碰到的有这么几类：Wright最早定义的基于固定指数的写法，Weir和Cockerham基于方差分量做的无偏估计，以及Hudson提出的基于等位基因频率的估计量。vcftools里默认的--weir-fst-pop走的是Weir-Cockerham路线，ANGSD的realSFS fst在低深度框架下走的是类似Hudson的思路，pixy这种带似然框架的工具又有自己的实现。

这些公式在数据质量好、样本量均衡时数值差异不大，但千万不要拿两个软件的结果直接做数值对比。我见过有人把vcftools算出来一个窗口的Fst是0.35，pixy算出来同样窗口是0.28，然后就断定某个工具出了问题——其实两边定义不同，可比的是趋势而不是绝对值。正确的习惯是：同一个项目从头到尾只用一套公式，报告里写明用了哪个工具、哪个参数，别人才能复现。

3.2 单点Fst不靠谱，窗口化才是标准操作

很多刚接触Fst的同学喜欢直接看单个SNP的Fst，找那些Fst接近1的位点。这个做法在逻辑上就走偏了。单个双等位基因位点的Fst在样本量不足时方差极大，尤其是频率极低或极高的位点，一个样本的分组差异就能让单点Fst从0.02跳到0.9。而且真实选择信号往往是作用在一段连锁区域上的，单点看等于管中窥豹。

标准做法是滑窗：把基因组按固定长度（比如100kb）切成窗口，每个窗口内对位点级统计量做平均或者加权合并，然后沿染色体滑动（比如步长20kb或50kb）。窗口大小是个权衡——小了噪声压不住，大了把真正的信号稀释掉。我自己一般在全基因组范围先用100kb窗口扫一遍出候选区域，再对候选区域用20kb甚至10kb的窗口细看信号轮廓。窗口的边界效应也值得留意，信号峰值落在两个窗口交界处时，建议把窗口起点平移半个窗口再跑一次。

3.3 dxy和Fst一起看，才能避免被多样性误导

一个实际的例子：两个群体在某一段区域都没有什么多态，所有位点都是同一种纯合状态，这时候Fst会因为缺乏变异而在数值上变得不稳定，甚至出现很高的值。看起来像是极端分化，但其实两个群体在这块区域几乎没有差异——这就是相对指标的天生局限。

在这个场景下，dxy的价值就出来了。dxy直接数两个群体间的平均每位点碱基差异，它不受群体内部多样性的牵制。判断规则我习惯用“两维坐标”：Fst高但dxy也高的区域，才是真正的分化热点；Fst高但dxy低的区域，多半是多样性塌陷造成的伪信号，需要去看是否对应选择性清除、拷贝数缺失或比对死角。反过来，dxy高但Fst不高的区域，说明两边各自积累了不同的变异，只是因为各自内部多样性太高把相对值拉平了，这种情况在分化时间不长的姐妹群体里很常见，单独看Fst很容易漏掉。

4. 选择信号检测：把统计量组合成证据链

4.1 不同统计量看到的是不同时间尺度的选择

选择信号检测从来不是“跑一个软件然后数峰”的游戏。选择发生的时间尺度不同，留下痕迹的统计量也不同。一个比较实用的分类方式是把检测策略按信号类型分成几组：

• 基于多样性降低：选择性清除会让局部区域核苷酸多样性π显著低于基因组背景，代表策略是π比值、CLR类方法。
• 基于频谱偏移：选择扫荡后低频变异尚未恢复，Tajima's D显著负值；平衡选择则让D显著正值。
• 基于群体分歧：局部适应会让特定群体在某个区域的等位基因频率与其他群体明显不同，代表是Fst窗口、PBS。
• 基于单倍型结构：近期选择扫荡造成长单倍型延长，代表是iHS、nSL、XP-EHH，这些方法通常需要先做单倍型相位分析。

一张表可以快速评估自己手头的数据适合用什么：

如果你的项目数据连单倍型相位都还没做，那iHS这类方法暂时不适用，老老实实从π、Tajima's D和Fst这些基础量入手，先建立证据链的下限。

4.2 滑窗参数和异常值阈值怎么定

滑窗参数的选择直接影响你能发现多窄的信号。全基因组扫描常用50kb到100kb窗口，步长取窗口的1/5到1/2。窗口太小时估计噪声大，太大时一个窗口可能跨越多个重组区段，把选择信号周围的背景多样性都平均进来。我一般会跑两个尺度的窗口，一个大窗口（100kb/20kb步长）做全基因组定位，一个小窗口（20kb/10kb步长）在候选区域里精确定位。

阈值方面，不要迷信某个固定数值，比如“Fst>0.4就算候选”。基因组不同区域的突变率、重组率、基因密度差异巨大，固定阈值会把全基因组大约等比例的窗口都标记成异常，没有任何统计意义。更合理的做法是取经验分布的尾部，比如每个统计量在全基因组窗口分布的前5%或后5%作为异常值候选，然后要求多个指标同时异常才进入下一步。比如我常用的一组标准是：Fst处于该染色体前5%、π比值（两个群体π的比值）处于后5%、Tajima's D显著负值，三者同时满足的窗口才列为强候选。

4.3 从候选峰到功能基因的最后一公里

统计上拿到一批候选窗口之后，下一步才是真正体现分析深度的地方。先把窗口和参考基因组的注释文件做overlap，找到落在窗口内或其上下游一定距离的基因。距离设置别太死板，选择作用可能落在调控区域，基因本体并不一定直接卡在信号峰正下方，上下游50kb以内的顺式调控区域都要考虑。

然后对窗口内的变异位点做功能效应预测，重点看非同义突变、剪接位点突变和启动子区域变异。若群体间存在明显频率差异的能力位点，再结合基因的已知功能去推测适应性场景。如果候选基因一批，可以做GO或通路富集，但样本量小的时候富集结果很可能只是噪声，我更倾向于逐一查阅候选基因的文献，确认它们是否在相应生理过程里有已知功能。

最后要记住：这一切都只是候选证据。真正的验证需要回到湿实验，比如功能实验、群体重测序的独立样本验证、或者至少用另一个不重叠的样本集重新扫描一遍，确认信号不是样本特异的。

5. 一套可以直接照抄的实战流程

5.1 跑之前先检查五件事

在敲任何统计命令之前，先花半小时检查输入数据。我踩过太多次“跑完才发现输入有问题”的坑了。清单如下：

1. BAM文件的染色体命名和参考基因组完全一致，别出现scaffold和chr的混用。
2. 每个样本的测序深度和覆盖度是否均衡，极端偏差样本最好直接剔除。
3. 亲缘关系：先做一遍个体间亲缘关系或主成分分析，得过近的样本（比如意外重复测序）会让群体统计量失真。
4. 参考基因组里重复区域和拼接gap是否做好mask，这直接关系到dxy和Fst的假信号。
5. 明确自己的分析单位：常染色体、性染色体、线粒体要分开处理。

样本量上，我不建议任何群体少于5个个体。少于这个数，SFS的低频bin基本是空壳，Fst窗口的方差会大到失去意义。自然群体和实验群体要求也不同，这里只说一般底线。

5.2 低覆盖度数据用ANGSD流程

假设你的BAM文件已经按群体列好，每个群体一个list文件，参考基因组和外群物种都准备好了。计算展开SFS的命令大致是这样：

# 对每个群体生成SAF文件 angsd -bam popA.bamlist -ref genome.fasta -anc outgroup.fasta \ -out popA -doSaf 1 -GL 1 -minMapQ 30 -minQ 20 \ -uniqueOnly 1 -remove_bads 1 -only_proper_pairs 1 \ -minInd 5 -maxDepth 200 -baq 1 # 得到一维SFS realSFS popA.saf.idx > popA.sfs # 两个群体的二维SFS realSFS popA.saf.idx popB.saf.idx > popA.popB.ml # 用二维SFS做窗口化Fst realSFS fst index popA.saf.idx popB.saf.idx \ -sfs popA.popB.ml -fstout popAB realSFS fst stats2 popAB.fst.idx -win 100000 -step 20000 > popAB.fst.txt

如果确实没有可靠外群，去掉-anc参数，后面realSFS加-fold 1生成折叠频谱即可。这里-minInd按你的样本量调整，低于这个数的位点会被当作缺失处理。

5.3 高覆盖度数据用VCF流程

覆盖度足够（建议30x以上），可以直接用VCF，流程会直观不少。vcftools的经典组合：

# 每个群体的样本名单，一行一个样本名 # Fst窗口 vcftools --vcf all.vcf.gz \ --weir-fst-pop popA.txt --weir-fst-pop popB.txt \ --fst-window-size 100000 --fst-window-step 20000 \ --out popAB # 每个群体的pi和Tajima's D vcftools --vcf all.vcf.gz --keep popA.txt \ --window-pi 100000 --window-pi-step 20000 --out popA_pi vcftools --vcf all.vcf.gz --keep popA.txt \ --TajimaD 100000 --out popA_td

dxy和更标准的Fst建议用pixy一起跑，它按位点处理基因型似然并输出pi、Fst和dxy三种统计量，可以避免用不同工具算出来的口径差异：

# populations文件：两列，样本名和群体名 pixy --stats pi fst dxy \ --vcf all.vcf.gz --populations populations.txt \ --window_size 100000 --n_cores 8 \ --bypass_invariant_check yes

pixy的--bypass_invariant_check建议打开，否则它只算多态位点，dxy会系统偏低。输出文件里的fst_window.txt和dxy_window.txt直接就能用于后续阈值筛选。

5.4 选择信号扫描的收尾

有了上面的pi、Fst、Tajima's D结果，就可以做交集筛选。我习惯把这些窗口结果按染色体分组，分别算百分位，避免用全基因组统一阈值掩盖染色体间异质性。一个窗口如果同时满足Fst在该染色体前5%、π比值后5%、Tajima's D为负值，就进入候选清单。

候选窗口清单出来后，用bedtools直接和基因注释文件做交集：

bedtools intersect -a candidates.bed -b genes.bed -wa -wb > candidate_genes.txt

然后再对候选基因做变异注释和文献排查。整个流程到这里就算闭环了，后面是湿实验的事，不在本文范围。

6. 我踩过的坑，逐个说给你听

6.1 MAF过滤之后算SFS，等于自己骗自己

这是新人最容易犯、也最伤的一个错误。为了“去掉测序错误”，很多人习惯先做MAF过滤（比如只保留minor allele frequency大于0.05的位点），然后再去画SFS。问题是，MAF过滤相当于把频谱最左边的若干bin全部清零，SFS的形状直接被削掉半边。后续任何依赖SFS的推断（Tajima's D、种群历史拟合）都会输出一个“假扩张”或“假选择”的结果。

处理测序错误的正确方式不是事后按频率过滤，而是在源头控制：比对质量、碱基质量、错误率模型、BAQ，这些参数该设的设到位。ANGSD的GL框架本身就包含错误率建模，不需要额外砍低频位点。

6.2 低深度硬call基因型，结果基本是噪声

低深度数据强行用GATK或者bcftools的genotype模块硬call一个基因型，然后拿这个VCF去算dxy和Fst，这个流程在5x深度下能保留多少真实信号，取决于运气。低频杂合位点最容易漏call，而漏掉杂合位点会直接抬高群体间差异的假象——因为零散的高频纯合错误看起来就像群体分化。

如果你手里的数据平均深度低于15x，我的建议是默认使用各级基因型似然方法（ANGSD、pixy二选一）。虽然学习成本高一点，但结果质量完全不是一个级别。

6.3 样本量不齐和缺失模式不同

样本量不齐是群体比较里的隐形杀手。一个群体10个样本，另一个群体3个样本，Fst对稀有等位基因的响应完全不同，窗口方差也不一样。更隐蔽的是缺失模式：如果A群体在某个区域普遍深度低、缺失率高，B群体深度正常，那这个区域的等位基因频率本身就已经被测序深度变相筛选过，算出来的Fst和π都有偏。

尽量让群体样本量接近，缺失标准也按群体分别检查。不要只用全数据集的minInd阈值，那样大样本群体可能没什么位点被过滤，小样本群体却被过滤掉一大半。

6.4 基因组里的“暗物质”会冒充选择信号

重复区域、着丝粒、拼接塌陷区域，这些地方比对质量天然差，reads堆出来一堆假多态，或者反过来几乎没reads覆盖。这两种情况都会制造虚假的极端统计量：覆盖塌陷区域产生的高dxy和高Fst尤其逼真，从数值看就像受了强选择。所以跑窗口之前一定要用repeat mask和低复杂度区域注释把基因组“暗物质”遮掉，或在最后删除落在这些区域的候选窗口。

另外，串联重复和拷贝数变异所在的窗口，我一般直接标灰不看。偶发拷贝数差异在个体间产生的read深度波动，会伪装成等位基因频率差异，pixy和vcftools都救不了这种数据。

6.5 峰和基因是两码事

拿到一个漂亮的选择信号峰，距离“这个基因在适应中起作用”还有很长的路。峰只告诉你这个区域有统计异常，异常的来源可能是选择、背景选择、群体结构、参考基因组缺陷，甚至只是统计方法的边界效应。候选基因列表里排在第一位的不代表就是因果基因，很可能只是它恰好落在窗口正中央。

我自己的做法是：把峰内所有高影响的变异位点列出来，看它们在两个群体间的频率差，再看注释类型，最后才结合文献判断哪些位点值得进入功能验证。这一步任何软件都替代不了人工判断。好工具给出的是候选方向，不是答案。

说到底，把我这几年和这些数字打交道的经验浓缩成一句话：每个统计量都是一块拼图，单独拎出来都可能骗人，但几个独立指标合在一起指向同一个方向的时候，结论往往就稳了。希望这篇能帮你少踩几个我踩过的坑。