diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 84109a9..1e4df1d 100644 --- a/.github/workflows/tests.yml +++ b/.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 diff --git a/CHANGELOG.md b/CHANGELOG.md index b16577a..3fe370e 100644 --- a/CHANGELOG.md +++ b/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 `_#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 diff --git a/bin/pangenome.py b/bin/pangenome.py new file mode 100644 index 0000000..b062d68 --- /dev/null +++ b/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 .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 .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 (.fasta), where pangene names it +#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: #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 : 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 .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): + """.AMAP.faa, next to .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 .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, .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()) diff --git a/bin/report/amap_report.js b/bin/report/amap_report.js index a1506a3..9ef4ce0 100644 --- a/bin/report/amap_report.js +++ b/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 }, diff --git a/bin/report/amap_report.py b/bin/report/amap_report.py index a8041fc..e2e12b5 100644 --- a/bin/report/amap_report.py +++ b/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 - _#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 diff --git a/docs/about/changelog.html b/docs/about/changelog.html index 9790964..2337032 100644 --- a/docs/about/changelog.html +++ b/docs/about/changelog.html @@ -96,7 +96,7 @@

[Unreleased]🔬 Correctness#

  • Paralog copies no longer take another gene's place: a protein with a paralog was placed at both loci, and which copy survived came down to insertion order. On H37Rv itself that lost Rv1150, Rv1945 and Rv2828c; on M. bovis it reported genes the strain lacks (the mce3 operon of RD7, the phiRv2 prophage of RD11, plcBC of RD5) as copies sitting on their paralogs. Each placement is now checked against where its protein belongs, predicted from the sample's own offset to the reference so that indels do not throw it off
  • Strand of a resolved conflict: taken from the protein's first placement anywhere in the genome instead of the one being resolved, so esxM (Rv1792), Rv3466 and two IS6110 copies came out at their own coordinates on the wrong strand; a Bakta call chosen over an alignment also took the alignment's strand
  • Reference pseudogenes keep their extent: PE21 (Rv2099c) and PPE48 (Rv3022c), which the reference annotates without a stop, were extended to the next one, and oxyR' (Rv2427A) and Rv3216 lost the partial final codon the reference gives them
  • The 1.0.1 fix for reference pseudogenes never fired: it looked for reference coordinates in a column that has not carried them since March 2026, and its tests wrote them in by hand. check_cds.py now reads the reference proteome instead
  • Alignments that stop short are completed: miniprot can miss a protein's N-terminus (Rv3611 on CP110619) or leave a repetitive C-terminus unaligned further from its stop than the 300 bp searched; the reference protein now says how much is missing
  • Ties between copies of a gene go to the one Bakta also called, then to the one where the protein belongs by synteny, not the one that comes first in the genome (E. coli K-12 b2635)
  • Draft assemblies: genes were compared on coordinates alone, and every contig starts at 1, so cut into 40 contigs CDC1551 had 87% of its proteins translated from the wrong contig. Every gene is now matched, translated and written on its own contig, and Bakta runs with --keep-contig-headers so that its contig names match the FASTA's
  • Case 4 is grouped strand by strand: about 8% of its groups (101 of 1255 on H37Rv) put an alignment on one strand together with a Bakta call on the other
  • A complete alignment is not displaced by a longer call in its frame: such a call ends at the same stop and only starts further upstream (plcD on CDC1551). Nor does a call that spans several reference genes win against them (Rv3022c and Rv3021c, one reading frame in CDC1551)
  • Reference genes are protected again: the check looked for coordinates in identifiers that no longer carry them, so any call of similar length could displace a reference gene (Rv1614)
  • Missing genes are recovered by what covers them: genes stored under Bakta locus tags looked missing by name, so a placement a codon off its gene went back in as a second copy. With the case 4 changes, ORFs annotated twice fall from 10 to 3 on CDC1551 and from 16 to 6 on M. bovis, where the reference has none
  • Real copies of a gene are kept: deduplicate_rvs.py deleted every copy of a gene but one, 34 genes on E. coli O157:H7 Sakai, which were then missing from the annotation and from the duplication analysis. Further copies get identifiers of their own
  • Genes that lost their stop no longer swallow the next gene: check_cds.py extended a gene without a stop to the next one in its frame, which after an indel is the next gene's, so that gene was annotated a second time in the same frame (M. bovis Rv0590 over Rv0590A, Rv3373 over Rv3374, Rv2293c over Rv2292c). Such a gene now keeps its extent as a pseudogene, unless the extension only restores residues its reference protein has. ORFs annotated twice fall from 3 to 2 on CDC1551 and from 6 to 3 on M. bovis, and none of those left comes from an extension
  • Small proteins: miniprot placed none of the 33 proteins of 8 to 29 residues in E. coli K-12, and put the toxin IbsC (b4665) only on its near-identical neighbour IbsD. Every protein without an exact placement is now looked for in the six-frame translation of the genome, which takes K-12 against itself from 4252 to all 4286 of its genes, its alternative-start isoforms counted with their gene, and Sakai from 198 to 220 of 220
  • Frameshifts and lost stops: a gene was extended to the next stop in its frame however far that was, and came out as a complete gene with a C-terminus its reference lacks; 18 in CDC1551 and 20 in M. bovis. It now keeps its extent and is flagged as a pseudogene when that stop lies more than ten residues past the end of the reference protein, as PGAP and Bakta treat a frameshift, which also covers a frame that runs into the next gene
  • Genes a contig end cuts: in a draft assembly they came out as pseudogenes, or as complete genes, without a word. A gene whose frame runs off the contig, or that opens at its edge without a start codon, is now marked partial at that end. In CDC1551 cut into 40 contigs, the larger piece of each of the 25 genes a cut runs through is marked so
  • The other pieces of a cut gene: miniprot reports only the placements within 5% of a protein's best, so only the larger piece of a gene a contig end cuts was annotated. What that placement lacks is now aligned again as a protein of its own, and its placement at the end of another contig is annotated under the gene's name and marked partial: 47 of the 50 pieces in the CDC1551 draft instead of 25, the other three being five or six codons long. A gene that opens at a contig's edge with a codon that can also start one, as one in six codons inside H37Rv genes can, is marked partial too when its protein begins inside its reference protein
  • Disrupted copies of a gene: deduplication kept further copies only when intact, so a whole copy that is a pseudogene was deleted with the pieces of split genes; it now stays when nearly as long as the reference protein
  • GFF3 references: each CDS line was translated on its own, so a CDS of several lines came out in fragments, and proteins under 11 residues were dropped to hide them. The lines are joined per CDS and translated with the bacterial code, reading transl_except, and E. coli K-12's GFF3 gives the same 4318 proteins as its GenBank record, 24 of which were missing
  • Genes longer than their reference keep its name: a Bakta call that won its locus over a shorter alignment in its frame kept no name, and the protein was missing from the table: Rv2024c in CDC1551 and M. bovis, mmpL6 (Rv1557, whole in M. bovis where TbD1 cut it in H37Rv), Rv1068c, Rv1917c and Rv3135 in M. bovis, fliD (b1923) and b0282 in Sakai. The call now takes the name when the alignment spans at least half of the protein, so a remnant does not name the gene it ends (the 67 residues of cyp141 that RD12 leaves in M. bovis)
  • Small proteins with a few substitutions: only exact copies were looked for, so Sakai's ten K-12 micro-proteins with one to four residues changed were not annotated. A protein of up to 50 residues with no exact copy is now looked for with up to a fifth of its residues changed; a stop codon in place of a residue counts, and the gene is then a pseudogene (b1229)
  • Synteny no longer hangs on the order of the rows: the anchors were sorted by start and searched by midpoint, so two anchors starting together sent a placement to whichever window the row order gave, and adding unrelated rows moved Sakai's b2127 24 kb and handed its locus to the paralog b1162. A locus several proteins are placed at exactly no longer anchors either: the single IS6110 of M. bovis made 34 anchors with reference positions all over the genome, and came out under all 17 IS6110 names of each of its genes; they are now Rv2814c and Rv2815c, its own
  • Genes without a COG category are counted once: count_annot.py filled in a missing category as Unassigned and split it into letters, as it splits JK into J and K, so every gene eggNOG gave no category was counted once in A, D, E, G, I and U and twice in N and S, in *-AMAP.COG.tsv and in the report's COG figures and statistics, where it also shrank every other category's share. A test-profile run with a stand-in for eggNOG that gave none of its 4036 genes a category reported 4036 genes in each of those six and 8072 in N and S. Such a gene now counts once, as NaN, which the report shows as None, and its percentage is of all the genome's genes like every other category's
  • COG categories R, W, X, Y and Z in the report: its COG figures took their categories from a list that ran from A to V without R, so a gene eggNOG put in R, W, X, Y or Z was plotted as NA in the overview and in none of the per-category figures, which stayed empty for those five, X - Mobilome (prophages and transposons) among them. The categories now come from the report's own list of figures, A to Z and None, so each is plotted in its own. The statistics table never used that list and is unchanged
  • The QC table counts the pseudogenes AMAP called: annotation_qc.py read the first of Is_Pseudogene, pseudogene and is_pseudogene that the master table has, and every master table has Mycobrowser's Is_Pseudogene, empty when Mycobrowser is off. So pseudogene_count in *-AMAP.qc.tsv and the report's QC table was 0 in every organism-agnostic run (E. coli Sakai against K-12, where AMAP called 129 pseudogenes), and in an MTBC run the number of the sample's genes named after one of the 13 that Mycobrowser marks as pseudogenes in H37Rv, whatever AMAP made of them: 13 for CDC1551, where AMAP calls 123 and finds 6 of those 13 intact. It now counts AMAP's own call, the master table's pseudogene column, the same genes the GFF3 marks pseudo=true and the GenBank /pseudo
  • The report's pseudogene chart counts AMAP's calls, sample by sample: its Pseudogene Detection section counted the master table's genes flagged in Mycobrowser's Is_Pseudogene or in AMAP's pseudogene column, lineage by lineage, as "Pseudogenes per Lineage (Mycobrowser annotation)". Is_Pseudogene says which genes are pseudogenes in H37Rv, not in the sample, so CDC1551 came out at 129 where AMAP calls 123 and the QC table says 123, the other 6 being genes Mycobrowser marks that CDC1551 has intact (Rv0947c, Rv1150, Rv2023A, Rv2098c, Rv3021c, Rv3324A), and without Mycobrowser the chart showed AMAP's calls under Mycobrowser's name. It also left out a lineage with none, and added up a lineage's samples, so a lineage with more samples looked as if it had more pseudogenes. REPORT_INPUTS now hands the report every master table row with the sample it is from, and the chart shows how many pseudogenes AMAP called in each sample, by lineage, a sample with none at 0, next to a table of them. The master table explorer shows the note that says why a gene is a pseudogene
  • A sample with no gene in a category counts 0 in the report: count_annot.py writes a row for each COG or Mycobrowser category a sample has genes in, and the report plotted the rows as they came, so a sample with no gene in a category was left out of that category's figures, of the Dunn tests between lineages run on them, of the overview and of the statistics table, instead of counting 0. Its lineage's box, mean and median were those of the samples that have some, a lineage none of whose samples has any dropped out of the comparison, and the comparison with the median never showed a sample that had lost a category. Each sample now counts 0 in every category another sample of the run has genes in, and the statistics table's Count is every sample of the lineage. A category no sample has keeps its empty figure, and *-AMAP.COG.tsv and *-AMAP.MycoBrowser.tsv are unchanged
  • The clean proteome leaves out the pseudogenes AMAP called: export_clean_proteome.py took them from the first of Is_Pseudogene, pseudogene and is_pseudogene that the master table has, which is Mycobrowser's, and looked each protein up by its whole FASTA ID, L4_CDC1551_Rv0001;+;0-1524 where the master table has L4_CDC1551_Rv0001, so it never found one and its log said 0 pseudogenes. Only the proteins with an internal stop codon were left out, and the pseudogenes without one stayed in *-AMAP.clean.faa, which the docs say has none: 42 of the 123 AMAP calls in CDC1551 (a frameshift or lost stop codon, no stop codon within reach, no stop codon as in the reference) and 27 of the 129 in E. coli Sakai against K-12. It now leaves out the genes whose pseudogene is yes, found by the ID before the ;, and on both genomes holds the CDS the GFF3 does not mark pseudo=true, 3891 and 5194, with the proteins the GenBank translates for them. A master table without that column stops it, where a warning let every pseudogene through

🐛 Pipeline#

-
  • Pseudogenes in the GFF3 and GenBank: the GFF3 never said which genes were pseudogenes, and the GenBank only for genes with an internal stop, though the master table knew. Both now carry AMAP's flag with its reason (pseudo=true and Note, /pseudo and /note), and mark the ends a contig cut off (partial=true with start_range or end_range, < or >); the master table and *-AMAP.tsv have them as their last columns
  • The default H37Rv reference under Docker: its files live inside the container and were handed to processes as paths, so Nextflow mounted /reference from the host, which Docker Desktop refuses and Docker on Linux mounts empty over the container's own. BUNDLED_REFERENCE copies them out first
  • -profile test: VALIDATE_INPUTS still required the sample sheet the test profile does without. 1.0.1 had only been checked with -preview, which runs no task
  • -resume kept results made by older code: Nextflow hashes a task's inputs and script, not the files it runs by path, so after an edit to the classifier or any script under bin/, -resume reused the old results. Every process that runs AMAP code now lists those files in ext code, which puts them in its cache key; edit bin/classify/cases.py and TYPE_CASES and FIX_CDS run again while PATHOTYPR and Bakta stay cached
  • PATHOTYPR no longer depends on the classifier: its script said pathotypr classify, and Nextflow adds any bin/ entry a script names as a bare word to the cache key, so a change in bin/classify/ re-ran PATHOTYPR, Bakta and everything after them. The subcommand is quoted
  • Nextflow 25.10 and later: the pipeline stopped before running anything, because the strict parser those versions use by default rejects the helper functions in nextflow.config, the statements outside the workflow in main.nf, the bare env(SAMPLE_NAME) output and publishDir paths that name a task input. Nextflow 26.10 also refuses the return that ended a --validate_only run. The resource caps are now the resourceLimits directive, which still takes --max_cpus, --max_memory and --max_time from the command line, a profile or a -c config. Nextflow 24.04, where that directive first appeared, is the new minimum
  • --skip_mycobrowser in organism-agnostic mode: a custom reference from another genus turned it on by reassigning the parameter, which current Nextflow ignores with a warning, so the run went ahead with Mycobrowser enabled. It is now worked out in the parameter's default
  • Tasks leave the pipeline untouched: the first time a copy of AMAP ran, Python cached the bytecode of AMAP's own modules inside it, in bin/classify/__pycache__, and with pathotypr classify naming that directory in PATHOTYPR's script, the first -resume re-ran PATHOTYPR, Bakta and everything after them. Tasks now run with PYTHONDONTWRITEBYTECODE=1
  • Samples with no Bakta call outside cases 2 and 3: when every Bakta call had an alignment at exactly its coordinates, which a small input can do, TYPE_CASES stopped with a KeyError, because the empty column of overlaps was read as a list of column names. Case 1 and case 4 now come out empty and the rest is classified as usual. An empty miniprot table stopped it the same way
  • Samples Bakta calls no CDS on: with no contig names from Bakta to compare, the contig check took them for different names and stopped the run, and a Bakta table with nothing but its header did not read. The genes miniprot places now come out as missing from Bakta, like any gene it did not call
  • --version: --help and the parameter docs have listed it since 1.0.0, but nothing read it, so nextflow run PathoGenOmics-Lab/AMAP --version started a run. It now prints the pipeline and Nextflow versions and exits before any task
  • Drafts named after one of their contigs: pathotypr writes a summary row for each contig that carries lineage markers, and PATHOTYPR took the sample's name from the first, so CDC1551 cut into 40 contigs, sample CDC1551d, came out as L4_contig_14, and so did every file after it. Two drafts whose first row named the same contig became one sample, the results of one overwriting the other's. A genome is now named <lineage>_<sampleId> from the sample sheet, the lineage of a draft being the one the markers of all its contigs add up to, and the contigs of a draft keep their names instead of the first one taking the sample's. A closed genome whose header is its sample ID comes out as before
  • Headers with a description: the name was the whole header line, so the tutorial's CDC1551 (NC_002755.2 Mycobacterium tuberculosis CDC1551, complete sequence) was named with spaces in it, Bakta failed on it three times and the run finished without it, and without failing. Such a genome now takes its sample ID like any other, and so does the test profile's: L4_H37Rv_test, where it was L4_CP110619_H37Rv
  • An input could be emptied: PATHOTYPR wrote the renamed genome next to the link Nextflow stages its input as, so when the new name was the file's own, sample G26689 in L4_G26689.fasta for instance, the input was overwritten with nothing and the run went on without it. The input is now staged in a directory of its own
  • Sample IDs with underscores: the steps after the correction told a sample's lineage from its sample ID by splitting its name at underscores, though typing only puts the lineage in front and sample IDs may hold underscores of their own. L4_H37Rv_test was counted as genome H37Rv, the tutorial's Beijing_1 and a Beijing_2 as one genome Beijing in the category counts, the report and the duplication network, an untyped H37Rv_test as lineage H37Rv, and each gene of an untyped sample with no underscore as a genome of its own in the network. PATHOTYPR now says which lineage it named a genome after and every step is handed the sample ID and the lineage. The case summary gains a Sample column and names an untyped sample as it is, not Untyped_<sampleId>, and the report no longer reads the first gene not in the reference as a header, nor stops when every gene has one
  • Genomes kept outside the run's folder: VALIDATE_INPUTS staged only the sample sheet and opened the paths written in it from inside its container, which sees only the folder that the work directory and the sheet share. Genomes anywhere else, /data/genomes while running from ~/project for instance, were "FASTA file not found" and the run stopped before anything ran, and a relative path, which the sheet may use, was read from the task directory with or without a container. The genomes are now inputs of the process, each in a directory of its own so that two with the same file name do not collide, and each is read where Nextflow staged it
  • Runs with functional annotation stopped at validation: VALIDATE_INPUTS looked for --eggnog_db in its own task directory, where the default eggnog_db, relative to the launch directory, never is, and took a database not built yet for an error, though EGGNOG_DB builds it. It now looks where EGGNOG_DB builds it, staged when it exists, and stops the run only when eggnog_db.ok is there without a complete database, since Nextflow then never runs EGGNOG_DB again
  • A Bakta database on the host: --bakta_db reached the tasks as a path the container does not mount. With -profile docker, validation called the database missing and stopped the run, and where nothing validated it, as under -profile test, every sample failed three times and was ignored, and the run ended with nothing annotated and exit 0. A database on the host is now staged into BAKTA_ANNOTATOR, PREPARE_REFERENCE and VALIDATE_INPUTS like any other input; the default inside the container is used as before
  • A warning on every run: validation expected a db.sqlite3 that no Bakta database has. It now checks for version.json and bakta.db, which Bakta cannot run without, so their absence stops the run, and names the database's type and version
  • --strict_validation: the documentation said validation warnings could be made to stop the run, but nothing passed the script's --strict on. The new parameter does
  • eggNOG's database under Docker: EGGNOG_DB built it in a folder it worked out for itself from --eggnog_db, and having no inputs, its container saw that folder only when Nextflow happened to mount it with the work directory and the pipeline's bin/, through a parent folder they share. Otherwise, as when running from /data/project with the pipeline in ~/.nextflow, the database was downloaded and built inside the container and lost with it, and only eggnog_db.ok reached the host: later runs skipped EGGNOG_DB, and ANNOT_EGGNOGMAP_BATCH found no bacteria.dmnd. A database already there without eggnog_db.ok was built again. main.nf now makes the folder on the host and hands it to EGGNOG_DB as an input. A folder left with only eggnog_db.ok needs that file deleted, as validation says
  • eggNOG's database in a folder of any name: ANNOT_EGGNOGMAP_BATCH asked eggNOG-mapper for eggnog_db/bacteria.dmnd and --data_dir eggnog_db, but Nextflow stages a folder under its own name, so with a --eggnog_db called anything else, as the HPC guide's /shared/databases/eggnog or the tutorial's ~/databases/eggnog, EGGNOG_DB built the database and the annotation stopped at "DIAMOND database eggnog_db/bacteria.dmnd not present", with or without a container. The folder is now staged under that one name, whatever it is called. The task's tag and first log line also gave a single sample's FAA size in bytes as the number of samples (batch_1519647_samples under -profile test)
  • Reference files on the host: BAKTA_ANNOTATOR gave Bakta --ref_annotation, or --bakta_proteins, as a path written into its script, which a container does not see, nor a task directory when it is relative, so the organism-agnostic mode failed under -profile docker: with K-12 as the reference for Sakai, Bakta called the "user proteins file" not valid three times, the sample was ignored and the run ended with exit 0. Likewise --pathotyper_snps and --pathotyper_ref on the host left every genome untyped, as NA_<sampleId>, --mycobrowser_ref left the Mycobrowser categories out, and one --miniprot_ref_* on the host while the others kept their defaults made BUNDLED_REFERENCE fail, so nothing was corrected. A file on the host is now staged into the tasks that read it, each in a directory of its own, and a default is read inside the container as before
  • Reference files that are nowhere: a typo in one of those paths, or a default on a run without the container, went unnoticed until the samples were ignored or left untyped. Validation now checks each reference file the run reads and that is not on the host, inside the container where the tasks read it, and stops the run when it is not there. A missing Mycobrowser table only warns, since GET_OUTPUTS goes on without it
  • The pangenome without lineage typing: PANGENE aligned each sample's proteins onto the genomes PATHOTYPR renamed, so with --skip_typing, which every custom reference from another genus needs, it had none and wrote empty files for every sample. It now aligns them onto the sample sheet's genomes, each copied to <sampleId>.fasta first, so that two files of the same name (assembly.fasta from two assemblies) stay apart and a .fna is not passed over
  • Genes not in the reference, sample by sample: the report counted only the samples that had some, so a sample whose genes all matched the reference dropped out of its lineage's box instead of counting 0, and when no sample had any (H37Rv against itself, as in -profile test) the section was an empty plot. Every sample is now counted, the section says so when none has any, and each point is drawn at its count, on an axis from 0, where it was jittered up and down by as much as 0.4
  • The report's pangenome section lost a genome: pangene names each genome <name>#0, and the report read the merged matrix with R's default comment character, so its header stopped at the first #. One genome column was dropped and the gene names became row names, so the core/accessory and frequency plots used a single column and the heatmap never appeared. The missing-gene and duplication tables, which carry the reference annotation's gene IDs as given, are now read without a comment character too
  • Genes across the origin of a circular sequence: Bakta reads a sequence as a complete replicon when its FASTA description says complete (NCBI's "complete genome" and "complete sequence") or circular=true (Unicycler), joins its ends and calls a gene across its origin from its start to its stop, start after end. TYPE_CASES stopped on such a gene: E. coli O157:H7 Sakai with NCBI's headers, whose stcE runs across the origin of plasmid pO157 (92563..2502), failed three times and was left out, and the run finished without failing. Past it, check_cds.py dropped such a gene as invalid coordinates, deduplication and the QC table gave it a negative length, the GFF3 wrote it backwards, the GenBank writer stopped on it and the Panaroo GFF3 spanned it over the whole sequence. It now goes through whole and is written as Bakta writes it: start after end in the tables, join() across the origin in a circular GenBank record, and in the GFF3 its end plus the sequence's length on a region marked Is_circular=true. A gene near the origin of such a sequence is read across it too, so Rop on plasmid pOSAK1, of which Bakta called the first 141 bp, comes out whole. Annotated against K-12, Sakai with NCBI's headers now has 5331 genes in common with the same genome under bare names, which Bakta reads as linear, and differs in stcE and Rop, whole where the linear reading has 804 of stcE's 886 residues and no Rop, and in two phage genes Bakta starts elsewhere on a complete genome
  • Typed closed genomes were annotated as linear: PATHOTYPR renamed a genome of one sequence to <lineage>_<sampleId> and nothing else, so the description Bakta reads a complete, circular replicon from never reached it, whatever the header said, and a gene across the origin of a closed genome came out as two pieces under one name, marked 5' and 3' partial. Only the sequence's ID is replaced now, and the description after it stays, as with --skip_typing. Typed from NCBI's headers, CDC1551, CP110619, AL123456.3 and M. bovis get the same genes as before, in place, name, flags and protein, though Bakta starts a few of its own calls elsewhere on a complete genome (dnaA at 1 instead of 34), so the gene IDs that are Bakta locus tags can be numbered differently. H37Rv turned so that the origin cuts gyrB, gyrA or pknB gives back all 4031 GenBank CDS, where it gave 4030 with the cut gene in two pieces. The GenBank record of such a genome is marked circular, with the description in its definition, and its GFF3 has a region marked Is_circular=true
  • The report's duplication summary: it printed the number of homologous community pairs as markdown from a chunk that did not pass its output on as markdown, so the report showed ## **48** homologous community pairs detected across all genomes. in a code box. It is now a sentence. The example report's duplication table had columns GET_ALLDUP never writes, so its Network Metrics tab said there were none; the demo data now has the columns GET_ALLDUP writes and the example is rendered again
  • Modularity in homologous_communities.tsv: Modularity1 and Modularity2 were computed with the whole community as its only group, which makes them 0 whatever the community, as on all 29 rows of CDC1551 against M. bovis. The two columns are gone
  • COG categories in homologous_communities.tsv: with numpy 2, as in the container, Elements1 and Elements2 read {'L': np.int64(4)}, and they counted each protein eggNOG gave no category in category + or -, since its ID keeps the strand FIX_CDS wrote. They are now plain counts, the largest first and ties in alphabetical order rather than in an order that changed from run to run, and a protein without a category counts as -, eggNOG's mark for none
  • The order of -profile: from Nextflow 25.10, whose strict parser applies the profiles in the order -profile names them and lets the last one win, the laptop caps of the container profiles took over from the others: -profile slurm,singularity held every task of a cluster run to the launching machine's CPUs and 16 GB, where -profile singularity,slurm gave it 48 CPUs and 128 GB, -profile test,docker ran the test with 16 GB and -profile highmem,slurm with 128 GB. No profile sets a cap now: they are set after the profiles, from whether highmem, test or a scheduler profile is on, so every order gives what Nextflow 24.04 gave. The parameter docs gave 24 CPUs and 50 GB as the defaults, which no run got
  • The duplication files are the same from one run to the next: homologous_communities.tsv listed its pairs of communities in the order of a set, which changes with Python's hash seed from one run to the next, and with it which community of a pair came first and the last digit of AvgClustering, an average over the community's proteins that networkx went through in the order of their set: 0.5833333333333333 or 0.5833333333333334 for the same pair of CDC1551 and M. bovis. The pairs now come in the order of the communities, the sample first by name on the left, and the metrics go through the proteins in alphabetical order. Each family in duplications_*_groups.tsv lists its proteins in alphabetical order, where they too came in the order of a set, and each sample's two network drawings no longer start from a random layout and a random partition, so every file both steps write is the same run after run
  • No web browser: both duplication steps wrote their networks with pyvis' show(), which also opens the file in a web browser, so a run without a container on a desktop opened a browser tab for every network. They are written as before, without being opened
  • A custom reference with a gene across the origin, or with several sequences: PREPARE_REFERENCE named each reference protein after where its pieces lie, so tagA on plasmid pO157 of E. coli O157:H7 Sakai, join(92527..92721,1..2502) in NCBI's GenBank record, became pO157p01_1-92721, the whole plasmid, and lost its stop; from NCBI's GFF3, which writes it 92527..95223, it came out as its first 65 of 898 residues. Bakta's stcE on the reference, 92563..2502, matched no gene and was listed among the false positives, and the calls were compared with the genes of every sequence at once: seven calls on pO157 and two on the chromosome counted as matched by genes of the other at the same numbers, and the GenBank's 1-92721 hid two more. With Sakai as its own reference, the sample with NCBI's headers came out with tagA as 1..2502, 65 codons short and not marked partial (GenBank), or from Bakta's later start (GFF3), and with bare headers the GFF3 lost the 2502 bp piece altogether. A gene across the origin now ends past the sequence's last base, as GFF3 writes it (pO157p01_92527-95223, with its stop), each protein names its sequence (contig=), Bakta keeps the FASTA's names on the reference, and each call is compared with the genes of its own sequence, on both sides of the origin; calls on a sequence the annotation has no gene on are not judged. Both formats now give the same 5155 proteins and the same 260 false positives, and the sample with NCBI's headers gets tagA exactly, 92527..2502. K-12's 142 false positives are unchanged
  • pipeline_info/params.yml is the run's own: with -resume it came back from the cache, so it held the first run's parameters whatever later runs were given, since its script reads params as a whole, which Nextflow 26.10's default parser leaves out of the task's hash. And all but the first parameter fell outside params:, with a stray YAML line at the end, because Nextflow strips a script's indentation once the script is filled in. DUMP_PARAMS now never comes from the cache, and writes each value as a JSON string under params:
  • Bakta's SLURM QoS only on SLURM: Bakta's jobs asked for --qos=medium under every scheduler, though only SLURM's sbatch takes it, and under -profile sge it took the place of the memory request the sge profile gives every task. The slurm profile now sets it for Bakta
  • 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
+
  • Pseudogenes in the GFF3 and GenBank: the GFF3 never said which genes were pseudogenes, and the GenBank only for genes with an internal stop, though the master table knew. Both now carry AMAP's flag with its reason (pseudo=true and Note, /pseudo and /note), and mark the ends a contig cut off (partial=true with start_range or end_range, < or >); the master table and *-AMAP.tsv have them as their last columns
  • The default H37Rv reference under Docker: its files live inside the container and were handed to processes as paths, so Nextflow mounted /reference from the host, which Docker Desktop refuses and Docker on Linux mounts empty over the container's own. BUNDLED_REFERENCE copies them out first
  • -profile test: VALIDATE_INPUTS still required the sample sheet the test profile does without. 1.0.1 had only been checked with -preview, which runs no task
  • -resume kept results made by older code: Nextflow hashes a task's inputs and script, not the files it runs by path, so after an edit to the classifier or any script under bin/, -resume reused the old results. Every process that runs AMAP code now lists those files in ext code, which puts them in its cache key; edit bin/classify/cases.py and TYPE_CASES and FIX_CDS run again while PATHOTYPR and Bakta stay cached
  • PATHOTYPR no longer depends on the classifier: its script said pathotypr classify, and Nextflow adds any bin/ entry a script names as a bare word to the cache key, so a change in bin/classify/ re-ran PATHOTYPR, Bakta and everything after them. The subcommand is quoted
  • Nextflow 25.10 and later: the pipeline stopped before running anything, because the strict parser those versions use by default rejects the helper functions in nextflow.config, the statements outside the workflow in main.nf, the bare env(SAMPLE_NAME) output and publishDir paths that name a task input. Nextflow 26.10 also refuses the return that ended a --validate_only run. The resource caps are now the resourceLimits directive, which still takes --max_cpus, --max_memory and --max_time from the command line, a profile or a -c config. Nextflow 24.04, where that directive first appeared, is the new minimum
  • --skip_mycobrowser in organism-agnostic mode: a custom reference from another genus turned it on by reassigning the parameter, which current Nextflow ignores with a warning, so the run went ahead with Mycobrowser enabled. It is now worked out in the parameter's default
  • Tasks leave the pipeline untouched: the first time a copy of AMAP ran, Python cached the bytecode of AMAP's own modules inside it, in bin/classify/__pycache__, and with pathotypr classify naming that directory in PATHOTYPR's script, the first -resume re-ran PATHOTYPR, Bakta and everything after them. Tasks now run with PYTHONDONTWRITEBYTECODE=1
  • Samples with no Bakta call outside cases 2 and 3: when every Bakta call had an alignment at exactly its coordinates, which a small input can do, TYPE_CASES stopped with a KeyError, because the empty column of overlaps was read as a list of column names. Case 1 and case 4 now come out empty and the rest is classified as usual. An empty miniprot table stopped it the same way
  • Samples Bakta calls no CDS on: with no contig names from Bakta to compare, the contig check took them for different names and stopped the run, and a Bakta table with nothing but its header did not read. The genes miniprot places now come out as missing from Bakta, like any gene it did not call
  • --version: --help and the parameter docs have listed it since 1.0.0, but nothing read it, so nextflow run PathoGenOmics-Lab/AMAP --version started a run. It now prints the pipeline and Nextflow versions and exits before any task
  • Drafts named after one of their contigs: pathotypr writes a summary row for each contig that carries lineage markers, and PATHOTYPR took the sample's name from the first, so CDC1551 cut into 40 contigs, sample CDC1551d, came out as L4_contig_14, and so did every file after it. Two drafts whose first row named the same contig became one sample, the results of one overwriting the other's. A genome is now named <lineage>_<sampleId> from the sample sheet, the lineage of a draft being the one the markers of all its contigs add up to, and the contigs of a draft keep their names instead of the first one taking the sample's. A closed genome whose header is its sample ID comes out as before
  • Headers with a description: the name was the whole header line, so the tutorial's CDC1551 (NC_002755.2 Mycobacterium tuberculosis CDC1551, complete sequence) was named with spaces in it, Bakta failed on it three times and the run finished without it, and without failing. Such a genome now takes its sample ID like any other, and so does the test profile's: L4_H37Rv_test, where it was L4_CP110619_H37Rv
  • An input could be emptied: PATHOTYPR wrote the renamed genome next to the link Nextflow stages its input as, so when the new name was the file's own, sample G26689 in L4_G26689.fasta for instance, the input was overwritten with nothing and the run went on without it. The input is now staged in a directory of its own
  • Sample IDs with underscores: the steps after the correction told a sample's lineage from its sample ID by splitting its name at underscores, though typing only puts the lineage in front and sample IDs may hold underscores of their own. L4_H37Rv_test was counted as genome H37Rv, the tutorial's Beijing_1 and a Beijing_2 as one genome Beijing in the category counts, the report and the duplication network, an untyped H37Rv_test as lineage H37Rv, and each gene of an untyped sample with no underscore as a genome of its own in the network. PATHOTYPR now says which lineage it named a genome after and every step is handed the sample ID and the lineage. The case summary gains a Sample column and names an untyped sample as it is, not Untyped_<sampleId>, and the report no longer reads the first gene not in the reference as a header, nor stops when every gene has one
  • Genomes kept outside the run's folder: VALIDATE_INPUTS staged only the sample sheet and opened the paths written in it from inside its container, which sees only the folder that the work directory and the sheet share. Genomes anywhere else, /data/genomes while running from ~/project for instance, were "FASTA file not found" and the run stopped before anything ran, and a relative path, which the sheet may use, was read from the task directory with or without a container. The genomes are now inputs of the process, each in a directory of its own so that two with the same file name do not collide, and each is read where Nextflow staged it
  • Runs with functional annotation stopped at validation: VALIDATE_INPUTS looked for --eggnog_db in its own task directory, where the default eggnog_db, relative to the launch directory, never is, and took a database not built yet for an error, though EGGNOG_DB builds it. It now looks where EGGNOG_DB builds it, staged when it exists, and stops the run only when eggnog_db.ok is there without a complete database, since Nextflow then never runs EGGNOG_DB again
  • A Bakta database on the host: --bakta_db reached the tasks as a path the container does not mount. With -profile docker, validation called the database missing and stopped the run, and where nothing validated it, as under -profile test, every sample failed three times and was ignored, and the run ended with nothing annotated and exit 0. A database on the host is now staged into BAKTA_ANNOTATOR, PREPARE_REFERENCE and VALIDATE_INPUTS like any other input; the default inside the container is used as before
  • A warning on every run: validation expected a db.sqlite3 that no Bakta database has. It now checks for version.json and bakta.db, which Bakta cannot run without, so their absence stops the run, and names the database's type and version
  • --strict_validation: the documentation said validation warnings could be made to stop the run, but nothing passed the script's --strict on. The new parameter does
  • eggNOG's database under Docker: EGGNOG_DB built it in a folder it worked out for itself from --eggnog_db, and having no inputs, its container saw that folder only when Nextflow happened to mount it with the work directory and the pipeline's bin/, through a parent folder they share. Otherwise, as when running from /data/project with the pipeline in ~/.nextflow, the database was downloaded and built inside the container and lost with it, and only eggnog_db.ok reached the host: later runs skipped EGGNOG_DB, and ANNOT_EGGNOGMAP_BATCH found no bacteria.dmnd. A database already there without eggnog_db.ok was built again. main.nf now makes the folder on the host and hands it to EGGNOG_DB as an input. A folder left with only eggnog_db.ok needs that file deleted, as validation says
  • eggNOG's database in a folder of any name: ANNOT_EGGNOGMAP_BATCH asked eggNOG-mapper for eggnog_db/bacteria.dmnd and --data_dir eggnog_db, but Nextflow stages a folder under its own name, so with a --eggnog_db called anything else, as the HPC guide's /shared/databases/eggnog or the tutorial's ~/databases/eggnog, EGGNOG_DB built the database and the annotation stopped at "DIAMOND database eggnog_db/bacteria.dmnd not present", with or without a container. The folder is now staged under that one name, whatever it is called. The task's tag and first log line also gave a single sample's FAA size in bytes as the number of samples (batch_1519647_samples under -profile test)
  • Reference files on the host: BAKTA_ANNOTATOR gave Bakta --ref_annotation, or --bakta_proteins, as a path written into its script, which a container does not see, nor a task directory when it is relative, so the organism-agnostic mode failed under -profile docker: with K-12 as the reference for Sakai, Bakta called the "user proteins file" not valid three times, the sample was ignored and the run ended with exit 0. Likewise --pathotyper_snps and --pathotyper_ref on the host left every genome untyped, as NA_<sampleId>, --mycobrowser_ref left the Mycobrowser categories out, and one --miniprot_ref_* on the host while the others kept their defaults made BUNDLED_REFERENCE fail, so nothing was corrected. A file on the host is now staged into the tasks that read it, each in a directory of its own, and a default is read inside the container as before
  • Reference files that are nowhere: a typo in one of those paths, or a default on a run without the container, went unnoticed until the samples were ignored or left untyped. Validation now checks each reference file the run reads and that is not on the host, inside the container where the tasks read it, and stops the run when it is not there. A missing Mycobrowser table only warns, since GET_OUTPUTS goes on without it
  • The pangenome without lineage typing: PANGENE aligned each sample's proteins onto the genomes PATHOTYPR renamed, so with --skip_typing, which every custom reference from another genus needs, it had none and wrote empty files for every sample. It now aligns them onto the sample sheet's genomes, each copied to <sampleId>.fasta first, so that two files of the same name (assembly.fasta from two assemblies) stay apart and a .fna is not passed over
  • Genes not in the reference, sample by sample: the report counted only the samples that had some, so a sample whose genes all matched the reference dropped out of its lineage's box instead of counting 0, and when no sample had any (H37Rv against itself, as in -profile test) the section was an empty plot. Every sample is now counted, the section says so when none has any, and each point is drawn at its count, on an axis from 0, where it was jittered up and down by as much as 0.4
  • The report's pangenome section lost a genome: pangene names each genome <name>#0, and the report read the merged matrix with R's default comment character, so its header stopped at the first #. One genome column was dropped and the gene names became row names, so the core/accessory and frequency plots used a single column and the heatmap never appeared. The missing-gene and duplication tables, which carry the reference annotation's gene IDs as given, are now read without a comment character too
  • Genes across the origin of a circular sequence: Bakta reads a sequence as a complete replicon when its FASTA description says complete (NCBI's "complete genome" and "complete sequence") or circular=true (Unicycler), joins its ends and calls a gene across its origin from its start to its stop, start after end. TYPE_CASES stopped on such a gene: E. coli O157:H7 Sakai with NCBI's headers, whose stcE runs across the origin of plasmid pO157 (92563..2502), failed three times and was left out, and the run finished without failing. Past it, check_cds.py dropped such a gene as invalid coordinates, deduplication and the QC table gave it a negative length, the GFF3 wrote it backwards, the GenBank writer stopped on it and the Panaroo GFF3 spanned it over the whole sequence. It now goes through whole and is written as Bakta writes it: start after end in the tables, join() across the origin in a circular GenBank record, and in the GFF3 its end plus the sequence's length on a region marked Is_circular=true. A gene near the origin of such a sequence is read across it too, so Rop on plasmid pOSAK1, of which Bakta called the first 141 bp, comes out whole. Annotated against K-12, Sakai with NCBI's headers now has 5331 genes in common with the same genome under bare names, which Bakta reads as linear, and differs in stcE and Rop, whole where the linear reading has 804 of stcE's 886 residues and no Rop, and in two phage genes Bakta starts elsewhere on a complete genome
  • Typed closed genomes were annotated as linear: PATHOTYPR renamed a genome of one sequence to <lineage>_<sampleId> and nothing else, so the description Bakta reads a complete, circular replicon from never reached it, whatever the header said, and a gene across the origin of a closed genome came out as two pieces under one name, marked 5' and 3' partial. Only the sequence's ID is replaced now, and the description after it stays, as with --skip_typing. Typed from NCBI's headers, CDC1551, CP110619, AL123456.3 and M. bovis get the same genes as before, in place, name, flags and protein, though Bakta starts a few of its own calls elsewhere on a complete genome (dnaA at 1 instead of 34), so the gene IDs that are Bakta locus tags can be numbered differently. H37Rv turned so that the origin cuts gyrB, gyrA or pknB gives back all 4031 GenBank CDS, where it gave 4030 with the cut gene in two pieces. The GenBank record of such a genome is marked circular, with the description in its definition, and its GFF3 has a region marked Is_circular=true
  • The report's duplication summary: it printed the number of homologous community pairs as markdown from a chunk that did not pass its output on as markdown, so the report showed ## **48** homologous community pairs detected across all genomes. in a code box. It is now a sentence. The example report's duplication table had columns GET_ALLDUP never writes, so its Network Metrics tab said there were none; the demo data now has the columns GET_ALLDUP writes and the example is rendered again
  • Modularity in homologous_communities.tsv: Modularity1 and Modularity2 were computed with the whole community as its only group, which makes them 0 whatever the community, as on all 29 rows of CDC1551 against M. bovis. The two columns are gone
  • COG categories in homologous_communities.tsv: with numpy 2, as in the container, Elements1 and Elements2 read {'L': np.int64(4)}, and they counted each protein eggNOG gave no category in category + or -, since its ID keeps the strand FIX_CDS wrote. They are now plain counts, the largest first and ties in alphabetical order rather than in an order that changed from run to run, and a protein without a category counts as -, eggNOG's mark for none
  • The order of -profile: from Nextflow 25.10, whose strict parser applies the profiles in the order -profile names them and lets the last one win, the laptop caps of the container profiles took over from the others: -profile slurm,singularity held every task of a cluster run to the launching machine's CPUs and 16 GB, where -profile singularity,slurm gave it 48 CPUs and 128 GB, -profile test,docker ran the test with 16 GB and -profile highmem,slurm with 128 GB. No profile sets a cap now: they are set after the profiles, from whether highmem, test or a scheduler profile is on, so every order gives what Nextflow 24.04 gave. The parameter docs gave 24 CPUs and 50 GB as the defaults, which no run got
  • The duplication files are the same from one run to the next: homologous_communities.tsv listed its pairs of communities in the order of a set, which changes with Python's hash seed from one run to the next, and with it which community of a pair came first and the last digit of AvgClustering, an average over the community's proteins that networkx went through in the order of their set: 0.5833333333333333 or 0.5833333333333334 for the same pair of CDC1551 and M. bovis. The pairs now come in the order of the communities, the sample first by name on the left, and the metrics go through the proteins in alphabetical order. Each family in duplications_*_groups.tsv lists its proteins in alphabetical order, where they too came in the order of a set, and each sample's two network drawings no longer start from a random layout and a random partition, so every file both steps write is the same run after run
  • No web browser: both duplication steps wrote their networks with pyvis' show(), which also opens the file in a web browser, so a run without a container on a desktop opened a browser tab for every network. They are written as before, without being opened
  • A custom reference with a gene across the origin, or with several sequences: PREPARE_REFERENCE named each reference protein after where its pieces lie, so tagA on plasmid pO157 of E. coli O157:H7 Sakai, join(92527..92721,1..2502) in NCBI's GenBank record, became pO157p01_1-92721, the whole plasmid, and lost its stop; from NCBI's GFF3, which writes it 92527..95223, it came out as its first 65 of 898 residues. Bakta's stcE on the reference, 92563..2502, matched no gene and was listed among the false positives, and the calls were compared with the genes of every sequence at once: seven calls on pO157 and two on the chromosome counted as matched by genes of the other at the same numbers, and the GenBank's 1-92721 hid two more. With Sakai as its own reference, the sample with NCBI's headers came out with tagA as 1..2502, 65 codons short and not marked partial (GenBank), or from Bakta's later start (GFF3), and with bare headers the GFF3 lost the 2502 bp piece altogether. A gene across the origin now ends past the sequence's last base, as GFF3 writes it (pO157p01_92527-95223, with its stop), each protein names its sequence (contig=), Bakta keeps the FASTA's names on the reference, and each call is compared with the genes of its own sequence, on both sides of the origin; calls on a sequence the annotation has no gene on are not judged. Both formats now give the same 5155 proteins and the same 260 false positives, and the sample with NCBI's headers gets tagA exactly, 92527..2502. K-12's 142 false positives are unchanged
  • pipeline_info/params.yml is the run's own: with -resume it came back from the cache, so it held the first run's parameters whatever later runs were given, since its script reads params as a whole, which Nextflow 26.10's default parser leaves out of the task's hash. And all but the first parameter fell outside params:, with a stray YAML line at the end, because Nextflow strips a script's indentation once the script is filled in. DUMP_PARAMS now never comes from the cache, and writes each value as a JSON string under params:
  • Bakta's SLURM QoS only on SLURM: Bakta's jobs asked for --qos=medium under every scheduler, though only SLURM's sbatch takes it, and under -profile sge it took the place of the memory request the sge profile gives every task. The slurm profile now sets it for Bakta
  • 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

🧹 Code#

@@ -104,7 +104,7 @@

🧹 Code📚 Documentation#

  • A website in plain HTML, drawn with real genes: MkDocs is gone. The front page starts a newcomer in three steps and builds the command for their machine, the pipeline is a metro map whose stations are main.nf's processes, and the four cases, the missing genes and the harder cases (paralogs, small proteins, contig ends, frameshifts, copies) are drawn on genes of CDC1551, M. bovis and E. coli from the benchmark. The tutorial and user guide stay markdown, readable on GitHub, and tools/build_docs.py turns them into pages of the same site, dark by default
  • The site is checked against the code: tools/check_site.py runs every drawn window through bin/, compares the thresholds on the pages with the constants in the code, the metro map with main.nf and every option with the pipeline's parameters, so a change that makes the site wrong fails the docs workflow. tools/figure_data.py redraws the windows from a benchmark run
  • What the docs said that the pipeline does not do: with --skip_functional_annotation a run still writes AMAP/, with the corrected table and sequences, where the outputs page said it did not; no step writes PATHOTYPR/renamed_fastas.tsv; and the 19-residue toxin the algorithm page names is IbsC (b4665), not an ldr toxin
  • A drawing on every page that explains something: the duplication family and communities of CDC1551 and M. bovis, from a real run of both duplication steps; the sample sheet with its tabs and the messages the validator prints; the ranges a genome is validated against; the container and what it holds; a run on a cluster; what happens to a task that fails; the report section by section; the five outcomes and the three inputs of a decision; the command builder on the quickstart. Most are built from the code or the config, so they change when it does
  • The benchmark on the front page: the seven genomes whose answer is known, and how many of their genes AMAP finds, checked against tests/benchmark/README.md. A shared link to any page shows a card drawn from the front page's genes
  • More that the docs got wrong: the duplication page listed ten network metrics, the table has four, and its columns ending in 1 and 2 are the two communities of a pair, not two levels; Nextflow needs Java 17, not 11; --skip_functional_annotation leaves no GFF3 or GenBank, which come after eggNOG; tasks are capped at 16 GB on one machine and at 128 GB under a scheduler profile, not 50 GB; and lineage typing is not skipped on its own for other species, only Mycobrowser is
  • The tutorial's numbers are today's: its case counts and QC line, and their copies on the front page, the outputs page and the custom-reference page, came from runs in July. They now come from the tutorial's own commands run again on 18 October 2026 with --skip_functional_annotation, the QC line from GET_OUTPUTS' own script run on that run's corrected genes: CDC1551 has 4014 genes and 123 pseudogenes as AMAP calls them, where the old line counted Mycobrowser's 13, and Sakai keeps 1387 genes K-12 lacks, not 1356. The tutorial also says where --validate_only writes its verdict and where the number of reference proteins extracted is, as neither reaches the console
  • The README in the site's style: the metro map and the five outcomes in light and dark, the benchmark's scores, and the links a newcomer needs first
  • The README opens with the logo, and its drawing plays: the picture at the top was a still of the front page's drawing, dark whatever GitHub's theme. It now plays the drawing's four steps, in light or dark as the theme is, made by tools/pictures.mjs from the front page's own figure. The README's metro map gains the interactive report's station
  • One DOI: the README, CITATION.cff and every page of the website cited the DOI of one version on Zenodo, and the manifest, main.nf and the interactive report the concept DOI, which always resolves to the latest version. All of them cite the concept DOI now
  • Bakta's options on a SLURM cluster: the HPC page says how to give Bakta a site's options, since for Bakta the QoS the slurm profile sets takes the place of -process.clusterOptions
  • The website is served as it is built: GitHub Pages serves docs/ from main, and by default runs Jekyll over it first, which would have made a page of each markdown file with front matter at the address of the page tools/build_docs.py builds from it. docs/.nojekyll turns Jekyll off, and tools/check_site.py fails without it

🧪 Testing#

-
  • tests/selftest/h37rv_selftest.sh annotates H37Rv with itself inside the container and fails unless every GenBank CDS comes back exactly
  • tests/test_cache_keys.py fails when a process runs AMAP code that is not in its cache key, or when its script names a bin/ entry as a bare word without declaring it
  • tests/benchmark/benchmark.sh scores the correction algorithm on seven genomes whose right answer is known (H37Rv twice, CDC1551 complete and in 40 contigs, M. bovis, E. coli K-12 and Sakai), downloading them and running Bakta and miniprot in the container the first time, and lists every gene that changed since an earlier run. A gene whose sequence changed counts as found when a gene with its name covers half of its own bases there, not half of the stretch it lands on, which an insertion can make far longer (63 kb with a prophage inside for b1004 in Sakai); a gene mostly deleted is not counted. Changed genes found: 50 of 50 in CP110619, 860 of 862 in CDC1551, 1515 of 1515 in M. bovis and 3363 of 3367 in Sakai, and the README says why each of the other six is not
  • tests/test_name_by_lineage.py runs PATHOTYPR's own script, rendered as Nextflow renders it, with a stand-in for pathotypr that writes the summaries pathotypr wrote for CDC1551 complete and in 40 contigs, and fails unless Bakta would read a renamed genome as circular exactly when it reads the original so
  • tests/test_sample_ids.py runs TYPE_CASES, GET_OUTPUTS, GET_COUNT, GET_DUPLICATIONS, REPORT_INPUTS and GET_ALLDUP as Nextflow renders them, chained as the pipeline chains them, for a typed and an untyped run whose sample IDs hold underscores, and then the report's own R code where R is installed, and FINAL_REPORT's own script on a run's tables and on one that leaves every table it can empty
  • tests/test_validate_inputs.py runs VALIDATE_INPUTS's own script, rendered as Nextflow renders it, in a task directory laid out as Nextflow lays one out: genomes in a folder the task cannot see, relative paths, two genomes with the same file name and a sheet of 500 samples; and checks both databases, eggNOG's in each state EGGNOG_DB can leave it in
  • 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_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
  • tests/test_qc_pseudogenes.py runs the steps after the correction as tests/test_sample_ids.py chains them, with FIX_CDS calling some genes pseudogenes and Mycobrowser marking another, with Mycobrowser on and off, and checks that *-AMAP.qc.tsv, REPORT_INPUTS and the report's own QC table count the pseudogenes AMAP called, as many as the GFF3 and the GenBank mark
  • tests/test_duplication_reproducible.py runs get_duplications.py and network_genomes.py under three hash seeds, on the Diamond tables of three genomes with a family whose average clustering changes with the order of its proteins, and fails unless every file they write is the same under all three; it also fails when either script opens a web browser, with hits or without, as a stand-in browser that records what it is asked to open tells, and when a network page asks a CDN for its drawing library
  • tests/test_custom_reference.py extracts a reference whose plasmid has a gene across its origin from a GenBank record and from GFF3 written either way, on either strand, compares Bakta's calls with it sequence by sequence, and runs PREPARE_REFERENCE's own script, rendered as Nextflow renders it, with a stand-in for Bakta that names sequences as Bakta does
  • tests/test_amap_report.py writes the report from four genomes of a real run cut to two stretches of the reference (tests/data/report/, rebuilt by its make_fixture.py), where M. bovis lacks RD4 and a contig end cuts dnaE2, plus a gene across the origin of a circular plasmid, and checks that every gene, and every Bakta call AMAP did not keep, comes back out of the page as it went in, with Python and, where node is installed, with the page's own JavaScript, along with the states, the flags, the correction's steps and what the process hands the script
  • tests/test_dump_params.py runs DUMP_PARAMS with Nextflow where it is installed, under both parsers, with values holding quotes, $, backticks, a backslash and a line break, and again with -resume and one value changed; everywhere, it fails if the step can come from the cache
  • tests/test_cluster_options.py reads which profile gives Bakta the QoS, and with Nextflow resolves Bakta's and every task's cluster options for each scheduler, in both orders with a container and under both parsers
  • tests/test_environment_pins.py fails when environment.yml pins a tool at a version docker/Dockerfile does not install
  • 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
+
  • tests/selftest/h37rv_selftest.sh annotates H37Rv with itself inside the container and fails unless every GenBank CDS comes back exactly
  • tests/test_cache_keys.py fails when a process runs AMAP code that is not in its cache key, or when its script names a bin/ entry as a bare word without declaring it
  • tests/benchmark/benchmark.sh scores the correction algorithm on seven genomes whose right answer is known (H37Rv twice, CDC1551 complete and in 40 contigs, M. bovis, E. coli K-12 and Sakai), downloading them and running Bakta and miniprot in the container the first time, and lists every gene that changed since an earlier run. A gene whose sequence changed counts as found when a gene with its name covers half of its own bases there, not half of the stretch it lands on, which an insertion can make far longer (63 kb with a prophage inside for b1004 in Sakai); a gene mostly deleted is not counted. Changed genes found: 50 of 50 in CP110619, 860 of 862 in CDC1551, 1515 of 1515 in M. bovis and 3363 of 3367 in Sakai, and the README says why each of the other six is not
  • tests/test_name_by_lineage.py runs PATHOTYPR's own script, rendered as Nextflow renders it, with a stand-in for pathotypr that writes the summaries pathotypr wrote for CDC1551 complete and in 40 contigs, and fails unless Bakta would read a renamed genome as circular exactly when it reads the original so
  • tests/test_sample_ids.py runs TYPE_CASES, GET_OUTPUTS, GET_COUNT, GET_DUPLICATIONS, REPORT_INPUTS and GET_ALLDUP as Nextflow renders them, chained as the pipeline chains them, for a typed and an untyped run whose sample IDs hold underscores, and then the report's own R code where R is installed, and FINAL_REPORT's own script on a run's tables and on one that leaves every table it can empty
  • tests/test_validate_inputs.py runs VALIDATE_INPUTS's own script, rendered as Nextflow renders it, in a task directory laid out as Nextflow lays one out: genomes in a folder the task cannot see, relative paths, two genomes with the same file name and a sheet of 500 samples; and checks both databases, eggNOG's in each state EGGNOG_DB can leave it in
  • 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 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
  • tests/test_qc_pseudogenes.py runs the steps after the correction as tests/test_sample_ids.py chains them, with FIX_CDS calling some genes pseudogenes and Mycobrowser marking another, with Mycobrowser on and off, and checks that *-AMAP.qc.tsv, REPORT_INPUTS and the report's own QC table count the pseudogenes AMAP called, as many as the GFF3 and the GenBank mark
  • tests/test_duplication_reproducible.py runs get_duplications.py and network_genomes.py under three hash seeds, on the Diamond tables of three genomes with a family whose average clustering changes with the order of its proteins, and fails unless every file they write is the same under all three; it also fails when either script opens a web browser, with hits or without, as a stand-in browser that records what it is asked to open tells, and when a network page asks a CDN for its drawing library
  • tests/test_custom_reference.py extracts a reference whose plasmid has a gene across its origin from a GenBank record and from GFF3 written either way, on either strand, compares Bakta's calls with it sequence by sequence, and runs PREPARE_REFERENCE's own script, rendered as Nextflow renders it, with a stand-in for Bakta that names sequences as Bakta does
  • tests/test_amap_report.py writes the report from four genomes of a real run cut to two stretches of the reference (tests/data/report/, rebuilt by its make_fixture.py), where M. bovis lacks RD4 and a contig end cuts dnaE2, plus a gene across the origin of a circular plasmid, and checks that every gene, and every Bakta call AMAP did not keep, comes back out of the page as it went in, with Python and, where node is installed, with the page's own JavaScript, along with the states, the flags, the correction's steps and what the process hands the script
  • tests/test_dump_params.py runs DUMP_PARAMS with Nextflow where it is installed, under both parsers, with values holding quotes, $, backticks, a backslash and a line break, and again with -resume and one value changed; everywhere, it fails if the step can come from the cache
  • tests/test_cluster_options.py reads which profile gives Bakta the QoS, and with Nextflow resolves Bakta's and every task's cluster options for each scheduler, in both orders with a container and under both parsers
  • tests/test_environment_pins.py fails when environment.yml pins a tool at a version docker/Dockerfile does not install
  • 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
  • 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#

Bug fixes, most of them found by running the pipeline on real genomes rather than by reading it. Several changed what the pipeline outputs, so results from 1.0.0 are not directly comparable.

🔬 Correctness#

diff --git a/docs/assets/figures.js b/docs/assets/figures.js index 92efa98..b2d6833 100644 --- a/docs/assets/figures.js +++ b/docs/assets/figures.js @@ -778,8 +778,8 @@ }, PANGENE: { name: "Pangenome", proc: "PANGENE", x: 930, y: 320, label: "below", lineOnly: "pan", - text: "Each sample's corrected proteins aligned with miniprot onto every other sample's genome, then a pangene graph and a presence/absence matrix.", - reads: "the corrected proteins, and every genome: as PATHOTYPR renamed it, or with --skip_typing as the sample sheet gives it, under its sample ID", writes: "PANGENE/{sample}/", flags: "", doc: "how-it-works/workflow.html#9-pangenome", exit: "PANGENE/", exitSide: "right", functional: true + text: "The reference proteome and the genes no reference gene matches, aligned with miniprot to every genome, its own included, then one pangene graph for the run and a presence/absence matrix with a column per sample.", + reads: "the reference proteome, the corrected genes, and every genome: as PATHOTYPR renamed it, or with --skip_typing as the sample sheet gives it, under its sample ID", writes: "PANGENE/pangenome.gfa, .bubbles.txt, .Rtab, .stat", flags: "", doc: "how-it-works/workflow.html#9-pangenome", exit: "PANGENE/", exitSide: "right", functional: true }, FINAL_REPORT: { name: "Report", proc: "FINAL_REPORT", x: 1050, y: 170, kind: "hub-end", label: "above", @@ -950,7 +950,7 @@ var text = state.custom && s.textCustom ? s.textCustom : s.text; var status = ""; if (stationOff(key, state)) status = '

Not run with the switches above.

'; - if (key === "PANGENE" && state.skipTyping && !state.skipFunctional) status = '

With --skip_typing it aligns onto the genomes of the sample sheet, each under its sample ID.

'; + if (key === "PANGENE" && state.skipTyping && !state.skipFunctional) status = '

With --skip_typing it searches the genomes of the sample sheet, each under its sample ID.

'; info.innerHTML = "

" + esc(s.name) + "" + esc(proc) + "

" + esc(text) + "

" + status + (s.doc ? '

More on this step

' : "") + "
" + diff --git a/docs/assets/search-index.js b/docs/assets/search-index.js index 9bea74a..a44967e 100644 --- a/docs/assets/search-index.js +++ b/docs/assets/search-index.js @@ -404,6 +404,12 @@ window.AMAP_DOCS_INDEX = [ }, { "page": "Outputs", +"section": "PANGENE/pangenome.Rtab", +"url": "user-guide/outputs.html#pangenepangenomertab", +"text": "The gene presence/absence matrix of the whole run, the one both reports draw the pangenome from: one row per gene, one column per sample, 1 where pangene found the gene in that sample's genome. Every genome is searched with the same proteins, the reference proteome and the genes no reference gene matches, so a gene with a 1 in every column is one all the samples share. A gene is named after its reference gene (Rv1506c) or, when no reference gene matches it, after one sample's copy of it (L4_CDC1551_01396, the ID in that sample's master table)." +}, +{ +"page": "Outputs", "section": "final_report.html", "url": "user-guide/outputs.html#final_reporthtml", "text": "Self-contained, so it travels by email. Covers per-sample quality, missing reference genes by lineage, PE/PPE family distribution, KEGG pathways, duplication communities and network topology, the pangenome (core, softcore, shell, cloud, frequency spectrum, presence/absence heatmaps), COG and Mycobrowser distributions, case statistics and the associated tests, plus a searchable master-table explorer." @@ -682,7 +688,7 @@ window.AMAP_DOCS_INDEX = [ "page": "The workflow", "section": "9. Pangenome", "url": "how-it-works/workflow.html#9-pangenome", -"text": "pangene builds a pangenome graph from the corrected annotations and produces a gene presence/absence matrix, which the report turns into core, softcore, shell and cloud classifications, a frequency spectrum and presence/absence heatmaps. Each sample's corrected proteins are aligned with miniprot onto every other sample's genome: the ones PATHOTYPR renamed or, with --skip_typing, the genomes of the sample sheet, each under its sample ID." +"text": "pangene builds one pangenome graph for the whole run and a gene presence/absence matrix, which the report turns into core, softcore, shell and cloud classifications, a frequency spectrum and presence/absence heatmaps. Every genome is searched with the same proteins: the reference proteome, and the genes of each sample that no reference gene matches. miniprot aligns them to every sample's genome, its own included: the ones PATHOTYPR renamed or, with --skip_typing, the genomes of the sample sheet, each under its sample ID. A row of the matrix is then a gene, named after its reference gene or, ou" }, { "page": "The workflow", diff --git a/docs/how-it-works/workflow.html b/docs/how-it-works/workflow.html index 2324d7f..b2e70c3 100644 --- a/docs/how-it-works/workflow.html +++ b/docs/how-it-works/workflow.html @@ -121,8 +121,8 @@

8. Duplication analysisDuplication metrics explains the numbers.

9. Pangenome#

-

pangene builds a pangenome graph from the corrected annotations and produces a gene presence/absence matrix, which the report turns into core, softcore, shell and cloud classifications, a frequency spectrum and presence/absence heatmaps.

-

Each sample's corrected proteins are aligned with miniprot onto every other sample's genome: the ones PATHOTYPR renamed or, with --skip_typing, the genomes of the sample sheet, each under its sample ID.

+

pangene builds one pangenome graph for the whole run and a gene presence/absence matrix, which the report turns into core, softcore, shell and cloud classifications, a frequency spectrum and presence/absence heatmaps.

+

Every genome is searched with the same proteins: the reference proteome, and the genes of each sample that no reference gene matches. miniprot aligns them to every sample's genome, its own included: the ones PATHOTYPR renamed or, with --skip_typing, the genomes of the sample sheet, each under its sample ID. A row of the matrix is then a gene, named after its reference gene or, outside the reference, after one sample's copy of it, and a column a sample, so the genes the samples share are the rows found in every column.

10. Report#

A single self-contained HTML file built from every output above. Sample quality, missing reference genes by lineage, PE/PPE families, KEGG pathways, duplication communities and topology, the pangenome, COG and Mycobrowser distributions, case statistics with the associated tests, and a searchable table of every gene.

Alongside it, pipeline_info/ collects the Nextflow execution report, timeline, DAG and trace, plus software_versions.yml and params.yml for your methods section.

diff --git a/docs/how-it-works/workflow.md b/docs/how-it-works/workflow.md index 9858095..4918013 100644 --- a/docs/how-it-works/workflow.md +++ b/docs/how-it-works/workflow.md @@ -113,14 +113,18 @@ community. ## 9. Pangenome -[pangene](https://github.com/lh3/pangene) builds a pangenome graph from the -corrected annotations and produces a gene presence/absence matrix, which the -report turns into core, softcore, shell and cloud classifications, a frequency -spectrum and presence/absence heatmaps. - -Each sample's corrected proteins are aligned with miniprot onto every other -sample's genome: the ones PATHOTYPR renamed or, with `--skip_typing`, the -genomes of the sample sheet, each under its sample ID. +[pangene](https://github.com/lh3/pangene) builds one pangenome graph for the +whole run and a gene presence/absence matrix, which the report turns into core, +softcore, shell and cloud classifications, a frequency spectrum and +presence/absence heatmaps. + +Every genome is searched with the same proteins: the reference proteome, and +the genes of each sample that no reference gene matches. miniprot aligns them +to every sample's genome, its own included: the ones PATHOTYPR renamed or, with +`--skip_typing`, the genomes of the sample sheet, each under its sample ID. A +row of the matrix is then a gene, named after its reference gene or, outside +the reference, after one sample's copy of it, and a column a sample, so the +genes the samples share are the rows found in every column. ## 10. Report diff --git a/docs/user-guide/outputs.html b/docs/user-guide/outputs.html index 8475c87..3604f47 100644 --- a/docs/user-guide/outputs.html +++ b/docs/user-guide/outputs.html @@ -133,11 +133,10 @@

The full treeWhat a healthy summary looks like
CDC1551 vs H37Rv4,077 Bakta calls
case 2 2,813case 3 24case 4 1,100case 1, kept 11missing, put back 97
Sakai vs K-125,367 Bakta calls
case 2 3,324case 3 31case 4 507case 1, kept 1,387missing, put back 133
The same table, drawn. The number that differs by two orders of magnitude is case 1: the prophages and the LEE island that O157:H7 carries and K-12 does not.

The one number that differs by two orders of magnitude is case 1, and that is biology: CDC1551 is the same species as the reference so almost nothing is genuinely new, while O157:H7 carries prophages and the LEE island that K-12 does not have.

What should worry you is the opposite: case 2 at or near zero. Every gene disagreeing with the reference by a few bases means something systematic is wrong with the reference, not with 4000 individual genes. The usual cause is a reference proteome whose coordinates do not describe the sequence its proteins were translated from.

+

PANGENE/pangenome.Rtab#

+

The gene presence/absence matrix of the whole run, the one both reports draw the pangenome from: one row per gene, one column per sample, 1 where pangene found the gene in that sample's genome. Every genome is searched with the same proteins, the reference proteome and the genes no reference gene matches, so a gene with a 1 in every column is one all the samples share. A gene is named after its reference gene (Rv1506c) or, when no reference gene matches it, after one sample's copy of it (L4_CDC1551_01396, the ID in that sample's master table).

final_report.html#

Self-contained, so it travels by email. Covers per-sample quality, missing reference genes by lineage, PE/PPE family distribution, KEGG pathways, duplication communities and network topology, the pangenome (core, softcore, shell, cloud, frequency spectrum, presence/absence heatmaps), COG and Mycobrowser distributions, case statistics and the associated tests, plus a searchable master-table explorer.

amap_report.html#

@@ -198,7 +199,7 @@

Next stepsEdit this page on GitHub

- +