gwseq-io 0.2.1

Rust library for processing bigWig, bigBed, BAM, CRAM and HiC files
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
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95

//! 500 random 1 kb loci, BAM against CRAM, as the review measured them.
use std::time::Instant;
fn main() {
    let args: Vec<String> = std::env::args().collect();
    let dir = "local/test_data";
    let reference = format!("{dir}/mm10.fa.gz");
    let n: usize = args.get(1).map(|s| s.parse().unwrap()).unwrap_or(500);
    let parallel: i64 = args.get(2).map(|s| s.parse().unwrap()).unwrap_or(4);

    // A deterministic spread of loci over chr1..chr19.
    let mut seed = 0x1234_5678_9abc_def0u64;
    let mut next = || {
        seed ^= seed >> 12;
        seed ^= seed << 25;
        seed ^= seed >> 27;
        seed.wrapping_mul(0x2545_F491_4F6C_DD1D)
    };
    let mut chrs = Vec::new();
    let mut starts = Vec::new();
    let mut ends = Vec::new();
    for _ in 0..n {
        let c = (next() % 19) + 1;
        let s = (next() % 60_000_000) as i64 + 3_000_000;
        chrs.push(format!("chr{c}"));
        starts.push(s);
        ends.push(s + 1000);
    }
    let locs = gwseq_io::genomic::Locs::spans(&chrs, &starts, &ends).unwrap();

    let run = |label: &str, open: &dyn Fn() -> gwseq_io::align::Alignments| {
        let reader = open();
        let request = gwseq_io::bam::EntriesRequest::new(locs.clone());
        let t = Instant::now();
        let out = reader.read_entries(&request).unwrap();
        let total: usize = out.iter().map(|v| v.len()).sum();
        println!(
            "{label:<34} {:>7.2} s   {total} records",
            t.elapsed().as_secs_f64()
        );
    };

    run("BAM (.bai)", &|| {
        gwseq_io::align::Alignments::Bam(
            gwseq_io::bam::BamReader::open(
                &format!("{dir}/AtTPax7_H3K4me.10M.bam"),
                None,
                parallel,
                None,
                None,
            )
            .unwrap(),
        )
    });
    run("CRAM (.crai, reference)", &|| {
        gwseq_io::align::Alignments::Cram(
            gwseq_io::cram::CramReader::open(
                &format!("{dir}/AtTPax7_H3K4me.10M.cram"),
                None,
                Some(&reference),
                parallel,
                None,
                None,
            )
            .unwrap(),
        )
    });
    run("CRAM (.crai, no reference)", &|| {
        gwseq_io::align::Alignments::Cram(
            gwseq_io::cram::CramReader::open(
                &format!("{dir}/AtTPax7_H3K4me.10M.cram"),
                None,
                Some("/dev/null"),
                parallel,
                None,
                None,
            )
            .unwrap(),
        )
    });
    // A named index that is not there falls back to building one from the
    // container headers, which is the path a file with no `.crai` takes.
    run("CRAM (built index, reference)", &|| {
        gwseq_io::align::Alignments::Cram(
            gwseq_io::cram::CramReader::open(
                &format!("{dir}/AtTPax7_H3K4me.10M.cram"),
                Some("/nonexistent.crai"),
                Some(&reference),
                parallel,
                None,
                None,
            )
            .unwrap(),
        )
    });
}