MegaSSR.sh

#!/bin/bash
#
# MegaSSR v3.0 - SSR Identification and Primer Design Pipeline
# Author: Zakaria Mahmoud
# Website: www.mahmoud.ma
# Date: March 01, 2026
#
# Clean rebuild with modular architecture
#

#=============================================================================
# Strict Error Handling
#=============================================================================
set -eo pipefail
# Temporarily disable error trap for debugging
# trap 'error_handler $? $LINENO' ERR

#=============================================================================
# Script Constants
#=============================================================================
readonly SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
readonly SCRIPT_NAME="$(basename "${BASH_SOURCE[0]}")"
readonly VERSION="3.0.0"
readonly TIMESTAMP="$(date +%Y%m%d_%H%M%S)"

# Load configuration
if [[ -f "${SCRIPT_DIR}/config/megassr.env" ]]; then
    source "${SCRIPT_DIR}/config/megassr.env"
fi

# Load utility functions
if [[ -f "${SCRIPT_DIR}/utils/functions.sh" ]]; then
    source "${SCRIPT_DIR}/utils/functions.sh"
fi

# Default paths
BIN_DIR="${SCRIPT_DIR}/bin/Script"
RESULTS_DIR="${SCRIPT_DIR}/results"
LOG_DIR="${SCRIPT_DIR}/logs"
TMP_DIR="${SCRIPT_DIR}/tmp"

#=============================================================================
# Global Variables
#=============================================================================
ANALYSIS_TYPE=1
FASTA_FILE=""
GFF_FILE=""
PROJECT_NAME="megassr_${TIMESTAMP}"
OUTPUT_DIR=""
EMAIL=""
THREADS="${NUM_THREADS:-4}"

# Additional parameters for compatibility
ALLELES="no"
MAX_ALLELE_LENGTH="1000"
PRIMER_IMAGE="50"

# SSR detection parameters (from config or defaults)
MONO_THRESHOLD="${MONO_THRESHOLD:-20}"
DI_THRESHOLD="${DI_THRESHOLD:-6}"
TRI_THRESHOLD="${TRI_THRESHOLD:-5}"
TETRA_THRESHOLD="${TETRA_THRESHOLD:-4}"
PENTA_THRESHOLD="${PENTA_THRESHOLD:-3}"
HEXA_THRESHOLD="${HEXA_THRESHOLD:-3}"
COMPOUND_DISTANCE="${COMPOUND_DISTANCE:-100}"

# Primer design parameters
PRIMER_MIN="${PRIMER_MIN_SIZE:-18}"
PRIMER_OPT="${PRIMER_OPT_SIZE:-20}"
PRIMER_MAX="${PRIMER_MAX_SIZE:-22}"
PCR_MIN="${PCR_MIN_SIZE:-250}"
PCR_MAX="${PCR_MAX_SIZE:-500}"

# Feature type for annotation
GENIC_FEATURE="gene"

# Checkpoint/Resume options
RESUME_MODE=false
RESTART_MODE=false
CHECKPOINT_FILE=""

# NCBI FTP download options
NCBI_FASTA_URL=""
NCBI_GFF_URL=""

#=============================================================================
# Usage Information
#=============================================================================
usage() {
    cat << EOF
MegaSSR v${VERSION} - SSR Identification and Primer Design Pipeline

USAGE:
    $SCRIPT_NAME --analysis <type> --fasta <file> --project <name> [options]

REQUIRED:
    --analysis, -A <1|2>    Analysis type:
                            1 = SSR detection + primer design (no GFF needed)
                                - Identifies all SSRs in the genome
                                - Designs primers for ALL SSRs
                                - Generates statistics and 3 visualization plots
                            2 = Full analysis with gene annotation (requires GFF)
                                - Identifies all SSRs in the genome
                                - Classifies SSRs as genic/intergenic using GFF
                                - Designs primers for both genic and intergenic SSRs
                                - Generates detailed statistics and 7 plots
    --fasta, -F <file>      Input FASTA file (genome sequence)
    --project, -P <name>    Project name

OPTIONAL:
    --gff, -G <file>        GFF annotation file (required for --analysis 2)
    --feature, -g <type>    Genic feature type (default: ${GENIC_FEATURE})
    --threads, -T <n>       Number of threads (default: ${THREADS})
    --email, -E <addr>      Email for notifications
    
    SSR Mining Parameters:
    --mono <n>              Mononucleotide threshold (default: ${MONO_THRESHOLD})
    --di <n>                Dinucleotide threshold (default: ${DI_THRESHOLD})
    --tri <n>               Trinucleotide threshold (default: ${TRI_THRESHOLD})
    --tetra <n>             Tetranucleotide threshold (default: ${TETRA_THRESHOLD})
    --penta <n>             Pentanucleotide threshold (default: ${PENTA_THRESHOLD})
    --hexa <n>              Hexanucleotide threshold (default: ${HEXA_THRESHOLD})
    --compound <n>          Max distance between 2 SSRs (default: ${COMPOUND_DISTANCE})
    
    Primer Design Options:
    --primer-min, -s <n>        Minimum primer size (default: ${PRIMER_MIN})
    --primer-opt, -o <n>        Optimal primer size (default: ${PRIMER_OPT})
    --primer-max, -S <n>        Maximum primer size (default: ${PRIMER_MAX})
    --product-min, -r <n>       Min product size range (default: ${PCR_MIN})
    --product-max, -R <n>       Max product size range (default: ${PCR_MAX})
    
    In Silico Validation:
    --validate, -v <yes|no>     Validate designed SSR primers (default: ${ALLELES})
    --allele-length, -l <n>     Max length for allele search (default: ${MAX_ALLELE_LENGTH})
    --primers-per-image, -i <n> Max primers per image (default: ${PRIMER_IMAGE})
    
    Checkpoint/Resume:
    --resume                    Resume from last completed phase
    --restart                   Restart from beginning (clear checkpoints)
    
    NCBI Download:
    --ncbi-fasta <url>      Download FASTA from NCBI FTP link
                            Example: https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/.../xxx_genomic.fna.gz
    --ncbi-gff <url>        Download GFF from NCBI FTP link
                            Example: https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/.../xxx_genomic.gff.gz
    
    -h, --help              Show this help message

EXAMPLES:
    # SSR identification with primer design (no GFF needed)
    $SCRIPT_NAME --analysis 1 --fasta genome.fa --project my_project
    
    # Full analysis with gene annotation and custom parameters
    $SCRIPT_NAME --analysis 2 --fasta genome.fa --gff annotation.gff --project my_project \\
        --mono 20 --di 6 --tri 5 --tetra 4 --penta 3 --hexa 3 \\
        --primer-min 18 --primer-opt 20 --primer-max 22 \\
        --product-min 250 --product-max 500 --feature gene
    
    # Resume from last checkpoint (if pipeline failed or was interrupted)
    $SCRIPT_NAME --analysis 2 --fasta genome.fa --gff annotation.gff --project my_project --resume
    
    # Download from NCBI FTP and run analysis
    $SCRIPT_NAME --analysis 2 --project my_project \\
        --ncbi-fasta https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/002/846/995/GCF_002846995.1_ASM284699v1/GCF_002846995.1_ASM284699v1_genomic.fna.gz \\
        --ncbi-gff https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/002/846/995/GCF_002846995.1_ASM284699v1/GCF_002846995.1_ASM284699v1_genomic.gff.gz
    
    # Mix local file and NCBI download
    $SCRIPT_NAME --analysis 2 --project my_project \\
        --fasta local_genome.fa \\
        --ncbi-gff https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/.../xxx_genomic.gff.gz
    
    # Restart fresh (ignore previous checkpoints)
    $SCRIPT_NAME --analysis 2 --fasta genome.fa --gff annotation.gff --project my_project --restart

CHECKPOINT SYSTEM:
    MegaSSR automatically saves progress after each phase. If the pipeline fails
    or is interrupted, use --resume to continue from the last completed phase.
    
    Phases: 2=SSR Detection, 3=Gene Annotation, 4=Genic Primers, 
            5=Intergenic Primers, 6=Post-processing, 7=Statistics, 8=Visualizations
    
    Checkpoint file: results/<project>/.megassr_checkpoint

EOF
    exit 0
}

#=============================================================================
# Parse Command Line Arguments
#=============================================================================
parse_arguments() {
    while [[ $# -gt 0 ]]; do
        case $1 in
            # Required arguments
            --analysis|-A) ANALYSIS_TYPE="$2"; shift 2 ;;
            --fasta|-F) FASTA_FILE="$2"; shift 2 ;;
            --gff|-G) GFF_FILE="$2"; shift 2 ;;
            --project|-P) PROJECT_NAME="$2"; shift 2 ;;
            
            # General options
            --feature|-g) GENIC_FEATURE="$2"; shift 2 ;;
            --threads|-T|-t) THREADS="$2"; shift 2 ;;
            --email|-e|-E) EMAIL="$2"; shift 2 ;;
            
            # SSR mining parameters
            --mono) MONO_THRESHOLD="$2"; shift 2 ;;
            --di) DI_THRESHOLD="$2"; shift 2 ;;
            --tri) TRI_THRESHOLD="$2"; shift 2 ;;
            --tetra) TETRA_THRESHOLD="$2"; shift 2 ;;
            --penta) PENTA_THRESHOLD="$2"; shift 2 ;;
            --hexa) HEXA_THRESHOLD="$2"; shift 2 ;;
            --compound|-C) COMPOUND_DISTANCE="$2"; shift 2 ;;
            
            # Primer design options
            --primer-min|-s) PRIMER_MIN="$2"; shift 2 ;;
            --primer-opt|-o|-O) PRIMER_OPT="$2"; shift 2 ;;
            --primer-max|-S) PRIMER_MAX="$2"; shift 2 ;;
            --product-min|-r) PCR_MIN="$2"; shift 2 ;;
            --product-max|-R) PCR_MAX="$2"; shift 2 ;;
            
            # In silico validation options
            --validate|-v|-B) ALLELES="$2"; shift 2 ;;
            --allele-length|-l|-L) MAX_ALLELE_LENGTH="$2"; shift 2 ;;
            --primers-per-image|-i|-I) PRIMER_IMAGE="$2"; shift 2 ;;
            
            # Checkpoint/Resume options
            --resume) RESUME_MODE=true; shift ;;
            --restart) RESTART_MODE=true; shift ;;
            
            # NCBI FTP download
            --ncbi-fasta) NCBI_FASTA_URL="$2"; shift 2 ;;
            --ncbi-gff) NCBI_GFF_URL="$2"; shift 2 ;;
            
            -h|--help) usage ;;
            *) log ERROR "Unknown option: $1"; usage ;;
        esac
    done
    
    # Validate required arguments
    # If NCBI URL is provided, local file is not required (will be downloaded)
    if [[ -z "$NCBI_FASTA_URL" ]] && [[ -z "$FASTA_FILE" ]]; then
        log ERROR "Missing required arguments: provide --fasta or --ncbi-fasta"
        usage
    fi
    
    if [[ -z "$PROJECT_NAME" ]]; then
        log ERROR "Missing required argument: --project"
        usage
    fi
    
    # For Analysis Type 2, GFF is required (from local file or NCBI)
    if [[ "$ANALYSIS_TYPE" == "2" ]] && [[ -z "$GFF_FILE" ]] && [[ -z "$NCBI_GFF_URL" ]]; then
        log ERROR "Analysis type 2 requires GFF file: use --gff or --ncbi-gff"
        exit 1
    fi
}

#=============================================================================
# NCBI FTP Download Functions
#=============================================================================

# Validate NCBI FTP URL format
validate_ncbi_url() {
    local url="$1"
    local file_type="$2"  # "FASTA" or "GFF"
    
    if [[ ! "$url" =~ ^https://ftp\.ncbi\.nlm\.nih\.gov/ ]]; then
        log ERROR "Invalid NCBI FTP URL for $file_type"
        log ERROR "URL must start with: https://ftp.ncbi.nlm.nih.gov/"
        log ERROR "Provided: $url"
        return 1
    fi
    return 0
}

# Download a file from NCBI FTP
download_ncbi_file() {
    local url="$1"
    local output_path="$2"
    local file_type="$3"  # "FASTA" or "GFF"
    
    log INFO "Downloading $file_type from NCBI FTP..."
    log INFO "URL: $url"
    
    # Validate URL
    validate_ncbi_url "$url" "$file_type" || exit 1
    
    # Download using wget or curl
    if command -v wget &> /dev/null; then
        wget -q --show-progress -O "$output_path" "$url" || {
            log ERROR "Failed to download $file_type file"
            exit 1
        }
    elif command -v curl &> /dev/null; then
        curl -# -L -o "$output_path" "$url" || {
            log ERROR "Failed to download $file_type file"
            exit 1
        }
    else
        log ERROR "Neither wget nor curl is available for downloading"
        exit 1
    fi
    
    log INFO "$file_type downloaded: $output_path"
}

# Download FASTA from NCBI if URL provided
download_ncbi_fasta() {
    local url="$1"
    local download_dir="$2"
    
    mkdir -p "$download_dir"
    
    local filename=$(basename "$url")
    local output_path="${download_dir}/${filename}"
    
    download_ncbi_file "$url" "$output_path" "FASTA"
    
    # Set global variable
    FASTA_FILE="$output_path"
}

# Download GFF from NCBI if URL provided
download_ncbi_gff() {
    local url="$1"
    local download_dir="$2"
    
    mkdir -p "$download_dir"
    
    local filename=$(basename "$url")
    local output_path="${download_dir}/${filename}"
    
    download_ncbi_file "$url" "$output_path" "GFF"
    
    # Set global variable
    GFF_FILE="$output_path"
}

#=============================================================================
# Checkpoint Functions
#=============================================================================

# Initialize checkpoint file path
init_checkpoint() {
    CHECKPOINT_FILE="${OUTPUT_DIR}/.megassr_checkpoint"
}

# Save checkpoint after completing a phase
# NOTE: Prevents duplicate entries by removing any existing entry for the same phase
save_checkpoint() {
    local phase="$1"
    local timestamp=$(date '+%Y-%m-%d %H:%M:%S')
    
    # Remove any existing entry for this phase to prevent duplicates
    if [[ -f "$CHECKPOINT_FILE" ]]; then
        local temp_file=$(mktemp)
        grep -v "^${phase}|" "$CHECKPOINT_FILE" > "$temp_file" 2>/dev/null || true
        mv "$temp_file" "$CHECKPOINT_FILE"
    fi
    
    echo "${phase}|${timestamp}" >> "$CHECKPOINT_FILE"
    log DEBUG "Checkpoint saved: Phase $phase completed at $timestamp"
}

# Get the last completed phase number
get_last_checkpoint() {
    if [[ -f "$CHECKPOINT_FILE" ]]; then
        local last_line=$(tail -1 "$CHECKPOINT_FILE" 2>/dev/null)
        local phase=$(echo "$last_line" | cut -d'|' -f1)
        echo "$phase"
    else
        echo "0"
    fi
}

# Check if a phase should be skipped (already completed)
should_skip_phase() {
    local phase="$1"
    local last_completed=$(get_last_checkpoint)
    
    if [[ "$RESUME_MODE" == "true" ]] && [[ "$last_completed" -ge "$phase" ]]; then
        return 0  # true - skip this phase
    fi
    return 1  # false - run this phase
}

# Clear all checkpoints
clear_checkpoints() {
    if [[ -f "$CHECKPOINT_FILE" ]]; then
        rm -f "$CHECKPOINT_FILE"
        log INFO "Checkpoints cleared - starting fresh"
    fi
}

# Display checkpoint status
show_checkpoint_status() {
    if [[ -f "$CHECKPOINT_FILE" ]]; then
        log INFO "Previous run detected. Checkpoint history:"
        while IFS='|' read -r phase timestamp; do
            log INFO "  Phase $phase completed at $timestamp"
        done < "$CHECKPOINT_FILE"
        local last=$(get_last_checkpoint)
        log INFO "Resuming from Phase $((last + 1))"
    fi
}

#=============================================================================
# Initialize Pipeline
#=============================================================================
initialize_pipeline() {
    log INFO "=========================================="
    log INFO "MegaSSR v${VERSION}"
    log INFO "=========================================="
    log INFO "Project: $PROJECT_NAME"
    log INFO "Analysis Type: $ANALYSIS_TYPE"
    
    # Create output directory first (needed for NCBI download)
    OUTPUT_DIR="${RESULTS_DIR}/${PROJECT_NAME}"
    mkdir -p "$OUTPUT_DIR"
    mkdir -p "${OUTPUT_DIR}/fasta"
    mkdir -p "${OUTPUT_DIR}/intermediate"
    mkdir -p "${LOG_DIR}"
    mkdir -p "${TMP_DIR}/${PROJECT_NAME}"
    
    # Download from NCBI if URLs provided
    local download_dir="${OUTPUT_DIR}/ncbi_download"
    
    if [[ -n "$NCBI_FASTA_URL" ]]; then
        log INFO "=========================================="
        log INFO "Downloading from NCBI FTP"
        log INFO "=========================================="
        download_ncbi_fasta "$NCBI_FASTA_URL" "$download_dir"
    fi
    
    if [[ -n "$NCBI_GFF_URL" ]]; then
        if [[ -z "$NCBI_FASTA_URL" ]]; then
            log INFO "=========================================="
            log INFO "Downloading from NCBI FTP"
            log INFO "=========================================="
        fi
        download_ncbi_gff "$NCBI_GFF_URL" "$download_dir"
    fi
    
    if [[ -n "$NCBI_FASTA_URL" ]] || [[ -n "$NCBI_GFF_URL" ]]; then
        log INFO "NCBI download complete!"
        log INFO "=========================================="
    fi
    
    log INFO "FASTA: $FASTA_FILE"
    [[ -n "$GFF_FILE" ]] && log INFO "GFF: $GFF_FILE"
    log INFO "Threads: $THREADS"
    log INFO "=========================================="
    
    # Initialize checkpoint system
    init_checkpoint
    
    # Handle restart mode - clear existing checkpoints
    if [[ "$RESTART_MODE" == "true" ]]; then
        clear_checkpoints
    fi
    
    # Handle resume mode - show status
    if [[ "$RESUME_MODE" == "true" ]]; then
        show_checkpoint_status
    fi
    
    # Validate input files
    validate_fasta "$FASTA_FILE" || exit 1
    if [[ -n "$GFF_FILE" ]]; then
        validate_gff "$GFF_FILE" || exit 1
        # Check chromosome ID consistency between FASTA and GFF
        validate_chromosome_consistency "$FASTA_FILE" "$GFF_FILE"
    fi
    
    # Activate conda environment (non-critical)
    log DEBUG "About to call activate_conda_env"
    set +e  # Disable exit on error
    activate_conda_env "$CONDA_ENV"
    local conda_result=$?
    set -e  # Re-enable exit on error
    log DEBUG "activate_conda_env returned: $conda_result"
    
    log INFO "Initialization complete"
}

#=============================================================================
# Phase 2: SSR Identification using MISA
#=============================================================================
run_ssr_identification() {
    log INFO "=========================================="
    log INFO "PHASE 2: SSR Identification (MISA)"
    log INFO "=========================================="
    
    local fasta_dir="${OUTPUT_DIR}/fasta"
    local genome_fasta="${fasta_dir}/${PROJECT_NAME}_genomic.fa"
    
    # Copy FASTA to working directory (decompress if gzipped)
    log DEBUG "Preparing genome FASTA file"
    if [[ "$FASTA_FILE" == *.gz ]]; then
        log DEBUG "Decompressing gzipped FASTA file"
        zcat "$FASTA_FILE" > "${genome_fasta}.tmp"
    else
        cp "$FASTA_FILE" "${genome_fasta}.tmp"
    fi
    
    # Normalize FASTA headers: replace spaces with underscores for consistency
    # MISA converts spaces to underscores in chromosome IDs, so we do the same in the FASTA
    # This ensures extraction scripts can match chromosome IDs correctly
    log DEBUG "Normalizing FASTA headers (spaces → underscores)"
    sed 's/^>\([^ ]*\) />\1_/; :a; s/^\(>[^ ]*\) /\1_/; ta' "${genome_fasta}.tmp" > "$genome_fasta"
    rm -f "${genome_fasta}.tmp"
    
    # Create MISA configuration file
    local misa_ini="${fasta_dir}/misa.ini"
    cat > "$misa_ini" << EOF
definition(unit_size,min_repeats): 1-${MONO_THRESHOLD} 2-${DI_THRESHOLD} 3-${TRI_THRESHOLD} 4-${TETRA_THRESHOLD} 5-${PENTA_THRESHOLD} 6-${HEXA_THRESHOLD}
interruptions(max_difference_between_2_SSRs): ${COMPOUND_DISTANCE}
EOF
    
    log DEBUG "MISA configuration:"
    log DEBUG "  Mono: ${MONO_THRESHOLD}, Di: ${DI_THRESHOLD}, Tri: ${TRI_THRESHOLD}"
    log DEBUG "  Tetra: ${TETRA_THRESHOLD}, Penta: ${PENTA_THRESHOLD}, Hexa: ${HEXA_THRESHOLD}"
    log DEBUG "  Compound distance: ${COMPOUND_DISTANCE}"
    
    # Run MISA
    log INFO "Running MISA for SSR detection"
    cd "$fasta_dir" || exit 1
    
    perl "${BIN_DIR}/misa.pl" "$genome_fasta" "$fasta_dir" 2>&1 | tee -a "${LOG_DIR}/${PROJECT_NAME}.log" || {
        log ERROR "MISA execution failed"
        exit 1
    }
    
    cd - > /dev/null
    
    # Check MISA output
    local misa_output="${genome_fasta}.misa"
    if [[ ! -f "$misa_output" ]] || [[ ! -s "$misa_output" ]]; then
        log ERROR "MISA did not produce output file"
        exit 1
    fi
    
    # Create formatted SSR list (add organism ID and convert type codes)
    local ssr_file="${OUTPUT_DIR}/intermediate/${PROJECT_NAME}_SSR.txt"
    awk -F'\t' -v org="$PROJECT_NAME" '{print "\t" org "\t" $1 "\t" $2 "\t" $3 "\t" $4 "\t" $5 "\t" $6 "\t" $7}' "$misa_output" | \
    sed 's/\tp1\t/\tMono\t/g' | \
    sed 's/\tp2\t/\tDi\t/g' | \
    sed 's/\tp3\t/\tTri\t/g' | \
    sed 's/\tp4\t/\tTetra\t/g' | \
    sed 's/\tp5\t/\tPenta\t/g' | \
    sed 's/\tp6\t/\tHexa\t/g' | \
    sed 's/\tc\t/\tCompound\t/g' | \
    sed 's/\tc\*\t/\tCompound\t/g' > "$ssr_file"
    
    local ssr_count=$(wc -l < "$ssr_file")
    log INFO "SSR identification complete: $ssr_count SSRs found"
    
    # Create SSR statistics summary
    local stats_file="${genome_fasta}.statistics"
    if [[ -f "$stats_file" ]]; then
        log DEBUG "MISA statistics generated"
        cp "$stats_file" "${OUTPUT_DIR}/${PROJECT_NAME}-MegaSSR_Results/SSR_statistics.txt" 2>/dev/null || true
    fi
    
    # For Analysis Type 1 (no GFF): treat ALL SSRs as intergenic
    # This matches the behavior of the old MegaSSR where -A 1 designed primers for all SSRs
    if [[ "$ANALYSIS_TYPE" == "1" ]]; then
        log INFO "Analysis Type 1: Treating all $ssr_count SSRs as intergenic (no gene annotation)"
        local intergenic_file="${OUTPUT_DIR}/intermediate/${PROJECT_NAME}_Accept_Intergenic_SSR.txt"
        cp "$ssr_file" "$intergenic_file"
        log DEBUG "Created intergenic SSR file from all SSRs"
    fi
    
    log INFO "SSR identification phase complete"
}

#=============================================================================
# Phase 3: Gene Annotation (Genic vs Intergenic Classification)
#=============================================================================
run_gene_annotation() {
    log INFO "=========================================="
    log INFO "PHASE 3: Gene Annotation"
    log INFO "=========================================="
    
    local intermediate_dir="${OUTPUT_DIR}/intermediate"
    local ssr_file="${intermediate_dir}/${PROJECT_NAME}_SSR.txt"
    local gff_features="${intermediate_dir}/${PROJECT_NAME}_gff_features.txt"
    
    # Extract specified features from GFF
    log INFO "Extracting '$GENIC_FEATURE' features from GFF"
    python3 "${BIN_DIR}/extract_gff_features.py" \
        "$GFF_FILE" \
        "$GENIC_FEATURE" \
        "$gff_features" || {
        log ERROR "Failed to extract GFF features"
        exit 1
    }
    
    local feature_count=$(wc -l < "$gff_features")
    log DEBUG "Extracted $feature_count $GENIC_FEATURE features"
    
    # Run genic SSR classification
    log INFO "Classifying SSRs as genic or intergenic"
    local genic_with_feature="${intermediate_dir}/${PROJECT_NAME}.Genic_SSR_with_feature.txt"
    local genic_for_db="${OUTPUT_DIR}/${PROJECT_NAME}-MegaSSR_Results/${PROJECT_NAME}.Genic_SSR_with_feature.txt"
    local genic_no_feature="${intermediate_dir}/${PROJECT_NAME}.Genic_SSR_none_feature.txt"
    
    mkdir -p "${OUTPUT_DIR}/${PROJECT_NAME}-MegaSSR_Results"
    
    perl "${BIN_DIR}/get-genic-SSR.pl" \
        "$ssr_file" \
        "$gff_features" \
        "$genic_with_feature" \
        "$genic_for_db" \
        "$genic_no_feature" 2>&1 | tee -a "${LOG_DIR}/${PROJECT_NAME}.log" || {
        log WARN "Genic annotation completed with warnings"
    }
    
    # Count results
    local genic_count=0
    [[ -f "$genic_with_feature" ]] && genic_count=$(wc -l < "$genic_with_feature")
    log INFO "Genic SSRs identified: $genic_count"
    
    # Create intergenic SSR file using grep to exclude genic SSRs
    # The genic_no_feature file contains SSRs that ARE in genes (despite the name)
    # It may contain duplicates (SSRs overlapping multiple genes), so deduplicate first
    # Use grep -F -x -v -f to find SSRs in original list NOT in genic list
    local intergenic_file="${intermediate_dir}/${PROJECT_NAME}_Accept_Intergenic_SSR.txt"
    
    if [[ -f "$genic_no_feature" ]] && [[ -s "$genic_no_feature" ]]; then
        log DEBUG "Finding intergenic SSRs using grep exclusion (after deduplication)"
        # Deduplicate genic file first to avoid grep issues with duplicate patterns
        sort -u "$genic_no_feature" > "${genic_no_feature}.uniq"
        grep -F -x -v -f "${genic_no_feature}.uniq" "$ssr_file" > "$intergenic_file" || {
            log WARN "No intergenic SSRs found (all SSRs are genic)"
            touch "$intergenic_file"
        }
        rm -f "${genic_no_feature}.uniq"
    else
        # If no genic SSRs found, all SSRs are intergenic
        log DEBUG "No genic SSRs found, all SSRs are intergenic"
        cp "$ssr_file" "$intergenic_file"
    fi
    
    local intergenic_count=0
    [[ -f "$intergenic_file" ]] && intergenic_count=$(wc -l < "$intergenic_file")
    log INFO "Intergenic SSRs identified: $intergenic_count"
    
    # Validation
    local total_ssr=$(wc -l < "$ssr_file")
    log INFO "Total SSRs: $total_ssr (Genic: $genic_count, Intergenic: $intergenic_count)"
    
    log INFO "Gene annotation phase complete"
}

#=============================================================================
# Phase 4: Genic Primer Design
#=============================================================================
run_genic_primer_design() {
    log INFO "=========================================="
    log INFO "PHASE 4: Genic Primer Design"
    log INFO "=========================================="
    
    local intermediate_dir="${OUTPUT_DIR}/intermediate"
    local results_dir="${OUTPUT_DIR}/${PROJECT_NAME}-MegaSSR_Results"
    local genic_with_feature="${results_dir}/${PROJECT_NAME}.Genic_SSR_with_feature.txt"
    
    # Check if genic SSR file exists
    if [[ ! -f "$genic_with_feature" ]] || [[ ! -s "$genic_with_feature" ]]; then
        log WARN "No genic SSRs found, skipping genic primer design"
        return 0
    fi
    
    local genic_count=$(wc -l < "$genic_with_feature")
    # Subtract 1 for header if present
    if head -1 "$genic_with_feature" | grep -q "Process Id"; then
        genic_count=$((genic_count - 1))
    fi
    log INFO "Designing primers for $genic_count genic SSRs"
    
    # Create formatted genic SSR file for primer design (like reference table11_Genic_misa_feature.txt)
    # Input format (16 fields): empty, project_id, chr_id, repeat_num, repeat_type, repeat_seq, repeat_len, start, end, gene, gene_start, gene_end, ., strand, ., annotation
    # Output format (14 fields): empty, project_id, chr_id, repeat_num, repeat_type, repeat_seq, repeat_len, start, end, gene, gene_start, gene_end, strand, annotation
    local genic_formatted="${intermediate_dir}/${PROJECT_NAME}_Genic_SSR_for_primer.txt"
    # Note: Input file has NO header at this point (header is added in Phase 6 post-processing)
    # Reformat 16-field input to 14-field output for primer design
    awk -F '\t' '{ print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10"\t"$11"\t"$12"\t"$14"\t"$16 }' "$genic_with_feature" > "$genic_formatted"
    log DEBUG "Created formatted genic SSR file for primer design ($(wc -l < "$genic_formatted") entries)"
    
    # Create working directories (clean first to avoid append issues)
    local gdesignprimer_dir="${OUTPUT_DIR}/gdesignprimer"
    local gdesignprimer_results="${OUTPUT_DIR}/gdesignprimerresults"
    rm -rf "$gdesignprimer_dir" "$gdesignprimer_results"
    mkdir -p "$gdesignprimer_dir" "$gdesignprimer_results"
    
    # Split formatted genic SSR file into chunks of 500 for parallel processing
    log DEBUG "Splitting genic SSRs into chunks of 500"
    cd "$gdesignprimer_dir" || exit 1
    split -d -l 500 "$genic_formatted"
    
    local split_count=$(ls -1 x* 2>/dev/null | wc -l)
    log DEBUG "Created $split_count split files for processing"
    
    if [[ $split_count -eq 0 ]]; then
        log ERROR "Failed to create split files"
        exit 1
    fi
    
    # Process each split file with gdesignprimer_threads.py
    log INFO "Extracting flanking sequences for genic SSRs"
    local genome_fasta="${OUTPUT_DIR}/fasta/${PROJECT_NAME}_genomic.fa"
    
    # gdesignprimer_threads.py expects: genome_fasta split_dir results_dir threads extractseq_script p3_in_script
    python3 "${BIN_DIR}/gdesignprimer_threads.py" \
        "$genome_fasta" \
        "$gdesignprimer_dir" \
        "$gdesignprimer_results" \
        "$THREADS" \
        "${BIN_DIR}/extractseq-id-start-end-genic.pl" \
        "${BIN_DIR}/modified_p3_in_genic.pl" || {
        log WARN "Warning during sequence extraction"
    }
    
    # Concatenate all extracted sequences
    local flanking_fa="${results_dir}/Genic_SSR_flanking_regions.fa"
    if ls "$gdesignprimer_results"/*.fa 1>/dev/null 2>&1; then
        cat "$gdesignprimer_results"/*.fa > "$flanking_fa"
        local flanking_count=$(grep -c "^>" "$flanking_fa" || echo "0")
        log INFO "Extracted flanking sequences: $flanking_count"
    else
        log WARN "No flanking sequences extracted"
        touch "$flanking_fa"
    fi
    
    # Run actual primer design using gdesignprimer2_threads.py
    log INFO "Running Primer3 for genic SSR primer design"
    python3 "${BIN_DIR}/gdesignprimer2_threads.py" \
        "$gdesignprimer_dir" \
        "$gdesignprimer_results" \
        "$THREADS" \
        "${BIN_DIR}/extractseq-id-start-end-genic.pl" \
        "${BIN_DIR}/modified_p3_in_genic.pl" \
        "$PROJECT_NAME" \
        "${BIN_DIR}/modified_p3_out-genic.pl" \
        "${BIN_DIR}/print-primers-line-genicCCC.pl" \
        "$intermediate_dir" || {
        log ERROR "Primer design failed"
        exit 1
    }
    
    # Deduplicate primers (by strand and GC%)
    local primers_raw="${intermediate_dir}/${PROJECT_NAME}.Extract-Genic-seq-out.fasta.results2.txt"
    local primers_dedup1="${intermediate_dir}/${PROJECT_NAME}.Extract-Genic-seq-out.fasta.results2_2.txt"
    local primers_final="${results_dir}/${PROJECT_NAME}.Genic-primers.txt"
    
    if [[ -f "$primers_raw" ]] && [[ -s "$primers_raw" ]]; then
        log DEBUG "Deduplicating primers by strand (field 15)"
        awk '!seen[$15]++' "$primers_raw" > "$primers_dedup1"
        
        log DEBUG "Deduplicating primers by GC% (field 19)"
        awk '!seen[$19]++' "$primers_dedup1" > "$primers_final"
        
        local final_count=$(wc -l < "$primers_final")
        log INFO "Genic primers designed: $final_count"
        
        # Create statistics file
        local stats_file="${intermediate_dir}/${PROJECT_NAME}.Extract-Genic-out.fasta.stat"
        echo -e "$final_count\t${PROJECT_NAME}.Genic-primers.txt" > "$stats_file"
    else
        log WARN "No primers generated"
        touch "$primers_final"
    fi
    
    log INFO "Genic primer design phase complete"
}

#=============================================================================
# Phase 5: Intergenic Primer Design
#=============================================================================
run_intergenic_primer_design() {
    log INFO "=========================================="
    log INFO "PHASE 5: Intergenic Primer Design"
    log INFO "=========================================="
    
    local intermediate_dir="${OUTPUT_DIR}/intermediate"
    local results_dir="${OUTPUT_DIR}/${PROJECT_NAME}-MegaSSR_Results"
    local intergenic_file="${intermediate_dir}/${PROJECT_NAME}_Accept_Intergenic_SSR.txt"
    
    # Ensure results directory exists (important for Analysis Type 1 which skips Phase 3-4)
    mkdir -p "$results_dir"
    
    # Check if intergenic SSR file exists
    if [[ ! -f "$intergenic_file" ]] || [[ ! -s "$intergenic_file" ]]; then
        log WARN "No intergenic SSRs found, skipping intergenic primer design"
        return 0
    fi
    
    local intergenic_count=$(wc -l < "$intergenic_file")
    log INFO "Designing primers for $intergenic_count intergenic SSRs"
    
    # Create working directories (clean first to avoid append issues)
    local designprimer_dir="${OUTPUT_DIR}/designprimer"
    local designprimer_results="${OUTPUT_DIR}/designprimerresults"
    rm -rf "$designprimer_dir" "$designprimer_results"
    mkdir -p "$designprimer_dir" "$designprimer_results"
    
    # Split intergenic SSR file into chunks of 500 for parallel processing
    log DEBUG "Splitting intergenic SSRs into chunks of 500"
    cd "$designprimer_dir" || exit 1
    split -d -l 500 "$intergenic_file"
    
    local split_count=$(ls -1 x* 2>/dev/null | wc -l)
    log DEBUG "Created $split_count split files for processing"
    
    if [[ $split_count -eq 0 ]]; then
        log ERROR "Failed to create split files"
        exit 1
    fi
    
    # Process split files with xdesignprimer_threads.py for sequence extraction
    log INFO "Extracting flanking sequences for intergenic SSRs"
    local genome_fasta="${OUTPUT_DIR}/fasta/${PROJECT_NAME}_genomic.fa"
    
    # xdesignprimer_threads.py expects: genome_fasta split_dir results_dir threads extractseq_script p3_in_script
    python3 "${BIN_DIR}/xdesignprimer_threads.py" \
        "$genome_fasta" \
        "$designprimer_dir" \
        "$designprimer_results" \
        "$THREADS" \
        "${BIN_DIR}/extractseq-id-start-end-intergenic.pl" \
        "${BIN_DIR}/modified_p3_in.pl" || {
        log WARN "Warning during sequence extraction"
    }
    
    # Concatenate all extracted sequences
    local flanking_fa="${results_dir}/Intergenic_SSR_flanking_regions.fa"
    if ls "$designprimer_results"/*.fa 1>/dev/null 2>&1; then
        cat "$designprimer_results"/*.fa > "$flanking_fa"
        local flanking_count=$(grep -c "^>" "$flanking_fa" || echo "0")
        log INFO "Extracted flanking sequences: $flanking_count"
    else
        log WARN "No flanking sequences extracted"
        touch "$flanking_fa"
    fi
    
    # Run actual primer design using designprimer3_threads.py
    log INFO "Running Primer3 for intergenic SSR primer design"
    python3 "${BIN_DIR}/designprimer3_threads.py" \
        "$designprimer_dir" \
        "$designprimer_results" \
        "$THREADS" \
        "${BIN_DIR}/extractseq-id-start-end-intergenic.pl" \
        "${BIN_DIR}/modified_p3_in.pl" \
        "$PROJECT_NAME" \
        "${BIN_DIR}/modified_p3_out-intergenic.pl" \
        "${BIN_DIR}/print-primers-line-nongenicCCC.pl" \
        "$intermediate_dir" || {
        log ERROR "Primer design failed"
        exit 1
    }
    
    # Deduplicate primers (by SSR position field 9 and GC% field 13)
    local primers_raw="${intermediate_dir}/${PROJECT_NAME}.interGenic-primers3.txt"
    local primers_dedup1="${intermediate_dir}/${PROJECT_NAME}.interGenic-primers3_2.txt"
    local primers_final="${results_dir}/${PROJECT_NAME}.interGenic-primers.txt"
    
    if [[ -f "$primers_raw" ]] && [[ -s "$primers_raw" ]]; then
        log DEBUG "Deduplicating primers by SSR position (field 9)"
        awk '!seen[$9]++' "$primers_raw" > "$primers_dedup1"
        
        log DEBUG "Deduplicating primers by GC% (field 13)"
        awk '!seen[$13]++' "$primers_dedup1" > "$primers_final"
        
        local final_count=$(wc -l < "$primers_final")
        log INFO "Intergenic primers designed: $final_count"
        
        # Create statistics file
        local stats_file="${intermediate_dir}/${PROJECT_NAME}.Extract-intergenic-out.fasta.stat"
        echo -e "$final_count\t${PROJECT_NAME}.interGenic-primers.txt" > "$stats_file"
    else
        log WARN "No primers generated"
        touch "$primers_final"
    fi
    
    log INFO "Intergenic primer design phase complete"
}

#=============================================================================
# Phase 6: Post-processing and Output Formatting
#=============================================================================
run_post_processing() {
    log INFO "=========================================="
    log INFO "PHASE 6: Post-processing"
    log INFO "=========================================="
    
    local intermediate_dir="${OUTPUT_DIR}/intermediate"
    local results_dir="${OUTPUT_DIR}/${PROJECT_NAME}-MegaSSR_Results"
    
    # Generate SSR statistics from MISA output
    generate_ssr_statistics
    
    # Create non-redundant SSR library if usearch available
    create_nonredundant_library
    
    # Export final CSVs with friendly names
    export_csv_files
    
    log INFO "Post-processing complete"
}

generate_ssr_statistics() {
    log INFO "Generating SSR statistics"
    
    local misa_stats="${OUTPUT_DIR}/fasta/${PROJECT_NAME}_genomic.fa.statistics"
    local results_dir="${OUTPUT_DIR}/${PROJECT_NAME}-MegaSSR_Results"
    
    if [[ ! -f "$misa_stats" ]]; then
        log WARN "MISA statistics file not found"
        return 1
    fi
    
    # Copy MISA statistics to results
    local stats_file="${results_dir}/SSR_statistics.txt"
    cp "$misa_stats" "$stats_file"
    
    log INFO "SSR statistics generated"
}

create_nonredundant_library() {
    log INFO "Creating non-redundant SSR library"
    
    local results_dir="${OUTPUT_DIR}/${PROJECT_NAME}-MegaSSR_Results"
    local genic_fa="${results_dir}/Genic_SSR_flanking_regions.fa"
    local intergenic_fa="${results_dir}/Intergenic_SSR_flanking_regions.fa"
    local combined_fa="${OUTPUT_DIR}/intermediate/${PROJECT_NAME}_combined_SSR_flanking.fa"
    local nonred_fa="${results_dir}/Non-redundant_SSR_library.fasta"
    local nonred_log="${results_dir}/SSR.non-redundant.log"
    
    # Check if both flanking region files exist
    if [[ ! -f "$genic_fa" ]] && [[ ! -f "$intergenic_fa" ]]; then
        log WARN "No flanking region files found, skipping non-redundant library"
        return 0
    fi
    
    # Combine flanking regions (both files may not exist)
    > "$combined_fa"  # Create empty file first
    [[ -f "$genic_fa" ]] && cat "$genic_fa" >> "$combined_fa"
    [[ -f "$intergenic_fa" ]] && cat "$intergenic_fa" >> "$combined_fa"
    
    local seq_count=$(grep -c "^>" "$combined_fa" || echo "0")
    log INFO "Combined $seq_count flanking sequences"
    
    # Define usearch binary path
    # Priority: 1) USEARCH_BIN env var, 2) bundled binary, 3) system usearch
    # NOTE: Bundled binary is i86linux32; set USEARCH_BIN for x86_64 systems if needed
    local usearch_bin="${USEARCH_BIN:-${BIN_DIR}/usearch11.0.667_i86linux32}"
    
    # Try to use usearch for clustering
    local usearch_cmd=""
    if [[ -n "${USEARCH_BIN:-}" ]] && [[ -x "$USEARCH_BIN" ]]; then
        usearch_cmd="$USEARCH_BIN"
        log INFO "Using custom usearch from USEARCH_BIN: $USEARCH_BIN"
    elif [[ -x "$usearch_bin" ]]; then
        # Test if bundled binary is architecture-compatible
        if "$usearch_bin" --version >/dev/null 2>&1; then
            usearch_cmd="$usearch_bin"
            log INFO "Using bundled usearch for clustering"
        else
            log WARN "Bundled usearch binary not compatible with this architecture"
        fi
    fi
    
    # Fallback to system usearch if bundled not available/compatible
    if [[ -z "$usearch_cmd" ]] && command -v usearch >/dev/null 2>&1; then
        usearch_cmd="usearch"
        log INFO "Using system usearch for clustering"
    fi
    
    if [[ -n "$usearch_cmd" ]]; then
        local temp_dir="${OUTPUT_DIR}/usearch_temp"
        mkdir -p "$temp_dir"
        
        # Step 1: Sort by length (as in original pipeline)
        local sorted_fa="${temp_dir}/SSR_with_flanking_regions_sorted.fa"
        "$usearch_cmd" -sortbylength "$combined_fa" \
            -fastaout "$sorted_fa" \
            -log "${temp_dir}/usearch_sort.log" 2>&1
        
        # Step 2: Cluster at 90% identity
        "$usearch_cmd" -cluster_fast "$sorted_fa" \
            -id 0.90 \
            -centroids "${temp_dir}/centroids.fa" \
            -uc "${temp_dir}/result.uc" \
            -consout "$nonred_fa" \
            -msaout "${temp_dir}/aligned.fasta" \
            -log "${temp_dir}/usearch_cluster.log" 2>&1
        
        if [[ $? -eq 0 ]] && [[ -f "$nonred_fa" ]]; then
            # Copy log and clean up sensitive lines
            sed '2d;3d;11d' "${temp_dir}/usearch_cluster.log" > "$nonred_log" 2>/dev/null || \
                cp "${temp_dir}/usearch_cluster.log" "$nonred_log"
            
            local cluster_count=$(grep -c "^>" "$nonred_fa" || echo "0")
            log INFO "Created non-redundant library: $cluster_count clusters (from $seq_count sequences)"
            rm -rf "$temp_dir"
            return 0
        else
            log WARN "usearch clustering failed, falling back to simple deduplication"
            rm -rf "$temp_dir"
        fi
    fi
    
    # Fallback: simple deduplication by sequence
    log INFO "Using simple deduplication (usearch not available or failed)"
    
    # Simple deduplication using awk - just copy as-is (no real dedup for now)
    # Real clustering requires usearch
    cp "$combined_fa" "$nonred_fa"
    echo "Copied $seq_count sequences (simple deduplication - install usearch for clustering)" > "$nonred_log" 2>&1
    
    local dedup_count=$(grep -c "^>" "$nonred_fa" || echo "0")
    log INFO "Created non-redundant library: $dedup_count unique sequences (from $seq_count total)"
}

export_csv_files() {
    log INFO "Exporting CSV files"
    
    local results_dir="${OUTPUT_DIR}/${PROJECT_NAME}-MegaSSR_Results"
    local intermediate_dir="${OUTPUT_DIR}/intermediate"
    
    # POSIX-safe sed -i replacement function (works on both Linux and macOS)
    # Usage: sed_inplace 'sed_expression' 'file'
    sed_inplace() {
        local expr="$1"
        local file="$2"
        local temp_file=$(mktemp)
        sed "$expr" "$file" > "$temp_file" && mv "$temp_file" "$file"
    }
    
    # Add headers to genic primers file (like reference)
    local genic_primers="${results_dir}/${PROJECT_NAME}.Genic-primers.txt"
    if [[ -f "$genic_primers" ]] && [[ -s "$genic_primers" ]]; then
        # Add header matching reference format (POSIX-safe)
        sed_inplace '1s/^/\tSequence ID\tProcess Id\tRepeat number\tRepeat Type\tRepeat Start\tRepeat End\tPrimer Start\tPrimer End\tRepeat Type\tRepeat Length\tGene\tGene Start\tGene End\tSequence Strand (Gene)\tForward Primer\tForward Tm\tForward Size (bp)\tForward GC\tReverse Primer\tReverse Tm\tReverse Size (bp)\tReverse GC\tProduct Size (bp)\tAnnotation\n/' "$genic_primers"
        
        # Also create renamed version like reference
        cp "$genic_primers" "${results_dir}/Desinged genic SSR primer.csv"
        log INFO "Exported genic primers CSV with header"
    fi
    
    # Add headers to intergenic primers file (like reference)
    local intergenic_primers="${results_dir}/${PROJECT_NAME}.interGenic-primers.txt"
    if [[ -f "$intergenic_primers" ]] && [[ -s "$intergenic_primers" ]]; then
        # Add header matching reference format (POSIX-safe)
        sed_inplace '1s/^/\tSequence ID\tProcess Id\tRepeat number\tRepeat Type\tRepeat Sequence\tRepeat Length\tRepeat Start\tRepeat End\tPrimer Start\tPrimer End\tForward Primer\tForward Tm\tForward Size (bp)\tForward GC\tReverse Primer\tReverse Tm\tReverse Size (bp)\tReverse GC\tProduct Size (bp)\n/' "$intergenic_primers"
        
        # Also create renamed version like reference
        cp "$intergenic_primers" "${results_dir}/Desinged non-genic SSR primers.csv"
        log INFO "Exported intergenic primers CSV with header"
    fi
    
    # Add headers to genic SSR annotations (like reference)
    local genic_ssr="${results_dir}/${PROJECT_NAME}.Genic_SSR_with_feature.txt"
    if [[ -f "$genic_ssr" ]] && [[ -s "$genic_ssr" ]]; then
        # Add header matching reference format (POSIX-safe)
        sed_inplace '1s/^/\tProcess Id\tSequence ID\tRepeat number\tRepeat Type\tRepeat Sequence\tRepeat Length\tRepeat Start\tRepeat End\tGene\tGene Start\tGene End\t\tSequence Strand (Gene)\t\tAnnotation\n/' "$genic_ssr"
        
        # Also create renamed version like reference
        cp "$genic_ssr" "${results_dir}/Genic SSR repeats with annotations.csv"
        log INFO "Exported genic SSR annotations CSV with header"
    fi
    
    # Add headers to intergenic SSR list (like reference)
    local intergenic_ssr="${intermediate_dir}/${PROJECT_NAME}_Accept_Intergenic_SSR.txt"
    if [[ -f "$intergenic_ssr" ]] && [[ -s "$intergenic_ssr" ]]; then
        # Create copy with header
        local intergenic_ssr_csv="${results_dir}/non-genic SSR repeats with annotations.csv"
        sed '1s/^/\tProcess Id\tSequence Id\tRepeat number\tRepeat Type\tRepeat Sequence\tRepeat Length\tRepeat Start\tRepeat End\n/' "$intergenic_ssr" > "$intergenic_ssr_csv"
        log INFO "Exported intergenic SSR list CSV with header"
    fi
    
    # Copy MISA output as CSV with header
    local misa_file="${OUTPUT_DIR}/fasta/${PROJECT_NAME}_genomic.fa.misa"
    if [[ -f "$misa_file" ]] && [[ -s "$misa_file" ]]; then
        local misa_csv="${results_dir}/Identified SSR motifs table.csv"
        sed '1s/^/\tProcess Id\tSequence Id\tRepeat number\tRepeat Type\tRepeat Sequence\tRepeat Length\tRepeat Start\tRepeat End\n/' "$misa_file" > "$misa_csv"
        log INFO "Exported SSR motifs table CSV with header"
    fi
    
    # Rename flanking regions to match reference
    if [[ -f "${results_dir}/Genic_SSR_flanking_regions.fa" ]]; then
        cp "${results_dir}/Genic_SSR_flanking_regions.fa" "${results_dir}/Genic SSR flanking regions.fa"
    fi
    if [[ -f "${results_dir}/Intergenic_SSR_flanking_regions.fa" ]]; then
        cp "${results_dir}/Intergenic_SSR_flanking_regions.fa" "${results_dir}/Non-genic SSR flanking regions.fa"
    fi
    
    # Rename non-redundant library to match reference
    if [[ -f "${results_dir}/Non-redundant_SSR_library.fasta" ]]; then
        cp "${results_dir}/Non-redundant_SSR_library.fasta" "${results_dir}/Non-redundant SSR library.fasta"
    fi
}

#=============================================================================
# Phase 7: Advanced Statistics
#=============================================================================
generate_advanced_statistics() {
    log INFO "=========================================="
    log INFO "PHASE 7: Advanced Statistics"
    log INFO "=========================================="
    
    local results_dir="${OUTPUT_DIR}/${PROJECT_NAME}-MegaSSR_Results"
    local intermediate_dir="${OUTPUT_DIR}/intermediate"
    local fasta_dir="${OUTPUT_DIR}/fasta"
    local stats_file="${fasta_dir}/${PROJECT_NAME}_genomic.fa.statistics"
    
    # Create temporary directory for intermediate files
    local temp_dir="${OUTPUT_DIR}/temp_stats"
    mkdir -p "$temp_dir"
    
    # Extract frequency tables from MISA statistics
    log INFO "Extracting frequency tables from MISA statistics"
    
    # 1. Distribution of different repeat type classes - Create proper CSV format
    # Note: MISA statistics counts compound SSRs by their constituent unit sizes (expands them)
    # This matches the reference hpc-results which shows 5811 total instead of 5446
    if [[ -f "$stats_file" ]]; then
        # Extract the distribution section and create CSV with header
        local dist_csv="${results_dir}/The distribution of the different SSR classes.csv"
        echo -e "\tProcess Id\tRepeat Type\tTotal No. of present" > "$dist_csv"
        
        # Parse MISA statistics for distribution - look for lines after "Unit size" header
        # Format: "1       333" (number tab number)
        awk -v org="$PROJECT_NAME" '
            /^Unit size/ { found=1; next }
            /^$/ && found { exit }
            found && /^[1-6][[:space:]]+[0-9]+/ {
                type = ""
                if ($1 == 1) type = "Mono"
                else if ($1 == 2) type = "Di"
                else if ($1 == 3) type = "Tri"
                else if ($1 == 4) type = "Tetra"
                else if ($1 == 5) type = "Penta"
                else if ($1 == 6) type = "Hexa"
                if (type != "") print "\t" org "\t" type "\t" $2
            }
        ' "$stats_file" >> "$dist_csv"
        
        # Check if parsing succeeded
        local dist_count=$(tail -n +2 "$dist_csv" | wc -l)
        if [[ "$dist_count" -eq 0 ]]; then
            log WARN "Failed to parse MISA statistics, using SSR file fallback"
            # Fallback: count from SSR file (but this won't expand compounds)
            local ssr_file="${intermediate_dir}/${PROJECT_NAME}_SSR.txt"
            if [[ -f "$ssr_file" ]]; then
                awk -F'\t' -v org="$PROJECT_NAME" '
                    { types[$5]++ }
                    END { for (t in types) print "\t" org "\t" t "\t" types[t] }
                ' "$ssr_file" >> "$dist_csv"
            fi
        fi
        
        log DEBUG "Created distribution CSV"
    fi
    
    # 2. Frequency of identified SSR motifs - Create simple format CSV
    local freq_csv="${results_dir}/Frequency of the identified SSR motifs.csv"
    local stats_file="${fasta_dir}/${PROJECT_NAME}_genomic.fa.statistics"
    
    if [[ -f "$stats_file" ]]; then
        echo -e "Repeats\ttotal" > "$freq_csv"
        
        # Extract frequency table from MISA statistics (between "Frequency of identified SSR motifs" and next section)
        # Use grep -Pzo for multi-line extraction, then awk to get motif and total
        grep -Pzo "(?s)(?<=Frequency of identified SSR motifs\n)(.*?)(?=(\nFrequency of classified repeat types \(considering sequence complementary\)|\Z))" "$stats_file" | \
            tr '\0' '\n' | \
            awk -F'\t' 'NF > 0 { print $1 "\t" $NF }' | \
            tail -n +3 | \
            grep -v '^$' >> "$freq_csv"
        
        log INFO "Created frequency of SSR motifs CSV"
    else
        log WARN "Statistics file not found: $stats_file"
    fi
    
    # 3. Frequency with complementarity - Extract from different section of MISA statistics
    local comp_csv="${results_dir}/Frequency of identified SSR motifs considering complementarity.csv"
    
    if [[ -f "$stats_file" ]]; then
        echo -e "Repeats\ttotal" > "$comp_csv"
        
        # Extract complementarity frequency table from MISA statistics
        grep -Pzo "(?s)(?<=Frequency of classified repeat types \(considering sequence complementary\)\n)(.*?)(?=(\nFrequency of identified SSR motifs|\Z))" "$stats_file" | \
            tr '\0' '\n' | \
            awk -F'\t' 'NF > 0 { print $1 "\t" $NF }' | \
            tail -n +3 | \
            sed '$d' >> "$comp_csv"
        
        log INFO "Created frequency with complementarity CSV"
    fi
    
    # 4. Generate comparison analyses (genic vs non-genic)
    log INFO "Generating genic vs non-genic comparison analyses"
    
    local genic_file="${intermediate_dir}/${PROJECT_NAME}.Genic_SSR_with_feature.txt"
    local intergenic_file="${intermediate_dir}/${PROJECT_NAME}_Accept_Intergenic_SSR.txt"
    
    if [[ -f "$genic_file" ]] && [[ -f "$intergenic_file" ]]; then
        # Extract BASE motifs matching reference method:
        # (AT)7 → AT, (AT)11gtgg(AT)6 → AT (first motif only)
        # sed 's/).*//' removes everything after first ), then sed 's/.*(//' removes up to last (
        # Genic file: motif is in column 5
        # Intergenic file: motif is in column 6
        awk -F'\t' '{
            if ($5 != "") {
                motif = $5
                # Remove everything after first )
                sub(/\).*/, "", motif)
                # Remove everything up to and including last (
                sub(/.*\(/, "", motif)
                if (motif != "") print motif
            }
        }' "$genic_file" | sort | uniq -c | awk '{print $2 "\t" $1}' > "$temp_dir/genic_motif_counts.txt"
        
        awk -F'\t' '{
            if ($6 != "") {
                motif = $6
                sub(/\).*/, "", motif)
                sub(/.*\(/, "", motif)
                if (motif != "") print motif
            }
        }' "$intergenic_file" | sort | uniq -c | awk '{print $2 "\t" $1}' > "$temp_dir/intergenic_motif_counts.txt"
        
        # Extract unique motifs
        cut -f1 "$temp_dir/genic_motif_counts.txt" | sort -u > "$temp_dir/genic_motifs.txt"
        cut -f1 "$temp_dir/intergenic_motif_counts.txt" | sort -u > "$temp_dir/intergenic_motifs.txt"
        
        # Find shared and unique motifs
        comm -12 "$temp_dir/genic_motifs.txt" "$temp_dir/intergenic_motifs.txt" > "$temp_dir/shared_motifs.txt"
        comm -23 "$temp_dir/genic_motifs.txt" "$temp_dir/intergenic_motifs.txt" > "$temp_dir/unique_genic_motifs.txt"
        comm -13 "$temp_dir/genic_motifs.txt" "$temp_dir/intergenic_motifs.txt" > "$temp_dir/unique_intergenic_motifs.txt"
        
        # Create shared motifs CSV with both counts (using join - much faster)
        local shared_csv="${results_dir}/Motif sequences shared between genic and non-genic regions.csv"
        echo -e "Repeat Sequence\tNO. of the present (Genic)\tNO. of the present (Non-genic)" > "$shared_csv"
        join -t $'\t' -1 1 -2 1 \
            <(sort -t $'\t' -k1,1 "$temp_dir/genic_motif_counts.txt") \
            <(sort -t $'\t' -k1,1 "$temp_dir/intergenic_motif_counts.txt") \
            >> "$shared_csv"
        
        # Create unique genic motifs CSV with reference format (tab-separated with Process ID)
        {
            echo -e "\tProcess Id\tNO. of the present (Genic)\tRepeat Sequence"
            while read -r motif; do
                [[ -z "$motif" ]] && continue
                count=$(grep "^${motif}	" "$temp_dir/genic_motif_counts.txt" | head -1 | cut -f2)
                [[ -z "$count" ]] && count=0
                echo -e "\t${PROJECT_NAME}\t${count}\t${motif}"
            done < "$temp_dir/unique_genic_motifs.txt"
        # NOTE: Filename "Repeate" is intentionally preserved for backward compatibility
        # with existing reference outputs and downstream tools
        } > "${results_dir}/Repeate sequences Unique to genic regions.csv"
        
        # Create unique intergenic motifs CSV with reference format (tab-separated with Process ID)
        # NOTE: Filename "Repeate" is intentionally preserved for backward compatibility
        {
            echo -e "\tProcess Id\tNO. of the present (Non-Genic)\tRepeat Sequence"
            while read -r motif; do
                [[ -z "$motif" ]] && continue
                count=$(grep "^${motif}	" "$temp_dir/intergenic_motif_counts.txt" | head -1 | cut -f2)
                [[ -z "$count" ]] && count=0
                echo -e "\t${PROJECT_NAME}\t${count}\t${motif}"
            done < "$temp_dir/unique_intergenic_motifs.txt"
        } > "${results_dir}/Repeate sequences Unique to non-genic regions.csv"
        
        log INFO "Created motif comparison CSVs"
    fi
    
    # 5. Primer statistics
    local genic_primers="${results_dir}/${PROJECT_NAME}.Genic-primers.txt"
    local intergenic_primers="${results_dir}/${PROJECT_NAME}.interGenic-primers.txt"
    
    if [[ -f "$genic_primers" ]] && [[ -f "$genic_file" ]]; then
        local total_genic=$(wc -l < "$genic_file")
        local designed_genic=$(wc -l < "$genic_primers")
        # Subtract 1 for header
        designed_genic=$((designed_genic - 1))
        local success_rate=$(awk "BEGIN {printf \"%.2f\", ($designed_genic * 100.0) / $total_genic}")
        
        echo "Metric,Value" > "${results_dir}/Genic SSR primer statistics.csv"
        echo "Total Genic SSRs,$total_genic" >> "${results_dir}/Genic SSR primer statistics.csv"
        echo "Primers Designed,$designed_genic" >> "${results_dir}/Genic SSR primer statistics.csv"
        echo "Success Rate (%),$success_rate" >> "${results_dir}/Genic SSR primer statistics.csv"
    fi
    
    # 6. Non-genic primer statistics
    if [[ -f "$intergenic_primers" ]] && [[ -f "$intergenic_file" ]]; then
        local total_intergenic=$(wc -l < "$intergenic_file")
        local designed_intergenic=$(wc -l < "$intergenic_primers")
        # Subtract 1 for header
        designed_intergenic=$((designed_intergenic - 1))
        local success_rate=$(awk "BEGIN {printf \"%.2f\", ($designed_intergenic * 100.0) / $total_intergenic}")
        
        echo "Metric,Value" > "${results_dir}/Statistics of non-genic SSR primers.csv"
        echo "Total Non-genic SSRs,$total_intergenic" >> "${results_dir}/Statistics of non-genic SSR primers.csv"
        echo "Primers Designed,$designed_intergenic" >> "${results_dir}/Statistics of non-genic SSR primers.csv"
        echo "Success Rate (%),$success_rate" >> "${results_dir}/Statistics of non-genic SSR primers.csv"
        
        log INFO "Created primer statistics CSVs"
    fi
    
    # 7. Summary of identified SSR motifs
    local summary_csv="${results_dir}/Summary of identified SSR motifs.csv"
    echo -e "Category\tCount" > "$summary_csv"
    
    local total_ssr=$(wc -l < "${intermediate_dir}/${PROJECT_NAME}_SSR.txt" 2>/dev/null || echo "0")
    local genic_count=$(wc -l < "$genic_file" 2>/dev/null || echo "0")
    local intergenic_count=$(wc -l < "$intergenic_file" 2>/dev/null || echo "0")
    local genic_primers_count=$(($(wc -l < "$genic_primers" 2>/dev/null || echo "1") - 1))
    local intergenic_primers_count=$(($(wc -l < "$intergenic_primers" 2>/dev/null || echo "1") - 1))
    
    echo -e "Total SSRs\t$total_ssr" >> "$summary_csv"
    echo -e "Genic SSRs\t$genic_count" >> "$summary_csv"
    echo -e "Non-genic SSRs\t$intergenic_count" >> "$summary_csv"
    echo -e "Genic Primers\t$genic_primers_count" >> "$summary_csv"
    echo -e "Non-genic Primers\t$intergenic_primers_count" >> "$summary_csv"
    
    log INFO "Created summary CSV"
    log INFO "Advanced statistics generation complete"
}

#=============================================================================
# Phase 8: Visualizations
#=============================================================================
generate_visualizations() {
    log INFO "=========================================="
    log INFO "PHASE 8: Visualizations"
    log INFO "=========================================="
    
    local results_dir="${OUTPUT_DIR}/${PROJECT_NAME}-MegaSSR_Results"
    local intermediate_dir="${OUTPUT_DIR}/intermediate"
    local temp_dir="${OUTPUT_DIR}/temp_stats"
    
    # Check if matplotlib is available
    if ! python3 -c "import matplotlib" 2>/dev/null; then
        log WARN "matplotlib not installed, skipping visualizations"
        log WARN "Install with: conda install matplotlib or pip install matplotlib"
        return 0
    fi
    
    log INFO "Generating comprehensive visualization suite"
    
    # Create visualization script inline to ensure all data is processed correctly
    python3 -W ignore - "$results_dir" "$intermediate_dir" "$PROJECT_NAME" "$temp_dir" <<'VISUALIZATION_SCRIPT'
#!/usr/bin/env python3
"""MegaSSR Comprehensive Visualization Suite"""

import sys
import os
import warnings
warnings.filterwarnings('ignore')

# Set non-interactive backend before importing pyplot
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np

# Optional imports
try:
    import seaborn as sns
    sns.set_style("whitegrid")
    HAS_SEABORN = True
except ImportError:
    HAS_SEABORN = False

try:
    from matplotlib_venn import venn2
    HAS_VENN = True
except ImportError:
    HAS_VENN = False

def parse_distribution_file(filepath):
    """Parse SSR class distribution data"""
    data = {}
    type_map = {'1': 'Mono', '2': 'Di', '3': 'Tri', '4': 'Tetra', '5': 'Penta', '6': 'Hexa', 'c': 'Compound'}
    
    try:
        with open(filepath, 'r') as f:
            for line in f:
                line = line.strip()
                if not line or line.startswith(('Process', '\t', '#')):
                    continue
                parts = line.split('\t')
                if len(parts) >= 2:
                    # Handle format: Type\tCount or OrganismID\tType\tCount
                    if len(parts) == 2:
                        type_name, count = parts[0], parts[1]
                    elif len(parts) >= 3:
                        type_name, count = parts[-2], parts[-1]
                    
                    # Map numeric types to names
                    type_name = type_map.get(type_name, type_name)
                    
                    try:
                        data[type_name] = int(count)
                    except ValueError:
                        pass
    except Exception as e:
        print(f"Warning: Could not parse distribution file: {e}")
    
    return data

def parse_frequency_file(filepath):
    """Parse motif frequency data - handles both CSV (Motif,Count) and tab format (\\tProcessId\\tCount\\tMotif)"""
    data = {}
    try:
        with open(filepath, 'r') as f:
            for line in f:
                line_stripped = line.strip()
                # Skip empty lines and headers (check stripped line for headers)
                if not line_stripped or line_stripped.startswith(('Repeat', 'motif', '#', 'Process', 'NO.')):
                    continue
                
                # Try tab-delimited first (reference format: \\tProcessId\\tCount\\tMotif)
                parts = line.rstrip('\n').split('\t')  # Don't strip leading tab
                if len(parts) == 4 and not parts[0]:  # Tab format with leading empty field
                    # Format: \\tProcessId\\tCount\\tMotif (fields: ['', 'ProcessId', 'Count', 'Motif'])
                    try:
                        motif = parts[3].strip()
                        count = int(parts[2])
                        if motif and count > 0:
                            data[motif] = count
                    except (ValueError, IndexError):
                        pass
                elif len(parts) >= 2:
                    # CSV format: Motif,Count or simple tab format
                    parts = line_stripped.replace(',', '\t').split('\t')
                    motif = parts[0].strip()
                    try:
                        count = int(parts[-1])
                        if motif and count > 0:
                            data[motif] = count
                    except ValueError:
                        pass
    except Exception as e:
        print(f"Warning: Could not parse frequency file: {e}")
    return data

def parse_comparison_file(filepath):
    """Parse genic vs non-genic comparison data"""
    genic = {}
    nongenic = {}
    try:
        with open(filepath, 'r') as f:
            for line in f:
                line = line.strip()
                if not line or line.startswith(('Repeat', 'motif', '#')):
                    continue
                parts = line.replace(',', '\t').split('\t')
                if len(parts) >= 3:
                    motif = parts[0].strip()
                    try:
                        genic[motif] = int(parts[1])
                        nongenic[motif] = int(parts[2])
                    except (ValueError, IndexError):
                        pass
    except Exception as e:
        print(f"Warning: Could not parse comparison file: {e}")
    return genic, nongenic

def plot_distribution(data, output_path):
    """Create SSR class distribution bar chart"""
    if not data:
        print("No distribution data available")
        return False
    
    # Define order
    order = ['Mono', 'Di', 'Tri', 'Tetra', 'Penta', 'Hexa', 'Compound']
    categories = [c for c in order if c in data]
    values = [data[c] for c in categories]
    
    if not categories:
        print("No valid distribution categories found")
        return False
    
    fig, ax = plt.subplots(figsize=(10, 6))
    colors = plt.cm.viridis(np.linspace(0.2, 0.8, len(categories)))
    bars = ax.bar(categories, values, color=colors, edgecolor='black', linewidth=0.8)
    
    # Add value labels on bars
    for bar, val in zip(bars, values):
        height = bar.get_height()
        ax.annotate(f'{val:,}',
                    xy=(bar.get_x() + bar.get_width() / 2, height),
                    xytext=(0, 3), textcoords="offset points",
                    ha='center', va='bottom', fontsize=10, fontweight='bold')
    
    ax.set_xlabel('SSR Type', fontsize=12, fontweight='bold')
    ax.set_ylabel('Count', fontsize=12, fontweight='bold')
    ax.set_title('Distribution of Different SSR Classes', fontsize=14, fontweight='bold')
    
    plt.tight_layout()
    plt.savefig(output_path, dpi=150, bbox_inches='tight', facecolor='white')
    plt.close()
    print(f"Created: {output_path}")
    return True

def plot_frequency(data, output_path, title="Frequency of Identified SSR Motifs", top_n=60):
    """Create motif frequency bar chart"""
    if not data:
        print("No frequency data available")
        return False
    
    # Sort by frequency and get top N
    sorted_data = sorted(data.items(), key=lambda x: x[1], reverse=True)[:top_n]
    motifs = [x[0] for x in sorted_data]
    counts = [x[1] for x in sorted_data]
    
    fig, ax = plt.subplots(figsize=(16, 8))
    colors = plt.cm.tab20c(np.linspace(0, 1, len(motifs)))
    bars = ax.bar(range(len(motifs)), counts, color=colors, edgecolor='black', linewidth=0.5)
    
    # Add value labels on top of bars
    for i, (bar, count) in enumerate(zip(bars, counts)):
        height = bar.get_height()
        ax.annotate(f'{count}',
                    xy=(bar.get_x() + bar.get_width() / 2, height),
                    xytext=(0, 2), textcoords="offset points",
                    ha='center', va='bottom', fontsize=7, fontweight='bold',
                    rotation=90)
    
    ax.set_xlabel('Repeat', fontsize=12, fontweight='bold')
    ax.set_ylabel('Count', fontsize=12, fontweight='bold')
    ax.set_title(title + f'\n*Only the highest {top_n} values are shown', 
                 fontsize=14, fontweight='bold')
    ax.set_xticks(range(len(motifs)))
    ax.set_xticklabels(motifs, rotation=90, fontsize=8, ha='center')
    ax.grid(axis='y', alpha=0.3, linestyle='--')
    
    plt.tight_layout()
    plt.savefig(output_path, dpi=150, bbox_inches='tight', facecolor='white')
    plt.close()
    print(f"Created: {output_path}")
    return True

def plot_comparison_grouped(genic_data, nongenic_data, output_path):
    """Create grouped bar chart for genic vs non-genic comparison"""
    if not genic_data and not nongenic_data:
        print("No comparison data available")
        return False
    
    # Get all motifs and sort by total count (use top 60 to match reference)
    all_motifs = set(genic_data.keys()) | set(nongenic_data.keys())
    motif_totals = {m: genic_data.get(m, 0) + nongenic_data.get(m, 0) for m in all_motifs}
    sorted_motifs = sorted(motif_totals.keys(), key=lambda x: motif_totals[x], reverse=True)[:60]
    
    genic_vals = [genic_data.get(m, 0) for m in sorted_motifs]
    nongenic_vals = [nongenic_data.get(m, 0) for m in sorted_motifs]
    
    x = np.arange(len(sorted_motifs))
    width = 0.38  # Optimal width for grouped bars - wide enough to be visible, narrow enough to not overlap
    
    fig, ax = plt.subplots(figsize=(20, 10))
    bars1 = ax.bar(x - width/2, genic_vals, width, label='Genic', color='#4472C4', edgecolor='black', linewidth=0.5)
    bars2 = ax.bar(x + width/2, nongenic_vals, width, label='Non Genic', color='#ED7D31', edgecolor='black', linewidth=0.5)
    
    # Add value labels on top of bars with 90 degree rotation
    for bar in bars1:
        height = bar.get_height()
        if height > 0:
            ax.annotate(f'{int(height)}',
                       xy=(bar.get_x() + bar.get_width() / 2, height),
                       xytext=(0, 3), textcoords="offset points",
                       ha='center', va='bottom', fontsize=7, fontweight='bold',
                       rotation=90)
    
    for bar in bars2:
        height = bar.get_height()
        if height > 0:
            ax.annotate(f'{int(height)}',
                       xy=(bar.get_x() + bar.get_width() / 2, height),
                       xytext=(0, 3), textcoords="offset points",
                       ha='center', va='bottom', fontsize=7, fontweight='bold',
                       rotation=90)
    
    ax.set_xlabel('Repeat', fontsize=12, fontweight='bold')
    ax.set_ylabel('Count', fontsize=12, fontweight='bold')
    ax.set_title('Shared repeats between genic and non-genic regions\n*Only the highest 60 values are shown', 
                 fontsize=14, fontweight='bold', color='red')
    ax.set_xticks(x)
    ax.set_xticklabels(sorted_motifs, rotation=90, ha='center', fontsize=9)
    ax.legend(loc='upper right', fontsize=11)
    ax.grid(axis='y', alpha=0.3, linestyle='--')
    
    plt.tight_layout()
    plt.savefig(output_path, dpi=150, bbox_inches='tight', facecolor='white')
    plt.close()
    print(f"Created: {output_path}")
    return True

def plot_venn(genic_only, nongenic_only, shared_count, output_path):
    """Create Venn diagram for genic vs non-genic SSRs
    
    Args:
        genic_only: Count of motifs unique to genic regions (already exclusive)
        nongenic_only: Count of motifs unique to non-genic regions (already exclusive)
        shared_count: Count of motifs shared between both regions
    """
    if not HAS_VENN:
        print("matplotlib_venn not available, skipping Venn diagram")
        return False
    
    fig, ax = plt.subplots(figsize=(10, 8))
    
    # Use counts directly - they are already unique (no overlap subtraction needed)
    venn = venn2(subsets=(genic_only, nongenic_only, shared_count),
                 set_labels=('Genic region', 'Non-genic region'),
                 set_colors=('#9370DB', '#5DADE2'),
                 alpha=0.7,
                 ax=ax)
    
    # Style the diagram
    for text in venn.set_labels:
        if text:
            text.set_fontsize(14)
            text.set_fontweight('bold')
    for text in venn.subset_labels:
        if text:
            text.set_fontsize(14)
            text.set_fontweight('bold')
    
    ax.set_title('Venn diagram showing a comparison of genic and non-genic SSR', 
                 fontsize=14, fontweight='normal')
    
    plt.tight_layout()
    plt.savefig(output_path, dpi=150, bbox_inches='tight', facecolor='white')
    plt.close()
    print(f"Created: {output_path}")
    return True

def plot_unique_motifs(data, output_path, title, color):
    """Create vertical bar chart for unique motifs"""
    if not data:
        print(f"No unique motif data for: {title}")
        return False
    
    # Sort and get top 60 (matching frequency plot limit)
    sorted_data = sorted(data.items(), key=lambda x: x[1], reverse=True)[:60]
    motifs = [x[0] for x in sorted_data]
    counts = [x[1] for x in sorted_data]
    
    # Adjust figure size for vertical bars
    fig, ax = plt.subplots(figsize=(16, 10))
    x_pos = np.arange(len(motifs))
    bars = ax.bar(x_pos, counts, color=color, edgecolor='black', linewidth=0.5)
    
    ax.set_xticks(x_pos)
    ax.set_xticklabels(motifs, fontsize=8, rotation=90)
    ax.set_ylabel('Count', fontsize=12, fontweight='bold')
    ax.set_xlabel('Repeat', fontsize=12, fontweight='bold')
    ax.set_title(title + '\n*Only the highest 60 values are shown', fontsize=14, fontweight='bold')
    
    # Add value labels on top of bars
    for bar, val in zip(bars, counts):
        ax.annotate(f'{val}', xy=(bar.get_x() + bar.get_width()/2, bar.get_height()),
                    xytext=(0, 3), textcoords="offset points",
                    ha='center', va='bottom', fontsize=9)
    
    plt.tight_layout()
    plt.savefig(output_path, dpi=150, bbox_inches='tight', facecolor='white')
    plt.close()
    print(f"Created: {output_path}")
    return True

def main():
    if len(sys.argv) < 5:
        print("Usage: script.py results_dir intermediate_dir project_name temp_dir")
        sys.exit(1)
    
    results_dir = sys.argv[1]
    intermediate_dir = sys.argv[2]
    project_name = sys.argv[3]
    temp_dir = sys.argv[4]
    
    print(f"Generating visualizations in: {results_dir}")
    
    # Pre-flight check: verify required CSV files exist
    plots_skipped = []
    
    # 1. Distribution plot
    dist_file = os.path.join(results_dir, "The distribution of the different SSR classes.csv")
    if os.path.exists(dist_file):
        data = parse_distribution_file(dist_file)
        if data:
            plot_distribution(data, os.path.join(results_dir, "Distribution of the different SSR classes.png"))
        else:
            plots_skipped.append(("Distribution plot", "Empty or invalid data"))
    else:
        plots_skipped.append(("Distribution plot", f"Missing: {os.path.basename(dist_file)}"))
    
    # 2. Frequency plot
    freq_file = os.path.join(results_dir, "Frequency of the identified SSR motifs.csv")
    if os.path.exists(freq_file):
        data = parse_frequency_file(freq_file)
        if data:
            plot_frequency(data, os.path.join(results_dir, "Frequency of the identified SSR motifs.png"))
        else:
            plots_skipped.append(("Frequency plot", "Empty or invalid data"))
    else:
        plots_skipped.append(("Frequency plot", f"Missing: {os.path.basename(freq_file)}"))
    
    # 3. Complementarity frequency plot
    comp_file = os.path.join(results_dir, "Frequency of identified SSR motifs considering complementarity.csv")
    if os.path.exists(comp_file):
        data = parse_frequency_file(comp_file)
        if data:
            plot_frequency(data, os.path.join(results_dir, "SSR distribution considering sequence complementarity.png"),
                          title="SSR Distribution Considering Sequence Complementarity")
        else:
            plots_skipped.append(("Complementarity plot", "Empty or invalid data"))
    else:
        plots_skipped.append(("Complementarity plot", f"Missing: {os.path.basename(comp_file)}"))
    
    # 4. Venn diagram - Use motif counts from CSV files (not SSR counts)
    # NOTE: Filenames preserve legacy naming ("Repeate" typo) for reference compatibility
    unique_genic_csv = os.path.join(results_dir, "Repeate sequences Unique to genic regions.csv")
    unique_nongenic_csv = os.path.join(results_dir, "Repeate sequences Unique to non-genic regions.csv")
    shared_csv = os.path.join(results_dir, "Motif sequences shared between genic and non-genic regions.csv")
    
    venn_files_exist = all(os.path.exists(f) for f in [unique_genic_csv, unique_nongenic_csv, shared_csv])
    if venn_files_exist:
        # Count unique motifs (subtract 1 for header)
        with open(unique_genic_csv, 'r') as f:
            genic_only = sum(1 for _ in f) - 1
        with open(unique_nongenic_csv, 'r') as f:
            nongenic_only = sum(1 for _ in f) - 1
        with open(shared_csv, 'r') as f:
            shared_motifs = sum(1 for _ in f) - 1
        
        if HAS_VENN:
            plot_venn(genic_only, nongenic_only, shared_motifs,
                     os.path.join(results_dir, "Venn diagram of genic and non-genic SSR.png"))
        else:
            plots_skipped.append(("Venn diagram", "matplotlib_venn not installed"))
    else:
        missing = [os.path.basename(f) for f in [unique_genic_csv, unique_nongenic_csv, shared_csv] if not os.path.exists(f)]
        plots_skipped.append(("Venn diagram", f"Missing: {', '.join(missing)}"))
    
    # 5. Common repeats (use shared motifs data for bar chart)
    shared_file = os.path.join(results_dir, "Motif sequences shared between genic and non-genic regions.csv")
    if os.path.exists(shared_file):
        genic_data, nongenic_data = parse_comparison_file(shared_file)
        if genic_data or nongenic_data:
            plot_comparison_grouped(genic_data, nongenic_data, 
                                   os.path.join(results_dir, "Common repeats between genic and non-genic regions.png"))
        else:
            plots_skipped.append(("Common repeats plot", "Empty or invalid data"))
    else:
        plots_skipped.append(("Common repeats plot", f"Missing: {os.path.basename(shared_file)}"))
    
    # 6. Unique genic motifs
    # NOTE: Filename preserves legacy naming ("Repeate" typo) for reference compatibility
    unique_genic_file = os.path.join(results_dir, "Repeate sequences Unique to genic regions.csv")
    if os.path.exists(unique_genic_file):
        data = parse_frequency_file(unique_genic_file)
        if data:
            plot_unique_motifs(data, os.path.join(results_dir, "Unique repeats in the genic region.png"),
                              "Unique Repeats in the Genic Region", '#3498db')
        else:
            plots_skipped.append(("Unique genic motifs plot", "Empty or invalid data"))
    else:
        plots_skipped.append(("Unique genic motifs plot", f"Missing: {os.path.basename(unique_genic_file)}"))
    
    # 7. Unique non-genic motifs
    # NOTE: Filename preserves legacy naming ("Repeate" typo) for reference compatibility
    unique_nongenic_file = os.path.join(results_dir, "Repeate sequences Unique to non-genic regions.csv")
    if os.path.exists(unique_nongenic_file):
        data = parse_frequency_file(unique_nongenic_file)
        if data:
            plot_unique_motifs(data, os.path.join(results_dir, "Unique repeats of the non-genic region.png"),
                              "Unique Repeats of the Non-genic Region", '#e74c3c')
        else:
            plots_skipped.append(("Unique non-genic motifs plot", "Empty or invalid data"))
    else:
        plots_skipped.append(("Unique non-genic motifs plot", f"Missing: {os.path.basename(unique_nongenic_file)}"))
    
    # Print summary of skipped plots
    if plots_skipped:
        print(f"\nWARNING: {len(plots_skipped)} plot(s) skipped due to missing data:")
        for plot_name, reason in plots_skipped:
            print(f"  - {plot_name}: {reason}")
    
    print("Visualization generation complete!")

if __name__ == "__main__":
    main()
VISUALIZATION_SCRIPT

    if [[ $? -eq 0 ]]; then
        log INFO "Visualization generation complete"
    else
        log WARN "Some visualizations may have failed - check output files"
    fi
    
    # Clean up temporary directory
    rm -rf "$temp_dir"
}

#=============================================================================
# Analysis Type 1 Specific Functions
# These functions handle the simplified pipeline without gene annotation
#=============================================================================

#=============================================================================
# Phase 6 (Type 1): Post-processing for SSR-only analysis
#=============================================================================
run_post_processing_type1() {
    log INFO "=========================================="
    log INFO "PHASE 6: Post-processing (Type 1)"
    log INFO "=========================================="
    
    local intermediate_dir="${OUTPUT_DIR}/intermediate"
    local results_dir="${OUTPUT_DIR}/${PROJECT_NAME}-MegaSSR_Results"
    
    # Generate SSR statistics from MISA output
    generate_ssr_statistics
    
    # Create non-redundant SSR library (intergenic only for Type 1)
    create_nonredundant_library_type1
    
    # Export final CSVs with friendly names
    export_csv_files_type1
    
    log INFO "Post-processing complete"
}

create_nonredundant_library_type1() {
    log INFO "Creating non-redundant SSR library (Type 1)"
    
    local results_dir="${OUTPUT_DIR}/${PROJECT_NAME}-MegaSSR_Results"
    local intergenic_fa="${results_dir}/Intergenic_SSR_flanking_regions.fa"
    local nonred_fa="${results_dir}/Non-redundant_SSR_library.fasta"
    local nonred_log="${results_dir}/SSR.non-redundant.log"
    
    # Check if flanking region file exists
    if [[ ! -f "$intergenic_fa" ]] || [[ ! -s "$intergenic_fa" ]]; then
        log WARN "No flanking region file found, skipping non-redundant library"
        return 0
    fi
    
    local seq_count=$(grep -c "^>" "$intergenic_fa" || echo "0")
    log INFO "Processing $seq_count flanking sequences"
    
    # Define usearch binary path
    local usearch_bin="${USEARCH_BIN:-${BIN_DIR}/usearch11.0.667_i86linux32}"
    
    # Try to use usearch for clustering
    local usearch_cmd=""
    if [[ -n "${USEARCH_BIN:-}" ]] && [[ -x "$USEARCH_BIN" ]]; then
        usearch_cmd="$USEARCH_BIN"
        log INFO "Using custom usearch from USEARCH_BIN: $USEARCH_BIN"
    elif [[ -x "$usearch_bin" ]]; then
        if "$usearch_bin" --version >/dev/null 2>&1; then
            usearch_cmd="$usearch_bin"
            log INFO "Using bundled usearch for clustering"
        else
            log WARN "Bundled usearch binary not compatible with this architecture"
        fi
    fi
    
    # Fallback to system usearch
    if [[ -z "$usearch_cmd" ]] && command -v usearch >/dev/null 2>&1; then
        usearch_cmd="usearch"
        log INFO "Using system usearch for clustering"
    fi
    
    if [[ -n "$usearch_cmd" ]]; then
        local temp_dir="${OUTPUT_DIR}/usearch_temp"
        mkdir -p "$temp_dir"
        
        # Step 1: Sort by length
        local sorted_fa="${temp_dir}/SSR_with_flanking_regions_sorted.fa"
        "$usearch_cmd" -sortbylength "$intergenic_fa" \
            -fastaout "$sorted_fa" \
            -log "${temp_dir}/usearch_sort.log" 2>&1
        
        # Step 2: Cluster at 90% identity
        "$usearch_cmd" -cluster_fast "$sorted_fa" \
            -id 0.90 \
            -centroids "${temp_dir}/centroids.fa" \
            -uc "${temp_dir}/result.uc" \
            -consout "$nonred_fa" \
            -msaout "${temp_dir}/aligned.fasta" \
            -log "${temp_dir}/usearch_cluster.log" 2>&1
        
        if [[ $? -eq 0 ]] && [[ -f "$nonred_fa" ]]; then
            sed '2d;3d;11d' "${temp_dir}/usearch_cluster.log" > "$nonred_log" 2>/dev/null || \
                cp "${temp_dir}/usearch_cluster.log" "$nonred_log"
            
            local cluster_count=$(grep -c "^>" "$nonred_fa" || echo "0")
            log INFO "Created non-redundant library: $cluster_count clusters (from $seq_count sequences)"
            rm -rf "$temp_dir"
            return 0
        else
            log WARN "usearch clustering failed, falling back to simple copy"
            rm -rf "$temp_dir"
        fi
    fi
    
    # Fallback: just copy the flanking sequences
    log INFO "Using flanking sequences as non-redundant library (usearch not available)"
    cp "$intergenic_fa" "$nonred_fa"
    echo "Non-redundant library: $seq_count sequences (no clustering applied)" > "$nonred_log"
}

export_csv_files_type1() {
    log INFO "Exporting CSV files (Type 1)"
    
    local results_dir="${OUTPUT_DIR}/${PROJECT_NAME}-MegaSSR_Results"
    local intermediate_dir="${OUTPUT_DIR}/intermediate"
    
    # POSIX-safe sed -i replacement function
    sed_inplace() {
        local expr="$1"
        local file="$2"
        local temp_file=$(mktemp)
        sed "$expr" "$file" > "$temp_file" && mv "$temp_file" "$file"
    }
    
    # Add headers to intergenic primers file (renamed to "SSR primers" for Type 1)
    local intergenic_primers="${results_dir}/${PROJECT_NAME}.interGenic-primers.txt"
    if [[ -f "$intergenic_primers" ]] && [[ -s "$intergenic_primers" ]]; then
        sed_inplace '1s/^/\tSequence ID\tProcess Id\tRepeat number\tRepeat Type\tRepeat Sequence\tRepeat Length\tRepeat Start\tRepeat End\tPrimer Start\tPrimer End\tForward Primer\tForward Tm\tForward Size (bp)\tForward GC\tReverse Primer\tReverse Tm\tReverse Size (bp)\tReverse GC\tProduct Size (bp)\n/' "$intergenic_primers"
        
        # Create renamed version matching old reference format
        cp "$intergenic_primers" "${results_dir}/Desinged SSR primer.csv"
        log INFO "Exported SSR primers CSV with header"
    fi
    
    # Copy SSR list as "Accept_Intergenic_SSR" with header (matches reference)
    local ssr_file="${intermediate_dir}/${PROJECT_NAME}_SSR.txt"
    if [[ -f "$ssr_file" ]] && [[ -s "$ssr_file" ]]; then
        local ssr_csv="${results_dir}/${PROJECT_NAME}.Accept_Intergenic_SSR.txt"
        sed '1s/^/\tProcess Id\tSequence Id\tRepeat number\tRepeat Type\tRepeat Sequence\tRepeat Length\tRepeat Start\tRepeat End\n/' "$ssr_file" > "$ssr_csv"
        log INFO "Exported SSR list CSV with header"
    fi
    
    # Rename flanking regions to match reference
    if [[ -f "${results_dir}/Intergenic_SSR_flanking_regions.fa" ]]; then
        cp "${results_dir}/Intergenic_SSR_flanking_regions.fa" "${results_dir}/SSR flanking regions.fa"
    fi
    
    # Rename non-redundant library to match reference
    if [[ -f "${results_dir}/Non-redundant_SSR_library.fasta" ]]; then
        cp "${results_dir}/Non-redundant_SSR_library.fasta" "${results_dir}/Non-redundant SSR library.fasta"
    fi
}

