Skip to content

Search every genome with the same proteins for the pangenome - #69

Merged
Paururo merged 2 commits into
mainfrom
pangene-shared-genes
Oct 24, 2026
Merged

Paururo merged 2 commits into
mainfrom
pangene-shared-genes

Conversation

@Paururo

@Paururo Paururo commented Oct 24, 2026 •

Copy link
Copy Markdown
Member

The bug

PANGENE aligned each sample's own corrected proteins to every other sample's genome and ran pangene once per sample, so each Rtab was (this sample's genes) x (the other genomes), its columns named <genome>_<sample>#0. REPORT_INPUTS then outer-joined the Rtabs on gene names that carry the sample's name and coordinates: A4_Mbovis_Rv0001;+;0-1524 and L4_CDC1551_Rv0001;+;0-1524 were two rows. With N samples the matrix had N(N-1) columns, a gene could be in at most N-1 of them, and no gene was ever core. A sample's own genome was never searched either, so a gene only one genome has never entered the matrix: the interactive report gave every genome 0 genes "in no other genome".

The fix

pangene is meant to be given one set of proteins aligned to every genome. PANGENE is now one task for the run:

  1. bin/pangenome.py proteins writes the set: the reference proteome (H37Rv or a custom one), each protein named <gene>:<protein> so that pangene calls the gene by its locus (Rv1506c, b0002), and each sample's genes that no reference gene matches (an empty rvs in <sample>.AMAP.tsv), named after the sample's gene. A sequence already in the set is not added again.
  2. miniprot aligns the set to every sample's genome, the sample's own included: N alignments where there were N(N-1).
  3. One pangene run, then gfa2matrix, call and gfatools stat.
  4. bin/pangenome.py matrix names each column after its sample (L4_CDC1551, no #0) and gives a genome none of the proteins was found in a column of zeros, which pangene leaves out.

Of the four designs measured on real data before coding (reference genes only, every sample's proteins pooled, reference plus the genes outside it, per-sample runs merged by reference name), this is the third: rows a reader can look up, genes outside the reference kept and matched across samples by pangene, no new tool.

  • Outputs are PANGENE/pangenome.{gfa,bubbles.txt,Rtab,stat}, one set for the run instead of PANGENE/{sample}/. .stat was never published before: no output declared it.
  • REPORT_INPUTS still passes the matrix through merge_rtabs.py, which leaves a single matrix as it is (byte for byte in every run below). final_report.html reads it unchanged; the pangenome card of amap_report.html now reads a column per genome.
  • PANGENE_GENOME (Build the pangenome without lineage typing聽#48) stays: under --skip_typing it still gives each genome of the sheet its own <sampleId>.fasta.
  • Docs (outputs, workflow, the website's PANGENE station) and CHANGELOG updated, pages rebuilt with tools/build_docs.py.

Before and after, real runs

Typed runs under -profile docker with Nextflow 26.04.6, eggNOG-mapper replaced by a stand-in (its database is about 45 GB, and nothing here reads it); main is 6351c0c.

Run main this branch
CDC1551 and M. bovis, first 2 Mb 3539 rows x 2 columns, each row in one; report "2 samples, 3539 gene families", 100% shell 1809 genes x A4_Mbovis, L4_CDC1551: 1768 core, 17 only in M. bovis, 24 only in CDC1551
the same and H37Rv 5379 rows x 6 columns, none in all six; report "6 samples, 5379 gene families", 100% shell, heatmap of 200 unnamed M. bovis rows 1810 genes x 3 samples; report "3 samples, 1810 gene families", 1760 core (97.2%), 50 shell, heatmap of the 50 by name
E. coli K-12 (first 400 kb) as custom reference, K-12 and Sakai (first 460 kb), --skip_typing, both files called assembly.fna 603 rows x 2 columns (Sakai_K12#0, K12_Sakai#0), none in both 465 genes x K12, Sakai: 296 in both, 63 only in K-12, 106 only in Sakai

The 50 variable genes of the three-genome run, by where they are missing:

  • from M. bovis only (24): RD4 (Rv1506c to Rv1516c), RD10 (Rv1255c to Rv1257c), Rv0095c, Rv0221, Rv0222, Rv0336, Rv1197, Rv2106, Rv3023c, Rv3325
  • from CDC1551 only (16): the phiRv1 prophage (Rv1573 to Rv1586c), Rv0793, Rv1199c
  • from H37Rv only (7): three genes of the RvD2 region that CDC1551 and M. bovis carry (A4_Mbovis_01833, A4_Mbovis_01834, L4_CDC1551_01839), L4_CDC1551_00745, and Rv1767, Rv1768, Rv3724A past the end of the H37Rv window
  • in one genome only: Rv1769 and Rv1770 in M. bovis (window edge), L4_CDC1551_01396 (an Rv1320c-like protein) in CDC1551

The pangenome card of amap_report.html, same run: on main every genome had 0 genes in no other genome and its own count in every other one (1761, 1762, 1764); now 1760 for all three, and M. bovis has 2 genes of its own and CDC1551 1.

PANGENE took 3.6 s and 382 MB for the three 2 Mb genomes. A reference gene can still be in no column: in the E. coli run, 15 of K-12's 374, six of them proteins of 10 to 24 residues, which miniprot does not place.

Tests

753 tests, up from 744.

  • tests/test_pangene_genomes.py is rewritten. It renders PANGENE_GENOME's and PANGENE's scripts as Nextflow does and runs them on tests/data/pangenome/: real genomes of CDC1551, M. bovis and H37Rv cut around RD4 and RvD2 out of the three-genome run above, rebuilt byte for byte by its make_fixture.py. Stand-ins for miniprot, pangene, k8 and gfatools hand names on as the real ones do. It checks that:

    • every genome, the sample's own included, is searched with the same proteins;
    • the set is the reference proteome plus the genes outside it, with a protein two samples share once;
    • there is a column per genome, typed or not, and a genome with no hit keeps a column of zeros;
    • one genome works, and so do as many genomes as jobs at once (the bash macOS has as /bin/bash calls an empty array unbound under set -u);
    • what PANGENE publishes is what it declares, and bin/pangenome.py behaves on its own.

    Where the real tools are installed, as in the container, RD4 must come out missing from M. bovis only and the three RvD2 genes in CDC1551 and M. bovis only, one row each.

  • tests/test_report_tables.py, tests/test_amap_report.py and tests/test_sample_ids.py give the reports the matrix in its new shape. Where R is installed, the report must count the genes all the samples share.

  • Each part broken on purpose fails 2 to 5 tests: no guard on the last wait, #0 kept, no column of zeros, a shared sequence added twice, the sample's own genome left out.

Ran:

  • Python 3.11 here: 711 passed, 42 skipped.
  • In the container (Python 3.9, R 4.4.3, the real tools): 713 passed, 40 skipped.
  • With R off the PATH, as in CI, the R parts skip.
  • nextflow lint (26.10.2): no errors. The only new warnings are the projectDir ones every module has, and main.nf has two fewer.
  • -preview is clean on 24.04.4 and 26.10.2, typed and with --skip_typing.
  • tools/build_docs.py --check and tools/check_site.py pass.

CI

The second commit runs the Python 3.9 job on ubuntu-24.04. GitHub's ubuntu-latest became Ubuntu 26.04, for which setup-python has no Python 3.9, so on this PR the job failed at setup before running a test, and it would on any other: main's last run, on 20 October, was still on 24.04. 3.9 is the container's Python; 3.11 stays on ubuntu-latest.

PANGENE aligned each sample's own proteins to the other samples' genomes and the report merged the matrices on gene names that carry the sample, so no gene was ever in two samples and nothing was core. It now aligns the reference proteome and each sample's genes that no reference gene matches to every genome, the sample's own included, runs pangene once and names each column after its sample. On the first 2 Mb of CDC1551, M. bovis and H37Rv the report goes from 5379 families in 6 columns, none core, to 1810 genes in 3 samples, 1760 of them core.
@Paururo Paururo added bug Something isn't working documentation Improvements or additions to documentation labels Oct 24, 2026
GitHub's ubuntu-latest became Ubuntu 26.04, for which setup-python has no Python 3.9, so the job stopped at setup before running a test. Python 3.9 is the container's, so its job now runs on 24.04 and 3.11 stays on the latest.
@Paururo
Paururo merged commit 1a1b395 into main Oct 24, 2026
3 checks passed
@Paururo
Paururo deleted the pangene-shared-genes branch October 24, 2026 09:23
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

bug Something isn't working documentation Improvements or additions to documentation

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant