rustalign_cli/

inspect.rs

//! rustalign-inspect - Inspect and query FM-index

use anyhow::Result;
use clap::{Parser, Subcommand};

/// RustAlign index inspection options
#[derive(Debug, Parser)]
pub struct InspectOptions {
    /// Index basename
    #[clap(required = true)]
    pub index: String,

    /// Inspection subcommand
    #[clap(subcommand)]
    pub command: InspectCommand,
}

/// Inspection subcommands
#[derive(Debug, Subcommand)]
pub enum InspectCommand {
    /// Print index summary
    Summary,

    /// Extract reference sequences
    Extract,

    /// Print index statistics
    Stats,

    /// Test exact_match functionality
    Test,

    /// Debug position resolution
    DebugPos {
        /// Genome offset to resolve
        #[clap(long)]
        offset: Option<u64>,

        /// Test specific known positions
        #[clap(long)]
        test: bool,
    },
}

/// Run the rustalign-inspect tool
pub fn run_with_options(opts: InspectOptions) -> Result<()> {
    run_with_opts(opts)
}

/// Run the rustalign-inspect tool (legacy - parses args)
pub fn run() -> Result<()> {
    let opts = InspectOptions::parse();
    run_with_options(opts)
}

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

fn run_with_opts(opts: InspectOptions) -> Result<()> {
    let index = rustalign_fmindex::Ebwt::load(&opts.index)?;
    let params = index.params();

    match opts.command {
        InspectCommand::Summary => {
            println!("Index: {}", opts.index);
            println!("Reference length: {}", params.len);
            println!("BWT length: {}", params.bwt_len);
            println!("Number of sides: {}", params.num_sides);
            println!("Side size: {} bytes", params.side_sz);
            println!("Total size: {} bytes", params.total_size());
        }
        InspectCommand::Extract => {
            println!("Extract not yet implemented");
        }
        InspectCommand::Stats => {
            println!("EbwtParams:");
            println!("  len: {}", params.len);
            println!("  bwt_len: {}", params.bwt_len);
            println!("  bwt_sz: {}", params.bwt_sz);
            println!("  line_sz: {}", params.line_sz);
            println!("  side_sz: {}", params.side_sz);
            println!("  side_bwt_sz: {}", params.side_bwt_sz);
            println!("  side_bwt_len: {}", params.side_bwt_len);
            println!("  num_sides: {}", params.num_sides);
            println!("  off_rate: {}", params.off_rate);
            println!("  ftab_chars: {}", params.ftab_chars);
            println!("  ftab_len: {}", params.ftab_len);
            println!("  ebwt_tot_len: {}", params.ebwt_tot_len);
        }
        InspectCommand::Test => {
            use rustalign_common::Nuc;
            use rustalign_fmindex::FmIndex;

            println!("Testing exact_match...");

            // Try a pattern from the first read: GATTGCTCGC (10 chars)
            let pattern: Vec<Nuc> = "GATTGCTCGC"
                .chars()
                .map(|c| Nuc::from_ascii(c as u8).unwrap())
                .collect();
            println!("\nSearching for: GATTGCTCGC");

            match index.exact_match(&pattern) {
                Ok(positions) => {
                    println!("Found {} matches", positions.len());
                    if !positions.is_empty() {
                        println!(
                            "First 10 positions: {:?}",
                            &positions[..positions.len().min(10)]
                        );
                    }
                }
                Err(e) => {
                    println!("Error: {:?}", e);
                }
            }

            // Try 10 A's
            let pattern: Vec<Nuc> = vec![Nuc::A; 10];
            println!("\nSearching for: AAAAAAAAAA");

            match index.exact_match(&pattern) {
                Ok(positions) => {
                    println!("Found {} matches", positions.len());
                }
                Err(e) => {
                    println!("Error: {:?}", e);
                }
            }
        }
        InspectCommand::DebugPos { offset, test } => {
            println!("=== Debug Position Resolution ===");
            println!();

            // Print chromosome info
            println!("Chromosomes:");
            for (i, name) in index.refnames().iter().enumerate() {
                let len = index.plen().get(i).copied().unwrap_or(0);
                println!("  [{}] {} ({} bp)", i, name, len);
            }
            println!();

            // Print rstarts array (first 10 entries)
            println!("rstarts array (first 10 fragments):");
            let rstarts = index.rstarts();
            for i in 0..10.min(index.n_frag() as usize) {
                let frag_start = rstarts[i * 3];
                let text_id = rstarts[i * 3 + 1];
                let frag_off = rstarts[i * 3 + 2];
                println!(
                    "  [{}] frag_start={}, text_id={}, frag_off={}",
                    i, frag_start, text_id, frag_off
                );
            }
            println!("  ... (showing first 10 of {} fragments)", index.n_frag());
            println!();

            // Also print the last few fragments
            let n_frag = index.n_frag() as usize;
            if n_frag > 5 {
                println!("rstarts array (last 5 fragments):");
                for i in (n_frag - 5)..n_frag {
                    let frag_start = rstarts[i * 3];
                    let text_id = rstarts[i * 3 + 1];
                    let frag_off = rstarts[i * 3 + 2];
                    println!(
                        "  [{}] frag_start={}, text_id={}, frag_off={}",
                        i, frag_start, text_id, frag_off
                    );
                }
                println!();
            }

            // Test specific offset if provided
            if let Some(off) = offset {
                println!("Resolving offset {}:", off);
                match index.joined_to_text_off(100, off, true) {
                    Ok((chr_idx, text_off, chr_len)) => {
                        println!(
                            "  Result: chr_idx={}, text_off={}, chr_len={}",
                            chr_idx, text_off, chr_len
                        );
                        if (chr_idx as usize) < index.refnames().len() {
                            println!("  Chromosome name: {}", index.refnames()[chr_idx as usize]);
                        }
                    }
                    Err(e) => {
                        println!("  Error: {:?}", e);
                    }
                }
            }

            // Run tests if requested
            if test {
                println!("\n=== Testing known positions ===");
                println!("Note: The 'off' parameter should be the BWT offset (without N's),");
                println!("not the genome offset (with N's).");
                println!();

                // Test cases from critical errors:
                // For BWT offsets, we need to subtract the cumulative N count
                // Let's test with the actual BWT length vs genome length
                let bwt_len = index.params().len;
                println!("BWT length (len parameter): {}", bwt_len);

                // Test specific BWT row resolutions
                // Test row 0 (should map to start of first chromosome)
                println!("\nTesting BWT row 0:");
                match index.get_offset(0) {
                    Ok(off) => {
                        println!("  get_offset(0) = {}", off);
                        match index.joined_to_text_off(100, off, true) {
                            Ok((chr_idx, text_off, _chr_len)) => {
                                let chr_name = if (chr_idx as usize) < index.refnames().len() {
                                    index.refnames()[chr_idx as usize].clone()
                                } else {
                                    format!("?{}", chr_idx)
                                };
                                println!(
                                    "  joined_to_text_off({}, {}) = {}@{}",
                                    100, off, chr_name, text_off
                                );
                            }
                            Err(e) => println!("  Error in joined_to_text_off: {:?}", e),
                        }
                    }
                    Err(e) => println!("  Error in get_offset: {:?}", e),
                }

                // Test the zOff position (should return offset 0)
                println!("\nzOff = {}", index.z_off());
                match index.get_offset(index.z_off() as u64) {
                    Ok(off) => println!("  get_offset(zOff) = {} (should be 0)", off),
                    Err(e) => println!("  Error: {:?}", e),
                }

                // Find which fragment contains a specific BWT offset
                // For Chr1 position 11789656, the BWT offset should be the same (no N's before it)
                println!("\n=== Binary search trace for offset 11789656 ===");
                let target_off = 11789656u64;
                let rstarts = index.rstarts();
                let n_frag = index.n_frag() as usize;

                let mut top = 0usize;
                let mut bot = n_frag;
                let mut steps = 0;

                while top < bot && steps < 20 {
                    let mid = top + ((bot - top) >> 1);
                    let lower = rstarts[mid * 3] as u64;
                    let upper = if mid + 1 >= n_frag {
                        index.params().len
                    } else {
                        rstarts[(mid + 1) * 3] as u64
                    };

                    steps += 1;
                    if steps <= 10 {
                        println!(
                            "  step {}: mid={}, lower={}, upper={}, target={}",
                            steps, mid, lower, upper, target_off
                        );
                    }

                    if lower <= target_off && upper > target_off {
                        let tidx = rstarts[mid * 3 + 1];
                        let frag_off_base = rstarts[mid * 3 + 2];
                        let frag_off = target_off - lower;
                        let text_off = frag_off + frag_off_base as u64;

                        println!("  FOUND: fragment {}", mid);
                        println!("    tidx (chr_idx) = {}", tidx);
                        println!("    frag_off_base = {}", frag_off_base);
                        println!("    frag_off = {} - {} = {}", target_off, lower, frag_off);
                        println!(
                            "    text_off = {} + {} = {}",
                            frag_off, frag_off_base, text_off
                        );
                        break;
                    } else if lower <= target_off {
                        top = mid + 1;
                    } else {
                        bot = mid;
                    }
                }
            }
        }
    }

    Ok(())
}