home *** CD-ROM | disk | FTP | other *** search
/ Crawly Crypt Collection 1 / crawlyvol1.bin / apps / science / clustalv / myers.c < prev    next >
C/C++ Source or Header  |  1993-04-11  |  11KB  |  539 lines
  1. #include <stdio.h>
  2. #include <string.h>
  3. #include <stdlib.h>
  4. #include "clustalv.h"
  5.  
  6. /*
  7. *    Prototypes
  8. */
  9.  
  10. extern Boolean read_tree(int *,double *,int *,int *,int *);
  11. extern void *ckalloc(size_t);
  12.  
  13. void init_myers(void);
  14. void myers(int);
  15.  
  16. static void group_gap(int,int,char *);
  17. static void add_ggaps(void);
  18. static void add(int);
  19. static int calc_weight(int,int,int,int);
  20. static void fill_pam(void);
  21. static int diff(int,int,int,int,int,int);
  22. static void do_align(int,int *);
  23. static int *alist;
  24. static int    ** naa1,**naa2,**naas;
  25.  
  26. /*
  27. *    Global variables
  28. */
  29.  
  30. extern int nseqs, next, profile1_nseqs;    /* DES */
  31. extern int *seqlen_array;
  32. extern char **seq_array;
  33. extern int pamo[];
  34. extern Boolean dnaflag,percent;
  35. extern int xover,big_pam;
  36. extern int nblocks;
  37. extern FILE *clustal_outfile,*tree;
  38. extern char treename[],seqname[],mtrxnam[],**names;
  39. extern char *matptr,idmat[],pam100mt[],pam250mt[];
  40. extern int ktup,wind_gap,window;
  41. extern int signif;
  42.  
  43.  
  44. int      gap_open,      gap_extend;
  45. int  dna_gap_open,  dna_gap_extend;
  46. int prot_gap_open, prot_gap_extend;
  47.  
  48. Boolean is_weight;
  49.  
  50. extern int *zza,*zzb,*zzc,*zzd;        /* allocated in show_pair.c     */
  51. static int *group;
  52. static int print_ptr,last_print,pos1,pos2;
  53. extern int * displ;                    /* allocated in show_pair.c        */
  54. static int pam[21][21];
  55.  int    weights[21][21];
  56. static int *fst_list,*snd_list;
  57.  
  58. #define gap(x)  ((x) <= 0 ? 0 : gap_open + gap_extend * (x))
  59.  
  60. void init_myers()
  61. {
  62.     register int i;
  63.     
  64.     group = (int *)ckalloc( (MAXN+1) * sizeof (int));
  65.  
  66.     alist = (int *)ckalloc( (MAXN+1) * sizeof (int));
  67.     fst_list = (int *)ckalloc( (MAXN+1) * sizeof (int));
  68.     snd_list = (int *)ckalloc( (MAXN+1) * sizeof (int));
  69.  
  70.     naa1 = (int **)ckalloc(21 * sizeof (int *) );
  71.     for(i=0;i<21;i++)
  72.         naa1[i]=(int *)ckalloc( (MAXLEN+1)*sizeof (int));
  73.     naa2 = (int **)ckalloc(21 * sizeof (int *) );
  74.     for(i=0;i<21;i++)
  75.         naa2[i]=(int *)ckalloc( (MAXLEN+1)*sizeof (int));
  76.     naas = (int **)ckalloc(21 * sizeof (int *) );
  77.     for(i=0;i<21;i++)
  78.         naas[i]=(int *)ckalloc( (MAXLEN+1)*sizeof (int));
  79.  
  80.     dna_gap_open    = 10;    /* Default gap penalties for DNA */
  81.     dna_gap_extend  = 10;
  82.  
  83.     prot_gap_open   = 10;    /* Default gap penalties for protein */
  84.     prot_gap_extend = 10;
  85.  
  86.     is_weight = TRUE;
  87. }
  88.  
  89.  
  90.  
  91. static void group_gap(int len,int sclass,char *seq)
  92. {
  93.     int i,j,k,xtra;
  94.     
  95.     for(i=1;i<=nseqs;++i)
  96.         if(group[i] == sclass) {
  97.             xtra = len - seqlen_array[i];
  98.             if(xtra>0)
  99.                 for(j=1;j<=xtra;++j)
  100.                     seq_array[i][seqlen_array[i]+j] = -1;
  101.             for(j=1;j<=len;++j)
  102.                 if(seq[j] == '-') {
  103.                     for(k=len;k>=j+1;--k)
  104.                     seq_array[i][k] = seq_array[i][k-1];
  105.                     seq_array[i][j] = -1;
  106.                 }
  107.             seqlen_array[i] = len;
  108.         }
  109. }
  110.  
  111.  
  112. static void add_ggaps()
  113. {
  114.     int i,j,k,pos,to_do,len;
  115.     char str1[MAXLEN+1],str2[MAXLEN+1];
  116.     
  117.     pos=1;
  118.     to_do=print_ptr-1;
  119.     
  120.     for(i=1;i<=to_do;++i) {
  121.         if(displ[i]==0) {
  122.             str1[pos]=str2[pos]='*';
  123.             ++pos;
  124.         }
  125.         else {
  126.             if((k=displ[i])>0) {
  127.                 for(j=0;j<=k-1;++j) {
  128.                     str2[pos+j]='*';
  129.                     str1[pos+j]='-';
  130.                 }
  131.                 pos += k;
  132.             }
  133.             else {
  134.                 k = (displ[i]<0) ? displ[i] * -1 : displ[i];
  135.                 for(j=0;j<=k-1;++j) {
  136.                     str1[pos+j]='*';
  137.                     str2[pos+j]='-';
  138.                 }
  139.             pos += k;
  140.             }
  141.         }
  142.     }
  143.     
  144.     len = --pos;
  145.     group_gap(len,1,str1);
  146.     group_gap(len,2,str2);
  147. }
  148.  
  149.  
  150. static void add(int v)
  151. {
  152.     
  153.     if(last_print<0) {
  154.         displ[print_ptr-1] = v;
  155.         displ[print_ptr++] = last_print;
  156.     }
  157.     else 
  158.         last_print = displ[print_ptr++] = v;
  159. }
  160.  
  161.  
  162. static int calc_weight(int iat,int jat,int v1,int v2)
  163. {
  164.     int sum,i,j,lookn,ret;
  165.     int ipos,jpos;
  166.     
  167.     sum=lookn=0;
  168.     ipos = v1 + iat -1;
  169.     jpos = v2 + jat -1;
  170.     
  171.     ret=0;
  172.     if(pos1>=pos2) {
  173.         for(i=1;i<=pos2;++i) {
  174.             j=seq_array[alist[i]][jpos];
  175.             if(j>0) {
  176.                 sum += naas[j][ipos];
  177.                 ++lookn;
  178.             }
  179.         }    
  180.     }
  181.     else {
  182.         for(i=1;i<=pos1;++i) {
  183.             j=seq_array[alist[i]][ipos];
  184.             if(j>0) {
  185.                 sum += naas[j][jpos];
  186.                 ++lookn;
  187.             }
  188.         }
  189.     }
  190.     
  191.     if(sum > 0 ) ret = sum / lookn;
  192.     return ret;
  193. }
  194.  
  195.  
  196. static void fill_pam()
  197. {
  198.     register int i,j,pos;
  199.     
  200.     pos=0;
  201.     
  202.     for(i=0;i<20;++i)
  203.         for(j=0;j<=i;++j)
  204.             pam[i][j]=pamo[pos++];
  205.     
  206.     for(i=0;i<20;++i)
  207.         for(j=0;j<=i;++j)
  208.             pam[j][i]=pam[i][j];
  209.     
  210.     if(dnaflag) {
  211.         xover=4;
  212.         big_pam=8;
  213.         for(i=0;i<5;++i)
  214.             for(j=0;j<5;++j) {
  215.                 if(i==j)
  216.                     weights[i][j]=0;
  217.                 else
  218.                     weights[j][i]=10;
  219.             }
  220.         if(is_weight) {
  221.             weights[1][3]=4;
  222.             weights[3][1]=4;
  223.             weights[2][4]=4;
  224.             weights[4][2]=4;
  225.         }
  226.     }
  227.     else {
  228. /*
  229.     fprintf(stdout,"\nxover = %d; big_pam = %d\n",xover,big_pam);
  230. */
  231.         for(i=1;i<21;++i)
  232.             for(j=1;j<21;++j) {
  233.                 weights[i][j] = big_pam - pam[i-1][j-1];
  234. /*
  235.                 fprintf(stdout,"\n%2d vs %2d:  %5d",i,j,weights[i][j]);
  236. */
  237.             }
  238.         for(i=0;i<21;++i) {
  239.             weights[0][i] = xover;
  240.             weights[i][0] = xover;
  241.         }
  242.     }
  243. }
  244.  
  245. static int diff(int v1,int v2,int v3,int v4,int st,int en)
  246. {
  247.     int ctrc,ctri,ctrj,i,j,k,l,m,n,p,q,flag;
  248.     
  249.     q  = gap_open + gap_extend;
  250.     
  251.     if(v4<=0)  {
  252.         if(v3>0) {
  253.             if(last_print<0)
  254.                 last_print = displ[print_ptr-1] -= v3;
  255.             else
  256.                 last_print = displ[print_ptr++] = -(v3);
  257.         }
  258.     
  259.                 return gap(v3);
  260.     }
  261.     
  262.     if(v3<=1) {
  263.         if(v3<=0) {
  264.             add(v4);
  265.                         return gap(v4);
  266.         }
  267.         if(st>en)
  268.             st=en;
  269.  
  270. /***************if(!v4)*********BUG********************************/
  271.  
  272.         ctrc = (st+gap_extend) + gap(v4);
  273.         ctrj = 0;
  274.         for(j=1;j<=v4;++j) {
  275.             k = calc_weight(1,j,v1,v2) + gap(v4-j) + gap(j-1);
  276.             if(k<ctrc) {
  277.                 ctrc = k;
  278.                 ctrj = j;
  279.             }
  280.         }
  281.  
  282.         if(!ctrj) {
  283.             add(v4);
  284.             if(last_print<0)
  285.                 last_print = displ[print_ptr-1] -= 1;
  286.             else
  287.                 last_print = displ[print_ptr++] = -1;
  288.         }
  289.         else {
  290.             if(ctrj>1)
  291.                 add(ctrj-1);
  292.             displ[print_ptr++] = last_print = 0;
  293.             if(ctrj<v4)
  294.                 add(v4-ctrj);
  295.         }
  296.         return ctrc;
  297.     }
  298.     
  299.     
  300.     ctri = v3 / 2;
  301.     zza[0] = 0;
  302.     p = gap_open;
  303.     for(j=1;j<=v4;++j) {
  304.         p += gap_extend;
  305.         zza[j] = p;
  306.         zzb[j] = p + gap_open;
  307.     }
  308.     
  309.     p=st;
  310.     for(i=1;i<=ctri;++i) {
  311.         n=zza[0];
  312.         p += gap_extend;
  313.         k = p;
  314.         zza[0]=k;
  315.         l = p+gap_open;
  316.         for(j=1;j<=v4;++j) {
  317.             k += q;
  318.             l += gap_extend;
  319.             if(k<l)
  320.                 l=k;
  321.             k = zza[j] + q;
  322.             m = zzb[j] + gap_extend;
  323.             if(k<m)
  324.                 m=k;
  325.             k = n + calc_weight(i,j,v1,v2);
  326.             if(l<k)
  327.                 k=l;
  328.             if(m<k)
  329.                 k=m;
  330.             n=zza[j];
  331.             zza[j]=k;
  332.             zzb[j]=m;
  333.         }
  334.     }
  335.     
  336.     zzb[0]=zza[0];
  337.     zzc[v4]=0;
  338.     p=gap_open;
  339.     for(j=v4-1;j>-1;--j) {
  340.         p += gap_extend;
  341.         zzc[j]=p;
  342.         zzd[j]=p+gap_open;
  343.     }
  344.     p=en;
  345.     for(i=v3-1;i>=ctri;--i) {
  346.         n=zzc[v4];
  347.         p += gap_extend;
  348.         k = p;
  349.         zzc[v4] = k;
  350.         l = p+gap_open;
  351.         for(j=v4-1;j>=0;--j) {
  352.             k += q;
  353.             l += gap_extend;
  354.             if(k<l)
  355.                 l=k;
  356.             k = zzc[j] + q;
  357.             m = zzd[j] + gap_extend;
  358.             if(k<m)
  359.                 m=k;
  360.             k = n + calc_weight(i+1,j+1,v1,v2);
  361.             if(l<k)
  362.                 k=l;
  363.             if(m<k)
  364.                 k=m;
  365.             n=zzc[j];
  366.             zzc[j]=k;
  367.             zzd[j]=m;
  368.         }
  369.     }
  370.  
  371.     zzd[v4]=zzc[v4];
  372.     ctrc=zza[0]+zzc[0];
  373.     ctrj=0;
  374.     flag=1;
  375.     for(j=0;j<=v4;++j) {
  376.         k = zza[j] + zzc[j];
  377.         if(k<=ctrc)
  378.             if(k<ctrc || ((zza[j]!=zzb[j]) && (zzc[j]==zzd[j]))) {
  379.                 ctrc=k;
  380.                 ctrj=j;
  381.             }
  382.     }
  383.  
  384.     for(j=v4;j>=0;--j) {
  385.         k = zzb[j] + zzd[j] - gap_open;
  386.         if(k<ctrc) {
  387.             ctrc=k;
  388.             ctrj=j;
  389.             flag=2;
  390.         }
  391.     }
  392.  
  393.     /* Conquer recursively around midpoint  */
  394.  
  395.     if(flag==1) {             /* Type 1 gaps  */
  396.         diff(v1,v2,ctri,ctrj,st,gap_open);
  397.         diff(v1+ctri,v2+ctrj,v3-ctri,v4-ctrj,gap_open,en);
  398.     }
  399.     else {
  400.         diff(v1,v2,ctri-1,ctrj,st,0);
  401.         if(last_print<0)                         /* Delete 2 */
  402.             last_print = displ[print_ptr-1] -= 2;
  403.         else
  404.             last_print = displ[print_ptr++] = -2;
  405.         diff(v1+ctri+1,v2+ctrj,v3-ctri-1,v4-ctrj,0,en);
  406.     }
  407.     return ctrc;       /* Return the score of the best alignment */
  408. }
  409.  
  410.  
  411.  
  412.  
  413.  
  414.  
  415. static void do_align(int v1,int *score)
  416. {
  417.     int ctrc,i,j,k,l1,l2,n;
  418.     int t_arr[21];
  419.  
  420.     l1=l2=pos1=pos2=0;
  421.  
  422.     for(i=1;i<=MAXLEN;++i) {
  423.         for(j=0;j<21;++j)
  424.             naa1[j][i]=naa2[j][i]=0;
  425.         for(j=1;j<21;++j)
  426.             naas[j][i]=0;
  427.     }
  428.  
  429.     for(i=1;i<=nseqs;++i) {
  430.         if(group[i]==1) {
  431.             fst_list[++pos1]=i;
  432.             for(j=1;j<=seqlen_array[i];++j)
  433.                 if(seq_array[i][j]>0) {
  434.                  ++naa1[seq_array[i][j]][j];
  435.                  ++naa1[0][j];
  436.                 }
  437.             if(seqlen_array[i]>l1)
  438.                 l1=seqlen_array[i];
  439.         }
  440.         else if(group[i]==2) {
  441.             snd_list[++pos2]=i;
  442.             for(j=1;j<=seqlen_array[i];++j)
  443.                 if(seq_array[i][j]>0) {
  444.                     ++naa2[seq_array[i][j]][j];
  445.                     ++naa2[0][j];
  446.                 }
  447.             if(seqlen_array[i]>l2) l2=seqlen_array[i];
  448.         }
  449.     }
  450.  
  451.     if(pos1>=pos2) {
  452.         for(i=1;i<=pos2;++i)
  453.             alist[i]=snd_list[i];
  454.         for(n=1;n<=l1;++n) {
  455.             for(i=1;i<21;++i)
  456.                 t_arr[i]=0;
  457.             for(i=1;i<21;++i)
  458.                 if(naa1[i][n]>0)
  459.                     for(j=1;j<21;++j)
  460.                         t_arr[j] += (weights[i][j]*naa1[i][n]);
  461.             k = naa1[0][n];
  462.             if(k>0)
  463.                 for(i=1;i<21;++i)
  464.                     naas[i][n]=t_arr[i]/k;
  465.         }
  466.     }
  467.     else {
  468.         for(i=1;i<=pos1;++i)
  469.             alist[i]=fst_list[i];
  470.         for(n=1;n<=l2;++n) {
  471.             for(i=1;i<21;++i)
  472.                 t_arr[i]=0;
  473.             for(i=1;i<21;++i)
  474.                 if(naa2[i][n]>0)
  475.                     for(j=1;j<21;++j)
  476.                         t_arr[j] += (weights[i][j]*naa2[i][n]);
  477.             k = naa2[0][n];
  478.             if(k>0)
  479.                 for(i=1;i<21;++i)
  480.                     naas[i][n]=t_arr[i]/k;
  481.         }
  482.     }
  483.  
  484.     *score=diff(1,1,l1,l2,v1,v1);    /* Myers and Miller alignment now */
  485. }
  486.  
  487.  
  488. void myers(int align_type)   /* align_type = 0 for full progressive alignment*/
  489.                  /* align_type = 1 for a profile alignment */
  490. {
  491. /*    static char *nbases = "XACGT";
  492.     static char seq1[MAXLEN+1];     */
  493.     char        temp[MAXLINE];
  494.     int val,i,j,k,a,b,len,set;
  495.     int sets,chunks,ident,lv1,pos,copt,flag,entries,idummy,score,ptr;
  496.     double dummy;
  497.  
  498.     nblocks=0;
  499.  
  500.     fprintf(stdout,"\nStart of Multiple Alignment\n");
  501.  
  502.     fill_pam();
  503.  
  504.     if(align_type == 0)  {        /* a full progressive alignment */
  505.         sets=0;
  506.         tree=fopen(treename,"r");
  507.         while(fgets(temp,MAXLINE,tree)!=NULL) ++sets;
  508.         fseek(tree,0,0);
  509.         fprintf(stdout,"There are %d sets\n",sets);
  510.     }
  511.     else                /* just one set (a profile alignment) */
  512.         sets = 1;
  513.  
  514.     fprintf(stdout,"Aligning...\n");
  515.  
  516.     for(set=1;set<=sets;++set) {
  517.         if(align_type == 0)
  518.             read_tree(&entries,&dummy,&idummy,&idummy,group);
  519.         else  {
  520.             for(i=1; i<=profile1_nseqs; ++i)
  521.                 group[i] = 1;
  522.             for(i=profile1_nseqs+1; i<=nseqs; ++i)
  523.                  group[i] = 2;
  524.              entries = nseqs;
  525.         }
  526.         last_print=0;
  527.         print_ptr=1;
  528.         do_align(gap_open,&score);
  529.         fprintf(stdout,"Set %d: Entries:%d    Score:%d\n",set,entries,score);
  530.         add_ggaps();
  531.     }
  532.     if(align_type == 0) fclose(tree);
  533.  
  534. /* make the rest (output stuff) into routine clustal_out in file amenu.c */
  535.  
  536. }
  537.  
  538.  
  539.