#!/bin/bash
# ============================================================================
# GlaciaBase Metagenome Assembly & 16S Analysis Pipeline
# ============================================================================
# Author: Alba Tull
# Date: March 2026
# Purpose: Downloads sequencing data from GlaciaBase/NCBI SRA, classifies
#          runs as WGS (shotgun) or AMPLICON (16S), then routes each to
#          the appropriate analysis pipeline:
#            - WGS → fastp QC → MEGAHIT assembly → MetaBAT2 binning →
#                    CheckM2 quality → GTDB-Tk taxonomy
#            - AMPLICON → fastp QC → QIIME2/DADA2 denoising →
#                         SILVA taxonomy → diversity analysis
#
# Requirements:
#   - SRA Toolkit (prefetch, fasterq-dump)
#   - fastp
#   - MEGAHIT
#   - MetaBAT2
#   - samtools, minimap2 (for binning coverage)
#   - CheckM2
#   - GTDB-Tk (with GTDB database downloaded)
#   - QIIME2 (conda environment)
#
# Usage:
#   ./glaciabase_metagenome_pipeline.sh <input_file> <output_dir>
#
#   <input_file> is a TSV with columns: SRR_ACCESSION  BIOPROJECT  ASSAY_TYPE  ORGANISM
#   You can export this from GlaciaBase's Extract Data tab or manually create it.
#
#   Example input_file (runs.tsv):
#     SRR3454536	SRP074055	AMPLICON	soil metagenome
#     SRR3454537	SRP074055	WGS	soil metagenome
#     SRR7654321	PRJNA901488	WGS	marine metagenome
#
# ============================================================================

set -euo pipefail

# ── Configuration ──────────────────────────────────────────────────────────
THREADS=${THREADS:-8}                        # CPU threads (override with env var)
MIN_CONTIG_LEN=${MIN_CONTIG_LEN:-1000}       # Minimum contig length for assembly
CHECKM2_DB=${CHECKM2_DB:-""}                 # Path to CheckM2 database
GTDBTK_DB=${GTDBTK_DB:-""}                   # Path to GTDB-Tk database
SILVA_CLASSIFIER=${SILVA_CLASSIFIER:-""}     # Path to QIIME2 SILVA classifier .qza
QIIME2_ENV=${QIIME2_ENV:-"qiime2-2024.5"}   # Name of QIIME2 conda environment

# Colors for terminal output
RED='\033[0;31m'
GREEN='\033[0;32m'
BLUE='\033[0;34m'
YELLOW='\033[1;33m'
NC='\033[0m' # No Color

