2026/9/19 7:58:06
FEATURED · 精选文章 
很多做GWAS的朋友应该都有过这种体验辛辛苦苦跑完关联分析拿到一堆显著位点结果注释到基因头上却发现大部分位点都在基因间区或者内含子里根本说不清楚到底通过什么机制影响表型。这时候如果手里刚好有一套eQTL数据就能往前推一步——看看这个位点是不是通过调控某个基因的表达量间接影响了你的性状。而SMRSummary-based Mendelian Randomization这个工具就是专门干这件事的它把GWAS数据和eQTL数据整合到一起用孟德尔随机化的思路判断基因表达与性状之间的因果关系是目前做共定位和因果推断最常用的手段之一。这份指南我会直接从实际应用角度出发完整走一遍SMR的分析流程包括方法原理、数据准备、命令行操作、结果解读和常见坑。无论你是刚接触生信的学生还是已经跑过几轮GWAS但没碰过SMR的研究人员照着我这个流程走基本能跑通整个分析。1. SMR方法的核心逻辑与适用场景1.1 为什么GWAS显著位点需要eQTL分析来补充GWAS的结果本质上只告诉你在某个位点上不同基因型的人群表型有显著差异。但大多数显著的SNP位点并不是直接改变蛋白质氨基酸序列的错义突变而是落在内含子、基因间区这些非编码区域。这种位点要解释生物学机制最常见的一条路就是看它是否影响基因表达——即这个位点是不是一个eQTLexpression Quantitative Trait Locus表达数量性状位点。如果SNP的基因型能显著影响某个基因的表达量而这个基因表达量的高低又和疾病风险相关那么这条“SNP→基因表达→表型”的因果链条就串起来了。SMR正是用来检验这条链条的它把GWAS的关联信号和eQTL的关联信号放在一起判断两者是否来自同一个因果变异。我举个例子你就明白了。假设你GWAS发现rs123456这个位点与冠心病显著关联但这个位点位于某个基因上游20kb的增强子区域内。你查了GTEx数据库的eQTL数据发现rs123456是某个基因的cis-eQTL携带风险等位基因的人该基因表达量显著升高。这时候用SMR计算如果仪器变量工具变量检验通过你就可以合理推测这个位点是通过上调该基因表达来增加冠心病风险的。这就比单纯报告“发现了一个关联位点”深入了一大步。1.2 SMR与MR分析的异同以及HEIDI检验的作用SMR的统计学框架本质上是一个基于汇总统计数据的孟德尔随机化分析。经典的MR分析通常需要个体水平的基因型-表型数据或者至少需要完整的SNP-暴露和SNP-结局的汇总统计量。SMR做的是类似的事情把基因表达量当作暴露把GWAS性状当作结局然后利用cis-eQTL的SNP作为工具变量来推断因果关系。不过SMR做了一些简化——它不需要两个独立样本的完整数据只要有一个包含SNP效应量的GWAS汇总数据和一个eQTL汇总数据就够了。这里必须提一下HEIDI检验Heterogeneity in Dependent Instruments这是SMR非常核心的一步。SMR分析得出显著关联后存在一种可能其实因果变异并不是同一个只是两个不同的连锁变异分别驱动了eQTL信号和GWAS信号导致SMR出现假阳性。HEIDI检验就是用来排除这种连锁不平衡造成的假阳性。一句话总结就是SMR负责检测关联HEIDI负责排除假阳性。我实际跑的时候通常要求p_SMR小于0.05且p_HEIDI大于0.01才算通过——SMR显著但HEIDI也显著比如p_HEIDI小于0.01的位点我会果断放弃因为那很可能是LD造成的假信号。1.3 SMR适合回答哪些科学问题SMR不是万能的但它有一套非常明确的应用场景。最常见的三个用途一是给GWAS显著位点做功能注释把“风险位点”升级为“风险位点→风险基因→机制”二是跨组织比较比如用GTEx的多个组织eQTL数据分别跑SMR看某个基因的表达效应是否具有组织特异性这往往能提示疾病的关键组织三是筛选药物靶点通过SMR找到与疾病显著关联的基因后可以进一步评估这些基因作为药物靶点的可行性。需要特别说明的是SMR和另一种常用的共定位方法coloc有区别。coloc回答的问题是“两个关联信号是否共享同一个因果变异”它的输出是一个后验概率并不直接给出因果方向而SMR则直接给出基因表达对性状的因果效应估计值beta和方向。两者经常结合使用SMR筛一遍、coloc验证一遍结论会更有说服力。2. 数据准备与工具链选择2.1 GWAS汇总数据的格式转换与质量控制跑SMR之前最花时间的准备工作其实是数据格式整理。SMR软件官方要求的GWAS输入格式叫 .ma 格式本质上是一个文本文件每一行代表一个SNP。但不同来源的GWAS数据字段名五花八门有的叫OR、有的叫Beta、有的直接给P值不给效应量。所以第一步是把你的GWAS汇总数据统一转换成SMR能识别的格式。SMR的.ma格式必须包含以下几列SNP、Chr、BP、A1效应等位基因也就是被估效应量的那个等位基因、A2另一个等位基因、FreqA1的频率、BetaA1的效应量、SE标准误、PP值。如果你的原始数据是logistic回归得到的OR值和95%置信区间记得先用ln(OR)转换成BetaSE可以用log(CI上限)减去log(OR)再除以1.96得到。这一步很容易出错建议写个小脚本统一处理不要手动改行数一多必然出问题。质量控制方面我比较关注三个点第一检查等位基因是否与参考基因组版本匹配如果GWAS用的dbSNP版本和参考面板不一致位点的rs号可能对不上第二过滤掉无效SNP——比如allele不是ATCG的、频率在0或1的、P值缺失的第三最好顺便记录一下原始GWAS的总样本量因为SMR计算时样本量会影响到变异标准误的估计不同位点样本量差异较大的话建议按样本量分层跑。2.2 eQTL数据来源GTEx、eQTLGen与特定组织选择eQTL数据的质量决定了SMR结果的下限。目前最常用的来源是GTExGenotype-Tissue Expression数据库它有V8版本覆盖49个组织的cis-eQTL数据。GTEx的SMR格式数据可以直接从SMR官方提供的资源页下载对应文件为GTEx_All_V8_eQTL_SNP_bgen不需要自己再处理。这个文件是bgen格式需要通过SMR的--eqtl-summarize命令转换成二进制格式。值得注意的是GTEx的样本量在不同组织中差异很大比如肌肉骨骼组织样本量可能有700多而一些小组织只有100出头样本量小的组织eQTL检出效率低SMR的统计效能也会跟着打折。如果是研究血液相关性状eQTLGen数据库是个更好的选择。eQTLGen整合了31,684例外周血样本的cis-eQTL数据统计效能远高于GTEx的全血数据。我在跑血液指标相关GWAS时基本优先用eQTLGen。另外还有一些特定数据库可以参考比如BraINeek脑组织、Metabrain脑组织、DICE免疫细胞等按你的研究方向选用即可。不管用哪个数据库都建议下载后先检查一下版本和基因组坐标是否跟你的GWAS数据一致。如果GWAS用的是GRCh37坐标eQTL数据是GRCh38坐标不经过坐标转换就硬跑结果会非常难看——大量SNP因为坐标对不上而缺失或者更糟的是错配到错误的位点。2.3 软件环境与参考面板准备SMR的软件本体是统一的命令行工具在Linux服务器上用起来最顺手。Windows也有编译好的exe版本但处理大规模数据时容易内存不足我还是强烈建议在Linux环境下跑。安装步骤不复杂去SMR官方下载对应系统版本的压缩包解压后就能直接运行。需要注意SMR依赖两个外部工具一是PLINK用于计算LD和进行一些数据操作二是bgen工具包比如bgenix如果eQTL数据是bgen格式需要先用它来做格式转换或者提取特定区域的基因型数据。参考面板我用的是1000 Genomes Phase 3的VCF文件这是SMR官方教程里推荐的。参考面板主要用于估计SNP之间的LD结构SMR分析HEIDI检验时需要用到这个信息。如果你的数据集有明显的种族来源差异比如做东亚人群的GWAS那最好准备相应的东亚参考面板。我自己做过一个教训很深的项目手里是东亚人群的GWAS参考面板却用了欧洲人群的1000G数据结果HEIDI检验的p值系统性偏小一大批本来靠谱的位点被判成假阳性排查了很久才发现是参考面板搞错了。3. SMR完整实操流程与关键命令解析3.1 建立eQTL二进制格式数据集eQTL原始文件下载下来后第一次使用需要建立二进制索引格式。这里以GTEx全血数据为例# 从SMR官方资源页下载GTEx cis-eQTL数据bgen格式 wget https://yanglab.westlake.edu.cn/data/smr/GTEx_All_V8_eQTL_SNP_bgen/Whole_Blood.v8.EUR.zip unzip Whole_Blood.v8.EUR.zip # 将bgen格式转换为SMR二进制格式 ./bgenix -g Whole_Blood.v8.EUR.bgen -list Whole_Blood.v8.EUR.snps # 使用SMR的--eqtl-summarize命令生成bti索引 ./smr --eqtl-summarize \ --eqtl-bgen Whole_Blood.v8.EUR.bgen \ --eqtl-info Whole_Blood.v8.EUR.samples \ --eqtl-summary whole_blood_eur.txt \ --out whole_blood_eur这一步跑完之后会生成一个.epi后缀的二进制文件和一个.epg后缀的索引文件后续跑SMR时通过--beqtl-summarize参数调用。这里有个细节SMR官方提供的GTEx数据默认是经过MAF过滤的但不同组织过滤阈值略有不同建议在分析前用--eqtl-summary输出文件确认一下各组织的信息完整度。如果你的eQTL数据不是bgen格式而是普通的文本格式比如从eQTLGen官网下载的txt.gz处理方式会不一样。eQTLGen提供的汇总文件包含SNP、Gene、beta、se、p等字段你需要先根据SMR文档里的说明整理成规定的列名格式然后用--eqtl-summarize参数配合--eqtl-text来创建二进制文件。这里注意文本格式的数据必须包含allele频率这一列Freq否则SMR会报错甚至静默跳过大量SNP。3.2 GWAS数据格式检查与标准化有了eQTL的二进制文件之后就需要把自己的GWAS数据整理成.ma格式并标准化。这一步我一般用一个小脚本来完成避免手动修改。假设你的原始数据是一个csv文件包含rsid、chromosome、position、effect_allele、other_allele、eaf、beta、se、pvalue这些列标准的处理逻辑是import pandas as pd import numpy as np df pd.read_csv(my_gwas.csv) df df.rename(columns{ rsid: SNP, chromosome: Chr, position: BP, effect_allele: A1, other_allele: A2, eaf: Freq, beta: Beta, se: SE, pvalue: P }) # 确保A1是效应等位基因这里假设原始数据已经是 # 过滤掉无效的allele df df[df[A1].isin([A, T, C, G]) df[A2].isin([A, T, C, G])] # 过滤掉频率为0/1的SNP df df[(df[Freq] 0) (df[Freq] 1)] # 如果原始数据是OR值需要转换 # df[Beta] np.log(df[OR]) df[[SNP, Chr, BP, A1, A2, Freq, Beta, SE, P]].to_csv(my_gwas.ma, sep\t, indexFalse)这里要特别强调等位基因方向的问题。SMR计算时要求eQTL的效应等位基因和GWAS的效应等位基因能正确对齐如果两个数据集中同一个SNP报告的是相反链上的等位基因Beta会完全反向。常见的错误出现在A/T、C/G这种互补等位基因上因为无法简单判断是否对齐。虽然有工具可以做变异链翻转但实际操作中我更推荐在整理GWAS数据时用PLINK的--flip命令检查与参考面板的一致性确保大部分SNP的等位基因方向和参考面板一致。3.3 运行SMR核心分析数据都准备好后运行SMR就非常简单了。基本命令如下./smr \ --bfile /path/to/1000G_EUR \ --gwas-summary my_gwas.ma \ --beqtl-summarize whole_blood_eur \ --out my_smr_result \ --thread-num 8 \ --peqtl-smr 0.05 \ --peqtl-heidi 0.05 \ --maf 0.01 \ --diff-freq 0.2这里逐个解释参数的用意--bfile参数用来指定PLINK格式的参考面板必须是二进制的bed/bim/fam前缀。SMR会自动计算相关SNP之间的LD矩阵用于HEIDI检验。--gwas-summary指定整理好的GWAS汇总文件。--beqtl-summarize指定eQTL二进制文件前缀。--out指定输出文件前缀。--thread-num设置多线程数大文件建议至少给8否则速度感人。--peqtl-smr控制在基因层面的显著性阈值一般0.05即可如果想更严格也可以设成0.01。--peqtl-heidi控制HEIDI检验的显著性阈值0.05是官方推荐值。之前的实践中我往往还会把阈值设成0.01来进一步过滤具体看后续验证需求。--maf过滤掉低于1%的SNP。--diff-freq控制eQTL和GWAS中效应等位基因频率差的阈值如果某个SNP在两个数据集中的频率相差超过0.2会怀疑是数据错配默认是0.2通常不需要改。运行过程中SMR会在屏幕上打印当前处理到哪个染色体哪个基因区间进程速度取决于参考面板的LD计算量和你eQTL数据的规模。以GTEx全血数据加一个500万SNP的GWAS文件为例在32核服务器上大概跑30~50分钟如果是eQTLGen这种百万级eQTL的数据可能需要1~2小时。跑完后输出一堆文件最核心的是.smr结尾的结果文件以及.txt结尾的详细记录文件。3.4 输出文件解读与条件筛选SMR的结果文件是纯文本的每一行代表一个基因和某个SNP的检验结果。核心列包括Gene基因名比如ENSG编号建议后续映射回symbol名便于阅读。ProbeChr、ProbeBP探针所在位置。topSNP与基因表达关联最强的SNP的rs号。topSNP.Chr、topSNP.BP这个SNP的位置。A1、A2这个位点的等位基因。FreqA1频率。b_GWAS、se_GWAS、p_GWAS该SNP在GWAS中的效应量、标准误、P值。b_eQTL、se_eQTL、p_eQTL该SNP在eQTL数据中的效应量、标准误、P值。b_SMR、se_SMR、p_SMRSMR计算得到的因果效应估计和P值这是最核心的一列。p_HEIDIHEIDI检验的P值。nsnp_HEIDIHEIDI检验使用的SNP数量。拿到结果后我通常按以下标准筛选显著位点p_SMR 0.05 / 检验基因总数Bonferroni校正如果检验了几千个基因阈值会非常苛刻。我自己一般先用0.05做个初步筛选再对候选位点做多重检验校正避免漏掉真正有信号的基因。p_HEIDI 0.01表示没有异质性证据即GWAS和eQTL的信号很可能共享同一个因果变异。如果p_HEIDI小于0.01会再检查nsnp_HEIDI如果太少比如小于3结果不够可靠直接把该基因放弃。筛选出来的基因再检查一下效应方向是否一致——b_GWAS和b_eQTL方向一致时b_SMR也为正表示基因表达越高风险越大方向相反则为保护效应。这种方向一致性检查能避免一些怪异的错误结果。4. 结果可视化与生物学解释4.1 用曼哈顿图展示SMR信号SMR跑完只有一串数字表格可读性太差。我习惯用R来做可视化最常用的是类似GWAS曼哈顿图的SMR曼哈顿图把每个基因的p_SMR画出来横轴是基因组位置纵轴是-log10(p_SMR)。这里推荐用CMplot这个R包一行代码就能出图library(CMplot) smr_result - read.delim(my_smr_result.txt, header TRUE) # 确保有Chr和BP信息这里以topSNP的位置为例 CMplot(smr_result, plot.type m, col c(grey30, skyblue), LOG10 TRUE, threshold list(0.05/nrow(smr_result)), threshold.lty 2, threshold.lwd 1, threshold.col red, amplify TRUE, chr.den.col NULL, file pdf, file.name smr_manhattan.pdf, dpi 300, width 12, height 6)如果你的数据里基因位置信息比较稀疏也可以先按染色体分组统计每条染色体上的基因数量再用ggplot2画。但CMplot曼哈顿图的好处是它会自动标出超过阈值线的点旁边还会附上对应的基因名对后续挑选候选基因非常方便。我看曼哈顿图时会特别关注那些颜色最突出的点然后回到表格中逐个核对p_SMR和p_HEIDI。4.2 基因表达与表型因果效应的方向判断SMR分析的最终产出除了显著基因的列表还有一个关键信息是b_SMR的方向和大小。b_SMR的单位解释起来要小心——它表示的是基因表达量每变化一个标准差SD时对GWAS性状log(OR)或线性表型的影响。如果GWAS是二分类性状比如疾病b_SMR大于0表示基因表达升高会提升疾病风险小于0则相反。如果GWAS是连续性状比如血压b_SMR表示基因表达每增加1个SD血压平均变化多少。举一个实际案例我曾有个项目研究某炎症因子与冠心病的关联SMR的结果显示b_SMR为0.35p_SMR最小达到了1.2×10⁻⁶p_HEIDI为0.23这意味着基因表达水平每提高1个标准差冠心病的log(OR)增加0.35大约对应OR为1.42。这个效应量在孟德尔随机化研究中算是比较强的信号后续也通过coloc验证两信号共定位的后验概率高达0.94进一步增强了结论的可信度。解读的时候注意一个坑eQTL数据里的beta是基因表达量标准化后的效应不同组织、不同数据库的标准差定义不同因此b_SMR不能跨数据库直接比较绝对值大小。我做跨组织比较时主要看方向一致性和p值显著性不看beta数值大小。4.3 与coloc结果联动验证SMR的结果虽然方便但审稿人往往会要求额外的共定位验证。我通常的做法是SMR筛出来的显著基因全部跑一遍coloc的共定位分析。SMR和coloc的关系可以这么理解SMR告诉你有因果关联coloc告诉你两个信号是否共享同一个因果变异。如果coloc的PP.H4两个信号共定位的后验概率大于0.8那这个位点的可信度就非常高。coloc分析的输入也是GWAS汇总数据和eQTL汇总数据但需要先提取特定基因区域内的SNP。实际操作时我会从SMR结果中挑出显著基因以topSNP为中心提取上下游500kb范围内的SNP效应量整理成coloc要求的格式然后用R的coloc.abf函数逐个分析。如果某个基因在SMR中显著但coloc显示PP.H3两个信号独立的后验概率很高基本可以确认是LD导致的假阳性我一般会把它从最终列表里剔除。5. 常见问题与排查技巧实录5.1 eQTL数据与GWAS数据坐标不一致这是最让我头疼的问题也是用户在高频提问中出现最多的情况。只要GWAS和eQTL数据来自不同的参考基因组版本GRCh37 vs GRCh38染色体坐标就会整体偏移。SMR官方虽然提供了部分数据的坐标转换说明但实际运行起来如果坐标不一致通常表现为结果文件中大量基因的nsnp_HEIDI为0或者topSNP根本没有出现在eQTL数据里。排查方法也很直接随机挑几个已知的位点检查它们在两个数据中的染色体位置是否一致。如果发现不一致建议先用Liftover工具把GWAS结果转换成eQTL数据的坐标版本或者反过来转换。这里必须提醒一下Liftover转换完后要再做一次等位基因方向校验因为坐标转换过程中有可能出现链方向变化导致A/T或者C/G互换。5.2 内存不足与运行时间过长跑SMR的时候内存不足也是一个高频问题。eQTLGen的bgen文件有几十个GBGTEx虽然小一些但如果一次跑全部组织内存消耗会快速飙升。如果服务器内存不够比如16GB以下建议分染色体运行SMR最后再合并结果。SMR命令中没有专门的单染色体参数但可以通过修改输入数据的方式实现——比如用PLINK按染色体拆分参考面板同时准备对应染色体上的GWAS子集。这样并行跑22条染色体最后用cat命令合并结果文件。我遇到过最极端的情况是eQTLGen全量数据在32G内存的机器上跑到一半直接报Killed错误。后面改用三号染色体中心区域单独跑三分钟就完成了。所以遇到内存问题优先考虑分块而不是升级硬件。5.3 HEIDI检验p值几乎都为0如果发现所有结果的p_HEIDI都接近0基本可以确定是参考面板出了问题。常见原因有参考面板的人群和你GWAS数据的种族不一致或者参考面板的LD结构与eQTL和GWAS样本不匹配。我之前做东亚人群数据时用欧洲参考面板跑HEIDI就是这个结果换回东亚参考面板后p_HEIDI恢复正常。另外如果GWAS样本量和eQTL样本量有著名的重叠比如同一个队列同时贡献了GWAS数据和eQTL数据也会导致HEIDI检验p值系统性偏小这种情况建议在论文里说明局限性或者尝试用独立数据集做验证。5.4 结果文件中大量基因缺失如果你发现输出的.smr结果文件里基因数量远少于预期多半是过滤阈值设置得太严。特别是--peqtl-smr如果设置成0.001很多基因连被检验的资格都没有。还有一种可能是你的GWAS数据格式里SNP命名方式和eQTL数据不一致——比如GWAS用的是chr:pos:ref:alt格式而eQTL是rs号SMR内部做ID匹配时如果对不上直接跳过。这种问题排查方法很简单取一个eQTL数据的top SNP在GWAS数据里搜索一下看看能不能找到、格式是否一致很快就能定位问题。5.5 多重检验校正与结果报告规范最后说一个报告层面的问题。SMR一次会检验成千上万个基因如果不做多重检验校正假阳性会爆炸。常用做法是FDR比如q值小于0.05或Bonferroni0.05除以检验基因数。我自己的习惯是用Bonferroni做严格筛选然后在补充材料里报告FDR结果给审稿人留出宽松的空间。论文方法部分需要写清楚使用的eQTL数据版本、样本量、过滤标准、参考面板版本以及SMR软件的版本号和关键参数这些细节直接决定结果的可重复性。写在最后的一点心得SMR这套流程我前前后后跑了不下几十遍。最开始只把它当个黑盒工具输两个文件出来一张表格后来踩的坑多了才逐渐明白每一步参数背后都是统计学假设和生物学常识的博弈。如果你现在正准备跑自己第一份SMR分析我的建议是不要急着上全基因组的大数据先挑一个明确的候选基因区域跑通整个流程看清楚每一步输出的中间文件和结果字段再扩展到全基因组范围。这样至少你能在半小时内完成一次完整的分析迭代而不是在等了两小时任务跑完之后才发现第一步的数据格式就错了。另外一个小技巧SMR和coloc是可以互相印证的但两者的前提假设并不完全一样建议在正式分析之前用自己领域内一篇已发表的高质量论文里的某个位点做一次“标定”——用同一套数据看看你能不能复现出对方的SMR信号。能复现出来你的数据整理和参数设置基本就没有问题了。技术细节说再多都不如自己亲手跑通一遍来得踏实。希望这份指南能帮你少走一些弯路。
SMR分析实战：从GWAS到eQTL的因果推断全流程指南
很多做GWAS的朋友应该都有过这种体验辛辛苦苦跑完关联分析拿到一堆显著位点结果注释到基因头上却发现大部分位点都在基因间区或者内含子里根本说不清楚到底通过什么机制影响表型。这时候如果手里刚好有一套eQTL数据就能往前推一步——看看这个位点是不是通过调控某个基因的表达量间接影响了你的性状。而SMRSummary-based Mendelian Randomization这个工具就是专门干这件事的它把GWAS数据和eQTL数据整合到一起用孟德尔随机化的思路判断基因表达与性状之间的因果关系是目前做共定位和因果推断最常用的手段之一。这份指南我会直接从实际应用角度出发完整走一遍SMR的分析流程包括方法原理、数据准备、命令行操作、结果解读和常见坑。无论你是刚接触生信的学生还是已经跑过几轮GWAS但没碰过SMR的研究人员照着我这个流程走基本能跑通整个分析。1. SMR方法的核心逻辑与适用场景1.1 为什么GWAS显著位点需要eQTL分析来补充GWAS的结果本质上只告诉你在某个位点上不同基因型的人群表型有显著差异。但大多数显著的SNP位点并不是直接改变蛋白质氨基酸序列的错义突变而是落在内含子、基因间区这些非编码区域。这种位点要解释生物学机制最常见的一条路就是看它是否影响基因表达——即这个位点是不是一个eQTLexpression Quantitative Trait Locus表达数量性状位点。如果SNP的基因型能显著影响某个基因的表达量而这个基因表达量的高低又和疾病风险相关那么这条“SNP→基因表达→表型”的因果链条就串起来了。SMR正是用来检验这条链条的它把GWAS的关联信号和eQTL的关联信号放在一起判断两者是否来自同一个因果变异。我举个例子你就明白了。假设你GWAS发现rs123456这个位点与冠心病显著关联但这个位点位于某个基因上游20kb的增强子区域内。你查了GTEx数据库的eQTL数据发现rs123456是某个基因的cis-eQTL携带风险等位基因的人该基因表达量显著升高。这时候用SMR计算如果仪器变量工具变量检验通过你就可以合理推测这个位点是通过上调该基因表达来增加冠心病风险的。这就比单纯报告“发现了一个关联位点”深入了一大步。1.2 SMR与MR分析的异同以及HEIDI检验的作用SMR的统计学框架本质上是一个基于汇总统计数据的孟德尔随机化分析。经典的MR分析通常需要个体水平的基因型-表型数据或者至少需要完整的SNP-暴露和SNP-结局的汇总统计量。SMR做的是类似的事情把基因表达量当作暴露把GWAS性状当作结局然后利用cis-eQTL的SNP作为工具变量来推断因果关系。不过SMR做了一些简化——它不需要两个独立样本的完整数据只要有一个包含SNP效应量的GWAS汇总数据和一个eQTL汇总数据就够了。这里必须提一下HEIDI检验Heterogeneity in Dependent Instruments这是SMR非常核心的一步。SMR分析得出显著关联后存在一种可能其实因果变异并不是同一个只是两个不同的连锁变异分别驱动了eQTL信号和GWAS信号导致SMR出现假阳性。HEIDI检验就是用来排除这种连锁不平衡造成的假阳性。一句话总结就是SMR负责检测关联HEIDI负责排除假阳性。我实际跑的时候通常要求p_SMR小于0.05且p_HEIDI大于0.01才算通过——SMR显著但HEIDI也显著比如p_HEIDI小于0.01的位点我会果断放弃因为那很可能是LD造成的假信号。1.3 SMR适合回答哪些科学问题SMR不是万能的但它有一套非常明确的应用场景。最常见的三个用途一是给GWAS显著位点做功能注释把“风险位点”升级为“风险位点→风险基因→机制”二是跨组织比较比如用GTEx的多个组织eQTL数据分别跑SMR看某个基因的表达效应是否具有组织特异性这往往能提示疾病的关键组织三是筛选药物靶点通过SMR找到与疾病显著关联的基因后可以进一步评估这些基因作为药物靶点的可行性。需要特别说明的是SMR和另一种常用的共定位方法coloc有区别。coloc回答的问题是“两个关联信号是否共享同一个因果变异”它的输出是一个后验概率并不直接给出因果方向而SMR则直接给出基因表达对性状的因果效应估计值beta和方向。两者经常结合使用SMR筛一遍、coloc验证一遍结论会更有说服力。2. 数据准备与工具链选择2.1 GWAS汇总数据的格式转换与质量控制跑SMR之前最花时间的准备工作其实是数据格式整理。SMR软件官方要求的GWAS输入格式叫 .ma 格式本质上是一个文本文件每一行代表一个SNP。但不同来源的GWAS数据字段名五花八门有的叫OR、有的叫Beta、有的直接给P值不给效应量。所以第一步是把你的GWAS汇总数据统一转换成SMR能识别的格式。SMR的.ma格式必须包含以下几列SNP、Chr、BP、A1效应等位基因也就是被估效应量的那个等位基因、A2另一个等位基因、FreqA1的频率、BetaA1的效应量、SE标准误、PP值。如果你的原始数据是logistic回归得到的OR值和95%置信区间记得先用ln(OR)转换成BetaSE可以用log(CI上限)减去log(OR)再除以1.96得到。这一步很容易出错建议写个小脚本统一处理不要手动改行数一多必然出问题。质量控制方面我比较关注三个点第一检查等位基因是否与参考基因组版本匹配如果GWAS用的dbSNP版本和参考面板不一致位点的rs号可能对不上第二过滤掉无效SNP——比如allele不是ATCG的、频率在0或1的、P值缺失的第三最好顺便记录一下原始GWAS的总样本量因为SMR计算时样本量会影响到变异标准误的估计不同位点样本量差异较大的话建议按样本量分层跑。2.2 eQTL数据来源GTEx、eQTLGen与特定组织选择eQTL数据的质量决定了SMR结果的下限。目前最常用的来源是GTExGenotype-Tissue Expression数据库它有V8版本覆盖49个组织的cis-eQTL数据。GTEx的SMR格式数据可以直接从SMR官方提供的资源页下载对应文件为GTEx_All_V8_eQTL_SNP_bgen不需要自己再处理。这个文件是bgen格式需要通过SMR的--eqtl-summarize命令转换成二进制格式。值得注意的是GTEx的样本量在不同组织中差异很大比如肌肉骨骼组织样本量可能有700多而一些小组织只有100出头样本量小的组织eQTL检出效率低SMR的统计效能也会跟着打折。如果是研究血液相关性状eQTLGen数据库是个更好的选择。eQTLGen整合了31,684例外周血样本的cis-eQTL数据统计效能远高于GTEx的全血数据。我在跑血液指标相关GWAS时基本优先用eQTLGen。另外还有一些特定数据库可以参考比如BraINeek脑组织、Metabrain脑组织、DICE免疫细胞等按你的研究方向选用即可。不管用哪个数据库都建议下载后先检查一下版本和基因组坐标是否跟你的GWAS数据一致。如果GWAS用的是GRCh37坐标eQTL数据是GRCh38坐标不经过坐标转换就硬跑结果会非常难看——大量SNP因为坐标对不上而缺失或者更糟的是错配到错误的位点。2.3 软件环境与参考面板准备SMR的软件本体是统一的命令行工具在Linux服务器上用起来最顺手。Windows也有编译好的exe版本但处理大规模数据时容易内存不足我还是强烈建议在Linux环境下跑。安装步骤不复杂去SMR官方下载对应系统版本的压缩包解压后就能直接运行。需要注意SMR依赖两个外部工具一是PLINK用于计算LD和进行一些数据操作二是bgen工具包比如bgenix如果eQTL数据是bgen格式需要先用它来做格式转换或者提取特定区域的基因型数据。参考面板我用的是1000 Genomes Phase 3的VCF文件这是SMR官方教程里推荐的。参考面板主要用于估计SNP之间的LD结构SMR分析HEIDI检验时需要用到这个信息。如果你的数据集有明显的种族来源差异比如做东亚人群的GWAS那最好准备相应的东亚参考面板。我自己做过一个教训很深的项目手里是东亚人群的GWAS参考面板却用了欧洲人群的1000G数据结果HEIDI检验的p值系统性偏小一大批本来靠谱的位点被判成假阳性排查了很久才发现是参考面板搞错了。3. SMR完整实操流程与关键命令解析3.1 建立eQTL二进制格式数据集eQTL原始文件下载下来后第一次使用需要建立二进制索引格式。这里以GTEx全血数据为例# 从SMR官方资源页下载GTEx cis-eQTL数据bgen格式 wget https://yanglab.westlake.edu.cn/data/smr/GTEx_All_V8_eQTL_SNP_bgen/Whole_Blood.v8.EUR.zip unzip Whole_Blood.v8.EUR.zip # 将bgen格式转换为SMR二进制格式 ./bgenix -g Whole_Blood.v8.EUR.bgen -list Whole_Blood.v8.EUR.snps # 使用SMR的--eqtl-summarize命令生成bti索引 ./smr --eqtl-summarize \ --eqtl-bgen Whole_Blood.v8.EUR.bgen \ --eqtl-info Whole_Blood.v8.EUR.samples \ --eqtl-summary whole_blood_eur.txt \ --out whole_blood_eur这一步跑完之后会生成一个.epi后缀的二进制文件和一个.epg后缀的索引文件后续跑SMR时通过--beqtl-summarize参数调用。这里有个细节SMR官方提供的GTEx数据默认是经过MAF过滤的但不同组织过滤阈值略有不同建议在分析前用--eqtl-summary输出文件确认一下各组织的信息完整度。如果你的eQTL数据不是bgen格式而是普通的文本格式比如从eQTLGen官网下载的txt.gz处理方式会不一样。eQTLGen提供的汇总文件包含SNP、Gene、beta、se、p等字段你需要先根据SMR文档里的说明整理成规定的列名格式然后用--eqtl-summarize参数配合--eqtl-text来创建二进制文件。这里注意文本格式的数据必须包含allele频率这一列Freq否则SMR会报错甚至静默跳过大量SNP。3.2 GWAS数据格式检查与标准化有了eQTL的二进制文件之后就需要把自己的GWAS数据整理成.ma格式并标准化。这一步我一般用一个小脚本来完成避免手动修改。假设你的原始数据是一个csv文件包含rsid、chromosome、position、effect_allele、other_allele、eaf、beta、se、pvalue这些列标准的处理逻辑是import pandas as pd import numpy as np df pd.read_csv(my_gwas.csv) df df.rename(columns{ rsid: SNP, chromosome: Chr, position: BP, effect_allele: A1, other_allele: A2, eaf: Freq, beta: Beta, se: SE, pvalue: P }) # 确保A1是效应等位基因这里假设原始数据已经是 # 过滤掉无效的allele df df[df[A1].isin([A, T, C, G]) df[A2].isin([A, T, C, G])] # 过滤掉频率为0/1的SNP df df[(df[Freq] 0) (df[Freq] 1)] # 如果原始数据是OR值需要转换 # df[Beta] np.log(df[OR]) df[[SNP, Chr, BP, A1, A2, Freq, Beta, SE, P]].to_csv(my_gwas.ma, sep\t, indexFalse)这里要特别强调等位基因方向的问题。SMR计算时要求eQTL的效应等位基因和GWAS的效应等位基因能正确对齐如果两个数据集中同一个SNP报告的是相反链上的等位基因Beta会完全反向。常见的错误出现在A/T、C/G这种互补等位基因上因为无法简单判断是否对齐。虽然有工具可以做变异链翻转但实际操作中我更推荐在整理GWAS数据时用PLINK的--flip命令检查与参考面板的一致性确保大部分SNP的等位基因方向和参考面板一致。3.3 运行SMR核心分析数据都准备好后运行SMR就非常简单了。基本命令如下./smr \ --bfile /path/to/1000G_EUR \ --gwas-summary my_gwas.ma \ --beqtl-summarize whole_blood_eur \ --out my_smr_result \ --thread-num 8 \ --peqtl-smr 0.05 \ --peqtl-heidi 0.05 \ --maf 0.01 \ --diff-freq 0.2这里逐个解释参数的用意--bfile参数用来指定PLINK格式的参考面板必须是二进制的bed/bim/fam前缀。SMR会自动计算相关SNP之间的LD矩阵用于HEIDI检验。--gwas-summary指定整理好的GWAS汇总文件。--beqtl-summarize指定eQTL二进制文件前缀。--out指定输出文件前缀。--thread-num设置多线程数大文件建议至少给8否则速度感人。--peqtl-smr控制在基因层面的显著性阈值一般0.05即可如果想更严格也可以设成0.01。--peqtl-heidi控制HEIDI检验的显著性阈值0.05是官方推荐值。之前的实践中我往往还会把阈值设成0.01来进一步过滤具体看后续验证需求。--maf过滤掉低于1%的SNP。--diff-freq控制eQTL和GWAS中效应等位基因频率差的阈值如果某个SNP在两个数据集中的频率相差超过0.2会怀疑是数据错配默认是0.2通常不需要改。运行过程中SMR会在屏幕上打印当前处理到哪个染色体哪个基因区间进程速度取决于参考面板的LD计算量和你eQTL数据的规模。以GTEx全血数据加一个500万SNP的GWAS文件为例在32核服务器上大概跑30~50分钟如果是eQTLGen这种百万级eQTL的数据可能需要1~2小时。跑完后输出一堆文件最核心的是.smr结尾的结果文件以及.txt结尾的详细记录文件。3.4 输出文件解读与条件筛选SMR的结果文件是纯文本的每一行代表一个基因和某个SNP的检验结果。核心列包括Gene基因名比如ENSG编号建议后续映射回symbol名便于阅读。ProbeChr、ProbeBP探针所在位置。topSNP与基因表达关联最强的SNP的rs号。topSNP.Chr、topSNP.BP这个SNP的位置。A1、A2这个位点的等位基因。FreqA1频率。b_GWAS、se_GWAS、p_GWAS该SNP在GWAS中的效应量、标准误、P值。b_eQTL、se_eQTL、p_eQTL该SNP在eQTL数据中的效应量、标准误、P值。b_SMR、se_SMR、p_SMRSMR计算得到的因果效应估计和P值这是最核心的一列。p_HEIDIHEIDI检验的P值。nsnp_HEIDIHEIDI检验使用的SNP数量。拿到结果后我通常按以下标准筛选显著位点p_SMR 0.05 / 检验基因总数Bonferroni校正如果检验了几千个基因阈值会非常苛刻。我自己一般先用0.05做个初步筛选再对候选位点做多重检验校正避免漏掉真正有信号的基因。p_HEIDI 0.01表示没有异质性证据即GWAS和eQTL的信号很可能共享同一个因果变异。如果p_HEIDI小于0.01会再检查nsnp_HEIDI如果太少比如小于3结果不够可靠直接把该基因放弃。筛选出来的基因再检查一下效应方向是否一致——b_GWAS和b_eQTL方向一致时b_SMR也为正表示基因表达越高风险越大方向相反则为保护效应。这种方向一致性检查能避免一些怪异的错误结果。4. 结果可视化与生物学解释4.1 用曼哈顿图展示SMR信号SMR跑完只有一串数字表格可读性太差。我习惯用R来做可视化最常用的是类似GWAS曼哈顿图的SMR曼哈顿图把每个基因的p_SMR画出来横轴是基因组位置纵轴是-log10(p_SMR)。这里推荐用CMplot这个R包一行代码就能出图library(CMplot) smr_result - read.delim(my_smr_result.txt, header TRUE) # 确保有Chr和BP信息这里以topSNP的位置为例 CMplot(smr_result, plot.type m, col c(grey30, skyblue), LOG10 TRUE, threshold list(0.05/nrow(smr_result)), threshold.lty 2, threshold.lwd 1, threshold.col red, amplify TRUE, chr.den.col NULL, file pdf, file.name smr_manhattan.pdf, dpi 300, width 12, height 6)如果你的数据里基因位置信息比较稀疏也可以先按染色体分组统计每条染色体上的基因数量再用ggplot2画。但CMplot曼哈顿图的好处是它会自动标出超过阈值线的点旁边还会附上对应的基因名对后续挑选候选基因非常方便。我看曼哈顿图时会特别关注那些颜色最突出的点然后回到表格中逐个核对p_SMR和p_HEIDI。4.2 基因表达与表型因果效应的方向判断SMR分析的最终产出除了显著基因的列表还有一个关键信息是b_SMR的方向和大小。b_SMR的单位解释起来要小心——它表示的是基因表达量每变化一个标准差SD时对GWAS性状log(OR)或线性表型的影响。如果GWAS是二分类性状比如疾病b_SMR大于0表示基因表达升高会提升疾病风险小于0则相反。如果GWAS是连续性状比如血压b_SMR表示基因表达每增加1个SD血压平均变化多少。举一个实际案例我曾有个项目研究某炎症因子与冠心病的关联SMR的结果显示b_SMR为0.35p_SMR最小达到了1.2×10⁻⁶p_HEIDI为0.23这意味着基因表达水平每提高1个标准差冠心病的log(OR)增加0.35大约对应OR为1.42。这个效应量在孟德尔随机化研究中算是比较强的信号后续也通过coloc验证两信号共定位的后验概率高达0.94进一步增强了结论的可信度。解读的时候注意一个坑eQTL数据里的beta是基因表达量标准化后的效应不同组织、不同数据库的标准差定义不同因此b_SMR不能跨数据库直接比较绝对值大小。我做跨组织比较时主要看方向一致性和p值显著性不看beta数值大小。4.3 与coloc结果联动验证SMR的结果虽然方便但审稿人往往会要求额外的共定位验证。我通常的做法是SMR筛出来的显著基因全部跑一遍coloc的共定位分析。SMR和coloc的关系可以这么理解SMR告诉你有因果关联coloc告诉你两个信号是否共享同一个因果变异。如果coloc的PP.H4两个信号共定位的后验概率大于0.8那这个位点的可信度就非常高。coloc分析的输入也是GWAS汇总数据和eQTL汇总数据但需要先提取特定基因区域内的SNP。实际操作时我会从SMR结果中挑出显著基因以topSNP为中心提取上下游500kb范围内的SNP效应量整理成coloc要求的格式然后用R的coloc.abf函数逐个分析。如果某个基因在SMR中显著但coloc显示PP.H3两个信号独立的后验概率很高基本可以确认是LD导致的假阳性我一般会把它从最终列表里剔除。5. 常见问题与排查技巧实录5.1 eQTL数据与GWAS数据坐标不一致这是最让我头疼的问题也是用户在高频提问中出现最多的情况。只要GWAS和eQTL数据来自不同的参考基因组版本GRCh37 vs GRCh38染色体坐标就会整体偏移。SMR官方虽然提供了部分数据的坐标转换说明但实际运行起来如果坐标不一致通常表现为结果文件中大量基因的nsnp_HEIDI为0或者topSNP根本没有出现在eQTL数据里。排查方法也很直接随机挑几个已知的位点检查它们在两个数据中的染色体位置是否一致。如果发现不一致建议先用Liftover工具把GWAS结果转换成eQTL数据的坐标版本或者反过来转换。这里必须提醒一下Liftover转换完后要再做一次等位基因方向校验因为坐标转换过程中有可能出现链方向变化导致A/T或者C/G互换。5.2 内存不足与运行时间过长跑SMR的时候内存不足也是一个高频问题。eQTLGen的bgen文件有几十个GBGTEx虽然小一些但如果一次跑全部组织内存消耗会快速飙升。如果服务器内存不够比如16GB以下建议分染色体运行SMR最后再合并结果。SMR命令中没有专门的单染色体参数但可以通过修改输入数据的方式实现——比如用PLINK按染色体拆分参考面板同时准备对应染色体上的GWAS子集。这样并行跑22条染色体最后用cat命令合并结果文件。我遇到过最极端的情况是eQTLGen全量数据在32G内存的机器上跑到一半直接报Killed错误。后面改用三号染色体中心区域单独跑三分钟就完成了。所以遇到内存问题优先考虑分块而不是升级硬件。5.3 HEIDI检验p值几乎都为0如果发现所有结果的p_HEIDI都接近0基本可以确定是参考面板出了问题。常见原因有参考面板的人群和你GWAS数据的种族不一致或者参考面板的LD结构与eQTL和GWAS样本不匹配。我之前做东亚人群数据时用欧洲参考面板跑HEIDI就是这个结果换回东亚参考面板后p_HEIDI恢复正常。另外如果GWAS样本量和eQTL样本量有著名的重叠比如同一个队列同时贡献了GWAS数据和eQTL数据也会导致HEIDI检验p值系统性偏小这种情况建议在论文里说明局限性或者尝试用独立数据集做验证。5.4 结果文件中大量基因缺失如果你发现输出的.smr结果文件里基因数量远少于预期多半是过滤阈值设置得太严。特别是--peqtl-smr如果设置成0.001很多基因连被检验的资格都没有。还有一种可能是你的GWAS数据格式里SNP命名方式和eQTL数据不一致——比如GWAS用的是chr:pos:ref:alt格式而eQTL是rs号SMR内部做ID匹配时如果对不上直接跳过。这种问题排查方法很简单取一个eQTL数据的top SNP在GWAS数据里搜索一下看看能不能找到、格式是否一致很快就能定位问题。5.5 多重检验校正与结果报告规范最后说一个报告层面的问题。SMR一次会检验成千上万个基因如果不做多重检验校正假阳性会爆炸。常用做法是FDR比如q值小于0.05或Bonferroni0.05除以检验基因数。我自己的习惯是用Bonferroni做严格筛选然后在补充材料里报告FDR结果给审稿人留出宽松的空间。论文方法部分需要写清楚使用的eQTL数据版本、样本量、过滤标准、参考面板版本以及SMR软件的版本号和关键参数这些细节直接决定结果的可重复性。写在最后的一点心得SMR这套流程我前前后后跑了不下几十遍。最开始只把它当个黑盒工具输两个文件出来一张表格后来踩的坑多了才逐渐明白每一步参数背后都是统计学假设和生物学常识的博弈。如果你现在正准备跑自己第一份SMR分析我的建议是不要急着上全基因组的大数据先挑一个明确的候选基因区域跑通整个流程看清楚每一步输出的中间文件和结果字段再扩展到全基因组范围。这样至少你能在半小时内完成一次完整的分析迭代而不是在等了两小时任务跑完之后才发现第一步的数据格式就错了。另外一个小技巧SMR和coloc是可以互相印证的但两者的前提假设并不完全一样建议在正式分析之前用自己领域内一篇已发表的高质量论文里的某个位点做一次“标定”——用同一套数据看看你能不能复现出对方的SMR信号。能复现出来你的数据整理和参数设置基本就没有问题了。技术细节说再多都不如自己亲手跑通一遍来得踏实。希望这份指南能帮你少走一些弯路。 
