mshowalign.c

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

    initseq(&seqc0, &seqc0a, &seqc1, &seqca, maxc);

    pvm_recv(pinums[bbp->seq->wrkr],ALN2TYPE);
    pvm_upkbyte(seqc0,maxc,1);
    if (m_msg.q_ann_flg) pvm_upkbyte(seqc0a,maxc,1);
    pvm_upkbyte(seqc1,maxc,1);
    pvm_upkbyte(seqca,maxc,1);
#endif	/* PVM_SRC */
#ifdef MPI_SRC
    MPI_Recv(int_msg_b,4,MPI_INT,bbp->seq->wrkr,ALN1TYPE,MPI_COMM_WORLD,
	     &mpi_status);
    nc = int_msg_b[0];
    lc = int_msg_b[1];
    maxc = int_msg_b[2];
    bbp->sw_score = int_msg_b[3];
    MPI_Recv(&percent,1,MPI_FLOAT,bbp->seq->wrkr,ALN2TYPE,MPI_COMM_WORLD,
	     &mpi_status);
    MPI_Recv(&gpercent,1,MPI_FLOAT,bbp->seq->wrkr,ALN2TYPE,MPI_COMM_WORLD,
	     &mpi_status);
    MPI_Recv(aln_p,sizeof(struct a_struct),MPI_BYTE,
	     bbp->seq->wrkr,ALN3TYPE,MPI_COMM_WORLD,&mpi_status);

    initseq(&seqc0, &seqc0a, &seqc1, &seqca, maxc);
    MPI_Recv(seqc0,maxc,MPI_BYTE,bbp->seq->wrkr,ALN2TYPE,MPI_COMM_WORLD,&mpi_status);
    if (m_msg.q_ann_flg)
      MPI_Recv(seqc0a,maxc,MPI_BYTE,bbp->seq->wrkr,ALN2TYPE,MPI_COMM_WORLD,&mpi_status);
    MPI_Recv(seqc1,maxc,MPI_BYTE,bbp->seq->wrkr,ALN3TYPE,MPI_COMM_WORLD,&mpi_status);
    MPI_Recv(seqca,maxc,MPI_BYTE,bbp->seq->wrkr,ALN3TYPE,MPI_COMM_WORLD,&mpi_status);
#endif	/* MPI_SRC */

    loffset = bbp->seq->l_offset-l_off;
    lsw_score = bbp->sw_score;
