|
| 1 | +//! ColumnKit cookbook: a custom per-locus analyzer that inherits the bounded |
| 2 | +//! memory contract for free. The only domain logic is `on_column` (~3 lines) — |
| 3 | +//! by running it through `run_bounded_whole_genome` it gets the SAME bounded |
| 4 | +//! per-contig whole-genome walk and the SAME working-set bound that `plan` / |
| 5 | +//! `variants --index --enforce` admit, with no contract code re-derived. |
| 6 | +//! |
| 7 | +//! Run with: `cargo run --example columnkit_coverage` |
| 8 | +
|
| 9 | +use std::collections::BTreeMap; |
| 10 | +use std::io::{self, Write}; |
| 11 | +use std::sync::Arc; |
| 12 | + |
| 13 | +use rosalind::call::{run_bounded_whole_genome, ColumnAnalyzer}; |
| 14 | +use rosalind::core::{AlignedRead, CigarOp, CigarOpKind, Position, SamFlags}; |
| 15 | +use rosalind::genomics::{GenomeIndex, IndexReader, IndexWriter}; |
| 16 | +use rosalind::{PileupColumn, PileupParams, SliceSource}; |
| 17 | + |
| 18 | +/// Emit a per-locus coverage track (contig, 1-based pos, depth). This is the |
| 19 | +/// entire analyzer a builder writes — bounded memory, determinism, and a |
| 20 | +/// verifiable receipt come from the driver, not from here. |
| 21 | +struct CoverageTrack; |
| 22 | + |
| 23 | +impl ColumnAnalyzer for CoverageTrack { |
| 24 | + fn header(&self) -> Option<String> { |
| 25 | + Some("#contig\tpos\tdepth\n".to_string()) |
| 26 | + } |
| 27 | + |
| 28 | + fn params(&self) -> BTreeMap<String, String> { |
| 29 | + BTreeMap::from([("analyzer".to_string(), "coverage".to_string())]) |
| 30 | + } |
| 31 | + |
| 32 | + fn on_column( |
| 33 | + &mut self, |
| 34 | + col: &PileupColumn, |
| 35 | + contig: &str, |
| 36 | + out: &mut dyn Write, |
| 37 | + ) -> io::Result<()> { |
| 38 | + writeln!(out, "{contig}\t{}\t{}", col.locus.pos.0 + 1, col.depth()) |
| 39 | + } |
| 40 | +} |
| 41 | + |
| 42 | +fn read(pos: u32, seq: &[u8]) -> AlignedRead { |
| 43 | + AlignedRead { |
| 44 | + contig: 0, |
| 45 | + pos: Position(pos), |
| 46 | + mapq: 60, |
| 47 | + flags: SamFlags(0), |
| 48 | + cigar: vec![CigarOp::new(CigarOpKind::Match, seq.len() as u32)], |
| 49 | + seq: Arc::from(seq.to_vec().into_boxed_slice()), |
| 50 | + qual: Arc::from(vec![40u8; seq.len()].into_boxed_slice()), |
| 51 | + } |
| 52 | +} |
| 53 | + |
| 54 | +fn main() -> Result<(), Box<dyn std::error::Error>> { |
| 55 | + // A tiny in-memory genome + a few reads. A real run would pass a persisted |
| 56 | + // index and a coordinate-sorted BAM (`StreamingBamSource`) — the driver and |
| 57 | + // the contract are identical either way. |
| 58 | + let dir = std::env::temp_dir().join("rosalind-columnkit-example"); |
| 59 | + std::fs::create_dir_all(&dir)?; |
| 60 | + let idx_path = dir.join("ref.idx"); |
| 61 | + let index = |
| 62 | + GenomeIndex::from_named_sequences(&[("chr1".to_string(), b"ACGTACGTACGTACGT".to_vec())])?; |
| 63 | + IndexWriter::create(&idx_path)?.write_genome_index(&index)?; |
| 64 | + let loaded = IndexReader::open(&idx_path)?; |
| 65 | + let ref_view = loaded.reference_view()?; |
| 66 | + let contigs = loaded.contigs(); |
| 67 | + |
| 68 | + let reads = vec![ |
| 69 | + read(0, b"ACGTACGT"), |
| 70 | + read(0, b"ACGTACGT"), |
| 71 | + read(4, b"ACGTACGT"), |
| 72 | + ]; |
| 73 | + |
| 74 | + let mut analyzer = CoverageTrack; |
| 75 | + let mut out = io::stdout().lock(); |
| 76 | + let (ws, _skips) = run_bounded_whole_genome( |
| 77 | + &mut analyzer, |
| 78 | + SliceSource::new(reads), |
| 79 | + &ref_view, |
| 80 | + contigs, |
| 81 | + PileupParams::default(), |
| 82 | + &mut out, |
| 83 | + )?; |
| 84 | + out.flush()?; |
| 85 | + |
| 86 | + eprintln!( |
| 87 | + "\n# inherited a bounded working set of {} bytes — coverage-bounded, \ |
| 88 | + independent of input size; this is the value `plan`/`--enforce` admit.", |
| 89 | + ws.bytes |
| 90 | + ); |
| 91 | + |
| 92 | + std::fs::remove_dir_all(&dir).ok(); |
| 93 | + Ok(()) |
| 94 | +} |
0 commit comments