use crate::alphabet::K;
pub const TMM: usize = 0;
pub const TMI: usize = 1;
pub const TMD: usize = 2;
pub const TIM: usize = 3;
pub const TII: usize = 4;
pub const TDM: usize = 5;
pub const TDD: usize = 6;
#[derive(Debug, Clone, Copy, Default)]
pub struct EvParams {
pub msv_mu: f64,
pub msv_lambda: f64,
pub vit_mu: f64,
pub vit_lambda: f64,
pub fwd_tau: f64,
pub fwd_lambda: f64,
}
#[derive(Debug, Clone, Copy, Default)]
pub struct Cutoffs {
pub ga: Option<(f64, f64)>,
pub tc: Option<(f64, f64)>,
pub nc: Option<(f64, f64)>,
}
#[derive(Debug, Clone)]
pub struct P7Hmm {
pub name: String,
pub acc: Option<String>,
pub desc: Option<String>,
pub m: usize,
pub mat: Vec<[f32; K]>,
pub ins: Vec<[f32; K]>,
pub t: Vec<[f32; 7]>,
pub compo: Option<[f32; K]>,
pub evparam: EvParams,
pub cutoffs: Cutoffs,
}
#[derive(Debug)]
pub struct ParseError(pub String);
impl std::fmt::Display for ParseError {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
write!(f, "hmm parse error: {}", self.0)
}
}
impl std::error::Error for ParseError {}
#[inline]
fn score_to_prob(tok: &str) -> f32 {
if tok == "*" {
0.0
} else {
(-tok.parse::<f32>().unwrap_or(f32::INFINITY)).exp()
}
}
impl P7Hmm {
pub fn read_one<I>(lines: &mut std::iter::Peekable<I>) -> Result<Option<Self>, ParseError>
where
I: Iterator<Item = std::io::Result<String>>,
{
let magic = loop {
match lines.next() {
None => return Ok(None),
Some(Ok(l)) if l.trim().is_empty() => continue,
Some(Ok(l)) => break l,
Some(Err(e)) => return Err(ParseError(format!("io: {e}"))),
}
};
if !magic.starts_with("HMMER3/") {
return Err(ParseError(format!("expected HMMER3/f header, got {magic:?}")));
}
let mut name = String::new();
let mut acc = None;
let mut desc = None;
let mut m: usize = 0;
let mut alph = String::new();
let mut ev = EvParams::default();
let mut cut = Cutoffs::default();
loop {
let line = next_line(lines)?
.ok_or_else(|| ParseError("EOF in header".into()))?;
let t = line.trim();
if t.starts_with("HMM ") || t == "HMM" {
break;
}
let p: Vec<&str> = t.split_whitespace().collect();
if p.is_empty() {
continue;
}
match p[0] {
"NAME" => name = p.get(1).unwrap_or(&"").to_string(),
"ACC" => acc = p.get(1).map(|s| s.to_string()),
"DESC" => desc = Some(t[p[0].len()..].trim().to_string()),
"LENG" => m = p.get(1).and_then(|s| s.parse().ok()).unwrap_or(0),
"ALPH" => alph = p.get(1).unwrap_or(&"").to_string(),
"STATS" if p.len() >= 5 && p[1] == "LOCAL" => {
let a: f64 = p[3].parse().unwrap_or(0.0);
let b: f64 = p[4].parse().unwrap_or(0.0);
match p[2] {
"MSV" => { ev.msv_mu = a; ev.msv_lambda = b; }
"VITERBI" => { ev.vit_mu = a; ev.vit_lambda = b; }
"FORWARD" => { ev.fwd_tau = a; ev.fwd_lambda = b; }
_ => {}
}
}
"GA" if p.len() >= 3 => cut.ga = Some((pf(p[1]), pf(p[2]))),
"TC" if p.len() >= 3 => cut.tc = Some((pf(p[1]), pf(p[2]))),
"NC" if p.len() >= 3 => cut.nc = Some((pf(p[1]), pf(p[2]))),
_ => {}
}
}
if !alph.eq_ignore_ascii_case("amino") {
return Err(ParseError(format!(
"rustyhmmer is amino-only; model {name:?} has ALPH {alph:?}"
)));
}
if m == 0 {
return Err(ParseError(format!("model {name:?} has LENG 0")));
}
let mut hmm = P7Hmm {
name,
acc,
desc,
m,
mat: vec![[0.0; K]; m + 1],
ins: vec![[0.0; K]; m + 1],
t: vec![[0.0; 7]; m + 1],
compo: None,
evparam: ev,
cutoffs: cut,
};
let _ = next_line(lines)?;
let first = next_line(lines)?
.ok_or_else(|| ParseError("EOF before COMPO/inserts".into()))?;
let fp: Vec<&str> = first.split_whitespace().collect();
let node0_ins_line = if fp.first() == Some(&"COMPO") {
let mut c = [0.0f32; K];
for i in 0..K {
c[i] = score_to_prob(fp.get(i + 1).copied().unwrap_or("*"));
}
hmm.compo = Some(c);
next_line(lines)?.ok_or_else(|| ParseError("EOF at node-0 inserts".into()))?
} else {
first
};
parse_emission(&node0_ins_line, &mut hmm.ins[0])?;
let node0_t = next_line(lines)?
.ok_or_else(|| ParseError("EOF at node-0 transitions".into()))?;
parse_transition(&node0_t, &mut hmm.t[0])?;
for k in 1..=m {
let mline = next_line(lines)?
.ok_or_else(|| ParseError(format!("EOF at node {k} match")))?;
if mline.trim_start().starts_with("//") {
return Err(ParseError(format!("model ended early at node {k}")));
}
let mp: Vec<&str> = mline.split_whitespace().collect();
for i in 0..K {
hmm.mat[k][i] = score_to_prob(mp.get(i + 1).copied().unwrap_or("*"));
}
let iline = next_line(lines)?
.ok_or_else(|| ParseError(format!("EOF at node {k} insert")))?;
parse_emission(&iline, &mut hmm.ins[k])?;
let tline = next_line(lines)?
.ok_or_else(|| ParseError(format!("EOF at node {k} transition")))?;
parse_transition(&tline, &mut hmm.t[k])?;
}
if let Some(Ok(l)) = lines.peek() {
if l.trim() == "//" {
let _ = lines.next();
}
}
Ok(Some(hmm))
}
pub fn read_all(path: &str) -> Result<Vec<Self>, ParseError> {
use std::io::BufRead;
let f = std::fs::File::open(path).map_err(|e| ParseError(format!("open {path}: {e}")))?;
let mut lines = std::io::BufReader::new(f).lines().peekable();
let mut out = Vec::new();
while let Some(hmm) = Self::read_one(&mut lines)? {
out.push(hmm);
}
Ok(out)
}
}
#[inline]
fn pf(tok: &str) -> f64 {
tok.parse().unwrap_or(0.0)
}
fn next_line<I>(lines: &mut std::iter::Peekable<I>) -> Result<Option<String>, ParseError>
where
I: Iterator<Item = std::io::Result<String>>,
{
match lines.next() {
None => Ok(None),
Some(Ok(l)) => Ok(Some(l)),
Some(Err(e)) => Err(ParseError(format!("io: {e}"))),
}
}
fn parse_emission(line: &str, out: &mut [f32; K]) -> Result<(), ParseError> {
let p: Vec<&str> = line.split_whitespace().collect();
if p.len() < K {
return Err(ParseError(format!("emission line has {} fields, need {K}", p.len())));
}
for i in 0..K {
out[i] = score_to_prob(p[i]);
}
Ok(())
}
fn parse_transition(line: &str, out: &mut [f32; 7]) -> Result<(), ParseError> {
let p: Vec<&str> = line.split_whitespace().collect();
if p.len() < 7 {
return Err(ParseError(format!("transition line has {} fields, need 7", p.len())));
}
for i in 0..7 {
out[i] = score_to_prob(p[i]);
}
Ok(())
}
#[cfg(test)]
mod tests {
use super::*;
fn td(name: &str) -> String {
format!("{}/testdata/{}", env!("CARGO_MANIFEST_DIR"), name)
}
#[test]
fn parse_globins4() {
let hmms = P7Hmm::read_all(&td("globins4.hmm")).unwrap();
assert_eq!(hmms.len(), 1);
let h = &hmms[0];
assert_eq!(h.name, "globins4");
assert_eq!(h.m, 149);
assert!((h.evparam.msv_mu - -9.9014).abs() < 1e-4);
assert!((h.evparam.msv_lambda - 0.70957).abs() < 1e-5);
assert!((h.evparam.fwd_tau - -4.1637).abs() < 1e-4);
assert!(h.cutoffs.ga.is_none());
assert!(h.compo.is_some());
for k in 1..=h.m {
let s: f32 = h.mat[k].iter().sum();
assert!((s - 1.0).abs() < 1e-2, "node {k} emission sum {s}");
}
}
#[test]
fn parse_fn3_cutoffs() {
let hmms = P7Hmm::read_all(&td("fn3.hmm")).unwrap();
let h = &hmms[0];
let (ga_seq, ga_dom) = h.cutoffs.ga.expect("fn3 has GA");
assert!((ga_seq - 8.00).abs() < 1e-6);
assert!((ga_dom - 7.20).abs() < 1e-6);
}
}