source: tags/arb-6.0/GDE/CLUSTALW/sequence.c

Last change on this file was 1754, checked in by westram, 21 years ago

updated to version 1.83
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 33.4 KB
Line 
1/********* Sequence input routines for CLUSTAL W *******************/
2/* DES was here.  FEB. 1994 */
3/* Now reads PILEUP/MSF and CLUSTAL alignment files */
4
5#include <stdio.h>
6#include <string.h>
7#include <ctype.h>
8#include <stdlib.h>
9#include "clustalw.h"   
10
11#define MIN(a,b) ((a)<(b)?(a):(b))
12
13
14
15/*
16*       Prototypes
17*/
18
19static char * get_seq(char *,sint *,char *);
20static char * get_clustal_seq(char *,sint *,char *,sint);
21static char * get_msf_seq(char *,sint *,char *,sint);
22static void check_infile(sint *);
23static void p_encode(char *, char *, sint);
24static void n_encode(char *, char *, sint);
25static sint res_index(char *,char);
26static Boolean check_dnaflag(char *, sint);
27static sint count_clustal_seqs(void);
28static sint count_pir_seqs(void);
29static sint count_msf_seqs(void);
30static sint count_rsf_seqs(void);
31static void get_swiss_feature(char *line,sint len);
32static void get_rsf_feature(char *line,sint len);
33static void get_swiss_mask(char *line,sint len);
34static void get_clustal_ss(sint length);
35static void get_embl_ss(sint length);
36static void get_rsf_ss(sint length);
37static void get_gde_ss(sint length);
38static Boolean cl_blankline(char *line);
39
40/*
41 *      Global variables
42 */
43extern sint max_names;
44FILE *fin;
45extern Boolean usemenu, dnaflag, explicit_dnaflag;
46extern Boolean interactive;
47extern char seqname[];
48extern sint nseqs;
49extern sint *seqlen_array;
50extern sint *output_index;
51extern char **names,**titles;
52extern char **seq_array;
53extern Boolean profile1_empty, profile2_empty;
54extern sint gap_pos2;
55extern sint max_aln_length;
56extern char *gap_penalty_mask, *sec_struct_mask;
57extern sint struct_penalties;
58extern char *ss_name;
59extern sint profile_no;
60extern sint debug;
61
62char *amino_acid_codes   =    "ABCDEFGHIKLMNPQRSTUVWXYZ-";  /* DES */
63static sint seqFormat;
64static char chartab[128];
65static char *formatNames[] = {"unknown","EMBL/Swiss-Prot","PIR",
66                              "Pearson","GDE","Clustal","Pileup/MSF","RSF","USER","PHYLIP","NEXUS"};
67
68void fill_chartab(void) /* Create translation and check table */
69{
70        register sint i;
71        register char c;
72       
73        for(i=0;i<128;chartab[i++]=0);
74        for(i=0;(c=amino_acid_codes[i]);i++)
75                chartab[(int)c]=chartab[tolower(c)]=c;
76}
77
78static char * get_msf_seq(char *sname,sint *len,char *tit,sint seqno)
79/* read the seqno_th. sequence from a PILEUP multiple alignment file */
80{
81        static char line[MAXLINE+1];
82        char *seq = NULL;
83        sint i,j,k;
84        unsigned char c;
85
86        fseek(fin,0,0);                 /* start at the beginning */
87
88        *len=0;                         /* initialise length to zero */
89        for(i=0;;i++) {
90                if(fgets(line,MAXLINE+1,fin)==NULL) return NULL; /* read the title*/
91                if(linetype(line,"//") ) break;             /* lines...ignore*/
92        }
93
94        while (fgets(line,MAXLINE+1,fin) != NULL) {
95                if(!blankline(line)) {
96
97                        for(i=1;i<seqno;i++) fgets(line,MAXLINE+1,fin);
98                        for(j=0;j<=strlen(line);j++) if(line[j] != ' ') break;
99                        for(k=j;k<=strlen(line);k++) if(line[k] == ' ') break;
100                        strncpy(sname,line+j,MIN(MAXNAMES,k-j)); 
101                        sname[MIN(MAXNAMES,k-j)]=EOS;
102                        rtrim(sname);
103                        blank_to_(sname);
104
105                        if(seq==NULL)
106                                seq=(char *)ckalloc((MAXLINE+2)*sizeof(char));
107                        else
108                                seq=(char *)ckrealloc(seq,((*len)+MAXLINE+2)*sizeof(char));
109                        for(i=k;i<=MAXLINE;i++) {
110                                c=line[i];
111                                if(c == '.' || c == '~' ) c = '-';
112                                if(c == '*') c = 'X';
113                                if(c == '\n' || c == EOS) break; /* EOL */
114                                c=chartab[c];
115                                if(c) seq[++(*len)]=c;
116                        }
117
118                        for(i=0;;i++) {
119                                if(fgets(line,MAXLINE+1,fin)==NULL) return seq;
120                                if(blankline(line)) break;
121                        }
122                }
123        }
124        return seq;
125}
126
127static Boolean cl_blankline(char *line)
128{
129        int i;
130
131        if (line[0] == '!') return TRUE;
132       
133        for(i=0;line[i]!='\n' && line[i]!=EOS;i++) {
134                if( isdigit(line[i]) ||
135                    isspace(line[i]) ||
136                    (line[i] == '*') ||
137                    (line[i] == ':') ||
138                    (line[i] == '.')) 
139                        ;
140                else
141                        return FALSE;
142        }
143        return TRUE;
144}
145
146static char * get_clustal_seq(char *sname,sint *len,char *tit,sint seqno)
147/* read the seqno_th. sequence from a clustal multiple alignment file */
148{
149        static char line[MAXLINE+1];
150        static char tseq[MAXLINE+1];
151        char *seq = NULL;
152        sint i,j;
153        unsigned char c;
154
155        fseek(fin,0,0);                 /* start at the beginning */
156
157        *len=0;                         /* initialise length to zero */
158        fgets(line,MAXLINE+1,fin);      /* read the title line...ignore it */
159
160        while (fgets(line,MAXLINE+1,fin) != NULL) {
161                if(!cl_blankline(line)) {
162
163                        for(i=1;i<seqno;i++) fgets(line,MAXLINE+1,fin);
164                        for(j=0;j<=strlen(line);j++) if(line[j] != ' ') break;
165
166                        sscanf(line,"%s%s",sname,tseq);
167                        for(j=0;j<MAXNAMES;j++) if(sname[j] == ' ') break;
168                        sname[j]=EOS;
169                        rtrim(sname);
170                        blank_to_(sname);
171
172                        if(seq==NULL)
173                                seq=(char *)ckalloc((MAXLINE+2)*sizeof(char));
174                        else
175                                seq=(char *)ckrealloc(seq,((*len)+MAXLINE+2)*sizeof(char));
176                        for(i=0;i<=MAXLINE;i++) {
177                                c=tseq[i];
178                                /*if(c == '\n' || c == EOS) break;*/ /* EOL */
179                                if(isspace(c) || c == EOS) break; /* EOL */
180                                c=chartab[c];
181                                if(c) seq[++(*len)]=c;
182                        }
183
184                        for(i=0;;i++) {
185                                if(fgets(line,MAXLINE+1,fin)==NULL) return seq;
186                                if(cl_blankline(line)) break;
187                        }
188                }
189        }
190
191        return seq;
192}
193
194static void get_clustal_ss(sint length)
195/* read the structure data from a clustal multiple alignment file */
196{
197        static char title[MAXLINE+1];
198        static char line[MAXLINE+1];
199        static char lin2[MAXLINE+1];
200        static char tseq[MAXLINE+1];
201        static char sname[MAXNAMES+1];
202        sint i,j,len,ix,struct_index=0;
203        char c;
204
205       
206        fseek(fin,0,0);                 /* start at the beginning */
207
208        len=0;                          /* initialise length to zero */
209        if (fgets(line,MAXLINE+1,fin) == NULL) return;  /* read the title line...ignore it */
210
211        if (fgets(line,MAXLINE+1,fin) == NULL) return;  /* read the next line... */
212/* skip any blank lines */
213        for (;;) {
214                if(fgets(line,MAXLINE+1,fin)==NULL) return;
215                if(!blankline(line)) break;
216        }
217
218/* look for structure table lines */
219        ix = -1;
220        for(;;) {
221                if(line[0] != '!') break;
222                if(strncmp(line,"!SS",3) == 0) {
223                        ix++;
224                        sscanf(line+4,"%s%s",sname,tseq);
225                        for(j=0;j<MAXNAMES;j++) if(sname[j] == ' ') break;
226                        sname[j]=EOS;
227                        rtrim(sname);
228                blank_to_(sname);
229                if (interactive) {
230                                strcpy(title,"Found secondary structure in alignment file: ");
231                                strcat(title,sname);
232                                (*lin2)=prompt_for_yes_no(title,"Use it to set local gap penalties ");
233                        }
234                        else (*lin2) = 'y';
235                        if ((*lin2 != 'n') && (*lin2 != 'N'))  {               
236                                struct_penalties = SECST;
237                                struct_index = ix;
238                                for (i=0;i<length;i++)
239                                {
240                                        sec_struct_mask[i] = '.';
241                                        gap_penalty_mask[i] = '.';
242                                }
243                                strcpy(ss_name,sname);
244                                for(i=0;len < length;i++) {
245                                        c = tseq[i];
246                                        if(c == '\n' || c == EOS) break; /* EOL */
247                                        if (!isspace(c)) sec_struct_mask[len++] = c;
248                                }
249                        }
250                }
251                else if(strncmp(line,"!GM",3) == 0) {
252                        ix++;
253                        sscanf(line+4,"%s%s",sname,tseq);
254                        for(j=0;j<MAXNAMES;j++) if(sname[j] == ' ') break;
255                        sname[j]=EOS;
256                        rtrim(sname);
257                blank_to_(sname);
258                if (interactive) {
259                                strcpy(title,"Found gap penalty mask in alignment file: ");
260                                strcat(title,sname);
261                                (*lin2)=prompt_for_yes_no(title,"Use it to set local gap penalties ");
262                        }
263                        else (*lin2) = 'y';
264                        if ((*lin2 != 'n') && (*lin2 != 'N'))  {               
265                                struct_penalties = GMASK;
266                                struct_index = ix;
267                                for (i=0;i<length;i++)
268                                        gap_penalty_mask[i] = '1';
269                                        strcpy(ss_name,sname);
270                                for(i=0;len < length;i++) {
271                                        c = tseq[i];
272                                        if(c == '\n' || c == EOS) break; /* EOL */
273                                        if (!isspace(c)) gap_penalty_mask[len++] = c;
274                                }
275                        }
276                }
277                if (struct_penalties != NONE) break;
278                if(fgets(line,MAXLINE+1,fin)==NULL) return;
279        }
280                       
281        if (struct_penalties == NONE) return;
282       
283/* skip any more comment lines */
284        while (line[0] == '!') {
285                if(fgets(line,MAXLINE+1,fin)==NULL) return;
286        }
287
288/* skip the sequence lines and any comments after the alignment */
289        for (;;) {
290                if(isspace(line[0])) break;
291                if(fgets(line,MAXLINE+1,fin)==NULL) return;
292        }
293                       
294
295/* read the rest of the alignment */
296       
297        for (;;) {
298/* skip any blank lines */
299                        for (;;) {
300                                if(!blankline(line)) break;
301                                if(fgets(line,MAXLINE+1,fin)==NULL) return;
302                        }
303/* get structure table line */
304                        for(ix=0;ix<struct_index;ix++) {
305                                if (line[0] != '!') {
306                                        if(struct_penalties == SECST)
307                                                error("bad secondary structure format");
308                                        else
309                                                error("bad gap penalty mask format");
310                                        struct_penalties = NONE;
311                                        return;
312                                }
313                                if(fgets(line,MAXLINE+1,fin)==NULL) return;
314                        }
315                        if(struct_penalties == SECST) {
316                                if (strncmp(line,"!SS",3) != 0) {
317                                        error("bad secondary structure format");
318                                        struct_penalties = NONE;
319                                        return;
320                                }
321                                sscanf(line+4,"%s%s",sname,tseq);
322                                for(i=0;len < length;i++) {
323                                        c = tseq[i];
324                                        if(c == '\n' || c == EOS) break; /* EOL */
325                                        if (!isspace(c)) sec_struct_mask[len++] = c;
326                                }                       
327                        }
328                        else if (struct_penalties == GMASK) {
329                                if (strncmp(line,"!GM",3) != 0) {
330                                        error("bad gap penalty mask format");
331                                        struct_penalties = NONE;
332                                        return;
333                                }
334                                sscanf(line+4,"%s%s",sname,tseq);
335                                for(i=0;len < length;i++) {
336                                        c = tseq[i];
337                                        if(c == '\n' || c == EOS) break; /* EOL */
338                                        if (!isspace(c)) gap_penalty_mask[len++] = c;
339                                }                       
340                        }
341
342/* skip any more comment lines */
343                while (line[0] == '!') {
344                        if(fgets(line,MAXLINE+1,fin)==NULL) return;
345                }
346
347/* skip the sequence lines */
348                for (;;) {
349                        if(isspace(line[0])) break;
350                        if(fgets(line,MAXLINE+1,fin)==NULL) return;
351                }
352        }
353}
354
355static void get_embl_ss(sint length)
356{
357        static char title[MAXLINE+1];
358        static char line[MAXLINE+1];
359        static char lin2[MAXLINE+1];
360        static char sname[MAXNAMES+1];
361        char feature[MAXLINE+1];
362        sint i;
363
364/* find the start of the sequence entry */
365        for (;;) {
366                while( !linetype(line,"ID") )
367                        if (fgets(line,MAXLINE+1,fin) == NULL) return;
368                       
369        for(i=5;i<=strlen(line);i++)  /* DES */
370                        if(line[i] != ' ') break;
371                strncpy(sname,line+i,MAXNAMES); /* remember entryname */
372                for(i=0;i<=strlen(sname);i++)
373                        if(sname[i] == ' ') {
374                                sname[i]=EOS;
375                                break;
376                        }
377                sname[MAXNAMES]=EOS;
378                rtrim(sname);
379        blank_to_(sname);
380               
381/* look for secondary structure feature table / gap penalty mask */
382                while(fgets(line,MAXLINE+1,fin) != NULL) {
383                        if (linetype(line,"FT")) {
384                                sscanf(line+2,"%s",feature);
385                                if (strcmp(feature,"HELIX") == 0 ||
386                                    strcmp(feature,"STRAND") == 0)
387                                {
388
389                                if (interactive) {
390                                        strcpy(title,"Found secondary structure in alignment file: ");
391                                        strcat(title,sname);
392                                        (*lin2)=prompt_for_yes_no(title,"Use it to set local gap penalties ");
393                                }
394                                else (*lin2) = 'y';
395                                if ((*lin2 != 'n') && (*lin2 != 'N'))  {               
396                                        struct_penalties = SECST;
397                                        for (i=0;i<length;i++)
398                                                sec_struct_mask[i] = '.';
399                                        do {
400                                                get_swiss_feature(&line[2],length);
401                                                fgets(line,MAXLINE+1,fin);
402                                        } while( linetype(line,"FT") );
403                                }
404                                else {
405                                        do {
406                                                fgets(line,MAXLINE+1,fin);
407                                        } while( linetype(line,"FT") );
408                                }
409                                strcpy(ss_name,sname);
410                                }
411                        }
412                        else if (linetype(line,"GM")) {
413                                if (interactive) {
414                                        strcpy(title,"Found gap penalty mask in alignment file: ");
415                                        strcat(title,sname);
416                                        (*lin2)=prompt_for_yes_no(title,"Use it to set local gap penalties ");
417                                }
418                                else (*lin2) = 'y';
419                                if ((*lin2 != 'n') && (*lin2 != 'N'))  {               
420                                        struct_penalties = GMASK;
421                                        for (i=0;i<length;i++)
422                                                gap_penalty_mask[i] = '1';
423                                        do {
424                                                get_swiss_mask(&line[2],length);
425                                                fgets(line,MAXLINE+1,fin);
426                                        } while( linetype(line,"GM") );
427                                }
428                                else {
429                                        do {
430                                                fgets(line,MAXLINE+1,fin);
431                                        } while( linetype(line,"GM") );
432                                }
433                                strcpy(ss_name,sname);
434                        }
435                        if (linetype(line,"SQ"))
436                                break; 
437
438                        if (struct_penalties != NONE) break;                   
439                }
440                                               
441        }
442                                               
443}
444
445static void get_rsf_ss(sint length)
446{
447        static char title[MAXLINE+1];
448        static char line[MAXLINE+1];
449        static char lin2[MAXLINE+1];
450        static char sname[MAXNAMES+1];
451        sint i;
452
453/* skip the comments */
454        while (fgets(line,MAXLINE+1,fin) != NULL) {
455                if(line[strlen(line)-2]=='.' &&
456                                 line[strlen(line)-3]=='.')
457                        break;
458        }
459
460/* find the start of the sequence entry */
461        for (;;) {
462                while (fgets(line,MAXLINE+1,fin) != NULL)
463                        if( *line == '{' ) break;
464
465                while( !keyword(line,"name") )
466                        if (fgets(line,MAXLINE+1,fin) == NULL) return;
467                       
468        for(i=5;i<=strlen(line);i++)  /* DES */
469                        if(line[i] != ' ') break;
470                strncpy(sname,line+i,MAXNAMES); /* remember entryname */
471                for(i=0;i<=strlen(sname);i++)
472                        if(sname[i] == ' ') {
473                                sname[i]=EOS;
474                                break;
475                        }
476                sname[MAXNAMES]=EOS;
477                rtrim(sname);
478        blank_to_(sname);
479               
480/* look for secondary structure feature table / gap penalty mask */
481                while(fgets(line,MAXLINE+1,fin) != NULL) {
482                        if (keyword(line,"feature")) {
483                                if (interactive) {
484                                        strcpy(title,"Found secondary structure in alignment file: ");
485                                        strcat(title,sname);
486                                        (*lin2)=prompt_for_yes_no(title,"Use it to set local gap penalties ");
487                                }
488                                else (*lin2) = 'y';
489                                if ((*lin2 != 'n') && (*lin2 != 'N'))  {               
490                                        struct_penalties = SECST;
491                                        for (i=0;i<length;i++)
492                                                sec_struct_mask[i] = '.';
493                                        do {
494                                                if(keyword(line,"feature"))
495                                                        get_rsf_feature(&line[7],length);
496                                                fgets(line,MAXLINE+1,fin);
497                                        } while( !keyword(line,"sequence") );
498                                }
499                                else {
500                                        do {
501                                                fgets(line,MAXLINE+1,fin);
502                                        } while( !keyword(line,"sequence") );
503                                }
504                                strcpy(ss_name,sname);
505                        }
506                        else if (keyword(line,"sequence"))
507                                break; 
508
509                        if (struct_penalties != NONE) break;                   
510                }
511                                               
512        }
513                                               
514}
515
516static void get_gde_ss(sint length)
517{
518        static char title[MAXLINE+1];
519        static char line[MAXLINE+1];
520        static char lin2[MAXLINE+1];
521        static char sname[MAXNAMES+1];
522        sint i, len, offset = 0;
523        unsigned char c;
524
525        for (;;) {
526                line[0] = '\0';
527/* search for the next comment line */
528                while(*line != '"')
529                        if (fgets(line,MAXLINE+1,fin) == NULL) return;
530
531/* is it a secondary structure entry? */
532                if (strncmp(&line[1],"SS_",3) == 0) {
533                        for (i=1;i<=MAXNAMES-3;i++) {
534                                if (line[i+3] == '(' || line[i+3] == '\n')
535                                                break;
536                                sname[i-1] = line[i+3];
537                        }
538                        i--;
539                        sname[i]=EOS;
540                        if (sname[i-1] == '(') sscanf(&line[i+3],"%d",&offset);
541                        else offset = 0;
542                        for(i--;i > 0;i--) 
543                                if(isspace(sname[i])) {
544                                        sname[i]=EOS;   
545                                }
546                                else break;             
547                        blank_to_(sname);
548
549                        if (interactive) {
550                                strcpy(title,"Found secondary structure in alignment file: ");
551                                strcat(title,sname);
552                                (*lin2)=prompt_for_yes_no(title,"Use it to set local gap penalties ");
553                        }
554                        else (*lin2) = 'y';
555                        if ((*lin2 != 'n') && (*lin2 != 'N'))  {               
556                                struct_penalties = SECST;
557                                for (i=0;i<length;i++)
558                                        sec_struct_mask[i] = '.';
559                                len = 0;
560                                while(fgets(line,MAXLINE+1,fin)) {
561                                        if(*line == '%' || *line == '#' || *line == '"') break;
562                                        for(i=offset;i < length;i++) {
563                                                c=line[i];
564                                                if(c == '\n' || c == EOS) 
565                                                        break;                  /* EOL */
566                                                sec_struct_mask[len++]=c;
567                                        }
568                                        if (len > length) break;
569                                }
570                                strcpy(ss_name,sname);
571                        }
572                }
573/* or is it a gap penalty mask entry? */
574                else if (strncmp(&line[1],"GM_",3) == 0) {
575                        for (i=1;i<=MAXNAMES-3;i++) {
576                                if (line[i+3] == '(' || line[i+3] == '\n')
577                                                break;
578                                sname[i-1] = line[i+3];
579                        }
580                        i--;
581                        sname[i]=EOS;
582                        if (sname[i-1] == '(') sscanf(&line[i+3],"%d",&offset);
583                        else offset = 0;
584                        for(i--;i > 0;i--) 
585                                if(isspace(sname[i])) {
586                                        sname[i]=EOS;   
587                                }
588                                else break;             
589                        blank_to_(sname);
590
591                        if (interactive) {
592                                strcpy(title,"Found gap penalty mask in alignment file: ");
593                                strcat(title,sname);
594                                (*lin2)=prompt_for_yes_no(title,"Use it to set local gap penalties ");
595                        }
596                        else (*lin2) = 'y';
597                        if ((*lin2 != 'n') && (*lin2 != 'N'))  {               
598                                struct_penalties = GMASK;
599                                for (i=0;i<length;i++)
600                                        gap_penalty_mask[i] = '1';
601                                len = 0;
602                                while(fgets(line,MAXLINE+1,fin)) {
603                                        if(*line == '%' || *line == '#' || *line == '"') break;
604                                        for(i=offset;i < length;i++) {
605                                                c=line[i];
606                                                if(c == '\n' || c == EOS) 
607                                                        break;                  /* EOL */
608                                                gap_penalty_mask[len++]=c;
609                                        }
610                                        if (len > length) break;
611                                }
612                                strcpy(ss_name,sname);
613                        }
614                }
615                if (struct_penalties != NONE) break;                   
616        }                       
617                       
618}
619
620static void get_swiss_feature(char *line, sint len)
621{
622        char c, s, feature[MAXLINE+1];
623        int  i, start_pos, end_pos;
624       
625        if (sscanf(line,"%s%d%d",feature,&start_pos,&end_pos) != 3) {
626                return;
627        }
628
629        if (strcmp(feature,"HELIX") == 0) {
630                c = 'A';
631                s = '$';
632        }
633        else if (strcmp(feature,"STRAND") == 0) {
634                c = 'B';
635                s = '%';
636        }
637        else
638                return;
639                       
640        if(start_pos >=len || end_pos>=len) return;
641
642        sec_struct_mask[start_pos-1] = s;
643        for (i=start_pos;i<end_pos-1;i++)
644                sec_struct_mask[i] = c;
645        sec_struct_mask[end_pos-1] = s;
646               
647}
648
649static void get_rsf_feature(char *line, sint len)
650{
651        char c, s;
652        char str1[MAXLINE+1],str2[MAXLINE+1],feature[MAXLINE+1];
653        int  i, tmp,start_pos, end_pos;
654       
655        if (sscanf(line,"%d%d%d%s%s%s",&start_pos,&end_pos,&tmp,str1,str2,feature) != 6) {
656                return;
657        }
658
659        if (strcmp(feature,"HELIX") == 0) {
660                c = 'A';
661                s = '$';
662        }
663        else if (strcmp(feature,"STRAND") == 0) {
664                c = 'B';
665                s = '%';
666        }
667        else
668                return;
669                       
670        if(start_pos>=len || end_pos >= len) return;
671        sec_struct_mask[start_pos-1] = s;
672        for (i=start_pos;i<end_pos-1;i++)
673                sec_struct_mask[i] = c;
674        sec_struct_mask[end_pos-1] = s;
675               
676}
677
678static void get_swiss_mask(char *line, sint len)
679{
680        int  i, value, start_pos, end_pos;
681       
682        if (sscanf(line,"%d%d%d",&value,&start_pos,&end_pos) != 3) {
683                return;
684        }
685
686        if (value < 1 || value > 9) return;
687       
688        if(start_pos>=len || end_pos >= len) return;
689        for (i=start_pos-1;i<end_pos;i++)
690                gap_penalty_mask[i] = value+'0';
691               
692}
693
694static char * get_seq(char *sname,sint *len,char *tit)
695{
696        static char line[MAXLINE+1];
697        char *seq = NULL;
698        sint i, offset = 0;
699        unsigned char c=EOS;
700        Boolean got_seq=FALSE;
701
702        switch(seqFormat) {
703
704/************************************/
705                case EMBLSWISS:
706                        while( !linetype(line,"ID") )
707                                if (fgets(line,MAXLINE+1,fin) == NULL) return NULL;
708                       
709                        for(i=5;i<=strlen(line);i++)  /* DES */
710                                if(line[i] != ' ') break;
711                        strncpy(sname,line+i,MAXNAMES); /* remember entryname */
712                        for(i=0;i<=strlen(sname);i++)
713                                if(sname[i] == ' ') {
714                                        sname[i]=EOS;
715                                        break;
716                                }
717
718                        sname[MAXNAMES]=EOS;
719                        rtrim(sname);
720                        blank_to_(sname);
721
722                                               
723                        while( !linetype(line,"SQ") )
724                                fgets(line,MAXLINE+1,fin);
725                               
726                        *len=0;
727                        while(fgets(line,MAXLINE+1,fin)) {
728                                if(got_seq && blankline(line)) break;
729                                if( strlen(line) > 2 && line[strlen(line)-2]=='.' && line[strlen(line)-3]=='.' ) 
730                                        continue;
731                                if(seq==NULL)
732                                        seq=(char *)ckalloc((MAXLINE+2)*sizeof(char));
733                                else
734                                        seq=(char *)ckrealloc(seq,((*len)+MAXLINE+2)*sizeof(char));
735                                for(i=0;i<=MAXLINE;i++) {
736                                        c=line[i];
737                                if(c == '\n' || c == EOS || c == '/')
738                                        break;                  /* EOL */
739                                c=chartab[c];
740                                if(c) {
741                                        got_seq=TRUE;
742                                        seq[++(*len)]=c;
743                                }
744                                }
745                        if(c == '/') break;
746                        }
747                break;
748               
749/************************************/
750                case PIR:
751                        while(*line != '>')
752                                fgets(line,MAXLINE+1,fin);                     
753                        for(i=4;i<=strlen(line);i++)  /* DES */
754                                if(line[i] != ' ') break;
755                        strncpy(sname,line+i,MAXNAMES); /* remember entryname */
756                        sname[MAXNAMES]=EOS;
757                        rtrim(sname);
758                        blank_to_(sname);
759
760                        fgets(line,MAXLINE+1,fin);
761                        strncpy(tit,line,MAXTITLES);
762                        tit[MAXTITLES]=EOS;
763                        i=strlen(tit);
764                        if(tit[i-1]=='\n') tit[i-1]=EOS;
765                       
766                        *len=0;
767                        while(fgets(line,MAXLINE+1,fin)) {
768                                if(seq==NULL)
769                                        seq=(char *)ckalloc((MAXLINE+2)*sizeof(char));
770                                else
771                                        seq=(char *)ckrealloc(seq,((*len)+MAXLINE+2)*sizeof(char));
772                                for(i=0;i<=MAXLINE;i++) {
773                                        c=line[i];
774                                if(c == '\n' || c == EOS || c == '*')
775                                        break;                  /* EOL */
776                       
777                                c=chartab[c];
778                                if(c) seq[++(*len)]=c;
779                                }
780                        if(c == '*') break;
781                        }
782                break;
783/***********************************************/
784                case PEARSON:
785                        while(*line != '>')
786                                fgets(line,MAXLINE+1,fin);
787                       
788                        for(i=1;i<=strlen(line);i++)  /* DES */
789                                if(line[i] != ' ') break;
790                        strncpy(sname,line+i,MAXNAMES); /* remember entryname */
791                        for(i=1;i<=strlen(sname);i++)  /* DES */
792                                if(sname[i] == ' ') break;
793                        sname[i]=EOS;
794                        rtrim(sname);
795                        blank_to_(sname);
796
797                        *tit=EOS;
798                       
799                        *len=0;
800                        while(fgets(line,MAXLINE+1,fin)) {
801                                if(seq==NULL)
802                                        seq=(char *)ckalloc((MAXLINE+2)*sizeof(char));
803                                else
804                                        seq=(char *)ckrealloc(seq,((*len)+MAXLINE+2)*sizeof(char));
805                                for(i=0;i<=MAXLINE;i++) {
806                                        c=line[i];
807                                if(c == '\n' || c == EOS || c == '>')
808                                        break;                  /* EOL */
809                       
810                                c=chartab[c];
811                                if(c) seq[++(*len)]=c;
812                        }
813                        if(c == '>') break;
814                        }
815                break;
816/**********************************************/
817                case GDE:
818                        if (dnaflag) {
819                                while(*line != '#')
820                                        fgets(line,MAXLINE+1,fin);
821                        }
822                        else {
823                                while(*line != '%')
824                                        fgets(line,MAXLINE+1,fin);
825                        }
826                       
827                        for (i=1;i<=MAXNAMES;i++) {
828                                if (line[i] == '(' || line[i] == '\n')
829                                    break;
830                                sname[i-1] = line[i];
831                        }
832                        i--;
833                        sname[i]=EOS;
834                        if (sname[i-1] == '(') sscanf(&line[i],"%d",&offset);
835                        else offset = 0;
836                        for(i--;i > 0;i--) 
837                                if(isspace(sname[i])) {
838                                        sname[i]=EOS;   
839                                }
840                                else break;             
841                        blank_to_(sname);
842
843                        *tit=EOS;
844                       
845                        *len=0;
846                        for (i=0;i<offset;i++) seq[++(*len)] = '-';
847                        while(fgets(line,MAXLINE+1,fin)) {
848                        if(*line == '%' || *line == '#' || *line == '"') break;
849                                if(seq==NULL)
850                                        seq=(char *)ckalloc((MAXLINE+2)*sizeof(char));
851                                else
852                                        seq=(char *)ckrealloc(seq,((*len)+MAXLINE+2)*sizeof(char));
853                                for(i=0;i<=MAXLINE;i++) {
854                                        c=line[i];
855                                if(c == '\n' || c == EOS) 
856                                        break;                  /* EOL */
857                       
858                                c=chartab[c];
859                                if(c) seq[++(*len)]=c;
860                                }
861                        }
862                break;
863/***********************************************/
864                case RSF:
865                        while(*line != '{')
866                                if (fgets(line,MAXLINE+1,fin) == NULL) return NULL;
867                       
868                        while( !keyword(line,"name") )
869                                if (fgets(line,MAXLINE+1,fin) == NULL) return NULL;
870                       
871                        for(i=5;i<=strlen(line);i++)  /* DES */
872                                if(line[i] != ' ') break;
873                        strncpy(sname,line+i,MAXNAMES); /* remember entryname */
874                        for(i=0;i<=strlen(sname);i++)
875                                if(sname[i] == ' ') {
876                                        sname[i]=EOS;
877                                        break;
878                                }
879
880                        sname[MAXNAMES]=EOS;
881                        rtrim(sname);
882                        blank_to_(sname);
883
884                                               
885                        while( !keyword(line,"sequence") )
886                                if (fgets(line,MAXLINE+1,fin) == NULL) return NULL;
887                               
888                        *len=0;
889                        while(fgets(line,MAXLINE+1,fin)) {
890                                if(seq==NULL)
891                                        seq=(char *)ckalloc((MAXLINE+2)*sizeof(char));
892                                else
893                                        seq=(char *)ckrealloc(seq,((*len)+MAXLINE+2)*sizeof(char));
894                                for(i=0;i<=MAXLINE;i++) {
895                                        c=line[i];
896                                        if(c == EOS || c == '}')
897                                                break;                  /* EOL */
898                                        if( c=='.')
899                                                seq[++(*len)]='-';
900                                        c=chartab[c];
901                                        if(c) seq[++(*len)]=c;
902                                }
903                                if(c == '}') break;
904                        }
905                break;
906/***********************************************/
907        }
908       
909        seq[*len+1]=EOS;
910
911        return seq;
912}
913
914
915sint readseqs(sint first_seq) /*first_seq is the #no. of the first seq. to read */
916{
917        char line[FILENAMELEN+1];
918        char fileName[FILENAMELEN+1];
919
920        static char *seq1,sname1[MAXNAMES+1],title[MAXTITLES+1];
921        sint i,j;
922        sint no_seqs;
923        static sint l1;
924        static Boolean dnaflag1;
925       
926        if(usemenu)
927                getstr("Enter the name of the sequence file",line);
928        else
929                strcpy(line,seqname);
930        if(*line == EOS) return -1;
931
932        if ((sscanf(line,"file://%s",fileName) == 1 )) {
933          strcpy(line,fileName);
934        }
935
936        if((fin=fopen(line,"r"))==NULL) {
937                error("Could not open sequence file (%s) ",line);
938                return -1;      /* DES -1 => file not found */
939        }
940        strcpy(seqname,line);
941        no_seqs=0;
942        check_infile(&no_seqs);
943        info("Sequence format is %s",formatNames[seqFormat]);
944        if(seqFormat==NEXUS)
945                error("Cannot read nexus format");
946
947/* DES DEBUG
948        fprintf(stdout,"\n\n File name = %s\n\n",seqname);
949*/
950        if(no_seqs == 0)
951                return 0;       /* return the number of seqs. (zero here)*/
952
953/*
954        if((no_seqs + first_seq -1) > MAXN) {
955                error("Too many sequences. Maximum is %d",(pint)MAXN);
956                return 0;
957        }
958*/
959
960/* DES */
961/*      if(seqFormat == CLUSTAL) {
962                info("no of sequences = %d",(pint)no_seqs);
963                return no_seqs;
964        }
965*/
966        max_aln_length = 0;
967
968/* if this is a multiple alignment, or profile 1 - free any memory used
969by previous alignments, then allocate memory for the new alignment */
970        if(first_seq == 1) {
971                max_names = 0;
972                free_aln(nseqs);
973                alloc_aln(no_seqs);
974        }
975/* otherwise, this is a profile 2, and we need to reallocate the arrays,
976leaving the data for profile 1 intact */
977        else realloc_aln(first_seq,no_seqs);
978
979        for(i=1;i<first_seq;i++)
980        {
981                if(seqlen_array[i]>max_aln_length)
982                        max_aln_length=seqlen_array[i];
983                if(strlen(names[i])>max_names)
984                        max_names=strlen(names[i]);
985        }
986
987        for(i=first_seq;i<=first_seq+no_seqs-1;i++) {    /* get the seqs now*/
988                output_index[i] = i;    /* default output order */
989                if(seqFormat == CLUSTAL) 
990                        seq1=get_clustal_seq(sname1,&l1,title,i-first_seq+1);
991                else if(seqFormat == MSF)
992                            seq1=get_msf_seq(sname1,&l1,title,i-first_seq+1);
993                else
994                        seq1=get_seq(sname1,&l1,title);
995
996                if(seq1==NULL) break;
997/* JULIE */
998/*  Set max length of dynamically allocated arrays in prfalign.c */
999                if (l1 > max_aln_length) max_aln_length = l1;
1000                seqlen_array[i]=l1;                   /* store the length */
1001                strcpy(names[i],sname1);              /*    "   "  name   */
1002                strcpy(titles[i],title);              /*    "   "  title  */
1003
1004                if(!explicit_dnaflag) {
1005                        dnaflag1 = check_dnaflag(seq1,l1); /* check DNA/Prot */
1006                        if(i == 1) dnaflag = dnaflag1;
1007                }                       /* type decided by first seq*/
1008                else
1009                        dnaflag1 = dnaflag;
1010
1011                alloc_seq(i,l1);
1012
1013                if(dnaflag)
1014                        n_encode(seq1,seq_array[i],l1); /* encode the sequence*/
1015                else                                    /* as ints  */
1016                        p_encode(seq1,seq_array[i],l1);
1017                if(seq1!=NULL) seq1=ckfree(seq1);
1018        }
1019
1020
1021        max_aln_length *= 2;
1022/*
1023   JULIE
1024   check sequence names are all different - otherwise phylip tree is
1025   confused.
1026*/
1027        for(i=1;i<=first_seq+no_seqs-1;i++) {
1028                for(j=i+1;j<=first_seq+no_seqs-1;j++) {
1029                        if (strncmp(names[i],names[j],MAXNAMES) == 0) {
1030                                error("Multiple sequences found with same name, %s (first %d chars are significant)", names[i],MAXNAMES);
1031                                return 0;
1032                        }
1033                }
1034        }
1035        for(i=first_seq;i<=first_seq+no_seqs-1;i++)
1036        {
1037                if(seqlen_array[i]>max_aln_length)
1038                        max_aln_length=seqlen_array[i];
1039        }
1040       
1041/* look for a feature table / gap penalty mask (only if this is a profile) */
1042        if (profile_no > 0) {
1043                rewind(fin);
1044                struct_penalties = NONE;
1045                gap_penalty_mask = (char *)ckalloc((max_aln_length+1) * sizeof (char));
1046                sec_struct_mask = (char *)ckalloc((max_aln_length+1) * sizeof (char));
1047                ss_name = (char *)ckalloc((MAXNAMES+1) * sizeof (char));
1048
1049                if (seqFormat == CLUSTAL) {
1050                        get_clustal_ss(max_aln_length);
1051                }
1052                else if (seqFormat == GDE) {
1053                        get_gde_ss(max_aln_length);
1054                }
1055                else if (seqFormat == EMBLSWISS) {
1056                        get_embl_ss(max_aln_length);
1057                }
1058                else if (seqFormat == RSF) {
1059                        get_rsf_ss(max_aln_length);
1060                }
1061        }
1062
1063        for(i=first_seq;i<=first_seq+no_seqs-1;i++)
1064        {
1065                if(strlen(names[i])>max_names)
1066                        max_names=strlen(names[i]);
1067        }
1068
1069        if(max_names<10) max_names=10;
1070
1071        fclose(fin);
1072                       
1073        return no_seqs;    /* return the number of seqs. read in this call */
1074}
1075
1076
1077static Boolean check_dnaflag(char *seq, sint slen)
1078/* check if DNA or Protein
1079   The decision is based on counting all A,C,G,T,U or N.
1080   If >= 85% of all characters (except -) are as above => DNA  */
1081{
1082        sint i, c, nresidues, nbases;
1083        float ratio;
1084        char *dna_codes="ACGTUN";
1085       
1086        nresidues = nbases = 0; 
1087        for(i=1; i <= slen; i++) {
1088                if(seq[i] != '-') {
1089                        nresidues++;
1090                        if(seq[i] == 'N')
1091                                nbases++;
1092                        else {
1093                                c = res_index(dna_codes, seq[i]);
1094                                if(c >= 0)
1095                                        nbases++;
1096                        }
1097                }
1098        }
1099        if( (nbases == 0) || (nresidues == 0) ) return FALSE;
1100        ratio = (float)nbases/(float)nresidues;
1101/* DES  fprintf(stdout,"\n nbases = %d, nresidues = %d, ratio = %f\n",
1102                (pint)nbases,(pint)nresidues,(pint)ratio); */
1103        if(ratio >= 0.85) 
1104                return TRUE;
1105        else
1106                return FALSE;
1107}
1108
1109
1110
1111static void check_infile(sint *nseqs)
1112{
1113        char line[MAXLINE+1];
1114        sint i; 
1115
1116        *nseqs=0;
1117        while (fgets(line,MAXLINE+1,fin) != NULL) {
1118                if(!blankline(line)) 
1119                        break;
1120        }
1121
1122        for(i=strlen(line)-1;i>=0;i--)
1123                if(isgraph(line[i])) break;
1124        line[i+1]=EOS;
1125       
1126        for(i=0;i<=6;i++) line[i] = toupper(line[i]);
1127
1128        if( linetype(line,"ID") ) {                                     /* EMBL/Swiss-Prot format ? */
1129                seqFormat=EMBLSWISS;
1130                (*nseqs)++;
1131        }
1132        else if( linetype(line,"CLUSTAL") ) {
1133                seqFormat=CLUSTAL;
1134        }
1135        else if( linetype(line,"PILEUP") ) {
1136                seqFormat = MSF;
1137        }
1138        else if( linetype(line,"!!AA_MULTIPLE_ALIGNMENT") ) {
1139                seqFormat = MSF;
1140                dnaflag = FALSE;
1141        }
1142        else if( linetype(line,"!!NA_MULTIPLE_ALIGNMENT") ) {
1143                seqFormat = MSF;
1144                dnaflag = TRUE;
1145        }
1146        else if( strstr(line,"MSF") && line[strlen(line)-1]=='.' &&
1147                                 line[strlen(line)-2]=='.' ) {
1148                seqFormat = MSF;
1149        }
1150        else if( linetype(line,"!!RICH_SEQUENCE") ) {
1151                seqFormat = RSF;
1152        }
1153        else if( linetype(line,"#NEXUS") ) {
1154                seqFormat=NEXUS;
1155                return;
1156        }
1157        else if(*line == '>') {                                         /* no */
1158                seqFormat=(line[3] == ';')?PIR:PEARSON; /* distinguish PIR and Pearson */
1159                (*nseqs)++;
1160        }
1161        else if((*line == '"') || (*line == '%') || (*line == '#')) {
1162                seqFormat=GDE; /* GDE format */
1163                if (*line == '%') {
1164                        (*nseqs)++;
1165                        dnaflag = FALSE;
1166                }
1167                else if (*line == '#') {
1168                        (*nseqs)++;
1169                        dnaflag = TRUE;
1170                }
1171        }
1172        else {
1173                seqFormat=UNKNOWN;
1174                return;
1175        }
1176
1177        while(fgets(line,MAXLINE+1,fin) != NULL) {
1178                switch(seqFormat) {
1179                        case EMBLSWISS:
1180                                if( linetype(line,"ID") )
1181                                        (*nseqs)++;
1182                                break;
1183                        case PIR:
1184                                *nseqs = count_pir_seqs();
1185                                fseek(fin,0,0);
1186                                return;
1187                        case PEARSON:
1188                                if( *line == '>' )
1189                                        (*nseqs)++;
1190                                break;
1191                        case GDE:
1192                                if(( *line == '%' ) && ( dnaflag == FALSE))
1193                                        (*nseqs)++;
1194                                else if (( *line == '#') && ( dnaflag == TRUE))
1195                                        (*nseqs)++;
1196                                break;
1197                        case CLUSTAL:
1198                                *nseqs = count_clustal_seqs();
1199/* DES */                       /* fprintf(stdout,"\nnseqs = %d\n",(pint)*nseqs); */
1200                                fseek(fin,0,0);
1201                                return;
1202                        case MSF:
1203                                *nseqs = count_msf_seqs();
1204                                fseek(fin,0,0);
1205                                return;
1206                        case RSF:
1207                                fseek(fin,0,0);
1208                                *nseqs = count_rsf_seqs();
1209                                fseek(fin,0,0);
1210                                return;
1211                        case USER:
1212                        default:
1213                                break;
1214                }
1215        }
1216        fseek(fin,0,0);
1217}
1218
1219
1220static sint count_pir_seqs(void)
1221/* count the number of sequences in a pir alignment file */
1222{
1223        char line[MAXLINE+1],c;
1224        sint  nseqs, i;
1225        Boolean seq_ok;
1226
1227        seq_ok = FALSE;
1228        while (fgets(line,MAXLINE+1,fin) != NULL) { /* Look for end of first seq */
1229                if(*line == '>') break;
1230                for(i=0;seq_ok == FALSE;i++) {
1231                        c=line[i];
1232                        if(c == '*') {
1233                                seq_ok = TRUE;  /* ok - end of sequence found */
1234                                break;
1235                        }                       /* EOL */
1236                        if(c == '\n' || c == EOS)
1237                                break;                  /* EOL */
1238                }
1239                if (seq_ok == TRUE)
1240                        break;
1241        }
1242        if (seq_ok == FALSE) {
1243                error("PIR format sequence end marker '*'\nmissing for one or more sequences.");
1244                return (sint)0; /* funny format*/
1245        }
1246       
1247       
1248        nseqs = 1;
1249       
1250        while (fgets(line,MAXLINE+1,fin) != NULL) {
1251                if(*line == '>') {              /* Look for start of next seq */
1252                        seq_ok = FALSE;
1253                        while (fgets(line,MAXLINE+1,fin) != NULL) { /* Look for end of seq */
1254                                if(*line == '>') {
1255                                        error("PIR format sequence end marker '*' missing for one or more sequences.");
1256                                        return (sint)0; /* funny format*/
1257                                }
1258                                for(i=0;seq_ok == FALSE;i++) {
1259                                        c=line[i];
1260                                        if(c == '*') {
1261                                                seq_ok = TRUE;  /* ok - sequence found */
1262                                                break;
1263                                        }                       /* EOL */
1264                                        if(c == '\n' || c == EOS)
1265                                                break;                  /* EOL */
1266                                }
1267                                if (seq_ok == TRUE) {
1268                                        nseqs++;
1269                                        break;
1270                                }
1271                        }
1272                }
1273        }
1274        return (sint)nseqs;
1275}
1276
1277
1278static sint count_clustal_seqs(void)
1279/* count the number of sequences in a clustal alignment file */
1280{
1281        char line[MAXLINE+1];
1282        sint  nseqs;
1283
1284        while (fgets(line,MAXLINE+1,fin) != NULL) {
1285                if(!cl_blankline(line)) break;          /* Look for next non- */
1286        }                                               /* blank line */
1287        nseqs = 1;
1288
1289        while (fgets(line,MAXLINE+1,fin) != NULL) {
1290                if(cl_blankline(line)) return nseqs;
1291                nseqs++;
1292        }
1293
1294        return (sint)0; /* if you got to here-funny format/no seqs.*/
1295}
1296
1297static sint count_msf_seqs(void)
1298{
1299/* count the number of sequences in a PILEUP alignment file */
1300
1301        char line[MAXLINE+1];
1302        sint  nseqs;
1303
1304        while (fgets(line,MAXLINE+1,fin) != NULL) {
1305                if(linetype(line,"//")) break;
1306        }
1307
1308        while (fgets(line,MAXLINE+1,fin) != NULL) {
1309                if(!blankline(line)) break;             /* Look for next non- */
1310        }                                               /* blank line */
1311        nseqs = 1;
1312
1313        while (fgets(line,MAXLINE+1,fin) != NULL) {
1314                if(blankline(line)) return nseqs;
1315                nseqs++;
1316        }
1317
1318        return (sint)0; /* if you got to here-funny format/no seqs.*/
1319}
1320
1321static sint count_rsf_seqs(void)
1322{
1323/* count the number of sequences in a GCG RSF alignment file */
1324
1325        char line[MAXLINE+1];
1326        sint  nseqs;
1327
1328        nseqs = 0;
1329/* skip the comments */
1330        while (fgets(line,MAXLINE+1,fin) != NULL) {
1331                if(line[strlen(line)-2]=='.' &&
1332                                 line[strlen(line)-3]=='.')
1333                        break;
1334        }
1335
1336        while (fgets(line,MAXLINE+1,fin) != NULL) {
1337                if( *line == '{' )
1338                      nseqs++;
1339        }
1340        return (sint)nseqs;
1341}
1342
1343static void p_encode(char *seq, char *naseq, sint l)
1344{                               /* code seq as ints .. use gap_pos2 for gap */
1345        register sint i;
1346/*      static char *aacids="CSTPAGNDEQHRKMILVFYW";*/
1347       
1348        for(i=1;i<=l;i++)
1349                if(seq[i] == '-')
1350                        naseq[i] = gap_pos2;
1351                else
1352                        naseq[i] = res_index(amino_acid_codes,seq[i]);
1353        naseq[i] = -3;
1354}
1355
1356static void n_encode(char *seq,char *naseq,sint l)
1357{                               /* code seq as ints .. use gap_pos2 for gap */
1358        register sint i;
1359/*      static char *nucs="ACGTU";      */
1360       
1361        for(i=1;i<=l;i++) {
1362        if(seq[i] == '-')                  /* if a gap character -> code = gap_pos2 */
1363                        naseq[i] = gap_pos2;   /* this is the code for a gap in */
1364                else {                     /* the input files */
1365                        naseq[i]=res_index(amino_acid_codes,seq[i]);
1366                }
1367        }
1368        naseq[i] = -3;
1369}
1370
1371static sint res_index(char *t,char c)
1372{
1373        register sint i;
1374       
1375        for(i=0;t[i] && t[i] != c;i++)
1376                ;
1377        if(t[i]) return(i);
1378        else return -1;
1379}
Note: See TracBrowser for help on using the repository browser.