Inference¶

The infer task is the main one: it runs the MCMC that estimates the effective population size shared by all loci together with the selection coefficient at each locus separately.

$ ./sweeplink infer --counts alleleCounts.txt --meta meta.txt --mu_a_A 1e-8 --out myrun

Only three things are required (--counts, --meta, --mu_a_A); everything else on this page has a default. For a worked walk-through rather than a list, see the tutorial.

Required input¶

--meta <file>: Which samples belong to which population and time point. Three tab-separated columns, no header: sample label, time in generations, population name. The labels are the column names of the counts file. Time runs forward, so 0 is the oldest time point. See file formats.
--counts <file>: Allele counts per locus. Columns are chromosome, position, then one column per sample label from the meta file. Each cell is c/n — derived allele count out of n sampled chromosomes. Missing data is 0/0.
--mu_a_A <rate>: Mutation rate from the ancestral allele a to the derived allele A, per generation. Mandatory — the run stops without it.
--mu_A_a <rate>: The reverse rate, A → a. Default 0, i.e. the derived allele can arise but never revert. Set it if you do not know which allele is ancestral — see the FAQ.

Warning

--mu is not a mutation rate. It is one of the linkage model’s transition parameters (kappa, nu, mu). The mutation rates are --mu_a_A and --mu_A_a.

Output¶

--out <prefix>: Prefix for every output file. Default sweeplink. See file formats for what each file contains.
--logFile <file>, --silent, --suppressWarnings: Redirect, suppress or quieten the console log.

Which per-locus files are written¶

The per-locus parameters — s, h and z — each have a trace, a state posteriors and a posterior mode file. They are only written when they would actually say something:

Parameter	Written by default?
`s`	Yes, unless `--s.update false` holds it fixed.
`h`	Only with `--h.update true`, since otherwise it is fixed.
`z`	No. `z` is the linkage model’s internal state and is computed from `s` and `h`, so it never holds anything those files do not.

This is not only about tidiness. A trace has one column per locus and one row per retained iteration, so on a whole chromosome a constant trace runs to gigabytes.

--<param>.writeFiles <bool>: Force the decision either way, e.g. --z.writeFiles true to inspect the linkage model’s own state, or --s.writeFiles false to skip the largest files when you only want the population size.

Population size¶

SweepLink samples log10N, not N — so read the population size off log10N, whose posterior is roughly symmetric, rather than off the strongly right-skewed N.

--N <value>: Starting value for the population size in every population. Default 10000. This sets where the chain begins, not the prior.
--log10N.alpha <shape>, --log10N.lambda <rate>: The prior, which is a Gamma on log10N in shape–rate form. Defaults 16.0 and 4.0, giving a mean log10N of 4 (i.e. N = 10,000) and a standard deviation of 1 — one order of magnitude, and clearly human-scaled. If you work on a species whose N is nowhere near 10,000, change this — see the FAQ for the arithmetic.
--log10N.update false: Hold the population size fixed at its starting value instead of sampling it. Used together with --N in the second round of the two-round workflow.

Note

N is the diploid effective population size: the drift term of the diffusion is x(1-x)/2N, so there are 2N chromosomes.

Selection¶

The grid on the selection coefficient is symmetric around zero: the negative half mirrors the positive half, which is what lets a deleterious variant be as reachable as a beneficial one.

--numStatesAbsS <M>: Number of positive values on the grid. Default 8. The full grid has 2 × (n+1) states — the mirrored negatives, two zeros, and the positives.
--abs_s_max <value>: Largest |s| when the grid is built regularly. Default 0.08. The grid is then evenly spaced: with the defaults, 0.01, 0.02, … 0.08.
--grid_abs_s <v1,v2,...>: An explicit, strictly increasing grid of positive values instead. Its length must equal --numStatesAbsS, and it takes precedence over --abs_s_max.
--s.update false: Stop sampling selection; every locus is then neutral unless you also pass --s <file>, a file of grid indices (one per line) fixing each locus.

Important

The s output files store grid indices, not values. With the default grid of 18 states, index 0 is the most negative and index 17 the most positive — and because the grid is symmetric, both index 8 and index 9 mean neutral. Any sum over the s states has to add the two together.

Dominance¶

--h.update true: Infer the dominance coefficient. Off by default — h is held fixed unless you ask for it, and its output files are then not written at all (see output files).
--numStatesH <n>: Number of dominance states. Default 3.
--grid_h <v1,v2,...>: An explicit grid on [0, 1] instead of evenly spaced values. Its length must equal --numStatesH.

Fitnesses are w_AA = 1+s, w_Aa = 1+hs, w_aa = 1, so s is the selection coefficient of the derived allele. A single h is shared across populations; s is per population.

Linkage¶

--linkage <bool>

