1 /*
2 * modelcodon.cpp
3 *
4 * Created on: May 24, 2013
5 * Author: minh
6 */
7
8 #include "modelcodon.h"
9 #include <string>
10
11
12 const double MIN_OMEGA_KAPPA = 0.001;
13 const double MAX_OMEGA_KAPPA = 50.0;
14
15 /* Empirical codon model restricted (Kosiol et al. 2007), source: http://www.ebi.ac.uk/goldman/ECM/ */
16 string model_ECMrest1 =
17 "11.192024 \
18 1.315610 0.010896 \
19 5.427076 4.756288 24.748755 \
20 1.658051 0.000000 0.000000 0.000000 \
21 0.000000 1.913571 0.000000 0.000000 13.889102 \
22 0.000000 0.000000 2.952332 0.000000 44.407955 13.681751 \
23 0.000000 0.000000 0.000000 8.126914 17.057443 65.097021 12.991861 \
24 6.610894 0.000000 0.000000 0.000000 2.206054 0.000000 0.000000 0.000000 \
25 0.000000 5.177930 0.000000 0.000000 0.000000 5.615472 0.000000 0.000000 19.942818 \
26 3.347364 0.000000 0.000000 0.000000 6.191481 0.000000 0.000000 0.000000 0.582084 0.000000 \
27 0.000000 1.558523 0.000000 0.000000 0.000000 9.339206 0.000000 0.000000 0.000000 0.144278 44.777964 \
28 0.000000 0.000000 0.000000 5.369644 0.000000 0.000000 0.000000 4.662001 0.000000 0.000000 0.677177 0.073268 \
29 2.090751 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 \
30 0.000000 2.266373 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 8.905484 \
31 0.000000 0.000000 75.752638 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 56.803876 7.811205 \
32 0.000000 0.000000 0.000000 20.877218 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 8.432339 22.078564 5.650116 \
33 0.000000 0.000000 0.000000 0.000000 1.769355 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.263838 0.000000 0.000000 0.000000 \
34 0.000000 0.000000 0.000000 0.000000 0.000000 2.704601 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.389735 0.000000 0.000000 17.461627 \
35 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 3.312811 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.393680 0.000000 35.480963 12.053827 \
36 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.303480 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.477616 8.407091 28.557939 11.295213 \
37 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 3.444964 0.000000 0.000000 0.000000 0.000000 1.583116 0.000000 0.000000 0.000000 1.021682 0.000000 0.000000 0.000000 \
38 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 7.087801 0.000000 0.000000 0.000000 0.000000 3.230751 0.000000 0.000000 0.000000 3.774544 0.000000 0.000000 28.086160 \
39 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 7.419058 0.000000 0.000000 0.000000 5.381868 0.000000 3.440380 1.918904 \
40 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.812540 0.000000 0.000000 0.000000 1.794388 1.086327 5.369463 14.959151 \
41 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.617091 0.000000 0.000000 0.779565 0.000000 0.000000 0.000000 0.334165 0.000000 0.000000 0.000000 3.019726 0.000000 0.000000 0.000000 \
42 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.632945 0.000000 0.000000 2.250770 0.000000 0.000000 0.000000 1.699302 0.000000 0.000000 0.000000 7.016899 0.000000 0.000000 14.603857 \
43 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 3.023939 0.000000 0.000000 0.000000 1.693662 0.000000 0.000000 0.000000 6.415757 0.000000 99.459951 14.930266 \
44 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 3.026086 0.000000 0.000000 0.000000 1.462945 0.000000 0.000000 0.000000 3.144296 0.000000 0.000000 0.000000 19.920977 30.804750 79.483730 13.919752 \
45 1.682029 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 4.301225 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 \
46 0.000000 0.786043 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 6.381841 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 10.140728 \
47 0.000000 0.000000 10.116588 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 5.134459 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 18.298900 4.623936 \
48 0.000000 0.000000 0.000000 7.911096 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.570123 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.281784 1.303951 2.082128 \
49 0.000000 0.000000 0.000000 0.000000 38.229100 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 6.578976 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.801564 0.000000 0.000000 0.000000 \
50 0.000000 0.000000 0.000000 0.000000 0.000000 15.793595 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.434550 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.231468 0.000000 0.000000 6.035740 \
51 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 6.033932 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.925575 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 4.962350 0.000000 28.307876 6.967655 \
52 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 17.103904 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.238450 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 8.155285 19.578982 38.414969 12.678802 \
53 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.245405 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 5.004762 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.501054 0.000000 0.000000 0.000000 11.715476 0.000000 0.000000 0.000000 \
54 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.228361 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 4.105602 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.292691 0.000000 0.000000 0.000000 2.134740 0.000000 0.000000 13.863648 \
55 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 6.404436 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.647620 0.000000 0.000000 0.000000 3.919360 0.000000 4.929483 0.366267 \
56 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.715692 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.975074 0.000000 0.000000 0.000000 5.869857 1.010212 0.982893 10.762877 \
57 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 4.719489 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 3.834666 0.000000 0.000000 0.000000 0.578118 0.000000 0.000000 0.000000 39.399322 0.000000 0.000000 0.000000 16.623529 0.000000 0.000000 0.000000 \
58 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.047654 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 5.033630 0.000000 0.000000 0.000000 0.437779 0.000000 0.000000 0.000000 21.337943 0.000000 0.000000 0.000000 7.784768 0.000000 0.000000 26.637668 \
59 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 92.372238 0.000000 0.000000 0.000000 1.903175 0.000000 0.000000 0.000000 0.754055 0.000000 0.000000 0.000000 8.423762 0.000000 1.792245 0.120900 \
60 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.825082 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 133.296291 0.000000 0.000000 0.000000 2.231662 0.000000 0.000000 0.000000 22.577271 0.000000 0.000000 0.000000 21.000358 3.324581 6.011970 36.292705 \
61 2.261813 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.473623 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 7.096281 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 \
62 0.000000 1.923392 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 5.914972 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 10.137337 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 6.669955 \
63 0.000000 0.000000 2.362720 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 3.737489 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 25.294298 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 26.045078 3.531461 \
64 0.000000 0.000000 0.000000 2.022101 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.164805 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.078444 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 8.901167 21.657664 11.898141 \
65 0.000000 0.000000 0.000000 0.000000 5.540052 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.159185 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 5.107629 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 3.682092 0.000000 0.000000 0.000000 \
66 0.000000 0.000000 0.000000 0.000000 0.000000 7.675838 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 3.120189 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.312255 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 4.308415 0.000000 0.000000 6.516319 \
67 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 9.880382 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.923972 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 3.064069 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 6.291148 0.000000 21.910225 5.090423 \
68 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 21.863158 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 6.034856 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 25.461549 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 6.166554 5.512586 20.715347 9.529141 \
69 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.367553 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.383706 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 6.091654 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.352915 0.000000 0.000000 0.000000 0.693026 0.000000 0.000000 0.000000 \
70 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.294702 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 3.006827 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 3.686074 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.208522 0.000000 0.000000 0.000000 1.866565 0.000000 0.000000 10.605899 \
71 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 4.485369 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.811398 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.277861 0.000000 0.000000 0.000000 2.774445 0.000000 2.710610 0.650088 \
72 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 7.686782 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.090641 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.476105 0.000000 0.000000 0.000000 9.441919 1.296294 3.779053 10.153570 \
73 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.104727 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.041150 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 10.590780 0.000000 0.000000 0.000000 0.503385 0.000000 0.000000 0.000000 1.541379 0.000000 0.000000 0.000000 1.042624 0.000000 0.000000 0.000000 \
74 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.552851 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.252470 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 4.285543 0.000000 0.000000 0.000000 0.542717 0.000000 0.000000 0.000000 2.303487 0.000000 0.000000 0.000000 1.561629 0.000000 0.000000 9.488520";
75 string model_ECMrest = model_ECMrest1 +
76 "0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.091041 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.432410 0.000000 0.000000 0.000000 0.702411 0.000000 0.000000 0.000000 2.985093 0.000000 0.000000 0.000000 0.874000 0.000000 20.518100 4.120953 \
77 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.810856 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 4.803738 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 5.388514 0.000000 0.000000 0.000000 0.302501 0.000000 0.000000 0.000000 6.644971 0.000000 0.000000 0.000000 1.393810 13.246936 18.064826 19.084271 \
78 \
79 0.022103 0.021383 0.016387 0.015425 0.011880 0.011131 0.009750 0.008956 0.015965 0.015782 0.006025 0.007029 0.011880 0.014467 0.017386 0.007600 0.028839 0.010007 0.010100 0.010642 0.011843 0.011097 0.011703 0.016076 0.020211 0.008311 0.014148 0.004800 0.007837 0.025576 0.023441 0.013551 0.020102 0.013424 0.020201 0.015528 0.012142 0.023006 0.020171 0.030001 0.026344 0.010142 0.011679 0.010372 0.008195 0.019047 0.018938 0.010901 0.022747 0.019005 0.028307 0.015908 0.018853 0.028198 0.024532 0.033223 0.031878 0.016852 0.022982 0.015796 0.010191 \
80 \
81 TTT TTC TTA TTG TCT TCC TCA TCG TAT TAC TGT TGC TGG CTT CTC CTA CTG CCT CCC CCA \
82 CCG CAT CAC CAA CAG CGT CGC CGA CGG ATT ATC ATA ATG ACT ACC ACA ACG AAT AAC AAA \
83 AAG AGT AGC AGA AGG GTT GTC GTA GTG GCT GCC GCA GCG GAT GAC GAA GAG GGT GGC GGA \
84 GGG";
85
86 /* Empirical codon model unrestricted (Kosiol et al. 2007), source: http://www.ebi.ac.uk/goldman/ECM/ */
87 string model_ECMunrest1 =
88 "16.011531 \
89 2.395822 0.151858 \
90 1.204356 0.675537 18.541946 \
91 0.773935 0.052602 0.249707 0.274990 \
92 0.030074 0.656004 0.011609 0.158873 23.655090 \
93 0.278090 0.056677 1.184813 0.611887 35.921779 15.982573 \
94 0.034137 0.198277 0.010188 0.694091 11.510965 35.359077 17.424222 \
95 4.317981 0.503397 0.798582 0.337279 0.688169 0.047115 0.341791 0.058136 \
96 0.481042 4.483501 0.033529 0.177833 0.069588 0.524116 0.070809 0.213967 24.177765 \
97 0.733587 0.076912 0.645571 0.395942 1.811753 0.343463 0.751980 0.143447 0.822999 0.054860 \
98 0.045951 0.561620 0.040012 0.240632 0.138244 1.323765 0.121937 0.493179 0.068342 0.628438 56.838378 \
99 0.786871 1.183337 0.271072 0.632947 0.069758 0.081312 0.195833 0.410046 1.140051 1.421996 0.264556 0.210115 \
100 2.016257 0.207692 12.035723 11.161511 0.277929 0.000186 0.000289 0.000000 0.485469 0.000299 0.543240 0.000674 0.010122 \
101 0.083684 2.306110 1.373823 5.651603 0.000085 0.342813 0.000096 0.000344 0.000116 0.622089 0.000466 0.674176 0.113701 15.874441 \
102 1.036474 0.198558 27.219895 16.560966 0.000678 0.000186 0.496046 0.000115 0.016650 0.011978 0.020649 0.021578 0.017106 21.437257 8.808275 \
103 0.073550 1.341144 1.045943 12.455337 0.000000 0.001022 0.000000 0.266943 0.004815 0.308859 0.002639 0.265948 0.504866 4.802017 15.484088 8.319767 \
104 0.324368 0.000141 0.001358 0.003499 2.846677 0.196358 0.544474 0.078776 0.337879 0.000479 0.239715 0.000270 0.061833 0.822643 0.036254 0.181411 0.014388 \
105 0.000140 0.285635 0.000000 0.000382 0.101204 2.487136 0.072352 0.432520 0.000116 0.310416 0.000000 0.215779 0.032564 0.026571 0.648769 0.040087 0.149771 23.496083 \
106 0.025217 0.006558 0.261069 0.005535 0.487542 0.138742 3.121656 0.151589 0.032140 0.025873 0.002795 0.010250 0.070308 0.065669 0.016609 1.073790 0.040917 40.922701 15.426733 \
107 0.004063 0.079161 0.000000 0.112999 0.021444 0.371063 0.064924 2.075226 0.004177 0.037372 0.000155 0.004585 0.215788 0.007978 0.118229 0.016442 0.495176 10.291826 33.453780 15.127582 \
108 0.638696 0.001312 0.026551 0.040275 1.253945 0.002137 0.128111 0.073730 3.088481 0.340541 0.634065 0.001483 0.195073 0.664866 0.057328 0.438648 0.044742 0.775254 0.091276 0.286252 0.054021 \
109 0.000467 0.761771 0.000123 0.002163 0.000593 1.144692 0.014470 0.114551 0.265766 3.193996 0.000155 0.483076 0.369273 0.058614 0.617694 0.059927 0.330036 0.061583 0.730306 0.089835 0.364129 38.685701 \
110 0.126320 0.016628 0.576476 0.007508 0.508308 0.080383 2.066955 0.002179 0.486281 0.079236 0.163174 0.032232 0.055163 0.529045 0.071794 1.205738 0.033372 0.435109 0.074846 1.052040 0.063366 2.473439 0.751904 \
111 0.009760 0.107218 0.000000 0.250748 0.049246 0.423382 0.002122 1.519211 0.092070 0.332396 0.057910 0.105597 0.247490 0.079119 0.422671 0.105449 0.703795 0.107434 0.529594 0.184327 0.715716 1.106179 2.503268 17.923045 \
112 0.143832 0.000094 0.000741 0.003054 0.660622 0.001208 0.000579 0.001720 0.534375 0.001377 0.726908 0.077815 0.019696 0.663877 0.068758 0.134394 0.015019 0.500433 0.124232 0.063413 0.044676 2.460976 0.277265 1.164262 0.340811 \
113 0.000000 0.200806 0.000000 0.000064 0.000000 0.685812 0.000000 0.032106 0.000116 0.604541 0.012886 0.516927 0.176476 0.016022 0.544828 0.005436 0.563956 0.002398 0.563799 0.001702 0.798346 0.170088 2.478358 0.148940 2.029914 27.244097 \
114 0.030121 0.016020 0.136647 0.001527 0.006103 0.004089 0.557015 0.003211 0.043917 0.051686 0.232728 0.166150 0.146501 0.424607 0.112395 0.918198 0.041969 0.352807 0.154017 0.626603 0.091073 1.353860 0.526904 4.725840 0.617320 39.595443 12.677657 \
115 0.000934 0.027355 0.000000 0.127696 0.000085 0.004832 0.000000 1.903571 0.003713 0.081931 0.023909 0.183143 1.135910 0.039428 0.640495 0.040902 0.794366 0.009880 0.897101 0.010300 1.164525 0.316372 2.208430 0.299978 4.718199 12.868484 35.563093 30.574631 \
116 1.119411 0.059956 2.130663 1.292935 0.172403 0.000000 0.000386 0.000000 0.352731 0.000180 0.431456 0.000405 0.078312 3.330793 0.184010 1.328581 0.089308 0.292855 0.000096 0.002597 0.000246 0.193328 0.000000 0.078926 0.003859 0.076434 0.000000 0.000416 0.000000 \
117 0.056038 1.006045 0.042112 0.478019 0.000000 0.115975 0.000096 0.000344 0.000116 0.255975 0.000311 0.309643 0.136849 0.390190 3.765697 0.203017 2.469249 0.000096 0.270274 0.000448 0.021723 0.000469 0.127899 0.010543 0.105885 0.000118 0.238839 0.001248 0.003064 13.609310 \
118 1.075187 0.064968 5.159075 1.065537 0.000424 0.000000 0.403435 0.000000 0.573013 0.025454 0.069555 0.012138 0.170041 1.260239 0.136148 4.400610 0.048882 0.002014 0.000000 0.480521 0.000000 0.040109 0.000272 0.390087 0.000048 0.000000 0.000000 0.121855 0.000000 16.415611 5.784672 \
119 0.679370 0.800602 1.418466 3.062807 0.093491 0.042282 0.246094 0.527005 0.294368 0.300354 0.298091 0.324613 0.321642 1.220020 1.434579 1.635281 2.236557 0.081631 0.008455 0.042006 0.193459 0.323588 0.163406 0.443617 0.834976 0.028736 0.029786 0.015596 0.408680 1.155098 1.428293 2.230691 \
120 0.497293 0.000141 0.012473 0.015652 4.693944 0.487317 2.297807 0.199748 0.599932 0.000599 1.089585 0.001483 0.035939 0.831215 0.000060 0.004348 0.000070 1.050363 0.053805 0.345545 0.011476 0.898794 0.000272 0.374419 0.029088 0.344601 0.000000 0.001040 0.000000 1.266654 0.075878 0.351882 0.419831 \
121 0.000093 0.371541 0.000062 0.002990 0.196983 3.829580 0.150107 2.395833 0.000116 0.393545 0.000776 0.806071 0.037822 0.000396 0.897490 0.000679 0.073657 0.044125 0.870967 0.027138 0.527094 0.001125 0.969387 0.066519 0.617464 0.001060 1.481721 0.002079 0.017774 0.034573 1.066285 0.016701 0.433759 13.991583 \
122 0.079948 0.010539 0.871568 0.011134 2.018483 0.409721 5.709933 0.349846 0.208041 0.053363 0.415775 0.061227 0.060421 0.000857 0.000119 0.944560 0.000070 0.368538 0.026230 1.018005 0.015001 0.082373 0.021920 1.660885 0.001302 0.000589 0.000000 0.360575 0.000000 0.296014 0.048902 2.260424 0.853779 31.915858 8.373639 \
123 0.008032 0.036489 0.000062 0.586818 0.361164 1.562591 0.516594 4.919174 0.042119 0.177757 0.053874 0.173298 0.362681 0.000264 0.000595 0.000408 0.733202 0.079712 0.435243 0.083475 0.962295 0.076938 0.103080 0.002854 1.766961 0.004593 0.014243 0.002911 2.985421 0.090674 0.311759 0.154441 1.376727 12.116657 28.470047 19.459275 \
124 0.263567 0.000094 0.271628 0.077878 1.773102 0.000929 0.872084 0.040706 0.747870 0.042762 0.360038 0.000135 0.074859 0.259380 0.000000 0.019568 0.000140 0.787340 0.000192 0.096104 0.002705 2.691226 0.188587 1.759732 0.206851 0.682254 0.000068 0.009981 0.000981 0.239388 0.014351 0.256283 0.208924 2.057449 0.067554 1.243753 0.224397 \
125 0.000093 0.143614 0.002840 0.060699 0.003390 1.118208 0.187826 0.725836 0.046586 0.487814 0.000932 0.325422 0.053908 0.000330 0.187880 0.014540 0.034916 0.000959 0.532092 0.074161 0.095418 0.460970 2.203539 0.377447 0.985145 0.003180 0.863191 0.016636 0.065212 0.025405 0.175491 0.020837 0.219170 0.296066 1.346385 0.259909 0.822133 17.634677 \
126 0.148268 0.005996 0.612660 0.004963 0.340991 0.020909 1.496628 0.000459 0.467136 0.042942 0.099519 0.004316 0.051005 0.060922 0.002143 0.433620 0.000035 0.247195 0.012394 1.042457 0.000410 0.899262 0.143479 3.215817 0.285384 1.769879 0.296358 3.065510 0.011032 0.221536 0.023020 0.667397 0.275355 0.878163 0.089476 1.523251 0.199589 2.075154 0.413957 \
127 0.013122 0.043609 0.000062 0.333780 0.156468 0.251650 0.004438 0.768607 0.072867 0.140864 0.011489 0.034794 0.105226 0.039362 0.022384 0.000679 0.133032 0.220816 0.303613 0.012002 0.561523 0.324525 0.571469 0.461383 2.285052 0.560831 2.721043 0.034519 3.832200 0.041440 0.116405 0.056267 0.497593 0.291009 0.623366 0.256174 1.144639 0.524647 1.038682 12.931524 \
128 0.225554 0.000047 0.010929 0.009289 12.169045 0.083636 5.964323 0.681575 0.506470 0.000180 2.007768 0.181794 0.139046 0.322807 0.000060 0.009512 0.000105 1.035015 0.000288 0.021048 0.001066 1.293228 0.000453 1.177718 0.083261 0.772349 0.018146 0.530881 0.025006 0.359183 0.018727 0.269862 0.306375 4.943439 0.275865 2.397415 0.563566 4.971507 0.586685 1.293860 0.389004 \
129 0.000093 0.166706 0.000432 0.005790 0.094762 10.892976 1.049877 9.818281 0.000116 0.346890 0.248099 1.372357 0.167138 0.000330 0.300513 0.002718 0.017300 0.001247 1.337629 0.005195 0.104681 0.005060 1.282522 0.232701 1.418383 0.387941 1.320875 0.354545 1.360141 0.031140 0.238533 0.021539 0.304581 0.622868 4.699375 0.441084 2.871848 0.643789 4.127466 0.334224 0.928876 28.579806 \
130 0.140516 0.012600 0.423774 0.001718 0.047890 0.002044 1.094736 0.000115 0.424321 0.060909 0.144388 0.030883 0.145245 0.004747 0.000060 0.409432 0.000000 0.011511 0.000384 0.493508 0.000000 0.411115 0.025544 1.242140 0.000289 17.450524 1.113671 31.949764 2.418859 0.116039 0.012331 0.780008 0.305714 0.432918 0.021698 1.316696 0.087905 0.936840 0.273855 5.815294 1.197614 1.644621 0.403913 \
131 0.083310 0.056771 0.000247 0.506841 0.063994 0.028715 0.009261 0.514392 0.200499 0.269510 0.122186 0.070533 0.496706 0.009560 0.000952 0.000951 0.040320 0.060432 0.033244 0.009225 0.273876 0.140943 0.169022 0.003786 1.148387 6.798629 4.087042 15.287419 18.531553 0.093542 0.062369 0.317466 0.905810 0.309656 0.140701 0.511968 0.765495 0.454347 0.600415 1.868194 7.316623 1.477696 1.286990 43.916187 \
132 0.863970 0.065905 0.748196 0.529619 0.563995 0.000186 0.002219 0.000115 0.571505 0.000359 1.598824 0.004316 0.038763 1.897150 0.072866 0.555104 0.011580 0.708395 0.000192 0.009225 0.000246 0.338207 0.000091 0.150804 0.009600 0.336828 0.000000 0.003743 0.000000 6.546163 0.575921 2.577578 0.430124 2.898791 0.003020 0.056250 0.004868 0.318595 0.000295 0.201480 0.072138 0.501456 0.000219 0.032116 0.051709 \
133 0.026338 0.946136 0.005990 0.140867 0.000085 0.396897 0.000096 0.001261 0.000116 0.582801 0.001087 1.273908 0.092044 0.105361 2.720153 0.043756 1.085170 0.000192 0.760090 0.000537 0.031642 0.000375 0.491034 0.006757 0.069320 0.000589 0.648523 0.001871 0.005026 0.458622 7.487816 0.154207 0.350473 0.001027 2.862216 0.002515 0.010298 0.000000 0.215492 0.001930 0.056145 0.000097 0.381287 0.000205 0.002062 10.956917 \
134 0.565566 0.047403 2.299543 0.762425 0.001356 0.000279 0.813720 0.000115 0.236179 0.060430 0.754233 0.086986 0.058616 0.509595 0.031730 2.047159 0.020529 0.001151 0.000192 0.732111 0.000082 0.019586 0.002808 0.606945 0.000145 0.000353 0.000000 0.255147 0.000000 2.581654 0.313779 11.271062 0.569076 0.009008 0.001957 3.879355 0.002528 0.008844 0.002362 0.708204 0.000524 0.006305 0.000657 0.560642 0.002062 25.313949 5.637509 \
135 0.068927 0.371072 0.024699 1.108802 0.000254 0.001394 0.000772 0.451899 0.009630 0.116309 0.126844 0.530278 0.277072 0.091844 0.477439 0.173122 2.262490 0.000384 0.001537 0.000717 0.602428 0.032237 0.044112 0.000466 0.475496 0.001413 0.003287 0.000832 0.701399 0.953353 3.487342 0.911193 1.399673 0.003635 0.089643 0.011433 3.092500 0.000628 0.013582 0.000193 0.311885 0.001358 0.004160 0.000103 0.314379 8.832621 18.744445 13.945647 \
136 0.483563 0.000234 0.008892 0.035630 3.994417 0.771771 1.825878 0.250545 0.370309 0.000419 1.899865 0.011733 0.035860 0.879741 0.000417 0.004077 0.000772 1.272426 0.073021 0.346978 0.013526 0.420580 0.000091 0.272950 0.040907 0.324227 0.000000 0.001871 0.000000 0.536332 0.000253 0.001795 0.451134 2.359362 0.120121 0.324162 0.108501 0.643160 0.001329 0.216244 0.268662 2.323091 0.005328 0.013852 0.045374 2.600428 0.131929 0.662578 0.152215 \
137 0.000093 0.512203 0.000000 0.001654 0.111968 3.452570 0.084218 1.978448 0.000058 0.380788 0.000466 1.514097 0.044178 0.000132 1.191514 0.000136 0.103766 0.022062 1.394892 0.013077 0.971148 0.000562 0.696016 0.018814 0.634927 0.000236 1.610248 0.000416 0.020471 0.000040 0.536362 0.000000 0.169264 0.038559 2.159328 0.019665 0.498504 0.000045 0.565968 0.005500 0.373948 0.000485 2.214735 0.000000 0.001768 0.046583 2.620833 0.028569 0.682579 9.612709 \
138 0.109975 0.016535 1.041312 0.019406 1.931350 0.558500 4.380679 0.505677 0.176829 0.034737 0.806554 0.297371 0.031466 0.002374 0.000357 0.741135 0.000351 0.335924 0.069754 1.236457 0.087384 0.039265 0.004982 1.274118 0.002219 0.000589 0.000068 0.286547 0.000123 0.002262 0.000589 0.983926 0.517896 0.381796 0.077341 2.735831 0.318574 0.084620 0.041435 0.923644 0.004382 0.126382 0.063353 0.461729 0.004420 0.718719 0.092405 3.415722 0.415718 24.400553 6.746560 \
139 0.005884 0.074851 0.000000 0.220908 0.103323 1.262618 0.150589 4.658653 0.027035 0.106187 0.028567 0.586111 0.446015 0.000066 0.000893 0.000000 1.524024 0.014101 0.417565 0.017824 1.950083 0.080124 0.190037 0.001165 1.544626 0.001531 0.083744 0.000624 3.409178 0.000081 0.004629 0.000078 0.837302 0.023862 0.728891 0.049848 2.866325 0.003771 0.068501 0.000482 0.759132 0.006402 0.200205 0.000000 0.187832 0.054049 0.968351 0.081861 2.211488 5.140068 19.373137 11.561124 \
140 0.064397 0.000000 0.042112 0.038557 1.120532 0.003717 0.348448 0.117533 0.223763 0.015452 0.099985 0.000135 0.028249 0.129492 0.000000 0.012366 0.000491 0.661776 0.000769 0.147873 0.031560 0.746792 0.046739 0.706782 0.130873 0.162525 0.000000 0.007070 0.000368 0.066966 0.000042 0.001171 0.059065 0.928969 0.000559 0.092988 0.042595 3.529593 0.371685 0.604859 0.188097 1.702817 0.012481 0.030474 0.015763 0.153418 0.007112 0.078381 0.011491 0.396521 0.015140 0.189090 0.043198 \
141 0.000000 0.055366 0.000062 0.006808 0.000254 1.023142 0.007428 0.670108 0.010037 0.184704 0.000000 0.071612 0.066384 0.000066 0.135255 0.001359 0.015686 0.000096 0.976175 0.003672 0.644235 0.100928 0.975727 0.121389 0.928319 0.000236 0.915505 0.009981 0.150527 0.000000 0.032447 0.000000 0.011379 0.000158 1.013424 0.003354 0.095207 0.167041 2.729647 0.053168 0.426684 0.000388 2.005334 0.000718 0.008986 0.004101 0.119062 0.006776 0.041280 0.018617 0.802516 0.027912 0.702594 14.214694 \
142 0.084945 0.006464 0.287373 0.005472 0.330481 0.085680 1.265487 0.002179 0.257122 0.043721 0.028878 0.003641 0.009966 0.039560 0.002679 0.313495 0.000140 0.184749 0.105112 0.890822 0.005410 0.452442 0.106069 3.081614 0.536567 0.034978 0.025678 0.440217 0.000858 0.038612 0.009174 0.361403 0.033994 0.251423 0.109664 1.164866 0.003464 0.975582 0.193544 2.258321 0.308851 0.832592 0.308372 0.668173 0.004420 0.276499 0.042565 0.469281 0.055025 0.502355 0.140546 0.905488 0.227527 2.738552 0.892903 \
143 0.010974 0.034428 0.000000 0.159955 0.042380 0.283432 0.001061 1.029128 0.042815 0.136432 0.014439 0.013216 0.137634 0.004220 0.010061 0.000136 0.176300 0.034437 0.294294 0.001791 0.990330 0.159217 0.566034 0.343314 3.036767 0.007891 0.528692 0.001040 2.171984 0.003312 0.031984 0.000078 0.262465 0.033581 0.360196 0.000838 1.447392 0.149578 0.372719 0.159248 1.563846 0.129098 0.822643 0.000410 1.195790 0.049842 0.245019 0.053017 0.362328 0.106257 0.938586 0.157605 1.251589 1.091224 3.195698 12.984714 \
144 0.164659 0.000141 0.000741 0.003881 0.976185 0.001951 0.011673 0.007109 0.130940 0.000120 0.420899 0.045044 0.039313 0.169777 0.000060 0.000272 0.000175 0.418802 0.000288 0.002508 0.001312 0.388156 0.000091 0.042812 0.003377 0.241197 0.004656 0.042005 0.011768 0.069995 0.000000 0.000156 0.027479 0.380374 0.000112 0.000534 0.000374 1.322234 0.005905 0.048730 0.021649 2.382451 0.326035 0.037657 0.047437 0.164143 0.016776 0.072521 0.024883 1.572808 0.086923 0.585071 0.083552 0.629243 0.035120 0.089148 0.030223 \
145 0.000000 0.172889 0.000000 0.000191 0.000085 0.880032 0.000289 0.356038 0.000058 0.127388 0.007608 0.309374 0.105305 0.000000 0.240505 0.000000 0.047268 0.000096 0.636916 0.000090 0.395771 0.000843 0.566759 0.016193 0.336277 0.021435 0.676049 0.008942 0.703728 0.000283 0.055425 0.000000 0.018603 0.000000 0.518903 0.000000 0.006459 0.001122 1.110726 0.002863 0.176224 0.054025 2.392606 0.000821 0.012227 0.002050 0.201477 0.001557 0.051048 0.022214 1.797671 0.027973 1.398079 0.037461 1.228004 0.030585 0.239725 13.935950";
146 string model_ECMunrest = model_ECMunrest1 +
147 "0.113991 0.018315 0.201112 0.001082 0.012121 0.001951 1.720919 0.001720 0.082323 0.029826 0.197641 0.061497 0.073682 0.000330 0.000060 0.165784 0.000070 0.003549 0.000384 0.556204 0.000164 0.097554 0.004982 0.551493 0.000289 0.015310 0.000753 0.247245 0.010419 0.000283 0.000084 0.194319 0.037724 0.002449 0.000112 0.466770 0.000187 0.909861 0.280400 0.713961 0.001760 1.179053 0.298738 0.938439 0.165587 0.080337 0.009773 0.324696 0.016839 0.658541 0.036022 1.693998 0.046588 0.375097 0.067431 0.639125 0.053748 21.171295 5.214689 \
148 0.018773 0.032039 0.000000 0.175861 0.002797 0.002974 0.003376 2.163175 0.007948 0.014314 0.105884 0.183952 0.381671 0.000066 0.000119 0.000000 0.185038 0.001918 0.001441 0.001254 0.703092 0.084060 0.053714 0.003029 0.634203 0.043222 0.097165 0.143481 0.590833 0.000081 0.000295 0.000078 0.410199 0.000553 0.000447 0.000610 0.716441 0.194964 0.293884 0.001158 0.744000 0.684968 1.149846 0.069567 1.558784 0.032177 0.064227 0.074536 0.276276 0.238907 0.496552 0.672077 1.526141 0.235747 0.403521 0.136937 0.968146 13.981617 18.675227 25.640860 \
149 \
150 0.021414 0.021349 0.016195 0.015717 0.011798 0.010761 0.010366 0.008721 0.017237 0.016697 0.006441 0.007415 0.012744 0.015167 0.016798 0.007359 0.028497 0.010425 0.010408 0.011165 0.012199 0.010671 0.011040 0.017168 0.020730 0.008491 0.014604 0.004809 0.008158 0.024759 0.023762 0.012814 0.021180 0.012656 0.017882 0.013120 0.010682 0.022276 0.020321 0.031090 0.026699 0.010310 0.013701 0.009746 0.006788 0.019020 0.018419 0.010921 0.022626 0.018907 0.026817 0.016516 0.018288 0.028590 0.025285 0.034527 0.030606 0.016883 0.023659 0.016386 0.010223 \
151 \
152 TTT TTC TTA TTG TCT TCC TCA TCG TAT TAC TGT TGC TGG CTT CTC CTA CTG CCT CCC CCA \
153 CCG CAT CAC CAA CAG CGT CGC CGA CGG ATT ATC ATA ATG ACT ACC ACA ACG AAT AAC AAA \
154 AAG AGT AGC AGA AGG GTT GTC GTA GTG GCT GCC GCA GCG GAT GAC GAA GAG GGT GGC GGA \
155 GGG";
156
157 /* empirical codon model of Schneider, Cannarozzi, Gonnet 2005
158 */
159
160 string model_ECM_Schneider05 =
161 "15594\
162 787 609\
163 717 391 22864\
164 476 0 89 133\
165 0 444 39 51 15656\
166 77 11 2230 662 15154 11009\
167 99 0 0 2501 21330 19335 36566\
168 2355 360 127 41 400 18 28 0\
169 300 1911 76 49 0 260 23 70 25720\
170 378 30 81 96 1004 0 238 90 1025 0\
171 60 275 46 53 37 766 83 157 0 623 29318\
172 102 74 24 357 16 23 83 277 110 80 186 164\
173 1304 0 4065 4641 278 0 0 0 138 18 147 21 9\
174 0 1293 3595 2691 0 214 0 0 1 128 9 152 9 15530\
175 112 106 26307 6745 0 0 154 340 68 0 59 0 36 10903 8429\
176 186 136 4180 13314 20 34 0 0 39 14 17 9 51 7205 7138 16270\
177 97 0 7 72 2744 0 0 0 43 3 101 12 0 937 34 0 0\
178 0 82 18 3 65 2535 0 0 0 47 0 49 0 0 706 0 0 16327\
179 0 12 99 87 0 0 2829 131 7 0 11 0 14 0 0 1474 137 15604 10490\
180 26 3 11 109 0 0 244 5851 4 22 44 31 68 0 0 0 935 15588 21512 27394\
181 172 22 80 27 315 22 76 0 2184 0 475 6 63 676 55 0 43 756 38 56 173\
182 36 166 5 20 43 211 53 112 137 1624 21 276 20 0 390 96 40 0 617 97 246 30148\
183 35 18 93 63 143 90 403 0 122 56 70 60 38 168 38 1166 51 126 99 1405 296 1995 1176\
184 11 13 50 95 55 80 109 288 42 58 43 35 97 77 92 46 290 88 161 235 1234 1086 1057 15317\
185 70 7 49 98 326 11 1 0 362 0 1302 0 14 685 0 0 0 712 0 0 0 4390 0 646 232\
186 0 56 40 0 0 177 0 105 4 139 0 908 40 1 486 0 34 0 404 8 121 0 2689 217 355 42123\
187 8 0 0 25 0 0 163 53 68 0 51 53 130 0 45 966 0 0 23 393 0 133 261 3956 0 40836 19396\
188 10 8 2 45 0 50 31 385 0 17 50 109 576 21 0 7 340 0 0 57 1276 191 317 0 2408 25507 19595 39593\
189 437 29 201 399 72 30 0 18 64 1 90 1 10 1848 131 185 283 68 10 24 0 84 4 0 20 27 0 64 0\
190 23 356 111 27 15 35 6 15 7 42 10 62 7 133 1745 0 351 0 55 0 29 0 36 24 7 17 20 0 15 16920\
191 157 33 3243 833 0 54 209 88 29 1 14 16 14 134 176 3786 421 49 13 97 79 51 4 93 20 43 3 5 0 15113 13126\
192 126 89 536 1607 75 39 153 229 32 21 29 22 45 346 369 726 1104 34 27 55 92 61 20 96 116 29 23 0 50 781 495 2443\
193 82 23 67 79 2679 170 489 362 61 17 194 41 0 233 0 0 25 990 52 0 15 212 6 166 30 318 42 0 0 1552 0 0 262\
194 7 74 3 39 279 2315 113 0 21 33 10 153 9 0 186 0 44 76 817 0 0 15 171 37 61 14 118 71 33 0 1028 0 176 16316\
195 23 11 295 92 567 81 3129 578 12 0 75 22 19 4 0 288 46 64 0 1015 88 93 27 341 96 39 0 0 26 0 0 2434 714 16829 10306\
196 51 5 49 306 0 682 232 6944 52 9 134 0 44 55 28 0 181 0 170 0 1953 0 35 0 307 0 116 322 172 166 134 0 1616 20814 19942 34110\
197 41 0 48 36 289 32 139 78 316 0 191 0 11 110 0 22 0 115 34 55 46 1623 55 311 168 429 0 45 27 202 0 50 68 1425 16 237 182\
198 8 34 4 22 87 277 79 40 0 196 4 116 7 0 64 19 7 14 98 17 47 9 1196 148 142 0 264 0 58 0 145 21 42 122 984 165 67 17900\
199 4 6 59 28 54 24 133 53 16 8 20 15 9 30 17 100 0 26 18 91 39 169 94 1364 137 334 337 908 410 16 0 292 80 114 54 650 0 580 305\
200 6 0 4 39 26 24 52 70 11 11 13 16 24 7 3 0 36 14 14 28 105 84 110 33 727 353 569 299 1249 20 21 22 189 27 85 94 795 298 269 12057\
201 49 0 0 47 938 51 528 366 157 0 981 105 27 71 0 26 4 178 0 76 37 605 0 198 84 1145 0 0 44 345 0 91 78 2941 0 402 387 4830 0 245 146\
202 10 34 33 0 312 983 214 697 10 79 169 792 19 23 68 0 21 15 213 25 88 41 454 133 123 0 740 72 22 0 228 27 67 421 2248 343 618 178 3430 158 131 21641\
203 10 2 87 13 54 26 79 12 0 20 52 7 23 32 0 9 13 5 17 0 45 118 55 623 32 6203 2803 22415 6405 0 0 415 93 11 31 743 0 178 142 4144 390 638 339\
204 12 0 51 36 33 12 4 72 31 0 36 20 264 0 20 0 22 17 0 39 5 62 95 60 271 2885 5913 6064 18118 29 15 106 258 40 59 49 927 145 139 572 3644 544 444 25118\
205 446 0 175 236 218 61 96 0 80 0 108 60 20 1661 28 0 11 258 0 0 65 68 1 82 28 96 15 0 3 6623 53 389 259 1227 72 65 141 100 20 37 3 233 69 31 0\
206 6 460 73 35 15 199 49 0 7 69 28 105 4 18 1603 120 42 0 199 36 69 15 75 29 25 18 52 0 36 79 6084 158 190 148 1080 1 0 21 86 13 11 16 170 5 38 17102\
207 83 13 3202 251 70 2 396 0 18 2 52 24 6 131 5 2164 7 0 63 394 0 23 36 276 31 0 0 86 0 125 0 11512 626 0 20 1694 36 55 26 153 16 31 87 249 0 15228 10722\
208 60 42 0 1017 46 19 85 271 9 17 54 16 27 0 24 16 827 27 21 58 154 23 18 10 62 12 20 0 58 743 928 1567 1219 155 66 206 995 12 12 10 54 28 10 0 106 10663 9288 21679\
209 89 10 43 68 2454 180 301 0 54 11 203 32 0 183 8 14 31 1111 63 28 0 111 5 60 52 121 0 0 62 233 32 42 108 3490 0 0 0 252 38 43 31 628 101 47 0 2567 0 14 37\
210 1 78 42 24 250 2338 201 324 11 28 22 197 15 30 160 0 31 79 965 7 145 0 119 22 88 29 133 26 0 0 217 0 115 0 2888 0 36 25 224 13 45 0 608 24 51 0 2241 0 80 12825\
211 39 0 293 82 303 127 3185 250 19 1 87 27 17 29 0 300 46 91 19 1186 101 70 20 307 61 38 0 72 23 73 0 484 282 0 0 3718 0 149 28 216 48 268 106 108 88 0 0 3240 591 14816 8169\
212 19 59 0 480 450 354 0 9766 0 68 96 134 32 16 40 0 261 0 206 0 4882 109 45 182 246 0 94 0 379 0 77 315 392 0 193 0 7718 81 161 62 171 178 258 49 171 0 0 0 2849 13866 14245 22755\
213 14 0 0 16 138 0 52 20 122 0 41 0 3 21 0 5 0 44 7 28 37 378 0 132 54 66 8 0 13 36 0 19 8 149 2 83 0 1636 4 56 44 488 30 27 19 187 0 24 5 384 0 59 20\
214 0 13 16 0 1 115 35 76 0 76 9 24 1 10 20 0 3 8 72 6 40 22 311 76 72 1 57 0 11 0 20 0 9 54 124 0 77 33 1249 58 35 38 395 22 11 0 141 19 0 0 322 39 182 13482\
215 8 2 32 46 42 14 118 53 16 9 13 0 5 13 6 52 0 29 17 71 35 123 28 1246 112 69 0 24 0 0 5 95 22 38 36 196 0 187 74 626 0 114 82 213 0 62 24 367 16 64 41 636 0 1758 1114\
216 3 4 20 18 36 42 43 84 16 6 3 11 13 7 12 0 27 19 40 51 81 69 71 66 739 0 44 0 104 21 1 17 63 52 35 76 199 106 80 0 429 101 86 0 162 4 26 38 161 70 91 118 916 1216 1466 11181\
217 32 0 44 0 293 18 60 85 62 0 515 0 0 50 0 9 9 161 16 29 0 266 0 0 43 688 0 117 0 46 0 0 23 349 91 46 0 829 98 44 32 3432 6 0 0 557 0 0 21 1264 0 171 0 1058 0 28 30\
218 9 23 0 15 52 224 28 240 0 38 38 372 7 14 38 0 6 14 114 31 129 12 165 15 40 1 452 0 78 10 42 18 16 50 250 6 62 91 589 13 28 0 2627 25 0 15 383 0 0 72 1050 0 357 0 795 0 50 17792\
219 15 1 23 0 80 59 198 0 9 0 34 0 4 13 12 70 0 16 11 78 26 26 16 225 8 0 0 176 0 16 0 85 31 3 12 214 51 121 31 137 11 151 134 1071 0 0 0 629 0 133 44 1071 0 50 86 798 0 13656 9267\
220 2 11 9 122 10 49 162 532 14 3 42 39 214 0 1 0 39 9 55 36 251 29 31 53 148 0 0 0 515 0 22 14 116 77 0 84 405 67 171 13 144 398 328 0 1567 0 0 86 380 35 151 360 2466 88 32 0 749 11092 11188 21307\
221 0.019065 0.019572 0.008521 0.014125 0.017115 0.015966 0.013340 0.004325 0.013196 0.015988 0.010982 0.011900 0.012462 0.015017 0.017340 0.008031 0.037328 0.017599 0.014989 0.017755 0.005879 0.011592 0.014372 0.013669 0.033679 0.005419 0.008712 0.006185 0.008488 0.017469 0.019956 0.009364 0.021871 0.014391 0.015964 0.016870 0.005890 0.018236 0.020614 0.028227 0.031883 0.013622 0.019058 0.013502 0.011845 0.013592 0.013760 0.008372 0.026506 0.019952 0.021260 0.017896 0.005964 0.025167 0.024586 0.031093 0.038868 0.011549 0.017608 0.018437 0.013269\
222 TTT TTC TTA TTG TCT TCC TCA TCG TAT TAC TGT TGC TGG CTT CTC CTA CTG CCT CCC CCA \
223 CCG CAT CAC CAA CAG CGT CGC CGA CGG ATT ATC ATA ATG ACT ACC ACA ACG AAT AAC AAA \
224 AAG AGT AGC AGA AGG GTT GTC GTA GTG GCT GCC GCA GCG GAT GAC GAA GAG GGT GGC GGA \
225 GGG";
226
227
ModelCodon(const char * model_name,string model_params,StateFreqType freq,string freq_params,PhyloTree * tree)228 ModelCodon::ModelCodon(const char *model_name, string model_params, StateFreqType freq, string freq_params,
229 PhyloTree *tree) : ModelMarkov(tree)
230 {
231 half_matrix = false;
232 omega = kappa = kappa2 = 1.0;
233 fix_omega = fix_kappa = false;
234 fix_kappa2 = true;
235 codon_freq_style = CF_TARGET_CODON;
236 codon_kappa_style = CK_ONE_KAPPA;
237 ntfreq = new double[12];
238 empirical_rates = NULL;
239 int nrates = getNumRateEntries();
240 delete [] rates;
241 rates = new double[nrates];
242 empirical_rates = new double [nrates];
243
244 rate_attr = NULL;
245 computeRateAttributes();
246
247 init(model_name, model_params, freq, freq_params);
248 }
249
~ModelCodon()250 ModelCodon::~ModelCodon() {
251 if (rate_attr) {
252 delete [] rate_attr;
253 rate_attr = NULL;
254 }
255 if (empirical_rates) {
256 delete [] empirical_rates;
257 empirical_rates = NULL;
258 }
259 if (ntfreq) {
260 delete [] ntfreq;
261 ntfreq = NULL;
262 }
263 }
264
startCheckpoint()265 void ModelCodon::startCheckpoint() {
266 checkpoint->startStruct("ModelCodon");
267 }
268
saveCheckpoint()269 void ModelCodon::saveCheckpoint() {
270 startCheckpoint();
271 // CKP_ARRAY_SAVE(12, ntfreq);
272 CKP_SAVE(omega);
273 // CKP_SAVE(fix_omega);
274 // int codon_kappa_style = this->codon_kappa_style;
275 // CKP_SAVE(codon_kappa_style);
276 CKP_SAVE(kappa);
277 // CKP_SAVE(fix_kappa);
278 CKP_SAVE(kappa2);
279 // CKP_SAVE(fix_kappa2);
280 // int codon_freq_style = this->codon_freq_style;
281 // CKP_SAVE(codon_freq_style);
282 endCheckpoint();
283
284 ModelMarkov::saveCheckpoint();
285 }
286
restoreCheckpoint()287 void ModelCodon::restoreCheckpoint() {
288 ModelMarkov::restoreCheckpoint();
289 startCheckpoint();
290 // CKP_ARRAY_RESTORE(12, ntfreq);
291 CKP_RESTORE(omega);
292 // CKP_RESTORE(fix_omega);
293 // int codon_kappa_style;
294 // CKP_RESTORE(codon_kappa_style);
295 // this->codon_kappa_style = (CodonKappaStyle)codon_kappa_style;
296 CKP_RESTORE(kappa);
297 // CKP_RESTORE(fix_kappa);
298 CKP_RESTORE(kappa2);
299 // CKP_RESTORE(fix_kappa2);
300 // int codon_freq_style;
301 // CKP_RESTORE(codon_freq_style);
302 // this->codon_freq_style = (CodonFreqStyle)codon_freq_style;
303 endCheckpoint();
304
305 decomposeRateMatrix();
306 if (phylo_tree)
307 phylo_tree->clearAllPartialLH();
308
309 }
310
initCodon(const char * model_name,StateFreqType freq,bool reset_params)311 StateFreqType ModelCodon::initCodon(const char *model_name, StateFreqType freq, bool reset_params) {
312 string name_upper = model_name;
313 for (string::iterator it = name_upper.begin(); it != name_upper.end(); it++)
314 (*it) = toupper(*it);
315
316 if (name_upper == "MG") {
317 return initMG94(true, freq, CK_ONE_KAPPA);
318 } else if (name_upper == "MGK") {
319 return initMG94(false, freq, CK_ONE_KAPPA);
320 } else if (name_upper == "MG1KTS" || name_upper == "MGKAP2") {
321 return initMG94(false, freq, CK_ONE_KAPPA_TS);
322 } else if (name_upper == "MG1KTV" || name_upper == "MGKAP3") {
323 return initMG94(false, freq, CK_ONE_KAPPA_TV);
324 } else if (name_upper == "MG2K" || name_upper == "MGKAP4") {
325 return initMG94(false, freq, CK_TWO_KAPPA);
326 } else if (name_upper == "GY") {
327 return initGY94(false, CK_ONE_KAPPA);
328 } else if (name_upper == "GY0K" || name_upper == "GYKAP1") {
329 return initGY94(true, CK_ONE_KAPPA);
330 } else if (name_upper == "GY1KTS" || name_upper == "GYKAP2") {
331 return initGY94(false, CK_ONE_KAPPA_TS);
332 } else if (name_upper == "GY1KTV" || name_upper == "GYKAP3") {
333 return initGY94(false, CK_ONE_KAPPA_TV);
334 } else if (name_upper == "GY2K" || name_upper == "GYKAP4") {
335 return initGY94(false, CK_TWO_KAPPA);
336 } else if (name_upper == "ECM" || name_upper == "KOSI07" || name_upper == "ECMK07") {
337 if (!phylo_tree->aln->isStandardGeneticCode())
338 outError("For ECMK07 a standard genetic code must be used");
339 readCodonModel(model_ECMunrest, reset_params);
340 return FREQ_USER_DEFINED;
341 } else if (name_upper == "ECMREST") {
342 if (!phylo_tree->aln->isStandardGeneticCode())
343 outError("For ECMREST a standard genetic code must be used");
344 readCodonModel(model_ECMrest, reset_params);
345 return FREQ_USER_DEFINED;
346 } else if (name_upper == "SCHN05" || name_upper == "ECMS05") {
347 if (!phylo_tree->aln->isStandardGeneticCode())
348 outError("For ECMS05 a standard genetic code must be used");
349 readCodonModel(model_ECM_Schneider05, reset_params);
350 return FREQ_USER_DEFINED;
351 } else {
352 //cout << "User-specified model "<< model_name << endl;
353 //readParameters(model_name);
354 //name += " (user-defined)";
355 readCodonModelFile(model_name, reset_params);
356 return FREQ_USER_DEFINED;
357 }
358
359 return FREQ_UNKNOWN;
360 }
361
init(const char * model_name,string model_params,StateFreqType freq,string freq_params)362 void ModelCodon::init(const char *model_name, string model_params, StateFreqType freq, string freq_params)
363 {
364 int i, j;
365 for (i = 0; i < 12; i++)
366 ntfreq[i] = 0.25;
367 // initialize empirical_rates
368 for (i = 0; i < num_states; i++) {
369 double *this_emp_rate = &empirical_rates[i*num_states];
370 int *this_rate_attr = &rate_attr[i*num_states];
371 if (phylo_tree->aln->isStopCodon(i)) {
372 memset(this_emp_rate, 0, num_states*sizeof(double));
373 continue;
374 }
375 for (j = 0; j < num_states; j++) {
376 int attr = this_rate_attr[j];
377 if (attr & (CA_STOP_CODON+CA_MULTI_NT)) { // stop codon or multiple nt substitutions
378 this_emp_rate[j] = 0.0;
379 } else {
380 this_emp_rate[j] = 1.0;
381 }
382 }
383 }
384
385 ignore_state_freq = false;
386
387 StateFreqType def_freq = FREQ_UNKNOWN;
388 name = full_name = model_name;
389 size_t pos;
390 if ((pos=name.find('_')) == string::npos) {
391 def_freq = initCodon(model_name, freq, true);
392 } else {
393 def_freq = initCodon(name.substr(0, pos).c_str(), freq, false);
394 if (def_freq != FREQ_USER_DEFINED)
395 outError("Invalid model " + name + ": first component must be an empirical model"); // first model must be empirical
396 def_freq = initCodon(name.substr(pos+1).c_str(), freq, false);
397 if (def_freq == FREQ_USER_DEFINED) // second model must be parametric
398 outError("Invalid model " + name + ": second component must be a mechanistic model");
399 // adjust the constraint
400 if (codon_freq_style==CF_TARGET_CODON)
401 def_freq = FREQ_USER_DEFINED;
402 }
403
404 num_params = (!fix_omega) + (!fix_kappa) + (!fix_kappa2);
405
406 if (freq_params != "") {
407 readStateFreq(freq_params);
408 }
409 if (model_params != "") {
410 readRates(model_params);
411 }
412
413 // if (freq == FREQ_UNKNOWN || def_freq == FREQ_EQUAL) freq = def_freq;
414 if (freq == FREQ_UNKNOWN) freq = def_freq;
415 if (freq == FREQ_CODON_1x4 || freq == FREQ_CODON_3x4 || freq == FREQ_CODON_3x4C) {
416 //ntfreq = new double[12];
417 phylo_tree->aln->computeCodonFreq(freq, state_freq, ntfreq);
418 }
419 ModelMarkov::init(freq);
420 }
421
initMG94(bool fix_kappa,StateFreqType freq,CodonKappaStyle kappa_style)422 StateFreqType ModelCodon::initMG94(bool fix_kappa, StateFreqType freq, CodonKappaStyle kappa_style) {
423 /* Muse-Gaut 1994 model with 1 parameters: omega */
424
425 fix_omega = false;
426 this->fix_kappa = fix_kappa;
427 if (fix_kappa)
428 kappa = 1.0;
429 fix_kappa2 = true;
430 codon_freq_style = CF_TARGET_NT;
431 this->codon_kappa_style = kappa_style;
432 if (kappa_style == CK_TWO_KAPPA)
433 fix_kappa2 = false;
434
435 if (freq == FREQ_UNKNOWN || freq == FREQ_USER_DEFINED)
436 freq = FREQ_CODON_3x4;
437
438 switch (freq) {
439 case FREQ_CODON_1x4:
440 case FREQ_CODON_3x4:
441 case FREQ_CODON_3x4C:
442 phylo_tree->aln->computeCodonFreq(freq, state_freq, ntfreq);
443 break;
444 case FREQ_EMPIRICAL:
445 case FREQ_ESTIMATE:
446 case FREQ_USER_DEFINED:
447 outError("Invalid state frequency type for MG model, please use +F1X4 or +F3X4 or +F3X4C");
448 break;
449 default:
450 break;
451 }
452
453 // ignote state_freq because ntfreq is already used
454 ignore_state_freq = true;
455 combineRateNTFreq();
456
457 return FREQ_CODON_3x4;
458 }
459
initGY94(bool fix_kappa,CodonKappaStyle kappa_style)460 StateFreqType ModelCodon::initGY94(bool fix_kappa, CodonKappaStyle kappa_style) {
461 fix_omega = false;
462 this->fix_kappa = fix_kappa;
463 if (fix_kappa)
464 kappa = 1.0;
465 fix_kappa2 = true;
466 codon_freq_style = CF_TARGET_CODON;
467 this->codon_kappa_style = kappa_style;
468 if (kappa_style == CK_TWO_KAPPA)
469 fix_kappa2 = false;
470
471 return FREQ_EMPIRICAL;
472 }
473
474
computeRateAttributes()475 void ModelCodon::computeRateAttributes() {
476 char symbols_protein[] = "ARNDCQEGHILKMFPSTWYVX"; // X for unknown AA
477 int i, j, ts, tv;
478 int nrates = getNumRateEntries();
479 if (!rate_attr) {
480 rate_attr = new int[nrates];
481 memset(rate_attr, 0, sizeof(int)*nrates);
482 }
483 char aa_cost_change[400];
484 memset(aa_cost_change, 10, sizeof(char)*400);
485 for (i = 0; i < num_states; i++) {
486 int *rate_attr_row = &rate_attr[i*num_states];
487 if (phylo_tree->aln->isStopCodon(i)) {
488 for (j = 0; j < num_states; j++)
489 rate_attr_row[j] = CA_STOP_CODON;
490 continue;
491 }
492 for (j = 0; j < num_states; j++) {
493 if (j == i || phylo_tree->aln->isStopCodon(j)) {
494 rate_attr_row[j] = CA_STOP_CODON;
495 continue;
496 }
497 int nuc1, nuc2;
498 int attr = 0;
499 ts = tv = 0;
500 int codoni = phylo_tree->aln->codon_table[i];
501 int codonj = phylo_tree->aln->codon_table[j];
502 if (phylo_tree->aln->genetic_code[codoni] == phylo_tree->aln->genetic_code[codonj])
503 attr |= CA_SYNONYMOUS;
504 else
505 attr |= CA_NONSYNONYMOUS;
506
507
508 int nt_changes = ((codoni/16) != (codonj/16)) + (((codoni%16)/4) != ((codonj%16)/4)) + ((codoni%4) != (codonj%4));
509 int aa1 = strchr(symbols_protein, phylo_tree->aln->genetic_code[codoni]) - symbols_protein;
510 int aa2 = strchr(symbols_protein, phylo_tree->aln->genetic_code[codonj]) - symbols_protein;
511 ASSERT(aa1 >= 0 && aa1 < 20 && aa2 >= 0 && aa2 < 20);
512 if (nt_changes < aa_cost_change[aa1*20+aa2]) {
513 aa_cost_change[aa1*20+aa2] = aa_cost_change[aa2*20+aa1] = nt_changes;
514 }
515
516 if ((nuc1=codoni/16) != (nuc2=codonj/16)) {
517 if (abs(nuc1-nuc2)==2) { // transition
518 attr |= CA_TRANSITION_1NT;
519 ts++;
520 } else { // transversion
521 attr |= CA_TRANSVERSION_1NT;
522 tv++;
523 }
524 }
525 if ((nuc1=(codoni%16)/4) != (nuc2=(codonj%16)/4)) {
526 if (abs(nuc1-nuc2)==2) { // transition
527 attr |= CA_TRANSITION_2NT;
528 ts++;
529 } else { // transversion
530 attr |= CA_TRANSVERSION_2NT;
531 tv++;
532 }
533 }
534 if ((nuc1=codoni%4) != (nuc2=codonj%4)) {
535 if (abs(nuc1-nuc2)==2) { // transition
536 attr |= CA_TRANSITION_3NT;
537 ts++;
538 } else { // transversion
539 attr |= CA_TRANSVERSION_3NT;
540 tv++;
541 }
542 }
543 if (ts+tv>1)
544 attr |= CA_MULTI_NT;
545 else if (ts==1)
546 attr |= CA_TRANSITION;
547 else if (tv==1)
548 attr |= CA_TRANSVERSION;
549
550 rate_attr_row[j] = attr;
551 }
552 }
553
554 if (verbose_mode >= VB_MAX) {
555
556 // make cost matrix fulfill triangular inequality
557 for (int k = 0; k < 20; k++)
558 for (i = 0; i < 20; i++)
559 for (j = 0; j < 20; j++)
560 if (aa_cost_change[i*20+j] > aa_cost_change[i*20+k] + aa_cost_change[k*20+j])
561 aa_cost_change[i*20+j] = aa_cost_change[i*20+k] + aa_cost_change[k*20+j];
562
563 cout << "cost matrix by number of nt changes for TNT use" << endl;
564 cout << "smatrix =1 (aa_nt_changes)";
565 for (i = 0; i < 19; i++)
566 for (j = i+1; j < 20; j++)
567 cout << " " << symbols_protein[i] << "/" << symbols_protein[j] << " " << (int)aa_cost_change[i*20+j];
568 cout << ";" << endl;
569 cout << 20 << endl;
570 for (i = 0; i < 20; i++) {
571 aa_cost_change[i*20+i] = 0;
572 for (j = 0; j < 20; j++)
573 cout << (int)aa_cost_change[i*20+j] << " ";
574 cout << endl;
575 }
576
577 }
578
579 }
580
combineRateNTFreq()581 void ModelCodon::combineRateNTFreq() {
582 int i, j;
583 for (i = 0; i < num_states; i++) {
584 if (phylo_tree->aln->isStopCodon(i))
585 continue;
586 double *this_rate = &empirical_rates[i*num_states];
587 for (j = 0; j < num_states; j++) {
588 if (this_rate[j] == 0.0)
589 continue;
590
591 int codoni = phylo_tree->aln->codon_table[i];
592 int codonj = phylo_tree->aln->codon_table[j];
593 int nuc1, nuc2;
594
595 if ((nuc1=codoni/16) != (nuc2=codonj/16)) {
596 this_rate[j] *= ntfreq[nuc2];
597 }
598 if ((nuc1=(codoni%16)/4) != (nuc2=(codonj%16)/4)) {
599 this_rate[j] *= ntfreq[nuc2+4];
600 }
601 if ((nuc1=codoni%4) != (nuc2=codonj%4)) {
602 this_rate[j] *= ntfreq[nuc2+8];
603 }
604 }
605 }
606
607 }
608
609
readCodonModel(istream & in,bool reset_params)610 void ModelCodon::readCodonModel(istream &in, bool reset_params) {
611 int nrates = getNumRateEntries();
612
613 int i, j;
614 int nscodons = phylo_tree->aln->getNumNonstopCodons();
615
616 double * q = new double[nscodons*nscodons];
617 double *f = new double[nscodons];
618 for (i = 1; i < nscodons; i++) {
619 for (j = 0; j < i; j++) {
620 in >> q[i*nscodons+j];
621 //q[j*num_states+i] = q[i*num_states+j];
622 if (verbose_mode >= VB_MAX) cout << " " << q[i*nscodons+j];
623 }
624 if (verbose_mode >= VB_MAX) cout << endl;
625 }
626 for (i = 0; i < nscodons; i++)
627 in >> f[i];
628 StrVector codons;
629 codons.resize(nscodons);
630 IntVector state_map;
631 state_map.resize(nscodons);
632 for (i = 0; i < nscodons; i++) {
633 in >> codons[i];
634 if (codons[i].length() != 3)
635 outError("Input model has wrong codon format ", codons[i]);
636 int nt1 = phylo_tree->aln->convertState(codons[i][0], SEQ_DNA);
637 int nt2 = phylo_tree->aln->convertState(codons[i][1], SEQ_DNA);
638 int nt3 = phylo_tree->aln->convertState(codons[i][2], SEQ_DNA);
639 if (nt1 > 3 || nt2 > 3 || nt3 > 3)
640 outError("Wrong codon triplet ", codons[i]);
641 state_map[i] = phylo_tree->aln->non_stop_codon[nt1*16+nt2*4+nt3];
642 if (phylo_tree->aln->isStopCodon(state_map[i]) || state_map[i] == STATE_INVALID)
643 outError("Stop codon encountered");
644 if (verbose_mode >= VB_MAX)
645 cout << " " << codons[i] << " " << state_map[i];
646 }
647 if (verbose_mode >= VB_MAX) cout << endl;
648
649 //int row = 0, col = 1;
650 // since rates for codons is stored in lower-triangle, special treatment is needed
651 memset(empirical_rates, 0, nrates*sizeof(double));
652 memset(rates, 0, nrates*sizeof(double));
653 for (i = 1; i < nscodons; i++) {
654 for (j = 0; j < i; j++) {
655 int row = state_map[i], col = state_map[j];
656 if (row < col) {
657 int tmp = row;
658 row = col;
659 col = tmp;
660 }
661 // int id = col*(2*num_states-col-1)/2 + (row-col-1);
662 double qentry = q[i*nscodons+j];
663 int id = row*num_states+col;
664 ASSERT(id < nrates && id >= 0);
665 empirical_rates[id] = rates[id] = qentry;
666 id = col*num_states+row;
667 ASSERT(id < nrates && id >= 0);
668 empirical_rates[id] = rates[id] = qentry;
669 }
670 }
671 memset(state_freq, 0, num_states*sizeof(double));
672 for (i = 0; i < num_states; i++)
673 state_freq[i] = Params::getInstance().min_state_freq;
674 for (i = 0; i < nscodons; i++)
675 state_freq[state_map[i]] = f[i]-(num_states-nscodons)*Params::getInstance().min_state_freq/nscodons;
676
677 if (reset_params) {
678 fix_omega = fix_kappa = fix_kappa2 = true;
679 omega = kappa = kappa2 = 1.0;
680 }
681 delete [] f;
682 delete [] q;
683 }
684
readCodonModel(string & str,bool reset_params)685 void ModelCodon::readCodonModel(string &str, bool reset_params) {
686 try {
687 istringstream in(str);
688 readCodonModel(in, reset_params);
689 }
690 catch (const char *str) {
691 outError(str);
692 }
693 }
694
readCodonModelFile(const char * filename,bool reset_params)695 void ModelCodon::readCodonModelFile(const char *filename, bool reset_params) {
696 try {
697 ifstream in;
698 // set the failbit and badbit
699 in.exceptions(ios::failbit | ios::badbit);
700 in.open(filename);
701 // remove the failbit
702 in.exceptions(ios::badbit);
703
704 readCodonModel(in, reset_params);
705
706 in.clear();
707 // set the failbit again
708 in.exceptions(ios::failbit | ios::badbit);
709 in.close();
710 }
711 catch (const char *str) {
712 outError(str);
713 }
714 catch (...) {
715 outError(ERR_READ_INPUT, filename);
716 }
717 }
718
decomposeRateMatrix()719 void ModelCodon::decomposeRateMatrix() {
720 computeCodonRateMatrix();
721 ModelMarkov::decomposeRateMatrix();
722 }
723
computeCodonRateMatrix()724 void ModelCodon::computeCodonRateMatrix() {
725 // if (num_params == 0)
726 // return; // do nothing for empirical codon model
727
728 switch (codon_kappa_style) {
729 case CK_ONE_KAPPA:
730 computeCodonRateMatrix_1KAPPA();
731 break;
732 case CK_ONE_KAPPA_TS:
733 computeCodonRateMatrix_1KAPPATS();
734 break;
735 case CK_ONE_KAPPA_TV:
736 computeCodonRateMatrix_1KAPPATV();
737 break;
738 case CK_TWO_KAPPA:
739 computeCodonRateMatrix_2KAPPA();
740 break;
741 }
742 }
743
computeCodonRateMatrix_1KAPPA()744 void ModelCodon::computeCodonRateMatrix_1KAPPA() {
745 int nrates = getNumRateEntries();
746 memcpy(rates, empirical_rates, nrates*sizeof(double));
747 if (omega == 1.0 && kappa == 1.0)
748 return; // do nothing
749
750 int i, j;
751 double omega_kappa = omega*kappa;
752
753 for (i = 0; i < num_states; i++) {
754 double *this_rate = &rates[i*num_states];
755 int *this_rate_attr = &rate_attr[i*num_states];
756 if (phylo_tree->aln->isStopCodon(i)) {
757 continue;
758 }
759 for (j = 0; j < num_states; j++) {
760 if (this_rate[j] == 0.0) continue;
761 int attr = this_rate_attr[j];
762 if (attr & CA_SYNONYMOUS) { // synonymous
763 if (attr & CA_TRANSITION) // transition
764 this_rate[j] *= kappa;
765 } else if (attr & CA_NONSYNONYMOUS) { // non-synomyous
766 if (attr & CA_TRANSITION) // transition
767 this_rate[j] *= omega_kappa;
768 else // transversion
769 this_rate[j] *= omega;
770 }
771 }
772 }
773 }
774
computeCodonRateMatrix_1KAPPATS()775 void ModelCodon::computeCodonRateMatrix_1KAPPATS() {
776 int nrates = getNumRateEntries();
777 memcpy(rates, empirical_rates, nrates*sizeof(double));
778
779 int i, j;
780 double kappa_pow[] = {1.0, kappa, kappa*kappa, kappa*kappa*kappa};
781 double omega_kappa_pow[] = {omega, omega*kappa, omega*kappa*kappa, omega*kappa*kappa*kappa};
782
783 for (i = 0; i < num_states; i++) {
784 double *this_rate = &rates[i*num_states];
785 int *this_rate_attr = &rate_attr[i*num_states];
786 if (phylo_tree->aln->isStopCodon(i)) {
787 continue;
788 }
789 for (j = 0; j < num_states; j++) {
790 int attr = this_rate_attr[j];
791 if (this_rate[j] == 0.0) continue;
792 if (attr & CA_SYNONYMOUS) { // synonymous
793 int num = ((attr & CA_TRANSITION_1NT) != 0) + ((attr & CA_TRANSITION_2NT) != 0) + ((attr & CA_TRANSITION_3NT) != 0);
794 this_rate[j] *= kappa_pow[num];
795 } else if (attr & CA_NONSYNONYMOUS) { // non-synomyous
796 int num = ((attr & CA_TRANSITION_1NT) != 0) + ((attr & CA_TRANSITION_2NT) != 0) + ((attr & CA_TRANSITION_3NT) != 0);
797 this_rate[j] *= omega_kappa_pow[num];
798 }
799 }
800 }
801 }
802
computeCodonRateMatrix_1KAPPATV()803 void ModelCodon::computeCodonRateMatrix_1KAPPATV() {
804 int nrates = getNumRateEntries();
805 memcpy(rates, empirical_rates, nrates*sizeof(double));
806
807 int i, j;
808 double kappa_pow[] = {1.0, kappa, kappa*kappa, kappa*kappa*kappa};
809 double omega_kappa_pow[] = {omega, omega*kappa, omega*kappa*kappa, omega*kappa*kappa*kappa};
810
811 for (i = 0; i < num_states; i++) {
812 double *this_rate = &rates[i*num_states];
813 int *this_rate_attr = &rate_attr[i*num_states];
814 if (phylo_tree->aln->isStopCodon(i)) {
815 continue;
816 }
817 for (j = 0; j < num_states; j++) {
818 int attr = this_rate_attr[j];
819 if (this_rate[j] == 0.0) continue;
820 if (attr & CA_SYNONYMOUS) { // synonymous
821 int num = ((attr & CA_TRANSVERSION_1NT) != 0) + ((attr & CA_TRANSVERSION_2NT) != 0) + ((attr & CA_TRANSVERSION_3NT) != 0);
822 this_rate[j] *= kappa_pow[num];
823 } else if (attr & CA_NONSYNONYMOUS) { // non-synomyous
824 int num = ((attr & CA_TRANSVERSION_1NT) != 0) + ((attr & CA_TRANSVERSION_2NT) != 0) + ((attr & CA_TRANSVERSION_3NT) != 0);
825 this_rate[j] *= omega_kappa_pow[num];
826 }
827 }
828 }
829 }
830
computeCodonRateMatrix_2KAPPA()831 void ModelCodon::computeCodonRateMatrix_2KAPPA() {
832 int nrates = getNumRateEntries();
833 memcpy(rates, empirical_rates, nrates*sizeof(double));
834
835 int i, j;
836 double kappa_pow[] = {1.0, kappa, kappa*kappa, kappa*kappa*kappa};
837 double omega_kappa_pow[] = {omega, omega*kappa, omega*kappa*kappa, omega*kappa*kappa*kappa};
838 double kappa2_pow[] = {1.0, kappa2, kappa2*kappa2, kappa2*kappa2*kappa2};
839
840 for (i = 0; i < num_states; i++) {
841 double *this_rate = &rates[i*num_states];
842 int *this_rate_attr = &rate_attr[i*num_states];
843 if (phylo_tree->aln->isStopCodon(i)) {
844 continue;
845 }
846 for (j = 0; j < num_states; j++) {
847 int attr = this_rate_attr[j];
848 if (this_rate[j] == 0.0) continue;
849 if (attr & CA_SYNONYMOUS) { // synonymous
850 int numts = ((attr & CA_TRANSITION_1NT) != 0) + ((attr & CA_TRANSITION_2NT) != 0) + ((attr & CA_TRANSITION_3NT) != 0);
851 int numtv = ((attr & CA_TRANSVERSION_1NT) != 0) + ((attr & CA_TRANSVERSION_2NT) != 0) + ((attr & CA_TRANSVERSION_3NT) != 0);
852 this_rate[j] *= kappa_pow[numts]*kappa2_pow[numtv];
853 } else if (attr & CA_NONSYNONYMOUS) { // non-synomyous
854 int numts = ((attr & CA_TRANSITION_1NT) != 0) + ((attr & CA_TRANSITION_2NT) != 0) + ((attr & CA_TRANSITION_3NT) != 0);
855 int numtv = ((attr & CA_TRANSVERSION_1NT) != 0) + ((attr & CA_TRANSVERSION_2NT) != 0) + ((attr & CA_TRANSVERSION_3NT) != 0);
856 this_rate[j] *= omega_kappa_pow[numts]*kappa2_pow[numtv];
857 }
858 }
859 }
860 }
861
computeEmpiricalOmega()862 double ModelCodon::computeEmpiricalOmega() {
863 double dn = 0.0, ds = 0.0;
864 int i, j;
865 if (ignore_state_freq) {
866 for (i = 0; i < num_states; i++) {
867 if (phylo_tree->aln->isStopCodon(i))
868 continue;
869 double *this_rate = &rates[i*num_states];
870 int *this_rate_attr = &rate_attr[i*num_states];
871 for (j = 0; j < num_states; j++)
872 if (this_rate_attr[j] & CA_NONSYNONYMOUS)
873 dn += state_freq[i]*this_rate[j];
874 else
875 ds += state_freq[i]*this_rate[j];
876 }
877 } else {
878 for (i = 0; i < num_states; i++) {
879 if (phylo_tree->aln->isStopCodon(i))
880 continue;
881 double *this_rate = &rates[i*num_states];
882 int *this_rate_attr = &rate_attr[i*num_states];
883 for (j = 0; j < num_states; j++)
884 if (this_rate_attr[j] & CA_NONSYNONYMOUS)
885 dn += state_freq[i]*state_freq[j]*this_rate[j];
886 else
887 ds += state_freq[i]*state_freq[j]*this_rate[j];
888 }
889 }
890 return (dn/ds)*(0.21/0.79);
891 }
892
893
894
getVariables(double * variables)895 bool ModelCodon::getVariables(double *variables) {
896 int j;
897 bool changed = false;
898 if (num_params > 0) {
899 j = 1;
900 if (!fix_omega) {
901 changed |= (omega != variables[j]);
902 omega = variables[j++];
903 }
904 if (!fix_kappa) {
905 changed |= (kappa != variables[j]);
906 kappa = variables[j++];
907 }
908 if (!fix_kappa2) {
909 changed |= (kappa2 != variables[j]);
910 kappa2 = variables[j++];
911 }
912 ASSERT(j == num_params+1);
913 }
914 if (freq_type == FREQ_ESTIMATE) {
915 // 2015-09-07: relax the sum of state_freq to be 1, this will be done at the end of optimization
916 int ndim = getNDim();
917 changed |= memcmpcpy(state_freq, variables+(ndim-num_states+2), (num_states-1)*sizeof(double));
918 // double sum = 0;
919 // for (i = 0; i < num_states-1; i++)
920 // sum += state_freq[i];
921 // state_freq[num_states-1] = 1.0 - sum;
922
923 // BUG FIX 2015.08.28
924 // int nrate = getNDim();
925 // if (freq_type == FREQ_ESTIMATE) nrate -= (num_states-1);
926 // double sum = 1.0;
927 //// int i, j;
928 // for (i = 1; i < num_states; i++)
929 // sum += variables[nrate+i];
930 // for (i = 0, j = 1; i < num_states; i++)
931 // if (i != highest_freq_state) {
932 // state_freq[i] = variables[nrate+j] / sum;
933 // j++;
934 // }
935 // state_freq[highest_freq_state] = 1.0/sum;
936 }
937 return changed;
938 }
939
setVariables(double * variables)940 void ModelCodon::setVariables(double *variables) {
941 int j;
942 if (num_params > 0) {
943 j = 1;
944 if (!fix_omega)
945 variables[j++] = omega;
946 if (!fix_kappa)
947 variables[j++] = kappa;
948 if (!fix_kappa2)
949 variables[j++] = kappa2;
950
951 ASSERT(j == num_params+1);
952 }
953 if (freq_type == FREQ_ESTIMATE) {
954 // 2015-09-07: relax the sum of state_freq to be 1, this will be done at the end of optimization
955 int ndim = getNDim();
956 memcpy(variables+(ndim-num_states+2), state_freq, (num_states-1)*sizeof(double));
957
958 // BUG FIX 2015.08.28
959 // int nrate = getNDim();
960 // if (freq_type == FREQ_ESTIMATE) nrate -= (num_states-1);
961 // int i, j;
962 // for (i = 0, j = 1; i < num_states; i++)
963 // if (i != highest_freq_state) {
964 // variables[nrate+j] = state_freq[i] / state_freq[highest_freq_state];
965 // j++;
966 // }
967 }
968 }
969
setBounds(double * lower_bound,double * upper_bound,bool * bound_check)970 void ModelCodon::setBounds(double *lower_bound, double *upper_bound, bool *bound_check) {
971 int i, ndim = getNDim();
972
973 for (i = 1; i <= ndim; i++) {
974 //cout << variables[i] << endl;
975 lower_bound[i] = MIN_OMEGA_KAPPA;
976 upper_bound[i] = MAX_OMEGA_KAPPA;
977 bound_check[i] = false;
978 }
979
980 if (freq_type == FREQ_ESTIMATE) {
981 for (i = ndim-num_states+2; i <= ndim; i++) {
982 // lower_bound[i] = MIN_FREQUENCY/state_freq[highest_freq_state];
983 // upper_bound[i] = state_freq[highest_freq_state]/MIN_FREQUENCY;
984 lower_bound[i] = Params::getInstance().min_state_freq;
985 // upper_bound[i] = 100.0;
986 upper_bound[i] = 1.0;
987 bound_check[i] = false;
988 }
989 }
990 }
991
optimizeParameters(double gradient_epsilon)992 double ModelCodon::optimizeParameters(double gradient_epsilon) {
993
994 if (fixed_parameters)
995 return 0.0;
996
997 int ndim = getNDim();
998
999 // return if nothing to be optimized
1000 if (ndim == 0) return 0.0;
1001
1002 if (verbose_mode >= VB_MAX)
1003 cout << "Optimizing " << name << " model parameters..." << endl;
1004
1005
1006 double *variables = new double[ndim+1];
1007 double *upper_bound = new double[ndim+1];
1008 double *lower_bound = new double[ndim+1];
1009 bool *bound_check = new bool[ndim+1];
1010 double score;
1011
1012 for (int i = 0; i < num_states; i++)
1013 if (state_freq[i] > state_freq[highest_freq_state])
1014 highest_freq_state = i;
1015
1016 // by BFGS algorithm
1017 setVariables(variables);
1018 setBounds(lower_bound, upper_bound, bound_check);
1019 if (phylo_tree->params->optimize_alg.find("BFGS-B") == string::npos)
1020 score = -minimizeMultiDimen(variables, ndim, lower_bound, upper_bound, bound_check, max(gradient_epsilon, TOL_RATE));
1021 else
1022 score = -L_BFGS_B(ndim, variables+1, lower_bound+1, upper_bound+1, max(gradient_epsilon, TOL_RATE));
1023
1024 bool changed = getVariables(variables);
1025 // BQM 2015-09-07: normalize state_freq
1026 if (freq_type == FREQ_ESTIMATE) {
1027 scaleStateFreq(true);
1028 changed = true;
1029 }
1030 if (changed) {
1031 decomposeRateMatrix();
1032 phylo_tree->clearAllPartialLH();
1033 score = phylo_tree->computeLikelihood();
1034 }
1035
1036 delete [] bound_check;
1037 delete [] lower_bound;
1038 delete [] upper_bound;
1039 delete [] variables;
1040
1041 return score;
1042 }
1043
1044
writeInfo(ostream & out)1045 void ModelCodon::writeInfo(ostream &out) {
1046 if (name.find('_') == string::npos)
1047 out << "Nonsynonymous/synonymous ratio (omega): " << omega << endl;
1048 else
1049 out << "Empirical nonsynonymous/synonymous ratio (omega_E): " << computeEmpiricalOmega() << endl;
1050 out << "Transition/transversion ratio (kappa): " << kappa << endl;
1051 if (codon_kappa_style == CK_TWO_KAPPA)
1052 out << "Transition/transversion ratio 2 (kappa2): " << kappa2 << endl;
1053 }
1054
1055