Needleman-Wusnch algorithm as specified on http://en.wikipedia.org/wiki/Needleman_wunsch

NeedlemanWunsh.cpp

/* Needleman-Wusnch algorithm as specified on http://en.wikipedia.org/wiki/Needleman_wunsch */

#include<algorithm>
#include<iostream>
#include<string>

using namespace std;

string calcSeqAlignment(int S[4][4], int d, string seqA, string seqB);

const int A = 0, G = 1, C = 2, T = 3;
int convToIndex['T'+1];

int main()
{
	convToIndex['A'] = A;
	convToIndex['G'] = G;
	convToIndex['C'] = C;
	convToIndex['T'] = T;

	int d;
	string seqA, seqB, seqAlignment;
	int S[4][4];

	cout<<"Enter the similarity matrix: "<<endl;
	for(int i=0; i<4; ++i)
		for(int j=0; j<4; ++j)
			cin>>S[i][j];

	cout<<"Enter the gap penalty: "<<endl;
	cin>>d;

	cout<<"Enter the sequence A: "<<endl;
	cin>>seqA;

	cout<<"Enter the sequence B: "<<endl;
	cin>>seqB;

	seqAlignment = calcSeqAlignment(S, d, seqA, seqB);
	cout<<"The alignment is: "<<seqAlignment;

	return 0;
}

string calcSeqAlignment(int S[4][4], int d, string seqA, string seqB)
{
	size_t m = seqA.length();
	size_t n = seqB.length();
	int params[3];
	string alignmentA = "", alignmentB = "", seqAlignment;

	int **F = new int*[m+1];
	for(size_t i=0; i<=m; ++i)
		F[i] = new int[n+1];

	for(size_t i=0; i<=m; ++i)
		F[i][0] = d*i;
	for(size_t j=0; j<=n; ++j)
		F[0][j] = d*j;

	for(size_t i=1; i<=m; ++i)
		for(size_t j=1; j<=n; ++j){
			params[0] = F[i-1][j-1] + S[convToIndex[(size_t) seqA[i-1]]][convToIndex[(size_t) seqB[j-1]]];
			params[1] = F[i-1][j] + d;
			params[2] = F[i][j-1] + d;

			F[i][j] = *(max_element(params, params+3));
		}
	cout<<"The score is: "<<F[m][n]<<endl;
	
	int i=m, j=n;
	int Score, ScoreDiag, ScoreUp, ScoreLeft;
	while(i>0 && j>0){
		Score = F[i][j];
		ScoreDiag = F[i-1][j-1];
		ScoreUp = F[i][j-1];
		ScoreLeft= F[i-1][j];

		if(Score == ScoreDiag + S[convToIndex[(size_t) seqA[i-1]]][convToIndex[(size_t) seqB[j-1]]]){
			alignmentA = seqA[i-1]+ alignmentA;
			alignmentB = seqB[j-1] + alignmentB;
			i--, j--;
		} else if(Score == ScoreLeft + d) {
			alignmentA = seqA[i-1] + alignmentA;
			alignmentB = "-" + alignmentB;
			i--;
		} else if(Score == ScoreUp + d) {
			alignmentA = "-" + alignmentA;
			alignmentB = seqB[j-1] + alignmentB;
			j--;
		} else {
		}
	}

	while(i>0){
		alignmentA = seqA[i-1] + alignmentA;
		alignmentB = "-" + alignmentB;
		i--;
	}

	while(j>0){
		alignmentA = "-" + alignmentA;
		alignmentB = seqB[j-1] + alignmentB;
		j--;
	}

	seqAlignment = "\n"+alignmentA+"\n"+alignmentB+"\n";
	return seqAlignment;
}