src/swsse2/fastalib.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 "fastalib.h"
  20
  21 #define READ_BUFFER_SIZE (128 * 1024)
  22 #define SEQ_NAME_SIZE    (128)
  23
  24 FASTA_LIB *openLib (char *file, int pad)
  25 {
  26     FILE *fp;
  27
  28     FASTA_LIB *lib;
  29
  30     if ((fp = fopen (file, "r")) == NULL) {
  31         fprintf (stderr, "Unable to open file %s\n", file);
  32         exit (-1);
  33     }
  34
  35     lib = (FASTA_LIB *) malloc (sizeof (FASTA_LIB));
  36     if (!lib) {
  37         fprintf (stderr, "Unable to allocate memory for library header\n");
  38         exit (-1);
  39     }
  40
  41     lib->readBuffer = (char *) malloc (READ_BUFFER_SIZE);
  42     if (!lib->readBuffer) {
  43         fprintf (stderr, "Unable to allocate memory for read buffer\n");
  44         exit (-1);
  45     }
  46
  47     lib->seqBuffer = (unsigned char *) malloc (MAX_SEQ_LENGTH);
  48     if (!lib->seqBuffer) {
  49         fprintf (stderr, "Unable to allocate memory for sequence\n");
  50         exit (-1);
  51     }
  52
  53     lib->seqName = (char *) malloc (SEQ_NAME_SIZE);
  54     if (!lib->seqName) {
  55         fprintf (stderr, "Unable to allocate memory for sequence name\n");
  56         exit (-1);
  57     }
  58
  59     lib->size = (int) fread (lib->readBuffer, sizeof (char), READ_BUFFER_SIZE, fp);
  60     if (lib->size == 0 && !feof (fp)) {
  61         fprintf (stderr, "Error (%d) reading fasta file\n", ferror (fp));
  62         exit (-1);
  63     }
  64
  65     lib->pos = 0;
  66
  67     lib->fp = fp;
  68
  69     lib->sequences = 0;
  70     lib->residues = 0;
  71
  72     lib->pad;
  73
  74     return lib;
  75 }
  76
  77 static int
  78 readNextBlock (FASTA_LIB *lib)
  79 {
  80     FILE *fp = lib->fp;
  81     size_t size;
  82
  83     size = fread (lib->readBuffer, sizeof (char), READ_BUFFER_SIZE, fp);
  84     if (lib->size == 0 && !feof (fp)) {
  85         fprintf (stderr, "Error (%d) reading fasta file\n", ferror (fp));
  86         exit (-1);
  87     }
  88
  89     lib->pos = 0;
  90     lib->size = (int) size;
  91
  92     return lib->size;
  93 }
  94
  95 unsigned char *
  96 nextSeq (FASTA_LIB *lib, int *length)
  97 {
  98     int inx;
  99     int size;
 100     int done;
 101     int len;
 102
 103     char *name = lib->seqName;
 104     unsigned char *seq = lib->seqBuffer;
 105
 106     /* check if we are at the end of the library */
 107     if (lib->size == 0) {
 108         *length = 0;
 109         return NULL;
 110     }
 111
 112     if (lib->pos == lib->size) {
 113         readNextBlock (lib);
 114     }
 115
 116     inx = lib->pos;
 117
 118     /* check for the start of a sequence */
 119     if (lib->readBuffer[inx] != '>') {
 120         fprintf (stderr, "Error parsing fasta file expecting > found %c\n",
 121             lib->readBuffer[inx]);
 122         exit (-1);
 123     }
 124
 125     ++inx;
 126
 127     /* read in the sequence name */
 128     len = 0;
 129     done = 0;
 130     do {
 131         if (inx >= lib->size) {
 132             size = readNextBlock (lib);
 133             if (size == 0) {
 134                 *length = 0;
 135                 return NULL;
 136             }
 137             inx = lib->pos;
 138         } else if (lib->readBuffer[inx] == '\n') {
 139             *name = '\0';
 140             done = 1;
 141         } else if (len < SEQ_NAME_SIZE - 1) {
 142             *name++ = lib->readBuffer[inx];
 143             len++;
 144         }
 145         ++inx;
 146     } while (!done);
 147
 148     lib->pos = inx;
 149
 150     /* read in the sequence */
 151     len = 0;
 152     done = 0;
 153     do {
 154         if (inx >= lib->size) {
 155             size = readNextBlock (lib);
 156             if (size == 0) {
 157                 *seq = '\0';
 158                 done = 1;
 159             }
 160             inx = 0;
 161         } else if (isspace(lib->readBuffer[inx])) {
 162             ++inx;
 163         } else if (lib->readBuffer[inx] == '>') {
 164             *seq = '\0';
 165             done = 1;
 166         } else if (len >= MAX_SEQ_LENGTH) {
 167             fprintf (stderr, "Sequence %s exceeds maximum length\n",
 168                 lib->seqName);
 169             exit (-1);
 170         } else {
 171             int value = AMINO_ACID_VALUE[lib->readBuffer[inx]];
 172             if (value == -1) {
 173                 fprintf (stderr, "Unknown amino acid %c in sequence %s\n",
 174                     lib->readBuffer[inx], lib->seqName);
 175                 exit (-1);
 176             }
 177             *seq++ = (char) value;
 178             inx++;
 179             len++;
 180         }
 181     } while (!done);
 182
 183     lib->pos = inx;
 184     *length = len;
 185
 186     lib->sequences++;
 187     lib->residues += len;
 188
 189     /*  check if we need to pad the sequence to a multiple of 16  */
 190     if (lib->pad) {
 191         inx = 16 - (len % 16);
 192         while (inx--) {
 193             *seq++ = ALPHA_SIZE;
 194         }
 195         *seq = '\0';
 196     }
 197
 198     return lib->seqBuffer;
 199 }
 200
 201 void closeLib (FASTA_LIB *lib)
 202 {
 203     fclose (lib->fp);
 204
 205     free (lib->readBuffer);
 206     free (lib->seqBuffer);
 207     free (lib->seqName);
 208
 209     free (lib);
 210 }