A simple implementation of the Needleman-Wunsch Algorithm in C.

needle.c

#include <stdio.h>
#include <stdlib.h>

void reverse(char** Pstring);
int match(char a, char b);

int main(int argc, char* argv[]) {
  if(argc < 2) {
    fprintf(stderr, "Usage: needle <sequence 1> <sequence 2>");
  }

  FILE* fp1   = NULL;
  FILE* fp2   = NULL;
  char* seq1  = NULL;
  char* seq2  = NULL;
  char* aln1  = NULL;
  char* aln2  = NULL;
  int   len1  = 0;
  int   len2  = 0;
  int   len3  = 0;
  int   len4  = 0;
  int** score = NULL;
  int** point = NULL;
  int   c, i, j;

  int gap = -1; // Gap penalty

  seq1 = (char*)malloc(sizeof(char));
  seq2 = (char*)malloc(sizeof(char));

  aln1 = (char*)malloc(sizeof(char));
  aln2 = (char*)malloc(sizeof(char));

  if((fp1 = fopen(argv[1], "r")) != NULL && (fp2 = fopen(argv[2], "r")) != NULL) {
    while((c = fgetc(fp1)) != EOF) {
      if(c < 0x41 || (0x5a < c && c < 0x61) || 0x7a < c) continue;
      seq1 = (char*)realloc(seq1, sizeof(char) * (len1 + 1));
      seq1[len1] = c;
      ++len1;
    }

    while((c = fgetc(fp2)) != EOF) {
      if(c < 0x41 || (0x5a < c && c < 0x61) || 0x7a < c) continue;
      seq2 = (char*)realloc(seq2, sizeof(char) * (len2 + 1));
      seq2[len2] = c;
      ++len2;
    }
    
    fclose(fp1);
    fclose(fp2);
  } else {
    while((c = argv[1][len1]) != '\0') {
      if(c < 0x41 || (0x5a < c && c < 0x61) || 0x7a < c) continue;
      seq1 = (char*)realloc(seq1, sizeof(char) * (len1 + 1));
      seq1[len1] = c;
      ++len1;
    }

    while((c = argv[2][len2]) != '\0') {
      if(c < 0x41 || (0x5a < c && c < 0x61) || 0x7a < c) continue;
      seq2 = (char*)realloc(seq2, sizeof(char) * (len2 + 1));
      seq2[len2] = c;
      ++len2;
    }
  }

  score = (int**)malloc(sizeof(int*) * (len1 + 1));
  point = (int**)malloc(sizeof(int*) * (len1 + 1));

  for(i = 0; i < len1 + 1; ++i) {
    score[i] = (int*)malloc(sizeof(int) * (len2 + 1));
    point[i] = (int*)malloc(sizeof(int) * (len2 + 1));
  }

  for(i = 0; i < len1 + 1; ++i) {
    score[i][0] = i * gap;
    point[i][0] = 0;
  }

  for(j = 0; j < len2 + 1; ++j) {
    score[0][j] = j * gap;
    point[0][j] = 2;
  }

  for(i = 1; i < len1 + 1; ++i) {
    for(j = 1; j < len2 + 1; ++j) {
      int m_score = score[i-1][j-1] + match(seq1[i-1], seq2[j-1]);
      int l_score = score[i-1][j] + gap;
      int u_score = score[i][j-1] + gap;

      if(m_score >= l_score && m_score >= u_score) {
        score[i][j] = m_score;
        point[i][j] = 1;
      } else if (l_score >= m_score && l_score >= u_score) {
        score[i][j] = l_score;
        point[i][j] = 0;
      } else {
        score[i][j] = u_score;
        point[i][j] = 2;
      }
    }
  }

  --i;
  --j;

  while(i || j) {
    int t = point[i][j];

    aln1 = (char*)realloc(aln1, sizeof(char) * (len3 + 1));
    aln2 = (char*)realloc(aln2, sizeof(char) * (len4 + 1));

    if(t == 0) {
      aln1[len3] = seq1[i-1];
      aln2[len4] = '-';
      --i;
    } else if(t == 1) {
      aln1[len3] = seq1[i-1];
      aln2[len4] = seq2[j-1];
      --i;
      --j;
    } else {
      aln1[len3] = '-';
      aln2[len4] = seq2[j-1];
      --j;
    }

    ++len3;
    ++len4;
  }

  reverse(&aln1);
  reverse(&aln2);

  printf("%s\n%s\n", aln1, aln2);

  free(seq1);
  free(seq2);

  free(aln1);
  free(aln2);

  for(i = 0; i < len1 + 1; ++i) free(score[i]);
  for(i = 0; i < len1 + 1; ++i) free(point[i]);

  free(score);
  free(point);
}

void reverse(char** Pstring) {
  int i, c, len = 0;
  char* string = *(Pstring);

  for(i = 0; string[i] != '\0'; ++i) {
    ++len;
  }

  for(i = 0; i < len / 2; ++i) {
    c = string[i];
    string[i] = string[len - i - 1];
    string[len - i - 1] = c;
  }
}

int match(char a, char b) {
  /*
  ** Returns simple match.
  */
  return a == b ? 1 : -1;
}