source: trunk/GDE/CLUSTALW/pairalign.c

/* Change int h to int gh everywhere  DES June 1994 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "clustalw.h"

#define MIN(a,b) ((a)<(b)?(a):(b))
#define MAX(a,b) ((a)>(b)?(a):(b))

#define gap(k)  ((k) <= 0 ? 0 : g + gh * (k))
#define tbgap(k)  ((k) <= 0 ? 0 : tb + gh * (k))
#define tegap(k)  ((k) <= 0 ? 0 : te + gh * (k))

/*
 *      Prototypes
 */
static void add(sint v);
static sint calc_score(sint iat, sint jat, sint v1, sint v2);
static float tracepath(sint tsb1,sint tsb2);
static void forward_pass(char *ia, char *ib, sint n, sint m);
static void reverse_pass(char *ia, char *ib);
static sint diff(sint A, sint B, sint M, sint N, sint tb, sint te);
static void del(sint k);

/*
 *   Global variables
 */
#ifdef MAC
#define pwint   short
#else
#define pwint   int
#endif
static sint             int_scale;

extern double   **tmat;
extern float    pw_go_penalty;
extern float    pw_ge_penalty;
extern float    transition_weight;
extern sint     nseqs;
extern sint     max_aa;
extern sint     gap_pos1,gap_pos2;
extern sint     max_aln_length;
extern sint     *seqlen_array;
extern sint     debug;
extern sint     mat_avscore;
extern short    blosum30mt[],pam350mt[],idmat[],pw_usermat[],pw_userdnamat[];
extern short    clustalvdnamt[],swgapdnamt[];
extern short    gon250mt[];
extern short    def_dna_xref[],def_aa_xref[],pw_dna_xref[],pw_aa_xref[];
extern Boolean  dnaflag;
extern char     **seq_array;
extern char     *amino_acid_codes;
extern char     pw_mtrxname[];
extern char     pw_dnamtrxname[];

static float    mm_score;
static sint     print_ptr,last_print;
static sint     *displ;
static pwint    *HH, *DD, *RR, *SS;
static sint     g, gh;
static sint     seq1, seq2;
static sint     matrix[NUMRES][NUMRES];
static pwint    maxscore;
static sint     sb1, sb2, se1, se2;


