bpp-phyl  2.1.0
 All Classes Namespaces Files Functions Variables Friends Pages
TN93.cpp
Go to the documentation of this file.
1 //
2 // File: TN93.cpp
3 // Created by: Julien Dutheil
4 // Created on: Thu Jan 22 10:26:51 2004
5 //
6 
7 /*
8 Copyright or © or Copr. Bio++ Developement Team, (November 16, 2004)
9 
10 This software is a computer program whose purpose is to provide classes
11 for phylogenetic data analysis.
12 
13 This software is governed by the CeCILL license under French law and
14 abiding by the rules of distribution of free software. You can use,
15 modify and/ or redistribute the software under the terms of the CeCILL
16 license as circulated by CEA, CNRS and INRIA at the following URL
17 "http://www.cecill.info".
18 
19 As a counterpart to the access to the source code and rights to copy,
20 modify and redistribute granted by the license, users are provided only
21 with a limited warranty and the software's author, the holder of the
22 economic rights, and the successive licensors have only limited
23 liability.
24 
25 In this respect, the user's attention is drawn to the risks associated
26 with loading, using, modifying and/or developing or reproducing the
27 software by the user in light of its specific status of free software,
28 that may mean that it is complicated to manipulate, and that also
29 therefore means that it is reserved for developers and experienced
30 professionals having in-depth computer knowledge. Users are therefore
31 encouraged to load and test the software's suitability as regards their
32 requirements in conditions enabling the security of their systems and/or
33 data to be ensured and, more generally, to use and operate it in the
34 same conditions as regards security.
35 
36 The fact that you are presently reading this means that you have had
37 knowledge of the CeCILL license and that you accept its terms.
38 */
39 
40 #include "TN93.h"
41 #include "../FrequenciesSet/NucleotideFrequenciesSet.h"
42 
43 #include <Bpp/Numeric/Matrix/MatrixTools.h>
44 #include <Bpp/Numeric/Matrix/EigenValue.h>
45 
46 // From SeqLib:
47 #include <Bpp/Seq/Container/SequenceContainerTools.h>
48 
49 // From the STL:
50 #include <cmath>
51 
52 using namespace std;
53 
54 using namespace bpp;
55 
56 /******************************************************************************/
57 
58 TN93::TN93(
59  const NucleicAlphabet* alpha,
60  double kappa1,
61  double kappa2,
62  double piA,
63  double piC,
64  double piG,
65  double piT):
66  AbstractParameterAliasable("TN93."),
67  AbstractSubstitutionModel(alpha, "TN93."),
68  AbstractReversibleSubstitutionModel(alpha, "TN93."),
69  kappa1_(kappa1), kappa2_(kappa2),
70  piA_(piA), piC_(piC), piG_(piG), piT_(piT), piY_(), piR_(),
71  r_(), k1_(), k2_(),
72  theta_(piG + piC), theta1_(piA / (1. - theta_)), theta2_(piG / theta_),
73  exp1_(), exp21_(), exp22_(), l_(), p_(size_, size_)
74 {
75  addParameter_(new Parameter("TN93.kappa1", kappa1, &Parameter::R_PLUS_STAR));
76  addParameter_(new Parameter("TN93.kappa2", kappa2, &Parameter::R_PLUS_STAR));
77  addParameter_(new Parameter("TN93.theta" , theta_ , &FrequenciesSet::FREQUENCE_CONSTRAINT_SMALL));
78  addParameter_(new Parameter("TN93.theta1", theta1_, &FrequenciesSet::FREQUENCE_CONSTRAINT_SMALL));
79  addParameter_(new Parameter("TN93.theta2", theta2_, &FrequenciesSet::FREQUENCE_CONSTRAINT_SMALL));
80  p_.resize(size_, size_);
81  updateMatrices();
82 }
83 
84 /******************************************************************************/
85 
86 void TN93::updateMatrices()
87 {
88  kappa1_ = getParameterValue("kappa1");
89  kappa2_ = getParameterValue("kappa2");
90  theta_ = getParameterValue("theta" );
91  theta1_ = getParameterValue("theta1");
92  theta2_ = getParameterValue("theta2");
93  piA_ = theta1_ * (1. - theta_);
94  piC_ = (1. - theta2_) * theta_;
95  piG_ = theta2_ * theta_;
96  piT_ = (1. - theta1_) * (1. - theta_);
97  piR_ = piA_ + piG_;
98  piY_ = piT_ + piC_;
99  r_ = 1. / (2. * (piA_ * piC_ + piC_ * piG_ + piA_ * piT_ + piG_ * piT_ + kappa2_ * piC_ * piT_ + kappa1_ * piA_ * piG_));
100  k1_ = kappa2_ * piY_ + piR_;
101  k2_ = kappa1_ * piR_ + piY_;
102 
103  freq_[0] = piA_;
104  freq_[1] = piC_;
105  freq_[2] = piG_;
106  freq_[3] = piT_;
107 
108  generator_(0, 0) = -( piC_ + kappa1_*piG_ + piT_);
109  generator_(1, 1) = -( piA_ + piG_ + kappa2_*piT_);
110  generator_(2, 2) = -(kappa1_*piA_ + piC_ + piT_);
111  generator_(3, 3) = -( piA_ + kappa2_*piC_ + piG_ );
112 
113  generator_(1, 0) = piA_;
114  generator_(3, 0) = piA_;
115  generator_(0, 1) = piC_;
116  generator_(2, 1) = piC_;
117  generator_(1, 2) = piG_;
118  generator_(3, 2) = piG_;
119  generator_(0, 3) = piT_;
120  generator_(2, 3) = piT_;
121 
122  generator_(2, 0) = kappa1_ * piA_;
123  generator_(3, 1) = kappa2_ * piC_;
124  generator_(0, 2) = kappa1_ * piG_;
125  generator_(1, 3) = kappa2_ * piT_;
126 
127  // Normalization:
128  MatrixTools::scale(generator_, r_);
129 
130  // Exchangeability:
131  exchangeability_(0,0) = generator_(0,0) / piA_;
132  exchangeability_(0,1) = generator_(0,1) / piC_;
133  exchangeability_(0,2) = generator_(0,2) / piG_;
134  exchangeability_(0,3) = generator_(0,3) / piT_;
135 
136  exchangeability_(1,0) = generator_(1,0) / piA_;
137  exchangeability_(1,1) = generator_(1,1) / piC_;
138  exchangeability_(1,2) = generator_(1,2) / piG_;
139  exchangeability_(1,3) = generator_(1,3) / piT_;
140 
141  exchangeability_(2,0) = generator_(2,0) / piA_;
142  exchangeability_(2,1) = generator_(2,1) / piC_;
143  exchangeability_(2,2) = generator_(2,2) / piG_;
144  exchangeability_(2,3) = generator_(2,3) / piT_;
145 
146  exchangeability_(3,0) = generator_(3,0) / piA_;
147  exchangeability_(3,1) = generator_(3,1) / piC_;
148  exchangeability_(3,2) = generator_(3,2) / piG_;
149  exchangeability_(3,3) = generator_(3,3) / piT_;
150 
151  // We are not sure that the values are computed in this order :p
152  // Eigen values:
153  //eigenValues_[0] = 0;
154  //eigenValues_[1] = -r * (kappa2 * piY + piR);
155  //eigenValues_[2] = -r * (kappa1 * piR + piY);
156  //eigenValues_[3] = -r;
157 
158  // Eigen vectors and values:
159  EigenValue<double> ev(generator_);
160  rightEigenVectors_ = ev.getV();
161  MatrixTools::inv(rightEigenVectors_, leftEigenVectors_);
162  eigenValues_ = ev.getRealEigenValues();
163 }
164 
165 /******************************************************************************/
166 
167 double TN93::Pij_t(int i, int j, double d) const
168 {
169  l_ = rate_ * r_ * d;
170  exp1_ = exp(-l_);
171  exp22_ = exp(-k2_ * l_);
172  exp21_ = exp(-k1_ * l_);
173 
174  switch(i)
175  {
176  //A
177  case 0 : {
178  switch(j) {
179  case 0 : return piA_ * (1. + (piY_/piR_) * exp1_) + (piG_/piR_) * exp22_; //A
180  case 1 : return piC_ * (1. - exp1_); //C
181  case 2 : return piG_ * (1. + (piY_/piR_) * exp1_) - (piG_/piR_) * exp22_; //G
182  case 3 : return piT_ * (1. - exp1_); //T, U
183  }
184  }
185  //C
186  case 1 : {
187  switch(j) {
188  case 0 : return piA_ * (1. - exp1_); //A
189  case 1 : return piC_ * (1. + (piR_/piY_) * exp1_) + (piT_/piY_) * exp21_; //C
190  case 2 : return piG_ * (1. - exp1_); //G
191  case 3 : return piT_ * (1. + (piR_/piY_) * exp1_) - (piT_/piY_) * exp21_; //T, U
192  }
193  }
194  //G
195  case 2 : {
196  switch(j) {
197  case 0 : return piA_ * (1. + (piY_/piR_) * exp1_) - (piA_/piR_) * exp22_; //A
198  case 1 : return piC_ * (1. - exp1_); //C
199  case 2 : return piG_ * (1. + (piY_/piR_) * exp1_) + (piA_/piR_) * exp22_; //G
200  case 3 : return piT_ * (1. - exp1_); //T, U
201  }
202  }
203  //T, U
204  case 3 : {
205  switch(j) {
206  case 0 : return piA_ * (1. - exp1_); //A
207  case 1 : return piC_ * (1. + (piR_/piY_) * exp1_) - (piC_/piY_) * exp21_; //C
208  case 2 : return piG_ * (1. - exp1_); //G
209  case 3 : return piT_ * (1. + (piR_/piY_) * exp1_) + (piC_/piY_) * exp21_; //T, U
210  }
211  }
212  }
213  return 0;
214 }
215 
216 /******************************************************************************/
217 
218 double TN93::dPij_dt(int i, int j, double d) const
219 {
220  l_ = rate_ * r_ * d;
221  exp1_ = exp(-l_);
222  exp22_ = exp(-k2_ * l_);
223  exp21_ = exp(-k1_ * l_);
224 
225  switch(i)
226  {
227  //A
228  case 0 : {
229  switch(j) {
230  case 0 : return rate_ * r_ * (piA_ * -(piY_/piR_) * exp1_ - (piG_/piR_) * k2_ * exp22_); //A
231  case 1 : return rate_ * r_ * (piC_ * exp1_); //C
232  case 2 : return rate_ * r_ * (piG_ * -(piY_/piR_) * exp1_ + (piG_/piR_) * k2_ * exp22_); //G
233  case 3 : return rate_ * r_ * (piT_ * exp1_); //T, U
234  }
235  }
236  //C
237  case 1 : {
238  switch(j) {
239  case 0 : return rate_ * r_ * (piA_ * exp1_); //A
240  case 1 : return rate_ * r_ * (piC_ * -(piR_/piY_) * exp1_ - (piT_/piY_) * k1_ * exp21_); //C
241  case 2 : return rate_ * r_ * (piG_ * exp1_); //G
242  case 3 : return rate_ * r_ * (piT_ * -(piR_/piY_) * exp1_ + (piT_/piY_) * k1_ * exp21_); //T, U
243  }
244  }
245  //G
246  case 2 : {
247  switch(j) {
248  case 0 : return rate_ * r_ * (piA_ * -(piY_/piR_) * exp1_ + (piA_/piR_) * k2_ * exp22_); //A
249  case 1 : return rate_ * r_ * (piC_ * exp1_); //C
250  case 2 : return rate_ * r_ * (piG_ * -(piY_/piR_) * exp1_ - (piA_/piR_) * k2_ * exp22_); //G
251  case 3 : return rate_ * r_ * (piT_ * exp1_); //T, U
252  }
253  }
254  //T, U
255  case 3 : {
256  switch(j) {
257  case 0 : return rate_ * r_ * (piA_ * exp1_); //A
258  case 1 : return rate_ * r_ * (piC_ * -(piR_/piY_) * exp1_ + (piC_/piY_) * k1_ * exp21_); //C
259  case 2 : return rate_ * r_ * (piG_ * exp1_); //G
260  case 3 : return rate_ * r_ * (piT_ * -(piR_/piY_) * exp1_ - (piC_/piY_) * k1_ * exp21_); //T, U
261  }
262  }
263  }
264  return 0;
265 }
266 
267 /******************************************************************************/
268 
269 double TN93::d2Pij_dt2(int i, int j, double d) const
270 {
271  double r_2 = rate_ * rate_ * r_ * r_;
272  l_ = rate_ * r_ * d;
273  double k1_2 = k1_ * k1_;
274  double k2_2 = k2_ * k2_;
275  exp1_ = exp(-l_);
276  exp22_ = exp(-k2_ * l_);
277  exp21_ = exp(-k1_ * l_);
278 
279  switch(i)
280  {
281  //A
282  case 0 : {
283  switch(j) {
284  case 0 : return r_2 * (piA_ * (piY_/piR_) * exp1_ + (piG_/piR_) * k2_2 * exp22_); //A
285  case 1 : return r_2 * (piC_ * - exp1_); //C
286  case 2 : return r_2 * (piG_ * (piY_/piR_) * exp1_ - (piG_/piR_) * k2_2 * exp22_); //G
287  case 3 : return r_2 * (piT_ * - exp1_); //T, U
288  }
289  }
290  //C
291  case 1 : {
292  switch(j) {
293  case 0 : return r_2 * (piA_ * - exp1_); //A
294  case 1 : return r_2 * (piC_ * (piR_/piY_) * exp1_ + (piT_/piY_) * k1_2 * exp21_); //C
295  case 2 : return r_2 * (piG_ * - exp1_); //G
296  case 3 : return r_2 * (piT_ * (piR_/piY_) * exp1_ - (piT_/piY_) * k1_2 * exp21_); //T, U
297  }
298  }
299  //G
300  case 2 : {
301  switch(j) {
302  case 0 : return r_2 * (piA_ * (piY_/piR_) * exp1_ - (piA_/piR_) * k2_2 * exp22_); //A
303  case 1 : return r_2 * (piC_ * - exp1_); //C
304  case 2 : return r_2 * (piG_ * (piY_/piR_) * exp1_ + (piA_/piR_) * k2_2 * exp22_); //G
305  case 3 : return r_2 * (piT_ * - exp1_); //T, U
306  }
307  }
308  //T, U
309  case 3 : {
310  switch(j) {
311  case 0 : return r_2 * (piA_ * - exp1_); //A
312  case 1 : return r_2 * (piC_ * (piR_/piY_) * exp1_ - (piC_/piY_) * k1_2 * exp21_); //C
313  case 2 : return r_2 * (piG_ * - exp1_); //G
314  case 3 : return r_2 * (piT_ * (piR_/piY_) * exp1_ + (piC_/piY_) * k1_2 * exp21_); //T, U
315  }
316  }
317  }
318  return 0;
319 }
320 
321 /******************************************************************************/
322 
323 const Matrix<double> & TN93::getPij_t(double d) const
324 {
325  l_ = rate_ * r_ * d;
326  exp1_ = exp(-l_);
327  exp22_ = exp(-k2_ * l_);
328  exp21_ = exp(-k1_ * l_);
329 
330  //A
331  p_(0, 0) = piA_ * (1. + (piY_/piR_) * exp1_) + (piG_/piR_) * exp22_; //A
332  p_(0, 1) = piC_ * (1. - exp1_); //C
333  p_(0, 2) = piG_ * (1. + (piY_/piR_) * exp1_) - (piG_/piR_) * exp22_; //G
334  p_(0, 3) = piT_ * (1. - exp1_); //T, U
335 
336  //C
337  p_(1, 0) = piA_ * (1. - exp1_); //A
338  p_(1, 1) = piC_ * (1. + (piR_/piY_) * exp1_) + (piT_/piY_) * exp21_; //C
339  p_(1, 2) = piG_ * (1. - exp1_); //G
340  p_(1, 3) = piT_ * (1. + (piR_/piY_) * exp1_) - (piT_/piY_) * exp21_; //T, U
341 
342  //G
343  p_(2, 0) = piA_ * (1. + (piY_/piR_) * exp1_) - (piA_/piR_) * exp22_; //A
344  p_(2, 1) = piC_ * (1. - exp1_); //C
345  p_(2, 2) = piG_ * (1. + (piY_/piR_) * exp1_) + (piA_/piR_) * exp22_; //G
346  p_(2, 3) = piT_ * (1. - exp1_); //T, U
347 
348  //T, U
349  p_(3, 0) = piA_ * (1. - exp1_); //A
350  p_(3, 1) = piC_ * (1. + (piR_/piY_) * exp1_) - (piC_/piY_) * exp21_; //C
351  p_(3, 2) = piG_ * (1. - exp1_); //G
352  p_(3, 3) = piT_ * (1. + (piR_/piY_) * exp1_) + (piC_/piY_) * exp21_; //T, U
353 
354  return p_;
355 }
356 
357 const Matrix<double> & TN93::getdPij_dt(double d) const
358 {
359  l_ = rate_ * r_ * d;
360  exp1_ = exp(-l_);
361  exp22_ = exp(-k2_ * l_);
362  exp21_ = exp(-k1_ * l_);
363 
364  //A
365  p_(0, 0) = rate_ * r_ * (piA_ * -(piY_/piR_) * exp1_ - (piG_/piR_) * k2_ * exp22_); //A
366  p_(0, 1) = rate_ * r_ * (piC_ * exp1_); //C
367  p_(0, 2) = rate_ * r_ * (piG_ * -(piY_/piR_) * exp1_ + (piG_/piR_) * k2_ * exp22_); //G
368  p_(0, 3) = rate_ * r_ * (piT_ * exp1_); //T, U
369 
370  //C
371  p_(1, 0) = rate_ * r_ * (piA_ * exp1_); //A
372  p_(1, 1) = rate_ * r_ * (piC_ * -(piR_/piY_) * exp1_ - (piT_/piY_) * k1_ * exp21_); //C
373  p_(1, 2) = rate_ * r_ * (piG_ * exp1_); //G
374  p_(1, 3) = rate_ * r_ * (piT_ * -(piR_/piY_) * exp1_ + (piT_/piY_) * k1_ * exp21_); //T, U
375 
376  //G
377  p_(2, 0) = rate_ * r_ * (piA_ * -(piY_/piR_) * exp1_ + (piA_/piR_) * k2_ * exp22_); //A
378  p_(2, 1) = rate_ * r_ * (piC_ * exp1_); //C
379  p_(2, 2) = rate_ * r_ * (piG_ * -(piY_/piR_) * exp1_ - (piA_/piR_) * k2_ * exp22_); //G
380  p_(2, 3) = rate_ * r_ * (piT_ * exp1_); //T, U
381 
382  //T, U
383  p_(3, 0) = rate_ * r_ * (piA_ * exp1_); //A
384  p_(3, 1) = rate_ * r_ * (piC_ * -(piR_/piY_) * exp1_ + (piC_/piY_) * k1_ * exp21_); //C
385  p_(3, 2) = rate_ * r_ * (piG_ * exp1_); //G
386  p_(3, 3) = rate_ * r_ * (piT_ * -(piR_/piY_) * exp1_ - (piC_/piY_) * k1_ * exp21_); //T, U
387 
388  return p_;
389 }
390 
391 const Matrix<double> & TN93::getd2Pij_dt2(double d) const
392 {
393  double r_2 = rate_ * rate_ * r_ * r_;
394  l_ = rate_ * r_ * d;
395  double k1_2 = k1_ * k1_;
396  double k2_2 = k2_ * k2_;
397  exp1_ = exp(-l_);
398  exp22_ = exp(-k2_ * l_);
399  exp21_ = exp(-k1_ * l_);
400 
401  //A
402  p_(0, 0) = r_2 * (piA_ * (piY_/piR_) * exp1_ + (piG_/piR_) * k2_2 * exp22_); //A
403  p_(0, 1) = r_2 * (piC_ * - exp1_); //C
404  p_(0, 2) = r_2 * (piG_ * (piY_/piR_) * exp1_ - (piG_/piR_) * k2_2 * exp22_); //G
405  p_(0, 3) = r_2 * (piT_ * - exp1_); //T, U
406 
407  //C
408  p_(1, 0) = r_2 * (piA_ * - exp1_); //A
409  p_(1, 1) = r_2 * (piC_ * (piR_/piY_) * exp1_ + (piT_/piY_) * k1_2 * exp21_); //C
410  p_(1, 2) = r_2 * (piG_ * - exp1_); //G
411  p_(1, 3) = r_2 * (piT_ * (piR_/piY_) * exp1_ - (piT_/piY_) * k1_2 * exp21_); //T, U
412 
413  //G
414  p_(2, 0) = r_2 * (piA_ * (piY_/piR_) * exp1_ - (piA_/piR_) * k2_2 * exp22_); //A
415  p_(2, 1) = r_2 * (piC_ * - exp1_); //C
416  p_(2, 2) = r_2 * (piG_ * (piY_/piR_) * exp1_ + (piA_/piR_) * k2_2 * exp22_); //G
417  p_(2, 3) = r_2 * (piT_ * - exp1_); //T, U
418 
419  //T, U
420  p_(3, 0) = r_2 * (piA_ * - exp1_); //A
421  p_(3, 1) = r_2 * (piC_ * (piR_/piY_) * exp1_ - (piC_/piY_) * k1_2 * exp21_); //C
422  p_(3, 2) = r_2 * (piG_ * - exp1_); //G
423  p_(3, 3) = r_2 * (piT_ * (piR_/piY_) * exp1_ + (piC_/piY_) * k1_2 * exp21_); //T, U
424 
425  return p_;
426 }
427 
428 /******************************************************************************/
429 
430 void TN93::setFreq(std::map<int, double>& freqs)
431 {
432  piA_ = freqs[0];
433  piC_ = freqs[1];
434  piG_ = freqs[2];
435  piT_ = freqs[3];
436  vector<string> thetas(3);
437  thetas[0] = getNamespace() + "theta";
438  thetas[1] = getNamespace() + "theta1";
439  thetas[2] = getNamespace() + "theta2";
440  ParameterList pl = getParameters().subList(thetas);
441  pl[0].setValue(piC_ + piG_);
442  pl[1].setValue(piA_ / (piA_ + piT_));
443  pl[2].setValue(piG_ / (piC_ + piG_));
444  setParametersValues(pl);
445 }
446 
447 /******************************************************************************/
448