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

//! How many allocations one slice of a CRAM costs, and where.
//!
//! A counting global allocator, a single-threaded read of one region, and the
//! totals divided by the records that came back. Sampling profilers attribute
//! `malloc` to whoever is on the stack, which for a tree this deep is nobody
//! useful; this counts.
use std::alloc::{GlobalAlloc, Layout, System};
use std::sync::atomic::{AtomicU64, Ordering};

static ALLOCS: AtomicU64 = AtomicU64::new(0);
static BYTES: AtomicU64 = AtomicU64::new(0);

struct Counting;

unsafe impl GlobalAlloc for Counting {
    unsafe fn alloc(&self, layout: Layout) -> *mut u8 {
        ALLOCS.fetch_add(1, Ordering::Relaxed);
        BYTES.fetch_add(layout.size() as u64, Ordering::Relaxed);
        System.alloc(layout)
    }
    unsafe fn dealloc(&self, ptr: *mut u8, layout: Layout) {
        System.dealloc(ptr, layout)
    }
    unsafe fn realloc(&self, ptr: *mut u8, layout: Layout, new: usize) -> *mut u8 {
        ALLOCS.fetch_add(1, Ordering::Relaxed);
        BYTES.fetch_add(new.saturating_sub(layout.size()) as u64, Ordering::Relaxed);
        System.realloc(ptr, layout, new)
    }
}

#[global_allocator]
static GLOBAL: Counting = Counting;

fn snapshot() -> (u64, u64) {
    (
        ALLOCS.load(Ordering::Relaxed),
        BYTES.load(Ordering::Relaxed),
    )
}

fn main() {
    let dir = "local/test_data";
    let args: Vec<String> = std::env::args().collect();
    let file = args
        .get(1)
        .cloned()
        .unwrap_or(format!("{dir}/AtTPax7_H3K4me.10M.cram"));
    let reference = format!("{dir}/mm10.fa.gz");
    let with_reference = args.get(2).map(|s| s != "noref").unwrap_or(true);

    let before = snapshot();
    let reader = gwseq_io::cram::CramReader::open(
        &file,
        None,
        with_reference.then_some(reference.as_str()),
        1,
        None,
        None,
    )
    .expect("opens");
    let after = snapshot();
    println!(
        "open                {:>10} allocations {:>12} bytes",
        after.0 - before.0,
        after.1 - before.1
    );

    // A region wide enough that the records returned are most of the records
    // decoded, so dividing by them means something. A narrow locus still
    // decodes a whole slice and would flatter or damn the count at random.
    let chr = reader.chr_sizes().iter().next().map(|e| e.id.clone());
    let chr = chr.unwrap_or_else(|| "chr1".to_string());
    let locs = gwseq_io::genomic::Locs::spans(&[chr], &[3_000_000], &[6_000_000]).expect("locs");
    let request = gwseq_io::bam::EntriesRequest::new(locs).filter(false);
    let before = snapshot();
    let out = reader.read_entries(&request).expect("reads");
    let after = snapshot();
    let records: usize = out.iter().map(|v| v.len()).sum();
    let allocs = after.0 - before.0;
    println!(
        "read                {:>10} allocations {:>12} bytes   ({records} records)",
        allocs,
        after.1 - before.1
    );
    if records > 0 {
        println!(
            "                    {:>10.2} allocations per record",
            allocs as f64 / records as f64
        );
    }
}