**************************************************************************
* This is a simple example for showing how to call the TM-score
* subroutine to calculate the TM-score superpositions.
**************************************************************************
program example
PARAMETER(nmax=3000)
dimension x1(nmax),y1(nmax),z1(nmax),n1(nmax)
dimension x2(nmax),y2(nmax),z2(nmax),n2(nmax)
********* read the structure1--------->
open(unit=10,file='pdb1',status='old')
read(10,*)L1
do i=1,L1
read(10,103)du,n1(i),du,x1(i),y1(i),z1(i)
enddo
close(10)
********* read the structure2--------->
open(unit=10,file='pdb2',status='old')
read(10,*)L2
do i=1,L2
read(10,103)du,n2(i),du,x2(i),y2(i),z2(i)
enddo
close(10)
103 format(A22,i4,A4,3F8.3)
********* calculate TM-score ---------->
call TMscore(L1,x1,y1,z1,n1,L2,x2,y2,z2,n2,TM,Rcomm,Lcomm)
do i=1,L1
write(*,*)i,n1(i),x1(i),y1(i),z1(i)
enddo
do i=1,L2
write(*,*)i,n2(i),x2(i),y2(i),z2(i)
enddo
write(*,*)'TMscore=',TM
write(*,*)'Number of residues in common=',Lcomm
write(*,*)'RMSD of the common residues=',Rcomm
stop
end
*************************************************************************
*************************************************************************
* This is a subroutine to compare two structures and find the
* superposition that has the maximum TM-score.
* Reference: Yang Zhang, Jeffrey Skolnick, Proteins 2004 57:702-10.
*
* Explanations:
* L1--Length of the first structure
* (x1(i),y1(i),z1(i))--coordinates of i'th residue at the first structure
* n1(i)--Residue sequence number of i'th residue at the first structure
* L2--Length of the second structure
* (x2(i),y2(i),z2(i))--coordinates of i'th residue at the second structure
* n2(i)--Residue sequence number of i'th residue at the second structure
* TM--TM-score of the comparison
* Rcomm--RMSD of two structures in the common aligned residues
* Lcomm--Length of the common aligned regions
*
* Note:
* 1, Always put native as the second structure, by which TM-score
* is normalized.
* 2, The returned (x1(i),y1(i),z1(i)) are the rotated structure after
* TM-score superposition.
*************************************************************************
*************************************************************************
subroutine TMscore(L1,x1,y1,z1,n1,L2,x2,y2,z2,n2,TM,Rcomm,Lcomm)
PARAMETER(nmax=3000)
common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax)
common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
common/para/d,d0
common/align/n_ali,iA(nmax),iB(nmax)
common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score
dimension k_ali(nmax),k_ali0(nmax)
dimension L_ini(100),iq(nmax)
common/scores/score
double precision score,score_max
dimension xa(nmax),ya(nmax),za(nmax)
dimension x1(nmax),y1(nmax),z1(nmax),n1(nmax)
dimension x2(nmax),y2(nmax),z2(nmax),n2(nmax)
ccc RMSD:
double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)
double precision u(3,3),t(3),rms,drms !armsd is real
data w /nmax*1.0/
ccc
********* convert input data ****************
nseqA=L1
do i=1,nseqA
xa(i)=x1(i)
ya(i)=y1(i)
za(i)=z1(i)
nresA(i)=n1(i)
enddo
nseqB=L2
do i=1,L2
xb(i)=x2(i)
yb(i)=y2(i)
zb(i)=z2(i)
nresB(i)=n2(i)
enddo
******************************************************************
* pickup the aligned residues:
******************************************************************
k=0
do i=1,nseqA
do j=1,nseqB
if(nresA(i).eq.nresB(j))then
k=k+1
iA(k)=i
iB(k)=j
goto 205
endif
enddo
205 continue
enddo
n_ali=k !number of aligned residues
Lcomm=n_ali
if(n_ali.lt.1)then
c write(*,*)'There is no common residues in the input structures'
TM=0
Rcomm=0
return
endif
************/////
* parameters:
*****************
*** d0------------->
d0=1.24*(nseqB-15)**(1.0/3.0)-1.8
if(d0.lt.0.5)d0=0.5
*** d0_search ----->
d0_search=d0
if(d0_search.gt.8)d0_search=8
if(d0_search.lt.4.5)d0_search=4.5
*** iterative parameters ----->
n_it=20 !maximum number of iterations
d_output=5 !for output alignment
n_init_max=6 !maximum number of L_init
n_init=0
L_ini_min=4
if(n_ali.lt.4)L_ini_min=n_ali
do i=1,n_init_max-1
n_init=n_init+1
L_ini(n_init)=n_ali/2**(n_init-1)
if(L_ini(n_init).le.L_ini_min)then
L_ini(n_init)=L_ini_min
goto 402
endif
enddo
n_init=n_init+1
L_ini(n_init)=L_ini_min
402 continue
******************************************************************
* find the maximum score starting from local structures superposition
******************************************************************
score_max=-1 !TM-score
do 333 i_init=1,n_init
L_init=L_ini(i_init)
iL_max=n_ali-L_init+1
do 300 iL=1,iL_max !on aligned residues, [1,nseqA]
LL=0
ka=0
do i=1,L_init
k=iL+i-1 ![1,n_ali] common aligned
r_1(1,i)=xa(iA(k))
r_1(2,i)=ya(iA(k))
r_1(3,i)=za(iA(k))
r_2(1,i)=xb(iB(k))
r_2(2,i)=yb(iB(k))
r_2(3,i)=zb(iB(k))
LL=LL+1
ka=ka+1
k_ali(ka)=k
enddo
call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2
if(i_init.eq.1)then !global superposition
armsd=dsqrt(rms/LL)
Rcomm=armsd
endif
do j=1,nseqA
xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
enddo
d=d0_search-1
call score_fun !init, get scores, n_cut+i_ali(i) for iteration
if(score_max.lt.score)then
score_max=score
ka0=ka
do i=1,ka0
k_ali0(i)=k_ali(i)
enddo
endif
*** iteration for extending ---------------------------------->
d=d0_search+1
do 301 it=1,n_it
LL=0
ka=0
do i=1,n_cut
m=i_ali(i) ![1,n_ali]
r_1(1,i)=xa(iA(m))
r_1(2,i)=ya(iA(m))
r_1(3,i)=za(iA(m))
r_2(1,i)=xb(iB(m))
r_2(2,i)=yb(iB(m))
r_2(3,i)=zb(iB(m))
ka=ka+1
k_ali(ka)=m
LL=LL+1
enddo
call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2
do j=1,nseqA
xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
enddo
call score_fun !get scores, n_cut+i_ali(i) for iteration
if(score_max.lt.score)then
score_max=score
ka0=ka
do i=1,ka
k_ali0(i)=k_ali(i)
enddo
endif
if(it.eq.n_it)goto 302
if(n_cut.eq.ka)then
neq=0
do i=1,n_cut
if(i_ali(i).eq.k_ali(i))neq=neq+1
enddo
if(n_cut.eq.neq)goto 302
endif
301 continue !for iteration
302 continue
300 continue !for shift
333 continue !for initial length, L_ali/M
******** return the final rotation ****************
LL=0
do i=1,ka0
m=k_ali0(i) !record of the best alignment
r_1(1,i)=xa(iA(m))
r_1(2,i)=ya(iA(m))
r_1(3,i)=za(iA(m))
r_2(1,i)=xb(iB(m))
r_2(2,i)=yb(iB(m))
r_2(3,i)=zb(iB(m))
LL=LL+1
enddo
call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2
do j=1,nseqA
x1(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
y1(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
z1(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
enddo
TM=score_max
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
return
END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c 1, collect those residues with dis<d;
c 2, calculate score_GDT, score_maxsub, score_TM
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine score_fun
PARAMETER(nmax=3000)
common/stru/xa(nmax),ya(nmax),za(nmax),xb(nmax),yb(nmax),zb(nmax)
common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
common/para/d,d0
common/align/n_ali,iA(nmax),iB(nmax)
common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score
common/scores/score
double precision score,score_max
d_tmp=d
21 n_cut=0 !number of residue-pairs dis<d, for iteration
score_sum=0 !TMscore
do k=1,n_ali
i=iA(k) ![1,nseqA] reoder number of structureA
j=iB(k) ![1,nseqB]
dis=sqrt((xa(i)-xb(j))**2+(ya(i)-yb(j))**2+(za(i)-zb(j))**2)
if(dis.lt.d_tmp)then
n_cut=n_cut+1
i_ali(n_cut)=k ![1,n_ali], mark the residue-pairs in dis<d
endif
score_sum=score_sum+1/(1+(dis/d0)**2)
enddo
if(n_cut.lt.3.and.n_ali.gt.3)then
d_tmp=d_tmp+.5
goto 21
endif
score=score_sum/float(nseqB) !TM-score
return
end
cccccccccccccccc Calculate sum of (r_d-r_m)^2 cccccccccccccccccccccccccc
c w - w(m) is weight for atom pair c m (given)
c x - x(i,m) are coordinates of atom c m in set x (given)
c y - y(i,m) are coordinates of atom c m in set y (given)
c n - n is number of atom pairs (given)
c mode - 0:calculate rms only (given)
c 1:calculate rms,u,t (takes longer)
c rms - sum of w*(ux+t-y)**2 over all atom pairs (result)
c u - u(i,j) is rotation matrix for best superposition (result)
c t - t(i) is translation vector for best superposition (result)
c ier - 0: a unique optimal superposition has been determined(result)
c -1: superposition is not unique but optimal
c -2: no result obtained because of negative weights w
c or all weights equal to zero.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine u3b(w, x, y, n, mode, rms, u, t, ier)
double precision w(*), x(3,*), y(3,*)
integer n, mode
double precision rms, u(3,3), t(3)
integer ier
integer i, j, k, l, m1, m
integer ip(9), ip2312(4)
double precision r(3,3), xc(3), yc(3), wc
double precision a(3,3), b(3,3), e(3), rr(6), ss(6)
double precision e0, d, spur, det, cof, h, g
double precision cth, sth, sqrth, p, sigma
double precision sqrt3, tol, zero
data sqrt3 / 1.73205080756888d+00 /
data tol / 1.0d-2 /
data zero / 0.0d+00 /
data ip / 1, 2, 4, 2, 3, 5, 4, 5, 6 /
data ip2312 / 2, 3, 1, 2 /
wc = zero
rms = zero
e0 = zero
do i=1, 3
xc(i) = zero
yc(i) = zero
t(i) = zero
do j=1, 3
r(i,j) = zero
u(i,j) = zero
a(i,j) = zero
if( i .eq. j ) then
u(i,j) = 1.0
a(i,j) = 1.0
end if
end do
end do
ier = -1
if( n .lt. 1 ) return
ier = -2
do m=1, n
if( w(m) .lt. 0.0 ) return
wc = wc + w(m)
do i=1, 3
xc(i) = xc(i) + w(m)*x(i,m)
yc(i) = yc(i) + w(m)*y(i,m)
end do
end do
if( wc .le. zero ) return
do i=1, 3
xc(i) = xc(i) / wc
yc(i) = yc(i) / wc
end do
do m=1, n
do i=1, 3
e0=e0+w(m)*((x(i,m)-xc(i))**2+(y(i,m)-yc(i))**2)
d = w(m) * ( y(i,m) - yc(i) )
do j=1, 3
r(i,j) = r(i,j) + d*( x(j,m) - xc(j) )
end do
end do
end do
det = r(1,1) * ( (r(2,2)*r(3,3)) - (r(2,3)*r(3,2)) )
& - r(1,2) * ( (r(2,1)*r(3,3)) - (r(2,3)*r(3,1)) )
& + r(1,3) * ( (r(2,1)*r(3,2)) - (r(2,2)*r(3,1)) )
sigma = det
m = 0;
do j=1, 3
do i=1, j
m = m+1
rr(m) = r(1,i)*r(1,j) + r(2,i)*r(2,j) + r(3,i)*r(3,j)
end do
end do
spur = (rr(1)+rr(3)+rr(6)) / 3.0
cof = (((((rr(3)*rr(6) - rr(5)*rr(5)) + rr(1)*rr(6))
& - rr(4)*rr(4)) + rr(1)*rr(3)) - rr(2)*rr(2)) / 3.0
det = det*det
do i=1, 3
e(i) = spur
end do
if( spur .le. zero ) goto 40
d = spur*spur
h = d - cof
g = (spur*cof - det)/2.0 - spur*h
if( h .le. zero ) then
if( mode .eq. 0 ) then
goto 50
else
goto 30
end if
end if
sqrth = dsqrt(h)
d = h*h*h - g*g
if( d .lt. zero ) d = zero
d = datan2( dsqrt(d), -g ) / 3.0
cth = sqrth * dcos(d)
sth = sqrth*sqrt3*dsin(d)
e(1) = (spur + cth) + cth
e(2) = (spur - cth) + sth
e(3) = (spur - cth) - sth
if( mode .eq. 0 ) then
goto 50
end if
do l=1, 3, 2
d = e(l)
ss(1) = (d-rr(3)) * (d-rr(6)) - rr(5)*rr(5)
ss(2) = (d-rr(6)) * rr(2) + rr(4)*rr(5)
ss(3) = (d-rr(1)) * (d-rr(6)) - rr(4)*rr(4)
ss(4) = (d-rr(3)) * rr(4) + rr(2)*rr(5)
ss(5) = (d-rr(1)) * rr(5) + rr(2)*rr(4)
ss(6) = (d-rr(1)) * (d-rr(3)) - rr(2)*rr(2)
if( dabs(ss(1)) .ge. dabs(ss(3)) ) then
j=1
if( dabs(ss(1)) .lt. dabs(ss(6)) ) j = 3
else if( dabs(ss(3)) .ge. dabs(ss(6)) ) then
j = 2
else
j = 3
end if
d = zero
j = 3 * (j - 1)
do i=1, 3
k = ip(i+j)
a(i,l) = ss(k)
d = d + ss(k)*ss(k)
end do
if( d .gt. zero ) d = 1.0 / dsqrt(d)
do i=1, 3
a(i,l) = a(i,l) * d
end do
end do
d = a(1,1)*a(1,3) + a(2,1)*a(2,3) + a(3,1)*a(3,3)
if ((e(1) - e(2)) .gt. (e(2) - e(3))) then
m1 = 3
m = 1
else
m1 = 1
m = 3
endif
p = zero
do i=1, 3
a(i,m1) = a(i,m1) - d*a(i,m)
p = p + a(i,m1)**2
end do
if( p .le. tol ) then
p = 1.0
do i=1, 3
if (p .lt. dabs(a(i,m))) cycle
p = dabs( a(i,m) )
j = i
end do
k = ip2312(j)
l = ip2312(j+1)
p = dsqrt( a(k,m)**2 + a(l,m)**2 )
if( p .le. tol ) goto 40
a(j,m1) = zero
a(k,m1) = -a(l,m)/p
a(l,m1) = a(k,m)/p
else
p = 1.0 / dsqrt(p)
do i=1, 3
a(i,m1) = a(i,m1)*p
end do
end if
a(1,2) = a(2,3)*a(3,1) - a(2,1)*a(3,3)
a(2,2) = a(3,3)*a(1,1) - a(3,1)*a(1,3)
a(3,2) = a(1,3)*a(2,1) - a(1,1)*a(2,3)
30 do l=1, 2
d = zero
do i=1, 3
b(i,l) = r(i,1)*a(1,l) + r(i,2)*a(2,l) + r(i,3)*a(3,l)
d = d + b(i,l)**2
end do
if( d .gt. zero ) d = 1.0 / dsqrt(d)
do i=1, 3
b(i,l) = b(i,l)*d
end do
end do
d = b(1,1)*b(1,2) + b(2,1)*b(2,2) + b(3,1)*b(3,2)
p = zero
do i=1, 3
b(i,2) = b(i,2) - d*b(i,1)
p = p + b(i,2)**2
end do
if( p .le. tol ) then
p = 1.0
do i=1, 3
if( p .lt. dabs(b(i,1)) ) cycle
p = dabs( b(i,1) )
j = i
end do
k = ip2312(j)
l = ip2312(j+1)
p = dsqrt( b(k,1)**2 + b(l,1)**2 )
if( p .le. tol ) goto 40
b(j,2) = zero
b(k,2) = -b(l,1)/p
b(l,2) = b(k,1)/p
else
p = 1.0 / dsqrt(p)
do i=1, 3
b(i,2) = b(i,2)*p
end do
end if
b(1,3) = b(2,1)*b(3,2) - b(2,2)*b(3,1)
b(2,3) = b(3,1)*b(1,2) - b(3,2)*b(1,1)
b(3,3) = b(1,1)*b(2,2) - b(1,2)*b(2,1)
do i=1, 3
do j=1, 3
u(i,j) = b(i,1)*a(j,1) + b(i,2)*a(j,2) + b(i,3)*a(j,3)
end do
end do
40 do i=1, 3
t(i) = ((yc(i) - u(i,1)*xc(1)) - u(i,2)*xc(2)) - u(i,3)*xc(3)
end do
50 do i=1, 3
if( e(i) .lt. zero ) e(i) = zero
e(i) = dsqrt( e(i) )
end do
ier = 0
if( e(2) .le. (e(1) * 1.0d-05) ) ier = -1
d = e(3)
if( sigma .lt. 0.0 ) then
d = - d
if( (e(2) - e(3)) .le. (e(1) * 1.0d-05) ) ier = -1
end if
d = (d + e(2)) + e(1)
rms = (e0 - d) - d
if( rms .lt. 0.0 ) rms = 0.0
return
end