rustalign_cli/

align.rs

//! ferroux-align - Main alignment tool
//!
//! This is the main RustAlign alignment tool that reads sequence reads,
//! aligns them to a reference index using the FM-index, and outputs
//! SAM format alignments.

use anyhow::Result;
use clap::Parser;
use rustalign_aligner::{AlignerParams, AlignerWithRef, PeConstraints, PePolicy};
use rustalign_common::{Nuc, Strand};
use rustalign_io::{FastqReader, SamConfig, SamOpt, SamRecord, SamWriter};

/// Ferroux aligner - Fast and sensitive read alignment
///
/// Ferroux is an ultrafast and memory-efficient tool for aligning
/// sequencing reads to long reference sequences.
#[derive(Debug, Parser)]
#[command(name = "ferroux-align")]
#[command(author = "Ben Langmead <langmea@cs.jhu.edu>")]
#[command(version)]
#[command(about = "Ultrafast and memory-efficient read aligner", long_about = None)]
pub struct AlignOptions {
    // ========== Required Arguments ==========
    /// Index filename prefix (minus trailing .X.rai)
    #[clap(short = 'x', long)]
    pub index: String,

    /// Files with unpaired reads
    #[clap(short = 'U', long, value_name = "r")]
    pub reads_unpaired: Vec<String>,

    /// Files with #1 mates, paired with files in <m2>
    #[clap(short = '1', long, value_name = "m1")]
    pub reads_1: Vec<String>,

    /// Files with #2 mates, paired with files in <m1>
    #[clap(short = '2', long, value_name = "m2")]
    pub reads_2: Vec<String>,

    /// Output SAM file (default: stdout)
    #[clap(short = 'S', long, value_name = "sam")]
    pub sam_file: Option<String>,

    // ========== Input Options ==========
    /// Query input files are FASTQ .fq/.fastq (default)
    #[clap(short = 'q', long, conflicts_with = "fasta")]
    pub fastq: bool,

    /// Query input files are (multi-)FASTA .fa/.mfa
    #[clap(short = 'f', long, conflicts_with = "fastq")]
    pub fasta: bool,

    /// <m1>, <m2>, <r> are sequences themselves, not files
    #[clap(short = 'c', long)]
    pub sequences: bool,

    /// Skip the first <int> reads/pairs in the input
    #[clap(short = 's', long, value_name = "int")]
    pub skip: Option<usize>,

    /// Stop after first <int> reads/pairs
    #[clap(short = 'u', long = "upto", value_name = "int")]
    pub upto: Option<usize>,

    /// Trim <int> bases from 5'/left end of reads
    #[clap(short = '5', long = "trim5", value_name = "int")]
    pub trim5: Option<usize>,

    /// Trim <int> bases from 3'/right end of reads
    #[clap(short = '3', long = "trim3", value_name = "int")]
    pub trim3: Option<usize>,

    // ========== Alignment Options ==========
    /// Max # mismatches in seed alignment; can be 0 or 1
    #[clap(short = 'N', long, value_name = "int")]
    pub seed_mismatches: Option<u8>,

    /// Length of seed substrings; must be >3, <32
    #[clap(short = 'L', long, value_name = "int", default_value = "22")]
    pub seed_length: usize,

    /// Entire read must align; no clipping
    #[clap(long, conflicts_with = "local")]
    pub end_to_end: bool,

    /// Local alignment; ends might be soft clipped
    #[clap(long, conflicts_with = "end_to_end")]
    pub local: bool,

    /// Do not align forward (original) version of read
    #[clap(long = "nofw")]
    pub no_forward: bool,

    /// Do not align reverse-complement version of read
    #[clap(long = "norc")]
    pub no_reverse: bool,

    // ========== Scoring Options ==========
    /// Match bonus (0 for --end-to-end, 2 for --local)
    #[clap(long = "ma", value_name = "int")]
    pub match_bonus: Option<i32>,

    /// Max penalty for mismatch
    #[clap(long = "mp", value_name = "int", default_value = "6")]
    pub mismatch_penalty: i32,

