Deep Genomic Characterization of Bacillus subtilis SNUY1 & 168 Using MAYA (Multistage Analysis of Prokaryotic Annotations): An Integrated Snakemake-Based Bioinformatics Pipeline
- 1Amrita School of Biotechnology, Amrita Vishwa Vidyapeetham, Amritapuri, Kollam 690525, Kerala, India
- 2Sivasakthi Science Foundation, Trivandrum 695003, Kerala, India
Abstract
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.
Related
- TUDCA was identified as a potential therapeutic drug for nonalcoholic steatohepatitis complicated with hypercholesterolemia targeting ADRA2C by integrated bioinformatics and experimental verification
- Integrated bioinformatics analysis and mouse validation identify C7 and MFAP4 as candidate genes associated with type 2 diabetes and atherosclerosis
- Identification of Endoplasmic reticulum Homeostasis-associated diagnostic biomarkers for Sarcopenia through integrated Bioinformatics and Experimental Validation
- Computational design and molecular dynamics for developing a multistage serodiagnostic biomarker for tuberculosis
- Structure and mechanistic features of the prokaryotic minimal RNase P
BibTeX
@article{narendran2026deep,
title = {{Deep Genomic Characterization of Bacillus subtilis SNUY1 & 168 Using MAYA (Multistage Analysis of Prokaryotic Annotations): An Integrated Snakemake-Based Bioinformatics Pipeline}},
author = {Sastha Kumar Narendran and Anandhu Presannan and Lekshmi Indulekha and Abhijith Shibu and Adithiyakrishnan Rajeev and Satheesh Krishna Chakkadath and Vidhu S. Vijay and Pradeesh Babu and Geetha Kumar and Bipin G. Nair and Sabarinath Subramaniam and Nidheesh Melethadathil},
journal = {International Journal of Molecular Sciences},
year = {2026},
month = oct,
doi = {10.3390/ijms27209188},
url = {https://www.mdpi.com/1422-0067/27/20/9188},
note = {CC BY 4.0},
}
RIS
TY - JOUR
TI - Deep Genomic Characterization of Bacillus subtilis SNUY1 & 168 Using MAYA (Multistage Analysis of Prokaryotic Annotations): An Integrated Snakemake-Based Bioinformatics Pipeline
AU - Sastha Kumar Narendran
AU - Anandhu Presannan
AU - Lekshmi Indulekha
AU - Abhijith Shibu
AU - Adithiyakrishnan Rajeev
AU - Satheesh Krishna Chakkadath
AU - Vidhu S. Vijay
AU - Pradeesh Babu
AU - Geetha Kumar
AU - Bipin G. Nair
AU - Sabarinath Subramaniam
AU - Nidheesh Melethadathil
PY - 2026
DA - 2026/10/01
JO - International Journal of Molecular Sciences
DO - 10.3390/ijms27209188
UR - https://www.mdpi.com/1422-0067/27/20/9188
AB - 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.
SN - 1422-0067
ER -