compacc.c

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

    }
  }

  if (ppst->zsflag_f>=0) {
    if (ntt.carry==0) {
      fprintf(fd, "%7ld residues in %5ld sequences", ntt.length, ntt.entries);
    }
    else {
      db_tt = (double)ntt.carry*(double)LONG_MAX + (double)ntt.length;
      fprintf(fd, "%.0f residues in %5ld library sequences", db_tt, ntt.entries);
    }
    fprintf(fd, "%s\n",lib_range);
    fprintf(fd,"Statistics: %s\n",hist.stat_info);

    if (MAX_STATS < hist.entries)
#ifdef SAMP_STATS
      fprintf(fd," statistics sampled from %d to %ld sequences\n",
	      MAX_STATS,hist.entries);
#else
      fprintf(fd," statistics extrapolated from %d to %ld sequences\n",
	      MAX_STATS,hist.entries);
#endif
    /*    summ_stats(stat_info); */
      if (!m_msg.nohist && cum_hl > 0) {
	fprintf(fd," Kolmogorov-Smirnov  statistic: %6.4f (N=%d) at %3d\n",
	      max_dev/(float)cum_hl, n_chi_sq,max_i*hist.histint+hist.min_hist);
      }
      if (m_msg.markx & MX_M10FORM) {
	while ((bp=strchr(hist.stat_info,'\n'))!=NULL) *bp=' ';
	if (cum_hl <= 0) cum_hl = -1;
	sprintf(info_hstring[0],"; mp_extrap: %d %ld\n; mp_stats: %s\n; mp_KS: %6.4f (N=%d) at %3d\n",
		MAX_STATS,hist.entries,hist.stat_info,max_dev/(float)cum_hl, 
		n_chi_sq,max_i*hist.histint+hist.min_hist);
      }
  }
  if (m_msg.markx & MX_M10FORM) {
    if ((bp = strchr(info_gstring2[1],'\n'))!=NULL) *bp = ' ';
    sprintf(info_hstring[1],"; mp_Algorithm: %s\n; mp_Parameters: %s\n",info_gstring2[0],info_gstring2[1]);
    if (bp != NULL ) *bp = '\n';
  }

  if (ppst->other_info != NULL) {
    fputs(ppst->other_info, fd);
  }

  fprintf(fd,"Algorithm: %s\nParameters: %s\n",info_gstring2[0],info_gstring2[1]);
  

  fflush(fd);
}

extern char prog_name[], *verstr;

void s_abort (char *p,  char *p1)
{
  int i;

  fprintf (stderr, "\n***[%s] %s%s***\n", prog_name, p, p1);
#ifdef PCOMPLIB
#ifdef PVM_SRC
  for (i=FIRSTNODE; i< nnodes; i++) pvm_kill(pinums[i]);
  pvm_exit();
#endif
#ifdef MPI_SRC
  MPI_Abort(MPI_COMM_WORLD,1);
  MPI_Finalize();
#endif
#endif
  exit (1);
}

#ifndef MPI_SRC
void w_abort (char *p, char *p1)
{
  fprintf (stderr, "\n***[%s] %s%s***\n\n", prog_name, p, p1);
  exit (1);
}
#endif

#ifndef PCOMPLIB
/* copies from from to to shuffling */

extern int nrand(int);

void
shuffle(unsigned char *from, unsigned char *to, int n)
{
  int i,j; unsigned char tmp;

  if (from != to) memcpy((void *)to,(void *)from,n);

  for (i=n; i>0; i--) {
    j = nrand(i);
    tmp = to[j];
    to[j] = to[i-1];
    to[i-1] = tmp;
  }
  to[n] = 0;
}

