rustyhmmer 0.1.2

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
116
117
118
119
120
121
122
123
124
125
126
127

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;
        // Parallelize over PROFILES (coarse-grained), not over sequences within
        // each profile. A profile DB (~10^4-10^5 models) offers far more
        // parallelism than a single genome's sequence set, and most profiles are
        // rejected cheaply by the MSV filter — so a per-profile `par_iter` over
        // sequences paid rayon fork/join overhead on ~10^5 near-empty dispatches
        // (measured ~35% core usage on 16 threads). One coarse dispatch over the
        // profile DB — the scheme HMMER/pyhmmer use — keeps the cores saturated.
        // `par_iter().map(..).collect::<Vec<_>>()` is order-preserving
        // (IndexedParallelIterator), so the concatenated output is byte-identical
        // to the previous serial (model-then-sequence) order.
        let per_model: Vec<Vec<HmmHit>> = self
            .models
            .par_iter()
            .map(|model| {
                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,
                };
                let mut local = Vec::new();
                if !matches!(self.cutoff, Cutoff::Evalue(_)) && seq_bit_t.is_none() {
                    return local;
                }

                let hits: Vec<Hit> = seqs.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];
                    local.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,
                    });
                }
                local
            })
            .collect();
        per_model.into_iter().flatten().collect()
    }
}