fibertools-rs 0.12.0

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

use crate::cli::FiberHmmOptions;
use crate::utils::basemods::M6A_TYPE;
use crate::*;
use anyhow::Error;
use bio::alphabets::dna::revcomp;

pub fn run_fiber_hmm(opts: &mut FiberHmmOptions) -> Result<(), Error> {
    let mut bam = opts.input.bam_reader();
    let mut out = opts.input.bam_writer(&opts.out);
    for fiber in opts.input.fibers(&mut bam) {
        // m6a positions on the forward (molecular) strand
        let _m6a: Vec<u32> = fiber
            .annotations
            .get_forward_coords(M6A_TYPE)
            .map(|v| v.into_iter().map(|(s, _)| s).collect())
            .unwrap_or_default();
        // forward strand sequence
        let mut seq = fiber.record.seq().as_bytes();
        if fiber.record.is_reverse() {
            seq = revcomp(seq);
        }
        log::trace!("Run HMM on fiber: {seq:?}");

        // TODO hmm stuff on the one read

        out.write(&fiber.record).unwrap();
    }

    Ok(())
}