B05
经典正选择信号在 1000 Genomes 30× 高覆盖重测序数据上的独立复现:以数据版本更替为自变量,量化低覆盖填补对 iHS/XP-EHH 扫描的影响
1 · 研究问题
基于 1000 Genomes Project phase 3 低覆盖(平均 7.4×)数据发表的经典正选择(positive selection)信号,在同一批个体的 30× 高覆盖重测序(NYGC,3,202 例)call set 上重跑相同的单体型扫描统计量后,能否重现?重现率与信号强度、等位基因频率谱、填补(imputation)依赖程度之间是什么关系——即"低覆盖 + 统计填补"这一数据生成方式给单体型类选择扫描带来了多大的系统偏差?
2 · 研究背景与空白
背景。 检测人类基因组上的近期正选择,主流方法依赖单体型长度:iHS(integrated haplotype score)、nSL、XP-EHH(cross-population extended haplotype homozygosity)都通过比较携带不同等位基因的单体型衰减速度来识别"频率升高但单体型尚未被重组打断"的位点。这类统计量对单体型定相(phasing)质量极其敏感——而定相质量又直接取决于测序覆盖度与填补算法。
已有工作到哪一步。 1000 Genomes phase 3(2015)的 2,504 例基于低覆盖 WGS(平均 7.4×)+ 高覆盖外显子组 + 芯片基因型的组合调用,是过去十年绝大多数人类选择扫描的底层数据;LCT/MCM6(欧洲乳糖耐受)、SLC24A5(欧洲肤色)、EDAR(东亚毛发与汗腺)等经典位点的信号强度都是在这套数据上标定的。2022 年 Byrska-Bishop 等(Cell)发布了同一批个体的 30× 高覆盖重测序(3,202 例,含 602 个三口之家),明确报告相对 phase 3 在罕见 SNV、INDEL 与结构变异上灵敏度与精确度显著提升。此外 2019 年已有工作把 phase 3 重新比对/调用到 GRCh38 上,使两版可在同一坐标系比较。
空白在于:数据版本更新了,但没有人系统检验过基于旧版数据的选择扫描结论在新版上是否成立。 Byrska-Bishop 等评估的是变异调用层面的灵敏度/精确度,不是下游群体遗传推断的稳定性;而选择扫描依赖的是单体型结构,恰恰是低覆盖填补最可能引入系统性伪迹的地方。这是研究手册中最稳的差异化轴——换样本(此处是换同一样本的新数据版本)做独立检验:方法完全公开可复现、算力门槛可控、结论正负都成立(重现则为经典结论提供一次独立支持;不重现则指出一整类文献的数据依赖性)。
3 · 可检验假设
- H1:在 30× 数据上重跑 iHS 后,phase 3 中排名前 100 的候选窗口(每群体)有 ≥ 80% 仍落在 30× 结果的前 1% 窗口内;但重现率随信号强度单调下降——phase 3 排名 500–1000 的窗口重现率 ≤ 50%,两档差值的 95% 自助法置信区间不跨 0。
- H2(机制假设):不重现的窗口在两版之间的定相开关错误率(switch error)代理指标(用 602 个三口之家的孟德尔一致性/单体型长度分布推得)显著高于重现窗口;且不重现窗口富集于低重组率区域与低等位基因频率区间(0.05 < MAF < 0.15)。若 H2 被否证(不重现窗口无系统特征),说明差异是随机噪声,经典结论的稳健性反而更强,这同样是有效结论。
4 · 量化验收标准
- 方法学校验(硬门槛):先在 phase 3 数据本身上重跑扫描,要求复现已发表的经典结果:(a) CEU(欧洲裔)中 LCT/MCM6 区域与 SLC24A5 区域落入 iHS |score| 的全基因组前 0.1% 窗口;(b) CHB/CHS(东亚)中 EDAR 区域落入前 0.1%;(c) 与已发表的 1000 Genomes Selection Browser / Voight 等公布的 iHS 值在共有 SNP 上的 Spearman 相关 ≥ 0.85。三项任一不达标,管线不可信,后续版本对比全部无效。
- 统计口径预先写死:(a) 扫描以 100 kb 非重叠窗口汇总(同时报 50 kb 与 200 kb 作敏感性分析);(b) 全基因组多重检验一律报经验分位数与 BH-FDR,不报裸 p 值;(c) "重现率"这类不平衡二分类判断(真正的选择窗口占全基因组比例 < 1%)报 PR-AUC 与固定召回处的精确率,不报 accuracy,原因是负类占 99% 以上,全判"非候选"即得 99% accuracy 而毫无信息;(d) 所有重现率给 1,000 次自助法(按染色体块重抽样,保留连锁不平衡结构)95% 置信区间;(e) 两版比较必须限制在两版共有的双等位 SNP 且 MAF ≥ 0.05 的子集上,共有位点数与各版特有位点数逐染色体报出。
- 混杂因素必须显式处理并写明:(a) 坐标系统——两版必须在 GRCh38 上比较(使用 1000G phase 3 的 GRCh38 重调用版本,而非自行 liftOver;若必须 liftOver,须报告映射失败率);(b) 遗传图谱——两版使用同一份 GRCh38 重组图谱,图谱来源与版本写入方法;(c) 样本集——30× 版含 698 个亲属样本,主分析必须限制到与 phase 3 完全相同的 2,504 例无关个体,亲属样本只用于 H2 的定相质量评估。
- 规模下限:≥ 5 条染色体(须含 chr2、chr4、chr11、chr15 这些携带经典靶点的染色体,外加 chr22 或 chr20 作为"无已知强信号"的阴性对照染色体);≥ 4 个群体(至少覆盖 EUR、EAS、AFR、SAS 各一);每染色体每群体的分析 SNP 数 ≥ 200,000。
- 可复现性:bcftools 过滤命令、selscan 参数、图谱文件校验和、全部脚本与结果表开源;提供 chr22 单染色体的降规模复现包(≤ 2 小时可跑完)。
5 · 数据与工具
| 用途 | 来源 / 工具 |
|---|---|
| 30× 高覆盖定相 call set | 1000 Genomes FTP:http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/(2,504 无关个体的定相 VCF,GRCh38)与 .../working/20220422_3202_phased_SNV_INDEL_SV/(3,202 例含亲属)。完全开放访问,无需任何 DAC 审批——1000 Genomes 的知情同意即为全球公开发布。单染色体压缩 VCF 约 1–6 GB,5 条染色体约 15–35 GB 磁盘 |
| phase 3 低覆盖定相 call set(GRCh38 重调用版) | 1000 Genomes FTP 的 GRCh38 phase 3 目录(2019 年 remapping 发布)。确切路径与文件命名需核实;若 GRCh38 版不可用,退回 GRCh37 版 + 官方 chain 文件 liftOver,并报告映射失败率 |
| 样本–群体对应表 | 1000 Genomes integrated_call_samples_v3.20130502.ALL.panel 及 30× 版对应的 pedigree 文件 |
| GRCh38 重组图谱 | Beagle 发布的 GRCh38 genetic maps(公开下载)或 HapMap GRCh37 图谱经官方 liftOver。两版分析必须用同一份图谱,来源写入方法 |
| VCF 处理 | bcftools(view/filter/annotate/isec),流式处理,内存占用与文件大小无关,MB 级 |
| 选择扫描 | selscan v2(C++,多线程;iHS/nSL/XP-EHH 均支持);交叉校验用 Python scikit-allel(allel.ihs)。单群体单体型矩阵:CEU n=99 → 198 单体型 × 1.5M SNP 的 int8 数组约 0.3 GB,16 GB 内存充裕;若一次性载入全部 2,504 例(5,008 单体型 × 1.5M)则约 7.5 GB,接近上限,因此管线一律按群体切片 |
| 算力口径 | 纯 CPU。selscan iHS:200 单体型 × 1M SNP、4 线程约 15–60 分钟/染色体/群体。5 染色体 × 4 群体 × 2 版本 = 40 次运行 ≈ 15–40 小时,可后台分夜跑。全基因组 26 群体两版全跑属超出范围(约 25 倍工作量 + 200 GB 以上磁盘),方案中不设此路径 |
| 已发表选择信号对照 | Voight 等 2006 与 1000 Genomes Selection Browser 公布的 iHS 结果表;Byrska-Bishop 等 2022 的变异调用评估表。仅用于校验与对比,不计入本项目的数据贡献 |
| 伦理 | 1000 Genomes 为完全去标识化、以群体为单位公开发布的开放数据,不涉及可识别个人、不涉及人体实验,无需伦理审批。方案中不得做任何个体层面的表型推断 |
6 · 方法路径
- 装环境(bcftools + selscan + scikit-allel),用 chr22 小片段跑通两条管线,先完成验收标准第 1 条的三项经典信号复现,把复现结果与文献值并列写入
validation/。 - 核实数据版本与坐标:确认 phase 3 GRCh38 重调用版的可用性与路径;用
bcftools isec求两版共有 SNP 集合,逐染色体报共有/特有位点数与 MAF 分布。此步是全课题最关键的前置核实,不得边做边发现。 - 按预先写死的规则做过滤(双等位 SNP、MAF ≥ 0.05、非缺失、共有位点集),按群体切片生成单体型输入;过滤前后位点数逐步记录。
- 在两版数据上分别跑 iHS 与 nSL(4 群体 × 5 染色体),汇总到 100 kb 窗口,生成排名表。
- 计算重现率:按 phase 3 排名分档(1–100、101–500、501–1000)统计落入 30× 前 1% 的比例,自助法给区间;同时报 PR-AUC 与固定召回处精确率。检验 H1。
- 检验 H2:用 602 个三口之家计算区域级定相一致性代理指标,把不重现窗口与重现窗口在重组率、MAF 谱、定相一致性上做对比(BH 校正)。
- 独立交叉校验:换一类不依赖单体型定相的统计量(群体分化 F_ST 与 PBS)在同两版数据上重算同一批窗口。若 F_ST/PBS 的两版一致性远高于 iHS,则差异确由定相/填补引起而非位点集差异——这是本课题最关键的一次对照。
7 · 新颖性边界
本课题不提出新的选择检测统计量,不声称发现任何新的选择位点,不声称 LCT/SLC24A5/EDAR 处于正选择为本项目发现(这些是二十年来的经典结论),不声称 30× 数据优于 phase 3(Byrska-Bishop 等已在变异调用层面报告)。丘奖 2025 年铜奖论文《Association Study of SNPs in X-Chromosome Inactivation Escape Regions…》做的是关联分析,与本课题的选择扫描/数据版本对比路径不同。
已有工作具体完成了什么:Voight 等与后续大量论文在 phase 3(及更早的 HapMap)上标定了 iHS/XP-EHH 的候选位点清单;Byrska-Bishop 等(Cell, 2022)报告了 30× 相对 phase 3 在 SNV/INDEL/SV 调用灵敏度与精确度上的提升,并发布了定相 panel。没有人把这两件事接起来,检验下游选择推断是否随数据版本改变。
本项目的贡献(且是主结论):把评价维度从"变异调用质量"换成"下游群体遗传推断的版本稳健性",给出分信号强度档的重现率曲线,以及"哪一类窗口最不稳健"的可操作特征刻画(H2)。这是研究手册所列最稳的差异化轴的一个变体:换的不是新群体,而是同一群体的新数据生成方式。
为什么有价值:过去十年数以百计的人类适应性进化结论建立在 phase 3 之上。若高排名信号重现率高,这些结论获得一次真正独立的技术支持;若中等排名信号大面积不重现,则文献中"新发现的选择位点"这一类主张需要重新评估。两种结果都构成有效结论。
风险声明:若重现率接近 100%(差异不显著),必须给出自助法误差棒证明本设计有能力分辨假设中的分档差异,否则"高度稳健"的结论无效。
8 · 决策门槛(go / no-go)
- 第 4 周末:确认 phase 3 GRCh38 重调用版可下载且与 30× 版坐标一致。若不可得,立即降级:其一,改用官方 chain 文件把 30× 版 liftOver 到 GRCh37,与原始 phase 3 比较,并把 liftOver 失败率作为一条限制条款报出;其二,把比较对象从"两个数据版本"改为"同一 30× 数据的两种定相方式"(如 SHAPEIT/Beagle 统计定相 vs 三口之家 pedigree 定相),此时研究问题变为"定相方法对选择扫描的影响",主结论框架(分档重现率 + 不重现窗口特征刻画 + F_ST 对照)完全保留。
- 第 8 周末(硬门槛):验收标准第 1 条必须通过。若与已发表 iHS 的相关 < 0.85,先排查图谱版本与祖先等位(ancestral allele)标注(这两项是 iHS 最常见的错误源);四周内仍不通过则降级:把校验基准从"数值相关"改为"经典位点的排名分位数复现"(更稳健,阈值定为三个经典位点均入前 1%),并在论文中说明数值不可比的原因。
- 第 14 周末:确认磁盘与运行时间。若 5 条染色体的下载+运行超出可承受范围,立即降级:染色体数降到 3(chr2 + chr15 + chr22 对照),群体数降到 2(CEU + CHB),并在结论中相应收窄推广范围。主结论框架不变。
- 第 28 周末:若 H1 与 H2 均未达阈值,按风险声明写"版本稳健"结论,并把 F_ST/PBS 对照与误差棒作为核心证据。
- 第 36 周结果冻结,第 44 周英文 PPT 初稿。
- 需提前核实而非边做边发现:(a) 30× 定相 VCF 的确切 FTP 路径与是否已按无关个体子集发布;(b) GRCh38 重组图谱的来源与是否覆盖全部目标染色体;(c) 祖先等位标注文件(Ensembl EPO 比对导出)在 GRCh38 上的可得性——iHS 的极性依赖此文件,缺失会使结果无意义。三项在第 6 周前查清。
- 选择前提:适合已熟悉命令行、有耐心处理大文件与坐标系统的学生。本路线的最大困难是"两版数据的位点集不同"这一根本性不可比问题——这个困难本身就是研究内容,处理它(共有位点子集 + F_ST 对照)就是本课题方法学的核心。
- 预算裁剪顺序:先砍 nSL(只留 iHS),再砍 H2 的定相一致性分析(保留重组率与 MAF 谱特征刻画),再砍群体数,最后砍染色体数。砍到只剩"chr2 + chr15、CEU + CHB、iHS 两版分档重现率 + F_ST 对照"时,主结论仍成立。