read_align.rs

/// Read alignment driver function
use crate::align::score::{AlignmentScorer, SpliceMotif};
use crate::align::seed::Seed;
use crate::align::stitch::{cluster_seeds, stitch_seeds_with_jdb_debug};
use crate::align::transcript::Transcript;
use crate::error::Error;
use crate::index::GenomeIndex;
use crate::params::{IntronMotifFilter, IntronStrandFilter, Parameters};
use crate::stats::UnmappedReason;
use rand::{SeedableRng, rngs::StdRng, seq::SliceRandom};
use std::hash::{DefaultHasher, Hash, Hasher};

/// Derive a deterministic per-read RNG seed from `run_rng_seed` + the read name.
///
/// STAR seeds `std::mt19937` once per chunk/thread (`runRNGseed*(iChunk+1)`),
/// then advances the state sequentially per read. ruSTAR parallelises per-read
/// via rayon, so we instead fold the read name into the seed — this keeps tie
/// breaks reproducible regardless of thread count while still honoring the
/// user's `--runRNGseed` value.
fn per_read_seed(run_rng_seed: u64, read_name: &str) -> u64 {
    let mut hasher = DefaultHasher::new();
    read_name.hash(&mut hasher);
    run_rng_seed.wrapping_mul(hasher.finish().wrapping_add(1))
}

/// Shuffle the prefix of `items` whose `score_fn` equals the first element's score.
///
/// Mirrors STAR's `ReadAlign_multMapSelect` / `funPrimaryAlignMark`: best-scoring
/// alignments are randomized so primary selection (index 0) is not biased by
/// upstream sort order. Non-tied elements are left alone.
fn shuffle_tied_prefix<T>(items: &mut [T], score_fn: impl Fn(&T) -> i32, seed: u64) {
    let Some(first) = items.first() else {
        return;
    };
    let best = score_fn(first);
    let tied = items.iter().take_while(|t| score_fn(t) == best).count();
    if tied < 2 {
        return;
    }
    items[..tied].shuffle(&mut StdRng::seed_from_u64(seed));
}

/// Result of aligning a single read: (transcripts, chimeric_alignments, n_for_mapq, unmapped_reason)
pub type AlignReadResult = (
    Vec<Transcript>,
    Vec<crate::chimeric::ChimericAlignment>,
    usize,
    Option<UnmappedReason>,
);

/// Paired-end alignment result
#[derive(Debug, Clone)]
pub struct PairedAlignment {
    /// Transcript for mate1
    pub mate1_transcript: Transcript,
    /// Transcript for mate2
    pub mate2_transcript: Transcript,
    /// Read positions for mate1 in transcript (start, end)
    pub mate1_region: (usize, usize),
    /// Read positions for mate2 in transcript (start, end)
    pub mate2_region: (usize, usize),
    /// Whether this is a proper pair (same chr, concordant orientation, distance)
    pub is_proper_pair: bool,
    /// Signed insert size (TLEN) - genomic distance between mate starts
    pub insert_size: i32,
    /// Combined pair score: sum of per-mate finalized scores (each includes genomic length penalty).
    /// Used for multi-mapper score-range ranking and mappedFilter quality check.
    pub combined_wt_score: i32,
    /// Combined coverage: sum of exon read spans from both mates.
    /// Mirrors STAR's nMatch check in mappedFilter.
    pub combined_n_match: u32,
}

/// Result of paired-end alignment, covering all mapping outcomes.
#[derive(Debug, Clone)]
pub enum PairedAlignmentResult {
    /// Both mates mapped and paired successfully
    BothMapped(Box<PairedAlignment>),
    /// Only one mate mapped; rescue failed or was not attempted for the other
    HalfMapped {
        /// Transcript of the mapped mate
        mapped_transcript: Transcript,
        /// true = mate1 is the mapped mate, false = mate2 is mapped
        mate1_is_mapped: bool,
    },
}