# ── Argument Parsing ───────────────────────────────────────────────────────
if [ $# -lt 2 ]; then
    echo -e "${RED}Usage: $0 <input_file.tsv> <output_dir>${NC}"
    echo ""
    echo "  input_file.tsv : Tab-separated file with columns:"
    echo "                   SRR_ACCESSION  BIOPROJECT  ASSAY_TYPE  ORGANISM"
    echo "  output_dir     : Directory for all output files"
    echo ""
    echo "  ASSAY_TYPE must be 'WGS' or 'AMPLICON'"
    echo ""
    echo "  Environment variables:"
    echo "    THREADS=8              Number of CPU threads"
    echo "    CHECKM2_DB=path        CheckM2 database path"
    echo "    GTDBTK_DB=path         GTDB-Tk database path"
    echo "    SILVA_CLASSIFIER=path  QIIME2 SILVA classifier .qza"
    echo "    QIIME2_ENV=name        Conda environment name for QIIME2"
    exit 1
fi

INPUT_FILE="$1"
OUTPUT_DIR="$2"

# ── Validation ─────────────────────────────────────────────────────────────
if [ ! -f "$INPUT_FILE" ]; then
    echo -e "${RED}Error: Input file '$INPUT_FILE' not found${NC}"
    exit 1
fi

# Check required tools
check_tool() {
    if ! command -v "$1" &> /dev/null; then
        echo -e "${YELLOW}Warning: '$1' not found in PATH. Install it before running the ${2:-$1} step.${NC}"
    else
        echo -e "${GREEN}  ✓ $1 found${NC}"
    fi
}

echo -e "${BLUE}Checking required tools...${NC}"
check_tool prefetch "SRA download"
check_tool fasterq-dump "SRA FASTQ extraction"
check_tool fastp "quality control"
check_tool megahit "metagenome assembly"
check_tool metabat2 "contig binning"
check_tool minimap2 "read mapping"
check_tool samtools "BAM processing"
check_tool checkm2 "MAG quality assessment"
check_tool gtdbtk "taxonomic classification"
echo ""

# ── Create Directory Structure ─────────────────────────────────────────────
echo -e "${BLUE}Creating output directories...${NC}"
mkdir -p "$OUTPUT_DIR"/{raw_fastq,qc_fastq,logs}
mkdir -p "$OUTPUT_DIR"/wgs/{assembly,binning,checkm2,gtdbtk}
mkdir -p "$OUTPUT_DIR"/amplicon/{import,denoise,taxonomy,diversity}
mkdir -p "$OUTPUT_DIR"/reports

# ── Classify Runs ──────────────────────────────────────────────────────────
# Split input file into WGS and AMPLICON run lists
WGS_RUNS="$OUTPUT_DIR/wgs_runs.tsv"
AMP_RUNS="$OUTPUT_DIR/amplicon_runs.tsv"

awk -F'\t' '$3 == "WGS" || $3 == "Metagenomics" || $3 == "WGA" || $3 == "WXS"' "$INPUT_FILE" > "$WGS_RUNS"
awk -F'\t' '$3 == "AMPLICON" || $3 == "16S" || $3 ~ /[Aa]mplicon/' "$INPUT_FILE" > "$AMP_RUNS"

WGS_COUNT=$(wc -l < "$WGS_RUNS")
AMP_COUNT=$(wc -l < "$AMP_RUNS")

echo -e "${GREEN}Run classification:${NC}"
echo -e "  WGS/shotgun runs:  ${WGS_COUNT}"
echo -e "  AMPLICON/16S runs: ${AMP_COUNT}"
echo ""

# ── Generate Annotation Report ─────────────────────────────────────────────
# This is the "annotate 16S vs shotgun" deliverable — a clear classification
# of every run with its data type, suitable for methods sections or PI review
echo -e "${BLUE}Generating annotation report...${NC}"

cat > "$OUTPUT_DIR/reports/run_annotation_report.tsv" << 'HEADER'
SRR_ACCESSION	BIOPROJECT	ASSAY_TYPE	DATA_CLASS	ORGANISM	ANALYSIS_PIPELINE	NOTES
HEADER

while IFS=$'\t' read -r srr bp assay organism rest; do
    case "$assay" in
        WGS|Metagenomics|WGA)
            echo -e "${srr}\t${bp}\t${assay}\tSHOTGUN\t${organism}\tMEGAHIT_assembly→MetaBAT2_binning→CheckM2→GTDB-Tk\tFull metagenome; functional + taxonomic analysis possible"
            ;;
        WXS)
            echo -e "${srr}\t${bp}\t${assay}\tSHOTGUN\t${organism}\tMEGAHIT_assembly→MetaBAT2_binning→CheckM2→GTDB-Tk\tExome capture; may have reduced genome coverage"
            ;;
        AMPLICON|16S)
            echo -e "${srr}\t${bp}\t${assay}\tAMPLICON_16S\t${organism}\tDADA2_denoise→SILVA_taxonomy→diversity_analysis\tTaxonomic census only; no functional gene content"
            ;;
        RNA-Seq)
            echo -e "${srr}\t${bp}\t${assay}\tTRANSCRIPTOME\t${organism}\tSKIPPED\tRNA-Seq requires separate metatranscriptome pipeline"
            ;;
        *)
            echo -e "${srr}\t${bp}\t${assay}\tUNKNOWN\t${organism}\tSKIPPED\tAssay type '${assay}' not recognized; manual review needed"
            ;;
    esac
