psmc-rs: fast PSMC inference that reproduces the reference implementation

Yufan Han

The pairwise sequentially Markovian coalescent (PSMC) infers the history of effective population size from a single diploid genome and is one of the most widely used methods of population genomics. Its reference implementation, the C program psmc, runs on one thread, which makes bootstraps and large cohorts slow. We present psmc-rs, a reimplementation of PSMC in Rust that fits the same model as C psmc. On 40 simulated and 55 real genomes, the two programs' estimates of population size differed by at most 0.15% over the time range each genome informs. On the simulated genomes, their histories were at least 1,500 times closer to each other than either was to the true history. On a simulated genome of 30 chromosomes of 100 Mb, psmc-rs was 14.1 times as fast as C psmc on one thread and 386 times as fast on 30 threads, and it used at most 148 MB of memory, against C psmc's 1,062 MB. Its results were byte-identical on any number of threads. Fitting all 2,504 unrelated high-coverage genomes of the 1000 Genomes Project took psmc-rs 257 core-hours, where C psmc would have needed about 3,700. psmc-rs also reads VCF and multihetsep files, accounts for partly callable bins, and runs PSMC's block bootstrap in parallel. It is open source under the MIT license, and the pipeline behind every result, with the inputs and fits of the 2,504 genomes, is openly available.