#=============================================================================
# Phase 7 (Type 1): Statistics for SSR-only analysis
#=============================================================================
generate_statistics_type1() {
    log INFO "=========================================="
    log INFO "PHASE 7: Statistics (Type 1)"
    log INFO "=========================================="
    
    local results_dir="${OUTPUT_DIR}/${PROJECT_NAME}-MegaSSR_Results"
    local intermediate_dir="${OUTPUT_DIR}/intermediate"
    local fasta_dir="${OUTPUT_DIR}/fasta"
    local stats_file="${fasta_dir}/${PROJECT_NAME}_genomic.fa.statistics"
    
    # 1. Distribution of different repeat type classes
    if [[ -f "$stats_file" ]]; then
        local dist_csv="${results_dir}/The distribution of the different SSR classes.csv"
        echo -e "\tProcess Id\tRepeat Type\tTotal No. of present" > "$dist_csv"
        
        awk -v org="$PROJECT_NAME" '
            /^Unit size/ { found=1; next }
            /^$/ && found { exit }
            found && /^[1-6][[:space:]]+[0-9]+/ {
                type = ""
                if ($1 == 1) type = "Mono"
                else if ($1 == 2) type = "Di"
                else if ($1 == 3) type = "Tri"
                else if ($1 == 4) type = "Tetra"
                else if ($1 == 5) type = "Penta"
                else if ($1 == 6) type = "Hexa"
                if (type != "") print "\t" org "\t" type "\t" $2
            }
        ' "$stats_file" >> "$dist_csv"
        
        log INFO "Created distribution CSV"
    fi
    
    # 2. Frequency of identified SSR motifs
    local freq_csv="${results_dir}/Frequency of the identified SSR motifs.csv"
    if [[ -f "$stats_file" ]]; then
        echo -e "Repeats\ttotal" > "$freq_csv"
        
        grep -Pzo "(?s)(?<=Frequency of identified SSR motifs\n)(.*?)(?=(\nFrequency of classified repeat types \(considering sequence complementary\)|\Z))" "$stats_file" | \
            tr '\0' '\n' | \
            awk -F'\t' 'NF > 0 { print $1 "\t" $NF }' | \
            tail -n +3 | \
            grep -v '^$' >> "$freq_csv"
        
        log INFO "Created frequency of SSR motifs CSV"
    fi
    
    # 3. Frequency with complementarity
    local comp_csv="${results_dir}/Frequency of identified SSR motifs considering complementarity.csv"
    if [[ -f "$stats_file" ]]; then
        echo -e "Repeats\ttotal" > "$comp_csv"
        
        grep -Pzo "(?s)(?<=Frequency of classified repeat types \(considering sequence complementary\)\n)(.*?)(?=(\nFrequency of identified SSR motifs|\Z))" "$stats_file" | \
            tr '\0' '\n' | \
            awk -F'\t' 'NF > 0 { print $1 "\t" $NF }' | \
            tail -n +3 | \
            sed '$d' >> "$comp_csv"
        
        log INFO "Created frequency with complementarity CSV"
    fi
    
    # 4. Summary of identified SSR motifs
    local summary_csv="${results_dir}/Summary of identified SSR motifs.csv"
    local ssr_file="${intermediate_dir}/${PROJECT_NAME}_SSR.txt"
    local primers_file="${results_dir}/${PROJECT_NAME}.interGenic-primers.txt"
    
    local total_ssr=$(wc -l < "$ssr_file" 2>/dev/null || echo "0")
    local primers_count=0
    if [[ -f "$primers_file" ]]; then
        primers_count=$(($(wc -l < "$primers_file") - 1))  # Subtract header
    fi
    
    echo -e "Category\tCount" > "$summary_csv"
    echo -e "Total SSRs\t$total_ssr" >> "$summary_csv"
    echo -e "SSR Primers Designed\t$primers_count" >> "$summary_csv"
    
    # 5. Primer statistics
    if [[ -f "$primers_file" ]] && [[ $total_ssr -gt 0 ]]; then
        local success_rate=$(awk "BEGIN {printf \"%.2f\", ($primers_count * 100.0) / $total_ssr}")
        
        local primer_stats="${results_dir}/SSR primer statistics.csv"
        echo "Metric,Value" > "$primer_stats"
        echo "Total SSRs,$total_ssr" >> "$primer_stats"
        echo "Primers Designed,$primers_count" >> "$primer_stats"
        echo "Success Rate (%),$success_rate" >> "$primer_stats"
        
        log INFO "Created primer statistics CSV"
    fi
    
    log INFO "Statistics generation complete"
}

