pairalign.c

是有关基因比对的经典算法的实现。这对于初学计算生物学的人是非常重要的算法。

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