mshowalign.c

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

/* copyright (c) 1996, 1997, 1998, 1999 William R. Pearson and the
   U. of Virginia */

/* $Name: fa35_03_06 $ - $Id: mshowalign.c,v 1.58 2008/01/11 14:49:39 wrp Exp $ */

/* mshowalign.c - show sequence alignments in pvcomplib */

/* 
   this is a merged version of showalign.c that works properly with
   both the comp_lib (serial, threaded) and PCOMPLIB parallel versions
   of the programs.

   In the serial and current threaded versions of the programs,
   showalign gets a list of high scoring sequences and must
   re_getlib() the sequence, do_walign(), and then calculate the
   alignment.

   In the PCOMPLIB parallel versions, the worker programs do the
   aligning, so showalign() must send them the appropriate messages to
   have the alignment done, and then collect the alignment results

*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "defs.h"
#include "msg.h"
#include "structs.h"
#include "param.h"

#ifdef PCOMPLIB
#ifdef PVM_SRC
#include "pvm3.h"
extern int pinums[];
#endif
#ifdef MPI_SRC
#include "mpi.h"
#endif
#include "p_mw.h"
#else
#include "mm_file.h"
#include "mw.h"
#endif

#ifndef PCOMPLIB

/* used to position the library sequence for re_getlib - also gets
   description */
#define RANLIB (m_fptr->ranlib)

extern struct lmf_str *
re_openlib(struct lmf_str *, int outtty);

int
re_getlib(unsigned char *aa1, int maxn, int maxt,
	  int loff, int cont, int term_code,
	  long *l_offset, long *l_off, 
	  struct lmf_str *m_fptr);

#include "drop_func.h"

#endif


extern void cal_coord(int n0, int n1, long qoffset, long loffset,
		      struct a_struct *aln);

void initseq(char **, char **, char **, char **, int);
void freeseq(char **, char **, char **, char **);

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);

void
do_lav(FILE *fp, struct a_struct *aln, char *seqc, float percent, int is_mirror);

extern void discons(FILE *fd, struct mngmsg m_msg, struct pstruct *ppst,
		    char *seqc0, char *seqc0a, char *seqc1, char *seqca,
		    int nc, 
		    int n0, int n1, char *name0, char *name1, int nml,
		    struct a_struct *aln);

extern void disgraph(FILE *fd, int n0, int n1,
		     float percent, int score,
		     int min0, int min1, int max0, int max1, long sq0off,
		     char *name0, char *name1, int nml, int llen, int markx);

extern double (*find_zp)(int score, double escore, int length, double comp,void *);
extern double zs_to_bit(double, int, int);
extern double s_to_bit(int score, int n0, int  n1, void *pu);
extern double bit_to_E (double bit, int n0, int n1, long db_size, void *pu);
extern double zs_to_E(double zs, int n1, int dnaseq, long db_size, struct db_str db);

extern void
do_url1(FILE *, struct mngmsg, struct pstruct *, char *, int,
	struct a_struct , long);

#ifndef A_MARK
#define A_MARK ">>"
#endif

static char l_name[200];	/* link name */

#ifdef PCOMPLIB
#define BBP_INFO(info) bbp->seq->info
#else
#define BBP_INFO(info) bbp->seq->info
#endif

/* this version does not check for m_msg.e_cut because nshow/nbest has
   already been set to limit on e_cut */

