mshowbest.c

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

  strncpy(l_name,bline,sizeof(l_name)); /* get rid of text after second "|" */
  l_name[sizeof(l_name)-1]='\0';
  if ((bp=strchr(l_name,' '))!=NULL) *bp=0;
  if ((bp=strchr(&l_name[3],'|'))!=NULL) *bp='\0';
  if (m_msg->nframe > 2) sprintf(&l_name[strlen(l_name)],"_%d",bbp->frame+1);
  else if (m_msg->nframe > 0 && bbp->frame == 1)
    strncat(l_name,"_r",sizeof(l_name));
  if (bbp->seq->cont-1 > 0) {
    sprintf(tmp_str,":%d",bbp->seq->cont-1);
    strncat(l_name,tmp_str,sizeof(l_name)-strlen(l_name));
  }


#ifndef PCOMPLIB
  aln_p = &(m_msg->aln);
  if (m_msg->stages>1 || m_msg->markx & MX_M9SUMM) {
    if (bbp->seq->aa1b == NULL) {

      /* get the sequence but don't save it */
      n1 = re_getlib(aa1save,maxn,m_msg->maxt3,
		     m_msg->l_overlap,bbp->seq->cont,m_msg->term_code,
		     &loffset,&l_off,bbp->seq->m_file_p);
      aa1 = aa1save;
    }
    else {
      n1 = bbp->seq->n1;
      aa1 = bbp->seq->aa1b;
      loffset = bbp->seq->l_offset;
      l_off = bbp->seq->l_off;
    }

    if (! m_msg->markx & MX_M9SUMM) {
      do_opt (aa0[bbp->frame], m_msg->n0, aa1, n1, bbp->frame, ppst, f_str[bbp->frame], &rst);
      bbp->rst.score[2]=rst.score[2];
    }
    else {
      do_walign(aa0[bbp->frame],m_msg->n0, aa1, n1, bbp->frame, 
		ppst, f_str[bbp->frame],
		&bbp->a_res, &bbp->have_ares);

      
      /* save the alignment encoding for future use */
      if (bbp->have_ares && ! (bbp->have_ares & 0x2)) {
	if ((tres = calloc(bbp->a_res.nres+1,sizeof(int)))!=NULL) {
	  memcpy(tres,bbp->a_res.res,sizeof(int)*bbp->a_res.nres);
	  bbp->a_res.res = tres;
	  bbp->have_ares |= 0x2;	/* set 0x2 if has local copy */
	}
	else {
	  bbp->have_ares = 0;
	}
      }

      aln_func_vals(bbp->frame, aln_p);

      /* we need this for offset information for calc_code, but it is
	 incomplete so we must do it again */
      cal_coord(m_msg->n0,BBP_INFO(n1),
		m_msg->q_offset + (m_msg->q_off-1) + (m_msg->sq0off-1),
		loffset + (l_off-1) + (m_msg->sq1off-1),
		aln_p);

      seqc_max = bbp->a_res.nres + 4*m_msg->aln.llen+4;
      seq_code = NULL;
      seq_code_len = 0;
      if (m_msg->show_code == SHOW_CODE_ALIGN) {
	seq_code=(char *)calloc(seqc_max,sizeof(char));
	/* if we have an annotation string, allocate space for the encoded annotation */
	if (m_msg->ann_arr[0] != '\0') {
	  /* the annotation encoding can be considerably longer than the alignment encoding */
	  annc_max = 4*seqc_max;
	  ann_code=(char *)calloc(annc_max,sizeof(char));
	}
	else {
	  ann_code = NULL;
	  annc_max = 0;
	}
	if (seq_code != NULL) {
	  lc=calc_code(aa0[bbp->frame],m_msg->n0,
		       aa1,n1, 
		       &m_msg->aln,&bbp->a_res,
		       ppst,seq_code,seqc_max,
		       m_msg->aa0a, m_msg->ann_arr,
		       ann_code, annc_max,
		       f_str[bbp->frame]);
	  seq_code_len = strlen(seq_code);
	  if (ann_code != NULL) ann_code_len = strlen(ann_code);
	}
      }
      else {
	lc=calc_id(aa0[bbp->frame],m_msg->n0,aa1,n1,
		   &m_msg->aln, &bbp->a_res,
		   ppst,f_str[bbp->frame]);
      }
      m_msg->aln.a_len = lc;

      nident = m_msg->aln.nident;
      if (lc > 0) percent = (100.0*(float)nident)/(float)lc;
      else percent = -1.00;

      ngap = m_msg->aln.ngap_q + m_msg->aln.ngap_l;
#ifndef SHOWSIM
      if (lc-ngap > 0) gpercent = (100.0*(float)nident)/(float)(lc-ngap);
      else gpercent = -1.00;
#else
      if (lc-ngap > 0) gpercent = (100.0*(float)m_msg->aln.nsim)/(float)(lc);
      else gpercent = -1.00;
#endif	/* SHOWSIM */

    }
  }
