前阵子帮一位做复杂疾病遗传学的同事看数据他的GWAS在某个基因间区落下了一个非常漂亮的关联峰曼哈顿图上闪闪发亮但整段区域几乎找不到错义突变也没有能直接解释功能的位点。审稿人一句“please provide functional interpretation”直接把他卡住了。这种反馈近两年越来越常见——光有显著性已经不够了你得回答“这个信号到底通过哪个基因、多大程度上影响了性状”。而SMRSummary-data-based Mendelian Randomization就是干这件事的把GWAS的关联信号和公共eQTL数据放在同一个孟德尔随机化框架里用汇总数据判断基因表达对性状的因果效应同时用HEIDI检验把“真因果”和“连锁导致的假关联”区分开。这篇文章不打算铺开讲数学推导重点是我从数据准备到结果解读完整跑过一遍之后的实操经验包括那些文档里没写的坑。1. SMR要回答的问题从一条GWAS信号到基因调控证据1.1 一条关联峰不等于一个因果基因很多刚接触SMR的人会有一个误区觉得GWAS某个区域最显著的SNP离哪个基因最近那个基因就是“致病基因”。但真实情况远没那么简单。我自己做过一个数据GWAS信号落在两个基因之间的非编码区按距离硬排最接近的是一个长链非编码RNA但后续用表达数据去看真正被这个信号调控的其实是下游30kb外的一个蛋白编码基因。这种情况在文献里太常见了有人统计过相当一部分风险位点并不会映射到物理距离最近的基因。这就是为什么需要功能基因组数据介入。eQTLexpression quantitative trait locus告诉我们某个SNP的基因型会改变某个基因的表达量。GTEx、eQTLGen这些公共数据库积累了不同组织里大量的cis-eQTL信息。但问题来了——GWAS说这个SNP和疾病相关eQTL说这个SNP和某个基因的表达相关两件事都是相关性怎么从“相关”走到“因果”SMR提供的就是这一步用孟德尔随机化的思路把基因表达当作暴露把疾病当作结局用cis-eQTL的SNP作为工具变量从两个独立GWAS汇总数据中估计“基因表达改变一个单位疾病风险改变多少”。因为等位基因在受精时是随机分配的相当于一次天然的随机分组实验混杂因素影响比普通观察性研究小得多。1.2 孟德尔随机化怎么套进汇总数据传统孟德尔随机化通常需要个体水平的基因型和表型数据拿到手之后跑两阶段回归。SMR的精妙之处在于它只需要汇总统计量就能完成类似的分析。核心逻辑可以简化为两步。第一步从eQTL数据中得到工具变量SNP对基因表达的效果值记为β_eQTL第二步从同一SNP的GWAS结果中得到它对性状的效果值记为β_GWAS。假如基因表达真的是中间机制那么这两个效应的比值就是基因表达对性状的因果效应即β_SMR β_GWAS / β_eQTL。显著性检验则用delta方法来估计这个比值的标准误。当然这个推导成立的前提是工具变量有效也就是cis-eQTL的SNP只能通过改变基因表达来影响疾病不能有水平多效性。这个假设在实践中很难完美满足所以SMR的作者才配套设计了HEIDI检验来兜底。1.3 HEIDI检验到底在筛什么HEIDIHeterogeneity in Dependent Instruments这个名字翻译过来是“依赖工具变量的异质性检验”听起来很学术思想其实很直白。如果一个因果模型是真的——某个基因表达导致疾病而且cis-eQTL区间里的SNP都是通过影响这个基因表达来起作用——那么不管你拿区间里哪个SNP去算β_SMR得到的因果效应估计都应该差不多落在随机误差范围内。但如果GWAS信号和eQTL信号只是因为连锁不平衡LD碰巧“共享”了同一个显著位点而不是同一个因果变异驱动那么不同SNP各自算出来的β_SMR就会散得很开产生显著异质性。HEIDI检验就是把这种异质性量化出来给一个p值。p_HEIDI显著比如小于0.05说明这些SNP的效应不一致你的阳性结果更可能是连锁造成的假信号p_HEIDI不显著说明没有检测到异质性结果支持基因表达对性状的因果效应。这就是SMR比单纯共定位分析强的地方不仅告诉你“是不是同一个信号”还顺带给了因果效应的估计值。1.4 SMR、共定位与经典MR一张表看明白我刚跑这个流程的时候最困惑的就是SMR、coloc、还有传统MR之间到底什么关系后来自己画了个对比表才彻底理清。方法输入数据核心问题输出结果典型应用场景SMR HEIDIGWAS汇总数据 eQTL/mQTL/pQTL汇总数据基因表达是否因果影响性状因果效应估计 异质性检验功能注释、药物靶点验证coloc共定位GWAS汇总数据 eQTL汇总数据两个信号是否共享同一个因果变异后验概率PP.H4等判断关联峰是否来自同一变异传统MR个体水平数据或汇总数据暴露是否因果影响结局因果效应估计危险因素因果推断通常暴露是性状而非分子实际项目里我会把SMR和coloc配合使用SMR给出效应方向和大小coloc验证共享变异的概率。两者结论一致时结果非常硬不一致时优先排查数据质量和LD参考面板的问题而不是急着下结论。2. 开跑前的三套数据GWAS汇总数据、eQTL数据与LD参考面板的选配2.1 GWAS汇总数据的最低要求SMR对GWAS数据的要求说高不高说低也不低。最低限度你得有这八列SNP ID、染色体、位置、效应等位基因A1、另一个等位基因A2、效应量Beta、标准误Se、P值。样本量N建议也要有虽然SMR核心计算不一定强依赖N但后续质控和审稿回复时没有会很被动。这里有个高频错误我必须单独拎出来说如果你下载的GWAS数据里给的是OR值一定要先转换成log(OR)再当Beta用不能直接塞进去。我见过有人拿OR直接跑出来的β_SMR大得离谱还以为发现了重要基因。另外A1必须是效应等位基因这个方向定错了后面所有结果符号都是反的。很多公开GWAS数据没有直接给出等位基因频率SMR的MA格式里Freq这一列可以填也可以不填但建议还是补上。频率信息有两个用途一是质控时用--diff-freq检查GWAS和eQTL数据里的等位基因是否一致二是遇到链方向问题时能多一个判断依据。如果你手里确实没有频率可以先不做这一步但后续遇到可疑结果时一定要回过头来补。2.2 eQTL数据选哪套要看组织和人群eQTL数据的选择直接决定你的分析有没有生物学意义。同样是跑SMR用血液eQTL解释神经精神疾病用肝脏eQTL解释血脂性状逻辑上说得通吗不完全是。我通常会按组织相关性选。数据源组织/来源样本规模适用场景eQTLGen全血约31,684人血液相关性状、免疫疾病、目前样本量最大的cis/trans-eQTLGTEx v849个组织838个供体多组织功能注释目标组织明确的场景BrainMeta / BRAINEAC脑组织约1,000神经精神疾病、认知性状PsychENCODE脑组织约1,000精神疾病、脑组织eQTLSMR官网预处理数据多组织视数据库而定直接下载BESD格式省去格式转换从实操角度SMR官网提供了大量已经处理成BESD格式的eQTL数据包括GTEx v7/v8、eQTLGen、BrainMeta等能下载就直接下载省掉很多格式转换的麻烦。但如果你的研究人群和这些数据库来源人群差异很大或者需要的组织没有覆盖那就只能自己拿原始eQTL文件转。人群匹配度也很重要。eQTLGen和GTEx主要都是欧洲裔样本用来注释欧洲人群GWAS结果很顺手如果是东亚人群的GWASLD参考面板选EASeQTL层面也尽量找东亚人群的数据库找不到时至少要在方法和局限性里诚实说明。2.3 LD参考面板最容易忽略却最致命的一环很多人第一次跑SMR会漏掉LD参考面板心想“我都给齐GWAS和eQTL了为什么还要一个PLINK格式的参考面板”因为HEIDI检验要算多个SNP之间的连锁不平衡结构没有LD信息这一步就成了空中楼阁。建议直接用千人基因组Phase 3样本量2504人覆盖五大洲人群分EUR、EAS、AFR、AMR、SAS等子面板已经是这个领域的默认配置。下载PLINK二进制格式bed/bim/fam按染色体拆分存储后续分析按需调用。需要注意LD参考面板的基因组版本必须和你的GWAS数据一致GWAS是GRCh37/hg19面板也得是hg19GWAS是GRCh38面板就得换成hg38。版本对不上轻则SNP匹配率低重则结果完全不可信。有些人会问能不能用自己项目的个体基因型做LD参考可以而且理论上更贴近你的研究人群只要样本量够至少几百人。但对多数公共数据研究来说千人基因组已经足够。3. 文件准备三步走MA格式、BESD格式与探针清单3.1 GWAS原始文件转MA格式SMR的GWAS输入是一种简单的文本格式官方叫MA格式本质上就是带表头的制表符分隔文件。列名必须是SNP Chr Pos A1 A2 Freq Beta Se P N这种固定命名SMR靠表头认列顺序可以打乱但名字不能改。从原始GWAS文件转MA格式最稳妥的方式是用awk原地处理。假设你的原始文件列依次是rsid、chr、pos、effect_allele、other_allele、beta、se、p、n命令就这样写awk BEGIN{OFS\t; print SNP, Chr, Pos, A1, A2, Freq, Beta, Se, P, N} NR1 {next} {print $1, $2, $3, $4, $5, NA, $6, $7, $8, $9} gwas_raw.txt gwas.ma这里我把Freq填了NA如果你的数据里有频率列就替换成对应字段。跑完之后一定要做两步检查第一wc -l看行数防止表头处理和换行符问题第二head看前几行确认没有列错位。这一步看着简单实际翻车率极高我见过不止一次因为原始文件列顺序和预想不一样导致最终结果全是垃圾的情况。另外一个容易踩的点原始GWAS文件里如果有重复的rsIDSMR会直接报错或者随机取一条最好先做去重。去重时保留P值最小的那一行是比较合理的策略。3.2 eQTL摘要数据转BESD格式如果SMR官网没有现成的BESD数据你需要自己转。BESD是SMR的二进制eQTL格式包含besd和besi两个文件转换命令非常简单smr --eqtl-summary eqtl_clean.txt --make-besd --out myeqtl真正费时间的是准备eqtl_clean.txt这个文本文件。SMR要求的eQTL文本格式包含这些列列名含义示例SNP变异IDrs123456Chr染色体1Pos位置bp1000000A1效应等位基因AA2另一个等位基因GFreqA1频率0.23Probe探针ID一般是Ensembl基因IDENSG00000123456Gene基因名BRCA1Ori基因方向 或 -BetaeQTL效应量0.15Se标准误0.03PP值1.2e-20这里最麻烦的是Ori列。GTEx和SMR官网预处理数据里的Ori是按转录本方向填的如果你拿到的eQTL原始数据没有这个信息可以统一填。但要注意如果真实方向是负链强行填可能导致部分SNP的效应方向被翻转最终结果里这些信号的β_SMR符号可能异常。稳妥起见如果数据源本身没有方向信息我会在转换后随机抽几个已知阳性信号做验证确认方向没被系统性翻转。转换前还有一个强烈建议先过滤。eQTL原始文件通常极大尤其eQTLGen的trans-eQTL部分动辄几千万行直接转BESD慢不说P值列为NA或0的坏行还会中断转换。我一般先按P 1e-5或P 5e-8筛一遍只保留显著信号再转换。SMR做的是cis-eQTL区间内的分析显著SNP之外的弱信号对结果贡献不大但能大幅缩小文件体积。3.3 探针清单控制分析范围SMR允许通过--probe参数指定一个探针清单文件只分析你关心的基因。一行一个探针ID可以是Ensembl ID也可以是eQTL数据里的原始探针ID比如Illumina的ILMN_开头必须和你BESD里的Probe列值能对上。我强烈建议你们不要图省事直接跑全基因组所有探针。原因有两个第一多组织数据全量跑一遍耗时很长尤其LD面板按全基因组读取每次都要加载几个G的文件第二分析几万个基因意味着几万次检验多重检验校正后阈值会严到离谱很多真实的弱信号会被一并砍掉。实际操作中先拿GWAS显著位点去refGene或UCSC注释一下把位点上下游1Mb范围内的基因都选出来建一个几十到几百个探针的清单精确打击。我自己常用的一种做法是先用--extract-snp限定GWAS显著位点及其LD伙伴附近的SNP再用--probe限定候选基因两个限制条件一起上速度能快一个数量级而且结果更容易解释。4. 跑通SMR与HEIDI命令行参数与批量实战4.1 基础命令逐项拆解当数据文件都准备好后正式分析只有一条命令。以欧洲人群GWAS加GTEx某个组织的BESD为例smr \ --bfile /path/to/1kg_phase3_eur \ --gwas-summary gwas.ma \ --beqtl-summary myeqtl \ --extract-snp sig_snps.snplist \ --probe candidate_probes.txt \ --out smr_result \ --thread-num 8每个参数我都解释一下方便你不看文档也能快速上手--bfileLD参考面板的PLINK前缀不含.bed后缀。这个文件提供SNP之间的LD结构HEIDI检验依赖它。--gwas-summary上一步做好的GWAS MA格式文件。--beqtl-summaryBESD格式的eQTL数据前缀也不需要写.besd后缀。--extract-snp只分析这个列表里的SNP。我一般把GWAS显著位点上下游1Mb窗口内的所有SNP都放进去既覆盖了cis-eQTL区间又避免了全基因组扫描的耗时。--probe探针清单限定要分析的基因。--out输出文件前缀。--thread-num线程数多核机器上能明显提速。运行过程中SMR会在屏幕打印每个探针的分析进度也会提示有多少SNP没能匹配LD面板。如果匹配率低于80%建议停下来检查SNP ID的版本是否统一而不是继续跑下去。4.2 几个值得设置的进阶参数除了基础参数有几个参数在实际项目中经常用到我整理成表格供参考。需要说明的是不同SMR版本的默认值略有差异跑之前最好smr --help再确认一遍。参数作用我的常用设置--heidi-method选择HEIDI检验的实现方法默认方法即可一般用0--heidi-min-nsnpHEIDI检验最少需要多少个SNP3低于这个值不输出HEIDI结果--peqtl-threshold筛选工具变量SNP的eQTL P阈值5e-8或者1e-5数据量大时用更严的阈值--diff-freqGWAS与eQTL之间等位基因频率差异上限0.1到0.2用于质控--diff-freq-prop允许有多少比例SNP超过频率差异阈值0.05--maf最小等位基因频率过滤0.01--peqtl-threshold这个参数值得多说一句。SMR默认会用eQTL数据里所有可用的SNP来构建工具变量但如果你拿到的eQTL数据保留了大量不显著的SNP自己转BESD时没做过滤就会出现这种情况它们不是好的工具变量反而会增加HEIDI检验的噪音。我会根据样本量调整eQTLGen这种大样本数据cis-eQTL可以压到5e-8GTEx某些样本量较小的组织压到1e-5就行太严会损失大量有效探针。4.3 多组织批量跑法与提速技巧做功能注释的时候几乎没人会只看一个组织。GTEx v8有49个组织挨个手敲命令不现实写个循环脚本就解决了for tis in Brain_Cortex Brain_Cerebellum Liver Lung Whole_Blood; do smr \ --bfile /path/to/1kg_phase3_eur \ --gwas-summary gwas.ma \ --beqtl-summary /path/to/GTEx_v8/${tis} \ --probe candidate_probes.txt \ --out smr_${tis} \ --thread-num 4 done wait几个提示后台并行任务不要开太多否则内存和IO都会成为瓶颈一般同时跑四到六个刚好组织名单里不要有空格不然循环变量会拆开报错每跑完一个组织看一眼.log文件确认没有大量SNP匹配失败。如果LD参考面板是全基因组体量的千人基因组每次加载都挺耗时。加速的关键是把LD面板按染色体拆分然后配合--extract-snp只加载目标区域的SNP。更极端的做法是用PLINK先按区域提取一份小的bed文件再把--bfile指向它速度能快上好几倍。5. 结果怎么读p_SMR、p_HEIDI与一张火山图5.1 .smr文件长什么样跑完后SMR会输出一个以.smr结尾的结果文件每行是一个探针的分析结果关键列如下列名含义probeIDeQTL探针IDgene基因名或Ensembl IDchr染色体probe_bp探针位置topSNP最显著的eQTL SNPnsnp参与分析的SNP数量FreqtopSNP的频率Beta_SMR基因表达对性状的因果效应估计SE_SMRBeta_SMR的标准误p_SMRSMR检验的P值Beta_HEIDIHEIDI检验的效应异质性估计p_HEIDIHEIDI检验的P值nsnp_HEIDIHEIDI检验使用的SNP数量举例来说某一行可能长这样ILMN_1234567 ENSG00000123456 1 1000000 rs123456 50 0.35 0.42 0.11 1.3e-7 0.08 0.23 12解读起来就是这个基因的eQTL信号和GWAS信号之间存在显著的SMR关联p_SMR 1.3e-7同时没有检测到明显的异质性p_HEIDI 0.23是一个值得关注的候选基因。5.2 什么才算可信的阳性结果拿到结果列表后不要直接抓p_SMR最小的基因就完事。我的筛序逻辑一般是这样的第一步看p_SMR是否达到多重检验校正后的阈值。Bonferroni最简单的做法是0.05 / 探针数比如分析了200个基因阈值就是2.5e-4。嫌严格可以上FDRq值小于0.05即可。第二步看p_HEIDI是否大于0.05。这里记住一个容易绕晕的点p_SMR要小p_HEIDI要大。p_HEIDI大说明没有异质性支持因果关联p_HEIDI小说明可能是连锁造成的假信号。第三步看效应方向。如果某个已经报道过的药物靶基因GWAS和eQTL的效应方向是反的说明你的基因型到表达再到表型这条链条中间可能出了问题哪怕统计显著也要打问号。第四步看nsnp_HEIDI。如果参与HEIDI检验的SNP太少小于3个这个p_HEIDI的可靠性就很有限结果只能作为候选不能当作确定因果。还有一个容易被忽视的点SMR的阳性结果里经常出现一个探针对应多个基因名的情况特别是同一个基因的不同转录本。如果某个基因有多个探针行建议取p_SMR最显著的那一行做代表其余标注为同一基因的重复信号避免后续统计把基因数重复计算。5.3 用R画一张能直接放进论文的火山图结果文件本身是文本肉眼扫几百行不现实。我习惯用R快速画一张火山图横轴是Beta_SMR纵轴是-log10(p_SMR)颜色表示HEIDI是否通过。代码很简单library(ggplot2) res - read.table(smr_result.smr, header TRUE, stringsAsFactors FALSE) res$pass_heidi - ifelse(res$p_HEIDI 0.05 !is.na(res$p_HEIDI), HEIDI pass, HEIDI fail) res$logP - -log10(res$p_SMR) sig - res[res$p_SMR 5e-4 res$pass_heidi HEIDI pass, ] ggplot(res, aes(x Beta_SMR, y logP, color pass_heidi)) geom_point(alpha 0.7, size 1.8) ggrepel::geom_text_repel(data sig, aes(label gene), size 3) scale_color_manual(values c(HEIDI pass #C0392B, HEIDI fail #7F8C8D)) geom_hline(yintercept -log10(5e-4), linetype dashed, color grey30) theme_bw(base_size 14) labs(x Beta_SMR, y -log10(p_SMR), color )右上角那些又红又高的点就是你最需要关注的候选基因。如果一张图上红色点稀少也别急着失望SMR本身就是偏保守的工具假阳性率控制得比较好阳性少是常态。6. 翻车记录格式错位、链方向与参考面板的隐性坑6.1 等位基因链方向不一致的典型症状与排查SMR虽然会自动尝试翻转等位基因实现对齐但A/T和C/G这种互补对称的SNP在翻转后模棱两可很容易出问题。我遇到过最典型的一次跑出来一批p_SMR很小但Beta_SMR符号全反的“阳性基因”一开始以为是发现了新的负向调控后来发现是原始GWAS文件用的负链和eQTL数据用的是正链两种数据在链方向没有对齐。排查链方向问题有几个信号结果里大量基因的Beta_SMR符号异常同一个基因在多个组织里符号不一致--diff-freq质控提示明显偏多的SNP频率差异超过阈值。处理办法是回到原始数据确认GWAS和eQTL文件各自用的参考基因组链和参考等位基因编码方式在清洗阶段统一为正链。这个工作必须在转换MA格式之前做跑完SMR再发现问题只能回头全部重来。6.2 基因组坐标build不一致看似能跑实则全废某个公共数据库的eQTL数据是GRCh38坐标而你手里的GWAS和千人基因组LD面板还是GRCh37这种情况下SMR也能“跑通”因为SNP ID很多是rsID能匹配上不少真正出问题的是一些依赖坐标的SNP会被错误配到别的区间。最要命的是这个问题不容易一眼发现结果看起来一切正常但HEIDI检验用的LD结构可能已经错得离谱。我的检查习惯是拿到任何一批数据先随机抽几个基因找到这些基因在GWAS和eQTL文件里的位置和UCSC/Ensembl上的基因组位置比对一下。如果整体偏差在几百bp内属于正常如果差了几百万甚至更多基本可以确定是build不一致用liftOver统一版本后再继续。6.3 rsID版本过老与染色体命名不一致还有一个高频问题rsID版本不匹配。GWAS用的是dbSNP 151的IDLD面板里是dbSNP 150的ID有些位点的rsID在两个版本间被合并或拆分导致SMR匹配不上。这种情况没有一劳永逸的解决办法能在数据准备阶段统一最好已经发生的话可以拿chr:pos:alleles格式做二次匹配SMR的老版本对这类变通支持不算好建议尽量用新版本。染色体命名也是一样。有的文件是chr1有的是1有的文件是1:1000000有的是rs123。最好在清洗阶段就统一成纯数字染色体加rsID或chr_pos_ref_alt的形式避免SMR内部解析出错。6.4 eQTL数据里的N多“NA”行和不显著行这个问题在自己转BESD文件时特别突出。一些公共eQTL数据会保留所有检测过的SNP哪怕基因型和表达量的关联检验根本没有通过P值也可能是一个很大的数或者直接是NA。SMR在转换过程中遇到这些行会报错或者悄悄丢弃结果就是你原本期待分析一万个基因转完BESD只剩三千个白白损失信息。所以我建议在转BESD之前就对原始eQTL文件做一次完整的QC剔除P值为NA或0的行剔除MAF过小的SNP按P 0.05或者更严格的阈值过滤。虽然理论上SMR的--peqtl-threshold参数也能控制进入分析的SNP阈值但文件清洗层面提前做一遍既能减小BESD体积又能避免转换中断。6.5 每次跑前都用一个“已知阳性”验证流程最后这点是我踩过几次坑后养成的习惯正式分析之前先挑一两个文献里已经反复验证过的基因当作阳性对照跑一遍SMR。比如做血脂相关分析就拿LDL受体或者PCSK9相关的位点试试做炎症性疾病就找IL-23R、TNFAIP3这类已知基因。如果连这些公认的阳性对照都跑不出来基本可以确定是上游数据或参数出了问题。此时先别急着改参数老老实实检查SNP匹配率、等位基因方向、坐标版本这三件事多半能找到罪魁祸首。用阳性对照验证流程就像装机后先跑一遍内存测试再投入使用省下的时间远远超过多花的这几分钟。对我个人来说凡是涉及公共数据的分析一定会先花半天时间把GWAS、eQTL、LD面板三者之间的SNP ID、等位基因和坐标版本彻底对齐再开始正式分析。这些准备工作虽然不产生任何漂亮结果但决定了后面所有输出是否可信。SMR这套流程跑通一次之后后续换数据、换组织、换疾病都很轻松真正慢的永远是第一次的数据清洗和格式统一。