minimap2_align.sh

#!/bin/bash
set -euo pipefail
export LC_ALL=C    # deterministic ASCII sort order across locales

# ═══════════════════════════════════════════════════════════════
# Usage:
#   minimap2_align.sh [options] <bait_assembly.gfa> <unclassified_reads.fastq.gz> [output_prefix] [threads] [min_mapq] [min_read_length]
#
# Arguments:
#   bait_assembly    Path to Flye GFA assembly graph (bait)
#   unclassified     Path to unclassified reads FASTQ (gzipped)
#   output_prefix    Prefix for output files (default: linear_alignments)
#   threads          Number of CPU threads (default: 128)
#   min_mapq         Minimum mapping quality for filtering (default: 10)
#   min_read_length  Minimum read length to keep (default: 200)
#
# Options:
#   -h, --help      Show this help message and exit
#
# Expected CWD: the rescue workspace directory created by prepare_workspace.sh
# ═══════════════════════════════════════════════════════════════

show_help() {
    cat << EOF
Usage: $(basename "$0") [options] <bait_assembly.gfa> <unclassified_reads.fastq.gz> [output_prefix] [threads] [min_mapq] [min_read_length]

Linear alignment of ONT unclassified reads against a linearised Flye assembly
using Minimap2, followed by MAPQ-based filtering and read ID extraction.

Reads are pre-filtered with seqkit seq to remove short reads and sanitise FASTQ
formatting before alignment (set min_read_length based on seqkit stats / NanoPlot).

Arguments:
  bait_assembly    Path to Flye GFA assembly graph (bait)
  unclassified     Path to unclassified reads FASTQ (gzipped)
  output_prefix    Prefix for output files (default: linear_alignments)
  threads          Number of CPU threads (default: 128)
  min_mapq         Minimum mapping quality for HQ filtering (default: 10)
  min_read_length  Minimum read length in bp (default: 200)

Options:
  -h, --help      Show this help message and exit

Output files:
  assembly_linear.fasta              Linearised contigs extracted from GFA
  <prefix>_sorted.bam                Sorted BAM alignment
  <prefix>_sorted.bam.csi            BAM index (CSI format, written by --write-index)
  <prefix>_flagstat.txt              samtools flagstat report
  <prefix>_all_ids.txt               All mapped read IDs (unfiltered)
  <prefix>_filtered_ids.txt          MAPQ-filtered mapped read IDs
  minimap2.log                       Minimap2 stderr log

Notes:
  - Run from inside the rescue workspace (created by prepare_workspace.sh)
  - Requires: minimap2, samtools, seqkit, awk
  - Reads are pre-filtered: seqkit seq -m <min_read_length> -g (length + whitespace)
  - The GFA is first linearised to FASTA (S-lines extracted)
  - Minimap2 uses the map-ont preset with --secondary=no

Example:
  $(basename "$0") assembly_graph.gfa unclassified_reads.fastq.gz linear_alignments 128 10 200
EOF
}

if [[ "${1:-}" == "-h" ]] || [[ "${1:-}" == "--help" ]]; then
    show_help
    exit 0
fi