    /// Penalty for non-A/C/G/Ts in read/ref
    #[clap(long = "np", value_name = "int", default_value = "1")]
    pub n_penalty: i32,

    /// Read gap open, extend penalties
    #[clap(long = "rdg", value_name = "int,int")]
    pub read_gap: Option<String>,

    /// Reference gap open, extend penalties
    #[clap(long = "rfg", value_name = "int,int")]
    pub ref_gap: Option<String>,

    // ========== Reporting Options ==========
    /// Report up to <int> alignments per read
    #[clap(short = 'k', long, value_name = "int")]
    pub report_alignments: Option<usize>,

    /// Report all alignments
    #[clap(short = 'a', long = "all")]
    pub report_all: bool,

    // ========== Effort Options ==========
    /// Give up extending after <int> failed extends
    #[clap(short = 'D', long, value_name = "int", default_value = "15")]
    pub failed_extends: usize,

    /// For reads with repetitive seeds, try <int> sets of seeds
    #[clap(short = 'R', long, value_name = "int", default_value = "2")]
    pub seed_sets: usize,

    // ========== Paired-end Options ==========
    /// Minimum fragment length
    #[clap(short = 'I', long = "minins", value_name = "int", default_value = "0")]
    pub min_insert: usize,

    /// Maximum fragment length
    #[clap(
        short = 'X',
        long = "maxins",
        value_name = "int",
        default_value = "500"
    )]
    pub max_insert: usize,

    /// -1, -2 mates align fw/rev, rev/fw, fw/fw
    #[clap(long)]
    pub orientation: Option<String>,

    // ========== Performance Options ==========
    /// Number of alignment threads
    #[clap(short = 'p', long, value_name = "int", default_value = "1")]
    pub threads: usize,

    /// Force SAM output order to match input
    #[clap(long = "reorder")]
    pub reorder: bool,

    // ========== Output Options ==========
    /// Print wall-clock time taken by search phases
    #[clap(short = 't', long = "time")]
    pub print_time: bool,

    /// Print nothing to stderr except serious errors
    #[clap(long)]
    pub quiet: bool,

    /// Suppress SAM records for unaligned reads
    #[clap(long = "no-unal")]
    pub no_unaligned: bool,

    /// Suppress header lines (@PG, @HD, etc.)
    #[clap(long = "no-head")]
    pub no_head: bool,

    /// Suppress @SQ header lines
    #[clap(long = "no-sq")]
    pub no_sq: bool,

    /// Use '='/'X', instead of 'M,' to specify matches/mismatches
    #[clap(long = "xeq")]
    pub xeq: bool,
}

/// Run the ferroux-align tool
pub fn run() -> Result<()> {
    let opts = AlignOptions::parse();
    run_with_opts(opts)
}

/// Run with explicit options (for subcommand support)
pub fn run_with_options(opts: AlignOptions) -> Result<()> {
    run_with_opts(opts)
}

#[allow(dead_code)]
fn main() -> Result<()> {
    run()
}

fn run_with_opts(opts: AlignOptions) -> Result<()> {
    // Load index to get reference names and lengths for SAM header
    let index_with_ref = rustalign_fmindex::EbwtWithRef::load(&opts.index)?;
    let ebwt = index_with_ref.ebwt();
    let ref_names = ebwt.refnames().to_vec();
    let ref_lengths = ebwt.plen().to_vec();

    // Use file output if -S is specified, otherwise stdout
    if let Some(ref sam_file) = opts.sam_file {
        let file = std::fs::File::create(sam_file)?;
        align_with_writer(&opts, file, &ref_names, &ref_lengths, index_with_ref)?;
    } else {
        let stdout = std::io::stdout();
        align_with_writer(&opts, stdout, &ref_names, &ref_lengths, index_with_ref)?;
    }

    Ok(())
}