#=============================================================================
# Phase 8 (Type 1): Visualizations (3 plots)
#=============================================================================
generate_visualizations_type1() {
    log INFO "=========================================="
    log INFO "PHASE 8: Visualizations (Type 1 - 3 plots)"
    log INFO "=========================================="
    
    local results_dir="${OUTPUT_DIR}/${PROJECT_NAME}-MegaSSR_Results"
    
    # Check if matplotlib is available
    if ! python3 -c "import matplotlib" 2>/dev/null; then
        log WARN "matplotlib not installed, skipping visualizations"
        log WARN "Install with: conda install matplotlib or pip install matplotlib"
        return 0
    fi
    
    log INFO "Generating visualizations for Type 1 analysis"
    
    # Create visualization script for Type 1 (3 plots only)
    python3 -W ignore - "$results_dir" "$PROJECT_NAME" <<'VISUALIZATION_TYPE1_SCRIPT'
#!/usr/bin/env python3
"""MegaSSR Type 1 Visualization Suite (3 plots)"""

import sys
import os
import warnings
warnings.filterwarnings('ignore')

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np

def parse_distribution_csv(filepath):
    """Parse SSR class distribution CSV"""
    data = {}
    try:
        with open(filepath, 'r') as f:
            for line in f:
                line = line.strip()
                if not line or 'Process Id' in line:
                    continue
                parts = line.split('\t')
                # Format after strip(): ProcessId\tType\tCount (3 fields)
                if len(parts) >= 3:
                    type_name = parts[1]
                    try:
                        count = int(parts[2])
                        data[type_name] = count
                    except ValueError:
                        pass
    except Exception as e:
        print(f"Warning: Could not parse distribution file: {e}")
    return data

def parse_frequency_csv(filepath):
    """Parse motif frequency CSV"""
    data = {}
    try:
        with open(filepath, 'r') as f:
            for line in f:
                line = line.strip()
                if not line or line.startswith('Repeats'):
                    continue
                parts = line.split('\t')
                if len(parts) >= 2:
                    motif = parts[0]
                    try:
                        count = int(parts[1])
                        if motif and count > 0:
                            data[motif] = count
                    except ValueError:
                        pass
    except Exception as e:
        print(f"Warning: Could not parse frequency file: {e}")
    return data

def plot_distribution(data, output_path):
    """Plot SSR class distribution"""
    if not data:
        print("Warning: No data for distribution plot")
        return
    
    # Order by standard SSR classes
    order = ['Mono', 'Di', 'Tri', 'Tetra', 'Penta', 'Hexa', 'Compound']
    labels = [k for k in order if k in data]
    values = [data[k] for k in labels]
    
    if not labels:
        print("Warning: No valid SSR classes found")
        return
    
    colors = ['#3498db', '#e74c3c', '#2ecc71', '#f39c12', '#9b59b6', '#1abc9c', '#34495e']
    
    fig, ax = plt.subplots(figsize=(10, 6))
    bars = ax.bar(labels, values, color=colors[:len(labels)], edgecolor='black', linewidth=1.2)
    
    ax.set_xlabel('SSR Type', fontsize=12, fontweight='bold')
    ax.set_ylabel('Count', fontsize=12, fontweight='bold')
    ax.set_title('Distribution of the Different SSR Classes', fontsize=14, fontweight='bold')
    
    # Add value labels on bars (no rotation for better visibility)
    max_val = max(values) if values else 1
    for bar, val in zip(bars, values):
        ax.text(bar.get_x() + bar.get_width()/2, bar.get_height(),
                str(val), ha='center', va='bottom', fontsize=11, fontweight='bold', rotation=0)
    
    # Add margin at top for labels
    ax.set_ylim(0, max_val * 1.1)
    
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    
    plt.tight_layout()
    plt.savefig(output_path, dpi=150, bbox_inches='tight')
    plt.close()
    print(f"Created: {os.path.basename(output_path)}")

def plot_frequency_bar(data, output_path, title, top_n=30):
    """Plot top N motif frequencies as bar chart"""
    if not data:
        print(f"Warning: No data for {title}")
        return
    
    # Sort by frequency and take top N
    sorted_data = sorted(data.items(), key=lambda x: x[1], reverse=True)[:top_n]
    
    if not sorted_data:
        return
    
    labels, values = zip(*sorted_data)
    
    fig, ax = plt.subplots(figsize=(14, 6))
    x = np.arange(len(labels))
    bars = ax.bar(x, values, color='#3498db', edgecolor='black', linewidth=0.5)
    
    # Add value labels on bars
    max_val = max(values) if values else 1
    for bar, val in zip(bars, values):
        height = bar.get_height()
        ax.text(bar.get_x() + bar.get_width()/2, height,
                str(int(val)), ha='center', va='bottom', fontsize=8, fontweight='bold')
    
    ax.set_xlabel('Repeat Motif', fontsize=12, fontweight='bold')
    ax.set_ylabel('Frequency', fontsize=12, fontweight='bold')
    ax.set_title(title, fontsize=14, fontweight='bold')
    ax.set_xticks(x)
    ax.set_xticklabels(labels, rotation=45, ha='right', fontsize=8)
    
    # Add margin at top for labels
    ax.set_ylim(0, max_val * 1.1)
    
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    
    plt.tight_layout()
    plt.savefig(output_path, dpi=150, bbox_inches='tight')
    plt.close()
    print(f"Created: {os.path.basename(output_path)}")

def main():
    results_dir = sys.argv[1]
    project_name = sys.argv[2]
    
    print(f"Generating Type 1 visualizations for {project_name}")
    
    # 1. Distribution of SSR classes
    dist_file = os.path.join(results_dir, "The distribution of the different SSR classes.csv")
    if os.path.exists(dist_file):
        data = parse_distribution_csv(dist_file)
        if data:
            plot_distribution(data, os.path.join(results_dir, "Distribution of the different SSR classes.png"))
    else:
        print(f"Warning: Missing {os.path.basename(dist_file)}")
    
    # 2. Frequency of identified SSR motifs
    freq_file = os.path.join(results_dir, "Frequency of the identified SSR motifs.csv")
    if os.path.exists(freq_file):
        data = parse_frequency_csv(freq_file)
        if data:
            plot_frequency_bar(data, os.path.join(results_dir, "Frequency of the identified SSR motifs.png"),
                              "Frequency of the Identified SSR Motifs")
    else:
        print(f"Warning: Missing {os.path.basename(freq_file)}")
    
    # 3. SSR distribution considering complementarity
    comp_file = os.path.join(results_dir, "Frequency of identified SSR motifs considering complementarity.csv")
    if os.path.exists(comp_file):
        data = parse_frequency_csv(comp_file)
        if data:
            plot_frequency_bar(data, os.path.join(results_dir, "SSR distribution considering sequence complementarity.png"),
                              "SSR Distribution Considering Sequence Complementarity")
    else:
        print(f"Warning: Missing {os.path.basename(comp_file)}")
    
    print("Type 1 visualization generation complete!")

if __name__ == "__main__":
    main()
VISUALIZATION_TYPE1_SCRIPT

    if [[ $? -eq 0 ]]; then
        log INFO "Type 1 visualization generation complete"
    else
        log WARN "Some visualizations may have failed - check output files"
    fi
}

