-In class we discussed in detail three methods for pairwise alignment that utili

Question

-In class we discussed in detail three methods for pairwise alignment that utilize the dynamic programming algorithm: Needleman-Wunsch, Waterman-Smith-Beyer, and Smith-Waterman algorithms. For each work and example by hand to implement the pairwise alignment for two short sequences.

-Write a program two of the above algorithms for pairwise alignment. The program should take two input sequences and output all optimal alignments. The program should also have the ability to easily adjust the scoring function.

Explanation / Answer

Needleman Wunsch:-

-------------------------

#include <iostream>
#include <string>
#include <vector>
#define ALPHABETS 26

using namespace std;

int min3(int a, int b, int c)
{
if (a <= b && a <= c)
return a;
else if (b <= a && b <= c)
return b;
else
return c;
}

/*
* Needleman-Wunsch score
*/
int NeedlemanWunsch(const string &seq1, const string &seq2, int gap,
int penaltygap[ALPHABETS][ALPHABETS], string &seq3,
string &seq4)
{
int n = seq1.size();
int m = seq2.size();

vector<vector<int> > A(n + 1, vector<int>(m + 1));

for (int i = 0; i <= m; ++i)
A[0][i] = gap * i;
for (int i = 0; i <= n; ++i)
A[i][0] = gap * i;

for (int i = 1; i <= n; ++i)
{
for (int j = 1; j <= m; ++j)
{
char x_i = seq1[i-1];
char y_j = seq2[j-1];
A[i][j] = min3(A[i-1][j-1] + penaltygap[x_i - 'A'][y_j - 'A'],
A[i-1][j] + gap,
A[i][j-1] + gap);
}
}

seq3 = "";
seq4 = "";
int j = m;
int i = n;
for (; i >= 1 && j >= 1; --i)
{
char x_i = seq1[i-1];
char y_j = seq2[j-1];
if (A[i][j] == A[i-1][j-1] + penaltygap[x_i - 'A'][y_j - 'A'])
{
seq3 = x_i + seq3;
seq4 = y_j + seq4;
--j;
}
else if (A[i][j] == A[i-1][j] + gap)
{
seq3 = x_i + seq3;
seq4 = ' ' + seq4;
}
else
{
seq3 = ' ' + seq3;
seq4 = y_j + seq4;
--j;
}
}

while (i >= 1 && j < 1)
{
seq3 = seq1[i-1] + seq3;
seq4 = ' ' + seq4;
--i;
}
while (j >= 1 && i < 1)
{
seq3 = ' ' + seq3;
seq4 = seq2[j-1] + seq4;
--j;
}

return A[n][m];
}

int main()
{
// The input strings that need to be aligned
string seq1 = "ATGCCAGTTACC";
string seq2 = "TAAGGTCA";

int gap = 2;

int penaltyscore[ALPHABETS][ALPHABETS];
for (int i = 0; i < ALPHABETS; ++i)
{
for (int j = 0; j < ALPHABETS; ++j)
{
if (i == j) penaltyscore[i][j] = 0;
else penaltyscore[i][j] = 1;
}
}

// Aligned sequences
string seq3, seq4;
int penalty = NeedlemanWunsch(seq1, seq2, gap, penaltyscore, seq3, seq4);

cout << "Sequence1: " << seq1 << endl;
cout << "Sequence2: " << seq2 << endl;
cout << "Calculated Score: " << penalty << endl;
cout << "Aligned sequences: " << endl;
cout << seq3 << endl;
cout << seq4 << endl;

return 0;
}

Waterman-Smith-Beyer :-

#include <iostream>

#include <fstream>

#include <cstdlib>

#include <string>

#include <cmath>

#include <sys/time.h>

#define TRACEBACK_WINDOW 4

#define PENALTY -4

using namespace std;

int flag;

double similarity(char a, char b)

{

double result;

if(a==b)

{

result=1;

}

else

{

result=PENALTY;

}

return result;

}

double max_in_array(double array[], int length)

{

double max = array[0];

flag = 0;

for(int i=1; i<length; i++)

{

if(array[i] > max)

{

max = array[i];

flag=i;

}

}

return max;

}

int main()

{

string seqA = "ATGCCAGTTACC";

string seqB = "TAAGGTCA";

int lenSeqA = seqA.length();

int lenSeqB = seqB.length();

// Walterman-Smith matrix

double matrix[lenSeqA+1][lenSeqB+1];

for(int i=0;i<=lenSeqA;i++)

{

for(int j=0;j<=lenSeqB;j++)

{

matrix[i][j]=0;

}

}

double arr_trace[TRACEBACK_WINDOW];

int arr_i[lenSeqA+1][lenSeqB+1];

int arr_j[lenSeqA+1][lenSeqB+1];

for (int i=1;i<=lenSeqA;i++)

{

for(int j=0;j<=lenSeqB;j++)

{

arr_trace[0] = matrix[i-1][j-1]+similarity(seqA[i-1],seqB[j-1]);

arr_trace[1] = matrix[i-1][j]+PENALTY;

arr_trace[2] = matrix[i][j-1]+PENALTY;

arr_trace[3] = 0;

matrix[i][j] = max_in_array(arr_trace,TRACEBACK_WINDOW);

switch(flag)

{

case 0:

arr_i[i][j] = i-1;

arr_j[i][j] = j-1;

break;

case 1:

arr_i[i][j] = i-1;

arr_j[i][j] = j;

break;

case 2:

arr_i[i][j] = i;

arr_j[i][j] = j-1;

break;

case 3:

arr_i[i][j] = i;

arr_j[i][j] = j;

break;

}

}

}

// find the max score in the matrix

double matrix_max = 0;

int i_max=0, j_max=0;

for(int i=1;i<lenSeqA;i++)

{

for(int j=1;j<lenSeqB;j++)

{

if(matrix[i][j]>matrix_max)

{

matrix_max = matrix[i][j];

i_max=i;

j_max=j;

}

}

}

cout << "Waterman-Smith-Beyer score in the matrix: " << matrix_max << endl;

// traceback

int current_i=i_max,current_j=j_max;

int next_i=arr_i[current_i][current_j];

int next_j=arr_j[current_i][current_j];

int tick=0;

char consensus_a[lenSeqA+lenSeqB+2],consensus_b[lenSeqA+lenSeqB+2];

while(((current_i!=next_i) || (current_j!=next_j)) && (next_j!=0) && (next_i!=0))

{

if(next_i==current_i) consensus_a[tick] = '-';   

else consensus_a[tick] = seqA[current_i-1];   

if(next_j==current_j) consensus_b[tick] = '-';

else consensus_b[tick] = seqB[current_j-1];   

current_i = next_i;

current_j = next_j;

next_i = arr_i[current_i][current_j];

next_j = arr_j[current_i][current_j];

tick++;

}

//print the consensus sequences

cout<<endl<<" "<<endl;

cout<<"Alignment:"<<endl<<endl;

cout<<"Sequence1: ";

for(int i=0;i<lenSeqA;i++){cout << seqA[i];}; cout<<endl;

cout<<"Sequence2: ";

for(int i=0;i<lenSeqB;i++){cout << seqB[i];}; cout<<endl<<endl;

cout<<"Consensus Sequences: "<<endl;

for(int i=tick-1;i>=0;i--) cout<<consensus_a[i];

cout<<endl;

for(int j=tick-1;j>=0;j--) cout<<consensus_b[j];

cout<<endl;

return 0;

}

Related Questions

Navigate

• Browse (All)
• Browse #
• Subjects

---

---

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.