我们的实现与导师代码的逐条对比

对象：F:\project\MR code\ 下导师的全部代码，与 mr-pipeline-r9 的实现。
一句话结论：MR 主体两套完全没有差异——我们的结果与他自己跑出的结果吻合到 1e-15；真正的差异只有两处（T1D 数据源、T2D 匹配键）。但下游三步（共定位计数、UpSet、PheWAS 方向）各有一处会改变结论的实现问题，「36 对 / 20 蛋白」与我们「19 对 / 13 蛋白」的落差已在代码层定位到根因。

〇、阅读覆盖范围（如实声明）

文件	行数	本次
`Whole story line.docx`（mrdata）/ `_v2.docx`（MR code）	各 88 段	逐字读完；两版文本完全相同
`1_r9_mr_diabetes.Rmd`	549	逐字读完
`2_diabete_validation.Rmd`	278	逐字读完
`4_r9_hyprcoloc_newdiabete.Rmd`	2067	逐字读完全部 2067 行（16 个逐组合重复块已逐块核对）
`significant_results_retinopathy.csv`	25	逐字读完
`hyprcoloc.R` vs `hyprcoloc_modify.R`	各 1582	完整 diff 读完（64 行差异）
`3_r9_upsetR.Rmd`	210	逐字读完
`5_r9_data_clean.Rmd`	340	逐字读完
`6_r9_phewas.R`	396	逐字读完
`7_network.Rmd`	514	逐字读完
`8_sensitivity_AMD.Rmd`	247	逐字读完
`R9_manifest.tsv`	2,4xx 行	按端点检索核对（非通读）

F:\project\MR code\ 下全部 8 个分析脚本已 100% 逐行读完。 未读的只剩 .Rhistory（交互历史，非流水线）与 result.pptx。下文凡标注行号者均为实读；未读部分不作论断。

一、★ 端点：五个并发症完全相同

此前根据 storyline 文档里的例数，我推断过"导师用了 FinnGen R10–R12、T2D 用的是 T2D_WIDE"。 读代码后确认这个推断是错的，已撤回。 他实际读的文件与写死的例数：

结局	他实际读的文件	他写的例数	我们	差异
糖网	`finngen_R9_DM_RETINOPATHY_EXMORE.gz`（F1 L72）	10,413（F1 L88）	10,413	无
黄斑病变	`finngen_R9_DM_MACULOPATHY_EXMORE.gz`（F1 L182）	3,572（F1 L198）	3,572	无
新生血管性青光眼	`finngen_R9_DM_NEOVASCULAR_GLAUCOMA.gz`（F1 L267）	1,100（F1 L282）	1,100	无
糖肾	`finngen_R9_DM_NEPHROPATHY_EXMORE.gz`（F1 L352）	4,111（F1 L379）	4,111	无
糖尿病神经病变	`finngen_R9_DM_NEUROPATHY.gz`（F1 L435）	2,843（F1 L461）	2,843	无
2 型糖尿病	`Mahajan.NatGenet2018b.T2D-noUKBB.European.txt`（F2 L81）	—	同一文件	无
1 型糖尿病	`GCST90475661.tsv.gz`（F2 L186）	—	`GCST90824163`	★有

文档里的例数为什么对不上

storyline 里的 12,681 / 1,353 / 4,526 / 49,101 都不是代码跑出来的：

F1 L87、L101 注释链接 https://r7.risteys.finngen.fi/phenocode/DM_RETINOPATHY，并写着"Retinopathy(12661)"——那是 R7 网页上另一个端点（非 EXMORE） 的数字，代码实际用的是 R9 EXMORE 的 10,413。
F2 L73–77 一整段注释：T1D (4526): r7.risteys.../T1D、 T2D wide(49101): r7.risteys.../T2D_WIDE——是早期 FinnGen 方案的遗留注释，实际代码在 L81 换成了 Mahajan。

所以：文档数字取自 R7 网页与遗留注释，代码跑的是 R9 文件加外部 T1D/T2D。这是文档与代码不一致，不是版本差异。 合稿时必须以代码为准重写 Table 1。

二、★★ 交叉验证：我们的糖网结果与他的逐条完全一致

他目录里留有自己跑出的 significant_results_retinopathy.csv（24 条）。逐条比对：

核对项	结果
显著蛋白个数	24 vs 24
蛋白集合	完全相同，无一进出
工具 SNP	24/24 为同一个 SNP
MR 效应量 beta	最大差 9.99×10⁻¹⁶
标准误 \|se\|	最大差 9.99×10⁻¹⁶

两套完全独立编写的代码，在相同输入上吻合到浮点精度。 这是我们并发症结果准确性最强的一条证据，也说明两边在 harmonise、Wald 比值、FDR 这三步上的实现是等价的。

三、真正的两处差异

差异一：T1D 数据源（★ 影响他的 T1D 全部结论）

他用 GCST90475661。我们在 2026-07-27 查明该数据集是 MVP 的 phecode「T1D」，实为 2 型糖尿病主导的表型：

TCF7L2（经典 T2D 位点）p = 9.5×10⁻⁶⁶
MHC 区（T1D 的主效区）仅 1.28 倍富集
与 T2D 的遗传相关 rg = 0.886

我们换成 GCST90824163（Nat Genet 2026，20,355 例 / 797,363 对照），换源后 T1D↔T2D 的 rg 从 0.886 降到 0.087，符合两者病因不同的常识。

含义：他那一支"T1D"的结果实际建立在一个标签错误的表型上，文中"T1D 与 T2D 共享三个蛋白"这类结论需要重做。详见外部 T1D 数据源与样本重叠。

差异二：T2D 的变异匹配键（★ 实测丢 42.3% 工具）

他不按 rsID 匹配，而是拼字符串（F2 L86–88、L96）：

ready_exposure_cis$clean_variant <- sub(":imp:v1.*", "", Variant.ID)   # CHROM:GENPOS:A0:A1
temp_t2d$clean_variant <- paste0(Chr, ":", Pos, ":", NEA, ":", EA)     # chr:pos:NEA:EA
selected_outcomes_diabetetp2 <- subset(temp_t2d, clean_variant %in% ready_exposure_cis$clean_variant)

这要求 UKB-PPP 的 A0:A1 顺序与 Mahajan 的 NEA:EA 顺序完全一致，顺序相反的变异直接匹配不上、静默丢失（不报错、不计数）。

实测（2026-07-31，用 1,955 个 cis 变异对 Mahajan 全库 21,508,698 个变异）：

匹配方式	命中
按他的键 `chr:pos:NEA:EA`	956
按相反顺序 `chr:pos:EA:NEA`	700
合计（真实可匹配数）	1,656

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

四、他代码里三个会影响结果的缺陷

缺陷 1　`se_mr` 没取绝对值 → 置信区间上下限颠倒

F1 L117（以及 L223、L307、L391、L473，F2 L141、L236，共七处同样写法）：

se_mr <- rows2$se.outcome[j] / rows1$beta.exposure[i]     # ← 分母是带符号的

正确写法是除以 abs(beta.exposure)。后果：beta.exposure < 0 时 se 为负，于是 lower_CI = ratio - 1.96*se 反而大于 upper_CI = ratio + 1.96*se。

在他自己那份糖网结果里，24 条有 13 条的 CI 上下限是颠倒的，例如：

蛋白	beta.exposure	se_mr	lower_CI	upper_CI
AGER	−0.170	−0.1021	−1.7548	−2.1548
AIF1	−0.111	−0.1310	−0.2180	−0.7315
APOE	−1.012	−0.0189	0.1672	0.0931

p 值不受影响（beta.exposure 在 z 里约掉了）。我们取了绝对值，CI 是对的。

缺陷 2　共定位与 MR 用了不同的糖尿病数据集

F4 L225–226（共定位步骤读入的结局）：

data_diabetetp2 <- read_delim('.../summary_stats_finngen_R9_T2D.gz')
data_diabetetp1 <- read_delim('.../summary_stats_finngen_R9_T1D.gz')

而 MR 步骤（F2）用的是 Mahajan 与 GCST90475661。 同一篇文章里，提名蛋白用的是一套数据，验证共定位用的是另一套。 共定位本应验证"MR 那个信号是不是同一个因果变异"，换了数据集就失去了验证意义。

缺陷 3　共定位的 SNP 匹配同样是等位顺序敏感的字符串

F4 L87–88 与 L110：

