tutorial.md

Tutorial

A complete run on a small dataset with known problems planted in it, so you can check that what Recore reports is what is actually there.

You need minimap2 on PATH and Recore installed (see Installation). Nothing is downloaded: the dataset is generated by a script in the repository.

1. Build the dataset

python3 benchmarks/make_tutorial_dataset.py --out tutorial_data

genomes            : 12
clusters           : 24
planted absences   : 11
written to         : tutorial_data

Twelve MTBC-style genomes across six lineages, with an L8 outgroup, and twenty-four gene clusters:

tutorial_data/
├── gene_presence_absence.csv   # what the pangenome pipeline thought it saw
├── pan_genome_reference.fa     # one representative sequence per cluster
├── genomes/                    # one assembly per genome
├── annotations/                # one GFF3 per genome
└── species.tree                # Newick, matching the lineage structure

Eleven cells are marked absent in the matrix that should not be. They are there on purpose, in four flavours:

Cluster	Genomes	What was planted
Rv0003	3	the full-length gene is in the assembly, the matrix says absent
Rv0012	3	the same, in a different cluster
Rv0006	2	the gene is in the assembly, but the GFF has no CDS there
Rv0009	3 (all of L2)	the gene is there with an internal stop codon

Plus three cases that must not be corrected: Rv3810 is genuinely absent everywhere but L8, Rv0018c has only a diverged copy in L4_AS00404, and Rv1180/Rv1181 are carried by L2 alone.

2. Run Recore

recore \
  -p tutorial_data/gene_presence_absence.csv \
  -r tutorial_data/pan_genome_reference.fa \
  -g tutorial_data/genomes \
  -a tutorial_data/annotations \
  --synteny \
  --multi-allele \
  --matrix-policy balanced \
  --reference-genome L4_AS00401 \
  --tree tutorial_data/species.tree \
  --write-html-report \
  -t 4 \
  -o tutorial_out/

The lines that matter in the log:

Gaps: 50 total, 39 after prevalence filter (≥10%), 11 skipped
Core genome (before): 14 / 24 clusters
Final rescues: 11 (DROPPED: 6, MISSING_CDS: 2, SPLIT: 0, FRAGMENTED: 0,
                   ALIGNMENT_SPLIT: 0, START_VARIANT: 0, PARALOG: 0, PSEUDOGENE: 3)
Core genome: 14 → 17 (+3)

Read that against what was planted:

• 50 gaps is every empty cell in the matrix, not just the planted ones.
• 11 skipped are the Rv3810 cells. It sits in one genome out of twelve, which is under the 10% --min-prevalence floor, so it is never searched.
• 11 rescues is exactly the 11 planted cells, split into DROPPED 6 (Rv0003 and Rv0012), MISSING_CDS 2 (Rv0006) and PSEUDOGENE 3 (Rv0009).
• PARALOG: 0: the diverged copy in L4_AS00404 is 85% identical, under the 90% --min-identity bar, so it produced nothing at all. It must not, and it did not.

3. Eleven found, eight corrected

grep -E "^Rescued:|^Matrix-corrected:" tutorial_out/summary.txt

Rescued:               11
Matrix-corrected:      8

Eleven cells had evidence; eight reached the binary matrix. The gap is the whole point of the balanced policy: a gene with an internal stop codon is evidence of a disrupted homologue, not a clean presence call, so Rv0009 stays out of the binary matrix and stays visible everywhere else.

orthology_evidence.tsv says so per cell:

awk -F'\t' 'NR==1{for(i=1;i<=NF;i++) h[$i]=i; next}
            {print $h["gene_cluster"], $h["rescue_type"],
                   $h["evidence_state"], $h["matrix_action"]}' \
    tutorial_out/orthology_evidence.tsv | sort | uniq -c

gene_cluster	rescue_type	evidence_state	matrix_action
Rv0003 (x3)	DROPPED	STRICT_PRESENT	correct_presence
Rv0012 (x3)	DROPPED	STRICT_PRESENT	correct_presence
Rv0006 (x2)	MISSING_CDS	STRICT_PRESENT	correct_presence
Rv0009 (x3)	PSEUDOGENE:1stops	PSEUDOGENE	audit_only_policy

