mafft_core/
external.rs

1/// External tool integration for RNA and structure-based alignment modes.
2///
3/// Q-INS-i: McCaskill base-pair probabilities via `mxscarnamod`
4/// X-INS-i: CONTRAfold structure predictions via `contrafold`
5/// SCARNA-like: Structural homology via `dash_client`
6///
7/// These modes follow the same pattern as C MAFFT: call external tools
8/// as subprocesses, parse their output, and convert to constraint tables
9/// for use in the alignment DP.
10
11use std::io::Write;
12use std::path::{Path, PathBuf};
13use std::process::Command;
14
15/// Base-pair probability: position i pairs with position j with given probability.
16#[derive(Debug, Clone)]
17pub struct BasePairProb {
18    pub left: usize,
19    pub right: usize,
20    pub prob: f64,
21}
22
23/// Per-sequence base-pair probability table.
24#[derive(Debug, Clone)]
25pub struct SequenceBpp {
26    pub seq_index: usize,
27    pub pairs: Vec<BasePairProb>,
28}
29
30/// Find an external tool in PATH or MAFFT_BINARIES directory.
31pub fn find_tool(name: &str) -> Option<PathBuf> {
32    // Check MAFFT_BINARIES environment variable first
33    if let Ok(bindir) = std::env::var("MAFFT_BINARIES") {
34        let path = Path::new(&bindir).join(name);
35        if path.exists() {
36            return Some(path);
37        }
38    }
39
40    // Check PATH
41    if let Ok(output) = Command::new("which").arg(name).output() {
42        if output.status.success() {
43            let path = String::from_utf8_lossy(&output.stdout).trim().to_string();
44            if !path.is_empty() {
45                return Some(PathBuf::from(path));
46            }
47        }
48    }
49
50    None
51}
52
53/// Run McCaskill base-pair probability prediction on a single sequence.
54///
55/// Calls: `mxscarnamod -m -writebpp`
56/// Input: FASTA (ungapped)
57/// Output: `left right probability` per line
58pub fn run_mccaskill(sequence: &[u8], tool_path: &Path) -> Result<Vec<BasePairProb>, String> {
59    let tmpdir = std::env::temp_dir();
60    let infile = tmpdir.join("_mafftrs_mccaskillin");
61    let outfile = tmpdir.join("_mafftrs_mccaskillout");
62
63    // Write input FASTA
64    let mut f = std::fs::File::create(&infile)
65        .map_err(|e| format!("Cannot create temp file: {e}"))?;
66    writeln!(f, ">seq").map_err(|e| format!("Write error: {e}"))?;
67    f.write_all(sequence).map_err(|e| format!("Write error: {e}"))?;
68    writeln!(f).map_err(|e| format!("Write error: {e}"))?;
69    drop(f);
70
71    // Run mxscarnamod
72    let output = Command::new(tool_path)
73        .args(["-m", "-writebpp"])
74        .stdin(std::fs::File::open(&infile).map_err(|e| format!("Cannot open temp: {e}"))?)
75        .output()
76        .map_err(|e| format!("Failed to run mxscarnamod: {e}"))?;
77
78    if !output.status.success() {
79        return Err(format!("mxscarnamod failed: {}", String::from_utf8_lossy(&output.stderr)));
80    }
81
82    // Write output to file for parsing
83    std::fs::write(&outfile, &output.stdout)
84        .map_err(|e| format!("Cannot write output: {e}"))?;
85
86    // Parse output: "left right probability" per line
87    let content = String::from_utf8_lossy(&output.stdout);
88    let pairs = parse_mccaskill_output(&content);
89
90    // Cleanup
91    let _ = std::fs::remove_file(&infile);
92    let _ = std::fs::remove_file(&outfile);
93
94    Ok(pairs)
95}
96
97/// Parse McCaskill output format: `left right probability` per line.
98fn parse_mccaskill_output(content: &str) -> Vec<BasePairProb> {
99    let mut pairs = Vec::new();
100    for line in content.lines() {
101        let parts: Vec<&str> = line.split_whitespace().collect();
102        if parts.len() >= 3 {
103            if let (Ok(left), Ok(right), Ok(prob)) = (
104                parts[0].parse::<usize>(),
105                parts[1].parse::<usize>(),
106                parts[2].parse::<f64>(),
107            ) {
108                if prob >= 0.01 {
109                    pairs.push(BasePairProb { left, right, prob });
110                }
111            }
112        }
113    }
114    pairs
115}
116
117/// Run CONTRAfold structure prediction on a single sequence.
118///
119/// Calls: `contrafold predict <infile> --posteriors 0.01 <outfile>`
120/// Input: FASTA (single sequence)
121/// Output: `pos pair1:prob pair2:prob ...` per line (1-indexed)
122pub fn run_contrafold(sequence: &[u8], tool_path: &Path) -> Result<Vec<BasePairProb>, String> {
123    let tmpdir = std::env::temp_dir();
124    let infile = tmpdir.join("_mafftrs_contrafoldin");
125    let outfile = tmpdir.join("_mafftrs_contrafoldout");
126
127    // Write input FASTA
128    let mut f = std::fs::File::create(&infile)
129        .map_err(|e| format!("Cannot create temp file: {e}"))?;
130    writeln!(f, ">seq").map_err(|e| format!("Write error: {e}"))?;
131    f.write_all(sequence).map_err(|e| format!("Write error: {e}"))?;
132    writeln!(f).map_err(|e| format!("Write error: {e}"))?;
133    drop(f);
134
135    // Run contrafold
136    let status = Command::new(tool_path)
137        .args(["predict", infile.to_str().unwrap(), "--posteriors", "0.01", outfile.to_str().unwrap()])
138        .status()
139        .map_err(|e| format!("Failed to run contrafold: {e}"))?;
140
141    if !status.success() {
142        return Err("contrafold failed".to_string());
143    }
144
145    // Parse output
146    let content = std::fs::read_to_string(&outfile)
147        .map_err(|e| format!("Cannot read contrafold output: {e}"))?;
148    let pairs = parse_contrafold_output(&content);
149
150    // Cleanup
151    let _ = std::fs::remove_file(&infile);
152    let _ = std::fs::remove_file(&outfile);
153
154    Ok(pairs)
155}
156
157/// Parse CONTRAfold output: `pos pair1:prob pair2:prob ...` (1-indexed).
158fn parse_contrafold_output(content: &str) -> Vec<BasePairProb> {
159    let mut pairs = Vec::new();
160    for line in content.lines() {
161        let parts: Vec<&str> = line.split_whitespace().collect();
162        if parts.is_empty() { continue; }
163        if let Ok(left) = parts[0].parse::<usize>() {
164            let left = left.saturating_sub(1); // Convert 1-indexed to 0-indexed
165            for &part in &parts[1..] {
166                if let Some(colon) = part.find(':') {
167                    if let (Ok(right), Ok(prob)) = (
168                        part[..colon].parse::<usize>(),
169                        part[colon + 1..].parse::<f64>(),
170                    ) {
171                        let right = right.saturating_sub(1);
172                        if prob >= 0.01 {
173                            pairs.push(BasePairProb { left, right, prob });
174                        }
175                    }
176                }
177            }
178        }
179    }
180    pairs
181}
182
183/// Run DASH structural alignment client.
184///
185/// Calls: `dash_client -url <server> -i <infile> -hat3 <outfile>`
186/// Input: FASTA (ungapped)
187/// Output: hat3 format constraint table
188pub fn run_dash(
189    sequences: &[Vec<u8>],
190    names: &[String],
191    server_url: &str,
192) -> Result<Vec<(usize, usize, f64, usize, usize, usize, usize)>, String> {
193    let tool_path = find_tool("dash_client")
194        .ok_or_else(|| "dash_client not found in PATH or MAFFT_BINARIES".to_string())?;
195
196    let tmpdir = std::env::temp_dir();
197    let infile = tmpdir.join("_mafftrs_dashin");
198    let outfile = tmpdir.join("_mafftrs_hat3seed");
199
200    // Write input FASTA (ungapped)
201    let mut f = std::fs::File::create(&infile)
202        .map_err(|e| format!("Cannot create temp file: {e}"))?;
203    for (name, seq) in names.iter().zip(sequences.iter()) {
204        writeln!(f, ">{name}").map_err(|e| format!("Write error: {e}"))?;
205        let ungapped: Vec<u8> = seq.iter().filter(|&&c| c != b'-').copied().collect();
206        f.write_all(&ungapped).map_err(|e| format!("Write error: {e}"))?;
207        writeln!(f).map_err(|e| format!("Write error: {e}"))?;
208    }
209    drop(f);
210
211    // Run dash_client
212    let status = Command::new(&tool_path)
213        .args([
214            "-url", server_url,
215            "-i", infile.to_str().unwrap(),
216            "-hat3", outfile.to_str().unwrap(),
217        ])
218        .status()
219        .map_err(|e| format!("Failed to run dash_client: {e}"))?;
220
221    if !status.success() {
222        return Err("dash_client failed".to_string());
223    }
224
225    // Parse hat3 output
226    let content = std::fs::read_to_string(&outfile)
227        .map_err(|e| format!("Cannot read DASH output: {e}"))?;
228    let constraints = parse_hat3(&content);
229
230    // Cleanup
231    let _ = std::fs::remove_file(&infile);
232    let _ = std::fs::remove_file(&outfile);
233
234    Ok(constraints)
235}
236
237/// Parse hat3 format: `i j overlapaa opt start1 end1 start2 end2`
238fn parse_hat3(content: &str) -> Vec<(usize, usize, f64, usize, usize, usize, usize)> {
239    let mut constraints = Vec::new();
240    for line in content.lines() {
241        let parts: Vec<&str> = line.split_whitespace().collect();
242        if parts.len() >= 8 {
243            if let (Ok(i), Ok(j), Ok(opt), Ok(s1), Ok(e1), Ok(s2), Ok(e2)) = (
244                parts[0].parse::<usize>(),
245                parts[1].parse::<usize>(),
246                parts[3].parse::<f64>(),
247                parts[4].parse::<usize>(),
248                parts[5].parse::<usize>(),
249                parts[6].parse::<usize>(),
250                parts[7].parse::<usize>(),
251            ) {
252                constraints.push((i, j, opt, s1, e1, s2, e2));
253            }
254        }
255    }
256    constraints
257}
258
259/// Compute base-pair probabilities for all sequences using McCaskill (Q-INS-i).
260pub fn compute_bpp_mccaskill(sequences: &[Vec<u8>]) -> Result<Vec<SequenceBpp>, String> {
261    let tool = find_tool("mxscarnamod")
262        .ok_or_else(|| {
263            "mxscarnamod not found. Q-INS-i requires the McCaskill base-pair probability program.\n\
264             Install it and ensure it is in your PATH or set MAFFT_BINARIES.\n\
265             See: https://mafft.cbrc.jp/alignment/software/source.html".to_string()
266        })?;
267
268    let mut results = Vec::with_capacity(sequences.len());
269    for (idx, seq) in sequences.iter().enumerate() {
270        // Strip gaps for BPP prediction
271        let ungapped: Vec<u8> = seq.iter().filter(|&&c| c != b'-').copied().collect();
272        let pairs = run_mccaskill(&ungapped, &tool)?;
273        results.push(SequenceBpp { seq_index: idx, pairs });
274    }
275    Ok(results)
276}
277
278/// Compute base-pair probabilities for all sequences using CONTRAfold (X-INS-i).
279pub fn compute_bpp_contrafold(sequences: &[Vec<u8>]) -> Result<Vec<SequenceBpp>, String> {
280    let tool = find_tool("contrafold")
281        .ok_or_else(|| {
282            "contrafold not found. X-INS-i requires CONTRAfold.\n\
283             Install CONTRAfold v2.02+ and ensure it is in your PATH or set MAFFT_BINARIES.\n\
284             See: https://mafft.cbrc.jp/alignment/software/source.html".to_string()
285        })?;
286
287    let mut results = Vec::with_capacity(sequences.len());
288    for (idx, seq) in sequences.iter().enumerate() {
289        let ungapped: Vec<u8> = seq.iter().filter(|&&c| c != b'-').copied().collect();
290        let pairs = run_contrafold(&ungapped, &tool)?;
291        results.push(SequenceBpp { seq_index: idx, pairs });
292    }
293    Ok(results)
294}
295
296#[cfg(test)]
297mod tests {
298    use super::*;
299
300    #[test]
301    fn parse_mccaskill_format() {
302        let input = "0 5 0.85\n1 4 0.72\n2 3 0.005\n";
303        let pairs = parse_mccaskill_output(input);
304        assert_eq!(pairs.len(), 2); // 0.005 filtered out
305        assert_eq!(pairs[0].left, 0);
306        assert_eq!(pairs[0].right, 5);
307        assert!((pairs[0].prob - 0.85).abs() < 1e-10);
308    }
309
310    #[test]
311    fn parse_contrafold_format() {
312        let input = "1 2:0.95 5:0.42\n3 4:0.31\n";
313        let pairs = parse_contrafold_output(input);
314        assert_eq!(pairs.len(), 3);
315        // 1-indexed → 0-indexed
316        assert_eq!(pairs[0].left, 0);
317        assert_eq!(pairs[0].right, 1);
318        assert!((pairs[0].prob - 0.95).abs() < 1e-10);
319    }
320
321    #[test]
322    fn parse_hat3_format() {
323        let input = "0 1 100 5.8 10 20 30 40 info\n2 3 50 2.9 15 25 35 45 info\n";
324        let constraints = parse_hat3(input);
325        assert_eq!(constraints.len(), 2);
326        assert_eq!(constraints[0].0, 0); // i
327        assert_eq!(constraints[0].1, 1); // j
328        assert!((constraints[0].2 - 5.8).abs() < 1e-10); // opt
329    }
330
331    #[test]
332    fn find_tool_nonexistent() {
333        assert!(find_tool("nonexistent_tool_xyz_12345").is_none());
334    }
335}
mafft_core/external.rs

mafft_core/
external.rs
