use crate::bias::BiasFilter;
use crate::dp::{build_local_profile, p7_domaindef_local, p7_flogsum};
use crate::forward::ForwardFilter;
use crate::hmmfile::P7Hmm;
use crate::msv::{null_one, MsvProfile};
use crate::seqio::Seq;
use crate::vitfilter::VitFilter;
const LOG2: f64 = std::f64::consts::LN_2;
const LOG2R: f64 = 1.0 / std::f64::consts::LN_2;
const OMEGA: f64 = 1.0 / 256.0;
#[inline]
fn exp_logsurv(x: f64, mu: f64, lambda: f64) -> f64 {
if x < mu {
0.0
} else {
-lambda * (x - mu)
}
}
pub struct DomHit {
pub ienv: i64,
pub jenv: i64,
pub bitscore: f32,
pub dombias: f32,
pub lnp: f64,
}
pub struct Hit {
pub name: String,
pub desc: String,
pub score: f32,
pub pre_score: f32,
pub lnp: f64,
pub nexpected: f32,
pub nregions: i32,
pub nclustered: i32,
pub noverlaps: i32,
pub nenvelopes: i32,
pub ndom: i32,
pub best: usize,
pub domains: Vec<DomHit>,
}
impl Hit {
pub fn bias(&self) -> f32 {
self.pre_score - self.score
}
}
pub struct Model {
pub hmm: P7Hmm,
pub msv: MsvProfile,
pub vit: VitFilter,
pub bias: Option<BiasFilter>,
pub fwd: ForwardFilter,
pub f1: f64,
pub f2: f64,
pub f3: f64,
pub do_biasfilter: bool,
pub do_null2: bool,
}
impl Model {
pub fn new(hmm: P7Hmm) -> Self {
let msv = MsvProfile::build(&hmm);
let vit = VitFilter::build(&hmm);
let bias = BiasFilter::build(&hmm);
let fwd = ForwardFilter::build(&hmm);
Model {
hmm,
msv,
vit,
bias,
fwd,
f1: 0.02,
f2: 1e-3,
f3: 1e-5,
do_biasfilter: true,
do_null2: true,
}
}
pub fn search_one(&self, seq: &Seq) -> Option<Hit> {
let dsq = &seq.dsq;
let l = seq.len();
if l == 0 {
return None;
}
let nullsc = null_one(l) as f64;
let usc = self
.msv
.ssv_score(dsq, l)
.unwrap_or_else(|| self.msv.msv_score(dsq, l));
if self.msv.pvalue(usc, nullsc as f32) > self.f1 {
return None;
}
let filtersc = if self.do_biasfilter {
let fs = self.bias.as_ref().unwrap().filter_score(dsq, l);
if self.msv.pvalue(usc, fs) > self.f1 {
return None;
}
fs
} else {
nullsc as f32
};
let msv_p = self.msv.pvalue(usc, filtersc);
if msv_p > self.f2 {
let vfsc = self.vit.vit_score(dsq, l);
if self.vit.pvalue(vfsc, filtersc, l) > self.f2 {
return None;
}
}
let fwdsc = self.fwd.score(dsq, l);
if self.fwd.pvalue(fwdsc, filtersc) > self.f3 {
return None;
}
let mut gm = build_local_profile(&self.hmm, l);
let dd = p7_domaindef_local(&self.fwd, &mut gm, dsq, l, self.do_null2);
let domains = &dd.domains;
if domains.is_empty() {
return None;
}
let fwdsc = fwdsc as f64;
let ln = l as f64;
let seqbias_fwd = if self.do_null2 {
let n2sc_sum = crate::dp::esl_vec_fsum(&dd.n2sc) as f64;
p7_flogsum(0.0, (OMEGA.ln() + n2sc_sum) as f32) as f64
} else {
0.0
};
let mut pre_score = (fwdsc - nullsc) / LOG2;
let mut seq_score = (fwdsc - (nullsc + seqbias_fwd)) / LOG2;
let mut sum_nats = 0.0f64;
let mut ld = 0i64;
let mut dombias_sum = 0.0f64;
for d in domains {
let significant = if self.do_null2 {
(d.envsc - d.domcorrection) as f64 > 0.0
} else {
d.envsc as f64 > 0.0
};
if significant {
sum_nats += d.envsc as f64;
ld += d.jenv - d.ienv + 1;
dombias_sum += d.domcorrection as f64;
}
}
let seqbias_sum = if self.do_null2 {
p7_flogsum(0.0, (OMEGA.ln() + dombias_sum) as f32) as f64
} else {
0.0
};
sum_nats += (ln - ld as f64) * (ln / (ln + 3.0)).ln();
let pre2_score = (sum_nats - nullsc) / LOG2;
let sum_score = (sum_nats - (nullsc + seqbias_sum)) / LOG2;
if ld > 0 && sum_score > seq_score {
seq_score = sum_score;
pre_score = pre2_score;
}
let lnp = exp_logsurv(seq_score, self.fwd.ftau, self.fwd.flambda);
let mut dom_hits = Vec::with_capacity(domains.len());
let mut best = 0usize;
let mut best_bits = f32::NEG_INFINITY;
for (di, d) in domains.iter().enumerate() {
let ld_d = d.jenv - d.ienv + 1;
let bit_nats = d.envsc as f64 + (ln - ld_d as f64) * (ln / (ln + 3.0)).ln();
let dombias = if self.do_null2 {
p7_flogsum(0.0, (OMEGA.ln() + d.domcorrection as f64) as f32) as f64
} else {
0.0
};
let bitscore = ((bit_nats - (nullsc + dombias)) / LOG2) as f32;
let dlnp = exp_logsurv(
bitscore as f64,
self.fwd.ftau,
self.fwd.flambda,
);
if bitscore > best_bits {
best_bits = bitscore;
best = di;
}
dom_hits.push(DomHit {
ienv: d.ienv,
jenv: d.jenv,
bitscore,
dombias: dombias as f32,
lnp: dlnp,
});
}
Some(Hit {
name: seq.name.clone(),
desc: seq.desc.clone(),
score: seq_score as f32,
pre_score: pre_score as f32,
lnp,
nexpected: dd.nexpected,
nregions: dd.nregions,
nclustered: dd.nclustered,
noverlaps: dd.noverlaps,
nenvelopes: dd.nenvelopes,
ndom: domains.len() as i32,
best,
domains: dom_hits,
})
}
}
pub fn dombias_bits(dombias_nats: f32) -> f32 {
(dombias_nats as f64 * LOG2R) as f32
}
#[cfg(test)]
mod tests {
use super::*;
use crate::seqio::read_fasta;
#[test]
fn globins_scores_match_golden() {
let hmm = P7Hmm::read_all(&format!("{}/testdata/globins4.hmm", env!("CARGO_MANIFEST_DIR")))
.unwrap()
.pop()
.unwrap();
let model = Model::new(hmm);
let seqs = read_fasta(&format!("{}/testdata/globins45.fa", env!("CARGO_MANIFEST_DIR")))
.unwrap();
let z = 45.0_f64;
let mut found = 0;
for s in &seqs {
if let Some(h) = model.search_one(s) {
found += 1;
if h.name == "MYG_ESCGI" {
let e = h.lnp.exp() * z;
eprintln!("MYG_ESCGI: score={:.1} bias={:.1} E={:.2e}", h.score, h.bias(), e);
assert!((h.score - 215.6).abs() < 0.15, "score {:.2} != 215.6", h.score);
assert!((h.bias() - 2.9).abs() < 0.15, "bias {:.2} != 2.9", h.bias());
assert!((e / 8.7e-67).ln().abs() < 0.05, "E {:.3e} != 8.7e-67", e);
}
}
}
assert_eq!(found, 45, "expected 45 hits, got {found}");
}
}