Demo entry 6658210

Sequence Alignment Using Dynamic Programming

   

Submitted by anonymous on Nov 06, 2017 at 14:34
Language: C++. Code size: 3.9 kB.

#include <iostream>
using namespace std;

void get_sq(char *s, int len);
void display_sq(char *s, int len);
int max(int a, int b, int c);
int cal_value(char *s1, char *s2, int **map, int i, int j);
void display_map(int **a,int m, int n);
int *find_max(int **a, int m, int n);
int **find_pre(int **a, int m, int n, int *cur);
void display_rs(int **a, int m, int n, int *cur);

char s1[12] = {' ','g','C','O','E','L','A','C','A','N','T','H'};
char s2[9] = {' ','g','P','E','L','I','C','A','N'};

int main()
{
	/*
	//get the length of sequence
	int len1, len2;
	cout << "Sequence1's length:";
	cin >> len1;
	cout << "Sequence2's length:";
	cin >> len2;
	
	char *s1 = new char [len1+1];
	char *s2 = new char [len2+1];
	
	cout << "Sequence1: ";
	char c=cin.get();
	get_sq(s1, len1);
	cout << "Sequence2: ";
	get_sq(s2, len2);
	*/
	
	//display_sq(s1, len1);
	//display_sq(s2, len2);
	int m=9, n=12;
	//init the map
	int ** map;
	map = new int *[m];
	for (int i=0; i<m; i++) map[i] = new int [n];
	
	int val = 0;
	for (int i=1; i<m; i++){
		map[i][1] = val;
		val--;
	} 
	
	val = 0;
	for (int j=1; j<n; j++){
		map[1][j] = val;
		val--;
	}
	
	//obtain the map
	for (int i=2; i<m; i++){
		for (int j=2; j<n; j++){
			map[i][j] = cal_value(s2,s1,map,i,j);
		}
	}
	
	//display the map
	display_map(map,m,n); 
	
	//find the max value in the map
	int *max_id; 
	max_id = new int[2];
	max_id = find_max(map, m, n);
	
	display_rs(map, m, n, max_id);
} 


//max
int max(int a, int b, int c)
{
	int tmp;
	if (a > b) tmp =a;
	else tmp = b;
	if (c > tmp) tmp =c;
	return tmp;
}

//calculate value
int cal_value(char *s1, char *s2, int **map, int i, int j)
{
	int a, b, c;
	
	a = map[i-1][j] -1;
	b = map[i][j-1] -1;
	
	if (s1[i] == s2[j]) c = map[i-1][j-1]+1;
	else c = map[i-1][j-1]-1;
	
	return max(a,b,c);
}

//display map
void display_map(int **a,int m, int n)
{
	for(int i=1; i<m; i++){
		for(int j=1; j<n; j++) cout << a[i][j] << '\t';
		cout<<endl;
	}
}

//find the max value in the map
int *find_max(int **a, int m, int n)
{
	int *b;
	b = new int[2];
	int max = a[1][1];
	for (int i=1; i<m; i++){
		for (int j=1; j<n; j++){
			if (a[i][j] >= max) {
				max = a[i][j];
				//cout<<i<<' '<< j<<endl;
				b[0] = i; 
				b[1] = j;
			} 
		}
	}
	return b;
}

//find pre-node
int **find_pre(int **a, int m, int n, int *cur)
{
	int **pre;
	pre = new int*[2];
	for (int i=0; i<3; i++) pre[i] = new int[2];
	
	
	for (int i=0; i<3; i++){
		for (int j=0; j<2; j++) {
			pre[i][j] = 0;
		}
	}
			
	if ((a[cur[0]-1][cur[1]]-a[cur[0]][cur[1]]) == 1) {
		pre[1][0] = cur[0]-1;
		pre[1][1] = cur[1];
	}
	
	
	if ((a[cur[0]][cur[1]-1]-a[cur[0]][cur[1]]) == 1){
		pre[2][0] = cur[0];
		pre[2][1] = cur[1]-1;
	}
	
	if ((a[cur[0]-1][cur[1]-1]-a[cur[0]][cur[1]]) == 1){
		pre[0][0] = cur[0]-1;
		pre[0][1] = cur[1]-1;
	}
	
	if ((a[cur[0]-1][cur[1]-1]-a[cur[0]][cur[1]]) == -1){
		pre[0][0] = cur[0]-1;
		pre[0][1] = cur[1]-1;
	}
	
	return pre;
}

//display the result(¿ÉÒÔͨ¹ý¹¹½¨¶à²æÊ÷À´±íʾ½á¹û£¬È»ºó±éÀúÏÔʾ£¬µ«ÊÇÍüÁËÔõô²Ù×÷ÁË)
void display_rs(int **a, int m, int n, int *cur)
{
	int **pre;
	char *rs;
	rs = new char[9];
	
	while(1){
		pre = find_pre(a, m, n, cur);		
	
		for (int i=0; i<3; i++){
			if (pre[i][0]!=0) {
				if (i == 2) rs[9-cur[0]] = '-';
				else rs[9-cur[0]] = s2[cur[0]];
				cur[0] = pre[i][0];
				cur[1] = pre[i][1];
			}
		}
		if(s2[cur[0]]=='g') break;
	}
	
	cout<<endl;
	cout << rs[8] <<' '; 
	for (int i=7; i>0; i--)
	{
			cout << rs[i] << " ";
	}
}

// read the sequences
void get_sq(char *s, int len)
{	
	char c;
	s[0] = 'g';
	int i = 1;
	
	while((c = cin.get()) != '\n')
	{
		s[i] = c;
		i++;
	}
}

// display the sequences
void display_sq(char *s, int len)
{
	for (int i=0;i<=len;i++) {
		cout << s[i] << ' ';
	}
	cout<<endl;
}

This snippet took 0.02 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).