Search every genome with the same proteins for the pangenome

.github/workflows/tests.yml

@@ -22,11 +22,18 @@ concurrency:
 jobs:
   test:
     name: Python ${{ matrix.python-version }}
-    runs-on: ubuntu-latest
+    runs-on: ${{ matrix.os }}
     strategy:
       fail-fast: false
       matrix:
-        python-version: ["3.9", "3.11"]
+        # 3.9 is the Python of the AMAP container. setup-python has no 3.9 for
+        # Ubuntu 26.04, which ubuntu-latest now is, and the job stopped at
+        # setup before running a test, so it runs on 24.04.
+        include:
+          - python-version: "3.9"
+            os: ubuntu-24.04
+          - python-version: "3.11"
+            os: ubuntu-latest
 
     steps:
       - uses: actions/checkout@v4

CHANGELOG.md

@@ -81,6 +81,7 @@ directly comparable with 1.0.1.
 - **The Mycobrowser table is checked whenever it is read**: validation checked `--mycobrowser_ref` only with the functional annotation, though the interactive report reads it in every run, so under `--skip_functional_annotation` a wrong path went unwarned and the report came out without Mycobrowser's categories
 - **The duplication networks draw offline**: both duplication steps loaded their networks' drawing library from a CDN, so without the internet the pages drew nothing. The library now goes inside each page, about 0.7 MB
 - **`environment.yml` pins eggNOG-mapper 2.1.13**, the version the container ships, where it said 2.1.12
+- **The pangenome shared no gene between samples**: PANGENE aligned each sample's own proteins to the other samples' genomes, ran pangene once per sample and named each column `<genome>_<sample>#0`, and REPORT_INPUTS merged those matrices on gene names that carry the sample's name and coordinates (`A4_Mbovis_Rv0001;+;0-1524`, `L4_CDC1551_Rv0001;+;0-1524`), so no row was ever in two samples and no gene was ever core. With N samples the matrix had N(N-1) columns: on the first 2 Mb of CDC1551, M. bovis and H37Rv the report showed 5379 gene families in 6 samples, all shell, with 200 genes of one genome as its heatmap, and as no genome was searched for its own genes, the interactive report found no gene in one genome alone. PANGENE now uses pangene as it is meant to be used: one set of proteins, the reference proteome and each sample's genes that no reference gene matches, aligned to every genome, the sample's own included, and one graph for the run, N alignments where there were N(N-1). A row is a reference gene (`Rv1506c`) or, outside the reference, one sample's copy of a gene (`L4_CDC1551_01396`), and a column a sample, so the same run shows 1810 genes, 1760 of them in all three, RD4 and RD10 missing from M. bovis only, the phiRv1 prophage from CDC1551 only, and three genes of the RvD2 region that CDC1551 and M. bovis carry and H37Rv lacks. Its files are now `PANGENE/pangenome.{gfa,bubbles.txt,Rtab,stat}`, one set for the run instead of one per sample, the graph statistics, which were never published, among them
 
 ### 📊 Report
 - **A report organised by question, `results/amap_report.html`**: next to `final_report.html`, a new process, `AMAP_REPORT`, writes one self-contained file from each genome's own outputs: a table of the genomes with the values far from comparable genomes flagged, a map of every reference gene in every genome where regions of difference show as blocks, a lookup of one gene in every genome, the correction from Bakta's calls to AMAP's genes with the genes behind each step, the pseudogenes by lineage, what lies outside the reference, and the functional profile, the categories that differ most between lineages first. It needs only what the correction leaves, so a run with `--skip_functional_annotation`, which gets no `final_report.html`, gets this one, and it names the genomes of the sample sheet that produced no annotation. Fifteen MTBC genomes make 0.5 MB (`final_report.html` is 4.4 MB for twelve); it draws everything with its own JavaScript and needs only Python's standard library to write
@@ -112,8 +113,8 @@ directly comparable with 1.0.1.
 - **`tests/test_eggnog_db.py`** runs `EGGNOG_DB`'s own script, rendered as Nextflow renders it, with its folder staged as Nextflow stages an input and stand-ins for eggNOG-mapper's download and build: the database is built in that folder whatever it is called, one already there is used as it is, and a short `bacteria.dmnd` is built again
 - **`tests/test_eggnog_annotation.py`** runs `ANNOT_EGGNOGMAP_BATCH`'s own script the same way, with a stand-in for `emapper.py` that looks for the database as eggNOG-mapper does: the database is found whatever its folder is called, a folder called `split`, where the script writes each sample's annotations, is left as it was, and the tag and log count a single sample as one
 - **`tests/test_reference_files.py`** runs the scripts of `BAKTA_ANNOTATOR`, `PATHOTYPR`, `GET_OUTPUTS` and `BUNDLED_REFERENCE`, rendered as Nextflow renders them, their own Groovy included, with stand-ins for Bakta and pathotypr and each reference file on the host, under a relative path or inside the container; and fails when a process reads a reference file's param other than as what it falls back on when nothing is staged; and checks the Mycobrowser table whenever a step reads it
