Step 2: Inference

This is round 1 of the genome-wide scan. sweeplink infer reads the two files from Step 1: Preparing the Input, solves the Wright-Fisher diffusion at every locus, and runs an MCMC that estimates the effective population size N shared across all loci together with the selection coefficient s at each locus separately.

Neighbouring loci are linked in the model, so the evidence for a sweep is pooled across a region instead of resting on a single variant.

New to MCMC? Click for a quick primer

SweepLink uses Markov chain Monte Carlo, so it does not return one “best” estimate. It draws thousands of samples from the posterior distribution of each parameter, which gives you a value and the uncertainty around it.

A run has two phases:

  • Burn-in — early rounds where the sampler tunes its proposals and moves toward the region of high probability. These samples are thrown away.

  • Main chain — once tuned, the sampler explores the posterior. These samples are the result.

Three things follow from this:

  • It is stochastic. Two runs are not identical unless you fix the seed.

  • Longer is safer. More iterations explore the posterior more thoroughly.

  • The output is a distribution. You read the collection of samples, not any single one.

Running

sweeplink infer \
    --counts doc_alleleCounts.txt \
    --meta doc_meta.txt \
    --mu_a_A 1e-8 \
    --numThreads 4 \
    --out doc
--counts / --meta

The two files from Step 1.

--mu_a_A 1e-8

Mutation rate from the ancestral allele a to the derived allele A, per generation. 1e-8 suits our human-like simulation; use a rate appropriate to your organism. The reverse rate --mu_A_a defaults to 0.

--numThreads 4

The likelihood is computed per locus, so it parallelises well.

Important

The default is one thread, not all available cores. That is deliberate: the parallel parts of the likelihood do not sum in a fixed order, so a multi-threaded run is not reproducible bit-for-bit even with a fixed seed. With one thread it is. Ask for threads explicitly when you want speed, and drop back to one when you need an exactly repeatable run.

--out doc

Prefix for every output file.

Everything else — the grids, the numerical scheme, the MCMC length — has a sensible default. See Inference for the full list.

What it prints

SweepLink first echoes its entire configuration, then runs burn-in rounds followed by the main chain. Trimmed to the parts worth reading:

   - Will write output files with prefix doc (argument 'out').
   - Will model linkage between neighbouring loci (use 'linkage false' to change).
   - Will use deltaT = 0.1 for the grid on the time points (argument 'deltaT').
   - Will use 50 grid points for the grid on the allele frequencies x (argument 'nGridPoints').
   - Will use Chang-Cooper (CC) numerical scheme to solve PDE (argument 'PDEmethod').
   ...
   - Reading allele counts from file 'doc_alleleCounts.txt'... done!
      - Read 413 loci on 2 chromosomes.
      - Read 11 samples (combinations of time and populations).
   - Will use the following grid with 18 states on s: [-0.07407407407407407, ...,
     -0.009900990099009901, 0, 0, 0.01, 0.02, ..., 0.08] (arguments 'grid_abs_s' or 'abs_s_max').
   - Will use mutation rates of mu(a->A) = 1e-08 and mu(A->a) = 0.
   - Running with a maximum of 4 threads (argument 'numThreads')
   - Will not infer h (specify argument '--h.update true' if you want h inferred).
   - Reading MCMC parameters:
      - Will run an MCMC for 100000 iterations. (parameter 'iterations')
      - Will run 10 burnins of 1000 iterations each.
      - Will write every 10th iteration. (parameter 'thinning')
   - Running MCMC inference:
      - Running 10 burnins:
      1)  Burnin number 1:
         - Running a burnin of 1000 iterations ... done (in 7 seconds)!
           -> Acceptance rate kappa = 0.454773
           -> Acceptance rate nu = 0.312844
           -> Acceptance rate mu = 0.697303
           -> Acceptance rate log10N = 0.609391
         - Adjusting proposal parameters ... done!
      ...
      - Running MCMC chain:
         - Running an MCMC chain of 100000 iterations ... done (in 12 minutes)!
- sweepLink terminated successfully in 13 minutes!

Worth checking before you walk away:

  • Loci and samples. 413 loci on 2 chromosomes and 11 samples should match Step 1. If they do not, the counts or meta file is not what you think.

  • Linkage is on. That is round 1’s whole point; it becomes relevant in Step 4: Refining, where we turn it off.

  • Acceptance rates. kappa, nu and mu govern the linkage model, log10N the population size. SweepLink tunes them between burn-in rounds, so you do not need to touch anything — but rates near 0 or near 1 after tuning are a sign the chain is struggling.

Note

This takes about 12–15 minutes on the example data with 4 threads. Runtime scales with loci, timepoints and iterations.

Tip

Only every 10th sample is written (--thinning), which keeps files small while still leaving 10,000 samples. The state posteriors are accumulated over all 100,000 iterations regardless — thinning only affects the trace files.

The output files

Everything shares the doc prefix. Three of them matter here:

doc_s_statePosteriors.txt

For each locus, the posterior probability of every state on the s grid. This is the main result — everything in Step 3: Results is computed from it.

doc_trace.txt

The MCMC samples of the global parameters, including the population size.

doc.parameters

Every argument the run used, including the ones you left at their defaults. A complete, re-runnable record — pass it straight back to repeat the run:

sweeplink doc.parameters

It is also how the companion scripts discover the run’s settings, which is why they need nothing but --prefix.

The rest — doc_meanVar.txt and doc_config.txt — are described in File formats.

Two things about these files regularly catch people out:

Warning

Read N off log10N, not off N. The MCMC samples log10N, so N is strongly right-skewed and its arithmetic mean sits far out in the upper tail. On this run meanVar reports N_Pop1 as 794,906, while \(10^{\text{mean}(\log_{10}N)}\) is 82,899 — an order of magnitude apart. sweeplink-extract in Step 3: Results does the right thing for you.

Important

The s files store grid indices, not values of s. Index 0 is the most negative coefficient, index 17 the most positive, as printed in the run’s grid line above.

The grid contains two zeros: it is symmetric, and the negative and positive halves each end at \(s = 0\), so with the default grid both index 8 and index 9 mean neutral. Any calculation over the s states has to add the two together — the scripts in Step 3: Results do, and report a single neutral state.

Note

You will not find doc_h_* or doc_z_* files. The dominance coefficient h is held fixed unless you pass --h.update true, so its files would be constant, and z — the linkage model’s internal state — is computed from s and h, so it never holds anything they do not. Neither is written by default; see Inference if you want them anyway.

Next

Plot the scan and turn it into a list of loci under selection in Step 3: Results.