source: trunk/GDE/SUPPORT/count.c

Last change on this file was 19480, checked in by westram, 17 months ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 6.6 KB
Line 
1/*
2*       Copyright 1991 Steven Smith at the Harvard Genome Lab.
3*       All rights reserved.
4*/
5#include <math.h>
6#include "Flatio.c"
7
8#define FALSE 0
9#define TRUE 1
10
11#define JUKES 0
12#define OLSEN 1
13#define NONE 2
14
15#define Min(a,b) (a)<(b)?(a):(b)
16
17int width,start,jump,usecase,sim,correction;
18int tbl,numseq,special;
19char argtyp[255],argval[255];
20
21float   acwt=1.0, agwt=1.0, auwt=1.0, ucwt=1.0, ugwt=1.0, gcwt=1.0;
22
23float dist[200][200];
24
25struct data_format data[10000];
26float parta[200], partc[200], partg[200], partu[200], setdist();
27
28void getarg(av,ndx,atype,aval)
29char **av,atype[],aval[];
30int ndx;
31{
32        int i,j;
33        char c;
34        for(j=0;(c=av[ndx][j])!=' ' && c!= '=' && c!= '\0';j++)
35                atype[j]=c;
36        if (c=='=')
37        {
38                atype[j++] = c;
39                atype[j] = '\0';
40        }
41        else
42        {
43                atype[j] = '\0';
44                j++;
45        }
46
47        if(c=='=')
48        {
49                for(i=0;(c=av[ndx][j]) != '\0' && c!= ' ';i++,j++)
50                        aval[i] = c;
51                aval[i] = '\0';
52        }
53}
54
55void SetPart()
56{
57        int a,c,g,u,tot,i,j;
58        char nuc;
59
60        for(j=0;j<numseq;j++)
61        {
62                a=0;
63                c=0;
64                g=0;
65                u=0;
66                tot=0;
67
68                for(i=0;i<data[j].length;i++)
69                {
70                        nuc=data[j].nuc[i] | 32;
71                        switch (nuc)
72                        {
73                        case 'a':
74                                a++;
75                                tot++;
76                                break;
77
78                        case 'c':
79                                c++;
80                                tot++;
81                                break;
82
83                        case 'g':
84                                g++;
85                                tot++;
86                                break;
87
88                        case 'u':
89                                u++;
90                                tot++;
91                                break;
92
93                        case 't':
94                                u++;
95                                tot++;
96                                break;
97                        };
98                }
99                parta[j] = (float)a / (float)tot;
100                partc[j] = (float)c / (float)tot;
101                partg[j] = (float)g / (float)tot;
102                partu[j] = (float)u / (float)tot;
103        }
104}
105
106void Compare(a,b,num,denom)
107int a,b,*num,*denom;
108{
109        int mn,i,j,casefix,match,blank;
110        float fnum = 0.0;
111        struct data_format *da,*db;
112        char ac,bc;
113
114        casefix = (usecase)? 0:32;
115        *num = 0;
116        *denom = 0;
117
118        da = &data[a];
119        db = &data[b];
120        mn = Min(da->length,db->length);
121
122        for(j=0;j<mn;j+=jump)
123        {
124                match = TRUE;
125                blank = TRUE;
126                for(i=0;i<width;i++)
127                {
128                        ac = da->nuc[j+i] | casefix;
129                        bc = db->nuc[j+i] | casefix;
130                        if(ac == 't')
131                                ac = 'u';
132                        if(ac == 'T')
133                                ac = 'U';
134                        if(bc == 't')
135                                bc = 'u';
136                        if(bc == 'T')
137                                bc = 'U';
138
139                        if((ac=='-') || (ac|32)=='n' || (ac==' ') ||
140                            (bc== '-') || (bc|32)=='n' || (bc==' '));
141
142                        else
143                        {
144                                blank = FALSE;
145                                if(ac != bc)
146                                {
147                                        match = FALSE;
148                                        switch(ac)
149                                        {
150                                        case 'a':
151                                                if (bc == 'c') fnum+=acwt;
152                                                else if(bc == 'g') fnum+=agwt;
153                                                else if(bc == 'u') fnum+=auwt;
154                                                break;
155
156                                        case 'c':
157                                                if (bc == 'a') fnum+=acwt;
158                                                else if(bc == 'g') fnum+=gcwt;
159                                                else if(bc == 'u') fnum+=ucwt;
160                                                break;
161
162                                        case 'g':
163                                                if (bc == 'a') fnum+=agwt;
164                                                else if(bc == 'c') fnum+=gcwt;
165                                                else if(bc == 'u') fnum+=ugwt;
166                                                break;
167
168                                        case 'u':
169                                                if (bc == 'a') fnum+=auwt;
170                                                else if(bc == 'c') fnum+=ucwt;
171                                                else if(bc == 'g') fnum+=ugwt;
172                                                break;
173
174                                        case 't':
175                                                if (bc == 'a') fnum+=auwt;
176                                                else if(bc == 'c') fnum+=ucwt;
177                                                else if(bc == 'g') fnum+=ugwt;
178                                                break;
179
180                                        default:
181                                                break;
182                                        };
183                                }
184                        }
185
186                        if((blank == FALSE) && match)
187                        {
188                                (*num) ++;
189                                (*denom) ++;
190                        }
191                        else if(!blank)
192                                (*denom) ++;
193                }
194        }
195        if(special)
196                (*num) = *denom - (int)fnum;
197}
198
199void Report()
200{
201        int i,jj,j;
202
203        if(tbl)
204                printf("#\n#-\n#-\n#-\n#-\n");
205        for(jj=0,j=0;j<numseq;j++)
206        {
207                if(tbl)
208                        printf("%2d: %-.15s|",jj+1,data[j].name);
209
210                for (i=0;i<j;i++)
211                {
212                        if(sim)
213                                printf("%6.1f",100 - dist[i][j]*100.0);
214                        else
215                                printf("%6.1f",dist[i][j]*100.0);
216                }
217                printf("\n");
218                jj++;
219        }
220}
221
222
223
224int main(ac,av)
225int ac;
226char **av;
227{
228        int j,k;
229        extern int ReadFlat();
230        FILE *file;
231
232        width = 1;
233        jump = 1;
234        if(ac==1)
235        {
236                fprintf(stderr,
237                    "usage: %s [-sim] [-case] [-c=<none,olsen,jukes>] ",av[0]);
238                fprintf(stderr,"[-t] alignment_flat_file\n");
239                return 1;
240        }
241        for(j=1;j<ac-1;j++)
242        {
243                getarg(av,j,argtyp,argval);
244                if(strcmp(argtyp,"-s=") == 0)
245                {
246                        j++;
247                        sscanf(argval,"%d",&start);
248                        start --;
249                }
250                else if(strcmp(argtyp,"-m=") == 0)
251                {
252                        j++;
253                        sscanf(argval,"%d",&width);
254                }
255                else if(strcmp(argtyp,"-j=") == 0)
256                {
257                        j++;
258                        sscanf(argval,"%d",&jump);
259                }
260                else if(strcmp(argtyp,"-case") == 0)
261                        usecase = TRUE;
262
263                else if(strcmp(argtyp,"-sim") == 0)
264                        sim = TRUE;
265
266                else if(strcmp(argtyp,"-c=") == 0)
267                {
268                        if(strcmp(argval,"olsen") == 0)
269                                correction = OLSEN;
270
271                        else if(strcmp(argval,"none") == 0)
272                                correction = NONE;
273
274                        else if(strcmp(argval,"jukes") == 0)
275                                correction = JUKES;
276
277                        else
278                                fprintf(stderr,"Correction type %s %s\n",
279                                    argval,"unknown, using JUKES");
280                }
281                else if(strcmp("-t",argtyp) == 0)
282                        tbl = TRUE;
283
284                else if(strcmp("-ac=",argtyp) == 0 || strcmp("-ca=",argtyp)== 0)
285                {
286                        j++;
287                        sscanf(argval,"%f",&acwt);
288                        special = TRUE;
289                }
290                else if(strcmp("-au=",argtyp) == 0 || strcmp("-ua=",argtyp)== 0)
291                {
292                        j++;
293                        sscanf(argval,"%f",&auwt);
294                        special = TRUE;
295                }
296                else if(strcmp("-ag=",argtyp) == 0 || strcmp("-ga=",argtyp)== 0)
297                {
298                        j++;
299                        sscanf(argval,"%f",&agwt);
300                        special = TRUE;
301                }
302                else if(strcmp("-uc=",argtyp) == 0 || strcmp("-cu=",argtyp)== 0)
303                {
304                        j++;
305                        sscanf(argval,"%f",&ucwt);
306                        special = TRUE;
307                }
308                else if(strcmp("-ug=",argtyp) == 0 || strcmp("-gu=",argtyp)== 0)
309                {
310                        j++;
311                        sscanf(argval,"%f",&ugwt);
312                        special = TRUE;
313                }
314                else if(strcmp("-gc=",argtyp) == 0 || strcmp("-cg=",argtyp)== 0)
315                {
316                        j++;
317                        sscanf(argval,"%f",&gcwt);
318                        special = TRUE;
319                }
320                else if(strcmp("-transition=",argtyp) == 0)
321                {
322                        j++;
323                        sscanf(argval,"%f",&ucwt);
324                        agwt = ucwt;
325                        special = TRUE;
326                }
327                else if(strcmp("-transversion=",argtyp) == 0)
328                {
329                        j++;
330                        sscanf(argval,"%f",&gcwt);
331                        ugwt = gcwt;
332                        acwt = gcwt;
333                        auwt = gcwt;
334                        special = TRUE;
335                }
336        }
337
338
339        file = fopen(av[ac-1],"r");
340        if ((file == NULL) || (ac == 1))
341        {
342                fprintf(stderr,"Error opening input file %s\n",av[ac-1]);
343                return 1;
344        }
345
346        numseq = ReadFlat(file,data,10000);
347
348        fclose(file);
349        SetPart();
350
351        for(j=0;j<numseq-1;j++)
352            for(k=j+1;k<numseq;k++)
353            {
354                int num, denom;
355                Compare(j,k,&num,&denom);
356                dist[j][k] = setdist(num,denom,j,k);
357            }
358
359        Report();
360        return 0;
361}
362
363
364
365float setdist(num,denom,a,b)
366int num,denom,a,b;
367{
368        float cor;
369        switch (correction)
370        {
371        case OLSEN:
372                cor = parta[a]*parta[b]+
373                    partc[a]*partc[b]+
374                    partg[a]*partg[b]+
375                    partu[a]*partu[b];
376                break;
377
378        case JUKES:
379                cor = 0.25;
380                break;
381
382        case NONE:
383                cor = 0.0;
384                break;
385
386        default:
387                cor = 0.0;
388                break;
389        };
390
391        if(correction == NONE)
392                return(1.0 - (float)num/(float)denom);
393        else
394                return( -(1.0-cor)*log(1.0 / (1.0-cor)*((float)num/(float)denom-cor)));
395}
396
397
398
399
400
401int find(b,a)
402char *a,*b;
403{
404        int flag,lenb,lena;
405        int i,j;
406
407        flag=0;
408        lenb=strlen(b);
409        lena=strlen(a);
410        for (i=0;((i<lena) && flag==0);i++)
411        {
412                for(j=0;(j<lenb) && (a[i+j]==b[j]);j++);
413                flag=((j==lenb)? 1:0);
414        }
415        return flag;
416}
417
Note: See TracBrowser for help on using the repository browser.