source: tags/initial/GDE/MOLPHY/seqstat.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: 12.0 KB
Line 
1/*
2 * seqstat.c   Adachi, J.   1995.03.18
3 * Copyright (C) 1993-1995 J. Adachi & M. Hasegawa, All rights reserved.
4 */
5
6#define STGRAPH 0
7#define RRTEST 0
8#define RRTESTNUM 1000 /* 100000 */
9
10#include "protst.h"
11
12
13void
14seqinsdel(seqchar, numspc, numsite, insdel, numnoinsdel)
15cmatrix seqchar;
16int numspc, numsite;
17ivector insdel;
18int *numnoinsdel;
19{
20        int i, k, ninsdel;
21
22        ninsdel = 0;
23        for (k = 0; k < numsite; k++) {
24                insdel[k] = FALSE;
25                for (i = 0; i < numspc; i++) {
26                        if (seqchar[i][k] == Tpmradix) {
27                                insdel[k] = TRUE;
28                                ninsdel++;
29                                break;
30                        }
31                }
32        }
33        *numnoinsdel = numsite - ninsdel;
34} /* seqinsdel */
35
36
37#if RRTEST
38void
39rrtest(seqchar, numspc, numsite, insdel, numnoinsdel)
40cmatrix seqchar;
41int numspc, numsite;
42ivector insdel;
43int numnoinsdel;
44{
45        int i, j, k, dif, nsite, nsp1, comdiff;
46        double coefrand;
47        cvector seqchi, seqchj;
48        ivector diff;
49        imatrix hist;
50
51        nsp1 = numspc - 1;
52        coefrand = (double)numsite / ((double)RANDOM_MAX + 1.0);
53        diff = new_ivector(numspc);
54        hist = new_imatrix(numspc, numsite);
55        for (i = 0; i < numspc; i++) {
56                for (j = 0; j < numsite; j++) hist[i][j] = 0;
57        }
58        seqchi = seqchar[numspc - 1];
59
60        /* putchar('\n'); */
61        for (i = 0; i < RRTESTNUM; i++) {
62                for (j = 0; j < nsp1; j++) {
63                        seqchj = seqchar[j];
64                        for (k = 0, dif = 0; k < numsite; k++) {
65                                nsite = (int)( coefrand * (double)rand() ); /* RANDOM */
66                                if (!insdel[nsite]) {
67                                        if (seqchi[nsite] != seqchj[nsite]) dif++;
68                                }
69                        }
70                        diff[j] = dif;
71                        hist[j][dif]++;
72                }
73                diff[numspc - 1] = 0;
74                /*
75                printf("%5d", i+1);
76                for (j = 0; j < nsp1; j++) printf(" %4d",diff[j]); putchar('\n');
77                */
78        }
79
80        putchar('\n');
81        for (i = 0; i < nsp1 - 1; i++) {
82                for (j = i + 1; j < nsp1; j++) {
83                        comdiff = 0;
84                        for (k = 0; k < numsite; k++) {
85                                if (hist[i][k] <= hist[j][k]) {
86                                        dif = hist[i][k];
87                                } else {
88                                        dif = hist[j][k];
89                                }
90                                comdiff += dif;
91                        }
92                        printf("%3d %3d %8d %9.5f\n",
93                                i+1, j+1, comdiff, (double)comdiff/RRTESTNUM);
94                }
95        }
96
97        printf("\n%5s", "diff");
98        for (i = 0; i < nsp1; i++) printf(" %4d",i+1); putchar('\n');
99        for (j = 0; j < numsite; j++) {
100                printf("%5d", j);
101                for (i = 0; i < nsp1; i++) printf(" %4d",hist[i][j]); putchar('\n');
102        }
103
104        free_ivector(diff);
105        free_imatrix(hist);
106        exit(1);
107} /* rrtest */
108#endif /*RRTEST*/
109
110
111void
112seqdiff(seqchar, numspc, numsite, insdel, numnoinsdel)
113cmatrix seqchar;
114int numspc, numsite;
115ivector insdel;
116int numnoinsdel;
117{
118        int i, j, k, x, y, z, dif1, dif2, maxk, numk;
119        cvector seqchi, seqchj;
120        imatrix diff;
121#if STGRAPH
122        double v, s, vsd2, ssd2;
123#endif
124
125        diff = new_imatrix(numspc, numspc);
126        for (i = 0; i < numspc-1; i++) {
127                seqchi = seqchar[i];
128                for (j = i+1; j < numspc; j++) {
129                        seqchj = seqchar[j];
130                /*      for (k=0;k<numsite;k++) printf("%2d",seqchj[k]); putchar('\n'); */
131                        for (k = 0, dif1 = 0, dif2 = 0; k < numsite; k++) {
132                                if (!insdel[k]) {
133                                        if ((x = seqchi[k]) != (y = seqchj[k])) {
134#ifndef NUC
135                                                dif1++;
136#else
137                                                z = x + y;
138                                                if ((z == 1) || (z == 5))
139                                                        dif1++;
140                                                else
141                                                        dif2++;
142#endif /* NUC */
143                                        }
144                                }
145                        }
146#ifndef NUC
147                        diff[i][j] = dif1;
148                        diff[j][i] = dif1;
149#else
150                        diff[i][j] = dif1;
151                        diff[j][i] = dif2;
152#endif /* NUC */
153                }
154                diff[i][i] = 0;
155        }
156        diff[numspc-1][numspc-1] = 0;
157 
158        if (!Align_optn) {
159                if (numnoinsdel != numsite)
160                        printf("\nexcluding ins/del sites: %d\n", numnoinsdel);
161                numk = Maxelement;
162                for (k = 0; maxk = k + numk, k < numspc; k += numk) {
163                        if (maxk > numspc) maxk = numspc;
164                        if (Tpmradix == NUMAMI)
165                                printf("\n%4s%6s", "Diff","");
166                        else
167                                printf("\n%4s%-6s", " ","Ts");
168                        for (j = k; j < maxk; j++) printf("%4d", j + 1);
169                        putchar('\n');
170                        if (Tpmradix == NUMAMI)
171                                printf("%4s%6s", "","");
172                        else
173                                printf("%-4s%6s", "Tv","");
174                        for (j = k; j < maxk; j++) printf("%4.3s", Identif[j]);
175                        putchar('\n');
176                        for (i = 0; i < numspc; i++) {
177                                printf("%-4d%-6.6s", i + 1, Identif[i]);
178                                for (j = k; j < maxk; j++) {
179                                        if (i != j)
180                                                printf("%4d", diff[i][j]);
181                                        else
182                                                printf("%4.3s", Identif[i]);
183                                }
184                                putchar('\n');
185                        }
186                }
187        } /* Align_optn */
188
189#if STGRAPH
190        putchar('\n');
191        for (i = 0; i < numspc - 1; i++) {
192                for (j = i + 1; j < numspc; j++) {
193                        v = (double)diff[j][i];
194                        s =     (double)diff[i][j];
195                        vsd2 = 2 * sqrt(v * (1 - v/numsite)) / numsite;
196                        ssd2 = 2 * sqrt(s * (1 - s/numsite)) / numsite;
197                        printf("%9.5f%9.5f%9.5f%9.5f%9.5f%9.5f\n",
198                                v/numsite, v/numsite - vsd2, v/numsite + vsd2,
199                                s/numsite, s/numsite - ssd2, s/numsite + ssd2);
200                }
201        }
202#endif /* STGRAPH */
203
204        free_imatrix(diff);
205
206#if RRTEST
207        rrtest(seqchar, numspc, numsite, insdel, numnoinsdel);
208#endif /* RRTEST */
209
210} /* seqdiff */
211
212
213void
214seqfreq(seqchar, numspc, numsite, insdel, numnoinsdel)
215cmatrix seqchar;
216int numspc, numsite;
217ivector insdel;
218int numnoinsdel;
219{
220        int i, j, k, l, numfreq, numk, maxk;
221        double temp, bias, skew, nid1, nid2;
222        cvector seqchi;
223        imatrix freqs;
224        ivector freq;
225        ivector freqall;
226
227        nid1 = 1.0 / numnoinsdel;
228        nid2 = 1.0 / (numnoinsdel * numnoinsdel);
229        freqs = new_imatrix(numspc, Tpmradix + 1);
230        freqall = new_ivector(Tpmradix + 1);
231        for (j = 0; j < Tpmradix + 1; j++) freqall[j] = 0;
232        for (i = 0; i < numspc; i++) {
233                seqchi = seqchar[i];
234                freq = freqs[i];
235                for (j = 0; j < Tpmradix + 1; j++) freq[j] = 0;
236                for (k = 0; k < numsite; k++) {
237                        if (!insdel[k])
238                                freq[seqchi[k]]++;
239                }
240                for (j = 0; j < Tpmradix + 1; j++) freqall[j] += freq[j];
241        }
242
243        if (Tpmradix == NUMAMI)
244                numfreq = NUMAMI / 2;
245        else
246                numfreq = NUMNUC;
247
248        if (!Align_optn) {
249                printf("\n%-10s", ""); /* Frequencies, Composition */
250                if (Tpmradix == NUMAMI) {
251                        for (j = 0; j < numfreq; j++)
252                                printf("%3s%4s", Cacid1[j], Cacid3[j]);
253                } else {
254                        for (j = 0; j < numfreq; j++) printf("%5s%2s", Cacid1[j], " ");
255                        printf("%7.4s%7.4s%7.4s%7.4s", "A+T ", "G+C ", "Bias", "Skew");
256                }
257                putchar('\n');
258                for (i = 0; i < numspc; i++) {
259                        freq = freqs[i];
260                        printf("%-4d%-6.6s", i + 1, Identif[i]);
261                        for (j = 0; j < numfreq ; j++) {
262                                if (freq[j] != 0)
263                                        printf("%7.3f", (double)freq[j] / numnoinsdel);
264                                else
265                                        printf("%5.1f  ", (double)freq[j] / numnoinsdel);
266                        }
267                        if (Tpmradix == NUMNUC) {
268                                printf("%7.3f", (double)(freq[0]+freq[2]) / numnoinsdel);
269                                printf("%7.3f", (double)(freq[1]+freq[3]) / numnoinsdel);
270                                for (l = 0, bias = 0.0; l < Tpmradix; l++) {
271                                        temp = (freq[l] * nid1 - 0.25);
272                                        bias += temp * temp;
273                                }
274                                bias = bias * 4 / 3;
275                                printf("%7.3f", bias);
276                                skew = fabs(freq[0] - freq[2]) + fabs(freq[1] - freq[3]);
277                                skew /= numnoinsdel;
278                                printf("%7.3f", skew);
279                        }
280                        putchar('\n');
281                }
282                printf("%4s%-6s", " ", "mean");
283                for (j = 0; j < numfreq ; j++) {
284                        if (freqall[j] != 0)
285                                printf("%7.3f", (double)freqall[j] / (numspc * numnoinsdel));
286                        else
287                                printf("%5.1f  ", (double)freqall[j] / (numspc * numnoinsdel));
288                }
289                if (Tpmradix == NUMNUC) {
290                        printf("%7.3f",
291                                (double)(freqall[0]+freqall[2])/(numspc*numnoinsdel));
292                        printf("%7.3f",
293                                (double)(freqall[1]+freqall[3])/(numspc*numnoinsdel));
294                        for (l = 0, bias = 0.0; l < Tpmradix; l++) {
295                                temp = (freqall[l] * nid1 / numspc - 0.25);
296                                bias += temp * temp;
297                        }
298                        bias = bias * 4 / 3;
299                        printf("%7.3f", bias);
300                        skew = fabs(freqall[0] - freqall[2])
301                                 + fabs(freqall[1] - freqall[3]);
302                        skew /= (numspc * numnoinsdel);
303                        printf("%7.3f", skew);
304                }
305                putchar('\n');
306       
307                if (Tpmradix == NUMAMI) {
308                printf("\n%-10s", ""); /* Frequencies, Composition */
309                if (Tpmradix == NUMAMI) {
310                        for (j = numfreq; j < Tpmradix; j++)
311                                printf("%3s%4s", Cacid1[j], Cacid3[j]);
312                } else {
313                        for (j = numfreq; j < Tpmradix; j++)
314                                printf("%5s%2s", Cacid1[j], " ");
315                }
316                putchar('\n');
317                for (i = 0; i < numspc; i++) {
318                        freq = freqs[i];
319                        printf("%-4d%-6.6s", i + 1, Identif[i]);
320                        for (j = numfreq; j < Tpmradix ; j++) {
321                                if (freq[j] != 0)
322                                        printf("%7.3f", (double)freq[j] / numnoinsdel);
323                                else
324                                        printf("%5.1f  ", (double)freq[j] / numnoinsdel);
325                        }
326                        putchar('\n');
327                }
328                printf("%4s%-6s", " ", "mean");
329                for (j = numfreq; j < Tpmradix ; j++) {
330                        if (freqall[j] != 0)
331                                printf("%7.3f", (double)freqall[j] / (numspc * numnoinsdel));
332                        else
333                                printf("%5.1f  ", (double)freqall[j] / (numspc * numnoinsdel));
334                } putchar('\n');
335                }
336        } /* Align_optn */
337
338        if (!Align_optn) {
339                if (Maxspc > 28) putchar('\f');
340                numk = Maxelement;
341                for (k = 0; maxk = k + numk, k < numspc; k += numk) {
342                        if (maxk > numspc) maxk = numspc;
343#if 0
344                        printf("\n%4s%6s", "Bias"," x1e5");
345#else
346                        printf("\n%4s%6s", "Bias"," x1e3");
347#endif
348                        for (j = k; j < maxk; j++) printf("%4d", j + 1); putchar('\n');
349                        printf("%4s%6s", "","");
350                        for (j = k; j < maxk; j++) printf("%4.3s", Identif[j]);
351                        putchar('\n');
352                        for (i = 0; i < numspc; i++) {
353                                printf("%-4d%-6.6s", i + 1, Identif[i]);
354                                for (j = k; j < maxk; j++) {
355                                        if (i != j) {
356#if 0
357                                                for (l = 0, bias = 0.0; l < Tpmradix; l++) {
358                                                        temp = (freqs[i][l] - freqs[j][l]);
359                                                        bias += temp * temp;
360                                                }
361                                                bias = bias * nid2 * 0.5;
362                                                printf("%4.0f", bias * 100000);
363#else
364                                                for (l = 0, bias = 0.0; l < Tpmradix; l++) {
365                                                        bias += abs(freqs[i][l] - freqs[j][l]);
366                                                }
367                                                bias = bias * 0.5 / numnoinsdel;
368                                                printf("%4.0f", bias * 1000);
369#endif
370                                        } else {
371                                                printf("%4.3s", Identif[i]);
372                                        }
373                                }
374                                putchar('\n');
375                        }
376                }
377#if 0
378                printf("\n%d\n", numspc);
379                for (i = 0; i < numspc; i++) {
380                        printf("%s\n", Identif[i]);
381                        for (j = 0; j < numspc; j++) {
382                                        for (l = 0, bias = 0.0; l < Tpmradix; l++) {
383                                                temp = (freqs[i][l] - freqs[j][l]);
384                                                bias += temp * temp;
385                                        }
386                                        bias = bias * nid2 * 0.5;
387                                        printf(" %7.3f", bias * 100000);
388                        }
389                        putchar('\n');
390                }
391#endif
392        } /* Align_optn */
393
394        free_ivector(freqall);
395        free_imatrix(freqs);
396} /* seqfreq */
397
398
399void
400seqtran(seqchar, numspc, numsite, insdel, numnoinsdel)
401cmatrix seqchar;
402int numspc, numsite;
403ivector insdel;
404int numnoinsdel;
405{
406        int i, j, k, numradix, ii, jj;
407        cvector seqchi, seqchj;
408        imatrix trans;
409        ivector tranall;
410
411        trans = new_imatrix(Tpmradix + 1, Tpmradix + 1);
412        tranall = new_ivector(Tpmradix + 1);
413        for (i = 0; i < Tpmradix + 1; i++) {
414                for (j = 0; j < Tpmradix + 1; j++) trans[i][j] = 0;
415                tranall[j] = 0;
416        }
417        for (i = 0; i < numspc - 1; i++) {
418                seqchi = seqchar[i];
419                for (j = i + 1; j < numspc; j++) {
420                        seqchj = seqchar[j];
421                        for (k = 0; k < numsite; k++) {
422                                if (!insdel[k]) {
423                                        trans[seqchi[k]][seqchj[k]]++;
424                                        trans[seqchj[k]][seqchi[k]]++;
425                                }
426                        }
427                }
428        }
429
430        if (Tpmradix == NUMAMI)
431                numradix = NUMAMI / 2;
432        else
433                numradix = NUMNUC;
434
435        printf("\n%7s", " ");
436        if (Tpmradix == NUMAMI) {
437                for (j = 0; j < numradix; j++) printf("%3s%4s", Cacid1[j], Cacid3[j]);
438        } else {
439                for (j = 0; j < numradix; j++) printf("%5s%2s", Cacid1[j], " ");
440        }
441        putchar('\n');
442        for (i = 0; i < Tpmradix; i++) {
443                if (Tpmradix == NUMAMI) {
444                        printf("%3s%4s", Cacid1[i], Cacid3[i]);
445                } else {
446                        printf("%5s%2s", Cacid1[i], " ");
447                }
448                for (j = 0; j < numradix ; j++) {
449                        printf("%7d", trans[i][j]);
450                }
451                putchar('\n');
452        }
453
454        if (Tpmradix == NUMAMI) {
455        printf("\n%7s", " ");
456        if (Tpmradix == NUMAMI) {
457                for (j = numradix; j < Tpmradix; j++)
458                        printf("%3s%4s", Cacid1[j], Cacid3[j]);
459        } else {
460                for (j = numradix; j < Tpmradix; j++)
461                        printf("%5s%2s", Cacid1[j], " ");
462        }
463        putchar('\n');
464        for (i = 0; i < Tpmradix; i++) {
465                if (Tpmradix == NUMAMI) {
466                        printf("%3s%4s", Cacid1[i], Cacid3[i]);
467                } else {
468                        printf("%5s%2s", Cacid1[i], " ");
469                }
470                for (j = numradix; j < Tpmradix ; j++) {
471                        printf("%7d", trans[i][j]);
472                }
473                putchar('\n');
474        }
475        }
476
477#if 0
478        for (i = 0; i < numspc - 1; i++) {
479                seqchi = seqchar[i];
480                for (j = i + 1; j < numspc; j++) {
481                        seqchj = seqchar[j];
482                        for (ii = 0; ii < Tpmradix + 1; ii++) {
483                                for (jj = 0; jj < Tpmradix + 1; jj++) trans[ii][jj] = 0;
484                        }
485                        for (k = 0; k < numsite; k++) {
486                                if (!insdel[k]) {
487                                        trans[seqchi[k]][seqchj[k]]++;
488                                        trans[seqchj[k]][seqchi[k]]++;
489                                }
490                        }
491
492                        printf("\n%-3d %-3d  %s  %s\n", i+1, j+1, Identif[i], Identif[j]);
493                        for (jj = 0; jj < Tpmradix; jj++) printf("%4s", Cacid1[jj]);
494                        putchar('\n');
495                        for (ii = 0; ii < Tpmradix; ii++) {
496                                printf("%1s", Cacid1[ii]);
497                                printf("%3d", trans[ii][0]);
498                                for (jj = 1; jj < Tpmradix ; jj++) {
499                                        printf("%4d", trans[ii][jj]);
500                                }
501                                putchar('\n');
502                        }
503
504                }
505        }
506#endif
507
508        free_ivector(tranall);
509        free_imatrix(trans);
510} /* seqtran */
Note: See TracBrowser for help on using the repository browser.