Taxonomic Profiling

Read: ATCGATCGATCGATCGATCGATCG...
         ↓
K-mers: [ATCGATCGATCGATC], [TCGATCGATCGATCG], ...
         ↓
Database lookup → LCA assignment → Confidence score

Paired-end:
  Sample1_R1.fastq.gz
  Sample1_R2.fastq.gz

Single-end:
  Sample1.fastq.gz

python -m micos.cli full-run \
  --input-dir data/raw_input \
  --results-dir results \
  --threads 16 \
  --kneaddata-db /path/to/kneaddata_db \
  --kraken2-db /path/to/kraken2_db

python -m micos.cli run taxonomic-profiling \
  --input-dir data/cleaned_reads \
  --output-dir results/taxonomy \
  --threads 16 \
  --kraken2-db /path/to/kraken2_db \
  --confidence 0.1

# Step 1: Run Kraken2
kraken2 --db /path/to/kraken2_db \
  --paired sample_R1.fastq sample_R2.fastq \
  --output sample.kraken \
  --report sample.report \
  --confidence 0.1 \
  --threads 16
 
# Step 2: Run Bracken
bracken -d /path/to/kraken2_db \
  -i sample.report \
  -o sample.bracken \
  -r 150 \
  -l S
 
# Step 3: Convert to BIOM
kraken-biom *.report --fmt hdf5 -o feature-table.biom
 
# Step 4: Generate Krona
ktImportTaxonomy -o sample.krona.html sample.report

taxonomic_profiling:
  kraken2:
    # Classification confidence (0.0 - 1.0)
    # Higher = more stringent, fewer classified reads
    confidence: 0.1
    
    # Minimum base quality for k-mer matching
    min_base_quality: 20
    
    # Minimum hit groups required
    min_hit_groups: 2
    
    # Use memory-mapped database (faster, more RAM)
    memory_mapping: true
    
    # Include scientific names in output
    use_names: true

taxonomic_profiling:
  bracken:
    enabled: true
    
    # Read length (matches your data)
    read_length: 150
    
    # Taxonomic level for abundance
    # D, P, C, O, F, G, S
    level: "S"
    
    # Threshold for Bracken re-estimation
    threshold: 10

results/taxonomic_profiling/
├── raw/
│   ├── sample1.kraken      # Raw classification output
│   └── sample1.report      # Taxonomic report
├── bracken/
│   └── sample1.bracken     # Corrected abundances
├── krona/
│   └── sample1.krona.html  # Interactive visualization
└── biom/
    └── feature-table.biom  # QIIME2-compatible table

# Calculate classification percentage
# Look for: >70% classified at species level
 
grep -E "^\s+[0-9]+\s+[0-9]+\s+[0-9]+\s+[URS]" sample.report | \
  awk '{sum+=$2} END {print "Classified reads:", sum}'

# Top genera
grep "^G" sample.report | sort -k1,1nr | head -10
 
# Top species
grep "^S" sample.report | sort -k1,1nr | head -10

# Build custom Kraken2 database
kraken2-build --download-taxonomy --db custom_db
kraken2-build --add-to-library custom_genome.fa --db custom_db
kraken2-build --build --threads 16 --db custom_db

# Use higher-confidence threshold for strain analysis
taxonomic_profiling:
  kraken2:
    confidence: 0.3
    use_names: true

# Run with metadata for longitudinal analysis
python -m micos.cli run diversity-analysis \
  --input-biom feature-table.biom \
  --metadata metadata.tsv \
  --output-dir results/diversity

# 1. Lower confidence threshold
taxonomic_profiling:
  kraken2:
    confidence: 0.05
 
# 2. Use larger database
# Switch from MiniKraken to Standard
 
# 3. Check data quality
# Run FastQC on input reads

# Use memory mapping
taxonomic_profiling:
  kraken2:
    memory_mapping: true
 
# Or use SSD for database storage

# Check Kraken2 report format
head -5 sample.report
 
# Regenerate Krona
ktImportTaxonomy sample.report -o sample.krona.html