fn align_with_writer<W: std::io::Write>(
    opts: &AlignOptions,
    writer: W,
    ref_names: &[String],
    ref_lengths: &[u32],
    index_with_ref: rustalign_fmindex::EbwtWithRef,
) -> Result<()> {
    // Set up SAM output
    let config = SamConfig::default();
    let mut writer = SamWriter::new(writer, config);
    writer.write_header(ref_names, ref_lengths)?;

    // Align reads
    if !opts.reads_unpaired.is_empty() {
        align_unpaired(opts, &mut writer, index_with_ref)?;
    } else if !opts.reads_1.is_empty() && !opts.reads_2.is_empty() {
        align_paired(opts, &mut writer, index_with_ref)?;
    }

    Ok(())
}

fn align_unpaired<W: std::io::Write>(
    opts: &AlignOptions,
    writer: &mut SamWriter<W>,
    index_with_ref: rustalign_fmindex::EbwtWithRef,
) -> Result<()> {
    let reads_file = &opts.reads_unpaired[0];
    let reader = FastqReader::from_path(reads_file)?;

    // Create aligner with reference
    let aligner_params = rustalign_aligner::AlignerParams::default();
    let aligner = rustalign_aligner::AlignerWithRef::new(index_with_ref, aligner_params);

    // Get a reference to the index for position resolution
    let index_for_resolution = aligner.index();

    for record in reader {
        let record = record?;

        // Sequence is already parsed as Vec<Nuc>
        let nucs = record.seq.clone();

        // Create reverse complement
        let nucs_rc: Vec<rustalign_common::Nuc> =
            nucs.iter().rev().map(|&n| n.complement()).collect();

        // Align the read with actual reference extraction
        let result = aligner.align(&nucs, &nucs_rc);

        // Create SAM record
        let mut sam_rec = rustalign_io::SamRecord::new(record.id.clone());
        sam_rec.seq = nucs.clone();
        sam_rec.qual = record.qual;

        // Check if we got successful alignments
        if !result.alignments.is_empty() {
            let read_len = nucs.len() as u64;
            let index_is_forward = true;

            // Resolve positions for ALL candidate alignments
            // This allows us to use chromosome as a tie-breaker
            let mut resolved: Vec<(usize, String, u64, Strand, String)> = Vec::new();

            for (idx, aln) in result.alignments.iter().enumerate() {
                match index_for_resolution.ebwt().resolve_position(
                    aln.bwt_row,
                    read_len,
                    index_is_forward,
                ) {
                    Ok((chr_name, pos, _chr_len)) => {
                        // Adjust position by subtracting the seed offset
                        let adjusted_pos = if pos > aln.seed_offset as u64 {
                            pos - aln.seed_offset as u64
                        } else {
                            1
                        };

                        resolved.push((idx, chr_name, adjusted_pos, aln.strand, aln.cigar.clone()));
                    }
                    Err(_) => continue,
                }
            }

            // Select the best alignment using tie-breakers
            // Order: score > organelle > hit_count > chromosome_order > position > strand > seed_offset
            fn chromosome_order(chr: &str) -> u32 {
                match chr {
                    "Mt" => 1, // Mitochondria first among organelles
                    "Pt" => 2, // Plastid second among organelles
                    "1" => 3,
                    "2" => 4,
                    "3" => 5,
                    "4" => 6,
                    "5" => 7,
                    _ => 99, // Unknown chromosomes come last
                }
            }

            // For same chromosome, use position (lower wins)
            // For different chromosomes with same organelle status (both nuclear or both organelle),
            // use hit_count to prefer more unique alignments
            fn is_organelle(chr: &str) -> bool {
                chr == "Mt" || chr == "Pt"
            }

            // Find the best alignment using:
            // 1. Score (higher is better)
            // 2. Organelle preference (Mt/Pt > nuclear for NUMT/NUPT cases)
            //    For reads matching NUMT/NUPT regions, the organelle genome (Mt/Pt) is
            //    the "true" source, while nuclear copies are derived.
            // 3. hit_count (lower is better - more unique, for cross-chromosome comparison)
            // 4. Chromosome order (lower chr number wins)
            // 5. Position (lower wins for same chromosome - deterministic)
            // 6. Strand: forward preferred (bowtie2 behavior)
            // 7. seed_offset (only used as final tie-breaker)
            if let Some((_idx, chr_name, adjusted_pos, strand, cigar)) =
                resolved.iter().min_by(|a, b| {
                    // Get original alignment for comparison
                    let aln_a = &result.alignments[a.0];
                    let aln_b = &result.alignments[b.0];

                    // Primary: score (higher is better)
                    match aln_b.score.cmp(&aln_a.score) {
                        std::cmp::Ordering::Equal => {
                            // Secondary: organelle preference (Mt/Pt > nuclear)
                            let a_organelle = is_organelle(&a.1);
                            let b_organelle = is_organelle(&b.1);
                            match (a_organelle, b_organelle) {
                                (true, false) => std::cmp::Ordering::Less,
                                (false, true) => std::cmp::Ordering::Greater,
                                _ => {
                                    // Tertiary: hit_count (lower is better - more unique)
                                    // This helps differentiate between chromosomes with equally good matches
                                    match aln_a.seed_hit_count.cmp(&aln_b.seed_hit_count) {
                                        std::cmp::Ordering::Equal => {
                                            // Quaternary: chromosome order (lower chr number wins)
                                            match chromosome_order(&a.1)
                                                .cmp(&chromosome_order(&b.1))
                                            {
                                                std::cmp::Ordering::Equal => {
                                                    // Quinary: position (lower wins)
                                                    match a.2.cmp(&b.2) {
                                                        std::cmp::Ordering::Equal => {
                                                            // Senary: forward strand preferred
                                                            match (aln_a.strand, aln_b.strand) {
                                                                (
                                                                    Strand::Forward,
                                                                    Strand::Reverse,
                                                                ) => std::cmp::Ordering::Less,
                                                                (
                                                                    Strand::Reverse,
                                                                    Strand::Forward,
                                                                ) => std::cmp::Ordering::Greater,
                                                                _ => aln_a
                                                                    .seed_offset
                                                                    .cmp(&aln_b.seed_offset),
                                                            }
                                                        }
                                                        other => other,
                                                    }
                                                }
                                                other => other,
                                            }
                                        }
                                        other => other,
                                    }
                                }
                            }
                        }
                        other => other,
                    }
                })
            {
                sam_rec.rname = chr_name.clone();
                sam_rec.pos = *adjusted_pos as i32;
                sam_rec.flag = if *strand == Strand::Reverse { 16 } else { 0 };
                sam_rec.mapq = 255;
                sam_rec.cigar = if cigar.is_empty() {
                    format!("{}M", nucs.len())
                } else {
                    cigar.clone()
                };
            } else {
                // All position resolutions failed
                if !opts.quiet {
                    eprintln!(
                        "DEBUG: Position resolution failed for all alignments of {}",
                        record.id
                    );
                }
            }
        } else {
            // Debug: output why no alignment
            if !opts.quiet {
                if let Some(err) = result.error {
                    eprintln!("DEBUG: No alignment for {}: {}", record.id, err);
                } else {
                    eprintln!("DEBUG: No alignment for {} (unknown reason)", record.id);
                }
            }
        }
        // If no alignment, record stays as unmapped (flag=0, rname=*, etc.)

        // Skip writing unmapped reads if --no-unal is specified
        if sam_rec.rname == "*" && opts.no_unaligned {
            continue;
        }

        writer.write(&sam_rec)?;
    }

    Ok(())
}