done < "$INPUT_FILE" >> "$OUTPUT_DIR/reports/run_annotation_report.tsv"

echo -e "${GREEN}  Annotation report: $OUTPUT_DIR/reports/run_annotation_report.tsv${NC}"
echo ""


# ============================================================================
# STEP 1: DOWNLOAD FASTQ FILES
# ============================================================================
download_fastq() {
    local srr="$1"
    local outdir="$OUTPUT_DIR/raw_fastq"

    if [ -f "$outdir/${srr}_1.fastq" ] || [ -f "$outdir/${srr}_1.fastq.gz" ]; then
        echo -e "  ${YELLOW}Skipping $srr (already downloaded)${NC}"
        return 0
    fi

    echo -e "  Downloading $srr..."
    prefetch "$srr" --output-directory "$outdir" 2>"$OUTPUT_DIR/logs/${srr}_prefetch.log" || {
        echo -e "  ${RED}prefetch failed for $srr, trying fasterq-dump directly...${NC}"
    }
    fasterq-dump "$srr" --split-files -e "$THREADS" -O "$outdir" \
        2>"$OUTPUT_DIR/logs/${srr}_fasterq.log" || {
        echo -e "  ${RED}FAILED: Could not download $srr${NC}"
        return 1
    }
    echo -e "  ${GREEN}✓ $srr downloaded${NC}"
}

echo -e "${BLUE}═══ STEP 1: Downloading FASTQ files ═══${NC}"
while IFS=$'\t' read -r srr bp assay organism rest; do
    download_fastq "$srr"
done < "$INPUT_FILE"
echo ""


# ============================================================================
# STEP 2: QUALITY CONTROL (both WGS and AMPLICON)
# ============================================================================
run_qc() {
    local srr="$1"
    local indir="$OUTPUT_DIR/raw_fastq"
    local outdir="$OUTPUT_DIR/qc_fastq"

    if [ -f "$outdir/${srr}_1.fastq.gz" ]; then
        echo -e "  ${YELLOW}Skipping QC for $srr (already done)${NC}"
        return 0
    fi

    # Detect paired vs single-end
    if [ -f "$indir/${srr}_1.fastq" ] && [ -f "$indir/${srr}_2.fastq" ]; then
        echo -e "  QC (paired-end): $srr"
        fastp \
            -i "$indir/${srr}_1.fastq" -I "$indir/${srr}_2.fastq" \
            -o "$outdir/${srr}_1.fastq.gz" -O "$outdir/${srr}_2.fastq.gz" \
            --thread "$THREADS" \
            --detect_adapter_for_pe \
            --cut_front --cut_tail --cut_window_size 4 --cut_mean_quality 20 \
            --length_required 50 \
            --json "$OUTPUT_DIR/logs/${srr}_fastp.json" \
            --html "$OUTPUT_DIR/logs/${srr}_fastp.html" \
            2>"$OUTPUT_DIR/logs/${srr}_fastp.log"
    elif [ -f "$indir/${srr}.fastq" ]; then
        echo -e "  QC (single-end): $srr"
        fastp \
            -i "$indir/${srr}.fastq" \
            -o "$outdir/${srr}.fastq.gz" \
            --thread "$THREADS" \
            --cut_front --cut_tail --cut_window_size 4 --cut_mean_quality 20 \
            --length_required 50 \
            --json "$OUTPUT_DIR/logs/${srr}_fastp.json" \
            --html "$OUTPUT_DIR/logs/${srr}_fastp.html" \
            2>"$OUTPUT_DIR/logs/${srr}_fastp.log"
    else
        echo -e "  ${RED}No FASTQ found for $srr${NC}"
        return 1
    fi
    echo -e "  ${GREEN}✓ $srr QC complete${NC}"
}

echo -e "${BLUE}═══ STEP 2: Quality Control (fastp) ═══${NC}"
while IFS=$'\t' read -r srr bp assay organism rest; do
    run_qc "$srr"
