nmgetlib.c

序列对齐 Compare a protein sequence to a protein sequence database or a DNA sequence to a DNA sequenc

  }
  */

  sel_acc_ptr->fmt_term = fmt_term;

  /* now allocate some space for the ACC's */

  abuff_siz = new_buff_siz = 640000;

  if ((acc_buff = (char *)calloc(abuff_siz*10, sizeof(char)))==NULL) {
    fprintf(stderr, "Cannot allocate acc buff %d\n",abuff_siz*10);
    free(sel_acc_ptr);
    return NULL;
  }

  /* now iteratively read and reallocate space for buffer until its all read */

  acc_buff_p = acc_buff;
  while ((buff_siz = fread(acc_buff_p, sizeof(char), new_buff_siz, libf))==new_buff_siz) {
    if ((new_buff = realloc(acc_buff, (size_t)(abuff_siz+new_buff_siz)))==NULL) {
      fprintf(stderr, " cannot reallocate for acc_buf[%d]\n",abuff_siz);
      break;
    }
    else {
      acc_buff = new_buff;
      acc_buff_p = acc_buff + abuff_siz;
      abuff_siz += new_buff_siz;
    }
  }
  fclose(libf);

  acc_buff_max = acc_buff_p + buff_siz;

  /* now convert all the ACC lines (with \n) to null-terminated and
     count the number of aacc's */

  acc_cnt = 0;
  acc_buff_p = acc_buff;
  while (acc_buff_p < acc_buff_max && (bp = strchr(acc_buff_p,'\n'))!=NULL) {
    *bp = '\0';
    /*  also remove '\r'); */
    if ((bp1=strchr(acc_buff_p,'\r'))!=NULL) {*bp1 = '\0';}
    acc_cnt++;
    acc_buff_p = bp+1;
  }

  /* allocate the acc_list */
  if ((acc_list=(char **)calloc(acc_cnt+1, sizeof(char **)))==NULL) {
    fprintf(stderr," cannot allocate acc_list[%d]\n",acc_cnt+1);
    free(sel_acc_ptr);
    return NULL;
  }

  /* now load acc_list[] */
  for (i=0, acc_buff_p=acc_buff; i<acc_cnt; i++) {
    acc_list[i] = acc_buff_p;
    acc_buff_p += strlen(acc_buff_p)+1;
  }

  /* finally put everything in the structure to be returned */
  sel_acc_ptr->acc_buff = acc_buff;
  sel_acc_ptr->acc_list = acc_list;
  sel_acc_ptr->curr_entry = 0;
  sel_acc_ptr->max_entry = acc_cnt;
  return (void *)sel_acc_ptr;
}

int sel_acc_libstr(char *libstr, int gi, void *ptr) {
  struct sel_acc_str *sel_acc_ptr;
  char *curr_acc;
  char acc[MAX_SSTR], *acc_p, *bp;

  sel_acc_ptr = (struct sel_acc_str *)ptr;
  
  if (sel_acc_ptr->curr_entry >= sel_acc_ptr->max_entry) return -1;

  if ((bp = strchr(libstr,sel_acc_ptr->fmt_term))!=NULL) {
    *bp = '\0';
  }

  curr_acc = sel_acc_ptr->acc_list[sel_acc_ptr->curr_entry];
  if (libstr[2] == curr_acc[2] && libstr[1] == curr_acc[1] && 
      strncmp(libstr,curr_acc,MAX_UID)==0) {
    sel_acc_ptr->curr_entry++;
    return 1;
  }
  else {
    return 0;
  }
}

