1LoadFunctionLibrary("../codon.bf");
2LoadFunctionLibrary("../DNA.bf");
3LoadFunctionLibrary("../parameters.bf");
4LoadFunctionLibrary("../frequencies.bf");
5LoadFunctionLibrary("../../UtilityFunctions.bf");
6LoadFunctionLibrary("MG_REV.bf");
7
8/** @module models.codon.MG_REV_TRIP */
9
10/**
11 * @name models.codon.MG_REV_TRIP.ModelDescription
12 * @param {String} type
13 * @param {String} code
14 */
15lfunction models.codon.MG_REV_TRIP.ModelDescription(type, code) {
16
17    // piggy-back on the standard MG_REV model for most of the code
18
19    mg_base = models.codon.MG_REV.ModelDescription (type, code);
20    mg_base[utility.getGlobalValue("terms.description")] = "The Muse-Gaut 94 codon-substitution model coupled with the general time reversible (GTR) model of nucleotide substitution, which allows for two-hit and three-hit substitutions";
21    mg_base[utility.getGlobalValue("terms.model.q_ij")] = "models.codon.MG_REV_TRIP._GenerateRate";
22
23    return mg_base;
24}
25
26
27lfunction models.codon.MG_REV_TRIP._GenerateRate(fromChar, toChar, namespace, model_type, model) {
28    return models.codon.MG_REV_TRIP._GenerateRate_generic (fromChar, toChar, namespace, model_type,
29    model[utility.getGlobalValue("terms.translation_table")],
30        "alpha", utility.getGlobalValue("terms.parameters.synonymous_rate"),
31        "beta", utility.getGlobalValue("terms.parameters.nonsynonymous_rate"),
32        "omega", utility.getGlobalValue("terms.parameters.omega_ratio"),
33        "delta", utility.getGlobalValue("terms.parameters.multiple_hit_rate"),
34        "psi", utility.getGlobalValue("terms.parameters.triple_hit_rate"),
35        "psi_syn", utility.getGlobalValue("terms.parameters.triple_hit_rate_syn")
36        );
37}
38
39/**
40 * @name models.codon.MG_REV_TRIP._GenerateRate
41 * @param {Number} fromChar
42 * @param {Number} toChar
43 * @param {String} namespace
44 * @param {String} model_type
45 * @param {Matrix} _tt - translation table
46 */
47
48
49lfunction models.codon.MG_REV_TRIP._GenerateRate_generic (fromChar, toChar, namespace, model_type, _tt, alpha, alpha_term, beta, beta_term, omega, omega_term, delta, delta_term, psi, psi_term, psi_s, psi_s_term) {
50
51    _GenerateRate.p = {};
52    _GenerateRate.diff = models.codon.diff.complete(fromChar, toChar);
53    diff_count = utility.Array1D (_GenerateRate.diff);
54
55    if (diff_count == 1 || diff_count == 2 || diff_count == 3) {
56
57        _GenerateRate.p[model_type] = {};
58        _GenerateRate.p[utility.getGlobalValue("terms.global")] = {};
59
60        nuc_rate = "";
61
62        for (i = 0; i < diff_count; i += 1) {
63            if ((_GenerateRate.diff[i])[utility.getGlobalValue("terms.diff.from")] > (_GenerateRate.diff[i])[utility.getGlobalValue("terms.diff.to")]) {
64                nuc_p = "theta_" + (_GenerateRate.diff[i])[utility.getGlobalValue("terms.diff.to")] + (_GenerateRate.diff[i])[utility.getGlobalValue("terms.diff.from")];
65            } else {
66                nuc_p = "theta_" + (_GenerateRate.diff[i])[utility.getGlobalValue("terms.diff.from")] +(_GenerateRate.diff[i])[utility.getGlobalValue("terms.diff.to")];
67            }
68            nuc_p = parameters.ApplyNameSpace(nuc_p, namespace);
69            (_GenerateRate.p[utility.getGlobalValue("terms.global")])[terms.nucleotideRateReversible((_GenerateRate.diff[i])[utility.getGlobalValue("terms.diff.from")], (_GenerateRate.diff[i])[utility.getGlobalValue("terms.diff.to")])] = nuc_p;
70
71            nuc_rate = parameters.AppendMultiplicativeTerm (nuc_rate, nuc_p);
72       }
73
74
75        rate_entry = nuc_rate;
76
77        if (_tt[fromChar] != _tt[toChar]) {
78            if (model_type == utility.getGlobalValue("terms.global")) {
79                aa_rate = parameters.ApplyNameSpace(omega, namespace);
80                (_GenerateRate.p[model_type])[omega_term] = aa_rate;
81            } else {
82                aa_rate = beta;
83                (_GenerateRate.p[model_type])[beta_term] = aa_rate;
84            }
85            rate_entry += "*" + aa_rate;
86        } else {
87            if (model_type == utility.getGlobalValue("terms.local")) {
88                (_GenerateRate.p[model_type])[alpha_term] = alpha;
89                rate_entry += "*" + alpha;
90            } else {
91                _GenerateRate.p[utility.getGlobalValue("terms.model.rate_entry")] = nuc_rate;
92            }
93
94
95        }
96
97        if (diff_count == 2) {
98            if (model_type == utility.getGlobalValue("terms.global")) {
99                delta_rate = parameters.ApplyNameSpace(delta, namespace);
100                (_GenerateRate.p[model_type])[delta_term] = delta_rate;
101                rate_entry += "*" + delta_rate;
102            } else {
103                (_GenerateRate.p[model_type])[delta_term] = delta;
104                rate_entry += "*" + delta;
105            }
106        }
107
108
109        if (diff_count == 3) {
110            if (_tt[fromChar] != _tt[toChar]) {
111                if (model_type == utility.getGlobalValue("terms.global")) {
112                    psi_rate = parameters.ApplyNameSpace(psi, namespace);
113                    (_GenerateRate.p[model_type])[psi_term] = psi_rate;
114                    rate_entry += "*" + psi_rate;
115                } else {
116                    (_GenerateRate.p[model_type])[psi_term] = psi;
117                    rate_entry += "*" + psi;
118                }
119            } else {
120               //console.log (fromChar + " <-> " + toChar + " (" + _tt[fromChar] + ")");
121               if (model_type == utility.getGlobalValue("terms.global")) {
122                    psi_rate = parameters.ApplyNameSpace(psi_s, namespace);
123                    (_GenerateRate.p[model_type])[psi_s_term] = psi_rate;
124                    rate_entry += "*" + psi_rate;
125                } else {
126                    (_GenerateRate.p[model_type])[psi_s_term] = psi_s;
127                    rate_entry += "*" + psi_s;
128                }
129
130            }
131        }
132
133
134        _GenerateRate.p[utility.getGlobalValue("terms.model.rate_entry")] = rate_entry;
135    }
136    return _GenerateRate.p;
137}
138