source: tags/initial/GDE/MOLPHY/getseq.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: 10.5 KB
Line 
1/*
2 * getseq.c   Adachi, J.   1994.04.28
3 * Copyright (C) 1992, 1993 J. Adachi & M. Hasegawa, All rights reserved.
4 */
5
6#include "protml.h"
7#include "prot_tml.h"
8
9void
10getsize(ifp, maxspc, numsite, commentp)
11FILE *ifp;
12int *maxspc;
13int *numsite;
14char **commentp; 
15{
16        char *cp, *np;
17        char line[BUFLINE];
18
19        if (fgets(line, BUFLINE, ifp) != NULL) {
20                if (sscanf(line, "%d %d", maxspc, numsite) == 2) {
21                        for (cp = line; isdigit(*cp) || isspace(*cp) /*&& *cp != '\0'*/; cp++)
22                                ;
23                        *commentp = new_cvector(strlen(cp) + 1);
24                        if (*cp != '\0') {
25                                for (np = *commentp; *cp != '\n' && *cp != '\0'; cp++)
26                                        *np++ = *cp;
27                                *np = '\0';
28                        } else {
29                                **commentp = '\0';
30                        }
31                        if (Debug)
32                                fprintf(stdout, "%d OTUs, %d sites,  %s\n\n",
33                                        *maxspc, *numsite, *commentp);
34                } else {
35                        fputs(line, stderr);
36                        fprintf(stderr, "\nBad format, first line of input file.\n");
37                        exit(1);
38                }
39        } else {
40                fprintf(stderr, "\nCan't read input file.\n");
41                exit(1);
42        }
43        return;
44} /*_ getsize */
45
46
47/* getid */
48void
49getid(ifp, identif, sciname, engname, nl, notu)
50FILE *ifp;
51char **identif;
52char **sciname;
53char **engname;
54int *nl;
55int *notu;
56{
57        int i;
58        char line[BUFLINE];
59        char idbuf[BUFLINE];
60        char scibuf[BUFLINE];
61        char engbuf[BUFLINE];
62        char *cp, *np;
63
64        while (fgets(line, BUFLINE, ifp) != NULL) {
65                (*nl)++; /* line counter */
66                if (line[0] == '#') /* comment line skip */
67                        continue;
68                for (cp = line; isspace(*cp); cp++) /* white line skip */
69                        ;
70                if (*cp == '\0')
71                        continue;
72                for (np = idbuf; !isspace(*cp); *np++ = *cp++) /* identifier */
73                        ;
74                *np = '\0';
75                for (i = 0; i < *notu; i++) { /* ID check */
76                        if (strcmp(idbuf, identif[i]) == 0) {
77                                fprintf(stderr,"Identifier \"%s\" is double defined\n",idbuf);
78                                exit(1);
79                        }
80                }
81                identif[*notu] =
82                        (char *)malloc((unsigned)(strlen(idbuf) + 1) * sizeof(char));
83                if (identif[*notu] == NULL) maerror("in getid, identif");
84                strcpy(identif[*notu], idbuf);
85
86                for ( ; isspace(*cp); cp++) /* white char skip */
87                        ;
88                /* science name */
89                if (*cp != '(' && *cp != '\0' && *cp != '#' && *cp != '\n') {
90                        for (np = scibuf; *cp!='(' && *cp!='#' && *cp!='\0' && *cp!='\n'; *np++ = *cp++)
91                                ;
92                        for ( ; isspace(*(np-1)); np--)
93                                ;
94                        *np = '\0';
95                } else {
96                        scibuf[0] = '\0';
97                }
98                sciname[*notu] =
99                        (char *)malloc((unsigned)(strlen(scibuf) + 1) * sizeof(char));
100                if (identif[*notu] == NULL) maerror("in getid, sciname");
101                strcpy(sciname[*notu], scibuf);
102
103                /* english name */
104                if (*cp != '\0' && *cp != '#' && *cp != '\n') {
105                        for (np = engbuf; *cp!='#' && *cp!='\0' && *cp!='\n'; *np++ = *cp++)
106                                ;
107                        for ( ; isspace(*(np-1)); np--)
108                                ;
109                        *np = '\0';
110                } else {
111                        engbuf[0] = '\0';
112                }
113                engname[*notu] =
114                        (char *)malloc((unsigned)(strlen(engbuf) + 1) * sizeof(char));
115                if (identif[*notu] == NULL) maerror("in getid, engname");
116                strcpy(engname[*notu], engbuf);
117
118                if (Debug) {
119                        printf("%3d OTU \"%s\"", *notu+1, identif[*notu]);
120                        if (sciname[*notu] != '\0') {
121                                fputs(" sci:", stdout);
122                                fputs(sciname[*notu], stdout);
123                        }
124                        if (engname[*notu] != '\0') {
125                                putchar(' ');
126                                fputs(" eng:", stdout);
127                                fputs(engname[*notu], stdout);
128                        }
129                        putchar('\n');
130                }
131                break;
132        }
133        return;
134} /*_ getid */
135
136
137/* getsites */
138void
139getsites(ifp, seqchar, numsite, nl, notu)
140FILE *ifp;
141cmatrix seqchar;
142int numsite, *nl, *notu;
143{
144        char line[BUFLINE];
145        char c;
146        int i, imax, nsite;
147
148        nsite = 0;
149        while (fgets(line, BUFLINE, ifp) != NULL) {
150                (*nl)++; /* line counter */
151                if (line[0] == '#') /* comment line skip */
152                        continue;
153                for (i = 0, imax = strlen(line); i < imax && nsite < numsite; i++) {
154                        c = toupper(line[i]);
155                        if (isacid(c)) {
156                                seqchar[*notu][nsite++] = acid2int(c);
157                        } else if (isspace(c)) {
158                                /* skip white char */
159                        } else {
160                                fputs(line, stderr);
161                                fprintf(stderr, "Expected %s Acid Character: '%c'\n",
162                                        Infomol, c);
163                                fprintf(stderr, "%d line %d column, %d OTU %d site\n",
164                                        *nl+1, i+1, *notu+1, nsite+1);
165                                exit(1);
166                        }
167                }
168                if (nsite >= numsite) {
169                        if (i != imax && !isspace(c = toupper(line[i])) && c != '\0') {
170                                fputs(line, stderr);
171                                fprintf(stderr, "Inconsistent size of sequence: '%c'\n", c);
172                                fprintf(stderr, "%d line %d column, %d OTU %d site\n",
173                                        *nl+1, i+1, *notu+1, nsite+1);
174                                exit(1);
175                        }
176                        break;
177                }
178        }
179        return;
180} /* getsites */
181
182
183void getseq(ifp, identif, sciname, engname, seqchar, maxspc, numsite)
184FILE *ifp;
185char **identif, **sciname, **engname;
186cmatrix seqchar;
187int maxspc, numsite;
188{
189        int notu;
190        int nl = 0;
191
192        for (notu = 0; notu < maxspc; notu++) {
193                getid(ifp, identif, sciname, engname, &nl, &notu);
194                getsites(ifp, seqchar, numsite, &nl, &notu);
195        }
196        return;
197} /*_ getseq */
198
199
200/* getidsites */
201void
202getidsites(ifp, identif, seqchar, numsite, nl, notu)
203FILE *ifp;
204cmatrix identif, seqchar;
205int numsite, *nl, *notu;
206{
207        char line[BUFLINE];
208        char idbuf[BUFLINE];
209        char c;
210        int i, imax, nsite;
211
212        nsite = 0;
213        while (fgets(line, BUFLINE, ifp) != NULL) {
214                (*nl)++; /* line counter */
215                if (line[0] == '#') /* comment line skip */
216                        continue;
217                if (nsite == 0) {
218                        strncpy(idbuf, line, NMLNGTH);
219                        i = NMLNGTH;
220                        while (isspace(idbuf[--i]))
221                                ;
222                        idbuf[++i] = '\0';
223                        for (i = 0; i < *notu; i++) { /* ID check */
224                                if (strcmp(idbuf, identif[i]) == 0) {
225                                        fprintf(stderr,
226                                                "Identifier \"%s\" is double defined\n",idbuf);
227                                        exit(1);
228                                }
229                        }
230                        identif[*notu] =
231                                (char *)malloc((unsigned)(strlen(idbuf) + 1) * sizeof(char));
232                        if (identif[*notu] == NULL) maerror("in getidsites, identif");
233                        strcpy(identif[*notu], idbuf);
234                        for (i=NMLNGTH, imax=strlen(line); i<imax && nsite<numsite; i++) {
235                                c = toupper(line[i]);
236                                if (isacid(c)) {
237                                        seqchar[*notu][nsite++] = acid2int(c);
238                                } else if (isspace(c)) {
239                                        /* skip white char */
240                                } else {
241                                        fputs(line, stderr);
242                                        fprintf(stderr, "Expected %s Acid Character: '%c'\n",
243                                                Infomol, c);
244                                        fprintf(stderr, "%d line %d column, %d OTU %d site\n",
245                                                *nl+1, i+1, *notu+1, nsite+1);
246                                        exit(1);
247                                }
248                        }
249                } else {
250                        for (i=0, imax=strlen(line); i < imax && nsite < numsite; i++) {
251                                c = toupper(line[i]);
252                                if (isacid(c)) {
253                                        seqchar[*notu][nsite++] = acid2int(c);
254                                } else if (isspace(c)) {
255                                        /* skip white char */
256                                } else {
257                                        fputs(line, stderr);
258                                        fprintf(stderr, "Expected %s Acid Character: '%c'\n",
259                                                Infomol, c);
260                                        fprintf(stderr, "%d line %d column, %d OTU %d site\n",
261                                                *nl+1, i+1, *notu+1, nsite+1);
262                                        exit(1);
263                                }
264                        }
265                }
266                if (nsite >= numsite) {
267                        if (i != imax && !isspace(c = toupper(line[i])) && c != '\0') {
268                                fputs(line, stderr);
269                                fprintf(stderr, "Inconsistent size of sequence: '%c'\n", c);
270                                fprintf(stderr, "%d line %d column, %d OTU %d site\n",
271                                        *nl+1, i+1, *notu+1, nsite+1);
272                                exit(1);
273                        }
274                        break;
275                }
276        }
277        return;
278}
279/*_ getidsites */
280
281
282void getseqs(ifp, identif, seqchar, maxspc, numsite)
283FILE *ifp;
284char **identif;
285cmatrix seqchar;
286int maxspc, numsite;
287{
288        int notu;
289        int nl = 0;
290
291        for (notu = 0; notu < maxspc; notu++)
292                getidsites(ifp, identif, seqchar, numsite, &nl, &notu);
293        return;
294} /*_ getseqs */
295
296
297void getseqi(ifp, identif, seqchar, maxspc, numsite)
298FILE *ifp;
299char **identif;
300cmatrix seqchar;
301int maxspc, numsite;
302{
303        int notu, nl, nsite, ns, i, imax;
304        char line[BUFLINE];
305        char idbuf[BUFLINE];
306        char c, *cp;
307
308        nsite = 0;
309        nl = 0;
310        while (nsite < numsite) {
311        for (notu = 0; notu < maxspc; notu++) {
312                while (fgets(line, BUFLINE, ifp) != NULL) {
313                        nl++; /* line counter */
314                        for (cp = line; isspace(*cp); cp++) /* white line skip */
315                                ;
316                        if (*cp == '\0')
317                                continue;
318                        if (nsite == 0) {
319                                strncpy(idbuf, line, NMLNGTH);
320                                i = NMLNGTH;
321                                while (isspace(idbuf[--i]))
322                                        ;
323                                idbuf[++i] = '\0';
324                                for (i = 0; i < notu; i++) { /* ID check */
325                                        if (strcmp(idbuf, identif[i]) == 0) {
326                                                fprintf(stderr,
327                                                        "Identifier \"%s\" is double defined\n",idbuf);
328                                                exit(1);
329                                        }
330                                }
331                                identif[notu] =
332                                        (char *)malloc((unsigned)(strlen(idbuf)+1) * sizeof(char));
333                                if (identif[notu] == NULL) maerror("in getidsites, identif");
334                                strcpy(identif[notu], idbuf);
335                                i = NMLNGTH;
336                        } else {
337                                i = 0;
338                        }
339                        for (ns=nsite, imax=strlen(line); i<imax && ns<numsite; i++) {
340                                c = toupper(line[i]);
341                                if (isacid(c)) {
342                                        seqchar[notu][ns++] = acid2int(c);
343                                } else if (isspace(c)) {
344                                        /* skip white char */
345                                } else {
346                                        fputs(line, stderr);
347                                        fprintf(stderr, "Expected %s Acid Character: '%c'\n",
348                                                Infomol, c);
349                                        fprintf(stderr,"%d line %d column, %d OTU %d site\n",
350                                                nl+1, i+1, notu+1, ns+1);
351                                        exit(1);
352                                }
353                        }
354                        if (notu == maxspc-1) nsite = ns;
355                        if (nsite >= numsite) {
356                                if (i != imax && !isspace(c = toupper(line[i])) && c != '\0') {
357                                        fputs(line, stderr);
358                                        fprintf(stderr, "Inconsistent size of sequence: '%c'\n", c);
359                                        fprintf(stderr, "%d line %d column, %d OTU %d site\n",
360                                                nl+1, i+1, notu+1, nsite+1);
361                                        exit(1);
362                                }
363                                break;
364                        }
365                        break;
366                }
367        }
368        }
369        return;
370
371} /*_ getseqi */
372
373
374/* put identifier of otu */
375void
376fputid(ofp, name, maxcolumn)
377FILE *ofp;
378char *name;
379int maxcolumn;
380{
381        int i, imax;
382
383        if ((imax = strlen(name)) < maxcolumn) {
384                fputs(name, ofp);
385                for (i = imax; i < maxcolumn; i++)
386                        fputc(' ', ofp);
387        } else {
388                for (i = 0; i < maxcolumn; i++)
389                        fputc(name[i], ofp);
390        }
391}
392
393
394/* print sequences */
395void
396prsequence(ofp, identif, seqchar, maxspc, maxsite)
397FILE *ofp;
398char **identif;
399char **seqchar;
400int maxspc, maxsite;
401{
402        int i, nst, maxst, nsp;
403        char ch;
404        ivector counti;
405        cvector consenseq;
406
407        counti = new_ivector(Tpmradix + 1);
408        consenseq = new_cvector(maxsite);
409        for (nst = 0; nst < maxsite; nst++) { /* make consensus sequence */
410                for (i = 0; i < Tpmradix + 1; i++) counti[i] = 0;
411                for (nsp = 0; nsp < maxspc; nsp++) { /* count char */
412                        counti[ seqchar[nsp][nst] ]++; 
413                }
414                consenseq[nst] = '.';
415                for (i = 0; i < Tpmradix; i++) { /* + 1 */
416                        if (counti[i] > maxspc / 2) {
417                                consenseq[nst] = int2acid(i);
418                                break;
419                        }
420                }
421        }
422
423        fputc('\n', ofp);
424        for (i = 0; i < maxsite; i += Linesites) {
425                maxst = (i + Linesites < maxsite) ? (i + Linesites) : maxsite;
426                fputid(ofp, "CONSENSUS", 10);
427                for (nst = i; nst < maxst; nst++) {
428                        if (nst % 10 == 0) fputc(' ', ofp);
429                        fputc(consenseq[nst], ofp);
430                }
431                fputc('\n', ofp);
432                for (nsp = 0; nsp < maxspc; nsp++) {
433                        fputid(ofp, identif[nsp], 10);
434                        for (nst = i; nst < maxst; nst++) {
435                                if (nst % 10 == 0) fputc(' ', ofp);
436                                ch = chint2acid(seqchar[nsp][nst]);
437                                if (ch != consenseq[nst])
438                                        fputc(ch, ofp);
439                                else
440                                        fputc('.', ofp);
441                        }
442                        fputc('\n', ofp);
443                }
444                fputid(ofp, "", 10); 
445                for (nst = i; nst < maxst; nst+=10) {
446                        if (nst+10 <= maxst) printf(" %10d", nst+10);
447                }
448                fputc('\n', ofp);
449        }
450        free_cvector(consenseq);
451        free_ivector(counti);
452        return;
453} /*_ prsequence */
Note: See TracBrowser for help on using the repository browser.