-- **`tests/test_pangene_genomes.py`** runs PANGENE_GENOME and PANGENE's scripts as Nextflow renders them, with stand-ins for miniprot and pangene, and fails unless each sample is aligned onto every other sample's genome and never onto its own, typed or not
-- **`tests/test_report_tables.py`** merges two Rtabs as PANGENE writes them and checks that the report reads every table a `#` can reach without a comment character; where R is installed, as in the container, it also runs the report's pangenome code on the merged matrix
+- **`tests/test_pangene_genomes.py`** runs PANGENE_GENOME and PANGENE's scripts as Nextflow renders them on real genomes of CDC1551, M. bovis and H37Rv cut around RD4 and RvD2 (`tests/data/pangenome/`, cut again by its `make_fixture.py`), with stand-ins for miniprot, pangene, k8 and gfatools that hand names on as the real ones do, and fails unless every genome, the sample's own included, is searched with the same proteins, each column is named after its genome, typed or not, a genome none of them is found in keeps a column of zeros, and a protein two samples share goes in once. 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`** merges the matrix PANGENE writes as REPORT_INPUTS does and checks that the report reads every table a `#` can reach without a comment character; where R is installed, as in the container, it also runs the report's pangenome code on it and counts the genes all its samples share
 - **`tests/test_cog_counts.py`** counts COG categories as `count_annot.py` does, for genes in none, in two and in lowercase, then runs the steps after the correction as `tests/test_sample_ids.py` chains them, with eggNOG giving some genes no category, then no gene one, then some samples no gene in a category another sample has genes in, and reads what the report's own overview, per-COG figures and statistics table show where R is installed: each category in a figure of its own, with every sample, at 0 where it has no gene. The Mycobrowser figures, statistics and ratios to the median are read the same way for a sample with no gene in one of their categories
 - **`tests/test_duplication_report.py`** runs GET_DUPLICATIONS and GET_ALLDUP as Nextflow renders them, chained as the pipeline chains them, with eggNOG giving one copy of a duplicated gene a category, another none and the last no row at all, and checks the columns and COG counts of `homologous_communities.tsv`; where R is installed, it knits the report's duplication summary on that table and on the one `GET_ALLDUP` writes when nothing is duplicated. It also fails when a chunk of the report prints markdown without passing it on as markdown
 - **`tests/test_profile_caps.py`** fails when a profile sets a cap, or when the order of `-profile` changes one: it reads `nextflow.config` and applies its profiles in every order, as both of Nextflow's parsers do, and with Nextflow on PATH also asks `nextflow config`
@@ -127,7 +128,8 @@ directly comparable with 1.0.1.
 - **`tests/test_citation.py`** fails when anything cites AMAP with a DOI other than the concept DOI of the manifest
 - **`tests/test_clean_proteome.py`** runs the steps after the correction as `tests/test_sample_ids.py` chains them, with FIX_CDS calling some genes pseudogenes, one with a stop codon in its protein and the others without, and Mycobrowser marking a gene AMAP finds intact, with Mycobrowser on and off, and checks that `*-AMAP.clean.faa` holds every protein FIX_CDS wrote but the pseudogenes', as it wrote them: the CDS the GFF3 does not mark `pseudo=true`, with the GenBank's translations
 - **`tests/test_report_pseudogenes.py`** runs the steps after the correction as `tests/test_sample_ids.py` chains them, with FIX_CDS calling a different number of pseudogenes in each sample, none in one, and Mycobrowser marking a gene that every sample has intact, with Mycobrowser on and off, and checks that REPORT_INPUTS hands over each master table row with its sample; where R is installed, it runs the report's own pseudogene section on what REPORT_INPUTS made and checks that it draws each sample at the number AMAP called, the QC table's, a sample with none at 0, under a title that does not name Mycobrowser
-- **744 tests**, up from 92
+- **CI's Python 3.9 job runs on Ubuntu 24.04**: 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. 3.9 is the container's Python, so it now runs on 24.04; 3.11 stays on the latest
+- **753 tests**, up from 92
 
 ## [1.0.1] - 2026-07-25
 

bin/pangenome.py