# 蛋白侧
mutate(newid = paste(CHROM, GENPOS, ALLELE0, ALLELE1, sep = "_"))
# 疾病侧
mutate(newid = paste(chrom, pos, ref, alt, sep = "_"))

与差异二同一类问题：顺序不一致就丢，且 L114 用 intersect 逐个疾病求交集，要求该 SNP 在列表里每一个疾病数据集中都存在，集合会被逐步压缩。

另外 L163 与 L180 用硬编码列号取 beta/se：

colnames(protein_snp)[c(15, 10, 11)] <- c("SNP", "beta.outcome", "se.outcome")
colnames(disease_snp)[c(14, 9, 10)]  <- c("SNP", "beta.outcome", "se.outcome")

列顺序一变就会把错误的列当成 beta/se，且不会报错。

五、★★ 共定位「36 对 / 20 蛋白」的根因（本次读完 F3+F4+F5 后定位）

我们之前只知道「他报的 20 个蛋白里有 11 个从未进过簇」，不知道为什么。读完 3_r9_upsetR.Rmd（210 行）与 4_r9_hyprcoloc_newdiabete.Rmd（2067 行）后，根因有三层，每一层都独立地把数字放大。

根因 1　UpSet 统计的是「被测的组合」，不是「聚到一起的性状」

3_r9_upsetR.Rmd L175–182：

bin <- df_sig %>%
  mutate(tmp = str_split(disease_combination, "\\+")) %>%   # ★ disease_combination
  tidyr::unnest_longer(tmp, values_to = "dis") %>%
  mutate(flag = 1) %>%
  select(protein, disease_combination, dis, flag, posterior_prob) %>%
  filter(dis %in% basic) %>%
  distinct() %>%
  tidyr::pivot_wider(names_from = dis, values_from = flag, values_fill = 0)

HyPrColoc 的输出有两列：

disease_combination（F4 L432 由 cbind 写入）＝这一次喂给 hyprcoloc 的疾病清单
traits＝hyprcoloc 判定真正聚到一起的性状

拆的是前者。只要一个疾病被放进这次分析，它就在 UpSet 里被记一票，不管 hyprcoloc 有没有把它和蛋白聚进同一个簇。

根因 2　同一段代码里 `posterior_prob` 从头到尾没有被用来过滤

同一文件 L162–167：

thr <- 0.9                  # 
max_facets <- 20            # 
...
df_sig <- hypco_sum         # ★ 直接赋值，没有任何 filter

thr 定义之后再未出现（全文件检索 thr 只有 L162 这一处）。 posterior_prob 只在 L179 被 select 保留，没有出现在任何 filter 里。

→ 按这段代码的字面执行，UpSet 与「20」这个数计的是参与过共定位的蛋白数，不是共定位成功的蛋白数。存放结果的目录名恰好就叫 hypcoloc 20。

这一条的边界

sum_plot.csv 本身没有随代码提供。如果它是在这段脚本之外先筛过 PP 再存的，那么「20」仍可能是筛后的数。我只能确认这段脚本没筛，不能确认上游有没有筛。这一点需要直接问导师，不要替他下结论。

根因 3　子组合穷举导致同一对被重复计数

F4 的设计是：先按「这个 SNP 在哪几个疾病里显著」把 SNP 分成互斥的组，再对每组跑 hyprcoloc。但组内不是只跑一次，而是穷举所有子集：

disease_combinations <- unlist(
  lapply(1:k, function(n) combn(names(disease_datasets), n, simplify = FALSE)),
  recursive = FALSE)                      # → 2^k − 1 种组合

k 是该组的疾病数（F4 L375 / L483 / L808 / L1342 / L1668 / L1993 分别是 1/1/2/3/4/5）。按代码里 16 个块的分组规模统计：

分组疾病数 k	组合数 2^k−1	块数	该层 hyprcoloc 运行次数
1	1	4 块（dia2 / dia1 / mac / ret）	83
2	3	5 块	24
3	7	3 块	28
4	15	3 块	45
5	31	1 块	31
合计		16 块	约 211 次

一个落在五病组的蛋白，会被跑 31 次；其中含黄斑的组合有 16 种。同一个「蛋白–黄斑」关系因此最多可以在 16 行里各出现一次。 按行计数 = 按「测了多少次」计数，不是按「发现了多少个关系」计数。

根因 4　多疾病组合里 hyprcoloc 可以返回「不含蛋白」的簇

