pairwise.c

Go to the documentation of this file.

/* Pairwise identity program.
 *
 * $Id
 *
 * Reads in a gapped alignment sequences in a FASTA file format.
 * Assume gaps are represented by '-' characters.
 *
 * For each pair calculates the sequence percentage identity, see comments
 * for pairwise() for details.
 *
 * Output is in the format for OC to read in, i.e. the number of sequences,
 * followed by the sequence ID's, followed by the pairwise comparisons.
 *
 * Thu Dec  5 14:38:43 GMT 2002 Jon
 * Added checking for eqilength sequences
 *
 * Thu Jul 24 12:14:04 BST 2003 Jon
 * Moved FASTA reading into read_fasta function
 * Checked into CVS
 *
 * */

#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <string.h>
#include "pairwise.h"

#define ARRAY 50
#define FILEBUF 1000

/* Macro's added by Patrick, mucho neato, aparently ifndef aren't POSIX
 * anymore */
#if !defined(MAX)
#define MAX(A, B)       ((A) > (B) ? (A) : (B))
#endif

#if !defined(MIN)
#define MIN(A, B)       ((A) < (B) ? (A) : (B))
#endif

/* Macro's added by Jon, allows checking of NULL pointers from memory
 * allocations without function calls, but how can you return values? */
#if !defined(MALLOC)
#define MALLOC(PTR, SIZE) do { PTR = malloc(SIZE); if (PTR == NULL) fatal_sys_error("malloc returned NULL"); } while (0)
#endif

#if !defined(REALLOC)
#define REALLOC(PTR, SIZE) do { PTR = realloc(PTR, SIZE); if (PTR == NULL && SIZE != 0) fatal_sys_error("realloc returned NULL"); } while (0)
#endif

void
fatal_error (char *message) {
        printf(message);
        exit(EXIT_FAILURE);
}

void
fatal_sys_error (char *message) {
        perror(message);
        exit(EXIT_FAILURE);
}

FILE *
xfopen (const char *path, const char *mode) {
        FILE *fh;
        fh = fopen(path, mode);
        if (fh == NULL)
                fatal_sys_error("fopen returned NULL");
        return fh;
}

void *
xmalloc (size_t size) {
        void *ptr;
        ptr = (void *) malloc(size);
        if (ptr == NULL)
                fatal_sys_error("malloc returned NULL");
        return ptr;
}

void *
xrealloc (void *ptr, size_t size) {
        ptr = (void *) realloc(ptr, size);
        if (ptr == NULL && size != 0)
                fatal_sys_error("realloc returned NULL");
        return ptr;
}

/* If the same resdidue return 1 else return 0 */
int 
diff (char a, char b) {
        if (a == b && a == '-')
                return 0;
        if (a == b)
                return 1;
        return 0;
}

/* Does the pairwise comparison
 * Compares the regions of the alignment that overlap, common gaps are
 * ignored, and the number of common residues divided by the length of
 * the shortest sequence is the identity.
 * */
float
pairwise (struct fasta *a, struct fasta *b) {
        float result;
        int start, end, numres, i, id = 0;
        
        /* If the sequences don't overlap then the seq ID is 0 */
        if (a->end < b->start || a->start > b->end)
                return 0;

        start = MIN( a->start, b->start );
        end = MAX( a->end, b->end );
        numres = MAX( a->numres, b->numres );

        for (i = start; i < end; i++)
                id += diff(a->seq[i], b->seq[i]);

        result = 100 * (float) id / (float) numres;

        return result;
}

/* "Populate" the rest of the fasta structure with information
 * Gets the start position of the sequence (first non-gap character)
 * Gets the end position (the last non-gap character)
 * Gets the number of residues in the sequence (that aren't gaps)
 * Gaps are represented by '-'
 * */
void
populate (struct fasta *a) {
        int i;
        int len = strlen(a->seq);

        a->numres = 0;
        
        for (i = 0; i <= len; i++) {
                if (a->seq[i] != '-') {
                        a->start = i;
                        break;
                }
        }
        for (i = len; i > 0; i--) {
                if (a->seq[i] != '-') {
                        a->end = i;
                        break;
                }
        }
        for (i = a->start; i < a->end; i++) {
                if (a->seq[i] != '-') {
                        a->numres++;
                }
        }
}

