| 1 | /* Phyle of filogenetic tree calculating functions for CLUSTAL V */ |
|---|
| 2 | #include <stdio.h> |
|---|
| 3 | #include <string.h> |
|---|
| 4 | #include <stdlib.h> |
|---|
| 5 | #include <math.h> |
|---|
| 6 | #include "clustalv.h" |
|---|
| 7 | |
|---|
| 8 | /* |
|---|
| 9 | * Prototypes |
|---|
| 10 | */ |
|---|
| 11 | |
|---|
| 12 | extern void *ckalloc(size_t); |
|---|
| 13 | extern int getint(char *, int, int, int); |
|---|
| 14 | extern void get_path(char *, char *); |
|---|
| 15 | extern FILE * open_output_file(char *, char *, char *, char *); |
|---|
| 16 | |
|---|
| 17 | |
|---|
| 18 | void init_trees(void); |
|---|
| 19 | void phylogenetic_tree(void); |
|---|
| 20 | void bootstrap_tree(void); |
|---|
| 21 | void compare_tree(char **, char **, int *, int); |
|---|
| 22 | void tree_gap_delete(void); /*flag all positions in alignment that have a gap */ |
|---|
| 23 | void dna_distance_matrix(FILE *); |
|---|
| 24 | void prot_distance_matrix(FILE *); |
|---|
| 25 | void nj_tree(char **, FILE *); |
|---|
| 26 | void print_tree(char **, FILE *, int *); |
|---|
| 27 | void root_tree(char **, int); |
|---|
| 28 | |
|---|
| 29 | Boolean transition(int,int); |
|---|
| 30 | |
|---|
| 31 | /* |
|---|
| 32 | * Global variables |
|---|
| 33 | */ |
|---|
| 34 | |
|---|
| 35 | |
|---|
| 36 | extern double **tmat; /* general nxn array of reals; allocated from main */ |
|---|
| 37 | /* this is used as a distance matrix */ |
|---|
| 38 | extern double **smat; |
|---|
| 39 | extern Boolean dnaflag; /* TRUE for DNA seqs; FALSE for proteins */ |
|---|
| 40 | extern Boolean tossgaps; /* Ignore places in align. where ANY seq. has a gap*/ |
|---|
| 41 | extern Boolean kimura; /* Use correction for multiple substitutions */ |
|---|
| 42 | extern Boolean empty; /* any sequences in memory? */ |
|---|
| 43 | extern Boolean usemenu; /* interactive (TRUE) or command line (FALSE) */ |
|---|
| 44 | extern void error(char *,...); /* error reporting */ |
|---|
| 45 | extern int nseqs; /* total no. of seqs. */ |
|---|
| 46 | extern int *seqlen_array; /* the lengths of the sequences */ |
|---|
| 47 | extern char **seq_array; /* the sequences */ |
|---|
| 48 | extern char **names; /* the seq. names */ |
|---|
| 49 | extern char seqname[]; /* name of input file */ |
|---|
| 50 | |
|---|
| 51 | static double *av; |
|---|
| 52 | static int *kill; |
|---|
| 53 | static int ran_factor; |
|---|
| 54 | int boot_ntrials; /* number of bootstrap trials */ |
|---|
| 55 | unsigned int boot_ran_seed; /* random number generator seed */ |
|---|
| 56 | static int root_sequence; |
|---|
| 57 | static int *boot_positions; |
|---|
| 58 | static int *boot_totals; |
|---|
| 59 | static char **standard_tree; |
|---|
| 60 | static char **sample_tree; |
|---|
| 61 | static FILE *phy_tree_file; |
|---|
| 62 | static char phy_tree_name[FILENAMELEN+1]; |
|---|
| 63 | static Boolean verbose; |
|---|
| 64 | static char *tree_gaps; /* array of weights; 1 for use this posn.; 0 don't */ |
|---|
| 65 | |
|---|
| 66 | |
|---|
| 67 | void init_trees() |
|---|
| 68 | { |
|---|
| 69 | register int i; |
|---|
| 70 | |
|---|
| 71 | tree_gaps = (char *)ckalloc( (MAXLEN+1) * sizeof (char) ); |
|---|
| 72 | |
|---|
| 73 | boot_positions = (int *)ckalloc( (MAXLEN+1) * sizeof (int) ); |
|---|
| 74 | boot_totals = (int *)ckalloc( (MAXN+1) * sizeof (int) ); |
|---|
| 75 | |
|---|
| 76 | kill = (int *) ckalloc( (MAXN+1) * sizeof (int) ); |
|---|
| 77 | av = (double *) ckalloc( (MAXN+1) * sizeof (double) ); |
|---|
| 78 | |
|---|
| 79 | standard_tree = (char **) ckalloc( (MAXN+1) * sizeof (char *) ); |
|---|
| 80 | for(i=0; i<MAXN+1; i++) |
|---|
| 81 | standard_tree[i] = (char *) ckalloc( (MAXN+1) * sizeof(char) ); |
|---|
| 82 | |
|---|
| 83 | sample_tree = (char **) ckalloc( (MAXN+1) * sizeof (char *) ); |
|---|
| 84 | for(i=0; i<MAXN+1; i++) |
|---|
| 85 | sample_tree[i] = (char *) ckalloc( (MAXN+1) * sizeof(char) ); |
|---|
| 86 | |
|---|
| 87 | boot_ntrials = 1000; |
|---|
| 88 | boot_ran_seed = 111; |
|---|
| 89 | kimura = FALSE; |
|---|
| 90 | tossgaps = FALSE; |
|---|
| 91 | } |
|---|
| 92 | |
|---|
| 93 | |
|---|
| 94 | |
|---|
| 95 | void phylogenetic_tree() |
|---|
| 96 | { char path[FILENAMELEN+1]; |
|---|
| 97 | int j; |
|---|
| 98 | |
|---|
| 99 | if(empty) { |
|---|
| 100 | error("You must load an alignment first"); |
|---|
| 101 | return; |
|---|
| 102 | } |
|---|
| 103 | |
|---|
| 104 | get_path(seqname,path); |
|---|
| 105 | |
|---|
| 106 | if((phy_tree_file = open_output_file( |
|---|
| 107 | "\nEnter name for tree output file ",path, |
|---|
| 108 | phy_tree_name,"nj")) == NULL) return; |
|---|
| 109 | |
|---|
| 110 | for(j=1; j<=seqlen_array[1]; ++j) |
|---|
| 111 | boot_positions[j] = j; |
|---|
| 112 | |
|---|
| 113 | verbose = TRUE; /* Turn on screen output */ |
|---|
| 114 | if(dnaflag) |
|---|
| 115 | dna_distance_matrix(phy_tree_file); |
|---|
| 116 | else |
|---|
| 117 | prot_distance_matrix(phy_tree_file); |
|---|
| 118 | |
|---|
| 119 | verbose = TRUE; /* Turn on output */ |
|---|
| 120 | nj_tree(standard_tree,phy_tree_file); |
|---|
| 121 | fclose(phy_tree_file); |
|---|
| 122 | /* |
|---|
| 123 | print_tree(standard_tree,phy_tree_file); |
|---|
| 124 | */ |
|---|
| 125 | fprintf(stdout,"\nPhylogenetic tree file created: [%s]",phy_tree_name); |
|---|
| 126 | } |
|---|
| 127 | |
|---|
| 128 | |
|---|
| 129 | |
|---|
| 130 | |
|---|
| 131 | |
|---|
| 132 | Boolean transition(int base1, int base2) /* TRUE if transition; else FALSE */ |
|---|
| 133 | /* |
|---|
| 134 | |
|---|
| 135 | assumes that the bases of DNA sequences have been translated as |
|---|
| 136 | a,A = 1; c,C = 2; g,G = 3; t,T,u,U = 4; X or N = 0; "-" < 0; |
|---|
| 137 | |
|---|
| 138 | A <--> G and T <--> C are transitions; all others are transversions. |
|---|
| 139 | |
|---|
| 140 | */ |
|---|
| 141 | { |
|---|
| 142 | if( ((base1 == 1) && (base2 == 3)) || ((base1 == 3) && (base2 == 1)) ) |
|---|
| 143 | return TRUE; /* A <--> G */ |
|---|
| 144 | if( ((base1 == 4) && (base2 == 2)) || ((base1 == 2) && (base2 == 4)) ) |
|---|
| 145 | return TRUE; /* T <--> C */ |
|---|
| 146 | return FALSE; |
|---|
| 147 | } |
|---|
| 148 | |
|---|
| 149 | |
|---|
| 150 | void tree_gap_delete() /* flag all positions in alignment that have a gap */ |
|---|
| 151 | { /* in ANY sequence */ |
|---|
| 152 | int seqn, posn; |
|---|
| 153 | |
|---|
| 154 | for(posn=1; posn<=seqlen_array[1]; ++posn) { |
|---|
| 155 | tree_gaps[posn] = 0; |
|---|
| 156 | for(seqn=1; seqn<=nseqs; ++seqn) { |
|---|
| 157 | if(seq_array[seqn][posn] <= 0) { |
|---|
| 158 | tree_gaps[posn] = 1; |
|---|
| 159 | break; |
|---|
| 160 | } |
|---|
| 161 | } |
|---|
| 162 | } |
|---|
| 163 | } |
|---|
| 164 | |
|---|
| 165 | |
|---|
| 166 | |
|---|
| 167 | void nj_tree(char **tree_description, FILE *tree) |
|---|
| 168 | { |
|---|
| 169 | register int i; |
|---|
| 170 | int l[4],nude,k; |
|---|
| 171 | int nc,mini,minj,j,j1,ii,jj; |
|---|
| 172 | double fnseqs,fnseqs2,sumd; |
|---|
| 173 | double diq,djq,dij,d2r,dr,dio,djo,da; |
|---|
| 174 | double tmin,total,dmin; |
|---|
| 175 | double bi,bj,b1,b2,b3,branch[4]; |
|---|
| 176 | int typei,typej,type[4]; /* 0 = node; 1 = OTU */ |
|---|
| 177 | |
|---|
| 178 | fnseqs = (double)nseqs; |
|---|
| 179 | |
|---|
| 180 | /*********************** First initialisation ***************************/ |
|---|
| 181 | |
|---|
| 182 | if(verbose) { |
|---|
| 183 | fprintf(tree,"\n\n\t\t\tNeighbor-joining Method\n"); |
|---|
| 184 | fprintf(tree,"\n Saitou, N. and Nei, M. (1987)"); |
|---|
| 185 | fprintf(tree," The Neighbor-joining Method:"); |
|---|
| 186 | fprintf(tree,"\n A New Method for Reconstructing Phylogenetic Trees."); |
|---|
| 187 | fprintf(tree,"\n Mol. Biol. Evol., 4(4), 406-425\n"); |
|---|
| 188 | fprintf(tree,"\n\n This is an UNROOTED tree\n"); |
|---|
| 189 | fprintf(tree,"\n Numbers in parentheses are branch lengths\n\n"); |
|---|
| 190 | } |
|---|
| 191 | |
|---|
| 192 | mini = minj = 0; |
|---|
| 193 | |
|---|
| 194 | for(i=1;i<=nseqs;++i) |
|---|
| 195 | { |
|---|
| 196 | tmat[i][i] = av[i] = 0.0; |
|---|
| 197 | kill[i] = 0; |
|---|
| 198 | } |
|---|
| 199 | |
|---|
| 200 | /*********************** Enter The Main Cycle ***************************/ |
|---|
| 201 | |
|---|
| 202 | /* for(nc=1; nc<=(nseqs-3); ++nc) { */ /**start main cycle**/ |
|---|
| 203 | for(nc=1; nc<=(nseqs-3); ++nc) { |
|---|
| 204 | sumd = 0.0; |
|---|
| 205 | for(j=2; j<=nseqs; ++j) |
|---|
| 206 | for(i=1; i<j; ++i) { |
|---|
| 207 | tmat[j][i] = tmat[i][j]; |
|---|
| 208 | sumd = sumd + tmat[i][j]; |
|---|
| 209 | smat[i][j] = smat[j][i] = 0.0; |
|---|
| 210 | } |
|---|
| 211 | |
|---|
| 212 | tmin = 99999.0; |
|---|
| 213 | |
|---|
| 214 | /*.................compute SMATij values and find the smallest one ........*/ |
|---|
| 215 | |
|---|
| 216 | for(jj=2; jj<=nseqs; ++jj) |
|---|
| 217 | if(kill[jj] != 1) |
|---|
| 218 | for(ii=1; ii<jj; ++ii) |
|---|
| 219 | if(kill[ii] != 1) { |
|---|
| 220 | diq = djq = 0.0; |
|---|
| 221 | |
|---|
| 222 | for(i=1; i<=nseqs; ++i) { |
|---|
| 223 | diq = diq + tmat[i][ii]; |
|---|
| 224 | djq = djq + tmat[i][jj]; |
|---|
| 225 | } |
|---|
| 226 | |
|---|
| 227 | dij = tmat[ii][jj]; |
|---|
| 228 | d2r = diq + djq - (2.0*dij); |
|---|
| 229 | dr = sumd - dij -d2r; |
|---|
| 230 | fnseqs2 = fnseqs - 2.0; |
|---|
| 231 | total= d2r+ fnseqs2*dij +dr*2.0; |
|---|
| 232 | total= total / (2.0*fnseqs2); |
|---|
| 233 | smat[ii][jj] = total; |
|---|
| 234 | |
|---|
| 235 | if(total < tmin) { |
|---|
| 236 | tmin = total; |
|---|
| 237 | mini = ii; |
|---|
| 238 | minj = jj; |
|---|
| 239 | } |
|---|
| 240 | } |
|---|
| 241 | |
|---|
| 242 | |
|---|
| 243 | /*.................compute branch lengths and print the results ........*/ |
|---|
| 244 | |
|---|
| 245 | |
|---|
| 246 | dio = djo = 0.0; |
|---|
| 247 | for(i=1; i<=nseqs; ++i) { |
|---|
| 248 | dio = dio + tmat[i][mini]; |
|---|
| 249 | djo = djo + tmat[i][minj]; |
|---|
| 250 | } |
|---|
| 251 | |
|---|
| 252 | dmin = tmat[mini][minj]; |
|---|
| 253 | dio = (dio - dmin) / fnseqs2; |
|---|
| 254 | djo = (djo - dmin) / fnseqs2; |
|---|
| 255 | bi = (dmin + dio - djo) * 0.5; |
|---|
| 256 | bj = dmin - bi; |
|---|
| 257 | bi = bi - av[mini]; |
|---|
| 258 | bj = bj - av[minj]; |
|---|
| 259 | |
|---|
| 260 | if( av[mini] > 0.0 ) |
|---|
| 261 | typei = 0; |
|---|
| 262 | else |
|---|
| 263 | typei = 1; |
|---|
| 264 | if( av[minj] > 0.0 ) |
|---|
| 265 | typej = 0; |
|---|
| 266 | else |
|---|
| 267 | typej = 1; |
|---|
| 268 | |
|---|
| 269 | if(verbose) { |
|---|
| 270 | fprintf(tree,"\n Cycle%4d = ",nc); |
|---|
| 271 | if(typei == 0) |
|---|
| 272 | fprintf(tree,"Node:%4d (%9.5f) joins ",mini,bi); |
|---|
| 273 | else |
|---|
| 274 | fprintf(tree," SEQ:%4d (%9.5f) joins ",mini,bi); |
|---|
| 275 | if(typej == 0) |
|---|
| 276 | fprintf(tree,"Node:%4d (%9.5f)",minj,bj); |
|---|
| 277 | else |
|---|
| 278 | fprintf(tree," SEQ:%4d (%9.5f)",minj,bj); |
|---|
| 279 | fprintf(tree,"\n"); |
|---|
| 280 | } |
|---|
| 281 | |
|---|
| 282 | for(i=1; i<=nseqs; i++) |
|---|
| 283 | tree_description[nc][i] = 0; |
|---|
| 284 | |
|---|
| 285 | if(typei == 0) { |
|---|
| 286 | for(i=nc-1; i>=1; i--) |
|---|
| 287 | if(tree_description[i][mini] == 1) { |
|---|
| 288 | for(j=1; j<=nseqs; j++) |
|---|
| 289 | if(tree_description[i][j] == 1) |
|---|
| 290 | tree_description[nc][j] = 1; |
|---|
| 291 | break; |
|---|
| 292 | } |
|---|
| 293 | } |
|---|
| 294 | else |
|---|
| 295 | tree_description[nc][mini] = 1; |
|---|
| 296 | |
|---|
| 297 | if(typej == 0) { |
|---|
| 298 | for(i=nc-1; i>=1; i--) |
|---|
| 299 | if(tree_description[i][minj] == 1) { |
|---|
| 300 | for(j=1; j<=nseqs; j++) |
|---|
| 301 | if(tree_description[i][j] == 1) |
|---|
| 302 | tree_description[nc][j] = 1; |
|---|
| 303 | break; |
|---|
| 304 | } |
|---|
| 305 | } |
|---|
| 306 | else |
|---|
| 307 | tree_description[nc][minj] = 1; |
|---|
| 308 | |
|---|
| 309 | if(dmin <= 0.0) dmin = 0.0001; |
|---|
| 310 | av[mini] = dmin * 0.5; |
|---|
| 311 | |
|---|
| 312 | /*........................Re-initialisation................................*/ |
|---|
| 313 | |
|---|
| 314 | fnseqs = fnseqs - 1.0; |
|---|
| 315 | kill[minj] = 1; |
|---|
| 316 | |
|---|
| 317 | for(j=1; j<=nseqs; ++j) |
|---|
| 318 | if( kill[j] != 1 ) { |
|---|
| 319 | da = ( tmat[mini][j] + tmat[minj][j] ) * 0.5; |
|---|
| 320 | if( (mini - j) < 0 ) |
|---|
| 321 | tmat[mini][j] = da; |
|---|
| 322 | if( (mini - j) > 0) |
|---|
| 323 | tmat[j][mini] = da; |
|---|
| 324 | } |
|---|
| 325 | |
|---|
| 326 | for(j=1; j<=nseqs; ++j) |
|---|
| 327 | tmat[minj][j] = tmat[j][minj] = 0.0; |
|---|
| 328 | |
|---|
| 329 | |
|---|
| 330 | /****/ } /**end main cycle**/ |
|---|
| 331 | |
|---|
| 332 | /******************************Last Cycle (3 Seqs. left)********************/ |
|---|
| 333 | |
|---|
| 334 | nude = 1; |
|---|
| 335 | |
|---|
| 336 | for(i=1; i<=nseqs; ++i) |
|---|
| 337 | if( kill[i] != 1 ) { |
|---|
| 338 | l[nude] = i; |
|---|
| 339 | nude = nude + 1; |
|---|
| 340 | } |
|---|
| 341 | |
|---|
| 342 | b1 = (tmat[l[1]][l[2]] + tmat[l[1]][l[3]] - tmat[l[2]][l[3]]) * 0.5; |
|---|
| 343 | b2 = tmat[l[1]][l[2]] - b1; |
|---|
| 344 | b3 = tmat[l[1]][l[3]] - b1; |
|---|
| 345 | branch[1] = b1 - av[l[1]]; |
|---|
| 346 | branch[2] = b2 - av[l[2]]; |
|---|
| 347 | branch[3] = b3 - av[l[3]]; |
|---|
| 348 | |
|---|
| 349 | |
|---|
| 350 | for(i=1; i<=nseqs; i++) |
|---|
| 351 | tree_description[nseqs-2][i] = 0; |
|---|
| 352 | |
|---|
| 353 | if(verbose) |
|---|
| 354 | fprintf(tree,"\n Cycle%4d (Last cycle, trichotomy):\n",nc); |
|---|
| 355 | |
|---|
| 356 | for(i=1; i<=3; ++i) { |
|---|
| 357 | if( av[l[i]] > 0.0) { |
|---|
| 358 | if(verbose) |
|---|
| 359 | fprintf(tree,"\n\t\t Node:%4d (%9.5f) ",l[i],branch[i]); |
|---|
| 360 | for(k=nseqs-3; k>=1; k--) |
|---|
| 361 | if(tree_description[k][l[i]] == 1) { |
|---|
| 362 | for(j=1; j<=nseqs; j++) |
|---|
| 363 | if(tree_description[k][j] == 1) |
|---|
| 364 | tree_description[nseqs-2][j] = i; |
|---|
| 365 | break; |
|---|
| 366 | } |
|---|
| 367 | } |
|---|
| 368 | else { |
|---|
| 369 | if(verbose) |
|---|
| 370 | fprintf(tree,"\n\t\t SEQ:%4d (%9.5f) ",l[i],branch[i]); |
|---|
| 371 | tree_description[nseqs-2][l[i]] = i; |
|---|
| 372 | } |
|---|
| 373 | if(i < 3) { |
|---|
| 374 | if(verbose) |
|---|
| 375 | fprintf(tree,"joins"); |
|---|
| 376 | } |
|---|
| 377 | } |
|---|
| 378 | |
|---|
| 379 | if(verbose) |
|---|
| 380 | fprintf(tree,"\n"); |
|---|
| 381 | |
|---|
| 382 | } |
|---|
| 383 | |
|---|
| 384 | |
|---|
| 385 | |
|---|
| 386 | |
|---|
| 387 | void bootstrap_tree() |
|---|
| 388 | { |
|---|
| 389 | int i,j,ranno; |
|---|
| 390 | int k,l,m,p; |
|---|
| 391 | unsigned int num; |
|---|
| 392 | char lin2[MAXLINE],path[MAXLINE+1]; |
|---|
| 393 | |
|---|
| 394 | if(empty) { |
|---|
| 395 | error("You must load an alignment first"); |
|---|
| 396 | return; |
|---|
| 397 | } |
|---|
| 398 | |
|---|
| 399 | get_path(seqname, path); |
|---|
| 400 | |
|---|
| 401 | if((phy_tree_file = open_output_file( |
|---|
| 402 | "\nEnter name for bootstrap output file ",path, |
|---|
| 403 | phy_tree_name,"njb")) == NULL) return; |
|---|
| 404 | |
|---|
| 405 | for(i=0;i<MAXN+1;i++) |
|---|
| 406 | boot_totals[i]=0; |
|---|
| 407 | |
|---|
| 408 | for(j=1; j<=seqlen_array[1]; ++j) /* First select all positions for */ |
|---|
| 409 | boot_positions[j] = j; /* the "standard" tree */ |
|---|
| 410 | |
|---|
| 411 | verbose = TRUE; /* Turn on screen output */ |
|---|
| 412 | if(dnaflag) |
|---|
| 413 | dna_distance_matrix(phy_tree_file); |
|---|
| 414 | else |
|---|
| 415 | prot_distance_matrix(phy_tree_file); |
|---|
| 416 | |
|---|
| 417 | verbose = TRUE; /* Turn on screen output */ |
|---|
| 418 | nj_tree(standard_tree, phy_tree_file); /* compute the standard tree */ |
|---|
| 419 | |
|---|
| 420 | fprintf(phy_tree_file,"\n\n\t\t\tBootstrap Confidence Limits\n\n"); |
|---|
| 421 | |
|---|
| 422 | ran_factor = RAND_MAX / seqlen_array[1]; |
|---|
| 423 | |
|---|
| 424 | if(usemenu) |
|---|
| 425 | boot_ran_seed = |
|---|
| 426 | getint("\n\nEnter seed no. for random number generator ",1,1000,boot_ran_seed); |
|---|
| 427 | |
|---|
| 428 | srand(boot_ran_seed); |
|---|
| 429 | fprintf(phy_tree_file,"\n Random number generator seed = %7u\n", |
|---|
| 430 | boot_ran_seed); |
|---|
| 431 | |
|---|
| 432 | if(usemenu) |
|---|
| 433 | boot_ntrials = |
|---|
| 434 | getint("\n\nEnter number of bootstrap trials ",1,10000,boot_ntrials); |
|---|
| 435 | |
|---|
| 436 | fprintf(phy_tree_file,"\n Number of bootstrap trials = %7d\n", |
|---|
| 437 | boot_ntrials); |
|---|
| 438 | |
|---|
| 439 | fprintf(phy_tree_file, |
|---|
| 440 | "\n\n Diagrammatic representation of the above tree: \n"); |
|---|
| 441 | fprintf(phy_tree_file,"\n Each row represents 1 tree cycle;"); |
|---|
| 442 | fprintf(phy_tree_file," defining 2 groups.\n"); |
|---|
| 443 | fprintf(phy_tree_file,"\n Each column is 1 sequence; "); |
|---|
| 444 | fprintf(phy_tree_file,"the stars in each line show 1 group; "); |
|---|
| 445 | fprintf(phy_tree_file,"\n the dots show the other\n"); |
|---|
| 446 | fprintf(phy_tree_file,"\n Numbers show occurrences in bootstrap samples."); |
|---|
| 447 | /* |
|---|
| 448 | print_tree(standard_tree, phy_tree_file, boot_totals); |
|---|
| 449 | */ |
|---|
| 450 | verbose = FALSE; /* Turn OFF screen output */ |
|---|
| 451 | |
|---|
| 452 | fprintf(stdout,"\n\nEach dot represents 10 trials\n\n"); |
|---|
| 453 | for(i=1; i<=boot_ntrials; ++i) { |
|---|
| 454 | for(j=1; j<=seqlen_array[1]; ++j) { /* select alignment */ |
|---|
| 455 | ranno = ( rand() / ran_factor ) + 1; /* positions for */ |
|---|
| 456 | boot_positions[j] = ranno; /* bootstrap sample */ |
|---|
| 457 | } |
|---|
| 458 | if(dnaflag) |
|---|
| 459 | dna_distance_matrix(phy_tree_file); |
|---|
| 460 | else |
|---|
| 461 | prot_distance_matrix(phy_tree_file); |
|---|
| 462 | nj_tree(sample_tree, phy_tree_file); /* compute 1 sample tree */ |
|---|
| 463 | compare_tree(standard_tree, sample_tree, boot_totals, nseqs); |
|---|
| 464 | if(i % 10 == 0) fprintf(stdout,"."); |
|---|
| 465 | if(i % 100 == 0) fprintf(stdout,"\n"); |
|---|
| 466 | } |
|---|
| 467 | |
|---|
| 468 | /* |
|---|
| 469 | fprintf(phy_tree_file,"\n\n Bootstrap totals for each group\n"); |
|---|
| 470 | */ |
|---|
| 471 | print_tree(standard_tree, phy_tree_file, boot_totals); |
|---|
| 472 | |
|---|
| 473 | fclose(phy_tree_file); |
|---|
| 474 | |
|---|
| 475 | fprintf(stdout,"\n\nBootstrap output file completed [%s]" |
|---|
| 476 | ,phy_tree_name); |
|---|
| 477 | } |
|---|
| 478 | |
|---|
| 479 | |
|---|
| 480 | void compare_tree(char **tree1, char **tree2, int *hits, int n) |
|---|
| 481 | { |
|---|
| 482 | int i,j,k,l; |
|---|
| 483 | int nhits1, nhits2; |
|---|
| 484 | |
|---|
| 485 | for(i=1; i<=n-3; i++) { |
|---|
| 486 | for(j=1; j<=n-3; j++) { |
|---|
| 487 | nhits1 = 0; |
|---|
| 488 | nhits2 = 0; |
|---|
| 489 | for(k=1; k<=n; k++) { |
|---|
| 490 | if(tree1[i][k] == tree2[j][k]) nhits1++; |
|---|
| 491 | if(tree1[i][k] != tree2[j][k]) nhits2++; |
|---|
| 492 | } |
|---|
| 493 | if((nhits1 == nseqs) || (nhits2 == nseqs)) hits[i]++; |
|---|
| 494 | } |
|---|
| 495 | } |
|---|
| 496 | } |
|---|
| 497 | |
|---|
| 498 | |
|---|
| 499 | |
|---|
| 500 | |
|---|
| 501 | void print_tree(char **tree_description, FILE *tree, int *totals) |
|---|
| 502 | { |
|---|
| 503 | int row,col; |
|---|
| 504 | |
|---|
| 505 | fprintf(tree,"\n"); |
|---|
| 506 | |
|---|
| 507 | for(row=1; row<=nseqs-3; row++) { |
|---|
| 508 | fprintf(tree," \n"); |
|---|
| 509 | for(col=1; col<=nseqs; col++) { |
|---|
| 510 | if(tree_description[row][col] == 0) |
|---|
| 511 | fprintf(tree,"*"); |
|---|
| 512 | else |
|---|
| 513 | fprintf(tree,"."); |
|---|
| 514 | } |
|---|
| 515 | if(totals[row] > 0) |
|---|
| 516 | fprintf(tree,"%7d",totals[row]); |
|---|
| 517 | } |
|---|
| 518 | fprintf(tree," \n"); |
|---|
| 519 | for(col=1; col<=nseqs; col++) |
|---|
| 520 | fprintf(tree,"%1d",tree_description[nseqs-2][col]); |
|---|
| 521 | fprintf(tree,"\n"); |
|---|
| 522 | } |
|---|
| 523 | |
|---|
| 524 | |
|---|
| 525 | |
|---|
| 526 | void dna_distance_matrix(FILE *tree) |
|---|
| 527 | { |
|---|
| 528 | int m,n,j,i; |
|---|
| 529 | int res1, res2; |
|---|
| 530 | double p,q,e,a,b,k; |
|---|
| 531 | |
|---|
| 532 | tree_gap_delete(); /* flag positions with gaps (tree_gaps[i] = 1 ) */ |
|---|
| 533 | |
|---|
| 534 | if(verbose) { |
|---|
| 535 | fprintf(tree,"\n"); |
|---|
| 536 | fprintf(tree,"\n DIST = percentage divergence (/100)"); |
|---|
| 537 | fprintf(tree,"\n p = rate of transition (A <-> G; C <-> T)"); |
|---|
| 538 | fprintf(tree,"\n q = rate of transversion"); |
|---|
| 539 | fprintf(tree,"\n Length = number of sites used in comparison"); |
|---|
| 540 | fprintf(tree,"\n"); |
|---|
| 541 | if(tossgaps) { |
|---|
| 542 | fprintf(tree,"\n All sites with gaps (in any sequence) deleted!"); |
|---|
| 543 | fprintf(tree,"\n"); |
|---|
| 544 | } |
|---|
| 545 | if(kimura) { |
|---|
| 546 | fprintf(tree,"\n Distances corrected by Kimura's 2 parameter model:"); |
|---|
| 547 | fprintf(tree,"\n\n Kimura, M. (1980)"); |
|---|
| 548 | fprintf(tree," A simple method for estimating evolutionary "); |
|---|
| 549 | fprintf(tree,"rates of base"); |
|---|
| 550 | fprintf(tree,"\n substitutions through comparative studies of "); |
|---|
| 551 | fprintf(tree,"nucleotide sequences."); |
|---|
| 552 | fprintf(tree,"\n J. Mol. Evol., 16, 111-120."); |
|---|
| 553 | fprintf(tree,"\n\n"); |
|---|
| 554 | } |
|---|
| 555 | } |
|---|
| 556 | |
|---|
| 557 | for(m=1; m<nseqs; ++m) /* for every pair of sequence */ |
|---|
| 558 | for(n=m+1; n<=nseqs; ++n) { |
|---|
| 559 | p = q = e = 0.0; |
|---|
| 560 | tmat[m][n] = tmat[n][m] = 0.0; |
|---|
| 561 | for(i=1; i<=seqlen_array[1]; ++i) { |
|---|
| 562 | j = boot_positions[i]; |
|---|
| 563 | if(tossgaps && (tree_gaps[j] > 0) ) |
|---|
| 564 | goto skip; /* gap position */ |
|---|
| 565 | res1 = seq_array[m][j]; |
|---|
| 566 | res2 = seq_array[n][j]; |
|---|
| 567 | if( (res1 < 1) || (res2 < 1) ) |
|---|
| 568 | goto skip; /* gap in a seq*/ |
|---|
| 569 | e = e + 1.0; |
|---|
| 570 | if(res1 != res2) { |
|---|
| 571 | if(transition(res1,res2)) |
|---|
| 572 | p = p + 1.0; |
|---|
| 573 | else |
|---|
| 574 | q = q + 1.0; |
|---|
| 575 | } |
|---|
| 576 | skip:; |
|---|
| 577 | } |
|---|
| 578 | |
|---|
| 579 | |
|---|
| 580 | /* Kimura's 2 parameter correction for multiple substitutions */ |
|---|
| 581 | |
|---|
| 582 | if(!kimura) { |
|---|
| 583 | k = (p+q)/e; |
|---|
| 584 | if(p > 0.0) |
|---|
| 585 | p = p/e; |
|---|
| 586 | else |
|---|
| 587 | p = 0.0; |
|---|
| 588 | if(q > 0.0) |
|---|
| 589 | q = q/e; |
|---|
| 590 | else |
|---|
| 591 | q = 0.0; |
|---|
| 592 | tmat[m][n] = tmat[n][m] = k; |
|---|
| 593 | if(verbose) /* if screen output */ |
|---|
| 594 | fprintf(tree, |
|---|
| 595 | "%4d vs.%4d: DIST = %7.4f; p = %6.4f; q = %6.4f; length = %6.0f\n" |
|---|
| 596 | ,m,n,k,p,q,e); |
|---|
| 597 | } |
|---|
| 598 | else { |
|---|
| 599 | if(p > 0.0) |
|---|
| 600 | p = p/e; |
|---|
| 601 | else |
|---|
| 602 | p = 0.0; |
|---|
| 603 | if(q > 0.0) |
|---|
| 604 | q = q/e; |
|---|
| 605 | else |
|---|
| 606 | q = 0.0; |
|---|
| 607 | a = 1.0/(1.0-2.0*p-q); |
|---|
| 608 | b = 1.0/(1.0-2.0*q); |
|---|
| 609 | k = 0.5*log(a) + 0.25*log(b); |
|---|
| 610 | tmat[m][n] = tmat[n][m] = k; |
|---|
| 611 | if(verbose) /* if screen output */ |
|---|
| 612 | fprintf(tree, |
|---|
| 613 | "%4d vs.%4d: DIST = %7.4f; p = %6.4f; q = %6.4f; length = %6.0f\n" |
|---|
| 614 | ,m,n,k,p,q,e); |
|---|
| 615 | |
|---|
| 616 | } |
|---|
| 617 | } |
|---|
| 618 | } |
|---|
| 619 | |
|---|
| 620 | |
|---|
| 621 | |
|---|
| 622 | void prot_distance_matrix(FILE *tree) |
|---|
| 623 | { |
|---|
| 624 | int m,n,j,i; |
|---|
| 625 | int res1, res2; |
|---|
| 626 | double p,e,a,b,k; |
|---|
| 627 | |
|---|
| 628 | tree_gap_delete(); /* flag positions with gaps (tree_gaps[i] = 1 ) */ |
|---|
| 629 | |
|---|
| 630 | if(verbose) { |
|---|
| 631 | fprintf(tree,"\n"); |
|---|
| 632 | fprintf(tree,"\n DIST = percentage divergence (/100)"); |
|---|
| 633 | fprintf(tree,"\n Length = number of sites used in comparison"); |
|---|
| 634 | fprintf(tree,"\n\n"); |
|---|
| 635 | if(tossgaps) { |
|---|
| 636 | fprintf(tree,"\n All sites with gaps (in any sequence) deleted"); |
|---|
| 637 | fprintf(tree,"\n"); |
|---|
| 638 | } |
|---|
| 639 | if(kimura) { |
|---|
| 640 | fprintf(tree,"\n Distances corrected by Kimura's empirical method:"); |
|---|
| 641 | fprintf(tree,"\n\n Kimura, M. (1983)"); |
|---|
| 642 | fprintf(tree," The Neutral Theory of Molecular Evolution."); |
|---|
| 643 | fprintf(tree,"\n Cambridge University Press, Cambridge, England."); |
|---|
| 644 | fprintf(tree,"\n\n"); |
|---|
| 645 | } |
|---|
| 646 | } |
|---|
| 647 | |
|---|
| 648 | for(m=1; m<nseqs; ++m) /* for every pair of sequence */ |
|---|
| 649 | for(n=m+1; n<=nseqs; ++n) { |
|---|
| 650 | p = e = 0.0; |
|---|
| 651 | tmat[m][n] = tmat[n][m] = 0.0; |
|---|
| 652 | for(i=1; i<=seqlen_array[1]; ++i) { |
|---|
| 653 | j = boot_positions[i]; |
|---|
| 654 | if(tossgaps && (tree_gaps[j] > 0) ) goto skip; /* gap position */ |
|---|
| 655 | res1 = seq_array[m][j]; |
|---|
| 656 | res2 = seq_array[n][j]; |
|---|
| 657 | if( (res1 < 1) || (res2 < 1) ) goto skip; /* gap in a seq*/ |
|---|
| 658 | e = e + 1.0; |
|---|
| 659 | if(res1 != res2) p = p + 1.0; |
|---|
| 660 | skip:; |
|---|
| 661 | } |
|---|
| 662 | |
|---|
| 663 | if(p <= 0.0) |
|---|
| 664 | k = 0.0; |
|---|
| 665 | else |
|---|
| 666 | k = p/e; |
|---|
| 667 | |
|---|
| 668 | if(kimura) |
|---|
| 669 | if(k > 0.0) k = - log(1.0 - k - (k * k/5.0) ); |
|---|
| 670 | |
|---|
| 671 | tmat[m][n] = tmat[n][m] = k; |
|---|
| 672 | if(verbose) /* if screen output */ |
|---|
| 673 | fprintf(tree, |
|---|
| 674 | "%4d vs.%4d DIST = %6.4f; length = %6.0f\n",m,n,k,e); |
|---|
| 675 | } |
|---|
| 676 | } |
|---|
| 677 | |
|---|
| 678 | |
|---|
