/*  This program counts the number of strongly polar (D,E,N,Q,K,R,H), anionic
*/
/*  (D,E), cationic (R,K,H), semipolar (G,S,T,C), hydroxy (S,T), total
*/
/*  hydrophobic (L,I,V,M,F,Y,W,A,P), amphipathic aromatic (Y,W), total
*/
/*  aromatic (Y,W,F), total aliphatic (L,I,V,M,A,P), branched chain aliphatic
*/
/*  (L,I,V).  The percentages are given as averages in a variable window.
*/
/*  8/25/94 st
*/

#include <stdio.h>
#include <string.h>

main(ac,av)
int ac;
char *av[];

{
        FILE *ALIGN,*ALLout,*AAFout;
        char code[6],seqarray[30][1000],str[200];
        int
ansarrleng,begin,c,e,end,i,length,numlines,numseq,pos,posn,algn_leng,row,start,s
top,wsize,x,y,z;
        float
Acnt,Rcnt,Ncnt,Dcnt,Gcnt,Ccnt,Qcnt,Ecnt,Hcnt,Icnt,Lcnt,Kcnt,Mcnt,Fcnt,Pcnt,Scnt,
Tcnt,Wcnt,Ycnt,Vcnt;
        float
str_polar,anionic,cationic,semipolar,hydroxy,tot_hydrophob,amphi_arom,tot_arom,t
ot_aliph,br_chain_aliph;
        float aafreq[29][1000],scorearray1[10][1000],scorearray2[10][1000],sum;

        wsize=atoi(av[1]);
        if(ac!=5 && wsize%2==0)
                {
                printf("chkall gathers statistical data on all the positions in
an alignment.\nThe amino acid groups searched for are\n\n");
                printf("  Strongly Polar            D E N Q K R H\n");
                printf("  Anionic                   D E\n");
                printf("  Cationic                  R K H\n");
                printf("  Semipolar                 G S T C\n");
                printf("  Hydroxy                   S T\n");
                printf("  Total Hydrophobic         L I V M F Y W A P\n");
                printf("  Amphipathic Aromatic      Y W\n");
                printf("  Total Aliphatic           L I V M A P\n");
                printf("  Branched Chain Aliphatic  L I V\n\n");
                printf("Usage: chkall ws rf wf1 wf2.\n");
                printf("    ws = window size-must be odd\n");
                printf("    rf = alignment in NEWAT format\n");
                printf("    wf1= save file for various group frequencies\n");
                printf("    wf2= save file frequency of amino acid
appearance\n");
                printf("    MAX length = 1000      MAX SEQ=30\n");
                exit (1);
                }
        if((ALIGN=fopen(av[2],"r"))==NULL)
                {
                printf("File %s does not exist.\n", av[2]);
                exit (1);
                }
        ALLout=fopen(av[3],"w");
        AAFout=fopen(av[4],"w");

/* Determine 2-D scorearray1 size */
        numlines=0;
        fgets(str,sizeof str,ALIGN);
        while (strstr(str,"*")==NULL)
                {
               fgets(str,sizeof str,ALIGN);
                numlines++;
                }
        length=(numlines)*30;

/* Read aligned sequences into the array seqarray */

        rewind(ALIGN);
        row=0;
        numseq=0;
        for(;;)
                {
                c=fgetc(ALIGN);
                if(c==EOF)
                        {
                        fclose(ALIGN);
                        break;
                        }
                else
                        {
                        ungetc(c,ALIGN);
                        fgets(str,sizeof str,ALIGN);
                        if(str[0]=='T')
                                {
                                numseq++;
                                y=0;
                                for(i=0;i<numlines;i++)
                                        {
                                        fgets(str,sizeof str,ALIGN);
                                        for(x=11;x<70;x=x+2)
                                                {
                                                seqarray[row][y]=str[x];
                                                y++;
                                                }
                                        }
                                row++;
                                }
                        }
                }
        fprintf(ALLout,"number of sequences is %d\n",numseq);

/* Determine protein's alignment length */

        x=0;
        while(seqarray[0][x]!='\*')
                x++;
        algn_leng=x;


/* Perform the tally and store data in scorearray */
/* 0th row = aliphatic percent */
/* 1st row = aromatic percent */
/* 2nd row = other */
/* scorearray1's horizontal positions coincide with alignment position */
/* aafreq stores the frequency of each amino acid's occurrence at an alignment
position */

        fprintf(ALLout,"aligned sequence length is %d\n",algn_leng);
        for(pos=0;pos<algn_leng;pos++)
                {
                str_polar=0;
                anionic=0;
                cationic=0;
                semipolar=0;
                hydroxy=0;
                tot_hydrophob=0;
                amphi_arom=0;
                tot_arom=0;
                tot_aliph=0;
                br_chain_aliph=0;
                Acnt=0;
                Rcnt=0;
                Ncnt=0;
                Dcnt=0;
                Gcnt=0;
                Ccnt=0;
                Qcnt=0;
                Ecnt=0;
                Hcnt=0;
                Icnt=0;
                Lcnt=0;
                Kcnt=0;
                Mcnt=0;
                Fcnt=0;
                Pcnt=0;
                Scnt=0;
                Tcnt=0;
                Wcnt=0;
                Ycnt=0;
                Vcnt=0;
                for(e=0;e<numseq;e++)
                        {
                        if(seqarray[e][pos]=='D' || seqarray[e][pos]=='E' ||
seqarray[e][pos]=='N' || seqarray[e][pos]=='Q' || seqarray[e][pos]=='K' ||
seqarray[e][pos]=='R' || seqarray[e][pos]=='H')
                                str_polar++;
                        if(seqarray[e][pos]=='D' || seqarray[e][pos]=='E')
                                anionic++;
                        if(seqarray[e][pos]=='R' || seqarray[e][pos]=='K' ||
seqarray[e][pos]=='H')
                                cationic++;
                        if(seqarray[e][pos]=='G' || seqarray[e][pos]=='S' ||
seqarray[e][pos]=='T' || seqarray[e][pos]=='C')
                                semipolar++;
                        if(seqarray[e][pos]=='S' || seqarray[e][pos]=='T')
                                hydroxy++;
                        if(seqarray[e][pos]=='L' || seqarray[e][pos]=='I' ||
seqarray[e][pos]=='V' || seqarray[e][pos]=='M' || seqarray[e][pos]=='F' ||
seqarray[e][pos]=='Y' || seqarray[e][pos]=='W' || seqarray[e][pos]=='A' ||
seqarray[e][pos]=='P')
                                tot_hydrophob++;
                        if(seqarray[e][pos]=='Y' || seqarray[e][pos]=='W')
                                amphi_arom++;
                        if(seqarray[e][pos]=='Y' || seqarray[e][pos]=='W' ||
seqarray[e][pos]=='F')
                                tot_arom++;
                        if(seqarray[e][pos]=='L' || seqarray[e][pos]=='I' ||
seqarray[e][pos]=='V' || seqarray[e][pos]=='M' || seqarray[e][pos]=='A' ||
seqarray[e][pos]=='P')
                                tot_aliph++;
                        if(seqarray[e][pos]=='L' || seqarray[e][pos]=='I' ||
seqarray[e][pos]=='V')
                                br_chain_aliph++;


                        if(seqarray[e][pos]=='A')
                                Acnt++;
                        else if(seqarray[e][pos]=='R')
                                Rcnt++;
                        else if(seqarray[e][pos]=='N')
                                Ncnt++;
                       else if(seqarray[e][pos]=='D')
                                Dcnt++;
                        else if(seqarray[e][pos]=='G')
                                Gcnt++;
                        else if(seqarray[e][pos]=='C')
                                Ccnt++;
                        else if(seqarray[e][pos]=='Q')
                                Qcnt++;
                        else if(seqarray[e][pos]=='E')
                                Ecnt++;
                        else if(seqarray[e][pos]=='H')
                                Hcnt++;
                        else if(seqarray[e][pos]=='I')
                                Icnt++;
                        else if(seqarray[e][pos]=='L')
                                Lcnt++;
                        else if(seqarray[e][pos]=='K')
                                Kcnt++;
                        else if(seqarray[e][pos]=='M')
                                Mcnt++;
                        else if(seqarray[e][pos]=='F')
                                Fcnt++;
                        else if(seqarray[e][pos]=='P')
                                Pcnt++;
                        else if(seqarray[e][pos]=='S')
                                Scnt++;
                        else if(seqarray[e][pos]=='T')
                                Tcnt++;
                        else if(seqarray[e][pos]=='W')
                                Wcnt++;
                        else if(seqarray[e][pos]=='Y')
                                Ycnt++;
                        else if(seqarray[e][pos]=='V')
                                Vcnt++;
                        }
                scorearray1[0][pos]=str_polar/numseq;
                scorearray1[1][pos]=anionic/numseq;
                scorearray1[2][pos]=cationic/numseq;
                scorearray1[3][pos]=semipolar/numseq;
                scorearray1[4][pos]=hydroxy/numseq;
                scorearray1[5][pos]=tot_hydrophob/numseq;
                scorearray1[6][pos]=amphi_arom/numseq;
                scorearray1[7][pos]=tot_arom/numseq;
                scorearray1[8][pos]=tot_aliph/numseq;
                scorearray1[9][pos]=br_chain_aliph/numseq;
                aafreq[0][pos]=Acnt/numseq;
                aafreq[1][pos]=Rcnt/numseq;
                aafreq[2][pos]=Ncnt/numseq;
                aafreq[3][pos]=Dcnt/numseq;
                aafreq[4][pos]=Gcnt/numseq;
                aafreq[5][pos]=Ccnt/numseq;
                aafreq[6][pos]=Qcnt/numseq;
                aafreq[7][pos]=Ecnt/numseq;
                aafreq[8][pos]=Hcnt/numseq;
                aafreq[9][pos]=Icnt/numseq;
                aafreq[10][pos]=Lcnt/numseq;
                aafreq[11][pos]=Kcnt/numseq;
                aafreq[12][pos]=Mcnt/numseq;
                aafreq[13][pos]=Fcnt/numseq;
                aafreq[14][pos]=Pcnt/numseq;
                aafreq[15][pos]=Scnt/numseq;
                aafreq[16][pos]=Tcnt/numseq;
                aafreq[17][pos]=Wcnt/numseq;
                aafreq[18][pos]=Ycnt/numseq;
                aafreq[19][pos]=Vcnt/numseq;
              }

/* Change counts as determined by window size */


        fprintf(ALLout,"window size is %d\n",wsize);
        ansarrleng=algn_leng-wsize+1;
        start=0;
        for(z=0;z<10;z++) /* Category counter */
                {
                for(x=0;x<ansarrleng;x++) /* pos'n counter */
                        {
                        sum=0;
                        for(y=start;y<start+wsize;y++) /* Averager */
                                sum=sum+scorearray1[z][y];
                        scorearray2[z][x]=(sum/wsize);
                        start++;
                        }
                start=0;
                }

/* Output aliphatic and aromatic data for Excel spreadsheet usage */

        stop=ansarrleng+5;
        fprintf(ALLout,"Data from 6 to %d\n",stop);
        fprintf(ALLout,"Position,Str
Polar,Anionic,Cationic,Semipolar,Hydroxy,Tot Hydrophob,Amphi Arom,Tot Arom,Tot
Aliph,Br Chain Aliph\n");
        begin=(wsize+1)/2;
        posn=begin;
        for (x=0;x<ansarrleng;x++)
                {

fprintf(ALLout,"%d,%1.2f,%1.2f,%1.2f,%1.2f,%1.2f,%1.2f,%1.2f,%1.2f,%1.2f,%1.2f\n
",posn,scorearray2[0][x],scorearray2[1][x],scorearray2[2][x],scorearray2[3][x],s
corearray2[4][x],scorearray2[5][x],scorearray2[6][x],scorearray2[7][x],scorearra
y2[8][x],scorearray2[9][x]);
                posn++;
                }

/* Output amino acid frequencies in Excel usable format */

        fprintf(AAFout,"number of sequences is %d\n",numseq);
        fprintf(AAFout,"alignment length is %d\n",algn_leng);
        fprintf(AAFout,"Position,A,R,N,D,G,C,Q,E,H,I,L,K,M,F,P,S,T,W,Y,V\n");
        posn=1;
        for (x=0;x<algn_leng;x++)
                {

fprintf(AAFout,"%d,%1.2f,%1.2f,%1.2f,%1.2f,%1.2f,%1.2f,%1.2f,%1.2f,%1.2f,%1.2f,%
1.2f,%1.2f,%1.2f,%1.2f,%1.2f,%1.2f,%1.2f,%1.2f,%1.2f,%1.2f\n",posn,aafreq[0][x],
aafreq[1][x],aafreq[2][x],aafreq[3][x],aafreq[4][x],aafreq[5][x],aafreq[6][x],aa
freq[7][x],aafreq[8][x],aafreq[9][x],aafreq[10][x],aafreq[11][x],aafreq[12][x],a
afreq[13][x],aafreq[14][x],aafreq[15][x],aafreq[16][x],aafreq[17][x],aafreq[18][
x],aafreq[19][x]);
                posn++;
                }


        fclose(ALLout);
        fclose(AAFout);
}