入行这几年我处理过不少GWAS数据也回答过很多次“这个显著位点到底是怎么影响性状的”这类问题。早期大家拿到一个显著SNP第一反应是找最近的基因然后把它当重点基因写进文章里。这个做法在部分场景下没问题但现实中GWAS信号往往落在非编码区最近基因和真正起调控作用的基因并不总是一回事。所以这几年我越来越依赖一个工具就是SMRSummary-data-based Mendelian Randomization。它能把GWAS的汇总统计数据和eQTL数据放在一起分析在基因水平上快速筛选出与性状可能存在因果关联的基因而且整个过程只需要公开的汇总数据不需要拿到个体的基因型和表达量。这篇文章我想把SMR从原理到实操完整地过一遍。内容包括SMR到底在做什么、HEIDI检验为什么重要、软件和参考数据怎么准备、GWAS数据怎么转成SMR要求的格式、完整运行命令怎么敲、输出文件怎么读、有哪些坑我替你们先踩过了。如果你手上有GWAS的summary数据想往下做功能注释或者找候选致病基因这套流程可以直接照着走。1. SMR到底在做什么原理与思路拆解很多第一次接触SMR的人容易被“孟德尔随机化”这个名字吓到觉得又是一个高深统计模型。其实剥开来看核心逻辑非常朴素一个SNP如果能同时影响基因表达又影响性状那么基因表达有可能是这个SNP影响性状的中间环节。SMR想检验的就是“基因表达的改变是否介导了SNP到性状的关联”也就是把eQTL和GWAS这两层信号缝合起来找到那个处于中间位置的基因。1.1 一句话版本用“基因表达”解释GWAS信号用一个生活化的类比。假设你在路上看到一辆车撞了护栏GWAS信号车旁边站着一个人候选基因你很难判断这个人是不是司机。要确认这个人和事故的关系你需要看关键证据比如监控拍到这个人当时在驾驶位上eQTL证据。SMR做的事情就是系统性地审查“这个人是否在驾驶位”而不是只看他站在车旁边。更具体一点SMR利用了这样一个链条某个SNP位点z它和基因表达x显著相关这就是cis-eQTL信号同时这个SNP又在GWAS中与性状y显著相关。如果我们假设“z通过改变x的表达进而影响y”那么基因表达x对性状y的效应量就可以用SNP对性状的效应量去除以SNP对表达的效应量来近似估计。这个比值在因果推断里叫Wald ratio。SMR检验的核心就是这个比值是否显著不为零。如果显著说明基因表达的差异确实和性状的变化存在关联这个基因值得重点关注。这里最妙的一点是整个计算只需要eQTL和GWAS各自的汇总统计量不需要同时拥有同一批个体的基因型和表达量。两个独立的大样本研究一个提供了SNP对表达的效应一个提供了SNP对性状的效应SMR就足够把它们缝合起来。1.2 SMR检验与HEIDI检验一对搭配使用的“侦探搭档”只有SMR检验是不够的因为存在一个很常见的假阳性来源——连锁不平衡LD。假设真正影响基因表达的是SNP A真正影响性状的是SNP B而A和B因为离得近天然存在LD那么SNP A既会在eQTL中显著也会在GWAS中显著看起来就像是表达介导了性状其实是两个不同因果变异靠LD“串台”了。为了对付这种干扰SMR配套了一个HEIDI检验Heterogeneity in Dependent Instruments。HEIDI的思路也很直接如果基因表达和性状确实由同一个因果变异驱动那么不管拿位点附近的哪个SNP来算得到的“表达对性状的效应估计值”应该都差不多彼此之间没有异质性如果杀掉top SNP后附近其他SNP算出来的效应量明显不一致那就说明这个区域里可能有两个不同的因果变异在各自起作用SMR信号很可能只是LD造成的假象。所以记住一个反直觉的判定标准HEIDI的p值越“不显著”反而越好。我见过不止一次新手把HEIDI p值0.05当成了显著结果留下方向恰恰搞反了。p_HEIDI很小说明存在异质性这个结果要扔掉。1.3 为什么用SMR而不是直接做共定位有朋友会问做共定位不是也行吗确实coloc等共定位方法也能判断eQTL和GWAS信号是否共享因果变异而且不少人会把SMR和coloc结合使用。但两者侧重点不同。coloc是贝叶斯思路输出的是“两个信号共定位”的后验概率它不直接给你一个效应量。SMR则是在孟德尔随机化框架下给出基因表达对性状的因果效应估计值还带标准误和p值天然适合全基因组范围的大规模扫描。实际分析中我的习惯是先用SMR做全基因组层面的初步筛选找出“表达影响性状”的候选基因然后对最感兴趣的几个位点再用coloc做一次交叉验证。这样既发挥了SMR计算快、适合筛库的优势又用coloc规避了单一方法可能存在的模型设定误差。2. 开工前准备软件、参考数据与格式约定SMR分析本身跑起来很快但前期的数据准备往往比运行更耗时。软件安装很简单真正的门槛在于把GWAS数据整理成SMR要求的格式以及准备合适的eQTL和LD参考数据。这一节我把需要的东西按照清单列清楚。2.1 软件安装与执行环境SMR是纯命令行的Linux工具在服务器上使用。安装方式先到SMR的官方页面下载最新版本的源码包解压后直接make编译几分钟就能完成。编译环境只需要系统装了gcc和zlib开发库一般服务器都自带。编译完会生成一个smr可执行文件运行./smr --version能出现版本号就算装好了。SMR依赖的底层计算不算特别重内存需求取决于你分析的范围。跑单条染色体的时候2到4G内存足够如果做全基因组扫描建议内存给到8G以上。CPU方面支持多线程--thread-num参数可以指定线程数我在集群上一般给到10到20个线程速度提升很明显。这里有一个小建议SMR更新不算频繁但不同版本对参数名和输出格式有细微差别。我建议直接装1.3.1以上版本并且每次换版本后先用一个小规模数据跑一遍确认输出文件格式再投入使用免得分析做完才发现列名变了。2.2 参考数据从哪来eQTL、LD面板、GWASSMR分析需要三类数据缺一不可第一类是你的核心输入——GWAS summary数据。这个数据来自你自己正在研究的性状的公共数据库或者你自己跑出来的全基因组关联分析结果。注意必须是汇总统计量级别SMR不需要个体数据。第二类是eQTL summary数据。SMR官方贴心地提供了多种格式的现成数据最常用的是GTEx项目的eQTL数据直接以BESD二进制格式提供去官网的Data Download页面按需下载对应组织即可。此外还有eQTLGen血液组织样本量很大、ROSMAP脑组织、CAGE皮肤组织等可选。选哪个组织取决于你的研究性状比如研究血脂代谢就优先用肝脏、脂肪组织研究神经精神性状就优先用大脑各分区。第三类是LD参考面板。SMR在做HEIDI检验和SMR估计时需要知道位点附近的LD结构。官方推荐使用1000 Genomes Phase 3的参考面板PLINK格式的bed/bim/fam文件。面板的人群要和你的GWAS人群匹配欧洲人群的GWAS就选EUR面板东亚人群就选EAS面板跨人群混用会导致LD估计偏差直接让HEIDI检验失真。2.3 数据的三个前提条件在动手准备之前先检查一遍你的GWAS数据是否满足三个基本条件不满足的话后续分析会非常别扭首先GWAS数据中必须包含效应量beta如果原始数据提供的是OR值需要先取对数转成log(OR)SMR格式里列名可以叫b但本质要求线性刻度下的效应量。其次必须包含样本量NSMR算标准误需要用到。最后SNP的rs号、等位基因A1/A2必须和参考面板一致至少需要是一个你能对得上的版本比如都基于dbSNP build 150或类似版本避免rs号新旧版本对不上。满足这三个条件后后面的事情就顺多了。3. 数据清洗与格式转换格式转换是整个流程里最容易出问题的一步也是我每次帮别人排查时最常见的报错来源。SMR对输入格式的要求其实不复杂但它死板得很字段名不认识就跳过等位基因方向不对就全部对不上。这一节我按实际操作顺序拆开讲。3.1 标准格式字段说明先看SMR标准GWAS文本格式通常以.ma后缀名结尾核心字段包括SNPrs号比如rs629301A1效应等位基因effect alleleA2另一个等位基因freqA1的频率bA1对应的效应量注意是beta不是ORse效应量的标准误pGWAS的p值n样本量如果你下载的公共GWAS数据列名不一样比如用的是effect_allele、other_allele、beta、p_value这种不要紧SMR官方提供了一个gwas2smr的转换脚本通过配置文件映射原始列名和目标列名。设好原始列位置它自动输出符合SMR规范的格式。脚本对新用户更友好我基本都用它。这里有一个必须强调的细节如果原始数据里没有直接给beta而是给的OR比如很多二分类性状的GWAS一定要先对OR取自然对数。我见过有同学直接把OR填进b列导致后面所有效应量的符号和大小都不对HEIDI检验也乱成一团。3.2 用脚本把原始GWAS转成SMR格式我自己常用的做法是先用R或者Python做清洗因为公共GWAS数据往往有大量冗余行和不规范字段。以R为例核心流程大致是读入原始数据选出需要的列检查rs号是否存在重复去掉重复SNP然后处理等位基因链方向问题。等位基因的方向是这步的重点。对于AT和CG这两种类型如果参考面板无法确定正负链分析时容易出问题稳妥的做法是直接把这类多态性高的ambiguous SNP剔除。频率接近0.5的AT/CG SNP无法通过频率判断是否翻转链保留下来会污染结果。我会在代码里先过滤掉MAF介于0.4到0.6之间的AT和CG位点。清洗完之后用gwas2smr脚本或者直接手动写表头输出成SMR格式。输出之前再检查一遍BETA列没有缺失SE全部大于0P值在0到1之间N列有值。这些基础检查每次都能揪出一批原始数据里的格式混乱。3.3 eQTL数据两种使用路线eQTL数据的准备有两条路。最省心的是直接从SMR官网下载已经转好的BESD格式文件下载后解压就能用文件是一套以.besd结尾加上配套索引文件的形式。SMR官方提供了GTEx v8所有组织的版本我强烈建议优先走这条路因为自己从GTEx原始数据转BESD非常繁琐要处理样本级表达量、协变量、基因注释很多新手耗在这里好几天。另一条路是如果你有自己的eQTL数据比如自己的RNA-seq和基因型数据跑出来的结果那需要整理成文本格式。格式要求列名包括SNP、Chr、BP、A1、A2、Freq、Probe、Gene、b、se、p、n其中Probe是探针IDGene是基因名。整理好后可以直接用--eqtl-summary参数读取文本格式也可以转成BESD再读取。文本格式路径适合数据量小的场景。如果eQTL数据非常大比如上千万个探针位点建议还是转成BESD因为BESD是按块压缩的SMR运行起来内存开销小、速度快。转BESD工具也集成在SMR软件包里--make-besd参数就能生成。4. 完整SMR运行流程数据和参考都齐了接下来就是最核心的操作环节。我会用一次完整的分析示例说明假设场景是我有一个欧洲人群的代谢性状GWAS想用肝脏组织的GTEx v8 eQTL数据做SMR找到影响该性状的候选基因。4.1 命令行基本操作整个SMR全基因组扫描的命令行其实非常简洁下面是我常用的一份./smr \ --bfile /ref/1000G_EUR \ --gwas-summary /data/trait.ma \ --beqtl-summary /ref/GTEx_v8_Liver \ --out /output/liver_smr_result \ --thread-num 20 \ --maf 0.01参数含义依次是--bfile指定LD参考面板的PLINK文件前缀--gwas-summary指定GWAS汇总数据--beqtl-summary指定BESD格式的eQTL数据前缀--out指定输出文件前缀--thread-num指定线程数--maf过滤低频SNP。这一步跑完会生成两个核心输出文件文件名是liver_smr_result.sml和liver_smr_result.sns这两个文件的差异和读取方式我放到第五节详细说。4.2 常用参数选择与计算取舍SMR常用参数里有几个需要动脑子思考而不是直接默认跑完不管。第一个是--cis-window。这个参数定义了多少距离内的eQTL SNP被认为是cis-eQTL默认是2000kb也就是上下游各2Mb。这个区间范围基本覆盖了大多数基因的调控区域。如果分析的是特定基因组区域可以调小一点全基因组扫描时保持默认就好。第二个是--peqtl和--pgwas。这两个参数是筛选阈值分别对应eQTL和GWAS位点的p值上限。默认值通常是5e-8但实际操作中可以放宽一点比如eQTL信号可以放宽到1e-4因为有些基因的eQTL效应本身比较温和。我一般会把--peqtl设成1e-4--pgwas设成5e-8在保证信号质量的同时尽量多纳入一些候选基因。第三个是--probe。如果你只关心某几个特定基因可以用这个参数指定探针列表SMR只分析这些基因速度飞快。全基因组扫描则不需要这个参数。第四个是--extract-snp。如果只想分析某个染色体区域用这个参数提取SNP范围。我在做“已知GWAS信号附近找基因”这种任务时通常先提取GWAS显著位点前后几Mb的SNP再配合指定的探针列表这样跑出来更精准内存和耗时都低得多。4.3 一次全基因组扫描的完整示例假设前面那份命令已经跑起来实际运行中控制台会不断打印当前在处理的探针显示探针ID、基因名、染色体位置和分析进度。整个全基因组扫描的速度相当快GTEx肝脏组织配1000G EUR参考20线程跑下来大概十几分钟到半小时具体取决于探针数量。跑完后我会立即做一件事检查生成的文件大小和内容头几行确认没有出现“all SNPs filtered out”这样的警告。出现这个警告通常说明某个文件里的SNP ID格式和参考面板不一致后面会专门讲这条坑。全基因组扫描输出的.sns文件里包含了所有探针的分析结果不管有没有达到显著阈值都记录在内。.sml则是经过筛选后的显著位点列表SMR默认会筛选出一个较宽松的候选集。我的习惯是先打开.sns文件看整体分布因为它能反映一个客观情况——到底有多少基因的表达和性状存在显著关联。如果一份结果里.sns有几百个基因达到p_SMR1e-4那这批数据质量很高后续筛选空间大如果只有零星几个那就要考虑是不是GWAS本身功效不够或者eQTL组织不合适。5. 结果解读与候选基因筛选命令跑完只是第一步怎么解读输出才是真正考验功力的时候。很多人拿到.sml文件看到p值很小就直接写进论文忽略HEIDI方向后来又推翻重来。这一节我讲讲完整的解读口径。5.1 输出文件怎么读.sml和.sns文件的列结构基本一致只是筛选条件不同。核心列包括Probe探针ID一般是基因表达的探针标识Gene基因名Chr、BP基因位置TopSNP该位点最显著的SNP rs号A1、A2、FreqTopSNP的等位基因信息和频率bSMR估计的基因表达对性状的效应量也就是基因表达每增加一个单位性状log-odds或标准化表型的变化seb的标准误pSMR检验的p值p_HEIDI异质性检验的p值nsnp_heidiHEIDI检验用的SNP数量我习惯在读取.sml后立刻用筛选条件把合格的候选基因挑出来。核心条件一般是两条p_SMR足够小比如小于经过Bonferroni校正后的阈值p_HEIDI足够大比如大于0.05或者保守一点的0.01。5.2 筛选标准单点P值HEIDI方向关于p_SMR的阈值很多人问到底多少算显著。这要看你做了多少次检验。如果全基因组扫描了约两万个基因Bonferroni校正后的阈值就是0.05除以分析的基因数大约在2.5e-6左右。如果只是分析特定基因组区域里的几百个基因阈值就可以放宽到0.05除以探针数约1e-4到1e-3量级。实际操作中我建议先用一个较宽松的阈值看结果比如p_SMR1e-4然后配合p_HEIDI0.05筛。如果初筛结果太多再收紧到p_SMR2.5e-6。最理想的状态是某个基因既满足严格的p_SMR又满足p_HEIDI0.05同时eQTL和GWAS信号都指向同一个top SNP那这个基因的可信度就非常高了。这里再提醒一次HEIDI方向p_HEIDI越大越好它表示“没有证据表明存在异质性”也就是支持单一因果变异驱动。p_HEIDI小于0.05说明结果要打问号可能是因为LD里藏着两个不同变异。不要把它当成普通的显著性p值去理解。5.3 可视化解读locus plotSMR自带了locus plot绘制功能只要在运行命令中加上--plot --plot-y-max 20这类参数就会在输出目录生成PDF文件每个显著的探针对应一张图。locus plot的左上角一般是eQTL的关联信号右上角是GWAS的关联信号中间一栏展示SMR分析估计的b_xy在所有SNP之间是否一致。这张图能快速判断一个信号质量如何如果各SNP的b_xy值差别很大呈现明显的发散状即使p_HEIDI勉强大于0.05我也会对这个结果保持警惕。如果b_xy的值基本聚集在同一个点附近那就比较可靠。我一般会先批量生成所有显著位点的locus plot用眼睛快速扫一遍重点看那些“图形漂亮”的位点——GWAS和eQTL信号峰形一致b_xy分布集中这种位点写进文章里底气足。5.4 定位到候选基因后的交叉验证SMR输出的候选基因本质上还是统计推断结果不能直接当成生物学结论。我通常会做几件事来交叉验证第一是检查该基因在多个eQTL组织中是否都显著。比如某基因在肝脏组织中SMR显著在脂肪组织中也显著我会更相信它不是组织特异性的偶然信号。第二是用coloc做一次共定位。SMR和coloc虽然机制不同但两者结论交叉一致时证据强度会高很多。SMR给出的是方向和效应量coloc给出的是共定位概率结合起来可以说是互为补充。第三是查一下该基因的已知生物学功能。就算统计上关联很强如果基因功能与性状毫无关系比如一个调控色素合成的基因关联上了血压那大概率是某种未被识别的混杂或者LD效应。生物学合理性永远是最后一道关卡。6. 常见问题与实战坑点这个章节是写给每一个准备动手跑SMR的人看的。我在实际分析中遇到过不少问题有的是数据格式问题有的是对方法理解的偏差这里挑最典型的几类列出来按出现频率排序。6.1 运行阶段常见报错第一个高频报错是“SNP IDs in your GWAS are not found in the reference panel”。遇到这个90%的原因是GWAS数据的rs号或者染色体编号格式和参考面板不一致。比如GWAS里的染色体是chr1格式而参考面板是1格式或者rs号带着版本后缀。解决方法是统一格式通常做一次字符串替换去掉chr前缀就行。第二个高频问题是等位基因方向不一致。GWAS报告里的A1是效应等位基因但不同数据库的allele编码可能基于正链、负链或者不同参考基因组版本。SMR会对等位基因做核对发现不一致时会自动翻转或剔除。如果你的GWAS和eQTL的SNP重叠率很差先检查是不是链方向的反了。遇到AT、CG这类ambiguous SNP前面说过直接过滤掉最省心。第三个问题是运行到一半崩溃报内存不足。这通常发生在LD矩阵计算环节。可以尝试减少--thread-num或者用--ld-ctrl-gb限制LD矩阵的内存控制参数如果还不行就把分析范围切到单条染色体分开跑。第四个问题是运行很快结束但.sns文件为空。这说明你的GWAS数据一个位点都没能匹配上参考面板。检查一下--gwas-summary的列名是不是符合SMR要求特别是SNP列名是否为SNPA1列是否为A1这些标准命名文件名大小写、分隔符是不是都正确。6.2 结果解读阶段的坑第一个结果层面的坑是我反复强调的HEIDI方向问题。P_HEIDI0.05的结果要谨慎处理意味着SMR信号可能是LD造成的。很多新手第一次跑出全基因组几百个显著基因兴高采烈准备写文章结果发现里面一半都是p_HEIDI0.05的假阳性全部得剔除。第二个坑是忽略了效应方向。SMR分析中b_xy的正负号代表基因表达对性状的风险或保护方向。如果一个基因的b_xy为正说明表达量越高性状值越大如果为负则反之。很多人只看p值不看b的符号结果后续在机制解释时发现方向和文献相悖不得不返工。第三个坑是样本重叠问题。SMR方法推导时假设GWAS和eQTL的样本是独立的。如果你的GWAS和eQTL数据来自同一队列样本重叠会导致检验统计量膨胀p值偏小假阳性增加。公共数据的GWAS和GTEx一般不存在这个问题但自己分析时要注意。第四个坑是用错的eQTL组织。这个在选题阶段就要想清楚。研究阿尔茨海默病用肝脏的eQTL结果当然大概率一片阴性或者有的没的信号。正确的做法是先确定靶组织再下载相应组织的BESD数据。6.3 我的几条使用习惯最后分享几个我自己总结的小技巧不算复杂但能避开不少麻烦。第一个习惯是拿到一个不熟悉的GWAS数据后先不急着全基因组跑而是先挑选几个已知的、经过文献验证的阳性位点做一次小范围测试确认SMR能检测到这些已知信号。如果已知位点一个都检测不到说明数据格式或者参考面板有问题再直接去全基因组扫描大概率也做不出好结果。第二个习惯是同一个分析我会同时保留宽松和严格的筛选结果。宽松结果给自己看用来评估整体信号格局严格结果给论文用用来列候选基因表。两者差异本身也很有价值能反映结果对不同阈值的稳定性。第三个习惯是跑完SMR后把显著的基因列表和区域内的已知基因做一次交叉比对看看有没有之前GWAS精细定位或者TWAS分析里反复出现的基因。多个方法互相印证的基因值得优先投入后续验证。最后的体会SMR这套流程我用了好几年最大的感受是它的“性价比”很高——只要你把数据格式整理干净运行本身非常快结果又很直观。但我也越来越清楚它给出的候选基因只是起点不是终点。统计上的因果关联不等于生物学上的因果机制最终还需要分子实验去验证基因表达变化对表型的影响。做分析的时候保持“一切统计信号都可能是假的”这种心态反而能避免不少弯路。如果后续你想把这个流程扩展到更深的层次可以考虑换组织、换eQTL来源比如eQTLGen、pQTL、mQTL也可以把SMR和TWAS、coloc的结果放在一起做多方法整合。数据和代码都不复杂难的是每到一个步骤都有清晰的判断标准。希望这篇指南能帮你少踩几个坑顺利跑出能够经得起推敲的结果。