【算法刷题】poj 1080 Human Gene Function

首先将要看到如何运用动态编程查找两个 DNA 序列的最长公共子序列（longest common subsequence，LCS）。发现了新的基因序列的生物学家通常想知道该基因序列与其他哪个序列最相似。查找 LCS 是计算两个序列相似程度的一种方法：LCS 越长，两个序列越相似。

子序列中的字符与子字符串中的字符不同，它们不需要是连续的。例如，ACE 是 ABCDE 的子序列，但不是它的子字符串。请看下面两个 DNA 序列：

• S1 = DE>GCCCTAGCGDE>
• S2 = DE>GCGCAATGDE>

这两个序列的 LCS 是 GCCAG。（请注意，这仅是一个 LCS，而不是唯一的 LCS，因为可能存在其他长度相同的公共子序列。这种最优化问题和其他最优化问题的解可能不止一个。）

LCS 算法

首先，考虑如何递归地计算 LCS。令：

• C1 是 S1 最右侧的字符
• C2 是 S2 最右侧的字符
• S1' 是 S1 中 “切掉” C1 的部分
• S2' 是 S2 中 “切掉” C2 的部分

有三个递归子问题：

• L1 = LCS(S1', S2)
• L2 = LCS(S1, S2')
• L3 = LCS(S1', S2')

结果表明（而且很容易使人相信）原始问题的解就是下面三个子序列中最长的一个：

• L1
• L2
• 如果 C1 等于 C2，则为 L3 后端加上 C1 ，如果 C1 不等于 C2，则为 L3。

（基线条件（base case）是 S1 或 S2 为长度为零的字符串的情形。在这种情况下，S1 和 S2 的 LCS 显然是长度为零的字符串。）

状态转移方程为 

         LCS(S1, S2)=max(LCS(S1, S2'),     LCS(S1', S2),     C1==C2 ? LCS(S1', S2')+1 : LCS(S1', S2'));

但是，就像计算斐波纳契数的递归过程一样，这个递归解需要多次计算相同的子问题。可以证明，这种递归解法需要耗费指数级的时间。相比之下，这一问题的动态编程解法的运行时间是 Θ(mn)，其中 m 和 n 分别是两个序列的长度。

poj1080提供了相似度计算方法，各字母匹配权重不一，需要对算法进行修改。

 题目大意如下：对给定的两个字符串（其中只包含ACGT四个字母），以向需要的地方添加空格space(即表中'-')的方式使两个字符等长，判定两个字符串的匹配最大程度。其匹配方式如图：

                                         

那么，类似地，其状态转移方程为

（其中State表示俩字符串或俩字符之间的相似度，俩字符之间相似度可以查阅上表）

State(S1, S2)=max(State(S1, S2')+State(C2,space),     State(S1', S2)+State(C1,space),    State(S1', S2')+State(C1,C2));

此时要注意边界的处理

State(S1, 0)=State(S1', S2)+State(C1,space)；

State(0, S2)=State(S1, S2')+State(space,C2)；

状态转移表

参考文档

http://blog.csdn.net/leohxj/article/details/6003430

1080 poj

http://hi.baidu.com/ch2009120504/item/6ac388f7f4826b58932af293

//
//  main.cpp
//  proj
//
//  Created by yi chen on 13-10-10.
//  Copyright (c) 2013年 yi chen. All rights reserved.
//POJ 1080 Human Gene Function

#include <iostream>
#include <string.h>

using namespace std;

struct gene {
    int length;
    char seq[101];
};

const int scoreTable[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},
};

const char space='-';

//将ACGT对应到scoreTable坐标
int translate(char & a){
    int res=0;
    switch (a) {
        case 'A' : res = 0;break;
        case 'C' : res = 1;break;
        case 'G' : res = 2;break;
        case 'T' : res = 3;break;
        case '-' : res = 4;break;
        default:
            res=-1;
            break;
    }
    return res;
}

//取scoreTable值
int calculateScore(char a, char b){
    int i=translate(a);
    int j=translate(b);
    return scoreTable[i][j];
}

//估算gene序列的相似度
int evaluateSimi(gene g1, gene g2){
    //申请内存 存储动规运算的中间值 避免重复运算
    int mapRows=g1.length+1,mapCols=g2.length+1;
    int ** map= new int * [mapRows];
    for (int i=0; i<mapRows; i++) {
        map[i]=new int [mapCols];
        memset(map[i], 0, sizeof(int)*mapCols);
    }
    
    //边界处理
    for (int i=1; i<mapRows; i++) {
        map[i][0]=map[i-1][0]+calculateScore(g1.seq[i-1], space);
    }
    for (int j=1; j<mapCols; j++) {
        map[0][j]=map[0][j-1]+calculateScore(space, g2.seq[j-1]);
    }
    
    //状态转移
    for (int i=0; i<g1.length; i++) {
        for (int j=0; j<g2.length; j++) {
            int score_char1_char2  = map[i][j]  +calculateScore(g1.seq[i], g2.seq[j]);
            
            int max=score_char1_char2;
            
            int score_char1_space2 = map[i][j+1]+calculateScore(g1.seq[i], '-');
            if (score_char1_space2>max) {
                max=score_char1_space2;
            }
            
            int score_space1_char2 = map[i+1][j]+calculateScore('-', g2.seq[j]);
            if (score_space1_char2>max) {
                max=score_space1_char2;
            }
            
            
            map[i+1][j+1]=max;
            
        }
    }
    //返回结果 释放内存
    int res=map[g1.length][g2.length];
    for (int i=0; i<mapRows; i++) {
        delete [] map[i];
    }
    delete [] map;
    
    return res;
}

int main()
{
    
    int numOfCases = 0;
    gene geneHeader[2];
    cin>>numOfCases;
    for (int i=0; i<numOfCases; i++) {
        cin>>geneHeader[0].length>>geneHeader[0].seq;
        cin>>geneHeader[1].length>>geneHeader[1].seq;
        int res=evaluateSimi(geneHeader[0],geneHeader[1]);
        cout<<res<<endl;
    }
    
    return 0;
}