This repository contains a fully reproducible computational pipeline for the analysis of soil microbial community responses to escalating antibiotic concentrations. The study employs a dual-marker amplicon sequencing approach, simultaneously profiling bacterial communities via the 16S rRNA gene (V3-V4 region) and fungal communities via the internal transcribed spacer region (ITS2), to characterize the taxonomic and ecological consequences of antibiotic contamination in agricultural soils planted with two crop species.
All raw sequencing data are publicly available from the European Nucleotide Archive under BioProject accession PRJEB95045.
The experimental framework follows a fully balanced 4 x 2 factorial design comprising 48 biological samples per marker gene (96 libraries total). Soil samples were collected from pots planted with Spinacia oleracea (spinach) or Raphanus sativus (radish) and amended with a mixture of five antibiotics at four concentration levels. Six biological replicates were maintained per treatment group.
| Factor | Levels |
|---|---|
| Plant species | Spinacia oleracea (spinach), Raphanus sativus (radish) |
| Antibiotic concentration | 0 (c0), 0.1 (c0.1), 1.0 (c1), 10.0 (c10) mg/kg soil dw |
| Biological replicates | 6 per treatment combination |
| Sequencing platform | Illumina NextSeq 2000, 2 x 301 bp paired-end |
The complete analytical workflow was implemented in QIIME 2, with computationally intensive steps executed via standalone shell scripts to manage memory constraints. The pipeline is organized into two parallel tracks, one for each marker gene, each documented in a dedicated Jupyter notebook and supported by modular execution scripts.
- Data Importation and Quality Assessment: Raw paired-end FASTQ files were imported using a fastq manifest and profiled for per-base sequence quality to empirically determine truncation parameters.
- Amplicon Sequence Variant (ASV) Inference: DADA2 was employed for denoising with the following empirically optimized parameters: forward reads truncated at 280 bp, reverse reads at 240 bp, with primer sequences removed by positional trimming (17 bp forward, 21 bp reverse). These parameters yield an effective overlap of approximately 22 bp for the V3-V4 insert, exceeding the 12 bp algorithmic minimum. Thread count was restricted to 2 to prevent out-of-memory failures on systems with 16 GB RAM.
- Taxonomic Classification: Taxonomic assignments were performed using a pre-trained Naive Bayes classifier against the SILVA 138.2 SSU NR99 reference database. Region-specific reference sequences (V3-V4, primers 341F/806R) were extracted via the RESCRIPt plugin to maximize classification accuracy.
- Phylogenetic Reconstruction: A prevalence filter (minimum 2 samples) was applied to remove transient singleton ASVs prior to multiple sequence alignment (MAFFT) and approximate maximum-likelihood tree inference (FastTree). This yielded a midpoint-rooted phylogenetic tree suitable for computing phylogenetically informed diversity metrics.
- Diversity and Statistical Analysis: Following contaminant excision (mitochondria and chloroplast sequences), core diversity metrics were computed at a rarefaction depth of 100,000 reads, retaining all 48 samples. Alpha diversity was assessed using Shannon entropy, Pielou evenness, observed features, and Faith phylogenetic diversity, with group significance evaluated via Kruskal-Wallis tests. Beta diversity was quantified using weighted and unweighted UniFrac, Bray-Curtis dissimilarity, and Jaccard distance, with treatment effects tested via PERMANOVA (999 permutations). Differential abundance analysis was performed at the genus level using ANCOM-BC.
- Primer Removal: Unlike the 16S pipeline, ITS primer sequences (ITS3ngs/ITS4ngs; Tedersoo et al., 2014) were removed using Cutadapt's adaptive alignment algorithm rather than positional trimming, as the conserved flanking regions exhibit length variation. The
--p-discard-untrimmedflag ensured that only reads containing a valid ITS primer pair were retained, effectively filtering non-target reads from the multiplexed library. - DADA2 Denoising: No fixed-length truncation was applied (
--p-trunc-len-f 0,--p-trunc-len-r 0) to preserve the natural length heterogeneity of the ITS2 region (~200 to 450 bp). Quality filtering relied exclusively on the maximum expected errors threshold (2.0). - Taxonomic Classification: ITS sequences were classified against the UNITE 10.0 reference database (release 2024-04-04, dynamic species hypothesis clustering, all eukaryotes) using the VSEARCH consensus taxonomy method. This alignment-based approach was selected over the Naive Bayes classifier for its reduced memory footprint while maintaining comparable classification accuracy.
- Non-Phylogenetic Diversity Analysis: Because the ITS region cannot be meaningfully aligned across distantly related fungal taxa, all diversity analyses were restricted to non-phylogenetic metrics (Bray-Curtis, Jaccard, Shannon, observed features, Pielou evenness). A rarefaction depth of 4,000 reads was selected, retaining 44 of 48 samples (92%). Alpha and beta group significance testing and ANCOM-BC differential abundance analysis were performed identically to the 16S pipeline.
qiime2_antibiotic_microbiome_analysis/
|-- README.md # This document
|-- license # MIT License
|-- .gitignore # Version control exclusion rules
|-- installation.sh # QIIME 2 environment setup script
|-- environment.yml # Conda environment specification
|-- data_download.ipynb # Notebook: raw data retrieval from ENA
|-- PRJEB95045_runinfo.tsv # ENA run accession metadata (96 runs)
|-- 01_metadata/ # Sample metadata and manifests
| |-- 16S_sample_metadata.tsv # 16S factorial design metadata
| |-- 16S_manifest.tsv # 16S QIIME 2 import manifest
| |-- ITS_sample_metadata.tsv # ITS factorial design metadata
| |-- ITS_manifest.tsv # ITS QIIME 2 import manifest
|-- 16S_raw/ # Raw 16S paired-end FASTQ files (48 samples)
| |-- SB_16S_*_1.fastq.gz # [IGNORED: re-download via data_download.ipynb]
| |-- SB_16S_*_2.fastq.gz
|-- ITS_raw/ # Raw ITS paired-end FASTQ files (48 samples)
| |-- SB_ITS_*_1.fastq.gz # [IGNORED: re-download via data_download.ipynb]
| |-- SB_ITS_*_2.fastq.gz
|-- 02_interim_artifacts/ # QIIME 2 intermediate processing artifacts
| |-- 16S/ # DADA2 outputs, filtered tables, summaries
| |-- ITS/ # DADA2 outputs, trimmed demux, summaries
|-- 03_phylogeny_and_taxonomy/ # Reference databases, trees, taxonomic assignments
| |-- 16S/ # SILVA 138.2, phylogenetic tree, 16S taxonomy
| | |-- tree_export/ # Exported Newick tree file
| |-- ITS/ # UNITE 10.0, ITS taxonomy, VSEARCH results
|-- 04_diversity_and_stats/ # Diversity metrics and statistical test results
| |-- 16S/ # Core metrics, alpha/beta significance, ANCOM-BC
| | |-- 16S-core-metrics/ # Rarefied table, distance matrices, PCoA
| |-- ITS/ # Core metrics, alpha/beta significance, ANCOM-BC
| | |-- ITS-core-metrics/ # Rarefied table, distance matrices, PCoA
|-- notebooks/ # Jupyter notebook documentation
| |-- 16S_pipeline/ # 16S analysis notebook and local artifacts
| |-- ITS_pipeline/ # ITS analysis notebook
|-- run_dada2_16s.sh # Standalone DADA2 denoising script
|-- run_phylogeny_16s.sh # Standalone phylogeny construction script
|-- run_taxonomy_16s.sh # Standalone taxonomic classification script
|-- run_diversity_16s.sh # Standalone diversity analysis script
|-- run_its_taxonomy.sh # Standalone ITS taxonomy pipeline script
| Software | Version | Purpose |
|---|---|---|
| QIIME 2 | 2026.1 | Core amplicon analysis framework |
| DADA2 | (bundled) | Amplicon sequence variant inference |
| MAFFT | (bundled) | Multiple sequence alignment |
| FastTree | (bundled) | Approximate maximum-likelihood phylogeny |
| Cutadapt | (bundled) | Adaptive primer removal (ITS) |
| RESCRIPt | (bundled) | Reference database acquisition |
| VSEARCH | (bundled) | Consensus taxonomic classification (ITS) |
| scikit-learn | 1.4.2 | Naive Bayes classifier (16S) |
| Python | 3.10+ | Scripting and notebook execution |
| pandas | (bundled) | Metadata manipulation |
The pipeline was developed and validated on a system with 16 GB RAM. Memory-intensive steps (DADA2 denoising, MAFFT alignment, Naive Bayes classification) are configured with thread-limited execution (--p-n-threads 2 or --p-n-jobs 2) to prevent out-of-memory failures. Approximately 15 GB of disk space is required for the complete set of intermediate and output artifacts (excluding raw FASTQ files).
# Clone the repository
git clone [https://github.com/Qasim-Hussain-Code/qiime2_antibiotic_microbiome_analysis.git](https://github.com/Qasim-Hussain-Code/qiime2_antibiotic_microbiome_analysis.git)
cd qiime2_antibiotic_microbiome_analysis
# Install the QIIME 2 conda environment
bash installation.sh
# Activate the environment
conda activate qiime2-amplicon-2026.1
# Verify the installation
qiime info
Raw sequencing data must be downloaded from the European Nucleotide Archive prior to analysis. The data_download.ipynb notebook automates this process using the ENA REST API. Alternatively, the raw FASTQ files can be retrieved manually using the accession numbers listed in PRJEB95045_runinfo.tsv.
The 16S pipeline is documented in notebooks/16S_pipeline/01_16S_rRNA_amplicon_analysis.ipynb. Due to the memory requirements of several processing steps, the following standalone scripts should be executed sequentially:
# DADA2 denoising (~45 minutes)
nohup bash run_dada2_16s.sh > dada2_16s.log 2>&1 &
# Phylogenetic tree construction (~30 minutes)
nohup bash run_phylogeny_16s.sh > phylogeny_16s.log 2>&1 &
# Taxonomic classification (~20 minutes)
nohup bash run_taxonomy_16s.sh > taxonomy_16s.log 2>&1 &
# Diversity analysis and statistical testing (~15 minutes)
nohup bash run_diversity_16s.sh > diversity_16s.log 2>&1 &
The ITS pipeline is documented in notebooks/ITS_pipeline/01_ITS_fungal_amplicon_analysis.ipynb. The notebook is designed to be executed interactively, with computationally intensive taxonomy steps optionally offloaded to:
nohup bash run_its_taxonomy.sh > its_taxonomy.log 2>&1 &
| Database | Version | Marker | Classification Method | Source |
|---|---|---|---|---|
| SILVA | 138.2 NR99 | 16S | Naive Bayes (scikit-learn) | SILVA Project |
| UNITE | 10.0 Dynamic | ITS | VSEARCH consensus taxonomy | UNITE Community |
Reference databases are acquired programmatically via the QIIME 2 RESCRIPt plugin during pipeline execution. Pre-trained classifiers and extracted reference sequences exceeding 100 MB are excluded from version control and are regenerated automatically by the pipeline scripts.
All analytical outputs are stored as QIIME 2 artifacts (.qza) and visualizations (.qzv). Visualization files can be inspected interactively using the QIIME 2 Viewer at https://view.qiime2.org without requiring a local QIIME 2 installation.
| Artifact | Description |
|---|---|
04_diversity_and_stats/16S/16S-core-metrics/ |
Rarefied tables, distance matrices, PCoA ordinations |
04_diversity_and_stats/16S/16S-alpha-*-significance.qzv |
Kruskal-Wallis alpha diversity tests |
04_diversity_and_stats/16S/16S-beta-*-*.qzv |
PERMANOVA beta diversity significance tests |
04_diversity_and_stats/16S/16S-ancombc-*.qzv |
ANCOM-BC differential abundance barplots |
04_diversity_and_stats/16S/16S-emperor-*.qzv |
Interactive PCoA Emperor plots |
04_diversity_and_stats/ITS/ITS-core-metrics/ |
Non-phylogenetic core diversity outputs |
04_diversity_and_stats/ITS/ITS-alpha-*-significance.qzv |
Kruskal-Wallis alpha diversity tests |
04_diversity_and_stats/ITS/ITS-beta-*-*.qzv |
PERMANOVA beta diversity significance tests |
04_diversity_and_stats/ITS/ITS-ancombc-*.qzv |
ANCOM-BC differential abundance barplots |
03_phylogeny_and_taxonomy/16S/16S-taxonomy.qzv |
16S taxonomic classification results |
03_phylogeny_and_taxonomy/ITS/ITS-taxonomy.qzv |
ITS taxonomic classification results |
03_phylogeny_and_taxonomy/ITS/ITS-taxa-barplot.qzv |
Interactive taxonomic composition barplot |
GitHub imposes a 100 MB per-file size limit. Files exceeding this threshold, including raw sequencing reads, demultiplexed QIIME 2 artifacts, and reference database files, are excluded from version control via .gitignore. Placeholder README.md files are placed in each directory containing ignored content to preserve the repository structure and document the expected files, their approximate sizes, and instructions for regeneration.
The following files are excluded (totaling approximately 6.5 GB):
| File | Size | Regeneration Method |
|---|---|---|
16S_raw/*.fastq.gz |
~2.4 GB | Re-download via data_download.ipynb |
ITS_raw/*.fastq.gz |
~3.1 GB | Re-download via data_download.ipynb |
02_interim_artifacts/16S/16S-paired-end-demux.qza |
~2.4 GB | Re-import via 16S notebook Section 3 |
02_interim_artifacts/ITS/ITS-paired-end-demux.qza |
~3.1 GB | Re-import via ITS notebook Section 3 |
02_interim_artifacts/ITS/ITS-trimmed-demux.qza |
~398 MB | Re-run ITS notebook Section 4 |
03_phylogeny_and_taxonomy/16S/silva-138.2-seqs.qza |
~193 MB | Re-run run_taxonomy_16s.sh Step 1 |
03_phylogeny_and_taxonomy/16S/silva-138.2-dna-seqs.qza |
~174 MB | Re-run run_taxonomy_16s.sh Step 1.5 |
03_phylogeny_and_taxonomy/16S/silva-138-99-nb-classifier.qza |
~209 MB | Re-download via run_taxonomy_16s.sh Step 3 |
Raw sequencing data are deposited in the European Nucleotide Archive (ENA) under BioProject accession PRJEB95045. The dataset comprises 96 paired-end sequencing runs (48 for 16S rRNA, 48 for ITS) generated on an Illumina NextSeq 2000 platform. Sample accessions, download URLs, and MD5 checksums are provided in PRJEB95045_runinfo.tsv.
If you use this pipeline or adapt any component of this workflow, please cite the original study:
Soil microbial and plant responses to increasing antibiotic concentration: a case study of five antibiotics. European Nucleotide Archive. (BioProject PRJEB95045)
This project is licensed under the MIT License. See the license file for details.