diff --git a/.gitignore b/.gitignore index eef27aea..e8db3752 100644 --- a/.gitignore +++ b/.gitignore @@ -9,6 +9,9 @@ __pycache__/ # C extensions *.so +# Databases +db/ + # Distribution / packaging .Python build/ diff --git a/Dockerfile b/Dockerfile index 888974a5..a1856c51 100644 --- a/Dockerfile +++ b/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 \ diff --git a/README.md b/README.md index 95582a26..1f3bf7ac 100644 --- a/README.md +++ b/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) diff --git a/bakta/config.py b/bakta/config.py index bfa1f9a5..bd01b959 100644 --- a/bakta/config.py +++ b/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 \ No newline at end of file + return tmp_path diff --git a/bakta/constants.py b/bakta/constants.py index 1796ae78..7cf3ac9f 100644 --- a/bakta/constants.py +++ b/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 diff --git a/bakta/db.py b/bakta/db.py index 5b4a2a41..0cf3e4ba 100644 --- a/bakta/db.py +++ b/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) diff --git a/bakta/features/cds.py b/bakta/features/cds.py index b1e034a6..1b4de6b1 100644 --- a/bakta/features/cds.py +++ b/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( diff --git a/bakta/features/nc_rna.py b/bakta/features/nc_rna.py index 16b94404..37066988 100644 --- a/bakta/features/nc_rna.py +++ b/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( diff --git a/bakta/features/nc_rna_region.py b/bakta/features/nc_rna_region.py index 3dfe88b9..78d40253 100644 --- a/bakta/features/nc_rna_region.py +++ b/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( diff --git a/bakta/features/orf.py b/bakta/features/orf.py index 251fd8b5..5ae7d555 100644 --- a/bakta/features/orf.py +++ b/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] diff --git a/bakta/features/ori_arch.py b/bakta/features/ori_arch.py new file mode 100644 index 00000000..65eea2fd --- /dev/null +++ b/bakta/features/ori_arch.py @@ -0,0 +1,83 @@ +import logging +import subprocess as sp + +from collections import OrderedDict +from pathlib import Path +from typing import Sequence + +import bakta.config as cfg +import bakta.constants as bc +import bakta.io.fasta as fasta +import bakta.utils as bu + + +log = logging.getLogger('ORI_ARCH') + + +def predict_orics(data: dict, sequences_path: Path) -> Sequence[dict]: + """Predict archaeal replication origins with Ori-Finder-Arch. + + Ori-Finder-Arch emits GFF3 records with type ``rep_origin``. Its own RIP + CDS records are deliberately not imported: Bakta's CDS caller remains the + sole source of CDS features. + """ + sequences = {seq['id']: seq for seq in data['sequences']} + orics = [] + # Ori-Finder-Arch accepts multi-record FASTA, but topology and assembly + # level are global options. Invoke it per replicon so mixed assemblies are + # never silently downgraded to a single linear/draft interpretation. + for index, sequence in enumerate(data['sequences'], start=1): + input_path = cfg.tmp_path.joinpath(f'ori-arch.{index}.fna') + output_path = cfg.tmp_path.joinpath(f'ori-arch.{index}.gff3') + fasta.export_sequences([sequence], input_path) + topology = sequence['topology'] + level = 'complete' if sequence.get('complete', data['genome']['complete']) else 'draft' + _run_ori_finder(input_path, output_path, topology, level) + orics.extend(_parse_orics(output_path, sequences)) + log.info('predicted=%i', len(orics)) + return orics + + +def _run_ori_finder(input_path: Path, output_path: Path, topology: str, level: str): + cmd = ['OriFinderArch', '-i', str(input_path), '-o', str(output_path), '-t', topology, '-l', level] + log.debug('cmd=%s', cmd) + proc = sp.run(cmd, cwd=str(cfg.tmp_path), env=cfg.env, stdout=sp.PIPE, stderr=sp.PIPE, universal_newlines=True) + if(proc.returncode != 0): + log.debug('stdout=%s, stderr=%s', proc.stdout, proc.stderr) + raise Exception(f'Ori-Finder-Arch error! exit_code={proc.returncode}; stderr={proc.stderr.strip()}') + + +def _parse_orics(output_path: Path, sequences: dict) -> Sequence[dict]: + if(not output_path.is_file()): + raise Exception('Ori-Finder-Arch completed without producing a GFF3 result') + orics = [] + with output_path.open() as fh: + for line in fh: + if(line.startswith('#')): + continue + fields = line.rstrip('\n').split('\t') + if(len(fields) != 9): + log.warning('skip malformed Ori-Finder-Arch GFF3 line: %s', line.rstrip()) + continue + sequence_id, source, feature_type, start, stop, score, strand, phase, attributes = fields + if(feature_type != 'rep_origin' or sequence_id not in sequences): + continue + ori = OrderedDict() + ori['type'] = bc.FEATURE_ORIC + ori['sequence'] = sequence_id + ori['start'] = int(start) + ori['stop'] = int(stop) + ori['strand'] = bc.STRAND_UNKNOWN + ori['product'] = 'origin of replication' + ori['inference'] = 'ab initio prediction:Ori-Finder-Arch:1.0.0' + ori['attributes'] = _parse_attributes(attributes) + if(score != '.'): + ori['score'] = float(score) + ori['nt'] = bu.extract_feature_sequence(ori, sequences[sequence_id]) + orics.append(ori) + log.info('oriC: seq=%s, start=%i, stop=%i', sequence_id, ori['start'], ori['stop']) + return orics + + +def _parse_attributes(raw: str) -> dict: + return {key: value.strip("'") for key, value in (field.split('=', 1) for field in raw.split(';') if '=' in field)} diff --git a/bakta/features/r_rna.py b/bakta/features/r_rna.py index b5f5cb92..aea227ce 100644 --- a/bakta/features/r_rna.py +++ b/bakta/features/r_rna.py @@ -16,11 +16,31 @@ log = logging.getLogger('R_RNA') +R_RNA_PROFILES = { + bc.DOMAIN_BACTERIA: { + 'database': 'rRNA', + 'models': { + 'RF00001': ('5S', 'rrf', 119, 'K01985', so.SO_RRNA_5S.id), + 'RF00177': ('16S', 'rrs', 1533, 'K01977', so.SO_RRNA_16S.id), + 'RF02541': ('23S', 'rrl', 2925, 'K01980', so.SO_RRNA_23S.id) + } + }, + bc.DOMAIN_ARCHAEA: { + 'database': 'rRNA-archaea', + 'models': { + 'RF00001': ('5S', 'rrf', 119, 'K01985', so.SO_RRNA_5S.id), + 'RF01959': ('16S', 'rrs', 1477, 'K01977', so.SO_RRNA_16S.id), + 'RF02540': ('23S', 'rrl', 2990, 'K01980', so.SO_RRNA_23S.id) + } + } +} + def predict_r_rnas(data: dict, sequences_path: Path): """Search for ribosomal RNA sequences.""" output_path = cfg.tmp_path.joinpath('rrna.tsv') + profile = R_RNA_PROFILES[cfg.domain or bc.DOMAIN_BACTERIA] cmd = [ 'cmscan', '--noali', @@ -34,7 +54,7 @@ def predict_r_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('rRNA'))) + cmd.append(str(cfg.db_path.joinpath(profile['database']))) cmd.append(str(sequences_path)) log.debug('cmd=%s', cmd) proc = sp.run( @@ -74,25 +94,15 @@ def predict_r_rnas(data: dict, sequences_path: Path): else: truncated = None - db_xrefs = [f'{bc.DB_XREF_GO}:0005840', f'{bc.DB_XREF_GO}:0003735'] - if(accession == 'RF00001'): - rrna_tag = '5S' - db_xrefs += [f'{bc.DB_XREF_RFAM}:RF00001', f'{bc.DB_XREF_KOFAM}:K01985', so.SO_RRNA_5S.id] - consensus_length = 119 - elif(accession == 'RF00177'): - rrna_tag = '16S' - db_xrefs += [f'{bc.DB_XREF_RFAM}:RF00177', f'{bc.DB_XREF_KOFAM}:K01977', so.SO_RRNA_16S.id] - consensus_length = 1533 - elif(accession == 'RF02541'): - rrna_tag = '23S' - db_xrefs += [f'{bc.DB_XREF_RFAM}:RF02541', f'{bc.DB_XREF_KOFAM}:K01980', so.SO_RRNA_23S.id] - consensus_length = 2925 - else: + model = profile['models'].get(accession) + if(model is None): log.warning( 'unknown rRNA detected! accession=%s, seq=%s, start=%i, stop=%i, strand=%s, length=%i, truncated=%s, score=%1.1f, evalue=%1.1e', accession, sequence_id, start, stop, strand, length, truncated, score, evalue ) continue + rrna_tag, rrna_gene, consensus_length, kofam_id, so_id = model + db_xrefs = [f'{bc.DB_XREF_GO}:0005840', f'{bc.DB_XREF_GO}:0003735', f'{bc.DB_XREF_RFAM}:{accession}', f'{bc.DB_XREF_KOFAM}:{kofam_id}', so_id] coverage = length / consensus_length if(coverage < HIT_COVERAGE_TRUNCATED): @@ -110,12 +120,7 @@ def predict_r_rnas(data: dict, sequences_path: Path): rrna['start'] = start rrna['stop'] = stop rrna['strand'] = bc.STRAND_FORWARD if strand == '+' else bc.STRAND_REVERSE - if(accession == 'RF00001'): - rrna['gene'] = 'rrf' - elif(accession == 'RF00177'): - rrna['gene'] = 'rrs' - elif(accession == 'RF02541'): - rrna['gene'] = 'rrl' + rrna['gene'] = rrna_gene rrna['product'] = f'{rrna_tag} ribosomal RNA' if(truncated): diff --git a/bakta/features/s_orf.py b/bakta/features/s_orf.py index 84c67065..1a63a802 100644 --- a/bakta/features/s_orf.py +++ b/bakta/features/s_orf.py @@ -25,6 +25,7 @@ def extract(data: dict): orfs = [] for seq in data['sequences']: nt_seq = Seq(seq['nt']) + start_codons = bc.DOMAIN_PROFILES[cfg.domain or bc.DOMAIN_BACTERIA]['sorf_start_codons'] for strand, strand_nt_seq in [(bc.STRAND_FORWARD, nt_seq), (bc.STRAND_REVERSE, nt_seq.reverse_complement())]: # strands +/- for frame in range(3): # frames 1/2/3 -> 0, 1, 2 seq_frame = strand_nt_seq[frame:] @@ -35,12 +36,15 @@ def extract(data: dict): seq_frame = seq_frame[:-residue] aa_seq = str(seq_frame.translate(table=cfg.translation_table, stop_symbol='*', to_stop=False, cds=False)) - aa_start = aa_seq.find('M') - aa_end = aa_seq.find('*', aa_start) - while aa_start > -1 and aa_end > -1: + start_positions = [pos for pos in range(0, len(seq_frame) - 2, 3) if str(seq_frame[pos:pos + 3]).upper() in start_codons] + for nt_start in start_positions: + aa_start = nt_start // 3 + aa_end = aa_seq.find('*', aa_start) + if(aa_end == -1): + continue orf_length = aa_end - aa_start if(orf_length >= bc.MIN_SORF_LENGTH and orf_length < bc.MAX_SORF_LENGTH): # get all CDS starts (M) - aa = aa_seq[aa_start:aa_end] + aa = 'M' + aa_seq[aa_start + 1:aa_end] (aa_digest, aa_hexdigest) = bu.calc_aa_hash(aa) if(strand == bc.STRAND_FORWARD): dna_start = aa_start * 3 + frame + 1 # +1: 0 based idx to 1 based @@ -71,9 +75,6 @@ def extract(data: dict): 'seq=%s, start=%i, stop=%i, strand=%s, frame=%i, aa-length=%i, aa=%s, nt=[%s..%s]', seq['id'], sorf['start'], sorf['stop'], strand, frame, len(aa), aa, nt[:10], nt[-10:] ) - aa_start = aa_seq.find('M', aa_start + 1) - if(aa_start > aa_end): - aa_end = aa_seq.find('*', aa_start) log.info('predicted=%i', len(orfs)) return orfs diff --git a/bakta/features/t_rna.py b/bakta/features/t_rna.py index 083fe28d..b75dfa3a 100644 --- a/bakta/features/t_rna.py +++ b/bakta/features/t_rna.py @@ -49,7 +49,7 @@ def predict_t_rnas(data: dict, sequences_path: Path): fasta_output_path = cfg.tmp_path.joinpath('trna.fasta') cmd = [ 'tRNAscan-SE', - '-B', + '-A' if cfg.domain == bc.DOMAIN_ARCHAEA else '-B', '--output', str(txt_output_path), '--fasta', str(fasta_output_path), '--thread', str(cfg.threads), @@ -100,6 +100,11 @@ def predict_t_rnas(data: dict, sequences_path: Path): trna[bc.PSEUDOGENE] = True trna['score'] = float(score) + # Archaeal tRNAs can contain introns. Preserve their coordinates so + # writers and downstream consumers do not silently discard them. + intron_begin, bounds_end = int(intron_begin), int(bounds_end) + if(intron_begin > 0 and bounds_end > 0): + trna['intron'] = {'start': min(intron_begin, bounds_end), 'stop': max(intron_begin, bounds_end)} nt = bu.extract_feature_sequence(trna, sequences[sequence_id]) # extract nt sequences trna['nt'] = nt diff --git a/bakta/features/tm_rna.py b/bakta/features/tm_rna.py index cd4d1436..ac25a4f5 100644 --- a/bakta/features/tm_rna.py +++ b/bakta/features/tm_rna.py @@ -71,8 +71,9 @@ def predict_tm_rnas(data: dict, sequences_path: Path): tmrna['start'] = start tmrna['stop'] = stop tmrna['strand'] = strand - tmrna['gene'] = 'ssrA' - tmrna['product'] = 'transfer-messenger RNA, SsrA' + profile = bc.DOMAIN_PROFILES[cfg.domain or bc.DOMAIN_BACTERIA] + tmrna['gene'] = profile['tmrna_gene'] + tmrna['product'] = profile['tmrna_product'] tmrna['db_xrefs'] = [so.SO_TMRNA.id] tmrna['tag'] = { 'start': start + tag_start - 1, diff --git a/bakta/io/gff.py b/bakta/io/gff.py index fa5d6a9f..c7fe77e3 100644 --- a/bakta/io/gff.py +++ b/bakta/io/gff.py @@ -87,6 +87,10 @@ def write_features(data: dict, features_by_sequence: Dict[str, dict], gff3_path: fh.write(f"{seq_id}\ttRNAscan-SE\tgene\t{start}\t{stop}\t.\t{feat['strand']}\t.\t{gene_annotations}\n") annotations = encode_annotations(annotations) fh.write(f"{seq_id}\ttRNAscan-SE\t{so.SO_TRNA.name}\t{start}\t{stop}\t.\t{feat['strand']}\t.\t{annotations}\n") + if('intron' in feat): + intron = feat['intron'] + intron_annotations = encode_annotations({'ID': f"{feat['locus']}.intron", 'Parent': feat['locus']}) + fh.write(f"{seq_id}\ttRNAscan-SE\tintron\t{intron['start']}\t{intron['stop']}\t.\t{feat['strand']}\t.\t{intron_annotations}\n") elif(feat['type'] == bc.FEATURE_TM_RNA): annotations = { 'ID': feat['locus'], @@ -328,14 +332,17 @@ def write_features(data: dict, features_by_sequence: Dict[str, dict], gff3_path: 'ID': feat['id'], 'Name': feat['product'] } + source = 'Ori-Finder-Arch' if feat.get('inference', '').startswith('ab initio prediction:Ori-Finder-Arch') else 'BLAST+' if(cfg.compliant): annotations['Note'] = feat['product'] else: annotations['product'] = feat['product'] - annotations['inference'] = 'similar to DNA sequence' + annotations['inference'] = feat.get('inference', 'similar to DNA sequence') + if('attributes' in feat): + annotations.update({f"ori_{key}": value for key, value in feat['attributes'].items()}) annotations = encode_annotations(annotations) feat_type = bc.INSDC_FEATURE_ORIGIN_REPLICATION if cfg.compliant else so.SO_ORIC.name - fh.write(f"{seq_id}\tBLAST+\t{feat_type}\t{start}\t{stop}\t.\t{feat['strand']}\t.\t{annotations}\n") + fh.write(f"{seq_id}\t{source}\t{feat_type}\t{start}\t{stop}\t.\t{feat['strand']}\t.\t{annotations}\n") elif(feat['type'] == bc.FEATURE_ORIV): annotations = { 'ID': feat['id'], diff --git a/bakta/io/insdc.py b/bakta/io/insdc.py index e2207750..bc6a2021 100644 --- a/bakta/io/insdc.py +++ b/bakta/io/insdc.py @@ -114,8 +114,10 @@ def build_biopython_sequence_list(data: dict, features: Sequence[dict]): elif(feature['type'] == bc.FEATURE_ORIC or feature['type'] == bc.FEATURE_ORIV): # TODO: Add fuzzy positions for oriC/oriV insdc_feature_type = bc.INSDC_FEATURE_ORIGIN_REPLICATION - qualifiers['inference'] = 'similar to DNA sequence' + qualifiers['inference'] = feature.get('inference', 'similar to DNA sequence') qualifiers['note'].append(feature['product']) + if('attributes' in feature): + qualifiers['note'].extend([f"Ori-Finder-Arch {key}: {value}" for key, value in feature['attributes'].items()]) if('product' in qualifiers): qualifiers['note'] = feature['product'] del qualifiers['product'] @@ -234,7 +236,18 @@ def build_biopython_sequence_list(data: dict, features: Sequence[dict]): start = feature['start'] - 1 stop = feature['stop'] - if('edge' in feature): + if(feature['type'] == bc.FEATURE_T_RNA and 'intron' in feature and 'edge' not in feature): + intron = feature['intron'] + intron_start = intron['start'] - 1 + intron_stop = intron['stop'] + if(start < intron_start < intron_stop < stop): + left_exon = FeatureLocation(start, intron_start, strand=strand) + right_exon = FeatureLocation(intron_stop, stop, strand=strand) + feature_location = CompoundLocation([right_exon, left_exon] if strand == -1 else [left_exon, right_exon]) + else: + log.warning('ignore invalid tRNA intron coordinates: locus=%s, intron=%s..%s', feature.get('locus'), intron['start'], intron['stop']) + feature_location = FeatureLocation(start, stop, strand=strand) + elif('edge' in feature): fl_1 = FeatureLocation(start, seq['length'], strand=strand) fl_2 = FeatureLocation(0, stop, strand=strand) if(feature['strand'] == bc.STRAND_REVERSE): diff --git a/bakta/main.py b/bakta/main.py index 9306140d..8a036b29 100755 --- a/bakta/main.py +++ b/bakta/main.py @@ -28,6 +28,7 @@ import bakta.features.s_orf as s_orf import bakta.features.gaps as gaps import bakta.features.ori as ori +import bakta.features.ori_arch as ori_arch import bakta.db as db import bakta.utils as bu import bakta.ups as ups @@ -55,7 +56,7 @@ def main(): # - test binary dependencies ############################################################################ cfg.setup(args) # check parameters and prepare global configuration - cfg.db_info = db.check(cfg.db_path) + cfg.db_info = db.check(cfg.db_path, cfg.domain) bu.test_dependencies() if(cfg.verbose): print(f'Bakta v{cfg.version}') @@ -68,6 +69,7 @@ def main(): if(cfg.user_proteins): print(f'\tuser proteins: {cfg.user_proteins}') if(cfg.user_hmms): print(f'\tuser hmms: {cfg.user_hmms}') print(f'\ttranslation table: {cfg.translation_table}') + print(f'\tdomain: {cfg.domain}') if(cfg.taxon): print(f'\ttaxon: {cfg.taxon}') if(cfg.plasmid): print(f'\tplasmid: {cfg.plasmid}') if(cfg.gram != '?'): print(f'\tgram: {cfg.gram}') @@ -142,6 +144,7 @@ def main(): 'taxon': cfg.taxon, 'complete': cfg.complete or complete_genome, 'gram': cfg.gram, + 'domain': cfg.domain, 'translation_table': cfg.translation_table }, 'stats': { @@ -297,11 +300,15 @@ def main(): print('\tconduct expert systems...') # conduct expert systems annotation cds_aa_path = cfg.tmp_path.joinpath('cds.expert.faa') orf.write_internal_faa(cdss, cds_aa_path) - log.debug('conduct expert system: amrfinder') - expert_amr_found = exp_amr.search(cdss, cds_aa_path) - print(f'\t\tamrfinder: {len(expert_amr_found)}') + if(cfg.domain == bc.DOMAIN_BACTERIA): + log.debug('conduct expert system: amrfinder') + expert_amr_found = exp_amr.search(cdss, cds_aa_path) + print(f'\t\tamrfinder: {len(expert_amr_found)}') + else: + print('\t\tskip bacterial AMRFinderPlus expert system for archaeal profile') log.debug('conduct expert system: aa seqs') - diamond_db_path = cfg.db_path.joinpath('expert-protein-sequences.dmnd') + expert_db_name = 'expert-protein-sequences-archaea.dmnd' if cfg.domain == bc.DOMAIN_ARCHAEA else 'expert-protein-sequences.dmnd' + diamond_db_path = cfg.db_path.joinpath(expert_db_name) expert_aa_found = exp_aa_seq.search(cdss, cds_aa_path, 'expert_proteins', diamond_db_path) print(f'\t\tprotein sequences: {len(expert_aa_found)}') @@ -436,15 +443,18 @@ def main(): else: print('detect oriCs/oriVs...') log.debug('detect oriC/V') - oriCs = ori.predict_oris(data, sequences_path, bc.FEATURE_ORIC) + oriCs = ori_arch.predict_orics(data, sequences_path) if cfg.domain == bc.DOMAIN_ARCHAEA else ori.predict_oris(data, sequences_path, bc.FEATURE_ORIC) data['features'].extend(oriCs) print(f'\tfound: {len(oriCs)}') - print('detect oriTs...') - log.debug('detect oriT') - oriTs = ori.predict_oris(data, sequences_path, bc.FEATURE_ORIT) - data['features'].extend(oriTs) - print(f'\tfound: {len(oriTs)}') + if(cfg.domain == bc.DOMAIN_ARCHAEA): + print('skip oriT detection for archaeal profile...') + else: + print('detect oriTs...') + log.debug('detect oriT') + oriTs = ori.predict_oris(data, sequences_path, bc.FEATURE_ORIT) + data['features'].extend(oriTs) + print(f'\tfound: {len(oriTs)}') ############################################################################ # Filter overlapping features diff --git a/bakta/proteins.py b/bakta/proteins.py index c40b1f86..cd5ecc33 100644 --- a/bakta/proteins.py +++ b/bakta/proteins.py @@ -41,6 +41,9 @@ def main(): arg_group_io.add_argument('--prefix', '-p', action='store', default=None, help='Prefix for output files') arg_group_io.add_argument('--force', '-f', action='store_true', help='Force overwriting existing output folder') + arg_group_domain = parser.add_argument_group('Domain') + arg_group_domain.add_argument('--domain', action='store', default=bc.DOMAIN_BACTERIA, choices=[bc.DOMAIN_BACTERIA, bc.DOMAIN_ARCHAEA], help='Annotation domain: bacteria/archaea (default = bacteria)') + arg_group_annotation = parser.add_argument_group('Annotation') arg_group_annotation.add_argument('--proteins', action='store', default=None, dest='proteins', help='Fasta file of trusted protein sequences') arg_group_annotation.add_argument('--hmms', action='store', default=None, dest='hmms', help='HMM file of trusted hidden markov models in HMMER format') @@ -83,8 +86,9 @@ def main(): sys.exit(f'ERROR: input proteins file ({args.input}) not valid!') log.info('input-path=%s', aa_path) + cfg.domain = args.domain cfg.check_db_path(args) - cfg.db_info = db.check(cfg.db_path) + cfg.db_info = db.check(cfg.db_path, cfg.domain) cfg.check_tmp_path(args) cfg.check_user_proteins(args) cfg.check_threads(args) @@ -251,12 +255,16 @@ def annotate_aa(aas: Sequence[dict]): print('\tconduct expert systems...') # conduct expert systems annotation aa_path = cfg.tmp_path.joinpath('aa.faa') orf.write_internal_faa(aas, aa_path) - log.debug('conduct expert system: amrfinder') - cfg.translation_table = 11 - expert_amr_found = exp_amr.search(aas, aa_path) - print(f'\t\tamrfinder: {len(expert_amr_found)}') + if(cfg.domain == bc.DOMAIN_BACTERIA): + log.debug('conduct expert system: amrfinder') + cfg.translation_table = 11 + expert_amr_found = exp_amr.search(aas, aa_path) + print(f'\t\tamrfinder: {len(expert_amr_found)}') + else: + print('\t\tskip bacterial AMRFinderPlus expert system for archaeal profile') log.debug('conduct expert system: aa seqs') - diamond_db_path = cfg.db_path.joinpath('expert-protein-sequences.dmnd') + expert_db_name = 'expert-protein-sequences-archaea.dmnd' if cfg.domain == bc.DOMAIN_ARCHAEA else 'expert-protein-sequences.dmnd' + diamond_db_path = cfg.db_path.joinpath(expert_db_name) expert_aa_found = exp_aa_seq.search(aas, aa_path, 'expert_proteins', diamond_db_path) print(f'\t\tprotein sequences: {len(expert_aa_found)}') if(cfg.user_proteins): diff --git a/bakta/utils.py b/bakta/utils.py index 37e9e5e2..d4a58905 100644 --- a/bakta/utils.py +++ b/bakta/utils.py @@ -78,6 +78,7 @@ def parse_arguments(): arg_group_io.add_argument('--force', '-f', action='store_true', help='Force overwriting existing output folder (except for current working directory)') arg_group_organism = parser.add_argument_group('Organism') + arg_group_organism.add_argument('--domain', action='store', default=bc.DOMAIN_BACTERIA, choices=[bc.DOMAIN_BACTERIA, bc.DOMAIN_ARCHAEA], help='Organism domain profile: bacteria/archaea (default = bacteria)') arg_group_organism.add_argument('--genus', action='store', default=None, help='Genus name') arg_group_organism.add_argument('--species', action='store', default=None, help='Species name') arg_group_organism.add_argument('--strain', action='store', default=None, help='Strain name') @@ -87,7 +88,7 @@ def parse_arguments(): arg_group_annotation.add_argument('--complete', action='store_true', help='All sequences are complete replicons (chromosome/plasmid[s])') arg_group_annotation.add_argument('--prodigal-tf', action='store', default=None, dest='prodigal_tf', help='Path to existing Prodigal training file to use for CDS prediction') arg_group_annotation.add_argument('--translation-table', action='store', type=int, default=11, choices=[11, 4, 25], dest='translation_table', help='Translation table: 11/4/25 (default = 11)') - arg_group_annotation.add_argument('--gram', action='store', default=bc.GRAM_UNKNOWN, choices=[bc.GRAM_POSITIVE, bc.GRAM_NEGATIVE, bc.GRAM_UNKNOWN], help=f'Gram type for signal peptide predictions: {bc.GRAM_POSITIVE}/{bc.GRAM_NEGATIVE}/{bc.GRAM_UNKNOWN} (default = {bc.GRAM_UNKNOWN})') + arg_group_annotation.add_argument('--gram', action='store', default=bc.GRAM_UNKNOWN, choices=[bc.GRAM_POSITIVE, bc.GRAM_NEGATIVE, bc.GRAM_UNKNOWN], help=f'Gram type for signal peptide predictions: {bc.GRAM_POSITIVE}/{bc.GRAM_NEGATIVE}/{bc.GRAM_UNKNOWN} (default = {bc.GRAM_UNKNOWN}; ignored for --domain archaea)') arg_group_annotation.add_argument('--locus', action='store', default=None, help="Locus prefix (default = 'contig')") arg_group_annotation.add_argument('--locus-tag', action='store', default=None, dest='locus_tag', help='Locus tag prefix (default = autogenerated)') arg_group_annotation.add_argument('--locus-tag-increment', action='store', type=int, default=1, choices=[1, 5, 10], dest='locus_tag_increment', help='Locus tag increment: 1/5/10 (default = 1)') @@ -238,27 +239,29 @@ def test_dependencies(): if(cfg.skip_cds is not None and cfg.skip_cds is False): test_dependency(DEPENDENCY_PYRODIGAL) - test_dependency(DEPENDENCY_AMRFINDERPLUS) + if(cfg.domain == bc.DOMAIN_BACTERIA): + test_dependency(DEPENDENCY_AMRFINDERPLUS) # test if AMRFinderPlus db is installed - amrfinderplus_db_path = cfg.db_path.joinpath('amrfinderplus-db') - amrfinderplus_db_latest_path = amrfinderplus_db_path.joinpath('latest') - process = sp.run( - [ - 'amrfinder', - '--debug', - '--database', str(amrfinderplus_db_latest_path) - ], capture_output=True) - if('No valid AMRFinder database found' in process.stderr.decode()): - log.error('AMRFinderPlus database not installed') - sys.exit(f"ERROR: AMRFinderPlus database not installed! Please, install AMRFinderPlus's internal database by executing: 'amrfinder_update --database {amrfinderplus_db_path}'. This must be done only once.") + if(cfg.domain == bc.DOMAIN_BACTERIA): + amrfinderplus_db_path = cfg.db_path.joinpath('amrfinderplus-db') + amrfinderplus_db_latest_path = amrfinderplus_db_path.joinpath('latest') + process = sp.run(['amrfinder', '--debug', '--database', str(amrfinderplus_db_latest_path)], capture_output=True) + if('No valid AMRFinder database found' in process.stderr.decode()): + log.error('AMRFinderPlus database not installed') + sys.exit(f"ERROR: AMRFinderPlus database not installed! Please, install AMRFinderPlus's internal database by executing: 'amrfinder_update --database {amrfinderplus_db_path}'. This must be done only once.") if((cfg.skip_cds is not None and cfg.skip_cds is False) or (cfg.skip_sorf is not None and cfg.skip_sorf is False)): test_dependency(DEPENDENCY_PYHMMER) test_dependency(DEPENDENCY_DIAMOND) if(cfg.skip_ori is not None and cfg.skip_ori is False): - test_dependency(DEPENDENCY_BLASTN) + if(cfg.domain == bc.DOMAIN_ARCHAEA): + for executable in ['OriFinderArch', 'fimo', 'hmmsearch']: + if(shutil.which(executable) is None): + sys.exit(f'ERROR: Archaeal oriC prediction requires {executable} on PATH. Install Ori-Finder-Arch, MEME Suite, and HMMER, or use --skip-ori.') + else: + test_dependency(DEPENDENCY_BLASTN) if(cfg.skip_plot is not None and cfg.skip_plot is False): test_dependency(DEPENDENCY_PYCIRCLIZE) @@ -531,4 +534,4 @@ def get_db_connection(): _local.connection = conn - return _local.connection \ No newline at end of file + return _local.connection diff --git a/db-scripts/annotate-swissprot.py b/db-scripts/annotate-swissprot.py index 5bb321f4..acbb7ad4 100644 --- a/db-scripts/annotate-swissprot.py +++ b/db-scripts/annotate-swissprot.py @@ -16,7 +16,9 @@ parser.add_argument('--taxonomy', action='store', help='Path to NCBI taxonomy node.dmp file.') parser.add_argument('--xml', action='store', help='Path to SwissProt xml file.') parser.add_argument('--db', action='store', help='Path to Bakta sqlite3 db file.') +parser.add_argument('--taxon', action='append', default=None, help='NCBI root taxon to include; repeat for multiple domains (default: 2, Bacteria).') args = parser.parse_args() +taxa = args.taxon or ['2'] taxonomy_path = Path(args.taxonomy).resolve() @@ -81,7 +83,7 @@ def is_taxon_child(child, LCA, taxonomy): tax_id = tax_property.get('id') if tax_property is not None else '1' org_name = elem_org.find('./{*}name[@type="scientific"]') org_name = org_name.text.lower() if org_name is not None else '' - if(is_taxon_child(tax_id, '2', taxonomy) or 'phage' in org_name): + if(any(is_taxon_child(tax_id, taxon, taxonomy) for taxon in taxa) or 'phage' in org_name): seq = elem.find('./{*}sequence').text.upper() seq_hash = hashlib.md5(seq.encode()) seq_hash_hexdigest = seq_hash.hexdigest() diff --git a/db-scripts/buid-db.sh b/db-scripts/buid-db.sh index 30e6f1fd..1375d529 100644 --- a/db-scripts/buid-db.sh +++ b/db-scripts/buid-db.sh @@ -1,104 +1,242 @@ -#!/bin/bash +#!/usr/bin/env bash -set -e +set -euo pipefail -mkdir db +# Build one database containing both domain profiles. The runtime +# --domain option selects its profile-specific assets and predictors. +BAKTA_TAXON_ARGS="--taxon 2 --taxon 2157" + +# Database helper scripts reside beside this entry point. This used to rely on +# an undocumented BAKTA_DB_SCRIPTS environment variable, making direct runs +# fail with paths such as /init-db.py. Resolve it from the script location so +# the build can be started from any working directory. +BAKTA_DB_SCRIPTS="$(cd -- "$(dirname -- "${BASH_SOURCE[0]}")" && pwd -P)" + +mkdir -p db cd db +# Guard against two invocations building into the same db/ directory at once: +# concurrent runs share unnamespaced intermediate files (e.g. rfam-genes.txt) +# and can silently interleave writes into them, corrupting the build. +exec 9>.build.lock +if ! flock -n 9; then + printf "Another buid-db.sh build is already running against %s -- aborting.\n" "$(pwd)" >&2 + exit 1 +fi + +# Every numbered step below is guarded by a .stepNN.done marker, touched only +# as the very last action of a successful step. A killed/crashed run can +# therefore be resumed by simply re-invoking this script from the same db/ +# directory: completed steps are skipped, and steps 1-4 additionally skip +# already-finished bacteria/archaea profiles within an interrupted step. +have_outputs() { + local f + for f in "$@"; do + [ -s "${f}" ] || return 1 + done + return 0 +} + +if have_outputs .build.done; then + printf "Bakta database already fully built in %s -- nothing to do.\n" "$(pwd)" + exit 0 +fi + printf "Create Bakta database\n" # download rRNA covariance models from Rfam printf "\n1/19: download rRNA covariance models from Rfam ...\n" -wget https://ftp.ebi.ac.uk/pub/databases/Rfam/CURRENT/Rfam.cm.gz -pigz -d Rfam.cm.gz -cmfetch Rfam.cm RF00001 > rRNA -cmfetch Rfam.cm RF00177 >> rRNA -cmfetch Rfam.cm RF02541 >> rRNA -cmpress rRNA -rm rRNA +if have_outputs .step01.done; then + printf " already built, skipping.\n" +else + if [ ! -s Rfam.cm ]; then + wget -c https://ftp.ebi.ac.uk/pub/databases/Rfam/CURRENT/Rfam.cm.gz + pigz -df Rfam.cm.gz + fi + for profile in bacteria archaea; do + if [ "${profile}" = "bacteria" ]; then + rRNA_name=rRNA + rRNA_ssu=RF00177 + rRNA_lsu=RF02541 + else + rRNA_name=rRNA-archaea + rRNA_ssu=RF01959 + rRNA_lsu=RF02540 + fi + if have_outputs "${rRNA_name}.i1m" "${rRNA_name}.i1i" "${rRNA_name}.i1f" "${rRNA_name}.i1p"; then + continue + fi + cmfetch Rfam.cm RF00001 > "${rRNA_name}" + cmfetch Rfam.cm "${rRNA_ssu}" >> "${rRNA_name}" + cmfetch Rfam.cm "${rRNA_lsu}" >> "${rRNA_name}" + cmpress -F "${rRNA_name}" + rm "${rRNA_name}" + done + touch .step01.done +fi # download and extract ncRNA gene covariance models from Rfam printf "\n2/19: download ncRNA gene covariance models from Rfam ...\n" -mysql --user rfamro --host mysql-rfam-public.ebi.ac.uk --port 4497 --database Rfam < ${BAKTA_DB_SCRIPTS}/ncRNA-genes.sql | tail -n +2 > rfam-genes.raw.txt -grep "antitoxin;" rfam-genes.raw.txt >> rfam-genes.txt -grep "antisense;" rfam-genes.raw.txt >> rfam-genes.txt -grep "ribozyme;" rfam-genes.raw.txt >> rfam-genes.txt -grep "sRNA;" rfam-genes.raw.txt >> rfam-genes.txt -cut -f1 ${BAKTA_DB_SCRIPTS}/ncRNA-genes.blocklist.txt > ncRNA-genes.blocklist -grep -e "Gene;[^ ]" rfam-genes.raw.txt | grep -v -f ncRNA-genes.blocklist >> rfam-genes.txt -cmfetch -o ncRNA-genes -f Rfam.cm rfam-genes.txt -cmpress ncRNA-genes -wget http://current.geneontology.org/ontology/external2go/rfam2go -awk -F ' ' '{print $1 "\t" $NF}' rfam2go > rfam-go.tsv -rm rfam-genes.raw.txt rfam-genes.txt ncRNA-genes.blocklist ncRNA-genes rfam2go +if have_outputs .step02.done; then + printf " already built, skipping.\n" +else + if [ ! -s Rfam.cm ]; then + wget -c https://ftp.ebi.ac.uk/pub/databases/Rfam/CURRENT/Rfam.cm.gz + pigz -df Rfam.cm.gz + fi + for profile in bacteria archaea; do + if [ "${profile}" = "bacteria" ]; then + RFAM_TAXON=Bacteria + NCRNA_GENES_NAME=ncRNA-genes + else + RFAM_TAXON=Archaea + NCRNA_GENES_NAME=ncRNA-genes-archaea + fi + if have_outputs "${NCRNA_GENES_NAME}.i1m" "${NCRNA_GENES_NAME}.i1i" "${NCRNA_GENES_NAME}.i1f" "${NCRNA_GENES_NAME}.i1p"; then + continue + fi + sed "s/LIKE 'Bacteria%'/LIKE '${RFAM_TAXON}%'/" ${BAKTA_DB_SCRIPTS}/ncRNA-genes.sql | mysql --user rfamro --host mysql-rfam-public.ebi.ac.uk --port 4497 --database Rfam | tail -n +2 > rfam-genes.raw.txt + rm -f rfam-genes.txt + # Not every domain has hits in every category (e.g. archaea has no + # antitoxin entries): grep exits 1 on zero matches, which is fatal + # under `set -e` unless guarded with `|| true`. + grep "antitoxin;" rfam-genes.raw.txt >> rfam-genes.txt || true + grep "antisense;" rfam-genes.raw.txt >> rfam-genes.txt || true + grep "ribozyme;" rfam-genes.raw.txt >> rfam-genes.txt || true + grep "sRNA;" rfam-genes.raw.txt >> rfam-genes.txt || true + cut -f1 ${BAKTA_DB_SCRIPTS}/ncRNA-genes.blocklist.txt > ncRNA-genes.blocklist + grep -e "Gene;[^ ]" rfam-genes.raw.txt | grep -v -f ncRNA-genes.blocklist >> rfam-genes.txt || true + sort -u rfam-genes.txt > rfam-genes.uniq.txt + cmfetch -o "${NCRNA_GENES_NAME}" -f Rfam.cm rfam-genes.uniq.txt + cmpress -F "${NCRNA_GENES_NAME}" + rm rfam-genes.raw.txt rfam-genes.txt rfam-genes.uniq.txt ncRNA-genes.blocklist "${NCRNA_GENES_NAME}" + done + if [ ! -s rfam-go.tsv ]; then + wget -c http://current.geneontology.org/ontology/external2go/rfam2go + awk -F ' ' '{print $1 "\t" $NF}' rfam2go > rfam-go.tsv + rm rfam2go + fi + touch .step02.done +fi # download and extract ncRNA regions (cis reg elements) covariance models from Rfam printf "\n3/19: download ncRNA region covariance models from Rfam ...\n" -mysql --user rfamro --host mysql-rfam-public.ebi.ac.uk --port 4497 --database Rfam < ${BAKTA_DB_SCRIPTS}/ncRNA-regions.sql | tail -n +2 > rfam-regions.raw.txt -grep "riboswitch;" rfam-regions.raw.txt >> rfam-regions.txt -grep "thermoregulator;" rfam-regions.raw.txt >> rfam-regions.txt -grep "leader;" rfam-regions.raw.txt >> rfam-regions.txt -grep "frameshift_element;" rfam-regions.raw.txt >> rfam-regions.txt -cut -f1 ${BAKTA_DB_SCRIPTS}/ncRNA-regions.blocklist.txt > ncRNA-regions.blocklist -grep -e "Cis-reg;[^ ]" rfam-regions.raw.txt | grep -v -f ncRNA-regions.blocklist >> rfam-regions.txt -cmfetch -o ncRNA-regions -f Rfam.cm rfam-regions.txt -cmpress ncRNA-regions -rm rfam-regions.raw.txt rfam-regions.txt ncRNA-regions.blocklist ncRNA-regions Rfam.cm +if have_outputs .step03.done; then + printf " already built, skipping.\n" +else + if [ ! -s Rfam.cm ]; then + wget -c https://ftp.ebi.ac.uk/pub/databases/Rfam/CURRENT/Rfam.cm.gz + pigz -df Rfam.cm.gz + fi + for profile in bacteria archaea; do + if [ "${profile}" = "bacteria" ]; then + RFAM_TAXON=Bacteria + NCRNA_REGIONS_NAME=ncRNA-regions + else + RFAM_TAXON=Archaea + NCRNA_REGIONS_NAME=ncRNA-regions-archaea + fi + if have_outputs "${NCRNA_REGIONS_NAME}.i1m" "${NCRNA_REGIONS_NAME}.i1i" "${NCRNA_REGIONS_NAME}.i1f" "${NCRNA_REGIONS_NAME}.i1p"; then + continue + fi + sed "s/LIKE 'Bacteria%'/LIKE '${RFAM_TAXON}%'/" ${BAKTA_DB_SCRIPTS}/ncRNA-regions.sql | mysql --user rfamro --host mysql-rfam-public.ebi.ac.uk --port 4497 --database Rfam | tail -n +2 > rfam-regions.raw.txt + rm -f rfam-regions.txt + grep "riboswitch;" rfam-regions.raw.txt >> rfam-regions.txt || true + grep "thermoregulator;" rfam-regions.raw.txt >> rfam-regions.txt || true + grep "leader;" rfam-regions.raw.txt >> rfam-regions.txt || true + grep "frameshift_element;" rfam-regions.raw.txt >> rfam-regions.txt || true + cut -f1 ${BAKTA_DB_SCRIPTS}/ncRNA-regions.blocklist.txt > ncRNA-regions.blocklist + grep -e "Cis-reg;[^ ]" rfam-regions.raw.txt | grep -v -f ncRNA-regions.blocklist >> rfam-regions.txt || true + sort -u rfam-regions.txt > rfam-regions.uniq.txt + cmfetch -o "${NCRNA_REGIONS_NAME}" -f Rfam.cm rfam-regions.uniq.txt + cmpress -F "${NCRNA_REGIONS_NAME}" + rm rfam-regions.raw.txt rfam-regions.txt rfam-regions.uniq.txt ncRNA-regions.blocklist "${NCRNA_REGIONS_NAME}" + done + touch .step03.done +fi +rm -f Rfam.cm # download and extract spurious ORF HMMs from AntiFam printf "\n4/19: download and extract spurious ORF HMMs from AntiFam ...\n" -mkdir antifam-dir -cd antifam-dir -wget https://ftp.ebi.ac.uk/pub/databases/Pfam/AntiFam/current/Antifam.tar.gz -tar -xzf Antifam.tar.gz -cd .. -mv antifam-dir/AntiFam_Bacteria.hmm antifam -hmmpress antifam -rm -r antifam antifam-dir/ - - -# download & extract oriT sequences +if have_outputs .step04.done; then + printf " already built, skipping.\n" +else + if [ ! -s antifam-dir/AntiFam_Bacteria.hmm ] || [ ! -s antifam-dir/AntiFam_Archaea.hmm ]; then + rm -rf antifam-dir + mkdir antifam-dir + (cd antifam-dir && wget -c https://ftp.ebi.ac.uk/pub/databases/Pfam/AntiFam/current/Antifam.tar.gz && tar -xzf Antifam.tar.gz) + fi + if ! have_outputs antifam.h3m antifam.h3i antifam.h3f antifam.h3p; then + cp antifam-dir/AntiFam_Bacteria.hmm antifam + hmmpress -f antifam + fi + if ! have_outputs antifam-archaea.h3m antifam-archaea.h3i antifam-archaea.h3f antifam-archaea.h3p; then + cp antifam-dir/AntiFam_Archaea.hmm antifam-archaea + hmmpress -f antifam-archaea + fi + rm -f antifam antifam-archaea + rm -rf antifam-dir/ + touch .step04.done +fi + + +# DoriC and MOB-suite records are retained for the bacterial profile only; +# archaeal annotation uses Ori-Finder-Arch at runtime. printf "\n5/19: download and extract oriT sequences from Mob-suite ...\n" -wget https://zenodo.org/records/10304948/files/data.tar.gz -tar -xvzf data.tar.gz -mv data/orit.fas ./orit.fna -rm -r data/ data.tar.gz +if have_outputs orit.fna; then + printf " already built, skipping.\n" +else + wget -c https://zenodo.org/records/10304948/files/data.tar.gz + tar -xvzf data.tar.gz + mv data/orit.fas ./orit.fna + rm -r data/ data.tar.gz +fi printf "\n5/19: download oriC/V sequences from DoriC ...\n" -curl 'https://tubic.org/doric/search/bacteria' \ - -H 'content-type: multipart/form-data; boundary=----WebKitFormBoundaryBDBZTWpS3orCjS0m' \ - --data-raw $'------WebKitFormBoundaryBDBZTWpS3orCjS0m\r\nContent-Disposition: form-data; name="assembly_level"\r\n\r\nComplete\r\n------WebKitFormBoundaryBDBZTWpS3orCjS0m\r\nContent-Disposition: form-data; name="topology"\r\n\r\nAll\r\n------WebKitFormBoundaryBDBZTWpS3orCjS0m\r\nContent-Disposition: form-data; name="chromosome_type"\r\n\r\nAll\r\n------WebKitFormBoundaryBDBZTWpS3orCjS0m\r\nContent-Disposition: form-data; name="oric_type"\r\n\r\nSingle\r\n------WebKitFormBoundaryBDBZTWpS3orCjS0m\r\nContent-Disposition: form-data; name="organism"\r\n\r\n\r\n------WebKitFormBoundaryBDBZTWpS3orCjS0m\r\nContent-Disposition: form-data; name="lineage"\r\n\r\n\r\n------WebKitFormBoundaryBDBZTWpS3orCjS0m\r\nContent-Disposition: form-data; name="download1"\r\n\r\nDownload\r\n------WebKitFormBoundaryBDBZTWpS3orCjS0m--\r\n' \ - --compressed > oric.csv -curl 'https://tubic.org/doric/search/plasmid' \ - -H 'content-type: multipart/form-data; boundary=----WebKitFormBoundaryMu32WgFUyqC7TO0d' \ - --data-raw $'------WebKitFormBoundaryMu32WgFUyqC7TO0d\r\nContent-Disposition: form-data; name="topology"\r\n\r\nAll\r\n------WebKitFormBoundaryMu32WgFUyqC7TO0d\r\nContent-Disposition: form-data; name="organism"\r\n\r\n\r\n------WebKitFormBoundaryMu32WgFUyqC7TO0d\r\nContent-Disposition: form-data; name="lineage"\r\n\r\n\r\n------WebKitFormBoundaryMu32WgFUyqC7TO0d\r\nContent-Disposition: form-data; name="download1"\r\n\r\nDownload\r\n------WebKitFormBoundaryMu32WgFUyqC7TO0d--\r\n' \ - --compressed > oriv.csv -python3 ${BAKTA_DB_SCRIPTS}/extract-ori.py --doric oric.csv --fasta ori.chromosome.fna -python3 ${BAKTA_DB_SCRIPTS}/extract-ori.py --doric oriv.csv --fasta ori.plasmid.fna -cat ori.chromosome.fna > oric.raw.fna -cat ori.plasmid.fna >> oric.raw.fna -cd-hit-est -i oric.raw.fna -o oric.fna -c 0.99 -s 0.99 -aS 0.99 -g 1 -r 1 -rm *.csv ori.*.fna oric.raw.fna oric.fna.clstr +if have_outputs oric.fna; then + printf " already built, skipping.\n" +else + curl 'https://tubic.org/doric/search/bacteria' \ + -H 'content-type: multipart/form-data; boundary=----WebKitFormBoundaryBDBZTWpS3orCjS0m' \ + --data-raw $'------WebKitFormBoundaryBDBZTWpS3orCjS0m\r\nContent-Disposition: form-data; name="assembly_level"\r\n\r\nComplete\r\n------WebKitFormBoundaryBDBZTWpS3orCjS0m\r\nContent-Disposition: form-data; name="topology"\r\n\r\nAll\r\n------WebKitFormBoundaryBDBZTWpS3orCjS0m\r\nContent-Disposition: form-data; name="chromosome_type"\r\n\r\nAll\r\n------WebKitFormBoundaryBDBZTWpS3orCjS0m\r\nContent-Disposition: form-data; name="oric_type"\r\n\r\nSingle\r\n------WebKitFormBoundaryBDBZTWpS3orCjS0m\r\nContent-Disposition: form-data; name="organism"\r\n\r\n\r\n------WebKitFormBoundaryBDBZTWpS3orCjS0m\r\nContent-Disposition: form-data; name="lineage"\r\n\r\n\r\n------WebKitFormBoundaryBDBZTWpS3orCjS0m\r\nContent-Disposition: form-data; name="download1"\r\n\r\nDownload\r\n------WebKitFormBoundaryBDBZTWpS3orCjS0m--\r\n' \ + --compressed > oric.csv + curl 'https://tubic.org/doric/search/plasmid' \ + -H 'content-type: multipart/form-data; boundary=----WebKitFormBoundaryMu32WgFUyqC7TO0d' \ + --data-raw $'------WebKitFormBoundaryMu32WgFUyqC7TO0d\r\nContent-Disposition: form-data; name="topology"\r\n\r\nAll\r\n------WebKitFormBoundaryMu32WgFUyqC7TO0d\r\nContent-Disposition: form-data; name="organism"\r\n\r\n\r\n------WebKitFormBoundaryMu32WgFUyqC7TO0d\r\nContent-Disposition: form-data; name="lineage"\r\n\r\n\r\n------WebKitFormBoundaryMu32WgFUyqC7TO0d\r\nContent-Disposition: form-data; name="download1"\r\n\r\nDownload\r\n------WebKitFormBoundaryMu32WgFUyqC7TO0d--\r\n' \ + --compressed > oriv.csv + python3 ${BAKTA_DB_SCRIPTS}/extract-ori.py --doric oric.csv --fasta ori.chromosome.fna + python3 ${BAKTA_DB_SCRIPTS}/extract-ori.py --doric oriv.csv --fasta ori.plasmid.fna + cat ori.chromosome.fna > oric.raw.fna + cat ori.plasmid.fna >> oric.raw.fna + cd-hit-est -i oric.raw.fna -o oric.fna -c 0.99 -s 0.99 -aS 0.99 -g 1 -r 1 + rm *.csv ori.*.fna oric.raw.fna oric.fna.clstr +fi # download NCBI Taxonomy DB printf "\n6/19: download NCBI Taxonomy DB ...\n" -mkdir taxonomy -cd taxonomy -wget https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz -tar -I pigz -xf taxdump.tar.gz -cd .. -mv taxonomy/nodes.dmp . -rm -rf taxonomy +if have_outputs nodes.dmp; then + printf " already built, skipping.\n" +else + rm -rf taxonomy + mkdir taxonomy + (cd taxonomy && wget -c https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz && tar -I pigz -xf taxdump.tar.gz) + mv taxonomy/nodes.dmp . + rm -rf taxonomy +fi ############################################################################ # Setup SQLite Bakta db ############################################################################ printf "\n7/19: setup SQLite Bakta db ...\n" -python3 ${BAKTA_DB_SCRIPTS}/init-db.py --db bakta.db +if have_outputs bakta.db; then + printf " already built, skipping.\n" +else + python3 ${BAKTA_DB_SCRIPTS}/init-db.py --db bakta.db +fi ############################################################################ @@ -108,24 +246,40 @@ python3 ${BAKTA_DB_SCRIPTS}/init-db.py --db bakta.db # - build PSCC Diamond db ############################################################################ printf "\n8/19: download UniProt UniRef50 ...\n" -wget https://ftp.expasy.org/databases/uniprot/current_release/uniref/uniref50/uniref50.xml.gz -for i in {1..200}; do - wget https://ftp.expasy.org/databases/uniprot/current_release/uniparc/fasta/active/uniparc_active_p${i}.fasta.gz - pigz -dc uniparc_active_p${i}.fasta.gz >> uniparc_active.fasta - rm uniparc_active_p${i}.fasta.gz -done -printf "\n8/19: read UniRef90 entries and build Protein Sequence Cluster sequence and information databases:\n" -python3 ${BAKTA_DB_SCRIPTS}/init-pscc.py --taxonomy nodes.dmp --uniref50 uniref50.xml.gz --uniparc uniparc_active.fasta --db bakta.db --pscc pscc.faa --pscc_sorf pscc_sorf.faa -printf "\n8/19: build PSCC Diamond db ...\n" -diamond makedb --in pscc.faa --db pscc -diamond makedb --in pscc_sorf.faa --db sorf -mkdir db-light -cp bakta.db db-light/ -mv pscc.dmnd sorf.dmnd db-light/ -cd db-light -python3 ${BAKTA_DB_SCRIPTS}/optimize-db.py --db bakta.db -cd .. -rm uniref50.xml.gz +if have_outputs .step08.done; then + printf " already built, skipping.\n" +else + if [ ! -s uniref50.xml.gz ]; then + wget -c https://ftp.expasy.org/databases/uniprot/current_release/uniref/uniref50/uniref50.xml.gz + fi + if ! have_outputs uniparc_active.fasta; then + printf "\n8/19: download UniParc active sequences (200 parts, resumable) ...\n" + for i in {1..200}; do + part="uniparc_active_p${i}.fasta.gz" + if [ ! -s "${part}" ]; then + wget -c "https://ftp.expasy.org/databases/uniprot/current_release/uniparc/fasta/active/${part}" + fi + done + printf "\n8/19: merging UniParc parts ...\n" + rm -f uniparc_active.fasta.partial + for i in {1..200}; do + pigz -dc "uniparc_active_p${i}.fasta.gz" >> uniparc_active.fasta.partial + done + mv uniparc_active.fasta.partial uniparc_active.fasta + rm -f uniparc_active_p*.fasta.gz + fi + printf "\n8/19: read UniRef90 entries and build Protein Sequence Cluster sequence and information databases:\n" + python3 ${BAKTA_DB_SCRIPTS}/init-pscc.py --taxonomy nodes.dmp --uniref50 uniref50.xml.gz --uniparc uniparc_active.fasta --db bakta.db --pscc pscc.faa --pscc_sorf pscc_sorf.faa ${BAKTA_TAXON_ARGS} + printf "\n8/19: build PSCC Diamond db ...\n" + diamond makedb --in pscc.faa --db pscc + diamond makedb --in pscc_sorf.faa --db sorf + mkdir -p db-light + cp bakta.db db-light/ + mv pscc.dmnd sorf.dmnd db-light/ + (cd db-light && python3 ${BAKTA_DB_SCRIPTS}/optimize-db.py --db bakta.db) + rm uniref50.xml.gz + touch .step08.done +fi ############################################################################ @@ -135,13 +289,20 @@ rm uniref50.xml.gz # - build PSC Diamond db ############################################################################ printf "\n9/19: download UniProt UniRef90 ...\n" -wget https://ftp.expasy.org/databases/uniprot/current_release/uniref/uniref90/uniref90.xml.gz -printf "\n9/19: read UniRef90 entries and build Protein Sequence Cluster sequence and information databases:\n" -python3 ${BAKTA_DB_SCRIPTS}/init-psc.py --taxonomy nodes.dmp --uniref90 uniref90.xml.gz --uniparc uniparc_active.fasta --db bakta.db --psc psc.faa --psc_sorf sorf.faa -printf "\n9/19: build PSC Diamond db ...\n" -diamond makedb --in psc.faa --db psc -diamond makedb --in sorf.faa --db sorf -rm uniref90.xml.gz +if have_outputs .step09.done; then + printf " already built, skipping.\n" +else + if [ ! -s uniref90.xml.gz ]; then + wget -c https://ftp.expasy.org/databases/uniprot/current_release/uniref/uniref90/uniref90.xml.gz + fi + printf "\n9/19: read UniRef90 entries and build Protein Sequence Cluster sequence and information databases:\n" + python3 ${BAKTA_DB_SCRIPTS}/init-psc.py --taxonomy nodes.dmp --uniref90 uniref90.xml.gz --uniparc uniparc_active.fasta --db bakta.db --psc psc.faa --psc_sorf sorf.faa ${BAKTA_TAXON_ARGS} + printf "\n9/19: build PSC Diamond db ...\n" + diamond makedb --in psc.faa --db psc + diamond makedb --in sorf.faa --db sorf + rm uniref90.xml.gz + touch .step09.done +fi ############################################################################ @@ -150,29 +311,57 @@ rm uniref90.xml.gz # - read, filter and transform UniRef100 entries and store to ips.db ############################################################################ printf "\n10/19: download UniProt UniRef100 ...\n" -wget https://ftp.expasy.org/databases/uniprot/current_release/uniref/uniref100/uniref100.xml.gz -printf "\n10/19: read, filter and store UniRef100 entries ...:\n" -python3 ${BAKTA_DB_SCRIPTS}/init-ups-ips.py --taxonomy nodes.dmp --uniref100 uniref100.xml.gz --uniparc uniparc_active.fasta --db bakta.db --ips ips.faa -rm uniref100.xml.gz uniparc_active.fasta.gz +if have_outputs .step10.done; then + printf " already built, skipping.\n" +else + if [ ! -s uniref100.xml.gz ]; then + wget -c https://ftp.expasy.org/databases/uniprot/current_release/uniref/uniref100/uniref100.xml.gz + fi + printf "\n10/19: read, filter and store UniRef100 entries ...:\n" + python3 ${BAKTA_DB_SCRIPTS}/init-ups-ips.py --taxonomy nodes.dmp --uniref100 uniref100.xml.gz --uniparc uniparc_active.fasta --db bakta.db --ips ips.faa ${BAKTA_TAXON_ARGS} + rm uniref100.xml.gz uniparc_active.fasta + touch .step10.done +fi ############################################################################ # Integrate NCBI nonredundant protein identifiers and COG db -# - download bacterial RefSeq nonredundant proteins and COG files +# - download bacterial and archaeal RefSeq nonredundant proteins and COG files # - annotate UPSs with NCBI nrp IDs (WP_*) # - annotate IPSs/PSCs with COG IDs, gene symbols, product descriptions (seq -> hash -> UniParc/WP_* -> UniRef100 -> UniRef90 -> PSC) ############################################################################ printf "\n11/19: download NCBI COG clusters and RefSeq nonredundant proteins ...\n" -wget https://ftp.ncbi.nih.gov/pub/COG/COG2024/data/cog-24.def.tab # COG IDs and functional class -wget https://ftp.ncbi.nih.gov/pub/COG/COG2024/data/cog-24.cog.csv # Mapping GenBank IDs -> COG IDs -for i in {1..652}; do - wget https://ftp.ncbi.nih.gov/refseq/release/bacteria/bacteria.wp_protein.${i}.protein.faa.gz - pigz -dc bacteria.wp_protein.${i}.protein.faa.gz | seqtk seq -CU >> refseq-bacteria-nrp.trimmed.faa - rm bacteria.wp_protein.${i}.protein.faa.gz -done -printf "\n11/19: annotate IPSs and PSCs ...\n" -python3 ${BAKTA_DB_SCRIPTS}/annotate-ncbi-nrp-cog.py --db bakta.db --nrp refseq-bacteria-nrp.trimmed.faa --cog-ids cog-24.def.tab --nrp-cog-mapping cog-24.cog.csv -rm refseq-bacteria-nrp.trimmed.faa +if have_outputs .step11.done; then + printf " already built, skipping.\n" +else + if [ ! -s cog-24.def.tab ]; then + wget -c https://ftp.ncbi.nih.gov/pub/COG/COG2024/data/cog-24.def.tab # COG IDs and functional class + fi + if [ ! -s cog-24.cog.csv ]; then + wget -c https://ftp.ncbi.nih.gov/pub/COG/COG2024/data/cog-24.cog.csv # Mapping GenBank IDs -> COG IDs + fi + mkdir -p refseq-nrp + for refseq_domain in bacteria archaea; do + refseq_url=https://ftp.ncbi.nih.gov/refseq/release/${refseq_domain} + wget -qO- "${refseq_url}/" | sed -n "s/.*href=\"\(${refseq_domain}\.wp_protein\.[0-9]*\.protein\.faa\.gz\)\".*/\1/p" > "refseq-nrp/${refseq_domain}.files.txt" + while read -r refseq_file; do + [ -s "refseq-nrp/${refseq_file}" ] || wget -c -O "refseq-nrp/${refseq_file}" "${refseq_url}/${refseq_file}" + done < "refseq-nrp/${refseq_domain}.files.txt" + done + printf "\n11/19: merging RefSeq nonredundant proteins ...\n" + rm -f refseq-nrp.trimmed.faa.partial + for refseq_domain in bacteria archaea; do + while read -r refseq_file; do + pigz -dc "refseq-nrp/${refseq_file}" | seqtk seq -CU >> refseq-nrp.trimmed.faa.partial + done < "refseq-nrp/${refseq_domain}.files.txt" + done + mv refseq-nrp.trimmed.faa.partial refseq-nrp.trimmed.faa + rm -rf refseq-nrp + printf "\n11/19: annotate IPSs and PSCs ...\n" + python3 ${BAKTA_DB_SCRIPTS}/annotate-ncbi-nrp-cog.py --db bakta.db --nrp refseq-nrp.trimmed.faa --cog-ids cog-24.def.tab --nrp-cog-mapping cog-24.cog.csv + rm refseq-nrp.trimmed.faa + touch .step11.done +fi ############################################################################ @@ -182,19 +371,30 @@ rm refseq-bacteria-nrp.trimmed.faa # - annotate PSCs ############################################################################ printf "\n12/19: download KEGG kofams HMM models...\n" -wget https://www.genome.jp/ftp/db/kofam/ko_list.gz -wget https://www.genome.jp/ftp/db/kofam/profiles.tar.gz -zcat ko_list.gz | grep full | awk '{ if($5>=0.77) print $0}' > hmms.kofam.selected.tsv -cut -f1 hmms.kofam.selected.tsv > hmms.ids.txt -tar -I pigz -xf profiles.tar.gz -for kofam in `cat profiles/prokaryote.hal`; do cat profiles/$kofam >> kofam-prok; done -hmmfetch -f -o kofams kofam-prok hmms.ids.txt -hmmpress kofams -printf "\n12/19: annotate PSCs...\n" -mkdir -p work/tblout work/domtblout -nextflow run ${BAKTA_DB_SCRIPTS}/hmmsearch.nf --in psc.faa --db kofams --no_tc --out hmmsearch.kofam.tblout -python3 ${BAKTA_DB_SCRIPTS}/annotate-kofams.py --db bakta.db --hmms hmms.kofam.selected.tsv --hmm-results hmmsearch.kofam.tblout -rm -rf profiles ko_list.gz kofam* hmmsearch.kofam.* hmms* +if have_outputs .step12.done; then + printf " already built, skipping.\n" +else + if [ ! -s ko_list.gz ]; then + wget -c https://www.genome.jp/ftp/db/kofam/ko_list.gz + fi + if [ ! -s profiles.tar.gz ]; then + wget -c https://www.genome.jp/ftp/db/kofam/profiles.tar.gz + fi + zcat ko_list.gz | grep full | awk '{ if($5>=0.77) print $0}' > hmms.kofam.selected.tsv + cut -f1 hmms.kofam.selected.tsv > hmms.ids.txt + rm -rf profiles + tar -I pigz -xf profiles.tar.gz + rm -f kofam-prok + for kofam in `cat profiles/prokaryote.hal`; do cat profiles/$kofam >> kofam-prok; done + hmmfetch -f -o kofams kofam-prok hmms.ids.txt + hmmpress -f kofams + printf "\n12/19: annotate PSCs...\n" + mkdir -p work/tblout work/domtblout + nextflow run ${BAKTA_DB_SCRIPTS}/hmmsearch.nf --in psc.faa --db kofams --no_tc --out hmmsearch.kofam.tblout + python3 ${BAKTA_DB_SCRIPTS}/annotate-kofams.py --db bakta.db --hmms hmms.kofam.selected.tsv --hmm-results hmmsearch.kofam.tblout + rm -rf profiles ko_list.gz profiles.tar.gz kofam* hmmsearch.kofam.* hmms* + touch .step12.done +fi ############################################################################ @@ -204,10 +404,17 @@ rm -rf profiles ko_list.gz kofam* hmmsearch.kofam.* hmms* # - annotate IPSs if IPS have no PSC UniRef90 identifier (seq -> hash -> UPS -> IPS) ############################################################################ printf "\n13/19: download UniProt/SwissProt ...\n" -wget https://ftp.expasy.org/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.xml.gz -printf "\n13/19: annotate IPSs and PSCs ...\n" -python3 ${BAKTA_DB_SCRIPTS}/annotate-swissprot.py --taxonomy nodes.dmp --xml uniprot_sprot.xml.gz --db bakta.db -rm uniprot_sprot.xml.gz +if have_outputs .step13.done; then + printf " already built, skipping.\n" +else + if [ ! -s uniprot_sprot.xml.gz ]; then + wget -c https://ftp.expasy.org/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.xml.gz + fi + printf "\n13/19: annotate IPSs and PSCs ...\n" + python3 ${BAKTA_DB_SCRIPTS}/annotate-swissprot.py --taxonomy nodes.dmp --xml uniprot_sprot.xml.gz --db bakta.db ${BAKTA_TAXON_ARGS} + rm uniprot_sprot.xml.gz + touch .step13.done +fi ############################################################################ @@ -216,19 +423,28 @@ rm uniprot_sprot.xml.gz # - annotate PSCs ############################################################################ printf "\n14/19: download NCBIfams HMM models...\n" -wget https://ftp.ncbi.nlm.nih.gov/hmm/current/hmm_PGAP.LIB -wget https://ftp.ncbi.nlm.nih.gov/hmm/current/hmm_PGAP.tsv -grep -v "(Provisional)" hmm_PGAP.tsv > hmms.non-prov.tsv -grep exception hmms.non-prov.tsv > hmms.ncbi.selected.tsv -grep equivalog hmms.non-prov.tsv >> hmms.ncbi.selected.tsv -sort hmms.ncbi.selected.tsv | uniq | cut -f1 > hmms.ids.txt -hmmfetch -f -o ncbifams hmm_PGAP.LIB hmms.ids.txt -hmmpress ncbifams -printf "\n14/19: annotate PSCs...\n" -mkdir -p work/tblout work/domtblout -nextflow run ${BAKTA_DB_SCRIPTS}/hmmsearch.nf --in psc.faa --db ncbifams --block 10000 --out hmmsearch.ncbifams.tblout -python3 ${BAKTA_DB_SCRIPTS}/annotate-ncbi-fams.py --db bakta.db --hmms hmms.ncbi.selected.tsv --hmm-results hmmsearch.ncbifams.tblout -rm ncbifams* hmms.* hmm_PGAP.* hmmsearch.ncbifams.tblout +if have_outputs .step14.done; then + printf " already built, skipping.\n" +else + if [ ! -s hmm_PGAP.LIB ]; then + wget -c https://ftp.ncbi.nlm.nih.gov/hmm/current/hmm_PGAP.LIB + fi + if [ ! -s hmm_PGAP.tsv ]; then + wget -c https://ftp.ncbi.nlm.nih.gov/hmm/current/hmm_PGAP.tsv + fi + grep -v "(Provisional)" hmm_PGAP.tsv > hmms.non-prov.tsv + grep exception hmms.non-prov.tsv > hmms.ncbi.selected.tsv || true + grep equivalog hmms.non-prov.tsv >> hmms.ncbi.selected.tsv || true + sort hmms.ncbi.selected.tsv | uniq | cut -f1 > hmms.ids.txt + hmmfetch -f -o ncbifams hmm_PGAP.LIB hmms.ids.txt + hmmpress -f ncbifams + printf "\n14/19: annotate PSCs...\n" + mkdir -p work/tblout work/domtblout + nextflow run ${BAKTA_DB_SCRIPTS}/hmmsearch.nf --in psc.faa --db ncbifams --block 10000 --out hmmsearch.ncbifams.tblout + python3 ${BAKTA_DB_SCRIPTS}/annotate-ncbi-fams.py --db bakta.db --hmms hmms.ncbi.selected.tsv --hmm-results hmmsearch.ncbifams.tblout + rm ncbifams* hmms.* hmm_PGAP.* hmmsearch.ncbifams.tblout + touch .step14.done +fi ############################################################################ @@ -238,17 +454,28 @@ rm ncbifams* hmms.* hmm_PGAP.* hmmsearch.ncbifams.tblout # - annotate PSCs ############################################################################ printf "\n15/19: download PHROGs ...\n" -wget https://phrogs.lmge.uca.fr/downloads_from_website/FAA_phrog.tar.gz -wget https://phrogs.lmge.uca.fr/downloads_from_website/phrog_annot_v4.tsv -tar -xzf FAA_phrog.tar.gz -cat FAA_phrog/*.faa >> phrogs-raw.faa -python3 ${BAKTA_DB_SCRIPTS}/extract-phrogs.py --annotation phrog_annot_v4.tsv --proteins phrogs-raw.faa --filtered-proteins phrogs.faa -diamond makedb --in phrogs.faa --db phrog -printf "\n15/19: annotate PSCs...\n" -python3 ${BAKTA_DB_SCRIPTS}/extract-hypotheticals.py --psc psc.faa --db bakta.db --hypotheticals hypotheticals.faa -nextflow run ${BAKTA_DB_SCRIPTS}/diamond.nf --in hypotheticals.faa --db phrog.dmnd --block 100000 --id 90 --qcov 80 --scov 80 --out diamond.phrog.psc.tsv -python3 ${BAKTA_DB_SCRIPTS}/annotate-phrogs.py --db bakta.db --annotation phrog_annot_v4.tsv --psc-alignments diamond.phrog.psc.tsv -rm -r FAA_phrog.tar.gz phrog_annot_v4.tsv FAA_phrog phrogs-raw.faa phrogs.faa phrog.dmnd hypotheticals.faa +if have_outputs .step15.done; then + printf " already built, skipping.\n" +else + if [ ! -s FAA_phrog.tar.gz ]; then + wget -c https://phrogs.lmge.uca.fr/downloads_from_website/FAA_phrog.tar.gz + fi + if [ ! -s phrog_annot_v4.tsv ]; then + wget -c https://phrogs.lmge.uca.fr/downloads_from_website/phrog_annot_v4.tsv + fi + rm -rf FAA_phrog + tar -xzf FAA_phrog.tar.gz + rm -f phrogs-raw.faa + cat FAA_phrog/*.faa >> phrogs-raw.faa + python3 ${BAKTA_DB_SCRIPTS}/extract-phrogs.py --annotation phrog_annot_v4.tsv --proteins phrogs-raw.faa --filtered-proteins phrogs.faa + diamond makedb --in phrogs.faa --db phrog + printf "\n15/19: annotate PSCs...\n" + python3 ${BAKTA_DB_SCRIPTS}/extract-hypotheticals.py --psc psc.faa --db bakta.db --hypotheticals hypotheticals.faa + nextflow run ${BAKTA_DB_SCRIPTS}/diamond.nf --in hypotheticals.faa --db phrog.dmnd --block 100000 --id 90 --qcov 80 --scov 80 --out diamond.phrog.psc.tsv + python3 ${BAKTA_DB_SCRIPTS}/annotate-phrogs.py --db bakta.db --annotation phrog_annot_v4.tsv --psc-alignments diamond.phrog.psc.tsv + rm -r FAA_phrog.tar.gz phrog_annot_v4.tsv FAA_phrog phrogs-raw.faa phrogs.faa phrog.dmnd hypotheticals.faa diamond.phrog.psc.tsv + touch .step15.done +fi ############################################################################ @@ -257,10 +484,17 @@ rm -r FAA_phrog.tar.gz phrog_annot_v4.tsv FAA_phrog phrogs-raw.faa phrogs.faa ph # - annotate IPSs with AMR info ############################################################################ printf "\n16/19: download AMR gene WP_* annotations from NCBI Pathogen AMR db ...\n" -wget https://ftp.ncbi.nlm.nih.gov/pathogen/Antimicrobial_resistance/AMRFinderPlus/database/latest/ReferenceGeneCatalog.txt -printf "\n16/19: annotate PSCs...\n" -python3 ${BAKTA_DB_SCRIPTS}/annotate-ncbi-amr.py --db bakta.db --genes ReferenceGeneCatalog.txt -rm ReferenceGeneCatalog.txt +if have_outputs .step16.done; then + printf " already built, skipping.\n" +else + if [ ! -s ReferenceGeneCatalog.txt ]; then + wget -c https://ftp.ncbi.nlm.nih.gov/pathogen/Antimicrobial_resistance/AMRFinderPlus/database/latest/ReferenceGeneCatalog.txt + fi + printf "\n16/19: annotate PSCs...\n" + python3 ${BAKTA_DB_SCRIPTS}/annotate-ncbi-amr.py --db bakta.db --genes ReferenceGeneCatalog.txt + rm ReferenceGeneCatalog.txt + touch .step16.done +fi ############################################################################ @@ -270,14 +504,21 @@ rm ReferenceGeneCatalog.txt # - annotate IPSs/PCSs with IS info ############################################################################ printf "\n17/19: download & extract ISfinder protein sequences ...\n" -wget https://github.com/oschwengers/ISfinder-sequences/raw/2e9162bd5e3448c86ec1549a55315e498bef72fc/IS.faa -python3 ${BAKTA_DB_SCRIPTS}/extract-is.py --input IS.faa --output is.transposase.faa -printf "\n17/19: annotate IPSs/PCSs ...\n" -diamond makedb --in is.transposase.faa --db is -nextflow run ${BAKTA_DB_SCRIPTS}/diamond.nf --in ips.faa --db is.dmnd --block 100000 --id 95 --qcov 90 --scov 90 --out diamond.is.ips.tsv -nextflow run ${BAKTA_DB_SCRIPTS}/diamond.nf --in psc.faa --db is.dmnd --block 100000 --id 90 --qcov 80 --scov 80 --out diamond.is.psc.tsv -python3 ${BAKTA_DB_SCRIPTS}/annotate-is.py --db bakta.db --ips-alignments diamond.is.ips.tsv --psc-alignments diamond.is.psc.tsv -rm is.transposase.faa is.dmnd diamond.is.ips.tsv diamond.is.psc.tsv +if have_outputs .step17.done; then + printf " already built, skipping.\n" +else + if [ ! -s IS.faa ]; then + wget -c https://github.com/oschwengers/ISfinder-sequences/raw/2e9162bd5e3448c86ec1549a55315e498bef72fc/IS.faa + fi + python3 ${BAKTA_DB_SCRIPTS}/extract-is.py --input IS.faa --output is.transposase.faa + printf "\n17/19: annotate IPSs/PCSs ...\n" + diamond makedb --in is.transposase.faa --db is + nextflow run ${BAKTA_DB_SCRIPTS}/diamond.nf --in ips.faa --db is.dmnd --block 100000 --id 95 --qcov 90 --scov 90 --out diamond.is.ips.tsv + nextflow run ${BAKTA_DB_SCRIPTS}/diamond.nf --in psc.faa --db is.dmnd --block 100000 --id 90 --qcov 80 --scov 80 --out diamond.is.psc.tsv + python3 ${BAKTA_DB_SCRIPTS}/annotate-is.py --db bakta.db --ips-alignments diamond.is.ips.tsv --psc-alignments diamond.is.psc.tsv + rm is.transposase.faa is.dmnd diamond.is.ips.tsv diamond.is.psc.tsv + touch .step17.done +fi ############################################################################ @@ -288,19 +529,33 @@ rm is.transposase.faa is.dmnd diamond.is.ips.tsv diamond.is.psc.tsv # - annotate hypothetical PSC via Pfam families ############################################################################ printf "\n18/19: download HMM models from Pfam ...\n" -wget https://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/Pfam-A.hmm.dat.gz -python3 ${BAKTA_DB_SCRIPTS}/extract-pfam.py --pfam Pfam-A.hmm.dat.gz --family pfam.families.tsv --non-family pfam.non-families.tsv -wget https://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/Pfam-A.hmm.gz -pigz -d Pfam-A.hmm.gz -hmmfetch -o pfam-families -f Pfam-A.hmm pfam.families.tsv -hmmpress pfam-families -hmmfetch -o pfam -f Pfam-A.hmm pfam.non-families.tsv -hmmpress pfam -python3 ${BAKTA_DB_SCRIPTS}/extract-hypotheticals.py --psc psc.faa --db bakta.db --hypotheticals hypotheticals.faa -mkdir -p work/tblout work/domtblout -nextflow run ${BAKTA_DB_SCRIPTS}/hmmsearch.nf --in hypotheticals.faa --db pfam-families --block 10000 --out hmmsearch.pfam-families.tblout -python3 ${BAKTA_DB_SCRIPTS}/annotate-pfam.py --db bakta.db --hmms pfam-families --hmm-results hmmsearch.pfam-families.tblout -rm pfam-families* pfam *.tsv Pfam* hmmsearch.pfam-families.tblout +if have_outputs .step18.done; then + printf " already built, skipping.\n" +else + if [ ! -s Pfam-A.hmm.dat.gz ]; then + wget -c https://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/Pfam-A.hmm.dat.gz + fi + python3 ${BAKTA_DB_SCRIPTS}/extract-pfam.py --pfam Pfam-A.hmm.dat.gz --family pfam.families.tsv --non-family pfam.non-families.tsv + if [ ! -s Pfam-A.hmm ]; then + if [ ! -s Pfam-A.hmm.gz ]; then + wget -c https://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/Pfam-A.hmm.gz + fi + pigz -df Pfam-A.hmm.gz + fi + hmmfetch -o pfam-families -f Pfam-A.hmm pfam.families.tsv + hmmpress -f pfam-families + hmmfetch -o pfam -f Pfam-A.hmm pfam.non-families.tsv + hmmpress -f pfam + python3 ${BAKTA_DB_SCRIPTS}/extract-hypotheticals.py --psc psc.faa --db bakta.db --hypotheticals hypotheticals.faa + mkdir -p work/tblout work/domtblout + nextflow run ${BAKTA_DB_SCRIPTS}/hmmsearch.nf --in hypotheticals.faa --db pfam-families --block 10000 --out hmmsearch.pfam-families.tblout + python3 ${BAKTA_DB_SCRIPTS}/annotate-pfam.py --db bakta.db --hmms pfam-families --hmm-results hmmsearch.pfam-families.tblout + # `pfam.h3*` (the non-family Pfam HMMs) are a runtime bakta asset and are + # kept; only the family HMMs used for this step's own hmmsearch, and the + # intermediate tsv/source files, are discarded here. + rm -f pfam-families* pfam pfam.families.tsv pfam.non-families.tsv Pfam-A.hmm.dat.gz Pfam-A.hmm hmmsearch.pfam-families.tblout hypotheticals.faa + touch .step18.done +fi ############################################################################ @@ -310,18 +565,38 @@ rm pfam-families* pfam *.tsv Pfam* hmmsearch.pfam-families.tblout # - import VFDB sequences ############################################################################ printf "\n19/19: download AA sequences for expert annotation system ...\n" -wget https://ftp.ncbi.nlm.nih.gov/pub/blastrules/4.2.2.tgz -tar -xzf 4.2.2.tgz -wget http://www.mgc.ac.cn/VFs/Down/VFDB_setA_pro.fas.gz -gunzip VFDB_setA_pro.fas.gz -python3 ${BAKTA_DB_SCRIPTS}/expert/setup-is.py --expert-sequence expert-protein-sequences.faa --proteins IS.faa -python3 ${BAKTA_DB_SCRIPTS}/expert/setup-ncbiblastrules.py --expert-sequence expert-protein-sequences.faa --ncbi-blastrule-tsv 4.2.2/data/blast-rules_4.2.2.tsv --proteins 4.2.2/data/proteins.fasta -python3 ${BAKTA_DB_SCRIPTS}/expert/setup-vfdb.py --expert-sequence expert-protein-sequences.faa --proteins VFDB_setA_pro.fas -diamond makedb --in expert-protein-sequences.faa --db expert-protein-sequences -rm -r 4.2.2/ 4.2.2.tgz IS.faa VFDB_setA_pro.fas expert-protein-sequences.faa +if have_outputs .step19.done; then + printf " already built, skipping.\n" +else + if [ ! -s 4.2.2.tgz ]; then + wget -c https://ftp.ncbi.nlm.nih.gov/pub/blastrules/4.2.2.tgz + fi + rm -rf 4.2.2 + tar -xzf 4.2.2.tgz + if [ ! -s VFDB_setA_pro.fas ]; then + if [ ! -s VFDB_setA_pro.fas.gz ]; then + wget -c http://www.mgc.ac.cn/VFs/Down/VFDB_setA_pro.fas.gz + fi + gunzip -f VFDB_setA_pro.fas.gz + fi + rm -f expert-protein-sequences.faa + python3 ${BAKTA_DB_SCRIPTS}/expert/setup-is.py --expert-sequence expert-protein-sequences.faa --proteins IS.faa + python3 ${BAKTA_DB_SCRIPTS}/expert/setup-ncbiblastrules.py --expert-sequence expert-protein-sequences.faa --ncbi-blastrule-tsv 4.2.2/data/blast-rules_4.2.2.tsv --proteins 4.2.2/data/proteins.fasta + python3 ${BAKTA_DB_SCRIPTS}/expert/setup-vfdb.py --expert-sequence expert-protein-sequences.faa --proteins VFDB_setA_pro.fas + diamond makedb --in expert-protein-sequences.faa --db expert-protein-sequences + rm -r 4.2.2/ 4.2.2.tgz IS.faa VFDB_setA_pro.fas expert-protein-sequences.faa + : "${BAKTA_ARCHAEAL_EXPERT_PROTEINS:?Set BAKTA_ARCHAEAL_EXPERT_PROTEINS to a curated archaeal protein FASTA.}" + diamond makedb --in "${BAKTA_ARCHAEAL_EXPERT_PROTEINS}" --db expert-protein-sequences-archaea + touch .step19.done +fi + +# The light database keeps its PSCC SQLite snapshot but requires the same +# profile-specific HMMs and expert databases as the full database. +cp antifam* ncRNA-genes* ncRNA-regions* rRNA* oric.fna orit.fna pfam* rfam-go.tsv expert-protein-sequences*.dmnd db-light/ # Cleanup ls -l bakta.db python3 ${BAKTA_DB_SCRIPTS}/optimize-db.py --db bakta.db --tmp /var/scratch/ ls -l bakta.db -rm psc.faa sorf.faa node.dmp +rm -f psc.faa sorf.faa nodes.dmp +touch .build.done diff --git a/db-scripts/init-psc.py b/db-scripts/init-psc.py index 29ed23e2..7dabd1e3 100644 --- a/db-scripts/init-psc.py +++ b/db-scripts/init-psc.py @@ -20,7 +20,9 @@ parser.add_argument('--db', action='store', help='Path to Bakta sqlite3 db file.') parser.add_argument('--psc', action='store', help='Path to PSC fasta file.') parser.add_argument('--psc_sorf', action='store', help='Path to sORF PSC fasta file.') +parser.add_argument('--taxon', action='append', default=None, help='NCBI root taxon to include; repeat for multiple domains (default: 2, Bacteria).') args = parser.parse_args() +taxa = args.taxon or ['2'] MAX_SORF_LENGTH = 30 @@ -98,7 +100,7 @@ def is_taxon_child(child, LCA, taxonomy): rep_member_tax_id = rep_member_dbref.find('./{*}property[@type="NCBI taxonomy"]') rep_member_tax_id = rep_member_tax_id.get('value') if rep_member_tax_id is not None else 1 - if(is_taxon_child(common_tax_id, '2', taxonomy) or is_taxon_child(rep_member_tax_id, '2', taxonomy) or 'phage' in rep_member_organism.lower()): + if(any(is_taxon_child(common_tax_id, taxon, taxonomy) or is_taxon_child(rep_member_tax_id, taxon, taxonomy) for taxon in taxa) or 'phage' in rep_member_organism.lower()): uniref90_id = elem.attrib['id'][9:] # remove 'UniRef90_' prefix product = rep_member_dbref.find('./{*}property[@type="protein name"]') diff --git a/db-scripts/init-pscc.py b/db-scripts/init-pscc.py index 82d48cf2..f5cacca5 100644 --- a/db-scripts/init-pscc.py +++ b/db-scripts/init-pscc.py @@ -20,7 +20,9 @@ parser.add_argument('--db', action='store', help='Path to Bakta sqlite3 db file.') parser.add_argument('--pscc', action='store', help='Path to PSCC fasta file.') parser.add_argument('--pscc_sorf', action='store', help='Path to sORF PSCC fasta file.') +parser.add_argument('--taxon', action='append', default=None, help='NCBI root taxon to include; repeat for multiple domains (default: 2, Bacteria).') args = parser.parse_args() +taxa = args.taxon or ['2'] PSCC_MIN_MEMBER_COUNT = 10 MAX_SORF_LENGTH = 30 @@ -99,7 +101,7 @@ def is_taxon_child(child, LCA, taxonomy): rep_member_organism = rep_member_organism.get('value') if rep_member_organism is not None else '' rep_member_tax_id = rep_member_dbref.find('./{*}property[@type="NCBI taxonomy"]') rep_member_tax_id = rep_member_tax_id.get('value') if rep_member_tax_id is not None else 1 - if(is_taxon_child(common_tax_id, '2', taxonomy) or is_taxon_child(rep_member_tax_id, '2', taxonomy) or 'phage' in rep_member_organism.lower()): + if(any(is_taxon_child(common_tax_id, taxon, taxonomy) or is_taxon_child(rep_member_tax_id, taxon, taxonomy) for taxon in taxa) or 'phage' in rep_member_organism.lower()): uniref50_id = elem.attrib['id'][9:] # remove 'UniRef50_' prefix product = rep_member_dbref.find('./{*}property[@type="protein name"]') if(product is not None): diff --git a/db-scripts/init-ups-ips.py b/db-scripts/init-ups-ips.py index 46db6c6c..c8f8b598 100644 --- a/db-scripts/init-ups-ips.py +++ b/db-scripts/init-ups-ips.py @@ -20,7 +20,9 @@ parser.add_argument('--uniparc', action='store', help='Path to UniParc fasta file.') parser.add_argument('--db', action='store', help='Path to Bakta sqlite3 db file.') parser.add_argument('--ips', action='store', help='Path to IPS fasta file.') +parser.add_argument('--taxon', action='append', default=None, help='NCBI root taxon to include; repeat for multiple domains (default: 2, Bacteria).') args = parser.parse_args() +taxa = args.taxon or ['2'] DISCARDED_PRODUCTS = [ 'hypothetical protein', @@ -98,7 +100,7 @@ def is_taxon_child(child, LCA, taxonomy): rep_member_tax_id = rep_member_tax_id.get('value') if rep_member_tax_id is not None else 1 # filter for bacterial or phage protein sequences - if(is_taxon_child(common_tax_id, '2', taxonomy) or is_taxon_child(rep_member_tax_id, '2', taxonomy) or 'phage' in rep_member_organism.lower()): + if(any(is_taxon_child(common_tax_id, taxon, taxonomy) or is_taxon_child(rep_member_tax_id, taxon, taxonomy) for taxon in taxa) or 'phage' in rep_member_organism.lower()): uniref100_id = elem.attrib['id'][10:] seq_representative = elem.find('./{*}representativeMember/{*}sequence') seq = seq_representative.text.upper() diff --git a/docs/archaeal-support-plan.md b/docs/archaeal-support-plan.md new file mode 100644 index 00000000..f3b2f03c --- /dev/null +++ b/docs/archaeal-support-plan.md @@ -0,0 +1,65 @@ +# Archaeal annotation support plan + +## Scope and compatibility + +Archaeal annotation is an explicit `--domain archaea` profile. `bacteria` +remains the default and must retain its current results and database layout. +One database build always creates both profiles; the profile is selected by the +user at runtime and is not inferred from sequence content. + +## Implemented + +- Domain-aware CLI configuration, JSON metadata, database checks, and protein + bulk-annotation configuration. +- Archaeal tRNA mode, rRNA models, and domain-selected Rfam ncRNA/regulatory + models. +- Domain-specific tmRNA labels, sORF initiation policy, and Sec/Pyl recoding + policy. +- Archaeal oriC prediction through Ori-Finder-Arch; bacterial oriC/oriT stays + on the existing BLAST/DoriC/MOB-suite path. +- UniRef/Swiss-Prot cluster builders accept multiple root taxa and include both + NCBI taxon `2` (Bacteria) and `2157` (Archaea). +- An online Ori-Finder-Arch installer plus MEME and HMMER environment support. +- Bacterial DoriC/MOB-suite, AMRFinderPlus, and ISfinder assets are retained + in the shared database, but are not used by the archaeal runtime profile. +- Ori-Finder-Arch provenance, RIP/origin attributes, scores, and tRNA intron + coordinates are retained in GFF and INSDC exports. + +## Remaining implementation work + +### Database and functional annotation + +1. Build, version, and publish dual-profile full and light database archives. + Required bacterial and archaeal assets are checked by `bakta.db`; no + production archive exists yet. +2. Curate a versioned archaeal expert-protein FASTA and supply it through + `BAKTA_ARCHAEAL_EXPERT_PROTEINS` during database construction. +3. Replace the bacterial AntiFam-only spurious-ORF filter with a validated + archaeal or domain-neutral model. +4. Review selected archaeal Rfam families and thresholds for precision; do not + assume that bacterial regulatory-RNA selection rules transfer unchanged. + +### Feature fidelity + +1. Validate tRNA intron coordinates and exports against curated archaeal + records. +2. Validate tmRNA calls and sORF start-codon policies across archaeal clades. +3. Validate Sec/Pyl recoding calls against known loci and add suitable Rfam + model coverage where required. +4. Validate Ori-Finder-Arch provenance, RIP type, and origin class in every + output format against downstream consumers. + +### Validation and release + +1. Add fixture databases and curated archaeal genomes for end-to-end tests. +2. Benchmark CDS boundaries, tRNAs, rRNAs, ncRNAs, pseudogenes, and oriCs; + report recall/precision by archaeal lineage and assembly quality. +3. Pin and checksum the upstream Ori-Finder-Arch binary release. +4. Update CLI/package wording, manuals, examples, and release notes. + +## Acceptance criteria + +The archaeal profile is release-ready only when a versioned database is +available, all profile dependencies are installable, end-to-end tests cover +complete and fragmented archaeal assemblies, and benchmark results meet the +project's agreed accuracy thresholds. diff --git a/environment.yml b/environment.yml index 6801d4eb..9c9b384f 100644 --- a/environment.yml +++ b/environment.yml @@ -16,7 +16,14 @@ dependencies: - infernal>=1.1.5 - piler-cr - pyhmmer>=0.12.0 + - hmmer>=3.4 + - meme>=5.5 - diamond==2.2.0 - blast>=2.17.0 - ncbi-amrfinderplus>=4.2.7 - pycirclize>=1.7.0 + # Test dependencies + - pytest>=8 + - pip + - pip: + - pytest-dependency>=0.6 diff --git a/scripts/install-ori-finder-arch.sh b/scripts/install-ori-finder-arch.sh new file mode 100755 index 00000000..29242c92 --- /dev/null +++ b/scripts/install-ori-finder-arch.sh @@ -0,0 +1,15 @@ +#!/usr/bin/env bash + +set -euo pipefail + +ORI_FINDER_ARCH_URL='https://tubic.tju.edu.cn/Ori-Finder-Arch/download/Ori-Finder-Arch.tar.gz' +ORI_FINDER_ARCH_SHA256='f6a29d446b82df93b28144c02d404eeee6a3ba83fd9e0a0f31f63ba6b606b703' +target_dir="${1:-${CONDA_PREFIX:?Provide an installation directory or activate the Conda environment}/bin}" +tmp_dir=$(mktemp -d) +trap 'rm -rf "${tmp_dir}"' EXIT + +mkdir -p "${target_dir}" +wget -q -O "${tmp_dir}/Ori-Finder-Arch.tar.gz" "${ORI_FINDER_ARCH_URL}" +echo "${ORI_FINDER_ARCH_SHA256} ${tmp_dir}/Ori-Finder-Arch.tar.gz" | sha256sum --check --status +tar -xzf "${tmp_dir}/Ori-Finder-Arch.tar.gz" -C "${tmp_dir}" +install -m 0755 "${tmp_dir}/OriFinderArch" "${target_dir}/OriFinderArch" diff --git a/setup.py b/setup.py index 9b52acb9..b253ea80 100644 --- a/setup.py +++ b/setup.py @@ -13,8 +13,8 @@ setup( name='bakta', version=bakta.__version__, - description='Bakta: rapid & standardized annotation of bacterial genomes, MAGs & plasmids', - keywords=['bioinformatics', 'annotation', 'bacteria', 'plasmids'], + description='Bakta: rapid & standardized annotation of prokaryotic genomes, MAGs & plasmids', + keywords=['bioinformatics', 'annotation', 'bacteria', 'archaea', 'plasmids'], long_description=long_description, long_description_content_type='text/markdown', license='GPLv3', diff --git a/test/test_organism_profiles.py b/test/test_organism_profiles.py new file mode 100644 index 00000000..bd08add6 --- /dev/null +++ b/test/test_organism_profiles.py @@ -0,0 +1,41 @@ +import sys + +import bakta.constants as bc +import bakta.config as cfg +import bakta.utils as bu +from bakta.features.r_rna import R_RNA_PROFILES +from bakta.features.ori_arch import _parse_attributes, _parse_orics + + +def test_archaeal_domain_argument_selects_archaeal_profile(monkeypatch): + monkeypatch.setattr(sys, 'argv', ['bakta', '--domain', 'archaea', 'genome.fna']) + assert bu.parse_arguments().domain == bc.DOMAIN_ARCHAEA + + +def test_archaeal_rrna_models_are_domain_specific(): + models = R_RNA_PROFILES[bc.DOMAIN_ARCHAEA]['models'] + assert 'RF01959' in models + assert 'RF02540' in models + assert 'RF00177' not in models + + +def test_archaeal_recoding_profile_includes_pyrrolysine(): + assert bc.DOMAIN_PROFILES[bc.DOMAIN_ARCHAEA]['recoding_codons']['TAG'] == 'pyrrolysine' + + +def test_archaeal_profile_ignores_gram_option(): + assert cfg.normalize_gram(bc.DOMAIN_ARCHAEA, bc.GRAM_POSITIVE) == bc.GRAM_UNKNOWN + assert cfg.normalize_gram(bc.DOMAIN_BACTERIA, bc.GRAM_POSITIVE) == bc.GRAM_POSITIVE + + +def test_ori_finder_arch_gff_attributes_are_preserved(): + assert _parse_attributes("type='Type I (RIP-adjacent)';RIP=CDC6") == {'type': 'Type I (RIP-adjacent)', 'RIP': 'CDC6'} + + +def test_ori_finder_arch_gff_is_mapped_to_bakta_feature(tmp_path): + result = tmp_path / 'ori.gff3' + result.write_text("seq1\tOri-Finder-Arch\trep_origin\t10\t20\t42.0\t.\t0\ttype='Type I';RIP=CDC6\n") + orics = _parse_orics(result, {'seq1': {'id': 'seq1', 'nt': 'A' * 100, 'length': 100}}) + assert orics[0]['type'] == bc.FEATURE_ORIC + assert orics[0]['score'] == 42.0 + assert orics[0]['attributes']['RIP'] == 'CDC6'