void showalign (FILE *fp,
#ifndef PCOMPLIB
		unsigned char **aa0, unsigned char *aa1save, int maxn,
#endif
		struct beststr **bptr, int nbest, int qlib, 
		struct mngmsg m_msg, struct pstruct *ppst, 
		char *info_gstring2
#ifndef PCOMPLIB
		, void **f_str
#endif
)
{
  unsigned char *aa1;
  char tmp_str[20];
  char info_str[200];
  char bline[2048], *bl_ptr, *bp, fmt[40];
  int tmp_len, l_llen;
  int t_have_ares;
  char name0[80], name0s[80], name1[200];
  int istart = 0, istop, i = 0, ib, nml;
  int  n1tot;
  struct beststr *bbp;
  struct a_res_str *cur_ares;
  int nc, lc, maxc, maxca, seqc_max, annc_max;
  double zscore, bits;
  struct a_struct *aln_p;
  float percent, gpercent;
  /* strings, lengths for conventional alignment */
  char *seqc0, *seqc0a, *seqc1, *seqca, *seqc;
  int seqc_len, seqca_len;
  /* strings, lengths, for encoded alignment for MX10 */
  char *seq_code=NULL, *ann_code=NULL;
  int seq_code_len=0, ann_code_len=0;
  long loffset, l_off;
  int lsw_score;
#ifdef PCOMPLIB
  struct stage2_str liblist;
  struct qmng_str qm_msg;
#ifdef MPI_SRC
  int int_msg_b[10];
  MPI_Status mpi_status;
#endif
#else
  int n1;
  struct lmf_str *m_fptr;
  int ngap;
#endif

#ifdef PCOMPLIB
  /* this function has its own copy of qm_msg, so we must fill it
     appropriately */
  qm_msg.n0 = m_msg.n0;
  strncpy(qm_msg.libstr,m_msg.qtitle,sizeof(qm_msg.libstr));
#endif

  aln_p = &m_msg.aln;
  /* set the name0,1 label length */
  if (m_msg.markx & MX_M10FORM) nml = 12;
  else nml = m_msg.nmlen;

  if (strlen(m_msg.qtitle) > 0) {
    if (m_msg.qtitle[0]=='>') strncpy(name0s,&m_msg.qtitle[1],sizeof(name0s));
    else strncpy(name0s,m_msg.qtitle,sizeof(name0s));
  }
  else {
    strncpy(name0s,m_msg.tname,sizeof(name0s));
  }
  name0s[sizeof(name0s)-1]='\0';

  if ((bp=strchr(name0s,' '))!=NULL) *bp='\0';

  if (m_msg.revcomp) name0[nml-1]='-';

  l_llen = m_msg.aln.llen;
  if ((m_msg.markx & MX_M9SUMM) && m_msg.show_code != SHOW_CODE_ID) {
    l_llen += 40;
    if (l_llen > 200) l_llen=200;
  }

  sprintf(fmt,"%s%%-%ds (%%d %s)\n",A_MARK,l_llen-5,m_msg.sqnam);

  if (!(m_msg.markx&MX_M10FORM)) fprintf (fp,"\n");

  if (m_msg.ashow < 0) m_msg.ashow = m_msg.nshow;
  istart = 0; istop = min(min(nbest,m_msg.ashow),m_msg.nshow);

  for (ib=istart; ib<istop; ib++) {
    bbp = bptr[ib];

#ifdef SHOWUN
    if (BBP_INFO(nsfnum) > 0 && sfn_cmp(m_msg.qsfnum,BBP_INFO(sfnum))) {
      istop = min(istop+1,nbest);
      continue;
    }
#endif
    if (m_msg.quiet==1 && ppst->zsflag>=0 
	&& bbp->rst.escore < m_msg.e_low) continue;

#ifndef PCOMPLIB
    /* get the alignment and score by re-aligning */

    if ((m_fptr=re_openlib(bbp->seq->m_file_p,!m_msg.quiet))==NULL)
      exit(1);

    /* get the description - do not "edit" it yet */

    if (!(m_msg.markx & MX_M10FORM)){
      if (m_msg.long_info) {tmp_len = sizeof(bline)-1;}
      else {tmp_len = l_llen-5;}
      RANLIB(bline,tmp_len,bbp->seq->lseek,bbp->seq->libstr,bbp->seq->m_file_p);
      bline[tmp_len]='\0';
    }
    else {
      RANLIB(bline,sizeof(bline),bbp->seq->lseek,bbp->seq->libstr,bbp->seq->m_file_p);
      bline[sizeof(bline)-1]='\0';
    }

    if (bbp->seq->aa1b == NULL) {

      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;
      l_off = bbp->seq->l_off;
      loffset = bbp->seq->l_offset;
    }

#ifdef DEBUG
    if (n1 != BBP_INFO(n1)) {
      fprintf(stderr," library sequence: %s lengths differ: %d != %d\n",
	      bline,BBP_INFO(n1), n1);
      fprintf(stderr, "offset is: %lld\n",bbp->seq->lseek);
    }
#endif

    if (!bbp->have_ares) {
      do_walign(aa0[bbp->frame],m_msg.n0, aa1, n1, bbp->frame, ppst,
		f_str[bbp->frame], &bbp->a_res, &t_have_ares);
      bbp->have_ares = t_have_ares;
    }
    else {
      pre_cons(aa1,n1,bbp->frame,f_str[bbp->frame]);
    }

    cur_ares = &bbp->a_res;

    aln_func_vals(bbp->frame, aln_p);

#else	/* PCOMPLIB - get the alignment information from a worker */

    /* we have a sequence that we need an alignment for -
       send a message to the appropriate worker to produce an alignment 
       qm_msg.slist == 1  -> one alignment
       qm_msg.s_func == DO_ALIGN_FLG -> use the alignment function
       send mngmsg (MSEQTYPE)
       then send number of sequence to be aligned
    */

    qm_msg.slist = 1;
    qm_msg.s_func = DO_ALIGN_FLG;

    liblist.seqnm = bbp->seq->seqnm;
    liblist.frame = bbp->frame;
#ifdef PVM_SRC
    pvm_initsend(PvmDataRaw);
    pvm_pkbyte((char *)&qm_msg,sizeof(struct qmng_str),1);
    pvm_send(pinums[bbp->seq->wrkr],MSEQTYPE);

    pvm_initsend(PvmDataRaw);
    pvm_pkbyte((char *)&liblist,sizeof(struct stage2_str),1);
    pvm_send(pinums[bbp->seq->wrkr],LISTTYPE);
#endif
#ifdef MPI_SRC
    MPI_Send(&qm_msg,sizeof(struct qmng_str),MPI_BYTE,bbp->seq->wrkr,
	     MSEQTYPE,MPI_COMM_WORLD);
    MPI_Send(&liblist,sizeof(struct stage2_str),MPI_BYTE,bbp->seq->wrkr,
	     LISTTYPE,MPI_COMM_WORLD);
#endif
    /* information should be sent */
    /* pick up description */
    strncpy(bline,bbp->seq->bline,l_llen-5);
    bline[l_llen-5]='\0';
#endif	/* PCOMPLIB */

    if (strlen(bline)==0) {
      bline[0]='>';
      strncpy(&bline[1],m_msg.lname,l_llen-5);
      bline[l_llen-5]='\0';
    }

    /* re-format bline */
    while ((bp=strchr(bline,'\n'))!=NULL) *bp=' ';
    if (m_msg.long_info) {
      tmp_len = strlen(bline);
      bl_ptr = bline;
      if (!(m_msg.markx & MX_M10FORM)) while (tmp_len > l_llen) {
	for (i=l_llen; i>10; i--)
	  if (bl_ptr[i]==' ') {
	    bl_ptr[i]='\n';
	    break;
	  }
	if (i <= 10) break;
	tmp_len -= i;
	bl_ptr += i;
      }
      bline[sizeof(bline)-1]='\0';
    }

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

    strncpy(name1,bline,sizeof(name1));
    name1[sizeof(name1)-1] = '\0';

    if ((!m_msg.markx & MX_M10FORM)) name1[nml]='\0';
    if ((bp = strchr(name1,' '))!=NULL) *bp = '\0';

  /* l_name is used to build an HTML link from the bestscore line to
     the alignment.  It can also be used to discriminate multiple hits
     from the same long sequence.  Text must match that in showbest.c */

    strncpy(name1,bline,sizeof(name1));
    name1[sizeof(name1)-1]='\0';
    if ((bp = strchr(name1,' '))!=NULL) *bp = '\0';
    strncpy(l_name,name1,sizeof(l_name));
    l_name[sizeof(l_name)-1]='\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.qframe >= 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));
    }

    if (!(m_msg.markx & MX_M10FORM)) name1[nml]='\0';

    /* print out score information; */

    if (m_msg.markx & MX_HTML ) {
      fprintf (fp,"<A name=%s>\n<tt><pre>\n",l_name);
    }
    strncpy(name0,name0s,nml);
    name0[nml]='\0';

    if (ppst->zsflag%10 == 6) {
      sprintf(info_str," comp: %.5f H: %.5f",bbp->rst.comp,bbp->rst.H);
    }
    else info_str[0]='\0';

    if (m_msg.markx & MX_M11OUT) {
      fprintf (fp, "s {\n   \"%s\" %ld %ld \n   \"%s\" %ld %ld\n}\n",
	       name0, m_msg.q_off,  m_msg.q_off+m_msg.n0-1, 
	       name1, bbp->seq->l_off, bbp->seq->l_off + bbp->seq->n1 - 1);
      fprintf (fp, "h {\n   \"%s\"\n   \"%s\"\n}\n", m_msg.qtitle, bline);
    }
    else if ((m_msg.markx & MX_ATYPE)!=7 && !(m_msg.markx & MX_M10FORM)) {
      fprintf (fp, fmt,bp=bline,n1tot);	 /* provides >>id  description (length) line */
    }
    else if (m_msg.markx & MX_M10FORM) {
      fprintf (fp,">>%s\n",bline);
    }

#ifdef PCOMPLIB
    /*  get the sw_score, alignment  information,  get seqc0, seqc1 */

#ifdef PVM_SRC
    /* get alignment lengths, percents */
    pvm_recv(pinums[bbp->seq->wrkr],ALN1TYPE);
    pvm_upkint(&nc,1,1);
    pvm_upkint(&lc,1,1);
    pvm_upkint(&maxc,1,1);

    pvm_upkfloat(&percent,1,1);
    pvm_upkfloat(&gpercent,1,1);

    pvm_upkint(&bbp->sw_score,1,1);
    pvm_upkbyte((char *)aln_p,sizeof(struct a_struct),1);