03 · 谐化与定向过滤

对应模块：R/05_harmonise.R参考：Hemani 2018（谐化）；Hemani 2017（Steiger） 实测核查：tools/verify_harmonise.R（2026-08-04 新增，见本页第六节） 文献依据：Burgess 2023 指南第 4 节 "Variant harmonization" 等四条原文，见第一·五节

零、先讲清楚基础：什么是「效应等位」，为什么必须对齐

这一节写给第一次接触 MR 的人

如果你已经熟悉 harmonisation，可以直接跳到第六节的实测结果。

1. 一个位点上有两种碱基

人类是二倍体，同一个位置上不同人可能带不同的碱基。比如某个位点 rs123，有人是 A、有人是 G。这两种就叫这个位点的两个等位基因（allele）。

2. 「效应等位」只是一个记账的参照物

GWAS 报告某个位点的效应时，必须说清楚「相对于谁」。于是约定挑一个等位作为效应等位（effect allele, EA），另一个叫非效应等位（other allele, OA），然后报告：

每多带一个效应等位，性状变化 β 个单位。

举例：rs123，EA=A，OA=G，β=+0.3 → 意思是「每多一个 A，蛋白水平高 0.3 个标准差」。

关键点：谁当效应等位是人为选的，不是天然属性。 同一个位点，A 库可能选 A 当 EA，B 库可能选 G 当 EA。

3. 换一个效应等位，β 就要变号

还是 rs123。如果改选 G 当效应等位：

	效应等位	β
原来	A	+0.3
换成 G 当 EA	G	−0.3

数值一样，符号相反。因为「多一个 A」和「多一个 G」在二倍体上是此消彼长的。

4. 不对齐会发生什么：方向彻底反过来

MR 的核心算式非常简单（单个工具时叫 Wald 比值）：

因果效应 = β(工具→结局) ÷ β(工具→暴露)

现在假设：

暴露库（UKB-PPP 测蛋白）：rs123，EA=A，β = +0.3（A 让蛋白升高）
结局库（FinnGen 测疾病）：rs123，EA=G，β = +0.2（G 让疾病风险升高）

如果不谐化直接相除：

0.2 ÷ 0.3 = +0.67   → 结论：蛋白升高 → 疾病风险升高（有害）

但这是错的。结局库的 EA 是 G，要先翻成 A：β 变成 −0.2。正确算法：

(−0.2) ÷ 0.3 = −0.67  → 结论：蛋白升高 → 疾病风险降低（保护）

后果

有害和保护，完全颠倒。 而且 p 值、置信区间宽度都不变——从统计量上看不出任何异常，图也画得出来，只有方向是反的。这就是为什么谐化是 MR 里最不能省、也最容易翻车的一步。

5. 「谐化」（harmonisation）就是干这件事

把结局侧的效应等位统一翻成和暴露侧一样，该变号的变号。翻完之后两边的 effect_allele 必须逐条相等。

6. 翻号时的三条纪律

翻号不是只改 β 一个数。三样必须一起动，两样绝不能动：

要跟着翻	不能翻
`effect_allele` / `other_allele` 标签	`se` —— 标准误与参照物无关，恒非负
`eaf` → `1 − eaf`	`p` —— 不变
`beta` → `−beta`

且只翻一侧。 Hemani 2018 原文： "switch the direction of the effect in either the exposure or outcome study" ——两侧都翻等于没翻。TwoSampleMR 的约定是翻结局侧。

本流水线实测（APOE rs429358，暴露 EA=C／结局 EA=T）：

	效应等位	eaf	beta
源文件 Mahajan	T	0.85	+0.093
谐化后我们存的	C	0.15	−0.093

三样一起翻了，且 EA_exp = C = EA_out，两侧对齐。

只翻 β 而忘了翻 eaf 会造出什么

会得到「β 相对 C、频率却是 T 的」这种内部不一致—— 这正是 T1D 源那个故障的形状（源文件出厂时 EA 与 eaf 就对不上，谐化再正确也救不回来，见第八节）。

7. 什么情况下翻号是错的

不是规则错，而是前提不成立——这时该做的是丢弃，不是翻号：

情形	为什么翻号无效	谁来挡
rsID 指向了不同变异（rsID 合并/拆分历史）	两边根本不是同一个位点	⚠️ TwoSampleMR 不检查，是残余风险
多等位位点	暴露的 G/T 与结局的 T/G 可能指不同的 alt	「等位集合必须一致」间接挡掉
回文 SNP 且 MAF ≈ 0.5	字母无法区分「顺序反」与「链反」	`action=2` 自动剔除
等位集合不同（如 G/C vs G/A）	无论怎么翻都对不上	`harmonise_data` 判 `mr_keep=FALSE`

本课题实测：EA 错位 0、mr_keep 后无残留模糊回文（见第六节）。

一句话

真正容易出错的不是「要不要翻」，而是 ① 翻的时候有没有把 eaf 一起翻、② 输入数据的 EA 与 eaf 本身是否自洽。本课题踩的是第 ② 个。

一、还有一类更麻烦的：回文 SNP

什么是回文

DNA 是双链互补的：A↔T、C↔G。

大多数位点靠等位字母就能判断是否需要翻转。比如暴露侧是 A/G、结局侧是 G/A——一眼看出是同一对等位、只是顺序反了，翻符号即可。

但如果某个位点的两个等位恰好是 A/T 或 C/G，问题就来了：

情况	暴露侧	结局侧	怎么解释
可能一	EA=A	EA=T	顺序反了 → 需翻符号
可能二	EA=A	EA=T	两个库用了互补链记录，其实是同一个等位 → 不该翻

光看字母无法区分这两种情况。 这类位点叫回文 SNP（palindromic SNP），也叫 ambiguous SNP。

怎么解决：用等位频率

如果这个等位在人群里的频率明显偏离 0.5（比如 0.2），就能判断：

暴露侧 EA=A，频率 0.2
结局侧 EA=T，频率 0.2 → 频率一致 → 说明是互补链记录的同一个等位，不该翻
结局侧 EA=T，频率 0.8 → 频率互补 → 说明是顺序反了，该翻

但频率接近 0.5 时就没救了

如果频率是 0.48 vs 0.52，两种解释都说得通——无法判断。这类叫模糊回文（ambiguous palindromic），标准做法是直接丢掉，不要冒险。

本流水线的策略（harmonise_action: 2）：能用频率判断的回文就判断，判断不了的显式剔除。

一·五、文献依据：以上做法出自哪里

前面两节讲的都不是本项目的发明。MR 领域的协议级指南里有专门一节讲这件事。

最权威的一条：Burgess 指南第 4 节

Burgess S, et al. Wellcome Open Res 2019;4:186（v3, 2023 update） 全文分十节，第 4 节标题就叫 "Variant harmonization"。原文：

"Genetic associations with exposures and outcomes are typically reported per additional copy of a particular allele. Hence, when combining summarized data on genetic associations, it is important to ensure that genetic associations are expressed per additional copy of the same allele."

"…one dataset may report the association per additional copy of the A allele, and another per additional copy of the T allele – but the same comparison is being made."

"Allele and strand information can be double-checked by comparing allele frequency information – if the allele frequencies are similar for the A and T alleles, then the researcher can be more confident that this is a strand mismatch."

"Additional care should be taken for palindromic variants – if the alleles were A and T (or C and G), then the same alleles would appear on both the positive and negative strands. In such a case, if the allele frequency is close to 50%, analysts may choose to drop the variant from the analysis if it is not possible to verify that the alleles have been correctly orientated."

指南自己给出的后果

"While this is a conservative policy, allele alignment problems have led to incorrect results in Mendelian randomization analyses, and retractions and corrections of manuscripts."

这是「为什么必须逐条核查」最硬的一句出处。

其余三条

