home *** CD-ROM | disk | FTP | other *** search
/ Crawly Crypt Collection 1 / crawlyvol1.bin / apps / science / clustalv / showpair.c < prev    next >
C/C++ Source or Header  |  1993-04-11  |  9KB  |  476 lines
  1. #include <stdio.h>
  2. #include <string.h>
  3. #include <stdlib.h>
  4. #include <math.h>
  5. #include "clustalv.h"    
  6.  
  7.  
  8. /*
  9. *    Prototypes
  10. */
  11.  
  12. extern void *ckalloc(size_t);
  13. extern void error(char *,...);
  14. extern int *seqlen_array;
  15. extern char **seq_array;
  16.  
  17. void init_show_pair(void);
  18. void show_pair(void);
  19. static void make_p_ptrs(int *,int *,int,int);
  20. static void make_n_ptrs(int *,int *,int,int);
  21. static void put_frag(int,int,int,int);
  22. static int frag_rel_pos(int,int,int,int);
  23. static void pair_align(int,int,int);
  24. static void des_quick_sort(int *, int *, int);
  25.  
  26.  
  27. /*
  28. *     Global variables
  29. */
  30.  
  31. extern int next,nseqs;
  32. extern Boolean dnaflag;
  33. extern double **tmat;
  34.  
  35. int ktup,window,wind_gap,signif;                  /* Pairwise aln. params */
  36. int  dna_ktup, dna_window, dna_wind_gap, dna_signif;  /* params for DNA */
  37. int prot_ktup,prot_window,prot_wind_gap,prot_signif;  /* params for prots */
  38.  
  39.  
  40. Boolean percent;
  41.  
  42. static int curr_frag,maxsf,vatend;
  43. static int **accum;
  44. int *displ;                        /* also used in myers.c     */
  45. int *zza, *zzb, *zzc, *zzd;        /* also used in myers.c     */
  46. static int *diag_index;
  47. static char *slopes;
  48.  
  49. void init_show_pair(void)
  50. {
  51.     register int i;
  52.  
  53.     accum = (int **)ckalloc( 5*sizeof (int *) );
  54.     for (i=0;i<5;i++)
  55.         accum[i] = (int *) ckalloc(FSIZE * sizeof (int) );
  56.  
  57.     displ      = (int *) ckalloc( (2*MAXLEN +1) * sizeof (int) );
  58.     slopes     = (char *)ckalloc( (2*MAXLEN +1) * sizeof (char));
  59.     diag_index = (int *) ckalloc( (2*MAXLEN +1) * sizeof (int) );
  60.  
  61.     zza = (int *)ckalloc( (MAXLEN+1) * sizeof (int) );
  62.     zzb = (int *)ckalloc( (MAXLEN+1) * sizeof (int) );
  63.  
  64.     zzc = (int *)ckalloc( (MAXLEN+1) * sizeof (int) );
  65.     zzd = (int *)ckalloc( (MAXLEN+1) * sizeof (int) );
  66.  
  67.     dna_ktup      = 2;   /* default parameters for DNA */
  68.     dna_wind_gap  = 5;
  69.     dna_signif    = 4;
  70.     dna_window    = 4;
  71.  
  72.     prot_ktup     = 1;   /* default parameters for proteins */
  73.     prot_wind_gap = 3;
  74.     prot_signif   = 5;
  75.     prot_window   = 5;
  76.  
  77.     percent=TRUE;
  78. }
  79.  
  80. static void make_p_ptrs(int *tptr,int *pl,int naseq,int l)
  81. {
  82.     static int a[]={ 0, 1, 20, 400 };
  83.     int i,j,limit,code,flag;
  84.     char residue;
  85.     
  86.     limit = (int) pow((double)20,(double)ktup);
  87.     for(i=1;i<=limit;++i)
  88.         pl[i]=0;
  89.     for(i=1;i<=l;++i)
  90.         tptr[i]=0;
  91.     
  92.     for(i=1;i<=(l-ktup+1);++i) {
  93.         code=0;
  94.         flag=FALSE;
  95.         for(j=1;j<=ktup;++j) {
  96.             residue = seq_array[naseq][i+j-1];
  97.             if(residue<=0) {
  98.                 flag=TRUE;
  99.                 break;
  100.             }
  101.             code += ((residue-1) * a[j]);
  102.         }
  103.         if(flag)
  104.             continue;
  105.         ++code;
  106.         if(pl[code]!=0)
  107.             tptr[i]=pl[code];
  108.         pl[code]=i;
  109.     }
  110. }
  111.  
  112.  
  113. static void make_n_ptrs(int *tptr,int *pl,int naseq,int len)
  114. {
  115.     static int pot[]={ 0, 1, 4, 16, 64, 256, 1024, 4096 };
  116.     int i,j,limit,code,flag;
  117.     char residue;
  118.     
  119.     limit = (int) pow((double)4,(double)ktup);
  120.     
  121.     for(i=1;i<=limit;++i)
  122.         pl[i]=0;
  123.     for(i=1;i<=len;++i)
  124.         tptr[i]=0;
  125.     
  126.     for(i=1;i<=len-ktup+1;++i) {
  127.         code=0;
  128.         flag=FALSE;
  129.         for(j=1;j<=ktup;++j) {
  130.             residue = seq_array[naseq][i+j-1];
  131.             if(residue<=0) {
  132.                 flag=TRUE;
  133.                 break;
  134.             }
  135.             code += ((residue-1) * pot[j]);
  136.         }
  137.         if(flag)
  138.             continue;
  139.         ++code;
  140.         if(pl[code]!=0)
  141.             tptr[i]=pl[code];
  142.         pl[code]=i;
  143.     }
  144. }
  145.  
  146.  
  147. static void put_frag(int fs,int v1,int v2,int flen)
  148. {
  149.     int end;
  150.     
  151.     accum[0][curr_frag]=fs;
  152.     accum[1][curr_frag]=v1;
  153.     accum[2][curr_frag]=v2;
  154.     accum[3][curr_frag]=flen;
  155.     
  156.     if(!maxsf) {
  157.         maxsf=1;
  158.         accum[4][curr_frag]=0;
  159.         return;
  160.     }
  161.     
  162.         if(fs >= accum[0][maxsf]) {
  163.         accum[4][curr_frag]=maxsf;
  164.         maxsf=curr_frag;
  165.         return;
  166.     }
  167.     else {
  168.         next=maxsf;
  169.         while(TRUE) {
  170.             end=next;
  171.             next=accum[4][next];
  172.             if(fs>=accum[0][next])
  173.                 break;
  174.         }
  175.         accum[4][curr_frag]=next;
  176.         accum[4][end]=curr_frag;
  177.     }
  178. }
  179.  
  180.  
  181. static int frag_rel_pos(int a1,int b1,int a2,int b2)
  182. {
  183.     int ret;
  184.     
  185.     ret=FALSE;
  186.     if(a1-b1==a2-b2) {
  187.         if(a2<a1)
  188.             ret=TRUE;
  189.     }
  190.     else {
  191.         if(a2+ktup-1<a1 && b2+ktup-1<b1)
  192.             ret=TRUE;
  193.     }
  194.     return ret;
  195. }
  196.  
  197.  
  198. static void des_quick_sort(int *array1, int *array2, int array_size)
  199. /*  */
  200. /* Quicksort routine, adapted from chapter 4, page 115 of software tools */
  201. /* by Kernighan and Plauger, (1986) */
  202. /* Sort the elements of array1 and sort the */
  203. /* elements of array2 accordingly */
  204. /*  */
  205. {
  206.     int temp1, temp2;
  207.     int p, pivlin;
  208.     int i, j;
  209.     int lst[50], ust[50];       /* the maximum no. of elements must be*/
  210.                                 /* < log(base2) of 50 */
  211.  
  212.     lst[1] = 1;
  213.     ust[1] = array_size;
  214.     p = 1;
  215.  
  216.     while(p > 0) {
  217.         if(lst[p] >= ust[p])
  218.             p--;
  219.         else {
  220.             i = lst[p] - 1;
  221.             j = ust[p];
  222.             pivlin = array1[j];
  223.             while(i < j) {
  224.                 for(i=i+1; array1[i] < pivlin; i++)
  225.                     ;
  226.                 for(j=j-1; j > i; j--)
  227.                     if(array1[j] <= pivlin) break;
  228.                 if(i < j) {
  229.                     temp1     = array1[i];
  230.                     array1[i] = array1[j];
  231.                     array1[j] = temp1;
  232.                     
  233.                     temp2     = array2[i];
  234.                     array2[i] = array2[j];
  235.                     array2[j] = temp2;
  236.                 }
  237.             }
  238.             
  239.             j = ust[p];
  240.  
  241.             temp1     = array1[i];
  242.             array1[i] = array1[j];
  243.             array1[j] = temp1;
  244.  
  245.             temp2     = array2[i];
  246.             array2[i] = array2[j];
  247.             array2[j] = temp2;
  248.  
  249.             if(i-lst[p] < ust[p] - i) {
  250.                 lst[p+1] = lst[p];
  251.                 ust[p+1] = i - 1;
  252.                 lst[p]   = i + 1;
  253.             }
  254.             else {
  255.                 lst[p+1] = i + 1;
  256.                 ust[p+1] = ust[p];
  257.                 ust[p]   = i - 1;
  258.             }
  259.             p = p + 1;
  260.         }
  261.     }
  262.     return;
  263.  
  264. }
  265.  
  266.  
  267.  
  268.  
  269.  
  270. static void pair_align(int seq_no,int l1,int l2)
  271. {
  272.     int pot[8],i,j,k,l,m,flag,limit,pos,tl1,vn1,vn2,flen,osptr,fs;
  273.     int tv1,tv2,encrypt,subt1,subt2,rmndr;
  274.     char residue;
  275.     
  276.     if(dnaflag) {
  277.         for(i=1;i<=ktup;++i)
  278.             pot[i] = (int) pow((double)4,(double)(i-1));
  279.         limit = (int) pow((double)4,(double)ktup);
  280.     }
  281.     else {
  282.         pot[1]=1;
  283.         pot[2]=20;
  284.         pot[3]=400;
  285.         limit = (int) pow(20.0,(double)ktup);
  286.     }
  287.     
  288.     tl1 = (l1+l2)-1;
  289.     
  290.     for(i=1;i<=tl1;++i) {
  291.         slopes[i]=displ[i]=0;
  292.         diag_index[i] = i;
  293.     }
  294.     
  295.  
  296. /* increment diagonal score for each k_tuple match */
  297.  
  298.     for(i=1;i<=limit;++i) {
  299.         vn1=zzc[i];
  300.         while(TRUE) {
  301.             if(!vn1) break;
  302.             vn2=zzd[i];
  303.             while(vn2 != 0) {
  304.                 osptr=vn1-vn2+l2;
  305.                 ++displ[osptr];
  306.                 vn2=zzb[vn2];
  307.             }
  308.             vn1=zza[vn1];
  309.         }
  310.     }
  311.  
  312. /* choose the top SIGNIF diagonals */
  313.  
  314.     des_quick_sort(displ, diag_index, tl1);
  315.  
  316.     j = tl1 - signif + 1;
  317.     if(j < 1) j = 1;
  318.  
  319. /* flag all diagonals within WINDOW of a top diagonal */
  320.  
  321.     for(i=tl1; i>=j; i--) 
  322.         if(displ[i] > 0) {
  323.             pos = diag_index[i];
  324.             l = (1  >pos-window) ? 1   : pos-window;
  325.             m = (tl1<pos+window) ? tl1 : pos+window;
  326.             for(; l <= m; l++) 
  327.                 slopes[l] = 1;
  328.         }
  329.  
  330.     for(i=1; i<=tl1; i++)  displ[i] = 0;
  331.  
  332.     
  333.     curr_frag=maxsf=0;
  334.     
  335.     for(i=1;i<=(l1-ktup+1);++i) {
  336.         encrypt=flag=0;
  337.         for(j=1;j<=ktup;++j) {
  338.             residue = seq_array[seq_no][i+j-1];
  339.             if(residue<=0) {
  340.                 flag=TRUE;
  341.                 break;
  342.             }
  343.             encrypt += ((residue-1)*pot[j]);
  344.         }
  345.         if(flag) continue;
  346.         ++encrypt;
  347.     
  348.         vn2=zzd[encrypt];
  349.     
  350.         flag=FALSE;
  351.         while(TRUE) {
  352.             if(!vn2) {
  353.                 flag=TRUE;
  354.                 break;
  355.             }
  356.             osptr=i-vn2+l2;
  357.             if(slopes[osptr]!=1) {
  358.                 vn2=zzb[vn2];
  359.                 continue;
  360.             }
  361.             flen=0;
  362.             fs=ktup;
  363.             next=maxsf;
  364.         
  365.         
  366.         /*
  367.         * A-loop
  368.         */
  369.         
  370.             while(TRUE) {
  371.                 if(!next) {
  372.                     ++curr_frag;
  373.                     if(curr_frag>=FSIZE) {
  374.                         fprintf(stdout,"(Partial alignment)");
  375.                         vatend=1;
  376.                         return;
  377.                     }
  378.                     displ[osptr]=curr_frag;
  379.                     put_frag(fs,i,vn2,flen);
  380.                 }
  381.                 else {
  382.                     tv1=accum[1][next];
  383.                     tv2=accum[2][next];
  384.                     if(frag_rel_pos(i,vn2,tv1,tv2)) {
  385.                         if(i-vn2==accum[1][next]-accum[2][next]) {
  386.                             if(i>accum[1][next]+(ktup-1))
  387.                                 fs=accum[0][next]+ktup;
  388.                             else {
  389.                                 rmndr=i-accum[1][next];
  390.                                 fs=accum[0][next]+rmndr;
  391.                             }
  392.                             flen=next;
  393.                             next=0;
  394.                             continue;
  395.                         }
  396.                         else {
  397.                             if(displ[osptr]==0)
  398.                                 subt1=ktup;
  399.                             else {
  400.                                 if(i>accum[1][displ[osptr]]+(ktup-1))
  401.                                     subt1=accum[0][displ[osptr]]+ktup;
  402.                                 else {
  403.                                     rmndr=i-accum[1][displ[osptr]];
  404.                                     subt1=accum[0][displ[osptr]]+rmndr;
  405.                                 }
  406.                             }
  407.                             subt2=accum[0][next]-wind_gap+ktup;
  408.                             if(subt2>subt1) {
  409.                                 flen=next;
  410.                                 fs=subt2;
  411.                             }
  412.                             else {
  413.                                 flen=displ[osptr];
  414.                                 fs=subt1;
  415.                             }
  416.                             next=0;
  417.                             continue;
  418.                         }
  419.                     }
  420.                     else {
  421.                         next=accum[4][next];
  422.                         continue;
  423.                     }
  424.                 }
  425.                 break;
  426.             }
  427.         /*
  428.         * End of Aloop
  429.         */
  430.         
  431.             vn2=zzb[vn2];
  432.         }
  433.     }
  434.     vatend=0;
  435. }         
  436.  
  437. void show_pair()
  438. {
  439.     int i,j,k,dsr;
  440.     double calc_score;
  441.     
  442.     fprintf(stdout,"\n\n");
  443.     
  444.     for(i=1;i<=nseqs;++i) {
  445.         if(dnaflag)
  446.             make_n_ptrs(zza,zzc,i,seqlen_array[i]);
  447.         else
  448.             make_p_ptrs(zza,zzc,i,seqlen_array[i]);
  449.         for(j=i+1;j<=nseqs;++j) {
  450.             if(dnaflag)
  451.                 make_n_ptrs(zzb,zzd,j,seqlen_array[j]);
  452.             else
  453.                 make_p_ptrs(zzb,zzd,j,seqlen_array[j]);
  454.             pair_align(i,seqlen_array[i],seqlen_array[j]);
  455.             if(!maxsf)
  456.                 calc_score=0.0;
  457.             else {
  458.                 calc_score=(double)accum[0][maxsf];
  459.                 if(percent) {
  460.                     dsr=(seqlen_array[i]<seqlen_array[j]) ?
  461.                             seqlen_array[i] : seqlen_array[j];
  462.                     calc_score = (calc_score/(double)dsr) * 100.0;
  463.                 }
  464.             }
  465.             tmat[i][j]=calc_score;
  466.             tmat[j][i]=calc_score;
  467.             if(calc_score>0.1)
  468.                 fprintf(stdout,"Sequences (%d:%d) Aligned. Score: %lg\n",
  469.                i,j,calc_score);
  470.             else
  471.                 fprintf(stdout,"Sequences (%d:%d) Not Aligned\n",i,j);
  472.         }
  473.     }
  474. }
  475.  
  476.