/// Align a read to the genome.
///
/// # Algorithm
/// 1. Find seeds (exact matches) using MMP search
/// 2. Cluster seeds by genomic proximity
/// 3. Stitch seeds within each cluster using DP
/// 4. Filter transcripts by quality thresholds
/// 5. Sort by score and limit to top N
/// 6. Detect chimeric alignments if enabled
///
/// # Arguments
/// * `read_seq` - Read sequence (encoded as 0=A, 1=C, 2=G, 3=T)
/// * `read_name` - Read name (needed for chimeric output)
/// * `index` - Genome index
/// * `params` - User parameters
///
/// # Returns
/// Tuple of (transcripts, chimeric alignments, n_for_mapq, unmapped_reason):
/// - transcripts: sorted by score (best first)
/// - chimeric alignments: sorted by score (best first)
/// - n_for_mapq: effective alignment count for MAPQ calculation (max of transcript count
///   and valid cluster count, to avoid undercounting from coordinate dedup on tandem repeats)
/// - unmapped_reason: `Some(reason)` if no alignments produced, `None` if mapped
pub fn align_read(
    read_seq: &[u8],
    read_name: &str,
    index: &GenomeIndex,
    params: &Parameters,
) -> Result<AlignReadResult, Error> {
    let debug_read = !params.read_name_filter.is_empty() && read_name == params.read_name_filter;

    // Step 1: Find seeds (seedMapMin from params)
    let min_seed_length = params.seed_map_min;
    let seeds = Seed::find_seeds(
        read_seq,
        index,
        min_seed_length,
        params,
        if debug_read { read_name } else { "" },
    )?;

    if debug_read {
        let total_positions: usize = seeds.iter().map(|s| s.sa_end - s.sa_start).sum();
        eprintln!(
            "[DEBUG {}] Seeds: {} seeds, {} total SA positions, read_len={}",
            read_name,
            seeds.len(),
            total_positions,
            read_seq.len()
        );
        let n_lr = seeds.iter().filter(|s| !s.search_rc).count();
        let n_rl = seeds.iter().filter(|s| s.search_rc).count();
        eprintln!(
            "  {} seeds total: {} L→R (sparse), {} R→L (sparse)",
            seeds.len(),
            n_lr,
            n_rl
        );
        for (i, s) in seeds.iter().enumerate() {
            let n_loci = s.sa_end - s.sa_start;
            eprintln!(
                "  seed[{}]: read_pos={}, len={}, n_loci={}, search_rc={}, sa=[{},{})",
                i, s.read_pos, s.length, n_loci, s.search_rc, s.sa_start, s.sa_end
            );
        }
    }

    if seeds.is_empty() {
        if debug_read {
            eprintln!("[DEBUG {}] No seeds found — unmapped", read_name);
        }
        return Ok((Vec::new(), Vec::new(), 0, Some(UnmappedReason::Other)));
    }

    // Step 2: Cluster seeds (STAR's bin-based windowing)
    // seed_per_window_nmax capacity eviction is handled inside cluster_seeds()
    let clusters = cluster_seeds(&seeds, index, params, read_seq.len());

    if debug_read {
        eprintln!(
            "[DEBUG {}] Clusters: {} clusters",
            read_name,
            clusters.len()
        );
        for (i, cluster) in clusters.iter().enumerate() {
            let chr_name = if cluster.chr_idx < index.genome.chr_name.len() {
                &index.genome.chr_name[cluster.chr_idx]
            } else {
                "unknown"
            };
            eprintln!(
                "  cluster[{}]: chr={}, is_reverse={}, seeds={}, anchor_bin={}",
                i,
                chr_name,
                cluster.is_reverse,
                cluster.alignments.len(),
                cluster.anchor_bin,
            );
            for (j, wa) in cluster.alignments.iter().enumerate() {
                let chr_pos = wa.genome_pos.saturating_sub(
                    if cluster.chr_idx < index.genome.chr_start.len() {
                        index.genome.chr_start[cluster.chr_idx]
                    } else {
                        0
                    },
                ) + 1; // 1-based
                eprintln!(
                    "    wa[{}]: read_pos={}, len={}, genome_pos={} ({}:{}), n_rep={}, is_anchor={}",
                    j,
                    wa.read_pos,
                    wa.length,
                    wa.genome_pos,
                    chr_name,
                    chr_pos,
                    wa.n_rep,
                    wa.is_anchor
                );
                if j >= 5 {
                    eprintln!(
                        "    ... ({} more WA entries)",
                        cluster.alignments.len() - j - 1
                    );
                    break;
                }
            }
        }
    }

    if clusters.is_empty() {
        if debug_read {
            eprintln!("[DEBUG {}] No clusters — unmapped", read_name);
        }
        return Ok((Vec::new(), Vec::new(), 0, Some(UnmappedReason::Other)));
    }

    // Cap total clusters (alignWindowsPerReadNmax)
    let mut clusters = clusters;
    clusters.truncate(params.align_windows_per_read_nmax);

    // NOTE: STAR's winReadCoverageRelativeMin filter is long-reads-only
    // (#ifdef COMPILE_FOR_LONG_READS in stitchPieces.cpp). Standard STAR
    // does NOT filter clusters by seed coverage. Removed to match STAR.

    // Step 2b: Detect chimeric alignments from multi-cluster seeds (Tier 2)
    let mut chimeric_alignments = Vec::new();
    if params.chim_segment_min > 0 && clusters.len() > 1 {
        use crate::chimeric::ChimericDetector;
        let detector = ChimericDetector::new(params);
        chimeric_alignments
            .extend(detector.detect_from_multi_clusters(&clusters, read_seq, read_name, index)?);
    }

    // Step 3: Stitch seeds within each cluster
    let scorer = AlignmentScorer::from_params(params);
    let mut transcripts = Vec::new();

    // Use junction DB for annotation-aware scoring if available
    let junction_db = if index.junction_db.is_empty() {
        None
    } else {
        Some(&index.junction_db)
    };

    for (ci, cluster) in clusters.iter().enumerate() {
        let debug_name = if debug_read { read_name } else { "" };
        let cluster_transcripts = stitch_seeds_with_jdb_debug(
            cluster,
            read_seq,
            index,
            &scorer,
            junction_db,
            params.align_transcripts_per_window_nmax,
            debug_name,
        )?;
        if debug_read {
            eprintln!(
                "[DEBUG {}] Cluster[{}]: {} transcripts from DP",
                read_name,
                ci,
                cluster_transcripts.len()
            );
            for (ti, t) in cluster_transcripts.iter().enumerate().take(5) {
                let chr_name = if t.chr_idx < index.genome.chr_name.len() {
                    &index.genome.chr_name[t.chr_idx]
                } else {
                    "unknown"
                };
                let cigar_str: String = t.cigar.iter().map(|op| format!("{}", op)).collect();
                eprintln!(
                    "  transcript[{}]: chr={}:{}-{} ({}) score={} mm={} junctions={} cigar={}",
                    ti,
                    chr_name,
                    t.genome_start,
                    t.genome_end,
                    if t.is_reverse { "-" } else { "+" },
                    t.score,
                    t.n_mismatch,
                    t.n_junction,
                    cigar_str
                );
            }
        }
        transcripts.extend(cluster_transcripts);
    }

    // Step 4a: Deduplicate and score-range filter — BEFORE quality filters.
    // STAR order: multMapSelect (score-range) → mappedFilter (quality gates).
    // Doing quality filters first is wrong: it can remove the high-scoring primary,
    // leaving a lower-scoring secondary that then passes as the "best" alignment.

    // Deduplicate transcripts with identical genomic coordinates AND CIGAR.
    transcripts.sort_by(|a, b| {
        (
            a.chr_idx,
            a.genome_start,
            a.genome_end,
            a.is_reverse,
            &a.cigar,
        )
            .cmp(&(
                b.chr_idx,
                b.genome_start,
                b.genome_end,
                b.is_reverse,
                &b.cigar,
            ))
            .then_with(|| b.score.cmp(&a.score))
    });
    transcripts.dedup_by(|a, b| {
        a.chr_idx == b.chr_idx
            && a.genome_start == b.genome_start
            && a.genome_end == b.genome_end
            && a.is_reverse == b.is_reverse
            && a.cigar == b.cigar
    });

    // Sort by score descending (deterministic tie-breaking).
    transcripts.sort_by(|a, b| {
        b.score
            .cmp(&a.score)
            .then_with(|| a.n_junction.cmp(&b.n_junction))
            .then_with(|| a.chr_idx.cmp(&b.chr_idx))
            .then_with(|| a.genome_start.cmp(&b.genome_start))
            .then_with(|| a.is_reverse.cmp(&b.is_reverse))
    });

    // Randomize primary among best-scoring ties (ReadAlign_multMapSelect.cpp:71-79).
    shuffle_tied_prefix(
        &mut transcripts,
        |t| t.score,
        per_read_seed(params.run_rng_seed, read_name),
    );

    // Score-range filter: keep only alignments within outFilterMultimapScoreRange of the best.
    // (STAR's multMapSelect step — must run before quality filters.)
    if !transcripts.is_empty() {
        let max_score = transcripts[0].score;
        let score_threshold = max_score - params.out_filter_multimap_score_range;
        transcripts.retain(|t| t.score >= score_threshold);
    }

    // Multimap count check: too many loci → unmapped.
    if transcripts.len() > params.out_filter_multimap_nmax as usize {
        let n_loci = transcripts.len();
        transcripts.clear();
        return Ok((
            transcripts,
            chimeric_alignments,
            n_loci,
            Some(UnmappedReason::TooManyLoci),
        ));
    }

    // Step 4: Quality filters (STAR's mappedFilter — runs after score-range selection).
    // STAR uses (Lread-1) for relative thresholds and casts to integer
    // (ReadAlign_mappedFilter.cpp lines 8-9)
    let read_length = read_seq.len() as f64;
    let lread_m1 = (read_seq.len() as f64) - 1.0;

    // Log filtering statistics
    let pre_filter_count = transcripts.len();
    let mut filter_reasons = std::collections::HashMap::new();

    transcripts.retain(|t| {
        // Absolute score threshold
        if t.score < params.out_filter_score_min {
            *filter_reasons.entry("score_min").or_insert(0) += 1;
            return false;
        }

        // Relative score threshold: STAR casts to intScore (i32)
        if t.score < (params.out_filter_score_min_over_lread * lread_m1) as i32 {
            *filter_reasons.entry("score_min_relative").or_insert(0) += 1;
            return false;
        }

        // Absolute mismatch count
        if t.n_mismatch > params.out_filter_mismatch_nmax {
            *filter_reasons.entry("mismatch_max").or_insert(0) += 1;
            log::debug!(
                "Filtered {}: {} mismatches > {} max (read_len={}, score={})",
                read_name,
                t.n_mismatch,
                params.out_filter_mismatch_nmax,
                read_length,
                t.score
            );
            return false;
        }

        // Relative mismatch count (mismatches / read_length)
        let mismatch_rate = t.n_mismatch as f64 / read_length;
        if mismatch_rate > params.out_filter_mismatch_nover_lmax {
            *filter_reasons.entry("mismatch_rate").or_insert(0) += 1;
            log::debug!(
                "Filtered {}: {:.1}% mismatch rate > {:.1}% max ({}/{} bases, score={})",
                read_name,
                mismatch_rate * 100.0,
                params.out_filter_mismatch_nover_lmax * 100.0,
                t.n_mismatch,
                read_length,
                t.score
            );
            return false;
        }

        // Absolute matched bases
        let n_matched = t.n_matched();
        if n_matched < params.out_filter_match_nmin {
            *filter_reasons.entry("match_min").or_insert(0) += 1;
            return false;
        }

        // Relative matched bases: STAR casts to uint (u32)
        if n_matched < (params.out_filter_match_nmin_over_lread * lread_m1) as u32 {
            *filter_reasons.entry("match_min_relative").or_insert(0) += 1;
            return false;
        }

        // Junction motif filtering
        match params.out_filter_intron_motifs {
            IntronMotifFilter::None => {
                // Accept all motifs
            }
            IntronMotifFilter::RemoveNoncanonical => {
                // Reject if any junction is non-canonical
                if t.junction_motifs.contains(&SpliceMotif::NonCanonical) {
                    *filter_reasons.entry("noncanonical_junction").or_insert(0) += 1;
                    return false;
                }
            }
            IntronMotifFilter::RemoveNoncanonicalUnannotated => {
                // Only reject if a non-canonical junction is NOT annotated in GTF
                if t.junction_motifs
                    .iter()
                    .zip(t.junction_annotated.iter())
                    .any(|(m, annotated)| *m == SpliceMotif::NonCanonical && !annotated)
                {
                    *filter_reasons
                        .entry("noncanonical_unannotated_junction")
                        .or_insert(0) += 1;
                    return false;
                }
            }
        }

        // Intron strand consistency filtering (outFilterIntronStrands)
        // STAR's RemoveInconsistentStrands removes transcripts that have junctions
        // with mixed intron strand (some imply + strand, some imply - strand).
        // This handles chimeric/impossible transcripts spanning both strands.
        // Note: a reverse-strand read CAN have + strand motifs (antisense reads
        // from + strand genes in unstranded RNA-seq) — this is valid and STAR
        // keeps such reads. Only mixed-strand within one transcript is filtered.
        if params.out_filter_intron_strands == IntronStrandFilter::RemoveInconsistentStrands {
            let mut has_plus = false;
            let mut has_minus = false;
            for motif in &t.junction_motifs {
                match motif.implied_strand() {
                    Some('+') => has_plus = true,
                    Some('-') => has_minus = true,
                    None => {}
                    _ => {}
                }
            }
            if has_plus && has_minus {
                *filter_reasons.entry("inconsistent_strand").or_insert(0) += 1;
                return false;
            }
        }

        true
    });

    // Log filtering summary if anything was filtered
    if pre_filter_count > transcripts.len() {
        let filtered = pre_filter_count - transcripts.len();
        log::debug!(
            "Read {}: Filtered {}/{} transcripts: {:?}",
            read_name,
            filtered,
            pre_filter_count,
            filter_reasons
        );
    }

    if debug_read {
        eprintln!(
            "[DEBUG {}] After quality filters: {}/{} transcripts remain (reasons: {:?})",
            read_name,
            transcripts.len(),
            pre_filter_count,
            filter_reasons
        );
    }

    // Step 3b: Detect chimeric alignments from soft-clips (Tier 1)
    if params.chim_segment_min > 0 {
        use crate::chimeric::ChimericDetector;
        let detector = ChimericDetector::new(params);
        for transcript in &transcripts {
            if let Some(chim) =
                detector.detect_from_soft_clips(transcript, read_seq, read_name, index)?
            {
                chimeric_alignments.push(chim);
            }
        }
    }

    // Note: STAR sometimes finds 2 equivalent indel placements in homopolymer runs
    // via its recursive stitcher's seed exploration (NH=2 instead of NH=1 for ~5 reads).
    // Generating equivalents post-hoc causes more harm than good (41 false NH=2 vs 5 fixed).
    // The root cause is jR scanning placing insertions at different positions — fixing that
    // would be a better approach than post-hoc enumeration.

    // Step 6: Filter chimeric alignments
    if params.chim_segment_min > 0 {
        chimeric_alignments.retain(|chim| {
            chim.meets_min_segment_length(params.chim_segment_min)
                && chim.meets_min_score(params.chim_score_min)
        });
    }

    // n_for_mapq = transcripts.len() after dedup and filtering.
    // Multi-transcript DP (Phase 16.10) produces multiple transcripts per window
    // for tandem repeats (e.g. rDNA), yielding correct NH → correct MAPQ.
    let n_for_mapq = transcripts.len();

    if debug_read {
        eprintln!(
            "[DEBUG {}] Final: {} transcripts, n_for_mapq={}",
            read_name,
            transcripts.len(),
            n_for_mapq
        );
        for (i, t) in transcripts.iter().enumerate() {
            let chr_name = if t.chr_idx < index.genome.chr_name.len() {
                &index.genome.chr_name[t.chr_idx]
            } else {
                "unknown"
            };
            let cigar_str: String = t.cigar.iter().map(|op| format!("{}", op)).collect();
            eprintln!(
                "  FINAL[{}]: chr={}:{}-{} ({}) score={} mm={} junctions={} cigar={}",
                i,
                chr_name,
                t.genome_start,
                t.genome_end,
                if t.is_reverse { "-" } else { "+" },
                t.score,
                t.n_mismatch,
                t.n_junction,
                cigar_str
            );
        }
    }

    let unmapped_reason = if transcripts.is_empty() {
        // Transcripts were generated by DP but all filtered out
        Some(UnmappedReason::TooShort)
    } else {
        None
    };

    Ok((
        transcripts,
        chimeric_alignments,
        n_for_mapq,
        unmapped_reason,
    ))
}