文献	原文要点
*Hemani et al. 2018 eLife* 7:e34408**（PMC5976434）"Harmonising exposure and outcome SNP effects" 节	顺序相反的情形 "are harmonised by flipping the sign* of the SNP-outcome effect"；回文 SNP "If reference strands are unknown, effect allele frequency can be used to resolve the ambiguity"*
TwoSampleMR harmonise vignette	`action=1` "Assume all alleles are presented on the forward strand"；`action=2` "Try to infer the forward strand alleles using allele frequency information"；`action=3` "…drop all palindromic SNPs"
*STROBE-MR（Skrivankova 2021 JAMA* 326:1614）** Item 4c 解释部分	必须报告 "the presence or absence and handling of strand alignment, and orientation of effect and non-effect alleles"
GWAS-SSF v1.1.0（规范 Table 1）	`effect_allele` = "Allele associated with the effect"（Mandatory）；`effect_allele_frequency` = "Frequency of the effect allele"（Mandatory）

本流水线每一步的依据对照

本流水线的做法	依据	性质
效应等位必须两侧一致	Burgess 2023 §4；Hemani 2018	✅ 指南
顺序相反 → 翻号保留，不是丢弃	Hemani 2018 原文	✅ 文献
`harmonise_action: 2`（用频率推链）	TwoSampleMR vignette；Burgess 2023 §4	✅ 指南＋包文档
回文且频率接近 0.5 → 丢弃	Burgess 2023 §4「close to 50%…drop」；vignette 的 MAF>0.42	✅ 指南
用等位频率交叉核对方向（`verify_harmonise.R`）	Burgess 2023 §4「double-checked by comparing allele frequency information」	✅ 指南明确推荐
以 1000G 为第三方逐位点裁决（`verify_input_sanity.R`）	同上思路的加强版	⚠️ 指南只说「比频率」，用外部面板做第三方是我们的加强
无 rsID 数据集用排序等位键匹配	—	⚠️ 无可引文献，是工程选择，依据为实测（1,656 vs 956）
判定「整列反转」的阈值（反转率 ≥50%、\|Δfreq\|<0.05）	—	⚠️ 经验值，无出处

为什么要专门列这张表

写 Methods 时，有指南出处的可以直接引，没出处的必须写明是本研究的操作定义。两者混为一谈是审稿人最容易抓的点。相关判定标准的完整分层见与导师代码逐条对比 · 第十一节。

二、转成 TwoSampleMR 格式

to_tsmr_exposure() / to_tsmr_outcome() 把内部标准 schema 喂给 TwoSampleMR::format_data()，分别生成 exposure / outcome 格式。结局只取工具 SNP 那几行。F 统计量、R² 这类非标准列会被 format_data 丢弃，所以代码额外 merge 回来供报告用。

三、谐化 + 回文处理（`harmonise_pair`）

调 TwoSampleMR::harmonise_data(action = ...)。harmonise_action 控制回文 SNP 策略：

`action`	策略	代价
1	假设所有等位已正向对齐	最激进，回文位点可能方向错
2（本课题默认）	用等位频率推断回文 SNP 的链方向	平衡；频率接近 0.5 的仍需另行剔除
3	直接丢弃所有回文 SNP	最保守，白白损失约 12% 的工具

谐化后再显式剔除残留的 palindromic & ambiguous（模糊回文）SNP，然后只保留 mr_keep == TRUE 的变异。

为什么用了 action=2 还要再剔一遍：harmonise_data 的 action=2 会给出 palindromic 与 ambiguous 两个标记，但不一定把它们从结果里删掉。多这一道显式剔除是保险，不依赖上游包的默认行为。

四、Steiger 定向过滤（`apply_steiger`）

MR 假设工具先影响暴露、再影响结局。若某 SNP 与结局的关联其实强于与暴露的关联，说明方向可能反了。TwoSampleMR::steiger_filtering() 比较每个 SNP 对暴露、对结局的 R²，剔除 steiger_dir == FALSE（方向错误）的 SNP。

Steiger 需要暴露、结局双方的 R²，所以必须在谐化之后（拿到配对数据）才能做。由 config 的 iv.steiger_filtering 开关（默认开）。

一个已知坑

二分类结局若不声明 units、ncase、prevalence，steiger_filtering() 会把对数几率当连续量算 R²，导致 rsq.outcome 全为 NA、定向过滤形同虚设。本流水线已在 to_tsmr_outcome() 里补齐这三项，并在 rsq.outcome 全 NA 时打警告。

五、单对编排（`build_harmonised`）

IV → exposure fmt / outcome fmt(仅工具SNP) → harmonise_data → 剔模糊回文 → mr_keep → Steiger

若结局与暴露无重叠 SNP，返回 NULL 并跳过该对（错误隔离，不影响其它对）。谐化后的 IV 数（nrow(h)）直接决定下一步选哪套 MR 方法。

六、本课题的实测核查结果

为什么要实测

「代码里调了 harmonise_data()」只能说明流程写对了，不能说明结果真的对齐了。方向错误在统计量上看不出来，所以必须直接检查产物。以下为 2026-08-04 实测。

6.1 效应等位是否逐条对齐

结局	工具数	EA 暴露 == 结局	EA 暴露 ≠ 结局	`b` == β结局/β暴露最大差
2 型糖尿病（Mahajan）	1,614	1,614	0	5.55×10⁻¹⁶
糖尿病视网膜病变	1,587	1,587	0	3.77×10⁻¹⁵

谐化后没有一条残留错位。 Wald 比值与「结局 β ÷ 暴露 β」吻合到浮点精度，说明 b 确实是用谐化后的 β 算的。

6.2 等位频率的独立佐证

若有等位没对齐，那些位点的 eaf 会变成 1 − eaf，相关系数会被明显拉低。

结局	eaf 暴露 ↔ 结局相关
2 型糖尿病	0.9991
糖尿病视网膜病变	0.9841

★2026-08-04 更新：T1D 支线未通过

上表只列了通过的结局。外部 T1D 源 GCST90824163 的等位频率相关为 −0.5931，审计报错。 根因是该文件的 effect_allele_frequency 列整列反转，详见第八节。 该问题尚未修复；修完前 T1D 相关的方向结论不可对外使用。

6.3 回文 SNP：最危险的那一类是零

结局	工具数	回文（A/T、C/G）	占比
2 型糖尿病	1,614	205	12.7%
糖尿病视网膜病变	1,587	195	12.3%
糖尿病黄斑病变	1,587	195	12.3%

保留下来的 205 个回文 SNP，其暴露与结局的 eaf 相关 0.9991、最大绝对差仅 0.034——说明它们的链方向是靠频率正确判定的，不是蒙对的。

结论：回文这条最容易翻车的路径，本课题没有留下任何模糊病例。

6.4 全数据源列映射审计（2026-08-05 新增）

前三项查的是「谐化后对不对」。这一项往上游一层，查每个数据源的哪一列被当成了效应等位—— 因为列认错了，后面全对齐也是对齐到错的东西上。

对 outcome_manifest.csv 全部 45 行逐个读表头、跑 detect_columns()，打印解析结果：

数据源族	表头	EA / OA / EAF 解析为	判定
FinnGen R9 / R13（36 行）	`ref alt … af_alt`	`alt` / `ref` / `af_alt`	✅ FinnGen 的 alt 即效应等位
GWAS-SSF 族（6 行，含 T1D、MVP 复制队列）	标准字段名	`effect_allele` / `other_allele` / `effect_allele_frequency`	✅ 名字对（T1D 的内容另见第八节）
Mahajan T2D（2 行）	`EA NEA EAF`	`EA` / `NEA` / `EAF`	✅
IAMDGC full	`ALLELE1 ALLELE0 A1FREQ`	`ALLELE1` / `ALLELE0` / `A1FREQ`	✅
Salem2019 糖肾（1 行）	`Allele1 Allele2 FreqA1`	别名表未覆盖 → OA/EAF 落空	⚠️ 但调用处有显式 `column_map`，不受影响