#else	/* not PCOMPLIB */

    while (cur_ares != NULL && cur_ares->nres > 0) {

      /* estimate space for alignment consensus */
      if (m_msg.aln.showall==1) {
	maxc = cur_ares->nres + max(cur_ares->min0,cur_ares->min1)+
	  max((m_msg.n0-cur_ares->max0),(n1-cur_ares->max1))+4;
      }
      else {
	maxc = cur_ares->nres + 4*m_msg.aln.llen+4;
      }

      /* get space to put the sequence alignment consensus */
      initseq(&seqc0, &seqc0a, &seqc1, &seqca, maxc);

#ifdef LALIGN
      if ((m_msg.markx & MX_M11OUT) == MX_M11OUT) {	/* lav output - skip lots o stuff */
	lsw_score = cur_ares->sw_score;
	zscore = (*find_zp)(lsw_score, 0.0, BBP_INFO(n1), 0.0, m_msg.pstat_void);
	bits = zs_to_bit(zscore, m_msg.n0, BBP_INFO(n1));


	lc=calc_code(aa0[bbp->frame],m_msg.n0,
		     aa1,n1, 
		     aln_p, cur_ares,
		     ppst,seqc0,maxc, m_msg.aa0a, m_msg.ann_arr,
		     seqc0a, maxc,
		     f_str[bbp->frame]);
	seqc_len = strlen(seqc0);
	seqca_len = strlen(seqc0a);

	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 (lc > 0) percent = (100.0*(float)aln_p->nident)/(float)lc;
	else percent = -1.00;

	fprintf (fp, "a {\n  s %d %.1f\n", lsw_score, bits);
	do_lav(fp, aln_p, seqc0, percent, 0);

	if (ppst->nseq == 1) {
	  fprintf (fp, "a {\n  s %d %.1f\n", lsw_score, bits);
	  do_lav(fp, aln_p, seqc0, percent, 1);
	}

	cur_ares = cur_ares->next;
	continue;
      }
#endif		/* ifdef LALIGN */

      /* build consensus from res, nres (done by workers if PCOMPLIB) */
      nc=calcons_a(aa0[bbp->frame],m_msg.n0,aa1,n1,
		   &lc,aln_p, cur_ares, ppst, 
		   seqc0, seqc1, seqca, 
		   m_msg.aa0a, seqc0a, m_msg.ann_arr,
		   f_str[bbp->frame]);

      /* PCOMPLIB workers return percent, gpercent, so calculate it here */
      if (lc > 0) percent = (100.0*(float)aln_p->nident)/(float)lc;
      else percent = -1.00;
      ngap = aln_p->ngap_q + aln_p->ngap_l;
#ifndef SHOWSIM
      if (lc-ngap> 0) gpercent =(100.0*(float)aln_p->nident)/(float)(lc-ngap);
#else
      if (lc > 0) gpercent =(100.0*(float)aln_p->nsim)/(float)lc;
#endif
      else gpercent = -1.00;

      lsw_score = cur_ares->sw_score;
#endif	/* not PCOMPLIB */

      if (max(strlen(seqc0),strlen(seqc1)) > nc) {
	fprintf(stderr," mshowalign: nc/maxc: %d/%d seqc0/1: %lu/%lu\n",
		nc,maxc,strlen(seqc0),strlen(seqc1));
      }

      /* here PCOMPLIB/comp_lib logic is the same */

#ifdef DEBUG
      /*
	if (lsw_score < bbp->rst.score[ppst->score_ix]) {
	fprintf(stderr," *** warning - SW score=%d < opt score=%d ***\n",
	lsw_score, bbp->rst.score[ppst->score_ix]);
	}
      */
#endif
      cal_coord(m_msg.n0,BBP_INFO(n1),
		m_msg.q_offset+(m_msg.q_off-1)+(m_msg.sq0off-1),
		bbp->seq->l_offset+(bbp->seq->l_off-1)+(m_msg.sq1off-1),
		aln_p);

      zscore = (*find_zp)(lsw_score, 0.0, BBP_INFO(n1), 0.0, m_msg.pstat_void);
      /*      bits = s_to_bit(lsw_score, m_msg.n0, BBP_INFO(n1), m_msg.pstat_void); */
      bits = zs_to_bit(zscore, m_msg.n0, BBP_INFO(n1));

      if ((m_msg.markx & MX_ATYPE)!=7 && !(m_msg.markx & MX_M10FORM)) {

	/* this code makes sense for library searches, but not for
	   multiple non-intersecting alignments */

#ifndef LALIGN
	if (m_msg.nframe > 2) 
	  fprintf (fp, "Frame: %d",bbp->frame+1);
	else if (m_msg.nframe > 1) 
	  fprintf (fp, "Frame: %c",(bbp->frame? 'r': 'f'));
	else if (m_msg.qframe >= 0 && bbp->frame > 0 ) {
	  fputs("rev-comp",fp);
	  name0[nml-1]='\0';
	  strcat(name0,"-");
	}

	if (m_msg.arelv > 0)
	  fprintf (fp, " %s: %3d", m_msg.alab[0],bbp->rst.score[0]);
	if (m_msg.arelv > 1)
	  fprintf (fp, " %s: %3d", m_msg.alab[1],bbp->rst.score[1]);
	if (m_msg.arelv > 2)
	  fprintf (fp, " %s: %3d", m_msg.alab[2],bbp->rst.score[2]);
	fprintf (fp,"%s",info_str);
	if (ppst->zsflag>=0) 
	  fprintf (fp, "  Z-score: %4.1f  bits: %3.1f E(): %4.2g", 
		   bbp->zscore,zs_to_bit(bbp->zscore,m_msg.n0,BBP_INFO(n1)),
		   bbp->rst.escore);
	fprintf (fp, "\n");
#else 
	if ((m_msg.markx & MX_M11OUT) == 0) {
	  fprintf (fp, " %s score: %d; ", m_msg.alabel, lsw_score);
	  fprintf (fp," %3.1f bits; E(%ld) <  %.2g\n", bits, ppst->zdb_size,
		   zs_to_E(zscore, BBP_INFO(n1), ppst->dnaseq, ppst->zdb_size, m_msg.db));
	}
#endif
      }
      else if (m_msg.markx & MX_M10FORM) {
#ifndef LALIGN
	if (m_msg.qframe > -1) {
	  if (m_msg.nframe > 2) {
	    fprintf(fp,"; %s_frame: %d\n",m_msg.f_id0,bbp->frame+1);
	  }
	  else {
	    fprintf(fp,"; %s_frame: %c\n",m_msg.f_id0,(bbp->frame > 0? 'r':'f'));
	  }
	}
	fprintf (fp, "; %s_%s: %3d\n", m_msg.f_id0,m_msg.alab[0],bbp->rst.score[0]);
	if (m_msg.arelv > 1)
	  fprintf (fp,"; %s_%s: %3d\n", m_msg.f_id0,m_msg.alab[1],bbp->rst.score[1]);
	if (m_msg.arelv > 2)
	  fprintf (fp,"; %s_%s: %3d\n", m_msg.f_id0,m_msg.alab[2],bbp->rst.score[2]);
	if (info_str[0]) fprintf(fp,"; %s_info: %s\n",m_msg.f_id0,info_str);
	if (ppst->zsflag>=0) 
	  fprintf (fp,"; %s_z-score: %4.1f\n; %s_bits: %3.1f\n; %s_expect: %6.2g\n",
		   m_msg.f_id0,bbp->zscore,
		   m_msg.f_id0,zs_to_bit(bbp->zscore, m_msg.n0, BBP_INFO(n1)),
		   m_msg.f_id0,bbp->rst.escore);
#else
	if ((m_msg.markx & MX_M11OUT) == 0) {
	  fprintf (fp,"; %s_%s score: %d\n", m_msg.f_id0, m_msg.alab[0], lsw_score);
	  fprintf (fp,"; %s_z-score: %4.1f\n; %s_bits %3.1f\n; %s_expect: %6.2g\n",
		   m_msg.f_id0, zscore, m_msg.f_id0, bits, m_msg.f_id0, 
		   zs_to_E(zscore, BBP_INFO(n1), ppst->dnaseq, ppst->zdb_size, m_msg.db));
	}
#endif	
      }

      do_show(fp, m_msg.n0, BBP_INFO(n1), lsw_score, name0, name1, nml,
	      m_msg, ppst, seqc0, seqc0a, seqc1, seqca,
	      nc, percent, gpercent, lc, aln_p);

      /* display the encoded alignment left over from showbest()*/

      if ((m_msg.markx & MX_M10FORM) &&
	  (m_msg.markx & MX_M9SUMM) && 
	  (m_msg.show_code == SHOW_CODE_ALIGN)) {

#ifdef PCOMPLIB
	seq_code = bbp->aln_code;
	seq_code_len = bbp->aln_code_n;
	ann_code = bbp->ann_code;
	ann_code_len = bbp->ann_code_n;
#else
	seq_code = cur_ares->aln_code;
	seq_code_len = cur_ares->aln_code_n;
	ann_code = cur_ares->ann_code;
	ann_code_len = cur_ares->ann_code_n;
#endif

	if (seq_code_len > 0 && seq_code != NULL) {
	  fprintf(fp,"; al_code: %s\n",seq_code);
	  free(seq_code);
	  if (ann_code_len > 0 && ann_code != NULL) {
	    fprintf(fp,"; al_code_ann: %s\n",ann_code);
	    free(ann_code);
	  }
	}
      }

      if (m_msg.markx & MX_HTML) fprintf(fp,"</pre></tt>\n<hr>\n");
      fflush(fp);

      freeseq(&seqc0,&seqc0a,&seqc1, &seqca);

#ifndef PCOMPLIB
    cur_ares = cur_ares->next;

    if (cur_ares != NULL) {
      if (m_msg.markx & MX_HTML) {
	fprintf (fp,"<tt><pre>\n");
      }
      else {
	fprintf(fp,">--\n");
      }
    }
    }
