Deep Genomic Characterization of Bacillus subtilis SNUY1 & 168 Using MAYA (Multistage Analysis of Prokaryotic Annotations): An Integrated Snakemake-Based Bioinformatics Pipeline

Bacillus subtilis is an extensively characterized Gram-positive bacterium with broad applications in agriculture, biotechnology, and biocontrol. Strain SNUY1, known for its antifungal properties and use as an edible coating, has been completely sequenced and deposited in NCBI (accession: CP199918), making it an important model for genomic and functional characterization. Here, we present the MAYA pipeline, a fully automated, rule-based Snakemake workflow that integrates eighteen bioinformatic tools into a single end-to-end framework. MAYA accepts either raw Illumina paired-end FASTQ reads or pre-assembled FASTA sequences, utilizing per-rule Conda environment isolation to ensure strict version control and reproducibility across all constituent tools from quality control and assembly through functional annotation, plasmid typing, and detection of biosynthetic and antimicrobial resistance gene clusters. To benchmark pipeline accuracy across both input modes, MAYA was evaluated using Bacillus subtilis subsp. subtilis str. 168(AL009126 OZ025638). The FASTA assembly yielded a circularized 4,215,606 bp chromosome (N50: 4.21 Mb, 43.51% GC) with 98.28% CheckM completeness and 100.0% BUSCO ortholog recovery. Prokka annotation predicted 4214 CDSs (86 tRNAs, 30 rRNAs) in FASTA mode compared to 4313 CDSs (106 tRNAs, 25 rRNAs) in FASTQ mode, reflecting expected read fragmentation in unassembled data. DIAMOND mapping against SwissProt yielded 4135 hits (99.90% mean identity) in FASTA mode and 4170 hits (99.75% mean identity) in FASTQ mode. Both modes identified identical ST1 Multilocus Sequence Typing (MLST) profiles and six core AMR determinants (mphK, vmlR, rphC, aadK, satA_Bsub, and tet(L)), while antiSMASH detected 15 biosynthetic gene cluster (BGC) regions (22 protoclusters) in FASTA mode versus 16 regions (23 protoclusters) in FASTQ mode. MOB-suite resolved two non-mobilizable plasmids in the FASTQ assembly (bioproject: PRJNA76): pAC670 and a novel 1666 bp plasmid. Furthermore, nucleotide BLASTN analysis resolved a 76-CDS discrepancy between published NCBI annotations and MAYA/Prokka, demonstrating 95.5% (4081/4272) stringent sequence match and confirming that differences stem from gene-model predictions rather than omitted loci. To demonstrate the FASTA-mode application workflow, the complete SNUY1 (NZ_CP199918) genome assembly was processed through MAYA, yielding a single, contiguous chromosome of 4,196,931 base pairs with a GC content of 43.69% (98.28% CheckM, 100% BUSCO completeness). Annotation identified 4196 coding sequences, while antiSMASH profiling revealed 13 candidate biosynthetic gene cluster (BGC) regions comprising 20 protoclusters (including NRPS, lanthipeptides, terpenes, and sactipeptides), aligning with the strain’s biocontrol capabilities. AMR profiling identified six chromosomally encoded resistance genes (mphK, vmlR, rphC, aadK, satA_Bsub, and tet(L)). In silico mobility analysis via MOB-suite demonstrated that these sequences lack functional conjugation machinery and are classified as non-mobilizable. Executing in 3–4 h on an Intel Xeon server (128 GB RAM, 29 CPU cores), MAYA delivers an efficient, audit-ready framework for end-to-end bacterial genomic profiling.