readmat.c

clustalx用来做基因序列分析

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include "clustalw.h"
#include "matrices.h"


/*
 *   Prototypes
 */
static Boolean commentline(char *line);


/*
 *   Global variables
 */

extern char 	*amino_acid_codes;
extern sint 	gap_pos1, gap_pos2;
extern sint 	max_aa;
extern short 	def_dna_xref[],def_aa_xref[];
extern sint 	mat_avscore;
extern sint 	debug;
extern Boolean  dnaflag;

extern Boolean user_series;
extern UserMatSeries matseries;
extern short usermatseries[MAXMAT][NUMRES][NUMRES];
extern short aa_xrefseries[MAXMAT][NUMRES+1];


void init_matrix(void)
{

   char c1,c2;
   short i, j, maxres;

   max_aa = strlen(amino_acid_codes)-2;
   gap_pos1 = NUMRES-2;          /* code for gaps inserted by clustalw */
   gap_pos2 = NUMRES-1;           /* code for gaps already in alignment */

/*
   set up cross-reference for default matrices hard-coded in matrices.h
*/
   for (i=0;i<NUMRES;i++) def_aa_xref[i] = -1;
   for (i=0;i<NUMRES;i++) def_dna_xref[i] = -1;

   maxres = 0;
   for (i=0;(c1=amino_acid_order[i]);i++)
     {
         for (j=0;(c2=amino_acid_codes[j]);j++)
          {
           if (c1 == c2)
               {
                  def_aa_xref[i] = j;
                  maxres++;
                  break;
               }
          }
         if ((def_aa_xref[i] == -1) && (amino_acid_order[i] != '*'))
            {
                error("residue %c in matrices.h is not recognised",
                                       amino_acid_order[i]);
            }
     }

   maxres = 0;
   for (i=0;(c1=nucleic_acid_order[i]);i++)
     {
         for (j=0;(c2=amino_acid_codes[j]);j++)
          {
           if (c1 == c2)
               {
                  def_dna_xref[i] = j;
                  maxres++;
                  break;
               }
          }
         if ((def_dna_xref[i] == -1) && (nucleic_acid_order[i] != '*'))
            {
                error("nucleic acid %c in matrices.h is not recognised",
                                       nucleic_acid_order[i]);
            }
     }
}

sint get_matrix(short *matptr, short *xref, sint matrix[NUMRES][NUMRES], Boolean neg_flag, sint scale)
{
   sint gg_score = 0;
   sint gr_score = 0;
   sint i, j, k, ix = 0;
   sint ti, tj;
   sint  maxres;
   sint av1,av2,av3,min, max;
/*
   default - set all scores to 0
*/
   for (i=0;i<=max_aa;i++)
      for (j=0;j<=max_aa;j++)
          matrix[i][j] = 0;

   ix = 0;
   maxres = 0;
   for (i=0;i<=max_aa;i++)
    {
      ti = xref[i];
      for (j=0;j<=i;j++)
       {
          tj = xref[j]; 
          if ((ti != -1) && (tj != -1))
            {
               k = matptr[ix];
               if (ti==tj)
                  {
                     matrix[ti][ti] = k * scale;
                     maxres++;
                  }
               else
                  {
                     matrix[ti][tj] = k * scale;
                     matrix[tj][ti] = k * scale;
                  }
               ix++;
            }
       }
    }

   --maxres;

   av1 = av2 = av3 = 0;
   for (i=0;i<=max_aa;i++)
    {
      for (j=0;j<=i;j++)
       {
           av1 += matrix[i][j];
           if (i==j)
              {
                 av2 += matrix[i][j];
              }
           else
              {
                 av3 += matrix[i][j];
              }
       }
    }

   av1 /= (maxres*maxres)/2;
   av2 /= maxres;
   av3 /= ((float)(maxres*maxres-maxres))/2;
  mat_avscore = -av3;

  min = max = matrix[0][0];
  for (i=0;i<=max_aa;i++)
    for (j=1;j<=i;j++)
      {
        if (matrix[i][j] < min) min = matrix[i][j];
        if (matrix[i][j] > max) max = matrix[i][j];
      }
if (debug>1) fprintf(stdout,"maxres %d\n",(pint)max_aa);
if (debug>1) fprintf(stdout,"average mismatch score %d\n",(pint)av3);
if (debug>1) fprintf(stdout,"average match score %d\n",(pint)av2);
if (debug>1) fprintf(stdout,"average score %d\n",(pint)av1);

/*
   if requested, make a positive matrix - add -(lowest score) to every entry
*/
  if (neg_flag == FALSE)
   {

if (debug>1) fprintf(stdout,"min %d max %d\n",(pint)min,(pint)max);
      if (min < 0)
        {
           for (i=0;i<=max_aa;i++)
            {
              ti = xref[i];
              if (ti != -1)
                {
                 for (j=0;j<=max_aa;j++)
                   {
                    tj = xref[j];
/*
                    if (tj != -1) matrix[ti][tj] -= (2*av3);
*/
                    if (tj != -1) matrix[ti][tj] -= min;
                   }
                }
            }
        }
/*
       gr_score = av3;
       gg_score = -av3;
*/

   }



  for (i=0;i<gap_pos1;i++)
   {
      matrix[i][gap_pos1] = gr_score;
      matrix[gap_pos1][i] = gr_score;
      matrix[i][gap_pos2] = gr_score;
      matrix[gap_pos2][i] = gr_score;
   }
  matrix[gap_pos1][gap_pos1] = gg_score;
  matrix[gap_pos2][gap_pos2] = gg_score;
  matrix[gap_pos2][gap_pos1] = gg_score;
  matrix[gap_pos1][gap_pos2] = gg_score;

  maxres += 2;

  return(maxres);
}