F4 L431–432：

protein_result <- data.frame(hyprcoloc_results$results)
protein_result <- cbind(protein_result, protein = protein_name, ...)   # ★ 无条件贴上蛋白名

当喂进去的是 {蛋白, 糖网, 黄斑} 时，hyprcoloc 完全可能返回一个 traits = "retinopathy, maculopathy" 的簇——糖网与黄斑之间的共定位（两者 rg=0.946，这个簇是真实且预期之中的）。但 L432 无条件把 protein=XXX 贴上去。 只看 protein 列 + posterior_prob，就会把一个疾病–疾病共定位读成蛋白–疾病共定位。

这正是我们此前查出的「20 个里 11 个从未进过簇」的产生机制——现在有代码依据了，不再是推断。

我们的口径

	导师（按代码字面）	我们
计数单位	hyprcoloc 输出行	去重后的 (蛋白, 疾病) 对
疾病来源	`disease_combination`（被测的）	`traits`（聚进簇的）
PP 过滤	这段脚本里没有	PP > 0.7
要求蛋白在簇内	否	是
结果	36 对 / 20 蛋白	19 对 / 13 蛋白

两个数不是同一个东西，不存在谁对谁错的比较；但论文里只能报后者，前者的定义无法支撑「蛋白与疾病共定位」这句话。我们的两张 UpSet（宽松/严格）与三种口径说明见共定位证据分档与 SuSiE 降级第四节。

附带查出：新生血管性青光眼被整条排除在共定位之外

3_r9_upsetR.Rmd L26 读入了 significant_results_neovascular_glaucoma.csv， 但 L41–46 只给 6 个结局加了二值列，L50–54 只合并了这 6 个，青光眼两处都不在其中：

significant_results_retinopathy$retinopathy <- 1
significant_results_neuropathy$neuropathy   <- 1
significant_results_nephropathy$nephropathy <- 1
significant_results_maculopathy$maculopathy <- 1
significant_results_diabetetp2$diabete2     <- 1
significant_results_diabetetp1$diabete1     <- 1
#   ← 没有 neovascular_glaucoma

后续所有 sameSNP_* 分组都从这 6 个集合派生，所以青光眼有 MR 结果、但共定位覆盖为 0，UpSet 图里也没有这条集合。读入了却没用，R 不会报错。

附带查出：`5_r9_data_clean.Rmd` 是上一轮（非 R9）的遗留脚本

	F4 输出	F5 读写
路径	`taoci/1/r9/hyprcoloc/result/cleaned_results/`	`taoci/1/hyprcoloc/hyprcoloc_Result/cleaned_results/`（无 r9）
对象名	`combined_protein_results_diabete2`、`dia1_mac_ret`…	`combined_protein_results_1_maculopathy_diabete1_nephropathy`、`combined_protein_results_18_diabete1`

两边的路径和对象名都对不上，F5 也从不自己读入这些对象（假定已在环境里）。文件名虽叫 5_r9_，内容指向的是 R9 之前那一轮。

而且 F5 的解析是按行号手写规则把 hyprcoloc 的文本结果切成列，例如 L304–305：

onetrait <- c(1,3:7,9,11,13:17)
twotrait <- c(2,8,10,12)

process_row 对不在这两个名单里的行号不返回任何值（隐式 NULL），下游 c(x, rep(NA, max_len - length(x))) 会把它变成整行 NA，不报错。对象名是 ..._18_diabete1（暗示 18 行），而规则只覆盖到第 17 行； retinopathy 块（L268–275）只覆盖 row_index <= 3； nephropathy_diabete2 块（L228–238）只覆盖 1–4 行。 这类越界在结果里表现为悄悄多出的空行，而不是报错。

六、★ PheWAS 的方向标注：标的是「疾病方向」，不是「性状方向」

6_r9_phewas.R L27–41 是全部方向逻辑：

if (rs_data$EA[i] == temp$effect_allele.outcome[j]) {
  combined_row <- cbind(rs_data[i, ], temp[j, ])
  combined_row$MatchType <- ifelse(combined_row$beta.outcome >= 0, 'positive', 'negative')
}
if (rs_data$NEA[i] == temp$effect_allele.outcome[j]) {
  combined_row <- cbind(rs_data[i, ], temp[j, ])
  combined_row$MatchType <- ifelse(combined_row$beta.outcome >= 0, 'negative', 'positive')
}

