系统学习下BOLT-LMM的软件手册,
1 概述
BOLT-LMM软件包目前由两种主要算法组成,即用于混合模型关联分析的BOLT-LMM算法和用于方差分量分析(即SNP遗传性的分区和遗传相关性的估计)的BOLT-REML算法。
我们推荐BOLT-LMM用于分析包含5,000多个样本的人类遗传数据集。BOLT-LMM中使用的算法依赖于仅在大样本量下成立的近似值,并且仅在人类数据集中进行了测试。对于少于5,000个样品的分析,我们建议使用GCTA或GEMMA软件。
我们还注意到,BOLT-LMM 关联测试统计量对于定量性状和(合理地)平衡的病例对照性状有效。 对于不平衡的病例控制特征,我们建议使用 SAIGE 软件(参见部分10[1]进行全面讨论)。
1.1 BOLT-LMM混合模型关联测试
BOLT-LMM 算法使用线性混合模型 (LMM) [1[2]] 计算用于分析表型和基因型之间关联的统计。默认情况下,BOLT-LMM 假设归因于被测试 SNP 以外的随机效应之前有一个正态值的贝叶斯混合效应。该模型推广了以前的混合模型关联方法(例如,EMMAX [2[3]],FaST-LMM [3[4],4[5],5[6],6[7]],GEMMA [7[8]],GRAMMAR-Gamma [8[9]],GCTA-LOCO [9[10]])使用的标准”无穷小”混合模型,为在控制误报的同时提高检测关联的能力提供了机会。此外,BOLT-LMM 应用算法改进来计算混合模型关联统计量的速度比基于特征分解的方法快得多,无论是在使用贝叶斯混合模型时,还是在专用于标准混合模型关联时。BOLT-LMM 在参考文献 [1[11]] 中进行了描述:
此外,参考文献[10[12]]探讨了BOLT-LMM在完整的英国生物样本库数据集上的性能:Loh P-R, Kichaev G, Gazal S, Schoech AP, and Price AL. Mixed model association for biobank-scale data sets. Nature Genetics, 2018.
1.2 BOLT-REML方差分量分析
BOLT-REML算法估计遗传力,由基因分型SNP和在同一组个体上的多个性状之间的遗传相关性来解释。与GCTA软件[11]一[13]样,BOLT-REML应用方差成分分析来执行这些任务,支持多成分建模以划分SNP遗传性和多性状建模以估计相关性。BOLT-REML在大样本量下应用蒙特卡罗算法,该算法比基于特征分解的方法更快地进行方差分量分析(例如GCTA)。BOLT-REML 在参考文献 [12[14]] 中进行了描述。
2 下载和安装
您可以在以下位置下载最新版本的 BOLT-LMM 软件:
http://data.broadinstitute.org/alkesgroup/BOLT-LMM/downloads/[15]
以前的版本也在old/子目录中提供。
2.1 更新日志
-
版本 2.3.6(2021 年 10 月 29 日): -
修复了在线性回归输出中缩放 BETA 和 SE 列时出现的错误。此错误仅影响 BOLT-LMM v2.3.5 在线性回归模式下为具有非单位方差的表型计算的效应大小;线性混合模型分析(–lmm/–lmmInfOnly/–lmmForceNonInf)的效应大小不受影响。 -
版本 2.3.5(2021 年 3 月 20 日): -
提高了表型/协变量文件处理效率(仅加载请求的列)。 -
添加了 BETA 和 SE 列,以便在 BOLT-LMM 以线性回归模式运行时输出。 -
版本 2.3.4(2019 年 8 月 10 日): -
阻止环境变量MKL_NUM_THREADS覆盖 –numThreads 参数。 -
版本 2.3.3(2019 年 8 月 3 日): -
添加了对 BGEN v1.2 数据中缺失值的支持。 -
完成模型拟合步骤后减少了内存使用量(通过在计算填充变异的关联测试期间释放不再需要的基因型)。 -
将 BLAS 库更新为英特尔 MKL 2019 Update 4,并修改了 bolt 可执行文件以动态链接英特尔线程库(libiomp5.so,现在随 BOLT-LMM 一起分发)。 -
改进的错误报告和文档;添加了本手册的常见问题解答部分。 -
版本 2.3.2(2018 年 3 月 10 日): -
在英国生物库v3填充发布中增加了对填充X染色体变异的支持(请参阅部分9.1[16]有关详细信息)。 -
添加了 –bgenSampleFileList 选项,以允许具有不同示例文件的多个 bgen/示例文件对。 -
将 –allowX 选项更改为始终设置(因此 –allowX 不再需要用于 X 染色体分析)。 -
添加了用于非人类分析的 Nautosomes 选项。 -
版本 2.3.1(2017 年 12 月 19 日): -
添加了在指定 –covarCol 或 –qCovarCol 时对 –covarFile 的检查。(以前,如果未指定 –covarFile,则这些参数将被静默忽略。 -
添加了用于 X 染色体分析的 –allowX 选项。 -
增加了有关病例控制特征分析的文档部分。 -
已将 BLAS 库更新到英特尔 MKL 2018 更新 1。 -
改进了错误报告和文档。 -
版本 2.3(2017 年 8 月 1 日): -
添加了对英国生物银行使用的BGEN v1.2填充文件格式的多线程支持。 -
添加了 –bgenMinINFO 参数。 -
对于先前四舍五入为零的极小 P 值,固定 P 值输出。 -
修复了由很长的BGEN等位基因名称引起的崩溃。 -
修复了 –noBgenIDcheck 标志中的错误。 -
改进了错误报告。 -
本文档末尾增加了一节,为N=50万英国生物样本库分析提供了建议。 -
版本 2.2(2015 年 11 月 13 日): -
添加了对以 BGEN 格式测试填充 SNP 的支持。 -
添加了按基本对坐标而不是 SNP 名称 (–LDscoresMatchBp) 查找 LD 分数的选项。 -
修复了hg19遗传图插值中的错误。 -
修复了QC过滤器中每个样本缺失率的错误。 -
改进了错误检查。 -
版本 2.1(2015 年 4 月 29 日): -
改进了对IMPUTE2文件的处理(大加速;INFO 输出列而不是F_MISS;MAF 过滤)。 -
修复了每个样本缺失筛选器 (–maxMissingPerIndiv) 中被忽略的错误。 -
对输入参数进行了细微的更改(删除了 –impute2CallThresh,这没有任何影响;添加了 –impute2MinMAF 和 –h2gGuess)。 -
版本 2.0(2015 年 3 月 13 日):添加了用于估计遗传性参数的 BOLT-REML 算法。修复了阻止 BOLT-LMM 在某些系统上运行的参数初始化错误。对参数检查进行了各种小改进。 -
(2014 年 12 月 8 日):GPLv3 下的许可源代码。 -
版本 1.2(2014 年 11 月 4 日):添加了对以 2 剂格式(Ricopili/plink2 format=2)测试填充 SNP 的支持。修复了导致nan遗传力估计的错误。 -
版本 1.1(2014 年 10 月 17 日):添加了对测试具有概率剂量的填充 SNP 的支持。 -
版本 1.0(2014 年 8 月 8 日):初始版本。
2.2 安装
这里补充下,conda大法好需要什么软件,可以conda上先搜索下有没有,免去折腾软件,咱是干活的,不是配置环境的:Anaconda.org[17]
# 新建个环境,防止冲突啦
conda create -n bolt-lmm
conda install -c bioconda bolt-lmm
如果你用conda安装的话,下边的内容就不用看了。
BOLT-LMM_vX.X.tar.gz下载包包含一个预编译的 64 位 Linux 可执行文件 bolt,我们已经在几个 Linux 系统上对其进行了测试。我们建议使用此可执行文件,因为它经过了良好的优化,不需要进一步安装。请注意,从 BOLT-LMM v2.3.3 开始,bolt 可执行文件动态链接 libiomp5.so 英特尔线程运行时库;此共享库在 BOLT-LMM 包的 lib/子目录中提供,并将由 bolt 可执行文件从该子目录中自动加载。
如果您希望从源代码(在 src/ 子目录中)编译自己的 BOLT-LMM 软件版本,则需要确保满足库依赖关系,并且需要对 Makefile 进行适当的修改:
-
库依赖项: -
BLAS/LAPACK 数值库。BOLT-LMM 软件的速度主要取决于与其相关的 BLAS/LAPACK 实现的效率。我们建议使用英特尔数学核心函数库 (MKL)(如果可用)(AMD 处理器除外);否则,ATLAS可能是一个很好的选择。 -
提升C++库。BOLT-LMM 链接指向 Boost program_options和 iostreams 库,这些库需要在下载并解压缩 Boost 后安装。 -
NLopt 数值优化库 [13[18]]。 -
生成文件:需要相应地修改库的路径。请注意,已发布版本的 Makefile 不会设置标志 -DUSE_MKL_MALLOC。此标志将打开英特尔 MKL 的快速内存管理器(用 mkl_malloc 替换对_mm_malloc的调用),这可能会提高内存性能,但我们观察到某些系统在使用mkl_malloc时崩溃。
作为参考,提供的 bolt 可执行文件是在哈佛医学院”Orchestra 2″研究计算集群上构建的,使用英特尔 icpc 16.0.2(带有 MKL 2019 Update 4)和 Boost 1.58.0 和 NLopt 2.4.2 库,方法是调用 make linking=static-except-glibc。
移植到 Windows。Remi Daviet创建了一个在Windows中编译的BOLT-LMM版本:http://remidaviet.com/software.php[19]
移植到 FreeBSD。BOLT-LMM 可以通过 FreeBSD ports system (pkg install bolt-lmm) 安装在 FreeBSD 上。请注意,此安装将仅使用高度可移植(并且可能不太快)的优化。
2.3 运行
要运行 bolt 可执行文件,只需在 Linux 命令行(在 BOLT-LMM 安装目录中)调用./bolt
,参数格式为--optionName=optionValue
example/子目录包含一个 bash 脚本run_example.sh
,该脚本演示了 BOLT-LMM 在小型示例数据集上的基本用法。同样,run_example_reml2.sh
演示了 BOLT-REML。
最小的 BOLT-LMM 调用如下所示:
--lmm --LDscoresFile=tables/LDSCORE.1000G_EUR.tab.gz
--statsFile=stats.tab
最小的 BOLT-REML 调用如下所示:
./bolt --bfile=geno --phenoFile=pheno.txt --phenoCol=phenoName --reml--modelSnps=modelSnps.txt
要执行多性状 BOLT-REML(即估计遗传相关性),请提供多个 –phenoCol=phenoName 参数。
2.5 帮助
要获取基本选项的列表,请运行:./bolt -h
要获取基本和高级选项的完整列表,请运行:./bolt --helpFull
3 计算要求
3.1 操作系统
目前,我们只在Linux计算环境中编译和测试了BOLT-LMM;但是,如果您希望尝试为其他操作系统编译 BOLT-LMM,则可以使用源代码。
3.2 内存
对于典型的数据集(M,N超过10,000),BOLT-LMM和BOLT-REML使用大约MN/4字节的存储器,其中M是SNP的数量,N是个体的数量。更准确地说:
-
M = bim 文件中满足所有条件的 SNP 的数量: -
未在任何 –exclude 文件中列出 -
通过 QC 过滤器筛选缺失 -
在 –modelSnps 文件中列出(如果已指定) -
N = fam文件中未在任何 –remove 文件中列出的个体数(但 QC 之前;即N 包括由于缺少基因型或协变量而被过滤的个体)
3.3 运行时间
在实践中,BOLT-LMM 和 BOLT-REML 的运行时间大致与MN 1.5 差不多。对完整的英国生物银行数据集(M ~ 700K SNP,N = 500K个体)的分析通常需要几天时间,使用单个计算节点的8个线程;有关更多详细信息,请参阅参考文献。 [1[20], 10[21]]。
3.3.1 多线程
在多核计算机上,通过使用--numThreads
选项调用多线程,可以减少运行时间。
4 输入/输出文件命名约定
4.1 自动 gzip [取消]压缩
BOLT-LMM 软件假定以 .gz结尾的输入文件经过 gzip 压缩,并即时自动解压缩(无需创建临时文件)。同样,BOLT-LMM 将 gzip 压缩的输出写入任何以.gz结尾的输出文件。
4.2 输入文件和协变量数组
顺序编号的输入文件和协变量的数组可以用速记 {i:j} 指定。例如,
data.chr{1:22}.bim
被解释为文件列表
data.chr1.bim, data.chr2.bim, ..., data.chr22.bim
5 输入
5.1 基因型
BOLT-LMM 软件采用 PLINK [14[22]] 二进制格式(fam/bim/fam)的基因型输入。对于一般的文件转换和数据操作,我们强烈推荐PLINK/PLINK2软件[15[23]]。
如果所有基因型都包含在具有相同文件前缀的单bed/bim/fam 文件三元组中,则只需使用命令行选项 –bfile=prefix 即可。基因型也可以通过使用多个 –bed 和 –bim 调用或使用上述文件数组速记(例如,–bim=data.chr{1:22}.bim),将基因型拆分为包含连续 SNP 集的多个bed和 bim 文件(例如,每个染色体一个床/bim 文件对)。
5.1.1 参考遗传图谱
BOLT-LMM 包包括参考图,如果您的 PLINK bim文件不包含遗传坐标(以摩根为单位),您可以使用这些参考图来插值来自 SNP 物理(碱基对)位置的遗传图谱坐标。(BOLT-LMM关联测试算法使用遗传位置来防止近端污染;BOLT-REML 不使用此信息。要使用参考map,请使用以下选项
--geneticMapFile=tables/genetic_map_hg*.txt.gz
选择与 bim 文件的物理坐标对应的版本(hg17、hg18、hg19 或 hg38)。您可以使用--geneticMapFile
选项,即使您的 PLINK bim 文件确实包含遗传坐标;在这种情况下,将忽略 bim 文件中的遗传坐标,而是使用插值坐标。
5.1.2 估算的SNP剂量
BOLT-LMM 关联测试算法支持使用基于检测的 PLINK 格式基因型子集(通常是直接基因型的子集)构建的混合模型,在任意数量的填充 SNP(具有实值”dosage”而不是检测基因型)下计算混合模型关联统计信息。(BOLT-REML 方差成分分析不支持剂量输入。)
在测试填充 SNP 时,BOLT-LMM 首先对 PLINK 格式的基因型(通过 –bfile 或 bed/bim/fam 提供)执行其通常的模型拟合,然后应用该模型扫描任何提供的填充SNP。第二步只需要适度的额外计算,不需要额外的RAM,因为它只是对BOLT-LMM在模型拟合期间计算的残余表型执行实值剂量SNP的基因组扫描(如GRAMMAR-Gamma [8[24]])。我们目前建议对~500K为检测基因型进行模型拟合;这种方法几乎不会牺牲任何统计能力,同时保持计算效率。
如果您手头只有填充的 SNP 数据,则需要预处理数据集,以便为 BOLT-LMM 创建 PLINK 格式的检测 SNP 子集。我们建议以下程序。
-
确定一组高可信度的 SNP(例如,基于 R2 或 INFO 分数),以创建初始检测基因集。 -
在这些 SNP 上以 PLINK 格式创建检测基因型。 -
使用 PLINK 将 LD 修剪至 ~500K SNP(通过 –indep-pairwise 50 5 r2 thresh 获得适当的 r2thresh)。 -
使用最终的检测 SNP 作为 –bfile(或 bed/bim/fam)参数运行 BOLT-LMM,使用以下格式之一将填充的 SNP 指定为附加关联测试 SNP。
填充SNP的剂量格式。此输入格式由一个或多个 –doseFile 参数组成,这些参数指定在填充 SNP 下包含实值基因型期望的文件。–doseFile 的每一行应按如下格式设置:
rsID chr pos allele1 allele0 [dosage = E[#allele1]] x N
缺失(uncalled)的剂量可以用 –9 指定。您还需要提供一个额外的 –doseFidIidFile,指定剂量对应的样品的 PLINK FID 和 IID。有关示例,请参阅example/子目录。
impute2格式。您也可以指定IMPUTE2 软件 [16[25]] 的输出。IMPUTE2 基因型文件格式如下:
snpID rsID pos allele1 allele0 [p(11) p(10) p(00)] x N
(BOLT-LMM 忽略 snpID 字段。在这里,每个基因型条目都包含个体是等位基因纯合子1,杂合子和等位基因纯合子0的个体概率,而不是剂量。这三个概率不需要总和为1,允许基因型不确定性;如果概率之和小于 –impute2CallThresh 参数,则 BOLT-LMM 将基因型视为缺失。
若要计算包含 IMPUTE2 SNP 的文件列表中的关联统计信息,可以列出--impute2FileList
文件中的文件。此文件的每一行应包含两个条目:一个染色体编号,后跟一个包含该染色体的 SNP 的 IMPUTE2 基因型文件。您还需要提供一个额外的 –impute2FidIidFile,指定 IMPUTE2 基因型所对应的样本的 PLINK FID 和 IID。有关示例,请参阅example/子目录。
填充SNP的2剂量格式您也可以指定为** Ricopili 流程和 plink2 –dose format=2 的输出**。此文件格式由文件对组成:(1)包含有关SNP位置信息的PLINK map文件;和(2)2剂量格式的基因型概率文件,由表头组成
SNP A1 A2 [FID IID] x N
后跟格式中每个 SNP 一行
rsID allele1 allele0 [p(11) p(10)] x N
每个条目的第三个基因型概率假定为 p(00)=1-p(11)-p(10)(与 IMPUTE2 格式不同)。要计算 2-dose 文件列表中 SNP 的关联统计信息,可以列出 –dose2FileList 文件中的文件。此文件的每一行应包含两个条目:一个 PLINK map文件,后跟相应的基因型文件,其中包含这些 SNP 的概率(像往常一样,如果任一文件以 .gz 结尾,则会自动解压缩;否则假定它是纯文本。有关示例,请参阅example/子目录。
BGEN 格式的填充 SNP。要计算一个或多个 BGEN 数据文件中 SNP 的关联统计信息,请使用 --bgenFile
指定 .bgen 文件,并使用--sampleFile
指定相应的.sample
文件。--bgenMinMAF
和--bgenMinINFO
选项允许将输出限制为通过最小等位基因频率和 INFO 阈值的 SNP。(注意:--bgenMinMAF
过滤应用于完整的 BGEN 文件(在任何样本排除之前),而 BOLT-LMM 输出中报告的 MAF 是在实际分析的样本子集中计算的。因此,某些 SNP 可能会通过--bgenMinMAF
过滤器,但在输出文件中报告的 MAF 较低;如果您希望排除此类 SNP,则需要对结果进行后处理。
请注意,从BOLT-LMM v2.3 开始,--bgenFile
选项允许多个 BGEN 文件。我们已经对英国生物银行 N=500K 版本中使用的BGEN v1.2 格式的文件实施了多线程处理,因此现在可以在单个作业中分析所有染色体的 BGEN v1.2 数据。对于 N=150K 版本中使用的 BGEN v1.1 数据的分析,我们建议染色体并行化以提高计算便捷度(使用来自每个作业中所有染色体的直接基因分型 PLINK 数据的完整 –bfile)。
此外,从 BOLT-LMM v2.3.2 开始,您也可以使用`–bgenSampleFileList 选项(而不是使用 –bgenFile 和 –sampleFile )指定空格分隔的 .bgen / .sample 文件对列表选项)。此选项可以分析不同BGEN文件具有不同样本集的数据集(例如,英国生物银行v3填充发布;部分9.1[26]).
警告:BGEN 格式包含几个子格式。我们仅实现了对英国生物样本库 N=150K 和 N=500K 版本中使用的版本(和特定数据布局)的支持。特别是对于BGEN v1.2,BOLT-LMM目前仅支持用于英国生物样本库N= 500K数据的8位编码。(从 BOLT-LMM v2.3.3 开始,现在允许BGEN v1.2数据中的缺失值。)
VCF格式的填充SNP,plink格式的外显子组测序获得SNP等。BOLT-LMM 不支持上面未列出的填充数据格式,因此我们建议使用 PLINK2 将其他数据格式转换为 BGEN v1.2。如上所述,您需要使用与英国生物银行相同的BGEN v1.2子格式:
-
通过以下方式指定 8 位编码: plink2 --export bgen-1.2 bits=8
-
如果数据包含相位信息,请在转换前擦除相位。 -
如果要分析X染色体数据,则需要创建一个BGEN v1.2文件,其中所有基因型都编码为二倍体。默认情况下,plink2 会将雄性编码为单倍体,但您可以通过在转换之前将所有个体的性别设置为雌性来强制它创建二倍体 X 染色体数据。
5.1.3 X染色体分析
从 v2.3.2 开始,BOLT-LMM 接受 X 染色体基因型,用于模型拟合(通过--bfile
或--bed/bim/fam
PLINK 格式输入)和对填充变异(例如,在 BGEN 文件中)的关联测试。雄性应该被编码为二倍体(就像PLINK对染色体代码23 = X非PAR所做的那样),使得雄性基因型被编码为0/2,雌性基因型被编码为0/1/2(对应于随机X失活模型)。没有必要将chrX分为PAR和非PAR;对于PLINK输入,您应该使用PLINK --merge-x
简单地将PAR和非PAR SNP合并成单个”23号染色体”。
填充的X染色体SNP也可以包括在BOLT-LMM关联测试中;同样,男性应以当前支持的格式之一(例如,BGEN v1.1或8位BGEN v1.2)编码为二倍体。(BGEN v1.2 包含一种数据格式,该格式以本机方式对单倍体和二倍体 SNP 混合进行编码,但 BOLT-LMM 目前不支持此格式。名为23,X,XY,PAR1和PAR2的染色体都是可以接受的。
5.2 表型
可以通过以下两种方式之一指定表型:
-
–phenoUseFam:此选项告诉 BOLT-LMM 和 BOLT-REML 使用 fam 文件的最后一列(第 6 列)作为表型。此列必须是数字,因此大小写对照表型应为 1,0 编码,缺失值应用 -9 表示。 -
–phenoFile 和 –phenoCol:或者,可以在单独的空格分隔文件(用 –phenoFile 指定)中提供表型,第一行包含列标题,后续行包含记录,每个单独一个。前两列必须是 FID 和 IID(个人的 PLINK 标识符)。任何数量的列都可以跟随;包含要分析的表型的列用 –phenoCol 指定。值 -9 和 NA 被解释为缺少的数据。列中的所有其他值都应为数字。标题行后面的行中的记录不需要按排序顺序排列,并且不需要与基因型数据(即fam文件)中的个体匹配;BOLT-LMM 和 BOLT-REML 将仅分析基因型和表型文件交集中的个体,如果这些集合不匹配,将输出警告。
5.3 协变量
协变量数据可以在文件(--covarFile
)中指定,其格式与上述替代表型文件相同。(如果对表型和协变量使用相同的文件;--phenoFile 和 --covarFile
仍必须同时指定。必须使用--covarCol
(对于分类协变量)或--qCovarCol
(对于定量协变量)选项指定要使用的每个协变量。分类协变量值允许是任何不包含空格的文本字符串;列中的每个唯一文本字符串都对应于一个类别。(为了防止用户意外地使用 –covarCol 而不是 –qCovarCol 指定定量协变量,如果分类协变量包含 10 个以上的不同值,BOLT-LMM 会引发错误;此上限可以使用 –covarMaxLevels 进行修改。定量协变量值必须是数值(NA 除外)。在任一情况下,值 -9 和 NA 都会被解释为缺少的数据。如果相同类型的协变量组按顺序编号,则可以使用数组速记来指定它们(例如,对于 PC1、PC2、…、PC10 列 ,–qCovarCol=PC{1:10})。
5.4 缺失数据处理
表型缺失的个体将被忽略。默认情况下,任何缺少协变量的个体也将被忽略。这种方法通常使用,称为”完整案例分析”。作为替代方案,我们还实现了”缺失指标方法”(通过--covarUseMissingIndic
选项),该选项添加了指标变量,将缺失状态划分为附加协变量。
plink数据(--bfile
或bed/bim/fam
)中缺失的基因型被替换为每SNP平均值。填充的基因型不应包含缺失的数据;标准填充软件总是生成基因型概率估计值,即使不确定性很高。
5.5 基因型QC
BOLT-LMM 和 BOLT-REML 可自动过滤缺失率超过阈值 0.1 的 SNP 和个体。可以使用 --maxMissingPerSnp
和--maxMissingPerIndiv
修改这些阈值。请注意,过滤不是基于次要等位基因频率或与Hardy-Weinberg平衡的偏差来执行。然而,每个SNP的等位基因频率和缺失度都包含在BOLT-LMM关联测试输出中,我们建议在跟踪显着关联时检查这些值和Hardy-Weinberg p值(使PLINK--hardy
可以轻松计算)。
5.6 用户指定的筛选
要从分析中删除的单个值可以在一个或多个 --remove
文件中指定,其中列出了 FID 和 IID(每行一个单独的)。类似地,要从分析中排除的 SNP 可以在一个或多个 “–exclude`列出 SNP ID(通常为 rs 编号)的文件中指定。
请注意,--exclude
过滤不适用于填充数据;作为后处理步骤,需要单独执行特定填充 SNP 的排除。
6 关联分析
6.1 混合模型关联测试
BOLT-LMM 计算了两个关联统计量,χ2 BOLT-LMM 和 χ2 BOLT-LMM-inf,在手稿 [1[27]] 中进行了详细描述。
-
BOLT-LMM: Association test on residuals from Bayesian modeling using a mixture-of-normals prior on SNP effect sizes. 这种方法可以拟合具有中度至大效应的位点的”非无穷小”性状,从而增加关联能力。 -
BOLT-LMM-inf: Standard (infinitesimal) mixed model association. 此统计量近似于基于特征分解的软件使用的标准方法。
6.2 BOLT-LMM 混合模型关联选项
BOLT-LMM 软件为混合模型分析提供了以下选项:
-
--lmm
:执行默认 BOLT-LMM 分析,其中包括 (1a) 估计遗传性参数,(1b) 计算 BOLT-LMM-inf 统计量,(2a) 估计高斯混合参数,以及 (2b) 仅在预期功率增加时计算 BOLT-LMM 统计量。如果 BOLT-LMM 基于交叉验证确定非无穷小模型可能不会产生功率增加,则不会计算 BOLT-LMM(贝叶斯)混合模型统计量。 -
--lmmInfOnly
:仅计算无穷小混合模型关联统计量(即步骤 1a 和 1b)。 -
--lmmForceNonInf
:计算 BOLT-LMM-inf 和 BOLT-LMM 统计数据,而不管后者的功率是否预期会增加。
6.2.1 参考LD分数表
需要一个参考LD评分表[17[28]]来校准BOLT-LMM统计数据。适用于欧洲血统样本分析的参考LD分数在tables/子目录中提供,可以使用选项指定
这里补充下我的使用过程
# 下载东亚人的,竟然不需要要过“长城“
wget https://storage.googleapis.com/broad-alkesgroup-public/LDSCORE/eas_ldscores.tar.bz2
bzip2 -d eas_ldscores.tar.bz2
# 合并数据
for a in eas_ldscores/*.ldscore.gz;do zcat $a|sed -n '1!p' >> eas.tab;done
# 然后加个表头,就可以啦
--LDscoresFile=tables/LDSCORE.1000G_EUR.tab.gz
对于非欧洲数据的分析,我们建议使用LDSC软件在1000个基因组样本的祖先匹配子集上计算LD评分。
默认情况下,表中的 LD 分数通过 rsID 与 PLINK 数据中的SNP 匹配。–LDscoresMatchBp 选项允许按基对坐标匹配 SNP。
6.2.2 限制混合模型中使用的SNP
如果从填充中可以获得数百万个 SNP,我们建议在执行关联分析时一次最多包含 100 万个 SNP(使用 –modelSnps 选项)。使用最多100万个SNP的LD修剪集应实现接近最佳的功率和校正混杂,同时降低计算成本并改善收敛性。请注意,即使指定了 –modelSnps 文件,基因型数据中的所有 SNP 仍会进行关联测试;只有混合模型中的随机效应被限制为 –modelSnps。另请注意,BOLT-LMM会自动执行一条染色体(LOCO)分析,从包含被测试SNP的染色体中省略SNP,以避免近端污染[4[29],9[30]]。
6.3 标准线性回归
设置--verboseStats
标志将在其他输出列中输出标准线性回归卡方统计数据和 p 值,CHISQ_LINREG和P_LINREG。请注意,与混合模型关联不同,线性回归容易受到总体分层的影响,因此在执行线性回归时,您可能希望将主成分(使用其他软件计算,例如,EIGENSOFT v6.0+ 中的 PLINK2 或 FastPCA [18[31]] )作为协变量。将PC作为协变量包括在内也将加速BOLT-LMM混合模型计算的收敛。
7 方差成分分析 (BOLT-REML)
使用 --reml
选项调用 BOLT-REML 算法来估计遗传性参数和遗传相关性。
7.1 多个方差成分
要将 SNP 分配给不同的方差分量,请指定一个--modelSnps
文件,其中每条以空格分隔的行都包含一个 SNP ID(通常为 rs 编号),后跟它所属的差异分量的名称。
7.2 多种特征
要执行多性状方差分量分析,请指定多个--phenoCol
参数值标志(对应于同一 --phenoFile
中的不同列)。BOLT-REML目前仅支持对一组个体的性状表型进行多性状分析,因此任何缺少至少一种表型的个体都将被忽略。对于 D 性状,BOLT-REML 估计每个方差分量的 D 遗传性参数和每个方差分量(包括残差方差分量)的 D(D-1)/2 相关性。
7.3 初始方差参数猜测
要指定一组用于启动 REML 迭代的方差参数(与 BOLT-REML 使用的默认过程相比,如果您有良好的初始猜测,这可能会节省时间),请使用 --remlGuessStr="string"
并采用以下格式。对于每个方差分量(从残差项开始,该项自动命名为 env/noise),请指定方差分量的名称,后跟初始猜测。例如,具有两个名为 vc1 和 vc2 的(非残差)方差分量(在 –modelSnps 文件中)的模型可能具有由以下指定方差参数猜测:
--remlGuessStr="env/noise 0.5 vc1 0.2 vc2 0.3"
请注意,估计值的总和必须等于 1;BOLT-REML将自动相应地使表型正常化。
对于 D 性状的多性状分析,--remlGuessStr
需要指定 D 方差比例的猜测和每个方差分量的 D(D-1)/2 成对相关性。将这些值视为上三角形矩阵的条目(对角线上有方差比例,对角线上有相关性),您应该在每个方差分量名称后指定这些 D(D+1)/2 值,方法是从左到右、从上到下读取它们。
7.4 损失一点准确性的速度
BOLT-REML 使用蒙特卡罗算法来提高 REML 优化速度 [12[32]]。默认情况下,BOLT-REML 使用 15 个蒙特卡罗试验执行初始优化,然后使用 100 个蒙特卡罗试验优化参数估计值。如果计算成本是一个问题(或执行探索性分析),则可以使用--remlNoRefine
标志(除了--reml
标志之外)跳过优化步骤。此选项通常提供 2–3 倍的加速,代价是标准误差高出约 1.03 倍。
8 输出
8.1 BOLT-LMM关联测试统计
BOLT-LMM 关联统计信息以制表符分隔的 --statsFile
文件输出,其中包含以下字段,每个 SNP 一行:
-
SNP:rs 编号或 ID 字符串 -
CHR:染色体 -
BP:物理(碱基对)位置 -
GENPOS:来自bim文件或从遗传图谱插值的遗传位置 -
等位基因1:bim文件中的第一个等位基因(通常是次要等位基因),用作效果等位基因 -
等位基因0:bim文件中的第二个等位基因,用作参考等位基因 -
A1FREQ:第一个等位基因的频率 -
F_MISS:该SNP基因型缺失的个体的比例 -
BETA:从 BOLT-LMM 近似到无穷小混合模型的效应大小 -
SE:效应大小的标准误差 -
P_BOLT_LMM_INF:无穷小混合模型关联检验 p 值 -
P_BOLT_LMM:非无穷小混合模型关联检验 p 值
可选附加输出。要输出所有关联检验的卡方统计数据,请设置--verboseStats
标志。
8.2 BOLT-REML 输出和日志记录
当分析完成时,BOLT-REML输出(即方差参数估计值和标准误差)只需打印到终端(标准输出)即可。随着分析的进行,BOLT-LMM 和 BOLT-REML 都将输出写入(stdout 和 stderr);我们建议保存此输出。如果您希望在命令行上同时查看此输出的同时保存它,则可以使用
./bolt [... list of options ...] 2>&1 | tee output.log
9.分析N=50万英国生物样本库数据的建议
BOLT-LMM的许多用户希望分析英国生物样本库数据。以下是计算N = 500K英国生物样本库样本的关联统计数据的一些提示(另见参考文献[10[33]]):
-
计算单个表型的全基因组关联测试统计数据通常需要几天到一周的多线程;我们建议使用** 8 个以上的线程**。 -
计算将需要高达100GB的内存,具体取决于模型中包含的直接基因分型SNP的数量(带有 --bed/--bim/--fam
)。如果计算成本是一个问题,则可以通过指定要在模型中使用 –modelSnps 的 SNP 子集(例如,通过对 MAF 或缺失进行过滤或通过 LD 修剪)来减少运行时间和 RAM。 -
还可以通过将主成分作为协变量包括在内来加快分析的模型拟合步骤(因为投影出顶部特征向量可以改善核矩阵的条件反射)。请注意,为了实现这种加速,主成分必须使用模型中使用的样本和SNP集进行计算(例如,使用 plink2 --pca approx
或EIGENSOFT v6.0 +快速模式)。在不同样本集和/或SNP(例如,英国生物银行提供的样本)上计算的预先计算的主成分往往不能提供太多的加速。 -
模型拟合后,BOLT-LMM 计算填充 SNP 的关联统计信息。对于 BGEN v1.2 数据,此计算现在是多线程的,并且应该足够快,可以将所有染色体包含在单个作业中,但是跨作业并行分析染色体子集当然也是允许的。设置` –bgenMinMAF/–bgenMinINFO 将减小输出文件大小并提高速度。 -
您将需要创建一个版本 --fam
文件,该文件的第 6 列中包含数值,并且还需要--remove
plink在数据中但不在填充数据中的个人。(如果不匹配,BOLT-LMM 将生成要--remove
的此类样本的列表。 -
下面是一个示例命令行:
> ./bolt '
> --bed=ukb_cal_chr{1:22}_v2.bed '
> --bim=ukb_snp_chr{1:22}_v2.bim '
> --fam=ukb1404_cal_chr1_v2_CURRENT.fixCol6.fam '
> --删除 =bolt.in_plink_but_not_imputed。FID_IID.976.txt '
> --remove=sampleQC/remove.nonWhite.FID_IID.txt '
> --exclude=snpQC/autosome_maf_lt_0.001.txt '
> --exclude=snpQC/autosome_missing_gt_0.1.txt '
> --phenoFile=ukb4777.phenotypes.tab '
> --phenoCol=height '
> --covarFile=ukb4777.covars.tab.gz '
> --covarCol=cov_ASSESS_CENTER '
> --covarCol=cov_GENO_ARRAY '
> --covarMaxLevels=30 '
> --qCovarCol=cov_AGE '
> --qCovarCol=cov_AGE_SQ '
> --qCovarCol=PC{1:20} '
> --LDscoresFile=tables/LDSCORE.1000G_EUR.tab.gz '
> --geneticMapFile=tables/genetic_map_hg19.txt.gz '
> --lmmForceNonInf '
> --numThreads=8 '
> --statsFile=bolt_460K_selfRepWhite.height.stats.gz '
> --bgenFile=ukb_imp_chr{1:22}_v2.bgen '
> --bgenMinMAF=1e-3 '
> --bgenMinINFO=0.3 '
> --sampleFile=ukb1404_imp_chr1_v2_s487406.sample '
> --statsFileBgenSnps=bolt_460K_selfRepWhite.height.bgen.stats.gz '
> --verboseStats
9.1 英国生物样本库 v3 填充发布
英国生物银行v3填充版本(2018年3月7日)中的BGEN文件具有与以前的v2版本相同的格式,但现在包括染色体X和XY(= PAR1 + PAR2)的文件。这些文件以与其他文件相同的BGEN子格式(8位编码,男性编码为二倍体)编码;但是,它们包含的样本略少于常染色体文件(对应于in。Phasing.Input.chrX and in.Phasing.Input.chrXY 字段的示例 QC 文件)。
如果您希望同时分析常染色体和X染色体数据,您可以执行以下任一操作:
-
通过使用 --bgenSampleFileList
选项指定以空格分隔的.bgen / .sample
文件对列表,一次对所有文件运行 BOLT-LMM。请注意,您需要--remove
任何 BGEN 文件中不存在的所有样本。(如果 BOLT-LMM 检测到缺少的样本,它将报告错误并编写此类样本的列表供您--remove
。 -
在两个单独的 BOLT-LMM 运行中分析常染色体和 chrX 变异(使用两次运行中的所有常染色体和 chrX 类型变异作为模型拟合的 PLINK 输入)。
第一种方法比第二种方法稍微方便一些,但代价是与第二种方法相比,去除了稍微多一点的样本(约1000个额外的样本)。
10 病例对照性状的关联分析
虽然BOLT-LMM的数学推导基于定量性状模型,但BOLT-LMM也可以应用于分析病例控制性状(只需将二元性状视为定量性状)。然而,需要注意的一个重要警告是,正如SAIGE论文[19[34]]中所述,BOLT-LMM测试统计数据可能会因不平衡的病例对照性状而变得校准错误(导致罕见SNP的假阳性关联)。
10.1 病例对照平衡指南
BOLT-LMM P值在多大程度上可能遭受二元性状的误校准是三个变量的函数:样本量,次要等位基因频率和病例对照比例。具体而言,当次要等位基因计数(MAC)乘以病例分数相对较小时,就会发生校准错误(对应于传统观点,即卡方检验统计量在预期计数较小时会分解)。(对于任何其他基于回归的卡方统计量中的 P 值也是如此。我们还注意到,如果数量性状具有极端异常值,则在分析数量性状时可能会出现类似的问题。
在我们的预印本的修订版中,我们探索了BOLT-LMM在英国生物样本库N= 500K数据上的性能[10[35]],我们包括了一套模拟,这些模拟改变了影响I型误差控制的三个关键参数(样本大小,次要等位基因频率和病例分数)。对于完整的英国生物样本库数据的分析,我们确定对于病例占比至少为10%的性状,BOLT-LMM测试统计数据对于具有MAF>0.1%的 SNP进行了良好校准。如果最小MAF增加,也可以容忍更极端的病例控制不平衡。这些模拟的完整结果在参考文献[10[36]]的补充表8中提供,我们建议咨询该表以确定BOLT-LMM是否适合特定的二元性状分析。
对于不适合BOLT-LMM分析的高度不平衡的病例控制设置,我们建议使用SAIGE [19[37]](通过使用鞍点近似克服了偏离渐近正态性的问题)。
10.2 优势比的估计
由于BOLT-LMM在分析病例控制特征时仍然使用线性模型(而不是逻辑模型),因此需要进行转换,以便将定量尺度上的SNP效应大小估计值(”betas”)转换为传统的比值比。一个合理的近似值是:
对数 OR = β / (μ * (1 – μ))),其中 μ = 大小写分数。
在应用上述变换以获得对数优势比时,SNP效应大小估计值的标准误差也应除以(μ * (1 – μ))。
或者,下面描述了更复杂的转换:http://cnsgenomics.com/shiny/LMOR/[38]
11 常见问题
用户最常问的问题是,当 BOLT-LMM 报告由接近 0 或 1 的遗传性估计值引起的错误时,该怎么办。旧版本的BOLT-LMM报告”错误:遗传性估计无效;无法继续分析”;较新版本试图澄清问题:
-
“错误:遗传力估计接近0;LMM 可能无法纠正混杂。相反,对不相关项使用 PC 校正的线性/逻辑回归。
当遗传性估计值达到0时,线性混合模型关联检验(包括BOLT-LMM和其他方法)全部退化为简单线性回归,因此出现错误消息。这种情况是危险的,因为**”混合模型”将不再纠正人群分层和相关性**。
当然,在这些情况下,您仍然可以运行关联检验,但您需要修剪到一组不相关的样本(如果您的样本集包含相关个体)并包括主成分协变量。
您可以使用 BOLT-LMM 执行线性回归,方法是在不带
--lmm
选项 (以及--verboseStats
选项)的情况下运行线性回归。 -
“错误:遗传力估计接近1;算法可能无法收敛。 由于样本量少或病例确定性低,分析可能不合适。
当样本数量较低时,最常出现此误差,导致估计的遗传性具有非常大的标准误差(甚至可能大于1),使得估计值可能在0到1范围内的任何位置,并可能达到其中一个边界。不建议将 BOLT-LMM 用于分析较小的样品;在这种情况下,我们建议尝试其他软件包,如GEMMA或GCTA。
12 网站和联系信息
软件更新将在此处和以下网站上发布:
http://www.hsph.harvard.edu/alkes-price/software/[39]
如果您对 BOLT-LMM 软件有任何意见或疑问,请联系 Po-Ru Loh,poruloh@broadinstitute.org。
13 许可证
BOLT-LMM 是 GNU 通用公共许可证 v3.0 (GPLv3) 下的自由软件。
参考资料
10: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#x1-5300010
[2][1: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xloh2015ng
[3][2: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xkang2010ng
[4][3: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xlippert2011nmeth
[5]4: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xlistgarten2012nmeth
[6]5: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xlistgarten2013ng
[7]6: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xlippert2013scirep
[8][7: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xzhou2012ng
[9][8: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xsvishcheva2012ng
[10][9: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xyang2014ng
[11][1: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xloh2015ng
[12][10: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xloh2018ng
[13][11]一: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xyang2011ajhg
[14][12: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xbolt_reml
[15]http://data.broadinstitute.org/alkesgroup/BOLT-LMM/downloads/: http://data.broadinstitute.org/alkesgroup/BOLT-LMM/downloads/
[16]9.1: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#x1-520009.1
[17]Anaconda.org: https://anaconda.org/
[18][13: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#XjohnsonNLopt
[19]http://remidaviet.com/software.php: http://remidaviet.com/software.php
[20][1: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xloh2015ng
[21]10: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xloh2018ng
[22][14: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xpurcell2007ajhg
[23][15: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xchang2015gigascience
[24][8: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xsvishcheva2012ng
[25][16: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xhowie2009pg
[26]9.1: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#x1-520009.1
[27][1: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xloh2015ng
[28][17: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xbuliksullivan2015ng
[29][4: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xlistgarten2012nmeth
[30]9: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xyang2014ng
[31][18: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xgalinsky2016ajhg
[32][12: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xbolt_reml
[33][10: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xloh2018ng
[34][19: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xzhou2018ng
[35][10: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xloh2018ng
[36][10: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xloh2018ng
[37][19: https://storage.googleapis.com/broad-alkesgroup-public/BOLT-LMM/BOLT-LMM_manual.html#Xzhou2018ng
[38]http://cnsgenomics.com/shiny/LMOR/: http://cnsgenomics.com/shiny/LMOR/
[39]http://www.hsph.harvard.edu/alkes-price/software/: http://www.hsph.harvard.edu/alkes-price/software/
本篇文章来源于微信公众号: 微因