Taxonomic Classification with Kraken2

conda create -n kraken2_env -c bioconda -c conda-forge \
    kraken2 bracken krakentools krona
conda activate kraken2_env

# Verify
kraken2 --version
bracken --version
ktImportTaxonomy --help | head -3

mkdir -p kraken2_db

# Step 1: Download NCBI taxonomy (required for all databases)
kraken2-build --download-taxonomy --db kraken2_db

# IMPORTANT: The script uses ftp.ncbi.nih.gov
# If your cluster blocks FTP, edit kraken2-build's download_taxonomy.sh:
# Change ftp:// to https://
# Example fix:
DB_DIR="${KRAKEN2_DEFAULT_DB:-kraken2_db}"
sed -i 's|ftp://ftp.ncbi.nih.gov|https://ftp.ncbi.nih.gov|g' \
    $(which kraken2-build)   # edit the installed script if needed

# Step 2: Add standard bacterial library
kraken2-build --download-library bacteria --db kraken2_db

# Step 3: Add viral library
kraken2-build --download-library viral --db kraken2_db

# Step 4: Add human library (to remove host contamination)
kraken2-build --download-library human --db kraken2_db

echo "Standard libraries downloaded. Directory size:"
du -sh kraken2_db/library/

mkdir -p custom_genomes

# Download specific bacterial genomes associated with nematodes
# Example: Wolbachia (common nematode endosymbiont)
# Taxon ID must be added to FASTA headers as: |kraken:taxid|XXXX

# Download Wolbachia genome (TaxID: 953)
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/008/025/GCF_000008025.1_ASM802v1/GCF_000008025.1_ASM802v1_genomic.fna.gz
gzip -d GCF_000008025.1_ASM802v1_genomic.fna.gz

# Add taxonomy ID to FASTA headers
# Format: >seq_name|kraken:taxid|953
sed 's/^>/&|kraken:taxid|953 /' GCF_000008025.1_ASM802v1_genomic.fna \
    > custom_genomes/Wolbachia_taxid953.fna

# Repeat for other organisms:
# Xanthomonas (TaxID: 338) - common soil bacterium
# Bacillus thuringiensis (TaxID: 1428) - nematode pathogen
# Microbotryum (TaxID: 5936) - plant fungal pathogen

# Add all custom genomes to the Kraken2 library
for genome in custom_genomes/*.fna; do
    echo "Adding: $genome"
    kraken2-build --add-to-library ${genome} --db kraken2_db
done

# Build the final database index
# --kmer-len 35  : k-mer length (default 35 for Kraken2)
# --minimizer-len 31 : minimizer length
# --minimizer-spaces 6 : spaces in minimizer (reduces memory)
kraken2-build --build \
    --db kraken2_db \
    --kmer-len 35 \
    --minimizer-len 31 \
    --minimizer-spaces 6 \
    --threads 24

echo "Database build complete!"
kraken2-build --db kraken2_db --report-minimizer-data

conda activate kraken2_env

DB="kraken2_db"
OUTDIR="kraken2_results"
mkdir -p ${OUTDIR}

# Classify single-end RNA-seq reads
for fq in data/*.fastq.gz; do
    SAMPLE=$(basename ${fq} .fastq.gz)
    echo "Classifying: ${SAMPLE}"

    kraken2 \
        --db ${DB} \
        --threads 16 \
        --output ${OUTDIR}/${SAMPLE}.kraken \
        --report ${OUTDIR}/${SAMPLE}.report \
        --gzip-compressed \
        ${fq}

    echo "  Report: ${OUTDIR}/${SAMPLE}.report"
done

# For paired-end reads:
# kraken2 --db ${DB} --threads 16 \
#   --paired reads_R1.fastq.gz reads_R2.fastq.gz \
#   --output sample.kraken --report sample.report

# Show top classifications from one sample
echo "Top 20 classifications:"
sort -rn -k3,3 ${OUTDIR}/${SAMPLE}.report | head -20

# Report format columns:
# 1: % reads covered by clade
# 2: # reads covered by clade (direct + children)
# 3: # reads directly assigned to taxon
# 4: taxonomic rank code (U,R,D,P,C,O,F,G,S)
# 5: taxonomy ID
# 6: scientific name

# Step 1: Build Bracken database (once per Kraken2 database)
bracken-build \
    -d kraken2_db \
    -t 16 \
    -l 150 \
    -k 35
# -l = read length (match your actual read length)
# -k = k-mer length used for Kraken2 database

# Step 2: Run Bracken on each sample report
for report in kraken2_results/*.report; do
    SAMPLE=$(basename ${report} .report)
    echo "Bracken: ${SAMPLE}"

    bracken \
        -d kraken2_db \
        -i ${report} \
        -o bracken_results/${SAMPLE}.bracken \
        -r 150 \
        -l S \
        -t 10
    # -l S = species level
    # -t 10 = minimum read threshold for re-estimation
done

# Combine all Bracken results into a single count table
python3 $(which combine_bracken_outputs.py) \
    --files bracken_results/*.bracken \
    --output combined_bracken.txt

echo "Combined table: combined_bracken.txt"
head -5 combined_bracken.txt

mkdir -p krona_plots

# Convert Kraken2 report format to Krona input
for report in kraken2_results/*.report; do
    SAMPLE=$(basename ${report} .report)

    # Convert Kraken2 report to Krona format
    kreport2krona.py \
        -r ${report} \
        -o krona_plots/${SAMPLE}.krona.txt

    # Generate interactive HTML Krona chart
    ktImportText \
        krona_plots/${SAMPLE}.krona.txt \
        -o krona_plots/${SAMPLE}.krona.html

    echo "Krona chart: krona_plots/${SAMPLE}.krona.html"
done

# Or combine all samples into one multi-sample Krona chart
# (allows switching between samples in the browser)
ALL_INPUTS=$(for r in kraken2_results/*.report; do
    SAMPLE=$(basename ${r} .report)
    kreport2krona.py -r ${r} -o /tmp/${SAMPLE}.krona.txt
    echo "/tmp/${SAMPLE}.krona.txt,${SAMPLE}"
done | tr '\n' ' ')

# Create multi-sample Krona HTML
ktImportText ${ALL_INPUTS} -o krona_all_samples.html
echo "Multi-sample Krona chart: krona_all_samples.html"
echo "Open in any browser to explore the interactive taxonomy wheel."