BOLT-LMM v2.3 User Manual

Po-Ru Loh

August 1, 2017

Contents

1 Overview

The BOLT-LMM software package currently consists of two main algorithms, the BOLT-LMM algorithm for mixed model association testing, and the BOLT-REML algorithm for variance components analysis (i.e., partitioning of SNP-heritability and estimation of genetic correlations).

We recommend BOLT-LMM for analyses of human genetic data sets containing more than 5,000 samples. The algorithms used in BOLT-LMM rely on approximations that hold only at large sample sizes and have been tested only in human data sets. For analyses of fewer than 5,000 samples, we recommend the GCTA software.

1.1 BOLT-LMM mixed model association testing

The BOLT-LMM algorithm computes statistics for testing association between phenotype and genotypes using a linear mixed model (LMM)  [1]. By default, BOLT-LMM assumes a Bayesian mixture-of-normals prior for the random effect attributed to SNPs other than the one being tested. This model generalizes the standard “infinitesimal” mixed model used by existing mixed model association methods (e.g., EMMAX  [2], FaST-LMM  [3456], GEMMA  [7], GRAMMAR-Gamma  [8], GCTA-LOCO  [9]), providing an opportunity for increased power to detect associations while controlling false positives. Additionally, BOLT-LMM applies algorithmic advances to compute mixed model association statistics much faster than existing methods, both when using the Bayesian mixture model and when specialized to standard mixed model association. BOLT-LMM is described in ref.  [1]:

Loh P-R, Tucker G, Bulik-Sullivan BK, Vilhjįlmsson BJ, Finucane HK, Salem RM, Chasman DI, Ridker PM, Neale BM, Berger B, Patterson N, and Price AL. Efficient Bayesian mixed model analysis increases association power in large cohorts. Nature Genetics, 2015.

1.2 BOLT-REML variance components analysis

The BOLT-REML algorithm estimates heritability explained by genotyped SNPs and genetic correlations among multiple traits measured on the same set of individuals. Like the GCTA software  [10], BOLT-REML applies variance components analysis to perform these tasks, supporting both multi-component modeling to partition SNP-heritability and multi-trait modeling to estimate correlations. BOLT-REML applies a Monte Carlo algorithm that is much faster than standard methods for variance components analysis (e.g., GCTA) at large sample sizes. BOLT-REML is described in ref.  [11]:

Loh P-R, Bhatia G, Gusev A, Finucane HK, Bulik-Sullivan BK, Pollack SJ, PGC-SCZ Working Group, de Candia TR, Lee SH, Wray NR, Kendler KS, O’Donovan MC, Neale BM, Patterson N, and Price AL. Contrasting genetic architectures of schizophrenia and other complex diseases using fast variance components analysis. Nature Genetics, 2015.

2 Download and installation

You can download the latest version of the BOLT-LMM software at:

http://data.broadinstitute.org/alkesgroup/BOLT-LMM/downloads/

Previous versions are also available in the old/ subdirectory.

2.1 Change log

2.2 Installation

The BOLT-LMM_vX.X.tar.gz download package contains a standalone (i.e., statically linked) 64-bit Linux executable, bolt, which we have tested on several Linux systems. We strongly recommend using this static executable because it is well-optimized and no further installation is required.

If you wish to compile your own version of the BOLT-LMM software from the source code (in the src/ subdirectory, provided in a separate download licensed under GNU GPLv3), you will need to ensure that library dependencies are fulfilled and will need to make appropriate modifications to the Makefile:

For reference, the provided bolt executable was created on the Harvard Medical School “Orchestra” research computing cluster using Intel Parallel Studio XE 2016 Composer (MKL 11.3) and the Boost 1.58 and NLopt 2.4.2 libraries by invoking make linking=static-except-glibc.

2.3 Running BOLT-LMM and BOLT-REML

To run the bolt executable, simply invoke ./bolt on the Linux command line (within the BOLT-LMM install directory) with parameters in the format --optionName=optionValue.

2.4 Examples

The example/ subdirectory contains a bash script run_example.sh that demonstrates basic use of BOLT-LMM on a small example data set. Likewise, run_example_reml2.sh demonstrates BOLT-REML.

2.5 Help

To get a list of basic options, run:

./bolt -h

To get a complete list of basic and advanced options, run:

./bolt --helpFull

3 Computing requirements

3.1 Operating system

At the current time we have only compiled and tested BOLT-LMM on Linux computing environments; however, the source code is available if you wish to try compiling BOLT-LMM for a different operating system.

3.2 Memory

