+++ /dev/null
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-#include <stdarg.h>
-
-#include "io_lib_header.h"
-#include "util_lib_header.h"
-#include "dp_lib_header.h"
-#include "define_header.h"
-
-static int first_seq, last_seq;
-
-static int nseqs;
-static char **names; /* the seq. names */
-
-static double **tmat;
-static double *av;
-static double *left_branch, *right_branch;
-static int *tkill;
-
-
-/*!
- * \file util_make_tree.c
- * \brief Source code tree algorithms
- */
-
-NT_node ** seq2cw_tree ( Sequence *S, char *tree)
-{
- Alignment *A;
- char *aln,command[1000];
- int tot_node;
-
-
-
- A=seq2clustalw_aln (S);
- aln=vtmpnam (NULL); if (!tree)tree=vtmpnam (NULL);
-
- output_fasta_aln (aln, A);
-
- sprintf ( command, "clustalw -infile=%s -newtree=%s -tree %s", aln, tree, TO_NULL_DEVICE);
- my_system (command);
- return read_tree (tree, &tot_node, S->nseq, S->name);
-}
-
-NT_node ** make_upgma_tree ( Alignment *A,int **distances,int gop, int gep, char **out_seq, char **out_seq_name, int out_nseq, char *tree_file, char *tree_mode)
- {
-
- if ( distances==NULL && A==NULL)
- {
- fprintf ( stderr, "\nError: make_nj_tree, must provide an alignment or a distance matrix [FATAL:%s]", PROGRAM);
- myexit (EXIT_FAILURE);
- }
- else if ( distances==NULL)
- {
-
- distances=get_dist_aln_array (A, "idmat");
- }
- return int_dist2upgma_tree (distances,A, A->nseq,tree_file);
- }
-NT_node ** make_nj_tree ( Alignment *A,int **distances,int gop, int gep, char **out_seq, char **out_seq_name, int out_nseq, char *tree_file, char *tree_mode)
- {
-
- if ( distances==NULL && A==NULL)
- {
- fprintf ( stderr, "\nError: make_nj_tree, must provide an alignment or a distance matrix [FATAL:%s]", PROGRAM);
- myexit (EXIT_FAILURE);
- }
- else if ( distances==NULL)
- {
- distances=get_dist_aln_array (A, "idmat");
- }
- return int_dist2nj_tree (distances,out_seq_name, out_nseq,tree_file);
- }
-
-NT_node ** int_dist2nj_tree (int **distances, char **out_seq_name, int out_nseq, char *tree_file)
- {
- int a, b;
- double **d;
- NT_node **T;
-
- d=declare_double( out_nseq, out_nseq);
- for ( a=0; a< out_nseq; a++)
- for ( b=0; b< out_nseq; b++)
- d[a][b]=distances[a][b];
-
- T=dist2nj_tree ( d, out_seq_name, out_nseq, tree_file);
-
- free_double (d, -1);
- return T;
- }
-NT_node ** float_dist2nj_tree (float **distances, char **out_seq_name, int out_nseq, char *tree_file)
- {
- int a, b;
- double **d;
- NT_node **T;
-
- d=declare_double( out_nseq, out_nseq);
- for ( a=0; a< out_nseq; a++)
- for ( b=0; b< out_nseq; b++)
- d[a][b]=distances[a][b];
- T=dist2nj_tree ( d, out_seq_name, out_nseq, tree_file);
- free_double (d, -1);
- return T;
- }
-NT_node ** dist2nj_tree (double **distances, char **out_seq_name, int out_nseq, char *tree_file)
- {
- int a, b;
- double **d_dist;
- int tot_node=0;
- NT_node **T;
-
- if ( !tree_file)tree_file=vtmpnam(NULL);
- d_dist=declare_double( out_nseq+1, out_nseq+1);
-
- for ( a=0; a<out_nseq; a++)
- for ( b=0; b< out_nseq; b++)
- {
- if ( a!=b)
- d_dist[a+1][b+1]=distances[a][b]/MAXID;
- else d_dist[a+1][b+1]=0;
- }
-
- guide_tree ( tree_file, d_dist, out_seq_name, out_nseq);
- free_double (d_dist,-1);
- T=read_tree (tree_file,&tot_node,out_nseq, out_seq_name);
-
- return T;
- }
-
-
-void guide_tree(char *fname, double **saga_tmat, char **saga_seq_name, int saga_nseq)
-/*
- Routine for producing unrooted NJ trees from seperately aligned
- pairwise distances. This produces the GUIDE DENDROGRAMS in
- PHYLIP format.
-*/
-{
-
- int i;
- FILE *fp;
- char **standard_tree;
-
-
-
- nseqs=saga_nseq;
- first_seq=1;
- last_seq=nseqs;
-
- names=saga_seq_name;
- tmat=saga_tmat;
-
- standard_tree = (char **) vmalloc( (nseqs+1) * sizeof (char *) );
- for(i=0; i<nseqs+1; i++)
- standard_tree[i] = (char *) vmalloc( (nseqs+1) * sizeof(char));
-
-
- nj_tree(standard_tree, nseqs);
-
- fp=fopen ( fname, "w");
- print_phylip_tree(standard_tree,fp,FALSE);
-
- vfree(left_branch);
- vfree(right_branch);
- vfree(tkill);
- vfree(av);
-
- for (i=0;i<nseqs+1;i++)
- vfree(standard_tree[i]);
-
-
- fclose(fp);
-
-
-}
-
-
-void nj_tree(char **tree_description, int nseq)
-{
- if ( nseq<100)
- slow_nj_tree (tree_description);
- else
- fast_nj_tree (tree_description);
-
-}
-
-
-
-void slow_nj_tree(char **tree_description)
-{
- int i;
- int l[4],nude,k;
- int nc,mini,minj,j,ii,jj;
- double fnseqs,fnseqs2=0,sumd;
- double diq,djq,dij,d2r,dr,dio,djo,da;
- double tmin,total,dmin;
- double bi,bj,b1,b2,b3,branch[4];
- int typei,typej; /* 0 = node; 1 = OTU */
-
- fnseqs = (double)nseqs;
-
-/*********************** First initialisation ***************************/
-
-
-
- mini = minj = 0;
-
- left_branch = (double *) vcalloc( (nseqs+2),sizeof (double) );
- right_branch = (double *) vcalloc( (nseqs+2),sizeof (double) );
- tkill = (int *) vcalloc( (nseqs+1),sizeof (int) );
- av = (double *) vcalloc( (nseqs+1),sizeof (double) );
-
- for(i=1;i<=nseqs;++i)
- {
- tmat[i][i] = av[i] = 0.0;
- tkill[i] = 0;
- }
-
-/*********************** Enter The Main Cycle ***************************/
-
- /**start main cycle**/
- for(nc=1; nc<=(nseqs-3); ++nc)
- {
-
- sumd = 0.0;
- for(j=2; j<=nseqs; ++j)
- for(i=1; i<j; ++i)
- {
- tmat[j][i] = tmat[i][j];
- sumd = sumd + tmat[i][j];
- }
- tmin = 99999.0;
-
-/*.................compute SMATij values and find the smallest one ........*/
-
- for(jj=2; jj<=nseqs; ++jj)
- if(tkill[jj] != 1)
- for(ii=1; ii<jj; ++ii)
- if(tkill[ii] != 1)
- {
- diq = djq = 0.0;
-
- for(i=1; i<=nseqs; ++i)
- {
- diq = diq + tmat[i][ii];
- djq = djq + tmat[i][jj];
- }
- dij = tmat[ii][jj];
- d2r = diq + djq - (2.0*dij);
- dr = sumd - dij -d2r;
- fnseqs2 = fnseqs - 2.0;
- total= d2r+ fnseqs2*dij +dr*2.0;
- total= total / (2.0*fnseqs2);
-
- /* commented out to have consistent results with CYGWIN: if(total < tmin)"*/
- if(total < tmin)
- {
- if ( tmin-total>MY_EPS)
- tmin = total;
- mini = ii;
- minj = jj;
- }
- }
-
-
-/*.................compute branch lengths and print the results ........*/
-
-
- dio = djo = 0.0;
- for(i=1; i<=nseqs; ++i) {
- dio = dio + tmat[i][mini];
- djo = djo + tmat[i][minj];
- }
-
- dmin = tmat[mini][minj];
- dio = (dio - dmin) / fnseqs2;
- djo = (djo - dmin) / fnseqs2;
- bi = (dmin + dio - djo) * 0.5;
- bj = dmin - bi;
- bi = bi - av[mini];
- bj = bj - av[minj];
-
- if( av[mini] > 0.0 )
- typei = 0;
- else
- typei = 1;
- if( av[minj] > 0.0 )
- typej = 0;
- else
- typej = 1;
-
-
-
-/*
- set negative branch lengths to zero. Also set any tiny positive
- branch lengths to zero.
-*/
- if( fabs(bi) < 0.0001) bi = 0.0;
- if( fabs(bj) < 0.0001) bj = 0.0;
-
-
-
-
- left_branch[nc] = bi;
- right_branch[nc] = bj;
-
- for(i=1; i<=nseqs; i++)
- tree_description[nc][i] = 0;
-
- if(typei == 0) {
- for(i=nc-1; i>=1; i--)
- if(tree_description[i][mini] == 1) {
- for(j=1; j<=nseqs; j++)
- if(tree_description[i][j] == 1)
- tree_description[nc][j] = 1;
- break;
- }
- }
- else
- tree_description[nc][mini] = 1;
-
- if(typej == 0) {
- for(i=nc-1; i>=1; i--)
- if(tree_description[i][minj] == 1) {
- for(j=1; j<=nseqs; j++)
- if(tree_description[i][j] == 1)
- tree_description[nc][j] = 1;
- break;
- }
- }
- else
- tree_description[nc][minj] = 1;
-
-
-/*
- Here is where the -0.00005 branch lengths come from for 3 or more
- identical seqs.
-*/
-/* if(dmin <= 0.0) dmin = 0.0001; */
- if(dmin <= 0.0) dmin = 0.000001;
- av[mini] = dmin * 0.5;
-
- /*........................Re-initialisation................................*/
-
- fnseqs = fnseqs - 1.0;
- tkill[minj] = 1;
-
- for(j=1; j<=nseqs; ++j)
- if( tkill[j] != 1 ) {
- da = ( tmat[mini][j] + tmat[minj][j] ) * 0.5;
- if( (mini - j) < 0 )
- tmat[mini][j] = da;
- if( (mini - j) > 0)
- tmat[j][mini] = da;
- }
-
- for(j=1; j<=nseqs; ++j)
- tmat[minj][j] = tmat[j][minj] = 0.0;
-
-
-
- }
- /*end main cycle**/
-
-/******************************Last Cycle (3 Seqs. left)********************/
-
- nude = 1;
-
-
- for(i=1; i<=nseqs; ++i)
- if( tkill[i] != 1 ) {
- l[nude] = i;
- nude = nude + 1;
- }
-
- b1 = (tmat[l[1]][l[2]] + tmat[l[1]][l[3]] - tmat[l[2]][l[3]]) * 0.5;
- b2 = tmat[l[1]][l[2]] - b1;
- b3 = tmat[l[1]][l[3]] - b1;
-
- branch[1] = b1 - av[l[1]];
- branch[2] = b2 - av[l[2]];
- branch[3] = b3 - av[l[3]];
-
-/* Reset tiny negative and positive branch lengths to zero */
- if( fabs(branch[1]) < 0.0001) branch[1] = 0.0;
- if( fabs(branch[2]) < 0.0001) branch[2] = 0.0;
- if( fabs(branch[3]) < 0.0001) branch[3] = 0.0;
-
- left_branch[nseqs-2] = branch[1];
- left_branch[nseqs-1] = branch[2];
- left_branch[nseqs] = branch[3];
-
- for(i=1; i<=nseqs; i++)
- tree_description[nseqs-2][i] = 0;
-
-
-
- for(i=1; i<=3; ++i) {
- if( av[l[i]] > 0.0) {
-
- for(k=nseqs-3; k>=1; k--)
- if(tree_description[k][l[i]] == 1) {
- for(j=1; j<=nseqs; j++)
- if(tree_description[k][j] == 1)
- tree_description[nseqs-2][j] = i;
- break;
- }
- }
- else {
-
- tree_description[nseqs-2][l[i]] = i;
- }
- if(i < 3) {
- }
- }
-}
-
-void print_phylip_tree(char **tree_description, FILE *tree, int bootstrap)
-{
-
- fprintf(tree,"(\n");
-
- two_way_split(tree_description, tree, nseqs-2,1,bootstrap);
- fprintf(tree,":%7.5f,\n",left_branch[nseqs-2]);
- two_way_split(tree_description, tree, nseqs-2,2,bootstrap);
- fprintf(tree,":%7.5f,\n",left_branch[nseqs-1]);
- two_way_split(tree_description, tree, nseqs-2,3,bootstrap);
-
- fprintf(tree,":%7.5f)",left_branch[nseqs]);
- if (bootstrap) fprintf(tree,"TRICHOTOMY");
- fprintf(tree,";\n");
-}
-
-
-void two_way_split
-(char **tree_description, FILE *tree, int start_row, int flag, int bootstrap)
-{
- int row, new_row, col, test_col=0;
- int single_seq;
-
- if(start_row != nseqs-2) fprintf(tree,"(\n");
-
- for(col=1; col<=nseqs; col++) {
- if(tree_description[start_row][col] == flag) {
- test_col = col;
- break;
- }
- }
-
- single_seq = TRUE;
- for(row=start_row-1; row>=1; row--)
- if(tree_description[row][test_col] == 1) {
- single_seq = FALSE;
- new_row = row;
- break;
- }
-
- if(single_seq) {
- tree_description[start_row][test_col] = 0;
- fprintf(tree,"%s",names[test_col+0-1]);
- }
- else {
- for(col=1; col<=nseqs; col++) {
- if((tree_description[start_row][col]==1)&&
- (tree_description[new_row][col]==1))
- tree_description[start_row][col] = 0;
- }
- two_way_split(tree_description, tree, new_row, (int)1, bootstrap);
- }
-
- if(start_row == nseqs-2) {
-/* if (bootstrap && (flag==1)) fprintf(tree,"[TRICHOTOMY]");
-*/
- return;
- }
-
- fprintf(tree,":%7.5f,\n",left_branch[start_row]);
-
- for(col=1; col<=nseqs; col++)
- if(tree_description[start_row][col] == flag) {
- test_col = col;
- break;
- }
-
- single_seq = TRUE;
- for(row=start_row-1; row>=1; row--)
- if(tree_description[row][test_col] == 1) {
- single_seq = FALSE;
- new_row = row;
- break;
- }
-
- if(single_seq) {
- tree_description[start_row][test_col] = 0;
- fprintf(tree,"%s",names[test_col+0-1]);
- }
- else {
- for(col=1; col<=nseqs; col++) {
- if((tree_description[start_row][col]==1)&&
- (tree_description[new_row][col]==1))
- tree_description[start_row][col] = 0;
- }
- two_way_split(tree_description, tree, new_row, (int)1, bootstrap);
- }
-
- fprintf(tree,":%7.5f)\n",right_branch[start_row]);
-
-
-}
-
-
-
-
-/****************************************************************************
- * [ Improvement ideas in fast_nj_tree() ] by DDBJ & FUJITSU Limited.
- * written by Tadashi Koike
- * (takoike@genes.nig.ac.jp)
- *******************
- * <IMPROVEMENT 1> : Store the value of sum of the score to temporary array,
- * and use again and again.
- *
- * In the main cycle, these are calculated again and again :
- * diq = sum of tmat[n][ii] (n:1 to last_seq-first_seq+1),
- * djq = sum of tmat[n][jj] (n:1 to last_seq-first_seq+1),
- * dio = sum of tmat[n][mini] (n:1 to last_seq-first_seq+1),
- * djq = sum of tmat[n][minj] (n:1 to last_seq-first_seq+1)
- * // 'last_seq' and 'first_seq' are both constant values //
- * and the result of above calculations is always same until
- * a best pair of neighbour nodes is joined.
- *
- * So, we change the logic to calculate the sum[i] (=sum of tmat[n][i]
- * (n:1 to last_seq-first_seq+1)) and store it to array, before
- * beginning to find a best pair of neighbour nodes, and after that
- * we use them again and again.
- *
- * tmat[i][j]
- * 1 2 3 4 5
- * +---+---+---+---+---+
- * 1 | | | | | |
- * +---+---+---+---+---+
- * 2 | | | | | | 1) calculate sum of tmat[n][i]
- * +---+---+---+---+---+ (n: 1 to last_seq-first_seq+1)
- * 3 | | | | | | 2) store that sum value to sum[i]
- * +---+---+---+---+---+
- * 4 | | | | | | 3) use sum[i] during finding a best
- * +---+---+---+---+---+ pair of neibour nodes.
- * 5 | | | | | |
- * +---+---+---+---+---+
- * | | | | |
- * V V V V V Calculate sum , and store it to sum[i]
- * +---+---+---+---+---+
- * sum[i] | | | | | |
- * +---+---+---+---+---+
- *
- * At this time, we thought that we use upper triangle of the matrix
- * because tmat[i][j] is equal to tmat[j][i] and tmat[i][i] is equal
- * to zero. Therefore, we prepared sum_rows[i] and sum_cols[i] instead
- * of sum[i] for storing the sum value.
- *
- * tmat[i][j]
- * 1 2 3 4 5 sum_cols[i]
- * +---+---+---+---+---+ +---+
- * 1 | # | # | # | # | --> | | ... sum of tmat[1][2..5]
- * + - +---+---+---+---+ +---+
- * 2 | # | # | # | --> | | ... sum of tmat[2][3..5]
- * + - + - +---+---+---+ +---+
- * 3 | # | # | --> | | ... sum of tmat[3][4..5]
- * + - + - + - +---+---+ +---+
- * 4 | # | --> | | ... sum of tmat[4][5]
- * + - + - + - + - +---+ +---+
- * 5 | --> | | ... zero
- * + - + - + - + - + - + +---+
- * | | | | |
- * V V V V V Calculate sum , sotre to sum[i]
- * +---+---+---+---+---+
- * sum_rows[i] | | | | | |
- * +---+---+---+---+---+
- * | | | | |
- * | | | | +----- sum of tmat[1..4][5]
- * | | | +--------- sum of tmat[1..3][4]
- * | | +------------- sum of tmat[1..2][3]
- * | +----------------- sum of tmat[1][2]
- * +--------------------- zero
- *
- * And we use (sum_rows[i] + sum_cols[i]) instead of sum[i].
- *
- *******************
- * <IMPROVEMENT 2> : We manage valid nodes with chain list, instead of
- * tkill[i] flag array.
- *
- * In original logic, invalid(killed?) nodes after nodes-joining
- * are managed with tkill[i] flag array (set to 1 when killed).
- * By this method, it is conspicuous to try next node but skip it
- * at the latter of finding a best pair of neighbor nodes.
- *
- * So, we thought that we managed valid nodes by using a chain list
- * as below:
- *
- * 1) declare the list structure.
- * struct {
- * int n; // entry number of node.
- * void *prev; // pointer to previous entry.
- * void *next; // pointer to next entry.
- * }
- * 2) construct a valid node list.
- *
- * +-----+ +-----+ +-----+ +-----+ +-----+
- * NULL<-|prev |<---|prev |<---|prev |<---|prev |<- - - -|prev |
- * | 0 | | 1 | | 2 | | 3 | | n |
- * | next|--->| next|--->| next|--->| next|- - - ->| next|->NULL
- * +-----+ +-----+ +-----+ +-----+ +-----+
- *
- * 3) when finding a best pair of neighbor nodes, we use
- * this chain list as loop counter.
- *
- * 4) If an entry was killed by node-joining, this chain list is
- * modified to remove that entry.
- *
- * EX) remove the entry No 2.
- * +-----+ +-----+ +-----+ +-----+
- * NULL<-|prev |<---|prev |<--------------|prev |<- - - -|prev |
- * | 0 | | 1 | | 3 | | n |
- * | next|--->| next|-------------->| next|- - - ->| next|->NULL
- * +-----+ +-----+ +-----+ +-----+
- * +-----+
- * NULL<-|prev |
- * | 2 |
- * | next|->NULL
- * +-----+
- *
- * By this method, speed is up at the latter of finding a best pair of
- * neighbor nodes.
- *
- *******************
- * <IMPROVEMENT 3> : Cut the frequency of division.
- *
- * At comparison between 'total' and 'tmin' in the main cycle, total is
- * divided by (2.0*fnseqs2) before comparison. If N nodes are available,
- * that division happen (N*(N-1))/2 order.
- *
- * We thought that the comparison relation between tmin and total/(2.0*fnseqs2)
- * is equal to the comparison relation between (tmin*2.0*fnseqs2) and total.
- * Calculation of (tmin*2.0*fnseqs2) is only one time. so we stop dividing
- * a total value and multiply tmin and (tmin*2.0*fnseqs2) instead.
- *
- *******************
- * <IMPROVEMENT 4> : some transformation of the equation (to cut operations).
- *
- * We transform an equation of calculating 'total' in the main cycle.
- *
- */
-
-
-void fast_nj_tree(char **tree_description)
-{
- register int i;
- int l[4],nude,k;
- int nc,mini,minj,j,ii,jj;
- double fnseqs,fnseqs2=0,sumd;
- double diq,djq,dij,dio,djo,da;
- double tmin,total,dmin;
- double bi,bj,b1,b2,b3,branch[4];
- int typei,typej; /* 0 = node; 1 = OTU */
-
- /* IMPROVEMENT 1, STEP 0 : declare variables */
- double *sum_cols, *sum_rows, *join;
-
- /* IMPROVEMENT 2, STEP 0 : declare variables */
- int loop_limit;
- typedef struct _ValidNodeID {
- int n;
- struct _ValidNodeID *prev;
- struct _ValidNodeID *next;
- } ValidNodeID;
- ValidNodeID *tvalid, *lpi, *lpj, *lpii, *lpjj, *lp_prev, *lp_next;
-
- /*
- * correspondence of the loop counter variables.
- * i .. lpi->n, ii .. lpii->n
- * j .. lpj->n, jj .. lpjj->n
- */
-
- fnseqs = (double)last_seq-first_seq+1;
-
-/*********************** First initialisation ***************************/
-
-
- if (fnseqs == 2) {
- return;
- }
-
- mini = minj = 0;
-
- left_branch = (double *) ckalloc( (nseqs+2) * sizeof (double) );
- right_branch = (double *) ckalloc( (nseqs+2) * sizeof (double) );
- tkill = (int *) ckalloc( (nseqs+1) * sizeof (int) );
- av = (double *) ckalloc( (nseqs+1) * sizeof (double) );
-
- /* IMPROVEMENT 1, STEP 1 : Allocate memory */
- sum_cols = (double *) ckalloc( (nseqs+1) * sizeof (double) );
- sum_rows = (double *) ckalloc( (nseqs+1) * sizeof (double) );
- join = (double *) ckalloc( (nseqs+1) * sizeof (double) );
-
- /* IMPROVEMENT 2, STEP 1 : Allocate memory */
- tvalid = (ValidNodeID *) ckalloc( (nseqs+1) * sizeof (ValidNodeID) );
- /* tvalid[0] is special entry in array. it points a header of valid entry list */
- tvalid[0].n = 0;
- tvalid[0].prev = NULL;
- tvalid[0].next = &tvalid[1];
-
- /* IMPROVEMENT 2, STEP 2 : Construct and initialize the entry chain list */
- for(i=1, loop_limit = last_seq-first_seq+1,
- lpi=&tvalid[1], lp_prev=&tvalid[0], lp_next=&tvalid[2] ;
- i<=loop_limit ;
- ++i, ++lpi, ++lp_prev, ++lp_next)
- {
- tmat[i][i] = av[i] = 0.0;
- tkill[i] = 0;
- lpi->n = i;
- lpi->prev = lp_prev;
- lpi->next = lp_next;
-
- /* IMPROVEMENT 1, STEP 2 : Initialize arrays */
- sum_cols[i] = sum_rows[i] = join[i] = 0.0;
- }
- tvalid[loop_limit].next = NULL;
-
- /*
- * IMPROVEMENT 1, STEP 3 : Calculate the sum of score value that
- * is sequence[i] to others.
- */
- sumd = 0.0;
- for (lpj=tvalid[0].next ; lpj!=NULL ; lpj = lpj->next) {
- double tmp_sum = 0.0;
- j = lpj->n;
- /* calculate sum_rows[j] */
- for (lpi=tvalid[0].next ; lpi->n < j ; lpi = lpi->next) {
- i = lpi->n;
- tmp_sum += tmat[i][j];
- /* tmat[j][i] = tmat[i][j]; */
- }
- sum_rows[j] = tmp_sum;
-
- tmp_sum = 0.0;
- /* Set lpi to that lpi->n is greater than j */
- if ((lpi != NULL) && (lpi->n == j)) {
- lpi = lpi->next;
- }
- /* calculate sum_cols[j] */
- for( ; lpi!=NULL ; lpi = lpi->next) {
- i = lpi->n;
- tmp_sum += tmat[j][i];
- /* tmat[i][j] = tmat[j][i]; */
- }
- sum_cols[j] = tmp_sum;
- }
-
-/*********************** Enter The Main Cycle ***************************/
-
- for(nc=1, loop_limit = (last_seq-first_seq+1-3); nc<=loop_limit; ++nc) {
-
- sumd = 0.0;
- /* IMPROVEMENT 1, STEP 4 : use sum value */
- for(lpj=tvalid[0].next ; lpj!=NULL ; lpj = lpj->next) {
- sumd += sum_cols[lpj->n];
- }
-
- /* IMPROVEMENT 3, STEP 0 : multiply tmin and 2*fnseqs2 */
- fnseqs2 = fnseqs - 2.0; /* Set fnseqs2 at this point. */
- tmin = 99999.0 * 2.0 * fnseqs2;
-
-
-/*.................compute SMATij values and find the smallest one ........*/
-
- mini = minj = 0;
-
- /* jj must starts at least 2 */
- if ((tvalid[0].next != NULL) && (tvalid[0].next->n == 1)) {
- lpjj = tvalid[0].next->next;
- } else {
- lpjj = tvalid[0].next;
- }
-
- for( ; lpjj != NULL; lpjj = lpjj->next) {
- jj = lpjj->n;
- for(lpii=tvalid[0].next ; lpii->n < jj ; lpii = lpii->next) {
- ii = lpii->n;
- diq = djq = 0.0;
-
- /* IMPROVEMENT 1, STEP 4 : use sum value */
- diq = sum_cols[ii] + sum_rows[ii];
- djq = sum_cols[jj] + sum_rows[jj];
- /*
- * always ii < jj in this point. Use upper
- * triangle of score matrix.
- */
- dij = tmat[ii][jj];
-
- /*
- * IMPROVEMENT 3, STEP 1 : fnseqs2 is
- * already calculated.
- */
- /* fnseqs2 = fnseqs - 2.0 */
-
- /* IMPROVEMENT 4 : transform the equation */
- /*-------------------------------------------------------------------*
- * OPTIMIZE of expression 'total = d2r + fnseqs2*dij + dr*2.0' *
- * total = d2r + fnseq2*dij + 2.0*dr *
- * = d2r + fnseq2*dij + 2(sumd - dij - d2r) *
- * = d2r + fnseq2*dij + 2*sumd - 2*dij - 2*d2r *
- * = fnseq2*dij + 2*sumd - 2*dij - 2*d2r + d2r *
- * = fnseq2*dij + 2*sumd - 2*dij - d2r *
- * = fnseq2*dij + 2*sumd - 2*dij - (diq + djq - 2*dij) *
- * = fnseq2*dij + 2*sumd - 2*dij - diq - djq + 2*dij *
- * = fnseq2*dij + 2*sumd - 2*dij + 2*dij - diq - djq *
- * = fnseq2*dij + 2*sumd - diq - djq *
- *-------------------------------------------------------------------*/
- total = fnseqs2*dij + 2.0*sumd - diq - djq;
-
- /*
- * IMPROVEMENT 3, STEP 2 : abbrevlate
- * the division on comparison between
- * total and tmin.
- */
- /* total = total / (2.0*fnseqs2); */
-
- if(total < tmin) {
- tmin = total;
- mini = ii;
- minj = jj;
- }
- }
- }
-
- /* MEMO: always ii < jj in avobe loop, so mini < minj */
-
-/*.................compute branch lengths and print the results ........*/
-
-
- dio = djo = 0.0;
-
- /* IMPROVEMENT 1, STEP 4 : use sum value */
- dio = sum_cols[mini] + sum_rows[mini];
- djo = sum_cols[minj] + sum_rows[minj];
-
- dmin = tmat[mini][minj];
- dio = (dio - dmin) / fnseqs2;
- djo = (djo - dmin) / fnseqs2;
- bi = (dmin + dio - djo) * 0.5;
- bj = dmin - bi;
- bi = bi - av[mini];
- bj = bj - av[minj];
-
- if( av[mini] > 0.0 )
- typei = 0;
- else
- typei = 1;
- if( av[minj] > 0.0 )
- typej = 0;
- else
- typej = 1;
-
-
-/*
- set negative branch lengths to zero. Also set any tiny positive
- branch lengths to zero.
-*/ if( fabs(bi) < 0.0001) bi = 0.0;
- if( fabs(bj) < 0.0001) bj = 0.0;
-
-
- left_branch[nc] = bi;
- right_branch[nc] = bj;
-
- for(i=1; i<=last_seq-first_seq+1; i++)
- tree_description[nc][i] = 0;
-
- if(typei == 0) {
- for(i=nc-1; i>=1; i--)
- if(tree_description[i][mini] == 1) {
- for(j=1; j<=last_seq-first_seq+1; j++)
- if(tree_description[i][j] == 1)
- tree_description[nc][j] = 1;
- break;
- }
- }
- else
- tree_description[nc][mini] = 1;
-
- if(typej == 0) {
- for(i=nc-1; i>=1; i--)
- if(tree_description[i][minj] == 1) {
- for(j=1; j<=last_seq-first_seq+1; j++)
- if(tree_description[i][j] == 1)
- tree_description[nc][j] = 1;
- break;
- }
- }
- else
- tree_description[nc][minj] = 1;
-
-
-/*
- Here is where the -0.00005 branch lengths come from for 3 or more
- identical seqs.
-*/
-/* if(dmin <= 0.0) dmin = 0.0001; */
- if(dmin <= 0.0) dmin = 0.000001;
- av[mini] = dmin * 0.5;
-
-/*........................Re-initialisation................................*/
-
- fnseqs = fnseqs - 1.0;
- tkill[minj] = 1;
-
- /* IMPROVEMENT 2, STEP 3 : Remove tvalid[minj] from chain list. */
- /* [ Before ]
- * +---------+ +---------+ +---------+
- * |prev |<-------|prev |<-------|prev |<---
- * | n | | n(=minj)| | n |
- * | next|------->| next|------->| next|----
- * +---------+ +---------+ +---------+
- *
- * [ After ]
- * +---------+ +---------+
- * |prev |<--------------------------|prev |<---
- * | n | | n |
- * | next|-------------------------->| next|----
- * +---------+ +---------+
- * +---------+
- * NULL---|prev |
- * | n(=minj)|
- * | next|---NULL
- * +---------+
- */
- (tvalid[minj].prev)->next = tvalid[minj].next;
- if (tvalid[minj].next != NULL) {
- (tvalid[minj].next)->prev = tvalid[minj].prev;
- }
- tvalid[minj].prev = tvalid[minj].next = NULL;
-
- /* IMPROVEMENT 1, STEP 5 : re-calculate sum values. */
- for(lpj=tvalid[0].next ; lpj != NULL ; lpj = lpj->next) {
- double tmp_di = 0.0;
- double tmp_dj = 0.0;
- j = lpj->n;
-
- /*
- * subtrace a score value related with 'minj' from
- * sum arrays .
- */
- if (j < minj) {
- tmp_dj = tmat[j][minj];
- sum_cols[j] -= tmp_dj;
- } else if (j > minj) {
- tmp_dj = tmat[minj][j];
- sum_rows[j] -= tmp_dj;
- } /* nothing to do when j is equal to minj. */
-
-
- /*
- * subtrace a score value related with 'mini' from
- * sum arrays .
- */
- if (j < mini) {
- tmp_di = tmat[j][mini];
- sum_cols[j] -= tmp_di;
- } else if (j > mini) {
- tmp_di = tmat[mini][j];
- sum_rows[j] -= tmp_di;
- } /* nothing to do when j is equal to mini. */
-
- /*
- * calculate a score value of the new inner node.
- * then, store it temporary to join[] array.
- */
- join[j] = (tmp_dj + tmp_di) * 0.5;
- }
-
- /*
- * 1)
- * Set the score values (stored in join[]) into the matrix,
- * row/column position is 'mini'.
- * 2)
- * Add a score value of the new inner node to sum arrays.
- */
- for(lpj=tvalid[0].next ; lpj != NULL; lpj = lpj->next) {
- j = lpj->n;
- if (j < mini) {
- tmat[j][mini] = join[j];
- sum_cols[j] += join[j];
- } else if (j > mini) {
- tmat[mini][j] = join[j];
- sum_rows[j] += join[j];
- } /* nothing to do when j is equal to mini. */
- }
-
- /* Re-calculate sum_rows[mini],sum_cols[mini]. */
- sum_cols[mini] = sum_rows[mini] = 0.0;
-
- /* calculate sum_rows[mini] */
- da = 0.0;
- for(lpj=tvalid[0].next ; lpj->n < mini ; lpj = lpj->next) {
- da += join[lpj->n];
- }
- sum_rows[mini] = da;
-
- /* skip if 'lpj->n' is equal to 'mini' */
- if ((lpj != NULL) && (lpj->n == mini)) {
- lpj = lpj->next;
- }
-
- /* calculate sum_cols[mini] */
- da = 0.0;
- for( ; lpj != NULL; lpj = lpj->next) {
- da += join[lpj->n];
- }
- sum_cols[mini] = da;
-
- /*
- * Clean up sum_rows[minj], sum_cols[minj] and score matrix
- * related with 'minj'.
- */
- sum_cols[minj] = sum_rows[minj] = 0.0;
- for(j=1; j<=last_seq-first_seq+1; ++j)
- tmat[minj][j] = tmat[j][minj] = join[j] = 0.0;
-
-
-/****/ } /*end main cycle**/
-
-/******************************Last Cycle (3 Seqs. left)********************/
-
- nude = 1;
-
- for(lpi=tvalid[0].next; lpi != NULL; lpi = lpi->next) {
- l[nude] = lpi->n;
- ++nude;
- }
-
- b1 = (tmat[l[1]][l[2]] + tmat[l[1]][l[3]] - tmat[l[2]][l[3]]) * 0.5;
- b2 = tmat[l[1]][l[2]] - b1;
- b3 = tmat[l[1]][l[3]] - b1;
-
- branch[1] = b1 - av[l[1]];
- branch[2] = b2 - av[l[2]];
- branch[3] = b3 - av[l[3]];
-
-/* Reset tiny negative and positive branch lengths to zero */
- if( fabs(branch[1]) < 0.0001) branch[1] = 0.0;
- if( fabs(branch[2]) < 0.0001) branch[2] = 0.0;
- if( fabs(branch[3]) < 0.0001) branch[3] = 0.0;
-
- left_branch[last_seq-first_seq+1-2] = branch[1];
- left_branch[last_seq-first_seq+1-1] = branch[2];
- left_branch[last_seq-first_seq+1] = branch[3];
-
- for(i=1; i<=last_seq-first_seq+1; i++)
- tree_description[last_seq-first_seq+1-2][i] = 0;
-
-
- for(i=1; i<=3; ++i) {
- if( av[l[i]] > 0.0) {
-
- for(k=last_seq-first_seq+1-3; k>=1; k--)
- if(tree_description[k][l[i]] == 1) {
- for(j=1; j<=last_seq-first_seq+1; j++)
- if(tree_description[k][j] == 1)
- tree_description[last_seq-first_seq+1-2][j] = i;
- break;
- }
- }
- else {
- tree_description[last_seq-first_seq+1-2][l[i]] = i;
- }
- if(i < 3) {;
- }
- }
- ckfree(sum_cols);
- ckfree(sum_rows);
- ckfree(join);
- ckfree(tvalid);
-}
-//////////////////////////////////////////////////////////////////////////////
-//
-// UPGMA_aln
-//////////////////////////////////////////////////////////////////////////////
-
-Alignment * upgma_merge_aln_rows (Alignment *A, int *ns, int **ls,int N, int**mat,int *used, int *n, Constraint_list *CL);
-int upgma_pair_wise (Alignment *A, int *ls0, int ns0, int *ls2, int ns2, Constraint_list *CL);
-
-
-Alignment * upgma_tree_aln ( Alignment*A, int nseq, Constraint_list *CL)
-{
- int a, b,n, *used;
- static int **mat;
- int **ls;
- int *ns;
- nseq=(CL->S)->nseq;
- mat=declare_int (nseq, nseq);
- ls=declare_int (nseq,nseq);
- ns=vcalloc (nseq,sizeof (int));
-
- for (a=0; a<nseq-1; a++)
- for (b=a+1; b<nseq; b++)
- {
- ns[a]=ns[b]=1;
- ls[a][0]=a;
- ls[b][0]=b;
- mat[a][b]=mat[b][a]=upgma_pair_wise(A,ls[a],ns[a],ls[b],ns[b],CL);
- HERE ("%d %d", a, b, mat[a][b]);
- }
-
- used=vcalloc (nseq, sizeof (int));
- n=nseq;
- while (n>1)
- {
- upgma_merge_aln_rows (A,ns, ls,nseq, mat, used, &n,CL);
- }
- print_aln (A);
- HERE ("finished");
- free_int ( mat, -1);
- free_int (ls, -1);
- vfree (ns);
- return A;
-}
-
-Alignment * upgma_merge_aln_rows (Alignment *A, int *ns, int **ls,int N, int**mat,int *used, int *n, Constraint_list *CL)
-{
-
- int a, b, w, best_a, best_b, set;
- float best_s;
-
- for (set=0,a=0; a<N-1; a++)
- {
- if (used[a])continue;
- for (b=a+1; b<N; b++)
- {
- if (used[b])continue;
- w=mat[a][b];
- if ( !set || w>best_s)
- {
- best_s=w;
- best_a=a;
- best_b=b;
- set=1;
- }
- }
- }
- used[best_b]=1;
-
- //merge a and b
- mat[best_a][best_b]=upgma_pair_wise (A, ls[best_a], ns[best_a], ls[best_b], ns[best_b], CL);
- for (a=0; a<ns[best_b]; a++)ls[best_a][ns[best_a]++]=ls[best_b][a];
- ns[best_b]=0;
-
- //update the a row
- for (a=0; a< A->nseq; a++)
- {
- if (a!=best_a && !used[a])
- mat[best_a][a]=mat[a][best_a]=upgma_pair_wise (A, ls[best_a], ns[best_a], ls[a], ns[a], CL);
- }
-
- HERE ("DONE");
-
- n[0]--;
- return A;
-}
-int upgma_pair_wise (Alignment *A, int *ls0, int ns0, int *ls1, int ns1, Constraint_list *CL)
-{
- static int **ls;
- static int *ns;
- static int *fl;
- int a, b, n;
-
- if ( !ls )
- {
- ls=vcalloc (2, sizeof (int*));
- ns=vcalloc (2, sizeof (int));
- fl=vcalloc ((CL->S)->nseq, sizeof (int));
- }
- ls[0]=ls0;
- ls[1]=ls1;
- ns[0]=ns0; ns[1]=ns1;
-
- fprintf ( stderr, "\n");
- for (a=0; a<ns0; a++)
- fprintf ( stderr, "%d", ls0[a]);
- fprintf ( stderr, "\n");
- for (a=0; a<ns1; a++)
- fprintf ( stderr, "%d", ls1[a]);
-
- ungap_sub_aln (A, ns[0], ls[0]);
- ungap_sub_aln (A, ns[1], ls[1]);
- HERE ("Align");
- pair_wise (A, ns, ls, CL);
- HERE ("Done");
-
- for (n=0,a=0;a<2; a++)
- for (b=0; b<ns[a]; b++)
- fl[n++]=ls[a][b];
-
- return sub_aln2ecl_raw_score (A,CL,n, fl);
-}
-
-//////////////////////////////////////////////////////////////////////////////
-//
-// UPGMA
-///////////////////////////////////////////////////////////////////////////////
-
-
-NT_node ** int_dist2upgma_tree (int **mat, Alignment *A, int nseq, char *fname)
-{
- NT_node *NL, T;
- int a, n, *used;
- int tot_node;
- if (upgma_node_heap (NULL))
- {
- printf_exit ( EXIT_FAILURE,stderr, "\nERROR: non empty heap in upgma [FATAL]");
- }
- NL=vcalloc (nseq, sizeof (NT_node));
-
- for (a=0; a<A->nseq; a++)
- {
- NL[a]=new_declare_tree_node ();
- upgma_node_heap (NL[a]);
- sprintf (NL[a]->name, "%s", A->name[a]);
- NL[a]->isseq=1;
- NL[a]->leaf=1;
- }
-
- used=vcalloc ( A->nseq, sizeof (int));
- n=A->nseq;
- while (n>1)
- {
- T=upgma_merge (mat, NL,used, &n, A->nseq);
- }
-
- vfree (used);
- vfclose (print_tree (T, "newick", vfopen (fname, "w")));
- upgma_node_heap (NULL);
- vfree (NL);
-
- return read_tree (fname,&tot_node,A->nseq, A->name);
-}
-NT_node upgma_merge (int **mat, NT_node *NL, int *used, int *n, int N)
-{
- NT_node P, LC, RC;
- int a, b, w, best_a, best_b, set;
- float best_s;
- P=new_declare_tree_node();
- upgma_node_heap (P);
-
- for (set=0,a=0; a<N-1; a++)
- {
- if (used[a])continue;
- for (b=a+1; b<N; b++)
- {
- if (used[b])continue;
- w=mat[a][b];
- if ( !set || w<best_s)
- {
- best_s=w;
- best_a=a;
- best_b=b;
- set=1;
- }
- }
- }
-
- for (a=0; a<N; a++)
- {
- if (!used[a])mat[best_a][a]=mat[a][best_a]=(mat[best_b][a]+mat[best_a][a])/2;
- }
- used[best_b]=1;
-
- LC=NL[best_a];
- RC=NL[best_b];
- best_s/=(float)100;
- LC->dist=RC->dist=best_s;
- LC->parent=RC->parent=P;
- P->left=LC;
- P->right=RC;
- NL[best_a]=P;
- n[0]--;
- return P;
-
-}
-
-
-//////////////////////////////////////////////////////////////////////////////
-//
-// SPLIT UPGMA
-///////////////////////////////////////////////////////////////////////////////
-
-int upgma_node_heap (NT_node X)
-{
- static int n;
- static NT_node *h;
- if ( X==NULL)
- {
- int a,r;
- if (n==0) return 0;
- for (a=0; a<n; a++)
- {
- free_tree_node (h[a]);
- }
- vfree (h);
- h=NULL;
- r=n;
- n=0;
- return r;
- }
- else
- {
- h=vrealloc (h, sizeof (NT_node)*(n+1));
- h[n++]=X;
- }
- return n;
-}
-NT_node split2upgma_tree (Split **S, Alignment *A, int nseq, char *fname)
-{
- NT_node *NL, T;
- int a, n, *used;
-
- NL=vcalloc (nseq, sizeof (NT_node));
- for (a=0; a<A->nseq; a++)
- {
-
- NL[a]=new_declare_tree_node ();
- NL[a]->lseq2=vcalloc (A->nseq+1, sizeof (int));
- NL[a]->lseq2[a]=1;
- sprintf (NL[a]->name, "%s", A->name[a]);
- NL[a]->isseq=1;
- NL[a]->leaf=1;
- NL[a]->dist=1;
- upgma_node_heap (NL[a]);
- }
- used=vcalloc ( A->nseq, sizeof (int));
- n=A->nseq;
- while (n>1)
- {
- T=split_upgma_merge (A,S, NL,used, &n, A->nseq);
- }
- vfree (used);
- fprintf ( stdout, "\n");
- vfclose (print_tree (T, "newick", vfopen (fname, "w")));
- upgma_node_heap (NULL);
- return T;
-}
-NT_node split_upgma_merge (Alignment *A, Split **S, NT_node *NL, int *used, int *n, int N)
-{
- NT_node P, LC, RC;
- int a, b, w, best_a, best_b, set;
- float best_s;
- static int **mat;
-
- if (!mat)
- {
- mat=declare_int (N, N);
- for (a=0; a<N-1; a++)
- for ( b=a+1; b<N; b++)
- {
- mat[a][b]=get_split_dist (A, NL[a], NL[b], S);
- }
- }
-
- P=new_declare_tree_node();
- upgma_node_heap (P);
- P->lseq2=vcalloc (N, sizeof (int));
- for (set=0,a=0; a<N-1; a++)
- {
- if (used[a])continue;
- for (b=a+1; b<N; b++)
- {
- if (used[b])continue;
- w=mat[a][b];
- if ( !set || w<best_s)
- {
- best_s=w;
- best_a=a;
- best_b=b;
- set=1;
- }
- }
- }
-
-
-
-
- LC=NL[best_a];
- RC=NL[best_b];
- best_s/=100;
-
- P->dist=1-best_s;
- LC->parent=RC->parent=P;
- P->left=LC;
- P->right=RC;
- P->bootstrap=best_s*100;
- used[best_b]=1;
-
- n[0]--;
-
- for (a=0; a<A->nseq; a++)
- {
- P->lseq2[a]=(LC->lseq2[a] || RC->lseq2[a])?1:0;
- }
-
- for (a=0; a<N; a++)
- {
- if (!used[a])mat[best_a][a]=mat[a][best_a]=(int)get_split_dist(A,P, NL[a], S);
- }
- NL[best_a]=P;
-
- return P;
-
-}
-float get_split_dist ( Alignment *A, NT_node L, NT_node R, Split **S)
-{
- static char *split;
-
-
- int n,a;
-
- if (!split)
- {
- split=vcalloc ( A->nseq+1, sizeof (char));
- }
-
-
-
- for ( a=0; a<A->nseq; a++)
- split[a]=((L->lseq2[a] || R->lseq2[a])?1:0)+'0';
-
- n=0;
- while (S[n])
- {
- float score;
- if ( strm (S[n]->split,split))
- {
- return score=100-S[n]->score;
- }
- n++;
- }
- return 100;
-}
-/******************************COPYRIGHT NOTICE*******************************/
-/*© Centro de Regulacio Genomica */
-/*and */
-/*Cedric Notredame */
-/*Fri Feb 18 08:27:45 CET 2011 - Revision 596. */
-/*All rights reserved.*/
-/*This file is part of T-COFFEE.*/
-/**/
-/* T-COFFEE is free software; you can redistribute it and/or modify*/
-/* it under the terms of the GNU General Public License as published by*/
-/* the Free Software Foundation; either version 2 of the License, or*/
-/* (at your option) any later version.*/
-/**/
-/* T-COFFEE is distributed in the hope that it will be useful,*/
-/* but WITHOUT ANY WARRANTY; without even the implied warranty of*/
-/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the*/
-/* GNU General Public License for more details.*/
-/**/
-/* You should have received a copy of the GNU General Public License*/
-/* along with Foobar; if not, write to the Free Software*/
-/* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA*/
-/*............................................... |*/
-/* If you need some more information*/
-/* cedric.notredame@europe.com*/
-/*............................................... |*/
-/**/
-/**/
-/* */
-/******************************COPYRIGHT NOTICE*******************************/
git://source.jalview.org / jabaws.git / blobdiff