配套的输入层自洽性检查（各取前 20 万行）：

检查	结果
Salem2019：`Beta` 是否等于 `log(OR)`	最大差 4.9×10⁻⁶，相关 1.000000 → `effect_type="beta"` 正确，未重复取对数
五个 FinnGen 端点：`β/se` 重算 p 对报告 p	中位 \|Δlog₁₀p\| ≈ 4.5×10⁻⁷
T1D_gcst	中位 1.9×10⁻⁵
T2D Mahajan	中位 3.8×10⁻³，最大 0.597 —— 系发表文件只保留 2–3 位有效数字所致（`Beta=0.071, SE=0.012`），非错误
暴露 `log10(p) (discovery)` 编码方向	取值全为正（10.8–5369.5）＝ −log₁₀p，与 `pval_encoding: neglog10` 一致；对照重算相关 0.9952

6.5 全结局 × 1000G 第三方裁决（2026-08-05 新增）

6.2 那张表比的是「我们自己的暴露 eaf ↔ 我们自己的结局 eaf」，属内部一致性。本项对每一个在用结局都以 1000G EUR 做外部裁决（tools/verify_eaf_all_outcomes.R，工具位点 1,875 个）：

release	结局	可比对	反转率	cor(eaf, 1000G)	判定
R9_dm	Retinopathy / Maculopathy / NeovascGlaucoma / Neuropathy / Nephropathy	1,210	2.2–2.3%	+0.983	✅
R9_dm / R13_dm	T1D_gcst	1,277	93.7%	−0.9922	❌ 整列反转
R13_dm	RetinopathyStrict / Maculopathy / NeovascGlaucoma / Neuropathy / Nephropathy	1,232	2.2–2.3%	+0.983	✅
R13_amd	AMD / WetAMD / DryAMD	1,232	2.2%	+0.983	✅
R13_reserve	AMD_iamdgc / AMD_iamdgc_full	1,252	1.4%	+0.994	✅
R13_repl	DR_mvp / Neuropathy_mvp	1,280	1.6%	+0.997	✅
R9_sens	T2D_mvp	1,280	1.6%	+0.997	✅
R9_sens	T1D_crouch（零重叠复制源）	1,030	2.0%	+0.996	✅
R13_dm	T2D_mahajan（无 rsID，按 chr:pos 另跑）	1,228	0.8%	+0.998	✅
R13_repl	Nephropathy_dkd Salem2019（METAL 列名，显式指列另跑）	1,293	1.6%	+0.998	✅
暴露	UKB-PPP	1,305	1.2%	+0.9928	✅

24 个数据源，只有 T1D_gcst 一个不通过（它在 R9_dm 与 R13_dm 各出现一次，是同一个文件）。其余的 1–2% 「反转」是人群频率差异与阈值边界造成的正常噪声。

查出的两处结构性隐患（当前未触发，但应修）

maf 被列进了 eaf 的别名表（config.yaml 的 column_aliases）。 MAF 是次等位频率，效应等位若为主等位就是错的——与第八节的 T1D 故障同一类型。当前所用文件都没有 maf 列，故未触发；换数据源就可能踩上。
pval 与 mlogp 同在 pval 别名表。FinnGen 两列都有，靠别名顺序（pval 在前）取对了。顺序一动就全错，且 p 值错了不会报任何错。

两条的共性：靠列表顺序侥幸对着，不是靠断言保证对着。

6.6 本轮未覆盖的部分（如实声明，不要当成已查）

未查项	影响什么	现有间接证据
多组学项目的视网膜 eQTL / sQTL 输入（独立目录 `multiomics-mr`）	ERMAP「视网膜方向与血浆相反」这一条	无。ERMAP 本就标注为「待验证」
PheWAS（OpenGWAS）输入	方向感知药靶 PheWAS	该分析两侧都用各自的 Wald 比值，不依赖单一 eaf 列
deCODE 的 eaf 列	—	该文件报的是 `ImpMAF`（次等位频率），本来就没有 EAF 可查；方向仅靠等位字母

另有一条未加断言的隐含假设：共定位 / 复制 / eQTL 脚本里的 flip 辅助函数（analysis/05_hyprcoloc.R、08_replicate_iamdgc.R、10_eqtl_coloc.R、26_coloc_susie.R、 28_decode_reselect.R、35_reverse_mr.R）只按等位字母配对，没有回文/链处理，等价于 action = 1（假设各源均为正向链）。

实测佐证：各结局的回文子集 eaf 与暴露高度一致（如 T2D n=205 相关 0.9991、最大差 0.034），支持「所有源均为正向链」
边界已量化：三个候选的工具全是非回文（IFNAR1 rs914142 A/G、ERMAP rs11210710 T/G、 APOL1 rs136168 A/G）；deCODE 重选的 14 个工具里只有 1 个回文（APOL1 rs377486200 A/T，MAF 0.082，且 APOL1 已因跨平台方向矛盾出局）
28_decode_reselect.R 是唯一自带回文保护的（palin && maf > 0.42 剔除），且 se <- abs(m$se_o / iv$beta) 取了绝对值

七、与导师实现的差别：这是「匹配」问题，不是「谐化」问题

这两步经常被混为一谈，必须分开：

步骤	做什么	做错的后果
匹配（matching）	在结局文件里找到这个变异	找不到 → 工具丢失 → 功效下降
谐化（harmonisation）	找到之后对齐效应等位	没对齐 → 方向反了

7.1 涉及的三个原始文件（可自行点开核对）

#	文件	大小	说明
1	`D:\mrdata\exposure\protein_info.csv`	4.06 MB	UKB-PPP cis 工具表（WSL 的 `data/exposure/protein_info.csv` 是指向它的软链）
2	`D:\mrdata\outcome\Mahajan.NatGenet2018b.T2D-noUKBB.European.txt`	1.12 GB	外部 T2D 结局，2150 万行
3	`F:\project\MR code\2_diabete_validation.Rmd`	15.9 KB	导师代码，相关在 L86 / L88 / L96

文件 2 太大不能用编辑器打开，看某一行用：

bash

sed -n '20272186p' "D:\mrdata\outcome\Mahajan.NatGenet2018b.T2D-noUKBB.European.txt"

7.2 三行代码逐行拆解

L86 —— 把暴露侧的变异 ID 削掉尾巴

ready_exposure_cis$clean_variant <- sub(":imp:v1.*", "",
    ready_exposure_cis$Variant.ID..CHROM.GENPOS..hg37..A0.A1.imp.v1.)

片段	意思
`ready_exposure_cis`	暴露表（UKB-PPP cis 工具）
`$Variant.ID..CHROM...`	取其中一列。这个怪名字是 R 自动改的——原列名 `Variant ID (CHROM:GENPOS (hg37):A0:A1:imp:v1)` 里的空格、括号、冒号全被 R 换成了点
`sub(a, b, x)`	在 x 里把第一处匹配 a 的地方替换成 b
`":imp:v1.*"`	正则式：`:imp:v1` 后面跟任意字符
`""`	替换成空＝删掉
`df$新列 <- ...`	新建一列，取名 `clean_variant`

效果：19:45411941:T:C:imp:v1 → 19:45411941:T:C

⚠️ 保留的是 UKB-PPP 原本的写法 chr:pos:A0:A1，即非效应等位在前。

L88 —— 给 Mahajan 也拼一个键

temp_t2d$clean_variant <- paste0(temp_t2d$Chr, ":", temp_t2d$Pos, ":",
    temp_t2d$NEA, ":", temp_t2d$EA)

片段	意思
`temp_t2d`	Mahajan 全库，21,508,698 行
`paste0(...)`	把几段字符串首尾相接，中间不加分隔符（`paste0` ＝ `paste` 且 `sep=""`）

效果：19 + : + 45411941 + : + C + : + T = 19:45411941:C:T

