binaries/src/fasta34/getseq.c

   1 /*      May, June 1987  - modified for rapid read of database
   2
   3         copyright (c) 1987,1988,1989,1992,1995,2000 William R. Pearson
   4
   5         This is one of three alternative files that can be used to
   6         read a database.  The three files are nxgetaa.c, nmgetaa.c, and
   7         mmgetaa.c.
   8
   9         nxgetaa.c contains the original code for reading databases, and
  10         is still used for Mac and PC versions of fasta33 (which do not
  11         use mmap).
  12
  13         nmgetaa.c and mmgetaa.c are used together.  nmgetaa.c provides
  14         the same functions as nxgetaa.c if memory mapping is not used,
  15         mmgetaa.c provides the database reading functions if memory
  16         mapping is used. The decision to use memory mapping is made on
  17         a file-by-file basis.
  18
  19         June 2, 1987 - added TFASTA
  20         March 30, 1988 - combined ffgetaa, fgetgb;
  21         April 8, 1988 - added PIRLIB format for unix
  22         Feb 4, 1989 - added universal subroutines for libraries
  23         December, 1995 - added range option file.name:1-1000
  24         Feb 22, 2002 - fix to allow "plain" text file queries
  25
  26         getnt.c associated subroutines for matching sequences */
  27
  28 /* $Name: fa_34_26_5 $ - $Id: getseq.c,v 1.13 2006/10/05 18:22:07 wrp Exp $ */
  29
  30 /*
  31         8-April-88
  32         The compile time #define PIRLIB allows this routine to be used
  33         to read protein and DNA sequence libraries in the NBRF/PIR
  34         VAX/VMS library format.  That is:
  35
  36         >P1;LCBO
  37         This is a line of description
  38         GTYH ... the sequence starts on this line
  39
  40         This may ease conversion from UWGCG format libraries. It
  41         has not been extensively tested.
  42
  43         In addition, sequence libraries with a '>' in the 4th position
  44         are recognized as NBRF format libraries for consistency with
  45         UWGCG
  46 */
  47
  48 /*      Nov 12, 1987    - this version checks to see if the sequence
  49         is DNA or protein by asking whether > 85% is A, C, G, T
  50
  51         May 5, 1988 - modify the DNA/PROTEIN checker by re-reading
  52         DNA sequences in order to check for 'U'.
  53 */
  54
  55 #include <stdio.h>
  56 #include <stdlib.h>
  57 #include <string.h>
  58
  59 #include "defs.h"
  60 #include "structs.h"
  61
  62 #ifndef SFCHAR
  63 #define SFCHAR ':'
  64 #endif
  65
  66 #ifdef VMS
  67 #define PIRLIB
  68 #endif
  69
  70 #define XTERNAL
  71 #include "uascii.h"
  72 #include "upam.h"
  73 #undef XTERNAL
  74
  75 #define YES 1
  76 #define NO 0
  77 #define MAXLINE 512
  78
  79 #ifndef min
  80 #define min(x,y) ((x) > (y) ? (y) : (x))
  81 #endif
  82
  83 #ifdef SUPERFAMNUM
  84 extern int nsfnum;      /* number of superfamily numbers */
  85 extern int sfnum[];     /* superfamily number from types 0 and 5 */
  86 extern int nsfnum_n;
  87 extern int sfnum_n[];
  88 #endif
  89
  90 #define NO_FORMAT 0
  91 #define FASTA_FORMAT 1
  92 #define GCG_FORMAT 2
  93
  94 static int seq_format=NO_FORMAT;
  95 static char seq_title[200];
  96
  97 int scanseq(unsigned char *, int, char *);
  98 void sf_sort(int *, int);
  99 extern void init_ascii(int is_ext, int *sascii, int is_dna);
 100
 101 /* getseq       - get a query sequence, possibly re-reading to set type
 102    returns      - length of query sequence or error = 0
 103
 104    char *filen  - name of file to be opened
 105    char *seq    - destination for query sequence
 106    int maxs     - maximum length of query
 107    char libstr[20] - short description (locus or acc)
 108    int *dnaseq  - -1 => use scanseq to determine sequence type
 109                    0 => must be protein
 110                    1 => must be DNA
 111    long *sq0off - offset into query specified by query_file:1001-2000
 112 */
 113
 114 int
 115 getseq(char *filen, int *qascii, unsigned char *seq, int maxs, char *libstr, long *sq0off)
 116 {
 117   FILE *fptr;
 118   char line[512],*bp, *bp1, *bpn, *tp;
 119   int i, rn, n;
 120   int ic;
 121   int sstart, sstop, sset=0;
 122   int llen, l_offset;
 123 #ifdef SUPERFAMNUM
 124   static char tline[MAXLINE];
 125 #endif
 126
 127   seq_title[0]='\0';
 128   libstr[0]='\0';
 129
 130   sstart = sstop = -1;
 131 #ifndef DOS
 132   if ((bp=strchr(filen,':'))!=NULL && *(bp+1)!='\0') {
 133 #else
 134   if ((bp=strchr(filen+3,':'))!=NULL && *(bp+1)!='\0') {
 135 #endif
 136     *bp='\0';
 137     if (*(bp+1)=='-') {
 138       sstart = 0;
 139       sscanf(bp+2,"%d",&sstop);
 140     }
 141     else {
 142       sscanf(bp+1,"%d-%d",&sstart,&sstop);
 143       sstart--;
 144       if (sstop <= 0 ) sstop = BIGNUM;
 145     }
 146     sset=1;
 147   }
 148   else {
 149     sstart = 0;
 150     sstop = BIGNUM;
 151   }
 152
 153   /* check for input from stdin */
 154   if (strcmp(filen,"-") && strcmp(filen,"@")) {
 155     if ((fptr=fopen(filen,"r"))==NULL) {
 156       fprintf(stderr," could not open %s\n",filen);
 157       return 0;
 158     }
 159   }
 160   else {
 161     fptr = stdin;
 162   }
 163   rn = n=0;
 164
 165   while(fgets(line,sizeof(line),fptr)!=NULL) {
 166 #ifdef PIRLIB
 167     if (line[0]=='>'&& (line[3]==';'||line[3]=='>'))
 168       fgets(line,sizeof(line),fptr);
 169     else
 170 #endif
 171     l_offset = 0;
 172     if (line[0]=='>') {
 173       seq_format = FASTA_FORMAT;
 174 #ifdef SUPERFAMNUM
 175       sfnum[nsfnum=0]= sfnum_n[nsfnum_n=0]=0;
 176       strncpy(tline,line+1,sizeof(tline));
 177       tline[sizeof(tline)-1]='\0';
 178
 179       if ((bp=strchr(tline,' ')) && (bp=strchr(bp+1,SFCHAR))) {
 180         if ((bp1=strchr(bp+1,SFCHAR))==NULL) {
 181           fprintf(stderr," second %c missing: %s\n",SFCHAR,tline);
 182         }
 183         else {
 184           if ((bpn=strchr(bp+1,NSFCHAR))!=NULL) *bpn = '\0';
 185           *bp1 = '\0';
 186           i = 0;
 187           if ((tp = strtok(bp+1," \t"))!=NULL)  {
 188             sfnum[i++] = atoi(tp);
 189             while ((tp = strtok((char *)NULL," \t")) != (char *)NULL) {
 190               if (isdigit(*tp)) sfnum[i++] = atoi(tp);
 191               if (i>=9) break;
 192             }
 193           }
 194           sfnum[nsfnum=i]= 0;
 195           if (nsfnum>1) sf_sort(sfnum,nsfnum);
 196           else {
 197             if (nsfnum < 1) fprintf(stderr," found | but no sfnum: %s\n",libstr);
 198           }
 199           if (bpn != NULL) {
 200             tp = strtok(bpn+1," \t");
 201             sfnum_n[0]=atoi(tp);
 202             i = 1;
 203             while ((tp=strtok(NULL," \t"))!=NULL) {
 204               sfnum_n[i++] = atoi(tp);
 205               if (i >= 10) {
 206                 fprintf(stderr,
 207                         " error - too many negative superfamilies: %d\n %s\n",
 208                         i,tline);
 209                 break;
 210               }
 211             }
 212             sfnum_n[nsfnum_n=i]=0;
 213             sf_sort(sfnum_n,nsfnum_n);
 214           }
 215         }
 216       }
 217       else {
 218         sfnum[nsfnum = 0] = 0;
 219         sfnum_n[nsfnum_n = 0] = 0;
 220       }
 221 #endif
 222       if ((bp=(char *)strchr(line,'\n'))!=NULL) *bp='\0';
 223       strncpy(seq_title,line+1,sizeof(seq_title));
 224       seq_title[sizeof(seq_title)-1]='\0';
 225       if ((bp=(char *)strchr(line,' '))!=NULL) *bp='\0';
 226       strncpy(libstr,line+1,12);
 227       libstr[12]='\0';
 228     }
 229     else if (seq_format==NO_FORMAT && strcmp(line,"..")==0) {
 230       seq_format = GCG_FORMAT;
 231 /*
 232       if (*dnaseq != 1) qascii['*'] = qascii['X'];
 233 */
 234       l_offset = 10;
 235       llen = strlen(line);
 236       while (strncmp(&line[llen-3],"..\n",(size_t)3) != 0) {
 237         if (fgets(line,sizeof(line),fptr)==NULL) return 0;
 238         llen = strlen(line);
 239       }
 240       bp = strtok(line," \t");
 241 /*
 242       if ((bp=(char *)strchr(line,' '))!=NULL) *bp='\0';
 243       else if ((bp=(char *)strchr(line,'\n'))!=NULL) *bp='\0';
 244 */
 245       if (bp!=NULL) strncpy(libstr,bp,12);
 246       else strncpy(libstr,filen,12);
 247       libstr[12]='\0';
 248       if (fgets(line,sizeof(line),fptr)==NULL) return 0;
 249     }
 250     else {
 251       if (libstr[0]=='\0') strncpy(libstr,filen,12);
 252       libstr[12]='\0';
 253     }
 254
 255     if (seq_format==GCG_FORMAT && strlen(line)<l_offset) continue;
 256
 257     if (line[0]!='>'&& line[0]!=';') {
 258       for (i=l_offset; (n<maxs && rn < sstop)&&
 259              ((ic=qascii[line[i]&AAMASK])<EL); i++)
 260         if (ic<NA && ++rn > sstart) seq[n++]= ic;
 261       if (ic == ES || rn > sstop) break;
 262     }
 263   }
 264
 265   if (n==maxs) {
 266     fprintf(stderr," sequence may be truncated %d %d\n",n,maxs);
 267     fflush(stderr);
 268   }
 269   if ((bp=strchr(libstr,'\n'))!=NULL) *bp = '\0';
 270   if ((bp=strchr(libstr,'\r'))!=NULL) *bp = '\0';
 271   seq[n]= EOSEQ;
 272
 273
 274   if (seq_format !=GCG_FORMAT)
 275     while(fgets(line,sizeof(line),fptr)!=NULL) {
 276 #ifdef PIRLIB
 277       if (line[0]=='>'&& (line[3]==';'||line[3]=='>'))
 278         fgets(line,sizeof(line),fptr);
 279       else
 280 #endif
 281         if (line[0]!='>'&& line[0]!=';') {
 282           for (i=0; (n<maxs && rn < sstop)&&
 283                  ((ic=qascii[line[i]&AAMASK])<EL); i++)
 284             if (ic<NA && ++rn > sstart ) seq[n++]= ic;
 285           if (ic == ES || rn > sstop) break;
 286         }
 287     }
 288   else {
 289     llen = strlen(line);
 290     while (strncmp(&line[llen-3],"..\n",(size_t)3) != 0) {
 291       if (fgets(line,sizeof(line),fptr)==NULL) return 0;
 292       llen = strlen(line);
 293     }
 294     while (fgets(line,sizeof(line),fptr)!=NULL) {
 295       if (strlen(line)<l_offset) continue;
 296       for (i=l_offset; (n<maxs && rn < sstop) &&
 297              ((ic=qascii[line[i]&AAMASK])<EL); i++)
 298         if (ic<NA && ++rn > sstart ) seq[n++]= ic;
 299       if (ic == ES || rn > sstop ) break;
 300     }
 301   }
 302
 303   if (n==maxs) {
 304     fprintf(stderr," sequence may be truncated %d %d\n",n,maxs);
 305     fflush(stderr);
 306   }
 307   seq[n]= EOSEQ;
 308
 309   if (fptr!=stdin) fclose(fptr);
 310
 311   if (sset==1) {
 312     sstart++;
 313     filen[strlen(filen)]=':';
 314     if (*sq0off==1 || sstart>=1) *sq0off = sstart;
 315   }
 316
 317   return n;
 318 }
 319
 320 int
 321 gettitle(char *filen, char *title, int len) {
 322   FILE *fptr;
 323   char line[512];
 324   char *bp;
 325   int sset;
 326 #ifdef WIN32
 327   char *strpbrk();
 328 #endif
 329
 330   sset = 0;
 331
 332   if (strncmp(filen,"-",1)==0 || strncmp(filen,"@",1)==0) {
 333     strncpy(title,seq_title,len);
 334     title[len-1]='\0';
 335     return (int)strlen(title);
 336   }
 337
 338   if ((bp=strchr(filen,':'))!=NULL) { *bp='\0'; sset=1;}
 339
 340
 341   if ((fptr=fopen(filen,"r"))==NULL) {
 342     fprintf(stderr," file %s was not found\n",filen);
 343     fflush(stderr);
 344     return 0;
 345   }
 346
 347   if (sset==1) filen[strlen(filen)]=':';
 348
 349   while(fgets(line,sizeof(line),fptr)!=NULL) {
 350     if (line[0]=='>'|| line[0]==';') goto found;
 351   }
 352   fclose(fptr);
 353   title[0]='\0';
 354   return 0;
 355
 356  found:
 357
 358 #ifdef PIRLIB
 359   if (line[0]=='>'&&(line[3]==';'||line[3]=='>')) {
 360     if ((bp = strchr(line,'\n'))!=NULL) *bp='\0';
 361     ll=strlen(line); line[ll++]=' '; line[ll]='\0';
 362     fgets(&line[ll],sizeof(line)-ll,fptr);
 363   }
 364 #endif
 365
 366 #ifdef WIN32
 367   bp = strpbrk(line,"\n\r");
 368 #else
 369   bp = strchr(line,'\n');
 370 #endif
 371   if (bp!=NULL) *bp = 0;
 372   strncpy(title,line,len);
 373   title[len-1]='\0';
 374   fclose(fptr);
 375   return strlen(title);
 376 }
 377