/// Align paired-end reads using per-mate independent seeding and pairwise cluster matching.
///
/// # Algorithm
/// 1. Seed mate1_seq and mate2_seq independently (correct Nstart per mate)
/// 2. Cluster each mate's seeds independently
/// 3. Stitch each mate's clusters independently via stitch_seeds_with_jdb
/// 4. Pair compatible transcripts (same chr, opposite strand, within alignMatesGapMax)
/// 5. Score-range filter → TooManyLoci check → dedup → quality filter
/// 6. If no pairs: rescue as HalfMapped using single-mate transcripts
///
/// # Arguments
/// * `mate1_seq` - First mate sequence (encoded)
/// * `mate2_seq` - Second mate sequence (encoded)
/// * `index` - Genome index
/// * `params` - Parameters (includes alignMatesGapMax)
///
/// # Returns
/// Tuple of (paired alignment results, n_for_mapq, unmapped_reason)
pub fn align_paired_read(
    mate1_seq: &[u8],
    mate2_seq: &[u8],
    read_name: &str,
    index: &GenomeIndex,
    params: &Parameters,
) -> Result<(Vec<PairedAlignmentResult>, usize, Option<UnmappedReason>), Error> {
    let len1 = mate1_seq.len();
    let len2 = mate2_seq.len();

    let debug_pe = !params.read_name_filter.is_empty() && read_name == params.read_name_filter;

    let scorer = AlignmentScorer::from_params(params);
    let junction_db = if index.junction_db.is_empty() {
        None
    } else {
        Some(&index.junction_db)
    };

    // Per-mate seeding: each mate seeded with its own Nstart positions.
    // mate2_seq (forward orientation) gives Nstart 0,37,74,112 in legitimate 5'-sequence,
    // avoiding adapter-RC contamination that RC(mate2) seeds at positions 0,37,74 would hit.
    // mate2 reverse-cluster transcripts (is_reverse=true) = mate2 on - strand (FR pairs).
    // mate2 forward-cluster transcripts (is_reverse=false) = mate2 on + strand (RF pairs).
    let s1 = Seed::find_seeds(mate1_seq, index, params.seed_map_min, params, "")?;
    let s2 = Seed::find_seeds(mate2_seq, index, params.seed_map_min, params, "")?;

    let mate1_clusters = cluster_seeds(&s1, index, params, len1);
    let mate2_clusters = cluster_seeds(&s2, index, params, len2);

    if debug_pe {
        eprintln!(
            "[DEBUG-PE] Per-mate: s1={} s2={} m1_clusters={} m2_clusters={}",
            s1.len(),
            s2.len(),
            mate1_clusters.len(),
            mate2_clusters.len()
        );
        for (i, c) in mate1_clusters.iter().enumerate() {
            eprintln!(
                "  m1_cluster[{}]: rev={} chr={} seeds={}",
                i,
                c.is_reverse,
                c.chr_idx,
                c.alignments.len()
            );
        }
        for (i, c) in mate2_clusters.iter().enumerate() {
            eprintln!(
                "  m2_cluster[{}]: rev={} chr={} seeds={}",
                i,
                c.is_reverse,
                c.chr_idx,
                c.alignments.len()
            );
        }
    }

    // Stitch each mate's clusters independently.
    // For mate2 reverse clusters: stitch_seeds_with_jdb internally uses RC(mate2_seq) as
    // stitch_read and sets transcript.is_reverse=true + transcript.read_seq=mate2_seq.
    // This correctly represents mate2 on - strand for FR pairs without explicit RC handling.
    let debug_name: &str = if debug_pe { read_name } else { "" };
    let mut mate1_transcripts: Vec<Transcript> = Vec::new();
    for cluster in mate1_clusters
        .iter()
        .take(params.align_windows_per_read_nmax)
    {
        let ts = stitch_seeds_with_jdb_debug(
            cluster,
            mate1_seq,
            index,
            &scorer,
            junction_db,
            params.align_transcripts_per_window_nmax,
            debug_name,
        )?;
        mate1_transcripts.extend(ts);
    }

    let mut mate2_transcripts: Vec<Transcript> = Vec::new();
    for cluster in mate2_clusters
        .iter()
        .take(params.align_windows_per_read_nmax)
    {
        let ts = stitch_seeds_with_jdb_debug(
            cluster,
            mate2_seq,
            index,
            &scorer,
            junction_db,
            params.align_transcripts_per_window_nmax,
            debug_name,
        )?;
        mate2_transcripts.extend(ts);
    }

    // Combined score threshold: use len1+len2 as denominator (matches STAR's combined Lread-1).
    let combined_score_threshold =
        (params.out_filter_score_min_over_lread * (len1 + len2) as f64) as i32;

    // Pair compatible transcripts (opposite strand, same chr, within alignMatesGapMax)
    let mut joint_pairs: Vec<PairedAlignment> = Vec::new();
    for t1 in &mate1_transcripts {
        for t2 in &mate2_transcripts {
            if let Some(pair) =
                try_pair_transcripts(t1, t2, len1, len2, params, combined_score_threshold)
            {
                joint_pairs.push(pair);
            }
        }
    }

    if debug_pe {
        eprintln!(
            "[DEBUG-PE] After pairing: pairs={} m1_ts={} m2_ts={}",
            joint_pairs.len(),
            mate1_transcripts.len(),
            mate2_transcripts.len()
        );
        for (i, pa) in joint_pairs.iter().enumerate() {
            eprintln!(
                "  pair[{}]: M1 chr={} pos={} rev={} score={} | M2 chr={} pos={} rev={} score={} | combined={}",
                i,
                pa.mate1_transcript.chr_idx,
                pa.mate1_transcript.genome_start,
                pa.mate1_transcript.is_reverse,
                pa.mate1_transcript.score,
                pa.mate2_transcript.chr_idx,
                pa.mate2_transcript.genome_start,
                pa.mate2_transcript.is_reverse,
                pa.mate2_transcript.score,
                pa.combined_wt_score
            );
        }
    }

    // --- Decision tree: score-filter, dedup, quality-filter, then half-mapped fallback ---
    // Step 1: score-range filter (STAR's multMapSelect).
    if !joint_pairs.is_empty() {
        let best_score = joint_pairs
            .iter()
            .map(|pa| pa.combined_wt_score)
            .max()
            .unwrap_or(0);
        let score_threshold = best_score - params.out_filter_multimap_score_range;
        joint_pairs.retain(|pa| pa.combined_wt_score >= score_threshold);
    }

    // Step 2: TooManyLoci check using pre-dedup count.
    if joint_pairs.len() > params.out_filter_multimap_nmax as usize {
        return Ok((Vec::new(), 0, Some(UnmappedReason::TooManyLoci)));
    }

    // Step 3: position dedup — remove exact (chr, mate1_pos, mate2_pos, strand, CIGAR) duplicates.
    joint_pairs.sort_by(|a, b| {
        let pos_cmp = (
            a.mate1_transcript.chr_idx,
            a.mate1_transcript.genome_start,
            a.mate1_transcript.is_reverse,
            a.mate2_transcript.genome_start,
            a.mate2_transcript.is_reverse,
        )
            .cmp(&(
                b.mate1_transcript.chr_idx,
                b.mate1_transcript.genome_start,
                b.mate1_transcript.is_reverse,
                b.mate2_transcript.genome_start,
                b.mate2_transcript.is_reverse,
            ));
        if pos_cmp != std::cmp::Ordering::Equal {
            return pos_cmp;
        }
        b.combined_wt_score
            .cmp(&a.combined_wt_score)
            .then_with(|| a.mate1_transcript.cigar.cmp(&b.mate1_transcript.cigar))
            .then_with(|| a.mate2_transcript.cigar.cmp(&b.mate2_transcript.cigar))
    });
    joint_pairs.dedup_by(|a, b| {
        a.mate1_transcript.chr_idx == b.mate1_transcript.chr_idx
            && a.mate1_transcript.genome_start == b.mate1_transcript.genome_start
            && a.mate1_transcript.is_reverse == b.mate1_transcript.is_reverse
            && a.mate1_transcript.cigar == b.mate1_transcript.cigar
            && a.mate2_transcript.genome_start == b.mate2_transcript.genome_start
            && a.mate2_transcript.is_reverse == b.mate2_transcript.is_reverse
            && a.mate2_transcript.cigar == b.mate2_transcript.cigar
    });

    // Post-finalization mate2-exon-subset dedup.
    {
        use crate::align::transcript::Exon;
        let exons_subset = |b: &[Exon], a: &[Exon]| -> bool {
            let total_b: u32 = b.iter().map(|e| (e.read_end - e.read_start) as u32).sum();
            if total_b == 0 {
                return false;
            }
            let mut covered = 0u32;
            for be in b {
                let b_diag = be.genome_start as i64 - be.read_start as i64;
                for ae in a {
                    let a_diag = ae.genome_start as i64 - ae.read_start as i64;
                    if a_diag == b_diag {
                        let r_start = be.read_start.max(ae.read_start);
                        let r_end = be.read_end.min(ae.read_end);
                        if r_start < r_end {
                            covered += (r_end - r_start) as u32;
                        }
                    }
                }
            }
            covered == total_b
        };
        let n = joint_pairs.len();
        let mut keep = vec![true; n];
        for i in 0..n {
            if !keep[i] {
                continue;
            }
            for j in 0..n {
                if i == j || !keep[j] {
                    continue;
                }
                let same_pos = joint_pairs[i].mate1_transcript.chr_idx
                    == joint_pairs[j].mate1_transcript.chr_idx
                    && joint_pairs[i].mate1_transcript.genome_start
                        == joint_pairs[j].mate1_transcript.genome_start
                    && joint_pairs[i].mate1_transcript.genome_end
                        == joint_pairs[j].mate1_transcript.genome_end
                    && joint_pairs[i].mate1_transcript.is_reverse
                        == joint_pairs[j].mate1_transcript.is_reverse
                    && joint_pairs[i].mate2_transcript.genome_start
                        == joint_pairs[j].mate2_transcript.genome_start
                    && joint_pairs[i].mate2_transcript.is_reverse
                        == joint_pairs[j].mate2_transcript.is_reverse;
                if same_pos
                    && joint_pairs[i].combined_wt_score < joint_pairs[j].combined_wt_score
                    && exons_subset(
                        &joint_pairs[i].mate2_transcript.exons,
                        &joint_pairs[j].mate2_transcript.exons,
                    )
                {
                    keep[i] = false;
                    break;
                }
            }
        }
        joint_pairs = joint_pairs
            .into_iter()
            .enumerate()
            .filter(|(i, _)| keep[*i])
            .map(|(_, p)| p)
            .collect();
    }

    joint_pairs.sort_by(|a, b| {
        b.combined_wt_score
            .cmp(&a.combined_wt_score)
            .then_with(|| a.mate1_transcript.chr_idx.cmp(&b.mate1_transcript.chr_idx))
            .then_with(|| {
                a.mate1_transcript
                    .genome_start
                    .cmp(&b.mate1_transcript.genome_start)
            })
            .then_with(|| {
                a.mate1_transcript
                    .is_reverse
                    .cmp(&b.mate1_transcript.is_reverse)
            })
    });

    // Randomize primary among best-scoring pairs (STAR's funPrimaryAlignMark).
    shuffle_tied_prefix(
        &mut joint_pairs,
        |pa| pa.combined_wt_score,
        per_read_seed(params.run_rng_seed, read_name),
    );

    // Step 4: quality filter (mappedFilter).
    filter_paired_transcripts(&mut joint_pairs, params);

    if !joint_pairs.is_empty() {
        let pe_mapq_n = joint_pairs.len().max(1);
        let results = joint_pairs
            .into_iter()
            .map(|pa| PairedAlignmentResult::BothMapped(Box::new(pa)))
            .collect();
        return Ok((results, pe_mapq_n, None));
    }

    // Half-mapped fallback: report the best-scoring single-mate transcript.
    // Apply per-mate quality threshold (outFilterScoreMinOverLread * (len - 1)).
    let thresh1 = ((params.out_filter_score_min_over_lread * (len1 as f64 - 1.0)) as i32)
        .max(params.out_filter_score_min);
    let thresh2 = ((params.out_filter_score_min_over_lread * (len2 as f64 - 1.0)) as i32)
        .max(params.out_filter_score_min);

    let best_m1 = mate1_transcripts
        .into_iter()
        .filter(|t| t.score >= thresh1)
        .max_by_key(|t| t.score);
    let best_m2 = mate2_transcripts
        .into_iter()
        .filter(|t| t.score >= thresh2)
        .max_by_key(|t| t.score);

    match (best_m1, best_m2) {
        (Some(t1), None) => Ok((
            vec![PairedAlignmentResult::HalfMapped {
                mapped_transcript: t1,
                mate1_is_mapped: true,
            }],
            1,
            None,
        )),
        (None, Some(t2)) => Ok((
            vec![PairedAlignmentResult::HalfMapped {
                mapped_transcript: t2,
                mate1_is_mapped: false,
            }],
            1,
            None,
        )),
        (Some(t1), Some(t2)) => {
            // Both have single-mate alignments but couldn't form a valid pair.
            // Report the higher-scoring mate as half-mapped.
            if t1.score >= t2.score {
                Ok((
                    vec![PairedAlignmentResult::HalfMapped {
                        mapped_transcript: t1,
                        mate1_is_mapped: true,
                    }],
                    1,
                    None,
                ))
            } else {
                Ok((
                    vec![PairedAlignmentResult::HalfMapped {
                        mapped_transcript: t2,
                        mate1_is_mapped: false,
                    }],
                    1,
                    None,
                ))
            }
        }
        (None, None) => Ok((Vec::new(), 0, Some(UnmappedReason::TooShort))),
    }
}