STRICT_PRESENT on Rv0006 is worth a second look. The gene was in the assembly with no CDS annotated over it, which on its own is only MISSING_CDS_UNVALIDATED. --synteny checked the flanking genes, found the locus where it should be, and that promoted the call.

4. What the policy and the evidence each decide

Run the same data six ways and the two levers separate cleanly:

for policy in strict balanced sensitive; do
  for flags in "--synteny" ""; do
    recore -p tutorial_data/gene_presence_absence.csv \
           -r tutorial_data/pan_genome_reference.fa \
           -g tutorial_data/genomes -a tutorial_data/annotations \
           $flags --multi-allele --matrix-policy "$policy" \
           -t 4 -o "out_${policy}" >/dev/null 2>&1
    printf "%-10s %-8s " "$policy" "${flags:---none}"
    grep "^Matrix-corrected:" "out_${policy}/summary.txt"
  done
done

Policy	With --synteny	Without
strict	8	0
balanced	8	6
sensitive	11	9

Three things fall out of that table:

• strict without annotation context corrects nothing. It only accepts same-locus evidence, and without synteny none of the eight correctable cells has any: six come through as SEQUENCE_ONLY and the two Rv0006 cells as MISSING_CDS_UNVALIDATED, and strict takes neither. (The other three are PSEUDOGENE either way.) Strict is not "balanced but fussier", it is a different question.
• balanced takes full-length sequence, but not an unconfirmed locus. The six Rv0003 and Rv0012 cells sit on annotated CDSs, so they are corrected either way; what changes is the provenance, STRICT_PRESENT with synteny and SEQUENCE_ONLY without. The two Rv0006 cells are a gene with no CDS annotated over it. Only synteny can say the locus is where the gene belongs, so without it balanced holds them back as audit-only.
• sensitive with synteny reaches 11. The extra three are the Rv0009 pseudogenes. That is the policy saying "a disrupted homologue counts as presence", which is a claim you have to be willing to make in writing. Without synteny it still takes them, and still holds back Rv0006.

5. Read the multistate matrix

The binary matrix is what downstream tools consume. The multistate matrix is what you read:

grep -E "^(Rv0003|Rv0009|Rv3810)," tutorial_out/gene_presence_absence_recore_states.csv \
  | cut -d, -f1,15-17

Rv0003,STRICT_PRESENT,STRICT_PRESENT,STRICT_PRESENT
Rv0009,PRESENT_INPUT,PRESENT_INPUT,PSEUDOGENE
Rv3810,NOT_SEARCHED_LOW_PREVALENCE,NOT_SEARCHED_LOW_PREVALENCE,NOT_SEARCHED_LOW_PREVALENCE

Those three columns are L1_AS00101, L1_AS00102 and L2_AS00201. Rv0003 shows STRICT_PRESENT where its false absences were planted; Rv0009 is untouched in the two L1 genomes and PSEUDOGENE in the first L2 one. Rv3810 was never searched, being under the prevalence floor, and says so instead of reading as an absence Recore had checked. Nothing was thrown away and nothing was quietly upgraded. The thirteen states are listed in Evidence states.

6. Open the report

Open tutorial_out/recore_report.html in a browser. It is one self-contained file, so double-clicking it works (checked in Chrome and in Safari's engine), and so does copying it off a cluster.

Things to try, in order:

1. Press Cmd+K and type Rv0009. The matrix jumps to that column. The three L2 genomes are the only ones not marked present, and the inspector shows the mechanism.
2. Set tip colour to lineage from the sliders icon in the phylogeny header. The L2 clade lights up as one block, which is what a clade-restricted signal looks like.
3. Find Rv1180. Only the three L2 tips carry it. This one is a genuine clade acquisition rather than an artefact, and the inspector says so.
4. Press S to open the synteny band. The flanking genes are stable across genomes, which is why the Rv0006 locus could be validated.
5. Press ? for the rest of the shortcuts.

The Run report page covers the viewer properly.

7. What changed, in one line

grep "Core genome" tutorial_out/summary.txt

Core genome (strict):  14 → 17 (+3)

Three clusters that were missing from at least one genome are now found in all twelve. They were always there; the matrix was wrong about them. That is the entire job.

Next

• Matrix policies, to decide which policy your own run should use.
• Output files, for the columns of everything in tutorial_out/.
• Inputs, for getting your own data to match up.