source: tags/initial/GDE/CLUSTALW/readmat.c

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


/*
 *   Prototypes
 */


/*
 *   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;

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 = 1;
   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++;
            }
       }
    }

   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];
              }
       }
    }

   --maxres;
   av1 /= (maxres*maxres)/2;
   av2 /= maxres;
   av3 /= ((float)(maxres*maxres-maxres))/2;
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);

  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 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;
*/

   }

  mat_avscore = -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_pos2][gap_pos1] = gg_score;

  maxres += 2;

  return(maxres);
}


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

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

        maxres = 0;
        while (fgets(inline1,1024,fd) != NULL)
          {
             if ((inline1[0] == '\0') || (inline1[0] == '#')) continue;
/*
   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 %s",filename);
                        return((sint)0);
                     }
               }
             codes[k] = '\0';
             break;
          }

        if (k == 0) 
          {
             error("wrong format in %s",filename);
             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 %s not recognised",
                                       codes[i],filename);
            }
     }


/*
   get the weights
*/

        ix = ix1 = 0;
        while (fgets(inline1,1024,fd) != NULL)
          {
             if (inline1[0] == '\n') continue;
             numargs = getargs(inline1, args, (int)(k+1));
             if (numargs == 0)
               {
                  error("wrong format in %s", filename);
                  return((sint)0);
               }
             if (isalpha(args[0][0])) farg=1;
             else farg=0;
             for (i=0;i<=ix;i++)
               {
                  if (xref[i] != -1)
                    {
                       f = atof(args[i+farg]);
                       usermat[ix1++] = (short)(f*10.0);
                    }
               }
             ix++;
          }
        if (ix != k+1)
          {
             error("wrong format in %s", filename);
             return((sint)0);
          }
     }

    fclose(fd);

  maxres += 2;

  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);
}