void *sel_acc_gi_init(FILE *libf, int *acc_off, char fmt_term) {
  struct sel_acc_str *sel_acc_ptr;
  char acc_line[MAX_STR];
  char *bp;
  int *gi_list;
  int *new_buff;	/* reallocated buffer size */
  int *acc_buff_p;
  int acc_cnt, i;
  int new_buff_siz;
  int abuff_siz;		/* allocated buffer size */
  int buff_siz;		/* fread buff_siz */

  if ((sel_acc_ptr = (struct sel_acc_str *)calloc(1,sizeof(struct sel_acc_str)))==NULL) {
    fprintf(stderr, "Cannot allocate struct sel_acc_str\n");
    return NULL;
  }

  /*
  if (fmt != NULL && *fmt != '\0') {
    sel_acc_ptr->fmt = (char *)calloc(strlen(fmt)+1,sizeof(char));
    strncpy(sel_acc_ptr->fmt, fmt, strlen(fmt)+1);
    sel_acc_ptr->fmt[strlen(fmt)] = '\0';
  }
  else {
    sel_acc_ptr->fmt = NULL;
  }
  */

  sel_acc_ptr->fmt_term = fmt_term;

  /* now allocate some space for the ACC's */

  abuff_siz = new_buff_siz = 64000;

  if ((gi_list = (int *)calloc(abuff_siz, sizeof(int)))==NULL) {
    fprintf(stderr, "Cannot allocate acc buff %d\n",abuff_siz);
    free(sel_acc_ptr);
    return NULL;
  }

  /* now iteratively read and reallocate space for buffer until its all read */

  acc_cnt = 0;
  while (fgets(acc_line, sizeof(acc_line), libf)!=NULL) {
    gi_list[acc_cnt++] = atoi(acc_line);
    if (acc_cnt >= abuff_siz) {
      if ((new_buff = realloc(gi_list,(abuff_siz + new_buff_siz)*sizeof(int)))==NULL) {
	fprintf(stderr,"cannot realloc gi_list[%d]\n",abuff_siz+new_buff_siz);
	break;
      }
      else {
	abuff_siz += new_buff_siz;
	gi_list = new_buff;
      }
    }
  }
  fclose(libf);

  /* finally put everything in the structure to be returned */
  sel_acc_ptr->gi_list = gi_list;
  sel_acc_ptr->curr_entry = 0;
  sel_acc_ptr->max_entry = acc_cnt;
  return (void *)sel_acc_ptr;
}

int sel_acc_gi(char *libstr, int gi, void *ptr) {
  struct sel_acc_str *sel_acc_ptr;
  char *bp;

  sel_acc_ptr = (struct sel_acc_str *)ptr;
  
  if (sel_acc_ptr->curr_entry >= sel_acc_ptr->max_entry) return -1;

  if (gi <= 0) {
    if (libstr) {
      /*
      if (sel_acc_ptr->fmt) {
	sscanf(libstr,sel_acc_ptr->fmt,&gi);
      }
      */
      gi = atoi(libstr);
    }
  }

  if (gi == sel_acc_ptr->gi_list[sel_acc_ptr->curr_entry]) {
    sel_acc_ptr->curr_entry++;
    return 1;
  }
  else {
    return 0;
  }
}

/* this version of the selection algorithm does not require sorted
   lists.  It hashes the initial list, and uses the hash table to
   lookup the library sequences.

*/

#define HASH_TABLE_MULT 8

void *sel_hacc_libstr_init(FILE *libf, int *acc_off, char fmt_term) {
  struct sel_acc_str *sel_acc_ptr;
  char acc_line[MAX_STR];
  char *bp, *bp1;
  char *acc_buff;
  char *acc_buff_max;	/* end of buffer */
  char *new_buff;	/* reallocated buffer size */
  char *acc_buff_p;
  char **acc_list;
  int acc_cnt, i;
  int new_buff_siz;
  int abuff_siz;		/* allocated buffer size */
  int buff_siz;		/* fread buff_siz */
  int hash_mask, hash_max;
  int hash_val, *acc_hash, *acc_hash_link;
  int link_save;

  if ((sel_acc_ptr = (struct sel_acc_str *)calloc(1,sizeof(struct sel_acc_str)))==NULL) {
    fprintf(stderr, "Cannot allocate struct sel_acc_str\n");
    return NULL;
  }

  /*
  if (fmt && *fmt != '\0') {
    sel_acc_ptr->fmt = (char *)calloc(strlen(fmt)+1,sizeof(char));
    strncpy(sel_acc_ptr->fmt, fmt, strlen(fmt)+1);
    sel_acc_ptr->fmt[strlen(fmt)] = '\0';
  }
  */

  sel_acc_ptr->fmt_term = fmt_term;

  /* now allocate some space for the ACC's */

  abuff_siz = new_buff_siz = 640000;

  if ((acc_buff = (char *)calloc(abuff_siz, sizeof(char)))==NULL) {
    fprintf(stderr, "Cannot allocate acc buff %d\n",abuff_siz);
    free(sel_acc_ptr);
    return NULL;
  }

  /* now iteratively read and reallocate space for buffer until its all read */

  acc_buff_p = acc_buff;
  while ((buff_siz = fread(acc_buff_p, sizeof(char), new_buff_siz, libf))==new_buff_siz) {
    if ((new_buff = realloc(acc_buff, (size_t)(abuff_siz+new_buff_siz)))==NULL) {
      fprintf(stderr, " cannot reallocate for acc_buf[%d]\n",abuff_siz);
      break;
    }
    else {
      acc_buff = new_buff;
      acc_buff_p = acc_buff + abuff_siz;
      abuff_siz += new_buff_siz;
    }
  }
  fclose(libf);

  acc_buff_max = acc_buff_p + buff_siz;

  /* now convert all the ACC lines (with \n) to null-terminated and
     count the number of aacc's */

  acc_cnt = 0;
  acc_buff_p = acc_buff;
  while (acc_buff_p < acc_buff_max && (bp = strchr(acc_buff_p,'\n'))!=NULL) {
    *bp = '\0';
    /*  also remove '\r'); */
    if ((bp1=strchr(acc_buff_p,'\r'))!=NULL) {*bp1 = '\0';}
    acc_cnt++;
    acc_buff_p = bp+1;
  }

  /* allocate the acc_list */
  if ((acc_list=(char **)calloc(acc_cnt+1, sizeof(char **)))==NULL) {
    fprintf(stderr," cannot allocate acc_list[%d]\n",acc_cnt+1);
    free(sel_acc_ptr);
    return NULL;
  }

  /* now load acc_list[] */
  for (i=0, acc_buff_p=acc_buff; i<acc_cnt; i++) {
    acc_list[i] = acc_buff_p;
    acc_buff_p += strlen(acc_buff_p)+1;
  }

  /* allocate the hash for the acc_list - we want a table that is
     about 4X acc_cnt and a power of 2 */
  for (hash_max = 4096; hash_max <= HASH_TABLE_MULT * acc_cnt; hash_max *= 2);
  hash_mask = hash_max - 1;

  if ((acc_hash = (int *)calloc(hash_max, sizeof(int)))==NULL) {
    fprintf(stderr, "cannot allocate acc_hash[%d]\n",hash_max*sizeof(int));
    return NULL;
  }

  /* allocate the acc_list link table */
  if ((acc_hash_link = (int *)calloc(acc_cnt+1,sizeof(char *)))==NULL) {
    fprintf(stderr, "cannot allocate acc_hash_link[%d]\n",acc_cnt*sizeof(char *));
    return NULL;
  }

  for (i=0; i<acc_cnt; i++) {
    hash_val = hash_func(acc_list[i]) & hash_mask;
    if ((link_save = acc_hash[hash_val]) != 0) {
      acc_hash_link[i+1]=link_save;
    }
    acc_hash[hash_val] = i+1;
  }

  /* finally put everything in the structure to be returned */
  sel_acc_ptr->acc_buff = acc_buff;
  sel_acc_ptr->acc_list = acc_list;
  sel_acc_ptr->curr_entry = 0;
  sel_acc_ptr->max_entry = acc_cnt;
  /* hash stuff */
  sel_acc_ptr->acc_hash = acc_hash;
  sel_acc_ptr->acc_hash_link = acc_hash_link;
  sel_acc_ptr->hash_mask = hash_mask;

  return (void *)sel_acc_ptr;
}

