initial commit to support archaeal genome annotations

.gitignore

@@ -9,6 +9,9 @@ __pycache__/
 # C extensions
 *.so
 
+# Databases
+db/
+
 # Distribution / packaging
 .Python
 build/

Dockerfile

@@ -28,6 +28,7 @@ COPY . /tmp/source/
 
 RUN source /opt/conda/bashrc && micromamba activate \
     && python3 -m pip install --no-cache /tmp/source/ \
+    && /tmp/source/scripts/install-ori-finder-arch.sh /opt/conda/bin \
     && echo '#!/bin/bash' > /entrypoint.sh \
     && echo 'bakta "$@"' >> /entrypoint.sh \
     && chmod +x /entrypoint.sh \

README.md

@@ -12,9 +12,13 @@
 [![Galaxy Toolshed - Tool Version](https://img.shields.io/galaxytoolshed/v/bakta/iuc/bakta?label=usegalaxy.eu)](https://usegalaxy.eu/root?tool_id=bakta)
 [![Static Badge](https://img.shields.io/badge/bio.tools-v1.12.1-blue?link=https%3A%2F%2Fbio.tools%2Fbakta)](https://bio.tools/bakta)
 
-# Bakta: rapid & standardized annotation of bacterial genomes, MAGs & plasmids
+# Bakta: rapid & standardized annotation of prokaryotic genomes, MAGs & plasmids
 
-Bakta is a tool for the rapid & standardized annotation of bacterial genomes and plasmids from both isolates and MAGs. It provides **dbxref**-rich, **sORF**-including and taxon-independent annotations in machine-readable `JSON` & bioinformatics standard file formats for automated downstream analysis.
+Bakta is a tool for the rapid & standardized annotation of bacterial and archaeal genomes and plasmids from both isolates and MAGs. It provides **dbxref**-rich, **sORF**-including and taxon-independent annotations in machine-readable `JSON` & bioinformatics standard file formats for automated downstream analysis.
+
+Archaeal annotation is selected with `--domain archaea`. The shared Bakta database contains both bacterial and archaeal profile assets; this option selects the archaeal tRNA/rRNA/ncRNA covariance models, archaeal/domain-neutral AntiFam and expert-protein references, and Ori-Finder-Arch plus MEME/FIMO for oriC prediction. AMRFinderPlus and oriT prediction are intentionally not run in this profile. Install the standalone Ori-Finder-Arch binary with `scripts/install-ori-finder-arch.sh` after activating the environment; the script downloads the upstream v1.0.0 release into `$CONDA_PREFIX/bin` and verifies its SHA-256 checksum.
+
+The remaining database curation and validation work is tracked in [the archaeal support plan](docs/archaeal-support-plan.md). Every database build requires curated inputs through `BAKTA_ARCHAEAL_EXPERT_PROTEINS` and `BAKTA_ARCHAEAL_ANTIFAM_HMM`; these prevent bacterial expert/AntiFam data from being used implicitly.
 
 ## Contents
 
@@ -63,8 +67,8 @@ Bakta annotates ncRNA cis-regulatory regions, oriC/oriV/oriT and assembly gaps a
 - **GFF3 & INSDC conform annotations**
 Bakta writes GFF3 and INSDC-compliant (Genbank & EMBL) annotation files ready for submission (checked via [GenomeTools GFF3Validator](http://genometools.org/cgi-bin/gff3validator.cgi), [table2asn_GFF](https://www.ncbi.nlm.nih.gov/genbank/genomes_gff/#run) and [ENA Webin-CLI](https://github.com/enasequence/webin-cli) for GFF3 and EMBL file formats, respectively for representative genomes of all ESKAPE species).
 
-- **Bacteria & plasmids only**
-Bakta was designed to annotate bacteria (isolates & MAGs) and plasmids, only. This decision by design has been made in order to tweak the annotation process regarding tools, preferences & databases and to streamline further development & maintenance of the software.
+- **Bacterial and archaeal profiles**
+The default bacterial profile supports isolates, MAGs, and plasmids. The archaeal profile is selected with `--domain archaea`; it uses domain-specific RNA models, protein references, recoding rules, and origin prediction while preserving the same standard output formats.
 
 - **Reasoning**
 By annotating bacterial genomes in a standardized, taxonomy-independent, high-throughput and local manner, Bakta aims at a well-balanced tradeoff between fully featured but computationally demanding pipelines like [PGAP](https://github.com/ncbi/pgap) and rapid highly customizable offline tools like [Prokka](https://github.com/tseemann/prokka). Indeed, Bakta is heavily inspired by Prokka (kudos to [Torsten Seemann](https://github.com/tseemann)) and many command line options are compatible for the sake of interoperability and user convenience. Hence, if Bakta does not fit your needs, please consider trying Prokka.
@@ -377,7 +381,7 @@ Exemplary annotation result files for several genomes (mostly ESKAPE species) ar
 
 ```bash
 usage: bakta [--db DB] [--min-contig-length MIN_CONTIG_LENGTH] [--prefix PREFIX] [--output OUTPUT] [--force]
-             [--genus GENUS] [--species SPECIES] [--strain STRAIN] [--plasmid PLASMID]
+             [--domain {bacteria,archaea}] [--genus GENUS] [--species SPECIES] [--strain STRAIN] [--plasmid PLASMID]
              [--complete] [--prodigal-tf PRODIGAL_TF] [--translation-table {11,4,25}] [--gram {+,-,?}]
              [--locus LOCUS] [--locus-tag LOCUS_TAG] [--locus-tag-increment {1,5,10}] [--keep-contig-headers] [--compliant]
              [--replicons REPLICONS] [--regions REGIONS] [--proteins PROTEINS] [--hmms HMMS] [--meta]
@@ -401,7 +405,9 @@ Input / Output:
                         Output directory (default = current working directory)
   --force, -f           Force overwriting existing output folder (except for current working directory)
 
-Organism:
+Domain:
+  --domain {bacteria,archaea}
+                        Annotation domain: bacteria/archaea (default = bacteria)
   --genus GENUS         Genus name
   --species SPECIES     Species name
   --strain STRAIN       Strain name
@@ -413,7 +419,7 @@ Annotation:
                         Path to existing Prodigal training file to use for CDS prediction
   --translation-table {11,4,25}
                         Translation table: 11/4/25 (default = 11)
-  --gram {+,-,?}        Gram type for signal peptide predictions: +/-/? (default = ?)
+  --gram {+,-,?}        Gram type for signal peptide predictions: +/-/? (default = ?; ignored for --domain archaea)
   --locus LOCUS         Locus prefix (default = 'contig')
   --locus-tag LOCUS_TAG
                         Locus tag prefix (default = autogenerated)

bakta/config.py

@@ -43,6 +43,7 @@
 strain = None
 plasmid = None
 taxon = None
+domain = None
 
 # annotation configuration
 complete = None
@@ -76,6 +77,11 @@
 skip_filter = None
 skip_plot = None
 
+
+def normalize_gram(domain: str, gram_value: str) -> str:
+    """Return the applicable Gram value for a domain profile."""
+    return bc.GRAM_UNKNOWN if domain == bc.DOMAIN_ARCHAEA else gram_value
+
 run_start = datetime.now()
 run_end = None
 
@@ -121,8 +127,10 @@ def setup(args):
     force = args.force
     log.info('force=%s', force)
 
-    # organism configurations
-    global genus, species, strain, plasmid, taxon
+    # domain configurations
+    global genus, species, strain, plasmid, taxon, domain
+    domain = args.domain
+    log.info('domain=%s', domain)
     genus = args.genus
     if(genus is not None):
         genus = genus.strip()
@@ -184,7 +192,9 @@ def setup(args):
     log.info('prodigal_tf=%s', prodigal_tf)
     translation_table = args.translation_table
     log.info('translation_table=%s', translation_table)
-    gram = args.gram
+    if(domain == bc.DOMAIN_ARCHAEA and args.gram != bc.GRAM_UNKNOWN):
+        log.info('ignore gram=%s for archaeal profile', args.gram)
+    gram = normalize_gram(domain, args.gram)
     log.info('gram=%s', gram)
     compliant = args.compliant
     log.info('compliant=%s', compliant)
@@ -437,4 +447,4 @@ def check_tmp_path(args: Namespace) -> Path:
     else:
         tmp_path = Path(tempfile.mkdtemp()).resolve()
     log.info('tmp-path=%s', tmp_path)
-    return tmp_path
+    return tmp_path

bakta/constants.py

@@ -193,6 +193,27 @@
 GRAM_NEGATIVE = '-'
 GRAM_UNKNOWN = '?'
 
+# Organism profiles
+DOMAIN_BACTERIA = 'bacteria'
+DOMAIN_ARCHAEA = 'archaea'
+
+# Domain-specific gene-calling and recoding policy. These are intentionally
+# explicit so that changes to one profile cannot alter the other.
+DOMAIN_PROFILES = {
+    DOMAIN_BACTERIA: {
+        'sorf_start_codons': ('ATG',),
+        'tmrna_gene': 'ssrA',
+        'tmrna_product': 'transfer-messenger RNA, SsrA',
+        'recoding_codons': {'TGA': 'selenocysteine'}
+    },
+    DOMAIN_ARCHAEA: {
+        'sorf_start_codons': ('ATG', 'GTG', 'TTG'),
+        'tmrna_gene': None,
+        'tmrna_product': 'transfer-messenger RNA',
+        'recoding_codons': {'TGA': 'selenocysteine', 'TAG': 'pyrrolysine'}
+    }
+}
+
 
 ############################################################################
 # Replicon types, length thresholds & topology

bakta/db.py

@@ -50,12 +50,31 @@
         'rRNA.i1m',
         'rRNA.i1p',
         'sorf.dmnd'
+]
+ARCHAEAL_FILE_NAMES = [
+        'antifam-archaea.h3f',
+        'antifam-archaea.h3i',
+        'antifam-archaea.h3m',
+        'antifam-archaea.h3p',
+        'ncRNA-genes-archaea.i1f',
+        'ncRNA-genes-archaea.i1i',
+        'ncRNA-genes-archaea.i1m',
+        'ncRNA-genes-archaea.i1p',
+        'ncRNA-regions-archaea.i1f',
+        'ncRNA-regions-archaea.i1i',
+        'ncRNA-regions-archaea.i1m',
+        'ncRNA-regions-archaea.i1p',
+        'rRNA-archaea.i1f',
+        'rRNA-archaea.i1i',
+        'rRNA-archaea.i1m',
+        'rRNA-archaea.i1p',
+        'expert-protein-sequences-archaea.dmnd'
     ]
 FILE_PERMISSIONS = stat.S_IRUSR | stat.S_IWUSR | stat.S_IRGRP | stat.S_IWGRP | stat.S_IROTH
 DIR_PERMISSIONS = stat.S_IRUSR | stat.S_IWUSR | stat.S_IXUSR | stat.S_IRGRP | stat.S_IWGRP | stat.S_IXGRP | stat.S_IROTH | stat.S_IXOTH
 
 
-def check(db_path: Path) -> dict:
+def check(db_path: Path, domain: str = bc.DOMAIN_BACTERIA) -> dict:
     """Check if database directory exists, is accessible and contains necessary files."""
 
     if(db_path is None):
@@ -91,7 +110,20 @@ def check(db_path: Path) -> dict:
         log.error('wrong database version detected! required=%i, detected=%i', bakta.__db_schema_version__, db_info['major'])
         sys.exit(f"ERROR: wrong database version detected!\nBakta version {cfg.version} requires database version {bakta.__db_schema_version__}.x, but {db_info['major']}.{db_info['minor']} was detected. Please, update Bakta or download a compatible database version from https://doi.org/10.5281/zenodo.4247253")
 
-    required_db_files = FILE_NAMES
+    # New database builds contain both profiles. Require only the selected
+    # profile here so existing bacterial database releases remain usable.
+    profile_files = {
+        'ncRNA-genes.i1f', 'ncRNA-genes.i1i', 'ncRNA-genes.i1m', 'ncRNA-genes.i1p',
+        'ncRNA-regions.i1f', 'ncRNA-regions.i1i', 'ncRNA-regions.i1m', 'ncRNA-regions.i1p',
+        'rRNA.i1f', 'rRNA.i1i', 'rRNA.i1m', 'rRNA.i1p'
+    }
+    required_db_files = [file_name for file_name in FILE_NAMES if file_name not in profile_files and file_name != 'expert-protein-sequences.dmnd']
+    if domain == bc.DOMAIN_ARCHAEA:
+        required_db_files = [file_name for file_name in required_db_files if file_name not in {'oric.fna', 'orit.fna'}]
+        required_db_files.extend(ARCHAEAL_FILE_NAMES)
+    else:
+        required_db_files.extend(profile_files)
+        required_db_files.append('expert-protein-sequences.dmnd')
     required_db_files.append('psc.dmnd' if db_info['type'] == 'full' else 'pscc.dmnd')
     for file_name in required_db_files:
         path = db_path.joinpath(file_name)

bakta/features/cds.py

@@ -489,16 +489,21 @@ def revise_translational_exceptions(data: dict, cdss: Sequence[dict]):
             if(
                 cds_a['strand'] == cds_b['strand'] and  # up- and downstream ORFs on the same strand
                 cds_a['frame'] == cds_b['frame'] and  # up- and downstream ORFs on the same frame
-                upstream_stop_codon == 'TGA' and  # tRNAScan-SE 2.0 only predicts tRNA-Sec with UCA anticodons, therefore we can only detect TGA stop codons
+                upstream_stop_codon in bc.DOMAIN_PROFILES[cfg.domain or bc.DOMAIN_BACTERIA]['recoding_codons'] and
                 (cds_b['start'] - cds_a['stop']) < 100):  # up- and downstream ORFs in close proximity
                 cds_pairs = cds_pairs_per_sequence[cds_a['sequence']]
-                cds_pairs.append((cds_a, cds_b))
+                cds_pairs.append((cds_a, cds_b, upstream_stop_codon))
 
     recoding_regions = [ncrna_region for ncrna_region in data['features'] if ncrna_region['type'] == bc.FEATURE_NC_RNA_REGION  and  ncrna_region['class'] == so.SO_CIS_REG_RECODING_STIMULATION_REGION]  #  Selenocysteine insertion sequences
     for recoding_region in recoding_regions:
-        if('selenocysteine' in recoding_region.get('product', '').lower()):
+        product = recoding_region.get('product', '').lower()
+        recoding_profile = bc.DOMAIN_PROFILES[cfg.domain or bc.DOMAIN_BACTERIA]
+        recoding_type = next((value for codon, value in recoding_profile['recoding_codons'].items() if value in product), None)
+        if(recoding_type is not None):
             cds_pairs = cds_pairs_per_sequence[recoding_region['sequence']]
-            for cds_a, cds_b in cds_pairs:  # find CDS pair around recoding region
+            for cds_a, cds_b, stop_codon in cds_pairs:  # find CDS pair around recoding region
+                if(recoding_profile['recoding_codons'][stop_codon] != recoding_type):
+                    continue
                 strand = cds_a['strand']
                 if(
                     strand == recoding_region['strand'] and  # everything is on the same strand
@@ -519,16 +524,16 @@ def revise_translational_exceptions(data: dict, cdss: Sequence[dict]):
                         aa[-1] == '*' and  # ends with stop *
                         aa[1:-1].count('*') == 1  # contains exactly 1 additional stop (*) somewhere in between
                         ):
-                        aa = aa.replace('*', 'U', 1)  # replace internal stop codon by U -> selenocysteine
+                        aa = aa.replace('*', 'U' if recoding_type == 'selenocysteine' else 'O', 1)
                         aa = aa[:-1]  # remove stop asterisk
                         seleno_cds['aa'] = aa
                         seleno_cds['aa_digest'], seleno_cds['aa_hexdigest'] = bu.calc_aa_hash(aa)
                         seleno_cds['exception'] = {
-                            'type': 'selenocysteine',
-                            'aa': 'Sec',
+                            'type': recoding_type,
+                            'aa': 'Sec' if recoding_type == 'selenocysteine' else 'Pyl',
                             'start': cds_a['stop'] - 2 if strand == bc.STRAND_FORWARD else cds_b['start'],
                             'stop': cds_a['stop'] if strand == bc.STRAND_FORWARD else cds_b['start'] + 2,
-                            'codon_position': aa.find('U') + 1
+                            'codon_position': aa.find('U' if recoding_type == 'selenocysteine' else 'O') + 1
                         }
                         cdss.append(seleno_cds)
                         log.info(

bakta/features/nc_rna.py

@@ -34,7 +34,8 @@ def predict_nc_rnas(data: dict, sequences_path: Path):
     if(data['stats']['size'] >= 1000000):
         cmd.append('-Z')
         cmd.append(str(2 * data['stats']['size'] // 1000000))
-    cmd.append(str(cfg.db_path.joinpath('ncRNA-genes')))
+    database = 'ncRNA-genes-archaea' if cfg.domain == bc.DOMAIN_ARCHAEA else 'ncRNA-genes'
+    cmd.append(str(cfg.db_path.joinpath(database)))
     cmd.append(str(sequences_path))
     log.debug('cmd=%s', cmd)
     proc = sp.run(

bakta/features/nc_rna_region.py

@@ -33,7 +33,8 @@ def predict_nc_rna_regions(data: dict, sequences_path: Path):
     if(data['stats']['size'] >= 1000000):
         cmd.append('-Z')
         cmd.append(str(2 * data['stats']['size'] // 1000000))
-    cmd.append(str(cfg.db_path.joinpath('ncRNA-regions')))
+    database = 'ncRNA-regions-archaea' if cfg.domain == bc.DOMAIN_ARCHAEA else 'ncRNA-regions'
+    cmd.append(str(cfg.db_path.joinpath(database)))
     cmd.append(str(sequences_path))
     log.debug('cmd=%s', cmd)
     proc = sp.run(

bakta/features/orf.py

@@ -20,7 +20,8 @@ def detect_spurious(orfs: Sequence[dict]):
     orf_by_aa_digest = get_orf_dictionary(orfs)
     alphabet: "AA" = pyhmmer.easel.Alphabet.amino()
     proteins: "DigitalSequenceBlock[AA]" = TextSequenceBlock(TextSequence(sequence=orf['aa'], name=get_orf_key(orf)) for orf in orfs).digitize(alphabet)
-    with pyhmmer.plan7.HMMFile(cfg.db_path.joinpath('antifam'), alphabet=alphabet) as hmm_fh:
+    antifam_name = 'antifam-archaea' if cfg.domain == bc.DOMAIN_ARCHAEA else 'antifam'
+    with pyhmmer.plan7.HMMFile(cfg.db_path.joinpath(antifam_name), alphabet=alphabet) as hmm_fh:
         for top_hits in pyhmmer.hmmsearch(hmm_fh, proteins, bit_cutoffs='gathering', cpus=cfg.threads):
             for hit in top_hits:
                 orf = orf_by_aa_digest[hit.name]