source: branches/profile/GDE/SUPPORT/count.c

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