Optimal Pair-wise Sequence Alignment
It implements Smith-Waterman with affine gap penalties. It requires at least one blank line between the two sequences. Ignores input lines with non-alphabetical characters
This program is provided as is with no warranty. The author is not responsible for any damage caused either directly or indirectly by using this program. Anybody is free to do whatever he/she wants with this program as long as this header section is preserved.
Created on 2005-02-01 by Roger Zhang ([email protected])

Last compiled under Linux with gcc-3

/*******************************************************
*     MYCPLUS Sample Code - https://www.mycplus.com     *
*                                                     *
*   This code is made available as a service to our   *
*      visitors and is provided strictly for the      *
*               purpose of illustration.              *
*                                                     *
* Please direct all inquiries to saqib at mycplus.com *
*******************************************************/

//====
// opsa.cpp
// - Optimal Pair-wise Sequence Alignment
// - implements Smith-Waterman with affine gap penalties
// - requires at least one blank line between the two sequences
// - ignores input lines with non-alphabetical characters
// Notes
// - This program is provided as is with no warranty.
// - The author is not responsible for any damage caused
//   either directly or indirectly by using this program.
// - Anybody is free to do whatever he/she wants with this
//   program as long as this header section is preserved.
// Created on 2005-02-01 by
// - Roger Zhang ([email protected])
// Modifications
// - Roger Zhang on 2005-02-18
//   - Hard coded the blosum matrix to remove external dependency
// -
// Last compiled under Linux with gcc-3
//====

#include 
#include 
#include 
#include "Array2D.hpp" // see http://cs.smu.ca/~r_zhang/a2h
#include "strUtil.hpp" // see http://cs.smu.ca/~r_zhang/suh

using namespace std;
using namespace RZ;

const double LARGE_NUMBER = 65536.;
const double GAP_OPENING_COST = 10.;
const double GAP_EXTENSION_COST = .1;
const double NEW_GAP_COST = GAP_OPENING_COST + GAP_EXTENSION_COST;

