source: branches/profile/GDE/PHYML/bionj.c

Last change on this file was 4073, checked in by westram, 13 years ago
  • phyml 2.4.5
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.9 KB
Line 
1/*
2
3PHYML :  a program that  computes maximum likelihood  phylogenies from
4DNA or AA homologous sequences
5
6Copyright (C) Stephane Guindon. Oct 2003 onward
7
8All parts of  the source except where indicated  are distributed under
9the GNU public licence.  See http://www.opensource.org for details.
10
11*/
12
13/*
14
15The code below is an implementation of the building tree algorithm
16described in "BIONJ: an improved version of the NJ algorithm based
17on a simple model of sequence data." (1997) O. Gascuel. Mol Biol Evol.
1814:685-95. 
19
20*/
21
22#include "bionj.h"
23
24void Bionj(matrix *mat)
25{
26  int x,y,i;
27  double vxy,lx,ly,lamda,score;
28
29  For(i,mat->tree->n_otu)
30    mat->tip_node[i] = mat->tree->noeud[i];
31
32
33  while(mat->r > 3)
34    {
35      x = y = 0;
36      vxy = .0;
37      score = .0;
38      Compute_Sx(mat);
39      Best_Pair(mat,&x,&y,&score);
40      vxy=Variance(mat,x,y);
41      lx=Br_Length(mat,x,y);   
42      ly=Br_Length(mat,y,x);
43      lamda=Lamda(mat,x,y,vxy); 
44      Update_Mat(mat,x,y,lx,ly,vxy,lamda);
45      Update_Tree(mat,x,y,lx,ly,score);     
46    }
47
48  Finish(mat);
49  i=0;
50
51  Init_Tree_Edges(mat->tree->noeud[0],
52                  mat->tree->noeud[0]->v[0],
53                  mat->tree,&i);
54}
55 
56/*********************************************************/
57
58void Bionj_Scores(matrix *mat)
59{
60  int i;
61
62  For(i,2*mat->n_otu-3)
63    {
64      if(!mat->tree->noeud[i]->tax)
65        {
66          mat->tree->noeud[i]->b[0]->nj_score = 
67            mat->tree->noeud[i]->score[0];
68        }
69    }
70}
71
72/*********************************************************/
73
74void Finish(matrix *mat)
75{
76  double dxy,dxz,dyz;
77  int x,y,z;
78  node *nx,*ny,*nz,*new;
79  int i;
80
81  dxy = dxz = dyz = -1.;
82  x = y = z = -1;
83
84  For(i,mat->n_otu)
85    {
86      if(mat->on_off[i])
87        {
88          if(x < 0) x=i;
89          else if(y < 0) y = i;
90          else if(z < 0) z = i;
91        }
92    }
93
94  dxy = Dist(mat,x,y);
95  dxz = Dist(mat,x,z);
96  dyz = Dist(mat,y,z);
97
98  nx = mat->tip_node[x];
99  ny = mat->tip_node[y];
100  nz = mat->tip_node[z];
101
102  new = mat->tree->noeud[mat->curr_int];
103  new->num = mat->curr_int;
104  new->v[0] = nx;
105  new->v[1] = ny;
106  new->v[2] = nz;
107
108
109  nx->v[0] = new;
110  ny->v[0] = new;
111  nz->v[0] = new;
112 
113  Make_Edge_Light(new,nx);
114  Make_Edge_Light(new,ny);
115  Make_Edge_Light(new,nz);
116 
117  nx->b[0]->l = .5*(dxy-dyz+dxz);
118  ny->b[0]->l = .5*(dyz-dxz+dxy);
119  nz->b[0]->l = .5*(dxz-dxy+dyz);
120   
121  new->b[0]->l = nx->b[0]->l;
122  new->b[1]->l = ny->b[0]->l;
123  new->b[2]->l = nz->b[0]->l;
124}
125
126/*********************************************************/
127
128void Update_Mat(matrix *mat, int x, int y, double lx, double ly, double vxy, double lamda)
129{
130  int i;
131  int a,b;
132 
133  a = b = -1;
134  For(i,mat->n_otu)
135    {
136      if((mat->on_off[i]) && (i != x) && (i != y))
137        {
138          if(x > i)
139            {
140              a=x;
141              b=i;
142            }
143          else
144            {
145              a=i;
146              b=x;
147            }
148          mat->dist[a][b]=Dist_Red(mat,x,lx,y,ly,i,lamda);
149          mat->dist[b][a]=Var_Red(mat,x,y,i,lamda,vxy);
150        }
151    }
152}
153
154/*********************************************************/
155
156void Update_Tree(matrix *mat, int x, int y, double lx, double ly, double score)
157{
158  node *new, *nx, *ny;
159
160  nx = mat->tip_node[x];
161  ny = mat->tip_node[y];
162  new = mat->tree->noeud[mat->curr_int];
163  nx->v[0] = new;
164  ny->v[0] = new;
165  new->v[1] = nx;
166  new->v[2] = ny;
167 
168  new->num = mat->curr_int;
169
170  Make_Edge_Light(new,nx);
171  Make_Edge_Light(new,ny);
172
173  nx->b[0]->l = lx;
174  ny->b[0]->l = ly;
175 
176  new->b[1]->l = lx;
177  new->b[2]->l = ly;
178  new->score[0] = score;
179
180  nx->l[0] = lx;
181  ny->l[0] = ly;
182 
183  new->l[1] = lx;
184  new->l[2] = ly;
185 
186  mat->tip_node[x] = new;
187  mat->on_off[y] = 0;
188  mat->curr_int++;
189  mat->r--;
190}
191
192/*********************************************************/
193
194void Best_Pair(matrix *mat, int *x, int *y,double *score)
195{
196  int i,j;
197  double Qij,Qmin,Qmin2;
198  double **t_Qij;
199
200  t_Qij = (double **)mCalloc(mat->n_otu,sizeof(double *));
201  For(i,mat->n_otu)
202    t_Qij[i] = (double *)mCalloc(mat->n_otu,sizeof(double));
203
204  Qmin = 1.e+10;
205  Qij = Qmin;
206
207  For(i,mat->n_otu)
208    {
209      if(mat->on_off[i])
210        {
211          for(j=0;j<i;j++)
212            {
213              if(mat->on_off[j])
214                {
215                  Qij = Q_Agglo(mat,i,j);
216                  t_Qij[i][j] = Qij;
217
218                  if(Qij < Qmin)
219                    {
220                      *x = i;
221                      *y = j;
222                      Qmin = Qij;
223                    }
224                }
225            }
226        }
227    }
228
229  Qmin2 = 1e+10;
230
231  For(i,mat->n_otu)
232    {
233      if((i != *y) && (i != *x) && (t_Qij[*x][i] < Qmin2)) Qmin2 = t_Qij[*x][i];
234    }
235
236  For(i,mat->n_otu)
237    {
238      if((i != *y) && (i != *x) && (t_Qij[i][*y] < Qmin2)) Qmin2 = t_Qij[i][*y];
239    }
240
241  *score = fabs(Qmin2 - Qmin)/fabs(Qmin);
242 
243  For(i,mat->n_otu) free(t_Qij[i]);
244  free(t_Qij);
245}
246/*********************************************************/
247
248void Compute_Sx(matrix *mat)
249{
250  int i,j;
251 
252  For(i,mat->n_otu)
253    {
254      mat->dist[i][i] = .0;
255      if(mat->on_off[i])
256        {
257          For(j,mat->n_otu)
258            {
259              if((i != j) && (mat->on_off[j]))
260                {
261                  mat->dist[i][i] += Dist(mat,i,j);
262                }
263            }
264        }
265    }
266}
267             
268/*********************************************************/
269
270double Sum_S(matrix *mat, int i)
271{
272  return mat->dist[i][i];
273}
274
275/*********************************************************/
276
277double Dist(matrix *mat, int x, int y)
278{
279    if(x > y)
280      return(mat->dist[x][y]);
281    else
282      return(mat->dist[y][x]);
283}
284
285/*********************************************************/
286
287double Variance(matrix *mat, int x, int y)
288{
289    if(x > y)
290      {
291        return(mat->dist[y][x]);
292      }
293    else
294      {
295        return(mat->dist[x][y]);
296      }
297}
298
299/*********************************************************/
300
301double Br_Length(matrix *mat, int x, int y)
302{
303    return .5*(Dist(mat,x,y)+
304                (Sum_S(mat,x)-Sum_S(mat,y))/(double)(mat->r-2.)); 
305}
306
307/*********************************************************/
308
309double Dist_Red(matrix *mat, int x, double lx, int y, double ly, int i, double lamda)
310{
311  double Dui;
312  Dui=lamda*(Dist(mat,x,i)-lx)
313     +(1.-lamda)*(Dist(mat,y,i)-ly);
314  return(Dui);
315}
316
317/*********************************************************/
318
319double Var_Red(matrix *mat, int x, int y, int i, double lamda, double vxy)
320{
321  double Vui;
322  Vui=lamda*(Variance(mat,x,i))
323     +(1.-lamda)*(Variance(mat,y,i))
324    -lamda*(1.-lamda)*vxy;
325  return(Vui);
326}
327
328/*********************************************************/
329
330double Lamda(matrix *mat, int x, int y, double vxy)
331{
332    double lamda=0.0;
333    int i;
334   
335    if(mat->method == 0) /* NJ (Saitou & Nei, 1987) */
336      lamda = 0.5;
337    else /* BioNJ (Gascuel, 1997) */
338      {
339        if(vxy==0.0)
340          lamda=0.5;
341        else
342          {
343            For(i,mat->n_otu)
344              {
345                if((x != i) && (y != i) && (mat->on_off[i]))
346                  lamda = lamda + Variance(mat,y,i) - Variance(mat,x,i);
347              }
348            lamda = 0.5 + lamda/(2.*(mat->r-2)*vxy);
349          }
350       
351        if(lamda > 1.0)
352          lamda = 0.5;/*1.0;*/
353        else if(lamda < 0.0)
354          lamda = 0.5;/*0.0;*/
355      }
356
357    return(lamda);
358}
359
360/*********************************************************/
361
362double Q_Agglo(matrix *mat, int x, int y)
363{
364  double Qxy;
365
366  Qxy = .0;
367
368  Qxy=(mat->r-2.)*Dist(mat,x,y)
369    -Sum_S(mat,x)
370    -Sum_S(mat,y); 
371  return(Qxy);                       
372}
373
374/*********************************************************/
375
376void Bionj_Br_Length(matrix *mat)
377{
378  int x;
379
380  x = Bionj_Br_Length_Post(mat->tree->noeud[0],
381                           mat->tree->noeud[0]->v[0],
382                           mat);
383  mat->tree->noeud[0]->b[0]->l = Dist(mat,0,x);
384}
385
386/*********************************************************/
387
388int Bionj_Br_Length_Post(node *a, node *d, matrix *mat)
389{
390  int i;
391
392  if(d->tax)
393    {
394      return d->num;
395    }
396  else
397    {
398      int d_v1, d_v2;
399      double lx, ly, vxy,lamda;
400      int x,y;
401
402      d_v1 = d_v2 = -1;
403      For(i,3)
404        if(d->v[i] != a) {(d_v1 < 0)?(d_v1 = i):(d_v2 = i);}
405     
406
407      x = Bionj_Br_Length_Post(d,d->v[d_v1],mat);
408      y = Bionj_Br_Length_Post(d,d->v[d_v2],mat);
409
410      vxy = .0;
411      Compute_Sx(mat);
412      vxy=Variance(mat,(x),(y));
413      lx=Br_Length(mat,(x),(y));   
414      ly=Br_Length(mat,(y),(x));
415      lamda=Lamda(mat,(x),(y),vxy); 
416      Update_Mat(mat,(x),(y),lx,ly,vxy,lamda);
417
418      d->b[d_v1]->l = lx;
419      d->b[d_v2]->l = ly;
420     
421      mat->on_off[y] = 0;
422      mat->r--;
423
424      return x;
425    }
426}
427
428/*********************************************************/
429
430
431
432
433
434
Note: See TracBrowser for help on using the repository browser.