⚠️ 这里的顺序是手写死的：NEA 在前、EA 在后。

L96 —— 用键做筛选

selected_outcomes_diabetetp2 <- subset(temp_t2d,
    clean_variant %in% ready_exposure_cis$clean_variant)

片段	意思
`%in%`	判断左边每个值在不在右边那个集合里，返回一串 TRUE/FALSE
`subset(df, 条件)`	只保留条件为 TRUE 的行

三行连起来的隐含要求

L86 保留 A0:A1，L88 强制写成 NEA:EA。两个字符串要相等，就必须

UKB-PPP 的 A0 恰好等于 Mahajan 的 NEA，且 A1 恰好等于 EA

代码里没有任何一处检查这个前提。

7.3 具体例子甲：APOE rs429358（被漏掉）

暴露 protein_info.csv 第 708 行（rs429358 全文出现 25 次，其余都是别的蛋白的 trans，只有此行 cis/trans = cis）：

列	值
`Variant ID (CHROM:GENPOS (hg37):A0:A1:imp:v1)`	`19:45411941:T:C:imp:v1`
`rsID`	rs429358
`Assay Target`	APOE
`A1FREQ (discovery)`	0.1557
`BETA (discovery, wrt. A1)`	−1.012

→ A0（非效应）= T，A1（效应）= C；BETA 列名自带 "wrt. A1"，即 −1.012 是相对 C 的。

结局 Mahajan 第 20,272,186 行：

SNP           Chr  Pos        EA  NEA  EAF   Beta   SE     Pvalue
19:45411941   19   45411941   T   C    0.85  0.093  0.012  2.5e-14

→ 效应等位是 T，+0.093 是相对 T 的。

两边各自挑了不同的等位当参照——这是常态，不是错误。

怎么知道效应等位是 T？——四层证据，只信第一层是不够的

这个问题对每一个数据源都要回答一遍，不能想当然。

① 列名自证（最直接，但最弱）

表头写的是 EA / NEA ＝ Effect Allele / Non-Effect Allele。这是文件自己的声明。

② 官方数据字典（最强的文档依据）

F:\project\Mahajan.et.al.2018b.European.GWAS.readme.v2.pdf（DIAMANTE 官方发布，2022-03-05）第 2 页原文：

In the summary files, for each SNP, we have provided the following information:
Chromosome and position (build 37, base-pairs).
Effect (EA) and non-effect allele (NEA), aligned to the forward strand.
Effect allele frequency (EAF).
Log-odds ratio for the effect allele (Beta) and the corresponding standard error (SE).
P-value for association (Pvalue).

一次说清四件事：

问题	readme 的答复
EA 是不是效应等位	是
`Beta` 相对谁	"Log-odds ratio for the effect allele" → 相对 EA
`EAF` 是谁的频率	effect allele frequency → EA 的
哪个在前	字段列举顺序与实际列序 `SNP Chr Pos EA NEA EAF Beta SE Pvalue` 一致 → EA 在前

另有一条我们原先只能靠实测推断的："aligned to the forward strand"——官方声明已对齐正向链。这解释了为什么本课题实测那 4 个对不上的位点里反向互补 0 个。

这份 readme 不在数据目录里

它在 F:\project\ 下，而数据在 D:\mrdata\outcome\； .zip 解包后只有一个 txt，不含 readme。 换机器/交接时极易丢失——建议把它复制一份到数据目录旁边。

即便如此，③④ 两层仍然要做：readme 只能证明作者的意图，不能证明文件内容与意图一致 （T1D 源 GCST90824163 声明遵循 GWAS-SSF，实际 effect_allele_frequency 整列反转，见第八节）。

③ 频率对外部参照

若 EA=T，则 EAF 应当是 T 的频率。用 1000G EUR 实测：

位点	Mahajan 报的 EA / EAF	1000G EUR 实测	判定
rs429358（本例）	T / 0.85	C=0.1551 → T=0.8449	✅ 对上
rs7903146	T / 0.30	T=0.3171	✅ 对上

全库层面已跑过（见 6.5）： Mahajan cor(eaf, f_1000G) = +0.9981、反转率 0.8%。

④ 哨兵位点方向硬验（最强的一层）

挑一个方向已成教科书共识的位点，看数据是否复现它。

TCF7L2 rs7903146 是 2 型糖尿病最著名的位点，T 是风险等位 （Grant 2006 Nat Genet 38:320 首报，此后被反复复制；文献 OR ≈ 1.35–1.40）。

Mahajan 第 15,393,081 行附近该位点原始行：

SNP            Chr  Pos         EA  NEA  EAF  Beta  SE     Pvalue
10:114758349   10   114758349   T   C    0.3  0.3   0.009  4.5e-238

EA = T、Beta = +0.30 → exp(0.30) = 1.35
与文献共识的 OR 方向一致、量级一致
p = 4.5e-238 —— 这是全库最强的信号之一，与 TCF7L2 的地位相符

④ 一次同时验证了三件事：EA 列的语义、Beta 的参照物是谁、以及符号约定。仅靠①（列名）无法排除「作者把 EA/NEA 写反了」这种情况；③④ 能。

这正是我们后来加进审计的东西

tools/verify_input_sanity.R 的 C 项做的就是③（频率对 1000G）。 ④（哨兵位点硬断言）目前尚未写成代码，是待补项—— 本页 8.6 记录的 T1D 故障，如果当时有④会更早暴露。

键的比较：

L86 暴露侧 → "19:45411941:T:C"
L88 结局侧 → "19:45411941:C:T"        (paste0(19,":",45411941,":",NEA="C",":",EA="T"))

"19:45411941:T:C" == "19:45411941:C:T"   →   FALSE

L96 的 %in% 判 FALSE → 这一行不进 selected_outcomes_diabetetp2。R 不报任何错， subset() 只是返回一个更小的表。

7.4 对照例子乙：MANSC4 rs34311349（能命中）

暴露：12:27925037:T:A（A0=T，A1=A，BETA wrt A1 = +0.938）结局 Mahajan 第 15,393,081 行：

12:27925037   12   27925037   A   T    0.22  -0.069  0.01   5.9e-12

L86 暴露侧 → "12:27925037:T:A"
L88 结局侧 → "12:27925037:T:A"        (NEA="T", EA="A")
                                       →   TRUE，命中

差别仅在于 Mahajan 这一条恰好把 A 当效应等位，与 UKB-PPP 的 A1 撞上了。

例子丙：PAM rs149802978，暴露 5:102346866:C:G，Mahajan 第 6,648,000 行EA=C, NEA=G → 结局键 5:102346866:G:C → 不相等，同样漏掉，而它在 T2D 端 p=4.6e-15。

7.5 实测统计

1,955 个 cis 变异 × Mahajan 全库 21,508,698 个变异，按 chr:pos 能对上的有 1,660 个：

情况	数量
顺序一致（`A0==NEA` 且 `A1==EA`），他能命中	956
顺序相反，他静默丢掉	700
等位不匹配 / 多等位	4

丢 700 个，占 700/1656 = 42.3%。 我们按 rsID 匹配，得到 1,614 个 T2D 工具。

7.6 同一个例子，讲清「效应等位必须一致」

还是 APOE rs429358。两边的话是：

暴露：每多一个 C，APOE 蛋白 −1.012 SD
结局：每多一个 T，T2D log-odds +0.093

参照物不同，不能直接相除。好比一个人说「往东 1 km」，另一个说「往西 3 km」，要先统一方向。

统一到 C（暴露侧的效应等位）：

「每多一个 T，+0.093」 ⟺ 「每多一个 C，−0.093」

二倍体上多一个 C 就等于少一个 T，所以符号翻转。现在两边都以 C 为参照：

Wald 比值 = (-0.093) / (-1.012) = +0.0919   →   OR = 1.096
解读：每升高 1 SD APOE，T2D 风险【升高】