#endif
  }
  if (fp!=stdout) fprintf(fp,"\n");
}

void do_show(FILE *fp, int n0,int n1, int score,
	     char *name0, char *name1, int nml,
	     struct mngmsg m_msg, struct pstruct *ppst,
	     char *seqc0, char *seqc0a,  char *seqc1, char *seqca, int nc,
	     float percent, float gpercent, int lc,
	     struct a_struct *aln)
{
  int tmp;

  if (m_msg.markx & MX_AMAP && (m_msg.markx & MX_ATYPE)==7)
    disgraph(fp, n0, n1, percent, score,
	     aln->amin0, aln->amin1, aln->amax0, aln->amax1, m_msg.sq0off,
	     name0, name1, nml, aln->llen, m_msg.markx);
  else if (m_msg.markx & MX_M10FORM) {
    if (ppst->sw_flag && m_msg.arelv>0)
      fprintf(fp,"; %s_score: %d\n",m_msg.f_id1,score);
    fprintf(fp,"; %s_ident: %5.3f\n",m_msg.f_id1,percent/100.0);
#ifndef SHOWSIM
    fprintf(fp,"; %s_gident: %5.3f\n",m_msg.f_id1,gpercent/100.0);
#else
    fprintf(fp,"; %s_sim: %5.3f\n",m_msg.f_id1,gpercent/100.0);
#endif

    fprintf(fp,"; %s_overlap: %d\n",m_msg.f_id1,lc);
    discons(fp, m_msg, ppst,
	    seqc0, seqc0a, seqc1, seqca, nc,
	    n0, n1, name0, name1, nml, aln);
  }
  else {
#ifndef LALIGN
    fprintf(fp,"%s score: %d; ",m_msg.alabel, score);
#endif
#ifndef SHOWSIM
    fprintf(fp,"%4.1f%% identity (%4.1f%% ungapped) in %d %s overlap (%ld-%ld:%ld-%ld)\n",
	    percent,gpercent,lc,m_msg.sqnam,aln->d_start0,aln->d_stop0,
	    aln->d_start1,aln->d_stop1);
#else
    fprintf(fp,"%4.1f%% identity (%4.1f%% similar) in %d %s overlap (%ld-%ld:%ld-%ld)\n",
	    percent,gpercent,lc,m_msg.sqnam,aln->d_start0,aln->d_stop0,
	    aln->d_start1,aln->d_stop1);
#endif

    if (m_msg.markx & MX_HTML) {
      do_url1(fp, m_msg, ppst, l_name,n1,*aln,aln->l_offset);
    }

    if (m_msg.markx & MX_AMAP && (m_msg.markx & MX_ATYPE)!=7) {
      fputc('\n',fp);
      tmp = n0;

      if (m_msg.qdnaseq == SEQT_DNA && m_msg.ldnaseq== SEQT_PROT)
	tmp /= 3;

      disgraph(fp, tmp, n1, percent, score,
	       aln->amin0, aln->amin1,
	       aln->amax0, aln->amax1,
	       m_msg.sq0off,
	       name0, name1, nml, aln->llen,m_msg.markx);
    }

    discons(fp, m_msg, ppst, 
	    seqc0, seqc0a, seqc1, seqca, nc,
	    n0, n1, name0, name1, nml, aln);

    fputc('\n',fp);

  }
}

