use crate::pipeline::{dombias_bits, Hit};
pub fn fmt_g2(x: f64) -> String {
if x == 0.0 {
return "0".to_string();
}
let neg = x < 0.0;
let ax = x.abs();
let sci = format!("{:.1e}", ax);
let (mant, exp_str) = sci.split_once('e').unwrap();
let exp: i32 = exp_str.parse().unwrap();
let body = if exp >= -4 && exp < 2 {
let dec = (1 - exp).max(0) as usize;
let f = format!("{:.*}", dec, ax);
if f.contains('.') {
f.trim_end_matches('0').trim_end_matches('.').to_string()
} else {
f
}
} else {
let mant = if mant.contains('.') {
mant.trim_end_matches('0').trim_end_matches('.')
} else {
mant
};
let (sign, digits) = match exp_str.strip_prefix('-') {
Some(d) => ("-", d),
None => ("+", exp_str.trim_start_matches('+')),
};
format!("{}e{}{:0>2}", mant, sign, digits)
};
if neg {
format!("-{body}")
} else {
body
}
}
#[derive(Clone, Copy)]
pub enum ReportThresh {
Evalue,
Bits { dom_t: f64 },
}
#[derive(Clone, Copy)]
pub struct Widths {
pub tnamew: usize,
pub taccw: usize,
pub qnamew: usize,
pub qaccw: usize,
}
impl Widths {
pub fn compute(qname: &str, qacc: &str, reported: &[&Hit]) -> Self {
let tnamew = 20.max(reported.iter().map(|h| h.name.len()).max().unwrap_or(0));
let taccw = 10;
let qnamew = 20.max(qname.len());
let qaccw = if qacc.is_empty() { 10 } else { 10.max(qacc.len()) };
Widths { tnamew, taccw, qnamew, qaccw }
}
}
pub fn header(w: &Widths) -> String {
let grp = w.tnamew + w.qnamew + w.taccw + w.qaccw + 2;
let mut s = String::new();
s.push('#');
s.push_str(&format!(
"{:grp$} {:>22} {:>22} {:>33}\n",
"", "--- full sequence ----", "--- best 1 domain ----", "--- domain number estimation ----",
grp = grp,
));
s.push('#');
s.push_str(&format!(
"{:<tn$} {:<ta$} {:<qn$} {:<qa$} {:>9} {:>6} {:>5} {:>9} {:>6} {:>5} {:>5} {:>3} {:>3} {:>3} {:>3} {:>3} {:>3} {:>3} {}\n",
" target name", "accession", "query name", "accession",
" E-value", " score", " bias", " E-value", " score", " bias",
"exp", "reg", "clu", " ov", "env", "dom", "rep", "inc", "description of target",
tn = w.tnamew - 1, ta = w.taccw, qn = w.qnamew, qa = w.qaccw,
));
s.push('#');
s.push_str(&format!(
"{:>tn$} {:>ta$} {:>qn$} {:>qa$} {:>9} {:>6} {:>5} {:>9} {:>6} {:>5} {:>5} {:>3} {:>3} {:>3} {:>3} {:>3} {:>3} {:>3} {}\n",
"-------------------", "----------", "--------------------", "----------",
"---------", "------", "-----", "---------", "------", "-----", "-----",
"---", "---", "---", "---", "---", "---", "---", "---------------------",
tn = w.tnamew - 1, ta = w.taccw, qn = w.qnamew, qa = w.qaccw,
));
s
}
pub fn format_row(hit: &Hit, qname: &str, qacc: &str, z: f64, domz: f64, thresh: ReportThresh, w: &Widths) -> String {
let e_full = hit.lnp.exp() * z;
let bd = &hit.domains[hit.best];
let e_dom = bd.lnp.exp() * z;
let (nrep, ninc) = match thresh {
ReportThresh::Evalue => {
let seq_included = e_full <= 0.01;
let nrep = hit.domains.iter().filter(|d| d.lnp.exp() * domz <= 10.0).count() as i32;
let ninc = if seq_included {
hit.domains.iter().filter(|d| d.lnp.exp() * domz <= 0.01).count() as i32
} else {
0
};
(nrep, ninc)
}
ReportThresh::Bits { dom_t } => {
let n = hit.domains.iter().filter(|d| d.bitscore as f64 >= dom_t).count() as i32;
(n, n)
}
};
let desc = if hit.desc.is_empty() { "-" } else { &hit.desc };
let qacc = if qacc.is_empty() { "-" } else { qacc };
format!(
"{tname:<tnamew$} {tacc:<taccw$} {qname:<qnamew$} {qacc:<qaccw$} {efull:>9} {sc:6.1} {bi:5.1} \
{edom:>9} {dsc:6.1} {dbi:5.1} {exp:5.1} {reg:3} {clu:3} {ov:3} {env:3} {dom:3} {rep:3} {inc:3} {desc}",
tname = hit.name,
tacc = "-",
qname = qname,
qacc = qacc,
tnamew = w.tnamew,
taccw = w.taccw,
qnamew = w.qnamew,
qaccw = w.qaccw,
efull = fmt_g2(e_full),
sc = hit.score,
bi = hit.bias(),
edom = fmt_g2(e_dom),
dsc = bd.bitscore,
dbi = dombias_bits(bd.dombias),
exp = hit.nexpected,
reg = hit.nregions,
clu = hit.nclustered,
ov = hit.noverlaps,
env = hit.nenvelopes,
dom = hit.ndom,
rep = nrep,
inc = ninc,
desc = desc,
)
}
pub fn sort_hits(hits: &mut [Hit]) {
hits.sort_by(|a, b| a.lnp.partial_cmp(&b.lnp).unwrap());
}
#[cfg(test)]
mod tests {
use super::*;
use crate::hmmfile::P7Hmm;
use crate::pipeline::Model;
use crate::seqio::read_fasta;
#[test]
fn globins_tblout_rows_byte_identical() {
let base = env!("CARGO_MANIFEST_DIR");
let hmm = P7Hmm::read_all(&format!("{base}/testdata/globins4.hmm"))
.unwrap()
.pop()
.unwrap();
let qname = hmm.name.clone();
let model = Model::new(hmm);
let seqs = read_fasta(&format!("{base}/testdata/globins45.fa")).unwrap();
let z = seqs.len() as f64;
let mut hits: Vec<Hit> = seqs.iter().filter_map(|s| model.search_one(s)).collect();
sort_hits(&mut hits);
let reported: Vec<&Hit> = hits.iter().collect();
let w = Widths::compute(&qname, "", &reported);
let mine: Vec<String> = hits
.iter()
.map(|h| format_row(h, &qname, "", z, z, ReportThresh::Evalue, &w))
.collect();
let golden = std::fs::read_to_string(format!("{base}/golden/globins4.tblout")).unwrap();
let gold_rows: Vec<String> = golden
.lines()
.filter(|l| !l.starts_with('#'))
.map(|l| l.trim_end().to_string())
.collect();
assert_eq!(mine.len(), gold_rows.len(), "row count");
let mut nbad = 0;
for (m, g) in mine.iter().zip(gold_rows.iter()) {
if m.trim_end() != g.trim_end() {
nbad += 1;
if nbad <= 5 {
eprintln!("MINE : {:?}", m.trim_end());
eprintln!("GOLD : {:?}", g);
}
}
}
assert_eq!(nbad, 0, "{nbad}/{} rows differ from golden", gold_rows.len());
}
}