BOLT-LMM：十万级样本GWAS高效关联分析的原理与实战-育师

如果你手头的数据已经大到需要为“跑完一次GWAS要几天”发愁，BOLT-LMM就是那种能把时间压缩到几小时的工具。它由Broad Institute团队开发，专门面向几十万样本规模的混合模型关联分析。我最早是在一个约35万样本的队列里遇到性能问题的，当时对比了GEMMA、FaST-LMM和BOLT-LMM三套方案，最后真正能在一晚上出全基因组结果的只有BOLT-LMM。这篇笔记是我通读其论文和官方文档之后，按自己的理解整理的原理摘要和安装使用记录，希望能帮正在选型或卡在编译这步的朋友。

1. 当GWAS样本量冲到几十万，传统工具先顶不住了

1.1 从简单回归到混合模型的演进

先说说为什么GWAS分析会走到混合模型这一步。早年的全基因组关联分析面对的数据规模不大，几千个样本、几十万个SNP，用PLINK跑logistic regression或者线性回归就够了。但数据规模一上来，问题就暴露了：人群分层（population stratification）和样本间的潜在亲缘关系会把关联信号整体抬高，直接导致大量假阳性。为了压制这种结构效应，有人把主成分（PCA）作为协变量放进去，但PCA只能吸收掉一部分祖先成分，并不能完全刻画个体间的复杂亲缘关系。

线性混合模型（LMM）就是把个体间的关系写进模型的解法，它用遗传关系矩阵（GRM，Genetic Relatedness Matrix）作为随机效应的协方差结构。大家熟悉的GEMMA、FaST-LMM、GCTA都是走这条路。数学形式很简单：对每个SNP，假设Y = Xβ + Sg + ε，其中g服从均值为0、方差为σg²K的正态分布，K是基因型矩阵计算出的亲缘关系矩阵。这里的关键是，加入随机效应项后，检验每个SNP时都要考虑整个K矩阵的逆，这在大样本下非常昂贵。

1.2 GEMMA们哪里不够用

如果样本只有一两万，GEMMA这类工具其实够用。问题在于，现在很多项目拿到的都是生物库级别的大队列，动辄十几万、几十万人，SNP数也有数百万甚至上千万。精确的LMM算法需要反复对N×N的矩阵做特征分解或者Cholesky分解，时间复杂度在O(MN²)附近，十万样本下基本就是天文数字。有人会退回到PCA校正的简单回归，但损失了混合模型对隐性亲缘关系的校正能力，对统计学审稿人也不好交代。

我梳理过一个简单的工具对比表格，方便后面选型参考：

工具	算法类型	大体量数据表现	适用场景
PLINK常规回归	固定效应模型	速度快但校准差	小样本、初筛
GEMMA	精确LMM	万级样本可接受	小样本高精度
FaST-LMM	近似LMM	中大规模	中规模数据折中
BOLT-LMM	近似LMM+LD Score回归	十万级样本依然高效	生物库级别GWAS

这张表不是官方给的，是我自己跑了几个数据集后总结出来的选型感觉。BOLT-LMM的最大卖点，就是把这个计算瓶颈解决得比较优雅。它由Broad Institute的Po-Ru Loh等人开发，最早论文发表于Nature Genetics 2015年。它的核心思路是：先利用LD Score回归估算模型参数，再用一种近似的变分推断方法拟合混合模型，最后通过低秩近似计算每个SNP的关联统计量。整套流程让原本O(MN²)的大矩阵计算变成若干个轻量级矩阵运算，这也是它敢说“数十万样本、全基因组扫描只要几小时”的底气。

1.3 我判断要不要用BOLT-LMM的场景标准

在决定是否引入BOLT-LMM之前，我当时列了一个很简单的判断清单，这里分享出来：

样本量：如果有效样本超过5万，BOLT-LMM的收益就很明显；低于1万的话，GEMMA之类的精确LMM完全够用。
数据规模：基因组范围的标记数越大，BOLT-LMM在时间上的优势越突出。
性状遗传结构：高度多基因的性状（身高、BMI、很多免疫指标）特别适合BOLT-LMM，因为它的统计量设计本身考虑了多基因背景。
计算资源：BOLT-LMM对内存仍有要求，但比GEMMA等更可控；没有高性能集群也能跑。

这个清单不是官方的，是我自己踩过几轮后总结出来的选型标准。如果你的场景落在前三项中的多数，直接往下看安装部分，大概率值得。

2. 拆开BOLT-LMM的黑盒：两步近似法到底在算些什么