done < "$INPUT_FILE"
echo ""


# ============================================================================
# STEP 3A: WGS METAGENOME ASSEMBLY (MEGAHIT)
# ============================================================================
if [ "$WGS_COUNT" -gt 0 ]; then
    echo -e "${BLUE}═══ STEP 3A: Metagenome Assembly (MEGAHIT) ═══${NC}"
    echo -e "  Assembling ${WGS_COUNT} WGS runs..."

    # Collect all QC'd WGS reads for co-assembly
    R1_FILES=""
    R2_FILES=""
    SE_FILES=""

    while IFS=$'\t' read -r srr bp assay organism rest; do
        if [ -f "$OUTPUT_DIR/qc_fastq/${srr}_1.fastq.gz" ]; then
            R1_FILES="${R1_FILES:+$R1_FILES,}$OUTPUT_DIR/qc_fastq/${srr}_1.fastq.gz"
            R2_FILES="${R2_FILES:+$R2_FILES,}$OUTPUT_DIR/qc_fastq/${srr}_2.fastq.gz"
        elif [ -f "$OUTPUT_DIR/qc_fastq/${srr}.fastq.gz" ]; then
            SE_FILES="${SE_FILES:+$SE_FILES,}$OUTPUT_DIR/qc_fastq/${srr}.fastq.gz"
        fi
    done < "$WGS_RUNS"

    # Build MEGAHIT command
    MEGAHIT_CMD="megahit -t $THREADS --min-contig-len $MIN_CONTIG_LEN -o $OUTPUT_DIR/wgs/assembly"
    if [ -n "$R1_FILES" ]; then
        MEGAHIT_CMD="$MEGAHIT_CMD -1 $R1_FILES -2 $R2_FILES"
    fi
    if [ -n "$SE_FILES" ]; then
        MEGAHIT_CMD="$MEGAHIT_CMD -r $SE_FILES"
    fi

    echo -e "  Running: megahit (this may take hours for large datasets)..."
    eval "$MEGAHIT_CMD" 2>"$OUTPUT_DIR/logs/megahit.log" || {
        echo -e "  ${RED}MEGAHIT assembly failed. Check $OUTPUT_DIR/logs/megahit.log${NC}"
    }

    CONTIGS="$OUTPUT_DIR/wgs/assembly/final.contigs.fa"
    if [ -f "$CONTIGS" ]; then
        CONTIG_COUNT=$(grep -c "^>" "$CONTIGS")
        echo -e "  ${GREEN}✓ Assembly complete: ${CONTIG_COUNT} contigs${NC}"
    fi
    echo ""


    # ════════════════════════════════════════════════════════════════════════
    # STEP 4A: CONTIG BINNING (MetaBAT2)
    # ════════════════════════════════════════════════════════════════════════
    echo -e "${BLUE}═══ STEP 4A: Contig Binning (MetaBAT2) ═══${NC}"

    if [ -f "$CONTIGS" ]; then
        # Map reads back to assembly for coverage information
        echo -e "  Mapping reads to contigs for coverage..."
        SORTED_BAM="$OUTPUT_DIR/wgs/binning/mapped_sorted.bam"

        if [ -n "$R1_FILES" ]; then
            # Use first pair of reads for mapping (or concatenate all)
            FIRST_R1=$(echo "$R1_FILES" | cut -d',' -f1)
            FIRST_R2=$(echo "$R2_FILES" | cut -d',' -f1)
            minimap2 -t "$THREADS" -a "$CONTIGS" "$FIRST_R1" "$FIRST_R2" 2>/dev/null | \
                samtools sort -@ "$THREADS" -o "$SORTED_BAM"
            samtools index "$SORTED_BAM"
        fi

        # Generate coverage depth file
        echo -e "  Calculating coverage depth..."
        jgi_summarize_bam_contig_depths --outputDepth "$OUTPUT_DIR/wgs/binning/depth.txt" "$SORTED_BAM" \
            2>"$OUTPUT_DIR/logs/depth.log"

        # Run MetaBAT2
        echo -e "  Running MetaBAT2..."
        metabat2 \
            -i "$CONTIGS" \
            -a "$OUTPUT_DIR/wgs/binning/depth.txt" \
            -o "$OUTPUT_DIR/wgs/binning/bins/bin" \
            -t "$THREADS" \
            --minContig 1500 \
            2>"$OUTPUT_DIR/logs/metabat2.log"

        BIN_COUNT=$(ls "$OUTPUT_DIR/wgs/binning/bins/"*.fa 2>/dev/null | wc -l)
        echo -e "  ${GREEN}✓ Binning complete: ${BIN_COUNT} bins (MAGs)${NC}"
    fi
    echo ""


    # ════════════════════════════════════════════════════════════════════════
    # STEP 5A: MAG QUALITY ASSESSMENT (CheckM2)
    # ════════════════════════════════════════════════════════════════════════
    echo -e "${BLUE}═══ STEP 5A: MAG Quality (CheckM2) ═══${NC}"

    if [ -n "$CHECKM2_DB" ] && [ "$BIN_COUNT" -gt 0 ]; then
        checkm2 predict \
            --input "$OUTPUT_DIR/wgs/binning/bins/" \
            --output-directory "$OUTPUT_DIR/wgs/checkm2" \
            --threads "$THREADS" \
            --database_path "$CHECKM2_DB" \
            2>"$OUTPUT_DIR/logs/checkm2.log"

        echo -e "  ${GREEN}✓ CheckM2 report: $OUTPUT_DIR/wgs/checkm2/quality_report.tsv${NC}"

        # Summarize: count high/medium/low quality MAGs by MIMAG standards
        if [ -f "$OUTPUT_DIR/wgs/checkm2/quality_report.tsv" ]; then
            echo -e "  MAG quality summary (MIMAG standards):"
            awk -F'\t' 'NR>1 {
                comp=$2; cont=$3;
                if (comp >= 90 && cont < 5) hq++;
                else if (comp >= 50 && cont < 10) mq++;
                else lq++;
            } END {
                printf "    High quality (>90%% complete, <5%% contamination): %d\n", hq;
                printf "    Medium quality (>50%% complete, <10%% contamination): %d\n", mq;
                printf "    Low quality: %d\n", lq;
            }' "$OUTPUT_DIR/wgs/checkm2/quality_report.tsv"
        fi
    else
        if [ -z "$CHECKM2_DB" ]; then
            echo -e "  ${YELLOW}Skipped: Set CHECKM2_DB environment variable to enable${NC}"
        fi
    fi
    echo ""


    # ════════════════════════════════════════════════════════════════════════
    # STEP 6A: TAXONOMIC CLASSIFICATION (GTDB-Tk)
    # ════════════════════════════════════════════════════════════════════════
    echo -e "${BLUE}═══ STEP 6A: Taxonomy (GTDB-Tk) ═══${NC}"

    if [ -n "$GTDBTK_DB" ] && [ "$BIN_COUNT" -gt 0 ]; then
        export GTDBTK_DATA_PATH="$GTDBTK_DB"
        gtdbtk classify_wf \
            --genome_dir "$OUTPUT_DIR/wgs/binning/bins/" \
            --out_dir "$OUTPUT_DIR/wgs/gtdbtk" \
            --cpus "$THREADS" \
            --extension fa \
            2>"$OUTPUT_DIR/logs/gtdbtk.log"

        echo -e "  ${GREEN}✓ GTDB-Tk results: $OUTPUT_DIR/wgs/gtdbtk/${NC}"
    else
        if [ -z "$GTDBTK_DB" ]; then
            echo -e "  ${YELLOW}Skipped: Set GTDBTK_DB environment variable to enable${NC}"
        fi
    fi
    echo ""
