Repository navigation
Search every genome with the same proteins for the pangenome - #69
Merged
Merged
Conversation
PANGENE aligned each sample's own proteins to the other samples' genomes and the report merged the matrices on gene names that carry the sample, so no gene was ever in two samples and nothing was core. It now aligns the reference proteome and each sample's genes that no reference gene matches to every genome, the sample's own included, runs pangene once and names each column after its sample. On the first 2 Mb of CDC1551, M. bovis and H37Rv the report goes from 5379 families in 6 columns, none core, to 1810 genes in 3 samples, 1760 of them core.
GitHub's ubuntu-latest became Ubuntu 26.04, for which setup-python has no Python 3.9, so the job stopped at setup before running a test. Python 3.9 is the container's, so its job now runs on 24.04 and 3.11 stays on the latest.
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
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.
The bug
PANGENE aligned each sample's own corrected proteins to every other sample's genome and ran pangene once per sample, so each Rtab was (this sample's genes) x (the other genomes), its columns named
<genome>_<sample>#0. REPORT_INPUTS then outer-joined the Rtabs on gene names that carry the sample's name and coordinates:A4_Mbovis_Rv0001;+;0-1524andL4_CDC1551_Rv0001;+;0-1524were two rows. With N samples the matrix had N(N-1) columns, a gene could be in at most N-1 of them, and no gene was ever core. A sample's own genome was never searched either, so a gene only one genome has never entered the matrix: the interactive report gave every genome 0 genes "in no other genome".The fix
pangene is meant to be given one set of proteins aligned to every genome. PANGENE is now one task for the run:
bin/pangenome.py proteinswrites the set: the reference proteome (H37Rv or a custom one), each protein named<gene>:<protein>so that pangene calls the gene by its locus (Rv1506c,b0002), and each sample's genes that no reference gene matches (an emptyrvsin<sample>.AMAP.tsv), named after the sample's gene. A sequence already in the set is not added again.gfa2matrix,callandgfatools stat.bin/pangenome.py matrixnames each column after its sample (L4_CDC1551, no#0) and gives a genome none of the proteins was found in a column of zeros, which pangene leaves out.Of the four designs measured on real data before coding (reference genes only, every sample's proteins pooled, reference plus the genes outside it, per-sample runs merged by reference name), this is the third: rows a reader can look up, genes outside the reference kept and matched across samples by pangene, no new tool.
PANGENE/pangenome.{gfa,bubbles.txt,Rtab,stat}, one set for the run instead ofPANGENE/{sample}/..statwas never published before: no output declared it.merge_rtabs.py, which leaves a single matrix as it is (byte for byte in every run below).final_report.htmlreads it unchanged; the pangenome card ofamap_report.htmlnow reads a column per genome.--skip_typingit still gives each genome of the sheet its own<sampleId>.fasta.tools/build_docs.py.Before and after, real runs
Typed runs under
-profile dockerwith Nextflow 26.04.6, eggNOG-mapper replaced by a stand-in (its database is about 45 GB, and nothing here reads it); main is 6351c0c.A4_Mbovis,L4_CDC1551: 1768 core, 17 only in M. bovis, 24 only in CDC1551--skip_typing, both files calledassembly.fnaSakai_K12#0,K12_Sakai#0), none in bothK12,Sakai: 296 in both, 63 only in K-12, 106 only in SakaiThe 50 variable genes of the three-genome run, by where they are missing:
A4_Mbovis_01833,A4_Mbovis_01834,L4_CDC1551_01839),L4_CDC1551_00745, and Rv1767, Rv1768, Rv3724A past the end of the H37Rv windowL4_CDC1551_01396(an Rv1320c-like protein) in CDC1551The pangenome card of
amap_report.html, same run: on main every genome had 0 genes in no other genome and its own count in every other one (1761, 1762, 1764); now 1760 for all three, and M. bovis has 2 genes of its own and CDC1551 1.PANGENE took 3.6 s and 382 MB for the three 2 Mb genomes. A reference gene can still be in no column: in the E. coli run, 15 of K-12's 374, six of them proteins of 10 to 24 residues, which miniprot does not place.
Tests
753 tests, up from 744.
tests/test_pangene_genomes.pyis rewritten. It renders PANGENE_GENOME's and PANGENE's scripts as Nextflow does and runs them ontests/data/pangenome/: real genomes of CDC1551, M. bovis and H37Rv cut around RD4 and RvD2 out of the three-genome run above, rebuilt byte for byte by itsmake_fixture.py. Stand-ins for miniprot, pangene, k8 and gfatools hand names on as the real ones do. It checks that:/bin/bashcalls an empty array unbound underset -u);bin/pangenome.pybehaves on its own.Where the real tools are installed, as in the container, RD4 must come out missing from M. bovis only and the three RvD2 genes in CDC1551 and M. bovis only, one row each.
tests/test_report_tables.py,tests/test_amap_report.pyandtests/test_sample_ids.pygive the reports the matrix in its new shape. Where R is installed, the report must count the genes all the samples share.Each part broken on purpose fails 2 to 5 tests: no guard on the last wait,
#0kept, no column of zeros, a shared sequence added twice, the sample's own genome left out.Ran:
nextflow lint(26.10.2): no errors. The only new warnings are theprojectDirones every module has, andmain.nfhas two fewer.-previewis clean on 24.04.4 and 26.10.2, typed and with--skip_typing.tools/build_docs.py --checkandtools/check_site.pypass.CI
The second commit runs the Python 3.9 job on
ubuntu-24.04. GitHub'subuntu-latestbecame Ubuntu 26.04, for which setup-python has no Python 3.9, so on this PR the job failed at setup before running a test, and it would on any other: main's last run, on 20 October, was still on 24.04. 3.9 is the container's Python; 3.11 stays onubuntu-latest.