if [[ $# -lt 2 ]]; then
    show_help
    exit 1
fi

# ───────────────────────────────────────────────────────────────
# Argument parsing
# ───────────────────────────────────────────────────────────────
bait_assembly="$1"
unclassified="$2"
output_prefix="${3:-linear_alignments}"
threads="${4:-128}"
min_mapq="${5:-10}"
min_read_length="${6:-200}"

sorted_bam="${output_prefix}_sorted.bam"
flagstat_file="${output_prefix}_flagstat.txt"
all_ids="${output_prefix}_all_ids.txt"
filtered_ids="${output_prefix}_filtered_ids.txt"

# ───────────────────────────────────────────────────────────────
# Logging function
# ───────────────────────────────────────────────────────────────
log() {
    echo "[$(date '+%Y-%m-%d %H:%M:%S')] $1"
}

# ═══════════════════════════════════════════════════════════════
# Validate Inputs
# ═══════════════════════════════════════════════════════════════
log "=== Validating Inputs ==="

if [[ ! -f "${bait_assembly}" ]]; then
    log "ERROR: Assembly graph not found: ${bait_assembly}"
    exit 1
fi

if [[ ! -f "${unclassified}" ]]; then
    log "ERROR: Unclassified reads not found: ${unclassified}"
    exit 1
fi

log "Bait assembly:    ${bait_assembly}"
log "Unclassified:     ${unclassified}"
log "Output prefix:    ${output_prefix}"
log "Threads:          ${threads}"
log "Min MAPQ:         ${min_mapq}"
log "Min read length:  ${min_read_length} bp"

# ═══════════════════════════════════════════════════════════════
# Dependency Check
# ═══════════════════════════════════════════════════════════════
for cmd in minimap2 samtools seqkit awk grep; do
    if ! command -v "${cmd}" &> /dev/null; then
        log "ERROR: ${cmd} could not be found. Please ensure it is installed and in your PATH."
        exit 1
    fi
done

# ═══════════════════════════════════════════════════════════════
# GFA to FASTA Conversion
# ═══════════════════════════════════════════════════════════════
log "=== Converting GFA to linearised FASTA ==="

awk '/^S/{print ">"$2"\n"$3}' "${bait_assembly}" > assembly_linear.fasta

# Validate conversion
contig_count=$(grep -c '^>' assembly_linear.fasta || true)
if [[ ${contig_count} -eq 0 ]]; then
    log "ERROR: No contigs extracted from GFA — check file format"
    exit 1
fi

log "Extracted ${contig_count} contigs:"
seqkit stats --quiet assembly_linear.fasta

# ═══════════════════════════════════════════════════════════════
# Minimap2 Alignment
# ═══════════════════════════════════════════════════════════════
log "=== Minimap2: Linear Mapping ==="
log "Pre-filter: seqkit seq -m ${min_read_length} -g (length ≥${min_read_length} bp + whitespace removal)"
log "Minimap2 parameters:"
log "  -ax map-ont      ONT read preset"
log "  --secondary=no   Suppress secondary alignments"
log "  -t ${threads}  Thread count"

time_start=$(date +%s)

seqkit seq -m "${min_read_length}" -g "${unclassified}" \
    | minimap2 \
        -ax map-ont \
        -t "${threads}" \
        --secondary=no \
        assembly_linear.fasta \
        - \
        2> minimap2.log \
    | samtools sort \
        -@ "${threads}" \
        -o "${sorted_bam}" \
        --write-index \
        -

time_end=$(date +%s)
elapsed=$(( time_end - time_start ))
log "Minimap2 + sort finished in ${elapsed}s"

# Validate output
if [[ ! -s "${sorted_bam}" ]]; then
    log "ERROR: Minimap2 produced no output — check minimap2.log"
    exit 1
fi

ls -lh "${sorted_bam}"

# Flagstat report
log "=== Alignment Statistics ==="
samtools flagstat -@ "${threads}" "${sorted_bam}" | tee "${flagstat_file}"

# ═══════════════════════════════════════════════════════════════
# Extract Read IDs
# ═══════════════════════════════════════════════════════════════
log "=== Extracting Mapped Read IDs ==="

# all_ids: -F 4 only (exclude unmapped); supplementary reads harmlessly deduplicate
#           by read name via sort -u, so including them costs nothing and keeps us
#           maximally inclusive for circular-chromosome spanning reads.
# filtered_ids: -F 2052 (= -F 4 + -F 2048); minimap2 assigns MAPQ=0 to supplementary
#           alignments, so MAPQ filtering should be applied to primary alignments only.
# Note: -F 256 (secondary) is omitted — --secondary=no already suppresses them in minimap2.
samtools view -F 4 "${sorted_bam}" | cut -f1 | sort -u > "${all_ids}"
all_count=$(grep -c . "${all_ids}" || true)

samtools view -F 2052 -q "${min_mapq}" "${sorted_bam}" | cut -f1 | sort -u > "${filtered_ids}"
filtered_count=$(grep -c . "${filtered_ids}" || true)

log "All mapped reads (primary):       ${all_count}"
log "MAPQ-filtered reads (≥${min_mapq}):       ${filtered_count}"

# ═══════════════════════════════════════════════════════════════
# Metrics
# ═══════════════════════════════════════════════════════════════
if [[ -f "rescue_metrics.csv" ]]; then
    echo "$(date -Iseconds),minimap2,runtime_seconds,${elapsed}"         >> rescue_metrics.csv
    echo "$(date -Iseconds),minimap2,all_mapped_reads,${all_count}"     >> rescue_metrics.csv
    echo "$(date -Iseconds),minimap2,hq_mapped_reads,${filtered_count}" >> rescue_metrics.csv
    log "Metrics appended to rescue_metrics.csv"
else
    log "NOTE: rescue_metrics.csv not found — skipping metrics append"
fi

# ═══════════════════════════════════════════════════════════════
# Summary
# ═══════════════════════════════════════════════════════════════
log "────────────────────────────────────────────────────────────"
log "Minimap2 alignment complete."
log "  Linearised FASTA: assembly_linear.fasta  (${contig_count} contigs)"
log "  Sorted BAM:       ${sorted_bam}"
log "  Flagstat:          ${flagstat_file}"
log "  All read IDs:      ${all_ids}  (${all_count} reads)"
log "  HQ read IDs:       ${filtered_ids}  (${filtered_count} reads)"
log "  Aligner log:       minimap2.log"
log ""
log "Next step: compare with GraphAligner results and extract rescued reads"
log "────────────────────────────────────────────────────────────"