 Technical note
 Open Access
 Open Peer Review
 Published:
Secondgeneration PLINK: rising to the challenge of larger and richer datasets
GigaScience volume 4, Article number: 7 (2015)
Abstract
Background
PLINK 1 is a widely used opensource C/C++ toolset for genomewide association studies (GWAS) and research in population genetics. However, the steady accumulation of data from imputation and wholegenome sequencing studies has exposed a strong need for faster and scalable implementations of key functions, such as logistic regression, linkage disequilibrium estimation, and genomic distance evaluation. In addition, GWAS and populationgenetic data now frequently contain genotype likelihoods, phase information, and/or multiallelic variants, none of which can be represented by PLINK 1’s primary data format.
Findings
To address these issues, we are developing a secondgeneration codebase for PLINK. The first major release from this codebase, PLINK 1.9, introduces extensive use of bitlevel parallelism, \(O\left (\sqrt {n}\right)\)time/constantspace HardyWeinberg equilibrium and Fisher’s exact tests, and many other algorithmic improvements. In combination, these changes accelerate most operations by 14 orders of magnitude, and allow the program to handle datasets too large to fit in RAM. We have also developed an extension to the data format which adds lowoverhead support for genotype likelihoods, phase, multiallelic variants, and reference vs. alternate alleles, which is the basis of our planned second release (PLINK 2.0).
Conclusions
The secondgeneration versions of PLINK will offer dramatic improvements in performance and compatibility. For the first time, users without access to highend computing resources can perform several essential analyses of the featurerich and very large genetic datasets coming into use.
Findings
Because of its broad functionality and efficient binary file format, PLINK is widely employed in dataprocessing pipelines that are established for genetrait mapping and populationgenetic studies. However, the five years since the final firstgeneration update (v1.07), however, have witnessed the introduction of new algorithms and analytical approaches, the growth in size of typical datasets, as well as wide deployment of multicore processors.
In response, we have developed PLINK 1.9, a comprehensive performance, scaling, and usability update. Our data indicate that its speedups frequently exceed two, and sometimes even three, orders of magnitude for several commonly used operations. PLINK 1.9’s core functional domains are unchanged from that of its predecessor—data management, summary statistics, population stratification, association analysis, identitybydescent estimation [1] —and it is usable as a dropin replacement in most cases, requiring no changes to existing scripts. To support easier interoperation with newer software, for example BEAGLE 4 [2], IMPUTE2 [3], GATK [4], VCFtools [5], BCFtools [6] and GCTA [7], features such as the import/export of VCF and Oxfordformat files and an efficient crossplatform genomic relationship matrix (GRM) calculator have been introduced. Most pipelines currently employing PLINK 1.07 can expect to benefit from upgrading to PLINK 1.9.
A major problem remains: PLINK’s core file format can only represent unphased, biallelic data; however we are developing a second update, PLINK 2.0, to address this.
Improvements in PLINK 1.9
Bitlevel parallelism
Modern ×86 processors are designed to operate on data in (usually 64bit) machine word or (≥ 128bit) vector chunks. The PLINK 1 binary file format supports this well: the format’s packed 2bit data elements can, with the use of bit arithmetic, easily be processed 32 or 64 at a time. However, most existing programs fail to exploit opportunities for bitlevel parallelism; instead their loops painstakingly extract and operate on a single data element at a time. Replacement of these loops with bitparallel logic is, by itself, enough to speed up numerous operations by more than one order of magnitude.
For example, when comparing two DNA segments, it is frequently useful to start by computing their Hamming distance. Formally, define two sequences {a_{1},a_{2},…,a_{ m }} and {b_{1},b_{2},…,b_{ m }} where each a_{ i } and b_{ i } has a value in {0,1,2,ϕ}, representing either the number of copies of the major allele or (ϕ) the absence of genotype data. Also define an intersection set I_{a,b}:={i:a_{ i }≠ϕ and b_{ i }≠ϕ}. The “identitybystate” measure computed by PLINK can then be expressed as
where I_{a,b} denotes the size of set I_{a,b}, while a_{ i }−b_{ i } is the absolute value of a_{ i } minus b_{ i }. The old calculation proceeded roughly as follows:
IBS0 := 0 IBS1 := 0 IBS2 := 0 For i∈{1,2,…,m}:

If a_{ i }=ϕ or b_{ i }=ϕ, skip

otherwise, if a_{ i }=b_{ i }, increment IBS2

otherwise, if (a_{ i }=2 and b_{ i }=0), or (a_{ i }=0 and b_{ i }=2), increment IBS0

otherwise, increment IBS1
Return \(\frac {0\cdot \text {IBS}0 + 1\cdot \text {IBS}1 + 2\cdot \text {IBS}2}{2\cdot (\text {IBS}0 + \text {IBS}1 + \text {IBS}2)}\)
We replaced this with roughly the following, based on bitwise operations on 960marker blocks:
Pad the ends of {a_{ i }} and {b_{ i }} with ϕs, if necessary A_{ i }:={01_{2} if a_{ i }=ϕ,00_{2} if a_{ i }=0,10_{2} if a_{ i }=1,11_{2} if a_{ i }=2} B_{ i }:={01_{2} if b_{ i }=ϕ,00_{2} if b_{ i }=0,10_{2} if b_{ i }=1,11_{2} if b_{ i }=2} C_{ i }:={00_{2} if a_{ i }=ϕ,11_{2} otherwise} D_{ i }:={00_{2} if b_{ i }=ϕ,11_{2} otherwise} diff := 0 obs := 0For i∈{1,961,1921,…,m^{′}−959}:

E:=A_{i..i+959}XORB_{i..i+959}

F:=C_{i..i+959}ANDD_{i..i+959}

diff := diff + popcount(EANDF)

