序列比对

最新推荐文章于 2025-05-19 21:24:39 发布

原创最新推荐文章于 2025-05-19 21:24:39 发布 · 925 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#n2 #ios

优化理论专栏收录该内容

4 篇文章

订阅专栏

本文介绍了一种基于动态规划的DNA序列比对算法，并详细展示了如何通过C++实现该算法来评估两个DNA序列的相似度。算法首先计算序列间的匹配得分，然后回溯找出最佳比对路径，最后输出比对结果及执行耗时。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

#include<stdio.h>
#include<time.h>
#include<iostream.h>
#include<fstream.h>
#include <iomanip.h>

void reverse(char *p);//对字符串进行旋转

int k=-2;

int a[100][100];
char ms[100];

char s[100]="ATGGTGCACCTGACTCCTGAGGAGAAGTCTGCCGTTACTGCCCTGTGGGGCAAGGTGAACGTGGATTAAGTTGGTGGTGAGGCCCTGGGCAG";                  //"acgtcact";
char t[100]="ATGCTGACTGCTGAGGAGAAGGCTGCCGTCACCGGCTTCTGGGGCAAGGTGAAAGTGGATGAAGTTGGTGCTGAGGCCCTGGGCAG";                                                                                                 //"atgtacta";

char as[100]="";
char at[100]="";

int count=0;

int getLength(char *p)
{
    int i=0;
    while(*p!='/0')
    {
        i++;
        p++;
    }
    return i;
}

int score(char * p1,char * p2)
{
    if(*p1==*p2)
        return 1;
    else
        return -1;
}
int maxint(int n1,int n2,int n3)
{
    int temp=(n1>=n2)?n1:n2;
    temp=(temp>=n3)?temp:n3;
    return temp;
}

void contact(char *p,char *q)
{
    char * origin=p;
    while(*p)p++;
    *p=*q;
}

int align(char *s,char *t,int m,int n)
{
    for(int i=0;i<=m;i++)
    {
        a[i][0]=k*i;
    }
    for(int j=0;j<=n;j++)
    {
        a[0][j]=k*j;
    }
    for(i=1;i<=m;i++)
    {
        for(j=1;j<=n;j++)
        {
        int temp=score(s+i-1,t+j-1);
        a[i][j]=maxint(a[i-1][j]+k,a[i][j-1]+k,a[i-1][j-1]+temp);
        }
    }
    return a[m][n];
}

void traceBack(int i,int j)
{
    if((i==0)&&(j==0))return;
    else if((i>0)&&(j>0)&&(a[i-1][j-1]+score(s+i-1,t+j-1)==a[i][j]))
    {
        contact(as,s+(i-1));
        contact(at,t+(j-1));
        count++;
        traceBack(i-1,j-1);
    }
    else if((i>0)&&(a[i-1][j]+k==a[i][j]))
    {
        contact(as,s+(i-1));
        contact(at,"-");
        count++;
        traceBack(i-1,j);
    }
    else
    {
        contact(as,"-");
        contact(at,t+(j-1));
        count++;
        traceBack(i,j-1);
    }
}

void main()
{
    //时间开始
    time_t   start=clock();

    int m=getLength(s);
    int n=getLength(t);
    fstream fout("d://visualC//AlignOut.txt",ios::out);
    fout<<"m:"<<m<<endl;
    fout<<"n:"<<n<<endl;

    fout<<align(s,t,m,n)<<endl;

    for(int i=0;i<=m;i++)
    {
        for(int j=0;j<=n;j++)
        {
            fout<<setw(8)<<a[i][j]<<" ";
        }
        fout<<endl;
    }

    traceBack(m,n);

    reverse(as);
    reverse(at);

    fout<<"as串为："<<as<<endl;
    fout<<"ts串为："<<at<<endl;
//时间结束
time_t end=clock();
double dur=static_cast<double>(end-start)/CLOCKS_PER_SEC*1000;
cout<<"/n您的程序执行所耗费的时间为:"<<dur<<"毫秒。"<<endl;
    /*cout<<"m:"<<m<<endl;
    cout<<"n:"<<n<<endl;

    cout<<align(s,t,m,n)<<endl;

    for(int i=0;i<=m;i++)
    {
        for(int j=0;j<=n;j++)
        {
            cout<<a[i][j]<<" ";
        }
        cout<<endl;
    }

    traceBack(m,n);

    reverse(as);
    reverse(at);

    cout<<"as串为："<<as<<endl;
    cout<<"ts串为："<<at<<endl;*/

}

    void reverse(char *p)
    {
        cout<<"********"<<p<<endl;
        int n=getLength(p);
        for(int i=0,j=n-1;i<=j;i++,j--)
        {
            char temp=*(p+i);
            *(p+i)=*(p+j);
            *(p+j)=temp;
        }
    }