@@ -0,0 +1,171 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+The proteins PANGENE searches every genome with, and the matrix it writes.
+
+pangene builds its graph from one set of proteins aligned to every genome, the
+way its documentation describes: a gene is then one row and each genome one
+column, its own included. AMAP's set is the reference proteome, each protein
+named after its gene (Rv0001 for Rv0001_1-1524), and the genes of each sample
+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.
+
+    pangenome.py proteins --reference REF.faa --out proteins.faa SAMPLE.AMAP.tsv [...]
+
+reads each sample's proteins from the <sample>.AMAP.faa next to its table.
+
+    pangenome.py matrix --out pangenome.Rtab GENES_BY_GENOME.Rtab GENOME.fasta [...]
+
+turns what `pangene.js gfa2matrix` writes into the matrix AMAP keeps: one
+column per genome, named after its file (<name>.fasta), where pangene names it
+<name>#0, and a column of zeros for a genome none of the proteins was found
+in, which pangene leaves out.
+"""
+
+import argparse
+import csv
+import re
+import sys
+from pathlib import Path
+
+# Rv0001_1-1524: a reference gene and where the reference has it
+COORDS = re.compile(r'_\d+-\d+$')
+# pangene names a genome after its PAF and the haplotype: <name>#0
+HAPLOTYPE = re.compile(r'#\d+$')
+
+
+def read_fasta(path):
+    """(first word of the header, sequence) for each record. The H37Rv
+    proteome in the container has Windows line ends."""
+    name, seq = None, []
+    with open(path) as fh:
+        for line in fh:
+            line = line.strip()
+            if line.startswith('>'):
+                if name is not None:
+                    yield name, ''.join(seq)
+                words = line[1:].split()
+                name, seq = (words[0] if words else ''), []
+            elif line:
+                seq.append(line)
+    if name is not None:
+        yield name, ''.join(seq)
+
+
+def gene_name(name):
+    """pangene reads a protein's name as <gene>:<protein> and takes the gene
+    up to the first colon, so a colon of the gene's own would cut it short."""
+    return name.replace(':', '_')
+
+
+def genes_outside_reference(table):
+    """IDs of the genes of a <sample>.AMAP.tsv that no reference gene matches."""
+    with open(table, newline='') as fh:
+        rows = csv.DictReader(fh, delimiter='\t')
+        if 'id' not in (rows.fieldnames or []) or 'rvs' not in rows.fieldnames:
+            sys.exit(f"{table}: no id and rvs columns, as FIX_CDS writes them")
+        return {row['id'] for row in rows if not (row['rvs'] or '').strip()}
+
+
+def protein_file(table):
+    """<sample>.AMAP.faa, next to <sample>.AMAP.tsv."""
+    table = Path(table)
+    return table.with_name(table.name[:-len('.tsv')] + '.faa')
+
+
+def write_proteins(reference, tables, out):
+    """Write the set; return (reference proteins, sample proteins, repeats left out)."""
+    seen, genes = set(), set()
+    n_ref = n_own = n_repeat = 0
+    with open(out, 'w') as fh:
+        for name, seq in read_fasta(reference):
+            seq = seq.rstrip('*')
+            if not seq:
+                continue
+            gene = gene_name(COORDS.sub('', name) or name)
+            fh.write(f">{gene}:{name}\n{seq}\n")
+            seen.add(seq)
+            genes.add(gene)
+            n_ref += 1
+        for table in sorted(tables):
+            outside = genes_outside_reference(table)
+            if not outside:
+                continue
+            for name, seq in read_fasta(protein_file(table)):
+                gene_id = name.split(';')[0]
+                seq = seq.rstrip('*')
+                if gene_id not in outside or not seq:
+                    continue
+                if seq in seen:
+                    n_repeat += 1
+                    continue
+                gene = gene_name(gene_id)
+                if gene in genes:
+                    # Never two genes under one name: pangene would take
+                    # them for forms of the same gene.
+                    k = 2
+                    while f"{gene}_{k}" in genes:
+                        k += 1
+                    gene = f"{gene}_{k}"
+                fh.write(f">{gene}:{name}\n{seq}\n")
+                seen.add(seq)
+                genes.add(gene)
+                n_own += 1
+    return n_ref, n_own, n_repeat
+
+
+def genome_name(path):
+    """A genome's name, as PANGENE names its alignments: the file's name
+    without .fasta."""
+    name = Path(path).name
+    return name[:-len('.fasta')] if name.endswith('.fasta') else name
+
+
+def write_matrix(raw, genomes, out):
+    """Write the matrix; return (genes, genomes, genes in every genome)."""
+    with open(raw) as fh:
+        header = fh.readline().rstrip('\r\n').split('\t')
+        rows = [line.rstrip('\r\n').split('\t') for line in fh if line.strip()]
+    # gfa2matrix with no genome at all writes "Gene" and a tab
+    cols = [HAPLOTYPE.sub('', c) for c in header[1:]]
+    names = sorted({genome_name(g) for g in genomes} | {c for c in cols if c})
+    where = {name: i for i, name in enumerate(names)}
+    n_core = 0
+    with open(out, 'w') as fh:
+        fh.write('\t'.join(['Gene'] + names) + '\n')
+        for row in rows:
+            values = [0] * len(names)
+            for col, value in zip(cols, row[1:]):
+                if col:
+                    i = where[col]
+                    values[i] = max(values[i], 1 if value.strip() not in ('', '0') else 0)
+            n_core += all(values)
+            fh.write('\t'.join([row[0]] + [str(v) for v in values]) + '\n')
+    return len(rows), len(names), n_core
+
+
+def main(argv=None):
+    parser = argparse.ArgumentParser(description=__doc__.split('\n\n')[0].strip())
+    sub = parser.add_subparsers(dest='command', required=True)
+    p = sub.add_parser('proteins', help="Write the proteins every genome is searched with.")
+    p.add_argument('--reference', required=True, help="The reference proteome.")
+    p.add_argument('--out', required=True, help="The protein set to write.")
+    p.add_argument('tables', nargs='*', help="Each sample's <sample>.AMAP.tsv.")
+    m = sub.add_parser('matrix', help="Name the matrix's columns after the genomes.")
+    m.add_argument('--out', required=True, help="The matrix to write.")
+    m.add_argument('raw', help="What `pangene.js gfa2matrix` wrote.")
+    m.add_argument('genomes', nargs='*', help="Every genome aligned to, <name>.fasta.")
+    args = parser.parse_args(argv)
+
+    if args.command == 'proteins':
+        n_ref, n_own, n_repeat = write_proteins(args.reference, args.tables, args.out)
+        print(f"{n_ref} reference proteins and {n_own} of genes outside the reference "
+              f"({n_repeat} more already in the set)")
+    else:
+        n_genes, n_genomes, n_core = write_matrix(args.raw, args.genomes, args.out)
+        print(f"{n_genes} genes x {n_genomes} genomes, {n_core} in every genome")
+    return 0
+
+
+if __name__ == '__main__':
+    sys.exit(main())

