/*  This program counts the number of aromatic and aliphatic residues at a */
/*  given position in an alignment.  The percentages are given as averages */
/*  in a variable window.  This program is different than chkALAR2_1 in    */
/*  that data is also gathered for each alignment position.     8/28/94 st */

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

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

{
        FILE *ALIGN,*ALARout,*AAFout;
        char code[6],seqarray[30][1000],str[200];
        int
ansarrleng,begin,c,e,end,i,length,numlines,numseq,pos,posn,rleng,row,start,stop,
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
ALcount,aafreq[29][1000],ARcount,Other,scorearray1[3][1000],scorearray2[3][1000]
,sum;

        wsize=atoi(av[1]);
        if(ac!=5 && wsize%2==0)
                {
                printf("chkA gathers statistical data on all the positions in
an alignment.\n");
                printf("Usage: chkA 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 aliphatic and aromatic
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);
                }
        ALARout=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(ALARout,"number of sequences is %d\n",numseq);

/* Determine protein's alignment length */

        x=0;
        while(seqarray[0][x]!='\*')
                x++;
        rleng=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(ALARout,"aligned sequence length is %d\n",rleng);
        for(pos=0;pos<rleng;pos++)
                {
                ARcount=0;
                ALcount=0;
                Other=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]=='A' || seqarray[e][pos]=='V' ||
seqarray[e][pos]=='I' || seqarray[e][pos]=='L' || seqarray[e][pos]=='G')
                        ALcount++;
                        else if(seqarray[e][pos]=='F'  || seqarray[e][pos]=='Y'
|| seqarray[e][pos]=='W')
                        ARcount++;
                        else Other++;

                        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]=ALcount/numseq;
                scorearray1[1][pos]=ARcount/numseq;
                scorearray1[2][pos]=Other/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(ALARout,"window size is %d\n",wsize);
        ansarrleng=rleng-wsize+1;
        start=0;
        for(z=0;z<3;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(ALARout,"Data from 6 to %d\n",stop);
        fprintf(ALARout,"Position,Aliphatic,Aromatic,Other\n");
        begin=(wsize+1)/2;
        posn=begin;
        for (x=0;x<ansarrleng;x++)
                {

fprintf(ALARout,"%d,%1.2f,%1.2f,%1.2f\n",posn,scorearray2[0][x],scorearray2[1][x
],scorearray2[2][x]);
                posn++;
                }

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

        fprintf(AAFout,"numseq is %d\n",numseq);
        fprintf(AAFout,"alignment length is %d\n",rleng);
        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<rleng;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(ALARout);
        fclose(AAFout);
}