PReFerSim

What it simulates

PReFerSim simulates allele frequencies under the Poisson Random Field model. New mutations arise each generation, their number drawn from a Poisson distribution. Each mutation then moves through the population by Wright-Fisher drift with selection, independently of every other mutation. Sites are unlinked, which is why one run of 10 million sites equals ten runs of 1 million sites (the manual says so). That independence makes it fast.

You give a mutation rate, a distribution of fitness effects (DFE) and a demographic history made of epochs. Output files report the number of segregating sites per generation, the site frequency spectrum (SFS), genetic load, fixed and lost mutations and the frequency of every segregating mutation at the end.

Key assumptions: independent sites, a dominance coefficient h of 0.5 (additive) unless you change it, and a derived allele with selection coefficient s. In the manual's fitness scheme the genotype A1A2 has fitness 1 - h(2s) and A2A2 has 1 - 2s, where A2 is the derived allele. A positive s therefore means a harmful derived allele.

When to use it

Use it to ask what SFS or load a demographic history and a DFE produce, or to follow the trajectory of alleles that end at a given frequency. Typical questions: how much does a bottleneck change genetic load, and what does a gamma DFE look like in an exome-sized sample.

Alternatives:

• simuPOP fits better when sites are linked and you want recombination.
• stdpopsim fits better when you want a published demographic model as the starting point.

Install it yourself

PReFerSim 1.0 is on bioconda, with builds listed for linux-64 and osx-64.

conda create -n prefersim -c conda-forge -c bioconda prefersim=1.0

A micromamba dry run did not finish (downloads from conda.anaconda.org timed out), so the version rests on the anaconda.org index. The manual also describes compiling PReFerSim.c with gcc -g -o PReFerSim PReFerSim.c -lm -lgsl -lgslcblas -O3, which needs GSL.

The program takes a parameter file and a replicate number. The random seed and generator come from GSL environment variables:

GSL_RNG_SEED=42 GSL_RNG_TYPE=mrg PReFerSim params.txt 1

A minimal params.txt (keys are Name: value, any capitalization):

MutationRate: 80
DFEType: point
DFEPointSelectionCoefficient: 0
DemographicHistory: demography.txt
PrintSegSiteInfo: 1
PrintSFS: 1

The demography file has one epoch per line, oldest first: <population size in chromosomes> <generations>.

Reading the output

The recipe's helper run_prefersim.py writes /tmp/demography.txt and /tmp/params.txt from the command-line options, switches on seven output options, sets GSL_RNG_SEED and GSL_RNG_TYPE=mrg, and runs PReFerSim /tmp/params.txt 1. Files get the prefix prefersim and the replicate number 1.

The example has two epochs: 10000 chromosomes for 1 generation, then 5000 for 100 generations. The log says 32 mutations arrived before epoch 1 began and 2109 arose in total. The manual gives a Poisson mean of 40 new mutations per generation for a mutation rate of 80. The run averaged about 21 per generation, probably because the input scales with the current size: 5000 chromosomes is half the ancestral 10000.

• fixed_sites_out.txt has 101 rows (one per generation). Columns are mutations lost, mutations fixed and the sum of s of fixed mutations. Summed over rows, 1925 mutations were lost and 0 fixed.
• 2109 - 1925 leaves 184 segregating mutations. full_out.txt has exactly 184 rows: replicate, frequency, s, generation of origin and ID. The first row is 1 2e-04 0 101 2108, a new mutation at frequency 1/5000, which shows the population size is counted in chromosomes.
• num_snp_out.txt, sum_daf_out.txt, gen_load_out.txt and s_out.txt each hold 101 tab-separated values. The last value of num_snp_out.txt is 4, and the last of sum_daf_out.txt is 0.2.
• sfs_out.txt has 19 counts. The sample size n defaults to 20 in the manual. The first count is 4 and the rest are 0, so the sample holds four singletons. Four sites at 1/20 sum to 0.2, matching sum_daf_out.txt.

With no burn-in, nearly all variation is under 100 generations old. This is why the SFS is so sparse.

Your own runs

Every run changes one or more options. Each used seed 42 and one run only.

• --s 0.01 (harmful derived allele): 2095 mutations, 145 segregating, 6 singletons in the sample.
• --s=-0.01 (beneficial): 2060 mutations, 188 segregating, 12 singletons.
• --dfe gamma, default shape 0.184, scale 319.8626 and scaling 1000: 2097 mutations, 161 segregating, 8 singletons.
• --ne-anc 1000 --gen-anc 8000 --ne-now 100 --gen-now 100 (the manual's Bottleneck.txt epochs): 319,342 mutations, 291 fixed, 105 segregating. The SFS has all 19 classes filled, 89 polymorphic sites in the sample, and the final sum of derived allele frequency is 41.35.

The first three are single seeds, so their differences cannot be separated from sampling noise. Run many seeds before reading selection effects from a 100-generation run.

Pitfalls

• Sizes are in chromosomes (2N), as the manual states. 5000 means 2500 diploids.
• The default ancestral epoch is 1 generation, so there is no burn-in. The manual suggests about 16N generations in the first epoch to reach mutation-selection balance.
• Negative values break the helper's option parser when written with a space. The recipe writes --s=-0.01; --s -0.01 exits with a usage error.
• Sign of s: positive means the derived allele is harmful.
• The manual advises keeping the mutation rate (theta) at or below 2000 and says a run cannot exceed 2,000,000,000 sites.
• Use a different GSL_RNG_SEED for every replicate. Reusing it repeats the run.

Further reading

• Source, examples and manual: https://github.com/LohmuellerLab/PReFerSim
• PReFerSim manual (PDF): https://raw.githubusercontent.com/LohmuellerLab/PReFerSim/master/PReFerSim_Manual_v2.pdf
• Ortega-Del Vecchyo D, Marsden CD, Lohmueller KE (2016). Bioinformatics 32:3516-3518. https://doi.org/10.1093/bioinformatics/btw478