#endif	/* PCOMPLIB */

  n1tot = (BBP_INFO(n1tot_p)) ? *BBP_INFO(n1tot_p) : BBP_INFO(n1);

  bp = bline;
  if ((m_msg->markx & MX_HTML) && !strncmp(bline,"gi|",3)) {
    bp = strchr(bline+4,'|')+1;
    *(bp-1) = 0;
    gi_num = atoi(bline+3);
  }

#ifndef SUPERFAMNUM
  bp[m_msg->aln.llen-r_margin]='\0';
#else
  bp[m_msg->aln.llen-r_margin-5]='\0';
#endif

  if (m_msg->nframe == -1) bp[m_msg->aln.llen-r_margin]='\0';
  else bp[m_msg->aln.llen-(r_margin+4)]='\0';

#ifndef SUPERFAMNUM
  fprintf (fp, fmt,bp,n1tot);
#else
  if (m_msg->nframe == -1) {
    fprintf (fp, fmt,bp,BBP_INFO(sfnum[0]),n1tot);
  }
  else {fprintf (fp, fmt,bp,n1tot);}
#endif

  if (m_msg->nframe > 2) fprintf (fp, " [%d]", bbp->frame+1);
  else if (m_msg->nframe >= 0) fprintf(fp," [%c]",(bbp->frame > 0 ?'r':'f'));

  if (m_msg->srelv == 1) fprintf (fp, " %4d", bbp->rst.score[ppst->score_ix]);
  else {
    if (m_msg->srelv-1 > 0) fprintf (fp, " %4d", bbp->rst.score[0]);
    if (m_msg->srelv-1 > 1 || m_msg->stages>1)
      fprintf (fp, " %4d", bbp->rst.score[1]);
    fprintf (fp, " %4d", bbp->rst.score[ppst->score_ix]);
  }

  if (ppst->zsflag>=0) { 
    if (m_msg->z_bits==1) {
      fprintf (fp, " %.1f %7.2g",zs_to_bit(bbp->zscore,m_msg->n0,BBP_INFO(n1)),bbp->rst.escore);
    }
    else fprintf (fp, " %.1f %7.2g",bbp->zscore,bbp->rst.escore);
  }
  show_aux(fp,bbp);

#ifdef PCOMPLIB
  n1 = BBP_INFO(n1);
  percent = bbp->percent;
  gpercent = bbp->gpercent;
  aln_p = &(bbp->aln_d);
  seq_code = bbp->aln_code;
  seq_code_len = bbp->aln_code_n;
  ann_code = bbp->ann_code;
  ann_code_len = bbp->ann_code_n;
  loffset = bbp->seq->l_offset;
  l_off = 0;