static int ieven = 0;
/* copies from from to from shuffling, ieven changed for threads */
void
wshuffle(unsigned char *from, unsigned char *to, int n, int wsiz)
{
  int i,j, k, mm; 
  unsigned char tmp, *top;

  memcpy((void *)to,(void *)from,n);
	
  mm = n%wsiz;

  if (ieven) {
    for (k=0; k<(n-wsiz); k += wsiz) {
      top = &to[k];
      for (i=wsiz; i>0; i--) {
	j = nrand(i);
	tmp = top[j];
	top[j] = top[i-1];
	top[i-1] = tmp;
      }
    }
    top = &to[n-mm];
    for (i=mm; i>0; i--) {
      j = nrand(i);
      tmp = top[j];
      top[j] = top[i-1];
      top[i-1] = tmp;
    }
    ieven = 0;
  }
  else {
    for (k=n; k>=wsiz; k -= wsiz) {
      top = &to[k-wsiz];
      for (i=wsiz; i>0; i--) {
	j = nrand(i);
	tmp = top[j];
	top[j] = top[i-1];
	top[i-1] = tmp;
      }
    }
    top = &to[0];
    for (i=mm; i>0; i--) {
      j = nrand(i);
      tmp = top[j];
      top[j] = top[i-1];
      top[i-1] = tmp;
    }
    ieven = 1;
  }
  to[n] = 0;
}

#endif

int
sfn_cmp(int *q, int *s)
{
  if (*q == *s) return *q;
  while (*q && *s) {
    if (*q == *s) return *q;
    else if (*q < *s) q++;
    else if (*q > *s) s++;
  }
  return 0;
}

#ifndef MPI_SRC

#define ESS 49

void
revcomp(unsigned char *seq, int n, int *c_nt)
{
  unsigned char tmp;
  int i, ni;

  for (i=0, ni = n-1; i< n/2; i++,ni--) {
    tmp = c_nt[seq[i]];
    seq[i] = c_nt[seq[ni]];
    seq[ni] = tmp;
  }
  if ((n%2)==1) {
    i = n/2;
    seq[i] = c_nt[seq[i]];
  }
  seq[n]=0;
}
#endif

#ifdef PCOMPLIB

/* init_stage2 sets up the data structures necessary to send a subset
   of sequences to the nodes, and then collects the results
*/

/* wstage2[] FIRSTNODE .. nnodes has the next sequence to be do_opt()/do_walign()ed */
/* wstage2p[] is a list of sequence numbers/frames, to be sent to workers */
/* wstage2b[] is a list of bptr's that shares the index with wstage2p[] */

static int  wstage2[MAXWRKR +1];	/* count of second stage scores */
static struct stage2_str  *wstage2p[MAXWRKR+1]; /* list of second stage sequences */
static int  wstage2i[MAXWRKR+1];	/* index into second stage sequences */
static struct beststr *bbptr,
     **wstage2b[MAXWRKR+1];	/* reverse pointers to bestr */

