Rfam ncrna - #441
Open
gbouras13 wants to merge 10 commits into
Open
Conversation
pharokka had no covariance-model search, so structured RNAs (riboswitches,
ribozymes, regulatory sRNAs, introns) were invisible. `--rfam` scans the
genome against Rfam 15.1 with Infernal cmscan.
Opt-in, because it roughly doubles runtime on a small phage genome (~30-60s
measured). Flags and defaults are driven by measurement, not guesswork:
- `--rfam` filter preset is mandatory: Infernal derives filter strictness
from database size, so on a phage-sized target the defaults go permissive
and a 58 kb genome takes >25 min instead of 32 s.
- `--fmt 2` + `--clanin` enables clan competition; parsing drops `olp == '='`
rows, without which riboswitch clans emit piles of redundant overlapping
calls for one locus.
- Rfam tRNA/tmRNA hits are dropped by default. RF00023 does not hit the tmRNA
in NC_051700 that ARAGORN finds, so Rfam is additive and must never displace
the specialised tools. Locked in by a regression test.
- `--cpu` is passed through but does not help a single contig: Infernal
parallelises over the sequence database, so 8 cpus measured slower than 1 on
one 42 kb genome. It does help in --meta mode.
New outputs: `{prefix}_ncrna.tsv`, `{prefix}_cmscan.tblout`, `ncRNA` features
in the GFF/GBK, and an `ncRNAs` count per contig in `_cds_functions.tsv`.
Also fixed two latent database bugs found along the way: the tarball filename
and the PHROGs version marker were both hardcoded to 1.8.0 rather than derived
from the database version, so any version bump would have silently used stale
names.
The v1.11.0 database URL and md5 are placeholders pending the Zenodo upload.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Infernal parallelises over the sequence database rather than the model database, so with a handful of contigs there is nothing to divide up and the coordination overhead makes it slower. Measured on one 42 kb phage contig: 8 cpus took 59.4 s versus 36.8 s at 1 cpu. With 100 contigs the same scan went 238.1 s -> 132.7 s, so threads are worth using once there are enough sequences. Clamped silently in cmscan_threads(). The user asked for N threads for pharokka as a whole and every other step still uses them, so there is nothing here they need to act on. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This reverts commit 8fdea31.
The claim that threads do not help a single contig was wrong, and the clamp
reverted in the previous commit was based on it.
The original measurement used `cmsearch` against an unpressed flatfile, where
the CMs are the query and parallelism comes from the sequence database - so one
contig meant no parallelism. pharokka runs `cmscan` against the cmpress'd
database, where cmscan divides the *model* database across threads. There are
4227 models, so there is always work to parallelise regardless of contig count.
Measured on 8 cores (M1 Pro, 6P+2E), full Rfam 15.1, 1 thread -> 8 threads:
NC_043029 7.6 kb 3.0 s -> 0.7 s 4.3x
NC_004617 42.7 kb 26.4 s -> 4.9 s 5.3x
NC_051700 58.8 kb 20.0 s -> 4.8 s 4.2x
SAOMS1 140.0 kb 75.3 s -> 13.5 s 5.6x
100 contigs 587.0 kb 253.1 s -> 63.2 s 4.0x
Scaling is essentially flat across contig counts from 1 to 100 (~4-5x at 8
threads throughout), so there is no contig-count threshold worth special
casing. The clamp would have made the commonest case - one phage genome -
about 5x slower.
Runtime figures in the docs were also measured with the slow method and are
corrected: a typical 40 kb phage is ~5 s on 8 threads, not 30-60 s.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The measured cost no longer justifies opt-in: ~5 s for a typical 40 kb phage and ~14 s for a 140 kb phage on 8 threads, against a pharokka run that already takes minutes. `--rfam` is replaced by `--skip_rfam`. Meta mode is the exception and defaults to skipping, because cmscan runtime scales with assembly size (~2 min/Mbp on 8 threads, so ~3 h for a 100 Mbp metagenome). `--meta_rfam` forces it on, mirroring how --meta disables PyHMMER unless --meta_hmm is given. Because ncRNA annotation now runs by default, v1.11.0 requires the v1.11.0 database and Infernal; running against an older database fails with an explanatory error rather than silently skipping. `--skip_rfam` restores the previous behaviour exactly. Golden cases all pass --skip_rfam so they stay byte-identical and need no regeneration - requiring the ~800 MB Rfam database in the golden suite would be disproportionate when tests/test_rfam.py already covers the path end to end against a committed five-model subset. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
cmscan reads the cmpress'd Rfam.cm.i1{f,i,m,p} files and only uses the
Rfam.cm path as a base name, so the 329 MB flatfile is not needed at runtime.
Verified by running cmscan against a directory containing only the pressed
files: same three hits on NC_004617. Dropping it takes the installed database
from 2.4 GB to 2.1 GB.
Also adds coverage for check_rfam_installation, which is load-bearing now that
ncRNA annotation is on by default. Note it never returns False in practice --
pharokka's logger.error sink exits first -- so the tests assert SystemExit.
md5 recorded for the built v1.11.0 tarball (735 MB). db_url still needs the
real Zenodo record.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…t it
An end-to-end run showed ncRNAs appearing in the gff, gbk and _ncrna.tsv but
missing entirely from pharokka.tbl - create_tbl() handled CDS, tRNA, CRISPR and
tmRNA only. The full 178-test suite passed with this broken, because:
- all six golden cases pass --skip_rfam, so no golden output contained an
ncRNA at all
- test_overall's happy-path cases have no assertions; exec_command only
raises if pharokka exits non-zero, so they verify "did not crash"
So the expensive tests were paying the whole Rfam cost while checking none of
its output. Adds a "rfam" golden case on NC_004617 (three Rfam sRNA families:
SprD, SprX, rli28) - the one case that does not skip Rfam - plus
pharokka_ncrna.tsv in KEY_FILES. It now locks the ncRNA rows in the gff, gbk,
tbl and _cds_functions.tsv. The raw _cmscan.tblout is excluded as it embeds the
command line and run date.
.tbl specifics: ncRNA_class is mandatory on the ncRNA feature key, and the tbl
qualifier is db_xref rather than the GFF3 spelling Dbxref.
Also writes _ncrna.tsv with quote_style="never" so a family with no clan gives
an empty field rather than a literal "".
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Previously the row only appeared when Rfam annotation ran, so downstream parsers had to cope with a schema that changed depending on flags. It is now written unconditionally, one row per contig, with a count of 0 when Rfam was skipped. Golden outputs regenerated. The only substantive change is one added ncRNAs line per contig; the .gbk diffs are the LOCUS date line, which the comparator already ignores. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #441 +/- ##
==========================================
+ Coverage 87.88% 88.16% +0.27%
==========================================
Files 17 18 +1
Lines 2658 2839 +181
==========================================
+ Hits 2336 2503 +167
- Misses 322 336 +14 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
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.
No description provided.