rsomics_align_score/
lib.rs

1use rsomics_align_core::{Op, ScoreParams, needleman_wunsch, smith_waterman};
2use rsomics_common::{Result, RsomicsError};
3use std::io::Write;
4use std::path::Path;
5
6pub enum AlignMode {
7    Global,
8    Local,
9}
10
11#[allow(clippy::cast_precision_loss)]
12fn identity(ops: &[Op]) -> f64 {
13    if ops.is_empty() {
14        return 0.0;
15    }
16    let matches = ops.iter().filter(|o| matches!(o, Op::Match)).count();
17    matches as f64 / ops.len() as f64
18}
19
20pub fn align_pairs(
21    fasta: &Path,
22    mode: &AlignMode,
23    params: &ScoreParams,
24    output: &mut dyn Write,
25) -> Result<u64> {
26    let mut reader = needletail::parse_fastx_file(fasta)
27        .map_err(|e| RsomicsError::InvalidInput(format!("{}: {e}", fasta.display())))?;
28
29    let mut seqs: Vec<(String, Vec<u8>)> = Vec::new();
30    while let Some(result) = reader.next() {
31        let record = result.map_err(|e| RsomicsError::InvalidInput(format!("read: {e}")))?;
32        let name = std::str::from_utf8(record.id())
33            .map_err(|e| RsomicsError::InvalidInput(format!("name: {e}")))?
34            .to_string();
35        seqs.push((name, record.seq().to_vec()));
36    }
37
38    let mut count = 0u64;
39    writeln!(output, "seq1\tseq2\tscore\tidentity").map_err(RsomicsError::Io)?;
40
41    for i in 0..seqs.len() {
42        for j in (i + 1)..seqs.len() {
43            let aln = match mode {
44                AlignMode::Global => needleman_wunsch(&seqs[i].1, &seqs[j].1, params),
45                AlignMode::Local => smith_waterman(&seqs[i].1, &seqs[j].1, params),
46            }
47            .map_err(|e| RsomicsError::InvalidInput(format!("align: {e}")))?;
48
49            writeln!(
50                output,
51                "{}\t{}\t{}\t{:.4}",
52                seqs[i].0,
53                seqs[j].0,
54                aln.score,
55                identity(&aln.ops)
56            )
57            .map_err(RsomicsError::Io)?;
58            count += 1;
59        }
60    }
61    Ok(count)
62}
rsomics_align_score/lib.rs

rsomics_align_score/
lib.rs