int sel_hacc_libstr(char *libstr, int gi, void *ptr) {
  int i, i1;
  struct sel_acc_str *sel_acc_ptr;
  char *bp,  **acc_list;
  int hash_val;

  sel_acc_ptr = (struct sel_acc_str *)ptr;
  
  if (sel_acc_ptr->curr_entry >= sel_acc_ptr->max_entry) return -1;

  if ((bp = strchr(libstr,sel_acc_ptr->fmt_term))!=NULL) {
    *bp = '\0';
  }

  hash_val = hash_func(libstr) & sel_acc_ptr->hash_mask;

  if (sel_acc_ptr->acc_hash[hash_val] == 0) return 0;

  acc_list = sel_acc_ptr->acc_list;

  for (i=sel_acc_ptr->acc_hash[hash_val]; i > 0;
       i = sel_acc_ptr->acc_hash_link[i]) {
    i1 = i-1;
    if (libstr[2] == acc_list[i1][2] && libstr[1] == acc_list[i1][1] && 
      strncmp(libstr,acc_list[i1],MAX_UID)==0) {
      return 1;
    }
  }
  return 0;
}

void *sel_hacc_gi_init(FILE *libf, int *acc_off, char fmt_term) {
  struct sel_acc_str *sel_acc_ptr;
  char acc_line[MAX_STR];
  char *bp;
  int *gi_list;
  int *new_buff;	/* reallocated buffer size */
  int *acc_buff_p;
  int acc_cnt, i;
  int new_buff_siz;
  int abuff_siz;		/* allocated buffer size */
  int buff_siz;		/* fread buff_siz */
  int hash_val, hash_max, hash_mask, link_save;
  int *gi_hash, *gi_hash_link;

  if ((sel_acc_ptr = (struct sel_acc_str *)calloc(1,sizeof(struct sel_acc_str)))==NULL) {
    fprintf(stderr, "Cannot allocate struct sel_acc_str\n");
    return NULL;
  }

  /*
  if (fmt != NULL && *fmt != '\0') {
    sel_acc_ptr->fmt = (char *)calloc(strlen(fmt)+1,sizeof(char));
    strncpy(sel_acc_ptr->fmt, fmt, strlen(fmt)+1);
    sel_acc_ptr->fmt[strlen(fmt)] = '\0';
  }
  else {
    sel_acc_ptr->fmt = NULL;
  }
  */

  sel_acc_ptr->fmt_term = fmt_term;

  /* now allocate some space for the ACC's */

  abuff_siz = new_buff_siz = 64000;

  if ((gi_list = (int *)calloc(abuff_siz, sizeof(int)))==NULL) {
    fprintf(stderr, "Cannot allocate acc buff %d\n",abuff_siz);
    free(sel_acc_ptr);
    return NULL;
  }

  /* now iteratively read and reallocate space for buffer until its all read */

  acc_cnt = 0;
  while (fgets(acc_line, sizeof(acc_line), libf)!=NULL) {
    gi_list[acc_cnt++] = atoi(acc_line);
    if (acc_cnt >= abuff_siz) {
      if ((new_buff = realloc(gi_list,(abuff_siz + new_buff_siz)*sizeof(int)))==NULL) {
	fprintf(stderr,"cannot realloc gi_list[%d]\n",abuff_siz+new_buff_siz);
	break;
      }
      else {
	abuff_siz += new_buff_siz;
	gi_list = new_buff;
      }
    }
  }
  fclose(libf);

  /* allocate the hash for the gi_list - we want a table that is
     about 4X acc_cnt and a power of 2 */
  for (hash_max = 4096; hash_max <= HASH_TABLE_MULT * acc_cnt; hash_max *= 2);
  hash_mask = hash_max - 1;

  if ((gi_hash = (int *)calloc(hash_max, sizeof(int)))==NULL) {
    fprintf(stderr, "cannot allocate gi_hash[%d]\n",hash_max*sizeof(int));
    return NULL;
  }

  /* allocate the gi_list link table */
  if ((gi_hash_link = (int *)calloc(acc_cnt+1,sizeof(char *)))==NULL) {
    fprintf(stderr, "cannot allocate gi_hash_link[%d]\n",acc_cnt*sizeof(char *));
    return NULL;
  }

  for (i=0; i<acc_cnt; i++) {
    hash_val = fast_hash32(gi_list[i]) & hash_mask;
    if (gi_hash[hash_val] != 0) {
      link_save = gi_hash[hash_val];
      gi_hash_link[i+1]=link_save;
    }
    gi_hash[hash_val] = i+1;
  }

  /* finally put everything in the structure to be returned */
  sel_acc_ptr->gi_list = gi_list;
  sel_acc_ptr->curr_entry = 0;
  sel_acc_ptr->max_entry = acc_cnt;

  sel_acc_ptr->acc_hash = gi_hash;
  sel_acc_ptr->acc_hash_link = gi_hash_link;
  sel_acc_ptr->hash_mask = hash_mask;

  return (void *)sel_acc_ptr;
}