/// Resolved alignment with position information
#[derive(Clone)]
#[allow(dead_code)]
struct ResolvedAlignment {
    /// Index into the original alignment vector
    idx: usize,
    /// Chromosome name
    chr_name: String,
    /// Chromosome index (for paired-end comparison)
    chr_idx: u32,
    /// Adjusted position (1-based)
    pos: u64,
    /// Strand
    strand: Strand,
    /// CIGAR string
    cigar: String,
    /// Original alignment
    aln: rustalign_aligner::Alignment,
}

fn align_paired<W: std::io::Write>(
    opts: &AlignOptions,
    writer: &mut SamWriter<W>,
    index_with_ref: rustalign_fmindex::EbwtWithRef,
) -> Result<()> {
    let r1_file = &opts.reads_1[0];
    let r2_file = &opts.reads_2[0];

    let reader1 = FastqReader::from_path(r1_file)?;
    let reader2 = FastqReader::from_path(r2_file)?;

    // Create aligner
    let aligner_params = AlignerParams::default();
    let aligner = AlignerWithRef::new(index_with_ref, aligner_params);

    // Get reference to index for position resolution
    let index_for_resolution = aligner.index();

    // Create paired-end constraints from options
    let pe_policy = match opts.orientation.as_deref() {
        Some("fr") | None => PePolicy::FR, // Default is FR
        Some("rf") => PePolicy::RF,
        Some("ff") => PePolicy::FF,
        Some("rr") => PePolicy::RR,
        _ => {
            return Err(anyhow::anyhow!(
                "Invalid orientation: {}. Use fr, rf, ff, or rr",
                opts.orientation.as_ref().unwrap()
            ));
        }
    };

    let pe_constraints =
        PeConstraints::new(opts.min_insert as u32, opts.max_insert as u32, pe_policy);

    // Iterate through paired reads
    for (record1, record2) in reader1.zip(reader2) {
        let record1 = record1?;
        let record2 = record2?;

        // Parse sequences
        let nucs1 = record1.seq.clone();
        let nucs2 = record2.seq.clone();

        // Create reverse complements
        let nucs1_rc: Vec<Nuc> = nucs1.iter().rev().map(|&n| n.complement()).collect();
        let nucs2_rc: Vec<Nuc> = nucs2.iter().rev().map(|&n| n.complement()).collect();

        // Align both mates
        let result1 = aligner.align(&nucs1, &nucs1_rc);
        let result2 = aligner.align(&nucs2, &nucs2_rc);

        // Resolve positions for R1 alignments
        let mut resolved1: Vec<ResolvedAlignment> = Vec::new();
        if !result1.alignments.is_empty() {
            for (idx, aln) in result1.alignments.iter().enumerate() {
                if let Ok((chr_name, pos, _chr_len)) = index_for_resolution.ebwt().resolve_position(
                    aln.bwt_row,
                    nucs1.len() as u64,
                    true,
                ) {
                    let adjusted_pos = if pos > aln.seed_offset as u64 {
                        pos - aln.seed_offset as u64
                    } else {
                        1
                    };

                    // Get chromosome index
                    let chr_idx = index_for_resolution
                        .ebwt()
                        .refnames()
                        .iter()
                        .position(|n| n == &chr_name)
                        .unwrap_or(0) as u32;

                    resolved1.push(ResolvedAlignment {
                        idx,
                        chr_name,
                        chr_idx,
                        pos: adjusted_pos,
                        strand: aln.strand,
                        cigar: aln.cigar.clone(),
                        aln: aln.clone(),
                    });
                }
            }
        }

        // Resolve positions for R2 alignments
        let mut resolved2: Vec<ResolvedAlignment> = Vec::new();
        if !result2.alignments.is_empty() {
            for (idx, aln) in result2.alignments.iter().enumerate() {
                if let Ok((chr_name, pos, _chr_len)) = index_for_resolution.ebwt().resolve_position(
                    aln.bwt_row,
                    nucs2.len() as u64,
                    true,
                ) {
                    let adjusted_pos = if pos > aln.seed_offset as u64 {
                        pos - aln.seed_offset as u64
                    } else {
                        1
                    };

                    let chr_idx = index_for_resolution
                        .ebwt()
                        .refnames()
                        .iter()
                        .position(|n| n == &chr_name)
                        .unwrap_or(0) as u32;

                    resolved2.push(ResolvedAlignment {
                        idx,
                        chr_name,
                        chr_idx,
                        pos: adjusted_pos,
                        strand: aln.strand,
                        cigar: aln.cigar.clone(),
                        aln: aln.clone(),
                    });
                }
            }
        }

        // Find best concordant pair
        let best_pair = find_best_pair(
            &resolved1,
            &resolved2,
            &pe_constraints,
            nucs1.len() as u32,
            nucs2.len() as u32,
        );

        // Write SAM records
        match best_pair {
            Some((r1, r2, insert_size, is_concordant)) => {
                // Write R1
                let mut sam1 = SamRecord::new(record1.id.clone());
                sam1.seq = nucs1.clone();
                sam1.qual = record1.qual.clone();
                sam1.rname = r1.chr_name.clone();
                sam1.pos = r1.pos as i32;
                sam1.flag = calculate_flag(r1.strand, r2.strand, true, true, is_concordant);
                sam1.mapq = 255;
                sam1.cigar = if r1.cigar.is_empty() {
                    format!("{}M", nucs1.len())
                } else {
                    r1.cigar.clone()
                };
                sam1.rnext = if r1.chr_name == r2.chr_name {
                    "=".to_string()
                } else {
                    r2.chr_name.clone()
                };
                sam1.pnext = r2.pos as i32;
                sam1.tlen = insert_size;
                sam1.add_opt(SamOpt::Int("NM".to_string(), r1.aln.edits as i64));
                writer.write(&sam1)?;

                // Write R2
                let mut sam2 = SamRecord::new(record2.id.clone());
                sam2.seq = nucs2.clone();
                sam2.qual = record2.qual.clone();
                sam2.rname = r2.chr_name.clone();
                sam2.pos = r2.pos as i32;
                sam2.flag = calculate_flag(r2.strand, r1.strand, false, true, is_concordant);
                sam2.mapq = 255;
                sam2.cigar = if r2.cigar.is_empty() {
                    format!("{}M", nucs2.len())
                } else {
                    r2.cigar.clone()
                };
                sam2.rnext = if r1.chr_name == r2.chr_name {
                    "=".to_string()
                } else {
                    r1.chr_name.clone()
                };
                sam2.pnext = r1.pos as i32;
                sam2.tlen = -insert_size; // R2 has negative TLEN
                sam2.add_opt(SamOpt::Int("NM".to_string(), r2.aln.edits as i64));
                writer.write(&sam2)?;
            }
            None => {
                // No concordant pair found - write as unpaired
                if !opts.no_unaligned {
                    // Try to write single-end alignments if available
                    if let Some(r1) = select_best_alignment(&resolved1) {
                        let mut sam1 = SamRecord::new(record1.id.clone());
                        sam1.seq = nucs1;
                        sam1.qual = record1.qual;
                        sam1.rname = r1.chr_name;
                        sam1.pos = r1.pos as i32;
                        sam1.flag = calculate_flag(r1.strand, Strand::Forward, true, false, false);
                        sam1.mapq = 255;
                        sam1.cigar = format!("{}M", sam1.seq.len());
                        sam1.rnext = "*".to_string();
                        sam1.pnext = 0;
                        sam1.tlen = 0;
                        writer.write(&sam1)?;
                    } else {
                        // R1 unmapped
                        let mut sam1 = SamRecord::new(record1.id.clone());
                        sam1.seq = nucs1;
                        sam1.qual = record1.qual;
                        sam1.flag = 0x4 | 0x40 | 0x8; // UNMAP | FIRST | MUNMAP
                        writer.write(&sam1)?;
                    }

                    if let Some(r2) = select_best_alignment(&resolved2) {
                        let mut sam2 = SamRecord::new(record2.id.clone());
                        sam2.seq = nucs2;
                        sam2.qual = record2.qual;
                        sam2.rname = r2.chr_name;
                        sam2.pos = r2.pos as i32;
                        sam2.flag = calculate_flag(r2.strand, Strand::Forward, false, false, false);
                        sam2.mapq = 255;
                        sam2.cigar = format!("{}M", sam2.seq.len());
                        sam2.rnext = "*".to_string();
                        sam2.pnext = 0;
                        sam2.tlen = 0;
                        writer.write(&sam2)?;
                    } else {
                        // R2 unmapped
                        let mut sam2 = SamRecord::new(record2.id.clone());
                        sam2.seq = nucs2;
                        sam2.qual = record2.qual;
                        sam2.flag = 0x4 | 0x80 | 0x8; // UNMAP | SECOND | MUNMAP
                        writer.write(&sam2)?;
                    }
                }
            }
        }
    }

    Ok(())
}

