sources/pairwise/pairwise.h

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <errno.h>
   4 #include <string.h>
   5 #include <stdarg.h>
   6
   7 #include <sys/types.h>
   8 #include <sys/stat.h>
   9
  10 #define ARRAY 50
  11 #define FILEBUF 1024
  12
  13 /* Macro's added by Patrick, mucho neato, aparently ifndef aren't POSIX
  14  * anymore */
  15 #if !defined(MAX)
  16 #define MAX(A, B)       ((A) > (B) ? (A) : (B))
  17 #endif
  18
  19 #if !defined(MIN)
  20 #define MIN(A, B)       ((A) < (B) ? (A) : (B))
  21 #endif
  22
  23 /* Macro's added by Jon, allows checking of NULL pointers from memory
  24  * allocations without function calls, but how can you return values? */
  25 #if !defined(MALLOC)
  26 #define MALLOC(PTR, SIZE) do { PTR = malloc(SIZE); if (PTR == NULL) fatal_sys_error("malloc returned NULL"); } while (0)
  27 #endif
  28
  29 #if !defined(REALLOC)
  30 #define REALLOC(PTR, SIZE) do { PTR = realloc(PTR, SIZE); if (PTR == NULL && SIZE != 0) fatal_sys_error("realloc returned NULL"); } while (0)
  31 #endif
  32
  33 /* Structure for sequences from FASTA files */
  34 struct fasta {
  35         char *id;
  36         char *seq;
  37         int start;
  38         int end;
  39         int numres;
  40 };
  41
  42 void fatal_error (char *message);
  43 FILE * xfopen (const char *path, const char *mode);
  44 void * xmalloc (size_t size);
  45 void * xrealloc (void *ptr, size_t size);
  46 int diff (char a, char b);
  47 float pairwise (struct fasta *a, struct fasta *b);
  48 void populate (struct fasta *a);
  49 void check_length (struct fasta **array);
  50 struct fasta ** read_fasta (FILE *fh);
  51 char ** do_pairwise (struct fasta **array);
  52 void free_fasta (struct fasta *);
  53 int jsprintf (char *str, const char *const fmt, ...);
  54
  55 void
  56 free_fasta (struct fasta *seq) {
  57         free(seq->id);
  58         free(seq->seq);
  59         free(seq);
  60 }
  61
  62 void
  63 fatal_error (char *message) {
  64         fprintf(stderr, "%s", message);
  65         exit(EXIT_FAILURE);
  66 }
  67
  68 void
  69 fatal_sys_error (char *message) {
  70         perror(message);
  71         exit(EXIT_FAILURE);
  72 }
  73
  74 FILE *
  75 xfopen (const char *path, const char *mode) {
  76         FILE *fh;
  77         fh = fopen(path, mode);
  78         if (fh == NULL)
  79                 fatal_sys_error("fopen returned NULL");
  80         return fh;
  81 }
  82
  83 void *
  84 xmalloc (size_t size) {
  85         void *ptr;
  86         ptr = (void *) malloc(size);
  87         if (ptr == NULL)
  88                 fatal_sys_error("malloc returned NULL");
  89         return ptr;
  90 }
  91
  92 void *
  93 xrealloc (void *ptr, size_t size) {
  94         ptr = (void *) realloc(ptr, size);
  95         if (ptr == NULL && size != 0)
  96                 fatal_sys_error("realloc returned NULL");
  97         return ptr;
  98 }
  99
 100 /* If the same resdidue return 1 else return 0 */
 101 int
 102 diff (char a, char b) {
 103         if (a == b && a == '-')
 104                 return 0;
 105         if (a == b)
 106                 return 1;
 107         return 0;
 108 }
 109
 110 /* Does the pairwise comparison
 111  * Compares the regions of the alignment that overlap, common gaps are
 112  * ignored, and the number of common residues divided by the length of
 113  * the shortest sequence is the identity.
 114  * */
 115 float
 116 pairwise (struct fasta *a, struct fasta *b) {
 117         float result;
 118         int start, end, numres, i, id = 0;
 119
 120         /* If the sequences don't overlap then the seq ID is 0 */
 121         if (a->end < b->start || a->start > b->end)
 122                 return 0;
 123
 124         start = MIN( a->start, b->start );
 125         end = MAX( a->end, b->end );
 126         numres = MAX( a->numres, b->numres );
 127
 128         for (i = start; i < end; i++)
 129                 id += diff(a->seq[i], b->seq[i]);
 130
 131         result = 100 * (float) id / (float) numres;
 132
 133         return result;
 134 }
 135
 136 /* "Populate" the rest of the fasta structure with information
 137  * Gets the start position of the sequence (first non-gap character)
 138  * Gets the end position (the last non-gap character)
 139  * Gets the number of residues in the sequence (that aren't gaps)
 140  * Gaps are represented by '-'
 141  * */
 142 void
 143 populate (struct fasta *a) {
 144         int i;
 145         int len = strlen(a->seq);
 146
 147         a->numres = 0;
 148
 149         for (i = 0; i <= len; i++) {
 150                 if (a->seq[i] != '-') {
 151                         a->start = i;
 152                         break;
 153                 }
 154         }
 155         for (i = len; i > 0; i--) {
 156                 if (a->seq[i] != '-') {
 157                         a->end = i;
 158                         break;
 159                 }
 160         }
 161         for (i = a->start; i < a->end; i++) {
 162                 if (a->seq[i] != '-') {
 163                         a->numres++;
 164                 }
 165         }
 166 }
 167
 168 /* Makes sure that all of the sequences are the same length */
 169 void
 170 check_length (struct fasta **array) {
 171         int i, length;
 172
 173         if (array[0] != NULL)
 174                 length = strlen(array[0]->seq);
 175         else {
 176                 fatal_error("check_length() not passed an array of fasta structs\n");
 177         }
 178
 179         for (i = 0; array[i] != NULL; i++) {
 180                 if (length != (int) strlen(array[i]->seq)) {
 181                         fatal_error("Not all of the sequences are the same length\n");
 182                 }
 183         }
 184 }
 185
 186 /* Reads in a FASTA file and returns an array of fasta structures */
 187 struct fasta **
 188 read_fasta (FILE *fh) {
 189         int i, j, k, c, filesize = 0;
 190         char *file;
 191         struct fasta **array;
 192         struct stat file_stat;
 193
 194         array = (struct fasta **) malloc(ARRAY * sizeof(struct fasta *));
 195
 196         /* Find the size of the file and allocate it */
 197         fstat(fileno(fh), &file_stat);
 198         if (! S_ISREG(file_stat.st_mode)) {
 199                 fatal_error("Input file is not a regular file");
 200         }
 201         file = (char *) malloc(sizeof(char) * (file_stat.st_size + 1));
 202
 203         /* Read in the file and null terminate it */
 204         while ((c = getc(fh)) != EOF) {
 205                 file[filesize++] = c;
 206         }
 207         file[filesize] = '\0';
 208
 209         /* Parse the FASTA file into an array of structures */
 210         for (i = 0, j = 0, k = 0; i < filesize; i++) {
 211
 212                 /* Begining of FASTA record */
 213                 if (file[i] == '>') {
 214                         if (j % ARRAY == 0)
 215                                 array = (struct fasta **) xrealloc(array, sizeof(struct fasta *) * (j + ARRAY));
 216
 217                         /* Initalise the structure */
 218                         array[j] = (struct fasta *) xmalloc(sizeof(struct fasta));
 219                         array[j]->id = (char *) xmalloc(sizeof(char));
 220                         array[j]->seq = (char *) xmalloc(sizeof(char));
 221                         array[j]->numres = 0;
 222
 223                         /* Read in the ID, terminated by a newline or NULL (end of string) */
 224                         i++;
 225                         while (file[i] != '\0' && file[i] != '\n') {
 226                                 if (k % FILEBUF == 0)
 227                                         array[j]->id = (char *) xrealloc(array[j]->id, sizeof(char) * (FILEBUF + k));
 228                                 array[j]->id[k++] = file[i++];
 229                         }
 230                         array[j]->id = (char *) xrealloc(array[j]->id, sizeof(char) * k + 1);
 231                         array[j]->id[k] = '\0';
 232                         k = 0;
 233
 234                         /* read in the sequence, terminated by a '>' */
 235                         for (k = 0; file[i] != '\0' && file[i] != '>'; i++) {
 236                                 if (file[i] != '\n') {
 237                                         if (k % FILEBUF == 0)
 238                                                 array[j]->seq = (char *) xrealloc(array[j]->seq, sizeof(char) * (FILEBUF + k));
 239                                         array[j]->seq[k++] = file[i];
 240                                 }
 241                         }
 242                         array[j]->seq = (char *) xrealloc(array[j]->seq, sizeof(char) * k + 1);
 243                         array[j]->seq[k] = '\0';
 244
 245                         k = 0;
 246                         i--;
 247                         j++;
 248                 }
 249         }
 250
 251         free(file);
 252         array = (struct fasta **) xrealloc(array, sizeof(struct fasta *) * (j + 1));
 253         array[j] = NULL;
 254
 255         /* find the start and end points for the alignments */
 256         for (i = 0; array[i] != NULL; i++)
 257                 populate(array[i]);
 258         return array;
 259 }
 260
 261 char **
 262 do_pairwise (struct fasta **array) {
 263         int j = 0, i = 0, offset = 0, size;
 264         char **output;
 265
 266         /* Find the size of the output and preallocate the space in the output
 267          * array */
 268         {
 269                 while (array[i] != NULL)
 270                         i++;
 271
 272                 /* Number of IDs, each ID, then the upper corner distances, then for
 273                  * NULL terminator */
 274                 size = 1 + i + (i * i - i) / 2 + 1;
 275
 276                 output = (char **) xmalloc(sizeof(char *) * size);
 277
 278                 for (i = 0; i < size - 1; i++) {
 279                         output[i] = (char *) xmalloc(sizeof(char) * FILEBUF);
 280                         memset(output[i], 0x00, sizeof(char) * FILEBUF);
 281                 }
 282                 output[size - 1] = NULL;
 283         }
 284
 285         /* The number of IDs */
 286         i = 0;
 287         while (array[i] != NULL)
 288                 i++;
 289         jsprintf(output[offset++], "%i", i);
 290
 291         /* If this statement is above the other one then it core dumps!? */
 292
 293         /* Each ID */
 294         for (i = 0; array[i] != NULL; i++)
 295                 jsprintf(output[offset++], "%s", array[i]->id);
 296
 297         /* Upper corner pairwise comparisons */
 298         for (i = 0; array[i] != NULL; i++) {
 299                 for (j = i + 1; array[j] != NULL; j++) {
 300                         jsprintf(output[offset++], "%f", pairwise(array[i], array[j]));
 301                 }
 302         }
 303
 304         return output;
 305 }
 306
 307 /* Grows buffer dynamically to make input fit inside it.
 308  * Requires malloc'd pointer to string.
 309  * Returns length of string.
 310  * */
 311 int
 312 jsprintf (char *str, const char *const fmt, ...) {
 313         int i, size;
 314         va_list ap;
 315         va_start(ap, fmt);
 316
 317         str = (char *) realloc(str, sizeof(char) * FILEBUF);
 318
 319         for (i = 1; (size = vsnprintf(str, sizeof(char) * FILEBUF * i, fmt, ap)) == -1;  i++) {
 320                 if (i == 100)
 321                         fatal_error("Reached 100 iterations of vsnprinft, it's probably blown up. Check input");
 322                 str = (char *) xrealloc(str, sizeof(char) * FILEBUF * i);
 323         }
 324
 325         str = (char *) xrealloc(str, sizeof(char) * size + 1);
 326
 327         va_end(ap);
 328         return size;
 329 }