rustyhmmer 0.1.1

Pure-Rust HMMER3 hmmsearch with byte-identical --tblout output
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
112
113
114
115

use crate::pipeline::{dombias_bits, Hit, Model};
use crate::seqio::Seq;
use rayon::prelude::*;

#[derive(Clone, Copy, Debug, PartialEq)]
pub enum Cutoff {
    GatheringGa,
    TrustedTc,
    NoiseNc,
    Evalue(f64),
}

#[derive(Clone, Debug)]
pub struct HmmHit {
    pub query_name: String,
    pub query_acc: String,
    pub model_len: usize,
    pub target_name: String,
    pub seq_score: f32,
    pub seq_bias: f32,
    pub seq_evalue: f64,
    pub best_dom_score: f32,
    pub best_dom_bias: f32,
    pub best_dom_evalue: f64,
    pub env_from: i64,
    pub env_to: i64,
    pub n_domains: i32,
}

pub struct HmmAnnotator {
    models: Vec<Model>,
    cutoff: Cutoff,
}

impl HmmAnnotator {
    pub fn from_hmm_file(path: &str) -> Result<Self, String> {
        let hmms = crate::hmmfile::P7Hmm::read_all(path).map_err(|e| e.to_string())?;
        if hmms.is_empty() {
            return Err(format!("no models in {path}"));
        }
        Ok(Self {
            models: hmms.into_iter().map(Model::new).collect(),
            cutoff: Cutoff::Evalue(10.0),
        })
    }

    pub fn from_models(models: Vec<Model>) -> Self {
        Self {
            models,
            cutoff: Cutoff::Evalue(10.0),
        }
    }

    pub fn with_cutoff(mut self, cutoff: Cutoff) -> Self {
        self.cutoff = cutoff;
        self
    }

    pub fn len(&self) -> usize {
        self.models.len()
    }
    pub fn is_empty(&self) -> bool {
        self.models.is_empty()
    }

    pub fn search(&self, seqs: &[Seq]) -> Vec<HmmHit> {
        let z = seqs.len() as f64;
        let mut out = Vec::new();
        for model in &self.models {
            let seq_bit_t: Option<f64> = match self.cutoff {
                Cutoff::GatheringGa => model.hmm.cutoffs.ga.map(|(t1, _)| t1),
                Cutoff::TrustedTc => model.hmm.cutoffs.tc.map(|(t1, _)| t1),
                Cutoff::NoiseNc => model.hmm.cutoffs.nc.map(|(t1, _)| t1),
                Cutoff::Evalue(_) => None,
            };
            if !matches!(self.cutoff, Cutoff::Evalue(_)) && seq_bit_t.is_none() {
                continue;
            }

            let hits: Vec<Hit> = seqs
                .par_iter()
                .filter_map(|s| model.search_one(s))
                .collect();

            let qname = &model.hmm.name;
            let qacc = model.hmm.acc.clone().unwrap_or_default();
            for h in &hits {
                let reportable = match self.cutoff {
                    Cutoff::Evalue(e) => h.lnp.exp() * z <= e,
                    _ => (h.score as f64) >= seq_bit_t.unwrap(),
                };
                if !reportable {
                    continue;
                }
                let bd = &h.domains[h.best];
                out.push(HmmHit {
                    query_name: qname.clone(),
                    query_acc: qacc.clone(),
                    model_len: model.hmm.m,
                    target_name: h.name.clone(),
                    seq_score: h.score,
                    seq_bias: h.bias(),
                    seq_evalue: h.lnp.exp() * z,
                    best_dom_score: bd.bitscore,
                    best_dom_bias: dombias_bits(bd.dombias),
                    best_dom_evalue: bd.lnp.exp() * z,
                    env_from: bd.ienv,
                    env_to: bd.jenv,
                    n_domains: h.ndom,
                });
            }
        }
        out
    }
}