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
|
#include "wYangModel.h"
#include "codon.h"
#include "readDatMatrix.h" // for the normalizeQ function.
wYangModel::wYangModel(const MDOUBLE inW, const MDOUBLE inK,bool globalW, codon * coAlph):
_w(inW),_k(inK),_globalW(globalW),_coAlpha(NULL){
_coAlpha = (codon*)(coAlph->clone());
codonUtility::initSubMatrices(*_coAlpha);
homogenousFreq();
_Q.resize(alphabetSize());
for (int z=0; z < _Q.size();++z) _Q[z].resize(alphabetSize(),0.0);
updateQ();
}
wYangModel::wYangModel(const MDOUBLE inW, const MDOUBLE inK, const Vdouble& freq,bool globalW, codon * coAlph):
_w(inW),_k(inK),_globalW(globalW),_freq(freq),_coAlpha(NULL){
_coAlpha = (codon*)(coAlph->clone());
_Q.resize(alphabetSize());
codonUtility::initSubMatrices(*_coAlpha);
for (int z=0; z < _Q.size();++z) _Q[z].resize(alphabetSize(),0.0);
updateQ();
}
wYangModel& wYangModel::operator=(const wYangModel &other) {
_w = other._w;
_k = other._k;
_q2pt = other._q2pt;
_Q = other._Q;
_globalW = other._globalW;
_freq = other._freq;
if (_coAlpha) delete _coAlpha;
if (other._coAlpha)
_coAlpha = (codon*)(other._coAlpha->clone());
else
_coAlpha = NULL;
return *this;
}
void wYangModel::updateQ() {
int i,j;
MDOUBLE sum=0.0;
for (i=0; i < _Q.size();++i) {
for (j=i+1; j < _Q.size();++j) {
MDOUBLE val;
if (codonUtility::codonReplacement(i,j) == codonUtility::non_synonymous) {
if (codonUtility::codonDiff(i,j) == codonUtility::tr) val = _k*_w;
else if (codonUtility::codonDiff(i,j) == codonUtility::tv) val = _w;
else val = 0;//more than one substitution.
}
else {//synonymous
if (codonUtility::codonDiff(i,j) == codonUtility::tr) val = _k;
else if (codonUtility::codonDiff(i,j) == codonUtility::tv) val = 1;
else val = 0;//more than one substitution.
}
_Q[i][j] = val * _freq[j];
_Q[j][i] = val * _freq[i];
}
_Q[i][i] = 0.0; //temporary value
}
// filling the diagonal
for (i=0; i < _Q.size(); ++i){
sum = 0.0;
for (j=0; j < _Q.size(); ++j) {
sum += _Q[i][j];
}
_Q[i][i] = -sum;
}
if (_globalW == true) // w is not distributed, only one Q matrix
normalizeQ(_Q,_freq);
_q2pt.fillFromRateMatrix(_freq,_Q);
}
void wYangModel::norm(MDOUBLE scale){
for (int i=0; i < _Q.size(); ++i) {
for (int j=0; j < _Q.size(); ++j) {
_Q[i][j] *=scale;
}
}
_q2pt.fillFromRateMatrix(_freq,_Q);
}
MDOUBLE wYangModel::sumPijQij(){
MDOUBLE sum=0.0;
for (int i=0; i < _Q.size(); ++i) {
sum -= (_Q[i][i])*_freq[i];
}
return sum;
}
|
