fastq_rs 0.0.2

Multi purpose fastq toolkit.
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

use crate::common::utils::error_to_phred;
use crate::common::{bio_fastq_reader, mean_error_and_phred, mean_len, write_json};

use anyhow::Result;
use log::error;
use rayon::prelude::*;
use serde::{Deserialize, Serialize};
use std::path::PathBuf;
use std::sync::atomic::AtomicUsize;
use std::sync::atomic::Ordering::Relaxed;
use std::sync::{Arc, Mutex};

#[derive(Serialize, Deserialize)]
pub struct FastqStats {
    pub num_reads: usize,
    pub num_bases: usize,
    pub mean_error: f64,
    pub mean_phred: u8,
    pub mean_len: usize,
    pub shortest: Option<Vec<usize>>,
    pub longest: Option<Vec<usize>>,
}

/// # Arguments
/// * `fastq` - Path to FASTQ file (optional, defaults to stdin).
/// * `outfile` - Where to write results (optional, defaults to stdout).
///
/// # Returns
/// * `Ok(FastqStats)` if successful.
/// * `Err` if not.
pub fn fastq_stats(fastq: Option<PathBuf>, outfile: Option<PathBuf>) -> Result<FastqStats> {
    let reader = bio_fastq_reader(fastq)?;

    // Initialize thread safe variables.
    let num_reads = AtomicUsize::new(0);
    let num_bases = AtomicUsize::new(0);
    let mean_errors: Arc<Mutex<Vec<f64>>> = Arc::new(Mutex::new(Vec::new()));
    let read_lengths: Arc<Mutex<Vec<usize>>> = Arc::new(Mutex::new(Vec::new()));

    reader.records().par_bridge().for_each(|record| {
        let record = match record {
            Ok(record) => record,
            Err(e) => {
                error!("{:?}", e);
                return;
            }
        };

        let record_len: usize = record.seq().len();

        // Read error.
        let (mean_error, _) = mean_error_and_phred(record.qual());

        num_reads.fetch_add(1, Relaxed);
        num_bases.fetch_add(record_len, Relaxed);

        // Aggregate mean error and read length per read.
        mean_errors.lock().unwrap().push(mean_error);
        read_lengths.lock().unwrap().push(record_len);
    });

    // NOTE that for performance reasons, we calculate the mean of the mean read error rates.
    // To get the true mean error, we'd have to store every single nucleotide error rate, sum
    // them up and divide by the total number of bases (unfeasible for large files).
    let mean_errors = Arc::try_unwrap(mean_errors).unwrap().into_inner().unwrap();
    let mean_mean_error = mean_errors.iter().sum::<f64>() / mean_errors.len() as f64;
    let mean_mean_phred = error_to_phred(mean_mean_error);

    // Mean read length.
    let mut read_lengths = Arc::try_unwrap(read_lengths).unwrap().into_inner().unwrap();
    read_lengths.par_sort();

    let mean_len = mean_len(read_lengths.as_slice());

    let fastq_stats = FastqStats {
        num_reads: num_reads.into_inner(),
        num_bases: num_bases.into_inner(),
        mean_error: mean_mean_error,
        mean_phred: mean_mean_phred,
        mean_len,
        shortest: read_lengths.first_chunk::<5>().map(|c| c.to_vec()),
        longest: read_lengths.last_chunk::<5>().map(|c| c.to_vec()),
    };

    // Write json to output file.
    write_json(outfile, &fastq_stats)?;

    Ok(fastq_stats)
}