若不谐化、直接拿原始 β 相除：

b = (+0.093) / (-1.012) = -0.0919   →   OR = 0.912
解读：每升高 1 SD APOE，T2D 风险【降低】

同一个位点、同一份数据，结论从有害翻成保护。 PAM 同理：正确 OR 0.872（保护），不谐化得 OR 1.146（有害），方向相反。 MANSC4 因为两边效应等位本来就都是 A，翻不翻号结果一样（OR 0.929）。

为什么这一步最容易翻车

方向错了以后，p 值、置信区间宽度、森林图全都一模一样——从统计量上看不出任何异常。这正是本流水线要写 tools/verify_harmonise.R 逐条核对的原因（见第六节）。

7.7 原文的注释实况

2_diabete_validation.Rmd 在这几行是有注释的，但只写了「做什么」：

 85| #generate new variable in exposure for MR
 87| #generate new vaiable in temp_t2d          ← 原文即拼错的 "vaiable"
 95| #select SNPs from the whole database

没有任何一条提到「两边等位顺序必须一致」这个前提。

更关键的是筛完之后没有任何计数检查。全文检索 nrow / length / dim：

L56 length(unique(ready_exposure_cis$SNP)) —— 数了筛选前的暴露工具数
L96 之后 一个数都没数

所以「进去 1,955 个、出来 956 个」这件事，在他的分析流程里从头到尾没有出现过任何数字。

7.8 要说公道话的地方

L118 他是做了谐化的：

118| harm_rt_diabetetp2 <- harmonise_data(exposure_dat = ready_exposure_cis,
         outcome_dat = selected_outcomes_diabetetp2, action = 2)

且 L108 / L115–116 把两边的 clean_variant 都改名成 SNP，让 harmonise_data 按这个复合键去配。

由此产生一个有意思的结果：既然 L96 已经把「顺序不一致」的行全筛掉了，活下来的 956 条两边等位顺序天然相同——harmonise_data 拿到手时根本没有需要翻号的行。谐化这一步在他的流程里等于空跑。

7.9 那我们是怎么做的

两条路径，取决于结局文件有没有 rsID。

	路径一	路径二
适用	有 rsID（R9 七个结局里的六个）	无 rsID（只有 Mahajan T2D）
键	rsID	`chr:pos:排序后的两个等位`
方向	交给 `harmonise_data(action=2)`	同左

路径一 · 按 rsID 匹配（默认）

yaml

# config.yaml
 20|  match_by: "rsid"       # rsid(默认, build-safe) | position(无rsID数据集,如Mahajan T2D)
 23|  require_rsid: true     # 强制两侧都有 rsID，否则报错(position 模式忽略)

rsID 是位点的身份标识，与等位书写顺序无关、与基因组版本也无关—— 根本不会出现导师那个问题。require_rsid: true 保证不会悄悄退化：哪一侧缺 rsID 就直接报错，而不是静默少匹配。

路径二 · position 模式＋排序等位键

# R/02_standardize.R
 47|   # 3b) 匹配键：position 模式用 chr:pos:排序等位（无 rsID 数据集，如 Mahajan T2D）
 48|   #     排序等位使键与效应等位方向无关；谐化再对齐方向 [比固定 NEA:EA 顺序更稳健]
 49|   if (identical(match_by, "position")) {
 52|     # 等位用小写：TwoSampleMR::format_data 会 tolower(SNP)，小写键才能两侧一致
 54|     d[, SNP := tolower(paste(chr, as.integer(pos),
 55|                              pmin(effect_allele, other_allele),
 56|                              pmax(effect_allele, other_allele), sep = ":"))]
 57|     d <- d[!duplicated(SNP)]
 58|   }

analysis/01_screen.R L37–45 用同一口径给暴露侧建 hg37key_srt，注释里写明了理由。

pmin / pmax 按字母序排，所以键里只剩客观信息——还是 APOE 那个位点：

暴露  19:45411941:T:C  ─┐
                        ├─→  排序后都是  19:45411941:c:t   ✅ 对上
结局  19:45411941:C:T  ─┘

对比导师的固定顺序键：19:45411941:T:C vs 19:45411941:C:T → ✗ 不等 → 丢弃。

效果实测

同一批变异（1,955 个 cis 工具 × Mahajan 全库）：

匹配键	命中
导师的 `chr:pos:NEA:EA`（方向敏感）	956
我们的排序等位键（方向无关）	1,656

排序键之后仍对不上的只有 4 个，实测分了类：

蛋白	暴露等位	结局等位	原因
RHOC	G/C	G/A	等位集合真的不同（同位置不同变异/多等位）
TINAGL1	A/G	T/A	同上
PLB1	G/C	G/A	同上
NEXN	AT/A	T/A	indel

这四个本来就不该匹配——它们不是同一个变异。

为什么方向不放进键里，而交给谐化

问题	谁负责	依据
这是不是同一个变异？	匹配键（只含客观信息）	—
两边参照一不一样？不一样怎么办？	`harmonise_data(action=2)`	Hemani 2018："harmonised by flipping the sign"

职责分开，两件事才都能做对。 合并成一个方向敏感的键，就是下一节（7.10 第三层：拿「主观约定」当成了「客观身份」）讲的那个错误。

排序等位键的残余盲区（如实声明）

排序键不处理链翻转。若同一个非回文变异，一边记 A/G、另一边记 T/C（反向互补），排序后是 a:g 与 c:t，对不上，会漏。

路径一（rsID）没有这个问题：action=2 会识别并翻转
路径二实测：Mahajan 那 4 个漏掉的里反向互补 0 个，未触发
Mahajan 官方 readme 亦声明 "aligned to the forward strand"，与实测一致

Methods 应写明：position 匹配仅用于 Mahajan，且实测无链不一致病例。

本节讲的是「我们怎么做」，上一节讲的是「他哪里错」

两节要对照着看：他的问题出在匹配层，解法也在匹配层，与谐化无关。

7.10 最终判定：他到底做错了没有

分四层回答，每一层都有实测支撑。

第一层 · 概念理解：没错

L88 把 NEA 写在前，不是文件的列序（readme 与实际列序都是 EA 在前）。他是主动反过来写的——为了对齐 UKB-PPP Variant ID 的 A0:A1 格式（A0 = 非效应等位在前）。

这说明他正确读懂了两个文件的列语义。 不是抄错列、不是不知道哪个是效应等位。

第二层 · 效应量方向：没错，且已实证

保留的 956 条方向正确，而且 L118 确实调了 harmonise_data(action = 2)。

关键问题是：被丢掉的 700 个，是不是"碰巧"是某一类？ 若是，就会引入选择偏倚。实测（丢失组 vs 保留组，先把结局 β 统一对齐到暴露侧效应等位再比）：

比较项	丢失(700)	保留(956)	检验	判定
对齐后结局 β 为正的比例	53.0%	54.4%	χ² p = 0.61	无差异
对齐后 \|β\| 中位数	0.0098	0.0110	Wilcoxon p = 0.18	无差异
结局侧 p<0.05 的比例	13.7%	11.2%	χ² p = 0.14	无差异

→ 丢失对结局效应是非差异性的（non-differential）→ 保留部分的估计值无偏。

这里我自己差点报错，记录下来

第一次比较时我直接比了未对齐的 Beta，得到 53.0% vs 46.6%、χ² p = 0.002，看起来"丢失是差异性的"。这个显著是伪的——两组的 Beta 参照等位本来就不同（保留组参照 A1，丢失组参照 A0），比的是两个不同的量。对齐后 p = 0.61。

教训：任何跨组比较效应量之前，必须先确认两组的参照等位一致。 这与本页第零节讲的是同一件事。

第三层 · 真正做错的是什么：拿「主观约定」当成了「客观身份」

一个变异身上有两种完全不同的信息

以 19:45411941 为例：

	客观事实（身份）	主观选择（方向标注）