2.1 线性混合模型在GWAS里到底干了什么

很多教程直接把混合模型公式一贴，读者看完了还是不知道它“校正了什么”。我自己读文献时，脑子里一直想找一个直观类比。后来想出一个还算凑合的比喻：简单回归相当于在人群里按每个SNP逐个“拉关系”，但人群里本来就有亲戚、同乡、同族这些复杂关系，如果你不管这些，看见某个等位基因在某个家庭里频率高，就误以为它跟某个性状有关系。加PCA就是把大家按“大致来自哪几个祖先群体”分了个类，但同一个类别内部还有细微的亲缘结构。

LMM相当于把“任意两个人的亲缘程度”作为一张完整的网络图铺开，在网络里做关联检验。这张网络图就是遗传关系矩阵K。理论上，只要K估计得准，SNP的检验统计量就不会再被个体间的亲缘关系污染。代价是，模型里多了一个σg²，每次检验都要跟这个大矩阵打交道。传统精确LMM之所以慢，就是因为每个SNP的检验都要基于完整的K矩阵重新做分解。

2.2 LD Score回归为什么被BOLT-LMM当底座

读BOLT-LMM的论文时，我觉得最巧妙的一步是引入LD Score回归。很多人听到LD Score第一反应是LDSC这个做遗传相关性估计的工具，确实，BOLT-LMM参考了同一套数学框架。LD Score回归最基础的原理是：一个SNP的卡方统计量，它的大小既取决于这个SNP与因果变异的关联，又取决于这个SNP的LD Score——即它“连累”了多少邻居位点。

因为LD Score可以在参考面板上提前算好，几乎不增加分析成本。BOLT-LMM的思路是先用LD Score回归在全基因组范围内估计混合模型的参数（包括遗传力占比σg²/σp²），得到参数后就直接进入第二步近似计算，而不是像传统LMM那样用EM算法反复迭代。这一步非常关键，它把“估计方差组分”这项传统LMM中代价最高的步骤，转化成了一个轻量级的回归问题。整个算法的时间开销由此降了一个数量级。

2.3 变分推断和低秩近似解决了什么

BOLT-LMM还用到了变分贝叶斯。简单说，变分推断就是用一族简单的分布去近似难解的后验分布，再用这个近似分布计算每个SNP的检验统计量。官方文档里的原话是“a variational Bayes approach to fit the Gaussian mixture model”，目的是把随机效应项的后验均值算出来，之后每个SNP的似然比检验只需要做一次低秩更新，不需要重新求解全模型。这也是BOLT-LMM能够做到每个SNP检验又快又稳的核心。

不过要提醒一句：这个“近似”是有适用条件的。论文中验证的主要是常见变异（MAF通常高于0.1%且频率分析时小心处理）和由常见变异解释的多基因性状。如果你的研究重点在低频变异或单基因病那样的极端效应结构，BOLT-LMM的近似效果可能会打折扣，这时候宁愿花时间跑精确LMM或采用专门的低频变异分析方法。

2.4 BOLT-LMM与BOLT-LMM-inf，选哪个统计量

跑BOLT-LMM时你会注意到输出文件里有两组统计量：BOLT-LMM和BOLT-LMM-inf。这是我刚开始使用时比较困惑的地方，后面读文档才搞清楚。BOLT-LMM是在“有限标记数”模型下计算的关联统计量；BOLT-LMM-inf则假设标记数量趋近于无穷，也就是说它把未观测到的因果变异也纳入模型，理论上在多基因性状上统计效力更高，对人群分层的控制也更严格。

实际分析中，两者可以同时出结果。官方建议和多数应用经验是：如果研究性状是典型多基因结构，BOLT-LMM-inf更准确；如果样本量不大或运行开销敏感，那就以BOLT-LMM为准。两种统计量的P值可以画在同一个QQ图里看总体验证情况。如果两者差异非常大，通常说明模型的某个前提没满足，比如LD Score文件人群不匹配或表型分布异常。

3. 快速安装：从依赖到可执行文件

3.1 最容易被忽略的编译依赖

BOLT-LMM的安装我自己踩过一次坑，原因是前几年在一台比较旧的CentOS服务器上，gcc版本太低，编译直接报了一堆模板错误。所以先列依赖清单，这条经验很重要：

操作系统：Linux（Ubuntu 18.04以上或CentOS 7以上都比较顺畅）
编译器：gcc/g++版本建议在5.0以上，越新越好。BOLT-LMM是C++写的，对C++11/14特性有依赖
make工具：系统一般自带
Boost库：编译时需要用boost头文件，尤其regex等组件
数学库：BLAS/LAPACK或OpenBLAS
压缩和下载相关：zlib、libcurl