sint pairalign(sint istart, sint iend, sint jstart, sint jend)
{
  short  *mat_xref;
  static sint    si, sj, i;
  static sint    n,m,len1,len2;
  static sint    maxres;
  static short    *matptr;
  static char   c;
  static float gscale,ghscale;

  displ = (sint *)ckalloc((2*max_aln_length+1) * sizeof(sint));
  HH = (pwint *)ckalloc((max_aln_length) * sizeof(pwint));
  DD = (pwint *)ckalloc((max_aln_length) * sizeof(pwint));
  RR = (pwint *)ckalloc((max_aln_length) * sizeof(pwint));
  SS = (pwint *)ckalloc((max_aln_length) * sizeof(pwint));
                
#ifdef MAC
  int_scale = 10;
#else
  int_scale = 100;
#endif
  gscale=ghscale=1.0;
  if (dnaflag)
    {
      if (debug>1) fprintf(stdout,"matrix %s\n",pw_dnamtrxname);
      if (strcmp(pw_dnamtrxname, "iub") == 0)
        { 
          matptr = swgapdnamt;
          mat_xref = def_dna_xref;
        }
      else if (strcmp(pw_dnamtrxname, "clustalw") == 0)
        { 
          matptr = clustalvdnamt;
          mat_xref = def_dna_xref;
          gscale=0.6667;
          ghscale=0.751;
        }
      else
        {
          matptr = pw_userdnamat;
          mat_xref = pw_dna_xref;
        }
      maxres = get_matrix(matptr, mat_xref, matrix, TRUE, int_scale);
      if (maxres == 0) return((sint)-1);

      matrix[0][4]=transition_weight*matrix[0][0];
      matrix[4][0]=transition_weight*matrix[0][0];
      matrix[2][11]=transition_weight*matrix[0][0];
      matrix[11][2]=transition_weight*matrix[0][0];
      matrix[2][12]=transition_weight*matrix[0][0];
      matrix[12][2]=transition_weight*matrix[0][0];
    }
  else
    {
      if (debug>1) fprintf(stdout,"matrix %s\n",pw_mtrxname);
      if (strcmp(pw_mtrxname, "blosum") == 0)
        {
          matptr = blosum30mt;
          mat_xref = def_aa_xref;
        }
      else if (strcmp(pw_mtrxname, "pam") == 0)
        {
          matptr = pam350mt;
          mat_xref = def_aa_xref;
        }
      else if (strcmp(pw_mtrxname, "gonnet") == 0)
        {
          matptr = gon250mt;
          int_scale /= 10;
          mat_xref = def_aa_xref;
        }
      else if (strcmp(pw_mtrxname, "id") == 0)
        {
          matptr = idmat;
          mat_xref = def_aa_xref;
        }
      else
        {
          matptr = pw_usermat;
          mat_xref = pw_aa_xref;
        }

      maxres = get_matrix(matptr, mat_xref, matrix, TRUE, int_scale);
      if (maxres == 0) return((sint)-1);
    }


  for (si=MAX(0,istart);si<nseqs && si<iend;si++)
    {
      n = seqlen_array[si+1];
      len1 = 0;
      for (i=1;i<=n;i++) {
        c = seq_array[si+1][i];
        if ((c!=gap_pos1) && (c != gap_pos2)) len1++;
      }

      for (sj=MAX(si+1,jstart+1);sj<nseqs && sj<jend;sj++)
        {
          m = seqlen_array[sj+1];
          if(n==0 || m==0) {
            tmat[si+1][sj+1]=1.0;
            tmat[sj+1][si+1]=1.0;
            continue;
          }
          len2 = 0;
          for (i=1;i<=m;i++) {
            c = seq_array[sj+1][i];
            if ((c!=gap_pos1) && (c != gap_pos2)) len2++;
          }

          if (dnaflag) {
            g = 2 * (float)pw_go_penalty * int_scale*gscale;
            gh = pw_ge_penalty * int_scale*ghscale;
          }
          else {
            if (mat_avscore <= 0)
              g = 2 * (float)(pw_go_penalty + log((double)(MIN(n,m))))*int_scale;
            else
              g = 2 * mat_avscore * (float)(pw_go_penalty +
                                            log((double)(MIN(n,m))))*gscale;
            gh = pw_ge_penalty * int_scale;
          }

          if (debug>1) fprintf(stdout,"go %d ge %d\n",(pint)g,(pint)gh);

          /*
            align the sequences
          */
          seq1 = si+1;
        seq2 = sj+1;

        forward_pass(&seq_array[seq1][0], &seq_array[seq2][0],
           n, m);

        reverse_pass(&seq_array[seq1][0], &seq_array[seq2][0]);

        last_print = 0;
        print_ptr = 1;
/*
        sb1 = sb2 = 1;
        se1 = n-1;
        se2 = m-1;
*/

/* use Myers and Miller to align two sequences */

        maxscore = diff(sb1-1, sb2-1, se1-sb1+1, se2-sb2+1, 
        (sint)0, (sint)0);
 
/* calculate percentage residue identity */

        mm_score = tracepath(sb1,sb2);

                if(len1==0 || len2==0) mm_score=0;
                else
                        mm_score /= (float)MIN(len1,len2);

        tmat[si+1][sj+1] = ((float)100.0 - mm_score)/(float)100.0;
        tmat[sj+1][si+1] = ((float)100.0 - mm_score)/(float)100.0;

if (debug>1)
{
        fprintf(stdout,"Sequences (%d:%d) Aligned. Score: %d CompScore:  %d\n",
                           (pint)si+1,(pint)sj+1, 
                           (pint)mm_score, 
                           (pint)maxscore/(MIN(len1,len2)*100));
}
else
{
        info("Sequences (%d:%d) Aligned. Score:  %d",
                                      (pint)si+1,(pint)sj+1, 
                                      (pint)mm_score);
}

   }
  }
   displ=ckfree((void *)displ);
   HH=ckfree((void *)HH);
   DD=ckfree((void *)DD);
   RR=ckfree((void *)RR);
   SS=ckfree((void *)SS);


  return((sint)1);
}

static void add(sint v)
{

        if(last_print<0) {
                displ[print_ptr-1] = v;
                displ[print_ptr++] = last_print;
        }
        else
                last_print = displ[print_ptr++] = v;
}

