/*  This program counts the number of residues determined by the user at a */
/*  given position in an alignment.  The percentages are given as averages */
/*  in a variable window. 8/24/94 st                                       */
#include <stdio.h>
#include <string.h>

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

{
        FILE *ALIGN,*INPout,*AAFout;
        char T,code[6],residues[12],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
Inputcnt,aafreq[29][1000],ARcount,Other,scorearray1[2][1000],scorearray2[2][1000
],sum;

        wsize=atoi(av[1]);
        if(ac!=5 && wsize%2==0)
                {
                printf("aachk gathers statistical data on all the positions in
an alignment.\nThe user is allowed to input a group (max. of 10) of amino
acids.\n");
                printf("Usage: aachk 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 input amino acid
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);
                }
        INPout=fopen(av[3],"w");
        AAFout=fopen(av[4],"w");

/* Get the amino acids from the user */

        do
                {
                printf("Enter the amino acids with a space between each one
(max. of 10): ");
                gets(residues);
                T=strlen(residues);
                }
        while(T%2==0 && T<20);
        fprintf(INPout,"input residues %s\n",residues);


/* Approximate the length of the sequences for use in the loop for filling
seqarray below */
        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.  This is the array of
sequences.  */

        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(INPout,"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 = input percent */
/* 1nd 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(INPout,"aligned sequence length is %d\n",algn_leng);
        for(pos=0;pos<algn_leng;pos++)
                {
                Inputcnt=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]==residues[0] ||
seqarray[e][pos]==residues[2] || seqarray[e][pos]==residues[4] ||
seqarray[e][pos]==residues[6] || seqarray[e][pos]==residues[8] ||
seqarray[e][pos]==residues[10] || seqarray[e][pos]==residues[12] ||
seqarray[e][pos]==residues[14] || seqarray[e][pos]==residues[16] ||
seqarray[e][pos]==residues[18])
                                Inputcnt++;
                        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]=Inputcnt/numseq;
                scorearray1[1][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(INPout,"window size is %d\n",wsize);
        ansarrleng=algn_leng-wsize+1;
        start=0;
        for(z=0;z<2;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;
                }

/* Format the input  data for Excel spreadsheet usage */

        stop=ansarrleng+6;
        fprintf(INPout,"Data from 7 to %d\n",stop);
        fprintf(INPout,"Position,Input,Other\n");
        begin=(wsize+1)/2;
        posn=begin;
        for (x=0;x<ansarrleng;x++)
                {

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

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

        fprintf(AAFout,"numseq 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(INPout);
        fclose(AAFout);
}