obs := obs + popcount(F)
Return \(\frac {\text {obs}  \text {diff}}{\text {obs}}\).
The idea is that ({C_{ i }}AND {D_{ i }}) yields a bit vector with two ones for every marker where genotype data is present for both samples, and two 0 s elsewhere, so 2I_{a,b} is equal to the number of ones in that bit vector; while (({A_{ i }}XOR {B_{ i }})AND {C_{ i }}AND {D_{ i }}) yields a bit vector with a 1 for every nucleotide difference. Refer to Additional file 1 [8] for more computational details. Our timing data (see “Performance comparisons” below) indicate that this algorithm takes less than twice as long to handle a 960marker block as PLINK 1.07 takes to handle a single marker.
Bit population count
The “popcount” function above, defined as the number of ones in a bit vector, merits further discussion. Post2008 x86 processors support a specialized instruction that directly evaluates this quantity. However, thanks to 50 years of work on the problem, algorithms exist which evaluate bit population count nearly as quickly as the hardware instruction while sticking to universally available operations. Since PLINK is still used on some older machines, we took one such algorithm (previously discussed and refined by [9]), and developed an improved SSE2based implementation. (Note that SSE2 vector instructions are supported by even the oldest x8664 processors).
The applications of bit population count extend further than might be obvious at first glance. As another example, consider computation of the correlation coefficient r between a pair of genetic variants, where some data may be missing. Formally, let n be the number of samples in the dataset, and {x_{1},x_{2},…,x_{ n }} and {y_{1},y_{2},…,y_{ n }} contain genotype data for the two variants, where each x_{ i } and y_{ i } has a value in {0,1,2,ϕ}. In addition, define
The correlation coefficient of interest can then be expressed as
Given PLINK 1 binary data, I_{x,y}, \(\overline {v}\), \(\overline {w}\), \(\overline {v^{2}}\), and \(\overline {w^{2}}\) can easily be expressed in terms of bit population counts. The dot product \(\sum _{i=1}^{n}v_{i}w_{i}\) is trickier; to evaluate it, we preprocess the data so that the genotype bit vectors X and Y encode homozygote minor calls as 00_{2}, heterozygote and missing calls as 01_{2}, and homozygote major calls as 10_{2}, and then proceed as follows:

1.
Set Z := (XORY) AND01010101… _{2}

2.
Evaluate

popcount2(((XXORY) AND (10101010… _{2}  Z)) ORZ),
where popcount2() sums 2bit quantities instead of counting set bits. (This is actually cheaper than PLINK’s regular population count; the first step of software popcount() is reduction to a popcount2() problem).


3.
Subtract the latter quantity from n.
The key insight behind this implementation is that each v_{ i }w_{ i } term is in {−1,0,1}, and can still be represented in 2 bits in an additionfriendly manner. (This is not strictly necessary for bitwise parallel processing—the partial sum lookup algorithm discussed later handles 3bit outputs by padding the raw input data to 3 bits per genotype call—but it allows for unusually high efficiency). The exact sequence of operations that we chose to evaluate the dotproduct terms in a bitwise parallel fashion is somewhat arbitrary.
We note that when computing a matrix of correlation coefficients between all pairs of variants, if no genotype data is absent, then I_{x,y} is invariant, \(\overline {v}\) and \(\overline {v^{2}}\) do not depend on y, and \(\overline {w}\) and \(\overline {w^{2}}\) do not depend on x. Thus, these five values would not need to be recomputed for each variant pair at O(m^{2}n) total time cost; they could instead be precomputed outside the main loop at a total cost of O(mn) time and O(m) space. PLINK 1.9 optimizes this common case.
See popcount_longs() in plink_common.c for our primary bit population count function, and plink_ld.c for several correlation coefficient evaluation functions.
Multicore and cluster parallelism
Modern x86 processors also contain increasing numbers of cores, and computational workloads in genetic studies tend to contain large “embarrassingly parallel” steps which can easily exploit additional cores. Therefore, PLINK 1.9 autodetects the number of cores present in the machine it is running on, and many of its heavyduty operations default to employing roughly that number of threads. (This behavior can be manually controlled with the –threads flag.) Most of PLINK 1.9’s multithreaded computations use a simple set of crossplatform C functions and macros, which compile to pthread library idioms on Linux and OS X, and OSspecific idioms like _beginthreadex() on Windows.
PLINK 1.9 also contains improved support for distributed computation: the –parallel flag makes it easy to split large matrix computations across a cluster, while –writevarranges simplifies splitting of pervariant computations.
Graphics processing units (GPUs) remain as a major unexploited computational resource. We have made the development of GPUspecific code a low priority since their installed base is much smaller than that of multicore processors, and the speedup factor over wellwritten multithreaded code running on similarcost, less specialized hardware is usually less than 10x [10,11]. However, we do plan to build out GPU support for the heaviestduty computations after most of our other PLINK 2 development goals are achieved.
Memory efficiency
To make it possible for PLINK 1.9 to handle the huge datasets that benefit the most from these speed improvements, the program core no longer keeps the main genomic data matrix in memory; instead, most of its functions only load data for a single variant, or a small window of variants, at a time. Sample × sample matrix computations still normally require additional memory proportional to the square of the sample size, but –parallel gets around this:
calculates 1/40th of the genomic relationship matrix per run, with correspondingly reduced memory requirements.
Other noteworthy algorithms
Partial sum lookup
Each entry of a weighted genomic distance matrix between pairs of individuals is a sum of permarker terms. Given PLINK 1 binary data, for any specific marker, there are seven distinct cases at most:

1.
Both genotypes are homozygous for the major allele.

2.
One is homozygous major, and the other is heterozygous.

3.
One is homozygous major, and the other is homozygous minor.

4.
Both are heterozygous.

5.
One is heterozygous, and the other is homozygous minor.

6.
Both are homozygous minor.

7.
At least one genotype is missing.
For example, the GCTA genomic relationship matrix is defined by the following permarker increments, where q is the minor allele frequency:

1.
\(\frac {(22q)(22q)}{2q(1q)}\)

2.
\(\frac {(22q)(12q)}{2q(1q)}\)

3.
\(\frac {(22q)(02q)}{2q(1q)}\)

4.
\(\frac {(12q)(12q)}{2q(1q)}\)

5.
\(\frac {(12q)(02q)}{2q(1q)}\)

6.
\(\frac {(02q)(02q)}{2q(1q)}\)

7.
0; subtract 1 from the final denominator instead, in another loop
This suggests the following matrix calculation algorithm, as a first draft:

1.
Initialize all distance/relationship partial sums to zero.

2.
For each marker, calculate and save the seven possible increments in a lookup table, and then refer to the table when updating partial sums. This replaces several floating point adds/multiplies in the inner loop with a single addition operation.
We can substantially improve on this by handling multiple markers at a time. Since seven cases can be distinguished by three bits, we can compose a sequence of operations which maps a pair of padded 2bit genotypes to seven different 3bit values in the appropriate manner. On 64bit machines, 20 3bit values can be packed into a machine word—for example, let bits 02 describe the relation at marker #0, bits 35 describe the relation at marker #1, and so forth, all the way up to bits 5759 describing the relation at marker #19—so this representation lets us instruct the processor to act on 20 markers simultaneously.
Then, we need to perform the update
where the x_{ i }’s are bit trios, and the f_{ i }’s map them to increments. This could be done with 20 table lookups and floating point addition operations. Or, the update could be restructured as
where x_{{0−4}} denotes the lowestorder 15 bits, and f_{{0−4}} maps them directly to f_{0}(x_{0})+f_{1}(x_{1})+f_{2}(x_{2})+f_{3}(x_{3})+f_{4}(x_{4}); similarly for f_{{5−9}}, f_{{10−14}}, and f_{{15−19}}. In exchange for some precomputation—four tables with 2^{15} entries each; total size 1 MB, which is not onerous for modern L2/L3 caches—this restructuring licenses the use of four table lookups and adds per update instead of twenty. See fill_weights_r() and incr_dists_r() in plink_calc.c for source code.
HardyWeinberg equilibrium and Fisher’s exact tests
Under some population genetic assumptions such as minimal inbreeding, genotype frequencies for a biallelic variant can be expected to follow the HardyWeinberg proportions
where p is the frequency of allele A_{1} and q=1−p is the frequency of allele A_{2} [12]. It is now common for bioinformaticians to use an exact test for deviation from HardyWeinberg equilibrium (HWE) to help detect genotyping error and major violations of the HardyWeinberg assumptions.
PLINK 1.0 used the SNPHWE algorithm in a paper by Wigginton et al. [13] to execute this test. SNPHWE exploits the fact that, while the absolute likelihood of a contingency table involves large factorials which are fairly expensive to evaluate, the ratios between its likelihood and that of adjacent tables are simple since the factorials almost entirely cancel out [14]. More precisely, given n diploid samples containing a total of n_{1} copies of allele A_{1} and n_{2} copies of allele A_{2} (so n_{1}+n_{2}=2n), there are \(\frac {(2n)!}{n_{1}!n_{2}!}\) distinct ways for the alleles to be distributed among the samples, and \(\frac {(2^{n_{12}})(n!)}{((n_{1}n_{12})/2)!n_{12}!((n_{2}n_{12})/2)!}\) of those ways correspond to exactly n_{12} heterozygotes when n_{12} has the same parity as n_{1} and n_{2}. Under HardyWeinberg equilibrium, each of these ways is equally likely. Thus, the ratio between the likelihoods of observing exactly n_{12}=k+2 heterozygotes and exactly n_{12}=k heterozygotes, under HardyWeinberg equilibrium and fixed n_{1} and n_{2}, is
SNPHWE also recognizes that it is unnecessary to start the computation with an accurate absolute likelihood for one table. Since the final pvalue is computed as
it is fine for all computed likelihoods to be relative values off by a shared constant factor, since that constant factor will cancel out. This eliminates the need for loggamma approximation.
While studying the software, we made two additional observations:

Its size O(n) memory allocation (where n is the sum of all contingency table entries) could be avoided by reordering the calculation; it is only necessary to track a few partial sums.

