fibertools-rs 0.8.2

Fiber-seq toolkit in rust
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

use fibertools_rs::utils::bio_io;
use fibertools_rs::utils::input_bam::FiberFilters;
use rust_htslib::bam::Read;

fn get_fiber_data_from_test_bam(bam_file: &str) -> Vec<fibertools_rs::fiber::FiberseqData> {
    let mut bam = bio_io::bam_reader(bam_file);
    bam.records()
        .map(|r| {
            let record = r.unwrap();

            fibertools_rs::fiber::FiberseqData::new(record, None, &FiberFilters::default())
        })
        .collect::<Vec<_>>()
}

#[test]
fn test_msp_extract() -> Result<(), Box<dyn std::error::Error>> {
    let expected_msp_starts: Vec<Option<i64>> = vec![
        141, 376, 789, 977, 1112, 1316, 1453, 1691, 1974, 2134, 2283, 2455, 2677, 2784, 2920, 3057,
        3205, 3508, 3729, 4157, 4318, 4583, 4891, 5109, 5322, 5568, 5719, 5901, 6066, 6217, 6383,
        6583, 6727, 6921, 7090, 7276,
    ]
    .into_iter()
    .map(Some)
    .collect();
    let fiber_records = get_fiber_data_from_test_bam("tests/data/msp_nuc.bam");
    assert!(fiber_records.len() == 1);
    let fiber_data = &fiber_records[0];
    assert_eq!(fiber_data.msp.option_starts(), expected_msp_starts);
    Ok(())
}

#[test]
fn test_many_msps() {
    let fiber_records = get_fiber_data_from_test_bam("tests/data/all.bam");
    for fiber_data in fiber_records {
        let m6a_starts = fiber_data.m6a.starts();
        let m6a = m6a_starts.iter().collect::<Vec<_>>();
        if m6a.is_empty() {
            continue;
        }

        let msps = fiber_data
            .msp
            .starts()
            .iter()
            .copied()
            .chain(fiber_data.msp.ends().iter().map(|&x| x - 1))
            .collect::<Vec<_>>();
        eprintln!("m6a: {m6a:?}");
        eprintln!("msp: {msps:?}");
        /*
        eprintln!(
            "strand: {:?}\t{:?}",
            fiber_data.record.strand(),
            fiber_data.record.seq_len()
        );
         */
        for msp in msps {
            assert!(
                m6a.contains(&&msp),
                "{}\nreverse:{}",
                msp,
                fiber_data.record.is_reverse()
            );
        }
    }
}