UKB-PPP	这个位置上有 C 和 T 两种碱基	我挑 C 当效应等位
Mahajan	这个位置上有 C 和 T 两种碱基	我挑 T 当效应等位

左边一列两家完全相同——这就是同一个变异，没有任何疑问。

右边一列不同——但这只是各自的记账习惯，跟这个变异本身是什么毫无关系。就像两个人记同一笔账，一个记「收入 100」，另一个记「支出 −100」，说的是同一件事。

他的错：把两列揉进了一个键

paste0(Chr, ":", Pos, ":", NEA, ":", EA)
#      └──── 客观 ────┘  └──── 主观 ────┘

这个键里，前半截是身份，后半截是记账习惯。

于是程序比 19:45411941:T:C 与 19:45411941:C:T 时，看到不相等，得出的结论是：

❌ 「这不是同一个变异」

而真相是：

✅ 「这是同一个变异，只是两家的记账习惯不同」

他让「记账习惯不同」冒充了「不是同一个东西」。这就是错误的全部。

为什么必须拆成两步

要回答的问题	该看什么
① 这是不是同一个变异？	只能看客观信息（位置＋有哪两个碱基）
② 两边的参照一不一样？	看主观标注，不一样就翻号

一个键同时干这两件事，问题 ① 的答案就会被问题 ② 污染——参照不一样，就被误判成不是同一个变异。

Hemani 2018 明确要求这种情形应当 翻号保留（问题 ②），而不是让它决定纳入（问题 ①）。

修法：一行的事

# 他写的（把主观顺序编进了键）
paste0(Chr, ":", Pos, ":", NEA, ":", EA)
#   暴露 19:45411941:T:C    结局 19:45411941:C:T    →  ✗ 不等

# 排序一下（键里只剩客观信息）
paste0(Chr, ":", Pos, ":", pmin(NEA, EA), ":", pmax(NEA, EA))
#   暴露 19:45411941:C:T    结局 19:45411941:C:T    →  ✓ 相等

pmin / pmax 按字母序排 → 无论哪家挑谁当效应等位，键都一样。方向留给后面的 harmonise_data() 去翻。

本流水线 R/02_standardize.R L54–56 写的就是这个。同一批变异： 他的键命中 956，排序键命中 1,656。

为什么这不只是「少拿了点数据」

丢弃的判据是等位书写顺序——这件事跟蛋白、跟糖尿病、跟任何生物学毫无关系。

好的一面：正因为无关，丢失是随机的（见第二层实测）→ 留下来的估计值没有偏
坏的一面：等于闭着眼睛扔掉 42% 的工具，而且扔的时候没有任何提示—— subset() 只是返回一个更小的表，L96 之后全文没有一个 nrow

结果就是：8 个最强信号里 7 个（PAM、APOE、TSPAN8、NOTCH2、ABO…）从来没进过分析，而分析报告上看不出任何异常。

一句话

他不是算错了，是用记账习惯当身份证，把同一个变异认成了两个。

第四层 · 后果有多严重：比「丢 42%」更具体

（1）丢失不是随机的

	丢失(700)	保留(956)	检验
对齐后 EAF 中位数	0.275	0.200	Wilcoxon p < 1e-4
Mahajan 的 EA 频率 >0.5 的比例	72.6%	17.9%	—

丢失与等位频率强相关。不影响方向，但意味着保留下来的工具集不是原集合的随机样本。

（2）8 个全基因组显著的工具里，7 个被丢掉

因为单 SNP Wald 下 MR 的 p 恒等于结局 p（见第九节），这 8 个正是他 T2D 支线最强的 8 个潜在结果：

蛋白	rsID	位点	状态	T2D p	OR
PAM	rs149802978	5:102346866	★被丢弃	4.6e-15	0.872
APOE	rs429358	19:45411941	★被丢弃	2.5e-14	1.096
NCR3LG1	rs12146443	11:17384498	★被丢弃	9.3e-13	0.845
MANSC4	rs34311349	12:27925037	保留	5.9e-12	0.929
HLA-DRA	rs36096565	6:32560025	★被丢弃	7.5e-11	1.120
TSPAN8	rs11178654	12:71537329	★被丢弃	4.4e-10	1.115
NOTCH2	rs2641348	1:120437884	★被丢弃	1.4e-09	0.628
ABO	rs505922	9:136149229	★被丢弃	1.5e-09	1.046

其中 TSPAN8、NOTCH2、ABO、HLA 区都是已知的 T2D 位点。

这个 p 值是事后检验，不作推断用

Fisher 精确检验 p = 0.012，但我是先看到 7:1 这个不平衡才去做的检验，属于事后（post-hoc）分析，且事件数只有 8 个。 该数字只作描述，不能当作"丢失富集于强信号"的统计证据。 描述性事实本身（7/8 被丢）不受影响。

结论

维度	判定
他是否搞错了效应等位	否
他保留的 956 条估计值是否有方向偏误	否（实测非差异性丢失）
他是否做错了	是 —— 用方向敏感的键做匹配，且无任何计数检查
阳性结果是否可信	可信
阴性结果是否可解释	不可解释 —— 42.3% 的工具（含 8 个最强信号里的 7 个）从未被检验过

定稿措辞（合稿请照抄，不要改写成其他版本）：

该实现使用方向敏感的字符串键（chr:pos:NEA:EA）匹配变异，要求两个数据集恰好指定同一等位为非效应等位。实测 1,656 个可匹配变异中 700 个（42.3%）因等位顺序相反而被静默丢弃，其中包括 8 个结局侧全基因组显著位点中的 7 个。丢失对结局效应为非差异性（对齐后符号 p=0.61、量级 p=0.18、显著性 p=0.14）， 故保留部分的效应估计无偏；但阴性结果不可解释。 本流水线改按 rsID 匹配（无 rsID 数据集用位置＋无序等位键），方向交由 harmonise_data(action=2) 处理，同一批变异得到 1,614 个工具。

但共定位那一步没有这层保护

同样的顺序敏感字符串匹配也出现在共定位脚本里（见与导师代码逐条对比第四节缺陷 3），而共定位之后没有谐化来兜底。coloc.abf 基于 Z²、对等位方向不敏感，所以影响的仍是覆盖而非方向；但若换成对方向敏感的共定位方法，这里就会出问题。

八、已发现的真实故障：T1D 源的等位频率列整列反转

状态：2026-08-05 已修复并重跑验证

发现于 2026-08-04，修复于 2026-08-05。 cor(eaf, 1000G) 由 −0.9922 → +0.9922，verify_harmonise.R 由 27 PASS/1 FAIL 变为 28 PASS / 0 FAIL。修复方式与结果见 8.7 与判定标准与 skill 合规打分 · 第八节。

本节保留故障的完整记录——它是「输入层错误在统计量上完全隐形」的教学案例。

8.1 官方文档怎么说（这决定了失效模式）

TwoSampleMR 的 harmonise vignette 原文：

"if the outcome effect allele (A) were on the forward strand we would expect it to have a low allele frequency, but given it has a high frequency (0.91) we infer that the outcome GWAS is presenting the effect on the reverse strand."

也就是说：action = 2 判断回文 SNP 链方向的唯一依据就是等位频率。频率给反了，它的结论必然反。文档同时说明，回文 SNP 的 MAF 超过 0.42 时无法判断，会被直接丢弃 ——这也解释了为什么我们的产物里「模糊回文」为 0：它们已经被 TwoSampleMR 丢掉了。

这一列该装什么，有强制性标准可查。 GWAS Catalog 的 GWAS-SSF v1.1.0 规范 Table 1 原文：

字段	官方定义（原文）	强制性
`effect_allele`	"Column 2: Allele associated with the effect"	Mandatory
`other_allele`	"Column 3: The non-effect allele"	Mandatory
`beta`	"Column 4: Effect size as beta"	Mandatory（四选一）
`effect_allele_frequency`	"Column 6: Frequency of the effect allele"	Mandatory

