Step 4: Refining¶

Steps 2 and 3 answered which loci are affected by selection. That first round is not the best setting in which to ask how strong that selection is, for two reasons:

N was still being estimated at the same time, from all 413 loci. The selection coefficients and the population size trade off against each other, so each was inferred while the other was moving.
The linkage model pulls each locus toward its neighbours. That is exactly what you want for detection — it is how the evidence from a whole region gets pooled — but it also means a locus’s estimate is not its own.

This step removes both. It runs a second, focused round of inference on the 45 shortlisted loci only, with N held fixed and linkage switched off.

Running it¶

sweeplink-extract already printed the command in Step 3: Results. Copy it, and add --numThreads if you want it faster:

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

Two arguments do the work:

--N 82899 --log10N.update false: Start N at the value the first round estimated, and stop sampling it. Both are needed: --N sets the value, --log10N.update false freezes it. SweepLink samples log10N rather than N, which is why the flag names the log parameter.
--linkage false: Switch off the genome-wide model. This is the important one. The 45 loci are a shortlist scattered across a chromosome; they are no longer neighbours, and the physical distances between them mean nothing. A model whose transition probabilities are built from those distances would be modelling an artefact of how the shortlist was made. With linkage off, every locus is inferred on its own.

   - Will not model linkage between loci: each locus is inferred independently
     (argument 'linkage'). Arguments 'maxDistanceLoci', 'kappa', 'nu' and 'mu'
     have no effect.
   ...
      - Read 45 loci on 1 chromosomes.
   - Setting the update of parameter log10N from the command line. (parameter log10N.update)
- sweepLink terminated successfully in 35 seconds!

35 seconds, against 13 minutes for the first round. 45 loci instead of 413, one parameter fewer, and no linkage model to update.

Note

The burn-in reports no acceptance rates at all, which looks alarming the first time:

1)  Burnin number 1:
   - Running a burnin of 1000 iterations ... done (in 0 seconds)!
   - Adjusting proposal parameters ... done!

Nothing is wrong. Acceptance rates are reported for parameters with a tunable proposal, and in this run every one of them is fixed: log10N by --log10N.update false, and kappa, nu and mu — the linkage model’s own parameters — automatically, because with linkage off there is nothing left to inform them. s is still being sampled, but it moves on a discrete grid, which has nothing to tune. doc_r2_config.txt records all four as not updated.

The final figure¶

The same plotting script, with the second run added:

sweeplink-plot --prefix doc --round2-prefix doc_r2

The same two-panel figure, now with second-round estimates overlaid and joined to the first-round estimates, both in red. — Both rounds. Hollow red dots are the first round, filled red are the second, and the thin line joins the two estimates of the same locus.¶

The top panel is unchanged — it is the first round’s evidence, and what selected the loci in the first place. The bottom panel now carries both estimates, so you can see how far each one moved.

What changed¶

	First round	Second round
mean \(E\|s\|\)	0.0410	0.0479
range	0.0221 – 0.0496	0.0320 – 0.0589
`N`	inferred, sd 0.67 on \(\log_{10}\)	fixed at 82,899

That is the summary; --compare plots every one of the 45 loci individually, first round against second, against the line where the two would agree:

sweeplink-plot --prefix doc --round2-prefix doc_r2 --compare

This writes the usual scan figure plus a second file, doc_selection_comparison.png:

Scatter plot of round 1 posterior mean |s| against round 2, one point per locus, with a dashed y=x line. Almost every point sits above the line. — Round 1 vs round 2, one dot per locus. Above the dashed line means the second round’s estimate is larger.¶

82% of the loci moved up, and that direction is not a coincidence — it is the predicted effect of removing the linkage model, and our experiments confirm it happens in practice. With linkage on, every locus is pulled toward the neutral attractor it shares with its neighbours, which shrinks estimates toward zero; the more a locus’s neighbours look neutral, the more it gets pulled down even if its own signal is strong. Remove that shrinkage and the per-locus estimates rise to what the data at each locus alone supports.

