pickref

Rank a set of bacterial genomes and pick the best reference genome for downstream analyses (core-SNP calling, read mapping, pangenomics).

pickref replicates the dRep d_choose.py scoring formula and adds a Mash all-vs-all centrality term: the mean Mash-ANI of a genome to every other genome in the input set. The highest-scoring genome is the most representative, highest-quality choice.

Why

When you have thousands of closely related genomes (e.g. a serotype, an ST, an outbreak cluster) you need one reference to map reads or call core SNPs against. Picking it by hand is arbitrary; picking the best assembly quality alone ignores whether it is central to the population. pickref combines:

• Assembly quality (CheckM completeness / contamination / strain heterogeneity, N50, total length)
• Population centrality (average Mash-ANI to all other genomes in the set)

Scoring formula

score = completeness              # × W_com
      − contamination             # × W_con (penalty)
      + W_str × contamination × (strain_het / 100)   # strain-heterogeneity reprieve
      + log10(N50)                # × W_n50
      + log10(Length)             # × W_size
      + (centrality − 0.99)       # × W_cent
      + extra

All weights default to dRep's values (1.0, S_ani = 0.99, extra = 0) and can be overridden from the command line.

Note on strain heterogeneity: the + W_str × contamination × (het/100) term reduces the contamination penalty when the contaminating markers come from the same strain (high heterogeneity). In many real datasets het is 0.0 for most genomes, so this term has no effect — contamination is then penalised in full.

Requirements

• Python ≥ 3.9 with numpy and pandas
• mash on PATH (or pass --mash)

Install Python deps:

pip install -r requirements.txt

Input format

Three tab-separated inputs:

File	Columns
--list	one FASTA path per line
--quality	assembly, N50, Length (extra columns allowed; assembly matches the FASTA basename or its GCA_/GCF_ prefix)
--checkm	CheckM output TSV: first column bin_id, plus Completeness, Contamination, Strain heterogeneity

pickref joins --quality to --checkm either by full bin_id or by GCA_/GCF_ short prefix, so the two tables do not need identical IDs.

Usage

# rank a set of genomes
python3 pickref.py \
  --list    genomes.list \
  --quality assembly_quality.tsv \
  --checkm  checkm.tsv \
  -o        pickref_scores.tsv

Output: ranked TSV sorted by score descending — take the first row as your reference genome.

Options

Option	Default	Description
--list FILE	—	text file with one FASTA path per line (required)
--quality FILE	—	quality TSV: assembly, N50, Length (required)
--checkm FILE	—	CheckM output TSV (required)
-o, --out FILE	pickref_scores.tsv	output ranked table
--mash PATH	mash from PATH	mash executable
--sketch-size N	10000	mash sketch size
-t, --threads N	8	threads per mash dist call
--sketch-dir DIR	<out>/../sketch	cache dir for .msh sketches
--pairwise-out FILE	—	also write per-genome mean-ANI table
--limit N	—	only score the first N genomes (testing)
dRep weights	see formula	--w-com, --w-con, --w-str, --w-n50, --w-size, --w-cent, --s-ani, --extra

Run python3 pickref.py --help for details.

Try it now (toy example)

# regenerate the 3 synthetic genomes + input tables
python3 examples/make_demo_genomes.py

# rank them
python3 pickref.py \
  --list    examples/genomes.list \
  --quality examples/assembly_quality.tsv \
  --checkm  examples/checkm.tsv \
  -o        /tmp/pickref_demo.tsv

genome_2 is a 5%-diverged copy of genome_1 with strain_het = 100%; genome_3 is unrelated. After the run, genome_2 ranks first thanks to its het reprieve.

How it works

1. Sketch every FASTA with mash sketch.
2. Compute all-vs-all distances with per-genome rounds of mash dist (memory-friendly: only one genome's distances are held at a time).
3. Centrality of genome i = mean of (1 − Mash distance) to all others.
4. Score by the dRep formula, rank descending, write TSV.

Caveats

• All-vs-all Mash distance needs O(n) mash dist round trips; for tens of thousands of genomes expect hours (parallelise by splitting the list, or raise --threads).
• --list FASTA basenames must map to CheckM bin_id (or share a GCA_/GCF_ prefix).

License

MIT — see LICENSE.

Credits

Inspired by dRep (Olm et al. 2017, PeerJ) — scoring formula replicated from d_choose.py.