使用 SMR 整合 GWAS 与 eQTL

SMR（summary-data-based Mendelian randomization）整合 GWAS 与分子 QTL summary data，检验由遗传工具预测的分子性状——例如 gene expression——是否与复杂性状相关。

常见问题是：某个基因的 cis-eQTL 信号是否与疾病 GWAS 信号呈现一致关联。SMR 提供关联检验，HEIDI 使用区域内多个 variants 检查这种模式是否更像单一共享变异或 pleiotropy，而不是两个处于 LD 中的不同 causal variants。

SMR + HEIDI 与 Bayesian colocalization 不是同一种方法。它们可以互相补充，但不能把 SMR 结果直接写成“已经证明共定位或因果”。

三类输入

一次基础分析需要：

输入	作用	常见格式
GWAS summary statistics	复杂性状的 SNP association	`.ma`
Molecular QTL summary data	expression、protein 等分子性状 association	BESD 前缀
LD reference	SNP 间相关结构	PLINK `.bed/.bim/.fam` 前缀

三者应使用一致或可以可靠协调的 genome build、variant identifiers 与 allele definitions。LD reference 的遗传祖源还应尽量匹配 GWAS 与 QTL 研究人群。

准备 GWAS `.ma` 文件

SMR 的 GWAS summary 文件需要明确列名和字段含义：

SNP  A1  A2  freq  b  se  p  N

以 FinnGen 风格数据为例，可以先在 R 中筛选并重命名：

library(data.table)
library(dplyr)

outcome_gwas <- fread(
  "Outcome_Gwas/finngen_R10_C3_SQUOMOUS_CELL_CARCINOMA_SKIN_EXALLC.gz",
  data.table = FALSE
)

gwas_ma <- outcome_gwas |>
  transmute(
    SNP = rsids,
    A1 = alt,
    A2 = ref,
    freq = af_alt,
    b = beta,
    se = sebeta,
    p = pval,
    N = 317724
  ) |>
  filter(
    !is.na(SNP),
    SNP != "",
    !grepl(",", SNP)
  )

fwrite(
  gwas_ma,
  file = "Outcome_ma/SCC.ma",
  sep = "\t",
  quote = FALSE,
  row.names = FALSE
)

这里的 A1 必须对应 b 的 effect allele，freq 也应是同一个 allele 的频率。不能只按列位置猜测含义。

写出后重新读入检查：

check_ma <- fread("Outcome_ma/SCC.ma")

names(check_ma)
head(check_ma)
summary(check_ma)

同时检查重复 SNP、非法 allele、频率范围、极端标准误和样本量是否按 variant 变化。

全区域运行 SMR

假设：

EUR/EUR 是 PLINK LD reference 前缀；
Whole_Blood.lite 是 SMR 格式的 eQTL summary 前缀；
Outcome_ma/SCC.ma 是 GWAS 文件。

smr-1.3.1-win.exe \
  --bfile EUR/EUR \
  --gwas-summary Outcome_ma/SCC.ma \
  --beqtl-summary Whole_Blood.lite \
  --out results/scc_whole_blood

--bfile 接受的是 PLINK 文件共同前缀，而不是目录名。运行前应确认 .bed、.bim 和 .fam 三个文件同时存在。

只分析指定 probes

将需要的 probe IDs 写入纯文本文件，每行一个 ID：

ENSG00000228794
ENSG00000188976
ENSG00000187961

从已有 BESD 中提取：

smr-1.3.1-win.exe \
  --beqtl-summary Whole_Blood.lite \
  --extract-probe myprobe.list \
  --make-besd \
  --out selected_probes

再运行分析：

smr-1.3.1-win.exe \
  --bfile EUR/EUR \
  --gwas-summary Outcome_ma/SCC.ma \
  --beqtl-summary selected_probes \
  --out results/scc_selected_probes

.list 文件应使用纯文本，每行只有一个 ID，不带表头、引号或多余空格。