Which of the two is “right” depends on what you are asking. The first round’s estimates are the better description of a region under selection, because they use the information that the region shares. The second round’s are the better estimate for an individual locus, because they are not borrowing from anywhere else.

Warning

The second round’s uncertainty is understated, because N is treated as known when the first round only pinned it to within a factor of about five. The point estimates are the useful output here; the credible intervals around them are narrower than the analysis as a whole justifies.

Ranking the shortlist¶

45 loci is still a list, not an answer — the next question is usually which one first. The second round’s per-locus estimates, freed from the linkage model’s shrinkage, are the right numbers to rank by.

There is no dedicated “rank” command: sweeplink-extract already computes a signed mean_s for every locus it is given, so pointing it at the second round itself, with a threshold loose enough to keep all 45, reuses it as one:

sweeplink-extract --prefix doc_r2 --out-prefix doc_r2_ranked --max-samples-against 100000

doc_r2_ranked_selected.txt now holds all 45 loci with mean_s computed from the second round’s own posterior. Sort it:

tail -n +2 doc_r2_ranked_selected.txt | sort -t$'\t' -k4,4 -g -r | head -6

  63107   positive        0.0589262       100000  0       0.99999 doc_r2
  50001   positive        0.0583156       100000  0       0.99999 doc_r2
  55120   positive        0.0561655       100000  0       0.99999 doc_r2
  57525   positive        0.0549678       100000  0       0.99999 doc_r2
  51471   positive        0.054   100000  0       0.99999 doc_r2
  71549   positive        0.0511908       100000  0       0.99999 doc_r2

(sort -k4,4 -g without -r gives the strongest negative-direction loci instead, at the other end of the table.)

The simulated sweep sits at chr2:50001 with a true \(s\) of 0.05 — and it comes out second of 45, with a refined estimate of 0.0583. That is the payoff of the whole two-round procedure: the locus you would actually want to follow up on is right there near the top, not buried in the shortlist.

Note

We can check two files to see a comparison of selection values for this locus from round 1 and round 2:

chr2:50001	round 1	round 2
file	`doc_r2_selected.txt` (Step 3: Results)	`doc_r2_ranked_selected.txt` (above)
`mean_s`	0.0478024	0.0583156

A +22% move for this locus — more than the +17% the mean shifted by. The same two files, read together, give you this comparison for any locus on the shortlist, not just this one.

Annotating the figure¶

Once you know which positions matter, the next step is usually to look them up — nearby genes, prior GWAS hits, whatever your own annotation source is — and put the name on the figure rather than a bare coordinate.

--labels takes a tab-separated file with chr, pos and label columns — keep the label short, since it has to share the figure with the data:

chr  pos     label
2    50001   true sweep

The two panels then say two different things about the same locus: the top gets your annotation, the bottom gets the locus’s own chr:pos, worked out for you.

Since this dataset is simulated, the annotation is the one thing a real analysis would not have in advance. Passing it in shows what one looks like on the figure, and doubles as a check that the procedure really does land on this locus:

sweeplink-plot --prefix doc --round2-prefix doc_r2 --labels doc_true_variant.txt

The same two-round figure, with the simulated sweep's true locus labelled: "true sweep" above the top panel, "chr2:50,001" next to its marker in the bottom panel, joined by a leader line through both. — The locus, labelled. The leader line spans both panels; the annotation sits above the top one, italic, the position sits next to the marker in the bottom one, bold.¶

With real data, --labels works the same way once you have your own annotation file — one row per locus you want named, built from whatever lookup you use (a gene database, a table of previously reported hits, …). A curated handful of peaks is what it is built for; see sweeplink-plot for the full option reference.

Where to go next¶

doc_r2_selected.txt — the loci, their direction, a signed posterior mean of s (mean_s), and their support.
doc_r2_s_statePosteriors.txt — the full refined posterior over s at each of them. Remember that the states are grid indices, and that indices 8 and 9 both mean neutral (see Step 2: Inference).
doc_r2.parameters — the complete record of the run. Feed it back to sweepLink to repeat it exactly.

For every option of the inference itself, see Inference.