rs_data ＝该 SNP 的 PheWAS 结果（含 EA/NEA/性状 beta）
temp ＝ sig_Data（significant_mr_all.csv）里该 SNP 的 MR/疾病记录
combined_row$beta.outcome 来自 temp，即疾病的 beta

等位对齐本身是对的（EA 对上就保号，对上 NEA 就翻号）。 但被对齐、被赋号的自始至终是疾病效应量；PheWAS 文件自己的 beta 一次都没被读。

后果：同一个 SNP–疾病对下的所有 PheWAS 性状，MatchType 全都是同一个符号。 而 L119/L126 把它当成性状方向画图：

Correlation = combined_results_signif$MatchType,
signed_logP = logP * ifelse(Correlation == "positive", 1, -1)

图标题写的是 Directional PheWAS Plot (Signed -log10 p-value)。 纵轴是 PheWAS 性状的 p 值，符号却是疾病的方向——两者拼在一起没有可解释的含义，也无法回答药靶问题（「压低这个蛋白，这个性状会往哪边走」）。

这个标注一路传到 7_network.Rmd：L370/L410 的 Correlation = MatchType、 L293–296 按 disease_sign 切 positive_targets / negative_targets，以及 L319–332 用它删基因，全部继承同一个符号。

我们 Task 1 做的正是这件事的正确版本： sign(wald_trait) == sign(b_index) —— 两侧都用各自的 Wald 比值，同号判为潜在获益、异号判为潜在风险。见方向感知药靶 PheWAS。

同一文件的另外三处

行	问题	后果
L20–25	内层 `temp <- sig_data %>% filter(SNP == snp_id)` 与 `i` 无关，却放在 `i` 循环里；一个 SNP 有 m 条 MR 记录、n 条 PheWAS 关联时产生 m×n 行	行数膨胀
L81	`bonferroni_threshold <- 0.05 / nrow(combined_results)` 用的正是上面膨胀后的行数	阈值偏保守（不是偏松），会漏而不会多
L229 / L296	两张图的阈值线写死 `yintercept = 6`，与 L81 算出的 Bonferroni 无关	图上的线不是所声称的阈值
L46	`result <- result[, -c(13:18)]` 写死列号	列序一变即静默错位

覆盖面上：他的 PheWAS 目录里是 20 个 SNP 的 csv（L54 注释 "all 20 snps"）；我们是 92 个 SNP、15,789 条关联（OpenGWAS，50,164 个数据集）。

七、★ AMD 敏感性分析：代码读的是一般人群 AMD，标签写的是「diabete AMD」

8_sensitivity_AMD.Rmd：

行	内容
L71 注释	`https://r9.risteys.finngen.fi/phenocode/DM_AMD`
L72 代码	`read_delim('.../summary_stats_finngen_R9_H7_AMD.gz')`
L85	`selected_outcomes_AMD$outcome <- "diabete AMD"`
L87 注释	例数来源写的是 `r7.risteys...`（又一处 R7 网页）
L88	`samplesize.outcome = 8913`

用他自己随包提供的 R9_manifest.tsv 检索：

H7_AMD    8913 例 / 348936 对照    Age-related macular degeneration (whether dry or wet)
DM_AMD    ← 在 R9 manifest 中不存在

8913 这个数字是对的（H7_AMD 的真实例数）；
但这个端点是一般人群的年龄相关性黄斑变性，不是糖尿病相关；
注释指向的 DM_AMD 在 R9 里根本不存在；
结局却被命名为 "diabete AMD"。

→ 这一节的结论只能表述为「H7_AMD（一般人群 AMD）」，凡按「糖尿病性 AMD」解读的句子都要改。这与 Table 1 例数问题是同一类（注释/文档来自 R7 网页，代码用的是 R9 文件）。

同文件另两处：L113 se_mr <- rows2$se.outcome[j] / rows1$beta.exposure[i]同样没取绝对值（CI 翻转问题在 AMD 这一节重现）； L43 ready_exposure <- exposure_combined[, c(8,20,9,11,13,14,36,35,12,22,24,25,15)] 仍是写死列号。

