动态规划算法_LIS_相似基因

题目背景

大家都知道，基因可以看作一个碱基对序列。它包含了 44 种核苷酸，简记作 A, C, G, T。生物学家正致力于寻找人类基因的功能，以利用于诊断疾病和发明药物。

在一个人类基因工作组的任务中，生物学家研究的是：两个基因的相似程度。因为这个研究对疾病的治疗有着非同寻常的作用。

题目描述

两个基因的相似度的计算方法如下：

对于两个已知基因，例如 AGTGATG 和 GTTAG，将它们的碱基互相对应。当然，中间可以加入一些空碱基 -，例如：

这样，两个基因之间的相似度就可以用碱基之间相似度的总和来描述，碱基之间的相似度如下表所示：

那么相似度就是：(−3)+5+5+(−2)+(−3)+5+(−3)+5=9(−3)+5+5+(−2)+(−3)+5+(−3)+5=9。因为两个基因的对应方法不唯一，例如又有：

相似度为：(−3)+5+5+(−2)+5+(−1)+5=14(−3)+5+5+(−2)+5+(−1)+5=14。规定两个基因的相似度为所有对应方法中，相似度最大的那个。

输入格式

共两行。每行首先是一个整数 𝑛n，表示基因序列的长度；隔一个空格后是一个基因序列，序列中只含 A,C,G,TA,C,G,T 四种字母。1≤𝑛≤1001≤n≤100。

输出格式

仅一行，即输入基因的相似度。

输入输出样例

输入

7 AGTGATG
5 GTTAG

输出

14

解题思路：

7 AGTGATG  第一个碱基序列 i,i<=7
5 GTTAG     第二个碱基序列 j,j<=5

碱基相似度d[i][j]：

（1）定义状态

f[i][j]：第一个碱基序列的前i位非空脱氧核苷酸与第二个碱基序列的前j位非空脱氧核苷酸最大相似度。

f[2][3]表示第一个碱基序列AG与第二个碱基序列GTT的最大相似度

（2）转移方程

每对碱基有3种可能的状态：

a）非空碱基  非空碱基

b）非空碱基  空碱基

c）空碱基  非空碱基

模型转移式：

//d(ai,b-)碱基相似度d[][]，ai表示第一个序列第i个位置字符，b-表示第二个序列，空碱基位置

f[i][j]= max(f[i-1][j-1]+d(ai,bj) ,f[i][j-1]+ d(a-,bj) +f[i-1][j]+ d(ai,b-))

递推过程：

f[1][1] 时，如下图三种状态，其中A和G相似度最大，f[1][1]=-2 。

A	-A	A-
G	G-	-G

f[1][2] 时，如下图三种状态，其中A和GT相似度最大，f[1][2]=-3 。

A-	-A
GT	GT

f[2][1] 时，如下图三种状态，其中AG和G相似度最大，f[2][1]=2。

AG	AG
G-	-G

f[3][1] 时，如下图三种状态，其中AGT和G相似度最大，f[3][1]=1。

AGT	AGT	AGT
G--	-G-	--G

f[2][2] 时，如下图三种状态，其中AG和GT相似度最大，f[2][2]=0 。

//d(GT)碱基相似度d[][]行列对应的字符 G T

即：max(f[1][1]+d(GT) ,f[2][1]+ d(-T) +f[1][2]+ d(G-))

      =max(-2+(-2),2+(-2),(-3)+(-1))

      =max(0,0,-4)

      =0

AG	AG-	A-G
GT	-GT	GT-

f[3][2] 时，如下图三种状态，其中AGT和GT相似度最大，f[3][2]=-5 。

即：max(f[2][1]+d(TT) ,f[3][1]+ d(-T) +f[2][2]+ d(T-))

    =max(2+5,1+(-1),0+(-1))

   =max(7,0,0)

   =7

AGT	AGT-	AGT
- GT	- G-T	GT-
f[2][1]+dTT	f[3][1]+ d-T	f[2][2]+ dT-

（3）递推顺序

  确定for 循环下标顺序，从小到大还是从大到小

（4）边界

确定边界值，如下标为0时的取值。

f[0][0]=0

f[i][0]= f[i-1][0]+d[i][5]

f[0][j]= f[0][j-1]+d[5][j]

（5）代码实现

#include <bits/stdc++.h>
using namespace std;
int a[101],b[101],f[101][101];
int d[6][6]=
{
	{0,0,0,0,0,0},
	{0,5,-1,-2,-1,-3},
	{0,-1,5,-3,-2,-4},
	{0,-2,-3,5,-2,-2},
	{0,-1,-2,-2,5,-1},
	{0,-3,-4,-2,-1,0}
	
} ;
int charVerInt(char c)
{
	if (c=='A') return 1;
	if (c=='C') return 2;
	if (c=='G') return 3;
	return 4;
}
int main(){
	int num1,num2;
	cin >> num1;
	for (int i=1;i<=num1; i++)
	{
		char t;
		cin >> t ;
		a[i]=charVerInt(t);
	}
		cin >> num2;
	for (int i=1;i<=num2; i++)
	{
		char t;
		cin >> t ;
		b[i]=charVerInt(t);
	}
	f[0][0]=0;
	for(int i=1;i<=num1;i++)
	{
		f[i][0]=f[i-1][0]+d[a[i]][5];
	}
	for(int j=1;j<=num2;j++)
	{
		f[0][j]=f[0][j-1]+d[5][b[j]];
	}
	for(int i=1;i<=num1;i++)
	{
		for(int j=1;j<=num2;j++)
		{
			f[i][j]=max(f[i-1][j-1]+d[a[i]][b[j]],max(f[i-1][j]+d[a[i]][5],f[i][j-1]+d[5][b[j]]));
		}
	}
	cout << f[num1][num2];
	return 0;
}