fi


# ============================================================================
# STEP 3B: 16S AMPLICON ANALYSIS (QIIME2 / DADA2)
# ============================================================================
if [ "$AMP_COUNT" -gt 0 ]; then
    echo -e "${BLUE}═══ STEP 3B: 16S Amplicon Analysis (QIIME2/DADA2) ═══${NC}"
    echo -e "  Processing ${AMP_COUNT} AMPLICON runs..."

    # Check if QIIME2 is available
    if conda env list 2>/dev/null | grep -q "$QIIME2_ENV"; then
        echo -e "  Activating QIIME2 environment: $QIIME2_ENV"

        # Create a manifest file for QIIME2 import
        MANIFEST="$OUTPUT_DIR/amplicon/import/manifest.tsv"
        echo -e "sample-id\tforward-read\treverse-read" > "$MANIFEST"

        while IFS=$'\t' read -r srr bp assay organism rest; do
            FWD="$OUTPUT_DIR/qc_fastq/${srr}_1.fastq.gz"
            REV="$OUTPUT_DIR/qc_fastq/${srr}_2.fastq.gz"
            if [ -f "$FWD" ] && [ -f "$REV" ]; then
                echo -e "${srr}\t${FWD}\t${REV}" >> "$MANIFEST"
            fi
        done < "$AMP_RUNS"

        # Run QIIME2 steps within the conda environment
        # Using conda run to avoid needing to source activate
        QIIME_CMD="conda run -n $QIIME2_ENV"

        # Import reads
        echo -e "  Importing reads into QIIME2..."
        $QIIME_CMD qiime tools import \
            --type 'SampleData[PairedEndSequencesWithQuality]' \
            --input-path "$MANIFEST" \
            --output-path "$OUTPUT_DIR/amplicon/import/demux.qza" \
            --input-format PairedEndFastqManifestPhred33V2 \
            2>"$OUTPUT_DIR/logs/qiime2_import.log"

        # DADA2 denoising
        echo -e "  Running DADA2 denoising..."
        $QIIME_CMD qiime dada2 denoise-paired \
            --i-demultiplexed-seqs "$OUTPUT_DIR/amplicon/import/demux.qza" \
            --p-trunc-len-f 250 \
            --p-trunc-len-r 200 \
            --p-n-threads "$THREADS" \
            --o-table "$OUTPUT_DIR/amplicon/denoise/table.qza" \
            --o-representative-sequences "$OUTPUT_DIR/amplicon/denoise/rep-seqs.qza" \
            --o-denoising-stats "$OUTPUT_DIR/amplicon/denoise/stats.qza" \
            2>"$OUTPUT_DIR/logs/qiime2_dada2.log"

        echo -e "  ${GREEN}✓ DADA2 denoising complete${NC}"

        # Taxonomy assignment with SILVA
        if [ -n "$SILVA_CLASSIFIER" ] && [ -f "$SILVA_CLASSIFIER" ]; then
            echo -e "  Assigning taxonomy with SILVA..."
            $QIIME_CMD qiime feature-classifier classify-sklearn \
                --i-classifier "$SILVA_CLASSIFIER" \
                --i-reads "$OUTPUT_DIR/amplicon/denoise/rep-seqs.qza" \
                --o-classification "$OUTPUT_DIR/amplicon/taxonomy/taxonomy.qza" \
                --p-n-jobs "$THREADS" \
                2>"$OUTPUT_DIR/logs/qiime2_taxonomy.log"

            # Export taxonomy to TSV for easy viewing
            $QIIME_CMD qiime tools export \
                --input-path "$OUTPUT_DIR/amplicon/taxonomy/taxonomy.qza" \
                --output-path "$OUTPUT_DIR/amplicon/taxonomy/"

            echo -e "  ${GREEN}✓ Taxonomy: $OUTPUT_DIR/amplicon/taxonomy/taxonomy.tsv${NC}"
        else
            echo -e "  ${YELLOW}Taxonomy skipped: Set SILVA_CLASSIFIER to your .qza classifier path${NC}"
        fi

        # Alpha and beta diversity
        echo -e "  Computing diversity metrics..."
        $QIIME_CMD qiime diversity core-metrics-phylogenetic \
            --i-table "$OUTPUT_DIR/amplicon/denoise/table.qza" \
            --i-phylogeny <($QIIME_CMD qiime phylogeny align-to-tree-mafft-fasttree \
                --i-sequences "$OUTPUT_DIR/amplicon/denoise/rep-seqs.qza" \
                --o-alignment /dev/null \
                --o-masked-alignment /dev/null \
                --o-tree /dev/stdout \
                --o-rooted-tree /dev/null 2>/dev/null) \
            --p-sampling-depth 1000 \
            --output-dir "$OUTPUT_DIR/amplicon/diversity/" \
            2>"$OUTPUT_DIR/logs/qiime2_diversity.log" || {
            echo -e "  ${YELLOW}Diversity analysis requires >1 sample; skipped for single-sample runs${NC}"
        }

    else
        echo -e "  ${YELLOW}QIIME2 not found. Falling back to standalone DADA2 approach.${NC}"
        echo -e "  To install: conda create -n qiime2-2024.5 --channel https://packages.qiime2.org/qiime2/2024.5/shotgun/released qiime2"
        echo ""
        echo -e "  Alternative: Run 16S analysis in R with DADA2 package:"
        echo -e "    library(dada2)"
        echo -e "    # See: https://benjjneb.github.io/dada2/tutorial.html"
    fi
    echo ""