/// Attempt to pair two per-mate transcripts into a PairedAlignment.
///
/// Returns `None` if the mates are incompatible (same strand, different chr, too far, etc.).
fn try_pair_transcripts(
    t1: &Transcript,
    t2: &Transcript,
    len1: usize,
    len2: usize,
    params: &Parameters,
    combined_score_threshold: i32,
) -> Option<PairedAlignment> {
    // Must be same chromosome
    if t1.chr_idx != t2.chr_idx {
        return None;
    }
    // Must be opposite strands (FR or RF)
    if t1.is_reverse == t2.is_reverse {
        return None;
    }

    // Determine left mate (smaller genome_start) and right mate for distance/consistency checks
    let (left, right) = if t1.genome_start <= t2.genome_start {
        (t1, t2)
    } else {
        (t2, t1)
    };

    // Reject degenerate pairs (right ends before left starts → negative insert)
    if right.genome_end <= left.genome_start {
        return None;
    }

    // Genomic span check: use alignMatesGapMax if set, else fall back to win_bin_window_dist
    // (STAR's effective limit when alignMatesGapMax=0 is the window distance ~589kb)
    let span = right.genome_end - left.genome_start;
    let max_span = if params.align_mates_gap_max > 0 {
        params.align_mates_gap_max as u64
    } else {
        params.win_bin_window_dist()
    };
    if span > max_span {
        return None;
    }

    // Combined score: sum of per-mate finalized scores (each already includes per-mate span penalty)
    let combined_wt_score = t1.score + t2.score;

    // SCORE-GATE: reject pairs where score is below the absolute floor
    if combined_wt_score + params.out_filter_multimap_score_range < combined_score_threshold {
        return None;
    }

    // Junction consistency in overlap region
    if !pe_junctions_consistent(left, right) {
        return None;
    }

    // Combined coverage: sum of exon read spans from both mates
    let combined_n_match: u32 = t1
        .exons
        .iter()
        .map(|e| (e.read_end - e.read_start) as u32)
        .sum::<u32>()
        + t2.exons
            .iter()
            .map(|e| (e.read_end - e.read_start) as u32)
            .sum::<u32>();

    let is_proper_pair = check_proper_pair(t1, t2, params);
    let insert_size = calculate_insert_size(t1, t2);

    Some(PairedAlignment {
        mate1_transcript: t1.clone(),
        mate2_transcript: t2.clone(),
        mate1_region: (0, len1),
        mate2_region: (0, len2),
        is_proper_pair,
        insert_size,
        combined_wt_score,
        combined_n_match,
    })
}