在Ubuntu系统上可以直接用apt装：

sudo apt-get update sudo apt-get install build-essential g++ make zlib1g-dev libcurl4-openssl-dev libopenblas-dev libboost-all-dev

如果是CentOS/RHEL：

sudo yum install gcc gcc-c++ make zlib-devel libcurl-devel openblas-devel boost-devel

这里有个细节：旧版CentOS的默认gcc可能只有4.8，而BOLT-LMM某些版本要求C++11标准支持完整，升级gcc这一件事就能解决后续很多编译报错。我自己后来干脆用Developer Toolset在新系统上编译，省了不少事。

3.2 下载源码与编译

BOLT-LMM目前通过GitHub发布源码，整个项目比较规整，没有复杂的autotools流程。下载后解压，进目录看README，会发现编译命令出奇地简单——就是make：

wget <安装包链接> tar -xzf BOLT-LMM*.tar.gz cd BOLT-LMM_v* make

这里注意，下载时要选对版本。新版本对BGEN格式支持更好，表型文件缺失值的容错也更强。make之后会在当前目录生成可执行文件BOLT-LMM，可以用ls -l确认一下权限和大小。

有一点值得说明：官方仓库里已经捆绑了一些必要的头文件和辅助工具，比如计算LD Score时会用到的一些脚本。所以装的时候不要把个别脚本单独拎出来跑，否则后续会找不到依赖路径。

3.3 运行自检：验证安装成功

编译完不要直接跑大数据，先用--help做一次自检：

./BOLT-LMM --help

正常会输出一大段参数说明。如果这里报错“error while loading shared libraries”，多半是某个动态库没找到，用ldd BOLT-LMM看一下缺哪个，再补装对应的库。更稳妥的办法是把源码目录下自带的example或者自测数据跑一遍。我自己验证安装时习惯先构造一个1000样本的小模拟数据，用--bed、--phenoFile跑一下，确认能输出结果文件再上有价值的数据集。

4. 跑一个真实的GWAS任务：输入文件、命令行与结果解读

4.1 五类核心输入文件

在BOLT-LMM里，输入文件比PLINK稍多一点，我盘点一下：

基因型文件：最常用的是PLINK的bed/bim/fam三件套，或者BGEN格式。如果数据来自imputation，直接给BGEN比较方便。BIM文件里的等位基因做统一朝向，A1通常是被检验的效应等位基因。
表型文件：文本格式，至少三列：FID、IID、表型值。支持多个表型一起放，运行时用--phenoCol指定。
协变量文件：FID、IID加协变量列。可以是连续变量（用--qCovarCol标记，q代表quantitative），也可以是分类变量（用--covarCol标记）。
LD Score文件：通常从官方提供的参考面板下载，或者用配套脚本基于1000 Genomes等参考面板计算。这个文件必须和样本的人群来源匹配。
遗传图谱文件：提供物理位置到遗传位置的映射，官方一般随安装包提供或单独下载。

4.2 一个可以直接复制修改的命令行

下面这个命令是我在本地Linux服务器上验证好的模板：

./BOLT-LMM \ --bed=ukb_chr1_22.bed \ --bim=ukb_chr1_22.bim \ --fam=ukb_chr1_22.fam \ --phenoFile=pheno_bmi.txt \ --phenoCol=bmi \ --covarFile=covars.txt \ --covarCol=sex \ --covarCol=age \ --qCovarCol=age \ --LDscoresFile=eur_ldscores_hm3.txt.gz \ --geneticMapFile=genetic_map_hg19.txt \ --numThreads=10 \ --maxMissingPerSnp=0.02 \ --minMAF=0.001 \ --statsFile=bolt_bmi.txt \ --verbose

参数含义不用全背，重点记住几个：

--bed/--bim/--fam：基因型输入
--phenoFile/--phenoCol：指定表型文件及列名
--covarFile/--covarCol/--qCovarCol：协变量；没有分类变量时可省略covarCol
--LDscoresFile：LD Score路径，压缩的.gz文件也能直接读
--geneticMapFile：遗传图谱，缺少会报错
--numThreads：多线程加速，建议设置为你机器物理核数的一半到全部
--maxMissingPerSnp和--minMAF：SNP质控阈值，跟PLINK里的MAF过滤概念一致
--statsFile：结果输出路径

