seqtk_rs/
size.rs

1use crate::io_utils::{FaReader, FqReader, Output};
2use rayon::slice::ParallelSliceMut;
3
4/// Parses FASTQ file and computes read statistics.
5/// Outputs the results to [`std::io::stdout()`].
6///The output columns are:
7/// - `#seq`: Number of reads
8/// - `#bases`: Total number of bases
9/// - `avg_size`: Average read length
10/// - `mini_size`: Minimum read length
11/// - `med_size`: Median read length
12/// - `max_size`: Maximum read length
13/// - `N50`: N50 read length
14///
15/// # Arguments
16///
17/// * `path` - FASTQ path
18///
19/// # Errors
20///
21/// Return an error if the operation cannot be completed.
22pub fn calc_fq_size(path: &str) -> Result<(), std::io::Error> {
23    let fq_iter = FqReader::new(path)?;
24    let mut seq_len: Vec<usize> = Vec::new();
25    for record in fq_iter.records() {
26        match record {
27            Ok(read) => {
28                seq_len.push(read.seq().len());
29            }
30            Err(e) => eprintln!("Error read fASTQ: {}", e),
31        }
32    }
33    seq_len.par_sort_unstable();
34    let result = get_result_str(&seq_len);
35    let mut output = Output::new();
36
37    output.write(result)?;
38    Ok(())
39}
40/// Parses FASTA file and computes sequence statistics.
41/// Outputs the results to [`std::io::stdout()`].
42///The output columns are:
43/// - `#seq`: Number of sequences
44/// - `#bases`: Total number of bases
45/// - `avg_size`: Average sequence length
46/// - `mini_size`: Minimum sequence length
47/// - `med_size`: Median sequence length
48/// - `max_size`: Maximum sequence length
49/// - `N50`: N50 sequence length
50///
51/// # Arguments
52///
53/// * `path` - FASTA path
54///
55/// # Errors
56///
57/// Return an error if the operation cannot be completed.
58pub fn calc_fa_size(path: &str) -> Result<(), std::io::Error> {
59    let fa_iter = FaReader::new(path)?;
60    let mut seq_len: Vec<usize> = Vec::new();
61    for record in fa_iter.records() {
62        match record {
63            Ok(read) => {
64                seq_len.push(read.seq().len());
65            }
66            Err(e) => eprintln!("Error read fASTA: {}", e),
67        }
68    }
69    seq_len.par_sort_unstable();
70    let result = get_result_str(&seq_len);
71    let mut output = Output::new();
72    output.write(result)?;
73    Ok(())
74}
75
76fn get_result_str(sorted_seq_len: &[usize]) -> String {
77    // #seq, #bases, avg_size, min_size, med_size, max_size, N50
78    let sum: usize = sorted_seq_len.iter().sum();
79    let size = sorted_seq_len.len();
80    let median = match size {
81        0 => f64::NAN,
82        _ => {
83            let mid = size / 2;
84            match size % 2 {
85                1 => sorted_seq_len[mid] as f64,
86                _ => (sorted_seq_len[mid - 1] + sorted_seq_len[mid]) as f64 / 2.0,
87            }
88        }
89    };
90    let mut n50: usize = 0;
91    let half: usize = sum / 2;
92    let mut acc: usize = 0;
93    for &cur in sorted_seq_len.iter().rev() {
94        acc += cur;
95        if acc >= half {
96            n50 = cur;
97            break;
98        }
99    }
100    let result = format!(
101        "{}\t{}\t{:2}\t{}\t{}\t{}\t{}\n",
102        sorted_seq_len.len(),
103        sum,
104        sum as f64 / size as f64,
105        sorted_seq_len.first().unwrap(),
106        median,
107        sorted_seq_len.last().unwrap(),
108        n50
109    );
110    result
111}
seqtk_rs/size.rs

seqtk_rs/
size.rs