Since likelihoods decay supergeometrically as one moves away from the most probable table, only \(O(\sqrt {n})\) of the likelihoods can meaningfully impact the partial sums; the sum of the remaining terms is too small to consistently affect even the 10th significant digit in the final pvalue. By terminating the calculation when all the partial sums stop changing (due to the newest term being too tiny to be tracked by IEEE754 doubleprecision numbers), computational complexity is reduced from O(n) to \(O(\sqrt {n})\) with no loss of precision. See Figure 1 for an example.
PLINK 1.0 also has association analysis and quality control routines which perform Fisher’s exact test on 2×2 and 2×3 tables, using the FEXACT network algorithm from Mehta et al. [15,16]. The 2×2 case has the same mathematical structure as the HardyWeinberg equilibrium exact test, so it was straightforward to modify the earlytermination SNPHWE algorithm to handle it. The 2×3 case is more complicated, but retains the property that only \(O(\sqrt {\mathrm {\# of tables}})\) relative likelihoods need to be evaluated, so we were able to develop a function to handle it in O(n) time; see Figure 2 for more details. Our timing data indicate that our new functions are consistently faster than both FEXACT and the update to the network algorithm by Requena et al. [17].
Standalone source code for earlytermination SNPHWE and Fisher’s 2×2/ 2×3 exact test is posted at [18]. Due to recent calls for use of midp adjustments in biostatistics [19,20], all of these functions have midp modes, and PLINK 1.9 exposes them.
We note that, while the HardyWeinberg equilibrium exact test is only of interest to geneticists, Fisher’s exact test has wider application. Thus, we are preparing another paper which discusses these algorithms in more detail, with proofs of numerical error bounds and a full explanation of how the Fisher’s exact test algorithm extends to larger tables.
Haplotype block estimation
It can be useful to divide the genome into blocks of variants which appear to be inherited together most of the time, since observed recombination patterns are substantially more “blocklike” than would be expected under a model of uniform recombination [21]. PLINK 1.0’s –blocks command implements a method of identifying these haplotype blocks by Gabriel et al. [22]. (More precisely, it is a restricted port of Haploview’s [23] implementation of the method).
This method is based on 90% confidence intervals (as defined by Wall and Pritchard [21]) for Lewontin’s D^{′} disequilibrium statistic for pairs of variants. Depending on the confidence interval’s boundaries, a pair of variants is classified as “strong linkage disequilibrium (LD)”, “strong evidence for historical recombination”, or “inconclusive”; then, contiguous groups of variants where “strong LD” pairs outnumber “recombination” pairs by more than 19 to 1 are greedily selected, starting with the longest basepair spans.
PLINK 1.9 accelerates this in several ways:

Estimation of diplotype frequencies and maximumlikelihood D^{′} has been streamlined. Bit population counts are used to fill the contingency table; then we use the analytic solution to Hill’s diplotype frequency cubic equation [24,25] and only compute and compare log likelihoods in this step when multiple solutions to the equation are in the valid range.

90% confidence intervals were originally estimated by computing relative likelihoods at 101 points (corresponding to D^{′}=0,D^{′}=0.01,…,D^{′}=1) and checking where the resulting cumulative distribution function (cdf) crossed 5% and 95%. However, the likelihood function rarely has more than one extreme point in (0,1) (and the full solution to the cubic equation reveals the presence of additional extrema); it is usually possible to exploit this unimodality to establish good bounds on key cdf values after evaluating just a few likelihoods. In particular, many confidence intervals can be classified as “recombination” after inspection of just two of the 101 points; see Figure 3.

Instead of saving the classification of every variant pair and looking up the resulting massive table at a later point, we just update a small number of “strong LD pairs within last k variants” and “recombination pairs within last k variants” counts while processing the data sequentially, saving only final haploblock candidates. This reduces the amount of time spent looking up outofcache memory, and also allows much larger datasets to be processed.

Since “strong LD” pairs must outnumber “recombination” pairs by 19 to 1, it does not take many “recombination” pairs in a window before one can prove no haploblock can contain that window. When this bound is crossed, we take the opportunity to entirely skip classification of many pairs of variants.
Most of these ideas are implemented in haploview_blocks_classify() and haploview_blocks() in plink_ld.c. The last two optimizations were previously implemented in Taliun’s “LDExplorer” R package [26].
Coordinatedescent LASSO
PLINK 1.9 includes a basic coordinatedescent LASSO implementation [27] (–lasso), which can be useful for phenotypic prediction and related applications. See Vattikuti et al. for discussion of its theoretical properties [28].
Newly integrated thirdparty software
PLINK 1.0 commands
Many teams have significantly improved upon PLINK 1.0’s implementations of various commands and made their work open source. In several cases, their innovations have been integrated into PLINK 1.9; examples include

Pahl et al.’s PERMORY algorithm for fast permutation testing [29],

Wan et al.’s BOOST software for fast epistasis testing [30],

Ueki, Cordell, and Howey’s –fastepistasis variance correction and jointeffects test [31,32],

Taliun, Gamper, and Pattaro’s optimizations to Gabriel et al.’s haplotype block identification algorithm (discussed above) [26], and

Pascal Pons’s winning submission to the GWAS Speedup logistic regression crowdsourcing contest [33]. (The contest was designed by PoRu Loh, run by Babbage Analytics & Innovation and TopCoder, and subsequent analysis and code preparation were performed by Andrew Hill, Ragu Bharadwaj, and Scott Jelinsky. A manuscript is in preparation by these authors and Iain Kilty, Kevin Boudreau, Karim Lakhani and Eva Guinan.)
In all such cases, PLINK’s citation instructions direct users of the affected functions to cite the original work.
Multithreaded gzip
For many purposes, compressed text files strike a good balance between ease of interpretation, loading speed, and resource consumption. However, the computational cost of generating them is fairly high; it is not uncommon for data compression to take longer than all other operations combined. To make a dent in this bottleneck, we have written a simple multithreaded compression library function based on Mark Adler’s excellent pigz program [34], and routed most of PLINK 1.9’s gzipping through it. See parallel_compress() in pigz.c for details.
Convenience features
Import and export of Variant Call Format (VCF) and Oxfordformatted data
PLINK 1.9 can import data from Variant Call Format (–vcf), binary VCF (–bcf), and Oxfordformat (–data, –bgen) files. However, since it cannot handle genotype likelihoods, phase information or variants with more than two alleles, the import process can be quite lossy. Specifically,

With Oxfordformat files, genotype likelihoods smaller than 0.9 are normally treated as missing calls, and the rest are treated as hard calls. –hardcallthreshold can be used to change the threshold, or request independent pseudorandom calls based on the likelihoods in the file.

Phase is discarded.

By default, when a VCF variant has more than one alternate allele, only the most common alternate is retained; all other alternate calls are converted to missing. –bialleliconly can be used to skip variants with multiple alternate alleles.
Export to these formats is also possible, via –recode vcf and –recode oxford.
Unplaced contig and nonhuman species support
When the –allowextrachr or –aec flag is used, PLINK 1.9 allows datasets to contain unplaced contigs or other arbitrary chromosome names, and most commands will handle them in a reasonable manner. Also, arbitrary nonhuman species (with haploid or diploid genomes) can now be specified with –chrset.
Commandline help
To improve the experience of using PLINK interactively, we have expanded the –help flag’s functionality. When invoked with no parameters, it now prints an entire minimanual. Given keyword(s), it instead searches for and prints minimanual entries associated with those keyword(s), and handles misspelled keywords and keyword prefixes in a reasonable manner.
A comment on withinfamily analysis
Most of our discussion has addressed computational issues. However, there is one methodological issue that deserves a brief comment. The online documentation of PLINK 1.07 weighed the pros and cons of its permutation procedure for withinfamily analysis of quantitative traits (QFAM) with respect to the standard quantitative transmission disequilibrium test (QTDT) [35]. It pointed out that likelihoodbased QTDT enjoyed the advantages of computational speed and increased statistical power. However, a comparison of statistical power is only meaningful if both procedures are anchored to the same Type 1 error rate with respect to the null hypothesis of no linkage with a causal variant, and Ewens et al. has shown that the QTDT is not robust against certain forms of confounding (population stratification) [36]. On the other hand, the validity of a permutation procedure such as QFAM only depends on the applicability of Mendel’s laws. When this nicety is combined with the vast speedup of permutation in PLINK 1.9, a given user may now decide to rate QFAM more highly relative to QTDT when considering available options for withinfamily analysis.
Performance comparisons
In the following tables, running times are collected from seven machines operating on three datasets.

“Mac2” denotes a MacBook Pro with a 2.8 Ghz Intel Core 2 Duo processor and 4GB RAM running OS X 10.6.8.

“Mac12” denotes a Mac Pro with two 2.93 Ghz Intel 6core Xeon processors and 64GB RAM running OS X 10.6.8.

“Linux322” denotes a machine with a 2.4 Ghz Intel Core 2 Duo E6600 processor and 1GB RAM running 32bit Ubuntu Linux.

“Linux328” denotes a machine with a 3.4 Ghz Intel Core i73770 processor (8 cores) and 8GB RAM running 32bit Ubuntu Linux.

“Linux64512” denotes a machine with sixtyfour AMD 8core Opteron 6282 SE processors and 512GB RAM running 64bit Linux.

“Win322” denotes a laptop with a 2.4 Ghz Intel Core i52430 M processor (2 cores) and 4GB RAM running 32bit Windows 7 SP1.

“Win642” denotes a machine with a 2.3 Ghz Intel Celeron G1610T processor (2 cores) and 8GB RAM running 64bit Windows 8.

“synth1” refers to a 1000 sample, 100000 variant synthetic dataset generated with HAPGEN2 [37], while “synth1p” refers to the same dataset after one round of –indeppairwise 50 5 0.5 pruning (with 76124 markers remaining). For case/control tests, PLINK 1.9’s –tailpheno 0 command was used to downcode the quantitative phenotype to case/control.

“synth2” refers to a 4000 case, 6000 control synthetic dataset with 88025 markers on chromosomes 1922 generated by resampling HapMap and 1000 Genomes data with simuRare [38] and then removing monomorphic loci. “synth2p” refers to the same dataset after one round of –indeppairwise 700 70 0.7 pruning (with 71307 markers remaining).

“1000g” refers to the entire 1092 sample, 39637448 variant 1000 Genomes project phase 1 dataset [39]. “chr1” refers to chromosome 1 from this dataset, with 3001739 variants. “chr1snp” refers to chromosome 1 after removal of all nonSNPs and one round of –indeppairwise 20000 2000 0.5 pruning (798703 markers remaining). Pedigree information was not added to these datasets before our tests.
All times are in seconds. To reduce diskcaching variance, timing runs are preceded by “warmup” commands like plink –freq. PLINK 1.07 was run with the –noweb flag. “nomem” indicates that the program ran out of memory and there was no lowmemory mode or other straightforward workaround. A tilde indicates that runtime was extrapolated from several smaller problem instances.
Initialization and basic I/O
Table 1 displays execution times for plink –freq, one of the simplest operations PLINK can perform. These timings reflect fixed initialization and I/O overhead. (Due to the use of warmup runs, they do not include disk latency).
Identitybystate matrices, complete linkage clustering
The PLINK 1.0 –cluster –matrix flag combination launches an identitybystate matrix calculation and writes the result to disk, and then performs complete linkage clustering on the data; when –ppc is added, a pairwise population concordance constraint is applied to the clustering process. As discussed earlier, PLINK 1.9 employs an XOR/bit population count algorithm which speeds up the matrix calculation by a large constant factor; the computational complexity of the clustering algorithm has also been reduced, from O(n^{3}) to O(n^{2} logn). (Further improvement of clustering complexity, to O(n^{2}), is possible in some cases [40].)
In Table 2, we compare PLINK 1.07 and PLINK 1.9 execution times under three scenarios: identitybystate (IBS) matrix calculation only (–cluster –matrix –K [sample count  1] in PLINK 1.07, –distance ibs square in PLINK 1.9), IBS matrix + standard clustering (–cluster –matrix for both versions), and identitybydescent (IBD) report generation (–Zgenome.)
(Note that newer algorithms such as BEAGLE’s fastIBD [41] generate more accurate IBD estimates than PLINK –Zgenome. However, the –Zgenome report contains other useful information.)
Genomic relationship matrices
GCTA’s –makegrmbin command (–makegrm in early versions) calculates the variancestandardized genomic relationship matrix used by many of its other commands. The latest implementation as of this writing (v1.24) is very fast, but only runs on 64bit Linux, uses single instead of doubleprecision arithmetic, and has a high memory requirement.
PLINK 1.9’s implementation of this calculation is designed to compensate for GCTA 1.24’s limitations—it is crossplatform, works in lowmemory environments, and uses doubleprecision arithmetic while remaining within a factor of 25 on speed. See Table 3 for timing data. The comparison is with GCTA 1.24 on 64bit Linux, and v1.02 elsewhere.
Linkage disequilibriumbased variant pruning
The PLINK 1.0 –indeppairwise command is frequently used in preparation for analyses which assume approximate linkage equilibrium. In Table 4, we compare PLINK 1.07 and PLINK 1.9 execution times for some reasonable parameter choices. The r^{2} threshold for “synth2” was chosen to make the “synth1p” and “synth2p” pruned datasets contain similar number of SNPs, so Tables 2 and 3 could clearly demonstrate scaling with respect to sample size.
Haplotype block estimation
Table 5 demonstrates the impact of our rewrite of –blocks. Due to a minor bug in PLINK 1.0’s handling of lowMAF variants, we pruned each dataset to contain only variants with MAF ≥0.05 before running –blocks. 95506 markers remained in the “synth1” dataset, and 554549 markers remained in “chr1”. A question mark indicates that the extrapolated runtime may not be valid since we suspect Haploview or PLINK 1.07 would have run out of memory before finishing.
Association analysis max(T) permutation tests
PLINK 1.0’s basic association analysis commands were quite flexible, but the powerful max(T) permutation test suffered from poor performance. PRESTO [42] and PERMORY introduced major algorithmic improvements (including bit population count) which largely solved the problem. Table 6 demonstrates that PLINK 1.9 successfully extends the PERMORY algorithm to the full range of PLINK 1.0’s association analyses, while making Fisher’s exact test practical to use in permutation tests. (There is no 64bit Windows PERMORY build, so the comparisons on the Win642 machine are between 64bit PLINK and 32bit PERMORY.)
PLINK 2.0 design
Despite its computational advances, we recognize that PLINK 1.9 can ultimately still be an unsatisfactory tool for working with imputed genomic data, due to the limitations of the PLINK 1 binary file format. To address this, we designed a new core file format capable of representing most of the information emitted by modern imputation tools, which is the cornerstone of our plans for PLINK 2.0.
Multiple data representations
As discussed earlier, PLINK 1 binary is inadequate in three ways: likelihoods strictly between 0 and 1 cannot be represented, phase information cannot be stored, and variants are limited to two alleles. This can be addressed by representing all calls probabilistically, and introducing a few other extensions. Unfortunately, this would make PLINK 2.0’s representation of PLINK 1format data so inefficient that it would amount to a serious downgrade from PLINK 1.9 for many purposes.
Therefore, our new format defines several data representations, one of which is equivalent to PLINK 1 binary, and allows different files, or even variants within a single file, to use different representations. To work with this, PLINK 2.0 will include a translation layer which allows individual functions to assume a specific representation is used. As with the rest of PLINK’s source code, this translation layer will be GPLv3licensed open source; and unlike most of the other source code, we are explicitly designing it to be usable as a standalone library. PLINK 2.0 will also be able to convert files/variants from one data representation to another, making it practical for thirdparty tools lacking access to the library to demand a specific representation.
Reference vs. alternate alleles
The nowubiquitous VCF file format requires reference alleles to be distinguished from alternate alleles, and an increasing number of software tools and pipelines do not tolerate scrambling of the two. This presents an interoperability problem for PLINK: while it was theoretically possible to handle binary data with PLINK 1.0 in a manner that preserved the reference vs. alternate allele distinction when it was originally present, with constant use of –keepalleleorder and related flags, doing so was inconvenient and errorprone, especially since the accompanying native.ped/.map and.tped/.tfam text formats had no place to store that information. PLINK 1.9’s –a2allele flag, which can import that information from a VCF file, provides limited relief, but it is still necessary for users to fight against the program’s major/minorallele based design.
We aim to solve this problem for good in PLINK 2.0. The file format explicitly defines reference vs. alternate alleles, and this information will be preserved across runs by default. In addition, the file format will include a flag distinguishing provisional reference allele assignments from those derived from an actual reference genome. When PLINK 2.0 operates on.ped/.map or similar data lacking a reference vs. alternate distinction, it will treat a highestfrequency allele as the reference, while flagging it as a provisional assignment. When a file with flaggedasprovisional reference alleles is merged with another file with unflagged reference alleles, the unflagged reference allele assignments take precedence. (Merges involving conflicting unflagged reference alleles will fail unless the user specifies which source file takes precedence.) It will also be straightforward to import real reference allele assignments with an analogue of –a2allele.
Data compression
PLINK 1.9 demonstrates the power of a weak form of compressive genomics [43]: by using bit arithmetic to perform computation directly on compressed genomic data, it frequently exhibits far better performance than programs which require an explicit decompression step. But its “compressed format” is merely a tight packing which does not support the holy grail of true sublinear analysis.
To do our part to make “strong” sublinear compressive genomics a reality, the PLINK 2 file format will introduce support for “deviations from most common value” storage of lowMAF variants. For datasets containing many samples, this captures much of the storage efficiency benefit of having real reference genomes available, without the drawback of forcing all programs operating on the data to have access to a library of references. Thanks to PLINK 2.0’s translation layer and file conversion facilities, programmers will be able to ignore this feature during initial development of a tool, and then work to exploit it after basic functionality is in place.
We note that LDbased compression of variant groups is also possible, and Sambo’s SNPack software [44] applies this to the PLINK 1 binary format. We do not plan to support this in PLINK 2.0 due to the additional software complexity required to handle probabilistic and multiallelic data, but we believe this is a promising avenue for development and look forward to integrating it in the future.
Remaining limitations
PLINK 2.0 is designed to meet the needs of tomorrow’s genomewide association studies and populationgenetics research; in both contexts, it is appropriate to apply a single genomic coordinate system across all samples, and preferred sample sizes are large enough to make computational efficiency a serious issue.
Wholeexome and wholegenome sequencing also enables detailed study of structural variations which defy clean representation under a single coordinate system; and the number of individuals in such studies is typically much smaller than the tens or even hundreds of thousands which are sometimes required for effective GWAS. There are no plans to make PLINK suitable for this type of analysis; we strongly recommend the use of another software package, such as PLINK/SEQ [45], which is explicitly designed for it. This is why the PLINK 2 file format will still be substantially less expressive than VCF.
An important consequence is that, despite its ability to import and export VCF files, PLINK should not be used for management of genomic data which will be subject to both types of analysis, because it discards all information which is not relevant for its preferred type. However, we will continue to extend PLINK’s ability to interpret VCFlike formats and interoperate with other popular software.
Availability and requirements

Project name: Secondgeneration PLINK

Project (source code) home page:https://www.coggenomics.org/plink2/(https://github.com/chrchang/plinkng)

Operating systems: Linux (32/64bit), OS X (64bit Intel), Windows (32/64bit)

Programming language: C, C++

Other requirements (when recompiling): GCC version 4, a few functions also require LAPACK 3.2

License: GNU General Public License version 3.0 (GPLv3)

Any restrictions to use by nonacademics: none
Availability of supporting data
The test data and the source code snapshots supporting the results of this article are available in the GigaScience repository, GigaDB [8].
Abbreviations
 PLINK:

The software toolset that is the main subject of this paper. The name was originally shorthand for “population linkage”
 BEAGLE:

A software package capable of highaccuracy haplotype phasing, genotype imputation, and identitybydescent estimation, developed by Browning [2]
 GCTA:

Genomewide Complex Trait Analysis. This refers to both the statistical method and the software implementation discussed in [7]
 VCF:

Variant Call Format [5]
 x86:

A family of backward compatible instruction set architectures based on the Intel 8086 CPU
 IBS:

Identitybystate. A simple measure of genomic similarity, equal to the number of identical alleles divided by the number of observations
 popcount:

Bit population count. The number of ’1’ bits in a bit vector
 XOR:

Exclusiveor. A binary logical operation that evaluates to true if exactly one of its arguments is true
 SSE2:

Streaming SIMD Extensions 2. A SIMD (single instruction, multiple data) processor supplementary instruction set first introduced by Intel with the initial version of the Pentium 4 in 2001
 GPU:

Graphics processing unit
 HWE:

HardyWeinberg equilibrium
 SNP:

Singlenucleotide polymorphism
 FEXACT:

A network algorithm for evaluating Fisher’s exact test pvalues, developed by Mehta et al. [15,16]
 LD:

Linkage disequilibrium
 PERMORY:

A software package designed to perform efficient permutation tests for largescale genetic data sets, developed by Pahl et al. [29]
 GWAS:

GenomeWide Association Study
 QFAM:

A familybased quantitative trait association analysis procedure, introduced by PLINK 1.0, which combines a simple linear regression of phenotype on genotype with a special permutation test which corrects for family structure
 QTDT:

Quantitative Transmission Disequilibrium Tests, developed primarily by Abecasis et al. [35]
 Ghz:

Gigahertz
 GB:

Gigabyte
 RAM:

Randomaccess memory
 I/O:

Input/output
 MAF:

Minor allele frequency. Frequency of the least common allele that is still present in a population
 GPLv3:

GNU General Public License version 3
References
 1
Purcell S, Neale B, ToddBrown K, Thomas L, Ferreira M, Bender D, et al. Plink: A tool set for wholegenome association and populationbased linkage analyses. Am J Hum Genet. 2007; 81:559–75.
 2
Browning B, Browning S. Improving the accuracy and efficiency of identity by descent detection in population data. Genetics. 2013; 194:459–71.
 3
Howie B, Donnelly P, Marchini J. A flexible and accurate genotype imputation method for the next generation of genomewide association studies. PLoS Genet. 2009; 5:1000529.
 4
McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, et al. The genome analysis toolkit: A mapreduce framework for analyzing nextgeneration dna sequencing data. Genome Res. 2010; 20:1297–303.
 5
Danecek P, Auton A, Abecasis G, Albers C, Banks E, DePristo M, et al. The variant call format and vcftools. Bioinformatics. 2011; 27:2156–8.
 6
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, 1000 Genome Project Data Processing Subgroup, et al. The sequence alignment/map format and samtools. Bioinformatics. 2009; 25:2078–9.
 7
Yang J, Lee S, Goddard M, Visscher P. Gcta: A tool for genomewide complex trait analysis. Am J Hum Genet. 2011; 88:76–82.
 8
Chang C, Chow C, Tellier L, Vattikuti S, Purcell S, Lee J. Software and Supporting Material for "Secondgeneration PLINK: Rising to the Challenge of Larger and Richer Datasets". GigaScience Database. http://dx.doi.org/10.5524/100116.
 9
Dalke A. Update: Faster Population Counts. http://www.dalkescientific.com/writings/diary/archive/2011/11/02/faster_popcount_update.html.
 10
Lee V, Kim C, Chhugani J, Deisher M, Kim D, Nguyen A, et al. Debunking the 100x gpu vs. cpu myth: an evaluation of throughput computing on cpu and gpu. In: Proceedings of the 37th Annual International Symposium on Computer Architecture: 1923 June 2010. SaintMalo, France,: ACM: 2010. p. 451–460.
 11
Haque I, Pande V, Walters W. Anatomy of highperformance 2d similarity calculations. J Chem Inf Model. 2011; 51:2345–51.
 12
Hardy H. Mendelian proportions in a mixed population. Science. 1908; 28:49–50.
 13
Wigginton J, Cutler D, Abecasis G. A note on exact tests of hardyweinberg equilibrium. Am J Hum Genet. 2005; 76:887–93.
 14
Guo S, Thompson E. Performing the exact test of hardyweinberg proportion for multiple alleles. Biometrics. 1992; 48:361–72.
 15
Mehta C, Patel N. Algorithm 643: Fexact: a fortran subroutine for fisher’s exact test on unordered r ×c contingency tables. ACM Trans Math Softw. 1986; 12:154–61.
 16
Clarkson D, Fan Y, Joe H. A remark on algorithm 643: Fexact: an algorithm for performing fisher’s exact test in r x c contingency tables. ACM Trans Math Softw. 1993; 19:484–8.
 17
Requena F, Martín Ciudad N. A major improvement to the network algorithm for fisher’s exact test in 2 ×c contingency tables. J Comp Stat & Data Anal. 2006; 51:490–8.
 18
Chang C. Standalone C/C++ Exact Statistical Test Functions. https://github.com/chrchang/stats.
 19
Lydersen S, Fagerland M, Laake P. Recommended tests for association in 2 ×2 tables. Statist Med. 2009; 28:1159–75.
 20
Graffelman J, Moreno V. The mid pvalue in exact tests for hardyweinberg equilibrium. Stat Appl Genet Mol Bio. 2013; 12:433–48.
 21
Wall J, Pritchard J. Assessing the performance of the haplotype block model of linkage disequilibrium. Am J Hum Genet. 2003; 73:502–15.
 22
Gabriel S, Schaffner S, Nguyen H, Moore J, Roy J, Blumenstiel B, et al. The structure of haplotype blocks in the human genome. Science. 2002; 296:2225–9.
 23
Barrett J, Fry B, Maller J, Daly M. Haploview: analysis and visualization of ld and haplotype maps. Bioinformatics. 2005; 21:263–5.
 24
Hill W. Estimation of linkage disequilibrium in randomly mating populations. Heredity. 1974; 33:229–39.
 25
Gaunt T, Rodríguez S, Day I. Cubic exact solutions for the estimation of pairwise haplotype frequencies: implications for linkage disequilibrium analyses and a web tool ’cubex’. BMC Bioinformatics. 2007; 8:428.
 26
Taliun D, Gamper J, Pattaro C. Efficient haplotype block recognition of very long and dense genetic sequences. BMC Bioinformatics. 2014; 15:10.
 27
Friedman J, Hastie T, Höfling H, Tibshirani R. Pathwise coordinate optimization. Ann Appl Stat. 2007; 1:302–32.
 28
Vattikuti S, Lee J, Chang C, Hsu S, Chow C. Applying compressed sensing to genomewide association studies. GigaScience. 2014; 3:10.
 29
Steiß V, Letschert T, Schäfer H, Pahl R. Permorympi: A program for highspeed parallel permutation testing in genomewide association studies. Bioinformatics. 2012; 28:1168–9.
 30
Wan X, Yang C, Yang Q, Xue H, Fan X, Tang N, et al. Boost: A fast approach to detecting genegene interactions in genomewide casecontrol studies. Am J Hum Genet. 2010; 87:325–40.
 31
Ueki M, Cordell H. Improved statistics for genomewide interaction analysis. PLoS Genet. 2012; 8:1002625.
 32
Howey R. CASSI: GenomeWide Interaction Analysis Software. http://www.staff.ncl.ac.uk/richard.howey/cassi.
 33
GWASSpeedup Problem Statement. http://community.topcoder.com/longcontest/?module=ViewProblemStatement&rd=15637&pm=12525.
 34
Adler M. Pigz: Parallel Gzip. http://zlib.net/pigz/.
 35
Abecasis G, Cardon L, Cookson W. A general test of association for quantitative traits in nuclear families. Am J Hum Genet. 2000; 66:279–92.
 36
Ewens W, Li M, Spielman R. A review of familybased tests for linkage disequilibrium between a quantitative trait and a genetic marker. PLoS Genet. 2008; 4:1000180.
 37
Su Z, Marchini J, Donnelly P. Hapgen2: Simulation of multiple disease snps. Bioinformatics. 2011; 27:2304–5.
 38
Xu Y, Wu Y, Song C, Zhang H. Simulating realistic genomic data with rare variants. Genet Epidemiol. 2013; 37:163–72.
 39
The 1000 Genomes Project Consortium. An integrated map of genetic variation from 1,092 human genomes. Nature. 2012; 491:56–65.
 40
Defays D. An efficient algorithm for a complete link method. Comput J. 1977; 20:364–6.
 41
Browning B, Browning S. A fast, powerful method for detecting identity by descent. Am J Hum Genet. 2011; 88:173–82.
 42
Browning B. Presto: rapid calculation of order statistic distributions and multipletesting adjusted pvalues via permutation for one and twostage genetic association studies. BMC Bioinformatics. 2008; 9:309.
 43
Loh P, Baym M, Berger B. Compressive genomics. Nat Biotechnol. 2012; 30:627–30.
 44
Sambo F, Di Camillo B, Toffolo G, Cobelli C. Compression and fast retrieval of snp data. Bioinformatics. 2014; 30:495.
 45
PLINK/SEQ: A Library for the Analysis of Genetic Variation Data. https://atgu.mgh.harvard.edu/plinkseq/.
Acknowledgements
We thank Stephen D.H. Hsu for helpful discussions. We also continue to be thankful to PLINK 1.9 users who perform additional testing of the program, report bugs, and make useful suggestions.
Christopher Chang and Laurent Tellier were supported by BGI Hong Kong and Shenzhen Municipal Government of China grant CXB201108250094A. Carson Chow and Shashaank Vattikuti were supported by the Intramural Research Program of the NIH, NIDDK.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
SMP and Ch C designed the software. Ch C drafted the manuscript and did most of the v1.9 C/C++ programming. Ca C, SV, and JJL drove early v1.9 feature development and wrote MATLAB prototype code. Ca C, LCAMT, SV, SMP, and JJL assisted with v1.9 software testing. All authors read and approved the final manuscript.
Additional file
Additional file 1
Detailed description of software bit population count, as applied to identitybystate computation.
Rights and permissions
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Chang, C.C., Chow, C.C., Tellier, L.C. et al. Secondgeneration PLINK: rising to the challenge of larger and richer datasets. GigaSci 4, 7 (2015) doi:10.1186/s1374201500478
Received:
Accepted:
Published:
Keywords
 GWAS
 Population genetics
 Wholegenome sequencing
 Highdensity SNP genotyping
 Computational statistics