fi


# ============================================================================
# FINAL REPORT
# ============================================================================
echo -e "${BLUE}═══ PIPELINE COMPLETE ═══${NC}"
echo ""
echo -e "Output directory: $OUTPUT_DIR"
echo ""
echo "Directory structure:"
echo "  $OUTPUT_DIR/"
echo "  ├── raw_fastq/           # Downloaded FASTQ files"
echo "  ├── qc_fastq/            # Quality-trimmed reads (fastp)"
echo "  ├── logs/                 # All tool logs + fastp HTML reports"
echo "  ├── wgs_runs.tsv          # Classified WGS runs"
echo "  ├── amplicon_runs.tsv     # Classified AMPLICON runs"
echo "  ├── wgs/"
echo "  │   ├── assembly/         # MEGAHIT contigs (final.contigs.fa)"
echo "  │   ├── binning/bins/     # MetaBAT2 MAGs (.fa files)"
echo "  │   ├── checkm2/          # MAG quality report"
echo "  │   └── gtdbtk/           # Taxonomic classification"
echo "  ├── amplicon/"
echo "  │   ├── import/           # QIIME2 imported reads"
echo "  │   ├── denoise/          # DADA2 ASV table + rep sequences"
echo "  │   ├── taxonomy/         # SILVA taxonomy assignments"
echo "  │   └── diversity/        # Alpha + beta diversity results"
echo "  └── reports/"
echo "      └── run_annotation_report.tsv  # 16S vs shotgun classification"
echo ""
echo -e "${GREEN}Key output files:${NC}"
echo "  Annotation:  $OUTPUT_DIR/reports/run_annotation_report.tsv"
if [ "$WGS_COUNT" -gt 0 ]; then
    echo "  Contigs:     $OUTPUT_DIR/wgs/assembly/final.contigs.fa"
    echo "  MAG bins:    $OUTPUT_DIR/wgs/binning/bins/*.fa"
    echo "  MAG quality: $OUTPUT_DIR/wgs/checkm2/quality_report.tsv"
    echo "  Taxonomy:    $OUTPUT_DIR/wgs/gtdbtk/"
fi
if [ "$AMP_COUNT" -gt 0 ]; then
    echo "  ASV table:   $OUTPUT_DIR/amplicon/denoise/table.qza"
    echo "  16S taxonomy:$OUTPUT_DIR/amplicon/taxonomy/taxonomy.tsv"
fi
echo ""
echo -e "${GREEN}Done.${NC}"