static sint calc_score(sint iat,sint jat,sint v1,sint v2)
{
        sint ipos,jpos;
                sint ret;

        ipos = v1 + iat;
        jpos = v2 + jat;

        ret=matrix[(int)seq_array[seq1][ipos]][(int)seq_array[seq2][jpos]];

        return(ret);
}


static float tracepath(sint tsb1,sint tsb2)
{
        char c1,c2;
    sint  i1,i2,r;
    sint i,k,pos,to_do;
        sint count;
        float score;
        char s1[600], s2[600];

        to_do=print_ptr-1;
        i1 = tsb1;
        i2 = tsb2;

        pos = 0;
        count = 0;
        for(i=1;i<=to_do;++i) {

          if (debug>1) fprintf(stdout,"%d ",(pint)displ[i]);
          if(displ[i]==0) {
            c1 = seq_array[seq1][i1];
            c2 = seq_array[seq2][i2];
            
            if (debug>0)
              {
                if (c1>max_aa) s1[pos] = '-';
                else s1[pos]=amino_acid_codes[c1];
                if (c2>max_aa) s2[pos] = '-';
                else s2[pos]=amino_acid_codes[c2];
              }
            
            if ((c1!=gap_pos1) && (c1 != gap_pos2) &&
                (c1 == c2)) count++;
            ++i1;
            ++i2;
            ++pos;
          }
          else {
            if((k=displ[i])>0) {
              
              if (debug>0)
                for (r=0;r<k;r++)
                  {
                    s1[pos+r]='-';
                    if (seq_array[seq2][i2+r]>max_aa) s2[pos+r] = '-';
                    else s2[pos+r]=amino_acid_codes[seq_array[seq2][i2+r]];
                  }
              
              i2 += k;
              pos += k;
            }
            else {
              
              if (debug>0)
                for (r=0;r<(-k);r++)
                  {
                    s2[pos+r]='-';
                    if (seq_array[seq1][i1+r]>max_aa) s1[pos+r] = '-';
                    else s1[pos+r]=amino_acid_codes[seq_array[seq1][i1+r]];
                  }
              
              i1 -= k;
              pos -= k;
            }
          }
        }
        if (debug>0) fprintf(stdout,"\n");
        if (debug>0) 
          {
            for (i=0;i<pos;i++) fprintf(stdout,"%c",s1[i]);
            fprintf(stdout,"\n");
            for (i=0;i<pos;i++) fprintf(stdout,"%c",s2[i]);
            fprintf(stdout,"\n");
          }
        /*
          if (count <= 0) count = 1;
        */
        score = 100.0 * (float)count;
        return(score);
}


static void forward_pass(char *ia, char *ib, sint n, sint m)
{

  sint i,j;
  pwint f,hh,p,t;

  maxscore = 0;
  se1 = se2 = 0;
  for (i=0;i<=m;i++)
    {
       HH[i] = 0;
       DD[i] = -g;
    }

  for (i=1;i<=n;i++)
     {
        hh = p = 0;
                f = -g;

        for (j=1;j<=m;j++)
           {

              f -= gh; 
              t = hh - g - gh;
              if (f<t) f = t;

              DD[j] -= gh;
              t = HH[j] - g - gh;
              if (DD[j]<t) DD[j] = t;

              hh = p + matrix[(int)ia[i]][(int)ib[j]];
              if (hh<f) hh = f;
              if (hh<DD[j]) hh = DD[j];
              if (hh<0) hh = 0;

              p = HH[j];
              HH[j] = hh;

              if (hh > maxscore)
                {
                   maxscore = hh;
                   se1 = i;
                   se2 = j;
                }
           }
     }

}


static void reverse_pass(char *ia, char *ib)
{

  sint i,j;
  pwint f,hh,p,t;
  pwint cost;

  cost = 0;
  sb1 = sb2 = 1;
  for (i=se2;i>0;i--)
    {
       HH[i] = -1;
       DD[i] = -1;
    }

  for (i=se1;i>0;i--)
     {
        hh = f = -1;
        if (i == se1) p = 0;
        else p = -1;

        for (j=se2;j>0;j--)
           {

              f -= gh; 
              t = hh - g - gh;
              if (f<t) f = t;

              DD[j] -= gh;
              t = HH[j] - g - gh;
              if (DD[j]<t) DD[j] = t;

              hh = p + matrix[(int)ia[i]][(int)ib[j]];
              if (hh<f) hh = f;
              if (hh<DD[j]) hh = DD[j];

              p = HH[j];
              HH[j] = hh;

              if (hh > cost)
                {
                   cost = hh;
                   sb1 = i;
                   sb2 = j;
                   if (cost >= maxscore) break;
                }
           }
        if (cost >= maxscore) break;
     }

}