/// Find the best concordant pair from R1 and R2 alignments
fn find_best_pair(
    r1_alignments: &[ResolvedAlignment],
    r2_alignments: &[ResolvedAlignment],
    constraints: &PeConstraints,
    r1_len: u32,
    r2_len: u32,
) -> Option<(ResolvedAlignment, ResolvedAlignment, i32, bool)> {
    let mut best_pair: Option<(ResolvedAlignment, ResolvedAlignment, i32, bool, i32)> = None;

    for r1 in r1_alignments {
        for r2 in r2_alignments {
            // Check if this pair is concordant
            let strand1 = r1.strand == Strand::Forward;
            let strand2 = r2.strand == Strand::Forward;

            if let Some(insert_size) = constraints.is_concordant(
                r1.pos, r1_len, strand1, r1.chr_idx, r2.pos, r2_len, strand2, r2.chr_idx,
            ) {
                // Calculate combined score
                let combined_score = r1.aln.score + r2.aln.score;

                // Keep the best concordant pair
                let is_better = match &best_pair {
                    None => true,
                    Some((_, _, _, _, best_score)) => combined_score > *best_score,
                };

                if is_better {
                    best_pair = Some((r1.clone(), r2.clone(), insert_size, true, combined_score));
                }
            }
        }
    }

    // If no concordant pair, try to find a discordant pair (same chr but wrong orientation/insert)
    if best_pair.is_none() {
        for r1 in r1_alignments {
            for r2 in r2_alignments {
                if r1.chr_idx == r2.chr_idx {
                    let _strand1 = r1.strand == Strand::Forward;
                    let _strand2 = r2.strand == Strand::Forward;

                    // Calculate insert size anyway
                    let left_pos = r1.pos.min(r2.pos);
                    let right_pos = (r2.pos + r2_len as u64).max(r1.pos + r1_len as u64);
                    let insert_size = (right_pos - left_pos) as i32;

                    let combined_score = r1.aln.score + r2.aln.score;

                    let is_better = match &best_pair {
                        None => true,
                        Some((_, _, _, _, best_score)) => combined_score > *best_score,
                    };

                    if is_better {
                        best_pair =
                            Some((r1.clone(), r2.clone(), insert_size, false, combined_score));
                    }
                }
            }
        }
    }

    best_pair.map(|(r1, r2, insert, concordant, _)| (r1, r2, insert, concordant))
}

