seqair 0.1.0

Pure-Rust BAM/SAM/CRAM/FASTA reader and pileup engine
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70

#![allow(unused, reason = "example code")]

use clap::Parser as _;
use seqair::{
    Readers,
    bam::{
        RecordStore,
        record_store::{CustomizeRecordStore, SlimRecord},
    },
    reader::SegmentOptions,
};

#[derive(Debug, clap::Parser)]
struct Cli {
    input: std::path::PathBuf,     // BAM file
    reference: std::path::PathBuf, // FASTA file
    #[clap(long)]
    region: seqair_types::RegionString,
    #[clap(long, default_value_t = 20)]
    min_mapq: u8,
}

fn main() -> anyhow::Result<()> {
    let args = Cli::parse();

    let mut readers = Readers::open_customized(
        &args.input,
        &args.reference,
        ReadInfoBuilder { min_mapq: args.min_mapq },
    )?;
    let max_len = std::num::NonZeroU32::new(100_000).expect("non-zero literal");
    let opts = SegmentOptions::new(max_len);
    let plan: Vec<_> = readers.segments(&args.region, opts)?.collect();

    for segment in &plan {
        let mut engine = readers.pileup(segment)?;

        while let Some(column) = engine.pileups() {
            for aln in column.alignments() {
                let info = aln.extra();

                // ...
            }
        }
        // `engine` drops here; the RecordStore goes back into `readers`.
    }

    Ok(())
}

struct ReadInfo {
    read_group: Option<seqair_types::SmolStr>,
}

#[derive(Debug, Clone)]
struct ReadInfoBuilder {
    min_mapq: u8,
}

impl CustomizeRecordStore for ReadInfoBuilder {
    type Extra = ReadInfo;

    fn filter(&mut self, rec: &SlimRecord, _: &RecordStore<ReadInfo>) -> bool {
        !rec.flags.is_unmapped()
    }

    fn compute(&mut self, rec: &SlimRecord, store: &RecordStore<ReadInfo>) -> ReadInfo {
        ReadInfo { read_group: rec.aux(store).ok().and_then(|aux| aux.get("RG").ok()) }
    }
}