#=============================================================================
# Main Execution
#=============================================================================
main() {
    parse_arguments "$@"
    initialize_pipeline
    
    # Phase 2: SSR Identification
    if should_skip_phase 2; then
        log INFO "Skipping Phase 2 (already completed)"
    else
        run_ssr_identification
        save_checkpoint 2
    fi
    
    # Analysis Type 1: SSR identification + primer design for ALL SSRs (treated as intergenic)
    # This matches the old MegaSSR behavior where -A 1 designed primers without needing a GFF
    if [[ "$ANALYSIS_TYPE" == "1" ]]; then
        # Skip Phase 3 (Gene Annotation) and Phase 4 (Genic Primer Design)
        # All SSRs are treated as intergenic in Phase 2
        
        # Phase 5: Intergenic Primer Design (for ALL SSRs)
        if should_skip_phase 5; then
            log INFO "Skipping Phase 5 (already completed)"
        else
            run_intergenic_primer_design
            save_checkpoint 5
        fi
        
        # Phase 6: Post-processing
        if should_skip_phase 6; then
            log INFO "Skipping Phase 6 (already completed)"
        else
            run_post_processing_type1
            save_checkpoint 6
        fi
        
        # Phase 7: Statistics (simplified for Type 1)
        if should_skip_phase 7; then
            log INFO "Skipping Phase 7 (already completed)"
        else
            generate_statistics_type1
            save_checkpoint 7
        fi
        
        # Phase 8: Visualizations (3 plots for Type 1)
        if should_skip_phase 8; then
            log INFO "Skipping Phase 8 (already completed)"
        else
            generate_visualizations_type1
            save_checkpoint 8
        fi
    fi
    
    # Analysis Type 2: Full pipeline with gene annotation
    if [[ "$ANALYSIS_TYPE" == "2" ]]; then
        # Phase 3: Gene Annotation
        if should_skip_phase 3; then
            log INFO "Skipping Phase 3 (already completed)"
        else
            run_gene_annotation
            save_checkpoint 3
        fi
        
        # Phase 4: Genic Primer Design
        if should_skip_phase 4; then
            log INFO "Skipping Phase 4 (already completed)"
        else
            run_genic_primer_design
            save_checkpoint 4
        fi
        
        # Phase 5: Intergenic Primer Design
        if should_skip_phase 5; then
            log INFO "Skipping Phase 5 (already completed)"
        else
            run_intergenic_primer_design
            save_checkpoint 5
        fi
        
        # Phase 6: Post-processing
        if should_skip_phase 6; then
            log INFO "Skipping Phase 6 (already completed)"
        else
            run_post_processing
            save_checkpoint 6
        fi
        
        # Phase 7: Advanced Statistics
        if should_skip_phase 7; then
            log INFO "Skipping Phase 7 (already completed)"
        else
            generate_advanced_statistics
            save_checkpoint 7
        fi
        
        # Phase 8: Visualizations
        if should_skip_phase 8; then
            log INFO "Skipping Phase 8 (already completed)"
        else
            generate_visualizations
            save_checkpoint 8
        fi
    fi
    
    log INFO "=========================================="
    log INFO "Pipeline Complete"
    log INFO "Results: $OUTPUT_DIR"
    log INFO "=========================================="
}

# Run main function
main "$@"