只保留指定 SNPs

同样可以准备 rsID 列表：

rs559807721
rs761651383
rs9726668

smr-1.3.1-win.exe \
  --beqtl-summary Whole_Blood.lite \
  --extract-snp mysnp.list \
  --make-besd \
  --out selected_snps

筛选 SNP 前应确认这些 rsIDs 存在于 eQTL、GWAS 和 LD reference 的可协调集合中。只在某一个输入中出现的 SNP 不会形成有效分析信息。

从 ESD 文件批量构建 BESD

自有 QTL summary data 可以先整理为 ESD。每个 probe 的文件包含：

Chr  SNP  Bp  A1  A2  Freq  Beta  se  p

在 R 中创建一个 probe 的 ESD：

esd <- exposure_gwas |>
  transmute(
    Chr = sub("^chr", "", chromosome),
    SNP = rsid,
    Bp = position,
    A1 = effect_allele,
    A2 = other_allele,
    Freq = effect_allele_frequency,
    Beta = beta,
    se = se,
    p = p_value
  ) |>
  filter(
    !is.na(SNP),
    SNP != "",
    !grepl(",", SNP)
  )

fwrite(
  esd,
  file = "esd_cis/ENSG00000188976.esd",
  sep = "\t",
  quote = FALSE,
  row.names = FALSE
)

再建立 .flist，描述每个 probe 及其 ESD 路径：

Chr  ProbeID  GeneticDistance  ProbeBp  Gene  Orientation  PathOfEsd

flist <- data.frame(
  Chr = 1,
  ProbeID = "ENSG00000188976",
  GeneticDistance = 0,
  ProbeBp = 1000000,
  Gene = "GENE1",
  Orientation = "+",
  PathOfEsd = "esd_cis/ENSG00000188976.esd"
)

fwrite(
  flist,
  file = "all_flist.flist",
  sep = "\t",
  quote = FALSE,
  row.names = FALSE
)

实际批量处理中，probe position、strand、gene symbol 与 cis window 应从明确版本的 gene annotation 获得。不要用 gene length 代替 genetic distance，也不要在没有记录 genome build 的情况下混合坐标。

构建 BESD：

smr-1.3.1-win.exe \
  --eqtl-flist all_flist.flist \
  --make-besd \
  --out custom_eqtl

然后作为 --beqtl-summary custom_eqtl 运行 SMR。

解释 SMR 与 HEIDI

结果中需要把两个问题分开：

SMR test：遗传预测的分子性状是否与复杂性状相关？
HEIDI test：区域内多个 variants 的效应模式是否表现出异质性？

较小的 SMR p-value 支持 association，但仍可能受到 linkage、horizontal pleiotropy、样本重叠或输入 summary statistics 偏倚影响。较小的 HEIDI p-value 通常提示异质性，与“单一共享变异解释”不一致；未拒绝 HEIDI 也不等于证明共享 causal variant。

筛选结果时应同时考虑：

分析了多少 probes，以及 multiple testing 如何控制；
HEIDI 实际使用了多少 SNPs；
top eQTL 的强度；
GWAS 与 eQTL 的 allele 是否一致；
LD reference 与研究祖源是否匹配；
MHC 等复杂 LD 区域是否需要单独处理；
组织或细胞类型是否与研究问题相关。

最小输出记录

为了复现分析，应保留：

SMR executable 版本；
GWAS 数据来源、版本、build 和样本量；
QTL 数据来源、组织、版本和 build；
LD reference 的来源与祖源；
.ma、BESD 和筛选列表的生成代码；
命令行参数；
SMR 与 HEIDI 的完整输出，而不只保留显著行。

SMR 的优势是能直接利用大规模 summary data 整合分子 QTL 与复杂性状，但结果仍是遗传关联证据。将其与独立 colocalization、fine-mapping、不同组织 QTL 和功能实验结合，才能逐步缩小机制解释空间。

参考：SMR software and documentation；SMR data resources。