蛋白质序列

 1 #include<bits/stdc++.h>
 2 using namespace std;
 3 const int N=1e3+10,w1=2;
 4 
 5 int h[N][N];
 6 int x1,y1;
 7 string pre1,pre2;
 8 void dfs(int x,int y,string s1,string s2)
 9 {
10     if(h[x][y]==0) return;
11     if(s2[x]==s1[y])
12     {
13         pre1+=s1[y],pre2+=s2[x];
14         dfs(x-1,y-1,s1,s2);
15     }
16     else{
17         int a=h[x-1][y-1],b=h[x-1][y],c=h[x][y-1];
18         int maxv=max(a,max(b,c));
19         bool zs=0,up=0,z=0;
20         if(a==maxv) zs=1;
21         if(b==maxv) up=1;
22         if(c==maxv) z=1;
23         if(zs)
24         {
25         pre1+=s1[y],pre2+=s2[x];
26         dfs(x-1,y-1,s1,s2);
27         }
28         else if(up)
29         {
30             pre1+=s1[y],pre2+='_';
31             dfs(x-1,y,s1,s2);
32         }
33         else if(z)
34         {
35             pre1+='_',pre2+=s2[x];
36             dfs(x,y-1,s1,s2);
37         }
38     }
39 }
40 
41 double nsw(string s1,string s2)
42 {
43     
44     int m=s1.size(),n=s2.size();
45     
46     cout<<n<<' '<<m<<endl;
47     s1=" "+s1,s2=" "+s2;
48     int res=0;
49     for(int i=1;i<=n;i++)
50     for(int j=1;j<=m;j++)
51     {
52         int &x=h[i][j];
53         int t;
54         if(s1[j]==s2[i]) t=3;
55         else t=-3;
56         x=max(h[i-1][j-1]+t,max(h[i-1][j]-w1,max(h[i][j-1]-w1,0)));
57         if(res<x)
58         {
59             res=x;
60             x1=i,y1=j;
61         }
62     }
63     pre1="",pre2="";
64     dfs(x1,y1,s1,s2);
65     
66     reverse(pre1.begin(),pre1.end());
67     reverse(pre2.begin(),pre2.end());
68     
69     int cnt=0;
70     for(int i=0;i<pre1.size();i++)
71     if(pre1[i]==pre2[i]) cnt++;
72     
73     double ans=(cnt*1.0)/(sqrt(m*1.0)+sqrt(n*1.0));
74     
75     return ans;
76 }
77 int main()
78 {
79     string s1,s2;
80     printf("请输入第一个蛋白质序列: \n");
81     cin>>s1;
82     printf("请输入第二个蛋白质序列: \n");
83     cin>>s2;
84     
85     double ans=nsw(s1,s2);
86     ans+=nsw(s2,s1);
87     ans=ans/2.0;
88     
89     printf("这两个蛋白质的相似性为：\n");
90     printf("%.2f",ans);
91     return 0;
92 }