align/
align.rs

1//! Example: align a small batch of paired-end reads using the YARA mapper.
2//!
3//! Usage:
4//!
5//!     cargo run --example align -- <index-prefix>
6//!
7//! The `index_prefix` should point to a pre-built YARA index (e.g.,
8//! `ref/hla.fasta` if the index files are `ref/hla.fasta.yara.*`).
9
10use yara_mapper::{MapperOptions, ReadBatch, ReadEnd, Sensitivity, YaraMapper};
11
12fn main() {
13    let args: Vec<String> = std::env::args().collect();
14    if args.len() < 2 {
15        eprintln!("Usage: {} <index_prefix>", args[0]);
16        std::process::exit(1);
17    }
18    let index_prefix = &args[1];
19
20    // Configure the mapper.
21    let options = MapperOptions {
22        sensitivity: Sensitivity::Full,
23        threads: 2,
24        verbose: 1,
25        ..Default::default()
26    };
27
28    // Open the index.
29    eprintln!("Opening index: {index_prefix}");
30    let mapper = YaraMapper::open(index_prefix, &options).unwrap_or_else(|e| {
31        eprintln!("Failed to open index: {e}");
32        std::process::exit(1);
33    });
34
35    eprintln!("Loaded {} contigs", mapper.contig_count());
36    for (i, (name, len)) in mapper.contig_names().iter().zip(mapper.contig_lengths()).enumerate() {
37        eprintln!("  contig {i}: {name} ({len} bp)");
38        if i >= 5 {
39            eprintln!("  ... ({} more)", mapper.contig_count() - i - 1);
40            break;
41        }
42    }
43
44    // Create a synthetic read pair (these won't align to anything real).
45    let mut batch = ReadBatch::new();
46    batch
47        .push(
48            "test_read_1",
49            ReadEnd {
50                seq: b"ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGT",
51                qual: b"IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII",
52            },
53            ReadEnd {
54                seq: b"TGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCATGCA",
55                qual: b"IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII",
56            },
57        )
58        .expect("failed to add read pair");
59
60    // Map the reads.
61    eprintln!("Mapping {} read pair(s)...", batch.len());
62    match mapper.map_paired(&batch) {
63        Ok(records) => {
64            eprintln!("Got {} alignment record(s)", records.len());
65            for rec in &records {
66                println!(
67                    "pair={} read1={} contig={} pos={} rev={} unmapped={} secondary={} mapq={} nm={} cigar={} flag={}",
68                    rec.read_pair_index,
69                    rec.is_read1,
70                    rec.contig_id,
71                    rec.pos,
72                    rec.is_reverse,
73                    rec.is_unmapped,
74                    rec.is_secondary,
75                    rec.mapq,
76                    rec.nm,
77                    rec.cigar_string(),
78                    rec.flag,
79                );
80            }
81        }
82        Err(e) => {
83            eprintln!("Mapping failed: {e}");
84            std::process::exit(1);
85        }
86    }
87}
align/align.rs

align/
align.rs
