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
104
105
106
107
108
109
110
111
112
113
|
/*
* modelcodonparametric.cpp
*
* Created on: May 29, 2013
* Author: minh
*/
#include "modelcodonparametric.h"
ModelCodonParametric::ModelCodonParametric(const char *model_name, string model_params,
StateFreqType freq, string freq_params, PhyloTree *tree, bool count_rates) :
ModelCodon(tree, count_rates)
{
init(model_name, model_params, freq, freq_params);
}
ModelCodonParametric::~ModelCodonParametric() {
}
void ModelCodonParametric::init(const char *model_name, string model_params, StateFreqType freq, string freq_params)
{
StateFreqType def_freq = FREQ_UNKNOWN;
name = full_name = model_name;
string name_upper = model_name;
for (string::iterator it = name_upper.begin(); it != name_upper.end(); it++)
(*it) = toupper(*it);
if (name_upper == "JCC") {
name = "JCC";
def_freq = FREQ_EQUAL;
full_name = "JC-like codon model";
} else if (name_upper == "MG") {
initMG94();
} else if (name_upper == "GY") {
initGY94();
} else {
//cout << "User-specified model "<< model_name << endl;
readParameters(model_name);
//name += " (user-defined)";
}
if (freq_params != "") {
readStateFreq(freq_params);
}
if (model_params != "") {
readRates(model_params);
}
if (freq == FREQ_UNKNOWN || def_freq == FREQ_EQUAL) freq = def_freq;
ModelCodon::init(freq);
}
void ModelCodonParametric::initMG94() {
/* Muse-Gaut 1994 model with 1 parameters: omega */
int i,j;
IntVector group;
for (i = 0; i < num_states-1; i++) {
for (j = i+1; j < num_states; j++) {
if (isMultipleSubst(i, j))
group.push_back(0); // multiple substitution
else if (isSynonymous(i, j))
group.push_back(1); // synonymous substitution
else
group.push_back(2); // non-synonymous substitution
}
}
setRateGroup(group);
// set zero rate for multiple substitution and 1 for synonymous substitution
setRateGroupConstraint("x0=0,x1=1");
}
void ModelCodonParametric::initGY94() {
/* Yang-Nielsen 1998 model (also known as Goldman-Yang 1994) with 2 parameters: omega and kappa */
int i,j;
IntVector group;
for (i = 0; i < num_states-1; i++) {
for (j = i+1; j < num_states; j++) {
if (isMultipleSubst(i, j))
group.push_back(0); // multiple substitution
else if (isSynonymous(i, j)) {
if (isTransversion(i, j))
group.push_back(1); // synonymous transversion
else
group.push_back(2); // synonymous transition
} else {
if (isTransversion(i, j))
group.push_back(3); // non-synonymous transversion
else
group.push_back(4); // non-synonymous transition
}
}
}
setRateGroup(group);
// set zero rate for multiple substitution
// 1 for synonymous transversion
// and kappa*omega for non-synonymous transition
setRateGroupConstraint("x0=0,x1=1,x4=x2*x3");
}
void ModelCodonParametric::writeInfo(ostream &out) {
double *variables = new double[getNDim()+1];
setVariables(variables);
if (name == "MG") {
out << "Nonsynonymous/synonymous ratio (omega): " << variables[1] << endl;
} else if (name == "GY") {
out << "Transition/transversion ratio (kappa): " << variables[1] << endl;
out << "Nonsynonymous/synonymous ratio (omega): " << variables[2] << endl;
}
delete [] variables;
}
|