/// Check if paired alignment is a proper pair
fn check_proper_pair(
    mate1_trans: &Transcript,
    mate2_trans: &Transcript,
    params: &Parameters,
) -> bool {
    // Proper pair criteria:
    // 1. Both mates mapped (checked by caller)
    // 2. Same chromosome (checked by caller)
    // 3. Distance within alignMatesGapMax

    if params.align_mates_gap_max == 0 {
        return true; // Auto mode = unlimited
    }

    // Calculate genomic distance
    let start = mate1_trans.genome_start.min(mate2_trans.genome_start);
    let end = mate1_trans.genome_end.max(mate2_trans.genome_end);
    let genomic_span = end - start;

    genomic_span <= params.align_mates_gap_max as u64
}

/// Calculate signed insert size (TLEN)
fn calculate_insert_size(mate1_trans: &Transcript, mate2_trans: &Transcript) -> i32 {
    // STAR outSAMtlen=1 (default): tlen is computed from the combined PE transcript span,
    // not from max/min of individual mate endpoints.
    //
    // Forward cluster (mate1 forward, mate2 reverse — FR pair):
    //   tlen = mate2.genome_end - mate1.genome_start
    //   mate1 (imate=0 in combined transcript) gets +tlen
    //
    // Reverse cluster (mate1 reverse, mate2 forward — RF pair):
    //   tlen = mate1.genome_end - mate2.genome_start
    //   mate2 (imate=0 in combined transcript) gets +tlen, so mate1 gets -tlen
    //
    // This correctly handles:
    //   - Same-start overlapping pairs (uses trailing mate's end, not max of both)
    //   - RF pairs at same position (mate2 gets +tlen, not mate1)
    if !mate1_trans.is_reverse {
        // Forward cluster: mate1 on left
        (mate2_trans.genome_end - mate1_trans.genome_start) as i32
    } else {
        // Reverse cluster: mate2 on left, so mate1 gets negative tlen
        -((mate1_trans.genome_end - mate2_trans.genome_start) as i32)
    }
}

/// Filter paired transcripts by quality thresholds.
/// STAR's mappedFilter (ReadAlign_mappedFilter.cpp) applies ALL quality checks to the combined
/// read as a single unit (Lread = len1+1+len2). There are no per-mate quality filters for PE.
fn filter_paired_transcripts(paired_alns: &mut Vec<PairedAlignment>, params: &Parameters) {
    // STAR's mappedFilter (ReadAlign_mappedFilter.cpp) applies ALL quality thresholds to trBest
    // (the highest-scoring transcript), NOT to each individual transcript. If trBest passes,
    // all transcripts in the score window are included (they affect NH/MAPQ). If trBest fails,
    // the read is unmapped.
    //
    // Step 1: find the best pair and check quality thresholds on it.
    let best_pa = paired_alns.iter().max_by_key(|pa| pa.combined_wt_score);
    if let Some(best) = best_pa {
        let mate1_len = (best.mate1_region.1 - best.mate1_region.0) as f64;
        let mate2_len = (best.mate2_region.1 - best.mate2_region.0) as f64;
        // Lread-1 for the combined read: (len1+1+len2) - 1 = len1 + len2
        let combined_lread_m1 = mate1_len + mate2_len;

        let combined_nm = best.mate1_transcript.n_mismatch + best.mate2_transcript.n_mismatch;

        // Score: trBest->maxScore < outFilterScoreMin || < outFilterScoreMinOverLread*(Lread-1)
        // STAR checks the combined WT score (ReadAlign_mappedFilter.cpp), not the sum of per-mate
        // finalized scores. Using t1.score + t2.score incorrectly double-applies the genomic
        // length penalty (once per mate), rejecting valid overlapping short-insert pairs where
        // per-mate spans are small (penalty ≈ −2 each → sum penalty −4 vs combined penalty −2).
        let combined_score = best.combined_wt_score;
        if combined_score < params.out_filter_score_min
            || combined_score < (params.out_filter_score_min_over_lread * combined_lread_m1) as i32
        {
            paired_alns.clear();
            return;
        }

        // Match bases: trBest->nMatch < outFilterMatchNmin || < outFilterMatchNminOverLread*(Lread-1)
        // STAR's mappedFilter applies nMatch against Lread-1 on the JOINT combined transcript.
        // For overlapping pairs (both mates at same genome region), STAR's joint transcript covers
        // only one mate's worth → nMatch ≈ 137 < 198.
        // Fix: use the pre-split combined-read coverage captured before split_working_transcript,
        // which directly mirrors STAR's joint transcript nMatch without extension inflation.
        let combined_match = best.combined_n_match;
        if combined_match < params.out_filter_match_nmin
            || combined_match < (params.out_filter_match_nmin_over_lread * combined_lread_m1) as u32
        {
            paired_alns.clear();
            return;
        }

        // Mismatches: nMM > outFilterMismatchNmaxTotal || nMM/rLength > outFilterMismatchNoverLmax
        // outFilterMismatchNmaxTotal = min(outFilterMismatchNmax, outFilterMismatchNoverReadLmax*(len1+len2))
        if combined_nm > params.out_filter_mismatch_nmax
            || (combined_nm as f64)
                > params.out_filter_mismatch_nover_lmax * (mate1_len + mate2_len)
        {
            paired_alns.clear();
            return;
        }
    }

    // Step 2: STAR's mappedFilter: if nTr > outFilterMultimapNmax → unmapped=multi (clear all).
    // Unlike SE which truncates for output, PE must REJECT the entire pair when over the limit.
    if paired_alns.len() > params.out_filter_multimap_nmax as usize {
        paired_alns.clear();
    }
}

