/* PHYML : a program that computes maximum likelihood phylogenies from DNA or AA homologous sequences Copyright (C) Stephane Guindon. Oct 2003 onward All parts of the source except where indicated are distributed under the GNU public licence. See http://www.opensource.org for details. */ #include "utilities.h" #include "lk.h" #include "optimiz.h" #include "models.h" #include "free.h" #include "bionj.h" #include "simu.h" /*********************************************************/ /* NUMERICAL RECIPES ROUTINES FOR COMPUTING C(n,k) */ /*********************************************************/ double bico(int n, int k) { return floor(0.5+exp(factln(n)-factln(k)-factln(n-k))); } double factln(int n) { static double a[101]; if (n < 0){ Exit("Err: negative factorial in routine FACTLN"); } if (n <= 1) return 0.0; if (n <= 100) return a[n] ? a[n] : (a[n]=gammln(n+1.0)); else return gammln(n+1.0); } double gammln(double xx) { double x,tmp,ser; static double cof[6]={76.18009173,-86.50532033,24.01409822, -1.231739516,0.120858003e-2,-0.536382e-5}; int j; x=xx-1.0; tmp=x+5.5; tmp -= (x+0.5)*log(tmp); ser=1.0; for (j=0;j<=5;j++) { x += 1.0; ser += cof[j]/x; } return -tmp+log(2.50662827465*ser); } /*********************************************************/ /* END OF NUMERICAL RECIPES ROUTINES */ /*********************************************************/ double Pbinom(int N, int ni, double p) { return bico(N,ni)*pow(p,ni)*pow(1-p,N-ni); } /*********************************************************/ void Plim_Binom(double pH0, int N, double *pinf, double *psup) { *pinf = pH0 - 1.64*sqrt(pH0*(1-pH0)/(double)N); if(*pinf < 0) *pinf = .0; *psup = pH0 + 1.64*sqrt(pH0*(1-pH0)/(double)N); } /*********************************************************/ double LnGamma (double alpha) { /* returns ln(gamma(alpha)) for alpha>0, accurate to 10 decimal places. Stirling's formula is used for the central polynomial part of the procedure. Pike MC & Hill ID (1966) Algorithm 291: Logarithm of the gamma function. Communications of the Association for Computing Machinery, 9:684 */ double x=alpha, f=0, z; if (x<7) { f=1; z=x-1; while (++z<7) f*=z; x=z; f=-log(f); } z = 1/(x*x); return f + (x-0.5)*log(x) - x + .918938533204673 + (((-.000595238095238*z+.000793650793651)*z-.002777777777778)*z +.083333333333333)/x; } /*********************************************************/ double IncompleteGamma (double x, double alpha, double ln_gamma_alpha) { /* returns the incomplete gamma ratio I(x,alpha) where x is the upper limit of the integration and alpha is the shape parameter. returns (-1) if in error ln_gamma_alpha = ln(Gamma(alpha)), is almost redundant. (1) series expansion if (alpha>x || x<=1) (2) continued fraction otherwise RATNEST FORTRAN by Bhattacharjee GP (1970) The incomplete gamma integral. Applied Statistics, 19: 285-287 (AS32) */ int i; double p=alpha, g=ln_gamma_alpha; double accurate=1e-8, overflow=1e30; double factor, gin=0, rn=0, a=0,b=0,an=0,dif=0, term=0, pn[6]; if (x==0) return (0); if (x<0 || p<=0) return (-1); factor=exp(p*log(x)-x-g); if (x>1 && x>=p) goto l30; /* (1) series expansion */ gin=1; term=1; rn=p; l20: rn++; term*=x/rn; gin+=term; if (term > accurate) goto l20; gin*=factor/p; goto l50; l30: /* (2) continued fraction */ a=1-p; b=a+x+1; term=0; pn[0]=1; pn[1]=x; pn[2]=x+1; pn[3]=x*b; gin=pn[2]/pn[3]; l32: a++; b+=2; term++; an=a*term; for (i=0; i<2; i++) pn[i+4]=b*pn[i+2]-an*pn[i]; if (pn[5] == 0) goto l35; rn=pn[4]/pn[5]; dif=fabs(gin-rn); if (dif>accurate) goto l34; if (dif<=accurate*rn) goto l42; l34: gin=rn; l35: for (i=0; i<4; i++) pn[i]=pn[i+2]; if (fabs(pn[4]) < overflow) goto l32; for (i=0; i<4; i++) pn[i]/=overflow; goto l32; l42: gin=1-factor*gin; l50: return (gin); } /*********************************************************/ double PointChi2 (double prob, double v) { /* returns z so that Prob{x.999998 || v<=0) return (-1); g = LnGamma (v/2); xx=v/2; c=xx-1; if (v >= -1.24*log(p)) goto l1; ch=pow((p*xx*exp(g+xx*aa)), 1/xx); if (ch-e<0) return (ch); goto l4; l1: if (v>.32) goto l3; ch=0.4; a=log(1-p); l2: q=ch; p1=1+ch*(4.67+ch); p2=ch*(6.73+ch*(6.66+ch)); t=-0.5+(4.67+2*ch)/p1 - (6.73+ch*(13.32+3*ch))/p2; ch-=(1-exp(a+g+.5*ch+c*aa)*p2/p1)/t; if (fabs(q/ch-1)-.01 <= 0) goto l4; else goto l2; l3: x=PointNormal (p); p1=0.222222/v; ch=v*pow((x*sqrt(p1)+1-p1), 3.0); if (ch>2.2*v+6) ch=-2*(log(1-p)-c*log(.5*ch)+g); l4: q=ch; p1=.5*ch; if ((t=IncompleteGamma (p1, xx, g))<0) { printf ("\nerr IncompleteGamma"); return (-1); } p2=p-t; t=p2*exp(xx*aa+g+p1-c*log(ch)); b=t/ch; a=0.5*t-b*c; s1=(210+a*(140+a*(105+a*(84+a*(70+60*a))))) / 420; s2=(420+a*(735+a*(966+a*(1141+1278*a))))/2520; s3=(210+a*(462+a*(707+932*a)))/2520; s4=(252+a*(672+1182*a)+c*(294+a*(889+1740*a)))/5040; s5=(84+264*a+c*(175+606*a))/2520; s6=(120+c*(346+127*c))/5040; ch+=t*(1+0.5*t*s1-b*c*(s1-b*(s2-b*(s3-b*(s4-b*(s5-b*s6)))))); if (fabs(q/ch-1) > e) goto l4; return (ch); } /*********************************************************/ double PointNormal (double prob) { /* returns z so that Prob{xn_otu = n_otu; Init_Tree(tree); tree->noeud[n_otu]=(node *)mCalloc(1,sizeof(node)); tree->noeud[n_otu]->v = NULL; Make_Node_Light(tree->noeud[n_otu]); tree->noeud[n_otu]->num=n_otu; tree->noeud[n_otu]->tax=0; subs = Sub_Trees(s_tree,°ree); Clean_Multifurcation(subs,degree,3); if(!strlen(subs[2])) Exit("\n. Err: unrooted tree is needed\n"); sub_tp=(char *)mCalloc((int)strlen(s_tree),sizeof(char)); tree->has_branch_lengths = 1; For(i,degree) { strcpy(sub_tp,subs[i]); strcat(sub_tp,":"); if(strstr(s_tree,sub_tp)) { tree->noeud[n_otu]->l[i] = atof((char *)strstr(s_tree,sub_tp)+strlen(subs[i])+1); } else tree->has_branch_lengths = 0; } Free(sub_tp); n_int = n_ext = 0; For(i,degree) R_rtree(subs[i],tree->noeud[n_otu],tree,&n_int,&n_ext); Make_All_Edges_Light(tree->noeud[0],tree->noeud[0]->v[0]); i = 0; Init_Tree_Edges(tree->noeud[0],tree->noeud[0]->v[0],tree,&i); For(i,NODE_DEG_MAX) Free(subs[i]); Free(subs); return tree; } /*********************************************************/ void Make_All_Edges_Light(node *a, node *d) { int i; Make_Edge_Light(a,d); if(d->tax) return; else { For(i,3) { if(d->v[i] != a) Make_All_Edges_Light(d,d->v[i]); } } } /*********************************************************/ void Make_All_Edges_Lk(node *a, node *d, arbre *tree) { int i; Make_Edge_Lk(a,d,tree); if(d->tax) return; else { For(i,3) { if(d->v[i] != a) Make_All_Edges_Lk(d,d->v[i],tree); } } } /*********************************************************/ void R_rtree(char *s_tree, node *pere, arbre *tree, int *n_int, int *n_ext) { int i; char *sub_tp; node *fils; int n_otu = tree->n_otu; if(strstr(s_tree," ")) Exit("\n Err : tree must not contain a ' ' character\n"); fils=(node *)mCalloc(1,sizeof(node)); fils->v = NULL; Make_Node_Light(fils); if(s_tree[0] == '(') { char **subs; int degree; (*n_int)+=1; tree->noeud[n_otu+*n_int]=fils; fils->num=n_otu+*n_int; fils->tax=0; if(s_tree[(int)strlen(s_tree)-1] == '*') { fils->check_branch = 1; s_tree[(int)strlen(s_tree)-1] = '\0'; } For(i,3) { if(!pere->v[i]) { pere->v[i]=fils; fils->l[0]=pere->l[i]; break; } } fils->v[0]=pere; subs=Sub_Trees(s_tree,°ree); Clean_Multifurcation(subs,degree,2); sub_tp = (char *)mCalloc(T_MAX_LINE,sizeof(char)); strcpy(sub_tp,subs[0]); strcat(sub_tp,":"); if(strstr(s_tree,sub_tp)) { fils->l[1] = atof((char *)strstr(s_tree,sub_tp) +(int)strlen(subs[0])+1); } strcpy(sub_tp,subs[1]); strcat(sub_tp,":"); if(strstr(s_tree,sub_tp)) { fils->l[2] = atof((char *)strstr(s_tree,sub_tp) +(int)strlen(subs[1])+1); } Free(sub_tp); R_rtree(subs[0],fils,tree,n_int,n_ext); R_rtree(subs[1],fils,tree,n_int,n_ext); For(i,NODE_DEG_MAX) Free(subs[i]); Free(subs); } else { tree->noeud[*n_ext]=fils; fils->tax=1; For(i,3) { if(!pere->v[i]) { pere->v[i]=fils; fils->l[0]=pere->l[i]; break; } } if(s_tree[(int)strlen(s_tree)-1] == '*') { fils->check_branch = 1; s_tree[(int)strlen(s_tree)-1] = '\0'; } fils->v[0]=pere; strcpy(fils->name,s_tree); fils->num=*n_ext; (*n_ext)+=1; } } /*********************************************************/ void Clean_Multifurcation(char **subtrees, int current_deg, int end_deg) { if(current_deg <= end_deg) return; else { char *s_tmp; int i; s_tmp = (char *)mCalloc(T_MAX_LINE,sizeof(char)); strcat(s_tmp,"(\0"); strcat(s_tmp,subtrees[0]); strcat(s_tmp,",\0"); strcat(s_tmp,subtrees[1]); strcat(s_tmp,")\0"); Free(subtrees[0]); subtrees[0] = s_tmp; for(i=1;inoeud[tree->n_otu+i]->v[0]) || (!tree->noeud[tree->n_otu+i]->v[1]) || (!tree->noeud[tree->n_otu+i]->v[2])) i++; R_wtree(tree->noeud[tree->n_otu+i],tree->noeud[tree->n_otu+i]->v[0],s,tree); R_wtree(tree->noeud[tree->n_otu+i],tree->noeud[tree->n_otu+i]->v[1],s,tree); R_wtree(tree->noeud[tree->n_otu+i],tree->noeud[tree->n_otu+i]->v[2],s,tree); s[(int)strlen(s)-1]=')'; s[(int)strlen(s)]=';'; return s; } /*********************************************************/ void R_wtree(node *pere, node *fils, char *s_tree, arbre *tree) { int i,p; p = -1; if(fils->tax) { strcat(s_tree,fils->name); if((fils->b[0]) && (fils->b[0]->l != -1)) { strcat(s_tree,":"); sprintf(s_tree+(int)strlen(s_tree),"%f",fils->b[0]->l); fflush(stdout); } sprintf(s_tree+(int)strlen(s_tree),","); } else { s_tree[(int)strlen(s_tree)]='('; For(i,3) { if(fils->v[i] != pere) R_wtree(fils,fils->v[i],s_tree,tree); else p=i; } s_tree[(int)strlen(s_tree)-1]=')'; if(fils->b[0]->l != -1) { if(tree->print_boot_val) sprintf(s_tree+(int)strlen(s_tree),"%d",fils->b[p]->bip_score); strcat(s_tree,":"); sprintf(s_tree+(int)strlen(s_tree),"%f,",fils->b[p]->l); fflush(stdout); } } } /*********************************************************/ void Init_Tree(arbre *tree) { int i; tree->noeud = (node **)mCalloc(2*tree->n_otu-2,sizeof(node *)); tree->t_edges = (edge **)mCalloc(2*tree->n_otu-3,sizeof(edge *)); For(i,2*tree->n_otu-3) tree->t_edges[i] = NULL; tree->best_tree = NULL; tree->old_tree = NULL; tree->has_bip = 0; tree->best_loglk = UNLIKELY; tree->tot_loglk = UNLIKELY; tree->n_swap = 0; tree->tot_dloglk = NULL; tree->tot_d2loglk = NULL; tree->n_pattern = -1; tree->print_boot_val = 0; } /*********************************************************/ void Make_Edge_Light(node *a, node *d) { edge *b; b = (edge *)mCalloc(1,sizeof(edge)); Init_Edge_Light(b); b->left = a; b->rght = d; if(a->tax) {b->rght = a; b->left = d;} /* root */ /* a tip is necessary on the right side of the edge */ Make_Edge_Dirs(b,a,d); b->l = a->l[b->l_r]; if(a->tax) b->l = a->l[b->r_l]; if(b->l < BL_MIN) b->l = BL_MIN; b->l_old = b->l; } /*********************************************************/ void Init_Edge_Light(edge *b) { b->bip_score = 0; b->nj_score = .0; b->site_dlk = -1.; b->site_d2lk = -1.; b->site_dlk_rr = NULL; b->site_d2lk_rr = NULL; b->p_lk_left = NULL; b->p_lk_rght = NULL; b->Pij_rr = NULL; b->dPij_rr = NULL; b->d2Pij_rr = NULL; } /*********************************************************/ void Make_Edge_Dirs(edge *b, node *a, node *d) { int i; b->l_r = b->r_l = -1; For(i,3) { if((a->v[i]) && (a->v[i] == d)) { b->l_r = i; a->b[i] = b; } if((d->v[i]) && (d->v[i] == a)) { b->r_l = i; d->b[i] = b; } } if(a->tax) {b->r_l = 0; For(i,3) if(d->v[i]==a) {b->l_r = i; break;}} b->l_v1 = b->l_v2 = b->r_v1 = b->r_v2 = -1; For(i,3) { if(b->left->v[i] != b->rght) { if(b->l_v1 < 0) b->l_v1 = i; else b->l_v2 = i; } if(b->rght->v[i] != b->left) { if(b->r_v1 < 0) b->r_v1 = i; else b->r_v2 = i; } } } /*********************************************************/ void Make_Edge_Lk(node *a, node *d, arbre *tree) { int i,j; edge *b; b = NULL; For(i,3) if((a->v[i]) && (a->v[i] == d)) {b = a->b[i]; break;} b->diff_lk = 0.0; b->l_old = b->l; if(!b->Pij_rr) { b->best_conf = 1; b->ql = (double *)mCalloc(3,sizeof(double)); b->Pij_rr = (double ***)mCalloc(tree->mod->n_catg,sizeof(double **)); b->dPij_rr = (double ***)mCalloc(tree->mod->n_catg,sizeof(double **)); b->d2Pij_rr = (double ***)mCalloc(tree->mod->n_catg,sizeof(double **)); For(i,tree->mod->n_catg) { b->Pij_rr[i] = (double **)mCalloc(tree->mod->ns,sizeof(double *)); b->dPij_rr[i] = (double **)mCalloc(tree->mod->ns,sizeof(double *)); b->d2Pij_rr[i] = (double **)mCalloc(tree->mod->ns,sizeof(double *)); For(j,tree->mod->ns) { b->Pij_rr[i][j] = (double *)mCalloc(tree->mod->ns,sizeof(double )); b->dPij_rr[i][j] = (double *)mCalloc(tree->mod->ns,sizeof(double )); b->d2Pij_rr[i][j] = (double *)mCalloc(tree->mod->ns,sizeof(double )); } } b->site_p_lk_left = (double **)mCalloc((int)tree->mod->n_catg,sizeof(double *)); For(i,tree->mod->n_catg) b->site_p_lk_left[i] = (double *)mCalloc(tree->mod->ns,sizeof(double)); b->site_p_lk_rght = (double **)mCalloc((int)tree->mod->n_catg,sizeof(double *)); For(i,tree->mod->n_catg) b->site_p_lk_rght[i] = (double *)mCalloc(tree->mod->ns,sizeof(double)); b->site_dlk_rr = (double *)mCalloc((int)tree->mod->n_catg,sizeof(double)); b->site_d2lk_rr = (double *)mCalloc((int)tree->mod->n_catg,sizeof(double)); b->scale_left = b->scale_rght = 0; b->p_lk_left = NULL; b->p_lk_rght = NULL; b->sum_scale_f_rght = NULL; b->sum_scale_f_left = NULL; b->get_p_lk_left = 0; b->get_p_lk_rght = 0; b->ud_p_lk_left = 0; b->ud_p_lk_rght = 0; } } /*********************************************************/ void Make_Node_Light(node *n) { if(n->v) return; n->v = (node **)mCalloc(3,sizeof(node *)); n->l = (double *)mCalloc(3,sizeof(double)); n->b = (edge **)mCalloc(3,sizeof(edge *)); n->name = (char *)mCalloc(T_MAX_NAME,sizeof(char)); n->score = (double *)mCalloc(3,sizeof(double)); Init_Node_Light(n); } /*********************************************************/ void Init_Node_Light(node *n) { int i; n->check_branch = 0; /* n->is_attach = 0; */ /* n->is_free = 0; */ For(i,3) { n->v[i]=NULL; n->b[i]=NULL; n->l[i]=-1; } } /*********************************************************/ seq **Get_Seq(option *input) { seq **data; int i,j; char **buff; int n_unkn,n_removed,pos; int *remove; /* rewind(fp_seq); */ if(input->interleaved) data = Read_Seq_Interleaved(input->fp_seq,&(input->mod->n_otu)); else data = Read_Seq_Sequential(input->fp_seq,&(input->mod->n_otu)); if(data) { buff = (char **)mCalloc(input->mod->n_otu,sizeof(char *)); For(i,input->mod->n_otu) buff[i] = (char *)mCalloc(data[0]->len,sizeof(char)); remove = (int *)mCalloc(data[0]->len,sizeof(int)); n_removed = 0; For(i,data[0]->len) { For(j,input->mod->n_otu) { if((data[j]->state[i] == '?') || (data[j]->state[i] == '-')) data[j]->state[i] = 'X'; if(data[j]->state[i] == 'U') data[j]->state[i] = 'T'; if((input->mod->datatype == NT) && (data[j]->state[i] == 'N')) data[j]->state[i] = 'X'; } n_unkn = 0; For(j,input->mod->n_otu) if(data[j]->state[i] == 'X') n_unkn++; if(n_unkn == input->mod->n_otu) { remove[i] = 1; n_removed++; } For(j,input->mod->n_otu) buff[j][i] = data[j]->state[i]; } if(n_removed > 0) { if(input->mod->datatype == NT) printf("\n. %d sites are made from completely undetermined states ('X', '-', '?' or 'N')...\n",n_removed); else printf("\n. %d sites are made from completely undetermined states ('X', '-', '?')...\n",n_removed); } pos = 0; For(i,data[0]->len) { /* if(!remove[i]) */ /* { */ For(j,input->mod->n_otu) data[j]->state[pos] = buff[j][i]; pos++; /* } */ } For(i,input->mod->n_otu) data[i]->len = pos; For(i,input->mod->n_otu) Free(buff[i]); Free(buff); Free(remove); } return data; } /*********************************************************/ seq **Read_Seq_Sequential(FILE *in, int *n_otu) { int i; char *line; int len,readok; seq **data; char c; char *format = (char *)mCalloc(20, sizeof(char)); line = (char *)mCalloc(T_MAX_LINE,sizeof(char)); readok = len = 0; do { if(fscanf(in,"%s",line) == EOF) { Free(line); return NULL; } else { if(strcmp(line,"\n") && strcmp(line,"\n") && strcmp(line,"\t")) { *n_otu = atoi(line); data = (seq **)mCalloc(*n_otu,sizeof(seq *)); if(*n_otu <= 0) Exit("\n. Problem with sequence format\n"); fscanf(in,"%s",line); len = atoi(line); if(len <= 0) Exit("\n. Problem with sequence format\n"); else readok = 1; } } }while(!readok); /* while((c=fgetc(in))!='\n'); */ while(((c=fgetc(in))!='\n') && (c != ' ') && (c != '\r') && (c != '\t')); For(i,*n_otu) { data[i] = (seq *)mCalloc(1,sizeof(seq)); data[i]->len = 0; data[i]->name = (char *)mCalloc(T_MAX_NAME,sizeof(char)); data[i]->state = (char *)mCalloc(T_MAX_SEQ,sizeof(char)); sprintf(format, "%%%ds", T_MAX_NAME); fscanf(in, format, data[i]->name); while(data[i]->len < len) Read_One_Line_Seq(&data,i,in); if(data[i]->len != len) { printf("\n. Err: Problem with species %s's sequence (check the format)\n", data[i]->name); Exit(""); } } /* fgets(line,T_MAX_LINE,in); */ /* inter data sets */ Free(format); Free(line); return data; } /*********************************************************/ seq **Read_Seq_Interleaved(FILE *in, int *n_otu) { int i,end,num_block; char *line; int len,readok; seq **data; char c; char *format; line = (char *)mCalloc(T_MAX_LINE,sizeof(char)); format = (char *)mCalloc(T_MAX_NAME, sizeof(char)); readok = len = 0; do { if(fscanf(in,"%s",line) == EOF) { Free(format); Free(line); return NULL; } else { if(strcmp(line,"\n") && strcmp(line,"\r") && strcmp(line,"\t")) { *n_otu = atoi(line); data = (seq **)mCalloc(*n_otu,sizeof(seq *)); if(*n_otu <= 0) Exit("\n. Problem with sequence format\n"); fscanf(in,"%s",line); len = atoi(line); if(len <= 0) Exit("\n. Problem with sequence format\n"); else readok = 1; } } }while(!readok); while(((c=fgetc(in))!='\n') && (c != ' ') && (c != '\r') && (c != '\t')); end = 0; For(i,*n_otu) { data[i] = (seq *)mCalloc(1,sizeof(seq)); data[i]->len = 0; data[i]->name = (char *)mCalloc(T_MAX_NAME,sizeof(char)); data[i]->state = (char *)mCalloc(T_MAX_SEQ,sizeof(char)); sprintf(format, "%%%ds", T_MAX_NAME); fscanf(in, format, data[i]->name); if(!Read_One_Line_Seq(&data,i,in)) { end = 1; if((i != *n_otu) && (i != *n_otu-1)) { printf("\n. Err: Problem with species %s's sequence\n",data[i]->name); Exit(""); } break; } } if(data[0]->len == len) end = 1; if(!end) { end = 0; num_block = 1; do { num_block++; /* interblock */ if(!fgets(line,T_MAX_LINE,in)) break; if(line[0] != 13 && line[0] != 10) { printf("\n. One or more missing sequences in block %d\n",num_block-1); Exit(""); } For(i,*n_otu) if(data[i]->len != len) break; if(i == *n_otu) break; For(i,*n_otu) { if(data[i]->len > len) { printf("\n. Err: Problem with species %s's sequence\n",data[i]->name); Exit(""); } else if(!Read_One_Line_Seq(&data,i,in)) { end = 1; if((i != *n_otu) && (i != *n_otu-1)) { printf("\n. Err: Problem with species %s's sequence\n",data[i]->name); Exit(""); } break; } } }while(!end); } For(i,*n_otu) { if(data[i]->len != len) { printf("\n. Check sequence '%s' length...\n",data[i]->name); Exit(""); } } Free(format); Free(line); return data; } /*********************************************************/ int Read_One_Line_Seq(seq ***data, int num_otu, FILE *in) { char c; c=' '; while(1) { /* if((c == EOF) || (c == '\n') || (c == '\r')) break; */ if((c == EOF) || (c == 13) || (c == 10)) break; else if((c==' ') || (c=='\t')) {c=(char)fgetc(in); continue;} Uppercase(&c); /*if(strchr("ACGTUMRWSYKBDHVNXO?-.",c) == NULL)*/ if (strchr("ABCDEFGHIKLMNOPQRSTUVWXYZ?-.", c) == NULL) { printf("\n. Err: bad symbol: \"%c\" at position %d of species %s\n", c,(*data)[num_otu]->len,(*data)[num_otu]->name); Exit(""); } if(c == '.') { c = (*data)[0]->state[(*data)[num_otu]->len]; if(!num_otu) Exit("\n. Err: Symbol \".\" should not appear in the first sequence\n"); } (*data)[num_otu]->state[(*data)[num_otu]->len]=c; (*data)[num_otu]->len++; /* if(c=='U') c='T'; */ c = (char)fgetc(in); } if(c == EOF) return 0; else return 1; } /*********************************************************/ void Uppercase(char *ch) { /* convert ch to upper case -- either ASCII or EBCDIC */ *ch = isupper((int)*ch) ? *ch : toupper((int)*ch); } /*********************************************************/ allseq *Compact_Seq(seq **data, option *input) { allseq *alldata; int i,j,k/*,diff*/,site; int n_patt,which_patt,n_invar; char **sp_names; int n_otu; n_otu = input->mod->n_otu; sp_names = (char **)mCalloc(n_otu,sizeof(char *)); For(i,n_otu) { sp_names[i] = (char *)mCalloc(T_MAX_NAME,sizeof(char)); strcpy(sp_names[i],data[i]->name); } alldata = Make_Seq(n_otu,data[0]->len,sp_names); For(i,n_otu) Free(sp_names[i]); Free(sp_names); n_patt = which_patt = 0; /* diff = -1;*/ if(data[0]->len%input->mod->stepsize) { printf("\n. Sequence length is not a multiple of %d\n",input->mod->stepsize); Exit(""); } Fors(site,data[0]->len,input->mod->stepsize) { Fors(k,n_patt,input->mod->stepsize) { For(j,n_otu) { if(/*!Compare_Two_States*/strncmp(alldata->c_seq[j]->state+k, data[j]->state+site, input->mod->stepsize)) break; } if(j == n_otu) { which_patt = k; break; } } /* k = n_patt; */ if(k == n_patt) { For(j,n_otu) Copy_One_State(data[j]->state+site, alldata->c_seq[j]->state+n_patt, input->mod->stepsize); for(j=0;jc_seq[j]->state+n_patt,input->mod->datatype,input->mod->stepsize) || */ /* (/\*!Compare_Two_States*\/strncmp(alldata->c_seq[j]->state+n_patt, */ /* alldata->c_seq[0]->state+n_patt, */ /* input->mod->stepsize)))) */ if(!(Are_Compatible(alldata->c_seq[j]->state+n_patt, alldata->c_seq[0]->state+n_patt, input->mod->stepsize, input->mod->datatype))) break; } if(j==n_otu) { For(j,n_otu) { alldata->invar[n_patt] = Assign_State(alldata->c_seq[j]->state+n_patt, input->mod->datatype, input->mod->stepsize); break; } } else alldata->invar[n_patt] = -1; /* Print_Site(alldata,k,n_otu,"\n",input->mod->stepsize); */ alldata->sitepatt[site] = n_patt; alldata->wght[n_patt] += 1.; n_patt+=input->mod->stepsize; } else { alldata->sitepatt[site] = which_patt; alldata->wght[which_patt] += 1.; } } alldata->init_len = data[0]->len; alldata->crunch_len = n_patt; For(i,n_otu) alldata->c_seq[i]->len = n_patt; /* fprintf(stderr,"%d patterns found\n",n_patt); */ /* For(site,alldata->crunch_len) printf("%1.0f",alldata->wght[site]); */ n_invar=0; For(i,alldata->crunch_len) if(alldata->invar[i]>-1) n_invar+=(int)alldata->wght[i]; if(input->mod->datatype == NT) Get_Base_Freqs(alldata); else Get_AA_Freqs(alldata); /* fprintf(stderr,"Average nucleotides frequencies : \n"); */ /* fprintf(stderr,"%f %f %f %f\n", */ /* alldata->b_frq[0], */ /* alldata->b_frq[1], */ /* alldata->b_frq[2], */ /* alldata->b_frq[3]); */ return alldata; } /*********************************************************/ allseq *Compact_CSeq(allseq *data, model *mod) { allseq *alldata; int i,j,k,site; int n_patt,which_patt; int n_otu; n_otu = data->n_otu; alldata = (allseq *)mCalloc(1,sizeof(allseq)); alldata->n_otu=n_otu; alldata->c_seq = (seq **)mCalloc(n_otu,sizeof(seq *)); alldata->wght = (double *)mCalloc(data->crunch_len,sizeof(double)); alldata->b_frq = (double *)mCalloc(mod->ns,sizeof(double)); alldata->ambigu = (int *)mCalloc(data->crunch_len,sizeof(int)); alldata->invar = (int *)mCalloc(data->crunch_len,sizeof(int)); alldata->crunch_len = alldata->init_len = -1; For(j,n_otu) { alldata->c_seq[j] = (seq *)mCalloc(1,sizeof(seq)); alldata->c_seq[j]->name = (char *)mCalloc(T_MAX_NAME,sizeof(char)); strcpy(alldata->c_seq[j]->name,data->c_seq[j]->name); alldata->c_seq[j]->state = (char *)mCalloc(data->crunch_len,sizeof(char)); alldata->c_seq[j]->state[0] = data->c_seq[j]->state[0]; } n_patt = which_patt = 0; Fors(site,data->crunch_len,mod->stepsize) { Fors(k,n_patt,mod->stepsize) { For(j,n_otu) { if(/*!Compare_Two_States*/strncmp(alldata->c_seq[j]->state+k, data->c_seq[j]->state+site, mod->stepsize)) break; } if(j == n_otu) { which_patt = k; break; } } if(k == n_patt) { For(j,n_otu) Copy_One_State(data->c_seq[j]->state+site, alldata->c_seq[j]->state+n_patt, mod->stepsize); for(j=1;jc_seq[j]->state+n_patt, alldata->c_seq[j-1]->state+n_patt, mod->stepsize)) break; if(j==n_otu) alldata->invar[n_patt] = 1; alldata->wght[n_patt] += data->wght[site]; n_patt+=mod->stepsize; } else alldata->wght[which_patt] += data->wght[site]; /* Print_Site(alldata,k,n_otu,"\n",mod->stepsize); */ } alldata->init_len = data->crunch_len; alldata->crunch_len = n_patt; For(i,n_otu) alldata->c_seq[i]->len = n_patt; (mod->datatype == NT)? (Get_Base_Freqs(alldata)): (Get_AA_Freqs(alldata)); return alldata; } /*********************************************************/ void Get_Base_Freqs(allseq *data) { int i,j,k; double A,C,G,T; double fA,fC,fG,fT; double w; fA = fC = fG = fT = .25; For(k,8) { A = C = G = T = .0; For(i,data->n_otu) { For(j,data->crunch_len) { w = data->wght[j]; if(w) { switch(data->c_seq[i]->state[j]){ case 'A' : A+=w; break; case 'C' : C+=w; break; case 'G' : G+=w; break; case 'T' : T+=w; break; case 'U' : T+=w; break; case 'M' : C+=w*fC/(fC+fA); A+=w*fA/(fA+fC); break; case 'R' : G+=w*fG/(fA+fG); A+=w*fA/(fA+fG); break; case 'W' : T+=w*fT/(fA+fT); A+=w*fA/(fA+fT); break; case 'S' : C+=w*fC/(fC+fG); G+=w*fG/(fC+fG); break; case 'Y' : C+=w*fC/(fC+fT); T+=w*fT/(fT+fC); break; case 'K' : G+=w*fG/(fG+fT); T+=w*fT/(fT+fG); break; case 'B' : C+=w*fC/(fC+fG+fT); G+=w*fG/(fC+fG+fT); T+=w*fT/(fC+fG+fT); break; case 'D' : A+=w*fA/(fA+fG+fT); G+=w*fG/(fA+fG+fT); T+=w*fT/(fA+fG+fT); break; case 'H' : A+=w*fA/(fA+fC+fT); C+=w*fC/(fA+fC+fT); T+=w*fT/(fA+fC+fT); break; case 'V' : A+=w*fA/(fA+fC+fG); C+=w*fC/(fA+fC+fG); G+=w*fG/(fA+fC+fG); break; case 'N' : case 'X' : case '?' : case 'O' : case '-' : A+=w*fA; C+=w*fC; G+=w*fG; T+=w*fT; break; default : break; } } } } fA = A/(A+C+G+T); fC = C/(A+C+G+T); fG = G/(A+C+G+T); fT = T/(A+C+G+T); } data->b_frq[0] = fA; data->b_frq[1] = fC; data->b_frq[2] = fG; data->b_frq[3] = fT; } /*********************************************************/ void Get_AA_Freqs(allseq *data) { int i,j,k; double A,C,D,E,F,G,H,I,K,L,M,N,P,Q,R,S,T,V,W,Y; double fA,fC,fD,fE,fF,fG,fH,fI,fK,fL,fM,fN,fP,fQ,fR,fS,fT,fV,fW,fY; double w; double sum; fA = fC = fD = fE = fF = fG = fH = fI = fK = fL = fM = fN = fP = fQ = fR = fS = fT = fV = fW = fY = 1./20.; For(k,8) { A = C = D = E = F = G = H = I = K = L = M = N = P = Q = R = S = T = V = W = Y = .0; For(i,data->n_otu) { For(j,data->crunch_len) { w = data->wght[j]; if(w) { switch(data->c_seq[i]->state[j]){ case 'A' : A+=w; break; case 'C' : C+=w; break; case 'D' : D+=w; break; case 'E' : E+=w; break; case 'F' : F+=w; break; case 'G' : G+=w; break; case 'H' : H+=w; break; case 'I' : I+=w; break; case 'K' : K+=w; break; case 'L' : L+=w; break; case 'M' : M+=w; break; case 'N' : N+=w; break; case 'P' : P+=w; break; case 'Q' : Q+=w; break; case 'R' : R+=w; break; case 'S' : S+=w; break; case 'T' : T+=w; break; case 'V' : V+=w; break; case 'W' : W+=w; break; case 'Y' : Y+=w; break; case 'Z' : Q+=w; break; case 'X' : case '?' : case 'O' : case '-' : A+=w*fA; C+=w*fC; D+=w*fD; E+=w*fE; F+=w*fF; G+=w*fG; H+=w*fH; I+=w*fI; K+=w*fK; L+=w*fL; M+=w*fM; N+=w*fN; P+=w*fP; Q+=w*fQ; R+=w*fR; S+=w*fS; T+=w*fT; V+=w*fV; W+=w*fW; Y+=w*fY; break; default : break; } } } } sum = (A+C+D+E+F+G+H+I+K+L+M+N+P+Q+R+S+T+V+W+Y); fA = A/sum; fC = C/sum; fD = D/sum; fE = E/sum; fF = F/sum; fG = G/sum; fH = H/sum; fI = I/sum; fK = K/sum; fL = L/sum; fM = M/sum; fN = N/sum; fP = P/sum; fQ = Q/sum; fR = R/sum; fS = S/sum; fT = T/sum; fV = V/sum; fW = W/sum; fY = Y/sum; } data->b_frq[0] = fA; data->b_frq[1] = fR; data->b_frq[2] = fN; data->b_frq[3] = fD; data->b_frq[4] = fC; data->b_frq[5] = fQ; data->b_frq[6] = fE; data->b_frq[7] = fG; data->b_frq[8] = fH; data->b_frq[9] = fI; data->b_frq[10] = fL; data->b_frq[11] = fK; data->b_frq[12] = fM; data->b_frq[13] = fF; data->b_frq[14] = fP; data->b_frq[15] = fS; data->b_frq[16] = fT; data->b_frq[17] = fW; data->b_frq[18] = fY; data->b_frq[19] = fV; } /*********************************************************/ arbre *Read_Tree_File(FILE *fp_input_tree) { char *line; arbre *tree; int i; char c; line = (char *)mCalloc(T_MAX_LINE,sizeof(char)); do c=fgetc(fp_input_tree); while((c != '(') && (c != EOF)); if(c==EOF) { Free(line); return NULL; } i=0; for(;;) { if((c == ' ') || (c == '\n')) { c=fgetc(fp_input_tree); if(c==EOF) break; else continue; } line[i]=c; i++; c=fgetc(fp_input_tree); if(c==EOF || c==';') break; } tree = Read_Tree(line); Free(line); return tree; } /*********************************************************/ void Init_Tree_Edges(node *a, node *d, arbre *tree, int *cur) { int i,dir_a_d; dir_a_d = -1; For(i,3) if(a->v[i] == d) {dir_a_d = i; break;} tree->t_edges[*cur] = a->b[dir_a_d]; tree->t_edges[*cur]->num = *cur; *cur = *cur + 1; if(d->tax) return; else { For(i,3) { if(d->v[i] != a) Init_Tree_Edges(d,d->v[i],tree,cur); } } } /*********************************************************/ void Exit(const char *message) { fprintf(stderr,"%s",message); exit(1); } /*********************************************************/ void *mCalloc(int nb, size_t size) { void *allocated; if((allocated = calloc((size_t)nb,(size_t)size)) != NULL) { return allocated; } else Exit("\n. Err: low memory\n"); return NULL; } /*********************************************************/ void *mRealloc(void *p,int nb, size_t size) { if((p = realloc(p,(size_t)nb*size)) != NULL) return p; else Exit("\n. Err: low memory\n"); return NULL; } /*********************************************************/ arbre *Make_Light_Tree_Struct(int n_otu) { arbre *tree; int i; tree = (arbre *)mCalloc(1,sizeof(arbre )); tree->t_edges = (edge **)mCalloc(2*n_otu-3,sizeof(edge *)); tree->noeud = (node **)mCalloc(2*n_otu-2,sizeof(node *)); tree->n_otu = n_otu; For(i,2*n_otu-3) { tree->t_edges[i] = (edge *)mCalloc(1,sizeof(edge)); Init_Edge_Light(tree->t_edges[i]); } For(i,2*n_otu-2) { tree->noeud[i] = (node *)mCalloc(1,sizeof(node)); tree->noeud[i]->v = NULL; Make_Node_Light(tree->noeud[i]); } return tree; } /*********************************************************/ int Sort_Double_Decrease(const void *a, const void *b) { if((*(double *)(a)) >= (*(double *)(b))) return -1; else return 1; } /*********************************************************/ void qksort(double* A, int ilo, int ihi) { double pivot; // pivot value for partitioning array int ulo, uhi; // indices at ends of unpartitioned region int ieq; // least index of array entry with value equal to pivot double tempEntry; // temporary entry used for swapping if (ilo >= ihi) { return; } // Select a pivot value. pivot = A[(ilo + ihi)/2]; // Initialize ends of unpartitioned region and least index of entry // with value equal to pivot. ieq = ulo = ilo; uhi = ihi; // While the unpartitioned region is not empty, try to reduce its size. while (ulo <= uhi) { if (A[uhi] > pivot) { // Here, we can reduce the size of the unpartitioned region and // try again. uhi--; } else { // Here, A[uhi] <= pivot, so swap entries at indices ulo and // uhi. tempEntry = A[ulo]; A[ulo] = A[uhi]; A[uhi] = tempEntry; // After the swap, A[ulo] <= pivot. if (A[ulo] < pivot) { // Swap entries at indices ieq and ulo. tempEntry = A[ieq]; A[ieq] = A[ulo]; A[ulo] = tempEntry; // After the swap, A[ieq] < pivot, so we need to change // ieq. ieq++; // We also need to change ulo, but we also need to do // that when A[ulo] = pivot, so we do it after this if // statement. } // Once again, we can reduce the size of the unpartitioned // region and try again. ulo++; } } // Now, all entries from index ilo to ieq - 1 are less than the pivot // and all entries from index uhi to ihi + 1 are greater than the // pivot. So we have two regions of the array that can be sorted // recursively to put all of the entries in order. qksort(A, ilo, ieq - 1); qksort(A, uhi + 1, ihi); } /********************************************************/ void Print_Site(allseq *alldata, int num, int n_otu, char *sep, int stepsize) { int i,j; For(i,n_otu) { printf("%s ",alldata->c_seq[i]->name); For(j,stepsize) printf("%c",alldata->c_seq[i]->state[num+j]); printf("%s",sep); } fprintf(stderr,"%s",sep); } /*********************************************************/ void Print_Site_Lk(arbre *tree, FILE *fp) { int site; int catg; char *s; s = (char *)mCalloc(T_MAX_LINE,sizeof(char)); fprintf(fp,"Note : P(D|M) is the probability of site D given the model M (i.e., the site likelihood)\n"); if(tree->mod->n_catg > 1 || tree->mod->invar) fprintf(fp,"P(D|M,rr[x]) is the probability of site D given the model M and the relative rate\nof evolution rr[x], where x is the class of rate to be considered.\nWe have P(D|M) = \\sum_x P(x) x P(D|M,rr[x]).\n"); fprintf(fp,"\n\n"); sprintf(s,"Site"); fprintf(fp, "%-7s",s); sprintf(s,"P(D|M)"); fprintf(fp,"%-16s",s); if(tree->mod->n_catg > 1) { For(catg,tree->mod->n_catg) { sprintf(s,"P(D|M,rr[%d]=%5.4f)",catg+1,tree->mod->rr[catg]); fprintf(fp,"%-22s",s); } } if(tree->mod->invar) { sprintf(s,"P(D|M,rr[0]=0)"); fprintf(fp,"%-16s",s); } fprintf(fp,"\n"); For(site,tree->data->init_len) { fprintf(fp,"%-7d",site+1); fprintf(fp,"%-16g",exp(tree->site_lk[tree->data->sitepatt[site]])); if(tree->mod->n_catg > 1) { For(catg,tree->mod->n_catg) fprintf(fp,"%-22g",exp(tree->log_site_lk_cat[catg][tree->data->sitepatt[site]])); } if(tree->mod->invar) { if ((double)tree->data->invar[site] > -0.5) fprintf(fp,"%-16g",tree->mod->pi[tree->data->invar[tree->data->sitepatt[site]]]); else fprintf(fp,"%-16g",0.0); } fprintf(fp,"\n"); } Free(s); } /*********************************************************/ void Print_Seq(seq **data, int n_otu) { int i,j; printf("%d\t%d\n",n_otu,data[0]->len); For(i,n_otu) { For(j,23) { if(j<(int)strlen(data[i]->name)) putchar(data[i]->name[j]); else putchar(' '); } For(j,data[i]->len) /*FLT uncommented*/ /* For(j,2000)*//*FLT commented*/ { printf("%c",data[i]->state[j]); } printf("\n"); } } /*********************************************************/ void Print_CSeq(FILE *fp, allseq *alldata) { int i,j,k; int n_otu; n_otu = alldata->n_otu; fprintf(fp,"%d\t%d\n",n_otu,alldata->init_len); For(i,n_otu) { For(j,23) { if(j<(int)strlen(alldata->c_seq[i]->name)) fputc(alldata->c_seq[i]->name[j],fp); else fputc(' ',fp); } For(j,alldata->crunch_len) { For(k,alldata->wght[j]) fprintf(fp,"%c",alldata->c_seq[i]->state[j]); } fprintf(fp,"\n"); } fprintf(fp,"\n"); /* printf("\t"); */ /* For(j,alldata->crunch_len) */ /* printf("%.0f ",alldata->wght[j]); */ /* printf("\n"); */ } /*********************************************************/ void Order_Tree_Seq(arbre *tree, seq **data) { int i,j,n_otu; seq *buff; n_otu = tree->n_otu; For(i,n_otu) { For(j,n_otu) { if(!strcmp(tree->noeud[i]->name,data[j]->name)) break; } buff = data[j]; data[j] = data[i]; data[i] = buff; } } /*********************************************************/ void Order_Tree_CSeq(arbre *tree, allseq *data) { int i,j,n_otu_tree,n_otu_seq; seq *buff; n_otu_tree = tree->n_otu; n_otu_seq = data->n_otu; if(n_otu_tree != n_otu_seq) { /* printf("%d(tree) != %d(seq) \n",n_otu_tree,n_otu_seq); */ Exit("\n. The number of tips in the tree is not the same as the number of sequences\n"); } For(i,MAX(n_otu_tree,n_otu_seq)) { For(j,MIN(n_otu_tree,n_otu_seq)) { if(!strcmp(tree->noeud[i]->name,data->c_seq[j]->name)) break; } if(j==MIN(n_otu_tree,n_otu_seq)) { printf("\n. Err: %s is not found in sequences data set\n", tree->noeud[i]->name); Exit(""); } buff = data->c_seq[j]; data->c_seq[j] = data->c_seq[i]; data->c_seq[i] = buff; } } /*********************************************************/ matrix *Make_Mat(int n_otu) { matrix *mat; int i; mat = (matrix *)mCalloc(1,sizeof(matrix)); mat->n_otu = n_otu; mat->P = (double **)mCalloc(n_otu,sizeof(double *)); mat->Q = (double **)mCalloc(n_otu,sizeof(double *)); mat->dist = (double **)mCalloc(n_otu,sizeof(double *)); mat->on_off = (int *)mCalloc(n_otu,sizeof(int)); mat->name = (char **)mCalloc(n_otu,sizeof(char *)); mat->tip_node = (node **)mCalloc(n_otu,sizeof(node *)); For(i,n_otu) { mat->P[i] = (double *)mCalloc(n_otu,sizeof(double)); mat->Q[i] = (double *)mCalloc(n_otu,sizeof(double)); mat->dist[i] = (double *)mCalloc(n_otu,sizeof(double)); mat->name[i] = (char *)mCalloc(T_MAX_NAME,sizeof(char)); } return mat; } /*********************************************************/ void Init_Mat(matrix *mat, allseq *data) { int i; mat->n_otu = data->n_otu; mat->r = mat->n_otu; mat->curr_int = mat->n_otu; mat->method = 1; For(i,data->n_otu) { strcpy(mat->name[i],data->c_seq[i]->name); mat->on_off[i] = 1; } } /*********************************************************/ arbre *Make_Tree(allseq *data) { arbre *tree; int i; tree = (arbre *)mCalloc(1,sizeof(arbre )); tree->n_otu = data->n_otu; Init_Tree(tree); For(i,2*tree->n_otu-2) { tree->noeud[i] = (node *)mCalloc(1,sizeof(node)); tree->noeud[i]->v = NULL; Make_Node_Light(tree->noeud[i]); } For(i,tree->n_otu) { strcpy(tree->noeud[i]->name,data->c_seq[i]->name); tree->noeud[i]->tax = 1; tree->noeud[i]->num = i; } return tree; } /*********************************************************/ void Print_Dist(matrix *mat) { int i,j; For(i,mat->n_otu) { printf("%s ",mat->name[i]); For(j,mat->n_otu) printf("%9.6f ",mat->dist[i][j]); printf("\n"); } } /*********************************************************/ void Print_Node(node *a, node *d, arbre *tree) { int i; printf("N %2d %2d ",a->num,d->num); For(i,3) if(a->v[i] == d) {printf("%2d %f\n", a->b[i]->num, /* a->b[i]->check_this_one, */ a->b[i]->nj_score); break;} if(d->tax) return; else For(i,3) if(d->v[i] != a) Print_Node(d,d->v[i],tree); } /*********************************************************/ void Share_Lk_Struct(arbre *t_full, arbre *t_empt) { int i,n_otu; edge *b_e,*b_f; n_otu = t_full->n_otu; t_empt->root = t_empt->noeud[0]; t_empt->tot_loglk_sorted = t_full->tot_loglk_sorted; t_empt->log_site_lk_cat = t_full->log_site_lk_cat; t_empt->site_lk = t_full->site_lk; t_empt->tot_dloglk = t_full->tot_dloglk; t_empt->tot_d2loglk = t_full->tot_d2loglk; /* t_empt->mod = t_full->mod; */ /* t_empt->data = t_full->data; */ /* t_empt->mod->s_opt = t_full->mod->s_opt; */ /* For(i,2*n_otu-2) */ /* t_empt->noeud[i]->n_ex_nodes = t_full->noeud[i]->n_ex_nodes; */ For(i,2*n_otu-3) { b_f = t_full->t_edges[i]; b_e = t_empt->t_edges[i]; b_e->Pij_rr = b_f->Pij_rr; b_e->dPij_rr = b_f->dPij_rr; b_e->d2Pij_rr = b_f->d2Pij_rr; b_e->p_lk_left = b_f->p_lk_left; b_e->p_lk_rght = b_f->p_lk_rght; b_e->sum_scale_f_left = b_f->sum_scale_f_left; b_e->sum_scale_f_rght = b_f->sum_scale_f_rght; b_e->site_p_lk_left = b_f->site_p_lk_left; b_e->site_p_lk_rght = b_f->site_p_lk_rght; b_e->site_dlk_rr = b_f->site_dlk_rr; b_e->site_d2lk_rr = b_f->site_d2lk_rr; b_e->ql = b_f->ql; } } /*********************************************************/ /* void Init_Constant() */ /* { */ /* NODE_DEG_MAX = 50; */ /* BRENT_ITMAX = 100; */ /* BRENT_CGOLD = 0.3819660; */ /* BRENT_ZEPS = 1.e-10; */ /* MNBRAK_GOLD = 1.618034; */ /* MNBRAK_GLIMIT = 100.0; */ /* MNBRAK_TINY = 1.e-20; */ /* ALPHA_MIN = 0.04; */ /* ALPHA_MAX = 100; */ /* BL_MIN = 1.e-10; */ /* BL_START = 1.e-03; */ /* BL_MAX = 1.e+05; */ /* MIN_DIFF_LK = 1.e-06; */ /* GOLDEN_R = 0.61803399; */ /* GOLDEN_C = (1.0-GOLDEN_R); */ /* T_MAX_FILE = 200; */ /* T_MAX_LINE = 100000; */ /* T_MAX_NAME = 100; */ /* T_MAX_SEQ = 1000000; */ /* N_MAX_INSERT = 20; */ /* N_MAX_OTU = 4000; */ /* UNLIKELY = -1.e10; */ /* NJ_SEUIL = 0.1; */ /* ROUND_MAX = 100; */ /* DIST_MAX = 2.00; */ /* AROUND_LK = 50.0; */ /* PROP_STEP = 1.0; */ /* T_MAX_ALPHABET = 100; */ /* MDBL_MIN = 2.225074E-308; */ /* MDBL_MAX = 1.797693E+308; */ /* POWELL_ITMAX = 200; */ /* LINMIN_TOL = 2.0E-04; */ /* LIM_SCALE = 3; */ /* LIM_SCALE_VAL = 1.E-50; */ /* /\* LIM_SCALE = 300; *\/ */ /* /\* LIM_SCALE_VAL = 1.E-500; *\/ */ /* } */ /*********************************************************/ void Print_Mat(matrix *mat) { int i,j; printf("%d",mat->n_otu); printf("\n"); For(i,mat->n_otu) { For(j,13) { if(j>=(int)strlen(mat->name[i])) putchar(' '); else putchar(mat->name[i][j]); } For(j,mat->n_otu) { if(mat->dist[i][j] == -1) printf(" - "); else printf("%7.8f ",mat->dist[i][j]); } printf("\n"); } } /*********************************************************/ int Sort_Edges_Diff_Lk(edge **sorted_edges, int n_elem) { int i,j; edge *buff; For(i,n_elem-1) { for(j=i+1;jdiff_lk < sorted_edges[i]->diff_lk) { buff = sorted_edges[j]; sorted_edges[j] = sorted_edges[i]; sorted_edges[i] = buff; } } } return 1; } /*********************************************************/ void NNI(arbre *tree, edge *b_fcus, int do_swap) { node *v2,*v3,*v4; double lk1, lk2, lk3; double lk1_init, lk2_init, lk3_init; double bl_init; double l1,l2,l3; double l_infa, l_infb, l_max; /* double lk_infa, lk_infb, lk_max; */ double lk_init; bl_init = b_fcus->l; lk_init = tree->tot_loglk; b_fcus->best_conf = 1; b_fcus->diff_lk = .0; lk1 = lk2 = lk3 = UNLIKELY; v2 = v3 = v4 = NULL; v2 = b_fcus->left->v[b_fcus->l_v2]; v3 = b_fcus->rght->v[b_fcus->r_v1]; v4 = b_fcus->rght->v[b_fcus->r_v2]; l1 = l2 = l3 = -1.; /***********/ Swap(v2,b_fcus->left,b_fcus->rght,v3); tree->mod->s_opt->opt_bl = 0; tree->both_sides = 1; lk2_init = Update_Lk_At_Given_Edge(b_fcus,tree); l_infa = 10.*b_fcus->l; l_max = b_fcus->l; l_infb = BL_MIN; lk2 = Br_Len_Brent(l_infa,l_max,l_infb, 1.e-6, &(b_fcus->l), b_fcus,tree,1000); if(lk2 < lk2_init - MIN_DIFF_LK) { printf("%f %f %f %f\n",l_infa,l_max,l_infb,b_fcus->l); printf("%f -- %f \n",lk2_init,lk2); printf("\n. Err. in Optimize_Br_Len_Serie\n"); } l2 = b_fcus->l; Swap(v3,b_fcus->left,b_fcus->rght,v2); /***********/ /***********/ Swap(v2,b_fcus->left,b_fcus->rght,v4); b_fcus->l = bl_init; tree->mod->s_opt->opt_bl = 0; tree->both_sides = 1; lk3_init = Update_Lk_At_Given_Edge(b_fcus,tree); l_infa = 10.*b_fcus->l; l_max = b_fcus->l; l_infb = BL_MIN; lk3 = Br_Len_Brent(l_infa,l_max,l_infb, 1.e-6, &(b_fcus->l), b_fcus,tree,1000); if(lk3 < lk3_init - MIN_DIFF_LK) { printf("%f %f %f %f\n",l_infa,l_max,l_infb,b_fcus->l); printf("%f -- %f \n",lk3_init,lk3); printf("\n. Err. in NNI (2)\n"); } l3 = b_fcus->l; Swap(v4,b_fcus->left,b_fcus->rght,v2); /***********/ /***********/ b_fcus->l = bl_init; tree->mod->s_opt->opt_bl = 0; tree->both_sides = 1; lk1_init = Update_Lk_At_Given_Edge(b_fcus,tree); if((lk1_init < lk_init - MIN_DIFF_LK) || (lk1_init > lk_init + MIN_DIFF_LK)) { printf("\n\n. lk_init = %E; lk = %E\n", lk_init, lk1_init); Exit("\n. Err. in NNI (3)\n"); } l_infa = 10.*b_fcus->l; l_max = b_fcus->l; l_infb = BL_MIN; lk1 = Br_Len_Brent(l_infa,l_max,l_infb, 1.e-6, &(b_fcus->l), b_fcus,tree,1000); if(lk1 < lk_init - MIN_DIFF_LK) { printf("\n\n%f %f %f %f\n",l_infa,l_max,l_infb,b_fcus->l); printf("%f -- %f \n",lk1_init,lk1); printf("\n. Err. in NNI (3)\n"); } l1 = b_fcus->l; /***********/ b_fcus->ql[0] = l1; b_fcus->ql[1] = l2; b_fcus->ql[2] = l3; b_fcus->diff_lk = lk1 - MAX(lk2,lk3); if(lk2 > lk3) b_fcus->best_conf = 2; else b_fcus->best_conf = 3; if((do_swap) && ((lk2 > lk1+MDBL_MIN) || (lk3 > lk1+MDBL_MIN))) { tree->n_swap++; printf("Swap edge %d -> %f\n",b_fcus->num,MAX(lk2,lk3)); fflush(stdout); if(lk2 > lk3) { tree->best_loglk = lk2; Swap(v2,b_fcus->left,b_fcus->rght,v3); b_fcus->l = l2; tree->both_sides = 1; Lk(tree,tree->data); } else { tree->best_loglk = lk3; Swap(v2,b_fcus->left,b_fcus->rght,v4); b_fcus->l = l3; tree->both_sides = 1; Lk(tree,tree->data); } } else { b_fcus->l = bl_init; Update_PMat_At_Given_Edge(b_fcus,tree); tree->tot_loglk = lk_init; } } /*********************************************************/ void Swap(node *a, node *b, node *c, node *d) { int ab, ba, cd, dc; int i; /* \ /d \ /a \ / \ / \b__...__c/ -> \b__...__c/ / \ / \ / \ / \ /a \ /d \ */ ab = ba = cd = dc = -1; For(i,3) if(a->v[i] == b) { ab = i; break; } For(i,3) if(b->v[i] == a) { ba = i; break; } For(i,3) if(c->v[i] == d) { cd = i; break; } For(i,3) if(d->v[i] == c) { dc = i; break; } a->v[ab] = c; d->v[dc] = b; b->v[ba] = d; c->v[cd] = a; b->b[ba] = d->b[dc]; c->b[cd] = a->b[ab]; (a->b[ab]->left == b)? (a->b[ab]->left = c): (a->b[ab]->rght = c); (d->b[dc]->left == c)? (d->b[dc]->left = b): (d->b[dc]->rght = b); For(i,3) { if(a->b[ab]->left->v[i] == a->b[ab]->rght) a->b[ab]->l_r = i; if(a->b[ab]->rght->v[i] == a->b[ab]->left) a->b[ab]->r_l = i; if(d->b[dc]->left->v[i] == d->b[dc]->rght) d->b[dc]->l_r = i; if(d->b[dc]->rght->v[i] == d->b[dc]->left) d->b[dc]->r_l = i; } a->b[ab]->l_v1 = a->b[ab]->l_v2 = a->b[ab]->r_v1 = a->b[ab]->r_v2 = d->b[dc]->l_v1 = d->b[dc]->l_v2 = d->b[dc]->r_v1 = d->b[dc]->r_v2 = -1; For(i,3) { if(i != a->b[ab]->l_r) { if(a->b[ab]->l_v1 < 0) a->b[ab]->l_v1 = i; else a->b[ab]->l_v2 = i; } if(i != a->b[ab]->r_l) { if(a->b[ab]->r_v1 < 0) a->b[ab]->r_v1 = i; else a->b[ab]->r_v2 = i; } if(i != d->b[dc]->l_r) { if(d->b[dc]->l_v1 < 0) d->b[dc]->l_v1 = i; else d->b[dc]->l_v2 = i; } if(i != d->b[dc]->r_l) { if(d->b[dc]->r_v1 < 0) d->b[dc]->r_v1 = i; else d->b[dc]->r_v2 = i; } } } /*********************************************************/ void Update_All_Partial_Lk(edge *b_fcus, arbre *tree) { Update_SubTree_Partial_Lk(b_fcus->left->b[b_fcus->l_v1], b_fcus->left, b_fcus->left->v[b_fcus->l_v1], tree); Update_SubTree_Partial_Lk(b_fcus->left->b[b_fcus->l_v2], b_fcus->left, b_fcus->left->v[b_fcus->l_v2], tree); Update_SubTree_Partial_Lk(b_fcus->rght->b[b_fcus->r_v1], b_fcus->rght, b_fcus->rght->v[b_fcus->r_v1], tree); Update_SubTree_Partial_Lk(b_fcus->rght->b[b_fcus->r_v2], b_fcus->rght, b_fcus->rght->v[b_fcus->r_v2], tree); tree->tot_loglk = Lk_At_Given_Edge(tree,b_fcus); } /*********************************************************/ void Update_SubTree_Partial_Lk(edge *b_fcus, node *a, node *d, arbre *tree) { int i; Update_P_Lk(tree,b_fcus,a); if(d->tax) return; else For(i,3) if(d->v[i] != a) Update_SubTree_Partial_Lk(d->b[i],d,d->v[i],tree); } /*********************************************************/ double Update_Lk_At_Given_Edge(edge *b_fcus, arbre *tree) { /* if(b_fcus->l < BL_MIN) b_fcus->l = BL_MIN; */ /* For(i,tree->mod->n_catg) */ /* { */ /* PMat(b_fcus->l*tree->mod->rr[i], */ /* tree->mod, */ /* &b_fcus->Pij_rr[i]); */ /* } */ /* Updating partial likelihood after branch swapping */ Update_P_Lk(tree,b_fcus,b_fcus->left); Update_P_Lk(tree,b_fcus,b_fcus->rght); tree->tot_loglk = Lk_At_Given_Edge(tree,b_fcus); return tree->tot_loglk; } /*********************************************************/ void Update_PMat_At_Given_Edge(edge *b_fcus, arbre *tree) { int i; double len; len = -1.0; For(i,tree->mod->n_catg) { len = b_fcus->l*tree->mod->rr[i]; if(len < BL_MIN) len = BL_MIN; PMat(len,tree->mod,&b_fcus->Pij_rr[i]); } } /*********************************************************/ allseq *Make_Seq(int n_otu, int len, char **sp_names) { allseq *alldata; int j; alldata = (allseq *)mCalloc(1,sizeof(allseq)); alldata->n_otu = n_otu; alldata->c_seq = (seq **)mCalloc(n_otu,sizeof(seq *)); alldata->wght = (double *)mCalloc(len,sizeof(double)); alldata->b_frq = (double *)mCalloc(T_MAX_ALPHABET,sizeof(double)); alldata->ambigu = (int *)mCalloc(len,sizeof(int)); alldata->sitepatt = (int *)mCalloc(len,sizeof(int )); alldata->invar = (int *)mCalloc(len,sizeof(int)); alldata->crunch_len = alldata->init_len = -1; For(j,n_otu) { alldata->c_seq[j] = (seq *)mCalloc(1,sizeof(seq)); alldata->c_seq[j]->name = (char *)mCalloc(T_MAX_NAME,sizeof(char)); strcpy(alldata->c_seq[j]->name,sp_names[j]); alldata->c_seq[j]->state = (char *)mCalloc(len+1,sizeof(char)); } return alldata; } /*********************************************************/ allseq *Copy_CData(allseq *ori, model *mod) { allseq *new; int i,j,n_otu; char **sp_names; n_otu = ori->n_otu; sp_names = (char **)mCalloc(n_otu,sizeof(char *)); For(i,n_otu) { sp_names[i] = (char *)mCalloc(T_MAX_NAME,sizeof(char)); strcpy(sp_names[i],ori->c_seq[i]->name); } new = Make_Seq(n_otu,ori->init_len,sp_names); For(i,n_otu) Free(sp_names[i]); Free(sp_names); For(i,ori->init_len) new->sitepatt[i] = ori->sitepatt[i]; For(j,ori->crunch_len) { For(i,ori->n_otu) new->c_seq[i]->state[j] = ori->c_seq[i]->state[j]; new->wght[j] = ori->wght[j]; new->ambigu[j] = ori->ambigu[j]; new->invar[j] = ori->invar[j]; } For(i,ori->n_otu) { new->c_seq[i]->len = ori->c_seq[i]->len; strcpy(new->c_seq[i]->name,ori->c_seq[i]->name); } new->init_len = ori->init_len; new->clean_len = ori->clean_len; new->crunch_len = ori->crunch_len; For(i,mod->ns) new->b_frq[i] = ori->b_frq[i]; new->n_otu = ori->n_otu; return new; } /*********************************************************/ optimiz *Alloc_Optimiz() { optimiz *s_opt; s_opt = (optimiz *)mCalloc(1,sizeof(optimiz)); return s_opt; } /*********************************************************/ void Init_Optimiz(optimiz *s_opt) { s_opt->print = 1; s_opt->last_opt = 1; s_opt->opt_alpha = 0; s_opt->opt_kappa = 0; s_opt->opt_bl = 0; s_opt->opt_pinvar = 0; s_opt->init_lk = UNLIKELY; s_opt->n_it_max = 1000; } /*********************************************************/ int Filexists(char *filename) { FILE *fp; fp =fopen(filename,"r"); if (fp) { fclose(fp); return 1; } else return 0; } /*********************************************************/ FILE *Openfile(char *filename, int mode) { /* mode = 0 -> read */ /* mode = 1 -> write */ /* mode = 2 -> append */ FILE *fp; char *s; int open_test=0; /* s = (char *)mCalloc(T_MAX_FILE,sizeof(char)); */ /* strcpy(s,filename); */ s = filename; fp = NULL; switch(mode) { case 0 : { while(!(fp = (FILE *)fopen(s,"r")) && ++open_test<10) { printf("\nCan't open file %s, enter a new name : ",s); Getstring_Stdin(s); fflush(stdout); } break; } case 1 : { fp = (FILE *)fopen(s,"w"); break; } case 2 : { fp = (FILE *)fopen(s,"a"); break; } default : break; } /* Free(s); */ return fp; } /*********************************************************/ void Print_Fp_Out(FILE *fp_out, time_t t_beg, time_t t_end, arbre *tree, option *input, int n_data_set) { (void)n_data_set; // dont warn about unused debug param char *s; div_t hour,min; fprintf(fp_out,". Sequence file : [%s]\n\n", input->seqfile); /* fprintf(fp_out,". Data set [#%d]\n",n_data_set); FLT*/ (tree->mod->datatype == NT)? (fprintf(fp_out,". Model of nucleotides substitution : %s\n\n",input->modelname)): (fprintf(fp_out,". Model of amino acids substitution : %s\n\n",input->modelname)); /*was after Sequence file ; moved here FLT*/ s = (char *)mCalloc(T_MAX_LINE,sizeof(char)); fprintf(fp_out,". Initial tree : [%s]\n\n", (!input->inputtree)?("BIONJ"): (strcat(strcat(strcat(s,"user tree ("),input->inputtreefile),")"))); Free(s); fprintf(fp_out,". Number of taxa : %d\n\n",tree->n_otu);/*added FLT*/ fprintf(fp_out,"\n"); fprintf(fp_out,". Likelihood : loglk = %.5f\n\n",tree->tot_loglk);/*was last ; moved here FLT*/ fprintf(fp_out,". Discrete gamma model : %s\n", (tree->mod->n_catg>1)?("Yes"):("No\n")); if(tree->mod->n_catg > 1) { fprintf(fp_out," - Number of categories : %d\n",tree->mod->n_catg); fprintf(fp_out," - Gamma shape parameter : %.3f\n\n",tree->mod->alpha); } if(tree->mod->invar) fprintf(fp_out,". Proportion of invariant : %.3f\n\n",tree->mod->pinvar); /*was before Discrete gamma model ; moved here FLT*/ if(tree->mod->whichmodel <= 5) { fprintf(fp_out,". Transition/transversion ratio : %.3f\n\n",tree->mod->kappa); } else if(tree->mod->whichmodel == 6) { fprintf(fp_out,". Transition/transversion ratio for purines : %.3f\n", tree->mod->kappa*2.*tree->mod->lambda/(1.+tree->mod->lambda)); fprintf(fp_out,". Transition/transversion ratio for pyrimidines : %.3f\n\n", tree->mod->kappa*2./(1.+tree->mod->lambda)); } if(tree->mod->datatype == NT) { fprintf(fp_out,". Nucleotides frequencies :\n\n"); fprintf(fp_out," - f(A)=%8.5f\n",tree->mod->pi[0]); fprintf(fp_out," - f(C)=%8.5f\n",tree->mod->pi[1]); fprintf(fp_out," - f(G)=%8.5f\n",tree->mod->pi[2]); fprintf(fp_out," - f(T)=%8.5f\n\n",tree->mod->pi[3]); } /*****************************************/ if((tree->mod->whichmodel == 7) || (tree->mod->whichmodel == 8)) { int i,j; printf("\n"); fprintf(fp_out,". GTR relative rate parameters : \n\n"); fprintf(fp_out,"A <-> C %8.5f\n",*(tree->mod->rr_param[0])); fprintf(fp_out,"A <-> G %8.5f\n",*(tree->mod->rr_param[1])); fprintf(fp_out,"A <-> T %8.5f\n",*(tree->mod->rr_param[2])); fprintf(fp_out,"C <-> G %8.5f\n",*(tree->mod->rr_param[3])); fprintf(fp_out,"C <-> T %8.5f\n",*(tree->mod->rr_param[4])); fprintf(fp_out,"G <-> T 1.0 (fixed)\n\n"); fprintf(fp_out,"\n. Instantaneous rate matrix : \n"); fprintf(fp_out,"\n[A---------C---------G---------T------]\n"); For(i,4) { For(j,4) fprintf(fp_out,"%8.5f ",tree->mod->mat_Q[i*4+j]); fprintf(fp_out,"\n"); } fprintf(fp_out,"\n"); fprintf(fp_out,"eg., the instantaneous rate of change from 'C' to 'A' is %8.5f x %8.5f = %8.5f\n\n", tree->mod->pi[0], *(tree->mod->rr_param[0]), tree->mod->mat_Q[1*4+0]); } /*****************************************/ hour = div(t_end-t_beg,3600); min = div(t_end-t_beg,60 ); min.quot -= hour.quot*60; fprintf(fp_out,". Time used %dh%dm%ds\n", hour.quot,min.quot,(int)(t_end-t_beg)%60); if(t_end-t_beg > 60) fprintf(fp_out,". -> %d seconds\n",(int)(t_end-t_beg)); fprintf(fp_out,"\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n\n"); fflush(fp_out); } /*********************************************************/ /*FLT wrote this function*/ void Print_Fp_Out_Lines(FILE *fp_out, time_t t_beg, time_t t_end, arbre *tree, option *input, int n_data_set) { (void)t_beg; // dont warn about unused debug params (void)t_end; char *s; /*div_t hour,min;*/ if (n_data_set==1) { fprintf(fp_out,". Sequence file : [%s]\n\n", input->seqfile); (tree->mod->datatype == NT)? (fprintf(fp_out,". Model of nucleotides substitution : %s\n\n",input->modelname)): (fprintf(fp_out,". Model of amino acids substitution : %s\n\n",input->modelname)); s = (char *)mCalloc(T_MAX_LINE,sizeof(char)); fprintf(fp_out,". Initial tree : [%s]\n\n", (!input->inputtree)?("BIONJ"): (strcat(strcat(strcat(s,"user tree ("),input->inputtreefile),")"))); Free(s); fprintf(fp_out,"\n"); /*headline 1*/ fprintf(fp_out, ". Data\t"); fprintf(fp_out,"Nb of \t"); fprintf(fp_out,"Likelihood\t"); fprintf(fp_out, "Discrete \t"); if(tree->mod->n_catg > 1) fprintf(fp_out, "Number of \tGamma shape\t"); fprintf(fp_out,"Proportion of\t"); if(tree->mod->whichmodel <= 6) fprintf(fp_out,"Transition/ \t"); fprintf(fp_out,"Nucleotides frequencies \t"); if(tree->mod->whichmodel == 7) fprintf(fp_out,"Instantaneous rate matrix \t"); /* fprintf(fp_out,"Time\t");*/ fprintf(fp_out, "\n"); /*headline 2*/ fprintf(fp_out, " set\t"); fprintf(fp_out,"taxa\t"); fprintf(fp_out,"loglk \t"); fprintf(fp_out, "gamma model\t"); if(tree->mod->n_catg > 1) fprintf(fp_out, "categories\tparameter \t"); fprintf(fp_out,"invariant \t"); if(tree->mod->whichmodel <= 6) fprintf(fp_out,"transversion\t"); fprintf(fp_out,"f(A) f(C) f(G) f(T) \t"); if(tree->mod->whichmodel == 7) fprintf(fp_out,"[A---------C---------G---------T------]\t"); /* fprintf(fp_out,"used\t");*/ fprintf(fp_out, "\n"); /*headline 3*/ if(tree->mod->whichmodel == 6) { fprintf(fp_out," \t \t \t \t"); if(tree->mod->n_catg > 1) fprintf(fp_out," \t \t"); fprintf(fp_out," \t"); fprintf(fp_out,"purines pyrimid.\t"); fprintf(fp_out, "\n"); } fprintf(fp_out, "\n"); } /*line items*/ fprintf(fp_out," #%d\t",n_data_set); fprintf(fp_out,"%d \t",tree->n_otu); fprintf(fp_out,"%.5f\t",tree->tot_loglk); fprintf(fp_out,"%s \t", (tree->mod->n_catg>1)?("Yes"):("No ")); if(tree->mod->n_catg > 1) { fprintf(fp_out,"%d \t",tree->mod->n_catg); fprintf(fp_out,"%.3f \t",tree->mod->alpha); } /*if(tree->mod->invar)*/ fprintf(fp_out,"%.3f \t",tree->mod->pinvar); if(tree->mod->whichmodel <= 5) { fprintf(fp_out,"%.3f \t",tree->mod->kappa); } else if(tree->mod->whichmodel == 6) { fprintf(fp_out,"%.3f ", tree->mod->kappa*2.*tree->mod->lambda/(1.+tree->mod->lambda)); fprintf(fp_out,"%.3f\t", tree->mod->kappa*2./(1.+tree->mod->lambda)); } if(tree->mod->datatype == NT) { fprintf(fp_out,"%8.5f ",tree->mod->pi[0]); fprintf(fp_out,"%8.5f ",tree->mod->pi[1]); fprintf(fp_out,"%8.5f ",tree->mod->pi[2]); fprintf(fp_out,"%8.5f\t",tree->mod->pi[3]); } /* hour = div(t_end-t_beg,3600); min = div(t_end-t_beg,60 ); min.quot -= hour.quot*60; fprintf(fp_out,"%dh%dm%ds\t", hour.quot,min.quot,(int)(t_end-t_beg)%60); if(t_end-t_beg > 60) fprintf(fp_out,". -> %d seconds\t",(int)(t_end-t_beg)); */ /*****************************************/ if((tree->mod->whichmodel == 7) || (tree->mod->whichmodel == 8)) { int i,j; For(i,4) { if (i!=0) { /*format*/ fprintf(fp_out," \t \t \t \t"); if(tree->mod->n_catg > 1) fprintf(fp_out," \t \t"); fprintf(fp_out," \t \t"); } For(j,4) fprintf(fp_out,"%8.5f ",tree->mod->mat_Q[i*4+j]); if (i<3) fprintf(fp_out,"\n"); } } /*****************************************/ fprintf(fp_out, "\n\n"); fflush(fp_out); } /*********************************************************/ void Alloc_All_P_Lk(arbre *tree) { int i,j,k; int nbytes; nbytes = 0; For(i,2*tree->n_otu-3) { tree->t_edges[i]->get_p_lk_left = 1; tree->t_edges[i]->get_p_lk_rght = 1; tree->t_edges[i]->p_lk_left = (double ***)mCalloc(tree->data->crunch_len,sizeof(double **)); nbytes += tree->data->crunch_len * sizeof(double **); tree->t_edges[i]->p_lk_rght = (double ***)mCalloc(tree->data->crunch_len,sizeof(double **)); nbytes += tree->data->clean_len * sizeof(double **); For(j,tree->data->crunch_len) { tree->t_edges[i]->p_lk_left[j] = (double **)mCalloc(tree->mod->n_catg,sizeof(double *)); nbytes += tree->mod->n_catg * sizeof(double *); For(k,tree->mod->n_catg) { tree->t_edges[i]->p_lk_left[j][k] = (double *)mCalloc(tree->mod->ns,sizeof(double )); nbytes += tree->mod->ns * sizeof(double); } tree->t_edges[i]->p_lk_rght[j] = (double **)mCalloc(tree->mod->n_catg,sizeof(double *)); nbytes += tree->mod->n_catg * sizeof(double *); For(k,tree->mod->n_catg) { tree->t_edges[i]->p_lk_rght[j][k] = (double *)mCalloc(tree->mod->ns,sizeof(double )); nbytes += tree->mod->ns * sizeof(double); } } tree->t_edges[i]->sum_scale_f_left = (double *)mCalloc(tree->data->crunch_len,sizeof(double )); tree->t_edges[i]->sum_scale_f_rght = (double *)mCalloc(tree->data->crunch_len,sizeof(double )); } printf("NBYTES = %d\n",nbytes); } /*********************************************************/ matrix *K2P_dist(allseq *data, double g_shape) { int i,j,k; int diff; matrix *mat; double **len; len = (double **)mCalloc(data->n_otu,sizeof(double *)); For(i,data->n_otu) len[i] = (double *)mCalloc(data->n_otu,sizeof(double)); mat = Make_Mat(data->n_otu); Init_Mat(mat,data); For(i,data->c_seq[0]->len) { For(j,data->n_otu-1) { for(k=j+1;kn_otu;k++) { if(((data->c_seq[j]->state[i] == 'A' || data->c_seq[j]->state[i] == 'G') && (data->c_seq[k]->state[i] == 'C' || data->c_seq[k]->state[i] == 'T'))|| ((data->c_seq[j]->state[i] == 'C' || data->c_seq[j]->state[i] == 'T') && (data->c_seq[k]->state[i] == 'A' || data->c_seq[k]->state[i] == 'G'))) { diff++; mat->Q[j][k]+=data->wght[i]; len[j][k]+=data->wght[i]; len[k][j]=len[j][k]; } else if(((data->c_seq[j]->state[i] == 'A' && data->c_seq[k]->state[i] == 'G') || (data->c_seq[j]->state[i] == 'G' && data->c_seq[k]->state[i] == 'A'))|| ((data->c_seq[j]->state[i] == 'C' && data->c_seq[k]->state[i] == 'T') || (data->c_seq[j]->state[i] == 'T' && data->c_seq[k]->state[i] == 'C'))) { diff++; mat->P[j][k]+=data->wght[i]; len[j][k]+=data->wght[i]; len[k][j]=len[j][k]; } else if((data->c_seq[j]->state[i] == 'A' || data->c_seq[j]->state[i] == 'C' || data->c_seq[j]->state[i] == 'G' || data->c_seq[j]->state[i] == 'T')&& (data->c_seq[k]->state[i] == 'A' || data->c_seq[k]->state[i] == 'C' || data->c_seq[k]->state[i] == 'G' || data->c_seq[k]->state[i] == 'T')) { len[j][k]+=data->wght[i]; len[k][j]=len[j][k]; } } } } For(i,data->n_otu-1) for(j=i+1;jn_otu;j++) { if(len[i][j]) { mat->P[i][j] /= len[i][j]; mat->Q[i][j] /= len[i][j]; } else { mat->P[i][j] = .5; mat->Q[i][j] = .5; } mat->P[j][i] = mat->P[i][j]; mat->Q[j][i] = mat->Q[i][j]; if((1-2*mat->P[i][j]-mat->Q[i][j] <= .0) || (1-2*mat->Q[i][j] <= .0)) { mat->dist[i][j] = -1.; mat->dist[j][i] = -1.; continue; } mat->dist[i][j] = (g_shape/2)* (pow(1-2*mat->P[i][j]-mat->Q[i][j],-1./g_shape) + 0.5*pow(1-2*mat->Q[i][j],-1./g_shape) - 1.5); if(mat->dist[i][j] > DIST_MAX) { mat->dist[i][j] = DIST_MAX; } mat->dist[j][i] = mat->dist[i][j]; } For(i,data->n_otu) free(len[i]); free(len); return mat; } /*********************************************************/ matrix *JC69_Dist(allseq *data, model *mod) { int site,i,j,k; matrix *mat; double **len; len = (double **)mCalloc(data->n_otu,sizeof(double *)); For(i,data->n_otu) len[i] = (double *)mCalloc(data->n_otu,sizeof(double)); mat = Make_Mat(data->n_otu); Init_Mat(mat,data); Fors(site,data->c_seq[0]->len,mod->stepsize) { For(j,data->n_otu-1) { for(k=j+1;kn_otu;k++) { if((!Is_Ambigu(data->c_seq[j]->state+site,mod->datatype,mod->stepsize)) && (!Is_Ambigu(data->c_seq[k]->state+site,mod->datatype,mod->stepsize))) { len[j][k]+=data->wght[site]; len[k][j]=len[j][k]; if(/*!Compare_Two_States*/strncmp(data->c_seq[j]->state+site, data->c_seq[k]->state+site, mod->stepsize)) mat->P[j][k]+=data->wght[site]; } } } } For(i,data->n_otu-1) for(j=i+1;jn_otu;j++) { if(len[i][j]) { mat->P[i][j] /= len[i][j]; } else { mat->P[i][j] = 1.; } mat->P[j][i] = mat->P[i][j]; if((1.-(mod->ns)/(mod->ns-1.)*mat->P[i][j]) < .0) { mat->dist[i][j] = DIST_MAX; } else mat->dist[i][j] = -(mod->ns-1.)/(mod->ns)*log(1.-(mod->ns)/(mod->ns-1.)*mat->P[i][j]); if(mat->dist[i][j] > DIST_MAX) { mat->dist[i][j] = DIST_MAX; } mat->dist[j][i] = mat->dist[i][j]; } For(i,data->n_otu) free(len[i]); free(len); return mat; } /*********************************************************/ matrix *Hamming_Dist(allseq *data, model *mod) { int i,j,k; matrix *mat; double **len; len = (double **)mCalloc(data->n_otu,sizeof(double *)); For(i,data->n_otu) len[i] = (double *)mCalloc(data->n_otu,sizeof(double)); mat = Make_Mat(data->n_otu); Init_Mat(mat,data); For(i,data->c_seq[0]->len) { For(j,data->n_otu-1) { for(k=j+1;kn_otu;k++) { if((!Is_Ambigu(data->c_seq[j]->state+i,mod->datatype,mod->stepsize)) && (!Is_Ambigu(data->c_seq[k]->state+i,mod->datatype,mod->stepsize))) { len[j][k]+=data->wght[i]; len[k][j]=len[j][k]; if(data->c_seq[j]->state[i] != data->c_seq[k]->state[i]) mat->P[j][k]+=data->wght[i]; } } } } For(i,data->n_otu-1) for(j=i+1;jn_otu;j++) { if(len[i][j]) { mat->P[i][j] /= len[i][j]; } else { mat->P[i][j] = 1.; } mat->P[j][i] = mat->P[i][j]; mat->dist[i][j] = mat->P[i][j]; if(mat->dist[i][j] > DIST_MAX) { mat->dist[i][j] = DIST_MAX; } mat->dist[j][i] = mat->dist[i][j]; } For(i,data->n_otu) free(len[i]); free(len); return mat; } /*********************************************************/ int Is_Ambigu(char *state, int datatype, int stepsize) { int i; if(datatype == NT) { For(i,stepsize) { if(strchr("MRWSYKBDHVNXO?-.",state[i])) return 1; } } else { if(strchr("X?-.",state[0])) return 1; } return 0; } /*********************************************************/ void Check_Ambiguities(allseq *data, int datatype, int stepsize) { int i,j; Fors(j,data->crunch_len,stepsize) For(i,data->n_otu) { if(Is_Ambigu(data->c_seq[i]->state+j, datatype, stepsize)) { data->ambigu[j] = 1; break; } } } /*********************************************************/ int Assign_State(char *c, int datatype, int stepsize) { int state[3]; int i; state[0] = -1; if(datatype == NT) { For(i,stepsize) { switch(c[i]) { case 'A' : state[i]=0; break; case 'C' : state[i]=1; break; case 'G' : state[i]=2; break; case 'T' : state[i]=3; break; case 'U' : state[i]=3; break; default : state[i]=-1; break; } } return (stepsize>1)?(state[0]*16+state[1]*4+state[2]):(state[0]); } else { switch(c[0]){ case 'A' : state[0]=0; break; case 'R' : state[0]=1; break; case 'N' : state[0]=2; break; case 'D' : state[0]=3; break; case 'C' : state[0]=4; break; case 'Q' : state[0]=5; break; case 'E' : state[0]=6; break; case 'G' : state[0]=7; break; case 'H' : state[0]=8; break; case 'I' : state[0]=9; break; case 'L' : state[0]=10; break; case 'K' : state[0]=11; break; case 'M' : state[0]=12; break; case 'F' : state[0]=13; break; case 'P' : state[0]=14; break; case 'S' : state[0]=15; break; case 'T' : state[0]=16; break; case 'W' : state[0]=17; break; case 'Y' : state[0]=18; break; case 'V' : state[0]=19; break; case 'B' : state[0] = 2; break; case 'Z' : state[0] = 5; break; default : state[0]=-1; break; } return state[0]; } return -1; } /*********************************************************/ void Bootstrap(arbre *tree) { int *site_num, n_site; int replicate,j,k; int position,init_len; double buff; allseq *boot_data; arbre *boot_tree; model *boot_mod; matrix *boot_mat; char *s; /* double rf; */ tree->mod->s_opt->last_opt = 0; tree->print_boot_val = 1; Alloc_Bip(tree); Get_Bip(tree->noeud[0], tree->noeud[0]->v[0], tree); site_num = (int *)mCalloc(tree->data->init_len,sizeof(int)); n_site = 0; For(j,tree->data->crunch_len) For(k,tree->data->wght[j]) { site_num[n_site] = j; n_site++; } boot_data = Copy_CData(tree->data,tree->mod); boot_tree = NULL; printf("\n\n. Non parametric bootstrap analysis \n\n"); printf(" ["); fflush(NULL); For(replicate,tree->mod->bootstrap) { For(j,boot_data->crunch_len) boot_data->wght[j] = .0; init_len = 0; For(j,boot_data->init_len) { buff = rand(); buff /= (RAND_MAX+1.); buff *= tree->data->init_len; position = (int)floor(buff); boot_data->wght[site_num[position]] += 1.; init_len++; } if(init_len != tree->data->init_len) Exit("\n. Pb in copying sequences\n"); (!tree->mod->datatype)? (Get_Base_Freqs(boot_data)): (Get_AA_Freqs(boot_data)); boot_mod = Copy_Model(tree->mod); Init_Model(boot_data,boot_mod); if(tree->input->inputtree) { rewind(tree->input->fp_input_tree); boot_tree = Read_Tree_File(tree->input->fp_input_tree); } else { boot_mat = ML_Dist(boot_data,boot_mod); boot_mat->tree = Make_Tree(boot_data); Bionj(boot_mat); boot_tree = boot_mat->tree; Free_Mat(boot_mat); } boot_tree->mod = boot_mod; boot_tree->input = tree->input; boot_tree->data = boot_data; boot_tree->both_sides = 1; boot_tree->n_pattern = boot_tree->data->crunch_len/ boot_tree->mod->stepsize; boot_tree->mod->s_opt->print = 0; Order_Tree_CSeq(boot_tree,boot_data); Share_Lk_Struct(tree,boot_tree); Init_P_Lk_Tips(boot_tree); if(boot_tree->mod->s_opt->opt_topo) Simu(boot_tree,1000); else { if(boot_tree->mod->s_opt->opt_free_param) Round_Optimize(boot_tree,boot_tree->data); else Lk(boot_tree,boot_data); } Alloc_Bip(boot_tree); Get_Bip(boot_tree->noeud[0], boot_tree->noeud[0]->v[0], boot_tree); Compare_Bip(tree,boot_tree); if(tree->input->print_boot_trees) { s = Write_Tree(boot_tree); fprintf(tree->input->fp_boot_tree,"%s\n",s); Free(s); Print_Fp_Out_Lines(tree->input->fp_boot_stats,0,0,boot_tree,tree->input,replicate+1); } /* rf = .0; */ /* For(j,2*tree->n_otu-3) */ /* rf += tree->t_edges[j]->bip_score; */ printf("."); if(!((replicate+1)%10)) { printf("] %d/%d\n ",replicate+1,tree->mod->bootstrap); if(replicate != tree->mod->bootstrap-1) printf("["); } fflush(NULL); Free_Tree(boot_tree); Free_Model(boot_mod); } if(((replicate)%10)) printf("] %d/%d\n ",replicate,tree->mod->bootstrap); if(tree->input->print_boot_trees) { fclose(tree->input->fp_boot_tree); fclose(tree->input->fp_boot_stats); } Free_Cseq(boot_data); Free(site_num); } /*********************************************************/ void Update_BrLen_Invar(arbre *tree) { int i; For(i,2*tree->n_otu-3) tree->t_edges[i]->l*=(1.0-tree->mod->pinvar); } /*********************************************************/ void Getstring_Stdin(char *file_name) { fgets(file_name,T_MAX_LINE,stdin); if (strchr(file_name, '\n') != NULL) *strchr(file_name, '\n') = '\0'; } /*********************************************************/ void Print_Freq(arbre *tree) { switch(tree->mod->datatype) { case NT: { printf("A : %f\n",tree->mod->pi[0]); printf("C : %f\n",tree->mod->pi[1]); printf("G : %f\n",tree->mod->pi[2]); printf("T : %f\n",tree->mod->pi[3]); printf("U : %f\n",tree->mod->pi[4]); printf("M : %f\n",tree->mod->pi[5]); printf("R : %f\n",tree->mod->pi[6]); printf("W : %f\n",tree->mod->pi[7]); printf("S : %f\n",tree->mod->pi[8]); printf("Y : %f\n",tree->mod->pi[9]); printf("K : %f\n",tree->mod->pi[10]); printf("B : %f\n",tree->mod->pi[11]); printf("D : %f\n",tree->mod->pi[12]); printf("H : %f\n",tree->mod->pi[13]); printf("V : %f\n",tree->mod->pi[14]); printf("N : %f\n",tree->mod->pi[15]); break; } case AA: { printf("A : %f\n",tree->mod->pi[0]); printf("R : %f\n",tree->mod->pi[1]); printf("N : %f\n",tree->mod->pi[2]); printf("D : %f\n",tree->mod->pi[3]); printf("C : %f\n",tree->mod->pi[4]); printf("Q : %f\n",tree->mod->pi[5]); printf("E : %f\n",tree->mod->pi[6]); printf("G : %f\n",tree->mod->pi[7]); printf("H : %f\n",tree->mod->pi[8]); printf("I : %f\n",tree->mod->pi[9]); printf("L : %f\n",tree->mod->pi[10]); printf("K : %f\n",tree->mod->pi[11]); printf("M : %f\n",tree->mod->pi[12]); printf("F : %f\n",tree->mod->pi[13]); printf("P : %f\n",tree->mod->pi[14]); printf("S : %f\n",tree->mod->pi[15]); printf("T : %f\n",tree->mod->pi[16]); printf("W : %f\n",tree->mod->pi[17]); printf("Y : %f\n",tree->mod->pi[18]); printf("V : %f\n",tree->mod->pi[19]); printf("N : %f\n",tree->mod->pi[20]); break; } default : {break;} } } /*********************************************************/ double Num_Derivatives_One_Param(double (*func)(arbre *tree), arbre *tree, double f0, double *param, double stepsize, double *err, int precise) { int i,j; double errt,fac,hh,**a,ans; int n_iter; a = (double **)mCalloc(11,sizeof(double *)); For(i,11) a[i] = (double *)mCalloc(11,sizeof(double)); n_iter = 10; /* */ ans = .0; if (stepsize == 0.0) Exit("\n. h must be nonzero in Dfridr."); hh=stepsize; if(!precise) { *param = *param+hh; a[0][0] = (*func)(tree); a[0][0] -= f0; a[0][0] /= hh; *param = *param-hh; ans = a[0][0]; } else { *param = *param+hh; a[0][0] = (*func)(tree); /* *param = *param-2*hh; */ /* a[0][0] -= (*func)(tree); */ /* a[0][0] /= (2.0*hh); */ /* *param = *param+hh; */ a[0][0] -= f0; a[0][0] /= hh; *param = *param-hh; *err=1e30; for(i=1;i= 2.0*(*err)) break; } } For(i,11) Free(a[i]); Free(a); return ans; } /*********************************************************/ void Num_Derivative_Several_Param(arbre *tree, double *param, int n_param, double stepsize, double (*func)(arbre *tree), double *derivatives) { int i; double err,f0; f0 = (*func)(tree); For(i,n_param) { derivatives[i] = Num_Derivatives_One_Param(func, tree, f0, param+i, stepsize, &err, 0 ); } } /*********************************************************/ int Compare_Two_States(char *state1, char *state2, int state_size) { /* 1 the two states are identical */ /* 0 the two states are different */ int i; For(i,state_size) if(state1[i] != state2[i]) break; return (i==state_size)?(1):(0); } /*********************************************************/ void Copy_One_State(char *from, char *to, int state_size) { int i; For(i,state_size) to[i] = from[i]; } /*********************************************************/ model *Make_Model_Basic() { model *mod; int i; mod = (model *)mCalloc(1,sizeof(model)); mod->custom_mod_string = (char *)mCalloc(T_MAX_LINE,sizeof(char)); mod->user_b_freq = (double *)mCalloc(100,sizeof(double)); mod->rr_param = (double **)mCalloc(6,sizeof(double *)); mod->rr_param_values = (double *)mCalloc(6,sizeof(double)); mod->rr_param_num = (int **)mCalloc(6,sizeof(int *)); mod->n_rr_param_per_cat = (int *)mCalloc(6,sizeof(int)); mod->s_opt = (optimiz *)Alloc_Optimiz(); For(i,6) mod->rr_param_num[i] = (int *)mCalloc(6,sizeof(int)); return mod; } /*********************************************************/ void Make_Model_Complete(model *mod) { int i,j; mod->pi = (double *)mCalloc(mod->ns,sizeof(double)); mod->r_proba = (double *)mCalloc(mod->n_catg,sizeof(double)); mod->rr = (double *)mCalloc(mod->n_catg,sizeof(double)); mod->Pij_rr = (double ***)mCalloc(mod->n_catg,sizeof(double **)); mod->dPij_rr = (double ***)mCalloc(mod->n_catg,sizeof(double **)); mod->d2Pij_rr = (double ***)mCalloc(mod->n_catg,sizeof(double **)); For(i,mod->n_catg) { mod->Pij_rr[i] = (double **)mCalloc(mod->ns,sizeof(double *)); mod->dPij_rr[i] = (double **)mCalloc(mod->ns,sizeof(double *)); mod->d2Pij_rr[i] = (double **)mCalloc(mod->ns,sizeof(double *)); For(j,mod->ns) { mod->Pij_rr[i][j] = (double *)mCalloc(mod->ns,sizeof(double)); mod->dPij_rr[i][j] = (double *)mCalloc(mod->ns,sizeof(double)); mod->d2Pij_rr[i][j] = (double *)mCalloc(mod->ns,sizeof(double)); } } mod->mat_Q = (double *)mCalloc(mod->ns*mod->ns,sizeof(double)); mod->mat_Vr = (double *)mCalloc(mod->ns*mod->ns,sizeof(double)); mod->mat_Vi = (double *)mCalloc(mod->ns*mod->ns,sizeof(double)); mod->vct_eDmr = (double *)mCalloc(mod->ns ,sizeof(double)); mod->vct_ev = (double *)mCalloc(mod->ns ,sizeof(double)); } /*********************************************************/ model *Copy_Model(model *ori) { int i,j; model *cpy; cpy = Make_Model_Basic(); Copy_Optimiz(ori->s_opt,cpy->s_opt); /* cpy->c_code = ori->c_code; */ cpy->ns = ori->ns; cpy->n_catg = ori->n_catg; Make_Model_Complete(cpy); cpy->datatype = ori->datatype; cpy->n_otu = ori->n_otu; cpy->alpha_old = ori->alpha_old; cpy->kappa_old = ori->alpha_old; cpy->lambda_old = ori->lambda_old; cpy->pinvar_old = ori->pinvar_old; cpy->whichmodel = ori->whichmodel; cpy->seq_len = ori->seq_len; cpy->update_eigen = ori->update_eigen; cpy->kappa = ori->kappa; cpy->alpha = ori->alpha; cpy->lambda = ori->lambda; cpy->bootstrap = ori->bootstrap; cpy->invar = ori->invar; cpy->pinvar = ori->pinvar; cpy->stepsize = ori->stepsize; cpy->n_diff_rr_param = ori->n_diff_rr_param; For(i,cpy->n_diff_rr_param) { cpy->n_rr_param_per_cat[i] = ori->n_rr_param_per_cat[i]; For(j,cpy->n_rr_param_per_cat[i]) { cpy->rr_param_num[i][j] = ori->rr_param_num[i][j]; } } For(i,6) { cpy->rr_param_values[i] = ori->rr_param_values[i]; cpy->rr_param[i] = cpy->rr_param_values+i; } For(i,cpy->ns) { cpy->pi[i] = ori->pi[i]; cpy->user_b_freq[i] = ori->user_b_freq[i]; } For(i,cpy->n_catg) { cpy->r_proba[i] = ori->r_proba[i]; cpy->rr[i] = ori->rr[i]; } return cpy; } /*********************************************************/ void Set_Defaults_Input(option* input) { input->mod->datatype = 0; strcpy(input->modelname,"HKY"); strcpy(input->nt_or_cd,"nucleotides"); input->n_data_sets = 1; input->interleaved = 1; input->inputtree = 0; input->tree = NULL; input->phyml_tree_file_open_mode = 1; input->phyml_stat_file_open_mode = 1; input->seq_len = -1; input->n_data_set_asked = -1; input->print_boot_trees = 1; } /*********************************************************/ void Set_Defaults_Model(model *mod) { int i; strcpy(mod->custom_mod_string,"000000"); mod->whichmodel = 4; mod->n_catg = 1; mod->kappa = 4.0; mod->alpha = 2.0; mod->lambda = 1.0; mod->bootstrap = 0; mod->invar = 0; mod->pinvar = 0.0; mod->stepsize = 1; mod->ns = 4; mod->n_diff_rr_param = 0; For(i,6) { mod->rr_param_values[i] = 1.0; mod->rr_param[i] = mod->rr_param_values+i; } For(i,4) mod->user_b_freq[i] = -1.; } /*********************************************************/ void Set_Defaults_Optimiz(optimiz *s_opt) { s_opt->opt_alpha = 0; s_opt->opt_kappa = 0; s_opt->opt_lambda = 0; s_opt->opt_bl = 0; s_opt->opt_pinvar = 0; s_opt->opt_rr_param = 0; s_opt->opt_topo = 1; s_opt->opt_free_param = 1; } /*********************************************************/ void Copy_Optimiz(optimiz *ori, optimiz *cpy) { cpy->print = ori->print; cpy->opt_alpha = ori->opt_alpha; cpy->opt_kappa = ori->opt_kappa; cpy->opt_lambda = ori->opt_lambda; cpy->opt_pinvar = ori->opt_pinvar; cpy->opt_rr_param = ori->opt_rr_param; cpy->opt_free_param = ori->opt_free_param; cpy->opt_bl = ori->opt_bl; cpy->init_lk = ori->init_lk; cpy->n_it_max = ori->n_it_max; cpy->opt_topo = ori->opt_topo; } /*********************************************************/ void Get_Bip(node *a, node *d, arbre *tree) { if(d->tax) { d->bip_node[0][0] = d; d->bip_size[0] = 1; return; } else { int i,j,k; int d_a; d_a = -1; For(i,3) { if(d->v[i] != a) Get_Bip(d,d->v[i],tree); else d_a = i; } d->bip_size[d_a] = 0; For(i,3) if(d->v[i] !=a ) { For(j,3) { if(d->v[i]->v[j] == d) { For(k,d->v[i]->bip_size[j]) { d->bip_node[d_a][d->bip_size[d_a]] = d->v[i]->bip_node[j][k]; strcpy(d->bip_name[d_a][d->bip_size[d_a]],d->v[i]->bip_node[j][k]->name); d->bip_size[d_a]++; } break; } } } qsort(d->bip_name[d_a],d->bip_size[d_a],sizeof(char *),Sort_String); For(i,3) if(a->v[i] == d) { a->bip_size[i] = 0; For(j,tree->n_otu) { For(k,d->bip_size[d_a]) { if(d->bip_node[d_a][k] == tree->noeud[j]) break; } if(k == d->bip_size[d_a]) { a->bip_node[i][a->bip_size[i]] = tree->noeud[j]; strcpy(a->bip_name[i][a->bip_size[i]],tree->noeud[j]->name); a->bip_size[i]++; } } qsort(a->bip_name[i],a->bip_size[i],sizeof(char *),Sort_String); if(a->bip_size[i] != tree->n_otu - d->bip_size[d_a]) { printf("%d %d \n",a->bip_size[i],tree->n_otu - d->bip_size[d_a]); Exit("\n. Problem in counting bipartitions \n"); } break; } } } /*********************************************************/ void Alloc_Bip(arbre *tree) { int i,j,k; tree->has_bip = 1; For(i,2*tree->n_otu-2) { tree->noeud[i]->bip_size = (int *)mCalloc(3,sizeof(int)); tree->noeud[i]->bip_node = (node ***)mCalloc(3,sizeof(node **)); tree->noeud[i]->bip_name = (char ***)mCalloc(3,sizeof(char **)); For(j,3) { tree->noeud[i]->bip_node[j] = (node **)mCalloc(tree->n_otu,sizeof(node *)); tree->noeud[i]->bip_name[j] = (char **)mCalloc(tree->n_otu,sizeof(char *)); For(k,tree->n_otu) tree->noeud[i]->bip_name[j][k] = (char *)mCalloc(T_MAX_NAME,sizeof(char )); } } } /*********************************************************/ int Sort_Double_Increase(const void *a, const void *b) { if((*(double *)(a)) <= (*(double *)(b))) return -1; else return 1; } /*********************************************************/ int Sort_String(const void *a, const void *b) { return(strcmp((*(const char **)(a)), (*(const char **)(b)))); } /*********************************************************/ void Compare_Bip(arbre *tree1, arbre *tree2) { int i,j,k; edge *b1,*b2; char **bip1,**bip2; int bip_size; For(i,2*tree1->n_otu-3) { if((!tree1->t_edges[i]->left->tax) && (!tree1->t_edges[i]->rght->tax)) { b1 = tree1->t_edges[i]; For(j,2*tree2->n_otu-3) { if((!tree2->t_edges[j]->left->tax) && (!tree2->t_edges[j]->rght->tax)) { b2 = tree2->t_edges[j]; if(MIN(b1->left->bip_size[b1->l_r],b1->rght->bip_size[b1->r_l]) == MIN(b2->left->bip_size[b2->l_r],b2->rght->bip_size[b2->r_l])) { bip_size = MIN(b1->left->bip_size[b1->l_r],b1->rght->bip_size[b1->r_l]); if(b1->left->bip_size[b1->l_r] == b1->rght->bip_size[b1->r_l]) { if(b1->left->bip_name[b1->l_r][0][0] < b1->rght->bip_name[b1->r_l][0][0]) { bip1 = b1->left->bip_name[b1->l_r]; } else { bip1 = b1->rght->bip_name[b1->r_l]; } } else if(b1->left->bip_size[b1->l_r] < b1->rght->bip_size[b1->r_l]) { bip1 = b1->left->bip_name[b1->l_r]; } else { bip1 = b1->rght->bip_name[b1->r_l]; } if(b2->left->bip_size[b2->l_r] == b2->rght->bip_size[b2->r_l]) { if(b2->left->bip_name[b2->l_r][0][0] < b2->rght->bip_name[b2->r_l][0][0]) { bip2 = b2->left->bip_name[b2->l_r]; } else { bip2 = b2->rght->bip_name[b2->r_l]; } } else if(b2->left->bip_size[b2->l_r] < b2->rght->bip_size[b2->r_l]) { bip2 = b2->left->bip_name[b2->l_r]; } else { bip2 = b2->rght->bip_name[b2->r_l]; } if(bip_size == 1) Exit("\n. Problem in Compare_Bip\n"); For(k,bip_size) { if(strcmp(bip1[k],bip2[k])) break; } if(k == bip_size) { b1->bip_score++; b2->bip_score++; break; } } } } } } } /*********************************************************/ void Test_Multiple_Data_Set_Format(option *input) { char *line; line = (char *)mCalloc(T_MAX_LINE,sizeof(char)); input->n_trees = 0; while(fgets(line,T_MAX_LINE,input->fp_input_tree)) if(strstr(line,";")) input->n_trees++; Free(line); /* if((input->n_trees != input->n_data_sets) && */ /* (input->n_data_sets > 1)) */ /* Exit("\n. The number of trees should be the same as\n the number of data sets\n\n"); */ if((input->mod->bootstrap > 1) && (input->n_trees > 1)) Exit("\n. Bootstrap option is not allowed with multiple trees\n"); rewind(input->fp_input_tree); return; } /*********************************************************/ int Are_Compatible(char *statea, char *stateb, int stepsize, int datatype) { int i,j; char a,b; if(datatype == NT) { For(i,stepsize) { a = statea[i]; For(j,stepsize) { b = stateb[j]; switch(a) { case 'A': { switch(b) { case 'A' : case 'M' : case 'R' : case 'W' : case 'D' : case 'H' : case 'V' : case 'X' : {b=b; break;} default : return 0; } break; } case 'G': { switch(b) { case 'G' : case 'R' : case 'S' : case 'K' : case 'B' : case 'D' : case 'V' : case 'X' : {b=b; break;} default : return 0; } break; } case 'C': { switch(b) { case 'C' : case 'M' : case 'S' : case 'Y' : case 'B' : case 'H' : case 'V' : case 'X' : {b=b; break;} default : return 0; } break; } case 'T': { switch(b) { case 'T' : case 'W' : case 'Y' : case 'K' : case 'B' : case 'D' : case 'H' : case 'X' : {b=b; break;} default : return 0; } break; } case 'M' : { switch(b) { case 'M' : case 'A' : case 'C' : case 'R' : case 'W' : case 'S' : case 'Y' : case 'B' : case 'D' : case 'H' : case 'V' : case 'X' : {b=b; break;} default : return 0; } break; } case 'R' : { switch(b) { case 'R' : case 'A' : case 'G' : case 'M' : case 'W' : case 'S' : case 'K' : case 'B' : case 'D' : case 'H' : case 'V' : case 'X' : {b=b; break;} default : return 0; } break; } case 'W' : { switch(b) { case 'W' : case 'A' : case 'T' : case 'M' : case 'R' : case 'Y' : case 'K' : case 'B' : case 'D' : case 'H' : case 'V' : case 'X' : {b=b; break;} default : return 0; } break; } case 'S' : { switch(b) { case 'S' : case 'C' : case 'G' : case 'M' : case 'R' : case 'Y' : case 'K' : case 'B' : case 'D' : case 'H' : case 'V' : case 'X' : {b=b; break;} default : return 0; } break; } case 'Y' : { switch(b) { case 'Y' : case 'C' : case 'T' : case 'M' : case 'W' : case 'S' : case 'K' : case 'B' : case 'D' : case 'H' : case 'V' : case 'X' : {b=b; break;} default : return 0; } break; } case 'K' : { switch(b) { case 'K' : case 'G' : case 'T' : case 'R' : case 'W' : case 'S' : case 'Y' : case 'B' : case 'D' : case 'H' : case 'V' : case 'X' : {b=b; break;} default : return 0; } break; } case 'B' : { switch(b) { case 'B' : case 'C' : case 'G' : case 'T' : case 'M' : case 'R' : case 'W' : case 'S' : case 'Y' : case 'K' : case 'D' : case 'H' : case 'V' : case 'X' : {b=b; break;} default : return 0; } break; } case 'D' : { switch(b) { case 'D' : case 'A' : case 'G' : case 'T' : case 'M' : case 'R' : case 'W' : case 'S' : case 'Y' : case 'K' : case 'B' : case 'H' : case 'V' : case 'X' : {b=b; break;} default : return 0; } break; } case 'H' : { switch(b) { case 'H' : case 'A' : case 'C' : case 'T' : case 'M' : case 'R' : case 'W' : case 'S' : case 'Y' : case 'K' : case 'B' : case 'D' : case 'V' : case 'X' : {b=b; break;} default : return 0; } break; } case 'V' : { switch(b) { case 'V' : case 'A' : case 'C' : case 'G' : case 'M' : case 'R' : case 'W' : case 'S' : case 'Y' : case 'K' : case 'B' : case 'D' : case 'H' : case 'X' : {b=b; break;} default : return 0; } break; } case 'X' : { switch(b) { case 'X' : case 'A' : case 'C' : case 'G' : case 'T' : case 'M' : case 'R' : case 'W' : case 'S' : case 'Y' : case 'K' : case 'B' : case 'D' : case 'H' : case 'V' : {b=b; break;} default : return 0; } break; } default : { printf("\n. Err. in Are_Compatible\n"); printf("\n. Please check that characters `%c` and `%c`\n",a,b); printf(" correspond to existing amino-acids.\n"); Exit("\n"); return 0; } } } } } else { a = statea[0]; b = stateb[0]; switch(a) { case 'A' : { switch(b) { case 'A' : case 'X' : {b=b; break;} default : return 0; } break; } case 'R' : { switch(b) { case 'R' : case 'X' : {b=b; break;} default : return 0; } break; } case 'N' : { switch(b) { case 'N' : case 'B' : case 'X' : {b=b; break;} default : return 0; } break; } case 'B' : { switch(b) { case 'N' : case 'B' : case 'X' : {b=b; break;} default : return 0; } break; } case 'D' : { switch(b) { case 'D' : case 'X' : {b=b; break;} default : return 0; } break; } case 'C' : { switch(b) { case 'C' : case 'X' : {b=b; break;} default : return 0; } break; } case 'Q' : { switch(b) { case 'Q' : case 'Z' : case 'X' : {b=b; break;} default : return 0; } break; } case 'Z' : { switch(b) { case 'Q' : case 'Z' : case 'X' : {b=b; break;} default : return 0; } break; } case 'E' : { switch(b) { case 'E' : case 'X' : {b=b; break;} default : return 0; } break; } case 'G' : { switch(b) { case 'G' : case 'X' : {b=b; break;} default : return 0; } break; } case 'H' : { switch(b) { case 'H' : case 'X' : {b=b; break;} default : return 0; } break; } case 'I' : { switch(b) { case 'I' : case 'X' : {b=b; break;} default : return 0; } break; } case 'L' : { switch(b) { case 'L' : case 'X' : {b=b; break;} default : return 0; } break; } case 'K' : { switch(b) { case 'K' : case 'X' : {b=b; break;} default : return 0; } break; } case 'M' : { switch(b) { case 'M' : case 'X' : {b=b; break;} default : return 0; } break; } case 'F' : { switch(b) { case 'F' : case 'X' : {b=b; break;} default : return 0; } break; } case 'P' : { switch(b) { case 'P' : case 'X' : {b=b; break;} default : return 0; } break; } case 'S' : { switch(b) { case 'S' : case 'X' : {b=b; break;} default : return 0; } break; } case 'T' : { switch(b) { case 'T' : case 'X' : {b=b; break;} default : return 0; } break; } case 'W' : { switch(b) { case 'W' : case 'X' : {b=b; break;} default : return 0; } break; } case 'Y' : { switch(b) { case 'Y' : case 'X' : {b=b; break;} default : return 0; } break; } case 'V' : { switch(b) { case 'V' : case 'X' : {b=b; break;} default : return 0; } break; } case 'X' : { switch(b) { case 'A':case 'R':case 'N' :case 'B' :case 'D' : case 'C':case 'Q':case 'Z' :case 'E' :case 'G' : case 'H':case 'I':case 'L' :case 'K' :case 'M' : case 'F':case 'P':case 'S' :case 'T' :case 'W' : case 'Y':case 'V': case 'X' : {b=b; break;} default : return 0; } break; } default : { printf("\n. Err. in Are_Compatible\n"); printf("\n. Please check that characters `%c` and `%c`\n",a,b); printf(" correspond to existing amino-acids.\n"); Exit("\n"); return 0; } } } return 1; } /*********************************************************/ void Hide_Ambiguities(allseq *data) { int i; For(i,data->crunch_len) { if(data->ambigu[i]) { data->wght[i] = 0.0; } } } /*********************************************************/ /*********************************************************/ /*********************************************************/ /*********************************************************/ /*********************************************************/