LDSC 部分（L212–233）本身没问题：res 与 res2 只差 sample.prev，而 rg 对 observed/liability 尺度不变，两者应给出同一个 rg，属冗余而非错误；只是 population.prev 全为 NA，h² 停留在观测尺度、跨性状不可直接比。

八、其他差异（不一定是错，但两套不同）

项	导师	我们
工具强度过滤	算了 F 值（F1 L93）但从未用于过滤	F ≥ 10（实测最小 43.8）
LD clump	无	plink 本地 clump
Wald 计算	手工双重循环（F1 L107–154），`rows1` 取自未协调的暴露表	`TwoSampleMR::mr()`
重复 SNP 处理	双重循环产生交叉配对（他自己 L175 注释说有 22 个重复 SNP）	每蛋白单哨兵 SNP
缺失值	`na.omit(整个结局表)`——任一列有 NA 就整行丢	按需列判断
共定位窗口	±500 kb（F4 L63–64）	见流水线配置
共定位阈值	PP ≥ 0.7，单方法	PP ≥ 0.8 且两法一致，另设中间档
HyPrColoc 源码	`hyprcoloc_modify.R` 把 32 处 `Rmpfr::asNumeric` 改成基础 `as.numeric`，且 source 顺序让改版覆盖原版	用官方包
PheWAS 库	GWAS Atlas	OpenGWAS，50,164 个数据集
糖尿病易感性对照	无	有
反向 MR / LDSC / 外部复制 / MHC 分离 / 共定位分档	无	有

asNumeric → as.numeric 这一改动去掉了 HyPrColoc 对近似贝叶斯因子求和的高精度保护（ABF 可以大到超出 double 范围）。这一条是风险提示，未实测其数值后果。

九、一个双方共有的性质，写作时要说明

单 SNP Wald 比值用一阶标准误时：

{se}_{M R} = \frac{{se}_{o u t c o m e}}{| β_{e x p o s u r e} |}, z = \frac{β_{o u t c o m e} / β_{e x p o s u r e}}{{se}_{o u t c o m e} / | β_{e x p o s u r e} |} = \frac{β_{o u t c o m e}}{{se}_{o u t c o m e}}

MR 的 p 值恒等于结局 GWAS 在该 SNP 的 p 值，暴露侧完全约掉。

实测：他的 p 与 2·pnorm(−|β_out/se_out|) 相对差 1.96×10⁻¹⁴；我们的 se 与 se_out/|β_exp| 差 5×10⁻¹⁶。两边都是如此，这是标准做法不是谁的错误，但意味着显著性完全由结局侧决定，暴露侧只提供方向与效应量尺度。审稿人可能会问，方法学部分应主动说明。

十、给合稿的行动清单

按「必须改 → 应当改 → 需澄清」排序。

必须改（不改会出现与数据不符的陈述）

Table 1 的例数按代码重写（10,413 / 3,572 / 1,100 / 4,111 / 2,843），删掉 R7 网页来源的 12,681 / 1,353 / 4,526 / 49,101。
T1D 换源重做——GCST90475661 不能用。
T2D 改按 rsID 匹配，找回 42.3% 的工具后重跑。
se_mr 加绝对值，13/24 的 CI 需要重出（MR 主表与 AMD 敏感性分析两处都要）。
共定位改用 traits 列计数、并要求蛋白在簇内、并设 PP 阈值—— 「36 对 / 20 蛋白」按现有脚本口径不能写成「共定位成功」。
AMD 一节改写为「H7_AMD（一般人群 AMD）」，删去「diabete AMD」的表述。
PheWAS 的方向标注重做——现在的 MatchType 是疾病方向，不能用作性状方向；positive_targets / negative_targets 的划分随之作废。

应当改

共定位与 MR 统一数据集（他共定位读 FinnGen R9 T1D/T2D，MR 读 Mahajan/GCST90475661）。
青光眼补进 UpSet 与共定位，否则五个并发症里有一个没有共定位证据却未说明。
PheWAS 阈值线画成实际算出的 Bonferroni，不要写死 6。
弃用 5_r9_data_clean.Rmd（路径与对象名指向 R9 之前那一轮，且按行号手写解析规则）。
写死列号（F4 L163/L180、F6 L46、F8 L43）全部改成按列名取。

需要问导师（不替他下结论）