/// Extract splice junctions from a Transcript's CIGAR as (donor, acceptor) pairs.
/// Junction coords are in genomic space (0-based). Junctions only exist where CigarOp::RefSkip is.
fn extract_junctions_from_cigar(t: &Transcript) -> Vec<(u64, u64)> {
    let mut junctions = Vec::new();
    let mut genome_pos = t.genome_start;
    for op in &t.cigar {
        use crate::align::transcript::CigarOp;
        match op {
            CigarOp::Match(n) | CigarOp::Equal(n) | CigarOp::Diff(n) | CigarOp::Del(n) => {
                genome_pos += *n as u64;
            }
            CigarOp::RefSkip(n) => {
                let donor = genome_pos;
                let acceptor = genome_pos + *n as u64;
                junctions.push((donor, acceptor));
                genome_pos = acceptor;
            }
            CigarOp::Ins(_) | CigarOp::SoftClip(_) | CigarOp::HardClip(_) => {}
        }
    }
    junctions
}

/// D5: Check junction consistency in the overlapping region of paired-end mates.
/// When mates overlap in the genome, every splice junction in the overlap from the
/// left mate must appear in the right mate, and vice versa.
/// Implements STAR stitchWindowAligns.cpp check after overlap detection.
///
/// `left` is the mate with the lower genome_start, `right` is the other mate.
pub(crate) fn pe_junctions_consistent(left: &Transcript, right: &Transcript) -> bool {
    // Overlapping region: [overlap_start, overlap_end)
    let overlap_start = left.genome_start.max(right.genome_start);
    let overlap_end = left.genome_end.min(right.genome_end);
    if overlap_start >= overlap_end {
        return true; // No overlap — nothing to check
    }

    let left_juncs = extract_junctions_from_cigar(left);
    let right_juncs = extract_junctions_from_cigar(right);

    // Every junction from left that falls within the overlap must be in right too
    for (donor, acceptor) in &left_juncs {
        if *donor >= overlap_start
            && *acceptor <= overlap_end
            && !right_juncs.iter().any(|(d, a)| d == donor && a == acceptor)
        {
            return false;
        }
    }
    // Every junction from right that falls within the overlap must be in left too
    for (donor, acceptor) in &right_juncs {
        if *donor >= overlap_start
            && *acceptor <= overlap_end
            && !left_juncs.iter().any(|(d, a)| d == donor && a == acceptor)
        {
            return false;
        }
    }
    true
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::genome::Genome;
    use crate::index::packed_array::PackedArray;
    use crate::index::sa_index::SaIndex;
    use crate::index::suffix_array::SuffixArray;
    use clap::Parser;

    fn make_test_params() -> Parameters {
        // Parse empty args to get default parameters
        Parameters::try_parse_from(vec!["ruSTAR"]).unwrap()
    }

    fn make_test_index() -> GenomeIndex {
        // Simple genome: ACGTACGTNN (10 bases)
        let seq = vec![0, 1, 2, 3, 0, 1, 2, 3, 4, 4];
        let n_genome = 64u64; // Padded
        let mut sequence = vec![5u8; (n_genome * 2) as usize];
        sequence[0..seq.len()].copy_from_slice(&seq);

        // Build reverse complement
        for i in 0..n_genome as usize {
            let base = sequence[i];
            let complement = if base < 4 { 3 - base } else { base };
            sequence[2 * n_genome as usize - 1 - i] = complement;
        }

        let genome = Genome {
            sequence,
            n_genome,
            n_chr_real: 1,
            chr_name: vec!["chr1".to_string()],
            chr_length: vec![10],
            chr_start: vec![0, n_genome],
        };

        // Create dummy SA and SAindex (would need real index for actual alignment)
        let gstrand_bit = 33;
        let suffix_array = SuffixArray {
            data: PackedArray::new(gstrand_bit, 0),
            gstrand_bit,
            gstrand_mask: (1u64 << gstrand_bit) - 1,
        };

        let word_length = gstrand_bit + 3;
        let sa_index = SaIndex {
            data: PackedArray::new(word_length, 0),
            nbases: 14,
            genome_sa_index_start: vec![0],
            word_length,
            gstrand_bit,
        };

        GenomeIndex {
            genome,
            suffix_array,
            sa_index,
            junction_db: crate::junction::SpliceJunctionDb::empty(),
        }
    }

    #[test]
    fn test_align_read_no_seeds() {
        let index = make_test_index();
        let params = make_test_params();

        // Read with all N's (no seeds possible)
        let read_seq = vec![4, 4, 4, 4, 4, 4, 4, 4, 4, 4];

        let result = align_read(&read_seq, "READ_001", &index, &params);
        assert!(result.is_ok());

        let (transcripts, chimeras, n_for_mapq, unmapped_reason) = result.unwrap();
        assert_eq!(transcripts.len(), 0); // No alignment
        assert_eq!(chimeras.len(), 0); // No chimeric alignments
        assert_eq!(n_for_mapq, 0);
        assert_eq!(unmapped_reason, Some(UnmappedReason::Other));
    }

    #[test]
    fn test_transcript_filtering_score() {
        let index = make_test_index();
        let mut params = make_test_params();
        params.out_filter_score_min = 50;

        // Would need actual seeds and alignment to test this properly
        // This test just verifies the function doesn't crash
        let read_seq = vec![0, 1, 2, 3]; // ACGT
        let result = align_read(&read_seq, "READ_002", &index, &params);
        assert!(result.is_ok());
    }

    #[test]
    fn test_transcript_filtering_mismatch() {
        let index = make_test_index();
        let mut params = make_test_params();
        params.out_filter_mismatch_nmax = 2;

        let read_seq = vec![0, 1, 2, 3]; // ACGT
        let result = align_read(&read_seq, "READ_003", &index, &params);
        assert!(result.is_ok());
    }

    #[test]
    fn test_transcript_multimap_limit() {
        let index = make_test_index();
        let mut params = make_test_params();
        params.out_filter_multimap_nmax = 5;

        let read_seq = vec![0, 1, 2, 3]; // ACGT
        let result = align_read(&read_seq, "READ_004", &index, &params);
        assert!(result.is_ok());

        let (transcripts, _chimeras, _n_for_mapq, _reason) = result.unwrap();
        assert!(transcripts.len() <= 5);
    }

    #[test]
    fn test_align_paired_read_no_seeds() {
        let index = make_test_index();
        let params = make_test_params();

        // Both mates with all N's
        let mate1 = vec![4, 4, 4, 4, 4, 4, 4, 4];
        let mate2 = vec![4, 4, 4, 4, 4, 4, 4, 4];

        let result = align_paired_read(&mate1, &mate2, "test", &index, &params);
        assert!(result.is_ok());
        let (paired_alns, n_for_mapq, unmapped_reason) = result.unwrap();
        assert_eq!(paired_alns.len(), 0);
        assert_eq!(n_for_mapq, 0);
        assert!(unmapped_reason.is_some());
    }

    #[test]
    fn test_check_proper_pair_distance() {
        use crate::align::transcript::{CigarOp, Exon};

        let params = make_test_params();

        // Create two transcripts on same chromosome
        let t1 = Transcript {
            chr_idx: 0,
            genome_start: 1000,
            genome_end: 1100,
            is_reverse: false,
            exons: vec![Exon {
                genome_start: 1000,
                genome_end: 1100,
                read_start: 0,
                read_end: 100,
            }],
            cigar: vec![CigarOp::Match(100)],
            score: 100,
            n_mismatch: 0,
            n_gap: 0,
            n_junction: 0,
            junction_motifs: vec![],
            junction_annotated: vec![],
            read_seq: vec![0; 100],
        };

        let t2 = Transcript {
            chr_idx: 0,
            genome_start: 1200,
            genome_end: 1300,
            is_reverse: true,
            exons: vec![Exon {
                genome_start: 1200,
                genome_end: 1300,
                read_start: 0,
                read_end: 100,
            }],
            cigar: vec![CigarOp::Match(100)],
            score: 100,
            n_mismatch: 0,
            n_gap: 0,
            n_junction: 0,
            junction_motifs: vec![],
            junction_annotated: vec![],
            read_seq: vec![0; 100],
        };

        // Distance = 300bp, within default limit (auto mode = unlimited)
        assert!(check_proper_pair(&t1, &t2, &params));
    }

    #[test]
    fn test_check_proper_pair_too_far() {
        use crate::align::transcript::{CigarOp, Exon};

        let mut params = make_test_params();
        params.align_mates_gap_max = 100;

        let t1 = Transcript {
            chr_idx: 0,
            genome_start: 1000,
            genome_end: 1100,
            is_reverse: false,
            exons: vec![Exon {
                genome_start: 1000,
                genome_end: 1100,
                read_start: 0,
                read_end: 100,
            }],
            cigar: vec![CigarOp::Match(100)],
            score: 100,
            n_mismatch: 0,
            n_gap: 0,
            n_junction: 0,
            junction_motifs: vec![],
            junction_annotated: vec![],
            read_seq: vec![0; 100],
        };

        let t2 = Transcript {
            chr_idx: 0,
            genome_start: 1300,
            genome_end: 1400,
            is_reverse: true,
            exons: vec![Exon {
                genome_start: 1300,
                genome_end: 1400,
                read_start: 0,
                read_end: 100,
            }],
            cigar: vec![CigarOp::Match(100)],
            score: 100,
            n_mismatch: 0,
            n_gap: 0,
            n_junction: 0,
            junction_motifs: vec![],
            junction_annotated: vec![],
            read_seq: vec![0; 100],
        };

        // Distance = 400bp, exceeds limit of 100bp
        assert!(!check_proper_pair(&t1, &t2, &params));
    }

    #[test]
    fn test_calculate_insert_size_positive() {
        use crate::align::transcript::{CigarOp, Exon};

        // Mate1 is leftmost
        let t1 = Transcript {
            chr_idx: 0,
            genome_start: 1000,
            genome_end: 1100,
            is_reverse: false,
            exons: vec![Exon {
                genome_start: 1000,
                genome_end: 1100,
                read_start: 0,
                read_end: 100,
            }],
            cigar: vec![CigarOp::Match(100)],
            score: 100,
            n_mismatch: 0,
            n_gap: 0,
            n_junction: 0,
            junction_motifs: vec![],
            junction_annotated: vec![],
            read_seq: vec![0; 100],
        };

        let t2 = Transcript {
            chr_idx: 0,
            genome_start: 1200,
            genome_end: 1300,
            is_reverse: true,
            exons: vec![Exon {
                genome_start: 1200,
                genome_end: 1300,
                read_start: 0,
                read_end: 100,
            }],
            cigar: vec![CigarOp::Match(100)],
            score: 100,
            n_mismatch: 0,
            n_gap: 0,
            n_junction: 0,
            junction_motifs: vec![],
            junction_annotated: vec![],
            read_seq: vec![0; 100],
        };

        let tlen = calculate_insert_size(&t1, &t2);
        assert_eq!(tlen, 300); // Positive because mate1 is leftmost
    }

    #[test]
    fn test_strand_consistency_filter() {
        use crate::align::transcript::{CigarOp, Exon, Transcript};
        use crate::params::IntronStrandFilter;

        // Create a transcript with conflicting strand motifs (mixed + and - within one transcript)
        let t_inconsistent = Transcript {
            chr_idx: 0,
            genome_start: 1000,
            genome_end: 1300,
            is_reverse: false,
            exons: vec![Exon {
                genome_start: 1000,
                genome_end: 1300,
                read_start: 0,
                read_end: 100,
            }],
            cigar: vec![CigarOp::Match(100)],
            score: 100,
            n_mismatch: 0,
            n_gap: 0,
            n_junction: 2,
            junction_motifs: vec![SpliceMotif::GtAg, SpliceMotif::CtAc], // +strand and -strand
            junction_annotated: vec![],
            read_seq: vec![0; 100],
        };

        // Create a transcript with consistent strand motifs (all + strand)
        let t_consistent = Transcript {
            chr_idx: 0,
            genome_start: 1000,
            genome_end: 1300,
            is_reverse: false,
            exons: vec![Exon {
                genome_start: 1000,
                genome_end: 1300,
                read_start: 0,
                read_end: 100,
            }],
            cigar: vec![CigarOp::Match(100)],
            score: 100,
            n_mismatch: 0,
            n_gap: 0,
            n_junction: 2,
            junction_motifs: vec![SpliceMotif::GtAg, SpliceMotif::GcAg], // both + strand
            junction_annotated: vec![],
            read_seq: vec![0; 100],
        };

        // Note: STAR's RemoveInconsistentStrands filters transcripts where
        // junctions have MIXED implied strand (both + and - within one transcript).
        // It does NOT compare junction strand vs alignment strand — a reverse-strand
        // read at a GT/AG junction (antisense of + strand gene) is valid and kept.

        // Verify mixed-strand is detected
        let mut has_plus = false;
        let mut has_minus = false;
        for motif in &t_inconsistent.junction_motifs {
            match motif.implied_strand() {
                Some('+') => has_plus = true,
                Some('-') => has_minus = true,
                _ => {}
            }
        }
        assert!(has_plus && has_minus); // Inconsistent (mixed strands)

        // Verify consistent transcript has no conflict
        has_plus = false;
        has_minus = false;
        for motif in &t_consistent.junction_motifs {
            match motif.implied_strand() {
                Some('+') => has_plus = true,
                Some('-') => has_minus = true,
                _ => {}
            }
        }
        assert!(has_plus && !has_minus); // Consistent (all +)

        // Also verify: single CT/AC on forward-strand or single GT/AG on reverse-strand
        // are NOT filtered (STAR keeps these — antisense reads are valid)
        let ctac_only = vec![SpliceMotif::CtAc];
        let (mut hp, mut hm) = (false, false);
        for m in &ctac_only {
            match m.implied_strand() {
                Some('+') => hp = true,
                Some('-') => hm = true,
                _ => {}
            }
        }
        assert!(!hp && hm); // Only minus → NOT mixed → NOT filtered

        // Verify the filter enum
        assert_ne!(
            IntronStrandFilter::None,
            IntronStrandFilter::RemoveInconsistentStrands
        );
    }

    #[test]
    fn test_calculate_insert_size_negative() {
        use crate::align::transcript::{CigarOp, Exon};

        // RF pair: mate1 reverse (right side), mate2 forward (left side).
        // STAR reverse cluster: tlen = mate1.genome_end - mate2.genome_start = 1300 - 1000 = 300.
        // mate1 gets -tlen (negative) because mate2 (imate=0 in reverse cluster) is the left mate.
        let t1 = Transcript {
            chr_idx: 0,
            genome_start: 1200,
            genome_end: 1300,
            is_reverse: true, // mate1 is reverse strand → reverse cluster
            exons: vec![Exon {
                genome_start: 1200,
                genome_end: 1300,
                read_start: 0,
                read_end: 100,
            }],
            cigar: vec![CigarOp::Match(100)],
            score: 100,
            n_mismatch: 0,
            n_gap: 0,
            n_junction: 0,
            junction_motifs: vec![],
            junction_annotated: vec![],
            read_seq: vec![0; 100],
        };

        let t2 = Transcript {
            chr_idx: 0,
            genome_start: 1000,
            genome_end: 1100,
            is_reverse: false,
            exons: vec![Exon {
                genome_start: 1000,
                genome_end: 1100,
                read_start: 0,
                read_end: 100,
            }],
            cigar: vec![CigarOp::Match(100)],
            score: 100,
            n_mismatch: 0,
            n_gap: 0,
            n_junction: 0,
            junction_motifs: vec![],
            junction_annotated: vec![],
            read_seq: vec![0; 100],
        };

        let tlen = calculate_insert_size(&t1, &t2);
        assert_eq!(tlen, -300); // Negative for mate1 in RF pair (mate2 is the left mate)
    }

    #[test]
    fn test_noncanonical_unannotated_filter() {
        use crate::align::score::SpliceMotif;
        use crate::align::transcript::{CigarOp, Exon, Transcript};

        // Helper: check if a transcript would be filtered by RemoveNoncanonicalUnannotated
        // (mirrors the logic in the retain closure)
        let would_filter = |t: &Transcript| -> bool {
            t.junction_motifs
                .iter()
                .zip(t.junction_annotated.iter())
                .any(|(m, annotated)| *m == SpliceMotif::NonCanonical && !annotated)
        };

        let base_transcript = || Transcript {
            chr_idx: 0,
            genome_start: 1000,
            genome_end: 1200,
            is_reverse: false,
            exons: vec![Exon {
                genome_start: 1000,
                genome_end: 1200,
                read_start: 0,
                read_end: 100,
            }],
            cigar: vec![CigarOp::Match(100)],
            score: 100,
            n_mismatch: 0,
            n_gap: 0,
            n_junction: 1,
            junction_motifs: vec![],
            junction_annotated: vec![],
            read_seq: vec![0; 100],
        };

        // Case 1: NonCanonical + unannotated → should be filtered
        let mut t1 = base_transcript();
        t1.junction_motifs = vec![SpliceMotif::NonCanonical];
        t1.junction_annotated = vec![false];
        assert!(
            would_filter(&t1),
            "NonCanonical + unannotated should be filtered"
        );

        // Case 2: NonCanonical + annotated → should be KEPT
        let mut t2 = base_transcript();
        t2.junction_motifs = vec![SpliceMotif::NonCanonical];
        t2.junction_annotated = vec![true];
        assert!(
            !would_filter(&t2),
            "NonCanonical + annotated should be kept"
        );

        // Case 3: Canonical + unannotated → should be KEPT
        let mut t3 = base_transcript();
        t3.junction_motifs = vec![SpliceMotif::GtAg];
        t3.junction_annotated = vec![false];
        assert!(!would_filter(&t3), "Canonical + unannotated should be kept");

        // Case 4: Mixed — one canonical + one non-canonical unannotated → filtered
        let mut t4 = base_transcript();
        t4.junction_motifs = vec![SpliceMotif::GtAg, SpliceMotif::NonCanonical];
        t4.junction_annotated = vec![true, false];
        assert!(
            would_filter(&t4),
            "Mixed with unannotated non-canonical should be filtered"
        );

        // Case 5: Mixed — one canonical + one non-canonical annotated → kept
        let mut t5 = base_transcript();
        t5.junction_motifs = vec![SpliceMotif::GtAg, SpliceMotif::NonCanonical];
        t5.junction_annotated = vec![false, true];
        assert!(
            !would_filter(&t5),
            "Mixed with annotated non-canonical should be kept"
        );
    }

    #[test]
    fn test_align_paired_both_unmapped() {
        // Both mates are all N's → both unmapped → empty Vec
        let index = make_test_index();
        let params = make_test_params();

        let mate1 = vec![4, 4, 4, 4, 4, 4, 4, 4];
        let mate2 = vec![4, 4, 4, 4, 4, 4, 4, 4];

        let (results, n_for_mapq, unmapped_reason) =
            align_paired_read(&mate1, &mate2, "test", &index, &params).unwrap();
        assert!(results.is_empty(), "Both unmapped should return empty Vec");
        assert_eq!(n_for_mapq, 0);
        assert!(unmapped_reason.is_some(), "Should have unmapped reason");
    }

    #[test]
    fn test_paired_alignment_result_enum_variants() {
        use crate::align::transcript::{CigarOp, Exon};

        let transcript = Transcript {
            chr_idx: 0,
            genome_start: 1000,
            genome_end: 1100,
            is_reverse: false,
            exons: vec![Exon {
                genome_start: 1000,
                genome_end: 1100,
                read_start: 0,
                read_end: 100,
            }],
            cigar: vec![CigarOp::Match(100)],
            score: 100,
            n_mismatch: 0,
            n_gap: 0,
            n_junction: 0,
            junction_motifs: vec![],
            junction_annotated: vec![],
            read_seq: vec![0; 100],
        };

        // Test BothMapped variant
        let both = PairedAlignmentResult::BothMapped(Box::new(PairedAlignment {
            mate1_transcript: transcript.clone(),
            mate2_transcript: transcript.clone(),
            mate1_region: (0, 100),
            mate2_region: (0, 100),
            is_proper_pair: true,
            insert_size: 200,
            combined_wt_score: 0,
            combined_n_match: 200,
        }));
        assert!(matches!(both, PairedAlignmentResult::BothMapped(_)));

        // Test HalfMapped variant
        let half = PairedAlignmentResult::HalfMapped {
            mapped_transcript: transcript,
            mate1_is_mapped: true,
        };
        assert!(matches!(half, PairedAlignmentResult::HalfMapped { .. }));

        // Verify mate1_is_mapped
        if let PairedAlignmentResult::HalfMapped {
            mate1_is_mapped, ..
        } = half
        {
            assert!(mate1_is_mapped);
        }
    }

    #[test]
    fn shuffle_tied_prefix_is_deterministic() {
        // Same seed + same input → same permutation on reruns.
        let items: Vec<(i32, u32)> = (0..8).map(|i| (100, i)).collect();
        let mut a = items.clone();
        let mut b = items.clone();
        shuffle_tied_prefix(&mut a, |t| t.0, 12345);
        shuffle_tied_prefix(&mut b, |t| t.0, 12345);
        assert_eq!(a, b);
    }

    #[test]
    fn shuffle_tied_prefix_respects_ties() {
        // Only the top-score prefix gets shuffled; lower-scored tail is left alone.
        let mut items = vec![(100, 0u32), (100, 1), (100, 2), (50, 3), (40, 4)];
        shuffle_tied_prefix(&mut items, |t| t.0, 777);
        // Last two elements (non-tied) stay in place.
        assert_eq!(items[3], (50, 3));
        assert_eq!(items[4], (40, 4));
        // Tied prefix contains the original three items in some order.
        let mut top: Vec<u32> = items[..3].iter().map(|t| t.1).collect();
        top.sort();
        assert_eq!(top, vec![0, 1, 2]);
    }

    #[test]
    fn shuffle_tied_prefix_different_seeds_can_diverge() {
        // Probabilistic: for a tied set of 8, at least two seeds should disagree
        // on the chosen primary. (Exhaustive over a small seed range is fine.)
        let base: Vec<(i32, u32)> = (0..8).map(|i| (100, i)).collect();
        let mut firsts = std::collections::HashSet::new();
        for seed in 0..32u64 {
            let mut v = base.clone();
            shuffle_tied_prefix(&mut v, |t| t.0, seed);
            firsts.insert(v[0].1);
        }
        assert!(
            firsts.len() >= 2,
            "expected different seeds to pick different primaries, got {:?}",
            firsts
        );
    }

    #[test]
    fn shuffle_tied_prefix_noop_when_no_ties() {
        let mut items = vec![(100, 0u32), (90, 1), (80, 2)];
        let before = items.clone();
        shuffle_tied_prefix(&mut items, |t| t.0, 42);
        assert_eq!(items, before);
    }
}