sint read_matrix_series(char *filename, short *usermat, short *xref)
{
   FILE *fd = NULL, *matfd = NULL;
   char mat_filename[FILENAMELEN];
   char inline1[1024];
   sint  maxres = 0;
   sint nmat;
   sint n,llimit,ulimit;

   if (filename[0] == '\0')
     {
        error("comparison matrix not specified");
        return((sint)0);
     }
   if ((fd=fopen(filename,"r"))==NULL) 
     {
        error("cannot open %s", filename);
        return((sint)0);
     }

/* check the first line to see if it's a series or a single matrix */
   while (fgets(inline1,1024,fd) != NULL)
     {
        if (commentline(inline1)) continue;
	if(linetype(inline1,"CLUSTAL_SERIES"))
		user_series=TRUE;
	else
		user_series=FALSE;
        break;
     }

/* it's a single matrix */
  if(user_series == FALSE)
    {
	fclose(fd);
   	maxres=read_user_matrix(filename,usermat,xref);
   	return(maxres);
    }

/* it's a series of matrices, find the next MATRIX line */
   nmat=0;
   matseries.nmat=0;
   while (fgets(inline1,1024,fd) != NULL)
     {
        if (commentline(inline1)) continue;
	if(linetype(inline1,"MATRIX"))
	{
		if(sscanf(inline1+6,"%d %d %s",&llimit,&ulimit,mat_filename)!=3)
		{
			error("Bad format in file %s\n",filename);
   			fclose(fd);
			return((sint)0);
		}
		if(llimit<0 || llimit > 100 || ulimit <0 || ulimit>100)
		{
			error("Bad format in file %s\n",filename);
   			fclose(fd);
			return((sint)0);
		}
		if(ulimit<=llimit)
		{
			error("in file %s: lower limit is greater than upper (%d-%d)\n",filename,llimit,ulimit);
   			fclose(fd);
			return((sint)0);
		}
   		n=read_user_matrix(mat_filename,&usermatseries[nmat][0][0],&aa_xrefseries[nmat][0]);
		if(n<=0)
		{
			error("Bad format in matrix file %s\n",mat_filename);
   			fclose(fd);
			return((sint)0);
		}
		matseries.mat[nmat].llimit=llimit;
		matseries.mat[nmat].ulimit=ulimit;
		matseries.mat[nmat].matptr=&usermatseries[nmat][0][0];
		matseries.mat[nmat].aa_xref=&aa_xrefseries[nmat][0];
		nmat++;
	}
    }
   fclose(fd);
   matseries.nmat=nmat;

   maxres=n;
   return(maxres);

}

