| 1 | /* |
|---|
| 2 | * tridist.c Adachi, J. 1995.02.11 |
|---|
| 3 | * Copyright (C) 1993-1995 J. Adachi & M. Hasegawa, All rights reserved. |
|---|
| 4 | */ |
|---|
| 5 | |
|---|
| 6 | #define MAIN_MODULE 1 |
|---|
| 7 | #include "tridist.h" |
|---|
| 8 | |
|---|
| 9 | |
|---|
| 10 | void |
|---|
| 11 | copyright() |
|---|
| 12 | { |
|---|
| 13 | |
|---|
| 14 | #ifdef SD |
|---|
| 15 | fprintf(stderr, "SDdist %s (%s) ", VERSION, DATE ); |
|---|
| 16 | fprintf(stderr, "Star Division Phylogeny from Distance Matrix\n"); |
|---|
| 17 | #else /* SD */ |
|---|
| 18 | #ifdef NJ |
|---|
| 19 | fprintf(stderr, "NJdist %s (%s) ", VERSION, DATE ); |
|---|
| 20 | fprintf(stderr, "Neighbor Joining Phylogeny from Distance Matrix\n"); |
|---|
| 21 | #else /* NJ */ |
|---|
| 22 | fprintf(stderr, "Tridist %s (%s) ", VERSION, DATE ); |
|---|
| 23 | fprintf(stderr, "Triad Distance Phylogeny from Distance Matrix\n"); |
|---|
| 24 | #endif /* NJ */ |
|---|
| 25 | #endif /* SD */ |
|---|
| 26 | |
|---|
| 27 | fprintf(stderr, "Copyright (C) 1993-1995 J. Adachi & M. Hasegawa. "); |
|---|
| 28 | fprintf(stderr, "All rights reserved.\n"); |
|---|
| 29 | |
|---|
| 30 | #ifdef NJ |
|---|
| 31 | fprintf(stderr, "Ref: N. Saitou & M. Nei 1987."); |
|---|
| 32 | fprintf(stderr, " Molecular Biology and Evolution 4:406-425\n"); |
|---|
| 33 | #endif /* NJ */ |
|---|
| 34 | } |
|---|
| 35 | |
|---|
| 36 | |
|---|
| 37 | void |
|---|
| 38 | usage() |
|---|
| 39 | { |
|---|
| 40 | copyright(); |
|---|
| 41 | fprintf(stderr, "Usage: %s [-%s] distance_matrix_file -h : help\n", |
|---|
| 42 | Prog_name, SWITCHES); |
|---|
| 43 | fprintf(stderr, " or : protml -D sequence_file | %s [-%s]\n", |
|---|
| 44 | Prog_name, SWITCHES); |
|---|
| 45 | } |
|---|
| 46 | |
|---|
| 47 | |
|---|
| 48 | void |
|---|
| 49 | helpinfo() |
|---|
| 50 | { |
|---|
| 51 | copyright(); |
|---|
| 52 | fprintf(stderr, |
|---|
| 53 | "Usage: %s [switches] distance_matrix_file\n",Prog_name); |
|---|
| 54 | fprintf(stderr, "Switches:\n"); |
|---|
| 55 | fprintf(stderr, "-u UPGMA\n"); |
|---|
| 56 | fprintf(stderr, "-w branch length\n"); |
|---|
| 57 | fprintf(stderr, "-l Least squares\n"); |
|---|
| 58 | fprintf(stderr, "-S Sequential input format (PHYLIP)\n"); |
|---|
| 59 | fprintf(stderr, "-o num branch number of Out group \n"); |
|---|
| 60 | fprintf(stderr, "-T str output Tree file name\n"); |
|---|
| 61 | fprintf(stderr, "-t str output Topology file name\n"); |
|---|
| 62 | } |
|---|
| 63 | |
|---|
| 64 | static void prologue(); |
|---|
| 65 | static void distmethod(); |
|---|
| 66 | |
|---|
| 67 | |
|---|
| 68 | main(argc, argv) |
|---|
| 69 | int argc; |
|---|
| 70 | char **argv; |
|---|
| 71 | { |
|---|
| 72 | FILE *distfp, *treefp, *tplgyfp; |
|---|
| 73 | int i, ch; |
|---|
| 74 | char **cpp; |
|---|
| 75 | char *treefname, *tplgyfname; |
|---|
| 76 | extern int Optindex; /* index of next argument */ |
|---|
| 77 | extern char *Optargp; /* pointer to argument of current option */ |
|---|
| 78 | |
|---|
| 79 | Ct0 = time(NULL); |
|---|
| 80 | if((Prog_name = strrchr(argv[0], DIR_CHAR)) != NULL ) |
|---|
| 81 | Prog_name++; |
|---|
| 82 | else |
|---|
| 83 | Prog_name = argv[0]; |
|---|
| 84 | Tplgy_optn = FALSE; |
|---|
| 85 | |
|---|
| 86 | while((ch = mygetopt(argc, argv, SWITCHES)) != -1 ) { |
|---|
| 87 | switch(ch) { |
|---|
| 88 | case 'u': Upgma_optn = TRUE; break; |
|---|
| 89 | case 'l': Least_optn = TRUE; break; |
|---|
| 90 | case 'i': Info_optn = TRUE; break; |
|---|
| 91 | case 'm': cpp = &Optargp; |
|---|
| 92 | Multi_optn = TRUE; |
|---|
| 93 | if (i = strtol(Optargp, cpp, 10)) Numexe = i; |
|---|
| 94 | break; |
|---|
| 95 | case 'v': Verbs_optn = TRUE; break; |
|---|
| 96 | case 'w': Write_optn = TRUE; break; |
|---|
| 97 | case 'S': Seque_optn = TRUE; break; |
|---|
| 98 | case 'o': Outgr_optn = TRUE; |
|---|
| 99 | cpp = &Optargp; |
|---|
| 100 | if (i = strtol(Optargp, cpp, 10)) Outgroup = i-1; |
|---|
| 101 | break; |
|---|
| 102 | case 'T': Tfile_optn = TRUE; |
|---|
| 103 | treefname = new_cvector(strlen(Optargp)+strlen(TREEFEXT)+1); |
|---|
| 104 | *treefname = '\0'; |
|---|
| 105 | strcat(treefname, Optargp); |
|---|
| 106 | strcat(treefname, TREEFEXT); |
|---|
| 107 | break; |
|---|
| 108 | case 't': Tplgy_optn = TRUE; |
|---|
| 109 | tplgyfname = new_cvector(strlen(Optargp)+strlen(TPLGYFEXT)+1); |
|---|
| 110 | *tplgyfname = '\0'; |
|---|
| 111 | strcat(tplgyfname, Optargp); |
|---|
| 112 | strcat(tplgyfname, TPLGYFEXT); |
|---|
| 113 | break; |
|---|
| 114 | case 'z': Debug_optn = TRUE; break; |
|---|
| 115 | case 'Z': Debug = TRUE; break; |
|---|
| 116 | case 'h': |
|---|
| 117 | case 'H': helpinfo(); exit(1); |
|---|
| 118 | default : usage(); exit(1); |
|---|
| 119 | } |
|---|
| 120 | } |
|---|
| 121 | |
|---|
| 122 | if (!Tplgy_optn) { |
|---|
| 123 | Tplgy_optn = TRUE; |
|---|
| 124 | tplgyfname = new_cvector(strlen(TPLGYFILE)+1); |
|---|
| 125 | *tplgyfname = '\0'; |
|---|
| 126 | strcat(tplgyfname, TPLGYFILE); |
|---|
| 127 | } |
|---|
| 128 | |
|---|
| 129 | #ifdef DEBUG |
|---|
| 130 | if (Debug) { |
|---|
| 131 | printf("argc = %d\n",argc); |
|---|
| 132 | for(i = 0; i < argc; i++) printf("argv[%d] = %s\n",i,argv[i]); |
|---|
| 133 | putchar('\n'); |
|---|
| 134 | printf("\nOptindex = %d\n",Optindex); |
|---|
| 135 | printf("Optargp = %s\n",Optargp); |
|---|
| 136 | } |
|---|
| 137 | #endif |
|---|
| 138 | if (Optindex == argc) { |
|---|
| 139 | distfp = stdin; |
|---|
| 140 | } else if (Optindex + 1 == argc) { |
|---|
| 141 | if ((distfp = fopen(argv[Optindex++], "r")) == NULL) { |
|---|
| 142 | fprintf(stderr,"%s: can't open %s\n",Prog_name,argv[--Optindex]); |
|---|
| 143 | exit(1); |
|---|
| 144 | } |
|---|
| 145 | } else { |
|---|
| 146 | fprintf(stderr, "%s: Inconsistent number of file in command line\n", |
|---|
| 147 | Prog_name); |
|---|
| 148 | usage(); |
|---|
| 149 | exit(1); |
|---|
| 150 | } |
|---|
| 151 | |
|---|
| 152 | Cnoexe = 0; |
|---|
| 153 | if (Verbs_optn) copyright(); |
|---|
| 154 | if (!Multi_optn) Numexe = 1; |
|---|
| 155 | for (Cnoexe = 0; Cnoexe < Numexe; Cnoexe++) { |
|---|
| 156 | prologue(distfp, stdout); |
|---|
| 157 | distmethod(distfp, stdout); |
|---|
| 158 | } |
|---|
| 159 | if (distfp != stdin) fclose(distfp); |
|---|
| 160 | |
|---|
| 161 | if (Tfile_optn) { |
|---|
| 162 | if ((treefp = fopen(treefname, "w")) == NULL) { |
|---|
| 163 | fprintf(stderr,"%s: can't open tree file: %s\n", |
|---|
| 164 | Prog_name, treefname); |
|---|
| 165 | exit(1); |
|---|
| 166 | } else { |
|---|
| 167 | fputcphylogeny(treefp, Ctree); |
|---|
| 168 | fclose(treefp); |
|---|
| 169 | } |
|---|
| 170 | } |
|---|
| 171 | if (Tplgy_optn) { |
|---|
| 172 | if ((tplgyfp = fopen(tplgyfname, "w")) == NULL) { |
|---|
| 173 | fprintf(stderr,"%s: can't open topology file: %s\n", |
|---|
| 174 | Prog_name, tplgyfname); |
|---|
| 175 | exit(1); |
|---|
| 176 | } else { |
|---|
| 177 | fputs("1 ", tplgyfp); |
|---|
| 178 | header(tplgyfp, &Numspc, &Comment); |
|---|
| 179 | fputctopology(tplgyfp, Ctree); |
|---|
| 180 | fclose(tplgyfp); |
|---|
| 181 | } |
|---|
| 182 | } |
|---|
| 183 | |
|---|
| 184 | Relia_optn = FALSE; |
|---|
| 185 | Epsfile = new_cvector((unsigned)strlen(EPSFILE) + 1); |
|---|
| 186 | *Epsfile = '\0'; |
|---|
| 187 | strcat(Epsfile, EPSFILE); |
|---|
| 188 | if (1) { |
|---|
| 189 | if ((Epsfp = fopen(Epsfile, "w")) == NULL) { |
|---|
| 190 | fprintf(stderr, |
|---|
| 191 | "%s: can't open eps tree file: %s\n",Prog_name,Epsfile); |
|---|
| 192 | exit(1); |
|---|
| 193 | } else { |
|---|
| 194 | pstree(Epsfp, Ctree); |
|---|
| 195 | fclose(Epsfp); |
|---|
| 196 | } |
|---|
| 197 | } |
|---|
| 198 | |
|---|
| 199 | free_cmatrix(Identif); |
|---|
| 200 | free_cmatrix(Sciname); |
|---|
| 201 | free_cmatrix(Engname); |
|---|
| 202 | free_cvector(Comment); |
|---|
| 203 | return 0; |
|---|
| 204 | } |
|---|
| 205 | |
|---|
| 206 | |
|---|
| 207 | static void |
|---|
| 208 | prologue(ifp, ofp) |
|---|
| 209 | FILE *ifp; |
|---|
| 210 | FILE *ofp; |
|---|
| 211 | { |
|---|
| 212 | |
|---|
| 213 | getsize(ifp, &Numspc, &Comment); |
|---|
| 214 | if (!Multi_optn) header(ofp, &Numspc, &Comment); |
|---|
| 215 | Identif = (char **)malloc((unsigned)Numspc * sizeof(char *)); |
|---|
| 216 | if (Identif == NULL) maerror("in prologue, Identif"); |
|---|
| 217 | Sciname = (char **)malloc((unsigned)Numspc * sizeof(char *)); |
|---|
| 218 | if (Sciname == NULL) maerror("in prologue, Sciname"); |
|---|
| 219 | Engname = (char **)malloc((unsigned)Numspc * sizeof(char *)); |
|---|
| 220 | if (Engname == NULL) maerror("in prologue, Engname"); |
|---|
| 221 | Distanmat = new_dmatrix(Numspc, Numspc); |
|---|
| 222 | if (Seque_optn) |
|---|
| 223 | getdatas(ifp, Identif, Distanmat, Numspc); |
|---|
| 224 | else |
|---|
| 225 | getdata(ifp, Identif, Sciname, Engname, Distanmat, Numspc); |
|---|
| 226 | Maxbrnch = 2 * Numspc - 3; |
|---|
| 227 | Numpair = (Numspc * (Numspc - 1)) / 2; |
|---|
| 228 | } |
|---|
| 229 | |
|---|
| 230 | |
|---|
| 231 | static void |
|---|
| 232 | distmethod(ifp, ofp) |
|---|
| 233 | FILE *ifp; |
|---|
| 234 | FILE *ofp; |
|---|
| 235 | { |
|---|
| 236 | int n; |
|---|
| 237 | |
|---|
| 238 | if (Write_optn && Info_optn) prdistanmat(Distanmat, Numspc); |
|---|
| 239 | if (Least_optn) Distanvec = new_dvector(Numpair); |
|---|
| 240 | if (Least_optn) changedistan(Distanmat, Distanvec, Numspc); |
|---|
| 241 | if (!Cnoexe) Ctree = (Tree *)newtree(Maxbrnch); |
|---|
| 242 | if (Cnoexe) { |
|---|
| 243 | for (n = 0; n < Maxbrnch; n++) { |
|---|
| 244 | Ctree->brnchp[n]->length = 0.0; |
|---|
| 245 | Ctree->brnchp[n]->kinp->length = 0.0; |
|---|
| 246 | } |
|---|
| 247 | } |
|---|
| 248 | getproportion(&Proportion, Distanmat, Numspc); |
|---|
| 249 | distantree(Ctree, Distanmat, Numspc); |
|---|
| 250 | free_dmatrix(Distanmat); |
|---|
| 251 | if (!Multi_optn) prtopology(Ctree); |
|---|
| 252 | if (Least_optn) pathing(Ctree); |
|---|
| 253 | if (Least_optn) Lengths = new_dvector(Numbrnch); |
|---|
| 254 | if (Least_optn) lslength(Ctree, Distanvec, Lengths); |
|---|
| 255 | if (!Multi_optn && Write_optn) resulttree(Ctree); |
|---|
| 256 | if (!Multi_optn && Debug_optn) putchar('\n'); |
|---|
| 257 | if (Debug_optn) fputctopology(stdout, Ctree); |
|---|
| 258 | if (Least_optn) copylengths(Ctree, Lengths, Numbrnch); |
|---|
| 259 | if (Debug_optn) putchar('\n'); |
|---|
| 260 | if (Debug_optn) fputcphylogeny(stdout, Ctree); |
|---|
| 261 | if (Least_optn) free_dvector(Lengths); |
|---|
| 262 | if (Least_optn) free_dvector(Distanvec); |
|---|
| 263 | } /* distmethod */ |
|---|
