deepbiop-bam 0.1.11

Deep Learning Processing Library for Bam Format
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

use anyhow::Result;
use noodles::{bam, bgzf};
use rayon::prelude::*;
use std::{fs::File, num::NonZeroUsize, path::Path, thread};

use noodles::fastq;

// FIXME: The function has a bug since seq != qual

pub fn bam2fq(bam: &Path, threads: Option<usize>) -> Result<Vec<fastq::Record>> {
    let file = File::open(bam)?;
    let worker_count = if let Some(threads) = threads {
        NonZeroUsize::new(threads)
            .unwrap()
            .min(thread::available_parallelism().unwrap_or(NonZeroUsize::MIN))
    } else {
        thread::available_parallelism().unwrap_or(NonZeroUsize::MIN)
    };

    let decoder = bgzf::MultithreadedReader::with_worker_count(worker_count, file);
    let mut reader = bam::io::Reader::from(decoder);
    let _header = reader.read_header()?;

    reader
        .records()
        .par_bridge()
        .map(|result| {
            let record = result.unwrap();

            let seq = record.sequence().as_ref().to_vec();
            let qual = record.quality_scores().as_ref().to_vec();

            if seq.len() != qual.len() {
                let name = String::from_utf8_lossy(record.name().unwrap().as_ref()).to_string();
                return Err(anyhow::anyhow!(
                    "{} seq and qual length are not equal",
                    name
                ));
            }

            let fq_record = fastq::Record::new(
                fastq::record::Definition::new(record.name().unwrap().to_vec(), ""),
                seq,
                qual,
            );
            Ok(fq_record)
        })
        .collect::<Result<Vec<fastq::Record>>>()
}