KMP算法的应用——碱基序列匹配

碱基序列匹配（ACM真题）

地理项目是IBM和国家地理学会的合作研究项目，从成千上万捐献的DNA分析地球上人类是如何繁衍的。作为一个IBM的研究人员，请你写一个程序找出给定的DNA片段之间的相同之处，使得对个体的调查相关联。一个DNA碱基序列是指把在分子中发现的氮基的序列给罗列出来。有四种氮基：腺嘌呤 (A)、胸腺嘧啶(T)、鸟嘌呤(G)和胞嘧啶(D)，例如，一个6碱基DNA序列可以表示为 TAGACC。给出一个DNA碱基序列的集合，确定在所有序列中都出现的最长的碱基序列。

输入格式:
输入的第一行给出了整数n，表示测试数据集合的数目。每个测试数据集合由下述两部分组成：一个正整数m(2≤m≤10)，给出数据集合中碱基序列的数目。m行，每行给出一个60碱基的碱基序列。

输出格式:
对于输入的每个测试数据集合的所有的碱基序列，输出最长的相同的碱基子序列。如果最长的相同的碱基子序列的长度小于3，则输出“no significant commonalities”来代替碱基子序列。如果相同最长长度的子序列有多个，则仅输出按字母排序的第一个。

输入样例

3
2
GATACCAGATACCAGATACCAGATACCAGATACCAGATACCAGATACCAGATACCAGATA
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
3
GATACCAGATACCAGATACCAGATACCAGATACCAGATACCAGATACCAGATACCAGATA
GATACTAGATACTAGATACTAGATACTAAAGGAAAGGGAAAAGGGGAAAAAGGGGGAAAA
GATACCAGATACCAGATACCAGATACCAAAGGAAAGGGAAAAGGGGAAAAAGGGGGAAAA
3
CATCATCATCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
ACATCATCATAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
AACATCATCATTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT

输出样例

no significant commonalities
AGATAC
CATCATCAT

思路：KMP算法遍历每一个子串
KMP基本思想：发生不匹配时，主串S的i不回溯，子串T的j回溯到next[j]对应位置的k上。
①用字符串数组存储每一组的碱基序列，数组的元素是一串碱基序列；
②从第一串碱基序列（字符串数组的第一个元素）入手，从最长的子串（碱基序列本身）开始与其他序列进行比对，直到找到共有的子串为止；
③暴力列举每一个子串进行比对，采用KMP算法进行比较。

代码如下：

#include<bits/stdc++.h>
using namespace std;

void getNext(string S,int* next)  //得到子串下面的数组
{
    int j,k;
    j=0;k=-1;
    next[0]=-1;  //子串0号元素下面数为-1
    while(j<(S.size()-1))  //对子串所有元素下面赋值
    {
        if(k==-1||S[j]==S[k])  //如果k回到了第一个元素或者第j个元素等于第k个元素
        {
            j++;k++;  //j++；k++；
            next[j]=k;  //子串第j个元素下面的数为k
        }
        else
            k=next[k];  //k为第子串第k个元素下面的数
    }
}

bool Compare(string T,string *S,int n)  //返回该子串是否是每一个序列的子串
{
    int a=T.size();  //得到子串T的长度
    int next[a];  //建立子串的数组下标
    getNext(T,next);  //给子串数组赋值
    int results[n];  //建立一个大小为n的数组判断子串是不是n-1个主串的公共子串
    for(int i=0;i<n;i++)
    {
        results[i]=0;  //给数组hhh全赋初值0
    }
    for(int l=1;l<n;l++)
    {
        int aa=S[l].size();  //得到主串的长度
        int i=0,j=0;
        while(i<aa)  //当主串下标没到达尾部时
        {
            if(j==-1||S[l][i]==T[j]){
            ++i;
            ++j;
        }
        else
            j=next[j];
        if(j==a){
            results[l]=1;
            break;
        }
        }
    }
    for(int i=1;i<n;i++)  //查看该子串是否为每一个主串下面的子串
    {
        if(results[i]!=1)
            return false;  //不是则返回false
    }
    return true;  //反之是则返回true
}

void Deal(string *aar,int n)
{
    string key="Z";  //假定最长公共序列key
    string try1;  //第一个碱基序列的每一个字串
    int w=0;
    for(int i=60;i>=3;i--)  //从最长字符串长度开始作为子串长度
    {
        if(w!=0&&i<w)
                {
                    cout<<key<<endl;return ;
                }
        for(int k=0;k<=60-i;k++)  //开始位置
        {
            try1=aar[0].substr(k,i);  //第一个碱基序列的一个字串
            //cout<<try1<<endl;
            if(Compare(try1,aar,n))  //查看是否为公共字串
            {
                w=i;
                if((try1.size()>=key.size())&&(try1<key))
                key=try1;
            }
        }
    }
    if(key.size()<3)
    cout<<"no significant commonalities"<<endl;
    else
        cout<<key<<endl;
}

int main()
{
    int N;
    cin>>N; //输入数据集合的数目N
    for(int z=0;z<N;z++)  //输入集合的每一组元素
    {
        int n;
        cin>>n;  //输入数据集合中碱基序列的数目n
        string jjsz[n];  //建立jjsz[n]数组存放每一组碱基序列
        for(int x=0;x<n;x++)
        {
            cin>>jjsz[x];  //存放每一个碱基序列
        }
        Deal(jjsz,n);  //开始处理
    }
    return 0;
}