source: tags/initial/GDE/PHYLIP/protdist.c

Last change on this file was 2, checked in by oldcode, 24 years ago

Initial revision
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 38.9 KB
Line 
1#include "phylip.h"
2
3/*
4 * version 3.52c. (c) Copyright 1993 by Joseph Felsenstein. Written by Joseph
5 * Felsenstein, Akiko Fuseki, Sean Lamont, and Andrew Keeffe. Permission is
6 * granted to copy and use this program provided no fee is charged for it and
7 * provided that this copyright notice is not removed.
8 */
9
10#define nmlngth         10      /* number of characters in species name     */
11#define epsilon         0.000001/* a small number */
12
13#define ibmpc0          false
14#define ansi0           true
15#define vt520           false
16
17
18typedef long   *steparray;
19typedef enum {
20        ala, arg, asn, asp, cys, gln, glu, gly, his, ileu, leu, lys, met, phe, pro,
21        ser, thr, trp, tyr, val, del, stop, asx, glx, unk, quest
22}               aas;
23typedef enum {
24        universal, ciliate, mito, vertmito, flymito, yeastmito
25}               codetype;
26typedef enum {
27        chemical, hall, george
28}               cattype;
29
30typedef double  matrix[20][20];
31
32
33FILE           *infile, *outfile;
34long            spp, chars, datasets, ith;
35/*
36 * spp = number of species chars = number of sites in actual sequences
37 */
38double          freqa, freqc, freqg, freqt, ttratio, xi, xv, ease, fracchange;
39boolean         weights, printdata, progress, mulsets, interleaved, ibmpc,
40                vt52, ansi, basesequal, usepam, kimura, firstset;
41codetype        whichcode;
42cattype         whichcat;
43steparray       weight;
44Char          **nayme;
45aas           **node;
46aas             trans[4][4][4];
47double          pi[20];
48long            cat[(long) val - (long) ala + 1], numaa[(long) val - (long) ala + 1];
49double          eig[20];
50matrix          prob, eigvecs;
51double        **d;
52
53/* Local variables for makedists, propagated globally for c version: */
54double          tt, p, dp, d2p, q, elambdat;
55
56
57static double   pameigs[] = {-0.022091252, -0.019297602, 0.000004760, -0.017477817,
58        -0.016575549, -0.015504543, -0.002112213, -0.002685727,
59        -0.002976402, -0.013440755, -0.012926992, -0.004293227,
60        -0.005356688, -0.011064786, -0.010480731, -0.008760449,
61        -0.007142318, -0.007381851, -0.007806557, -0.008127024
62
63};
64static double   pamprobs[20][20] =
65{
66        {-0.01522976, -0.00746819, -0.13934468, 0.11755315, -0.00212101,
67                0.01558456, -0.07408235, -0.00322387, 0.01375826, 0.00448826,
68                0.00154174, 0.02013313, -0.00159183, -0.00069275, -0.00399898,
69        0.08414055, -0.01188178, -0.00029870, 0.00220371, 0.00042546},
70        {-0.07765582, -0.00712634, -0.03683209, -0.08065755, -0.00462872,
71                -0.03791039, 0.10642147, -0.00912185, 0.01436308, -0.00133243,
72                0.00166346, 0.00624657, -0.00003363, -0.00128729, -0.00690319,
73        0.17442028, -0.05373336, -0.00078751, -0.00038151, 0.01382718},
74        {-0.08810973, -0.04081786, -0.04066004, -0.04736004, -0.03275406,
75                -0.03761164, -0.05047487, -0.09086213, -0.03269598, -0.03558015,
76                -0.08407966, -0.07970977, -0.01504743, -0.04011920, -0.05182232,
77        -0.07026991, -0.05846931, -0.01016998, -0.03047472, -0.06280511},
78        {0.02513756, -0.00578333, 0.09865453, 0.01322314, -0.00310665,
79                0.05880899, -0.09252443, -0.02986539, -0.03127460, 0.01007539,
80                -0.00360119, -0.01995024, 0.00094940, -0.00145868, -0.01388816,
81        0.11358341, -0.12127513, -0.00054696, -0.00055627, 0.00417284},
82        {0.16517316, -0.00254742, -0.03318745, -0.01984173, 0.00031890,
83                -0.02817810, 0.02661678, -0.01761215, 0.01665112, 0.10513343,
84                -0.00545026, 0.01827470, -0.00207616, -0.00763758, -0.01322808,
85        -0.02202576, -0.07434204, 0.00020593, 0.00119979, -0.10827873},
86        {0.16088826, 0.00056313, -0.02579303, -0.00319655, 0.00037228,
87                -0.03193150, 0.01655305, -0.03028640, 0.01367746, -0.11248153,
88                0.00778371, 0.02675579, 0.00243718, 0.00895470, -0.01729803,
89        -0.02686964, -0.08262584, 0.00011794, -0.00225134, 0.09415650},
90        {-0.01739295, 0.00572017, -0.00712592, -0.01100922, -0.00870113,
91                -0.00663461, -0.01153857, -0.02248432, -0.00382264, -0.00358612,
92                -0.00139345, -0.00971460, -0.00133312, 0.01927783, -0.01053838,
93        -0.00911362, -0.01010908, 0.09417598, 0.01763850, -0.00955454},
94        {0.01728888, 0.01344211, 0.01200836, 0.01857259, -0.17088517,
95                0.01457592, 0.01997839, 0.02844884, 0.00839403, 0.00196862,
96                0.01391984, 0.03270465, 0.00347173, -0.01940984, 0.01233979,
97        0.00542887, 0.01008836, 0.00126491, -0.02863042, 0.00449764},
98        {-0.02881366, -0.02184155, -0.01566086, -0.02593764, -0.04050907,
99                -0.01539603, -0.02576729, -0.05089606, -0.00597430, 0.02181643,
100                0.09835597, -0.04040940, 0.00873512, 0.12139434, -0.02427882,
101        -0.02945238, -0.01566867, -0.01606503, 0.09475319, 0.02238670},
102        {0.04080274, -0.02869626, -0.05191093, -0.08435843, 0.00021141,
103                0.13043842, 0.00871530, 0.00496058, -0.02797641, -0.00636933,
104                0.02243277, 0.03640362, -0.05735517, 0.00196918, -0.02218934,
105        -0.00608972, 0.02872922, 0.00047619, 0.00151285, 0.00883489},
106        {-0.02623824, 0.00331152, 0.03640692, 0.04260231, -0.00038223,
107                -0.07480340, -0.01022492, -0.00426473, 0.01448116, 0.01456847,
108                0.05786680, 0.03368691, -0.10126924, -0.00147454, 0.01275395,
109        0.00017574, -0.01585206, -0.00015767, -0.00231848, 0.02310137},
110        {-0.00846258, -0.01508106, -0.01967505, -0.02772004, 0.01248253,
111                -0.01331243, -0.02569382, -0.04461524, -0.02207075, 0.04663443,
112                0.19347923, -0.02745691, 0.02288515, -0.04883849, -0.01084597,
113        -0.01947187, -0.00081675, 0.00516540, -0.07815919, 0.08035585},
114        {-0.06553111, 0.09756831, 0.00524326, -0.00885098, 0.00756653,
115                0.02783099, -0.00427042, -0.16680359, 0.03951331, -0.00490540,
116                0.01719610, 0.15018204, 0.00882722, -0.00423197, -0.01919217,
117        -0.02963619, -0.01831342, -0.00524338, 0.00011379, -0.02566864},
118        {-0.07494341, -0.11348850, 0.00241343, -0.00803016, 0.00492438,
119                0.00711909, -0.00829147, 0.05793337, 0.02734209, 0.02059759,
120                -0.02770280, 0.14128338, 0.01532479, 0.00364307, 0.05968116,
121        -0.06497960, -0.08113941, 0.00319445, -0.00104222, 0.03553497},
122        {0.05948223, -0.08959930, 0.03269977, -0.03272374, -0.00365667,
123                -0.03423294, -0.06418925, -0.05902138, 0.05746317, -0.02580596,
124                0.01259572, 0.05848832, 0.00672666, 0.00233355, -0.05145149,
125        0.07348503, 0.11427955, 0.00142592, -0.01030651, -0.04862799},
126        {-0.01606880, 0.05200845, -0.01212967, -0.06824429, -0.00234304,
127                0.01094203, -0.07375538, 0.08808629, 0.12394822, 0.02231351,
128                -0.03608265, -0.06978045, -0.00618360, 0.00274747, -0.01921876,
129        -0.01541969, -0.02223856, -0.00107603, -0.01251777, 0.05412534},
130        {0.01688843, 0.05784728, -0.02256966, -0.07072251, -0.00422551,
131                -0.06261233, -0.08502830, 0.08925346, -0.08529597, 0.01519343,
132                -0.05008258, 0.10931873, 0.00521033, 0.02593305, -0.00717855,
133        0.02291527, 0.02527388, -0.00266188, -0.00871160, 0.02708135},
134        {-0.04233344, 0.00076379, 0.01571257, 0.04003092, 0.00901468,
135                0.00670577, 0.03459487, 0.12420216, -0.00067366, -0.01515094,
136                0.05306642, 0.04338407, 0.00511287, 0.01036639, -0.17867462,
137        -0.02289440, -0.03213205, 0.00017924, -0.01187362, -0.03933874},
138        {0.01284817, -0.01685622, 0.00724363, 0.01687952, -0.00882070,
139                -0.00555957, 0.01676246, -0.05560456, -0.00966893, 0.06197684,
140                -0.09058758, 0.00880607, 0.00108629, -0.08308956, -0.08056832,
141        -0.00413297, 0.02973107, 0.00092948, 0.07010111, 0.13007418},
142        {0.00700223, -0.01347574, 0.00691332, 0.03122905, 0.00310308,
143                0.00946862, 0.03455040, -0.06712536, -0.00304506, 0.04267941,
144                -0.10422292, -0.01127831, -0.00549798, 0.11680505, -0.03352701,
145        -0.00084536, 0.01631369, 0.00095063, -0.09570217, 0.06480321}
146};
147
148openfile(fp, filename, mode, application, perm)
149        FILE          **fp;
150        char           *filename;
151        char           *mode;
152        char           *application;
153        char           *perm;
154{
155        FILE           *of;
156        char            file[100];
157        strcpy(file, filename);
158        while (1) {
159                of = fopen(file, mode);
160                if (of)
161                        break;
162                else {
163                        switch (*mode) {
164                        case 'r':
165                                printf("%s:  can't read %s\n", application, file);
166                                file[0] = '\0';
167                                while (file[0] == '\0') {
168                                        printf("Please enter a new filename>");
169                                        gets(file);
170                                }
171                                break;
172                        case 'w':
173                                printf("%s: can't write %s\n", application, file);
174                                file[0] = '\0';
175                                while (file[0] == '\0') {
176                                        printf("Please enter a new filename>");
177                                        gets(file);
178                                }
179                                break;
180                        }
181                }
182        }
183        *fp = of;
184        if (perm != NULL)
185                strcpy(perm, file);
186}
187
188
189
190void 
191uppercase(ch)
192        Char           *ch;
193{
194        (*ch) = (isupper(*ch) ? (*ch) : toupper(*ch));
195}                               /* uppercase */
196
197
198void 
199inputnumbers()
200{
201        /* input the numbers of species and of characters */
202        long            i;
203
204        fscanf(infile, "%ld%ld", &spp, &chars);
205
206        if (printdata)
207                fprintf(outfile, "%2ld species, %3ld  sites\n\n", spp, chars);
208        node = (aas **) Malloc(spp * sizeof(aas *));
209        if (firstset) {
210                for (i = 0; i < spp; i++)
211                        node[i] = (aas *) Malloc(chars * sizeof(aas));
212        }
213        weight = (steparray) Malloc(chars * sizeof(long));
214        d = (double **) Malloc(spp * sizeof(double *));
215        nayme = (char **) Malloc(spp * sizeof(char *));
216
217        for (i = 0; i < spp; ++i) {
218                d[i] = (double *) Malloc(spp * sizeof(double));
219                nayme[i] = (char *) malloc(nmlngth * sizeof(char));
220        }
221
222}                               /* inputnumbers */
223
224void 
225getoptions()
226{
227        /* interactively set options */
228        Char            ch;
229        Char            in[100];
230        boolean         done, done1;
231
232        if (printdata)
233                fprintf(outfile, "\nProtein distance algorithm, version 3.5\n\n");
234        putchar('\n');
235        weights = false;
236        printdata = false;
237        progress = true;
238        interleaved = true;
239        ttratio = 2.0;
240        whichcode = universal;
241        whichcat = george;
242        basesequal = true;
243        freqa = 0.25;
244        freqc = 0.25;
245        freqg = 0.25;
246        freqt = 0.25;
247        usepam = true;
248        kimura = false;
249        ease = 0.457;
250        do {
251                printf(ansi ? "\033[2J\033[H" :
252                       vt52 ? "\033E\033H" : "\n");
253                printf("\nProtein distance algorithm, version %s\n\n", VERSION);
254                printf("Settings for this run:\n");
255                printf("  P  Use PAM, Kimura or categories model?  %s\n",
256                       usepam ? "Dayhoff PAM matrix" :
257                       kimura ? "Kimura formula" : "Categories model");
258                if (!(usepam || kimura)) {
259                        printf("  C               Use which genetic code?  %s\n",
260                               (whichcode == universal) ? "Universal" :
261                               (whichcode == ciliate) ? "Ciliate" :
262                           (whichcode == mito) ? "Universal mitochondrial" :
263                               (whichcode == vertmito) ? "Vertebrate mitochondrial" :
264                            (whichcode == flymito) ? "Fly mitochondrial\n" :
265                               (whichcode == yeastmito) ? "Yeast mitochondrial" : "");
266                        printf("  A  Which categorization of amino acids?  %s\n",
267                               (whichcat == chemical) ? "Chemical" :
268                               (whichcat == george) ? "George/Hunt/Barker" : "Hall");
269
270                        printf("  E      Prob change category (1.0=easy):%8.4f\n", ease);
271                        printf("  T        Transition/transversion ratio:%7.3f\n", ttratio);
272                        printf("  F                     Base Frequencies:");
273                        if (basesequal)
274                                printf("  Equal\n");
275                        else
276                                printf("%7.3f%6.3f%6.3f%6.3f\n", freqa, freqc, freqg, freqt);
277                }
278                printf("  M           Analyze multiple data sets?");
279                if (mulsets)
280                        printf("  Yes, %2ld sets\n", datasets);
281                else
282                        printf("  No\n");
283                printf("  I          Input sequences interleaved?  %s\n",
284                       (interleaved ? "Yes" : "No, sequential"));
285                printf("  0   Terminal type (IBM PC, VT52, ANSI)?  %s\n",
286                       ibmpc ? "IBM PC" :
287                       ansi ? "ANSI" :
288                       vt52 ? "VT52" : "(none)");
289                printf("  1    Print out the data at start of run  %s\n",
290                       (printdata ? "Yes" : "No"));
291                printf("  2  Print indications of progress of run  %s\n",
292                       progress ? "Yes" : "No");
293                printf("\nAre these settings correct? (type Y or the letter for one to change)\n");
294                in[0] = '\0';
295                gets(in);
296                ch = in[0];
297                if (ch == '\n')
298                        ch = ' ';
299                uppercase(&ch);
300                done = (ch == 'Y');
301                if (!done) {
302                        if (ch == 'C' || ch == 'A' || ch == 'P' || ch == 'E' || ch == 'T' ||
303                            ch == 'F' || ch == 'M' || ch == 'I' || ch == '1' || ch == '2' ||
304                            ch == '0') {
305                                switch (ch) {
306
307                                case 'C':
308                                        printf("Which genetic code?\n");
309                                        printf(" type         for\n\n");
310                                        printf("   U           Universal\n");
311                                        printf("   M           Mitochondrial\n");
312                                        printf("   V           Vertebrate mitochondrial\n");
313                                        printf("   F           Fly mitochondrial\n");
314                                        printf("   Y           Yeast mitochondrial\n\n");
315                                        do {
316                                                printf("type U, M, V, F, or Y\n");
317                                                scanf("%c%*[^\n]", &ch);
318                                                getchar();
319                                                if (ch == '\n')
320                                                        ch = ' ';
321                                                uppercase(&ch);
322                                        } while (ch != 'U' && ch != 'M' && ch != 'V' && ch != 'F' && ch != 'Y');
323                                        switch (ch) {
324
325                                        case 'U':
326                                                whichcode = universal;
327                                                break;
328
329                                        case 'M':
330                                                whichcode = mito;
331                                                break;
332
333                                        case 'V':
334                                                whichcode = vertmito;
335                                                break;
336
337                                        case 'F':
338                                                whichcode = flymito;
339                                                break;
340
341                                        case 'Y':
342                                                whichcode = yeastmito;
343                                                break;
344                                        }
345                                        break;
346
347                                case 'A':
348                                        printf(
349                                               "Which of these categorizations of amino acids do you want to use:\n\n");
350                                        printf(
351                                               " all have groups: (Glu Gln Asp Asn), (Lys Arg His), (Phe Tyr Trp)\n");
352                                        printf(" plus:\n");
353                                        printf("George/Hunt/Barker:");
354                                        printf(" (Cys), (Met   Val  Leu  Ileu), (Gly  Ala  Ser  Thr    Pro)\n");
355                                        printf("Chemical:          ");
356                                        printf(" (Cys   Met), (Val  Leu  Ileu    Gly  Ala  Ser  Thr), (Pro)\n");
357                                        printf("Hall:              ");
358                                        printf(" (Cys), (Met   Val  Leu  Ileu), (Gly  Ala  Ser  Thr), (Pro)\n\n");
359                                        printf("Which do you want to use (type C, H, or G)\n");
360                                        do {
361                                                scanf("%c%*[^\n]", &ch);
362                                                getchar();
363                                                if (ch == '\n')
364                                                        ch = ' ';
365                                                uppercase(&ch);
366                                        } while (ch != 'C' && ch != 'H' && ch != 'G');
367                                        switch (ch) {
368
369                                        case 'C':
370                                                whichcat = chemical;
371                                                break;
372
373                                        case 'H':
374                                                whichcat = hall;
375                                                break;
376
377                                        case 'G':
378                                                whichcat = george;
379                                                break;
380                                        }
381                                        break;
382
383                                case 'P':
384                                        if (usepam) {
385                                                usepam = false;
386                                                kimura = true;
387                                        } else {
388                                                if (kimura)
389                                                        kimura = false;
390                                                else
391                                                        usepam = true;
392                                        }
393                                        break;
394
395                                case 'E':
396                                        printf("Ease of changing category of amino acid?\n");
397                                        do {
398                                                printf(" (1.0 if no difficulty of changing,\n");
399                                                printf(" less if less easy. Can't be negative\n");
400                                                scanf("%lf%*[^\n]", &ease);
401                                                getchar();
402                                        } while (ease > 1.0 || ease < 0.0);
403                                        break;
404
405                                case 'T':
406                                        do {
407                                                printf("Transition/transversion ratio?\n");
408                                                scanf("%lf%*[^\n]", &ttratio);
409                                                getchar();
410                                        } while (ttratio < 0.0);
411                                        break;
412
413                                case 'F':
414                                        do {
415                                                basesequal = false;
416                                                printf("Frequencies of bases A,C,G,T ?\n");
417                                                scanf("%lf%lf%lf%lf%*[^\n]", &freqa, &freqc, &freqg, &freqt);
418                                                getchar();
419                                                if (fabs(freqa + freqc + freqg + freqt - 1.0) >= 1.0e-3)
420                                                        printf("FREQUENCIES MUST SUM TO 1\n");
421                                        } while (fabs(freqa + freqc + freqg + freqt - 1.0) >= 1.0e-3);
422                                        break;
423
424                                case 'M':
425                                        mulsets = !mulsets;
426                                        if (mulsets) {
427                                                done1 = false;
428                                                do {
429                                                        printf("How many data sets?\n");
430                                                        scanf("%ld%*[^\n]", &datasets);
431                                                        getchar();
432                                                        done1 = (datasets >= 1);
433                                                        if (!done1)
434                                                                printf("BAD NUMBER OF DATA SETS:  it must be greater than 1\n");
435                                                } while (done1 != true);
436                                        }
437                                        break;
438
439                                case 'I':
440                                        interleaved = !interleaved;
441                                        break;
442
443                                case '0':
444                                        if (ibmpc) {
445                                                ibmpc = false;
446                                                vt52 = true;
447                                        } else {
448                                                if (vt52) {
449                                                        vt52 = false;
450                                                        ansi = true;
451                                                } else if (ansi)
452                                                        ansi = false;
453                                                else
454                                                        ibmpc = true;
455                                        }
456                                        break;
457
458                                case '1':
459                                        printdata = !printdata;
460                                        break;
461
462                                case '2':
463                                        progress = !progress;
464                                        break;
465                                }
466                        } else
467                                printf("Not a possible option!\n");
468                }
469        } while (!done);
470}                               /* getoptions */
471
472void 
473transition()
474{
475        /* calculations related to transition-transversion ratio */
476        double          aa, bb, freqr, freqy, freqgr, freqty;
477
478        freqr = freqa + freqg;
479        freqy = freqc + freqt;
480        freqgr = freqg / freqr;
481        freqty = freqt / freqy;
482        aa = ttratio * freqr * freqy - freqa * freqg - freqc * freqt;
483        bb = freqa * freqgr + freqc * freqty;
484        xi = aa / (aa + bb);
485        xv = 1.0 - xi;
486        if (xi <= 0.0 && xi >= -epsilon)
487                xi = 0.0;
488        if (xi < 0.0) {
489                printf("THIS TRANSITION-TRANSVERSION RATIO IS IMPOSSIBLE WITH");
490                printf(" THESE BASE FREQUENCIES\n");
491                exit(-1);
492        }
493}                               /* transition */
494
495
496void 
497doinit()
498{
499        /* initializes variables */
500        inputnumbers();
501        getoptions();
502        transition();
503}                               /* doinit */
504
505
506
507void 
508inputweights()
509{
510        /* input the character weights, 0-9 and A-Z for weights 10 - 35 */
511        Char            ch;
512        long            i;
513
514        for (i = 1; i < nmlngth; i++) {
515                ch = getc(infile);
516                if (ch == '\n')
517                        ch = ' ';
518        }
519        for (i = 0; i < chars; i++) {
520                do {
521                        if (eoln(infile)) {
522                                fscanf(infile, "%*[^\n]");
523                                getc(infile);
524                        }
525                        ch = getc(infile);
526                        if (ch == '\n')
527                                ch = ' ';
528                } while (ch == ' ');
529                weight[i] = 1;
530                if (isdigit(ch))
531                        weight[i] = ch - '0';
532                else if (isalpha(ch)) {
533                        uppercase(&ch);
534                        if (ch >= 'A' && ch <= 'I')
535                                weight[i] = ch - 55;
536                        else if (ch >= 'J' && ch <= 'R')
537                                weight[i] = ch - 55;
538                        else
539                                weight[i] = ch - 55;
540                } else {
541                        printf("BAD WEIGHT CHARACTER: %c\n", ch);
542                        exit(-1);
543                }
544        }
545        fscanf(infile, "%*[^\n]");
546        getc(infile);
547        weights = true;
548}                               /* inputweights */
549
550void 
551printweights()
552{
553        /* print out the weights of sites */
554        long            i, j, k;
555
556        fprintf(outfile, "    Sites are weighted as follows:\n");
557        fprintf(outfile, "        ");
558        for (i = 0; i <= 9; i++)
559                fprintf(outfile, "%3ld", i);
560        fprintf(outfile, "\n     *---------------------------------\n");
561        for (j = 0; j <= (chars / 10); j++) {
562                fprintf(outfile, "%5ld!  ", j * 10);
563                for (i = 0; i <= 9; i++) {
564                        k = j * 10 + i;
565                        if (k > 0 && k <= chars)
566                                fprintf(outfile, "%3ld", weight[k - 1]);
567                        else
568                                fprintf(outfile, "   ");
569                }
570                putc('\n', outfile);
571        }
572        putc('\n', outfile);
573}                               /* printweights */
574
575void 
576inputoptions()
577{
578        /* input the information on the options */
579        Char            ch;
580        long            extranum, i, cursp, curchs;
581
582        if (!firstset) {
583                if (eoln(infile)) {
584                        fscanf(infile, "%*[^\n]");
585                        getc(infile);
586                }
587                fscanf(infile, "%ld%ld", &cursp, &curchs);
588                if (cursp != spp) {
589                        printf("\nERROR: INCONSISTENT NUMBER OF SPECIES IN DATA SET %4ld\n",
590                               ith + 1);
591                        exit(-1);
592                }
593                chars = curchs;
594        }
595        extranum = 0;
596        while (!(eoln(infile))) {
597                ch = getc(infile);
598                if (ch == '\n')
599                        ch = ' ';
600                uppercase(&ch);
601                if (ch == 'W')
602                        extranum++;
603                else if (ch != ' ') {
604                        printf("BAD OPTION CHARACTER: %c\n", ch);
605                        exit(-1);
606                }
607        }
608        fscanf(infile, "%*[^\n]");
609        getc(infile);
610        for (i = 0; i < chars; i++)
611                weight[i] = 1;
612        for (i = 1; i <= extranum; i++) {
613                ch = getc(infile);
614                if (ch == '\n')
615                        ch = ' ';
616                uppercase(&ch);
617                if (ch == 'W')
618                        inputweights();
619                else {
620                        printf("ERROR: INCORRECT AUXILIARY OPTIONS LINE WHICH STARTS WITH %c\n",
621                               ch);
622                        exit(-1);
623                }
624        }
625        if (weights)
626                printweights();
627}                               /* inputoptions */
628
629void 
630inputdata()
631{
632        /* input the names and sequences for each species */
633        long            i, j, k, l, aasread, aasnew;
634        Char            charstate;
635        boolean         allread, done;
636        aas             aa;     /* temporary amino acid for input */
637
638        if (progress)
639                putchar('\n');
640        j = nmlngth + (chars + (chars - 1) / 10) / 2 - 5;
641        if (j < nmlngth - 1)
642                j = nmlngth - 1;
643        if (j > 37)
644                j = 37;
645        if (printdata) {
646                fprintf(outfile, "Name");
647                for (i = 1; i <= j; i++)
648                        putc(' ', outfile);
649                fprintf(outfile, "Sequences\n");
650                fprintf(outfile, "----");
651                for (i = 1; i <= j; i++)
652                        putc(' ', outfile);
653                fprintf(outfile, "---------\n\n");
654        }
655        aasread = 0;
656        allread = false;
657        while (!(allread)) {
658                allread = true;
659                if (eoln(infile)) {
660                        fscanf(infile, "%*[^\n]");
661                        getc(infile);
662                }
663                i = 1;
664                while (i <= spp) {
665                        if ((interleaved && aasread == 0) || !interleaved) {
666                                for (j = 0; j < nmlngth; j++) {
667                                        nayme[i - 1][j] = getc(infile);
668                                        if (nayme[i - 1][j] == '\n')
669                                                nayme[i - 1][j] = ' ';
670                                        if (eof(infile) || eoln(infile)) {
671                                                printf("ERROR: END-OF-LINE OR END-OF-FILE");
672                                                printf(" IN THE MIDDLE OF A SPECIES NAME\n");
673                                                exit(-1);
674                                        }
675                                }
676                        }
677                        if (interleaved)
678                                j = aasread;
679                        else
680                                j = 0;
681                        done = false;
682                        while (((!done) && (!(eoln(infile) | eof(infile))))) {
683                                if (interleaved)
684                                        done = true;
685                                while (((j < chars) & (!(eoln(infile) | eof(infile))))) {
686                                        charstate = getc(infile);
687                                        if (charstate == '\n')
688                                                charstate = ' ';
689                                        if (charstate == ' ' || (charstate >= '0' && charstate <= '9'))
690                                                continue;
691                                        uppercase(&charstate);
692                                        if ((!isalpha(charstate) && charstate != '.' && charstate != '?' &&
693                                             charstate != '-' && charstate != '*') || charstate == 'J' ||
694                                            charstate == 'O' || charstate == 'U') {
695                                                printf("WARNING -- BAD AMINO ACID:%c AT POSITION%5ld OF SPECIES %3ld\n",
696                                                       charstate, j + 1, i);
697                                                exit(-1);
698                                        }
699                                        j++;
700                                        if (charstate == '.') {
701                                                node[i - 1][j - 1] = node[0][j - 1];
702                                                continue;
703                                        }
704                                        switch (charstate) {
705
706                                        case 'A':
707                                                aa = ala;
708                                                break;
709
710                                        case 'B':
711                                                aa = asx;
712                                                break;
713
714                                        case 'C':
715                                                aa = cys;
716                                                break;
717
718                                        case 'D':
719                                                aa = asp;
720                                                break;
721
722                                        case 'E':
723                                                aa = glu;
724                                                break;
725
726                                        case 'F':
727                                                aa = phe;
728                                                break;
729
730                                        case 'G':
731                                                aa = gly;
732                                                break;
733
734                                        case 'H':
735                                                aa = his;
736                                                break;
737
738                                        case 'I':
739                                                aa = ileu;
740                                                break;
741
742                                        case 'K':
743                                                aa = lys;
744                                                break;
745
746                                        case 'L':
747                                                aa = leu;
748                                                break;
749
750                                        case 'M':
751                                                aa = met;
752                                                break;
753
754                                        case 'N':
755                                                aa = asn;
756                                                break;
757
758                                        case 'P':
759                                                aa = pro;
760                                                break;
761
762                                        case 'Q':
763                                                aa = gln;
764                                                break;
765
766                                        case 'R':
767                                                aa = arg;
768                                                break;
769
770                                        case 'S':
771                                                aa = ser;
772                                                break;
773
774                                        case 'T':
775                                                aa = thr;
776                                                break;
777
778                                        case 'V':
779                                                aa = val;
780                                                break;
781
782                                        case 'W':
783                                                aa = trp;
784                                                break;
785
786                                        case 'X':
787                                                aa = unk;
788                                                break;
789
790                                        case 'Y':
791                                                aa = tyr;
792                                                break;
793
794                                        case 'Z':
795                                                aa = glx;
796                                                break;
797
798                                        case '*':
799                                                aa = stop;
800                                                break;
801
802                                        case '?':
803                                                aa = quest;
804                                                break;
805
806                                        case '-':
807                                                aa = del;
808                                                break;
809                                        }
810                                        node[i - 1][j - 1] = aa;
811                                }
812                                if (interleaved)
813                                        continue;
814                                if (j < chars) {
815                                        fscanf(infile, "%*[^\n]");
816                                        getc(infile);
817                                } else if (j == chars)
818                                        done = true;
819                        }
820                        if (interleaved && i == 1)
821                                aasnew = j;
822                        fscanf(infile, "%*[^\n]");
823                        getc(infile);
824                        if ((interleaved && j != aasnew) || ((!interleaved) && j != chars)) {
825                                printf("ERROR: SEQUENCES OUT OF ALIGNMENT\n");
826                                exit(-1);
827                        }
828                        i++;
829                }
830                if (interleaved) {
831                        aasread = aasnew;
832                        allread = (aasread == chars);
833                } else
834                        allread = (i > spp);
835        }
836        if (printdata) {
837                for (i = 1; i <= ((chars - 1) / 60 + 1); i++) {
838                        for (j = 1; j <= spp; j++) {
839                                for (k = 0; k < nmlngth; k++)
840                                        putc(nayme[j - 1][k], outfile);
841                                fprintf(outfile, "   ");
842                                l = i * 60;
843                                if (l > chars)
844                                        l = chars;
845                                for (k = (i - 1) * 60 + 1; k <= l; k++) {
846                                        if (j > 1 && node[j - 1][k - 1] == node[0][k - 1])
847                                                charstate = '.';
848                                        else {
849                                                switch (node[j - 1][k - 1]) {
850
851                                                case ala:
852                                                        charstate = 'A';
853                                                        break;
854
855                                                case asx:
856                                                        charstate = 'B';
857                                                        break;
858
859                                                case cys:
860                                                        charstate = 'C';
861                                                        break;
862
863                                                case asp:
864                                                        charstate = 'D';
865                                                        break;
866
867                                                case glu:
868                                                        charstate = 'E';
869                                                        break;
870
871                                                case phe:
872                                                        charstate = 'F';
873                                                        break;
874
875                                                case gly:
876                                                        charstate = 'G';
877                                                        break;
878
879                                                case his:
880                                                        charstate = 'H';
881                                                        break;
882
883                                                case ileu:
884                                                        charstate = 'I';
885                                                        break;
886
887                                                case lys:
888                                                        charstate = 'K';
889                                                        break;
890
891                                                case leu:
892                                                        charstate = 'L';
893                                                        break;
894
895                                                case met:
896                                                        charstate = 'M';
897                                                        break;
898
899                                                case asn:
900                                                        charstate = 'N';
901                                                        break;
902
903                                                case pro:
904                                                        charstate = 'P';
905                                                        break;
906
907                                                case gln:
908                                                        charstate = 'Q';
909                                                        break;
910
911                                                case arg:
912                                                        charstate = 'R';
913                                                        break;
914
915                                                case ser:
916                                                        charstate = 'S';
917                                                        break;
918
919                                                case thr:
920                                                        charstate = 'T';
921                                                        break;
922
923                                                case val:
924                                                        charstate = 'V';
925                                                        break;
926
927                                                case trp:
928                                                        charstate = 'W';
929                                                        break;
930
931                                                case tyr:
932                                                        charstate = 'Y';
933                                                        break;
934
935                                                case glx:
936                                                        charstate = 'Z';
937                                                        break;
938
939                                                case del:
940                                                        charstate = '-';
941                                                        break;
942
943                                                case stop:
944                                                        charstate = '*';
945                                                        break;
946
947                                                case unk:
948                                                        charstate = 'X';
949                                                        break;
950
951                                                case quest:
952                                                        charstate = '?';
953                                                        break;
954                                                }
955                                        }
956                                        putc(charstate, outfile);
957                                        if (k % 10 == 0 && k % 60 != 0)
958                                                putc(' ', outfile);
959                                }
960                                putc('\n', outfile);
961                        }
962                        putc('\n', outfile);
963                }
964                putc('\n', outfile);
965        }
966        if (printdata)
967                putc('\n', outfile);
968}                               /* inputdata */
969
970
971void 
972doinput()
973{
974        /* reads the input data */
975        inputoptions();
976        inputdata();
977}                               /* doinput */
978
979
980void 
981code()
982{
983        /* make up table of the code 1 = u, 2 = c, 3 = a, 4 = g */
984        long            n;
985        aas             b;
986
987        trans[0][0][0] = phe;
988        trans[0][0][1] = phe;
989        trans[0][0][2] = leu;
990        trans[0][0][3] = leu;
991        trans[0][1][0] = ser;
992        trans[0][1][1] = ser;
993        trans[0][1][2] = ser;
994        trans[0][1][3] = ser;
995        trans[0][2][0] = tyr;
996        trans[0][2][1] = tyr;
997        trans[0][2][2] = stop;
998        trans[0][2][3] = stop;
999        trans[0][3][0] = cys;
1000        trans[0][3][1] = cys;
1001        trans[0][3][2] = stop;
1002        trans[0][3][3] = trp;
1003        trans[1][0][0] = leu;
1004        trans[1][0][1] = leu;
1005        trans[1][0][2] = leu;
1006        trans[1][0][3] = leu;
1007        trans[1][1][0] = pro;
1008        trans[1][1][1] = pro;
1009        trans[1][1][2] = pro;
1010        trans[1][1][3] = pro;
1011        trans[1][2][0] = his;
1012        trans[1][2][1] = his;
1013        trans[1][2][2] = gln;
1014        trans[1][2][3] = gln;
1015        trans[1][3][0] = arg;
1016        trans[1][3][1] = arg;
1017        trans[1][3][2] = arg;
1018        trans[1][3][3] = arg;
1019        trans[2][0][0] = ileu;
1020        trans[2][0][1] = ileu;
1021        trans[2][0][2] = ileu;
1022        trans[2][0][3] = met;
1023        trans[2][1][0] = thr;
1024        trans[2][1][1] = thr;
1025        trans[2][1][2] = thr;
1026        trans[2][1][3] = thr;
1027        trans[2][2][0] = asn;
1028        trans[2][2][1] = asn;
1029        trans[2][2][2] = lys;
1030        trans[2][2][3] = lys;
1031        trans[2][3][0] = ser;
1032        trans[2][3][1] = ser;
1033        trans[2][3][2] = arg;
1034        trans[2][3][3] = arg;
1035        trans[3][0][0] = val;
1036        trans[3][0][1] = val;
1037        trans[3][0][2] = val;
1038        trans[3][0][3] = val;
1039        trans[3][1][0] = ala;
1040        trans[3][1][1] = ala;
1041        trans[3][1][2] = ala;
1042        trans[3][1][3] = ala;
1043        trans[3][2][0] = asp;
1044        trans[3][2][1] = asp;
1045        trans[3][2][2] = glu;
1046        trans[3][2][3] = glu;
1047        trans[3][3][0] = gly;
1048        trans[3][3][1] = gly;
1049        trans[3][3][2] = gly;
1050        trans[3][3][3] = gly;
1051        if (whichcode == mito)
1052                trans[0][3][2] = trp;
1053        if (whichcode == vertmito) {
1054                trans[0][3][2] = trp;
1055                trans[2][3][2] = stop;
1056                trans[2][3][3] = stop;
1057                trans[2][0][2] = met;
1058        }
1059        if (whichcode == flymito) {
1060                trans[0][3][2] = trp;
1061                trans[2][0][2] = met;
1062                trans[2][3][2] = ser;
1063        }
1064        if (whichcode == yeastmito) {
1065                trans[0][3][2] = trp;
1066                trans[1][0][2] = thr;
1067                trans[2][0][2] = met;
1068        }
1069        n = 0;
1070        for (b = ala; (long) b <= (long) val; b = (aas) ((long) b + 1)) {
1071                n++;
1072                numaa[(long) b - (long) ala] = n;
1073        }
1074}                               /* code */
1075
1076
1077Static Void
1078cats()
1079{
1080        /* define categories of amino acids */
1081        aas             b;
1082
1083        /* fundamental subgroups */
1084        cat[(long) cys - (long) ala] = 1;
1085        cat[(long) met - (long) ala] = 2;
1086        cat[(long) val - (long) ala] = 3;
1087        cat[(long) leu - (long) ala] = 3;
1088        cat[(long) ileu - (long) ala] = 3;
1089        cat[(long) gly - (long) ala] = 4;
1090        cat[0] = 4;
1091        cat[(long) ser - (long) ala] = 4;
1092        cat[(long) thr - (long) ala] = 4;
1093        cat[(long) pro - (long) ala] = 5;
1094        cat[(long) phe - (long) ala] = 6;
1095        cat[(long) tyr - (long) ala] = 6;
1096        cat[(long) trp - (long) ala] = 6;
1097        cat[(long) glu - (long) ala] = 7;
1098        cat[(long) gln - (long) ala] = 7;
1099        cat[(long) asp - (long) ala] = 7;
1100        cat[(long) asn - (long) ala] = 7;
1101        cat[(long) lys - (long) ala] = 8;
1102        cat[(long) arg - (long) ala] = 8;
1103        cat[(long) his - (long) ala] = 8;
1104        if (whichcat == george) {
1105                /*
1106                 * George, Hunt and Barker: sulfhydryl, small hydrophobic,
1107                 * small hydrophilic, aromatic, acid/acid-amide/hydrophilic,
1108                 * basic
1109                 */
1110                for (b = ala; (long) b <= (long) val; b = (aas) ((long) b + 1)) {
1111                        if (cat[(long) b - (long) ala] == 3)
1112                                cat[(long) b - (long) ala] = 2;
1113                        if (cat[(long) b - (long) ala] == 5)
1114                                cat[(long) b - (long) ala] = 4;
1115                }
1116        }
1117        if (whichcat == chemical) {
1118                /*
1119                 * Conn and Stumpf:  monoamino, aliphatic, heterocyclic,
1120                 * aromatic, dicarboxylic, basic
1121                 */
1122                for (b = ala; (long) b <= (long) val; b = (aas) ((long) b + 1)) {
1123                        if (cat[(long) b - (long) ala] == 2)
1124                                cat[(long) b - (long) ala] = 1;
1125                        if (cat[(long) b - (long) ala] == 4)
1126                                cat[(long) b - (long) ala] = 3;
1127                }
1128        }
1129        /* Ben Hall's personal opinion */
1130        if (whichcat != hall)
1131                return;
1132        for (b = ala; (long) b <= (long) val; b = (aas) ((long) b + 1)) {
1133                if (cat[(long) b - (long) ala] == 3)
1134                        cat[(long) b - (long) ala] = 2;
1135        }
1136}                               /* cats */
1137
1138
1139void 
1140maketrans()
1141{
1142        /*
1143         * Make up transition probability matrix from code and category
1144         * tables
1145         */
1146        long            i, j, k, m, n, s, nb1, nb2;
1147        double          x, sum;
1148        long            sub[3], newsub[3];
1149        double          f[4], g[4];
1150        aas             b1, b2;
1151        double          TEMP, TEMP1, TEMP2, TEMP3;
1152
1153        for (i = 0; i <= 19; i++) {
1154                pi[i] = 0.0;
1155                for (j = 0; j <= 19; j++)
1156                        prob[i][j] = 0.0;
1157        }
1158        f[0] = freqt;
1159        f[1] = freqc;
1160        f[2] = freqa;
1161        f[3] = freqg;
1162        g[0] = freqc + freqt;
1163        g[1] = freqc + freqt;
1164        g[2] = freqa + freqg;
1165        g[3] = freqa + freqg;
1166        TEMP = f[0];
1167        TEMP1 = f[1];
1168        TEMP2 = f[2];
1169        TEMP3 = f[3];
1170        fracchange = xi * (2 * f[0] * f[1] / g[0] + 2 * f[2] * f[3] / g[2]) +
1171                xv * (1 - TEMP * TEMP - TEMP1 * TEMP1 - TEMP2 * TEMP2 - TEMP3 * TEMP3);
1172        for (i = 0; i <= 3; i++) {
1173                for (j = 0; j <= 3; j++) {
1174                        for (k = 0; k <= 3; k++) {
1175                                if (trans[i][j][k] != stop)
1176                                        sum += f[i] * f[j] * f[k];
1177                        }
1178                }
1179        }
1180        for (i = 0; i <= 3; i++) {
1181                sub[0] = i + 1;
1182                for (j = 0; j <= 3; j++) {
1183                        sub[1] = j + 1;
1184                        for (k = 0; k <= 3; k++) {
1185                                sub[2] = k + 1;
1186                                b1 = trans[i][j][k];
1187                                for (m = 0; m <= 2; m++) {
1188                                        s = sub[m];
1189                                        for (n = 1; n <= 4; n++) {
1190                                                memcpy(newsub, sub, sizeof(long) * 3L);
1191                                                newsub[m] = n;
1192                                                x = f[i] * f[j] * f[k] / (3.0 * sum);
1193                                                if (((s == 1 || s == 2) && (n == 3 || n == 4)) ||
1194                                                    ((n == 1 || n == 2) && (s == 3 || s == 4)))
1195                                                        x *= xv * f[n - 1];
1196                                                else
1197                                                        x *= xi * f[n - 1] / g[n - 1] + xv * f[n - 1];
1198                                                b2 = trans[newsub[0] - 1][newsub[1] - 1][newsub[2] - 1];
1199                                                if (b1 != stop) {
1200                                                        nb1 = numaa[(long) b1 - (long) ala];
1201                                                        pi[nb1 - 1] += x;
1202                                                        if (b2 != stop) {
1203                                                                nb2 = numaa[(long) b2 - (long) ala];
1204                                                                if (cat[(long) b1 - (long) ala] != cat[(long) b2 - (long) ala]) {
1205                                                                        prob[nb1 - 1][nb2 - 1] += x * ease;
1206                                                                        prob[nb1 - 1][nb1 - 1] += x * (1.0 - ease);
1207                                                                } else
1208                                                                        prob[nb1 - 1][nb2 - 1] += x;
1209                                                        } else
1210                                                                prob[nb1 - 1][nb1 - 1] += x;
1211                                                }
1212                                        }
1213                                }
1214                        }
1215                }
1216        }
1217        for (i = 0; i <= 19; i++)
1218                prob[i][i] -= pi[i];
1219        for (i = 0; i <= 19; i++) {
1220                for (j = 0; j <= 19; j++)
1221                        prob[i][j] /= sqrt(pi[i] * pi[j]);
1222        }
1223        /* computes pi^(1/2)*B*pi^(-1/2)  */
1224}                               /* maketrans */
1225
1226
1227
1228void 
1229givens(a, i, j, n, ctheta, stheta, left)
1230        double          (*a)[20];
1231        long            i, j, n;
1232        double          ctheta, stheta;
1233        boolean         left;
1234{
1235        /* Givens transform at i,j for 1..n with angle theta */
1236        long            k;
1237        double          d;
1238
1239        for (k = 0; k < n; k++) {
1240                if (left) {
1241                        d = ctheta * a[i - 1][k] + stheta * a[j - 1][k];
1242                        a[j - 1][k] = ctheta * a[j - 1][k] - stheta * a[i - 1][k];
1243                        a[i - 1][k] = d;
1244                } else {
1245                        d = ctheta * a[k][i - 1] + stheta * a[k][j - 1];
1246                        a[k][j - 1] = ctheta * a[k][j - 1] - stheta * a[k][i - 1];
1247                        a[k][i - 1] = d;
1248                }
1249        }
1250}                               /* givens */
1251
1252void 
1253coeffs(x, y, c, s, accuracy)
1254        double          x, y, *c, *s;
1255        double          accuracy;
1256{
1257        /* compute cosine and sine of theta */
1258        double          root;
1259
1260        root = sqrt(x * x + y * y);
1261        if (root < accuracy) {
1262                *c = 1.0;
1263                *s = 0.0;
1264        } else {
1265                *c = x / root;
1266                *s = y / root;
1267        }
1268}                               /* coeffs */
1269
1270void 
1271tridiag(a, n, accuracy)
1272        double          (*a)[20];
1273        long            n;
1274        double          accuracy;
1275{
1276        /* Givens tridiagonalization */
1277        long            i, j;
1278        double          s, c;
1279
1280        for (i = 2; i < n; i++) {
1281                for (j = i + 1; j <= n; j++) {
1282                        coeffs(a[i - 2][i - 1], a[i - 2][j - 1], &c, &s, accuracy);
1283                        givens(a, i, j, n, c, s, true);
1284                        givens(a, i, j, n, c, s, false);
1285                        givens(eigvecs, i, j, n, c, s, true);
1286                }
1287        }
1288}                               /* tridiag */
1289
1290void 
1291shiftqr(a, n, accuracy)
1292        double          (*a)[20];
1293        long            n;
1294        double          accuracy;
1295{
1296        /* QR eigenvalue-finder */
1297        long            i, j;
1298        double          approx, s, c, d, TEMP, TEMP1;
1299
1300        for (i = n; i >= 2; i--) {
1301                do {
1302                        TEMP = a[i - 2][i - 2] - a[i - 1][i - 1];
1303                        TEMP1 = a[i - 1][i - 2];
1304                        d = sqrt(TEMP * TEMP + TEMP1 * TEMP1);
1305                        approx = a[i - 2][i - 2] + a[i - 1][i - 1];
1306                        if (a[i - 1][i - 1] < a[i - 2][i - 2])
1307                                approx = (approx - d) / 2.0;
1308                        else
1309                                approx = (approx + d) / 2.0;
1310                        for (j = 0; j < i; j++)
1311                                a[j][j] -= approx;
1312                        for (j = 1; j < i; j++) {
1313                                coeffs(a[j - 1][j - 1], a[j][j - 1], &c, &s, accuracy);
1314                                givens(a, j, j + 1, i, c, s, true);
1315                                givens(a, j, j + 1, i, c, s, false);
1316                                givens(eigvecs, j, j + 1, n, c, s, true);
1317                        }
1318                        for (j = 0; j < i; j++)
1319                                a[j][j] += approx;
1320                } while (fabs(a[i - 1][i - 2]) > accuracy);
1321        }
1322}                               /* shiftqr */
1323
1324
1325void 
1326qreigen(prob, n)
1327        double          (*prob)[20];
1328        long            n;
1329{
1330        /* QR eigenvector/eigenvalue method for symmetric matrix */
1331        double          accuracy;
1332        long            i, j;
1333
1334        accuracy = 1.0e-6;
1335        for (i = 0; i < n; i++) {
1336                for (j = 0; j < n; j++)
1337                        eigvecs[i][j] = 0.0;
1338                eigvecs[i][i] = 1.0;
1339        }
1340        tridiag(prob, n, accuracy);
1341        shiftqr(prob, n, accuracy);
1342        for (i = 0; i < n; i++)
1343                eig[i] = prob[i][i];
1344        for (i = 0; i <= 19; i++) {
1345                for (j = 0; j <= 19; j++)
1346                        prob[i][j] = sqrt(pi[j]) * eigvecs[i][j];
1347        }
1348        /* prob[i][j] is the value of U' times pi^(1/2) */
1349}                               /* qreigen */
1350
1351
1352void 
1353pameigen()
1354{
1355        /* eigenanalysis for PAM matrix, precomputed */
1356        memcpy(prob, pamprobs, sizeof(pamprobs));
1357        memcpy(eig, pameigs, sizeof(pameigs));
1358        fracchange = 0.01;
1359}                               /* pameigen */
1360
1361
1362
1363void 
1364predict(nb1, nb2)
1365        long            nb1, nb2;
1366{
1367        /* make contribution to prediction of this aa pair */
1368        long            m;
1369        double          TEMP;
1370
1371        for (m = 0; m <= 19; m++) {
1372                elambdat = exp(tt * eig[m]);
1373                q = prob[m][nb1 - 1] * prob[m][nb2 - 1] * elambdat;
1374                p += q;
1375                dp += eig[m] * q;
1376                TEMP = eig[m];
1377                d2p += TEMP * TEMP * q;
1378        }
1379}                               /* predict */
1380
1381
1382void 
1383makedists()
1384{
1385        /* compute the distances */
1386        long            i, j, k, m, n, iterations, nb1, nb2;
1387        double          delta, lnlike, slope, curv;
1388        boolean         neginfinity, inf;
1389        aas             b1, b2;
1390
1391        fprintf(outfile, "%5ld\n", spp);
1392        if (progress)
1393                printf("Computing distances:\n");
1394        for (i = 1; i <= spp; i++) {
1395                if (progress)
1396                        printf("  ");
1397                if (progress) {
1398                        for (j = 0; j < nmlngth; j++)
1399                                putchar(nayme[i - 1][j]);
1400                }
1401                if (progress)
1402                        printf("   ");
1403                d[i - 1][i - 1] = 0.0;
1404                for (j = 0; j <= i - 2; j++) {
1405                        if (!kimura) {
1406                                if (usepam)
1407                                        tt = 10.0;
1408                                else
1409                                        tt = 1.0;
1410                                delta = tt / 2.0;
1411                                iterations = 0;
1412                                inf = false;
1413                                do {
1414                                        lnlike = 0.0;
1415                                        slope = 0.0;
1416                                        curv = 0.0;
1417                                        neginfinity = false;
1418                                        for (k = 0; k < chars; k++) {
1419                                                b1 = node[i - 1][k];
1420                                                b2 = node[j][k];
1421                                                if (b1 != stop && b1 != del && b1 != quest && b1 != unk &&
1422                                                    b2 != stop && b2 != del && b2 != quest && b2 != unk) {
1423                                                        p = 0.0;
1424                                                        dp = 0.0;
1425                                                        d2p = 0.0;
1426                                                        nb1 = numaa[(long) b1 - (long) ala];
1427                                                        nb2 = numaa[(long) b2 - (long) ala];
1428                                                        if (b1 != asx && b1 != glx && b2 != asx && b2 != glx)
1429                                                                predict(nb1, nb2);
1430                                                        else {
1431                                                                if (b1 == asx) {
1432                                                                        if (b2 == asx) {
1433                                                                                predict(3L, 3L);
1434                                                                                predict(3L, 4L);
1435                                                                                predict(4L, 3L);
1436                                                                                predict(4L, 4L);
1437                                                                        } else {
1438                                                                                if (b2 == glx) {
1439                                                                                        predict(3L, 6L);
1440                                                                                        predict(3L, 7L);
1441                                                                                        predict(4L, 6L);
1442                                                                                        predict(4L, 7L);
1443                                                                                } else {
1444                                                                                        predict(3L, nb2);
1445                                                                                        predict(4L, nb2);
1446                                                                                }
1447                                                                        }
1448                                                                } else {
1449                                                                        if (b1 == glx) {
1450                                                                                if (b2 == asx) {
1451                                                                                        predict(6L, 3L);
1452                                                                                        predict(6L, 4L);
1453                                                                                        predict(7L, 3L);
1454                                                                                        predict(7L, 4L);
1455                                                                                } else {
1456                                                                                        if (b2 == glx) {
1457                                                                                                predict(6L, 6L);
1458                                                                                                predict(6L, 7L);
1459                                                                                                predict(7L, 6L);
1460                                                                                                predict(7L, 7L);
1461                                                                                        } else {
1462                                                                                                predict(6L, nb2);
1463                                                                                                predict(7L, nb2);
1464                                                                                        }
1465                                                                                }
1466                                                                        } else {
1467                                                                                if (b2 == asx) {
1468                                                                                        predict(nb1, 3L);
1469                                                                                        predict(nb1, 4L);
1470                                                                                        predict(nb1, 3L);
1471                                                                                        predict(nb1, 4L);
1472                                                                                } else if (b2 == glx) {
1473                                                                                        predict(nb1, 6L);
1474                                                                                        predict(nb1, 7L);
1475                                                                                        predict(nb1, 6L);
1476                                                                                        predict(nb1, 7L);
1477                                                                                }
1478                                                                        }
1479                                                                }
1480                                                        }
1481                                                        if (p <= 0.0)
1482                                                                neginfinity = true;
1483                                                        else {
1484                                                                lnlike += log(p);
1485                                                                slope += dp / p;
1486                                                                curv += d2p / p - dp * dp / (p * p);
1487                                                        }
1488                                                }
1489                                        }
1490                                        iterations++;
1491                                        if (!neginfinity) {
1492                                                if (curv < 0.0) {
1493                                                        tt -= slope / curv;
1494                                                        if (tt > 10000.0) {
1495                                                                printf("\nWARNING: INFINITE DISTANCE BETWEEN SPECIES %ld AND %ld; -1.0 WAS WRITTEN\n", i, j);
1496                                                                tt = -1.0 / fracchange;
1497                                                                inf = true;
1498                                                                iterations = 20;
1499                                                        }
1500                                                } else {
1501                                                        if ((slope > 0.0 && delta < 0.0) || (slope < 0.0 && delta > 0.0))
1502                                                                delta /= -2;
1503                                                        tt += delta;
1504                                                }
1505                                        } else {
1506                                                delta /= -2;
1507                                                tt += delta;
1508                                        }
1509                                        if (tt < epsilon && !inf)
1510                                                tt = epsilon;
1511                                } while (iterations != 20);
1512                        } else {
1513                                m = 0;
1514                                n = 0;
1515                                for (k = 0; k < chars; k++) {
1516                                        b1 = node[i - 1][k];
1517                                        b2 = node[j][k];
1518                                        if ((long) b1 <= (long) val && (long) b2 <= (long) val) {
1519                                                if (b1 == b2)
1520                                                        m++;
1521                                                n++;
1522                                        }
1523                                }
1524                                p = 1 - (double) m / n;
1525                                dp = 1.0 - p - 0.2 * p * p;
1526                                if (dp < 0.0) {
1527                                        printf(
1528                                               "\nDISTANCE BETWEEN SEQUENCES %3ld AND %3ld IS TOO LARGE FOR KIMURA FORMULA\n",
1529                                               i, j + 1);
1530                                        tt = -1.0;
1531                                } else
1532                                        tt = -log(dp);
1533                        }
1534                        d[i - 1][j] = fracchange * tt;
1535                        d[j][i - 1] = d[i - 1][j];
1536                        if (progress)
1537                                putchar('.');
1538                }
1539                if (progress)
1540                        putchar('\n');
1541        }
1542        for (i = 0; i < spp; i++) {
1543                for (j = 0; j < nmlngth; j++)
1544                        putc(nayme[i][j], outfile);
1545                fprintf(outfile, "   ");
1546                for (j = 0; j < spp; j++)
1547                        fprintf(outfile, "%9.5f", d[i][j]);
1548                putc('\n', outfile);
1549        }
1550        if (progress)
1551                printf("\nOutput written to output file\n\n");
1552}                               /* makedists */
1553
1554
1555main(argc, argv)
1556        int             argc;
1557        Char           *argv[];
1558{                               /* ML Protein distances by PAM or categories
1559                                 * model */
1560        char            infilename[100], outfilename[100];
1561#ifdef MAC
1562        macsetup("Protdist", "");
1563#endif
1564        openfile(&infile, INFILE, "r", argv[0], infilename);
1565        openfile(&outfile, OUTFILE, "w", argv[0], outfilename);
1566
1567        ibmpc = ibmpc0;
1568        ansi = ansi0;
1569        vt52 = vt520;
1570        mulsets = false;
1571        datasets = 1;
1572        firstset = true;
1573        doinit();
1574        if (!kimura)
1575                code();
1576        if (!(usepam || kimura)) {
1577                cats();
1578                maketrans();
1579                qreigen(prob, 20L);
1580        } else {
1581                if (kimura)
1582                        fracchange = 1.0;
1583                else
1584                        pameigen();
1585        }
1586        for (ith = 1; ith <= datasets; ith++) {
1587                doinput();
1588                if (ith == 1)
1589                        firstset = false;
1590                if ((datasets > 1) && progress)
1591                        printf("\nData set # %ld:\n\n", ith);
1592                makedists();
1593        }
1594        FClose(outfile);
1595        FClose(infile);
1596#ifdef MAC
1597        fixmacfile(outfilename);
1598#endif
1599        exit(0);
1600}                               /* Protein distances */
1601
1602int 
1603eof(f)
1604        FILE           *f;
1605{
1606        register int    ch;
1607
1608        if (feof(f))
1609                return 1;
1610        ch = getc(f);
1611        if (ch == EOF)
1612                return 1;
1613        ungetc(ch, f);
1614        return 0;
1615}
1616
1617
1618int 
1619eoln(f)
1620        FILE           *f;
1621{
1622        register int    ch;
1623
1624        ch = getc(f);
1625        if (ch == EOF)
1626                return 1;
1627        ungetc(ch, f);
1628        return (ch == '\n');
1629}
1630
1631void 
1632memerror()
1633{
1634        printf("Error allocating memory\n");
1635        exit(-1);
1636}
1637
1638MALLOCRETURN   *
1639mymalloc(x)
1640        long            x;
1641{
1642        MALLOCRETURN   *mem;
1643        mem = (MALLOCRETURN *) malloc(x);
1644        if (!mem)
1645                memerror();
1646        else
1647                return (MALLOCRETURN *) mem;
1648}
Note: See TracBrowser for help on using the repository browser.