void 
do_lav(FILE *fp, struct a_struct *aln_p, char *seqc,
       float percent, int is_mirror) {
  int cur_b0, cur_b1, cur_e0, cur_e1;
  int ipercent;
  long len;
  char *seqc_p, *num_e;

  ipercent = (int)(percent+0.5);

  cur_b0 = aln_p->d_start0;
  cur_b1 = aln_p->d_start1;
  cur_e0 = aln_p->d_stop0;
  cur_e1 = aln_p->d_stop1;

  if (!is_mirror) {
    fprintf (fp, "  b %d %d\n  e %d %d\n", 
	     cur_b0, cur_b1, cur_e0, cur_e1);
  }
  else  {
    fprintf (fp, "  b %d %d\n  e %d %d\n", 
	     cur_b1, cur_b0, cur_e1, cur_e0);
  }

  seqc_p = seqc;
  
  while (*seqc_p) {
    if (*seqc_p == '=') {
      len = strtol(seqc_p+1, &num_e, 10);
      cur_e0 = cur_b0 + len - 1;
      cur_e1 = cur_b1 + len - 1;
      if (!is_mirror) {
	fprintf(fp, "  l %d %d %d %d %d\n",
		cur_b0, cur_b1, cur_e0, cur_e1,
		ipercent);
      }
      else  {
	fprintf(fp, "  l %d %d %d %d %d\n",
		cur_b1, cur_b0, cur_e1, cur_e0,
		ipercent);
      }
      cur_b0 = cur_e0 + 1;
      cur_b1 = cur_e1 + 1;
    }
    else if (*seqc_p == '+') {
      len = strtol(seqc_p+1, &num_e, 10);
      cur_b0 += len;
    }
    else {
      len = strtol(seqc_p+1, &num_e, 10);
      cur_b1 += len;
    }
    seqc_p = num_e;
  }

  fprintf (fp, "}\n");
}


#ifndef MPI_SRC
void	/* initialize consensus arrays */
initseq(char **seqc0, char **seqc0a, char **seqc1, char **seqca, int seqsiz)
{
  *seqc0=(char *)calloc((size_t)seqsiz*4,sizeof(char));
  if (*seqc0==NULL)
    {fprintf(stderr,"cannot allocate consensus arrays %d\n",seqsiz);
     exit(1);}
  *seqc0a=*seqc0 + seqsiz;
  *seqc1=*seqc0a + seqsiz;
  *seqca=*seqc1 + seqsiz;
}

void freeseq(char **seqc0, char **seqc0a, char **seqc1, char **seqca)
{
  free(*seqc0);
}
#endif