int sel_hacc_gi(char *libstr, int gi, void *ptr) {
  struct sel_acc_str *sel_acc_ptr;
  int hash_val;
  int *gi_list, i;
  char *bp;

  sel_acc_ptr = (struct sel_acc_str *)ptr;
  
  if (sel_acc_ptr->curr_entry >= sel_acc_ptr->max_entry) return -1;

  if (gi <= 0) {
    if (libstr) {
      /*
      if (sel_acc_ptr->fmt) {
	sscanf(libstr,sel_acc_ptr->fmt,&gi);
      }
      */
      gi = atoi(libstr);
    }
  }

  hash_val = fast_hash32(gi) & sel_acc_ptr->hash_mask;
  if (sel_acc_ptr->acc_hash[hash_val]==0) return 0;

  gi_list = sel_acc_ptr->gi_list;

  for (i=sel_acc_ptr->acc_hash[hash_val]; i > 0;
       i = sel_acc_ptr->acc_hash_link[i]) {
    if (gi_list[i-1] == gi) return 1;
  }
  return 0;
}


/* adapted from	http://burtleburtle.net/bob/hash/doobs.html */
unsigned int
hash_func(char *key)
{
  unsigned int hash;

  hash = 0;

  while (*key) {
    hash += *key++;
    hash += (hash << 10);
    hash ^= (hash >> 6);
  }

  hash += (hash << 3);
  hash ^= (hash >> 11);
  hash += (hash << 15);

  return hash;
}

unsigned int
fast_hash32 (unsigned int data) {
  int tmp, hash;


  hash  = data >> 16;
  tmp    = ((data & 0xFFFF) << 11) ^ hash;
  hash   = (hash << 16) ^ tmp;
  hash  += hash >> 11;

  /* Force "avalanching" of final 127 bits */
  hash ^= hash << 3;
  hash += hash >> 5;
  hash ^= hash << 4;
  hash += hash >> 17;
  hash ^= hash << 25;
  hash += hash >> 6;

  return hash;
}