/
pairhmm.rs
103 lines (80 loc) · 2.52 KB
/
pairhmm.rs
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
#![feature(test)]
extern crate test;
use bio::stats::pairhmm::*;
use bio::stats::{LogProb, Prob};
use test::Bencher;
static TEXT: &[u8] = b"GATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCTCTCCATGC\
ATTTGGTATTTTCGTCTGGGGGGTATGCACGCGATAGCATTGCGAGACGCTGGAGCCGGAGCACCCTATGTCGCAGTAT\
CTGTCTTTGATTCCTGCCTCATCCTATTATTTATCGCACCTACGTTCAATATTACAGGCGAACATACTTACTAAAGTGT";
static PATTERN: &[u8] = b"GGGTATGCACGCGATAGCATTGCGAGACGCTGGAGCCGGAGCACCCTATGTCGC";
// Single base insertion and deletion rates for R1 according to Schirmer et al.
// BMC Bioinformatics 2016, 10.1186/s12859-016-0976-y
static PROB_ILLUMINA_INS: Prob = Prob(2.8e-6);
static PROB_ILLUMINA_DEL: Prob = Prob(5.1e-6);
static PROB_ILLUMINA_SUBST: Prob = Prob(0.0021);
fn prob_emit_x_or_y() -> LogProb {
LogProb::from(Prob(1.0) - PROB_ILLUMINA_SUBST)
}
pub struct TestEmissionParams {
x: &'static [u8],
y: &'static [u8],
}
impl EmissionParameters for TestEmissionParams {
fn prob_emit_xy(&self, i: usize, j: usize) -> XYEmission {
if self.x[i] == self.y[j] {
XYEmission::Match(LogProb::from(Prob(1.0) - PROB_ILLUMINA_SUBST))
} else {
XYEmission::Mismatch(LogProb::from(PROB_ILLUMINA_SUBST / Prob(3.0)))
}
}
fn prob_emit_x(&self, _: usize) -> LogProb {
prob_emit_x_or_y()
}
fn prob_emit_y(&self, _: usize) -> LogProb {
prob_emit_x_or_y()
}
fn len_x(&self) -> usize {
self.x.len()
}
fn len_y(&self) -> usize {
self.y.len()
}
}
pub struct SemiglobalGapParams;
impl GapParameters for SemiglobalGapParams {
fn prob_gap_x(&self) -> LogProb {
LogProb::from(PROB_ILLUMINA_INS)
}
fn prob_gap_y(&self) -> LogProb {
LogProb::from(PROB_ILLUMINA_DEL)
}
fn prob_gap_x_extend(&self) -> LogProb {
LogProb::ln_zero()
}
fn prob_gap_y_extend(&self) -> LogProb {
LogProb::ln_zero()
}
}
pub struct SemiglobalAlignment;
impl StartEndGapParameters for SemiglobalAlignment {
fn free_start_gap_x(&self) -> bool {
true
}
fn free_end_gap_x(&self) -> bool {
true
}
}
#[bench]
fn pairhmm_semiglobal(b: &mut Bencher) {
let emission_params = TestEmissionParams {
x: TEXT,
y: PATTERN,
};
let gap_params = SemiglobalGapParams;
let mut pair_hmm = PairHMM::new(&gap_params);
pair_hmm.prob_related(&emission_params, &SemiglobalAlignment, Some(4));
b.iter(|| {
let p = pair_hmm.prob_related(&emission_params, &SemiglobalAlignment, Some(4));
assert!(*p <= 0.0);
});
}