/// Select the best alignment from a list
fn select_best_alignment(alignments: &[ResolvedAlignment]) -> Option<ResolvedAlignment> {
    alignments
        .iter()
        .max_by(|a, b| {
            // First compare by score
            match a.aln.score.cmp(&b.aln.score) {
                std::cmp::Ordering::Equal => {
                    // Then by hit count (lower is better)
                    b.aln.seed_hit_count.cmp(&a.aln.seed_hit_count)
                }
                other => other,
            }
        })
        .cloned()
}

/// Calculate SAM flag for paired-end alignment
fn calculate_flag(
    strand: Strand,
    mate_strand: Strand,
    is_first: bool,
    both_mapped: bool,
    is_concordant: bool,
) -> u16 {
    let mut flag: u16 = 0;

    // PAIRED (0x1)
    flag |= 0x1;

    // PROPER_PAIR (0x2) - set if concordant
    if is_concordant {
        flag |= 0x2;
    }

    // UNMAP (0x4) - not set since we have alignment

    // MUNMAP (0x8) - mate unmapped
    if !both_mapped {
        flag |= 0x8;
    }

    // REVERSE (0x10)
    if strand == Strand::Reverse {
        flag |= 0x10;
    }

    // MREVERSE (0x20)
    if mate_strand == Strand::Reverse {
        flag |= 0x20;
    }

    // FIRST_IN_PAIR (0x40) or SECOND_IN_PAIR (0x80)
    if is_first {
        flag |= 0x40;
    } else {
        flag |= 0x80;
    }

    flag
}