#endif

  if (m_msg->markx & MX_M9SUMM) {
    if (m_msg->show_code != SHOW_CODE_ID) {
      /* we need the coordinates for annotated SHOW_CODE_ALIGN */
      cal_coord(m_msg->n0,BBP_INFO(n1),
		m_msg->q_offset + (m_msg->q_off-1) + (m_msg->sq0off-1),
		loffset + (l_off-1) + (m_msg->sq1off-1),
		aln_p);

      if (m_msg->markx & MX_HTML) fprintf(fp,"<!-- ");
      /*            %_id  %_sim s-w alen an0  ax0  pn0  px0  an1  ax1  pn1  px1 gapq gapl fs  */
      /*                    alignment    min  max            min  max */
      /*                    sequence coordinate    min  max            min  max */
      fprintf(fp,"\t%5.3f %5.3f %4d %4d %4ld %4ld %4ld %4ld %4ld %4ld %4ld %4ld %3d %3d %3d",
	      percent/100.0,gpercent/100.0, 
#ifndef PCOMPLIB
	      bbp->a_res.sw_score,
#else
	      bbp->sw_score,
#endif
	      aln_p->a_len,
	      aln_p->d_start0,aln_p->d_stop0,
	      aln_p->q_offset+1, aln_p->q_offset+m_msg->n0,
	      aln_p->d_start1,aln_p->d_stop1,
	      aln_p->l_offset+1, aln_p->l_offset+BBP_INFO(n1),
	      aln_p->ngap_q,aln_p->ngap_l,aln_p->nfs);
      if (m_msg->show_code == SHOW_CODE_ALIGN
	  && seq_code_len > 0 && seq_code != NULL) {
	fprintf(fp,"\t%s",seq_code);
	if (ann_code_len > 0 && ann_code != NULL) {
	  fprintf(fp,"\t%s",ann_code);
	}
	/*      fprintf(fp," [%2d:%d]",bbp->wrkr,bbp->seqnm); */


	/* if we are doing MX_M10FORM and -m 9c, then we want to keep
	   the alignment code string for the alignment output - otherwise, we
	   can free() it

	   If PCOMPLIB, then it is stored in bbp->ann_code*, otherwise, it's in a_res.
	*/

#ifndef PCOMPLIB
	if (m_msg->markx & MX_M10FORM)  {
	  /* save encoded alignments in a_res */
	  if ((bbp->a_res.aln_code = (char *)calloc(seq_code_len+1,sizeof(char)))!=NULL) {
	    strncpy(bbp->a_res.aln_code,seq_code,seq_code_len+1);
	    bbp->a_res.aln_code[seq_code_len] = '\0';
	    bbp->a_res.aln_code_n = seq_code_len;
	  }
	}
	/* always free the originally allocated encoding */
	free(seq_code);
	seq_code = NULL;
	seq_code_len = 0;
#else
	/* only free it if not to be used */
	if (!(m_msg->markx & MX_M10FORM))  {
	  free(bbp->aln_code);
	  bbp->aln_code_n = 0;
	}
#endif

	/* also clean up ann_code(_n) */

	if (ann_code_len > 0 && ann_code != NULL) {
#ifndef PCOMPLIB
	  if (m_msg->markx & MX_M10FORM)  {
	    /* save encoded annotations in a_res */
	    if ((bbp->a_res.ann_code = (char *)calloc(ann_code_len+1,sizeof(char)))!=NULL) {
	      strncpy(bbp->a_res.ann_code,ann_code,ann_code_len+1);
	      bbp->a_res.ann_code[ann_code_len] = '\0';
	      bbp->a_res.ann_code_n = ann_code_len;
	    }
	    else {
	      bbp->a_res.ann_code = NULL;
	      bbp->a_res.ann_code_n = 0;
	    }
	  }
	  free(ann_code);
	  ann_code = NULL;
	  ann_code_len = 0;
#else
	  if (!(m_msg->markx & MX_M10FORM))  {
	    free(bbp->ann_code);
	    bbp->ann_code_n = 0;
	  }
#endif
	}
      }
      if (m_msg->markx & MX_HTML) fprintf(fp," -->");
    }
    else {
#ifdef SHOWSIM
      fprintf(fp," %5.3f %5.3f %4d", percent/100.0,(float)aln_p->nsim/(float)aln_p->a_len,aln_p->a_len);
#else
      fprintf(fp," %5.3f %4d", percent/100.0,aln_p->a_len);
#endif
    }
  }
  if (m_msg->markx & MX_HTML) fprintf(fp," <A HREF=\"#%s\">align</A>",l_name);
  fprintf (fp, "\n");
  fflush(fp);
  }

  if (quiet==0) {
    printf(" More scores? [0] ");
    fflush(stdout);
    if (fgets(rline,20,stdin)==NULL) exit(0);
    ntmp = 0;
    if (rline[0]!='\n' && rline[0]!=0) sscanf(rline,"%d",&ntmp);
    if (ntmp<=0) ntmp = 0;
    if (ntmp>0) {
      istart = istop;
      nshow += ntmp;
      goto l1;
    }
  }
  else if (quiet == 1)
    if (ib < nbest && (ppst->zsflag>=0 && bbp->rst.escore < m_msg->e_cut)) {
      if (m_msg->mshow_flg && istop >= m_msg->mshow) goto done;
      istart=istop;
      nshow += 10;
      goto l1;
    }

 done:
  m_msg->nshow = nshow;

  if (m_msg->markx & MX_HTML) fprintf(fp,"</pre></tt><p><hr><p>\n");
  if (fp!=stdout) fprintf(fp,"\n");
}

/*
  q[] has one set of sfnums, 0 terminated
  s[] has second
  return first match or 0
*/