Repository navigation
Call gene conversion only from reads that can carry it - #6
Merged
Merged
Conversation
freebayes-parallel ships inside the freebayes package, and its dependencies do not. Read out of the published image, it pipes GNU parallel into vcffirstheader, vcfstreamsort and vcfuniq, and none of those four are in the image: --freebayes_parallel is a documented option, so turning it on failed inside the container with a command not found. Found by listing the image's own bin/ against every binary the pipeline invokes. It is off by default, so no run has hit it, and the image needs a rebuild before the fix is real.
The first rebuild in two months failed, and neither cause was in this file's own logic. Debian bullseye left LTS on 2026-08-31: its security pool is being pruned while the index still advertises what was in it, so apt resolved curl, git and glibc to versions whose .deb files return 404. Base moved to bookworm, where all nine packages this image installs are present. The second would have failed on the next line. micromamba.snakepit.net is an alias served the *.mamba.pm wildcard certificate, which does not cover it, and the one it gets expired on 2026-09-15. Micromamba now comes from the project's own release channel, which takes a third-party CDN out of the build.
The first image in two months, and the first that contains what the pipeline actually invokes: its bin/ was listed against every binary the modules call, including the four that freebayes-parallel runs and that the previous image shipped without. 752 binaries against the old 607.
The digest is repeated in three files and grep found the third only after the first two were changed, which is the argument for a test that checks they agree rather than a habit of remembering.
The pinned digest is written down in six files and only one of them does anything: nextflow.config is what every run pulls, the rest are documentation. A stale copy there is worse than no mention, because it tells a reader they are running an image they are not. Updating them by grep found the last three only after the first three had been changed, which is the argument for a test rather than for remembering. The test compares every digest written anywhere in the docs against the one nextflow.config pins, accepting the abbreviated forms the prose uses, and a second test guards the guard: one that finds nothing to compare passes for the wrong reason.
The pinned reference carries the `docker://` scheme. Singularity requires it, Docker refuses it, and the documented laptop command `-profile local,docker` died with `invalid reference format` before pulling anything. Found while running the fixture cohort through the rebuilt image, which needed the scheme stripped by hand to run at all. The default keeps the scheme, because that is the path every cluster run takes, and the docker profile strips it. Dropping it from the default would fix a laptop and break every cluster, which is the more expensive direction. Verified by running the fixture cohort end to end under -profile local,docker with no override: 67 tasks, none failed. The test asserts both halves and that the two forms name the same image, since a profile that quietly resolved a different tag would satisfy either half alone.
Garnatxa's own documentation submits the Nextflow driver with sbatch rather than leaving it on the login node, where a long run is both against the usage rules and liable to be killed. This is that script for this pipeline, next to the site config it pairs with. The --time is explicit and that is the point: every QoS on that cluster defaults to six hours, and a driver killed at its limit orphans whatever it was waiting on. medium allows seven days. Per-task resources stay in conf/garnatxa.config. The script only carries what a launch decides: the samplesheet, the output directory, the revision, and whether gene conversion runs. The Kraken2 database needs no flag, because the profile already points at the shared copy; it needs a taxId column instead.
/storage/pgo is where the group's sequencing data is, and the Garnatxa profile did not bind it. That does not fail at submission. It fails inside the task, as a file plainly present on the login node and absent to the tool opening it, which is a slow way to be told about a mount. --extra_binds covers the next path nobody anticipated without editing a tracked config. It is declared in nextflow.config rather than beside the profile that consumes it, because main.nf validates flags against that block and rejected it as unrecognised when it was declared next to its use. Whether a command-line param reaches a GString in a profile was checked with a standalone reproducer rather than assumed: it does.
nucmer records the absolute paths of its two inputs in the delta's first line, and show-snps opens them to report the bases around each difference. Under a scheduler that gives each task its own scratch, that path belonged to another task on another node and is gone, so PARALOG_MAP died on a delta that was perfectly intact: show-snps failed (1): ERROR: Could not open file /scr/.../nxf.XXXX/reference.fa Found on a real cluster run and reproduced here by deleting the directory the alignment ran in. It had never been exercised: the chain tests stand MUMmer in with committed text, so no test ever handed show-snps a delta it had to open a file for. paralog_map.py now rewrites that one line to a file that is present, and the reference travels with the delta for that purpose. Everything after the header is offsets and is copied untouched. On the reference that hit this, the stage then finds 264 paralogous pairs and 6610 diagnostic sites.
The three Zenodo files live inside the container image, so declaring them as `path` inputs told Nextflow to stage them from the host and it added a bind for /opt/pathotypr, which Singularity refused because that directory does not exist on the cluster. As values they are interpolated into the script and resolved inside the container, where they do exist. A static test now derives the in-image parameters from the config and follows each call site to the input it lands on, so the same mistake on a future bundled asset fails before a run does.
… used Catalogue v1.0.0 gave each variant a single drug inherited from its gene and took the grade from the first catalogue row it found, so amikacin never appeared at all and 15,969 rows carried a grade belonging to another drug. v1.0.2 grades per variant-drug pair; the ancestor coordinate frame is the one bundled, because all of its REF alleles match the reference beside it and the H37Rv file mismatches 609 of them at the same coordinates. The report now expands a composite label into every drug it names and marks a detected-but-ungraded call instead of rendering it as a clean cell, the build verifies the downloads by checksum with --strict, and DUMP_VERSIONS records the catalogue version and its SHA-256 so the question this started from is answerable from the results next time.
The rebuilt image was verified by running it: /opt/pathotypr/VERSION names v1.0.2, the marker file hashes to the ancestor-frame file published on Zenodo, amikacin is present, the composite labels are there and no row carries the ungraded value any more. Six places write the digest down and four of them abbreviate it, which the consistency test caught after the two full copies were already updated.
FastP is run with --merge, so an overlapping pair leaves as one merged fragment in the orphan stream. MAP_READS has taken all three read streams from the start and LINEAGE_TYPING took two, so pathotypr typed whatever stayed unmerged: 8-19 percent of the reads in the cohort that surfaced this, and zero on a fully overlapping library, which reports Unclassified exactly as a genuinely unclassifiable sample does. The streams are concatenated and typed as one sample, since pathotypr rejects three -i arguments either way. That multi-member gzip is read whole, and that --paired only groups files, were both measured against the pinned binary. The join fails loudly on a missing partner, because the failure it replaces was silent.
Two sequential gzip calls closed the task, and gzip is single-threaded: 72 seconds per sample compressing 480 MB per mate while eleven of twelve reserved cores idled. bgzip with the task's threads does it in 2.8 seconds, produces valid gzip with identical bytes, and fastp reads the result whole, which is the check BGZF needs since it is a multi-member file. The 12-CPU and 80 GB requests came down to 8 and 56 GB, measured from a run that averaged 102 percent CPU and peaked at 39.7 GB, and Garnatxa tasks now declare a walltime instead of inheriting the six-hour QoS default that kept every job out of SLURM's backfill.
Three processes asked for a longer QoS. Two of them were measured at about twenty seconds each and never needed it, and the third only did because Kraken was paging a 133 GB database. Both queues now default to short and the long steps ask for one hour instead of two. The parameters live in the root config, not in the profile: main.nf validates command-line parameters against that block, so one declared only in a profile is refused. Verified by running a reproducer that prints task.clusterOptions, because the session config and what SLURM actually receives are resolved separately.
It was the last fixed memory directive in the pipeline, and on a 185-sample cohort it left 61 jobs queued behind QOSMaxMemoryPerUser while the cluster had capacity. With --memory-mapping the resident set follows the database: 34.7-39.7 GB measured against a 133 GB index, under 20 GB against the capped 16 GB one that answers the same question. --kraken_memory now sets it and it doubles on each retry, matching the eight other memory directives here, so too small a value costs a retry and not the run.
…M from inside it samtools sort -m is per thread and MERGE_AND_MARKDUP runs two sorts at once, each sized at 80 percent of the task's memory: 13.1 GB asked of an 8 GB request. A single run's BAM never fills those buffers, so it hid until the first merged sample streamed 2.8 GB through it. The kill was also unreadable from outside. A stage dies, its stdout closes mid-BGZF block, the next stage reports a short read and exits 1, and the OOM retry never fires. PIPESTATUS is now inspected and a signal death re-raised as 137, verified against a pipe whose middle stage kills itself, with a genuine non-zero exit still failing fast.
submitRateLimit binds whenever the tasks are shorter than the interval between submissions, and a 185-sample cohort is mostly short tasks: 740 of its ~4,300 are annotation steps of two or three seconds. At 50 a minute those spent a quarter of an hour on sbatch calls alone. Raised to 200 a minute. It is a courtesy limit toward the controller, not a correctness one, so the comment says to lower it again if submissions are refused.
A GFF built from a table writes pandas' NaN as four ordinary characters, and the attribute cascade only asked whether a key was present. On a real reference 3,001 CDS lines carry Name=nan next to a usable ID=Rv0001_1-1524 and none carry locus_tag, so Name won and every tract in a 185-sample cohort was labelled nan, with amino-acid changes reported as nan:P1443A. A placeholder value now counts as absent and the cascade falls through to the ID. The first guess at the cause was wrong: str() on a float NaN in the cohort writer would have shown up in every column, and only this one was affected, which is what pointed at the parser instead.
22 samples of a 185-sample cohort were annotated as one lineage, were really another, and were routed to the wrong reference. They carried 2,000 SNPs where their line mates carried 5 and produced 507 of the run's 522 gene-conversion tracts, because a diverged genotype carries the donor base at every paralogous locus by inheritance. Nothing failed; the metrics were correct about the wrong comparison. LINEAGE_MISMATCH flags a sample whose k-mer lineage disagrees with the others sharing its reference, which is independent evidence since k-mer typing ignores the reference. divergent_sample refuses to read tracts as conversion when one sample calls them in more than --max-locus-frac of its loci, measured at 0.20-0.39 for the mislabelled samples and at most 0.08 for every correct one. Both floors are absolute as well as proportional: a ratio over three loci says nothing, and firing there broke 61 existing tests before the floors went in.
… can fire The flag shipped working and never fired. qc_report computed the per-reference majority from summ[sid]["reference"], a key the summary has never had, so the majority was always empty; on the cohort it was written for it flagged nothing while 18 samples sat on the wrong reference. The reference now comes from the samplesheet's refId. The unit tests passed the majority in by hand and so never exercised the lookup; the new ones read it the way qc_report does, and one fails if the majority goes back to depending on the summary.
…t in The report took each Kraken2 report as a sample of its own, so a cohort of 185 samples sequenced over several runs gave 225 rows, and it took the top species as a sample's primary taxon. Kraken2 leaves most M. tuberculosis reads on the complex node, where k-mers cannot tell its members apart, so the species held about 5% of the reads and every clean culture fell under the 90% line: the report called 225 of 185 samples possibly contaminated, and the cultures that really were another organism looked no different. The reports of one sample are now summed by their reads. Each sample is placed in its dominant clade, found by following the most abundant child from the root while it holds at least half of its parent's reads: that stops at the complex for a TB culture and never goes below a species. The cohort's target is the clade most samples are dominated by, and target_pct, the share of a sample's classified reads inside it, is what contamination is now judged on.
…ed to The summary's date column is DAT_OUT, stamped with the day the pipeline processed the sample. The report read it as a collection date, so a whole cohort was sampled in the year it was run and the temporal panel reported a span of zero years. The date now comes only from a samplesheet column that names one (collection_date, sampling_date, date, year...), and a cohort without such a column has no temporal panel. Each sample now also carries its refId and the lineage the other samples on that reference agree on. That is what LINEAGE_MISMATCH compares, so the report can say what disagreed instead of only naming the flag.
On a 185-sample cohort the report was one scroll of 27 panels, about 150,000 pixels long. The SNP dynamics drew 800 trajectory cards on their own, the resistance panel listed all 4,986 calls, the table of all samples had 48 columns with a bar in every cell, and every flag was a code. Worse, several read-outs said things the data did not: 176 "mixed-lineage" samples above a table that counted none, "reference-bias suspects" that were passing samples mapped to their own ancestor, lineage markers counted as resistance, and Mapped % against Unmapped % as the headline correlation. The report now reads as seven pages, listed in the sidebar in the order a reader needs them, with a badge where a page needs attention: - Summary: what the run shows, in sentences with their numbers. Which samples to exclude and why, which are another organism or look mapped to the wrong reference, the resistance mutations beyond each lineage's own, the alleles that swept, the gene conversion called. Each finding links to the page that holds its evidence and says what it cannot prove. - Sample QC: the verdicts, the live thresholds, the flagged samples with the value behind each flag (Low depth 2.2x, needs 10x), the exclusion list, the table of all samples, lineages, contamination, distributions. - Genome & genes, Variants over time, Resistance, Gene conversion, and Diagnostics for the views that never flag a sample on their own. The read-outs are fixed at the source: mixed samples are counted from the MIXED flag; flags the page cannot recompute, such as LINEAGE_MISMATCH, are kept when it re-flags instead of disappearing; the reference-bias screen switches itself off when the cohort sits on its reference; a mutation that at least 90% of a lineage carries is listed as that lineage's marker; pairs that are one quantity measured twice never headline a correlation; and divergent_sample, the most common gene-conversion verdict on the cohort that introduced it, has a chip instead of reading "verdict not recognised". The table of all samples opens on the 14 metrics a verdict is made of and colours only the cells that tripped a flag. Resistance is grouped by mutation and opens on grades 1-2, contamination is one row per sample with the flagged ones first, and the dynamics draw 24 cards at a time. The scatter no longer draws its top values above the axis, the genome landscape key matches its track, and a printed report lays out every page without the mobile sidebar's backdrop over it.
The QC report page described one dashboard of panels. It now walks the seven pages, what each answers, and how the results read: flags in words with the codes kept in every export, pipeline flags that survive re-flagging, contamination measured against the cohort's target clade, lineage markers set apart from the resistance a cohort acquired, and the big panels opening small. The samplesheet's collection-date column gets a row of its own, LINEAGE_MISMATCH joins the flag codes, the exclusion basket is the exclusion list everywhere, and the demo report is re-rendered from its own payload so the documentation shows what a run now produces.
Open the QC report on what the run found
On a 185-sample cohort 1,520 of the 1,579 gene-conversion calls were rows the samples themselves had left ambiguous, promoted because another sample called the same stretch outright. They carried the donor's bases in 2 to 5% of the reads and were fitted at the model's 20% floor, and the calls vouching for them came from samples at 0.2x to 2x depth whose tracts were read by one to four molecules. A tract most of whose sites are under the depth floor, or whose reads carry the donor's bases below half the thinnest tract the model fits (10%), is now ambiguous, with a reason that quotes the reads. The cohort pass applies the same check, so neither kind vouches for another sample or is promoted by one, and re-running only that step corrects the per-sample files of an earlier run. A corroborated row says what it fell short of on its own, the Bayes factor, the read fraction or both, instead of always quoting the Bayes factor, and the report says its "Carried by" column cannot read below 20%.
The ubiquity rule divided by every sample in the run, so on a cohort mapped against two references an artefact of one could never reach the 90% that marks a reference artifact: 113 of 185 samples is 61%. The fraction is now taken over the samples mapped to the same reference, which the per-locus files say, leaving out the samples diverged from it. Those carry the donor's base at every paralogous locus by inheritance, so they would join every event on their reference and inflate it. The pooled locus statistics were keyed on the pair number alone, which mixed pair 61 of one reference with pair 61 of the other; they are keyed on the contig as well. The share of loci with a call that marks a sample divergent is counted in distinct stretches of its own reference rather than in pairs, because a gene family reports one stretch through a pair per relative and one conversion seen through eleven relatives put a clean sample on the threshold.
The summary card and the page badge counted every called event alike. An event found only in samples the QC fails, mixed or contaminated cultures where the donor's bases have a simpler explanation, is now counted apart and those samples are named, and the headline and badge count the events of the samples the QC keeps. A WARN sample counts with them: a sample to review is not a sample to exclude. The breakpoint support is counted in events like the headline. Counted in calls, one event carried by two samples read as more backed events than there were events.
divergent_sample was keyed on the sample alone, so a sample mapped against two references and diverged from one lost its tracts on the other as well. The share of loci was taken over only the contigs the sample had tracts on: five tracts on a plasmid made a sample with a clean chromosome divergent. The rows of each per-sample file now carry the sample and reference the file was written for, the share is taken over every contig of that reference, and the cohort count leaves a sample out only on the reference it is diverged from. The cohort pass demoted a gene_conversion or ambiguous row resting on too few reads, but not a reciprocal_exchange one, which the per-sample caller does demote. It now does, and keeps the file's own reason, with its copy-number and breakpoint notes, when that reason already says as much. The summary listed, among the events only failed samples carry, a failed sample whose event a sample the QC keeps also carries.
5 of 6 tasks
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
What changed and why
Gene conversion was called from reads that could not carry it, and the cohort step multiplied those calls. On the 185-sample cohort it was checked against, 1,520 of the 1,579 calls were rows the samples themselves had left
ambiguous, promoted because another sample had called the same stretch outright. The outright calls came from cultures at 0.2x to 2x that the QC fails, read by one to four molecules. With this branch the calls drop from 38 events to 4, two of them in the mixed cultures the QC also fails.--gconv_min_depth, or whose reads carry the donor's bases below 10%, is nowambiguous, with a reason that quotes the reads, and it takes no part in corroboration. The cohort step applies the same check, reciprocal exchanges included, so re-running only that step corrects an existing run.reference_artifact(113 of 185 samples is 61%). The pooled locus statistics are keyed on the contig as well as the pair number.divergent_sampleis judged per sample and reference, over every contig of that reference. A sample diverged from one of two references lost its tracts on the other as well, and five tracts on a plasmid made a sample with a clean chromosome divergent.This also brings to
mainwhat the branch already holds:LINEAGE_MISMATCHand thedivergent_sampleverdict, two guards against a sample mapped to a genome it does not belong to;nanis no longer read as a gene name.Each has its
CHANGELOG.mdentry.Outputs and defaults
<samplesheet>_gene_conversion.tsv: many rows that readgene_conversionare nowambiguous, with areasonthat quotes the reads.qc_flags.tsvcan carryLINEAGE_MISMATCH(WARN).--kraken_memory(24 GB),--qosand--qos_heavy(short),--extra_binds. The container digest moves to the image with catalogue v1.0.2.Two PRs are stacked on this one: #7 (the SNP matrix's depth), then #8 (the cohort analyses).
How it was verified
Each commit was also tested on its own. CI passes on the branch head, every job green, the stub run included on Nextflow 24.04.2 and the latest stable release: run 35875974618.
Checklist
tests/run_tests.shpasses (unit, js and lint locally; the pipeline leg in CI)ruff check .is clean, and I did not restyle code the change does not touchmodules/has astub:block (not applicable: no new process)tests/data/was not edited by hand (fixtures come frommake_test_data.py)CHANGELOG.mdentry and, where relevant, a docs update