Step 3: Results

Step 2’s inference has finished. Two things to do with it: look at the signal along the genome, and turn it into a list of loci under selection.

Note

What counts as “under selection”? By default, every one of the 100,000 MCMC samples must agree on the sign of s — no sample fell in a neutral state. That is deliberately strict; later on this page shows how to relax it and why the default is what it is.

The evidence and the landscape of selection

sweeplink-plot --prefix doc

That writes doc_selection.png. Nothing else is needed — the script reads the run’s settings out of doc.parameters.

Our scan. Chromosome 1 is flat, i.e. is neutral; everything is on chromosome 2, which is where the sweep was simulated.

Our scan. Chromosome 1 is flat, i.e. is neutral; everything is on chromosome 2, which is where the sweep was simulated.

The top panel is the evidence. Background loci are coloured only by which chromosome they are on (alternating grey and blue — no other meaning to that colour). The threshold is a red dashed line, all loci that are equal or above it, should be considered as “affected by selection”.

The bottom panel is the effect size of selection. The shaded area is the posterior mean of \(|s|\) at every locus — the magnitude of selection, tinted to match each chromosome’s colour above. Red dots mark the loci that passed the threshold.

Read together: the top panel says where there is something, the bottom says how strong it is.

Why the evidence axis is not a plain probability

The natural thing to plot would be the posterior probability itself. It does not work, because the posterior saturates:

50th percentile   0.00618
75th percentile   0.85172
90th percentile   0.99999
95th percentile   0.99999

The top 10% of loci all sit at the same value. With 100,000 MCMC samples, a locus where every sample agrees cannot be distinguished from another where every sample agrees — the resolution stops at one sample in 100,000. On a plain 0–1 axis they would all be crushed into the top sliver.

So the axis plots \(-\log_{10}(1-P)\) directly, which spreads the interesting range out into an ordinary, readable scale instead of crushing it into a sliver.

This is also why the reported posterior is never exactly 1. One pseudo-count is held on the neutral state, so a locus supported by all 100,000 samples reads \(100000/100001 = 0.99999\), not 1 — which keeps the axis finite instead of needing to plot infinity.

Tip

--region 2 zooms both panels to one chromosome, and --show opens the figure in an interactive window instead of writing a file. --out doc.svg gives a vector figure — best for a dataset this size, and what the figures on this page are. On a full genome, use PNG with --dpi 300 instead: a vector file stores one element per locus, so it grows with the data. For the rest — themes, direction markers — see sweeplink-plot.

The population size

We can also plot a posterior distribution of the inferred population size:

sweeplink-plot-posteriors --prefix doc

N:
  Pop1         82,899  [10**mean(log10 N)]
    median 58,653   IQR 27,781 - 175,700   95% 10,884 - 3,816,212
    range  4,821 - 265,842,000   sd 0.6618 (on log10)   n = 10,000

The posterior of N. Log10(N) by default — pass --linear for a raw N axis, though the tail above makes that much harder to read.

The posterior of population size N. It is in log-scale by default — pass --linear for a linear axis, though the tail above makes that much harder to read.

82,899 — the mean on the natural scale, 10^mean(log10 N) — is what the second round below fixes N at.

Extract loci under selection

sweeplink-extract --prefix doc --out-prefix doc_r2

Read 413 loci from 'doc_s_statePosteriors.txt' (100000 MCMC samples each;
selecting where at most 0 disagree, posterior >= 0.99999).

Selected 45 of 413 loci (18 positive, 27 negative).
  doc_r2_alleleCounts.txt
  doc_r2_meta.txt
  doc_r2_selected.txt

Second round: N fixed at the estimate of the first run, linkage off.
  Pop1: N = 82899 (from 10000 pooled samples)

  sweeplink infer \
      --counts doc_r2_alleleCounts.txt \
      --meta doc_r2_meta.txt \
      --mu_a_A 1e-8 \
      --N 82899 --log10N.update false \
      --linkage false \
      --out doc_r2

45 of 413 loci, and every one of them on chromosome 2 — which is where the sweep was simulated. Not one locus on chromosome 1 passed.

doc_r2_selected.txt lists them, including a posterior mean of s itself (signed — negative for the negative-direction loci):

chr  pos     direction       mean_s  samples_supporting      samples_against posterior       run
2    8126    negative        -0.0220707      100000  0       0.99999 doc
2    11896   negative        -0.0302666      100000  0       0.99999 doc
2    12639   negative        -0.0312808      100000  0       0.99999 doc
...

Note

Both directions show up around a sweep, and that is expected rather than a problem. As the selected haplotype rises it drags linked variants with it; a variant that happens to sit on the other haplotype is driven down, and shows up as negative selection. That is why 27 of the 45 are negative. The two-round procedure sharpens the estimates but does not try to separate the causal variant from its hitchhikers.

What counts as “under selection”

The default is strict: every MCMC sample must agree on the sign of s. No sample fell in a neutral state. Our experiments found this gives the best accuracy, so it is what you get with no arguments.

To relax it, say how many samples may disagree:

sweeplink-extract --prefix doc --out-prefix doc_loose --max-samples-against 5000

That keeps 85 loci instead of 45. The same threshold can also be written as a probability (--min-posterior 0.95) or in the evidence panel’s own units (--min-evidence 1.3) — all three select exactly the same 85 loci here. Pass whichever you used to sweeplink-plot as well, so the line in the figure matches the shortlist.

sweeplink-extract covers all three forms, the --direction filter, and the traps in --min-posterior.

Why N is fixed at 82,899

Note the line about N:

Pop1: N = 82899 (from 10000 pooled samples)

That value is \(10^{\text{mean}(\log_{10}N)}\) over the run’s MCMC samples, not the mean of N itself. The MCMC samples log10N, so N is strongly right-skewed and its arithmetic mean sits far out in the tail — doc_meanVar.txt reports N_Pop1 as 794,906, nearly ten times the 82,899 used here.

Warning

The second round treats N as a known, fixed value. The first round’s own posterior for N is genuinely uncertain — see Step 4: Refining for how much — so the second round’s estimates of s look more precise than the analysis as a whole justifies. Keep that in mind when reading them.

We already saw the shape of that skew above: the tail reaching past 10,000,000 is exactly what drags the mean so far from where almost all of the posterior actually sits.

sweeplink-plot-posteriors becomes more useful once you have several runs to compare — pass --prefix (or --dirs) more than one value, as you would for a dataset split by chromosome, and add --per-chromosome to see each run’s own violin alongside the pooled one, which is a quick way to check whether they actually agree. --params kappa nu mu (or betaMean betaESS) adds the linkage model’s own parameters to the figure — off by default, and skipped automatically wherever the run’s settings held them fixed (as they are here, with linkage on and a uniform prior).

Next

The three files written by sweeplink-extract are the input for Step 4: Refining.