sint read_user_matrix(char *filename, short *usermat, short *xref)
{
   double f;
   FILE *fd;
   sint  numargs,farg;
   sint i, j, k = 0;
   char codes[NUMRES];
   char inline1[1024];
   char *args[NUMRES+4];
   char c1,c2;
   sint ix1, ix = 0;
   sint  maxres = 0;
   float scale;

   if (filename[0] == '\0')
     {
        error("comparison matrix not specified");
       	return((sint)0);
     }

   if ((fd=fopen(filename,"r"))==NULL) 
   {
       	error("cannot open %s", filename);
       	return((sint)0);
   }
   maxres = 0;
   while (fgets(inline1,1024,fd) != NULL)
     {
        if (commentline(inline1)) continue;
	if(linetype(inline1,"CLUSTAL_SERIES"))
   	{
       		error("in %s - single matrix expected.", filename);
		fclose(fd);
       		return((sint)0);
   	}
/*
   read residue characters.
*/
        k = 0;
        for (j=0;j<strlen(inline1);j++)
          {
             if (isalpha((int)inline1[j])) codes[k++] = inline1[j];
             if (k>NUMRES)
                {
                   error("too many entries in matrix %s",filename);
		   fclose(fd);
                   return((sint)0);
                }
          }
        codes[k] = '\0';
        break;
    }

   if (k == 0) 
     {
        error("wrong format in matrix %s",filename);
  	fclose(fd);
        return((sint)0);
     }

/*
   cross-reference the residues
*/
   for (i=0;i<NUMRES;i++) xref[i] = -1;

   maxres = 0;
   for (i=0;(c1=codes[i]);i++)
     {
         for (j=0;(c2=amino_acid_codes[j]);j++)
           if (c1 == c2)
               {
                  xref[i] = j;
                  maxres++;
                  break;
               }
         if ((xref[i] == -1) && (codes[i] != '*'))
            {
                warning("residue %c in matrix %s not recognised",
                                       codes[i],filename);
            }
     }


/*
   get the weights
*/

   ix = ix1 = 0;
   while (fgets(inline1,1024,fd) != NULL)
     {
        if (inline1[0] == '\n') continue;
        if (inline1[0] == '#' ||
            inline1[0] == '!') break;
        numargs = getargs(inline1, args, (int)(k+1));
        if (numargs < maxres)
          {
             error("wrong format in matrix %s",filename);
  	     fclose(fd);
             return((sint)0);
          }
        if (isalpha(args[0][0])) farg=1;
        else farg=0;

/* decide whether the matrix values are float or decimal */
	scale=1.0;
	for(i=0;i<strlen(args[farg]);i++)
		if(args[farg][i]=='.')
		{
/* we've found a float value */
			scale=10.0;
			break;
		}

        for (i=0;i<=ix;i++)
          {
             if (xref[i] != -1)
               {
                  f = atof(args[i+farg]);
                  usermat[ix1++] = (short)(f*scale);
               }
          }
        ix++;
     }
   if (ix != k+1)
     {
        error("wrong format in matrix %s",filename);
  	fclose(fd);
        return((sint)0);
     }


  maxres += 2;
  fclose(fd);

  return(maxres);
}

int getargs(char *inline1,char *args[],int max)
{

	char	*inptr;
/*
#ifndef MAC
	char	*strtok(char *s1, const char *s2);
#endif
*/
	int	i;

	inptr=inline1;
	for (i=0;i<=max;i++)
	{
		if ((args[i]=strtok(inptr," \t\n"))==NULL)
			break;
		inptr=NULL;
	}

	return(i);
}


static Boolean commentline(char *line)
{
        int i;
 
        if(line[0] == '#') return TRUE;
        for(i=0;line[i]!='\n' && line[i]!=EOS;i++) {
                if(!isspace(line[i]))
			return FALSE;
        }
        return TRUE;
}