seqtk-rs 0.1.0

Pure rust version of seqtk
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
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111

use crate::io_utils::{FaReader, FqReader, Output};
use rayon::slice::ParallelSliceMut;

/// Parses FASTQ file and computes read statistics.
/// Outputs the results to [`std::io::stdout()`].
///The output columns are:
/// - `#seq`: Number of reads
/// - `#bases`: Total number of bases
/// - `avg_size`: Average read length
/// - `mini_size`: Minimum read length
/// - `med_size`: Median read length
/// - `max_size`: Maximum read length
/// - `N50`: N50 read length
///
/// # Arguments
///
/// * `path` - FASTQ path
///
/// # Errors
///
/// Return an error if the operation cannot be completed.
pub fn calc_fq_size(path: &str) -> Result<(), std::io::Error> {
    let fq_iter = FqReader::new(path)?;
    let mut seq_len: Vec<usize> = Vec::new();
    for record in fq_iter.records() {
        match record {
            Ok(read) => {
                seq_len.push(read.seq().len());
            }
            Err(e) => eprintln!("Error read fASTQ: {}", e),
        }
    }
    seq_len.par_sort_unstable();
    let result = get_result_str(&seq_len);
    let mut output = Output::new();

    output.write(result)?;
    Ok(())
}
/// Parses FASTA file and computes sequence statistics.
/// Outputs the results to [`std::io::stdout()`].
///The output columns are:
/// - `#seq`: Number of sequences
/// - `#bases`: Total number of bases
/// - `avg_size`: Average sequence length
/// - `mini_size`: Minimum sequence length
/// - `med_size`: Median sequence length
/// - `max_size`: Maximum sequence length
/// - `N50`: N50 sequence length
///
/// # Arguments
///
/// * `path` - FASTA path
///
/// # Errors
///
/// Return an error if the operation cannot be completed.
pub fn calc_fa_size(path: &str) -> Result<(), std::io::Error> {
    let fa_iter = FaReader::new(path)?;
    let mut seq_len: Vec<usize> = Vec::new();
    for record in fa_iter.records() {
        match record {
            Ok(read) => {
                seq_len.push(read.seq().len());
            }
            Err(e) => eprintln!("Error read fASTA: {}", e),
        }
    }
    seq_len.par_sort_unstable();
    let result = get_result_str(&seq_len);
    let mut output = Output::new();
    output.write(result)?;
    Ok(())
}

fn get_result_str(sorted_seq_len: &[usize]) -> String {
    // #seq, #bases, avg_size, min_size, med_size, max_size, N50
    let sum: usize = sorted_seq_len.iter().sum();
    let size = sorted_seq_len.len();
    let median = match size {
        0 => f64::NAN,
        _ => {
            let mid = size / 2;
            match size % 2 {
                1 => sorted_seq_len[mid] as f64,
                _ => (sorted_seq_len[mid - 1] + sorted_seq_len[mid]) as f64 / 2.0,
            }
        }
    };
    let mut n50: usize = 0;
    let half: usize = sum / 2;
    let mut acc: usize = 0;
    for &cur in sorted_seq_len.iter().rev() {
        acc += cur;
        if acc >= half {
            n50 = cur;
            break;
        }
    }
    let result = format!(
        "{}\t{}\t{:2}\t{}\t{}\t{}\t{}\n",
        sorted_seq_len.len(),
        sum,
        sum as f64 / size as f64,
        sorted_seq_len.first().unwrap(),
        median,
        sorted_seq_len.last().unwrap(),
        n50
    );
    result
}