而该文件随附的 GCST90824163.h.tsv.gz-meta.yaml 自己声明：

yaml

file_type: GWAS-SSF v1.0
is_harmonised: true
minor_allele_freq_lower_limit: 0.001

→ 这不是「口径不同」，是文件违反了它自己声明遵循的强制字段定义。 写 Methods 时可以这样表述，而不必用模糊说法。

8.2 实测：整列反转，93.7%

tools/verify_t1d_eaf.R 以 1000 Genomes EUR 为第三方，对全部工具位点逐个比对。只用非回文位点——这类位点的等位身份靠字母即可确定，不存在链歧义，因此「同一个等位在两个数据源里的频率」可以直接比。

判定	位点数	占比
★反转	1,196	93.7%
无法判定	59	4.6%
正确	21	1.6%

相关系数 EAF_T1D vs f_1000G = −0.9922
相关系数 EAF_T1D vs (1 − f_1000G) = +0.9922

那 21 个判为「正确」的全是频率贴近 0.5 的位点（如 0.5128 vs 0.5368），反转前后差别小于阈值、判不开而已，不是真的正确。59 个「无法判定」同理。

抽查三例与 1000G 对照：

SNP	1000G EUR	我们暴露侧（UKB-PPP）	T1D 源文件
rs9267797	T = 0.0189	T = 0.015 ✅	T = 0.985 ❌
rs641153	A = 0.0865	A = 0.096 ✅	A = 0.902 ❌
rs36096565	G = 0.1819	G = 0.212 ✅	G = 0.789 ❌

对照组：同一方法验证暴露侧（2026-08-05 补做）

只报「T1D 反了」不足以排除「是我们的比对方法有问题」。用完全相同的脚本逻辑验证暴露侧 UKB-PPP 的 A1FREQ (discovery)：

判定	位点数	占比
正确	1,254	96.1%
无法判定	35	2.7%
★反转	16	1.2%

cor(EAF_暴露, f_1000G) = +0.9928，与 T1D 的 −0.9922 恰成镜像。

→ 方法本身没问题；暴露侧干净，问题只在 T1D 那一列。 这一条同时是整条因果链最关键的一环的独立佐证：暴露的 BETA (discovery, wrt. A1) 与 A1FREQ (discovery) 均以 A1 为参照，而 A1 正是我们取的 effect_allele。

8.3 反转来自哪里：不是 GWAS Catalog 的谐化造成的

该文件是 GWAS Catalog 谐化过的 .h.tsv.gz，带 hm_code 列。官方 hm_code 表的相关两行：

代码	官方含义（原文）
10	"Forward strand; Alleles correct"
11	"Forward strand; Flipped alleles"

实测本课题 1,278 个 T1D 工具位点的 hm_code：

hm_code	位点数	★反转	无法判定	正确
10	1,277	1,197	59	21
11	1	1	0	0

几乎全部是 10 —— 即谐化流程判定「正向链、等位无需改动」，没有做过任何翻转。 既然没翻转过，反转就不可能是谐化引入的，只能来自作者提交的原始文件。

推论（对使用同一数据集的人有用）：GWAS Catalog 的谐化只核对等位身份是否与参考 VCF 一致， 不校验 effect_allele_frequency 是否真的对应 effect allele。作者侧的频率/等位错配会原样通过，并被盖上「Alleles correct」的戳。

8.4 造成了什么：一个位点看得最清楚

rs2273804（WARS 的 T1D 工具，回文 G/C）：

来源	效应等位	该等位频率	beta
1000G EUR	G	0.2247	—
暴露 UKB-PPP	G	0.2599 ✅	−0.258
T1D 源文件	G	0.7402 ❌	−0.0629
我们谐化后	G	0.2598	+0.0629 ← 符号被翻了

原始文件与我们存的效应等位都是 G，却翻了符号——同一个等位不该翻。 rs13406632、rs3815079、rs4823082 是同一模式。

8.5 影响面（已逐项实测）

不受影响

四个并发症端点（糖网 / 黄斑 / 糖肾 / 神经）：EA 零错位、Wald 比值差 < 10⁻¹⁴
2 型糖尿病（Mahajan）：eaf 相关 0.9991
T1D 的显著蛋白集合、p 值、FDR——符号翻转不改变 |β|、se、p，57 个仍是 57 个
T1D 的 1,502 个非回文工具：靠字母对齐，本来就对
三个候选 IFNAR1 / ERMAP / APOL1
共定位、跨平台、赢者诅咒、反向 MR、LDSC（rg 不使用 eaf）
糖尿病对照臂的**「并发症增强 11 / 并发症特异 10」**——这两类的判据用绝对值或 p 值，与方向无关

这里我预判错了一条，记录下来

本节原写「糖尿病驱动 30 也不受影响」。实测后是 35。

原因：判据本身确实与方向无关，但分类的输入变了——那 5 条原判「方向背离」的关联翻正后符号变成同号，于是从「方向背离」桶落进了「糖尿病驱动」桶。我只检查了判据公式，没考虑上游类别之间会发生迁移。

教训：判断「某个下游数字是否受影响」时，不能只看它的计算公式是否用到被改的量，还要看它的输入集合是否会因为别处的重分类而变化。

受影响

T1D 的 220 个回文工具（12.8%）方向反了
★「方向背离」6 条里有 5 条被重判（2026-08-05 实测已验证）：TRIM40（黄斑、糖网）、SIGLEC5（黄斑、糖网）、 WARS（糖网）；GALNT3（糖网）那条不受影响
★37% / 63% 这个主句数字可能微调——那 5 条翻正方向后会被重新归类：落进「并发症增强」会从 63% 挪到 37%，落进「糖尿病驱动」则不变。须跑完才知道。

8.6 为什么以前没查出来

以前的审计只查「数字有没有抄错、结论有没有过时」，没有任何一条检查输入数据本身是否自洽。

而方向错误不改变 p 值、置信区间宽度和图形，从统计量上完全看不出异常—— 只能靠外部参照（1000G）逐位点比对才能发现。这个洞从换源日 2026-07-29 存在至今。

8.7 修复：做了什么、结果如何

五步全部执行完毕（2026-08-05）：

#	做了什么	落点
1	manifest 增列 `eaf_orientation`（`effect`/`other`），不硬编码翻转	`outcome_manifest.csv`、`R/adapters/read_tsv.R`
2	启动期断言：每个结局读入后抽非回文位点对 1000G 比频率，反转率 >50% 直接 `stop()`	新增 `R/adapters/assert_eaf.R`；开关 `config.yaml` 的 `qc.assert_eaf`
3	重跑 R9_dm 主筛查＋糖尿病对照臂＋靶点表＋出图；同步 AMD 仓库并重跑 R13_dm	—
4	同步 PPT、中英文结果稿、docs 的数字	见下表
5	两个核查脚本纳入常规审计，并新增哨兵位点硬断言	`tools/verify_input_sanity.R` D 项

第 2 步之外还必须做的一件事

断点缓存指纹从 v1 升到 v2 并纳入 eaf_orientation。若不改指纹，改完 manifest 后缓存仍然命中，修复会静默不生效—— 这是「断点续跑掩盖上游修改」的经典坑，本项目此前已在别处踩过一次。

修复前后

指标	修复前	修复后
T1D `cor(eaf, 1000G EUR)`	−0.9922	+0.9922
T1D 反转率	93.7%	1.5%
`verify_harmonise.R`	27 PASS / 1 FAIL	28 PASS / 0 FAIL
各结局显著数	24/19/0/5/9/57/28	完全不变
糖尿病驱动	30（53%）	35（61%）
方向背离	6	1（仅剩 GALNT3×糖网）
散点异号	13	8
37% / 63% 主句	—	不变（5 条在同一桶内移动）
三个 A 级候选分层	—	不变