static int diff(sint A,sint B,sint M,sint N,sint tb,sint te)
{
  sint type;
  sint midi,midj,i,j;
  int midh;
  static pwint f, hh, e, s, t;
  
  if(N<=0)  {
    if(M>0) {
      del(M);
    }
    
    return(-(int)tbgap(M));
  }
  
  if(M<=1) {
    if(M<=0) {
      add(N);
      return(-(int)tbgap(N));
    }
    
    midh = -(tb+gh) - tegap(N);
    hh = -(te+gh) - tbgap(N);
    if (hh>midh) midh = hh;
    midj = 0;
    for(j=1;j<=N;j++) {
      hh = calc_score(1,j,A,B)
        - tegap(N-j) - tbgap(j-1);
      if(hh>midh) {
        midh = hh;
        midj = j;
      }
    }
    
    if(midj==0) {
      del(1);
      add(N);
    }
    else {
      if(midj>1)
        add(midj-1);
      displ[print_ptr++] = last_print = 0;
      if(midj<N)
        add(N-midj);
    }
    return midh;
  }
  
/* Divide: Find optimum midpoint (midi,midj) of cost midh */
  
  midi = M / 2;
  HH[0] = 0.0;
  t = -tb;
  for(j=1;j<=N;j++) {
    HH[j] = t = t-gh;
    DD[j] = t-g;
  }
  
  t = -tb;
  for(i=1;i<=midi;i++) {
    s=HH[0];
    HH[0] = hh = t = t-gh;
    f = t-g;
    for(j=1;j<=N;j++) {
      if ((hh=hh-g-gh) > (f=f-gh)) f=hh;
      if ((hh=HH[j]-g-gh) > (e=DD[j]-gh)) e=hh;
      hh = s + calc_score(i,j,A,B);
      if (f>hh) hh = f;
      if (e>hh) hh = e;
      
      s = HH[j];
      HH[j] = hh;
      DD[j] = e;
    }
  }
  
  DD[0]=HH[0];
  
  RR[N]=0;
  t = -te;
  for(j=N-1;j>=0;j--) {
    RR[j] = t = t-gh;
    SS[j] = t-g;
  }
  
  t = -te;
  for(i=M-1;i>=midi;i--) {
    s = RR[N];
    RR[N] = hh = t = t-gh;
    f = t-g;
    
    for(j=N-1;j>=0;j--) {
      
      if ((hh=hh-g-gh) > (f=f-gh)) f=hh;
      if ((hh=RR[j]-g-gh) > (e=SS[j]-gh)) e=hh;
      hh = s + calc_score(i+1,j+1,A,B);
      if (f>hh) hh = f;
      if (e>hh) hh = e;
      
      s = RR[j];
      RR[j] = hh;
      SS[j] = e;
      
    }
  }
  
  SS[N]=RR[N];
  
  midh=HH[0]+RR[0];
  midj=0;
  type=1;
  for(j=0;j<=N;j++) {
    hh = HH[j] + RR[j];
    if(hh>=midh)
      if(hh>midh || (HH[j]!=DD[j] && RR[j]==SS[j])) {
        midh=hh;
        midj=j;
      }
  }
  
  for(j=N;j>=0;j--) {
    hh = DD[j] + SS[j] + g;
    if(hh>midh) {
      midh=hh;
      midj=j;
      type=2;
    }
  }
  
  /* Conquer recursively around midpoint  */
  
  
  if(type==1) {             /* Type 1 gaps  */
    diff(A,B,midi,midj,tb,g);
    diff(A+midi,B+midj,M-midi,N-midj,g,te);
  }
  else {
    diff(A,B,midi-1,midj,tb,0.0);
    del(2);
    diff(A+midi+1,B+midj,M-midi-1,N-midj,0.0,te);
  }
  
  return midh;       /* Return the score of the best alignment */
}

static void del(sint k)
{
  if(last_print<0)
    last_print = displ[print_ptr-1] -= k;
  else
    last_print = displ[print_ptr++] = -(k);
}