sum_plot.csv 是否在 3_r9_upsetR.Rmd 之外先按 PP 筛过？ —— 决定「20」是「共定位成功的蛋白数」还是「参与共定位的蛋白数」。
storyline 正文几处内部不一致：同一句里黄斑出现 5 与 6 两个数；正文「45 个 protein–trait associations」与 Fig 1 标「45 proteins」（关联数 ≠ 蛋白数）；正文共定位「36 unique associations」与 Fig 2 标「20 proteins」；「Drug target」与「Sensitivity analysis」两节为空标题，Conclusion 亦为空。
hyprcoloc_modify.R 把 32 处 Rmpfr::asNumeric 换成基础 as.numeric 是有意为之（可能是为了绕开 Rmpfr 安装问题）还是临时改动？影响未实测。

十一、每条判断的官方依据（2026-08-05 补）

本页对导师代码的每一处判断，都能落到官方文档、标准规范或包源码上。以下逐条给出可核查的出处与原文，不用「一般认为」这类措辞。

11.1 效应等位（effect allele）—— 本页最核心的依据

定义

TwoSampleMR 官方 exposure vignette 原文：

effect_allele — "The allele of the SNP which has the effect marked in beta"

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（beta/OR/HR/z 四选一）
`effect_allele_frequency`	"Column 6: Frequency of the effect allele"	Mandatory

→ EA、beta、eaf 三者必须指向同一个等位，这是有强制性标准的，不是习惯。

「等位顺序相反」是常态，官方要求翻号处理而不是丢弃

Hemani et al., eLife 2018;7:e34408（MR-Base 平台论文， PMC5976434）， "Harmonising exposure and outcome SNP effects" 一节原文：

"A SNP with (for example) effect/non-effect alleles G/T for the exposure and T/G for the outcome are harmonised by flipping the sign of the SNP-outcome effect."

"To ensure that the effect sizes for the SNP reflect the same allele it is therefore necessary to switch the direction of the effect in either the exposure or outcome study."

→ 直接支撑本页差异二与缺陷 3。 导师的 chr:pos:NEA:EA 与 chr_pos_ref_alt 字符串键，把官方明确要求「翻号保留」的这批变异 静默丢弃了（实测 700/1,656 = 42.3%）。

回文 SNP 与链方向

同一篇 Hemani 2018 原文：

"SNPs with A/T or G/C alleles are known as palindromic SNPs, because their alleles are represented by the same pair of letters on the forward and reverse strands, which can introduce ambiguity into the identity of the effect allele... If reference strands are unknown, effect allele frequency can be used to resolve the ambiguity."

TwoSampleMR harmonise vignette 对三种 action 的原文定义：

`action`	官方原文
1	"Assume all alleles are presented on the forward strand"
2（默认，我们用的）	"Try to infer the forward strand alleles using allele frequency information"
3	"Correct the strand for non-palindromic SNPs, but drop all palindromic SNPs"

报告规范也把它列为必报项

STROBE-MR（Skrivankova et al. 2021， E&E 版 PMC8546498） Item 4c "Describe measurement, quality control and selection of genetic variants" 的解释部分要求报告：

"the presence or absence and handling of strand alignment, and orientation of effect and non-effect alleles"

术语表（Table 2）对 "Strand alignment" 的定义：

"ensures that the alleles in the exposure GWAS and the outcome GWAS are measured on the same DNA strand"

三条落地推论

effect_allele_frequency 必须是效应等位的频率——GWAS-SSF 强制字段定义。我们的外部 T1D 源违反了它自己 meta.yaml 里声明的 file_type: GWAS-SSF v1.0 （详见谐化与定向过滤 · 第八节）。
action = 2 定链的唯一依据就是频率 → 频率反 = 回文 SNP 方向反。
等位顺序相反是常态，官方要求翻号保留，不是匹配不上就丢。

11.2 `se_mr` 必须除以 `|β_exposure|` —— 有包源码为证

本页缺陷 1 的依据不是惯例，是官方实现。本机 TwoSampleMR（/home/research/R/x86_64-pc-linux-gnu-library/4.3/TwoSampleMR）源码：

mr_wald_ratio <- function(b_exp, b_out, se_exp, se_out, parameters) {
  ...
  b    <- b_out / b_exp
  se   <- se_out / abs(b_exp)        # ← 官方实现取绝对值
  pval <- stats::pnorm(abs(b)/se, lower.tail = FALSE) * 2
  return(list(b = b, se = se, pval = pval, nsnp = 1))
}

