src/swsse2/matrix.c

   1 /******************************************************************
   2   Copyright 2006 by Michael Farrar.  All rights reserved.
   3   This program may not be sold or incorporated into a commercial product,
   4   in whole or in part, without written consent of Michael Farrar.  For
   5   further information regarding permission for use or reproduction, please
   6   contact: Michael Farrar at farrar.michael@gmail.com.
   7 *******************************************************************/
   8
   9 /*
  10   Written by Michael Farrar, 2006.
  11   Please send bug reports and/or suggestions to farrar.michael@gmail.com.
  12 */
  13
  14 #include <stdio.h>
  15 #include <stdlib.h>
  16 #include <ctype.h>
  17
  18 #include "swsse2.h"
  19 #include "matrix.h"
  20
  21 #define BUF_SIZE 512
  22
  23 char *skipSpaces (char *line)
  24 {
  25     while (isspace (*line)) {
  26         ++line;
  27     }
  28
  29     return line;
  30 }
  31
  32 signed char *readMatrix (char *file)
  33 {
  34     FILE *fp;
  35     char line[BUF_SIZE];
  36
  37     signed char *matrix;
  38
  39     int mark[ALPHA_SIZE];
  40     int order[ALPHA_SIZE];
  41
  42     int done;
  43     int i;
  44
  45     int errors = 0;
  46     int count = 0;
  47
  48     if ((fp = fopen (file, "r")) == NULL) {
  49         fprintf (stderr, "Unable to open file %s\n", file);
  50         exit (-1);
  51     }
  52
  53     matrix = (signed char *) malloc (ALPHA_SIZE * ALPHA_SIZE);
  54     if (!matrix) {
  55         fprintf (stderr, "Unable to allocate memory for scoring matrix\n");
  56         exit (-1);
  57     }
  58
  59     /* initialize the order and mark arrays */
  60     for (i = 0; i < ALPHA_SIZE; ++i) {
  61         order[i] = -1;
  62         mark[i] = -1;
  63     }
  64
  65     /* read the first line of the matrix giving the amino acid order */
  66     done = 0;
  67     while (!done && fgets (line, BUF_SIZE, fp) != NULL) {
  68         char *ptr = skipSpaces (line);
  69         if (*ptr && *ptr != '#') {
  70
  71             while (*ptr && *ptr != '#') {
  72                 int inx = AMINO_ACID_VALUE[*ptr];
  73
  74                 if (inx == -1) {
  75                     fprintf (stderr, "Unknown amino acid %c in %s\n", *ptr, file);
  76                     ++errors;
  77                 } else if (mark[inx] != -1) {
  78                     fprintf (stderr, "Amino acid %c defined twice\n", *ptr);
  79                     ++errors;
  80                 } else if (count >= ALPHA_SIZE) {
  81                     /* this should not happen, but we will be safe */
  82                     fprintf (stderr, "Too many amino acids %d\n", count);
  83                     ++errors;
  84                 } else {
  85                     order[count++] = inx;
  86                     mark[inx] = inx;
  87                 }
  88                 ptr = skipSpaces (ptr + 1);
  89             }
  90
  91             done = 1;
  92         }
  93     }
  94
  95     /* make sure all amino acids are defined */
  96     for (i = 0; i < ALPHA_SIZE; ++i) {
  97         if (order[i] < 0) {
  98             fprintf (stderr, "Missing column for amino acid %c\n",
  99                 AMINO_ACIDS[i]);
 100             ++errors;
 101         }
 102         mark[i] = -1;
 103     }
 104
 105     if (errors > 0) {
 106         fprintf (stderr, "Terminating due to errors in matrix file\n");
 107         exit (-1);
 108     }
 109
 110     /* read the scores for the amino acids */
 111     while (fgets (line, BUF_SIZE, fp) != NULL) {
 112         signed char *row;
 113         char *ptr = skipSpaces (line);
 114         if (*ptr && *ptr != '#') {
 115             char aminoAcid = *ptr;
 116             int inx = AMINO_ACID_VALUE[*ptr];
 117             if (inx == -1) {
 118                 fprintf (stderr, "Unknown amino acid %c in matrix\n", *ptr);
 119                 ++errors;
 120             } else if (mark[inx] != -1) {
 121                 fprintf (stderr, "Row %c defined twice\n", *ptr);
 122                 ++errors;
 123             }
 124
 125             row = &matrix[inx * ALPHA_SIZE];
 126
 127             for (i = 0; i < ALPHA_SIZE; ++i) {
 128                 int sign = 1;
 129                 int num = 0;
 130
 131                 ptr = skipSpaces (ptr + 1);
 132
 133                 /* check the sign */
 134                 if (*ptr == '-') {
 135                     sign = -1;
 136                     ++ptr;
 137                 }
 138
 139                 do {
 140                     if (*ptr >= '0' && *ptr <= '9') {
 141                         num = num * 10 + (*ptr - '0');
 142                         ptr++;
 143                     } else {
 144                         char name[16];
 145                         char *pName;
 146                         if (isspace (*ptr)) {
 147                             pName = "space";
 148                         } else if (*ptr == 0) {
 149                             pName = "end of line";
 150                         } else {
 151                             name[0] = *ptr;
 152                             name[1] = 0;
 153                             pName = name;
 154                         }
 155                         fprintf (stderr, "Row %c Expecting digit found %s\n",
 156                             aminoAcid, pName);
 157                         exit (-1);
 158                     }
 159                 } while (*ptr && !isspace (*ptr));
 160
 161                 num = num * sign;
 162
 163                 if (num < -128 || num > 127) {
 164                     fprintf (stderr, "Weight %d out of range row %c\n",
 165                         aminoAcid, num);
 166                     ++errors;
 167                     num = 0;
 168                 }
 169
 170                 row[order[i]] = (char) num;
 171             }
 172
 173             if (i < ALPHA_SIZE) {
 174                 fprintf (stderr, "Amino acid row %c incomplete\n", aminoAcid);
 175                 ++errors;
 176             }
 177
 178             mark[inx] = 1;
 179         }
 180     }
 181
 182     /* make sure all amino acids are defined */
 183     for (i = 0; i < ALPHA_SIZE; ++i) {
 184         if (mark[i] < 0) {
 185             fprintf (stderr, "Missing row for amino acid %c\n",
 186                 AMINO_ACIDS[i]);
 187             ++errors;
 188         }
 189     }
 190
 191     if (errors) {
 192         fprintf (stderr, "Terminating due to errors in matrix %s\n", file);
 193     }
 194
 195     fclose (fp);
 196
 197     return matrix;
 198 }