HRC-1000G-check-bim.pl William Rayner (牛津大学) 写的 GWAS 基因型数据质控/预处理工具

• https://genepi.github.io/michigan-imputationserver/prepare-your-data/#download-tool-and-sites
• https://github.com/isglobal-brge/imputeInversion/blob/master/HRC-1000G-check-bim.pl

Read , lines 1 to 400

Read , lines 400 to 800

Read , lines 800 to 1093

这个脚本是 William Rayner (牛津大学) 写的 GWAS 基因型数据质控/预处理工具,用于在把 PLINK 格式的数据拿去做 HRC (Haplotype Reference Consortium) 或 1000 Genomes 参考面板的填补 (imputation) 之前,对 .bim 文件逐位点与参考面板进行比对和修正。

核心目的

把科研数据集的坐标、等位基因、链方向、ref/alt 都统一成参考面板的格式,使后续 imputation 服务器(如 Michigan Imputation Server / TOPMed)能够正确接受数据。

具体做了什么(逐条)

1. 读取参考面板

   • -h 用 HRC 参考文件(默认 HRC.r1.GRCh37.autosomes.mac5.sites.tab)
   • -g 用 1000G 参考文件(默认 1000GP_Phase3_combined.legend)
   • 把每个位点的「染色体-位置」映射到 id、ref/alt 等位基因、alt 等位基因频率

2. 逐行比对 .bim 文件里的每个 SNP,输出若干修正清单文件:

   • Exclude-*.txt —— 需要删除的 SNP(不在参考面板中的、重复的、Indel、无法匹配等位基因的、XY/Y/MT 染色体的)
   • Chromosome-*.txt —— 需要改染色体的 SNP
   • Position-*.txt —— 需要改位置的 SNP(错配的 SNP)
   • ID-*.txt —— 需要改 rs ID 的 SNP
   • Strand-Flip-*.txt —— 需要做链翻转 (flip) 的 SNP
   • Force-Allele1-*.txt —— 需要强制指定 ref 等位基因的 SNP
   • LOG-*.txt —— 详细日志
   • FreqPlot-*.txt —— 频率对比数据
   • Run-plink.sh —— 自动生成的 PLINK 命令脚本,把这些修正依次执行

3. 链方向 (strand) 判定:判断数据集的等位基因与参考面板是正链一致、反链一致,还是 ref/alt 反了,分 4 种情况处理(该 flip 的 flip、该换 ref/alt 的换)

4. 频率差质控:

   • 删除参考面板中 MAF > 40% 的回文 SNP(A/T、G/C 型),因为这类 SNP 链方向无法可靠判断
   • 删除数据集与参考面板等位基因频率差 > 阈值(默认 0.2,可用 -t 调)的 SNP

5. 处理特殊位点:

   • Indel(插入/缺失)—— 默认直接排除
   • 同一位置不同名字 / 位置错配的 SNP —— 按 rs ID 纠正位置或染色体
   • 重复位点 —— 排除
   • 23 号(X 染色体)之后的 XY、Y、MT —— 跳过/排除

6. 生成 PLINK 执行脚本 Run-plink.sh,按顺序执行:排除 SNP → 改染色体 → 改位置 → 翻转链 → (可选)改 ID → 强制 ref 等位基因 → 按染色体拆分,最终产出与参考面板对齐的、可以直接提交 imputation 的 PLINK 数据。

运行要求

• 命令行:perl HRC-1000G-check-bim.pl -b <bim> -f <frq> -r <参考面板> -h(HRC)或 -g -p <群体>(1000G)
• 输入:PLINK 的 .bim 文件 + 由 plink --freq 生成的 .frq 频率文件
• 约需 ~20GB 内存

简单说:它就是 imputation 前必做的"strand/坐标/等位基因/频率对齐检查器",输出一组 plink 更新清单和一条可执行的 plink 命令脚本。

为什么最终生成的 bim 中存在位点 ilmnseq_1:754198在HRC.r1.GRCh37.autosomes.mac5.sites.tab 中不存在，这样正常吗

Searched for text ilmnseq_1:754198, no results

Searched for files matching **/*.bim, no matches

Searched for files matching **/*sites.tab*, no matches

Searched for files matching **/*Exclude*, no matches

Searched for files matching **/*ID-*, no matches

Read

正常情况下,这是预期行为,不一定是错误。关键在于这个脚本的匹配逻辑:优先按「染色体:位置」匹配,而不是按 ID 匹配。原因分两步解释:

1. 为什么它没被排除

脚本主循环里,每个 SNP 只有三种归宿(见 while (<IN>) 部分):

ilmnseq_1:754198 这种 ID 在 HRC 里确实不存在(HRC 用的是 rs... 或 .)。但只要 HRC 在 chr1:754198 这个位置上有变异位点(ID 是 rs... 或 .),脚本就走第一分支:

• $id{'1-754198'} 命中 → 位置匹配
• $id{'1-754198'} ≠ ilmnseq_1:754198 → 判定为「ID 不一致」
• 又因为 ilmnseq_1:754198 不是 rs ID,$rs{$temp[1]} 不存在,不走「错配 SNP 纠正」分支
• 于是进入 $idmismatch++ 分支:写一行到 ID-*.txt(ilmnseq_1:754198 → HRC 的 ID),但 SNP 本身被保留

2. 为什么最终 bim 里还是旧 ID

脚本默认关闭了 ID 更新这一步。在生成 Run-plink.sh 的代码里,--update-name 这 4 行是注释掉的:

#update ids
#remove the following 4 lines if you want don't want to update the SNP identifiers to match the HRC
#print SH "$plink --bfile $tempfile --update-map $idfile --update-name --make-bed --out ";

所以后续 plink 流程只做了:排除 → 改染色体 → 改位置 → 翻转链 → 强制 ref 等位基因,而没有把 ID 换成 HRC 的 ID。最终 bim 里保留的仍是原始的 ilmnseq_1:754198 这个名字,尽管该位点在 HRC 里对应的是另一个 ID。

结论

• 正常。说明这个位点「位置」在 HRC 里存在(只是 ID 不同),等位基因/链方向也通过了检查,所以被保留了。
• 对后续 imputation 来说没有影响——imputation 服务器主要是按 chr:pos:ref:alt 匹配,ID 不一致通常不阻碍。

怎么确认(如果要)

1. 看 ID-<你的文件名>-HRC.txt:如果里面有 ilmnseq_1:754198 → rsXXXX 这样的行,就证明是上面这种情况(位置命中、ID 不同)。
2. 看 Exclude-<文件名>-HRC.txt:如果它出现在这里却还在最终 bim 里,那才说明 Run-plink.sh 的 --exclude 步骤没真正执行/执行顺序有问题。
3. 如果你确实希望最终 ID 与 HRC 完全一致,把 Run-plink.sh 里那 4 行注释打开(或手动加一步 plink --bfile ... --update-map ID-*.txt --update-name)。