?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
modified: SpatialOmicsCoord.py
[GalaxyCodeBases.git] / c_cpp / etc / KaKs_Calculator / src / GY94.cpp
blob7470b283c0c6d2aafbe9a46af20d289131b53539
1 /************************************************************
2 * Copyright (c) 2005, BGI of Chinese Academy of Sciences
3 * All rights reserved.
4
5 * Filename: GY94.cpp
6 * Abstract: Definition of GY94 class.
7
8 * Version: 1.0
9 * Author: Zhang Zhang (zhanghzhang@genomics.org.cn)
10 * Date: Oct.2, 2005
11
12 * Note: Source codes are taken from codeml.c in PAML.
13
14 References:
15 Goldman N, Yang Z (1994) A codon-based model of nucleotide
16 substitution for protein-coding DNA sequences. Mol. Biol.
17 Evol. 11:725-736.
18 *************************************************************/
19 #include "GY94.h"
20 //#include <malloc.h>
21 //#include<MALLOC>
22
23 int GeneticCode[][64] =
24 {{13,13,10,10,15,15,15,15,18,18,-1,-1, 4, 4,-1,17,
25 10,10,10,10,14,14,14,14, 8, 8, 5, 5, 1, 1, 1, 1,
26 9, 9, 9,12,16,16,16,16, 2, 2,11,11,15,15, 1, 1,
27 19,19,19,19, 0, 0, 0, 0, 3, 3, 6, 6, 7, 7, 7, 7}, /* 0:universal */
28
29 {13,13,10,10,15,15,15,15,18,18,-1,-1, 4, 4,17,17,
30 10,10,10,10,14,14,14,14, 8, 8, 5, 5, 1, 1, 1, 1,
31 9, 9,12,12,16,16,16,16, 2, 2,11,11,15,15,-1,-1,
32 19,19,19,19, 0, 0, 0, 0, 3, 3, 6, 6, 7, 7, 7, 7}, /* 1:vertebrate mt.*/
33
34 {13,13,10,10,15,15,15,15,18,18,-1,-1, 4, 4,17,17,
35 16,16,16,16,14,14,14,14, 8, 8, 5, 5, 1, 1, 1, 1,
36 9, 9,12,12,16,16,16,16, 2, 2,11,11,15,15, 1, 1,
37 19,19,19,19, 0, 0, 0, 0, 3, 3, 6, 6, 7, 7, 7, 7}, /* 2:yeast mt. */
38
39 {13,13,10,10,15,15,15,15,18,18,-1,-1, 4, 4,17,17,
40 10,10,10,10,14,14,14,14, 8, 8, 5, 5, 1, 1, 1, 1,
41 9, 9, 9,12,16,16,16,16, 2, 2,11,11,15,15, 1, 1,
42 19,19,19,19, 0, 0, 0, 0, 3, 3, 6, 6, 7, 7, 7, 7}, /* 3:mold mt. */
43
44 {13,13,10,10,15,15,15,15,18,18,-1,-1, 4, 4,17,17,
45 10,10,10,10,14,14,14,14, 8, 8, 5, 5, 1, 1, 1, 1,
46 9, 9,12,12,16,16,16,16, 2, 2,11,11,15,15,15,15,
47 19,19,19,19, 0, 0, 0, 0, 3, 3, 6, 6, 7, 7, 7, 7}, /* 4:invertebrate mt. */
48
49 {13,13,10,10,15,15,15,15,18,18, 5, 5, 4, 4,-1,17,
50 10,10,10,10,14,14,14,14, 8, 8, 5, 5, 1, 1, 1, 1,
51 9, 9, 9,12,16,16,16,16, 2, 2,11,11,15,15, 1, 1,
52 19,19,19,19, 0, 0, 0, 0, 3, 3, 6, 6, 7, 7, 7, 7}, /* 5:ciliate nuclear*/
53
54 {13,13,10,10,15,15,15,15,18,18,-1,-1, 4, 4,17,17,
55 10,10,10,10,14,14,14,14, 8, 8, 5, 5, 1, 1, 1, 1,
56 9, 9, 9,12,16,16,16,16, 2, 2, 2,11,15,15,15,15,
57 19,19,19,19, 0, 0, 0, 0, 3, 3, 6, 6, 7, 7, 7, 7}, /* 6:echinoderm mt.*/
58
59 {13,13,10,10,15,15,15,15,18,18,-1,-1, 4, 4, 4,17,
60 10,10,10,10,14,14,14,14, 8, 8, 5, 5, 1, 1, 1, 1,
61 9, 9, 9,12,16,16,16,16, 2, 2,11,11,15,15, 1, 1,
62 19,19,19,19, 0, 0, 0, 0, 3, 3, 6, 6, 7, 7, 7, 7}, /* 7:euplotid mt. */
63
64 {13,13,10,10,15,15,15,15,18,18,-1,-1, 4, 4,-1,17,
65 10,10,10,15,14,14,14,14, 8, 8, 5, 5, 1, 1, 1, 1,
66 9, 9, 9,12,16,16,16,16, 2, 2,11,11,15,15, 1, 1,
67 19,19,19,19, 0, 0, 0, 0, 3, 3, 6, 6, 7, 7, 7, 7},
68 /* 8:alternative yeast nu.*/
69
70 {13,13,10,10,15,15,15,15,18,18,-1,-1, 4, 4,17,17,
71 10,10,10,10,14,14,14,14, 8, 8, 5, 5, 1, 1, 1, 1,
72 9, 9,12,12,16,16,16,16, 2, 2,11,11,15,15, 7, 7,
73 19,19,19,19, 0, 0, 0, 0, 3, 3, 6, 6, 7, 7, 7, 7}, /* 9:ascidian mt. */
74
75 {13,13,10,10,15,15,15,15,18,18,-1, 5, 4, 4,-1,17,
76 10,10,10,10,14,14,14,14, 8, 8, 5, 5, 1, 1, 1, 1,
77 9, 9, 9,12,16,16,16,16, 2, 2,11,11,15,15, 1, 1,
78 19,19,19,19, 0, 0, 0, 0, 3, 3, 6, 6, 7, 7, 7, 7}, /* 10:blepharisma nu.*/
79
80 { 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4,
81 5, 5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8,
82 9, 9, 9, 9,10,10,10,10,11,11,11,11,12,12,12,12,
83 13,13,13,13,14,14,14,14,15,15,15,15,16,16,16,16} /* 11:Ziheng's regular code */
84 }; /* GeneticCode[icode][#codon] */
85
86
87 /********************************************
88 * Function: construction function
89 * Input Parameter: a candidate model
90 * Output:
91 * Return Value:
92
93 * Note: (Equal, Unequal: substitution rates)
94 JC, F81: rTC==rAG =rTA==rCG==rTG==rCA
95 K2P, HKY: rTC==rAG!=rTA==rCG==rTG==rCA
96 TNEF,TN: rTC!=rAG!=rTA==rCG==rTG==rCA
97 K3P, K3PUF: rTC==rAG!=rTA==rCG!=rTG==rCA
98 TIMEF, TIM: rTC!=rAG!=rTA==rCG!=rTG==rCA
99 TVMEF, TVM: rTC==rAG!=rTA!=rCG!=rTG!=rCA
100 SYM, GTR: rTC!=rAG!=rTA!=rCG!=rTG!=rCA
101 *********************************************/
102
103 //Constructor
104 GY94::GY94() {
105 }
106
107 GY94::GY94(string NulModel) {
108
109 name = "GY-"+NulModel;
110 Iround=0;
111 SIZEp=0;
112 Small_Diff=1e-6;
113 w_rndu=123456757;
114 com.ns = 2;
115 lnL = 0.0;
116
117 com.icode = genetic_code-1;
118 if (com.icode>11) com.icode = 0;
119 Nsensecodon = getNumNonsense(com.icode);
120 com.ncode = 64 - Nsensecodon;
121
122 //NulModel is set among a set of candidate models.
123 model = NulModel;
124
125 if (model=="JC" || model=="F81") com.nkappa = 0;
126 else if (model=="K2P" || model=="HKY") com.nkappa = 1;
127 else if (model=="TNEF"|| model=="TN") com.nkappa = 2;
128 else if (model=="K3P" || model=="K3PUF") com.nkappa = 2;
129 else if (model=="TIMEF" || model=="TIM") com.nkappa = 3;
130 else if (model=="TVMEF" || model=="TVM") com.nkappa = 4;
131 else if (model=="SYM" || model=="GTR") com.nkappa = 5;
132
133 //parameters' number: plus another two parameters (Ka/Ks and t)
134 com.np = 2 + com.nkappa;
135
136 }
137
138 GY94::~GY94() {
139
140 }
141
142
143 double GY94::rndu(void) {
144 w_rndu = w_rndu*69069+1;
145 return ldexp((double)w_rndu, -32);
146 }
147
148
149 /* x[i]=x0[i] + t*p[i] */
150 double GY94::fun_ls (double t, double x0[], double p[], double x[], int n) {
151 int i;
152 for (i=0; i<n; i++) x[i]=x0[i] + t*p[i];
153 return( lfun2dSdN(x,n) );
154 }
155
156
157
158 int GY94::PatternWeight () {
159
160 int *fpatti, lst=com.ls,h,ht,j,k=-1;
161 int gap=3;
162 int nb=1;
163 char *zt[2], b1[2]; /* b[][] data at site h */
164 double nc = 65;
165
166 //fpatti=(int*)malloc(lst*sizeof(int));
167 fpatti= new int[lst];
168
169 for(h=0; h<lst; h++) fpatti[h]=0;
170
171 for(j=0;j<com.ns;j++) {
172 zt[j]=(char*)malloc(lst*sizeof(char));
173 //zt[j]= new char[lst];
174 }
175
176 com.npatt=0;
177
178 for(h=0; h<com.ls; h++) {
179
180 b1[0]=com.z[0][h];
181 b1[1]=com.z[1][h];
182
183 //cout<<(char)b1[1]<<endl;
184
185 for(ht=0; ht<com.npatt; ht++) {
186 if(b1[0]==zt[0][ht] && b1[1]==zt[1][ht])
187 break;
188 }
189 if(ht==com.npatt) {
190 zt[0][com.npatt]=b1[0];
191 zt[1][com.npatt]=b1[1];
192
193 ht=com.npatt++;
194 }
195
196 fpatti[ht]++;
197 } //for (h)
198
199
200 //com.fpatt=(double*)malloc(com.npatt*sizeof(double));
201
202 for(h=0;h<com.npatt;h++) {
203 com.fpatt[h]=fpatti[h];
204 }
205
206 //free(fpatti);
207 delete []fpatti;
208
209 for(j=0; j<com.ns; j++) {
210 com.z[j]=(char*)realloc(zt[j],com.npatt*sizeof(char));
211 }
212
213 return 0;
214 }
215
216 //Preprocess in preparation for estimation
217 int GY94::preProcess(const char* seq1, const char* seq2) {
218
219 int i;
220
221 // com.space=NULL;
222 com.kappa=2; com.omega = 0.4;
223
224 setmark_61_64();
225
226 com.ls = strlen(seq1);
227
228 for(i=0, snp=0; i<com.ls; i++)
229 if (seq1[i]!=seq2[i]) snp++;
230
231 for(i=0; i<com.ns; i++) {
232 com.z[i]=new char[com.ls+1];
233 }
234
235 strcpy(com.z[0], seq1);
236 strcpy(com.z[1], seq2);
237
238 com.ls /=3;
239
240 EncodeSeqs();
241 PatternWeight();
242
243 return 0;
244 }
245
246 int GY94::setmark_61_64 (void) {
247
248 int i, *code=GeneticCode[com.icode];
249 //int c[3],aa0,aa1, by[3]={16,4,1};
250 //double nSilent, nStop, nRepl;
251
252 Nsensecodon=0;
253 for (i=0; i<64; i++) {
254 if (code[i]==-1) {
255 FROM64[i]=-1;
256 }
257 else {
258 FROM61[Nsensecodon] = i;
259 FROM64[i] = Nsensecodon++;
260 }
261 }
262
263 com.ncode=Nsensecodon;
264
265 return 0;
266 }
267
268 //Convert T,C,A,G to 0,1,2,3
269 int GY94::transform(char *z, int ls) {
270
271 int i;
272 char *p;
273
274 for (i=0,p=z; i<ls; i++,p++) {
275 *p=(char)convertChar(*p);
276 }
277
278 return 1;
279 }
280
281 void GY94::EncodeSeqs(void) {
282
283 int j,h,k, b[3];
284
285 //T C A G to 0 1 2 3
286 for(j=0; j<com.ns; j++) {
287 transform(com.z[j],com.ls*3);
288 }
289
290 //encode to codon 0-64
291 for(j=0; j<com.ns; j++) {
292 for(h=0; h<com.ls; h++) {
293 b[0]=com.z[j][h*3]; b[1]=com.z[j][h*3+1]; b[2]=com.z[j][h*3+2];
294 k=b[0]*16+b[1]*4+b[2];
295 com.z[j][h]=(char)FROM64[k];
296 }
297 }
298 }
299
300 int GY94::GetCodonFreqs(double pi[]) {
301 int n=com.ncode, i,j,ic,b[3];
302 double fb3x4[12], fb4[4];
303 int flag[CODON];
304
305 for(i=0, initArray(flag, 64); i<str1.length(); i+=3) {
306 j = getID(str1.substr(i,3));
307 if (getAminoAcid(j)!='!') flag[j] = 1;
308
309 j = getID(str2.substr(i,3));
310 if (getAminoAcid(j)!='!') flag[j] = 1;
311 }
312
313 //Whether sequences are long enough
314 if (sumArray(flag, CODON)==com.ncode) {
315 return 0;
316 }
317
318 //Codon frequency is estimated from nucleotide frequency(fb3x4) at the three positions
319 for(i=0,initArray(fb3x4,12),initArray(fb4,4); i<n; i++) {
320 ic=FROM61[i];
321 b[0]=ic/16;
322 b[1]=(ic/4)%4;
323 b[2]=ic%4;
324 for(j=0; j<3; j++) {
325 fb3x4[j*4+b[j]] += pi[i];
326 fb4[b[j]] += pi[i]/3.;
327 }
328 }
329
330 //use nul frequencies to get codon frequency
331 for (i=0; i<n; i++) {
332 ic=FROM61[i]; b[0]=ic/16; b[1]=(ic/4)%4; b[2]=ic%4;
333 pi[i]=fb3x4[b[0]]*fb3x4[4+b[1]]*fb3x4[8+b[2]];
334 }
335
336 scaleArray(1./sumArray(pi,n), pi, n);
337
338 return 0;
339 }
340
341 int GY94::gradientB(int n, double x[], double f0, double g[], double space[], int xmark[]) {
342 /* f0=fun(x) is always provided.
343 xmark=0: central; 1: upper; -1: down
344 */
345 int i,j;
346 double *x0=space, *x1=space+n, eh0=Small_Diff, eh; /* eh0=1e-6 || 1e-7 */
347
348 for(i=0; i<n; i++) {
349 eh=eh0*(fabs(x[i])+1);
350 if (xmark[i]==0 && SIZEp<1) { //central
351 for (j=0; j<n; j++) x0[j]=x1[j]=x[j];
352 eh=pow(eh,.67); x0[i]-=eh; x1[i]+=eh;
353 g[i]=(lfun2dSdN(x1,n)-(lfun2dSdN)(x0,n))/(eh*2.0);
354 }
355 else {//forward or backward
356 for(j=0; j<n; j++) x1[j]=x[j];
357 if (xmark[i]) eh*=-xmark[i];
358 x1[i] += eh;
359 g[i]=(lfun2dSdN(x1,n)-f0)/eh;
360 }
361 }
362 return(0);
363 }
364
365
366 double GY94::LineSearch2 (double *f, double x0[], double p[], double step, double limit, double e, double space[], int n) {
367 /* linear search using quadratic interpolation
368 from x0[] in the direction of p[],
369 x = x0 + a*p a ~(0,limit)
370 returns (a). *f: f(x0) for input and f(x) for output
371
372 x0[n] x[n] p[n] space[n]
373
374 adapted from Wolfe M. A. 1978. Numerical methods for unconstrained
375 optimization: An introduction. Van Nostrand Reinhold Company, New York.
376 pp. 62-73.
377 step is used to find the bracket and is increased or reduced as necessary,
378 and is not terribly important.
379 */
380 int ii=0, maxround=10, status, nsymb=0;
381 double *x=space, factor=4, small=1e-6, smallgapa=0.2;
382 double a0,a1,a2,a3,a4=-1,a5,a6, f0,f1,f2,f3,f4=-1,f5,f6;
383
384 a0=a1=0; f1=f0=*f;
385 a2=a0+step;
386 f2=fun_ls(a2, x0,p,x,n);
387 if (f2>f1) {
388 for (; ;) {
389 step/=factor;
390 if (step<small) return (0);
391 a3=a2; f3=f2;
392 a2=a0+step; f2=fun_ls(a2, x0,p,x,n);
393 if (f2<=f1) break;
394 }
395 }
396 else {
397 for (; ;) {
398 step*=factor;
399 if (step>limit) step=limit;
400 a3=a0+step; f3=fun_ls(a3, x0,p,x,n);
401
402 //obtain a bracket
403 if (f3>=f2) { //a1<a2<a3 and f1>f2<f3
404 break;
405 }
406
407 a1=a2; f1=f2; a2=a3; f2=f3;
408 if (step>=limit) {
409 *f=f3; return(a3);
410 }
411 }
412 }
413
414 /* iteration by quadratic interpolation, fig 2.2.9-10 (pp 71-71) */
415 for (ii=0; ii<maxround; ii++) { //a4 is the minimum from the parabola over (a1,a2,a3)
416
417 //2.2.4
418 a4 = (a2-a3)*f1+(a3-a1)*f2+(a1-a2)*f3;
419 if(fabs(a4)>1e-100) {
420 //2.2.3
421 a4 = ((a2*a2-a3*a3)*f1+(a3*a3-a1*a1)*f2+(a1*a1-a2*a2)*f3)/(2*a4);
422 }
423 if (a4>a3 || a4<a1) {//out of range: whether a1<a4<a3
424 a4=(a1+a2)/2;
425 }
426 else {
427 if((a4<=a2 && a2-a4>smallgapa*(a2-a1)) || (a4>a2 && a4-a2>smallgapa*(a3-a2)))
428 status='Y';
429 else
430 status='C';
431 }
432 f4 = fun_ls(a4, x0,p,x,n);
433
434 if (fabs(f2-f4)<e*(1+fabs(f2))) {
435 break;
436 }
437
438 if (a4<=a2) { /* fig 2.2.10 */
439 if (a2-a4>smallgapa*(a2-a1)) {
440 if (f4<=f2) { a3=a2; a2=a4; f3=f2; f2=f4; }
441 else { a1=a4; f1=f4; }
442 }
443 else {
444 if (f4>f2) {
445 a5=(a2+a3)/2; f5=fun_ls(a5, x0,p,x,n);
446 if (f5>f2) { a1=a4; a3=a5; f1=f4; f3=f5; }
447 else { a1=a2; a2=a5; f1=f2; f2=f5; }
448 }
449 else {
450 a5=(a1+a4)/2; f5=fun_ls(a5, x0,p,x,n);
451 if (f5>=f4)
452 { a3=a2; a2=a4; a1=a5; f3=f2; f2=f4; f1=f5; }
453 else {
454 a6=(a1+a5)/2; f6=fun_ls(a6, x0,p,x,n);
455 if (f6>f5)
456 { a1=a6; a2=a5; a3=a4; f1=f6; f2=f5; f3=f4; }
457 else { a2=a6; a3=a5; f2=f6; f3=f5; }
458 }
459 }
460 }
461 }
462 else { /* fig 2.2.9 */
463 if (a4-a2>smallgapa*(a3-a2)) {
464 if (f2>=f4) { a1=a2; a2=a4; f1=f2; f2=f4; }
465 else { a3=a4; f3=f4; }
466 }
467 else {
468 if (f4>f2) {
469 a5=(a1+a2)/2; f5=fun_ls(a5, x0,p,x,n);
470 if (f5>f2) { a1=a5; a3=a4; f1=f5; f3=f4; }
471 else { a3=a2; a2=a5; f3=f2; f2=f5; }
472 }
473 else {
474 a5=(a3+a4)/2; f5=fun_ls(a5, x0,p,x,n);
475 if (f5>=f4)
476 { a1=a2; a2=a4; a3=a5; f1=f2; f2=f4; f3=f5; }
477 else {
478 a6=(a3+a5)/2; f6=fun_ls(a6, x0,p,x,n);
479 if (f6>f5)
480 { a1=a4; a2=a5; a3=a6; f1=f4; f2=f5; f3=f6; }
481 else { a1=a5; a2=a6; f1=f5; f2=f6; }
482 }
483 }
484 }
485 }
486 }
487
488 if (f2>f0 && f4>f0) a4=0;
489 if (f2<=f4) { *f=f2; a4=a2; }
490 else *f=f4;
491
492 return a4;
493 }
494
495 double GY94::distance(double x[], double y[], int n) {
496
497 int i;
498 double t=0;
499
500 for (i=0; i<n; i++) t += square(x[i]-y[i]);
501
502 return sqrt(t);
503 }
504
505 int GY94::H_end (double x0[], double x1[], double f0, double f1, double e1, double e2, int n) {
506 /* Himmelblau termination rule. return 1 for stop, 0 otherwise.
507 */
508 double r;
509
510 if((r=norm(x0,n))<e2) r=1;
511 r*=e1;
512 if(distance(x1,x0,n)>=r)
513 return 0;
514
515 r=fabs(f0);
516 if(r<e2) r=1;
517 r*=e1;
518 if (fabs(f1-f0)>=r)
519 return 0;
520
521 return 1;
522 }
523
524 int GY94::ming2 (double *f, double x[], double xb[][2], double space[], double e, int n) {
525 /* n-variate minimization with bounds using the BFGS algorithm
526 g0[n] g[n] p[n] x0[n] y[n] s[n] z[n] H[n*n] C[n*n] tv[2*n]
527 xmark[n],ix[n]
528 Size of space should be (check carefully?)
529 #define spaceming2(n) ((n)*((n)*2+9+2)*sizeof(double))
530 nfree: # free variables
531 xmark[i]=0 for inside space; -1 for lower boundary; 1 for upper boundary.
532 x[] has initial values at input and returns the estimates in return.
533 ix[i] specifies the i-th free parameter
534 */
535 int i,j, i1,i2,it, maxround=2000, fail=0, *xmark, *ix, nfree;
536 int Ngoodtimes=2, goodtimes=0;
537 double small=1e-6, sizep0=0;
538 double f0, *g0, *g, *p, *x0, *y, *s, *z, *H, *C, *tv;
539 double w,v, alpha, am, h, maxstep=8;
540
541 if(n==0) return 0;
542
543 g0=space; g=g0+n; p=g+n; x0=p+n;
544 y=x0+n; s=y+n; z=s+n; H=z+n; C=H+n*n; tv=C+n*n;
545
546 xmark = (int*)(tv+2*n);
547 ix = xmark+n;
548
549 for(i=0; i<n; i++) {
550 xmark[i]=0;
551 ix[i]=i;
552 }
553
554 for(i=0,nfree=0;i<n;i++) {
555 if(x[i]<=xb[i][0]) {
556 x[i]=xb[i][0];
557 xmark[i]=-1;
558 continue;
559 }
560 if(x[i]>=xb[i][1]) {
561 x[i]=xb[i][1];
562 xmark[i]= 1;
563 continue;
564 }
565
566 ix[nfree++] = i;
567 }
568
569 f0=*f=lfun2dSdN(x,n);
570
571 copyArray(x,x0,n);
572 SIZEp=999;
573
574 gradientB (n, x0, f0, g0, tv, xmark);
575
576 initIdentityMatrix(H,nfree);
577
578 for(Iround=0; Iround<maxround; Iround++) {
579 for(i=0,initArray(p,n); i<nfree; i++) {
580 for(j=0; j<nfree; j++) p[ix[i]] -= H[i*nfree+j]*g0[ix[j]];
581 }
582
583 sizep0 = SIZEp;
584 SIZEp=norm(p,n); /* check this */
585
586 for (i=0,am=maxstep; i<n; i++) { /* max step length */
587 if (p[i]>0 && (xb[i][1]-x0[i])/p[i]<am)
588 am=(xb[i][1]-x0[i])/p[i];
589 else if (p[i]<0 && (xb[i][0]-x0[i])/p[i]<am)
590 am=(xb[i][0]-x0[i])/p[i];
591 }
592
593 if (Iround==0) {
594 h=fabs(2*f0*.01/innerp(g0,p,n)); /* check this?? */
595 h=min2(h,am/2000);
596 }
597 else {
598 h=norm(s,nfree)/SIZEp;
599 h=max2(h,am/500);
600 }
601 h=max2(h,1e-5); h=min2(h,am/5);
602 *f=f0;
603 alpha = LineSearch2(f,x0,p,h,am, min2(1e-3,e), tv,n); /* n or nfree? */
604
605 fail=0;
606 for(i=0; i<n; i++) x[i]=x0[i]+alpha*p[i];
607
608 w=min2(2,e*1000);
609 if(e<1e-4 && e>1e-6)
610 w=0.01;
611
612 if(Iround==0 || SIZEp<sizep0 || (SIZEp<.001 && sizep0<.001))
613 goodtimes++;
614 else
615 goodtimes=0;
616 if((n==1||goodtimes>=Ngoodtimes) && SIZEp<(e>1e-5?.05:.001) && (f0- *f<0.001) && H_end(x0,x,f0,*f,e,e,n))
617 break;
618
619 gradientB (n, x, *f, g, tv, xmark);
620
621 /* modify the working set */
622 for (i=0; i<n; i++) { /* add constraints, reduce H */
623
624 if (xmark[i])
625 continue;
626
627 if (fabs(x[i]-xb[i][0])<1e-6 && -g[i]<0)
628 xmark[i]=-1;
629 else if (fabs(x[i]-xb[i][1])<1e-6 && -g[i]>0)
630 xmark[i]=1;
631
632 if (xmark[i]==0)
633 continue;
634
635 copyArray (H, C, nfree*nfree);
636 for (it=0; it<nfree; it++)
637 if (ix[it]==i)
638 break;
639 for (i1=it; i1<nfree-1; i1++)
640 ix[i1]=ix[i1+1];
641 for (i1=0,nfree--; i1<nfree; i1++)
642 for (i2=0; i2<nfree; i2++)
643 H[i1*nfree+i2]=C[(i1+(i1>=it))*(nfree+1) + i2+(i2>=it)];
644 }
645
646 for (i=0,it=0,w=0; i<n; i++) { /* delete a constraint, enlarge H */
647 if (xmark[i]==-1 && -g[i]>w) {
648 it=i; w=-g[i];
649 }
650 else if (xmark[i]==1 && -g[i]<-w) {
651 it=i; w=g[i];
652 }
653 }
654
655 if (w>10*SIZEp/nfree) {
656 copyArray (H, C, nfree*nfree);
657
658 for (i1=0; i1<nfree; i1++) {
659 for (i2=0; i2<nfree; i2++) {
660 H[i1*(nfree+1)+i2]=C[i1*nfree+i2];
661 }
662 }
663
664 for(i1=0; i1<nfree+1; i1++) {
665 H[i1*(nfree+1)+nfree]=H[nfree*(nfree+1)+i1]=0;
666 }
667
668 H[(nfree+1)*(nfree+1)-1]=1;
669 xmark[it]=0;
670 ix[nfree++]=it;
671 }
672
673 for(i=0,f0=*f; i<nfree; i++) {
674 y[i]=g[ix[i]]-g0[ix[i]];
675 s[i]=x[ix[i]]-x0[ix[i]];
676 }
677 for (i=0; i<n; i++) {
678 g0[i]=g[i];
679 x0[i]=x[i];
680 }
681
682 //reduce loop
683 for (i=0,w=v=0.; i<nfree; i++) {
684 for (j=0,z[i]=0.; j<nfree; j++) {
685 z[i]+=H[i*nfree+j]*y[j];
686 }
687 w+=y[i]*z[i];
688 v+=y[i]*s[i];
689 }
690 if (fabs(v)<small) {
691 initIdentityMatrix(H,nfree);
692 fail=1;
693 continue;
694 }
695 for(i=0; i<nfree; i++) {
696 for(j=0; j<nfree; j++) {
697 H[i*nfree+j] += ((1+w/v)*s[i]*s[j]-z[i]*s[j]-s[i]*z[j])/v;
698 }
699 }
700
701 }//end of for(Iround...)
702
703 *f=lfun2dSdN(x,n);
704
705 if(nfree==n) {
706 copyArray(H, space, n*n);
707 return 1;
708 }
709 return 0;
710 }
711
712
713 void GY94::HouseholderRealSym(double a[], int n, double d[], double e[]) {
714
715 int m,k,j,i;
716 double scale,hh,h,g,f;
717
718 for (i=n-1;i>=1;i--) {
719 m=i-1;
720 h=scale=0;
721 if (m > 0) {
722 for (k=0;k<=m;k++)
723 scale += fabs(a[i*n+k]);
724 if (scale == 0)
725 e[i]=a[i*n+m];
726 else {
727 for (k=0;k<=m;k++) {
728 a[i*n+k] /= scale;
729 h += a[i*n+k]*a[i*n+k];
730 }
731 f=a[i*n+m];
732 g=(f >= 0 ? -sqrt(h) : sqrt(h));
733 e[i]=scale*g;
734 h -= f*g;
735 a[i*n+m]=f-g;
736 f=0;
737 for (j=0;j<=m;j++) {
738 a[j*n+i]=a[i*n+j]/h;
739 g=0;
740 for (k=0;k<=j;k++)
741 g += a[j*n+k]*a[i*n+k];
742 for (k=j+1;k<=m;k++)
743 g += a[k*n+j]*a[i*n+k];
744 e[j]=g/h;
745 f += e[j]*a[i*n+j];
746 }
747 hh=f/(h*2);
748 for (j=0;j<=m;j++) {
749 f=a[i*n+j];
750 e[j]=g=e[j]-hh*f;
751 for (k=0;k<=j;k++)
752 a[j*n+k] -= (f*e[k]+g*a[i*n+k]);
753 }
754 }
755 }
756 else
757 e[i]=a[i*n+m];
758 d[i]=h;
759 }
760 d[0]=e[0]=0;
761
762 /* Get eigenvectors */
763 for (i=0;i<n;i++) {
764 m=i-1;
765 if (d[i]) {
766 for (j=0;j<=m;j++) {
767 g=0;
768 for (k=0;k<=m;k++)
769 g += a[i*n+k]*a[k*n+j];
770 for (k=0;k<=m;k++)
771 a[k*n+j] -= g*a[k*n+i];
772 }
773 }
774 d[i]=a[i*n+i];
775 a[i*n+i]=1;
776 for (j=0;j<=m;j++) a[j*n+i]=a[i*n+j]=0;
777 }
778 }
779
780 int GY94::EigenTridagQLImplicit(double d[], double e[], int n, double z[]) {
781
782 int m,j,iter,niter=30, status=0, i,k;
783 double s,r,p,g,f,dd,c,b, aa,bb;
784
785 for (i=1;i<n;i++) e[i-1]=e[i]; e[n-1]=0;
786 for (j=0;j<n;j++) {
787 iter=0;
788 do {
789 for (m=j;m<n-1;m++) {
790 dd=fabs(d[m])+fabs(d[m+1]);
791 if (fabs(e[m])+dd == dd) break; /* ??? */
792 }
793 if (m != j) {
794 if (iter++ == niter) {
795 status=-1;
796 break;
797 }
798 g=(d[j+1]-d[j])/(2*e[j]);
799
800 /* r=pythag(g,1); */
801
802 if((aa=fabs(g))>1) r=aa*sqrt(1+1/(g*g));
803 else r=sqrt(1+g*g);
804
805 g=d[m]-d[j]+e[j]/(g+SIGN(r,g));
806 s=c=1;
807 p=0;
808 for (i=m-1;i>=j;i--) {
809 f=s*e[i];
810 b=c*e[i];
811
812 /* r=pythag(f,g); */
813 aa=fabs(f); bb=fabs(g);
814 if(aa>bb) { bb/=aa; r=aa*sqrt(1+bb*bb); }
815 else if(bb==0) r=0;
816 else { aa/=bb; r=bb*sqrt(1+aa*aa); }
817
818 e[i+1]=r;
819 if (r == 0) {
820 d[i+1] -= p;
821 e[m]=0;
822 break;
823 }
824 s=f/r;
825 c=g/r;
826 g=d[i+1]-p;
827 r=(d[i]-g)*s+2*c*b;
828 d[i+1]=g+(p=s*r);
829 g=c*r-b;
830 for (k=0;k<n;k++) {
831 f=z[k*n+i+1];
832 z[k*n+i+1]=s*z[k*n+i]+c*f;
833 z[k*n+i]=c*z[k*n+i]-s*f;
834 }
835 }
836 if (r == 0 && i >= j) continue;
837 d[j]-=p; e[j]=g; e[m]=0;
838 }
839 } while (m != j);
840 }
841 return status;
842 }
843
844 void GY94::EigenSort(double d[], double U[], int n) {
845
846 int k,j,i;
847 double p;
848
849 for (i=0;i<n-1;i++) {
850 p=d[k=i];
851 for (j=i+1;j<n;j++) {
852 if (d[j] >= p) p=d[k=j];
853 }
854 if (k != i) {
855 d[k]=d[i];
856 d[i]=p;
857 for (j=0;j<n;j++) {
858 p=U[j*n+i];
859 U[j*n+i]=U[j*n+k];
860 U[j*n+k]=p;
861 }
862 }
863 }
864 }
865
866 int GY94::eigenRealSym(double A[], int n, double Root[], double work[]) {
867
868 int status=0;
869 HouseholderRealSym(A, n, Root, work);
870 status=EigenTridagQLImplicit(Root, work, n, A);
871 EigenSort(Root, A, n);
872
873 return(status);
874 }
875
876
877 int GY94::eigenQREV (double Q[], double pi[], int n, double Root[], double U[], double V[], double spacesqrtpi[]) {
878
879 int i,j, inew, jnew, nnew, status;
880 double *pi_sqrt=spacesqrtpi, small=1e-6;
881
882 for(j=0,nnew=0; j<n; j++)
883 if(pi[j]>small)
884 pi_sqrt[nnew++]=sqrt(pi[j]);
885
886 /* store in U the symmetrical matrix S = sqrt(D) * Q * sqrt(-D) */
887
888 if(nnew==n) {
889 for(i=0; i<n; i++)
890 for(j=0,U[i*n+i] = Q[i*n+i]; j<i; j++)
891 U[i*n+j] = U[j*n+i] = (Q[i*n+j] * pi_sqrt[i]/pi_sqrt[j]);
892
893 status=eigenRealSym(U, n, Root, V);
894 for(i=0;i<n;i++) for(j=0;j<n;j++) V[i*n+j] = U[j*n+i] * pi_sqrt[j];
895 for(i=0;i<n;i++) for(j=0;j<n;j++) U[i*n+j] /= pi_sqrt[i];
896 }
897 else {
898 for(i=0,inew=0; i<n; i++) {
899 if(pi[i]>small) {
900 for(j=0,jnew=0; j<i; j++)
901 if(pi[j]>small) {
902 U[inew*nnew+jnew] = U[jnew*nnew+inew]
903 = Q[i*n+j] * pi_sqrt[inew]/pi_sqrt[jnew];
904 jnew++;
905 }
906 U[inew*nnew+inew] = Q[i*n+i];
907 inew++;
908 }
909 }
910
911 status=eigenRealSym(U, nnew, Root, V);
912
913 for(i=n-1,inew=nnew-1; i>=0; i--) /* construct Root */
914 Root[i] = (pi[i]>small ? Root[inew--] : 0);
915 for(i=n-1,inew=nnew-1; i>=0; i--) { /* construct V */
916 if(pi[i]>small) {
917 for(j=n-1,jnew=nnew-1; j>=0; j--)
918 if(pi[j]>small) {
919 V[i*n+j] = U[jnew*nnew+inew]*pi_sqrt[jnew];
920 jnew--;
921 }
922 else
923 V[i*n+j] = (i==j);
924 inew--;
925 }
926 else
927 for(j=0; j<n; j++) V[i*n+j] = (i==j);
928 }
929 for(i=n-1,inew=nnew-1; i>=0; i--) { /* construct U */
930 if(pi[i]>small) {
931 for(j=n-1,jnew=nnew-1;j>=0;j--)
932 if(pi[j]>small) {
933 U[i*n+j] = U[inew*nnew+jnew]/pi_sqrt[inew];
934 jnew--;
935 }
936 else
937 U[i*n+j] = (i==j);
938 inew--;
939 }
940 else
941 for(j=0;j<n;j++)
942 U[i*n+j] = (i==j);
943 }
944 }
945
946 return status;
947 }
948
949 /********************************************
950 * Function: parseSubRates
951 * Input Parameter: string, array of double
952 * Output: Parse substitution rates according to the given model
953 * Return Value: int
954 *********************************************/
955 int GY94::parseSubRates(string model, double kappa[]) {
956
957 /*
958 JC, F81: rTC==rAG =rTA==rCG==rTG==rCA
959 K2P, HKY: rTC==rAG!=rTA==rCG==rTG==rCA
960 TNEF,TN: rTC!=rAG!=rTA==rCG==rTG==rCA
961 K3P, K3PUF: rTC==rAG!=rTA==rCG!=rTG==rCA
962 TIMEF, TIM: rTC!=rAG!=rTA==rCG!=rTG==rCA
963 TVMEF, TVM: rTC==rAG!=rTA!=rCG!=rTG!=rCA
964 SYM, GTR: rTC!=rAG!=rTA!=rCG!=rTG!=rCA
965 */
966
967 kappa[5]=1.0; //Substitution rate between C and A
968
969 if(model=="JC" || model=="F81") {
970 //Q[i*n+j] = kappa[0] = 1;
971 kappa[0]=kappa[1]=kappa[2]=kappa[3]=kappa[4]=kappa[5];
972 }
973 else if(model=="K2P" || model=="HKY") {//one para.
974 //if ((b1+b2)==1 || (b1+b2)==5) Q[i*n+j] = kappa[0];
975 kappa[1]=kappa[0];
976 kappa[2]=kappa[3]=kappa[4]=kappa[5];
977 }
978 else if(model=="TNEF" || model=="TN") {//two para.
979 //if ((b1+b2)==1) Q[i*n+j] = kappa[0];
980 //else if ((b1+b2)==5) Q[i*n+j] = kappa[1];
981 kappa[2]=kappa[3]=kappa[4]=kappa[5];
982 }
983 else if(model=="K3P" || model=="K3PUF") {//two para.
984 //if ((b1+b2)==1 || (b1+b2)==5) Q[i*n+j] = kappa[0];
985 //else if ((b1+b2)==2 || (b1+b2)==4) Q[i*n+j] = kappa[1];
986 kappa[4]=kappa[5];
987 kappa[2]=kappa[3]=kappa[1];
988 kappa[1] = kappa[0];
989 }
990 else if(model=="TIMEF" || model=="TIM") {//three para.
991 //if ((b1+b2)==1) Q[i*n+j] = kappa[0];
992 //else if ((b1+b2)==5) Q[i*n+j] = kappa[1];
993 //else if ((b1+b2)==2 || (b1+b2)==4) Q[i*n+j] = kappa[2];
994 kappa[4]=kappa[5];
995 kappa[3]=kappa[2];
996 }
997 else if(model=="TVMEF" || model=="TVM") {//four para.
998 //if ((b1+b2)==1 || (b1+b2)==5) Q[i*n+j] = kappa[0];
999 //else if ((b1+b2)==2) Q[i*n+j] = kappa[1];
1000 //else if ((b1+b2)==4) Q[i*n+j] = kappa[2];
1001 //else if ((b1==0) && (b2==3)) Q[i*n+j] = kappa[3];
1002
1003 //TVMEF, TVM: rTC==rAG!=rTA!=rCG!=rTG!=rCA
1004 //SYM, GTR: rTC!=rAG!=rTA!=rCG!=rTG!=rCA
1005 kappa[4]=kappa[3];
1006 kappa[3]=kappa[2];
1007 kappa[2]=kappa[1];
1008 kappa[1]=kappa[0];
1009 }
1010 else if(model=="SYM" || model=="GTR") {//five para.
1011 }
1012 else {
1013 return 0;
1014 }
1015
1016 return 1;
1017 }
1018
1019 /* Construct Q: transition probability matrix 64*64 */
1020 int GY94::EigenQc (int getstats, double blength, double *S, double *dS, double *dN, double Root[], double U[], double V[], double kappa[], double omega, double Q[]) {
1021 /* This contructs the rate matrix Q for codon substitution and get the eigen
1022 values and vectors if getstats==0, or get statistics (dS & dN) if
1023 getstats==1.
1024 The routine is also called by Qcodon2aa for mechanistic amino acid
1025 substitution models.
1026 Input parameters are kappa, omega and com.pi (or com.fb61).
1027
1028 Statistics calculated include S, dS & dN.
1029 c0[0,1,2] and c[0,1,2] are rates for the 3 codon positions before and after
1030 selection. c4 is for 4-fold rates. ts[3] and tv[3] are transition/
1031 transversion rates, not calculated.
1032
1033 Under NSsites or other site-class models, this function does not scale
1034 Q but calculates the Qfactor_NS.
1035 DetailOutput() uses this function to calculate dS & dN; for each omega
1036 component, dS and dN are not scaled here but are scaled in DetailOutput()
1037 after Qfactor_NS is calculated.
1038
1039 aaDist=FIT1 & FIT2: ap,p*,av,v*, (and w0 for FIT2)
1040 The argument omega is used only if the model assumes one omega. For
1041 AAClasses, com.pomega is used instead.
1042 */
1043 int n=Nsensecodon, i,j,k, ic1,ic2,aa1,aa2;
1044 int ndiff,pos=0,from[3],to[3];
1045 double mr, rs0,ra0,rs,ra; /* rho's */
1046 double d4=0, d0[3],d[3],ts[3],tv[3]; /* rates at positions and 4-fold sites */
1047 double *pi=com.pi, w=-1, pijQij;
1048 double space_pisqrt[CODON];
1049
1050 for(i=0; i<3; d[i]=d0[i]=ts[i]=tv[i]=0, i++);
1051 initArray(Q, n*n);
1052
1053 //Equal codon frequency
1054 if(model=="JC" || model=="K2P" || model=="TNEF" || model=="K3P"||
1055 model=="TIMEF" ||model=="TVMEF"|| model=="SYM") {
1056 initArray(com.pi, 64, 1.0/(64-getNumNonsense(genetic_code)));
1057 }
1058
1059 //Parse substitution rates according to the given model
1060 parseSubRates(model, kappa);
1061
1062 //Construct Q: transition probability matrix 64*64
1063 for (i=0, rs0=ra0=rs=ra=0; i<n; i++) {
1064
1065 //codon i
1066 ic1=FROM61[i]; from[0]=ic1/16; from[1]=(ic1/4)%4; from[2]=ic1%4;
1067
1068 //codon j
1069 for(j=0; j<i; j++) {
1070 ic2=FROM61[j]; to[0]=ic2/16; to[1]=(ic2/4)%4; to[2]=ic2%4;
1071 for(k=0,ndiff=0; k<3; k++) {
1072 if(from[k]!=to[k]) {
1073 ndiff++; pos=k;
1074 }
1075 }
1076
1077 //consider only one difference between two codons
1078 if(ndiff!=1) continue;
1079
1080 int b1 = min2(from[pos],to[pos]);
1081 int b2 = max2(from[pos],to[pos]);
1082
1083 /* 01 23 02 13 03 12
1084 JC, F81: rTC==rAG =rTA==rCG==rTG==rCA
1085 K2P, HKY: rTC==rAG!=rTA==rCG==rTG==rCA
1086 TNEF,TN: rTC!=rAG!=rTA==rCG==rTG==rCA
1087 K3P, K3PUF: rTC==rAG!=rTA==rCG!=rTG==rCA
1088 TIMEF, TIM: rTC!=rAG!=rTA==rCG!=rTG==rCA
1089 TVMEF, TVM: rTC==rAG!=rTA!=rCG!=rTG!=rCA
1090 SYM, GTR: rTC!=rAG!=rTA!=rCG!=rTG!=rCA
1091 */
1092 if (b1==0 && b2==1) Q[i*n+j]=kappa[0]; /* TC */
1093 else if(b1==2 && b2==3) Q[i*n+j]=kappa[1]; /* AG */
1094 else if(b1==0 && b2==2) Q[i*n+j]=kappa[2]; /* TA */
1095 else if(b1==1 && b2==3) Q[i*n+j]=kappa[3]; /* CG */
1096 else if(b1==0 && b2==3) Q[i*n+j]=kappa[4]; /* TG */
1097 else if(b1==1 && b2==2) Q[i*n+j]=kappa[5]; /* CA=1 */
1098
1099 Q[j*n+i]=Q[i*n+j];
1100 Q[i*n+j]*=com.pi[j];
1101 Q[j*n+i]*=com.pi[i];
1102
1103 //probability
1104 pijQij=2*pi[i]*Q[i*n+j];
1105
1106 aa1=GeneticCode[com.icode][ic1];
1107 aa2=GeneticCode[com.icode][ic2];
1108 if(aa1==aa2) {//synonymous
1109 rs+=pijQij;
1110 }
1111 else {//nonsynonymous
1112 ra0+=pijQij;
1113 w = omega;
1114
1115 Q[i*n+j]*=w; Q[j*n+i]*=w;
1116 ra+=pijQij*w;
1117 }
1118
1119 } /* for (j) */
1120 } /* for (i) */
1121
1122 mr=rs+ra;
1123 if(getstats) {
1124 rs0=rs;
1125 w=(rs0+ra0); rs0/=w; ra0/=w; *S=rs0*3*com.ls;
1126 if(blength>=0) { /* calculates dS & dN */
1127 if(blength==0) *dS = *dN = 0;
1128 rs/=mr;
1129 ra/=mr;
1130 *dS=blength*rs/(3*rs0);
1131 *dN=blength*ra/(3*ra0);
1132 }
1133
1134 }
1135 else {
1136 for (i=0; i<n; i++) Q[i*n+i]=-sumArray(Q+i*n,n);
1137
1138 eigenQREV(Q, com.pi, n, Root, U, V, space_pisqrt);
1139
1140 for(i=0; i<n; i++) Root[i]/=mr;
1141
1142 }
1143 return 0;
1144 }
1145
1146 /* Return maximum-likelihood score */
1147 double GY94::lfun2dSdN(double x[], int np) {
1148 /* likelihood function for calculating dS and dN between 2 sequences,
1149 com.z[0] & com.z[1]:
1150 prob(i,j) = PI_i * p(i,j,t)
1151
1152 Data are clean and coded.
1153 Transition probability pijt is calculated for observed patterns only.
1154 */
1155 int n=com.ncode, h,k, ik, z0,z1;
1156 double fh,expt[CODON], lnL1=0;
1157 double *pkappa=com.KAPPA;
1158
1159 //cout<<name.c_str()<<": "<<com.ncode<<"\t"<<Nsensecodon<<endl;
1160
1161 k=1, ik=0;
1162 com.kappa=x[k];
1163 for (ik=0; ik<com.nkappa; ik++)
1164 pkappa[ik]=x[k++];
1165
1166 com.omega=x[1+com.nkappa];
1167
1168 EigenQc(0,-1,NULL,NULL,NULL, Root, U, V, pkappa, com.omega, PMat);
1169
1170 //t = x[0], exp(Qt)
1171 for(k=0; k<n; k++) {
1172 expt[k] = exp(x[0]*Root[k]);
1173 }
1174
1175 //com.npatt = number of patterns
1176 for (h=0; h<com.npatt; h++) {
1177
1178 if(com.fpatt[h]<Small_Diff) {
1179 continue;
1180 }
1181
1182 z0=com.z[0][h];
1183 z1=com.z[1][h];
1184
1185 for(k=0,fh=0;k<n;k++) {
1186 fh += U[z0*n+k]*expt[k]*V[k*n+z1];
1187 }
1188
1189 fh*=com.pi[z0];
1190
1191 lnL1-=log(fh)*com.fpatt[h];
1192 }
1193
1194 return lnL1;
1195 }
1196
1197 /* Main function to calculate Ka and Ks, called by "Run" */
1198 int GY94::PairwiseCodon(double space[]) {
1199
1200 char *pz0[2];
1201 int npatt0=com.npatt;
1202 double *fpatt0, ls0=com.ls;
1203 double fp[CODON*CODON];
1204 int n=com.ncode, j,k,h;
1205 double x[10]={.3,1,.5,.5,.5,.5,.3}, xb[10][2]={{1e-6,3.}};
1206 double kappab[2]={.01,30}, omegab[2]={0.001, 50.};
1207 double e=1e-6, dS,dN;
1208 double *pkappa=com.KAPPA;
1209
1210 //fpatt0=(double*)malloc(npatt0*3*sizeof(double));
1211 fpatt0= new double[npatt0*3];
1212 for(k=0; k<com.ns; k++) pz0[k]=com.z[k];
1213 com.z[0]=(char*)(fpatt0+npatt0);
1214 com.z[1]=com.z[0]+npatt0;
1215
1216 copyArray(com.fpatt, fpatt0, npatt0);
1217
1218 //t > snp/length
1219 if ((snp/length) > 1e-6) xb[0][0] = 3*snp/length;
1220
1221 //com.nkappa = 2;
1222 for(j=0; j<com.nkappa; j++) {
1223 xb[1+j][0]=kappab[0];
1224 xb[1+j][1]=kappab[1];
1225 }
1226
1227 k = 1+com.nkappa;
1228 xb[k][0] = omegab[0];
1229 xb[k][1] = omegab[1];
1230
1231 //fp: codon -> codon frequency
1232 initArray(fp, CODON*CODON);
1233 for(h=0; h<npatt0; h++) {
1234 j = max2(pz0[0][h], pz0[1][h]);
1235 k = min2(pz0[0][h], pz0[1][h]);
1236 fp[j*n+k] += fpatt0[h];
1237 }
1238
1239 //fp -> com.fpatt
1240 for(j=0,com.npatt=0;j<n;j++) {
1241 for(k=0;k<j+1;k++)
1242 if(fp[j*n+k]) {
1243 com.z[0][com.npatt]=(char)j;
1244 com.z[1][com.npatt]=(char)k;
1245 com.fpatt[com.npatt++]=fp[j*n+k];
1246 }
1247 }
1248
1249 //com.pi: sense codon's frequencies
1250 for(j=0,initArray(com.pi,n); j<com.npatt; j++) {
1251 com.pi[(int)com.z[0][j]]+=com.fpatt[j]/(2.*com.ls);
1252 com.pi[(int)com.z[1][j]]+=com.fpatt[j]/(2.*com.ls);
1253 }
1254
1255 GetCodonFreqs(com.pi);
1256
1257 /* initial values and bounds */
1258 //divergence time t
1259 x[0]=-1;
1260 if(x[0]>3)
1261 x[0]=1.5+rndu();
1262 if(x[0]<1e-6)
1263 x[0]=.5*rndu();
1264
1265 //kappas
1266 for(j=0; j<com.nkappa; j++)
1267 x[1+j]=.2+.4*rndu();
1268
1269 //Ka/Ks
1270 k=1+com.nkappa;
1271 x[k]=(3*x[k]+0.6*rndu())/4;
1272 x[k]=max2(x[k],0.01);
1273 x[k]=min2(x[k],2);
1274
1275 if ((snp/length) > 1e-6) x[0] = 3.0*snp/length;
1276
1277 ming2(&lnL, x, xb, space, e, com.np);
1278
1279 EigenQc(1, x[0], &S, &dS, &dN, NULL, NULL, NULL, pkappa, com.omega, PMat);
1280
1281 Ka = dN;
1282 Ks = dS;
1283 N = com.ls*3-S;
1284 Sd = Ks*S;
1285 Nd = Ka*N;
1286 double scale=(Sd+Nd)/snp;
1287 Sd /= scale;
1288 Nd /= scale;
1289
1290 lnL = -lnL;
1291 //AICc = -2log(lnL) + 2K + 2K(K+1)/(n-K-1), K=parameters' number, n=sample size
1292 AICc = -2*lnL + 2.*(com.nkappa+2)*(length/3)/((length/3)-(com.nkappa+2)-1.);
1293
1294 t = x[0]/3;
1295 //kappa = com.kappa;
1296
1297 parseSubRates(model, pkappa);
1298 copyArray(pkappa, KAPPA, NUMBER_OF_RATES);
1299
1300 k=com.np-1;
1301
1302 com.ls=(int)ls0;
1303 for(k=0; k<com.ns; k++) com.z[k]=pz0[k];
1304 com.npatt=npatt0;
1305
1306 for(h=0; h<npatt0; h++) com.fpatt[h]=fpatt0[h];
1307 //free(fpatt0);
1308 delete []fpatt0;
1309
1310 return 0;
1311 }
1312
1313 void GY94::FreeMemPUVR(void){
1314 //free(PMat);
1315 }
1316
1317 /* Main fuction for GY method */
1318 string GY94::Run(const char *seq1, const char *seq2) {
1319
1320 int i;
1321
1322 str1 = seq1;
1323 str2 = seq2;
1324
1325 preProcess(seq1, seq2);
1326
1327 com.sspace = max2(800000, 3*com.ncode*com.ncode*(int)sizeof(double));
1328 //com.space = (double*)realloc(com.space, com.sspace);
1329 com.space = new double[com.sspace];
1330
1331 PairwiseCodon(com.space);
1332
1333 FreeMemPUVR();
1334
1335 for(i=0; i<com.ns; i++) delete []com.z[i];
1336
1337 //free(com.fpatt);
1338 //free(com.space);
1339 delete []com.space;
1340
1341 return parseOutput();
1342 }
1343
A place to share codes with mates