rendered paste body/*Secondary Structure Prediction Program SIMPA96 */
/* Author Jonathan Levin E-Mail jon@glycine.jouy.inra.fr*/
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#define mxl 200
#define totres 150000
#define mxres 1500
#define mxseq 500
#define seg 15
#define cut 13
#define conf "Hb."
#define h 4
#define b 2
#define nchar 60
char dbs[] ="dbks.seq";
char dbo[] ="dbks.obs";
int readseq (char *, char [],char [][70], int [], int []);
void find (char [],char [],int [],int ,int [][3],
char [][70],char *, int *[],int ,int []);
void regul(float (*)[3],int, char *);
float big(float *);
int ibig(float *);
float max(float,float);
void fmax(char *, int *, float *);
void consensus (char [][mxseq],int ,int ,float (*)[3],int *,
char *,char *);
void HH (char [][mxseq],int,int,float *);
void pseudo(char [][mxseq],int,int,int *, char *);
void sipred (char *,char *,int,char *,int *, int, FILE *);
int mat[676];
float dc[3]={.86,1.10,.81};
int main (void)
{
FILE *fi1,*fi2,*fi3;
int dird=0,i,j,k,kk,dnseq,dnres[mxseq],dcres[mxseq],tnseq,tnres[mxseq],
tcres[mxseq],tusf,ali,l,ii,ncl,*ptseq[mxres],chk[mxres][3],mc,
nrts[mxseq],scor[mxres],scs[mxres];
float (*prop)[3],tprop[mxres][3],uprop[mxres][3],tmp[3];
char dseq[totres],dobs[totres],dnames[mxseq][70],usf[70],
tseq[totres],tnames[mxseq][70],yes[3],ucf[70],line[mxl],
clus[mxres][mxseq],prs[mxres],rpr[mxres*mxseq],mpr[mxres][mxseq],
pr[mxres],lpr[mxres],file3[70],uof[70];
char file1[70]="./";
char file2[70]="/usr/local/lib/";
/* data entry */
if((fi1=fopen("./mat.data","r"))== NULL)
{
if((fi1=fopen("/usr/local/lib/mat.data","r"))== NULL)
{
printf(" can't open mat.data file from . or /usr/local/lib\n");
exit(1);
}
dird=1;
}
for(i=0;i<26;i++)
{
for(k=0;k<26;k++)fscanf(fi1,"%2d",&mat[i*26+k]);
while(getc(fi1)!='\n');
}
fclose(fi1);
if(dird)
{
strcpy(file3,file2);
strcat(file3,dbs);
dnseq=readseq(file3,dseq,dnames,dnres,dcres);
strcpy(file3,file2);
strcat(file3,dbo);
dnseq=readseq(file3,dobs,dnames,dnres,dcres);
}
else
{
strcpy(file3,file1);
strcat(file3,dbs);
dnseq=readseq(file3,dseq,dnames,dnres,dcres);
strcpy(file3,file1);
strcat(file3,dbo);
dnseq=readseq(file3,dobs,dnames,dnres,dcres);
}
/* enter user files */
printf(" enter sequence file name (PIR,fasta format)\n");
gets(usf);
tnseq=readseq(usf,tseq,tnames,tnres,tcres);
if(tnseq<1)
{
printf(" no sequences found in file %s\n",usf);
exit(1);
}
else
for(i=0;i<tnseq;i++)
if(tnres[i]<20)
{
printf(" sequence %s contains less than 20 AA\n",tnames[i]);
exit(1);
}
printf(" %d sequences found in file %s\n",tnseq,usf);
tusf=1;
while(tusf)
{
printf(" is this correct? (y/n)\n");
gets(yes);
if(yes[0]=='y')tusf=0;
else if(yes[0]=='n')exit(1);
}
yes[0]='\0';
printf(" is there a clustal alignment file? (y/n)\n");
gets(yes);
ali=0;
if(yes[0]=='y')
{
printf(" enter clustal file name \n");
gets(ucf);
ali=1;
if((fi2=fopen(ucf,"r"))== NULL)
{
printf(" can't open clustal file: %s\n",ucf);
exit(1);
}
/* reads clustal file */
j=0;k=0;
for(i=0;i<199;i++)line[i]=' ';
while((line[17]<'A'||line[17]>'Z')&&line[17]!='-')
if(fgets(line,mxl,fi2)==NULL)
{
printf("Error reading clustal file %s\n",ucf);
printf("Are you sure this file is ");
printf("in CLUSTAL W(1.60) format ?\n");
exit(1);
}
else k++;
while(line[17]>='A'&&line[17]<='Z'||line[17]=='-')
{
line[17]=' ';
if(fgets(line,mxl,fi2)==NULL)
{
printf("Error reading clustal file %s\n",ucf);
printf("Are you sure this file is ");
printf("in CLUSTAL W(1.60) format ?\n");
exit(1);
}
j++;
}
if(j!=tnseq)
{
printf("ERROR no of seqs in seq file ");
printf("%d no of seq clustal file %d\n",tnseq,j);
exit(1);
}
rewind(fi2);
for(i=0;i<k;i++)
if(fgets(line,mxl,fi2)==NULL)
{
printf("Error reading clustal file %s\n",ucf);
printf("Are you sure this file is ");
printf("in CLUSTAL W(1.60) format ?\n");
exit(1);
}
l=0;ii=0;
while(l==0)
{
for(i=0;i<tnseq;i++)
{
j=ii*60;
k=16;
while(line[k]!='\n'){clus[j+k-16][i]=line[k];k++;}
if(fgets(line,mxl,fi2)==NULL)
{
printf("Error reading clustal file %s\n",ucf);
printf("Are you sure this file is ");
printf("in CLUSTAL W(1.60) format ?\n");
exit(1);
}
}
if(fgets(line,mxl,fi2)==NULL)l=1;
line[16]=' ';
while((line[16]<'A'||line[16]>'Z')&&line[16]!='-')
if(fgets(line,mxl,fi2)==NULL){l=1;line[16]='-';}
ii++;
ncl=j+k-16;
}
fclose(fi2);
/* check to see if same number of residues in seq file and clus file */
for(i=0;i<tnseq;i++)
{
k=tcres[i];
for(j=0;j<ncl;j++)
if(clus[j][i]!='-')
if(clus[j][i]==tseq[k])k++;
else
{
printf("In file %s sequence no %3d res no. %4d '%c'\nis ",
usf,i+1,k-tcres[i]+1,tseq[k]);
printf("different from\nfile");
printf(" %s sequence no %3d res no. %4d '%c'\n",ucf,i+1,j+1,
clus[j][i]);
exit(1);
}
}
}
printf("enter prediction ouput file name\n");
gets(uof);
if((fi3=fopen(uof,"w"))== NULL)
{
printf("can't open output file: %s\n",uof);
exit(1);
}
else printf(" opening prediction ouput file: %s\n", uof);
fprintf(fi3,"SIMPA96 secondary structure prediction program\n\n");
fprintf(fi3,"H = helix b = strand . = other\n\n");
fprintf(fi3,"Scale indicates prediction strength: ");
fprintf(fi3,"0-weak, 2-medium, 5-strong 9-very strong\n\n");
/* do prediction */
prop=calloc(tcres[tnseq]*3,sizeof(float));
for(ii=0;ii<tnseq;ii++)
{
for(i=0;i<tnres[ii];i++)ptseq[i]=&mat[(tseq[tcres[ii]+i]-65)*26];
find (dseq,dobs,dnres,dnseq,chk,dnames,tnames[ii],ptseq,
tnres[ii],dcres);
for(i=0;i<tnres[ii];i++)
for(k=0;k<3;k++)prop[tcres[ii]+i][k]=chk[i][k]*dc[k];
printf("%s\n",tnames[ii]);
}
/* regularize structure or do multiple alignment prediction */
if(ali)
{
kk=0,mc=0;
for(i=0;i<tnseq;i++)
{
for(k=0;k<tnres[i];k++)for(j=0;j<3;j++)
tprop[k][j]=prop[mc+k][j];
regul(tprop,tnres[i],prs);
for(k=0;k<tnres[i];k++)rpr[mc+k]=prs[k];
mc+=tnres[i];
}
for(i=0;i<ncl;i++)for(k=0;k<3;k++)uprop[i][k]=0;
nrts[0]=0;k=0;
for(i=0;i<tnseq;i++)
{
k+=tnres[i];
nrts[i+1]=k;
}
for(i=0;i<ncl;i++)
{
for(j=0;j<tnseq;j++)
if(clus[i][j]!='-')
{
for(k=0;k<3;k++)uprop[i][k]+=prop[nrts[j]][k];
mpr[i][j]=rpr[nrts[j]];
nrts[j]++;
}
else mpr[i][j]='-';
for(k=0;k<3;k++)tmp[k]=uprop[i][k];
fmax(&pr[i],&scor[i],tmp);
}
consensus(clus,tnseq,ncl,prop,tnres,pr,lpr);
for(i=0;i<ncl;i++)
if(clus[i][0]!='-')
{
prs[kk]=pr[i];
scs[kk]=scor[i];
kk++;
}
for(i=0;i<tnseq;i++)
for(j=0;j<25;j++)
if(tnames[i][j]=='\n'||tnames[i][j]=='\0')tnames[i][j]=' ';
fprintf(fi3,"multiple sequence prediction\n\n");
kk=ncl/nchar;
for(k=0;k<kk;k++)
{
for(i=0;i<19;i++)putc(' ',fi3);
for(i=k*nchar+10;i<=k*nchar+nchar;i=i+10)fprintf(fi3,"%10d",i);
putc('\n',fi3);
for(j=0;j<16;j++)putc(tnames[0][j],fi3);fprintf(fi3," ");
for(j=k*nchar;j<k*nchar+nchar;j++)putc(clus[j][0],fi3);
putc('\n',fi3);
for(i=0;i<tnseq;i++)
{
for(j=0;j<16;j++)putc(tnames[i][j],fi3);
fprintf(fi3," ");
for(j=k*nchar;j<k*nchar+nchar;j++)putc(mpr[j][i],fi3);
putc('\n',fi3);
}
fprintf(fi3,"consensus pred ");
for(j=k*nchar;j<k*nchar+nchar;j++)putc(pr[j],fi3);
putc('\n',fi3);
fprintf(fi3,"scale ");
for(j=k*nchar;j<k*nchar+nchar;j++)fprintf(fi3,"%1d",scor[j]);
putc('\n',fi3);putc('\n',fi3);putc('\n',fi3);
}
if(kk*nchar!=ncl)
{
for(i=0;i<19;i++)putc(' ',fi3);
for(i=kk*nchar+10;i<=ncl;i=i+10)fprintf(fi3,"%10d",i);
putc('\n',fi3);
for(j=0;j<16;j++)putc(tnames[0][j],fi3);fprintf(fi3," ");
for(j=kk*nchar;j<ncl;j++)putc(clus[j][0],fi3);putc('\n',fi3);
for(i=0;i<tnseq;i++)
{
for(j=0;j<16;j++)putc(tnames[i][j],fi3);fprintf(fi3," ");
for(j=kk*nchar;j<ncl;j++)putc(mpr[j][i],fi3);
putc('\n',fi3);
}
fprintf(fi3,"consensus pred ");
for(j=kk*nchar;j<ncl;j++)putc(pr[j],fi3);putc('\n',fi3);
fprintf(fi3,"scale ");
for(j=kk*nchar;j<ncl;j++)fprintf(fi3,"%1d",scor[j]);
putc('\n',fi3);putc('\n',fi3);putc('\n',fi3);
fprintf(fi3,"consensus prediction\n");
}
sipred(tnames[0],tseq,tcres[0],lpr,scs,tnres[0],fi3);
}
else
{
ii=0;mc=0;
fprintf(fi3,"secondary structure prediction without homologues\n\n");
for(i=0;i<tnseq;i++)
{
for(k=0;k<tnres[i];k++)
{
for(j=0;j<3;j++)tmp[j]=prop[ii][j];
ii++;
fmax(&prs[k],&scs[k],tmp);
}
for(k=0;k<tnres[i];k++)
for(j=0;j<3;j++)tprop[k][j]=prop[mc+k][j];
mc+=tnres[i];
regul(tprop,tnres[i],prs);
sipred(tnames[i],tseq,tcres[i],prs,scs,tnres[i],fi3);
}
}
fclose(fi3);
}
/* function sipred prints output seq,pred,scale */
void sipred (char *name,char *seq,int cres,char *pr,int *scale,int nres,
FILE *fi3)
{
int j,k,kk;
fprintf(fi3,"\n%s\n\n",name);
kk=nres/nchar;
for(k=0;k<kk;k++)
{
for(j=0;j<12;j++)putc(' ',fi3);
for(j=k*nchar+10;j<=k*nchar+nchar;j=j+10)fprintf(fi3,"%10d",j);
putc('\n',fi3);
fprintf(fi3,"sequence ");
for(j=k*nchar;j<k*nchar+nchar;j++)putc(seq[j+cres],fi3);
putc('\n',fi3);
fprintf(fi3,"prediction ");
for(j=k*nchar;j<k*nchar+nchar;j++)putc(pr[j],fi3);
putc('\n',fi3);
fprintf(fi3,"scale ");
for(j=k*nchar;j<k*nchar+nchar;j++)fprintf(fi3,"%1d",scale[j]);
putc('\n',fi3);putc('\n',fi3);putc('\n',fi3);
}
if(kk*nchar!=nres)
{
for(j=0;j<12;j++)putc(' ',fi3);
for(j=kk*nchar+10;j<=nres;j=j+10)fprintf(fi3,"%10d",j);
putc('\n',fi3);
fprintf(fi3,"sequence ");
for(j=kk*nchar;j<nres;j++)putc(seq[j+cres],fi3);
putc('\n',fi3);
fprintf(fi3,"prediction ");
for(j=kk*nchar;j<nres;j++)putc(pr[j],fi3);putc('\n',fi3);
fprintf(fi3,"scale ");
for(j=kk*nchar;j<nres;j++)fprintf(fi3,"%1d",scale[j]);
putc('\n',fi3);putc('\n',fi3);
}
}
/* function consensus tries to make a consensus prediction */
void consensus (char clus[][mxseq],int nseq,int ncl,float prop[][3],
int *nres,char *pr,char *lpr)
{
char mpr[mxres],psu[mxres];
float tprop[mxres][3],mprop[mxres][3],tmp[3],x,wei[mxseq];
int i,j,k,l,nrts[mxres],maxp[mxres][3],norm[mxres],ici[26],t;
HH(clus,nseq,ncl,wei);
pseudo(clus,nseq,ncl,nres,psu);
k=0;nrts[0]=0;
for(i=0;i<nseq;i++)
{
k+=nres[i];
nrts[i+1]=k;
}
for(i=0;i<ncl;i++)for(j=0;j<3;j++)tprop[i][j]=0;
for(i=0;i<ncl;i++)for(j=0;j<3;j++)maxp[i][j]=0;
for(i=0;i<ncl;i++)norm[i]=0;
for(i=0;i<ncl;i++)
{
for(j=0;j<26;j++)ici[j]=0;t=0;
for(j=0;j<nseq;j++)
if(clus[i][j]!='-')ici[clus[i][j]-65]++;
for(j=0;j<26;j++)if(ici[j]!=0)t++;
for(j=0;j<nseq;j++)
if(clus[i][j]!='-')
{
x=0;
for(k=0;k<3;k++)tmp[k]=prop[nrts[j]][k];
k=ibig(tmp);
maxp[i][k]++;
for(k=0;k<3;k++)x+=tmp[k];
for(k=0;k<3;k++)tmp[k]=tmp[k]/x;
tmp[1]=tmp[1]*1.04;
for(k=0;k<3;k++)tprop[i][k]+=tmp[k]*wei[j];
nrts[j]++;
norm[i]++;
}
}
for(i=0;i<ncl;i++)for(j=0;j<3;j++)
tprop[i][j]=tprop[i][j]/(float)norm[i];
j=0;
for(i=0;i<ncl;i++)
if(clus[i][0]!='-')
{
for(k=0;k<3;k++)mprop[j][k]=tprop[i][k];
j++;
}
regul(tprop, ncl, pr);
regul(mprop, nres[0], lpr);
}
void HH (char clus[][mxseq],int nseq,int ncl,float *wei)
{
int i,j,ici[26],t;
float u;
for(i=0;i<ncl;i++)wei[i]=0;
for(i=0;i<ncl;i++)
{
for(j=0;j<26;j++)ici[j]=0;t=0;
for(j=0;j<nseq;j++)
if(clus[i][j]!='-')ici[clus[i][j]-65]++;
for(j=0;j<26;j++)if(ici[j]!=0)t++;
for(j=0;j<nseq;j++)
if(clus[i][j]!='-')wei[j]+=1/(ici[clus[i][j]-65]*(float)t);
}
u=0;
for(i=0;i<nseq;i++)u+=wei[i];
for(i=0;i<nseq;i++)wei[i]=wei[i]/u;
}
void pseudo(char clus[][mxseq],int nseq,int ncl,int *nres, char *psu)
{
int i, j, nco[mxseq];
for(i=0;i<nseq;i++)nco[i]=0;
for(i=0;i<ncl;i++)psu[i]='X';
for(i=0;i<ncl;i++)
{
for(j=0;j<nseq;j++)
{
if(clus[i][j]=='-')
{
if(nco[j]!=0&&nco[j]<nres[j])psu[i]='-';
}
else nco[j]++;
}
}
}
/* calculates scale values and chosen conformation */
void fmax(char *c, int *sc, float tmp[])
{
int i,j;
float t=0,u=0;
for (i=0;i<3;i++)
if(tmp[i]>=t)
{
t=tmp[i];
j=i;
}
if(j==0)*c='H';
else if(j==1)*c='b';
else if(j==2)*c='.';
tmp[j]=0;
for (i=0;i<3;i++)
if(tmp[i]>=u)
{
u=tmp[i];
j=i;
}
*sc=(1-u/t)*10;
if(*sc>9)*sc=9;
}
/* regularisation using dynamic programming */
void regul(float (*prop)[3],int nres, char pr[])
{
float dyn[mxres][3],H,B;
int i,j,k,l;
for(i=0;i<nres;i++)for(j=0;j<3;j++)dyn[i][j]=0;
/* start filling dyn first and last residues will be coil*/
for(i=1;i<h+1;i++)dyn[i][0]=0;
for(i=1;i<b+1;i++)dyn[i][1]=0;
for(i=1;i<h+1;i++)dyn[h][0]+=prop[i][0];
for(i=1;i<b+1;i++)dyn[b][1]+=prop[i][1];
dyn[1][2]=prop[1][2];
if(h>=b)
{
for(i=2;i<b+1;i++)dyn[i][2]=dyn[i-1][2]+prop[i][2];
for(i=b+1;i<h+1;i++)
{
B=0;
for(j=i-b+1;j<=i;j++)B+=prop[j][1];
dyn[i][1]=max(big(&dyn[i-b][0])+B,dyn[i-1][1]+prop[i][1]);
dyn[i][2]=big(&dyn[i-1][0])+prop[i][2];
}
}
else
{
for(i=2;i<h+1;i++)dyn[i][2]=dyn[i-1][2]+prop[i][2];
for(i=h+1;i<b+1;i++)
{
H=0;
for(j=i-h+1;j<=i;j++)H+=prop[j][0];
dyn[i][0]=max(big(&dyn[i-h][0])+H,dyn[i-1][0]+prop[i][0]);
dyn[i][2]=big(&dyn[i-1][0])+prop[i][2];
}
}
/* fill up dyn from max(h,b) to end */
l=(h>=b)?h:b;
for(i=l+1;i<nres-1;i++)
{
H=0;
for(j=i-h+1;j<=i;j++)H+=prop[j][0];
dyn[i][0]=max(big(&dyn[i-h][0])+H,dyn[i-1][0]+prop[i][0]);
B=0;
for(j=i-b+1;j<=i;j++)B+=prop[j][1];
dyn[i][1]=max(big(&dyn[i-b][0])+B,dyn[i-1][1]+prop[i][1]);
dyn[i][2]=big(&dyn[i-1][0])+prop[i][2];
}
/* go back through dyn to choose best path */
i=nres-2;
k=ibig(&dyn[i][0]);
pr[i]=conf[k];
while(i>0)
{
if(i<0)
{
printf(" error at end of dyn **** i<0 ");
exit(1);
}
i--;
j=ibig(&dyn[i][0]);
if(j==k)pr[i]=conf[j];
else if(k==0)
{
for(l=i;l>i-h+1;l--)pr[l]=conf[k];
i=i-h+1;
k=ibig(&dyn[i][0]);
pr[i]=conf[k];
}
else if(k==1)
{
for(l=i;l>i-b+1;l--)pr[l]=conf[k];
i=i-b+1;
k=ibig(&dyn[i][0]);
pr[i]=conf[k];
}
else if(k==2)
{
k=ibig(&dyn[i][0]);
pr[i]=conf[k];
}
}
pr[0]='.';
pr[nres-1]='.';
}
/* function big returns largest value in row of dyn */
float big(float dyn[])
{
int i;
float a=0;
for(i=0;i<3;i++)if(a<dyn[i])a=dyn[i];
return a;
}
/* function ibig returns index of biggest value */
int ibig(float dyn[])
{
int i,j=0;
float a=0;
for(i=0;i<3;i++)
if(a<dyn[i]){a=dyn[i];j=i;}
return j;
}
/* function max returns biggest value */
float max(float x,float z)
{
return(x<z)?z:x;
}
/* fills up propensity matrix */
void find (char dseq[],char dobs[],int dnres[],int dnseq,int chk[][3],
char dnames[][70],char *tname, int *tseq[],int tres,int dcres[])
{
int ii,i,j,ik,k,iseq[mxres],iobs[mxres],tsiz,dsiz,ijk,it;
for(i=0;i<mxres;i++)for(j=0;j<3;j++)chk[i][j]=0;
for(ii=0;ii<dnseq;ii++)
{
for(i=0;i<dnres[ii];i++)
{
iseq[i]=dseq[dcres[ii]+i]-65;
if(dobs[dcres[ii]+i]=='H')iobs[i]=0;
else if(dobs[dcres[ii]+i]=='E')iobs[i]=1;
else if(dobs[dcres[ii]+i]=='C')iobs[i]=2;
}
if(tres<=dnres[ii])
{
tsiz=tres-seg-1;
dsiz=dnres[ii]-seg-1;
for(ik=0;ik<=dsiz;ik++)
{
tsiz++;
if(tsiz>=dsiz)tsiz=dsiz;
ijk=0;it=0;
for(k=0;k<seg;k++)ijk+=*(tseq[(it+k)]+iseq[ik+k]);
if(ijk>=cut)
for(k=0;k<seg;k++)chk[it+k][iobs[ik+k]]+=ijk;
for(j=ik;j<tsiz;j++)
{
ijk=ijk-*(tseq[it]+iseq[j])
+*(tseq[it+seg]+iseq[j+seg]);
it++;
if(ijk>=cut)
for(k=0;k<seg;k++)chk[it+k][iobs[j+1+k]]+=ijk;
}
}
tsiz=tres-seg;
for(ik=1;ik<tsiz;ik++)
{
ijk=0;it=0;
for(k=0;k<seg;k++)ijk+=*(tseq[(ik+k)]+iseq[it+k]);
if(ijk>=cut)
for(k=0;k<seg;k++)chk[ik+k][iobs[it+k]]+=ijk;
for(j=ik;j<tsiz;j++)
{
ijk=ijk-*(tseq[j]+iseq[it])
+*(tseq[j+seg]+iseq[it+seg]);
it++;
if(ijk>=cut)
for(k=0;k<seg;k++)chk[j+1+k][iobs[it+k]]+=ijk;
}
}
}
else
{
tsiz=tres-seg-1;
dsiz=dnres[ii]-seg-1;
for(ik=0;ik<=tsiz;ik++)
{
dsiz++;
if(dsiz>=tsiz)dsiz=tsiz;
ijk=0;it=0;
for(k=0;k<seg;k++)ijk+=*(tseq[(ik+k)]+iseq[it+k]);
if(ijk>=cut)
for(k=0;k<seg;k++)chk[ik+k][iobs[it+k]]+=ijk;
for(j=ik;j<dsiz;j++)
{
ijk=ijk-*(tseq[j]+iseq[it])
+*(tseq[j+seg]+iseq[it+seg]);
it++;
if(ijk>=cut)
for(k=0;k<seg;k++)chk[j+1+k][iobs[it+k]]+=ijk;
}
}
dsiz=dnres[ii]-seg;
for(ik=1;ik<dsiz;ik++)
{
ijk=0;it=0;
for(k=0;k<seg;k++)ijk+=*(tseq[(it+k)]+iseq[ik+k]);
if(ijk>=cut)
for(k=0;k<seg;k++)chk[it+k][iobs[ik+k]]+=ijk;
for(j=ik;j<dsiz;j++)
{
ijk=ijk-*(tseq[it]+iseq[j])
+*(tseq[it+seg]+iseq[j+seg]);
it++;
if(ijk>=cut)
for(k=0;k<seg;k++)chk[it+k][iobs[j+1+k]]+=ijk;
}
}
}
}
return;
}
/* reads seq (or obs) files in pir or my format */
int readseq (char *file, char seq[],char names[][70],int nres[],int cres[])
{
FILE *fi1;
char line[mxl];
int i,j,k,nseq=0,cr=0;
cres[0]=0;
if((fi1=fopen(file,"r"))== NULL)
{
printf(" cant open file %s\n",file);
exit(1);
}
while (line[0]!='>'&&line[0]!='!')
if(fgets(line,mxl-1,fi1)==NULL)
{
printf("file %s is not in FastA format\n",file);
exit(1);
}
if(line[0]=='!')
{
rewind(fi1);
while (fgets(line,mxl-1,fi1)!=NULL)
{
if(line[0]=='!')
{
i=0;j=0,k=0;
while(i<68&&line[i+1]!='\n'){names[nseq][i]=line[i+1];i++;}
names[nseq][i]='\n';names[nseq][i+1]='\0';
while(k==0)
{
fgets(line,mxl-1,fi1);
i=0;
while(line[i]!='\n'&&i<199)
if(line[i]>='A'&&line[i]<='Z')
{
seq[cr+j++]=line[i];
i++;
}
else if (line[i]=='@')
{
nres[nseq]=j;
cr+=j;
cres[nseq+1]=cr;
nseq++;
k=1;
i=200;
}
else i++;
}
}
}
fclose(fi1);
return nseq;
}
else
{
rewind(fi1);
while (fgets(line,mxl-1,fi1)!=NULL)
{
if(line[0]=='>')
{
i=0;j=0,k=0;
while(i<68&&line[i+1]!='\n'){names[nseq][i]=line[i+1];i++;}
names[nseq][i]='\n';names[nseq][i+1]='\0';
while(k==0)
{
if(fgets(line,mxl-1,fi1)!=NULL)
{
if (line[0]=='>')
{
nres[nseq]=j;
cr+=j;
cres[nseq+1]=cr;
nseq++;
i=0;j=0,k=0;
while(i<68&&line[i+1]!='\n')
{names[nseq][i]=line[i+1];i++;}
names[nseq][i]='\n';names[nseq][i+1]='\0';
}
else
{
i=0;
while(line[i]!='\n'&&i<199)
if(line[i]>='A'&&line[i]<='Z')
{
seq[cr+j++]=line[i];
i++;
}
else i++;
}
}
else
{
nres[nseq]=j;
cr+=j;
cres[nseq+1]=cr;
nseq++;
k=1;
fclose(fi1);
return nseq;
}
}
}
}
}
}