算法-动态规划-基因检测-最大相似度问题

题目描述

题目描述

大家都知道，基因可以看作一个碱基对序列。它包含了 4 种核苷酸，简记作 A,C,G,T 。生物学家正致力于寻找人类基因的功能，以利用于诊断疾病和发明药物。
在一个人类基因工作组的任务中，生物学家研究的是：两个基因的相似程度。因为这个研究对疾病的治疗有着非同寻常的作用。
两个基因的相似度的计算方法如下：
对于两个已知基因，例如 AGTGATG 和 GTTAG ，将它们的碱基互相对应。当然，中间可以加入一些空碱基 -，例如：

A	G	T	G	A	T	-	G
-	G	T	-	-	T	A	G

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

	A	C	G	T	-
A	5	-1	-2	-1	-3
C	-1	5	-3	-2	-4
G	-2	-3	5	-2	-2
T	-1	-2	-2	5	-1
-	-3	-4	-2	-1	*

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

A	G	T	G	A	T	G
-	G	T	T	A	-	G

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

输入格式

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

输出格式

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

输入样例

7 AGTGATG
5 GTTAG

输出样例

14

算法思路

1.该问题具有最优子结构性质和子问题重叠性质，故可使用动态规划算法求解。
2.递归结构：设有两个碱基序列
$$X=x_1{x}_2\cdot \cdot \cdot x_m，Y=y_1{y}_2\cdot \cdot \cdot y_n，$$其中xi, yj为字符A, G, C或T （1 ≤ i ≤ m，1 ≤ j ≤ n）。
设 dp[ i ][ j ] 为碱基序列X和Y的子序列$$\begin{gathered} Xi=x_1{x}_2\cdot \cdot \cdot x_i，Yj=y_1{y}_2\cdot \cdot \cdot y_j， \\ （X_0，Y_0表示空碱基序列） \end{gathered}$$对应的两段基因的最大相似度。
设 S[ x ][ y ] 是碱基x与y的相似度（若为‘-’则表示空碱基）。
当 i 或 j 为0时，dp[ i ][ j ] = 0。当 i 与 j 都不为0时，有以下三种情况：
1.序列 Xi-1 与 Yj-1 的最大相似度加上S[ xi ][ yj ]。
2.序列 Xi-1 与 Yj 的最大相似度加上S[ xi ][ - ]。
3.序列 Xi 与 Yj-1 的最大相似度加上S[ yj ][ - ]。
dp[ i ][ j ]的结果应为上述三种情况的最大值。即可按以下方法递归：
$$dp[i][j]=max\{\begin{array}{l}dp[i-1][j]+S[x_i][-]\\ dp[i][j-1]+S[y_j][-]\\ dp[i-1][j-1]+S[x_i][y_j]\end{array}$$

c语言代码

#include <stdio.h>
#include <string.h>
#define MAX_N 100

// 碱基对应相似度矩阵
const int similarity_matrix[5][5] = {
    { 5, -1, -2, -1, -3 },
    { -1, 5, -3, -2, -4 },
    { -2, -3, 5, -2, -2 },
    { -1, -2, -2, 5, -1 },
    { -3, -4, -2, -1, 0 }
};

// 获得不同碱基在矩阵中的位置 
int iGene(char a) {
	if(a == 'A') return 0;
	else if(a == 'C') return 1;
	else if(a == 'G') return 2;
	else if(a == 'T') return 3;
	else  return 4; 
}

int max(int a, int b, int c) {
    int max_value = a > b ? a : b;
    return max_value > c ? max_value : c;
}

int calculate_similarity(char* gene1, char* gene2, int n1, int n2) {
    // 声明并初始化dp数组
    int dp[MAX_N + 1][MAX_N + 1];
    memset(dp, 0, sizeof(dp));
    // 边界状态
    for (int i = 1; i <= n1; i++) {
        dp[i][0] = dp[i - 1][0] + similarity_matrix[iGene(gene1[i - 1])][4];
        
    }
    for (int i = 1; i <= n2; i++) {
        dp[0][i] = dp[0][i - 1] + similarity_matrix[iGene(gene2[i - 1])][4];
    }
    // 状态转移
    for (int i = 1; i <= n1; i++) {
        for (int j = 1; j <= n2; j++) {
            int case1 = dp[i - 1][j] + similarity_matrix[iGene(gene1[i - 1])][4];
            int case2 = dp[i][j - 1] + similarity_matrix[iGene(gene2[j - 1])][4];
            int case3 = dp[i - 1][j - 1] + similarity_matrix[iGene(gene1[i - 1])][iGene(gene2[j - 1])];
            dp[i][j] = max(case1, case2, case3);
        }
    }
    // 返回相似度最大值
    return dp[n1][n2];
}

int main() {
    int n1, n2;
    char gene1[MAX_N + 1], gene2[MAX_N + 1];
    scanf("%d %s %d %s", &n1, gene1, &n2, gene2);
    int similarity1 = calculate_similarity(gene1, gene2, n1, n2);
    int similarity2 = calculate_similarity(gene2, gene1, n2, n1);
    int similarity = similarity1 > similarity2 ? similarity1 : similarity2;
    printf("%d\n", similarity);
    return 0;
}