API tools faq
paste
Advertisement
Guest User

al2co

a guest
Sep 30th, 2011
364
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 57.93 KB | None | 0 0
raw download clone embed print report
  1. /*   al2co.c
  2. * ===========================================================================
  3. *
  4. *                            PUBLIC DOMAIN NOTICE
  5. *                   Department of Biochemistry  
  6. *       University of Texas Southwestern Medical Center at Dallas
  7. *
  8. *  This software is freely available to the public for use.We have not placed
  9. *  any restriction on its use or reproduction.
  10. *
  11. *  Although all reasonable efforts have been taken to ensure the accuracy
  12. *  and reliability of the software and data, the University of Texas
  13. *  Southwestern Medical Center does not and cannot warrant the performance
  14. *  or results that may be obtained by using this software or data. The
  15. *  University of Texas Southwestern Medical Center disclaims all warranties,
  16. *  express or implied, including warranties of performance, merchantability
  17. *  or fitness for any particular purpose.
  18. *
  19. *  Please cite the author in any work or product based on this material.
  20. *
  21. * ===========================================================================
  22. *
  23. * File Name:  al2co.c
  24. *
  25. * Author:  Jimin Pei, Nick Grishin
  26.  
  27. * Version Creation Date:  12/23/2000
  28. *
  29. * File Description:
  30. *       A program to calculate positional conservation in a sequence alignment
  31. *
  32. * Last Updated: 06/03/2001 (add -e option)
  33.  
  34. *  Last Updated: 08/13/2002 (MAXSEQNUM)
  35.  
  36.  version 1.2
  37.  Last updated: 05/06/03
  38.  
  39.  change 04.28.2011
  40.  
  41. *
  42. * --------------------------------------------------------------------------*/
  43.  
  44. #include <stdio.h>
  45. #include <stdlib.h>
  46. #include <math.h>
  47. #include <ctype.h>
  48. #include <string.h>
  49. //#include <malloc.h>
  50. #include <stddef.h>
  51.  
  52.  
  53. #define SQUARE(a) ((a)*(a))
  54. #define NUM_METHOD 9
  55. #define MAX_WINDOW 20
  56. #define MAX_DELTASITE 20
  57. #define MAXSTR   10001
  58. #define INDI -100
  59. #define MAXSEQNUM 50000
  60.  
  61. char *digit="0123456789";
  62. void nrerror(char error_text[]);
  63. char *cvector(long nl, long nh);
  64. int *ivector(long nl, long nh);
  65. double *dvector(long nl, long nh);
  66. int **imatrix(long nrl, long nrh, long ncl, long nch);
  67. double **dmatrix(long nrl, long nrh, long ncl, long nch);
  68. double ***d3tensor(long nrl,long nrh,long ncl,long nch,long ndl,long ndh);
  69.  
  70. int a3let2num(char *let);
  71. int am2num_c(int c);
  72. int am2num(int c);
  73. int am2numBZX(int c);
  74.  
  75. static void *mymalloc(int size);
  76. char *strsave(char *str);
  77. char *strnsave(char *str, int l);
  78. static char **incbuf(int n, char **was);
  79. static int *incibuf(int n, int *was);
  80.  
  81. void err_readali(int err_num);
  82. void readali(char *filename);
  83. static void printali(char *argt, int chunk, int n, int len, char **name, char **seq, int *start,int *csv);
  84. int **ali_char2int(char **aseq,int start_num, int start_seq);
  85. int **read_alignment2int(char *filename,int start_num,int start_seq);
  86.  
  87. void counter(int b);
  88. double effective_number(int **ali, int *marks, int n, int start, int end);
  89. double effective_number_nogaps(int **ali, int *marks, int n, int start, int end);
  90. double effective_number_nogaps_expos(int **ali, int *marks, int n, int start, int end, int pos);
  91.  
  92. void **freq(int **ali,double **f,int *num_gaps,int *effindiarr,double gap_threshold);
  93. double *overall_freq(int **ali, int startp, int endp, int *mark);
  94. double *overall_freq_wgt(int **ali,int startp,int endp,int *mark,double *wgt);
  95. double *h_weight(int **ali, int ip);
  96. double **h_freq(int **ali, double **f, double **hfr);
  97. double *entro_conv(double **f, int **ali, double *econv);
  98. double **ic_freq(int **ali, double **f, double **icf);
  99. double *variance_conv(double **f, int **ali, double **oaf, double *vconv);
  100. double *pairs_conv(double **f,int **ali,int **matrix1,int indx,double *pconv);
  101.  
  102.  
  103.  
  104. typedef struct _conv_info{
  105.         double **fq, **hfq, **icfq;
  106.         char *alifilename;
  107.         int alignlen;
  108.     int nali;
  109.     int *ngap;
  110.     int gapless50;
  111.     double eff_num_seq;
  112.     double *over_all_frq;
  113.     int  *eff_indi_arr;
  114.     double *avc,*csi;
  115.         double ***conv;
  116.             } conv_info;
  117. char **aname, **aseq;
  118. int nal, alilen, *astart, *alen;
  119. int **alignment;
  120. double **u_oaf,**h_oaf;
  121. char *am="-WFYMLIVACGPTSNQDEHRKBZX*.wfymlivacgptsnqdehrkbzx";
  122. char *am3[]={
  123. "---",
  124. "TRP",
  125. "PHE",
  126. "TYR",
  127. "MET",
  128. "LEU",
  129. "ILE",
  130. "VAL",
  131. "ALA",
  132. "CYS",
  133. "GLY",
  134. "PRO",
  135. "THR",
  136. "SER",
  137. "ASN",
  138. "GLN",
  139. "ASP",
  140. "GLU",
  141. "HIS",
  142. "ARG",
  143. "LYS",
  144. "ASX",
  145. "GLX",
  146. "UNK",
  147. "***"
  148. "...",
  149. };
  150.  
  151. double *nmlconv(double *conv, conv_info cvf, int wn,int mi);
  152. int  **read_aa_imatrix(FILE *fmat);
  153. int  **identity_imat(long n);
  154. void argument();
  155. void print_parameters(FILE *outfile,char *argi,char *argo,int nt,char *argt,int argb,char *args,int argm,int argf,int argc, int argw, char *argn,char *arga, char *arge, double argg, char *argp,char *argd);
  156.  
  157. main(int argc, char *argv[])
  158. {
  159.     FILE *fout, *fpdb,*matrixfile,*fpdbout,*fp,*ft;
  160.     conv_info convinfo;
  161.     char fstr[MAXSTR+1];
  162.     int i,j,k,l,nt=0;
  163.     double sumofconv,heteroatom;
  164.     int *ptoc1;
  165.     char *ptoc,tmpstr[200],*convstr;
  166.     int fcount=0, fi=0;
  167.     int **smatrix;
  168.     double *consv;
  169.     double **consvall;
  170.     int markali[MAXSEQNUM];
  171.     char ARG_I[500],ARG_O[500],ARG_T[500],ARG_P[500],ARG_D[500],ARG_S[500],ARG_N[500],ARG_A[500],ARG_E[500];
  172.     ARG_I[0]='\0';ARG_O[0]='\0';ARG_T[0]='\0';ARG_P[0]='\0';ARG_D[0]='\0';ARG_S[0]='\0';ARG_N[0]='\0';ARG_A[0]='\0';ARG_E[0]='\0';
  173.     int ARG_F=2,ARG_C=0,ARG_W=1,ARG_M=0,ARG_B=60;
  174.     double ARG_G=0.5;
  175.     double max_index, mini_index;
  176.     int *csv_index;
  177.     int NAL;
  178.     int **ALIGNMENT;
  179.  
  180.  
  181.     /*read input arguments */
  182.         if(argc<=2) { argument(); exit(0);}
  183.     for(i=1;i<argc;i++) {
  184.         if(strcmp(argv[i],"-i")==0) {strcpy(ARG_I,argv[i+1]);i++;continue;}
  185.         if(strcmp(argv[i],"-o")==0) {strcpy(ARG_O,argv[i+1]);i++;continue;}
  186.         if(strcmp(argv[i],"-t")==0) {strcpy(ARG_T,argv[i+1]);i++;continue;}
  187.         if(strcmp(argv[i],"-p")==0) {strcpy(ARG_P,argv[i+1]);i++;continue;}
  188.         if(strcmp(argv[i],"-d")==0) {strcpy(ARG_D,argv[i+1]);i++;continue;}
  189.         if(strcmp(argv[i],"-s")==0) {strcpy(ARG_S,argv[i+1]);i++;continue;}
  190.         if(strcmp(argv[i],"-n")==0) {strcpy(ARG_N,argv[i+1]);i++;continue;}
  191.         if(strcmp(argv[i],"-a")==0) {strcpy(ARG_A,argv[i+1]);i++;continue;}
  192.         if(strcmp(argv[i],"-e")==0) {strcpy(ARG_E,argv[i+1]);i++;continue;}
  193.         if(strcmp(argv[i],"-f")==0) {sscanf(argv[i+1],"%d",&ARG_F);i++;continue;}
  194.         if(strcmp(argv[i],"-b")==0) {sscanf(argv[i+1],"%d",&ARG_B);i++;continue;}
  195.         if(strcmp(argv[i],"-c")==0) {sscanf(argv[i+1],"%d",&ARG_C);i++;continue;}
  196.         if(strcmp(argv[i],"-w")==0) {sscanf(argv[i+1],"%d",&ARG_W);i++;continue;}
  197.         if(strcmp(argv[i],"-m")==0) {sscanf(argv[i+1],"%d",&ARG_M);i++;continue;}
  198.         if(strcmp(argv[i],"-g")==0) {sscanf(argv[i+1],"%lf",&ARG_G);i++;continue;}
  199.                 }
  200.    
  201.         if((ARG_F>2)||(ARG_F<0)){fprintf(stderr,"frequency calculation method(-f): \n0, unweighted; 1, Henikoff weight; 2, independent count\n");
  202.                     exit(0);}
  203.         if((ARG_C>2)||(ARG_C<0)){fprintf(stderr,"conservation calculation strategy(-c):\n0,entropy;1,variance;2,sumofpairs\n");
  204.                     exit(0);}
  205.         if((ARG_M>2)||(ARG_M<0)){fprintf(stderr,"matrix transform(-m):\n0, no transform;1,normalization;2,adjustment\n");
  206.                     exit(0);}
  207.     if((ARG_G>1.0)||(ARG_G<=0)){fprintf(stderr,"gap percentage(-g) to suppress calculation must be no more than 1 and more than 0 \n");
  208.             exit(0);}
  209.    
  210.     /* smatrix:  identity matrix if reading from input failed */
  211.         if((matrixfile=fopen(ARG_S,"r"))==NULL){
  212.                 if(strlen(ARG_S)!=0)fprintf(stderr, "Warning: not readable  matrixfile: %s; default: identity matrix\n",ARG_S);
  213.                 smatrix = identity_imat(24);   }
  214.         else {smatrix=read_aa_imatrix(matrixfile);
  215.           /*for(i=1;i<=20;i++){
  216.         for(j=1;j<=20;j++){
  217.             fprintf(stderr, "%d ",smatrix[i][j]);
  218.            
  219.                   }
  220.         fprintf(stderr,"\n");
  221.                 }*/
  222.            
  223.           fclose(matrixfile); }
  224.  
  225.     /* read alignment */
  226.     ALIGNMENT = read_alignment2int(ARG_I,1,1);
  227.     if( (strcmp(ARG_E, "T") == 0) || (strcmp(ARG_E, "t")==0) ) {
  228.         alignment = read_alignment2int(ARG_I,1,1);
  229.         alignment = alignment + 1;
  230.         NAL = nal;
  231.         nal = nal -1;
  232.     }
  233.     else {
  234.         alignment=read_alignment2int(ARG_I,1,1);
  235.         NAL = nal;
  236.     }
  237.  
  238.     if(nal>=MAXSEQNUM) {
  239.         fprintf(stderr, "Error: Number of sequences exceeds %d\n", MAXSEQNUM);
  240.         exit(0);
  241.     }
  242.     if(alignment==NULL){
  243.         fprintf(stderr, "alignment file not readable\n");
  244.                }
  245.  
  246.     /* memory allocation for the elements in convinfo*/
  247.     convinfo.alignlen=alilen;
  248.     convinfo.nali=nal;
  249.     convinfo.ngap=ivector(0,alilen);
  250.     convinfo.eff_indi_arr=ivector(0,alilen+1);
  251.     convinfo.alifilename=cvector(0,strlen(ARG_I));
  252.     strcpy(convinfo.alifilename,ARG_I);
  253.     convinfo.over_all_frq=dvector(0,20);
  254.     convinfo.fq=dmatrix(0,20,0,alilen);
  255.     convinfo.hfq=dmatrix(0,20,0,alilen);
  256.     convinfo.icfq=dmatrix(0,20,0,alilen);
  257.     convinfo.conv=d3tensor(0,MAX_WINDOW,0,NUM_METHOD,0,alilen);
  258.     convinfo.avc=dvector(0,NUM_METHOD);
  259.     convinfo.csi=dvector(0,NUM_METHOD);
  260.     consv=dvector(0,alilen);
  261.     consvall=dmatrix(0,NUM_METHOD,0,alilen);/* for all 9 methods */
  262.     csv_index=ivector(0,alilen);
  263.     for(i=1;i<=MAX_WINDOW;i++)
  264.         for(j=1;j<=NUM_METHOD;j++)
  265.             for(k=0;k<=alilen;k++)
  266.                 convinfo.conv[i][j][k]=INDI;
  267.  
  268.     /* get the frequencies */
  269.     freq(alignment,convinfo.fq,convinfo.ngap,convinfo.eff_indi_arr,ARG_G);
  270.     h_freq(alignment,convinfo.fq,convinfo.hfq);
  271.     ic_freq(alignment,convinfo.fq,convinfo.icfq);
  272.  
  273.     convinfo.gapless50=0;
  274.     for(i=1;i<=alilen;i++){
  275.         if(convinfo.ngap[i]<=0.5*convinfo.nali)
  276.         convinfo.gapless50++;
  277.                   }
  278.     for(i=1;i<=nal;i++){markali[i]=1;}
  279.     convinfo.over_all_frq=overall_freq(alignment,1,alilen,markali);
  280.  
  281.  
  282.     /* get the conservation indeces for method[1..NUM_METHOD] */
  283.     entro_conv(convinfo.fq,alignment,convinfo.conv[1][1]); 
  284.     entro_conv(convinfo.hfq,alignment,convinfo.conv[1][2]);
  285.     entro_conv(convinfo.icfq,alignment,convinfo.conv[1][3]);
  286.     variance_conv(convinfo.fq,alignment,u_oaf,convinfo.conv[1][4]);
  287.     variance_conv(convinfo.hfq,alignment,h_oaf,convinfo.conv[1][5]);
  288.     variance_conv(convinfo.icfq,alignment,u_oaf,convinfo.conv[1][6]);
  289.     pairs_conv(convinfo.fq,alignment,smatrix,ARG_M,convinfo.conv[1][7]);
  290.     pairs_conv(convinfo.hfq,alignment,smatrix,ARG_M,convinfo.conv[1][8]);
  291.     pairs_conv(convinfo.icfq,alignment,smatrix,ARG_M,convinfo.conv[1][9]);
  292.  
  293.     /* average over window size ARG_W */
  294.     for(i=1;i<=9;i++){
  295.         sumofconv=0;
  296.         k=0;
  297.         for(j=1;j<=alilen;j++){
  298.            if(convinfo.conv[1][i][j]==INDI)continue;
  299.            sumofconv+=convinfo.conv[1][i][j];
  300.            k++;
  301.                     }
  302.         convinfo.avc[i]=sumofconv/k;
  303.              }
  304.     if(ARG_W>1){
  305.         ptoc1=convinfo.eff_indi_arr;
  306.         for(i=1;i<=NUM_METHOD;i++){
  307.         j=ARG_W;
  308.                 if(j>*ptoc1){fprintf(stderr,"window size too big!\nARG_W:%d;   ptoc1:%d\n",j,*ptoc1);exit(0);}
  309.                 for(k=(j+1)/2;k<=convinfo.eff_indi_arr[0]-j/2;k++){
  310.                         sumofconv=0;
  311.                         for(l=k-(j-1)/2;l<=k+j/2;l++){
  312.                                 sumofconv+=convinfo.conv[1][i][*(ptoc1+l)];
  313.                                                      }
  314.                         convinfo.conv[j][i][*(ptoc1+k)]=sumofconv/j;
  315.                                                                      }
  316.         for(k=1;k<(j+1)/2;k++){
  317.             sumofconv=0;
  318.             for(l=1;l<=2*k-1;l++){
  319.                 sumofconv+=convinfo.conv[1][i][*(ptoc1+l)];
  320.                          }
  321.             convinfo.conv[j][i][*(ptoc1+k)]=convinfo.avc[i]+(sumofconv/(2*k-1)-convinfo.avc[i])*sqrt(1.0*(2*k-1)/ARG_W);
  322.                       }
  323.         for(k=0;k<=j/2-1;k++){
  324.             sumofconv=0;
  325.             for(l=1;l<=2*k+1;l++){
  326.                 sumofconv+=convinfo.conv[1][i][*(ptoc1+convinfo.eff_indi_arr[0]-l+1)];
  327.                          }
  328.             convinfo.conv[j][i][*(ptoc1+convinfo.eff_indi_arr[0]-k)]=convinfo.avc[i]+(sumofconv/(2*k+1)-convinfo.avc[i])*sqrt(1.0*(2*k+1)/ARG_W);
  329.                    }
  330.                                   }
  331.             }
  332.  
  333.  
  334.     /* normalize the conservation indices */
  335.     nmlconv(consv, convinfo, ARG_W, ARG_C*3+ARG_F+1);
  336.     if((strcmp(ARG_N,"F")==0)||(strcmp(ARG_N,"f")==0)){
  337.         for(i=1;i<=convinfo.alignlen;i++){
  338.         if(convinfo.conv[ARG_W][ARG_C*3+ARG_F+1][i]==INDI){
  339.             consv[i]=0-convinfo.csi[ARG_C*3+ARG_F+1]+convinfo.avc[ARG_C*3+ARG_F+1];}
  340.         else consv[i]=convinfo.conv[ARG_W][ARG_C*3+ARG_F+1][i];
  341.                       }
  342.                             }
  343.  
  344.     /* find the largest and the smallest conservation indices */
  345.     max_index=mini_index=consv[1];
  346.     for(i=2;i<=alilen;i++) {
  347.         if(consv[i]>max_index) max_index=consv[i];
  348.         if(consv[i]<mini_index) mini_index=consv[i];
  349.                 }
  350.     for(i=1;i<=alilen;i++) {
  351.         csv_index[i]=(int)(9.99*(consv[i]-mini_index)/(max_index-mini_index));
  352.                 }
  353.  
  354.         /* print file with conservation indices mapped to the alignment */
  355.         if((ft=fopen(ARG_T,"w"))==NULL){;}
  356.         else{
  357.                 nt=1;
  358.                 fclose(ft);
  359.                 if(ARG_B){;}
  360.                 else{ARG_B=60;}
  361.                 /*printali(ARG_T, ARG_B, nal, alilen, aname, aseq, astart, csv_index);*/
  362.         printali(ARG_T, ARG_B, NAL, alilen, aname, aseq, astart, csv_index);
  363.             }
  364.  
  365.     /* change b-factor to conservation index */
  366. if((fpdb=fopen(ARG_P,"r"))==NULL){
  367.     if(strlen(ARG_P)!=0){fprintf(stderr,"please give the right pdb file\n");
  368.                 }
  369.                  }
  370. else{
  371.         if((fpdbout=fopen(ARG_D,"w"))==NULL){
  372.                 if(strlen(ARG_D)!=0){fprintf(stderr,"output pdb file using default name\n");}
  373.                 fpdbout=stdout;}
  374.  
  375.     convstr=cvector(0,200);
  376.     while(fgets(fstr,MAXSTR,fpdb)!=NULL){
  377.         if(strncmp(fstr,"ATOM   ",7)==0)break;
  378.                         }
  379.     if(feof(fpdb)){fprintf(stderr,"end of the file reached: not a good pdb file");exit(0);}
  380.  
  381.     ptoc=fstr+23;
  382.     j=atoi(ptoc);i=1;
  383.     for(i=1;ALIGNMENT[1][i]==0;i++);
  384.     ptoc=fstr+17;
  385.     if(strncmp(ptoc,am3[ALIGNMENT[1][i]],3)!=0){
  386.         fprintf(stderr,"The first sequence does not match the pdb file\n");
  387.         exit(0);            }
  388.     while(!feof(fpdb)){
  389.         ptoc=fstr+66;
  390.         strcpy(tmpstr, ptoc);
  391.         *(fstr+61)='\0';
  392.         /*if(strncmp(fstr,"HETATM ",7)==0){
  393.             heteroatom=-1.0;
  394.             sprintf(convstr,"%5.2f",heteroatom);
  395.             strcat(fstr,convstr);
  396.             strcat(fstr,tmpstr);
  397.             fprintf(fpdbout,"%s",fstr);}
  398.         else{*/
  399.         if(consv[i]>-10){sprintf(convstr,"%5.2f",consv[i]);}
  400.         else{consv[i]=-9.9;sprintf(convstr,"%5.2f",consv[i]);}
  401.         strcat(fstr,convstr);
  402.         strcat(fstr,tmpstr);
  403.         fprintf(fpdbout,"%s",fstr);
  404.             /*}*/
  405.  
  406.         fgets(fstr,MAXSTR,fpdb);
  407.         if(strncmp(fstr,"ATOM   ",7)!=0){
  408.            if(strncmp(fstr,"HETATM ",7)==0){;}
  409.            else{fprintf(fpdbout,"END\n");break;}
  410.                         }
  411.         ptoc=fstr+23;
  412.         if(j!=atoi(ptoc)) {
  413.             /*fprintf(stderr,"i=%d  ptoc=%d\n", i, atoi(ptoc));*/
  414.  
  415.             /* change 04.28.2011 */
  416.             /*i=i+atoi(ptoc)-j;j=atoi(ptoc);*/
  417.             i++;j=atoi(ptoc);
  418.  
  419.             for(;ALIGNMENT[1][i]==0;i++);
  420.         }
  421.         ptoc=fstr+17;
  422.         if(strncmp(ptoc,am3[ALIGNMENT[1][i]],3)!=0){
  423.             if(strncmp(fstr,"HETATM ",7)==0){;}
  424.             if(i>alilen) {fprintf(fpdbout,"END\n");break;}
  425.             if(strncmp(ptoc,am3[ALIGNMENT[1][i]],3)!=0){
  426.                if(strncmp(fstr,"HETATM ",7)!=0){
  427.                 fprintf(stderr,"The first sequence does not match the pdb file\n");
  428.                 fprintf(stderr,"%s %s", am3[ALIGNMENT[1][i]], ptoc);
  429.                 exit(0);
  430.                                }
  431.                                   }
  432.                             }
  433.                }   
  434.     fclose(fpdb);fclose(fpdbout);
  435.     }
  436.  
  437.         /* print conservation indices */
  438.     if( (strcmp(ARG_A,"T")==0)||(strcmp(ARG_A,"t")==0) ) {
  439.     /* print results from all the 9 methods */
  440.     for(j=1;j<=NUM_METHOD;j++) {
  441.            nmlconv(consvall[j], convinfo, ARG_W, j);
  442.            if((strcmp(ARG_N,"F")==0)||(strcmp(ARG_N,"f")==0)){
  443.             for(i=1;i<=convinfo.alignlen;i++){
  444.                 if(convinfo.conv[ARG_W][j][i]==INDI)
  445.                         consvall[j][i]=0-convinfo.csi[j]+convinfo.avc[j];
  446.                 else consvall[j][i]=convinfo.conv[ARG_W][j][i];
  447.             }                              
  448.            }
  449.     }
  450.  
  451.     if((fout=fopen(ARG_O,"w"))==NULL) fout = stdout;
  452.         for(i=1;i<=convinfo.alignlen;i++){
  453.         fprintf(fout, "%-5d %c      ",i,am[alignment[1][i]]);
  454.         for(j=1;j<=NUM_METHOD;j++) {
  455.         fprintf(fout, "%8.3f", consvall[j][i]);
  456.         }
  457.         if(convinfo.conv[ARG_W][1][i]==INDI){
  458.         fprintf(fout,"     *\n");
  459.         } else fprintf(fout, "\n");
  460.         }
  461.         fprintf(fout,"* gap fraction no less than %5.2f; conservation set to M-S\n", ARG_G);
  462.         fprintf(fout,"  M: mean;  S: standard deviation\n");
  463.     print_parameters(fout,ARG_I,ARG_O,nt,ARG_T,ARG_B,ARG_S,ARG_M,ARG_F,ARG_C,ARG_W,ARG_N,ARG_A,ARG_E,ARG_G,ARG_P,ARG_D);
  464.     fclose(fout);
  465.  
  466.     }
  467.     else {
  468.         if((fout=fopen(ARG_O,"w"))==NULL){
  469.         for(i=1;i<=convinfo.alignlen;i++){
  470.                 if(convinfo.conv[ARG_W][ARG_C*3+ARG_F+1][i]==INDI){
  471.                 fprintf(stdout, "%-5d %c      %6.3f    *\n", i,am[alignment[1][i]], consv[i]);
  472.                                                                   }
  473.                 else{
  474.                 fprintf(stdout, "%-5d %c      %6.3f\n", i,am[alignment[1][i]], consv[i]);
  475.                     }
  476.                                          }
  477.  
  478.         fprintf(stdout,"* gap fraction no less than %5.2f; conservation set to M-S\n", ARG_G);
  479.     fprintf(stdout,"  M: mean;  S: standard deviation\n");
  480.         print_parameters(stdout,ARG_I,ARG_O,nt,ARG_T,ARG_B,ARG_S,ARG_M,ARG_F,ARG_C,ARG_W,ARG_N,ARG_A,ARG_E,ARG_G,ARG_P,ARG_D);
  481.                                          }
  482.         else{
  483.         for(i=1;i<=convinfo.alignlen;i++){
  484.                 if(convinfo.conv[ARG_W][ARG_C*3+ARG_F+1][i]==INDI){
  485.                 fprintf(fout, "%-5d %c      %6.3f *\n", i,am[alignment[1][i]], consv[i]);
  486.                                                                   }
  487.                 else{
  488.                 fprintf(fout, "%-5d %c      %6.3f\n", i,am[alignment[1][i]], consv[i]);
  489.                     }
  490.                                          }
  491.         fprintf(fout,"* gap fraction no less than %5.2f; conservation set to M-S\n", ARG_G);
  492.     fprintf(fout,"  M: mean;  S: standard deviation\n");
  493.         print_parameters(fout,ARG_I,ARG_O,nt,ARG_T,ARG_B,ARG_S,ARG_M,ARG_F,ARG_C,ARG_W,ARG_N,ARG_A,ARG_E,ARG_G,ARG_P,ARG_D);
  494.         fclose(fout);
  495.             }
  496.     }
  497.     /*print_parameters(stderr,ARG_I,ARG_O,nt,ARG_T,ARG_B,ARG_S,ARG_M,ARG_F,ARG_C,ARG_W,ARG_N,ARG_G,ARG_P,ARG_D);*/
  498.     exit(0);
  499.  
  500. }          
  501.  
  502. double *nmlconv(double *conv, conv_info cvf, int wn,int mi)
  503. {
  504.         double mean=0,vc=0;
  505.         int i,cps=0;
  506.  
  507.         for(i=1;i<=cvf.alignlen;i++){
  508.                 if(cvf.conv[wn][mi][i]==INDI) continue;
  509.                 cps++;
  510.                 mean+=cvf.conv[wn][mi][i];
  511.                               }
  512.         if(cps<=1) {fprintf(stderr,"less than one position\n");exit(0);}
  513.         mean=mean/cps;
  514.     /*fprintf(stderr, "mean value is: %f\n",mean);*/
  515.         cvf.avc[mi]=mean;
  516.         for(i=1;i<=cvf.alignlen;i++){
  517.                 if(cvf.conv[wn][mi][i]==INDI) continue;
  518.                 vc+=(cvf.conv[wn][mi][i]-mean)*(cvf.conv[wn][mi][i]-mean);
  519.                               }
  520.         vc=sqrt(vc/(cps-1));
  521.     /*fprintf(stderr, "sigma is: %f\n",vc);*/
  522.         cvf.csi[mi]=vc;
  523.  
  524.         if(vc==0){fprintf(stderr,"variance is zero\n");exit(0);}
  525.         for(i=1;i<=cvf.alignlen;i++){
  526.                 if(cvf.conv[wn][mi][i]==INDI) {conv[i]=-1.0;continue;}
  527.                 if(vc!=0)conv[i]=(cvf.conv[wn][mi][i]-mean)/vc;
  528.                 else conv[i]=0;
  529.                               }
  530. }
  531.  
  532. int  **read_aa_imatrix(FILE *fmat){
  533.  
  534. /* read matrix from file *fmat */
  535.  
  536. int i,ri,rj,c,flag,j,t;
  537. int col[31],row[31];
  538. char stri[11];
  539. int **mat;
  540.  
  541. mat=imatrix(0,25,0,25);
  542. for(i=0;i<=25;++i)for(j=0;j<=25;++j)mat[i][j]=0;
  543.  
  544. i=0;
  545. ri=0;
  546. rj=0;
  547.  
  548. flag=0;
  549.  
  550.  
  551. while( (c=getc(fmat)) != EOF){
  552.  
  553. if(flag==0 && c=='#'){flag=-1;continue;}
  554. else if(flag==-1 && c=='\n'){flag=0;continue;}
  555. else if(flag==-1){continue;}
  556. else if(flag==0 && c==' '){flag=1;continue;}
  557. else if(flag==1 && c=='\n'){flag=0;continue;}
  558. else if(flag==1 && c==' '){continue;}
  559. else if(flag==1){
  560.                 ++i;
  561.                 if(i>=25){nrerror("matrix has more than 25 columns: FATAL");exit(0);}
  562.                 col[i]=am2numBZX(c);
  563.                 continue;
  564.                 }
  565. else if(flag==0 && c!=' ' && c!='#'){
  566.                 ri=0;
  567.                 ++rj;
  568.                 if(rj>=25){nrerror("matrix has more than 25 rows: FATAL");exit(0);}
  569.                 row[rj]=am2numBZX(c);
  570.                 flag=2;
  571.                 continue;
  572.                 }
  573. else if (flag==2 && c==' '){for(i=0;i<=10;++i){stri[i]=' ';}j=0;continue;}
  574. else if (flag==2 && c=='\n'){flag=0;continue;}
  575. else if (flag==2){flag=3;stri[++j]=c;if(j>10){nrerror("string too long:FATAL");exit(0);}continue;}
  576. else if (flag==3 && c==' ' || flag==3 && c=='\n'){
  577.                         j=0;
  578.                         ++ri;
  579.                         t=atoi(stri);
  580.                         mat[row[rj]][col[ri]]=t;
  581.                         if (c=='\n')flag=0;else flag=2;
  582.                         continue;
  583.                         }
  584. else if (flag==3){stri[++j]=c;continue;}
  585.  
  586. }
  587.  
  588. return mat;
  589. }
  590.  
  591. int    **identity_imat(long n){
  592. /* allocates square integer identity matrix of length n+1: m[0.n][0.n] */
  593.  
  594. int i,j;
  595. int **m;
  596. m=imatrix(0,n,0,n);
  597. for(i=0;i<=n;++i)for(j=0;j<=n;++j){if(i==j)m[i][j]=1;else m[i][j]=0;}
  598. return m;
  599. }
  600.  
  601. void argument()
  602. {
  603. fprintf(stderr,"      al2co   arguments:\n");
  604. fprintf(stderr,"\n");
  605. fprintf(stderr,"  -i    Input alignment file [File in]\n");
  606. fprintf(stderr,"        Format: ClustalW or simple alignment format\n");
  607. fprintf(stderr,"        The title (first line) should begin with \"CLUSTAL\", or\n");
  608. fprintf(stderr,"        the title line should be deleted.\n");
  609. fprintf(stderr,"\n");
  610. fprintf(stderr,"  -o    Output file with conservation index for each position in the\n");
  611. fprintf(stderr,"        alignment [File out] Optional\n");
  612. fprintf(stderr,"        Default = STDOUT\n");
  613. fprintf(stderr,"\n");
  614. fprintf(stderr,"  -t    Output file with conservation index mapped to the alignment\n");
  615. fprintf(stderr,"        [File out] Optional\n");
  616. fprintf(stderr,"        Conservation indices are linearly rescaled to be from 0\n");
  617. fprintf(stderr,"        to 9.99. C'=9.99*(C-MIN)/(MAX-MIN), where C and C' are the\n");
  618. fprintf(stderr,"        the indices before and after rescaling respectively, MAX and\n");
  619. fprintf(stderr,"        MIN are the highest index and lowest index before rescaling\n");
  620. fprintf(stderr,"        respectively. The integer part of each rescaled index is\n");
  621. fprintf(stderr,"        written out along with the sequence alignment.\n");
  622. fprintf(stderr,"        Default = no output\n");
  623. fprintf(stderr,"\n");
  624. fprintf(stderr,"  -b    Block size of the output alignment file with conservation\n");
  625. fprintf(stderr,"        [Integer] Optional\n");
  626. fprintf(stderr,"        Default = 60\n");
  627. fprintf(stderr,"\n");
  628. fprintf(stderr,"  -s    Input file with the scoring matrix [File in] Optional\n");
  629. fprintf(stderr,"        Format: NCBI\n");
  630. fprintf(stderr,"        Notice: Scoring matrix is only used for sum-of-pairs measure\n");
  631. fprintf(stderr,"        with option -c  2.\n");
  632. fprintf(stderr,"        Default = identity matrix\n");
  633. fprintf(stderr,"\n");
  634. fprintf(stderr,"  -m    Scoring matrix transformation [Integer] Optional\n");
  635. fprintf(stderr,"        Options:\n");
  636. fprintf(stderr,"        0=no transformation,\n");
  637. fprintf(stderr,"        1=normalization S'(a,b)=S(a,b)/sqrt[S(a,a)*S(b,b)],\n");
  638. fprintf(stderr,"        2=adjustment S\"(a,b)=2*S(a,b)-(S(a,a)+S(b,b))/2\n");
  639. fprintf(stderr,"        Default = 0\n");
  640. fprintf(stderr,"\n");
  641. fprintf(stderr,"  -f    Weighting scheme for amino acid frequency estimation [Integer] Optional\n");
  642. fprintf(stderr,"        Options:\n");
  643. fprintf(stderr,"        0=unweighted,\n");
  644. fprintf(stderr,"        1=weighted by the modified method of Henikoff & Henikoff (2,3),\n");
  645. fprintf(stderr,"        2=independent-count based (1)\n");
  646. fprintf(stderr,"        Default = 2\n");
  647. fprintf(stderr,"\n");
  648. fprintf(stderr,"  -c    Conservation calculation method [Integer] Optional\n");
  649. fprintf(stderr,"        Options:\n");
  650. fprintf(stderr,"        0=entropy-based    C(i)=sum_{a=1}^{20}f_a(i)*ln[f_a(i)], where f_a(i)\n");
  651. fprintf(stderr,"          is the frequency of amino acid a at position i,\n");
  652. fprintf(stderr,"        1=variance-based   C(i)=sqrt[sum_{a=1}^{20}(f_a(i)-f_a)^2], where f_a\n");
  653. fprintf(stderr,"          is the overall frequency of amino acid a,\n");
  654. fprintf(stderr,"        2=sum-of-pairs measure   C(i)=sum_{a=1}^{20}sum_{b=1}^{20}f_a(i)*f_b(i)*S_{ab},\n");
  655. fprintf(stderr,"          where S_{ab} is the element of a scoring matrix for amino acids a and b\n");
  656. fprintf(stderr,"        Default = 0\n");
  657. fprintf(stderr,"\n");
  658. fprintf(stderr,"  -w    Window size used for averaging [Integer] Optional\n");
  659. fprintf(stderr,"        Default = 1\n");
  660. fprintf(stderr,"        Recommended value for motif analysis: 3\n");
  661. fprintf(stderr,"\n");
  662. fprintf(stderr,"  -n    Normalization option [T/F] Optional\n");
  663. fprintf(stderr,"        Subtract the mean from each conservation index and divide by the\n");
  664. fprintf(stderr,"        standard deviation.\n");
  665. fprintf(stderr,"        Default = T\n");
  666. fprintf(stderr,"\n");
  667. fprintf(stderr,"  -a    All methods option [T/F] Optional\n");
  668. fprintf(stderr,"        If set to true, the results of all 9 methods will be output.\n");
  669. fprintf(stderr,"        1. unweighted entropy measure; 2. Henikoff entropy measure;\n");
  670. fprintf(stderr,"        3. independent count entropy measure;\n");
  671. fprintf(stderr,"        4. unweighted variance measure; 5. Henikoff variance measure;\n");
  672. fprintf(stderr,"        6. independent count variance measure;\n");
  673. fprintf(stderr,"        7. unweighted matrix-based sum-of-pairs measure;\n");
  674. fprintf(stderr,"        8. Henikoff matrix-based sum-of-pairs measure;\n");
  675. fprintf(stderr,"        9. independent count matrix-based sum-of-pairs measure;\n");
  676. fprintf(stderr,"        Default = F\n");
  677. fprintf(stderr,"\n");
  678. fprintf(stderr,"  -e    Excluding the first sequence from calculation [T/F] Optional\n");
  679. fprintf(stderr,"        If set to true, the first sequence in the alignment will not\n");
  680. fprintf(stderr,"    be included in the conservation calculation.\n");
  681. fprintf(stderr,"        Default = F\n\n");
  682. fprintf(stderr,"  -g    Gap fraction to suppress conservation calculation [Real] Optional\n");
  683. fprintf(stderr,"        The value should be more than 0 and no more than 1. Conservation\n");
  684. fprintf(stderr,"        indices are calculated only for positions with gap fraction less\n");
  685. fprintf(stderr,"        than the specified value. Otherwise, conservation indices will\n");
  686. fprintf(stderr,"        be set to M-S, where M is the mean conservation value and S is\n");
  687. fprintf(stderr,"        the standard deviation.\n");
  688. fprintf(stderr,"        Default = 0.5\n");
  689. fprintf(stderr,"\n");
  690. fprintf(stderr,"  -p    Input pdb file [File in] Optional\n");
  691. fprintf(stderr,"        The sequence in the pdb file should match exactly the first sequence of\n");
  692. fprintf(stderr,"        the alignment.\n");
  693. fprintf(stderr,"\n");
  694. fprintf(stderr,"  -d    Output pdb file [File Out] Optional\n");
  695. fprintf(stderr,"        The B-factors are replaced by the conservation indices.\n");
  696. fprintf(stderr,"        Default = STDOUT\n");
  697. fprintf(stderr,"\n");
  698.  
  699.  
  700. }
  701.  
  702. void print_parameters(FILE *outfile,char *argi,char *argo,int nt,char *argt,int argb,char *args,int argm,int argf,int argc, int argw, char *argn,char *arga, char *arge, double argg, char *argp,char *argd)
  703. {
  704.     FILE *fp;
  705.     char *mt[]={
  706.         "no transformation",
  707.         "normalization S'(a,b)=S(a,b)/sqrt[S(a,a)*S(b,b)]",
  708.         "adjustment S\"(a,b)=2*S(a,b)-(S(a,a)+S(b,b))/2",
  709.            };
  710.     char *ws[]={
  711.         "unweighted",
  712.         "weighted by the modified method of Henikoff & Henikoff",
  713.         "independent-count based",
  714.            };
  715.     char *cm[]={
  716.         "entropy-based",
  717.         "variance-based",
  718.         "sum-of-pairs measure",
  719.            };
  720.  
  721.     fprintf(outfile, "\nal2co - The parameters are:\n\n");
  722.     fprintf(outfile, "Input alignment file - %s\n", argi);
  723.  
  724.     if((fp=fopen(argo,"r"))==NULL){
  725.         fprintf(outfile, "Output conservation - STDOUT\n");
  726.                      }
  727.     else{   fprintf(outfile, "Output conservation file - %s\n", argo); fclose(fp);}
  728.  
  729.     if(nt==0){;}
  730.     else if(nt==1)  {  
  731.         fprintf(outfile, "Output alignment file with index - %s;", argt);
  732.         fprintf(outfile, " Block size - %d\n",argb);
  733.             }
  734.  
  735.     if((fp=fopen(args,"r"))==NULL){;}
  736.     else {
  737.         if(argc==2) {
  738.             fprintf(outfile,"Input matrix file - %s\n",args);
  739.             fprintf(outfile,"Matrix transformation -  %s\n",mt[argm]);
  740.                 }
  741.         else {
  742.             fprintf(outfile,"Input matrix file - %s\n",args);
  743.             if( (strcmp(arga,"T")==0) || (strcmp(arga, "t")==0) ) {;}
  744.             else fprintf(outfile,"Warning - Matrix not used: matrix is for sum-of-pairs measure, -c 2\n"); }
  745.         fclose(fp);
  746.           }
  747.  
  748.     if( (strcmp(arga,"T")==0) || (strcmp(arga, "t")==0) ) {
  749.          fprintf(outfile, "All 9 methods are used - true \n");
  750.          fprintf(outfile, "   1. unweighted entropy measure; 2. Henikoff entropy measure;\n");
  751.          fprintf(outfile, "   3. independent count entropy measure;\n");
  752.          fprintf(outfile, "   4. unweighted variance measure; 5. Henikoff variance measure;\n");
  753.          fprintf(outfile, "   6. independent count variance measure;\n");
  754.          fprintf(outfile, "   7. unweighted matrix-based sum-of-pairs measure;\n");
  755.          fprintf(outfile, "   8. Henikoff matrix-based sum-of-pairs measure;\n");
  756.          fprintf(outfile, "   9. independent count matrix-based sum-of-pairs measure;\n");
  757.     }
  758.     else {
  759.       fprintf(outfile, "Weighting scheme - %s\n", ws[argf]);
  760.       fprintf(outfile, "Conservation calculation method - %s\n",cm[argc]);
  761.     }
  762.  
  763.     if( (strcmp(arge, "T")==0) || (strcmp(arge, "t")==0 ) ) {
  764.         fprintf(outfile, "The first sequence is not used in calculation\n");
  765.     }
  766.  
  767.     fprintf(outfile, "Window size - %d\n", argw);
  768.      
  769.     if((strcmp(argn,"f")==0)||(strcmp(argn,"F")==0)){
  770.         fprintf(outfile, "Conservation not normalized\n");}
  771.     else fprintf(outfile, "Conservation normalized to zero mean and unity variance\n");
  772.  
  773.     fprintf(outfile, "Gap fraction to suppress calculation - %5.2f\n",argg);
  774.     if((fp=fopen(argp,"r"))==NULL){;}
  775.     else{
  776.         fprintf(outfile,"Input pdb file - %s\n", argp);
  777.         fclose(fp);
  778.         if((fp=fopen(argd,"r"))==NULL){
  779.             fprintf(outfile, "Output pdb file - STDOUT\n");
  780.                         }
  781.         else {  fprintf(outfile, "Output pdb file - %s\n", argd);
  782.             fclose(fp);
  783.             }
  784.         }
  785.    
  786. }
  787. #define NR_END 1
  788. #define FREE_ARG char*
  789.  
  790.  
  791. void nrerror(char error_text[]){
  792. fprintf(stderr,"%s\n",error_text);
  793. fprintf(stderr,"FATAL - execution terminated\n");
  794. exit(1);
  795. }
  796.  
  797.  
  798. char *cvector(long nl, long nh){
  799. char *v;
  800. v=(char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
  801. if (!v) nrerror("allocation failure in ivector()");
  802. return v-nl+NR_END;
  803. }
  804.  
  805.  
  806. int *ivector(long nl, long nh){
  807. int *v;
  808. v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
  809. if (!v) nrerror("allocation failure in ivector()");
  810. return v-nl+NR_END;
  811. }
  812.  
  813. long *lvector(long nl, long nh){
  814. long int *v;
  815. v=(long int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(long int)));
  816. if (!v) nrerror("allocation failure in lvector()");
  817. return v-nl+NR_END;
  818. }
  819.  
  820. double *dvector(long nl, long nh){
  821. double *v;
  822. v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
  823. if (!v) nrerror("allocation failure in dvector()");
  824. return v-nl+NR_END;
  825. }
  826.  
  827.  
  828. int **imatrix(long nrl, long nrh, long ncl, long nch){
  829. long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
  830. int **m;
  831. m=(int **)malloc((size_t)((nrow+NR_END)*sizeof(int*)));
  832. if (!m) nrerror("allocation failure 1 in imatrix()");
  833. m += NR_END;
  834. m -= nrl;
  835.  
  836. m[nrl]=(int *)malloc((size_t)((nrow*ncol+NR_END)*sizeof(int)));
  837. if (!m[nrl]) nrerror("allocation failure 2 in imatrix()");
  838. m[nrl] += NR_END;
  839. m[nrl] -= ncl;
  840.  
  841. for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
  842.  
  843. return m;
  844.  
  845. }
  846.  
  847.  
  848. double **dmatrix(long nrl, long nrh, long ncl, long nch){
  849. long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
  850. double **m;
  851. m=(double **)malloc((size_t)((nrow+NR_END)*sizeof(double*)));
  852. if (!m) nrerror("allocation failure 1 in dmatrix()");
  853. m += NR_END;
  854. m -= nrl;
  855.  
  856. m[nrl]=(double *)malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
  857. if (!m[nrl]) nrerror("allocation failure 2 in dmatrix()");
  858. m[nrl] += NR_END;
  859. m[nrl] -= ncl;
  860.  
  861. for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
  862.  
  863. return m;
  864.  
  865. }
  866.  
  867.  
  868.  
  869. double ***d3tensor(long nrl,long nrh,long ncl,long nch,long ndl,long ndh){
  870. long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
  871. double ***t;
  872.  
  873. t=(double ***) malloc((size_t)((nrow+NR_END)*sizeof(double**)));
  874. if(!t)nrerror("allocation failure 1 in d3tensor()");
  875. t += NR_END;
  876. t -= nrl;
  877.  
  878. t[nrl]=(double **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double*)));
  879. if(!t[nrl])nrerror("allocation failure 2 in d3tensor()");
  880. t[nrl] += NR_END;
  881. t[nrl] -= ncl;
  882.  
  883. t[nrl][ncl]=(double *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(double)));
  884. if(!t[nrl][ncl])nrerror("allocation failure 3 in d3tensor()");
  885. t[nrl][ncl] += NR_END;
  886. t[nrl][ncl] -= ndl;
  887.  
  888. for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
  889. for(i=nrl+1;i<=nrh;i++){
  890.     t[i]=t[i-1]+ncol;
  891.     t[i][ncl]=t[i-1][ncl]+ncol*ndep;
  892.     for(j=ncl+1;j<=nch;j++)t[i][j]=t[i][j-1]+ndep;
  893.     }
  894. return t;
  895. }
  896.  
  897. int am2num_c(c)
  898. {
  899. switch (c) {
  900.              case 'W':
  901.                     c=1; break;
  902.          case 'w':
  903.             c=26; break;
  904.              case 'F':
  905.                     c=2; break;
  906.                  case 'f':
  907.             c=27; break;
  908.              case 'Y':
  909.                     c=3; break;
  910.          case 'y':
  911.                         c=28; break;
  912.              case 'M':
  913.                     c=4; break;
  914.          case 'm':
  915.                         c=29; break;
  916.              case 'L':
  917.                     c=5; break;
  918.          case 'l':
  919.                         c=30; break;
  920.              case 'I':
  921.                 c=6; break;
  922.          case 'i':
  923.                         c=31; break;
  924.              case 'V':
  925.                 c=7; break;
  926.          case 'v':
  927.                         c=32; break;
  928.              case 'A':
  929.             c=8; break;
  930.          case 'a':
  931.                         c=33; break;
  932.              case 'C':
  933.                     c=9; break;
  934.          case 'c':
  935.                         c=34; break;
  936.          case 'G':
  937.             c=10; break;
  938.          case 'g':
  939.                         c=35; break;
  940.              case 'P':
  941.                     c=11; break;
  942.          case 'p':
  943.                         c=36; break;
  944.              case 'T':
  945.             c=12; break;
  946.          case 't':
  947.                         c=37; break;
  948.              case 'S':
  949.             c=13; break;
  950.          case 's':
  951.                         c=38; break;
  952.              case 'N':
  953.                     c=14; break;
  954.          case 'n':
  955.                         c=39; break;
  956.              case 'Q':
  957.                     c=15; break;
  958.          case 'q':
  959.                         c=40; break;
  960.              case 'D':
  961.                     c=16; break;
  962.          case 'd':
  963.                         c=41; break;
  964.              case 'E':
  965.                     c=17; break;
  966.          case 'e':
  967.                         c=42; break;
  968.              case 'H':
  969.                     c=18; break;
  970.          case 'h':
  971.                         c=43; break;
  972.              case 'R':
  973.                     c=19; break;
  974.          case 'r':
  975.                         c=44; break;
  976.              case 'K':
  977.                     c=20; break;
  978.          case 'k':
  979.                         c=45; break;
  980.              default :
  981.             c=0;
  982.         }
  983. return (c);
  984. }
  985.  
  986.  
  987. int am2num(c)
  988. {
  989. switch (c) {
  990.              case 'W': case 'w':
  991.                     c=1; break;
  992.              case 'F': case 'f':
  993.                     c=2; break;
  994.              case 'Y': case 'y':
  995.                     c=3; break;
  996.              case 'M': case 'm':
  997.                     c=4; break;
  998.              case 'L': case 'l':
  999.                     c=5; break;
  1000.              case 'I': case 'i':
  1001.                 c=6; break;
  1002.              case 'V': case 'v':
  1003.                 c=7; break;
  1004.              case 'A': case 'a':
  1005.             c=8; break;
  1006.              case 'C': case 'c':
  1007.                     c=9; break;
  1008.          case 'G': case 'g':
  1009.             c=10; break;
  1010.              case 'P': case 'p':
  1011.                     c=11; break;
  1012.              case 'T': case 't':
  1013.             c=12; break;
  1014.              case 'S': case 's':
  1015.             c=13; break;
  1016.              case 'N': case 'n':
  1017.                     c=14; break;
  1018.              case 'Q': case 'q':
  1019.                     c=15; break;
  1020.              case 'D': case 'd':
  1021.                     c=16; break;
  1022.              case 'E': case 'e':
  1023.                     c=17; break;
  1024.              case 'H': case 'h':
  1025.                     c=18; break;
  1026.              case 'R': case 'r':
  1027.                     c=19; break;
  1028.              case 'K': case 'k':
  1029.                     c=20; break;
  1030.              default :
  1031.             c=0;
  1032.         }
  1033. return (c);
  1034. }
  1035.  
  1036. int am2numBZX(c)
  1037. {
  1038. switch (c) {
  1039.                  case 'W': case 'w':
  1040.                         c=1; break;
  1041.                  case 'F': case 'f':
  1042.                         c=2; break;
  1043.                  case 'Y': case 'y':
  1044.                         c=3; break;
  1045.                  case 'M': case 'm':
  1046.                         c=4; break;
  1047.                  case 'L': case 'l':
  1048.                         c=5; break;
  1049.                  case 'I': case 'i':
  1050.                         c=6; break;
  1051.                  case 'V': case 'v':
  1052.                         c=7; break;
  1053.                  case 'A': case 'a':
  1054.                         c=8; break;
  1055.                  case 'C': case 'c':
  1056.                         c=9; break;
  1057.                  case 'G': case 'g':
  1058.                         c=10; break;
  1059.                  case 'P': case 'p':
  1060.                         c=11; break;
  1061.                  case 'T': case 't':
  1062.                         c=12; break;
  1063.                  case 'S': case 's':
  1064.                         c=13; break;
  1065.                  case 'N': case 'n':
  1066.                         c=14; break;
  1067.                  case 'Q': case 'q':
  1068.                         c=15; break;
  1069.                  case 'D': case 'd':
  1070.                         c=16; break;
  1071.                  case 'E': case 'e':
  1072.                         c=17; break;
  1073.                  case 'H': case 'h':
  1074.                         c=18; break;
  1075.                  case 'R': case 'r':
  1076.                         c=19; break;
  1077.                  case 'K': case 'k':
  1078.                         c=20; break;
  1079.                  case 'B': case 'b':
  1080.                         c=21; break;
  1081.                  case 'Z': case 'z':
  1082.                         c=22; break;
  1083.                  case 'X': case 'x':
  1084.                         c=23; break;
  1085.                  case '*':
  1086.                         c=24; break;
  1087.                  default :
  1088.                         c=0;
  1089.                 }
  1090. return (c);
  1091. }
  1092.  
  1093.  
  1094. static char str[MAXSTR+1];
  1095.  
  1096. /*char **aname, **aseq;
  1097. int nal, alilen, *astart, *alen;
  1098. int **alignment;
  1099. */
  1100.  
  1101.  
  1102.  
  1103. static void *mymalloc(int size);
  1104. char *strsave(char *str);
  1105. char *strnsave(char *str, int l);
  1106. static char **incbuf(int n, char **was);
  1107. static int *incibuf(int n, int *was);
  1108.  
  1109. void readali(char *filename);
  1110. int **ali_char2int(char **aseq,int start_num, int start_seq);
  1111. int **read_alignment2int(char *filename,int start_num,int start_seq);
  1112.  
  1113. void counter(int b);
  1114. double effective_number(int **ali, int *marks, int n, int start, int end);
  1115. double effective_number_nogaps(int **ali, int *marks, int n, int start, int end);
  1116. double effective_number_nogaps_expos(int **ali, int *marks, int n, int start, int end, int pos);
  1117.  
  1118.  
  1119.  
  1120. static void *mymalloc(size)
  1121. int size;
  1122. {
  1123.     void *buf;
  1124.  
  1125.     if ((buf = malloc(size)) == NULL) {
  1126.         fprintf(stderr, "Not enough memory: %d\n", size);
  1127.         exit(1);
  1128.     }
  1129.     return buf;
  1130. }
  1131.  
  1132. char *strsave(str)
  1133. char *str;
  1134. {
  1135.     char *buf;
  1136.     int l;
  1137.  
  1138.     l = strlen(str);
  1139.     buf = mymalloc(l + 1);
  1140.     strcpy(buf, str);
  1141.     return buf;
  1142. }
  1143.  
  1144. char *strnsave(str, l)
  1145. char *str;
  1146. int l;
  1147. {
  1148.     char *buf;
  1149.  
  1150.     buf = mymalloc(l + 1);
  1151.     memcpy(buf, str, l);
  1152.     buf[l] = '\0';
  1153.     return buf;
  1154. }
  1155.  
  1156. static char **incbuf(n, was)
  1157. int n;
  1158. char **was;
  1159. {
  1160.     char **buf;
  1161.     char *aaa;
  1162.  
  1163.     buf = mymalloc((n+1) * sizeof(buf[0]));
  1164.     if (n > 0) {
  1165.         memcpy(buf, was, n * sizeof(was[0]));
  1166.         free(was);
  1167.     }
  1168.     buf[n] = NULL;
  1169.     return buf;
  1170. }
  1171.  
  1172. static int *incibuf(n, was)
  1173. int n, *was;
  1174. {
  1175.     int *ibuf;
  1176.  
  1177.     ibuf = mymalloc((n+1) * sizeof(ibuf[0]));
  1178.     if (n > 0) {
  1179.         memcpy(ibuf, was, n * sizeof(was[0]));
  1180.         free(was);
  1181.     }
  1182.     ibuf[n] = 0;
  1183.     return ibuf;
  1184. }
  1185. void err_readali(int err_num)
  1186. {
  1187.     fprintf(stderr,"Error with reading alignment: %d\n",err_num);
  1188. }
  1189.  
  1190. void readali(filename)
  1191. char *filename;
  1192. {
  1193.     FILE *fp;
  1194.     char *s, *ss, *seqbuf;
  1195.     int n, l, len, len0;
  1196.     int ii,mark=1;
  1197.  
  1198.     if ((fp = fopen(filename, "r")) == NULL) {
  1199.         fprintf(stderr, "No such file: \"%s\"\n", filename);
  1200.         err_readali(1);
  1201.         ;exit(1);
  1202.     }
  1203.     alilen = 0;
  1204.     nal = 0;
  1205.     n = 0;
  1206.     if(fgets(str, MAXSTR, fp) != NULL) {
  1207.         if(strncmp(str,"CLUSTAL",7)!=0){rewind(fp);}
  1208.                        }
  1209.     while (fgets(str, MAXSTR, fp) != NULL) {
  1210.         for (ss = str; isspace(*ss); ss++) ;
  1211.         if ((ii<=ss-str)&&(mark==0)) {continue;}
  1212.         if (*ss == '\0') {
  1213.             if (n == 0) {
  1214.                 continue;
  1215.             }
  1216.             if (nal == 0) {
  1217.                 if (n == 0) {
  1218.                     fprintf(stderr, "No alignments read\n");
  1219.                     err_readali(2);
  1220.                     exit(1);
  1221.                 }
  1222.                 nal = n;
  1223.             } else if (n != nal) {
  1224.                 fprintf(stderr, "Wrong nal, was: %d, now: %d\n", nal, n);
  1225.                 err_readali(3); exit(1);
  1226.             }
  1227.             n = 0;
  1228.             continue;
  1229.         }
  1230.         for (s = ss; *s != '\0' && !isspace(*s); s++) ;
  1231.         *s++ = '\0';
  1232.         if (nal == 0) {
  1233.             astart = incibuf(n, astart);
  1234.             alen = incibuf(n, alen);
  1235.             aseq = incbuf(n, aseq);
  1236.             aname = incbuf(n, aname);
  1237.             aname[n] = strsave(ss);
  1238.         } else {
  1239.             if (n < 0 || n >= nal) {
  1240.                 fprintf(stderr, "Bad sequence number: %d of %d\n", n, nal);
  1241.                 err_readali(4);  exit(1);
  1242.             }
  1243.             if (strcmp(ss, aname[n]) != 0) {
  1244.                 fprintf(stderr, "Names do not match");
  1245.                 fprintf(stderr, ", was: %s, now: %s\n", aname[n], ss);
  1246.                 err_readali(5);  exit(1);
  1247.             }
  1248.         }
  1249.         for (ss = s; isspace(*ss); ss++);
  1250.         if(mark==1){
  1251.         ii = ss-str;
  1252.         mark=0;}
  1253.        
  1254.         for (s = ss; isdigit(*s); s++) ;
  1255.         if (isspace(*s)) {
  1256.             if (nal == 0) {
  1257.                 astart[n] = atoi(ss);
  1258.             }
  1259.             for (ss = s; isspace(*ss); ss++);
  1260.         }
  1261.         for (s = ss, l = 0; *s != '\0' && !isspace(*s); s++) {
  1262.             if (isalpha(*s)) {
  1263.                 l++;
  1264.             }
  1265.         }
  1266.         len = s - ss;
  1267.         if (n == 0) {
  1268.             len0 = len;
  1269.             alilen += len;
  1270.         } else if (len != len0) {
  1271.             fprintf(stderr, "wrong len for %s", aname[n]);
  1272.             fprintf(stderr, ", was: %d, now: %d\n", len0, len);
  1273.             err_readali(6); exit(1);
  1274.         }
  1275.         alen[n] += l;
  1276.         if (aseq[n] == NULL) {
  1277.             aseq[n] = strnsave(ss, len);
  1278.         } else {
  1279.             seqbuf = mymalloc(alilen+1);
  1280.             memcpy(seqbuf, aseq[n], alilen-len);
  1281.             free(aseq[n]);
  1282.             aseq[n] = seqbuf;
  1283.             memcpy(seqbuf+alilen-len, ss, len);
  1284.             seqbuf[alilen] = '\0';
  1285.         }
  1286.         n++;
  1287.     }
  1288.     if (nal == 0) {
  1289.         if (n == 0) {
  1290.             fprintf(stderr, "No alignments read\n");
  1291.             err_readali(7);  exit(1);
  1292.         }
  1293.         nal = n;
  1294.     } else if (n != 0 && n != nal) {
  1295.         fprintf(stderr, "Wrong nal, was: %d, now: %d\n", nal, n);
  1296.         err_readali(8);  exit(1);
  1297.     }
  1298.     fclose(fp);
  1299. }
  1300.  
  1301. static void printali(char *argt, int chunk, int n, int len, char **name, char **seq, int *start,int *csv)
  1302. {
  1303.         int i, j, k, jj, mlen, sta, *pos;
  1304.     char csv_str[100], arg_t[500];
  1305.         char *sq;
  1306.     int *isq;
  1307.     FILE *fpp;
  1308.  
  1309.     strcpy(arg_t,argt);
  1310.     fpp=fopen(arg_t,"w");
  1311.     strcpy(csv_str, "Conservation:");
  1312.         pos = mymalloc(n * sizeof(pos[0]));
  1313.         memcpy(pos, start, n * sizeof(pos[0]));
  1314.         for (i=0; i < n && start[i] == 0; i++) ;
  1315.         sta = (i < n);
  1316.         for (i=1, mlen=strlen(name[0]); i < n; i++) {
  1317.                 if (mlen < strlen(name[i])) {
  1318.                         mlen = strlen(name[i]);
  1319.                 }
  1320.         }
  1321.         jj = 0;
  1322.         do {
  1323.                 if (jj > 0) {
  1324.                         fprintf(fpp, "\n");
  1325.                 }
  1326.                 for (i=0; i < n; i++) {
  1327.             if(i==0) {
  1328.                 fprintf(fpp,"\n\n");
  1329.                 for(k=0;k<mlen+4;k++){
  1330.                     if(k<strlen(csv_str)){
  1331.                     fprintf(fpp,"%c",csv_str[k]);}
  1332.                     else fprintf(fpp," ");
  1333.                            }
  1334.                 if (sta) {
  1335.                     fprintf(fpp,"%4d ", pos[i]);
  1336.                 }
  1337.                 isq = csv + jj;
  1338.                 for (j=0; j+jj < len && j < chunk; j++) {
  1339.                     if (isdigit(isq[j+1])) {
  1340.                         pos[i]++;
  1341.                     }
  1342.                     fprintf(fpp,"%d", isq[j+1]);
  1343.                 }
  1344.                 if (sta) {
  1345.                     fprintf(fpp, " %4d", pos[i]-1);
  1346.                 }
  1347.                 fprintf(fpp, "\n");
  1348.                   }
  1349.  
  1350.             for(k=0;k<mlen+4;k++){
  1351.             if(k<strlen(name[i])){
  1352.             fprintf(fpp,"%c",name[i][k]);}
  1353.             else fprintf(fpp," ");
  1354.                        }           
  1355.                         if (sta) {
  1356.                                 fprintf(fpp, "%4d ", pos[i]);
  1357.                         }
  1358.                         sq = seq[i] + jj;
  1359.                         for (j=0; j+jj < len && j < chunk; j++) {
  1360.                                 if (isalpha(sq[j])) {
  1361.                                         pos[i]++;
  1362.                                 }
  1363.                                 fprintf(fpp, "%c", sq[j]);
  1364.                         }
  1365.                         if (sta) {
  1366.                                 fprintf(fpp, " %4d", pos[i]-1);
  1367.                         }
  1368.                         fprintf(fpp, "\n");
  1369.                 }
  1370.                 jj += j;
  1371.         } while (jj < len);
  1372.         free(pos);
  1373.     fclose(fpp);
  1374. }
  1375.  
  1376. int **ali_char2int(char **aseq, int start_num, int start_seq){
  1377. /* fills the alignment ali[start_num..start_num+nal-1][start_seq..start_seq+alilen-1]
  1378. convetring charater to integer from aseq[0..nal-1][0..alilen-1]
  1379. */
  1380.  
  1381. int i,j,end_num,end_seq;
  1382. int **ali;
  1383. end_num=start_num+nal-1;
  1384. end_seq=start_seq+alilen-1;
  1385. ali=imatrix(start_num,end_num,start_seq,end_seq);
  1386. for(i=start_num;i<=end_num;++i)for(j=start_seq;j<=end_seq;++j)ali[i][j]=am2num(aseq[i-start_num][j-start_seq]);
  1387. return ali;
  1388. }
  1389.  
  1390. int **read_alignment2int(char *filename,int start_num,int start_seq){
  1391. int **ali;
  1392. readali(filename);
  1393. ali=ali_char2int(aseq,start_num,start_seq);
  1394. return ali;
  1395. }
  1396.  
  1397.  
  1398.  
  1399. double effective_number(int **ali, int *marks, int n, int start, int end){
  1400.  
  1401. /* from the alignment of n sequences ali[1..n][1..l]
  1402. calculates effective number of sequences that are marked by 1 in mark[1..n]
  1403. for the segment of positions ali[][start..end]
  1404. Neff=ln(1-0.05*N-of-different-letters-per-site)/ln(0.95)
  1405. */
  1406.  
  1407. int i,k,a,flag,ngc;
  1408. int *amco,lettercount=0,sitecount=0,gsitecount=0;
  1409. int *nogaps;
  1410. double letpersite=0,neff,randomcount=0.0;
  1411. amco=ivector(0,20);
  1412. nogaps=ivector(0,n);
  1413. for(i=0;i<=n;++i)nogaps[i]=0;
  1414. for(k=start;k<=end;++k){
  1415.     for(a=0;a<=20;++a)amco[a]=0;
  1416.     for(i=1;i<=n;++i)if(marks[i]==1)amco[ali[i][k]]++;
  1417.     flag=0;for(a=1;a<=20;++a)if(amco[a]>0){flag=1;lettercount++;}
  1418.     if(flag==1)sitecount++;
  1419.                }
  1420. letpersite=1.0*lettercount/sitecount;
  1421. for(k=start;k<=end;++k){
  1422.     ngc=0;for(i=1;i<=n;++i)if(marks[i]==1)if(ali[i][k]!=0)ngc++;
  1423.     if(ngc!=0)gsitecount++;
  1424.     nogaps[ngc]++;
  1425.                }
  1426. for(i=0;i<=n;++i)fprintf(stderr,"%3d %4d\n",i,nogaps[i]);
  1427.  
  1428. if(gsitecount!=sitecount){fprintf(stderr,"Counts didn't match in \"effective_number\": FATAL\n");exit(0);}
  1429.  
  1430. for(i=1;i<=n;++i)randomcount+=nogaps[i]*20.0*(1.0-exp(-0.05129329438755*i));
  1431. randomcount=randomcount/gsitecount;
  1432. letpersite=letpersite*20.0*(1.0-exp(-0.05129329438755*n))/randomcount;
  1433.  
  1434.  
  1435.  
  1436. neff=-log(1.0-0.05*letpersite)/0.05129329438755;
  1437. return neff;
  1438. }
  1439.  
  1440.  
  1441.  
  1442. double effective_number_nogaps(int **ali, int *marks, int n, int start, int end){
  1443.  
  1444. /* from the alignment of n sequences ali[1..n][1..l]
  1445. calculates effective number of sequences that are marked by 1 in mark[1..n]
  1446. for the segment of positions ali[][start..end]
  1447. Neff=ln(1-0.05*N-of-different-letters-per-site)/ln(0.95)
  1448. */
  1449.  
  1450. int i,k,a,flag;
  1451. int *amco,lettercount=0,sitecount=0;
  1452. double letpersite=0,neff;
  1453. amco=ivector(0,20);
  1454. for(k=start;k<=end;++k){
  1455.     flag=0;for(i=1;i<=n;++i)if(marks[i]==1 && ali[i][k]==0)flag=1;
  1456.     if(flag==1)continue;
  1457.     for(a=0;a<=20;++a)amco[a]=0;
  1458.     for(i=1;i<=n;++i)if(marks[i]==1)amco[ali[i][k]]++;
  1459.     flag=0;for(a=1;a<=20;++a)if(amco[a]>0){flag=1;lettercount++;}
  1460.     if(flag==1)sitecount++;
  1461.                }
  1462. if(sitecount==0)letpersite=0;
  1463. else letpersite=1.0*lettercount/sitecount;
  1464.  
  1465.  
  1466. neff=-log(1.0-0.05*letpersite)/0.05129329438755;
  1467. return neff;
  1468. }
  1469.  
  1470. int *letters;
  1471.  
  1472. void **freq(int **ali,double **f,int *num_gaps,int *effindiarr,double gap_threshold)
  1473. {
  1474.     int i,j,k,effnumind;
  1475.     int count[21];
  1476.    
  1477.     letters = ivector(0, alilen+1);
  1478.     letters[0]=1;
  1479.  
  1480.     /* find the number of frequences at each position */
  1481.     effnumind=0;
  1482.     for(j=1;j<=alilen;j++){
  1483.         for(i=0;i<=20;i++) count[i]=0;
  1484.         for(i=1;i<=nal;i++) {
  1485.             if(ali[i][j]<=20) {count[ali[i][j]]++;}
  1486.             else {if(ali[i][j]>25&&ali[i][j]<=45)
  1487.                 {count[ali[i][j]-25]++;}
  1488.                   else {
  1489.                 fprintf(stderr,"not good number for AA\n");
  1490.                 exit(0);
  1491.                    }
  1492.                  }
  1493.                     }
  1494.         num_gaps[j] = count[0];
  1495.         f[0][j] = count[0]*1.0/nal;
  1496.         if(f[0][j]>1) {
  1497.             fprintf(stderr,"gap number>total number\n");
  1498.             exit(0);
  1499.                   }
  1500.         if(f[0][j]>=gap_threshold) {   /* ignore the case where gaps occur >= gap_threshold(percentage of gaps)  */
  1501.             f[0][j]=INDI;
  1502.             continue;
  1503.                 }
  1504.         effnumind++;
  1505.         effindiarr[effnumind]=j;
  1506.         count[0]=nal-count[0];
  1507.         letters[j] = count[0];
  1508.         if(count[0]<=0){
  1509.             fprintf(stderr, "count[0] less than 0: %d\n",count[0]);
  1510.             exit(0);
  1511.                    }
  1512.         for(k=1;k<=20;k++){
  1513.             f[k][j]=count[k]*1.0/count[0];
  1514.                   }
  1515.     }
  1516.     effindiarr[effnumind+1]=INDI;/*set the last element negative*/
  1517.     effindiarr[0]=effnumind;
  1518. }          
  1519.  
  1520. /* calculating the standard error of frequencies */
  1521. double **sigmaf(double **f)
  1522. {
  1523.     double **sigma;
  1524.     int i,j;
  1525.  
  1526.     sigma = dmatrix(0,20,0,alilen);
  1527.     if(letters[0]!=1) {
  1528.         fprintf(stderr, "letters not defined previously\n");
  1529.         exit(0);
  1530.               }
  1531.     for(j=1;j<=alilen;j++) {
  1532.         if(letters[j]==0) {
  1533.             fprintf(stderr,"The column contains no letters:%d\n",j);
  1534.             for(i=1;i<=20;i++){sigma[i][j]=0;}
  1535.             continue;
  1536.                   }
  1537.         for(i=1;i<=20;i++) {
  1538.             if(f[i][j]==0) {sigma[i][j]=0;continue;}
  1539.             sigma[i][j]=sqrt(f[i][j]*(1-f[i][j])/letters[j]);
  1540.                    }
  1541.                 }
  1542.     return sigma;      
  1543. }
  1544.  
  1545. int totalw;
  1546. double *oaf_ip, *h_oaf_ip; /*unweighted,henikoff frequency at position ip */
  1547. double *overall_freq(int **ali, int startp, int endp, int *mark);
  1548. double *overall_freq_wgt(int **ali,int startp,int endp,int *mark,double *wgt);
  1549. double *h_weight(int **ali, int ip)
  1550. {
  1551.     double *hwt;
  1552.     int count[21];
  1553.         int amtypes;
  1554.     int i,j,k;
  1555.     int maxstart, maxs, miniend,minie;
  1556.     int *mark;
  1557.     int gapcount,totalcount,wsign;
  1558.  
  1559.     hwt = dvector(0, nal);
  1560.     mark = ivector(0, nal);
  1561.     oaf_ip = dvector(0,20);
  1562.     h_oaf_ip = dvector(0,20);
  1563.     for(i=1;i<=20;i++) h_oaf_ip[i] = 0;
  1564.  
  1565.     /* mark the sequences with gaps */
  1566.     for(i=1;i<=nal;i++) {
  1567.         // NEW: fix a typo 05/06/03
  1568.         if((ali[i][ip]>0&&ali[i][ip]<=20)||(ali[i][ip]>25&&ali[i][ip]<=45)) mark[i]=1;
  1569.         else if (ali[i][ip]==0)mark[i]=0;
  1570.                 }
  1571.    
  1572.     /* find the maxstart and miniend positions */
  1573.     maxstart = 1;
  1574.     for(i=1;i<=nal;i++){
  1575.         if(mark[i]==0) continue;
  1576.         maxs = 1;
  1577.             for(j=1;j<=alilen;j++) {
  1578.             if(ali[i][j]==0) maxs++;
  1579.             if(ali[i][j]>0) break;
  1580.                     }
  1581.         if(maxstart<maxs) maxstart = maxs;
  1582.                }
  1583.     miniend = alilen;
  1584.     for(i=1;i<=nal;i++){
  1585.         if(!mark[i]) continue;
  1586.         minie = alilen;
  1587.         for(j=alilen;j>0;j--) {
  1588.             if(ali[i][j]==0) minie--;
  1589.             if(ali[i][j]>0) break;
  1590.                       }
  1591.         if(miniend>minie) miniend = minie;
  1592.                }
  1593.     /* NEW: 05/06/03 */
  1594.     if(maxstart > miniend) maxstart = miniend;
  1595.  
  1596. /* special case where only one site is not gap */
  1597.         j=0;
  1598.         for(i=1;i<=nal;i++) {if(mark[i]>0) j++;}
  1599.         if(j==1) {
  1600.                 for(i=1;i<=nal;i++) hwt[i]=mark[i];
  1601.                 h_oaf_ip = overall_freq_wgt(ali,maxstart,miniend,mark,hwt);
  1602.                 oaf_ip = overall_freq(ali,1,alilen,mark);
  1603.                 return hwt;
  1604.                  }
  1605.  
  1606.     for(i=1;i<=nal;i++) hwt[i]=0;
  1607.     for(j=maxstart;j<=miniend;j++){
  1608.        
  1609.         amtypes = 0;
  1610.         for(i=0;i<=20;i++) count[i]=0;
  1611.         for(i=1;i<=nal;i++){
  1612.             if(mark[i]==0) continue;
  1613.             if(ali[i][j]>=0&&ali[i][j]<=20) count[ali[i][j]]++;
  1614.             else if(ali[i][j]>25&&ali[i][j]<=45) count[ali[i][j]-25]++;
  1615.                    }
  1616.         for(i=0;i<=20;i++) {
  1617.             if(count[i]>0) amtypes++;
  1618.                    }
  1619.         if(amtypes==1) continue; /* identical positions are excluded */
  1620.         gapcount=totalcount=0;
  1621.         for(i=1;i<=nal;i++){
  1622.             if(mark[i]==0) continue;
  1623.             if(ali[i][j]==0) gapcount++;
  1624.             totalcount++;
  1625.                    }
  1626.         if(gapcount>totalcount*0.5) continue;/*gap>50% excluded*/
  1627.  
  1628.         for(i=1;i<=nal;i++){
  1629.             if(mark[i]==0) continue;
  1630.             if(ali[i][j]>=0&&ali[i][j]<=20) {
  1631.                 if(count[ali[i][j]]>0)  {
  1632.                     hwt[i]+=1.0/(count[ali[i][j]]*amtypes);
  1633.                             }
  1634.                                }
  1635.             if(ali[i][j]>25&&ali[i][j]<=45) {
  1636.                 if(count[ali[i][j]-25]>0)  {
  1637.                     hwt[i]+=1.0/(count[ali[i][j]-25]*amtypes);
  1638.                                }
  1639.                             }
  1640.                    }
  1641.                 }
  1642.     /* test if all henikoff weights are zero for all sequences */
  1643.     wsign=0;
  1644.     for(i=1;i<=nal;i++) {
  1645.         if(hwt[i]>0) { wsign=1;break;}
  1646.                 }
  1647.     gapcount=0;
  1648.     if(wsign==0) {
  1649.         for(i=1;i<=nal;i++){
  1650.             if(ali[i][ip]==0) gapcount++;
  1651.                    }
  1652.         if(gapcount==nal){fprintf(stderr, "This position contains only gaps\n");exit(0);}
  1653.         for(i=1;i<=nal;i++){
  1654.             if(ali[i][ip]!=0) hwt[i]=1.0/(nal-gapcount);
  1655.                    }
  1656.              }
  1657.  
  1658.     h_oaf_ip = overall_freq_wgt(ali,maxstart,miniend,mark,hwt);
  1659.     oaf_ip = overall_freq(ali,1,alilen,mark);
  1660.     return hwt;
  1661. }
  1662.    
  1663. double **u_oaf,**h_oaf; /* unweighted and henikoff overall frequency */
  1664. double **h_freq(int **ali, double **f, double **hfr)
  1665. {
  1666.     int i,j;   
  1667.     double *hwt;
  1668.     double sumofweight;
  1669.  
  1670.     u_oaf = dmatrix(0,20,0,alilen);
  1671.     h_oaf = dmatrix(0,20,0,alilen);
  1672.    
  1673.     for(j=1;j<=alilen;j++) {
  1674.         if(f[0][j]==INDI) {hfr[0][j]=INDI;continue;}
  1675.         hwt = h_weight(ali, j); /* assign position specific weight */
  1676.         for(i=1;i<=20;i++) {
  1677.             u_oaf[i][j]=oaf_ip[i];
  1678.             h_oaf[i][j]=h_oaf_ip[i];
  1679.                    }
  1680.  
  1681.         sumofweight = 0;
  1682.         for(i=1;i<=nal;i++) {
  1683.             if(ali[i][j]>0) sumofweight+=hwt[i];
  1684.                     }
  1685.         for(i=1;i<=nal;i++) {
  1686.             if(ali[i][j]>0&&ali[i][j]<=20) {
  1687.                 hfr[ali[i][j]][j]+=hwt[i]/sumofweight;
  1688.                             }
  1689.             else if(ali[i][j]>25&&ali[i][j]<=45) {
  1690.                 hfr[ali[i][j]-25][j]+=hwt[i]/sumofweight;
  1691.                                  }
  1692.                     }
  1693.                 }
  1694. }
  1695.  
  1696.  
  1697. double *entro_conv(double **f, int **ali, double *econv)
  1698. {
  1699.     int i,j;
  1700.  
  1701.     for (j=1;j<=alilen;j++){
  1702.         econv[j]=0;
  1703.         if(f[0][j]==INDI) {econv[j]=INDI;continue;}
  1704.         for(i=1;i<=20;i++) {
  1705.             if(f[i][j]==0) continue;
  1706.             econv[j]+=f[i][j]*log(f[i][j]);
  1707.                    }
  1708.                 }
  1709. }
  1710.        
  1711. double *overall_freq(int **ali, int startp, int endp, int *mark)
  1712. {
  1713.     double *oaf;
  1714.     int gapcount,totalcount;
  1715.     int total=0;
  1716.     int i,j;
  1717.  
  1718.     oaf = dvector(0,20);
  1719.     for(i=1;i<=20;i++) oaf[i]=0;
  1720.     if(startp>endp) {
  1721.         fprintf(stderr, "start position larger than ending position\n");
  1722.         exit(1);
  1723.             }
  1724.     for(j=startp;j<=endp;j++) {
  1725.         /* excluding those that have >50% gaps*/
  1726.         gapcount=totalcount=0;
  1727.         for(i=1;i<=nal;i++){
  1728.             if(mark[i]==0) continue;
  1729.             totalcount++;
  1730.             if(ali[i][j]==0) gapcount++;
  1731.                    }
  1732.         if(gapcount>totalcount*0.5) continue;
  1733.  
  1734.         for(i=1;i<=nal;i++) {
  1735.             if(mark[i]==0) continue;
  1736.             if(ali[i][j]>0&&ali[i][j]<=20) {
  1737.                 oaf[ali[i][j]]+=1;
  1738.                 total++;
  1739.                             }
  1740.             if(ali[i][j]>25&&ali[i][j]<=45) {
  1741.                 oaf[ali[i][j]-25]+=1;
  1742.                 total++;
  1743.                             }
  1744.                   }
  1745.                 }
  1746.     for(i=1;i<=20;i++) {
  1747.         oaf[i]=oaf[i]/total;
  1748.                }       
  1749.     return oaf;
  1750. }
  1751.  
  1752. double *overall_freq_wgt(int **ali,int startp,int endp,int *mark,double *wgt)
  1753. {
  1754.         double *oaf;
  1755.         double total=0;
  1756.     int totalcount,gapcount;
  1757.         int i,j;
  1758.  
  1759.         oaf = dvector(0,20);
  1760.         for(i=1;i<=20;i++) oaf[i]=0;
  1761.         if(startp>endp) {
  1762.                 fprintf(stderr, "start position larger than ending position\n");                exit(1);
  1763.                         }
  1764.         for(j=startp;j<=endp;j++) {
  1765.         /* excluding those that have >50% gaps*/
  1766.         gapcount=totalcount=0;
  1767.         for(i=1;i<=nal;i++){
  1768.             if(mark[i]==0) continue;
  1769.             totalcount++;
  1770.             if(ali[i][j]==0) gapcount++;
  1771.                    }
  1772.         if(gapcount>totalcount*0.5) continue;
  1773.                 for(i=1;i<=nal;i++) {
  1774.                         if(mark[i]==0) continue;
  1775.                         if(ali[i][j]>0&&ali[i][j]<=20) {
  1776.                                 oaf[ali[i][j]]+=wgt[i];
  1777.                                 total+=wgt[i];
  1778.                                                         }
  1779.                         if(ali[i][j]>25&&ali[i][j]<=45) {
  1780.                                 oaf[ali[i][j]-25]+=wgt[i];
  1781.                                 total+=wgt[i];
  1782.                                                         }
  1783.                                   }
  1784.     if(total==0) fprintf(stderr,"total=0\n");
  1785.                                 }
  1786.         for(i=1;i<=20;i++) {
  1787.                 oaf[i]=oaf[i]/total;
  1788.                            }
  1789.         return oaf;
  1790. }
  1791.  
  1792. double **ic_freq(int **ali, double **f, double **icf)
  1793. {
  1794.         int i,j,k;
  1795.         int ele;
  1796.         double  *effnu;
  1797.         int *mark;
  1798.  
  1799.         mark = ivector(0,nal+10);
  1800.         effnu = dvector(0,20);
  1801.         for(i=0;i<=20;i++)
  1802.                 for(j=1;j<=alilen;j++)
  1803.                         icf[i][j]=0;
  1804.         for(i=0;i<=nal;i++) mark[i]=0;
  1805.  
  1806.         for(j=1;j<=alilen;j++){
  1807.                 if(f[0][j]==INDI) {icf[0][j]=INDI;continue;}
  1808.                 for(k=0;k<=20;++k)effnu[k]=0;
  1809.                 for(k=1;k<=20;++k){
  1810.                         for(i=1;i<=nal;++i){
  1811.                                 mark[i]=0;
  1812.                                 ele=ali[i][j];
  1813.                                 if(ele==k)mark[i]=1;
  1814.                                 ele=ali[i][j]-25;
  1815.                                 if(ele==k) mark[i]=1;
  1816.                                          }
  1817.                         effnu[k]=effective_number_nogaps(ali,mark,nal,1,alilen);
  1818.                         effnu[0]+=effnu[k];
  1819.                                   }
  1820.                 if(effnu[0]==0){fprintf(stderr,"all counts are zeros at the column %d: FATAL\n",j);exit(0);}
  1821.                 for(k=1;k<=20;k++) effnu[k]=effnu[k]/effnu[0];
  1822.                 for(k=1;k<=20;k++) icf[k][j] = effnu[k];
  1823.                                }
  1824. }
  1825.  
  1826. double *variance_conv(double **f, int **ali, double **oaf, double *vconv)
  1827. {
  1828.     int i,j;
  1829.  
  1830.     for(i=1;i<=alilen;i++) vconv[i]=0;
  1831.     for(j=1;j<=alilen;j++) {
  1832.         if(f[0][j]==INDI) {vconv[j]=INDI;continue;}
  1833.         for(i=1;i<=20;i++) {
  1834.             vconv[j]+=(f[i][j]-oaf[i][j])*(f[i][j]-oaf[i][j]);
  1835.                    }
  1836.         vconv[j]=sqrt(vconv[j]);
  1837.                    }
  1838. }
  1839.  
  1840. double *pairs_conv(double **f,int **ali,int **matrix1,int indx,double *pconv)
  1841. {
  1842.     int i,j,k;
  1843.     double **matrix2, **matrix3;
  1844.     FILE *fp;
  1845.  
  1846.     matrix2 = dmatrix(0,25,0,25);
  1847.     matrix3 = dmatrix(0,25,0,25);
  1848.     for(i=0;i<=25;i++){
  1849.         for(j=0;j<=25;j++) {
  1850.             matrix2[i][j]=0;
  1851.             matrix3[i][j]=0;
  1852.                    }
  1853.               }
  1854.  
  1855.     /* get the matrices */
  1856.     for(i=1;i<=24;i++){
  1857.         if(matrix1[i][i]==0) {
  1858.             fprintf(stderr, "diagonal elements zero: %d\n",i);
  1859.             exit(0);
  1860.                     }
  1861.         for(j=1;j<=24;j++)  {
  1862.             matrix2[i][j]=matrix1[i][j]*1.0/sqrt(matrix1[i][i]*matrix1[j][j]);
  1863.             matrix3[i][j]=matrix1[i][j]*2-0.5*(matrix1[i][i]+matrix1[j][j]);
  1864.                     }
  1865.                       }
  1866.    
  1867.     for(j=1;j<=alilen;j++) {
  1868.         if(f[0][j]==INDI){pconv[j]=INDI; continue;}
  1869.         pconv[j]=0;
  1870.         for(i=1;i<=20;i++) {
  1871.             for(k=1;k<=20;k++) {
  1872.                if(indx==0){
  1873.                 pconv[j]+=(f[k][j]*f[i][j])*matrix1[i][k];
  1874.                 continue;}
  1875.                if(indx==1){
  1876.                 pconv[j]+=(f[k][j]*f[i][j])*matrix2[i][k];
  1877.                 continue;}
  1878.                if(indx==2){
  1879.                 pconv[j]+=(f[k][j]*f[i][j])*matrix3[i][k];
  1880.                     continue;}
  1881.                fprintf(stderr,"not good index number of matrix\n");
  1882.                exit(0);
  1883.                        }
  1884.                    }
  1885.                 }
  1886. }
  1887.  
  1888.  
  1889.        
  1890.  
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!