数学上同样只能如此：标准误是标准差的估计，恒非负； delta 法一阶近似给出 ${se}_{M R} = {se}_{Y} / | β_{X} |$ （Burgess, Small & Thompson, Stat Methods Med Res 2017;26:2333–2355）。

导师七处写成 rows2$se.outcome[j] / rows1$beta.exposure[i]（分母带符号）， β_exposure < 0 时 se 为负 → lower_CI > upper_CI。 他自己那份糖网结果 24 条里 13 条如此（见第四节表）。

补充：他的数据对象是 TwoSampleMR 格式（beta.exposure / se.outcome），也就是说标准实现就在同一个包里，只是没有调用。

11.3 复现这些依据

bash

# 官方 Wald 比值实现
Rscript -e 'print(TwoSampleMR::mr_wald_ratio)'

# GWAS-SSF 规范（Table 1 字段定义）
curl -sLO https://raw.githubusercontent.com/EBISPOT/gwas-summary-statistics-standard/master/gwas-ssf_v1.1.0.pdf
pdftotext -layout gwas-ssf_v1.1.0.pdf - | grep -A2 "effect allele frequency"

尚未取得官方出处的两条

缺陷 2（共定位与 MR 用了不同数据集）与根因 1（disease_combination vs traits）依据的是读代码 + 我们自己的 hyprcoloc 产物，不是引文。 hyprcoloc 官方文档（?hyprcoloc）只说返回 "clusters of colocalized traits"，未逐列说明两个字段的区别。这一条的措辞应保持在「按代码字面执行会……」，不要写成「官方规定应……」。
asNumeric → as.numeric 的数值后果仍未实测，维持风险提示，不作论断。

我们的实现与导师代码的逐条对比 ​

〇、阅读覆盖范围（如实声明） ​

一、★ 端点：五个并发症完全相同 ​

文档里的例数为什么对不上 ​

二、★★ 交叉验证：我们的糖网结果与他的逐条完全一致 ​

三、真正的两处差异 ​

差异一：T1D 数据源（★ 影响他的 T1D 全部结论） ​

差异二：T2D 的变异匹配键（★ 实测丢 42.3% 工具） ​

四、他代码里三个会影响结果的缺陷 ​

缺陷 1 se_mr 没取绝对值 → 置信区间上下限颠倒 ​

缺陷 2 共定位与 MR 用了不同的糖尿病数据集 ​

缺陷 3 共定位的 SNP 匹配同样是等位顺序敏感的字符串 ​

五、★★ 共定位「36 对 / 20 蛋白」的根因（本次读完 F3+F4+F5 后定位） ​

根因 1 UpSet 统计的是「被测的组合」，不是「聚到一起的性状」 ​

根因 2 同一段代码里 posterior_prob 从头到尾没有被用来过滤 ​

根因 3 子组合穷举导致同一对被重复计数 ​

根因 4 多疾病组合里 hyprcoloc 可以返回「不含蛋白」的簇 ​

我们的口径 ​

附带查出：新生血管性青光眼被整条排除在共定位之外 ​

附带查出：5_r9_data_clean.Rmd 是上一轮（非 R9）的遗留脚本 ​

六、★ PheWAS 的方向标注：标的是「疾病方向」，不是「性状方向」 ​

同一文件的另外三处 ​

七、★ AMD 敏感性分析：代码读的是一般人群 AMD，标签写的是「diabete AMD」 ​

八、其他差异（不一定是错，但两套不同） ​

九、一个双方共有的性质，写作时要说明 ​

十、给合稿的行动清单 ​

必须改（不改会出现与数据不符的陈述） ​

应当改 ​

需要问导师（不替他下结论） ​

十一、每条判断的官方依据（2026-08-05 补） ​

11.1 效应等位（effect allele）—— 本页最核心的依据 ​

定义 ​

「等位顺序相反」是常态，官方要求翻号处理而不是丢弃 ​

回文 SNP 与链方向 ​

报告规范也把它列为必报项 ​

11.2 se_mr 必须除以 |β_exposure| —— 有包源码为证 ​

11.3 复现这些依据 ​

相关页面 ​