群体基因组学这几年是真的火但不少刚接触的朋友跟我聊的时候总感觉概念都眼熟真拿到自己的数据又不知道从哪里下手。等位基因频谱、群体分化、选择信号检测这三块是群体基因组分析的三大支柱也是从“拿到VCF文件”到“讲出生物学故事”之间最关键的几步。这篇东西我打算按自己做项目时的思路来写把每个概念的来龙去脉、背后的统计逻辑、跑的每一步命令还有踩过的坑都摊开说清楚希望能帮你把这条分析链路完整打通。我先把结论放在这儿这三块内容有很强的递进关系——等位基因频谱是理解数据的基本盘群体分化是描述群体关系的度量衡选择信号检测则是真正回答“哪个基因在进化中很重要”的核心武器。想做好分析缺了任何一环都会让结论站不住脚。接下来我从原理、实操到选型一步步拆开讲。1. 核心思路为什么把SFS、FST和选择检测放在同一条流水线上很多初学者会把这几个概念当成独立的模块分别去学这个方向就错了。真实项目里它们是环环相扣的先从测序数据里拿到最基础的多态位点信息然后通过等位基因频谱观察整体变异规律再用群体分化指标精细刻画亚群之间的关系最后基于这些背景信息去寻找偏离中性的位点。我习惯把这套流程称为“描述—比较—检测”三段式。1.1 一条完整的分析链路是怎样串起来的拿到重测序数据之后第一步是比对、变异检测、过滤最终得到全基因组的VCF文件。这个阶段不考虑生物学阐释只保证数据质量。接下来就要算等位基因频谱也就是把所有多态位点按照“次等位基因出现的频次”分类统计。这样做的目的是理解数据整体长什么样有什么样的测序深度、参考基因组方向是不是正确、有没有明显的缺失模式。紧接着用FST固定指数这一类指标去群体间做比较。为什么要这一步因为任何“选择信号”都是相对于背景来定义的。如果两个群体本身分化就非常剧烈某一个位点显示出的高分化很可能只是群体历史结构造成的假象而非自然选择的作用。最后才做选择信号检测。检测手段本质上是去寻找偏离中性模型预期的异常值。在已知群体分化背景的前提下你可以判断这些异常到底是真的受选择还是群体结构、基因流的自然结果。1.2 频繁踩坑的“视角混淆”问题我发现很多人在做选择信号检测的时候会犯一个经典错误直接在未做群体分层的样本里跑全基因组扫描。假如你有两个生态型存在明显分化却把它们混在一起当一个大群体去计算中性检验那么得到的“信号”往往是群体分层的噪音根本不是选择。正确的做法是“先描述、再比较、后检测”每一步的结果都在为下一步筛选提供背景信息。所以这篇文章的编排逻辑就是按这个流程来的先讲SFS怎么正确生成和解读再讲FST和PBS这类指标的生物学含义与参数选择最后讲选择信号检测中不同方法各自的适用场景以及怎么把多种证据串成一个让人信服的故事。2. 等位基因频谱SFS生成逻辑、中性预期与关键误区等位基因频谱Site Frequency SpectrumSFS是所有群体遗传分析的基石。简单来说它统计的是在一个群体中不同频次等位位点的数量分布。比如有10条染色体出现1次的突变有多少个、出现2次的有多少个、直到出现9次的有多少个。2.1 从VCF到SFS的三种姿势实际项目中生成SFS有三种常见方案按需求选择就好直接基于VCF文件用vcftools统计这种方法适合长期过滤后质量较高的基因型数据速度很快但也容易丢失低频变异的检出能力。基于BAM文件的ANGSDSAMtools流程直接用测序比对结果推算等位基因频谱能充分利用低频信息还能引入测序深度和质量过滤参数。这种方案对测序深度不高的数据特别友好重测序项目我基本都推荐。基于变异注释的衍生SFS如果你想区分同义与非同义位点就需要结合功能注释信息在得到基础SFS后分类汇总。我个人的建议是除非你的数据覆盖度和基因型检出率都很理想否则尽量直接用BAM做。原因很现实VCF里低频变异的基因型叫出率通常不高一个位点如果只有一半个体被成功分型那你的SFS从来源上就已经失真了。2.2 中性模型预期为什么SFS是“倒J型”在标准中性模型下一个样本量为n的二倍体群体中等位基因频谱的预期是每个频率类别的位点数量正比于1/f。这意味着低频变异应该是最多的。如果你的SFS显示出大量高频变异那么就需要警惕是否受到了选择、群体扩张或数据偏误的影响。这也是理解很多检验统计量的基础。比如中性偏离检验D值常称Tajimas D这个经典指标它的逻辑就是比较两种对遗传变异的估计量——基于平均成对差异的π值和基于多态位点数目的Watterson θ值。如果群体正在经历收缩或者说有正选择造成“搭车效应”那么高频变异会比中性预期偏多D值会显著偏离零。我在实际项目中见过很多用户对D值正负号的理解有问题。简单记法D为负说明低频变异过多常见于群体扩张、正选择或数据的低深度个体混入D为正说明高频变异过多常见于群体收缩、平衡选择或混合样本中存在群体结构。具体问题要结合SFS的形状一起判断别单独看一个值。2.3 计算SFS时的“致命细节”生成SFS看似只是跑一条命令实则细节决定成败尤其要注意以下几项极性SFS分“折叠”和“非折叠”两种。非折叠SFS需要指定一个外群或祖先状态来确定哪个是衍生等位基因如果参考基因组跟你的物种演化距离适中可以直接把参考基因组等位基因作为祖先状态。没有外群时只能做折叠SFS相当于把频率w和1-w合并灵敏度会明显下降。测序深度深度过低会导致低频等位基因被当成测序错误过滤掉深度过高又容易引入参考偏差。我一般按平均深度30X左右的数据将最低基因型质量设为20最低位点深度设为8到10。位点过滤顺序做SFS之前过滤做多了会把真实低频变异删掉做少了又会混入大量测序噪声。我摸索出来的经验是先过滤样本与位点质量再计算SFS最后在确认整体形态没问题之后再对SFS做进一步校正。注意SFS是所有后续分析的基础。如果你发现SFS异常前面做得再漂亮的FST和选择扫描结果都不可信因为基础分布已经不对了。3. 群体分化从单点FST到窗口PBS的完整推导群体分化通常是我们要回答的第一个生物学问题不同生态型或地理群体之间是否存在显著的遗传差异哪些区域差异最大。这里最常用的就是FST但市面上的FST版本五花八门算出来差别很大新手经常被搞晕。3.1 FST的本质不复杂但定义要选对FST这个概念可以理解为群体间遗传差异占总体遗传差异的比例。如果所有群体等位基因频率完全相同FST就是0如果各群体固定了完全不同的等位基因FST就接近1。但实际计算时有多种估计方式常见的有以下几种Wright的原始定义需要知道群体真实等位基因频率很少直接用于真实数据。Weir和Cockerham基于方差分析开发的估计量它把群体间和群体内的方差分量做矩估计是目前软件中最常用的算法。Hudson基于等位基因频率简单估计的FST计算速度快常用于滑动窗口扫描。基于成对序列差异的FST估计专门用于测序深度不齐或不依赖基因型时的估算。我自己的经验是做窗口扫描用Hudson的估计方式比较稳定做单点分析、候选位点筛查时可以结合Weir和Cockerham的估计值实际报告中一定要写清楚你用的是什么方法和参数。3.2 窗口化扫描的参数如何决定FST滑动窗口扫描中窗口大小和步长的选择会直接改变信号分辨率。窗口太小单点噪声大窗口太大真实信号被平均稀释。这个参数没有绝对标准但可以按以下逻辑去确定如果重测序序列全长约500Mb左右我通常起步用50kb窗口、10kb步长。如果是芯片数据或低密度标记直接用单点FST并按物理距离绘制可以绕开窗口平均的问题。如果已经有候选区间可以直接在候选区段内用5kb窗口精确定位。关键技巧是“先粗扫、再精扫”。先用大窗口找出信号区间再缩小窗口在区间内确定候选基因范围。如果一步就直接用小窗口全基因组范围的假阳性会多到你怀疑人生。3.3 PBS把三群体故事讲清楚的分支统计量FST可以告诉你哪些区域在群体间分化明显但它说不清楚分化具体发生在哪条演化分支上。这时候就要用分支统计量PBS群体分支统计量。它的逻辑是把三个群体的FST值做-log转换然后通过代数推算将每个群体分支上的独特分化量分离出来。计算方式上三群体的FST两两转换后可以得到一个值衡量目标群体自与其他两个群体分开之后积累的分化量。实际操作时注意三个群体的选择非常关键要按演化关系排好并且尽量选关系清楚的人群组合这样结果的解释才直接。PBS只能做三群体组合如果有四个以上群体可以组合多种三群体方案并取交集。PBS受群体历史影响很大种群瓶颈会让全基因组背景PBS升高所以不能只看绝对PBS值要看相对全基因组分布的百分位。实操时我用R脚本或Python并行计算PBS时有个习惯先生成所有位点的三群体FST再转成两两d值矩阵最后解出各分支的PBS。如果你用vcftools的滑动窗口FST输出记得保留位点位置因为PBS计算需要按窗口对应的区间合并。3.4 当群体历史比较复杂时怎么办现实中很多数据并不满足“三群体独立分化”的理想假设。比如存在广泛基因流、近期群体混合、种群瓶颈等。此时PBS的绝对值会受到严重干扰。我是这么处理的先用主成分分析、结构分析这些手段确认群体遗传结构与演化树基本一致。如果结构分析和PBS设定冲突最好放弃PBS改用基于扩展单倍型纯合度的方法或者用更稳健的群体分化统计量做两两比较。记住任何统计指标的生物学解读都依赖于群体历史模型模型不对指数再漂亮也是自欺欺人。4. 选择信号检测从频谱偏离到单倍型痕迹的综合研判选择信号检测是群体基因组学里最激动人心、也最容易翻车的一步。目前成熟的方法大致分三类基于SFS的方法、基于群体分化的方法、基于单倍型结构的方法。三者检测的目标不同单独用哪一种都有明显局限靠谱的研究基本都靠“多方法取交集”。4.1 基于SFS的检验适合什么场景这类方法包括中性偏离检验D值常称Tajimas D、D常称Fay和Wu的H等。它们通过比较真实SFS与中性预期SFS的偏离来判断是否存在选择。中性偏离检验D值对群体历史很敏感瓶颈和扩张都会造成假阳性所以通常不作为单一证据。对于新近完全扫荡形成的选择信号D值会显著为负但它的分辨率很低很难定位到具体位点。实际操作中我倾向于把这类统计量作为“数据质量与历史事件”的初步探针而不是直接用它出候选基因。如果整个群体的D值全基因组普遍异常那大概率是该群体经历过扩张或收缩应该先去联合历史动态模拟而不是急着找受选择基因。4.2 基于分化的方法要会挑“极端尾部”FST和PBS这类方法在当前数据里用得最广因为它们能直接给出受选择区域的候选位点。这类方法的假设是局部受选择会让两个亚群在该区域的等位基因频率发生显著分化远大于基因组背景的中性漂移。实际操作上我通常会先全基因组扫描获得所有窗口的FST或PBS值再取前1%或5%作为top候选。为什么不能直接设一个绝对值阈值因为不同群体的遗传背景、样本量、标记密度差异很大绝对阈值没有普适性。做选择信号检测时许多流程还会用到交叉群体扩展单倍型纯合度方法常称XP-EHH或群体分化方法常称XP-CLR。这两种方法的原理分别基于单倍型长度和群体分化模式。XP-EHH对仍在传播中的不完全扫荡敏感XP-CLR对已固定的扫荡和长期选择更有效。它们虽然原理不同但共享一个思想相对于背景中性模式受选择区域的等位基因频率分化程度或者单倍型纯合度会有突出的信号。4.3 基于单倍型的检验iHS与XP-EHH怎么选基于单倍型的方法利用的是“受选择位点附近的单倍型往往很长且纯合度高”这一特性。核心指标包括扩展单倍型纯合度EHH、内部单倍型纯合度得分iHS和跨群体扩展单倍型纯合度XP-EHH。iHS在单个群体内检测正在发生或者近期发生的正向选择比较同一核心位点上衍生与祖先等位基因携带的单倍型长度。iHS的绝对值大说明核心位点附近的单倍型异常长。XP-EHH在两个群体之间比较当目标群体的受选择单倍型长度和频率都显著高于另一个群体时信号会很强。它是当前应用最广的跨群体方法之一。以我的项目经验如果要找不同生态型对特定环境适应产生的信号XP-EHH往往比iHS效果好很多因为这里的对比本身就是生态型之间。而同一种群内部的驯化选择信号iHS更有优势。另外XP-EHH对重组率变化的鲁棒性比iHS好但需要两个群体质量都比较高。selscan这类工具可以同时计算和归一化这些指标跑的时候记得要提供遗传图谱或重组率参数别直接拿物理距离硬算。4.4 如何组合证据避免“单指标自嗨”单看任何一项指标都可能被假阳性欺骗多方法交集是最稳妥的做法。我常用的组合策略第一层全基因组XP-CLR或PBS扫描锁定前1%的极端窗口。第二层在候选窗口内计算FST、XP-EHH要求至少两项指标显著。第三层在最终候选区域检查SFS形态和单倍型网络确认不是组装错误或极端缺失导致的假信号。第四层结合功能注释筛选出非同义突变、启动子变异等有生物学意义的候选位点。5. 实战工作流从VCF到候选位点注释的一次性跑通理论部分讲得再多不如给你一套可以直接上手的实操流程。我按自己做真实项目的习惯整理了一份从VCF到候选位点注释的常用命令和参数希望能成为你的起步模板。5.1 数据准备与质控过滤第一步是对VCF做质控。过滤条件基于你的数据质量设定我常用以下参数# 保留双等位位点、过滤缺失过多、低QUAL位点 vcftools --vcf raw.vcf --remove-indels --min-alleles 2 --max-alleles 2 \ --max-missing 0.7 --minQ 30 --recode --out clean这个环节注意两点。一是别一上来就按极端的次要等位基因频率过滤比如只保留5%以上的位点那后面低频变异相关的分析就没法做了。我会在生成SFS之前保留所有多态位点只在做单倍型分析时考虑低频过滤。二是染色体名称和个体ID最好在项目开始时就统一好不然后面每跑一步都在跟命名错误斗智斗勇。5.2 SFS和中性检验的计算如果你有BAM文件且希望更灵敏地捕捉低频变异建议用以下思路生成SFS并计算θ和D值# 使用ANGSD生成等位基因频谱 angsd -bam bamlist.txt -ref genome.fa -out sample \ -doSaf 1 -doThetas 1 -GL 1 -minMapQ 30 -minQ 20 -anc genome.fa # 依据SFS计算核苷酸多样性与D值等指标 realSFS sample.saf.idx sample.sfs如果只有VCF且质量较好直接走vcftools的窗口D值更快vcftools --vcf clean.vcf --TajimaD 50000 --out window_tD注意ANGSD的祖先等位基因设置直接决定你的SFS是否是折叠的。若没有可靠外群建议用折叠SFS以免把参考基因组等位基因误当祖先状态导致衍生等位基因频谱失真。5.3 FST和PBS的计算流程先准备好两个群体ID列表比如pop1.txt、pop2.txt每个一行一个样本名。# 成对FST窗口扫描 vcftools --vcf clean.vcf --weir-fst-pop pop1.txt --weir-fst-pop pop2.txt \ --fst-window-size 50000 --fst-window-step 10000 --out pairwise_fst得到窗口FST后在三群体情境下继续算PBS。这一步通常要写个小脚本。每个窗口的三群体两两d值先算出来再代入公式解出各分支PBS。我的建议是直接转成每个SNP的FST再在窗口内取平均这样跟其他方法的比较更方便。5.4 单倍型扫描的标准化操作selscan是计算iHS和XP-EHH的常用工具。假设你的VCF已经完成样本和位点过滤再准备图谱文件染色体、物理位置、遗传位置# 计算iHS selscan --ihh --vcf clean.vcf --map genetic.map --threads 8 --out output # 计算XP-EHH selscan --xpehh --vcf clean.vcf --map genetic.map --vcf-ref popA.txt --vcf-target popB.txt --threads 8 --out xpehh_out无论哪种方法跑完之后一定要归一化将结果按等位基因频率分成若干箱在每个频率箱内标准化否则低频位点的统计量方差极大阈值法会失效。归一化后通常取|Z-score|大于2或3的区域为显著信号。5.5 候选位点的注释与可视化拿到候选区间后用功能注释工具做注释把变异位置映射到基因和转录本上。这里不做广告只用你熟悉的开源注释工具就行。注释结果里我非常看重几个字段是否非同义突变、是否落在启动子区域、该基因在哪些组织中高表达如果有转录组数据。可视化方面强烈建议画三样东西全基因组FST/PBS的曼哈顿图、显著窗口的局部LD热图、候选位点在群体中的等位基因频率条形图。不要怕画图麻烦评审人和合作者最信任的不是你跑出了多炫的统计量而是你能把信号讲清楚。6. 常见问题与排查技巧实录6.1 为什么不同软件算出的FST差那么大这是几乎每个项目都会遇到的问题。原因有几种软件默认过滤标准不同比如是否排除单倍型缺失、是否过滤低质量位点。这会导致参与计算的位点集合完全不同。部分软件用单倍型频率估计FST部分用基因型频率二者假设不同。随机抽样导致的方差特别是样本量很小时单点FST方差极大。我的建议是写清楚“用什么软件、什么参数、基于多少SNP”并且在同一研究中固定同一版本。不要混用不同方法的输出做趋势比较。6.2 SFS出现“高频变异过多”是怎么回事高频变异很多的SFS可能意味着群体结构混入样本里有明显亚群但没有分开。参考基因组偏差参考等位基因由于组装错误出现的虚假高频等位基因。选择的真实信号例如平衡选择维持多个高频等位基因。数据过滤不当把多等位基因位点的次要等位基因合并方式不对导致某些等位基因频率虚高。排查顺序上先做样本聚类看是否存在结构再对比各染色体SFS看是否存在局部异常。如果所有染色体高频变异都偏多优先怀疑群体结构与数据问题只有某一小段区域高频变异特别多才值得考虑平衡选择。6.3 选择信号扫描得到太多候选位点怎么办全基因组扫描通常会得到几百甚至上千个候选区域这再正常不过了。问题是很多假阳性来自重复区域、拷贝数变异或参考序列组装错误。我的处理步骤是用重复序列注释文件过滤掉落在重复区域的候选位点。用序列比对检查候选区是否存在异常高覆盖排除拷贝数变异干扰。如果数据中包含多个亲缘关系较近的个体考虑先做亲缘关系过滤否则群体结构会制造大量假信号。最后结合功能注释和文献信息做优先级排序。优先关注同义与非同义突变比值显著、且落在已知生态适应相关基因上下游的位点。6.4 运行速度和资源优化建议全基因组重测序数据算XP-EHH或全基因组SFS时计算量确实不小。我建议按染色体并行每条染色体单独跑扫描最后合并结果。这不仅能利用多核资源还能防止一条染色体上的故障拖垮全局。VCF文件先做双等位位点过滤、剔除高缺失位点数据量至少能降一半。单倍型扫描工具对内存占用不小采用分块运行或者每个染色体逐一扫描可以明显减少内存峰值。优先用BGZF压缩的VCF格式而不是纯文本VCF文件体积和读取速度差异非常大。7. 几点实操后的真心话做了这么多年群体基因组分析我最大的感触是方法永远在更新但分析的逻辑骨架没变——先理解数据怎么来的再理解统计量在测什么最后再谈生物学解释。如果你刚开始接触这套流程不用急着把每一种新方法都跑一遍。老老实实把SFS做对、把FST的窗口参数解释清楚、把XP-EHH和FST的联合使用练熟就能解决大多数实际项目里七八成的问题。我见过太多人栽在第一步数据质控没做好SFS极性搞反FST用了两种不同定义却在同一张图里混着画最后得出一堆互相矛盾的“候选基因”。另外凡是涉及全基因组扫描的分析一定要留好版本记录和分析日志。统计结果的微小差异经常来源于过滤样本数或位点数不一致没有日志几乎不可能排查。最后再安利一个小技巧拿到一套候选位点之后先别急着做功能富集先去看这些候选位点在染色体上的分布是否均匀。如果集中挤在几条染色体或者某个着丝粒附近大概率是结构变异或重组抑制区域的假信号。把这一步排除干净再进入功能注释你会省掉非常多的返工时间。 
长河君 · 2026/10/13 5:14:09 · 阅读 1.2k 等位基因频谱、群体分化与选择信号检测：群体基因组分析全攻略
群体基因组学这几年是真的火但不少刚接触的朋友跟我聊的时候总感觉概念都眼熟真拿到自己的数据又不知道从哪里下手。等位基因频谱、群体分化、选择信号检测这三块是群体基因组分析的三大支柱也是从“拿到VCF文件”到“讲出生物学故事”之间最关键的几步。这篇东西我打算按自己做项目时的思路来写把每个概念的来龙去脉、背后的统计逻辑、跑的每一步命令还有踩过的坑都摊开说清楚希望能帮你把这条分析链路完整打通。我先把结论放在这儿这三块内容有很强的递进关系——等位基因频谱是理解数据的基本盘群体分化是描述群体关系的度量衡选择信号检测则是真正回答“哪个基因在进化中很重要”的核心武器。想做好分析缺了任何一环都会让结论站不住脚。接下来我从原理、实操到选型一步步拆开讲。1. 核心思路为什么把SFS、FST和选择检测放在同一条流水线上很多初学者会把这几个概念当成独立的模块分别去学这个方向就错了。真实项目里它们是环环相扣的先从测序数据里拿到最基础的多态位点信息然后通过等位基因频谱观察整体变异规律再用群体分化指标精细刻画亚群之间的关系最后基于这些背景信息去寻找偏离中性的位点。我习惯把这套流程称为“描述—比较—检测”三段式。1.1 一条完整的分析链路是怎样串起来的拿到重测序数据之后第一步是比对、变异检测、过滤最终得到全基因组的VCF文件。这个阶段不考虑生物学阐释只保证数据质量。接下来就要算等位基因频谱也就是把所有多态位点按照“次等位基因出现的频次”分类统计。这样做的目的是理解数据整体长什么样有什么样的测序深度、参考基因组方向是不是正确、有没有明显的缺失模式。紧接着用FST固定指数这一类指标去群体间做比较。为什么要这一步因为任何“选择信号”都是相对于背景来定义的。如果两个群体本身分化就非常剧烈某一个位点显示出的高分化很可能只是群体历史结构造成的假象而非自然选择的作用。最后才做选择信号检测。检测手段本质上是去寻找偏离中性模型预期的异常值。在已知群体分化背景的前提下你可以判断这些异常到底是真的受选择还是群体结构、基因流的自然结果。1.2 频繁踩坑的“视角混淆”问题我发现很多人在做选择信号检测的时候会犯一个经典错误直接在未做群体分层的样本里跑全基因组扫描。假如你有两个生态型存在明显分化却把它们混在一起当一个大群体去计算中性检验那么得到的“信号”往往是群体分层的噪音根本不是选择。正确的做法是“先描述、再比较、后检测”每一步的结果都在为下一步筛选提供背景信息。所以这篇文章的编排逻辑就是按这个流程来的先讲SFS怎么正确生成和解读再讲FST和PBS这类指标的生物学含义与参数选择最后讲选择信号检测中不同方法各自的适用场景以及怎么把多种证据串成一个让人信服的故事。2. 等位基因频谱SFS生成逻辑、中性预期与关键误区等位基因频谱Site Frequency SpectrumSFS是所有群体遗传分析的基石。简单来说它统计的是在一个群体中不同频次等位位点的数量分布。比如有10条染色体出现1次的突变有多少个、出现2次的有多少个、直到出现9次的有多少个。2.1 从VCF到SFS的三种姿势实际项目中生成SFS有三种常见方案按需求选择就好直接基于VCF文件用vcftools统计这种方法适合长期过滤后质量较高的基因型数据速度很快但也容易丢失低频变异的检出能力。基于BAM文件的ANGSDSAMtools流程直接用测序比对结果推算等位基因频谱能充分利用低频信息还能引入测序深度和质量过滤参数。这种方案对测序深度不高的数据特别友好重测序项目我基本都推荐。基于变异注释的衍生SFS如果你想区分同义与非同义位点就需要结合功能注释信息在得到基础SFS后分类汇总。我个人的建议是除非你的数据覆盖度和基因型检出率都很理想否则尽量直接用BAM做。原因很现实VCF里低频变异的基因型叫出率通常不高一个位点如果只有一半个体被成功分型那你的SFS从来源上就已经失真了。2.2 中性模型预期为什么SFS是“倒J型”在标准中性模型下一个样本量为n的二倍体群体中等位基因频谱的预期是每个频率类别的位点数量正比于1/f。这意味着低频变异应该是最多的。如果你的SFS显示出大量高频变异那么就需要警惕是否受到了选择、群体扩张或数据偏误的影响。这也是理解很多检验统计量的基础。比如中性偏离检验D值常称Tajimas D这个经典指标它的逻辑就是比较两种对遗传变异的估计量——基于平均成对差异的π值和基于多态位点数目的Watterson θ值。如果群体正在经历收缩或者说有正选择造成“搭车效应”那么高频变异会比中性预期偏多D值会显著偏离零。我在实际项目中见过很多用户对D值正负号的理解有问题。简单记法D为负说明低频变异过多常见于群体扩张、正选择或数据的低深度个体混入D为正说明高频变异过多常见于群体收缩、平衡选择或混合样本中存在群体结构。具体问题要结合SFS的形状一起判断别单独看一个值。2.3 计算SFS时的“致命细节”生成SFS看似只是跑一条命令实则细节决定成败尤其要注意以下几项极性SFS分“折叠”和“非折叠”两种。非折叠SFS需要指定一个外群或祖先状态来确定哪个是衍生等位基因如果参考基因组跟你的物种演化距离适中可以直接把参考基因组等位基因作为祖先状态。没有外群时只能做折叠SFS相当于把频率w和1-w合并灵敏度会明显下降。测序深度深度过低会导致低频等位基因被当成测序错误过滤掉深度过高又容易引入参考偏差。我一般按平均深度30X左右的数据将最低基因型质量设为20最低位点深度设为8到10。位点过滤顺序做SFS之前过滤做多了会把真实低频变异删掉做少了又会混入大量测序噪声。我摸索出来的经验是先过滤样本与位点质量再计算SFS最后在确认整体形态没问题之后再对SFS做进一步校正。注意SFS是所有后续分析的基础。如果你发现SFS异常前面做得再漂亮的FST和选择扫描结果都不可信因为基础分布已经不对了。3. 群体分化从单点FST到窗口PBS的完整推导群体分化通常是我们要回答的第一个生物学问题不同生态型或地理群体之间是否存在显著的遗传差异哪些区域差异最大。这里最常用的就是FST但市面上的FST版本五花八门算出来差别很大新手经常被搞晕。3.1 FST的本质不复杂但定义要选对FST这个概念可以理解为群体间遗传差异占总体遗传差异的比例。如果所有群体等位基因频率完全相同FST就是0如果各群体固定了完全不同的等位基因FST就接近1。但实际计算时有多种估计方式常见的有以下几种Wright的原始定义需要知道群体真实等位基因频率很少直接用于真实数据。Weir和Cockerham基于方差分析开发的估计量它把群体间和群体内的方差分量做矩估计是目前软件中最常用的算法。Hudson基于等位基因频率简单估计的FST计算速度快常用于滑动窗口扫描。基于成对序列差异的FST估计专门用于测序深度不齐或不依赖基因型时的估算。我自己的经验是做窗口扫描用Hudson的估计方式比较稳定做单点分析、候选位点筛查时可以结合Weir和Cockerham的估计值实际报告中一定要写清楚你用的是什么方法和参数。3.2 窗口化扫描的参数如何决定FST滑动窗口扫描中窗口大小和步长的选择会直接改变信号分辨率。窗口太小单点噪声大窗口太大真实信号被平均稀释。这个参数没有绝对标准但可以按以下逻辑去确定如果重测序序列全长约500Mb左右我通常起步用50kb窗口、10kb步长。如果是芯片数据或低密度标记直接用单点FST并按物理距离绘制可以绕开窗口平均的问题。如果已经有候选区间可以直接在候选区段内用5kb窗口精确定位。关键技巧是“先粗扫、再精扫”。先用大窗口找出信号区间再缩小窗口在区间内确定候选基因范围。如果一步就直接用小窗口全基因组范围的假阳性会多到你怀疑人生。3.3 PBS把三群体故事讲清楚的分支统计量FST可以告诉你哪些区域在群体间分化明显但它说不清楚分化具体发生在哪条演化分支上。这时候就要用分支统计量PBS群体分支统计量。它的逻辑是把三个群体的FST值做-log转换然后通过代数推算将每个群体分支上的独特分化量分离出来。计算方式上三群体的FST两两转换后可以得到一个值衡量目标群体自与其他两个群体分开之后积累的分化量。实际操作时注意三个群体的选择非常关键要按演化关系排好并且尽量选关系清楚的人群组合这样结果的解释才直接。PBS只能做三群体组合如果有四个以上群体可以组合多种三群体方案并取交集。PBS受群体历史影响很大种群瓶颈会让全基因组背景PBS升高所以不能只看绝对PBS值要看相对全基因组分布的百分位。实操时我用R脚本或Python并行计算PBS时有个习惯先生成所有位点的三群体FST再转成两两d值矩阵最后解出各分支的PBS。如果你用vcftools的滑动窗口FST输出记得保留位点位置因为PBS计算需要按窗口对应的区间合并。3.4 当群体历史比较复杂时怎么办现实中很多数据并不满足“三群体独立分化”的理想假设。比如存在广泛基因流、近期群体混合、种群瓶颈等。此时PBS的绝对值会受到严重干扰。我是这么处理的先用主成分分析、结构分析这些手段确认群体遗传结构与演化树基本一致。如果结构分析和PBS设定冲突最好放弃PBS改用基于扩展单倍型纯合度的方法或者用更稳健的群体分化统计量做两两比较。记住任何统计指标的生物学解读都依赖于群体历史模型模型不对指数再漂亮也是自欺欺人。4. 选择信号检测从频谱偏离到单倍型痕迹的综合研判选择信号检测是群体基因组学里最激动人心、也最容易翻车的一步。目前成熟的方法大致分三类基于SFS的方法、基于群体分化的方法、基于单倍型结构的方法。三者检测的目标不同单独用哪一种都有明显局限靠谱的研究基本都靠“多方法取交集”。4.1 基于SFS的检验适合什么场景这类方法包括中性偏离检验D值常称Tajimas D、D常称Fay和Wu的H等。它们通过比较真实SFS与中性预期SFS的偏离来判断是否存在选择。中性偏离检验D值对群体历史很敏感瓶颈和扩张都会造成假阳性所以通常不作为单一证据。对于新近完全扫荡形成的选择信号D值会显著为负但它的分辨率很低很难定位到具体位点。实际操作中我倾向于把这类统计量作为“数据质量与历史事件”的初步探针而不是直接用它出候选基因。如果整个群体的D值全基因组普遍异常那大概率是该群体经历过扩张或收缩应该先去联合历史动态模拟而不是急着找受选择基因。4.2 基于分化的方法要会挑“极端尾部”FST和PBS这类方法在当前数据里用得最广因为它们能直接给出受选择区域的候选位点。这类方法的假设是局部受选择会让两个亚群在该区域的等位基因频率发生显著分化远大于基因组背景的中性漂移。实际操作上我通常会先全基因组扫描获得所有窗口的FST或PBS值再取前1%或5%作为top候选。为什么不能直接设一个绝对值阈值因为不同群体的遗传背景、样本量、标记密度差异很大绝对阈值没有普适性。做选择信号检测时许多流程还会用到交叉群体扩展单倍型纯合度方法常称XP-EHH或群体分化方法常称XP-CLR。这两种方法的原理分别基于单倍型长度和群体分化模式。XP-EHH对仍在传播中的不完全扫荡敏感XP-CLR对已固定的扫荡和长期选择更有效。它们虽然原理不同但共享一个思想相对于背景中性模式受选择区域的等位基因频率分化程度或者单倍型纯合度会有突出的信号。4.3 基于单倍型的检验iHS与XP-EHH怎么选基于单倍型的方法利用的是“受选择位点附近的单倍型往往很长且纯合度高”这一特性。核心指标包括扩展单倍型纯合度EHH、内部单倍型纯合度得分iHS和跨群体扩展单倍型纯合度XP-EHH。iHS在单个群体内检测正在发生或者近期发生的正向选择比较同一核心位点上衍生与祖先等位基因携带的单倍型长度。iHS的绝对值大说明核心位点附近的单倍型异常长。XP-EHH在两个群体之间比较当目标群体的受选择单倍型长度和频率都显著高于另一个群体时信号会很强。它是当前应用最广的跨群体方法之一。以我的项目经验如果要找不同生态型对特定环境适应产生的信号XP-EHH往往比iHS效果好很多因为这里的对比本身就是生态型之间。而同一种群内部的驯化选择信号iHS更有优势。另外XP-EHH对重组率变化的鲁棒性比iHS好但需要两个群体质量都比较高。selscan这类工具可以同时计算和归一化这些指标跑的时候记得要提供遗传图谱或重组率参数别直接拿物理距离硬算。4.4 如何组合证据避免“单指标自嗨”单看任何一项指标都可能被假阳性欺骗多方法交集是最稳妥的做法。我常用的组合策略第一层全基因组XP-CLR或PBS扫描锁定前1%的极端窗口。第二层在候选窗口内计算FST、XP-EHH要求至少两项指标显著。第三层在最终候选区域检查SFS形态和单倍型网络确认不是组装错误或极端缺失导致的假信号。第四层结合功能注释筛选出非同义突变、启动子变异等有生物学意义的候选位点。5. 实战工作流从VCF到候选位点注释的一次性跑通理论部分讲得再多不如给你一套可以直接上手的实操流程。我按自己做真实项目的习惯整理了一份从VCF到候选位点注释的常用命令和参数希望能成为你的起步模板。5.1 数据准备与质控过滤第一步是对VCF做质控。过滤条件基于你的数据质量设定我常用以下参数# 保留双等位位点、过滤缺失过多、低QUAL位点 vcftools --vcf raw.vcf --remove-indels --min-alleles 2 --max-alleles 2 \ --max-missing 0.7 --minQ 30 --recode --out clean这个环节注意两点。一是别一上来就按极端的次要等位基因频率过滤比如只保留5%以上的位点那后面低频变异相关的分析就没法做了。我会在生成SFS之前保留所有多态位点只在做单倍型分析时考虑低频过滤。二是染色体名称和个体ID最好在项目开始时就统一好不然后面每跑一步都在跟命名错误斗智斗勇。5.2 SFS和中性检验的计算如果你有BAM文件且希望更灵敏地捕捉低频变异建议用以下思路生成SFS并计算θ和D值# 使用ANGSD生成等位基因频谱 angsd -bam bamlist.txt -ref genome.fa -out sample \ -doSaf 1 -doThetas 1 -GL 1 -minMapQ 30 -minQ 20 -anc genome.fa # 依据SFS计算核苷酸多样性与D值等指标 realSFS sample.saf.idx sample.sfs如果只有VCF且质量较好直接走vcftools的窗口D值更快vcftools --vcf clean.vcf --TajimaD 50000 --out window_tD注意ANGSD的祖先等位基因设置直接决定你的SFS是否是折叠的。若没有可靠外群建议用折叠SFS以免把参考基因组等位基因误当祖先状态导致衍生等位基因频谱失真。5.3 FST和PBS的计算流程先准备好两个群体ID列表比如pop1.txt、pop2.txt每个一行一个样本名。# 成对FST窗口扫描 vcftools --vcf clean.vcf --weir-fst-pop pop1.txt --weir-fst-pop pop2.txt \ --fst-window-size 50000 --fst-window-step 10000 --out pairwise_fst得到窗口FST后在三群体情境下继续算PBS。这一步通常要写个小脚本。每个窗口的三群体两两d值先算出来再代入公式解出各分支PBS。我的建议是直接转成每个SNP的FST再在窗口内取平均这样跟其他方法的比较更方便。5.4 单倍型扫描的标准化操作selscan是计算iHS和XP-EHH的常用工具。假设你的VCF已经完成样本和位点过滤再准备图谱文件染色体、物理位置、遗传位置# 计算iHS selscan --ihh --vcf clean.vcf --map genetic.map --threads 8 --out output # 计算XP-EHH selscan --xpehh --vcf clean.vcf --map genetic.map --vcf-ref popA.txt --vcf-target popB.txt --threads 8 --out xpehh_out无论哪种方法跑完之后一定要归一化将结果按等位基因频率分成若干箱在每个频率箱内标准化否则低频位点的统计量方差极大阈值法会失效。归一化后通常取|Z-score|大于2或3的区域为显著信号。5.5 候选位点的注释与可视化拿到候选区间后用功能注释工具做注释把变异位置映射到基因和转录本上。这里不做广告只用你熟悉的开源注释工具就行。注释结果里我非常看重几个字段是否非同义突变、是否落在启动子区域、该基因在哪些组织中高表达如果有转录组数据。可视化方面强烈建议画三样东西全基因组FST/PBS的曼哈顿图、显著窗口的局部LD热图、候选位点在群体中的等位基因频率条形图。不要怕画图麻烦评审人和合作者最信任的不是你跑出了多炫的统计量而是你能把信号讲清楚。6. 常见问题与排查技巧实录6.1 为什么不同软件算出的FST差那么大这是几乎每个项目都会遇到的问题。原因有几种软件默认过滤标准不同比如是否排除单倍型缺失、是否过滤低质量位点。这会导致参与计算的位点集合完全不同。部分软件用单倍型频率估计FST部分用基因型频率二者假设不同。随机抽样导致的方差特别是样本量很小时单点FST方差极大。我的建议是写清楚“用什么软件、什么参数、基于多少SNP”并且在同一研究中固定同一版本。不要混用不同方法的输出做趋势比较。6.2 SFS出现“高频变异过多”是怎么回事高频变异很多的SFS可能意味着群体结构混入样本里有明显亚群但没有分开。参考基因组偏差参考等位基因由于组装错误出现的虚假高频等位基因。选择的真实信号例如平衡选择维持多个高频等位基因。数据过滤不当把多等位基因位点的次要等位基因合并方式不对导致某些等位基因频率虚高。排查顺序上先做样本聚类看是否存在结构再对比各染色体SFS看是否存在局部异常。如果所有染色体高频变异都偏多优先怀疑群体结构与数据问题只有某一小段区域高频变异特别多才值得考虑平衡选择。6.3 选择信号扫描得到太多候选位点怎么办全基因组扫描通常会得到几百甚至上千个候选区域这再正常不过了。问题是很多假阳性来自重复区域、拷贝数变异或参考序列组装错误。我的处理步骤是用重复序列注释文件过滤掉落在重复区域的候选位点。用序列比对检查候选区是否存在异常高覆盖排除拷贝数变异干扰。如果数据中包含多个亲缘关系较近的个体考虑先做亲缘关系过滤否则群体结构会制造大量假信号。最后结合功能注释和文献信息做优先级排序。优先关注同义与非同义突变比值显著、且落在已知生态适应相关基因上下游的位点。6.4 运行速度和资源优化建议全基因组重测序数据算XP-EHH或全基因组SFS时计算量确实不小。我建议按染色体并行每条染色体单独跑扫描最后合并结果。这不仅能利用多核资源还能防止一条染色体上的故障拖垮全局。VCF文件先做双等位位点过滤、剔除高缺失位点数据量至少能降一半。单倍型扫描工具对内存占用不小采用分块运行或者每个染色体逐一扫描可以明显减少内存峰值。优先用BGZF压缩的VCF格式而不是纯文本VCF文件体积和读取速度差异非常大。7. 几点实操后的真心话做了这么多年群体基因组分析我最大的感触是方法永远在更新但分析的逻辑骨架没变——先理解数据怎么来的再理解统计量在测什么最后再谈生物学解释。如果你刚开始接触这套流程不用急着把每一种新方法都跑一遍。老老实实把SFS做对、把FST的窗口参数解释清楚、把XP-EHH和FST的联合使用练熟就能解决大多数实际项目里七八成的问题。我见过太多人栽在第一步数据质控没做好SFS极性搞反FST用了两种不同定义却在同一张图里混着画最后得出一堆互相矛盾的“候选基因”。另外凡是涉及全基因组扫描的分析一定要留好版本记录和分析日志。统计结果的微小差异经常来源于过滤样本数或位点数不一致没有日志几乎不可能排查。最后再安利一个小技巧拿到一套候选位点之后先别急着做功能富集先去看这些候选位点在染色体上的分布是否均匀。如果集中挤在几条染色体或者某个着丝粒附近大概率是结构变异或重组抑制区域的假信号。把这一步排除干净再进入功能注释你会省掉非常多的返工时间。