1 *************************************************************************
2 *************************************************************************
3 * This is a subroutine to compare two structures and find the
4 * superposition that has the maximum TM-score.
5 * Reference: Yang Zhang, Jeffrey Skolnick, Proteins 2004 57:702-10.
8 * L1--Length of the first structure
9 * (x1(i),y1(i),z1(i))--coordinates of i'th residue at the first structure
10 * n1(i)--Residue sequence number of i'th residue at the first structure
11 * L2--Length of the second structure
12 * (x2(i),y2(i),z2(i))--coordinates of i'th residue at the second structure
13 * n2(i)--Residue sequence number of i'th residue at the second structure
14 * TM--TM-score of the comparison
15 * Rcomm--RMSD of two structures in the common aligned residues
16 * Lcomm--Length of the common aligned regions
19 * 1, Always put native as the second structure, by which TM-score
21 * 2, The returned (x1(i),y1(i),z1(i)) are the rotated structure after
22 * TM-score superposition.
23 *************************************************************************
24 *************************************************************************
25 subroutine TMscore(L1,x1,y1,z1,n1,L2,x2,y2,z2,n2,TM,Rcomm,Lcomm)
27 PARAMETER(nmax=maxres)
28 common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax)
29 common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
31 common/align/n_ali,iA(nmax),iB(nmax)
32 common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score
33 dimension k_ali(nmax),k_ali0(nmax)
34 dimension L_ini(100),iq(nmax)
36 double precision score,score_max
37 dimension xa(nmax),ya(nmax),za(nmax)
39 dimension x1(nmax),y1(nmax),z1(nmax),n1(nmax)
40 dimension x2(nmax),y2(nmax),z2(nmax),n2(nmax)
43 double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)
44 double precision u(3,3),t(3),rms,drms !armsd is real
48 ********* convert input data ****************
64 ******************************************************************
65 * pickup the aligned residues:
66 ******************************************************************
70 if(nresA(i).eq.nresB(j))then
79 n_ali=k !number of aligned residues
82 c write(*,*)'There is no common residues in the input structures'
92 d0=1.24*(nseqB-15)**(1.0/3.0)-1.8
96 if(d0_search.gt.8)d0_search=8
97 if(d0_search.lt.4.5)d0_search=4.5
98 *** iterative parameters ----->
99 n_it=20 !maximum number of iterations
100 d_output=5 !for output alignment
101 n_init_max=6 !maximum number of L_init
104 if(n_ali.lt.4)L_ini_min=n_ali
107 L_ini(n_init)=n_ali/2**(n_init-1)
108 if(L_ini(n_init).le.L_ini_min)then
109 L_ini(n_init)=L_ini_min
114 L_ini(n_init)=L_ini_min
117 ******************************************************************
118 * find the maximum score starting from local structures superposition
119 ******************************************************************
120 score_max=-1 !TM-score
121 do 333 i_init=1,n_init
123 iL_max=n_ali-L_init+1
124 do 300 iL=1,iL_max !on aligned residues, [1,nseqA]
128 k=iL+i-1 ![1,n_ali] common aligned
139 call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2
140 if(i_init.eq.1)then !global superposition
145 xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
146 yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
147 zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
150 call score_fun !init, get scores, n_cut+i_ali(i) for iteration
151 if(score_max.lt.score)then
158 *** iteration for extending ---------------------------------->
164 m=i_ali(i) ![1,n_ali]
175 call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2
177 xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
178 yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
179 zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
181 call score_fun !get scores, n_cut+i_ali(i) for iteration
182 if(score_max.lt.score)then
189 if(it.eq.n_it)goto 302
193 if(i_ali(i).eq.k_ali(i))neq=neq+1
195 if(n_cut.eq.neq)goto 302
197 301 continue !for iteration
199 300 continue !for shift
200 333 continue !for initial length, L_ali/M
202 ******** return the final rotation ****************
205 m=k_ali0(i) !record of the best alignment
214 call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2
216 x1(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
217 y1(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
218 z1(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
222 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
226 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
227 c 1, collect those residues with dis<d;
228 c 2, calculate score_GDT, score_maxsub, score_TM
229 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
232 PARAMETER(nmax=maxres)
234 common/stru/xa(nmax),ya(nmax),za(nmax),xb(nmax),yb(nmax),zb(nmax)
235 common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
237 common/align/n_ali,iA(nmax),iB(nmax)
238 common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score
240 double precision score,score_max
243 21 n_cut=0 !number of residue-pairs dis<d, for iteration
246 i=iA(k) ![1,nseqA] reoder number of structureA
248 dis=sqrt((xa(i)-xb(j))**2+(ya(i)-yb(j))**2+(za(i)-zb(j))**2)
251 i_ali(n_cut)=k ![1,n_ali], mark the residue-pairs in dis<d
253 score_sum=score_sum+1/(1+(dis/d0)**2)
255 if(n_cut.lt.3.and.n_ali.gt.3)then
259 score=score_sum/float(nseqB) !TM-score
264 cccccccccccccccc Calculate sum of (r_d-r_m)^2 cccccccccccccccccccccccccc
265 c w - w(m) is weight for atom pair c m (given)
266 c x - x(i,m) are coordinates of atom c m in set x (given)
267 c y - y(i,m) are coordinates of atom c m in set y (given)
268 c n - n is number of atom pairs (given)
269 c mode - 0:calculate rms only (given)
270 c 1:calculate rms,u,t (takes longer)
271 c rms - sum of w*(ux+t-y)**2 over all atom pairs (result)
272 c u - u(i,j) is rotation matrix for best superposition (result)
273 c t - t(i) is translation vector for best superposition (result)
274 c ier - 0: a unique optimal superposition has been determined(result)
275 c -1: superposition is not unique but optimal
276 c -2: no result obtained because of negative weights w
277 c or all weights equal to zero.
278 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
279 subroutine u3b(w, x, y, n, mode, rms, u, t, ier)
280 double precision w(*), x(3,*), y(3,*)
283 double precision rms, u(3,3), t(3)
286 integer i, j, k, l, m1, m
287 integer ip(9), ip2312(4)
288 double precision r(3,3), xc(3), yc(3), wc
289 double precision a(3,3), b(3,3), e(3), rr(6), ss(6)
290 double precision e0, d, spur, det, cof, h, g
291 double precision cth, sth, sqrth, p, sigma
293 double precision sqrt3, tol, zero
295 data sqrt3 / 1.73205080756888d+00 /
297 data zero / 0.0d+00 /
298 data ip / 1, 2, 4, 2, 3, 5, 4, 5, 6 /
299 data ip2312 / 2, 3, 1, 2 /
321 if( n .lt. 1 ) return
324 if( w(m) .lt. 0.0 ) return
327 xc(i) = xc(i) + w(m)*x(i,m)
328 yc(i) = yc(i) + w(m)*y(i,m)
331 if( wc .le. zero ) return
339 e0=e0+w(m)*((x(i,m)-xc(i))**2+(y(i,m)-yc(i))**2)
340 d = w(m) * ( y(i,m) - yc(i) )
342 r(i,j) = r(i,j) + d*( x(j,m) - xc(j) )
347 det = r(1,1) * ( (r(2,2)*r(3,3)) - (r(2,3)*r(3,2)) )
348 & - r(1,2) * ( (r(2,1)*r(3,3)) - (r(2,3)*r(3,1)) )
349 & + r(1,3) * ( (r(2,1)*r(3,2)) - (r(2,2)*r(3,1)) )
357 rr(m) = r(1,i)*r(1,j) + r(2,i)*r(2,j) + r(3,i)*r(3,j)
361 spur = (rr(1)+rr(3)+rr(6)) / 3.0
362 cof = (((((rr(3)*rr(6) - rr(5)*rr(5)) + rr(1)*rr(6))
363 & - rr(4)*rr(4)) + rr(1)*rr(3)) - rr(2)*rr(2)) / 3.0
369 if( spur .le. zero ) goto 40
372 g = (spur*cof - det)/2.0 - spur*h
373 if( h .le. zero ) then
374 if( mode .eq. 0 ) then
382 if( d .lt. zero ) d = zero
383 d = datan2( dsqrt(d), -g ) / 3.0
384 cth = sqrth * dcos(d)
385 sth = sqrth*sqrt3*dsin(d)
386 e(1) = (spur + cth) + cth
387 e(2) = (spur - cth) + sth
388 e(3) = (spur - cth) - sth
390 if( mode .eq. 0 ) then
396 ss(1) = (d-rr(3)) * (d-rr(6)) - rr(5)*rr(5)
397 ss(2) = (d-rr(6)) * rr(2) + rr(4)*rr(5)
398 ss(3) = (d-rr(1)) * (d-rr(6)) - rr(4)*rr(4)
399 ss(4) = (d-rr(3)) * rr(4) + rr(2)*rr(5)
400 ss(5) = (d-rr(1)) * rr(5) + rr(2)*rr(4)
401 ss(6) = (d-rr(1)) * (d-rr(3)) - rr(2)*rr(2)
403 if( dabs(ss(1)) .ge. dabs(ss(3)) ) then
405 if( dabs(ss(1)) .lt. dabs(ss(6)) ) j = 3
406 else if( dabs(ss(3)) .ge. dabs(ss(6)) ) then
420 if( d .gt. zero ) d = 1.0 / dsqrt(d)
426 d = a(1,1)*a(1,3) + a(2,1)*a(2,3) + a(3,1)*a(3,3)
427 if ((e(1) - e(2)) .gt. (e(2) - e(3))) then
437 a(i,m1) = a(i,m1) - d*a(i,m)
440 if( p .le. tol ) then
443 if (p .lt. dabs(a(i,m))) cycle
449 p = dsqrt( a(k,m)**2 + a(l,m)**2 )
450 if( p .le. tol ) goto 40
461 a(1,2) = a(2,3)*a(3,1) - a(2,1)*a(3,3)
462 a(2,2) = a(3,3)*a(1,1) - a(3,1)*a(1,3)
463 a(3,2) = a(1,3)*a(2,1) - a(1,1)*a(2,3)
468 b(i,l) = r(i,1)*a(1,l) + r(i,2)*a(2,l) + r(i,3)*a(3,l)
471 if( d .gt. zero ) d = 1.0 / dsqrt(d)
476 d = b(1,1)*b(1,2) + b(2,1)*b(2,2) + b(3,1)*b(3,2)
480 b(i,2) = b(i,2) - d*b(i,1)
483 if( p .le. tol ) then
486 if( p .lt. dabs(b(i,1)) ) cycle
492 p = dsqrt( b(k,1)**2 + b(l,1)**2 )
493 if( p .le. tol ) goto 40
504 b(1,3) = b(2,1)*b(3,2) - b(2,2)*b(3,1)
505 b(2,3) = b(3,1)*b(1,2) - b(3,2)*b(1,1)
506 b(3,3) = b(1,1)*b(2,2) - b(1,2)*b(2,1)
510 u(i,j) = b(i,1)*a(j,1) + b(i,2)*a(j,2) + b(i,3)*a(j,3)
515 t(i) = ((yc(i) - u(i,1)*xc(1)) - u(i,2)*xc(2)) - u(i,3)*xc(3)
518 if( e(i) .lt. zero ) e(i) = zero
523 if( e(2) .le. (e(1) * 1.0d-05) ) ier = -1
526 if( sigma .lt. 0.0 ) then
528 if( (e(2) - e(3)) .le. (e(1) * 1.0d-05) ) ier = -1
530 d = (d + e(2)) + e(1)
533 if( rms .lt. 0.0 ) rms = 0.0