source: tags/ms_r16q2/PROBE/pt_split.h

Last change on this file was 13970, checked in by westram, 7 years ago
File size: 3.9 KB
Line 
1// ================================================================ //
2//                                                                  //
3//   File      : pt_split.h                                         //
4//   Purpose   :                                                    //
5//                                                                  //
6//   Coded by Ralf Westram (coder@reallysoft.de) in December 2012   //
7//   Institute of Microbiology (Technical University Munich)        //
8//   http://www.arb-home.de/                                        //
9//                                                                  //
10// ================================================================ //
11
12#ifndef PT_SPLIT_H
13#define PT_SPLIT_H
14
15#ifndef PROBE_H
16#include "probe.h"
17#endif
18#ifndef PT_COMPLEMENT_H
19#include "PT_complement.h"
20#endif
21
22inline double get_bond_val(const PT_bond *bond, int probe, int seq) {
23    pt_assert(is_std_base(probe));
24    pt_assert(is_std_base(seq));
25    int idx = ((probe-int(PT_A))*4) + seq-(int)PT_A;
26    pt_assert(idx >= 0 && idx < 16);
27    return bond[idx].val;
28}
29
30class MaxBond {
31    double max_bond[PT_BASES];
32protected:
33    MaxBond() {}
34public:
35    MaxBond(const PT_bond *bond) {
36        for (int base = PT_A; base < PT_BASES; ++base) {
37            pt_assert(is_std_base(base));
38            max_bond[base] = get_bond_val(bond, get_complement(base), base);
39        }
40    }
41
42    double get_max_bond(int base) const {
43        pt_assert(is_std_base(base));
44        return max_bond[base];
45    }
46};
47
48class MismatchWeights : public MaxBond {
49    double weight[PT_BASES][PT_BASES];
50
51    double get_simple_wmismatch(const PT_bond *bonds, char probe, char seq) {
52        double max_bind = get_max_bond(probe);
53        double new_bind = get_bond_val(bonds, get_complement(probe), seq);
54        return max_bind - new_bind;
55    }
56
57    void init(const PT_bond *bonds) {
58        for (int probe = PT_A; probe < PT_BASES; ++probe) {
59            double sum = 0.0;
60            for (int seq = PT_A; seq < PT_BASES; ++seq) {
61                sum += weight[probe][seq] = get_simple_wmismatch(bonds, probe, seq);
62            }
63            weight[probe][PT_N] = sum/4.0;
64        }
65        for (int seq = PT_N; seq < PT_BASES; ++seq) {  // UNCHECKED_LOOP (doesnt matter)
66            double sum = 0.0;
67            for (int probe = PT_A; probe < PT_BASES; ++probe) {
68                sum += weight[probe][seq];
69            }
70            weight[PT_N][seq] = sum/4.0;
71        }
72
73        for (int i = PT_N; i<PT_BASES; ++i) {
74            weight[PT_QU][i] = weight[PT_N][i];
75            weight[i][PT_QU] = weight[i][PT_N];
76        }
77        weight[PT_QU][PT_QU] = weight[PT_N][PT_N];
78    }
79
80public:
81    MismatchWeights(const PT_bond *bonds) : MaxBond(bonds) { init(bonds); }
82    double get(int probe, int seq) const {
83        pt_assert(is_valid_base(probe) && is_valid_base(seq));
84        return weight[probe][seq];
85    }
86
87};
88
89
90class Splits : public MaxBond {
91    bool valid;
92    double split[PT_BASES][PT_BASES];
93
94    double calc_split(const PT_local *locs, char base, char ref) {
95        double max_bind = get_max_bond(base);
96        double new_bind = get_bond_val(locs->bond, get_complement(base), ref);
97
98        if (new_bind < locs->split)
99            return new_bind-max_bind; // negative values indicate split
100        // this mismatch splits the probe in several domains
101        return max_bind-new_bind;
102    }
103
104public:
105    Splits() : valid(false) {}
106    Splits(const PT_local *locs) : MaxBond(locs->bond), valid(true) {
107        for (int base = PT_A; base < PT_BASES; ++base) {
108            for (int ref = PT_A; ref < PT_BASES; ++ref) {
109                split[base][ref] = calc_split(locs, base, ref);
110            }
111        }
112    }
113
114    double check(char base, char ref) const {
115        pt_assert(valid);
116
117        pt_assert(is_std_base(base));
118        pt_assert(is_std_base(ref));
119
120        return split[safeCharIndex(base)][safeCharIndex(ref)];
121    }
122};
123
124#else
125#error pt_split.h included twice
126#endif // PT_SPLIT_H
Note: See TracBrowser for help on using the repository browser.