Default true. A hidden Markov model runs along the chromosome with transition probabilities that depend on the physical distance between consecutive loci; this is what pools the evidence for a sweep across a region.

Pass --linkage false to infer every locus independently. Do this when the loci are not neighbours on the genome — most importantly when re-analysing a shortlist extracted from an earlier run, where the distances between the survivors no longer mean anything. With linkage off, --maxDistanceLoci and the transition parameters kappa, nu and mu have no effect and are held fixed.

--maxDistanceLoci <bp>

Largest inter-locus distance used when binning distances for the HMM. Default 1000000, which yields about 20 distance groups. No effect with --linkage false.

Note

With linkage off the only prior left on the selection coefficient is flat over the grid. Since that grid is symmetric and therefore contains two neutral states, neutrality carries twice the prior weight of any single non-neutral value. The effect is small and conservative, but it is there.

The MCMC chain¶

--iterations <n>: Iterations of the main chain, after burn-in. Default 100000.
--burnin <n>: Iterations per burn-in round. Default 1000.
--numBurnin <n>: Number of burn-in rounds. Default 10. Proposals are retuned between rounds.
--thinning <n>: Write every n-th sample to the trace files. Default 10.

Tip

Thinning affects only the trace files. The state posteriors are accumulated over every iteration, so with the defaults they are counts out of 100,000 even though the traces hold 10,000 samples.

--writeBurnin: Also write traces during burn-in (off by default).
--thinningStateFile <n>: How often the state file is written. Default 1000 iterations. See resuming an interrupted run below.
--noWriteState: Do not write the state file at all. It is written by default, so an interrupted run can be resumed.
--writeDensities: Write the trace of the densities.
--config <file>, --initVals <file>: Start from a parameter-configuration table (same columns as <out>_config.txt) or from a file of initial values and jump sizes.

The diffusion grid and PDE solver¶

--deltaT <value>: Time step of the PDE grid, in generations. Default 0.1.
--nGridPoints <n>: Number of allele-frequency grid points, including the boundaries 0 and 1. Default 50.

The spacing of those points defaults to a quadratic grid, which puts more of them near 0 and 1 where allele frequency trajectories spend most of their time. Three alternatives are available; pass at most one:

--uniformGrid: Evenly spaced points. Simplest, but resolves the boundaries poorly.
--logisticGrid (with --scaleLogisticGrid <a>, default 10): Points placed by a logistic transform; a controls how strongly they concentrate at the boundaries.
--chebychevGrid: Chebychev spacing.
--PDEmethod <CC|CN>: Numerical scheme for the Wright-Fisher diffusion. CC (Chang-Cooper) is the default for the forward pass; CN (Crank-Nicolson) is the default when --backward is given. In our experiments CN is the less unstable of the two.
--backward: Evaluate the likelihood with a backward pass instead of a forward one. The forward priors below are then not used at all.
--guessNumLoci <n>: How many loci to pre-allocate for when reading the counts file. Default 100000. Performance only — it does not change results, and is deliberately not recorded in .parameters.

Priors on the initial allele frequency¶

The forward pass needs a prior on the allele frequency at the first time point, and it is set separately for two kinds of locus — a variant already segregating when sampling started is not in the same situation as one that appeared later.

--forwardPrior <name>: Applies to standing variants, i.e. loci already segregating at the first time point. Default uniform.
--newMutPrior <name>: Applies to new mutations, i.e. loci that start segregating later. Default neutral.

Both accept the same four values:

Value	Meaning
`uniform`	Flat over allele frequency.
`neutral`	The standard neutral density, proportional to 1/x.
`spike`	All mass at the lowest frequency — a variant entering the population.
`beta`	A Beta distribution whose mean and effective sample size are estimated along with everything else.

A locus counts as a standing variant if its first non-missing time point is polymorphic (0 < c < n); otherwise it is treated as a new mutation. Note that this classification depends on which allele you called ancestral.

--xMin <value>: Lower cutoff of the neutral and spike priors. Default 1e-4. Only read when one of the two priors is in use. The 1/x density has no finite integral down to zero, so it has to be truncated somewhere; fixing the cutoff rather than tying it to the grid means that refining --nGridPoints does not silently change what the prior means.

When neither prior is beta there is nothing to estimate and betaMean / betaESS are held fixed. When one of them is beta:

--betaMean <value>, --betaESS <value>: Set the Beta prior explicitly instead of estimating it from the data. Give both or neither — the run stops if only one is present.
--betaMean.update <bool>, --betaESS.update <bool>: Keep either one fixed during the MCMC.
--betaESS.lambda <rate>: Rate of the exponential prior on betaESS. Default 0.1.
--updatePriorEveryNthIter <n>: Update the Beta parameters only every n-th iteration. Default 5.
--initialPriorEstimation <name>: How the initial Beta prior is estimated from the data. One of allAtT0, not0AtT0 (default), not0AtAll, allAtAll — all* uses every time point, not0* skips non-segregating sites.

