binaries/src/fasta34/ncbl_lib.c

   1 /*      ncbl_lib.c      functions to read ncbi-blast format files from
   2                         setdb (blastp 1.3.2) format files
   3
   4                 copyright (c) 1992 William R. Pearson
   5 */
   6
   7 /* $Name: fa_34_26_5 $ - $Id: ncbl_lib.c,v 1.9 2006/10/05 18:22:07 wrp Exp $ */
   8
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <string.h>
  12
  13 #ifndef WIN32
  14 #define RBSTR "r"
  15 #else
  16 #define RBSTR "rb"
  17 #endif
  18
  19 #define XTERNAL
  20 #include "uascii.h"
  21 #include "upam.h"
  22 #include "ncbl_head.h"
  23 #include "mm_file.h"
  24
  25 int ncbl_getliba(char *, int, char *, int, fseek_t *, int *, struct lmf_str *, long *);
  26 int ncbl_getlibn(char *, int, char *, int, fseek_t *, int *, struct lmf_str *, long *);
  27
  28 void src_ulong_read();
  29
  30 #ifndef NCBL13_ONLY
  31 static void src_char_read();
  32 static void newname(char *, char *, char *, int);
  33 #else
  34 void src_char_read();
  35 void newname(char *, char *, char *, int);
  36 #endif
  37
  38 /* nt_btoa maps  from blast 2bit  format to ascii  characters */
  39 static char nt_btoa[5] = {"ACGT"};
  40
  41 static char aa_btoa[27]= {"-ARNDCQEGHILKMFPSTWYVBZX*"};
  42 static int aa_btof[32]; /* maps to fasta alphabet */
  43
  44 static FILE *tfile=NULL,        /* table of offsets, also DB info */
  45             *hfile=NULL,        /* description lines */
  46             *sfile=NULL;        /* binary sequence data */
  47
  48 static unsigned long lib_cnt, max_cnt, totlen, mxlen, dbline_len;
  49 static unsigned long *seq_beg, *hdr_beg;
  50 static unsigned char *ambiguity_ray;
  51 static long seq_format, dbtype, dbformat;
  52 static char dline[512];
  53
  54 #define NCBIBL13 11
  55
  56 struct lmf_str *
  57 ncbl_openlib(char *name, int ldnaseq)
  58 {
  59   char hname[256];
  60   char sname[256];
  61   char tname[256];
  62   long title_len;
  63   char *title_str;
  64   int rdtmp;
  65   int i;
  66   unsigned long line_len, c_len, clean_count;
  67
  68   if (ldnaseq!=1) {
  69     newname(tname,name,AA_TABLE_EXT,(int)sizeof(tname));
  70     if ((tfile = fopen(tname,RBSTR))==NULL) {
  71       fprintf(stderr," cannot open %s (%s.%s) table file\n",
  72               name,tname,NT_TABLE_EXT);
  73       return (-1);
  74     }
  75     seq_format = AAFORMAT;
  76   }
  77   else {
  78     newname(tname,name,NT_TABLE_EXT,(int)sizeof(tname));
  79     if ((tfile = fopen(tname,RBSTR))==NULL) {
  80       fprintf(stderr," cannot open %s (%s.%s) table file\n",
  81               name,tname,NT_TABLE_EXT);
  82       return (-1);
  83     }
  84     seq_format = NTFORMAT;
  85   }
  86
  87   src_ulong_read(tfile,&dbtype);
  88   src_ulong_read(tfile,&dbformat);
  89
  90   if (seq_format == AAFORMAT && (dbformat != seq_format || dbtype !=
  91                                  DB_TYPE_PRO)) {
  92     fprintf(stderr,"error - %s wrong type (%ld/%d) or format (%ld/%ld)\n",
  93             tname,dbtype,DB_TYPE_PRO,dbformat,seq_format);
  94     return (-1);
  95   }
  96   else if (seq_format == NTFORMAT && (dbformat != seq_format || dbtype !=
  97                                  DB_TYPE_NUC)) {
  98     fprintf(stderr,"error - %s wrong type (%ld/%d) or format (%ld/%ld)\n",
  99             tname,dbtype,DB_TYPE_NUC,dbformat,seq_format);
 100     return (-1);
 101   }
 102
 103   if (seq_format == AAFORMAT) {
 104     newname(hname,name,AA_HEADER_EXT,(int)sizeof(hname));
 105     if ((hfile = fopen(hname,RBSTR))==NULL) {
 106       fprintf(stderr," cannot open %s header file\n",hname);
 107       return (-1);
 108     }
 109     newname(sname,name,AA_SEARCHSEQ_EXT,(int)sizeof(sname));
 110     if ((sfile = fopen(sname,RBSTR))==NULL) {
 111       fprintf(stderr," cannot open %s sequence file\n",sname);
 112       return (-1);
 113     }
 114   }
 115   else {
 116     newname(hname,name,NT_HEADER_EXT,(int)sizeof(hname));
 117     if ((hfile = fopen(hname,RBSTR))==NULL) {
 118       fprintf(stderr," cannot open %s header file\n",hname);
 119       return (-1);
 120     }
 121     newname(sname,name,NT_SEARCHSEQ_EXT,(int)sizeof(sname));
 122     if ((sfile = fopen(sname,RBSTR))==NULL) {
 123       fprintf(stderr," cannot open %s sequence file\n",sname);
 124       return (-1);
 125     }
 126   }
 127
 128 /* all files should be open */
 129
 130   src_ulong_read(tfile,&title_len);
 131   rdtmp = title_len + ((title_len%4 !=0 ) ? 4-(title_len%4) : 0);
 132   if ((title_str = calloc((size_t)rdtmp,sizeof(char)))==NULL) {
 133     fprintf(stderr," cannot allocate title string (%d)\n",rdtmp);
 134     return(-1);
 135   }
 136   fread(title_str,(size_t)1,(size_t)rdtmp,tfile);
 137
 138   lib_cnt = 0;
 139   if (seq_format == AAFORMAT) {
 140     src_ulong_read(tfile,&max_cnt);
 141     src_ulong_read(tfile,&totlen);
 142     src_ulong_read(tfile,&mxlen);
 143
 144     /* fprintf(stderr," max_cnt: %d, totlen: %d\n",max_cnt,totlen); */
 145
 146     if ((seq_beg=(unsigned long *)calloc((size_t)max_cnt+1,sizeof(long)))==NULL) {
 147       fprintf(stderr," cannot allocate sequence pointers\n");
 148       return -1;
 149     }
 150     if ((hdr_beg=(unsigned long *)calloc((size_t)max_cnt+1,sizeof(long)))==NULL) {
 151       fprintf(stderr," cannot allocate header pointers\n");
 152       return -1;
 153     }
 154     for (i=0; i<max_cnt+1; i++) src_ulong_read(tfile,&seq_beg[i]);
 155     for (i=0; i<max_cnt+1; i++) src_ulong_read(tfile,&hdr_beg[i]);
 156
 157     for (i=0; i<sizeof(aa_btoa); i++) {
 158       if ((rdtmp=aascii[aa_btoa[i]])<NA) aa_btof[i]=rdtmp;
 159       else aa_btof[i]=aascii['X'];
 160     }
 161   }
 162   else if (seq_format == NTFORMAT) {
 163     src_ulong_read(tfile,&dbline_len);  /* length of uncompress DB lines */
 164     src_ulong_read(tfile,&max_cnt);     /* number of entries */
 165     src_ulong_read(tfile,&mxlen);       /* maximum length sequence */
 166     src_ulong_read(tfile,&totlen);      /* total count */
 167     src_ulong_read(tfile,&c_len);       /* compressed db length */
 168     src_ulong_read(tfile,&clean_count); /* count of nt's cleaned */
 169
 170     fseek(tfile,(size_t)((clean_count)*4),1);
 171                                          /* seek over clean_count */
 172     if ((seq_beg=(unsigned long *)calloc((size_t)max_cnt+1,sizeof(long)))==NULL) {
 173       fprintf(stderr," cannot allocate sequence pointers\n");
 174       return -1;
 175     }
 176     if ((hdr_beg=(unsigned long *)calloc((size_t)max_cnt+1,sizeof(long)))==NULL) {
 177       fprintf(stderr," cannot allocate header pointers\n");
 178       return -1;
 179     }
 180     if ((ambiguity_ray=
 181          (unsigned char *)calloc((size_t)max_cnt/CHAR_BIT+1,sizeof(char)))==NULL) {
 182       fprintf(stderr," cannot allocate ambiguity_ray\n");
 183       return -1;
 184     }
 185
 186     for (i=0; i<max_cnt+1; i++) src_ulong_read(tfile,&seq_beg[i]);
 187     fseek(tfile,(size_t)((max_cnt+1)*4),1);
 188                                          /* seek over seq_beg */
 189     for (i=0; i<max_cnt+1; i++) src_ulong_read(tfile,&hdr_beg[i]);
 190     for (i=0; i<max_cnt/CHAR_BIT+1; i++)
 191       src_char_read(tfile,&ambiguity_ray[i]);
 192   }
 193   return 1;
 194 }
 195
 196 void ncbl_closelib()
 197 {
 198   if (tfile !=NULL ) {fclose(tfile); tfile=NULL;}
 199   if (hfile !=NULL ) {fclose(hfile); hfile=NULL;}
 200   if (sfile !=NULL ) {fclose(sfile); sfile=NULL;}
 201 }
 202
 203 int
 204 ncbl_getliba(char *seq, int maxs,
 205              char *libstr, int n_libstr,
 206              fseek_t *libpos,
 207              int lcont)
 208 {
 209   register char *sptr;
 210   long seqcnt;
 211   long tmp;
 212   char ch;
 213   static long seq_len;
 214
 215   *libpos = lib_cnt;
 216   if (*lcont==0) {
 217     if (lib_cnt >= max_cnt) return -1;
 218     seq_len = seq_beg[lib_cnt+1] - seq_beg[lib_cnt] -1;
 219     tmp=(long)fgetc(sfile);     /* skip the null byte */
 220     if (tmp!=NULLB)
 221       fprintf(stderr," phase error: %ld:%ld found\n",lib_cnt,tmp);
 222     libstr[0]='\0';
 223     }
 224
 225   if (seq_len < maxs) {
 226     if ((tmp=fread(seq,(size_t)1,(size_t)seq_len,sfile))!=(size_t)seq_len) {
 227       fprintf(stderr," could not read sequence record: %ld %ld != %ld\n",
 228               *libpos,tmp,seq_len);
 229       goto error;
 230     }
 231     if (aa_btoa[seq[seq_len-1]]=='*') seqcnt = seq_len-1;
 232     else seqcnt=seq_len;
 233     lib_cnt++;
 234     *lcont = 0;
 235   }
 236   else {
 237     if (fread(seq,(size_t)1,(size_t)(maxs-1),sfile)!=(size_t)(maxs-1)) {
 238       fprintf(stderr," could not read sequence record: %ld %ld\n",
 239               *libpos,seq_len);
 240       goto error;
 241     }
 242     (*lcont)++;
 243     seqcnt = maxs-1;
 244     seq_len -= seqcnt;
 245   }
 246   sptr = seq+seqcnt;
 247
 248   while (--sptr >= seq) *sptr = aa_btof[*sptr];
 249
 250   seq[seqcnt]= EOSEQ;
 251   return (seqcnt);
 252
 253 error:  fprintf(stderr," error reading %ld at %ld\n",libstr,*libpos);
 254   fflush(stderr);
 255   return (-1);
 256 }
 257
 258 int
 259 ncbl_getlibn(char *seq, int maxs,
 260              char *libstr, int n_libstr,
 261              fseek_t *libpos, int *lcont)
 262 {
 263   register char *sptr, *tptr, stmp;
 264   long seqcnt;
 265   long tmp;
 266   char ch;
 267   static long seq_len;
 268   static int c_len,c_pad;
 269
 270   *libpos = lib_cnt;
 271   if (*lcont==0) {
 272     if (lib_cnt >= max_cnt) return -1;
 273     c_len = seq_beg[lib_cnt+1]/(CHAR_BIT/NBPN)
 274                 - seq_beg[lib_cnt]/(CHAR_BIT/NBPN);
 275     c_len -= NSENTINELS;
 276
 277     seq_len = c_len*(CHAR_BIT/NBPN);
 278     c_pad = seq_beg[lib_cnt] & ((CHAR_BIT/NBPN)-1);
 279     if (c_pad != 0) seq_len -= ((CHAR_BIT/NBPN) - c_pad);
 280
 281     tmp=fgetc(sfile);   /* skip the null byte */
 282     if (tmp!=NT_MAGIC_BYTE) {
 283       fprintf(stderr," phase error: %ld:%ld (%ld/%d) found\n",
 284               lib_cnt,seq_len,tmp,NT_MAGIC_BYTE);
 285       goto error;
 286     }
 287     libstr[0]='\0';
 288   }
 289
 290   if (seq_len < maxs-3) {
 291     seqcnt=(seq_len+3)/4;
 292     if (seqcnt==0) seqcnt++;
 293     if ((tmp=fread(seq,(size_t)1,(size_t)seqcnt,sfile))
 294         !=(size_t)seqcnt) {
 295       fprintf(stderr,
 296               " could not read sequence record: %s %ld %ld != %ld: %d\n",
 297               libstr,*libpos,tmp,seqcnt,*seq);
 298       goto error;
 299     }
 300     tmp=fgetc(sfile);   /* skip the null byte */
 301     if (tmp!=(unsigned char)NT_MAGIC_BYTE) {
 302       fprintf(stderr," phase2 error: %ld:%ld (%ld/%d) next ",
 303               lib_cnt,seqcnt,tmp,NT_MAGIC_BYTE);
 304
 305       goto error;
 306     }
 307     *lcont = 0;
 308     lib_cnt++;
 309   }
 310   else {
 311     seqcnt = ((maxs+3)/4)-1;
 312     if (fread(seq,(size_t)1,(size_t)(seqcnt),sfile)!=(size_t)(seqcnt)) {
 313       fprintf(stderr," could not read sequence record: %s %ld %ld\n",
 314               libstr,*libpos,seqcnt);
 315       goto error;
 316     }
 317     (*lcont)++;
 318   }
 319
 320   /* point to the last packed byte and to the end of the array
 321      seqcnt is the exact number of bytes read
 322      tptr points to the destination, use multiple of 4 to simplify math
 323      sptr points to the source, note that the last byte will be read 4 cycles
 324      before it is written
 325      */
 326
 327   sptr = seq + seqcnt;
 328   tptr = seq + 4*seqcnt;
 329   while (sptr>seq) {
 330     stmp = *--sptr;
 331     *--tptr = (stmp&3) +1;
 332     *--tptr = ((stmp >>= 2)&3)+1;
 333     *--tptr = ((stmp >>= 2)&3)+1;
 334     *--tptr = ((stmp >>= 2)&3)+1;
 335   }
 336   /*
 337     for (sptr=seq; sptr < seq+seq_len; sptr++) {
 338     printf("%c",nt[*sptr]);
 339     if ((int)(sptr-seq) % 60 == 59) printf("\n");
 340     }
 341     printf("\n");
 342     */
 343   if (seqcnt*4 >= seq_len) {    /* there was enough room */
 344     seq[seq_len]= EOSEQ;
 345     /* printf("%d\n",seq_len); */
 346     return seq_len;
 347   }
 348   else {                                /* not enough room */
 349     seq[seqcnt*4]=EOSEQ;
 350     seq_len -= 4*seqcnt;
 351     return (4*seqcnt);
 352   }
 353
 354 error:  fprintf(stderr," error reading %ld at %ld\n",libstr,*libpos);
 355   fflush(stderr);
 356   return (-1);
 357 }
 358
 359 void
 360 ncbl_ranlib(str,cnt,libpos)
 361         char *str; int cnt;
 362         long libpos;
 363 {
 364   char hline[256], *bp, *bp0;
 365   int llen;
 366   long spos;
 367
 368   lib_cnt = libpos;
 369   llen = hdr_beg[lib_cnt+1]-hdr_beg[lib_cnt];
 370   if (llen > sizeof(hline)) llen = sizeof(hline);
 371   fseek(hfile,hdr_beg[lib_cnt]+1,0);
 372
 373   fread(hline,(size_t)1,(size_t)(llen-1),hfile);
 374   hline[llen-1]='\0';
 375
 376   if (hline[9]=='|' || hline[10]=='|') {
 377     bp0 = strchr(hline+3,'|');
 378     if ((bp=strchr(bp0+1,' '))!=NULL) *bp='\0';
 379     if (dbformat == NTFORMAT &&
 380         (ambiguity_ray[lib_cnt/CHAR_BIT]&(1<<lib_cnt%CHAR_BIT))) {
 381       sprintf(str,"*%-9s ",bp0+1);
 382     }
 383     else sprintf(str,"%-10s ",bp0+1);
 384     strncat(str+11,bp+1,cnt-strlen(str));
 385   }
 386   else {
 387     if (dbformat == NTFORMAT &&
 388         (ambiguity_ray[lib_cnt/CHAR_BIT]&(1<<lib_cnt%CHAR_BIT))) {
 389       str[0]='*';
 390       strncpy(str+1,hline,cnt-1);
 391     }
 392     else strncpy(str,hline,cnt);
 393   }
 394   str[cnt-1]='\0';
 395
 396   if (dbformat == AAFORMAT)
 397     fseek(sfile,seq_beg[lib_cnt]-1,0);
 398   else {
 399     spos = (seq_beg[lib_cnt])/(CHAR_BIT/NBPN);
 400     fseek(sfile,spos-1,0);
 401   }
 402 }
 403
 404 void src_ulong_read(fd, val)
 405      FILE *fd;
 406      unsigned long *val;
 407 {
 408 #ifdef IS_BIG_ENDIAN
 409   fread((char *)val,(size_t)4,(size_t)1,fd);
 410 #else
 411   unsigned char b[4];
 412
 413   fread((char *)&b[0],(size_t)1,(size_t)4,fd);
 414   *val = 0;
 415   *val = (unsigned long)((unsigned long)((unsigned long)(b[0]<<8) +
 416          (unsigned long)b[1]<<8) + (unsigned long)b[2]<<8)+(unsigned long)b[3];
 417 #endif
 418 }
 419
 420 void src_long_read(fd,val)
 421      FILE *fd;
 422      long *val;
 423 {
 424 #ifdef IS_BIG_ENDIAN
 425   fread((char *)val,(size_t)4,(size_t)1,fd);
 426 #else
 427   unsigned char b[4];
 428
 429   fread((char *)&b[0],(size_t)1,(size_t)4,fd);
 430   *val = 0;
 431   *val = (long)((long)((long)(b[0]<<8)+(long)b[1]<<8)+(long)b[2]<<8)
 432           +(long)b[3];
 433 #endif
 434 }
 435
 436 #ifndef NCBL13_ONLY
 437 static void
 438 #else
 439 void
 440 #endif
 441 src_char_read(fd, val)
 442      FILE *fd;
 443      char *val;
 444 {
 445   fread(val,(size_t)1,(size_t)1,fd);
 446 }
 447
 448 #ifndef NCBL13_ONLY
 449 static void
 450 #else
 451 void
 452 #endif
 453 src_fstr_read(fd, val, slen)
 454      FILE *fd;
 455      char *val;
 456      long slen;
 457 {
 458   fread(val,(size_t)slen,(size_t)1,fd);
 459 }
 460
 461 #ifndef NCBL13_ONLY
 462 static void
 463 #else
 464 void
 465 #endif
 466 newname(char *nname, char *oname, char *suff, int maxn)
 467 {
 468   char *tptr;
 469
 470   if (oname[0]=='@') strncpy(nname,&oname[1],maxn);
 471   else strncpy(nname,oname,maxn);
 472   for (tptr=nname; *tptr=='.' && *tptr; tptr++);
 473   for (; *tptr!='.'&& *tptr; tptr++); /* get to '.' or EOS */
 474   *tptr++='.'; *tptr='\0';
 475   strncat(nname,suff,maxn);
 476 }
 477