翻正的 5 条：TRIM40（黄斑、糖网）、SIGLEC5（黄斑、糖网）、WARS（糖网）——与修复前的预测完全一致。

给合稿与投稿的提示

Methods 里应如实写明：我们发现该数据集的等位频率列方向异常并做了校正。这对后续使用同一数据集的人是有价值的信息。

对应 config 字段

yaml

iv:
  harmonise_action: 2        # 1=假设正向 2=按MAF推断 3=丢弃回文
  steiger_filtering: true    # Steiger 定向过滤

复现本页的实测

bash

cd /home/research/mr-pipeline-r9
Rscript tools/verify_harmonise.R

该脚本会对全部结局逐个检查上述三项，任何一项不通过即以非零码退出。

判定「算不算问题」的分级标准、阈值出处分层、以及按外部 skill 清单给本课题的打分，见判定标准与 skill 合规打分。

下一步：谐化后的工具变量进入 MR 估计方法。IV 数（1 / 2 / ≥3）决定用哪套方法。

03 · 谐化与定向过滤 ​

零、先讲清楚基础：什么是「效应等位」，为什么必须对齐 ​

1. 一个位点上有两种碱基 ​

2. 「效应等位」只是一个记账的参照物 ​

3. 换一个效应等位，β 就要变号 ​

4. 不对齐会发生什么：方向彻底反过来 ​

5. 「谐化」（harmonisation）就是干这件事 ​

6. 翻号时的三条纪律 ​

7. 什么情况下翻号是错的 ​

一、还有一类更麻烦的：回文 SNP ​

什么是回文 ​

怎么解决：用等位频率 ​

但频率接近 0.5 时就没救了 ​

一·五、文献依据：以上做法出自哪里 ​

最权威的一条：Burgess 指南第 4 节 ​

其余三条 ​

本流水线每一步的依据对照 ​

二、转成 TwoSampleMR 格式 ​

三、谐化 + 回文处理（harmonise_pair） ​

四、Steiger 定向过滤（apply_steiger） ​

五、单对编排（build_harmonised） ​

六、本课题的实测核查结果 ​

6.1 效应等位是否逐条对齐 ​

6.2 等位频率的独立佐证 ​

6.3 回文 SNP：最危险的那一类是零 ​

6.4 全数据源列映射审计（2026-08-05 新增） ​

6.5 全结局 × 1000G 第三方裁决（2026-08-05 新增） ​

6.6 本轮未覆盖的部分（如实声明，不要当成已查） ​

七、与导师实现的差别：这是「匹配」问题，不是「谐化」问题 ​

7.1 涉及的三个原始文件（可自行点开核对） ​

7.2 三行代码逐行拆解 ​

L86 —— 把暴露侧的变异 ID 削掉尾巴 ​

L88 —— 给 Mahajan 也拼一个键 ​

L96 —— 用键做筛选 ​

三行连起来的隐含要求 ​

7.3 具体例子甲：APOE rs429358（被漏掉） ​

怎么知道效应等位是 T？——四层证据，只信第一层是不够的 ​

7.4 对照例子乙：MANSC4 rs34311349（能命中） ​

7.5 实测统计 ​

7.6 同一个例子，讲清「效应等位必须一致」 ​

7.7 原文的注释实况 ​

7.8 要说公道话的地方 ​

7.9 那我们是怎么做的 ​

路径一 · 按 rsID 匹配（默认） ​

路径二 · position 模式 ＋ 排序等位键 ​

效果实测 ​

为什么方向不放进键里，而交给谐化 ​

7.10 最终判定：他到底做错了没有 ​

第一层 · 概念理解：没错 ​

第二层 · 效应量方向：没错，且已实证 ​

第三层 · 真正做错的是什么：拿「主观约定」当成了「客观身份」 ​

一个变异身上有两种完全不同的信息 ​

他的错：把两列揉进了一个键 ​

为什么必须拆成两步 ​

修法：一行的事 ​

为什么这不只是「少拿了点数据」 ​

第四层 · 后果有多严重：比「丢 42%」更具体 ​

结论 ​

但共定位那一步没有这层保护 ​

八、已发现的真实故障：T1D 源的等位频率列整列反转 ​

8.1 官方文档怎么说（这决定了失效模式） ​

8.2 实测：整列反转，93.7% ​

对照组：同一方法验证暴露侧（2026-08-05 补做） ​

8.3 反转来自哪里：不是 GWAS Catalog 的谐化造成的 ​

8.4 造成了什么：一个位点看得最清楚 ​

8.5 影响面（已逐项实测） ​

8.6 为什么以前没查出来 ​

8.7 修复：做了什么、结果如何 ​

对应 config 字段 ​

复现本页的实测 ​

03 · 谐化与定向过滤

零、先讲清楚基础：什么是「效应等位」，为什么必须对齐

1. 一个位点上有两种碱基

2. 「效应等位」只是一个记账的参照物

3. 换一个效应等位，β 就要变号

4. 不对齐会发生什么：方向彻底反过来

5. 「谐化」（harmonisation）就是干这件事

6. 翻号时的三条纪律

7. 什么情况下翻号是错的

一、还有一类更麻烦的：回文 SNP

什么是回文

怎么解决：用等位频率

但频率接近 0.5 时就没救了

一·五、文献依据：以上做法出自哪里

最权威的一条：Burgess 指南第 4 节

其余三条

本流水线每一步的依据对照

二、转成 TwoSampleMR 格式

三、谐化 + 回文处理（`harmonise_pair`）

四、Steiger 定向过滤（`apply_steiger`）

五、单对编排（`build_harmonised`）

六、本课题的实测核查结果

6.1 效应等位是否逐条对齐

6.2 等位频率的独立佐证

6.3 回文 SNP：最危险的那一类是零

6.4 全数据源列映射审计（2026-08-05 新增）

6.5 全结局 × 1000G 第三方裁决（2026-08-05 新增）

6.6 本轮未覆盖的部分（如实声明，不要当成已查）

七、与导师实现的差别：这是「匹配」问题，不是「谐化」问题

7.1 涉及的三个原始文件（可自行点开核对）

7.2 三行代码逐行拆解

L86 —— 把暴露侧的变异 ID 削掉尾巴

L88 —— 给 Mahajan 也拼一个键

L96 —— 用键做筛选

三行连起来的隐含要求

7.3 具体例子甲：APOE rs429358（被漏掉）

怎么知道效应等位是 T？——四层证据，只信第一层是不够的

7.4 对照例子乙：MANSC4 rs34311349（能命中）

7.5 实测统计

7.6 同一个例子，讲清「效应等位必须一致」

7.7 原文的注释实况

7.8 要说公道话的地方

7.9 那我们是怎么做的

路径一 · 按 rsID 匹配（默认）

路径二 · position 模式＋排序等位键

效果实测

为什么方向不放进键里，而交给谐化

7.10 最终判定：他到底做错了没有

第一层 · 概念理解：没错

第二层 · 效应量方向：没错，且已实证

第三层 · 真正做错的是什么：拿「主观约定」当成了「客观身份」

一个变异身上有两种完全不同的信息

他的错：把两列揉进了一个键

为什么必须拆成两步

修法：一行的事

为什么这不只是「少拿了点数据」

第四层 · 后果有多严重：比「丢 42%」更具体

结论

但共定位那一步没有这层保护

八、已发现的真实故障：T1D 源的等位频率列整列反转

8.1 官方文档怎么说（这决定了失效模式）

8.2 实测：整列反转，93.7%

对照组：同一方法验证暴露侧（2026-08-05 补做）

8.3 反转来自哪里：不是 GWAS Catalog 的谐化造成的

8.4 造成了什么：一个位点看得最清楚

8.5 影响面（已逐项实测）

8.6 为什么以前没查出来

8.7 修复：做了什么、结果如何

对应 config 字段

复现本页的实测