How it works

This page explains the model well enough to make the options make sense. It is not a substitute for the manuscript, and you do not need it to run the tutorial.

Figure 1 of the manuscript. A — the two layers. B — what the spatial layer buys you: isolated single-locus signals are suppressed, and a genuine sweep shows up as a contiguous region.

SweepLink models the data in two directions at once, and almost every option belongs to one of them:

• the temporal layer, down the figure — how the frequency of one variant changes through time, and how the counts you observe arise from it;

• the spatial layer, across the figure — how neighbouring loci along a chromosome inform each other.

Single-locus methods have only the first. The second is what the name refers to.

The temporal layer

For one locus, the quantity of interest is the frequency of the derived allele, \(x\), at each sampled time point: \(X_0, X_1, \dots X_K\) in the figure.

How the frequency moves. Between two time points the frequency evolves under the Wright-Fisher diffusion, written \(p(\Delta t, N, m)\) in the figure. It depends on the gap in generations, the effective population size \(N\) (larger \(N\) = less drift, since the drift variance is \(\frac{x(1-x)}{2N}\)), the migration matrix \(m\) if there is more than one population, and — this is the point — the selection coefficient \(s\) at that locus.

SweepLink does not simulate this; it solves the diffusion equation numerically on a grid of allele frequencies. That is where the numerical options come from: --nGridPoints and the grid-shape flags set how finely x is discretised, --deltaT sets the time step, and --PDEmethod chooses the solver.

How the counts arise. You never observe \(x\), only a sample: \(c\) copies of the derived allele out of \(n\) chromosomes. That step is a binomial draw, and it is the only place your data enter. Two consequences worth knowing:

• Sample sizes may differ freely between time points, and \(n = 0\) simply means “no observation here” — which is why missing data is written 0/0.

• Nothing about haplotypes or phasing is used. The model never sees more than a count per locus per time point per population. It does, though, assume that \(c\) counts the derived allele — the priors below and the sign of \(s\) both depend on it.

--emissionModel controls whether the binomial density is evaluated at a grid point or averaged over the frequency interval that grid point represents.

The spatial layer

Now the part that distinguishes sweepLink. Each locus carries a hidden state — \(s_1 \dots s_L\) in the figure, and z in the output files — that encodes how strongly it is being selected. These states form a Markov chain along the chromosome, with transition probabilities \(Q(d)\) that depend on the physical distance \(d\) between consecutive loci.

The effect is that neighbouring loci are not independent. A locus close to a strongly selected one is a priori more likely to be selected itself, and the closer it is the stronger that pull. This is the hitchhiking effect written into the prior.

Two things follow, both visible in panel B of the figure:

• A sweep is detected from a region, not a single variant. Evidence that would be unconvincing at any one locus adds up across the region, so weak sweeps are recoverable.

• Isolated blips are suppressed. A single locus that drifted to an extreme frequency, with unremarkable neighbours, is pulled back toward neutrality — so it does not become a false positive.

Transition matrix \(Q(d)\) is parameterised by kappa, nu and mu, which are estimated along with everything else. Distances are binned, up to --maxDistanceLoci.

Important

This layer is only meaningful when the loci really are neighbours. If you re-analyse a shortlist scattered across a chromosome, the distances between the survivors are an artefact of how the shortlist was made, and the model would be pooling evidence across gaps that are not there. That is what --linkage false is for, and why the tutorial’s second round uses it.

Note also that mu here is a transition parameter of this chain — not a mutation rate. The mutation rates are --mu_a_A and --mu_A_a.

The grid on s

Selection is not sampled continuously. It lives on a discrete grid, which is what makes the state posteriors in the output files a table of probabilities per state rather than a density.

The grid is symmetric around zero: the positive values you set with --abs_s_max or --grid_abs_s, and their mirrored negatives. Symmetry is deliberate — it makes a deleterious variant exactly as reachable as a beneficial one, so the method is not biased toward finding positive selection.

Important

Two consequences that catch people out:

1. The output files store grid indices, not values of s. Map them back with grid_abs_s from the .parameters file.

2. The grid contains two zeros — the negative half and the positive half each end at s = 0. With the default 18 states, both index 8 and index 9 mean neutral, and any sum over states must add the two together.

Fitness is \(w_{AA} = 1+s\), \(w_{Aa} = 1+hs\), \(w_{aa} = 1\), so s is the selection coefficient of the derived allele and h its dominance. h is shared across populations and is not inferred unless you pass --h.update true; s is per locus and per population.

Priors on the starting frequency

Solving the diffusion forward in time needs somewhere to start: a prior on the allele frequency at the first time point. A variant that was already segregating when sampling began is not in the same situation as one that appeared later, so the two get separate priors — --forwardPrior for standing variants, --newMutPrior for new mutations.

A locus counts as standing if its first non-missing time point is polymorphic.

What comes out

The MCMC returns full posterior distributions, not point estimates:

• for each locus, a posterior over the states of s — the main result, and what the “probability of selection” in the figures is computed from;

• for each population, a posterior over log10N (and migration, if there is more than one population);

• the parameters of the spatial layer, kappa, nu and mu.

Because the answer is a distribution, the natural question is not “what is s here” but “how much of the posterior agrees that s is non-zero” — which is exactly the criterion sweeplink-extract applies.

Where the likelihood tasks fit

Everything above describes infer, which samples from the posterior. The four calculateLL* tasks do something narrower: they evaluate the likelihood directly, at fixed parameter values, without any MCMC.

That is useful when you want to see the shape of the likelihood rather than a posterior summary — most commonly to check where your data put N before committing to a prior for it.