bin/report/amap_report.js

@@ -2060,7 +2060,7 @@
           }
           return { smp, total, all, none, some: total - all - none, others: p.others };
         });
-        el4 = card('Pangenome', 'What PANGENE found: each genome\'s genes aligned to every other genome. A gene in all of them is core to this run; one in none is this genome\'s own.', null,
+        el4 = card('Pangenome', 'What PANGENE found: the reference proteome and the genes no reference gene matches, searched for in every genome. A gene in every other genome is core to this run; one in none is this genome\'s own.', null,
           table({
             columns: [
               { key: 'genome', label: 'Genome', value: (r) => r.smp.id },

bin/report/amap_report.py

@@ -1068,28 +1068,26 @@ def build(args):
 
 
 def pangene_summary(rtab, samples):
-    """What PANGENE found, per genome: how many of its genes each other
-    genome carries. The merged matrix has one row per gene of every genome
-    and one column per genome a genome's genes were aligned to, named
-    <that genome>_<this genome>#0, so each genome's rows are read against
-    its own columns only."""
+    """What PANGENE found, per genome: how many of the genes it carries each
+    other genome carries too. The matrix has one row per gene of the run and
+    one column per genome, named after it, so a genome's genes are the rows
+    with a 1 in its column, and the other columns say where else they are."""
     if not rtab:
         return None
     names = [s['name'] for s in samples if not s.get('failed')]
+    columns = rtab['columns']
     out = {}
     for name in names:
-        cols = [i for i, c in enumerate(rtab['columns']) if c.split('#')[0].endswith('_' + name)]
-        if not cols:
+        if name not in columns:
             continue
+        own = columns.index(name)
+        others = [i for i in range(len(columns)) if i != own]
         hist = Counter()
-        n_genes = 0
-        for gene, values in rtab['rows']:
-            if not gene.startswith(name + '_'):
-                continue
-            n_genes += 1
-            hist[sum(values[i] for i in cols if i < len(values))] += 1
-        if n_genes:
-            out[name] = {'others': len(cols), 'hist': {str(k): v for k, v in sorted(hist.items())}}
+        for _gene, values in rtab['rows']:
+            if own < len(values) and values[own]:
+                hist[sum(values[i] for i in others if i < len(values))] += 1
+        if hist:
+            out[name] = {'others': len(others), 'hist': {str(k): v for k, v in sorted(hist.items())}}
     return out or None