gapsmith_align/
diamond.rs

1//! Diamond aligner (protein-protein).
2//!
3//! Command template mirrors `src/gapseq_find.sh:592–600`:
4//!
5//! ```text
6//! diamond makedb --in TARGET -d orgdb
7//! diamond blastp -d orgdb.dmnd -q QUERY --threads N --out alignments.tsv
8//!                --outfmt 6 qseqid pident evalue bitscore qcovhsp stitle sstart send
9//!                --query-cover 75
10//! ```
11//!
12//! By default we pass `--more-sensitive`, matching gapseq's default `aliArgs`
13//! for diamond.
14
15use crate::blast::which;
16use crate::error::{io_err, AlignError};
17use crate::hit::Hit;
18use crate::tsv::parse_tsv;
19use crate::{AlignOpts, Aligner};
20use std::fs::File;
21use std::io::BufReader;
22use std::path::Path;
23use std::process::Command;
24
25/// Columns emitted by our diamond invocation. Note `qcovhsp` (per-HSP) rather
26/// than `qcovs` — diamond doesn't support `qcovs`. Both are 0–100.
27const COLUMNS: &[&str] = &[
28    "qseqid", "pident", "evalue", "bitscore", "qcovhsp", "stitle", "sstart", "send",
29];
30
31pub struct DiamondAligner {
32    /// Pass `--more-sensitive` (default: true). Matches gapseq's aliArgs
33    /// default when the user doesn't override `-R`.
34    pub more_sensitive: bool,
35}
36
37impl DiamondAligner {
38    pub fn new() -> Self {
39        Self { more_sensitive: true }
40    }
41}
42impl Default for DiamondAligner {
43    fn default() -> Self {
44        Self::new()
45    }
46}
47
48impl Aligner for DiamondAligner {
49    fn name(&self) -> &'static str {
50        "diamond"
51    }
52
53    fn align(
54        &self,
55        query_fasta: &Path,
56        target_fasta: &Path,
57        opts: &AlignOpts,
58    ) -> Result<Vec<Hit>, AlignError> {
59        require("diamond")?;
60        let tmp = tempfile::tempdir().map_err(|e| io_err(Path::new("/tmp"), e))?;
61        let db_prefix = tmp.path().join("orgdb");
62        let out_tsv = tmp.path().join("aln.tsv");
63
64        // makedb
65        let out = Command::new("diamond")
66            .arg("makedb")
67            .arg("--in")
68            .arg(target_fasta)
69            .arg("-d")
70            .arg(&db_prefix)
71            .arg("--quiet")
72            .output()
73            .map_err(|e| io_err(Path::new("diamond"), e))?;
74        if !out.status.success() {
75            return Err(AlignError::ToolFailed {
76                tool: "diamond",
77                status: out.status,
78                stderr: String::from_utf8_lossy(&out.stderr).to_string(),
79            });
80        }
81
82        // blastp
83        let dmnd_db = {
84            let mut p = db_prefix.clone();
85            p.set_extension("dmnd");
86            p
87        };
88        let mut cmd = Command::new("diamond");
89        cmd.arg("blastp")
90            .arg("-d")
91            .arg(&dmnd_db)
92            .arg("-q")
93            .arg(query_fasta)
94            .arg("--threads")
95            .arg(opts.threads.to_string())
96            .arg("--out")
97            .arg(&out_tsv)
98            .arg("--outfmt")
99            .arg("6");
100        for c in COLUMNS {
101            cmd.arg(c);
102        }
103        cmd.arg("--query-cover").arg(opts.coverage_pct.to_string());
104        if self.more_sensitive {
105            cmd.arg("--more-sensitive");
106        }
107        if let Some(e) = opts.evalue {
108            cmd.arg("--evalue").arg(e.to_string());
109        }
110        if opts.quiet {
111            cmd.arg("--quiet");
112        }
113        for a in &opts.extra_args {
114            cmd.arg(a);
115        }
116
117        tracing::debug!(?cmd, "running diamond");
118        let out = cmd
119            .output()
120            .map_err(|e| io_err(Path::new("diamond"), e))?;
121        if !out.status.success() {
122            return Err(AlignError::ToolFailed {
123                tool: "diamond",
124                status: out.status,
125                stderr: String::from_utf8_lossy(&out.stderr).to_string(),
126            });
127        }
128        if !opts.quiet && !out.stderr.is_empty() {
129            eprintln!("{}", String::from_utf8_lossy(&out.stderr));
130        }
131
132        let f = File::open(&out_tsv).map_err(|e| io_err(&out_tsv, e))?;
133        parse_tsv(BufReader::new(f), false)
134    }
135}
136
137fn require(tool: &'static str) -> Result<(), AlignError> {
138    if which(tool).is_some() {
139        Ok(())
140    } else {
141        Err(AlignError::ToolMissing { tool })
142    }
143}
gapsmith_align/diamond.rs

gapsmith_align/
diamond.rs