For typical data sets (M, N exceeding 10,000), BOLT-LMM and BOLT-REML use approximately MN/4 bytes of memory, where M is the number of SNPs and N is the number of individuals. More precisely:

3.3 Running time

In practice, BOLT-LMM and BOLT-REML have running times that scale roughly with MN1.5. Our largest analyses of real data (M = 600K SNPs, N = 60K individuals) took ~1 day using a single computational core. We have also tested BOLT-LMM on simulated data sets containing up to N = 480K individuals; for more details, please see the BOLT-LMM manuscript  [1].

3.3.1 Multi-threading

On multi-core machines, running time can be reduced by invoking multi-threading using the --numThreads option.

4 Input/output file naming conventions

4.1 Automatic gzip [de]compression

The BOLT-LMM software assumes that input files ending in .gz are gzip-compressed and automatically decompresses them on-the-fly (i.e., without creating a temporary file). Similarly, BOLT-LMM writes gzip-compressed output to any output file ending in .gz.

4.2 Arrays of input files and covariates

Arrays of sequentially-numbered input files and covariates can be specified by the shorthand {i:j}. For example,

data.chr{1:22}.bim

is interpreted as the list of files

data.chr1.bim, data.chr2.bim, ..., data.chr22.bim

5 Input

5.1 Genotypes

The BOLT-LMM software takes genotype input in PLINK  [13] binary format (bed/bim/fam). For file conversion and data manipulation in general, we highly recommend the PLINK/PLINK2 software  [14].

If all genotypes are contained in a single bed/bim/fam file triple with the same file prefix, you may simply use the command line option --bfile=prefix. Genotypes may also be split into multiple bed and bim files containing consecutive sets of SNPs (e.g., one bed/bim file pair per chromosome) either by using multiple --bed and --bim invocations or by using the file array shorthand described above (e.g., --bim=data.chr{1:22}.bim).

5.1.1 Reference genetic maps

The BOLT-LMM package includes reference maps that you can use to interpolate genetic map coordinates from SNP physical (base pair) positions in the event that your PLINK bim file does not contain genetic coordinates (in units of morgans). (The BOLT-LMM association testing algorithm uses genetic positions to prevent proximal contamination; BOLT-REML does not use this information.) To use a reference map, use the option

--geneticMapFile=tables/genetic_map_hg##.txt.gz

selecting the build (hg17, hg18, or hg19) corresponding to the physical coordinates of your bim file. You may use the --geneticMapFile option even if your PLINK bim file does contain genetic coordinates; in this case, the genetic coordinates in the bim file will be ignored, and interpolated coordinates will be used instead.

5.1.2 Imputed SNP dosages

As of version 1.1, the BOLT-LMM association testing algorithm supports computation of mixed model association statistics at an arbitrary number of imputed SNPs (with real-valued “dosages” rather than hard-called genotypes) using a mixed model built on a subset of hard-called genotypes. (BOLT-REML variance components analysis does not support dosage input.) This approach requires only a trivial amount of additional computation and no additional RAM, as it simply applies a genome scan (as in GRAMMAR-Gamma  [8]) of real-valued dosage SNPs against the residual phenotypes that BOLT-LMM already computes. We expect that using a mixed model built on only a subset of ~500K hard-called genotypes should sacrifice almost no power while retaining the computational efficiency of BOLT-LMM.

If you have only imputed SNP data on hand, you will need to pre-process your data set to create a subset of hard-called SNPs in PLINK format for BOLT-LMM. We suggest the following procedure.

  1. Determine a high-confidence set of SNPs (e.g., based on IMPUTE2 INFO score) at which to create an initial hard-call set.
  2. Create hard-called genotypes at these SNPs in PLINK format.
  3. Use PLINK to LD prune to ~500K SNPs (via --indep-pairwise 50 5 r2thresh for an appropriate r2thresh).
  4. Run BOLT-LMM using the final hard-called SNPs as the --bfile (or bed/bim/fam) argument, specifying the imputed SNPs as additional association test SNPs using one of the formats below.

Imputed SNPs in dosage format. This input format consists of one or more --dosageFile parameters specifying files that contain real-valued genotype expectations at imputed SNPs. Each line of a --dosageFile should be formatted as follows:

rsID   chr   pos   allele1   allele0   [dosage = E[#allele1]] x N

Missing (i.e., uncalled) dosages can be specified with –9. You will also need to provide one additional --dosageFidIidFile specifying the PLINK FIDs and IIDs of samples that the dosages correspond to. See the example/ subdirectory for an example.

Imputed SNPs in IMPUTE2 format. You may also specify imputed SNPs as output by the IMPUTE2 software  [15]. The IMPUTE2 genotype file format is as follows:

snpID   rsID   pos   allele1   allele0   [p(11) p(10) p(00)] x N

(BOLT-LMM ignores the snpID field.) Here, instead of dosages, each genotype entry contains individual probabilities of the individual being homozygous for allele1, heterozygous, and homozygous for allele0. The three probabilities need not sum to 1, allowing for genotype uncertainty; if the sum of the probabilities is less than the --impute2CallThresh parameter, BOLT-LMM treats the genotype as missing.

To compute association statistics at a list of files containing IMPUTE2 SNPs, you may list the files within a --impute2FileList file. Each line of this file should contain two entries: a chromosome number followed by an IMPUTE2 genotype file containing SNPs from that chromosome. You will also need to provide one additional --impute2FidIidFile specifying the PLINK FIDs and IIDs of samples that the IMPUTE2 genotypes correspond to. See the example/ subdirectory for an example.

Imputed SNPs in 2-dosage format. You may also specify imputed SNPs as output by the Ricopili pipeline and plink2 --dosage format=2. This file format consists of file pairs: (1) PLINK map files containing information about SNP locations; and (2) genotype probability files in the 2-dosage format, which consists of a header line

SNP   A1   A2   [FID IID] x N

followed by one line per SNP in the format

rsID   allele1   allele0   [p(11) p(10)] x N

The third genotype probability for each entry is assumed to be p(00)=1-p(11)-p(10) (unlike with the IMPUTE2 format).

To compute association statistics at SNPs in a list of 2-dosage files, you may list the files within a --dosage2FileList file. Each line of this file should contain two entries: a PLINK map file followed by the corresponding genotype file containing probabilities for those SNPs. (As usual, if either file ends with .gz, it is automatically unzipped; otherwise it is assumed to be plain text.) See the example/ subdirectory for an example.

Imputed SNPs in BGEN format. To compute association statistics at SNPs in one or more BGEN data files, specify the .bgen file(s) with --bgenFile and the corresponding .sample file with --sampleFile. The --bgenMinMAF and --bgenMinINFO options allows limiting output to SNPs passing minimum allele frequency and INFO thresholds.

Note that starting with BOLT-LMM v2.3, the --bgenFile argument allows multiple BGEN files. We have implemented multi-threaded processing for files in the BGEN v1.2 format used in the UK Biobank N=500K release, so analyzing BGEN v1.2 data for all chromosomes within a single job should be feasible. For analyses of BGEN v1.1 data used in the N=150K release, we recommend parallelizing across chromosomes for computational convenience (using the full --bfile of directly genotyped PLINK data from all chromosomes in each job).

WARNING: The BGEN format comprises a few sub-formats; we have only implemented support for the versions (and specific data layouts) used in the UK Biobank N=150K and N=500K releases.

5.2 Phenotypes

Phenotypes may be specified in either of two ways:

5.3 Covariates

Covariate data may be specified in a file (--covarFile) with the same format as the alternate phenotype file described above. (The same file may be used for both phenotypes and covariates.) Each covariate to be used must be specified using either a --covarCol (for categorical covariates) or a --qCovarCol (for quantitative covariates) option. Categorical covariate values are allowed to be any text strings not containing whitespace; each unique text string in a column corresponds to a category. (To guard against users accidentally specifying quantitative covariates with --covarCol instead of --qCovarCol, BOLT-LMM throws an error if a categorical covariate contains more than 10 distinct values; this upper bound can be modified with --covarMaxLevels.) Quantitative covariate values must be numeric (with the exception of NA). In either case, values of -9 and -NA are interpreted as missing data. If groups of covariates of the same type are numbered sequentially, they may be specified using array shorthand (e.g., --qCovarCol=PC{1:10} for columns PC1, PC2, ..., PC10).

5.4 Missing data treatment

Individuals with missing phenotypes are ignored. By default, individuals with any missing covariates are also ignored; this approach is commonly used and referred to as “complete case analysis.” As an alternative, we have also implemented the “missing indicator method” (via the --covarUseMissingIndic option), which adds indicator variables demarcating missing status as additional covariates. Missing genotypes are replaced with per-SNP averages.

5.5 Genotype QC

BOLT-LMM and BOLT-REML automatically filter SNPs and individuals with missing rates exceeding thresholds of 0.1. These thresholds may be modified using --maxMissingPerSnp and --maxMissingPerIndiv. Note that filtering is not performed based on minor allele frequency or deviation from Hardy-Weinberg equilibrium. Allele frequency and missingness of each SNP are included in the BOLT-LMM association test output, however, and we recommend checking these values and Hardy-Weinberg p-values (which are easily computed using PLINK --hardy) when following up on significant associations.

5.6 User-specified filters

Individuals to remove from the analysis may be specified in one or more --remove files listing FID and IIDs (one individual per line). Similarly, SNPs to exclude from the analysis may be specified in one ore more --exclude files listing SNP IDs (typically rs numbers).

6 Association analysis (BOLT-LMM)

6.1 Mixed model association tests

BOLT-LMM computes two association statistics, χ2 BOLT-LMM and χ2 BOLT-LMM-inf, described in detail in our manuscript  [1].

6.2 BOLT-LMM mixed model association options

The BOLT-LMM software offers the following options for mixed model analysis:

6.2.1 Reference LD score tables

A table of reference LD scores  [16] is needed to calibrate the BOLT-LMM statistic. Reference LD scores appropriate for analyses of European-ancestry samples are provided in the tables/ subdirectory and can be specified using the option

--LDscoresFile=tables/LDSCORE.1000G_EUR.tab.gz

For analyses of non-European data, we recommend computing LD scores using the LDSC software on an ancestry-matched subset of the 1000 Genomes samples.

By default, LD scores in the table are matched to SNPs in the PLINK data by rsID. The --LDscoresMatchBp option allows matching SNPs by base pair coordinate.

6.2.2 Restricting SNPs used in the mixed model

If millions of SNPs are available from imputation, we suggest including at most 1 million SNPs at a time in the mixed model (using the --modelSnps option) when performing association analysis. Using an LD pruned set of at most 1 million SNPs should achieve near-optimal power and correction for confounding while reducing computational cost and improving convergence. Note that even when a file of --modelSnps is specified, all SNPs in the genotype data are still tested for association; only the random effects in the mixed model are restricted to the --modelSnps. Also note that BOLT-LMM automatically performs leave-one-chromosome-out (LOCO) analysis, leaving out SNPs from the chromosome containing the SNP being tested in order to avoid proximal contamination  [49].

6.3 Standard linear regression

Setting the --verboseStats flag will output standard linear regression chi-square statistics and p-values in additional output columns CHISQ_LINREG and P_LINREG. Note that unlike mixed model association, linear regression is susceptible to population stratification, so you may wish to include principal components (computed using other software, e.g., PLINK2 or FastPCA  [17] in EIGENSOFT v6.0+) as covariates when performing linear regression. Including PCs as covariates will also speed up convergence of BOLT-LMM’s mixed model computations.

7 Variance components analysis (BOLT-REML)

Using the --reml option invokes the BOLT-REML algorithm for estimating heritability parameters and genetic correlations.

7.1 Multiple variance components

To assign SNPs to different variance components, specify a --modelSnps file in which each whitespace-delimited line contains a SNP ID (typically an rs number) followed by the name of the variance component to which it belongs.

7.2 Multiple traits

To perform multi-trait variance components analysis, specify multiple --phenoCol parameter-value flags (corresponding to different columns in the same --phenoFile). BOLT-REML currently only supports multi-trait analysis of traits phenotyped on a single set of individuals, so any individuals with at least one missing phenotype will be ignored. For D traits, BOLT-REML estimates D heritability parameters per variance component and D(D-1)/2 correlations per variance component (including the residual variance component).

7.3 Initial variance parameter guesses

To specify a set of variance parameters at which to start REML iteration (which may save time compared to the default procedure used by BOLT-REML if you have good initial guesses), use --remlGuessStr="string" with the following format. For each variance component (starting with the residual term, which is automatically named env/noise), specify the name of the variance component followed by the initial guess. For instance, a model with two (non-residual) variance components named vc1 and vc2 (in the --modelSnps file) could have variance parameter guesses specified by:

--remlGuessStr="env/noise 0.5 vc1 0.2 vc2 0.3"

Note that the sum of the estimates must equal 1; BOLT-REML will automatically normalize the phenotype accordingly.

For multi-trait analysis of D traits, the --remlGuessStr needs to specify both guesses of D variance proportions and D(D-1)/2 pairwise correlations per variance component. Viewing these values as entries of an upper-triangular matrix (with variance proportions on the diagonal and correlations above the diagonal), you should specify these D(D+1)/2 values after each variance component name by reading them off left-to-right, top-to-bottom.

7.4 Trading a little accuracy for speed

BOLT-REML uses a Monte Carlo algorithm to increase REML optimization speed  [11]. By default, BOLT-REML performs an initial optimization using 15 Monte Carlo trials and then refines parameter estimates using 100 Monte Carlo trials. If computational cost is a concern (or to perform exploratory analyses), you can skip the refinement step using the --remlNoRefine flag (in addition to the --reml flag). This option typically gives 2–3x speedup at the cost of ~1.03x higher standard errors.

8 Output

8.1 BOLT-LMM association test statistics

BOLT-LMM association statistics are output in a tab-delimited --statsFile file with the following fields, one line per SNP:

Optional additional output. To output chi-square statistics for all association tests, set the --verboseStats flag.

8.2 BOLT-REML output and logging

BOLT-REML output (i.e., variance parameter estimates and standard errors) is simply printed to the terminal (stdout) when analysis finishes. Both BOLT-LMM and BOLT-REML write output to (stdout and stderr) as analysis proceeds; we recommend saving this output. If you wish to save this output while simultaneously viewing it on the command line, you may do so using

./bolt [... list of options ...] 2>&1 | tee output.log

9 Recommendations for analyzing N=500K UK Biobank data

Many users of BOLT-LMM wish to analyze UK Biobank data. Here are a few tips for computing association statistics on N=500K UK Biobank samples:

10 Website and contact info

Software updates will be posted at the following permanent website:

http://www.hsph.harvard.edu/alkes-price/software/

If you have comments or questions about the BOLT-LMM software, please contact Po-Ru Loh, loh@hsph.harvard.edu.

11 License

BOLT-LMM is free software under the GNU General Public License v3.0 (GPLv3).

References

1.   Loh, P.-R. et al. Efficient Bayesian mixed model analysis increases association power in large cohorts. Nature Genetics 47, 284–290 (2015).

2.   Kang, H. M. et al. Variance component model to account for sample structure in genome-wide association studies. Nature Genetics 42, 348–354 (2010).

3.   Lippert, C. et al. FaST linear mixed models for genome-wide association studies. Nature Methods 8, 833–835 (2011).

4.   Listgarten, J. et al. Improved linear mixed models for genome-wide association studies. Nature Methods 9, 525–526 (2012).

5.   Listgarten, J., Lippert, C. & Heckerman, D. FaST-LMM-Select for addressing confounding from spatial structure and rare variants. Nature Genetics 45, 470–471 (2013).

6.   Lippert, C. et al. The benefits of selecting phenotype-specific variants for applications of mixed models in genomics. Scientific Reports 3 (2013).

7.   Zhou, X. & Stephens, M. Genome-wide efficient mixed-model analysis for association studies. Nature Genetics 44, 821–824 (2012).

8.   Svishcheva, G. R., Axenovich, T. I., Belonogova, N. M., van Duijn, C. M. & Aulchenko, Y. S. Rapid variance components-based method for whole-genome association analysis. Nature Genetics 44, 1166–1170 (2012).

9.   Yang, J., Zaitlen, N. A., Goddard, M. E., Visscher, P. M. & Price, A. L. Advantages and pitfalls in the application of mixed-model association methods. Nature Genetics 46, 100–106 (2014).

10.   Yang, J., Lee, S. H., Goddard, M. E. & Visscher, P. M. GCTA: a tool for genome-wide complex trait analysis. American Journal of Human Genetics 88, 76–82 (2011).

11.   Loh, P.-R. et al. Contrasting genetic architectures of schizophrenia and other complex diseases using fast variance components analysis. Nature Genetics 47, 1385–1392 (2015).

12.   Johnson, S. G. The NLopt nonlinear-optimization package. URL http://ab-initio.mit.edu/nlopt.

13.   Purcell, S. et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. American Journal of Human Genetics 81, 559–575 (2007).

14.   Chang, C. C. et al. Second-generation PLINK: rising to the challenge of larger and richer datasets. GigaScience 4, 1–16 (2015).

15.   Howie, B. N., Donnelly, P. & Marchini, J. A flexible and accurate genotype imputation method for the next generation of genome-wide association studies. PLOS Genetics 5, e1000529 (2009).

16.   Bulik-Sullivan, B. et al. LD Score regression distinguishes confounding from polygenicity in genome-wide association studies. Nature Genetics 47, 291–295 (2015).

17.   Galinsky, K. J. et al. Fast principal-component analysis reveals convergent evolution of ADH1B in Europe and East Asia. American Journal of Human Genetics 98, 456–472 (2016).