Emission model¶

--emissionModel <cellAverage|pointwise>: How the probability of the observed counts is evaluated at a grid point. cellAverage (default) averages the binomial density over the allele-frequency interval the grid point represents; pointwise evaluates it at the grid point itself.

Trimming¶

--trimDataPoints <bool> or <bool,bool>: Drop leading and trailing non-segregating time points from each locus’s window. Default true,true. One value sets both; two set the early and late ends separately.

Migration and multiple populations¶

Populations are defined entirely by the third column of the meta file. s is estimated per population; h is shared.

--m <v1,v2,...>: The migration matrix, row-major, numPop² values. Each row is normalised to sum to 1, so the entries are forward migration proportions: m(i,j) is the fraction of population i that came from population j, and m(i,i) is the resident fraction. Default is the identity matrix, i.e. no migration.
--m.update false: Hold migration fixed instead of sampling it.

Note

There is no symmetric-migration option. If you want symmetry you have to impose it yourself, by passing a symmetric matrix with --m and freezing it with --m.update false. With a single population migration is inert.

Cost grows quickly with the number of populations: the hidden state is the joint frequency vector, so the state space is (nGridPoints-2)^numPop.

Performance and reproducibility¶

--numThreads <n|all|allButOne>: Default 1. The likelihood parallelises across loci, so more threads make a run substantially faster.

Warning

The parallel sums do not accumulate in a fixed order, so a multi-threaded run is not reproducible bit-for-bit, even with --fixedSeed. A single-threaded run is. Ask for threads when you want speed, and drop back to one when you need an exactly repeatable run.

--fixedSeed <n>, --addToSeed <n>: Fix the random seed, or offset it. The seed used is always recorded in <out>.parameters.

Setting any parameter directly¶

Every model parameter — m, log10N, s, h, betaMean, betaESS, kappa, nu, mu — accepts the same family of arguments:

Argument	Effect
`--<param> <value\|list\|file>`	Initial value(s). A filename is read as one column or one row.
`--<param>.update <bool>`	Whether it is sampled at all.
`--<param>.jumpSize <value>`	Initial proposal jump size.
`--<param>.propKernel <name>`	Proposal kernel: `normal`, `uniform`, `scaleLogNormal`, `integer`, `randomInteger`, `boolean`.
`--<param>.priorParameters <...>`	Hyperparameters of its prior (this is what `--log10N.alpha` and `--log10N.lambda` feed).

--N.update and --z.update have no effect: N is derived from log10N, and z is the internal state of the linkage model.

Resuming an interrupted run¶

Runs are long, and jobs get killed. Every run therefore writes <out>_state.txt, a snapshot of the current value and jump size of every parameter, rewritten in place every --thinningStateFile iterations. It is small — a few kilobytes per thousand loci, so a few MB genome-wide, against traces that run to gigabytes.

To carry on from where a dead run stopped, feed that file back as the initial values:

$ ./sweeplink infer --counts counts.txt --meta meta.txt --mu_a_A 1e-8 \
      --initVals dead_state.txt --numBurnin 0 --out rerun

--numBurnin 0: Burn-in can be skipped, because the state file restores the tuned jump sizes as well as the values. The chain therefore starts already converged and already tuned, which is the whole point of burn-in.

Tip

If <out>_state.txt is 0 bytes, the run died during burn-in — the file is created up front but only written once the main chain starts. There is nothing to resume from; start again from scratch. (The exception is --writeBurnin, which also writes state during burn-in; the header of the state file gives the iteration it came from.)

What is and is not carried over¶

	On resume
Parameter values, jump sizes	Continue where they stopped
Iteration counter	Restarts at 0
State posteriors	Accumulate from scratch

This is a warm restart, not a checkpoint: it saves the time spent converging, not the time spent sampling.

Warning

A killed run’s _statePosteriors.txt, _posteriorMode.txt and _meanVar.txt are 0 bytes — they are only written when the run finishes. Only the traces survive, because those are written as the chain goes. So there is no partial posterior to salvage, and the resumed run needs the full number of iterations you want the posterior built from.

Give the resumed run a different --out, or it will overwrite the dead run’s surviving traces.

Reproducing a run¶

Every run writes <out>.parameters listing every argument it used, including the ones left at their defaults and the random seed. Feeding that file back repeats the run:

$ ./sweeplink myrun.parameters

With --numThreads 1 the repeat is bitwise identical. Two arguments are deliberately not recorded because they do not affect results: --maxDistanceLoci and --guessNumLoci.