如果你有显式的亲缘关系矩阵需要强制校正，还可以加--GRM文件参数；大多数常见场景下不手动指定，BOLT-LMM会基于样本SNP自动完成计算。

4.3 结果文件怎么看

跑完后，--statsFile指定的文件就是核心结果。它一般包含下面这些列：

CHR、SNP、BP、GENPOS：染色体、SNP名、物理位置、遗传位置
ALLELE1/ALLELE0：效应等位基因/另一个等位基因
A1FREQ：效应等位基因频率
BETA、SE：效应量和标准误
CHISQ：卡方统计量
P_BOLT_LMM_INF和/或P_BOLT_LMM：两种统计量对应的P值

我拿到结果后，第一件事不是画曼哈顿图，而是计算基因组膨胀因子lambda。最快捷的方法是用R读入结果，取P值列转成卡方值，再除以卡方分布0.5分位数。理想情况lambda在1.0附近，如果超过1.1，说明统计量整体偏高，通常要先怀疑人群分层未校正干净或LD Score不匹配；低于0.9则可能你的质控过滤过严或样本量太小。

5. 实战中躲不开的坑与我的处理思路

5.1 表型文件里的“隐形炸弹”

BOLT-LMM读表型文件时对格式的要求有时比较严格，我自己被绊倒过好多次。首先是分隔符，官方支持空格或制表符，但不支持逗号；如果你从Excel直接导出CSV再改名，很容易在这里报错。其次是缺失值，官方推荐用NA表示，不要留空白或者写0，0会被当成真实的表型值参与分析。还有一点容易忽略：表型文件里的FID和IID必须与fam文件完全一致，顺序无所谓，但ID不能多不能少。如果样本ID对不上，BOLT-LMM会直接报错退出，不会自动帮你对齐。

实际操作中，我习惯先写一小段R代码统一检查表型文件和fam文件的ID交集，确认没有差异再提交任务，这一步能避免大量无效排队。

5.2 LD Score与参考面板人群不匹配

BOLT-LMM对LD Score文件的要求是：官方提供的文件本身按人群区分，比如欧洲人群、非洲人群、东亚人群等。如果你的样本是混合人群或来自中国南方某地队列，直接套用欧洲人群的LD Score，最典型的表现就是统计量膨胀或者紧缩。我的处理办法是：先用PCA看样本的祖先成分，确定最接近的参考人群，再选择对应LD Score；如果实在没有匹配的，可以自己基于参考面板基因型计算。这一步我一般放在正式全基因组扫描前，因为返工成本太大。

症状判断上有一个很实用的小技巧：如果结果里几乎所有SNP的P值都偏小，优先怀疑LD Score不匹配；如果只是个别区域膨胀，那更可能是结构变异或拷贝数区域的影响。

5.3 内存、线程与运行时间管理

BOLT-LMM虽然比很多工具省内存，但几十万样本下依然要准备充足的RAM。我遇到过的经验值差不多是：10万样本、全基因组扫描，每个线程占用几个GB内存，且随SNP数和线程数增长。为了稳妥，我的建议是先把--numThreads设小跑一小段，观察内存使用，再决定要不要加线程。内存不足时优先减少线程数，其次考虑分染色体跑，最后再用--memEstimate参数让程序先估算内存需求。这里务实地说，BOLT-LMM自己的内存估算功能挺好用，分配节点前跑一次能省很多冤枉时间。

5.4 迭代不收敛的排查路径

BOLT-LMM在运行时会输出一系列迭代日志，偶尔会提示模型不收敛或者方差组分估计异常。我遇到过的情况主要有三种：

表型分布严重偏离正态：比如原始计数数据没有做逆正态变换，BOLT-LMM对这类表型拟合会比较吃力，建议先做rank-based inverse normal transformation。
遗传力接近0：如果性状几乎不受遗传影响，σg²的估计会非常不稳定，运行时间反而变长。
样本间亲缘关系过密：数据里包含大量一级亲属时，GRM中会出现很大的块结构，导致低秩近似失效，建议先做亲缘关系剪枝，保留无亲缘关系样本。

遇到不收敛时，不要急着调参，先把表型分布和样本亲缘关系这两个基础问题排查掉，80%的情况都能解决。

我在实际项目中反复使用BOLT-LMM之后，一个比较深的体会是：工具的快速安装只是第一步，真正让分析结果站得住脚的，是你对模型假设的理解和对输入数据质量的把控。上面这些坑，大多不是从官方手册里直接能看到的，而是要在真实数据集上反复试错才能积累下来。希望这篇笔记能帮你少走一点弯路，把时间花在更有价值的生物学解读上。