const signed char BLOSUM[][25] = { // the blosum 62 scoring matrix
   {4, 0, 0, -2, -1, -2, 0, -2, -1, 0, -1, -1, // A
    -1, -2, 0, -1, -1, -1, 1, 0, 0, 0, -3, 0, -2},
   {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
   {0, 0, 9, -3, -4, -2, -3, -3, -1, 0, -3, -1, // C
    -1, -3, 0, -3, -3, -3, -1, -1, 0, -1, -2, 0, -2},
   {-2, 0, -3, 6, 2, -3, -1, -1, -3, 0, -1, -4, // D
    -3, 1, 0, -1, 0, -2, 0, -1, 0, -3, -4, 0, -3},
   {-1, 0, -4, 2, 5, -3, -2, 0, -3, 0, 1, -3, // E
    -2, 0, 0, -1, 2, 0, 0, -1, 0, -2, -3, 0, -2},
   {-2, 0, -2, -3, -3, 6, -3, -1, 0, 0, -3, 0, // F
    0, -3, 0, -4, -3, -3, -2, -2, 0, -1, 1, 0, 3},
   {0, 0, -3, -1, -2, -3, 6, -2, -4, 0, -2, -4, // G
    -3, 0, 0, -2, -2, -2, 0, -2, 0, -3, -2, 0, -3},
   {-2, 0, -3, -1, 0, -1, -2, 8, -3, 0, -1, -3, // H
    -2, 1, 0, -2, 0, 0, -1, -2, 0, -3, -2, 0, 2},
   {-1, 0, -1, -3, -3, 0, -4, -3, 4, 0, -3, 2, // I
    1, -3, 0, -3, -3, -3, -2, -1, 0, 3, -3, 0, -1},
   {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
   {-1, 0, -3, -1, 1, -3, -2, -1, -3, 0, 5, -2, // K
    -1, 0, 0, -1, 1, 2, 0, -1, 0, -2, -3, 0, -2},
   {-1, 0, -1, -4, -3, 0, -4, -3, 2, 0, -2, 4, // L
    2, -3, 0, -3, -2, -2, -2, -1, 0, 1, -2, 0, -1},
   {-1, 0, -1, -3, -2, 0, -3, -2, 1, 0, -1, 2, // M
    5, -2, 0, -2, 0, -1, -1, -1, 0, 1, -1, 0, -1},
   {-2, 0, -3, 1, 0, -3, 0, 1, -3, 0, 0, -3, // N
    -2, 6, 0, -2, 0, 0, 1, 0, 0, -3, -4, 0, -2},
   {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
   {-1, 0, -3, -1, -1, -4, -2, -2, -3, 0, -1, -3, // P
    -2, -2, 0, 7, -1, -2, -1, -1, 0, -2, -4, 0, -3},
   {-1, 0, -3, 0, 2, -3, -2, 0, -3, 0, 1, -2, // Q
    0, 0, 0, -1, 5, 1, 0, -1, 0, -2, -2, 0, -1},
   {-1, 0, -3, -2, 0, -3, -2, 0, -3, 0, 2, -2, // R
    -1, 0, 0, -2, 1, 5, -1, -1, 0, -3, -3, 0, -2},
   {1, 0, -1, 0, 0, -2, 0, -1, -2, 0, 0, -2, // S
    -1, 1, 0, -1, 0, -1, 4, 1, 0, -2, -3, 0, -2},
   {0, 0, -1, -1, -1, -2, -2, -2, -1, 0, -1, -1, // T
    -1, 0, 0, -1, -1, -1, 1, 5, 0, 0, -2, 0, -2},
   {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
   {0, 0, -1, -3, -2, -1, -3, -3, 3, 0, -2, 1, // V
    1, -3, 0, -2, -2, -3, -2, 0, 0, 4, -3, 0, -1},
   {-3, 0, -2, -4, -3, 1, -2, -2, -3, 0, -3, -2, // W
    -1, -4, 0, -4, -2, -3, -3, -2, 0, -3, 11, 0, 2},
   {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
   {-2, 0, -2, -3, -2, 3, -3, 2, -1, 0, -2, -1, // Y
    -1, -2, 0, -3, -1, -2, -2, -2, 0, -1, 2, 0, 7}
};

void read(string &sequence)
{
   string line;

   while (getline(cin, line)) {
      trim(line);
      if (line.empty()) {
         return;
      }
      for (int i = 0, n = line.length(); i < n; i++) {
         if (!isalpha(line[i] = toupper(line[i]))) {
            if (!sequence.empty()) {
               return;
            }
            line = "";
            break;
         }
      }
      sequence += line;
   }
}

double max(double x, double y)
{
   return x > y ? x : y;
}

double max(double x, double y, double z)
{
   return x > y ? max(x, z) : max(y, z);
}

double alignment(string &s1, string &s2)
{
   int n = s1.length() + 1, m = s2.length() + 1, i, j;

   Array2D r(n, m), t(n, m), s(n, m);

   //====
   // initialization

   r[0][0] = t[0][0] = s[0][0] = 0;

   for (i = 1; i < n; i++) {
      r[i][0] = -LARGE_NUMBER;
      s[i][0] = t[i][0] = -GAP_OPENING_COST - i * GAP_EXTENSION_COST;
   }

   for (j = 1; j < m; j++) {
      t[0][j] = -LARGE_NUMBER;
      s[0][j] = r[0][j] = -GAP_OPENING_COST - j * GAP_EXTENSION_COST;
   }

   //====
   // Smith-Waterman with affine gap costs

   for (i = 1; i < n; i++) {
      for (j = 1; j < m; j++) {
         r[i][j] =
            max(r[i][j - 1] - GAP_EXTENSION_COST, s[i][j - 1] - NEW_GAP_COST);
         t[i][j] =
            max(t[i - 1][j] - GAP_EXTENSION_COST, s[i - 1][j] - NEW_GAP_COST);
         s[i][j] =
            max(
               s[i - 1][j - 1] + BLOSUM[s1[i - 1] - 'A'][s2[j - 1] - 'A'],
               r[i][j], t[i][j]
            );
      }
   }

   //====
   // back tracking

   i = n - 1, j = m - 1;

   while (i > 0 || j > 0) {
      if (s[i][j] == r[i][j]) {
         s1.insert(i, 1, '-');
         j--;
      } else if (s[i][j] == t[i][j]) {
         s2.insert(j, 1, '-');
         i--;
      } else {
         i--, j--;
      }
   }

   //====
   // final score

   return s[s.height() - 1][s.width() - 1];
}

int main()
{
   string sequence1, sequence2;

   read(sequence1), read(sequence2);

   double score = alignment(sequence1, sequence2);

   cout << sequence1 << "\n\n" << sequence2 << "\n\nScore: " << score << endl;

   return 0;
}