+++ /dev/null
-/*****************************************************************
- * SQUID - a library of functions for biological sequence analysis
- * Copyright (C) 1992-2002 Washington University School of Medicine
- *
- * This source code is freely distributed under the terms of the
- * GNU General Public License. See the files COPYRIGHT and LICENSE
- * for details.
- *****************************************************************/
-
-/* dayhoff.c
- *
- * Routines for dealing with PAM matrices.
- *
- * Includes:
- * ParsePAMFile() -- read a PAM matrix from disk.
- *
- *
- * SRE - Fri Apr 2 11:23:45 1993
- * RCS $Id: dayhoff.c 217 2011-03-19 10:27:10Z andreas $ (Original squid RCS Id: dayhoff.c,v 1.5 2002/07/03 15:03:39 eddy Exp)
- */
-
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-#include <math.h>
-#include <ctype.h>
-#include "squid.h"
-
-/* Function: ParsePAMFile()
- *
- * Purpose: Given a pointer to an open file containing a PAM matrix,
- * parse the file and allocate and fill a 2D array of
- * floats containing the matrix. The PAM file is
- * assumed to be in the format that NCBI distributes
- * with BLAST. BLOSUM matrices also work fine, as
- * produced by Henikoff's program "MATBLAS".
- *
- * Parses both old format and new format BLAST matrices.
- * Old format just had rows of integers.
- * New format includes a leading character on each row.
- *
- * The PAM matrix is a 27x27 matrix, 0=A..25=Z,26=*.
- * Note that it's not a 20x20 matrix as you might expect;
- * this is for speed of indexing as well as the ability
- * to deal with ambiguous characters.
- *
- * Args: fp - open PAM file
- * ret_pam - RETURN: pam matrix, integers
- * ret_scale - RETURN: scale factor for converting
- * to real Sij. For instance, PAM120 is
- * given in units of ln(2)/2. This may
- * be passed as NULL if the caller
- * doesn't care.
- *
- * Returns: 1 on success; 0 on failure and sets squid_errno to
- * indicate the cause. ret_pam is allocated here and
- * must be freed by the caller (use FreePAM).
- */
-int
-ParsePAMFile(FILE *fp, int ***ret_pam, float *ret_scale)
-{
- int **pam;
- char buffer[512]; /* input buffer from fp */
- int order[27]; /* order of fields, obtained from header */
- int nsymbols; /* total number of symbols in matrix */
- char *sptr;
- int idx;
- int row, col;
- float scale;
- int gotscale = FALSE;
-
- scale = 0.0; /* just to silence gcc uninit warnings */
- if (fp == NULL) { squid_errno = SQERR_NODATA; return 0; }
-
- /* Look at the first non-blank, non-comment line in the file.
- * It gives single-letter codes in the order the PAM matrix
- * is arrayed in the file.
- */
- do {
- if (fgets(buffer, 512, fp) == NULL)
- { squid_errno = SQERR_NODATA; return 0; }
-
- /* Get the scale factor from the header.
- * For BLOSUM files, we assume the line looks like:
- * BLOSUM Clustered Scoring Matrix in 1/2 Bit Units
- * and we assume that the fraction is always 1/x;
- *
- * For PAM files, we assume the line looks like:
- * PAM 120 substitution matrix, scale = ln(2)/2 = 0.346574
- * and we assume that the number following the final '=' is our scale
- */
- if (strstr(buffer, "BLOSUM Clustered Scoring Matrix") != NULL &&
- (sptr = strchr(buffer, '/')) != NULL)
- {
- sptr++;
- if (! isdigit((int) (*sptr))) { squid_errno = SQERR_FORMAT; return 0; }
- scale = (float) (log(2.0) / atof(sptr));
- gotscale = TRUE;
- }
- else if (strstr(buffer, "substitution matrix,") != NULL)
- {
- while ((sptr = strrchr(buffer, '=')) != NULL) {
- sptr += 2;
- if (IsReal(sptr)) {
- scale = atof(sptr);
- gotscale = TRUE;
- break;
- }
- }
- }
- } while ((sptr = strtok(buffer, " \t\n")) == NULL || *sptr == '#');
-
- idx = 0;
- do {
- order[idx] = (int) *sptr - (int) 'A';
- if (order[idx] < 0 || order[idx] > 25) order[idx] = 26;
- idx++;
- } while ((sptr = strtok(NULL, " \t\n")) != NULL);
- nsymbols = idx;
-
- /* Allocate a pam matrix. For speed of indexing, we use
- * a 27x27 matrix so we can do lookups using the ASCII codes
- * of amino acid single-letter representations, plus one
- * extra field to deal with the "*" (terminators).
- */
- if ((pam = (int **) calloc (27, sizeof(int *))) == NULL)
- Die("calloc failed");
- for (idx = 0; idx < 27; idx++)
- if ((pam[idx] = (int *) calloc (27, sizeof(int))) == NULL)
- Die("calloc failed");
-
- /* Parse the rest of the file.
- */
- for (row = 0; row < nsymbols; row++)
- {
- if (fgets(buffer, 512, fp) == NULL)
- { squid_errno = SQERR_NODATA; return 0; }
-
- if ((sptr = strtok(buffer, " \t\n")) == NULL)
- { squid_errno = SQERR_NODATA; return 0; }
- for (col = 0; col < nsymbols; col++)
- {
- if (sptr == NULL) { squid_errno = SQERR_NODATA; return 0; }
-
- /* Watch out for new BLAST format, with leading characters
- */
- if (*sptr == '*' || isalpha((int) *sptr))
- col--; /* hack hack */
- else
- pam [order[row]] [order[col]] = atoi(sptr);
-
- sptr = strtok(NULL, " \t\n");
- }
- }
-
- /* Return
- */
- if (ret_scale != NULL)
- {
- if (gotscale) *ret_scale = scale;
- else
- {
-#ifdef CLUSTALO
- Warning("Failed to parse PAM matrix scale factor. Defaulting to ln(2)/2!");
-#else
- Warn("Failed to parse PAM matrix scale factor. Defaulting to ln(2)/2!");
-#endif
- *ret_scale = log(2.0) / 2.0;
- }
- }
- *ret_pam = pam;
- return 1;
-}
git://source.jalview.org / jabaws.git / blobdiff