主题
★ 07 · 代码审计与发表标准
回答两个问题："其他分析的代码有没有像孟德尔随机化那样逐一核对过？够不够发表高分 SCI？" 2026-07-24 做了三件事：① 逐脚本静态审查多组学 MR + 功能验证代码；② 对多组学 MR 头号命中做了独立数值复算（脱离流水线重算，证明算对了）；③ 核实了两处此前被我标错/标重的问题。准确性要求见 [[feedback_accuracy]]。
一、三种核对，现在各到什么程度
| 核对方式 | 主 MR 流水线 | 多组学 MR | 功能验证(dr-amd-targets) | 对接/MD |
|---|---|---|---|---|
| 静态代码审查 | ✅ | ✅ | ✅ 本次已做(7个核心脚本) | ❌(调成熟引擎,核参数) |
| 独立数值复算 | ✅ | ✅ 本次已做(见三) | ❌ | ❌ |
多组学 MR 现在和主 MR 同级：既过静态审查，也被独立复算证明。
二、多组学 MR 静态审查发现（含两处更正）
| 级别 | 位置 | 问题 | 状态 |
|---|---|---|---|
| 🔴→🟢 已修（2026-08-01） | 03d_parse_smr.R + 04_integrate.R L81 | ★HEIDI 判据把 p_HEIDI=NA 当成"通过"（is.na(p_HEIDI) | p_HEIDI>0.05）。算不出来＝无法排除连锁，与排除了连锁相反。这是本审计里唯一改动了已发布结论的缺陷：AGER 与 APOE 的甲基化层证据被撤下；全基因组稳健命中旧口径虚高 DR 43.7%、AMD 19.6–25.0% | 两处均改三分类 + 重跑 03d/04 + 重出矩阵与热图；"不可评估"的那批留档在 smr_heidi_not_evaluable.csv。CASP10/TNFRSF10A 三组学收敛不受影响 |
| 🟢 已修 | 04_integrate.R L81 | 热图 meQTL 格错标 "MR only"(应 "SMR+HEIDI") | 已改+重跑+换图 |
03c US_Blood | 出处已查明=Understanding Society 血 mQTL(见下) | 解决 | |
03c L27 | --diff-freq-prop 0.5 —— 我先前说"默认0.2、降低质控"是错的(见下) | 更正,降级 | |
| 🟠 应说明 | 01b L67 | eQTL FDR 把 AMD 4 相关结局 + 复制队列 IAMDGC 并进同一 BH 族 | 分族或论文注明 |
| 🟡 小限制 | extract_finngen | 逗号连多 rsID 行被精确匹配漏掉(多为罕见变异) | 方法里提一句 |
| 🟡 小限制 | run_coloc | quant 暴露未传 sdY,由 coloc 自估 | 标准做法,注明 |
| ⚪ 定位 | 04_integrate.R | 0/1/2 打分是启发式决策矩阵,非合并显著性 | 论文别当"综合显著" |
更正 1：US_Blood 出处已查明（不再是待核实）
.flist 内部路径写着 /mnt/data1/EPICQC/UnderstandingSociety/MatrixEQTL/...，探针是 EPIC 芯片 CpG(cg 开头,12.6 万探针)。US_Blood = Understanding Society(英国家庭纵向研究 UKHLS)全血 mQTL，SMR 官方格式，同目录还有配套的 Hannon et al. 血 mQTL 数据集(Aberdeen_Blood 等)。__MACOSX + ._ 资源叉文件证明是 Mac 压缩包解压 → 由用户下载而非本机生成。发表时按其官方下载页正式引用即可(数据身份已确定；主引用文献建议从下载出处核对)。
更正 2：SMR --diff-freq-prop 0.5 —— 我先前判重了
官方文档(已核)：--diff-freq 默认 0.2(逐 SNP：频率差>0.2 剔除)；--diff-freq-prop 默认 0.05(当被剔 SNP 占比超阈值就报错中止)。
- 我先前说"默认 0.2、放宽质控让坏 SNP 通过"——错。逐 SNP 的
--diff-freq 0.2没动、仍在剔坏 SNP(证据：三个*.snp_failed_freq_ck.list每次仍剔了约 1000 个 SNP)。 --diff-freq-prop 0.5只是把"因剔除过多而中止"的保险丝从 5% 抬到 50%，不改变用哪些 SNP。英国 UKHLS mQTL vs 芬兰/IAMDGC GWAS 频率天然有差，抬高它是让分析跑完的正常用法，不是偷工减料。发表时报告被剔 SNP 数即可。
三、独立数值复算：多组学 MR 头号命中已被证明
脚本 tools/verify_multiomics.R（不 source mr_funcs.R，重打 Z→β 公式、FinnGen 提取 awk、等位对齐、IVW 回归、独立重跑 coloc），对 CASP10 / TNFRSF10A × AMD 从原始 eQTLGen + FinnGen 重算，与流水线产出比：
| 量 | CASP10×AMD | TNFRSF10A×AMD |
|---|---|---|
| IVW β 相对差 | 2.2×10⁻¹⁵ | 2.9×10⁻¹⁵ |
| IVW se 相对差 | 2.1×10⁻¹⁵ | 0(逐位相同) |
| 工具数 | 13 = 13 ✓ | 23 = 23 ✓ |
| coloc PP.H4 | 0.9745 vs 0.9695 | 0.3548 vs 0.3511 |
IVW 估计吻合到浮点精度、工具数精确一致 → 多组学 MR 的 MR 主体计算被证明算对了，与主 MR 的 manual_verify 同级。coloc H4 差 ~0.005（独立提取的区域 SNP 数略不同：459/281 vs 480/288），结论方向完全一致（CASP10≈0.97 强共定位、TNFRSF10A≈0.35 弱）。
3.2 meQTL SMR 层：也已证明（tools/verify_layers.R）
独立重推 SMR 公式（b_SMR=b_GWAS/b_eQTL；T=z_E²z_G²/(z_E²+z_G²)；p=χ²₁；se=|b|/√T），并把 b_GWAS 对原始 FinnGen 核：
| 靶点(CpG top SNP) | b_SMR 相对差 | se 相对差 | p 相对差 | b_GWAS vs 原始 FinnGen |
|---|---|---|---|---|
| APOE (rs10414043) | 9.7×10⁻⁷ | 5.6×10⁻⁷ | 1.5×10⁻⁵ | 完全一致 |
| TNFRSF10A (rs13278062) | 2.2×10⁻⁶ | 1.1×10⁻⁶ | 9.2×10⁻⁶ | 完全一致(仅等位定向符号) |
| CASP10 (rs6435069) | 4.8×10⁻⁸ | 4.7×10⁻⁷ | 3.3×10⁻⁸ | 完全一致 |
残差 ~10⁻⁶ 仅因抄组分时取 6 位有效数字；SMR 公式与 GWAS 读取都证明对了。
3.3 代谢物层：复现到 ~1%（非逐位相同，如实）
独立重算 IVW（重打回归公式 + 独立 clump/等位对齐）：
| 代谢物×AMD | IVW b 相对差 | 工具数(独立 vs 流水线) |
|---|---|---|
| Total_TG | 3.4×10⁻³ | 225 vs 220 |
| MUFA | 8.6×10⁻³ | 183 vs 181 |
| Omega_3 | 1.5×10⁻² | 154 vs 152 |
为什么不是浮点一致：我的独立等位对齐比 TwoSampleMR harmonise(action=2) 多留了 2–5 个回文 SNP → IVW 估计差 ~0.3–1.5%。这说明结果对个别回文 SNP 不敏感、稳健，但不像 eQTL/meQTL 那样精确到浮点。要做到逐位一致需完全复刻 TwoSampleMR 的回文处理，价值不大。
3.4 视网膜 eQTL 层：也已证明（tools/verify_retina.R）
独立重解析 Advani 2024 视网膜 eQTL（重打变异解析 + eaf 对齐 + plink clump + Wald/IVW），对头号命中重算：
| 靶点×AMD | 方法 | b 相对差 | se 相对差 | 工具数 |
|---|---|---|---|---|
| CFH | Wald | 1.1×10⁻¹⁵ | 5.8×10⁻¹⁶ | 1 = 1 |
| MERTK | IVW | 6.5×10⁻¹⁶ | 1.1×10⁻¹⁶ | 2 = 2 |
浮点级吻合、工具数精确一致 → 视网膜层也证明算对了。（关键是补上流水线 MR 前的 clump 步骤：CFH 4 个显著 SNP→1、MERTK 110→2。）
四、功能验证代码静态审查（dr-amd-targets，本次审 10 个核心脚本）
| 脚本 | 方法 | 结论 |
|---|---|---|
13_amd_pseudobulk.R | 逐样本伪 bulk CPM + Wilcoxon/Welch | ✅ 正确，规避细胞级伪重复(Squair 2021)——正确推翻了假阳性 |
31_gct_nsr_vs_rpe.R | limma-voom + edgeR RPE vs NSR | ✅ RNA-seq DE 金标准，严谨 |
11_dr_scrna.R | Seurat + Harmony 整合、模块打分注释 | ✅ 标准流程；自动 argmax 注释偏粗，发表应人工核标记 |
12_amd_scrna_rpe.R | 同上 + 细胞级 FindMarkers | ✅ 标准；细胞级 FindMarkers 即被 13 伪 bulk 推翻的假阳性，结论未采用它 |
25_dr_slingshot.R | Slingshot 锚定 homeostatic、Spearman | ✅ 轨迹设置正确；伪时序本属探索性(无 p 值) |
21_dr_cellchat.R | CellChat triMean 通讯概率 | ✅ 标准；通讯概率非因果,探索性 |
02_wgcna_gsea.R | WGCNA 模块-性状 + 单基因 gseGO | ✅ 样本<12 自动跳过(DR 9 样本)，诚实不过度解读 |
03_gsea_fgsea.R | fgsea GO-BP，eps=0 精确 p | ✅ 标准 GSEA(与 02 gseGO 互为交叉) |
30_eqtl_coloc.R | 视网膜 eQTL(EGA METR)×疾病 coloc | ✅ 等位对齐正确、nsnp<5 返 NA |
01_bulk_analysis.R | 微阵列 DEG，base R t 检验 | 🟡 逻辑/符号对，但用基础 t 检验而非 limma moderated-t(微阵列发表标准)，功效偏低 |
功能验证唯一实质短板：01(微阵列 bulk DEG)应改 limma。AMD 侧微阵列本就阴性(且伪 bulk 佐证)、DR 阳性有单细胞交叉支持，影响有限；但发表前建议对微阵列也上 limma。未审的 04(msigdb GSEA 变体)/05(出图)/23(monocle,已被 slingshot 取代) 不承载独立结论。
五、够不够发表高分 SCI？
方法学够格：cis-MR + coloc.abf + SMR/HEIDI + 多效性全套 + BH-FDR + limma-voom + 伪 bulk，都是领域标准且实现正确；多组学头号命中已被独立证明。没有发现会翻转主结论的错误。
投稿前收尾清单（2026-07-25 已全部处理）：
- ✅ SMR 被剔 SNP 数已列（AMD 1,023 / DR 982 / IAMDGC 534，见 04·甲基化SMR）；US_Blood=Understanding Society，按下载页引用。
- ✅
01微阵列 DEG 已换 limma moderated-t：AMD 侧仍全部不显著(伪 bulk 佐证)、DR 侧显著性维持(MERTK FDR 0.018→0.012、SPRY2 →0.0008、TNFRSF10A →0.003、CFH →0.0001)，结论不变，靶点证据总表 已更新。 - ✅ eQTL FDR 已把 IAMDGC 单列为 replication 族：discovery(FinnGen) 与 replication(IAMDGC) 分别 BH；CASP10 0.0186→0.0163、TNFRSF10A 0.00139→0.00125，仍强显著。
- ✅ 参考文献已联网核对：eQTLGen(Võsa 2021 Nat Genet)、UKB-PPP(Sun 2023 Nature)、IAMDGC(Gorski/Grunin IOVS 2025)、Chen(Nat Genet 2023 PMID 36635386)、Nightingale 619k(Tambets R, et al. Nature 2026;655(8124):971-978, PMID 42162431) 均确认；仅 Z→β 近似式无单一公认出处（建议引 Zhu 2016 SMR 补充材料，公式本身已复算验证正确）。
- novelty 风险(TNFRSF10A 等被 2024–25 覆盖)靠 trans/多组学差异化 —— 属选题策略，非代码问题。
六、独立复算覆盖度（截至 2026-07-25）
| 层 | 静态审查 | 独立复算 |
|---|---|---|
| eQTL(CASP10/TNFRSF10A) | ✅ | ✅ 浮点级(1e-15) |
| meQTL SMR(APOE/TNFRSF10A/CASP10) | ✅ | ✅ 公式级(1e-6)+ GWAS 对原始核 |
| 代谢物(Total_TG/MUFA/Omega_3) | ✅ | ✅ ~1%(工具集差 2–5 SNP) |
| 视网膜 eQTL(CFH/MERTK) | ✅ | ✅ 浮点级(1e-15) |
| 功能验证 13 个脚本(全核心 + 04/05/23) | ✅ | —(非 MR 型，看方法正确性) |
多组学四层(eQTL/meQTL/代谢物/视网膜)独立复算已全覆盖。 功能验证 04(msigdb GSEA 变体)/05(GSEA 出图，已自带 PNG+PDF)/23(monocle3，Slingshot 的第二方法交叉验证) 已补审，标准合格；剩余仅装包/助手脚本不承载结论。
主 MR 共定位/LDSC 独立复算（2026-07-25，见 方法学审计 · 十三）：
- 共定位跨方法交叉验证 ✅：coloc.abf 对 TNFRSF10A pQTL×AMD 得 PP.H4=0.9994（共享 SNP 5612 vs hyprcoloc 5609），印证 hyprcoloc 支持判定。
- LDSC ⚠️→✅ 发现并已修 rg 的 se/z/p 被低估（delta 法 0.024 vs GenomicSEM jackknife 0.1594，rg 近 1 时低估）：
04_ldsc.R已改用原生 jackknife se、两家族 rg 表已重出（AMD↔WetAMD z 39.22→5.84），点估计/h2 本就无误。
脚本：tools/verify_multiomics.R(eQTL+coloc)、verify_layers.R(meQTL+代谢物)、verify_retina.R(视网膜)、mr-pipeline/tools/verify_coloc_pqtl.R、verify_ldsc.R、fix_ldsc_rg_se.R。