/* Makes sure that all of the sequences are the same length */
void
check_length (struct fasta **array) {
        int i, length;

        if (array[0] != NULL)
                length = strlen(array[0]->seq);
        else {
                fprintf(stderr, "check_length() not passed an array of fasta structs\n");
                return;
        }

        for (i = 0; array[i] != NULL; i++) {
                if (length != strlen(array[i]->seq)) {
                        fatal_error("Not all of the sequences are the same length\n");
                }
        }
}

/* Reads in a FASTA file and returns an array of fasta structures */
struct fasta **
read_fasta (FILE *fh) {
        int i, j, k, c, filesize = 1;
        char *file;
        struct fasta **array;

        array = (struct fasta **) xmalloc(ARRAY * sizeof(struct fasta *));
        file = (char *) xmalloc(sizeof(char));

        /* Allocate initial space for the file */
        file = xrealloc(file, sizeof(char) * filesize);
        file[filesize] = '\0';

        /* Read in the file */
        while ((c = getc(fh)) != EOF) {
                if (filesize % FILEBUF)
                        file = xrealloc(file, sizeof(char) * (FILEBUF + filesize));

                file[filesize] = c;
                filesize++;
        }

        /* Parse the FASTA file into an array of structures */
        for (i = 0, j = 0, k = 0; i < filesize; i++) {
                if (file[i] == '>') {
                        if (j % ARRAY == 0)
                                array = (struct fasta **) xrealloc(array, sizeof(struct fasta *) * (j + ARRAY));

                        array[j] = (struct fasta *) xmalloc(sizeof(struct fasta));
                        array[j]->id = (char *) xmalloc(sizeof(char));
                        array[j]->seq = (char *) xmalloc(sizeof(char));
                        array[j]->numres = 0;

                        i++;
                        while (file[i] != '\0' && file[i] != '\n') {
                                if (k % ARRAY == 0)
                                        array[j]->id = (char *) xrealloc(array[j]->id, sizeof(char) * (ARRAY + k));
                                array[j]->id[k] = file[i];
                                k++; i++;
                        }
                        array[j]->id = (char *) xrealloc(array[j]->id, sizeof(char) * (ARRAY + k));
                        array[j]->id[k] = '\0';
                        k = 0;

                        while (file[i] != '\0' && file[i] != '>') {
                                if (file[i] == '\n') {
                                        i++;
                                        continue;
                                }
                                if (k % ARRAY == 0)
                                        array[j]->seq = (char *) xrealloc(array[j]->seq, sizeof(char) * (ARRAY + k));
                                array[j]->seq[k] = file[i];
                                k++; i++;
                        }
                        array[j]->seq = (char *) xrealloc(array[j]->seq, sizeof(char) * (ARRAY + k));
                        array[j]->seq[k] = '\0';
                        k = 0;
                        i--;
                        j++;
                }
        }

        free(file);
        array = (struct fasta **) xrealloc(array, sizeof(struct fasta *) * j);
        array[j] = NULL;

        /* find the start and end points for the alignments */
        for (i = 0; array[i] != NULL; i++) {
                populate(array[i]);
        }

        return array;
}

int
main (int argc, char **argv) {
        FILE *fh;
        struct fasta **array;
        int i = 0, j;

        /* Read in the FASTA file */
        if (argc == 2) 
                fh = xfopen(argv[1], "r");
        else
                fh = stdin;

        array = read_fasta(fh);
        fclose(fh);

        check_length(array);
        
        /* start the OC output */
        while (array[i] != NULL)
                i++;

        fprintf(stdout, "%i\n", i);

        for (i = 0; array[i] != NULL; i++)
                fprintf(stdout, "%s\n", array[i]->id);
        
        /* do the pairwise comparison */
        for (i = 0; array[i] != NULL; i++) {
                for (j = i + 1; array[j] != NULL; j++) {
                        printf("%f\n", pairwise(array[i], array[j]));
                }
        }

        return EXIT_SUCCESS;
}

---