void
do_stage2(struct beststr **bptr, int nbest, struct mngmsg m_msg0,
	  int s_func, struct qmng_str *qm_msp) {

  int i, is, ib, iw, nres;
  int node, snode, node_done;
  int bufid, numt, tid;
  char errstr[120];
  struct comstr2 bestr2[BFR2+1];	/* temporary structure array */
  char *seqc_buff, *seqc;
  int seqc_buff_len, aln_code_n;
#ifdef MPI_SRC
  MPI_Status mpi_status;
#endif

  /* initialize the counter for each worker to 0 */
  for (iw = FIRSTNODE; iw < nnodes; iw++) wstage2[iw] = 0;

  /* for each result, bump the counter for the worker that has
     the sequence */
  for (ib = 0; ib < nbest; ib++ ) { wstage2[bptr[ib]->seq->wrkr]++;  }

  /* now allocate enough space to send each worker a 
     list of its sequences stage2_str {seqnm, frame} */
  for (iw = FIRSTNODE; iw < nnodes; iw++) {
    if (wstage2[iw]>0) {
      if ((wstage2p[iw]=
	   (struct stage2_str *)
	   calloc(wstage2[iw],sizeof(struct stage2_str)))==NULL) {
	sprintf(errstr," cannot allocate sequence listp %d %d",
		iw,wstage2[iw]);
	s_abort(errstr,"");
      }

      /* allocate space to remember the bptr's for each result */
      if ((wstage2b[iw]=(struct beststr **)
	   calloc(wstage2[iw],sizeof(struct beststr *)))==NULL) {
	sprintf(errstr," cannot allocate sequence listb %d %d",
		iw,wstage2[iw]);
	s_abort(errstr,"");
      }
      wstage2i[iw]=0;
    }
    else {
      wstage2p[iw] = NULL;
      wstage2b[iw] = NULL;
    }
  }

  /* for each result, set wstage2p[worker][result_index_in_worker] */
  for (is = 0; is < nbest; is++) {
    iw=bptr[is]->seq->wrkr;
    wstage2p[iw][wstage2i[iw]].seqnm = bptr[is]->seq->seqnm;
    wstage2p[iw][wstage2i[iw]].frame = bptr[is]->frame;
    wstage2b[iw][wstage2i[iw]] = bptr[is];
    wstage2i[iw]++;
  }


  /* at this point, wstage2i[iw] should equal wstage2[iw] */
  node_done = 0;
  for (node = FIRSTNODE; node < nnodes; node++) {

    /*    fprintf(stderr,"node: %d stage2: %d\n",node,wstage2[node]); */

    /* if a worker has no results, move on */
    if (wstage2[node]<=0) { node_done++; continue;}

    qm_msp->slist = wstage2[node];	/* set number of results to return */
    qm_msp->s_func = s_func;		/* set s_funct for do_opt/do_walign */
#ifdef PVM_SRC
    pvm_initsend(PvmDataRaw);
    pvm_pkbyte((char *)qm_msp,sizeof(struct qmng_str),1);
    pvm_send(pinums[node],MSEQTYPE);	/* send qm_msp */
    pvm_initsend(PvmDataRaw);	/* send the list of seqnm/frame */
    pvm_pkbyte((char *)wstage2p[node],wstage2[node]*sizeof(struct stage2_str),1);
    pvm_send(pinums[node],LISTTYPE);
#endif
#ifdef MPI_SRC
    MPI_Send(qm_msp,sizeof(struct qmng_str),MPI_BYTE,node,MSEQTYPE,
	     MPI_COMM_WORLD);
    MPI_Send((char *)wstage2p[node],wstage2[node]*
	     sizeof(struct stage2_str),MPI_BYTE,node,LISTTYPE,
	     MPI_COMM_WORLD);
#endif
  }
	
  /* all the workers have their list of sequences */
  /* reset the index of results to obtain */
  for (iw = 0; iw < nnodes; iw++) wstage2i[iw]=0;
	
  while (node_done < nnodes-FIRSTNODE) {
#ifdef PVM_SRC
    bufid = pvm_recv(-1,LISTRTYPE);	/* wait for results */
    pvm_bufinfo(bufid,NULL,NULL,&tid);
    /* get a chunk of comstr2 results */
    pvm_upkbyte((char *)&bestr2[0],sizeof(struct comstr2)*(BFR2+1),1);
    snode = (iw=tidtonode(tid));
    pvm_freebuf(bufid);
#endif
#ifdef MPI_SRC
    MPI_Recv((char *)&bestr2[0],sizeof(struct comstr2)*(BFR2+1),
	     MPI_BYTE,MPI_ANY_SOURCE,LISTRTYPE,MPI_COMM_WORLD,
	     &mpi_status);
    snode = mpi_status.MPI_SOURCE;
    iw = snode;
#endif

    seqc_buff = NULL;
    if (s_func == DO_OPT_FLG && m_msg0.show_code==SHOW_CODE_ALIGN) {
#ifdef PVM_SRC
      bufid = pvm_recv(tid,CODERTYPE);
      pvm_upkint(&seqc_buff_len,1,1);	/* get the code string length */
#endif
#ifdef MPI_SRC
      MPI_Recv((char *)&seqc_buff_len,1,MPI_INT, snode,
	       CODERTYPE,MPI_COMM_WORLD, &mpi_status);
#endif

      seqc=seqc_buff = NULL;
      if (seqc_buff_len > 0) {		/* allocate space for it */
	if ((seqc=seqc_buff=calloc(seqc_buff_len,sizeof(char)))==NULL) {
	  fprintf(stderr,"Cannot allocate seqc_buff: %d\n",seqc_buff_len);
	  seqc_buff_len=0;
	}
	else {
#ifdef PVM_SRC
	  pvm_upkbyte(seqc_buff,seqc_buff_len*sizeof(char),1);
#endif
#ifdef MPI_SRC
	  MPI_Recv((char *)seqc_buff,seqc_buff_len*sizeof(char),
		   MPI_BYTE,snode,CODERTYPE,MPI_COMM_WORLD, &mpi_status);
#endif
	}
      }
#ifdef PVM_SRC
      pvm_freebuf(bufid);
#endif
    }

    /* get number of results in this message */
    nres = bestr2[BFR2].seqnm & ~FINISHED;
    /* check to see if finished */
    if (bestr2[BFR2].seqnm&FINISHED) {node_done++;}

    seqc = seqc_buff;

    /* count through results from a specific worker */
    for (i=0,is=wstage2i[iw]; i < nres; i++,is++) {

      /* get the (saved) bptr for this result */
      bbptr=wstage2b[iw][is];
      /* consistency check seqnm's must agree */
      if (wstage2p[iw][is].seqnm ==  bbptr->seq->seqnm) {
	if (s_func == DO_CALC_FLG && m_msg0.last_calc_flg) {
	  bbptr->rst.score[0] = bestr2[i].rst.score[0];
	  bbptr->rst.score[1] = bestr2[i].rst.score[1];
	  bbptr->rst.score[2] = bestr2[i].rst.score[2];
	  bbptr->rst.escore = bestr2[i].rst.escore;
	  bbptr->rst.segnum = bestr2[i].rst.segnum;
	  bbptr->rst.seglen = bestr2[i].rst.seglen;
	}
	else if (m_msg0.stages > 1) {
	  bbptr->rst.score[0] = bestr2[i].rst.score[0];
	  bbptr->rst.score[1] = bestr2[i].rst.score[1];
	  bbptr->rst.score[2] = bestr2[i].rst.score[2];
	}

	if (s_func == DO_OPT_FLG && m_msg0.markx & MX_M9SUMM) {
	  /* get score, alignment information, percents */
#ifndef PCOMPLIB
	  bbptr->a_res.sw_score = bestr2[i].sw_score;
#else
	  bbptr->sw_score = bestr2[i].sw_score;
#endif
	  memcpy(&(bbptr->aln_d),&bestr2[i].aln_d,sizeof(struct a_struct));
	  bbptr->percent = bestr2[i].percent;
	  bbptr->gpercent = bestr2[i].gpercent;

	  if (m_msg0.show_code == 2) {	/* if show code */
	    /* length of encoding */
	    aln_code_n = bbptr->aln_code_n = bestr2[i].aln_code_n;
	    bbptr->aln_code = NULL;
	    /* right now we don't return ann encodings */
	    bbptr->ann_code = NULL;
	    bbptr->ann_code_n = 0;

	    if (aln_code_n > 0) {
	      if ((bbptr->aln_code = 
		   (char *)calloc(aln_code_n+1,sizeof(char)))==NULL) {
		fprintf(stderr,"cannot allocate seq_code[%d:%d]: %d\n",
			bbptr->seq->wrkr,bbptr->seq->seqnm,aln_code_n);
		seqc += aln_code_n+1;
		bbptr->aln_code_n = 0;
	      }
	      else {
		strncpy(bbptr->aln_code,seqc,aln_code_n);
		bbptr->aln_code[aln_code_n]='\0';
		seqc += aln_code_n+1;
	      }
	    }
	    else {
	      fprintf(stderr," aln_code_n <=0: %d\n",aln_code_n);
	    }
	  }
	}
      }
      else fprintf(stderr,"phase error in phase II return %d %d", iw,i);
    }
    if (seqc_buff != NULL) {
      free(seqc_buff);
      seqc_buff = NULL;
    }
    wstage2i[iw] += nres;
  }

  for (iw=FIRSTNODE; iw < nnodes; iw++) {
    if ((void *)wstage2p[iw]!=NULL) free((void *)wstage2p[iw]);
    if ((void *)wstage2b[iw]!=NULL) free((void *)wstage2b[iw]);
  }
}

#endif