source/unres/src_CSA/TMscore_subroutine.f

   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.
   6 *
   7 *     Explanations:
   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
  17 *
  18 *     Note:
  19 *     1, Always put native as the second structure, by which TM-score
  20 *        is normalized.
  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)
  26       include 'DIMENSIONS'
  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
  30       common/para/d,d0
  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)
  35       common/scores/score
  36       double precision score,score_max
  37       dimension xa(nmax),ya(nmax),za(nmax)
  38
  39       dimension x1(nmax),y1(nmax),z1(nmax),n1(nmax)
  40       dimension x2(nmax),y2(nmax),z2(nmax),n2(nmax)
  41
  42 ccc   RMSD:
  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
  45       data w /nmax*1.0/
  46 ccc
  47
  48 ********* convert input data ****************
  49       nseqA=L1
  50       do i=1,nseqA
  51          xa(i)=x1(i)
  52          ya(i)=y1(i)
  53          za(i)=z1(i)
  54          nresA(i)=n1(i)
  55       enddo
  56       nseqB=L2
  57       do i=1,L2
  58          xb(i)=x2(i)
  59          yb(i)=y2(i)
  60          zb(i)=z2(i)
  61          nresB(i)=n2(i)
  62       enddo
  63
  64 ******************************************************************
  65 *     pickup the aligned residues:
  66 ******************************************************************
  67       k=0
  68       do i=1,nseqA
  69          do j=1,nseqB
  70             if(nresA(i).eq.nresB(j))then
  71                k=k+1
  72                iA(k)=i
  73                iB(k)=j
  74                goto 205
  75             endif
  76          enddo
  77  205     continue
  78       enddo
  79       n_ali=k                   !number of aligned residues
  80       Lcomm=n_ali
  81       if(n_ali.lt.1)then
  82 c         write(*,*)'There is no common residues in the input structures'
  83          TM=0
  84          Rcomm=0
  85          return
  86       endif
  87
  88 ************/////
  89 *     parameters:
  90 *****************
  91 ***   d0------------->
  92       d0=1.24*(nseqB-15)**(1.0/3.0)-1.8
  93       if(d0.lt.0.5)d0=0.5
  94 ***   d0_search ----->
  95       d0_search=d0
  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
 102       n_init=0
 103       L_ini_min=4
 104       if(n_ali.lt.4)L_ini_min=n_ali
 105       do i=1,n_init_max-1
 106          n_init=n_init+1
 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
 110             goto 402
 111          endif
 112       enddo
 113       n_init=n_init+1
 114       L_ini(n_init)=L_ini_min
 115  402  continue
 116
 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
 122         L_init=L_ini(i_init)
 123         iL_max=n_ali-L_init+1
 124         do 300 iL=1,iL_max      !on aligned residues, [1,nseqA]
 125            LL=0
 126            ka=0
 127            do i=1,L_init
 128               k=iL+i-1          ![1,n_ali] common aligned
 129               r_1(1,i)=xa(iA(k))
 130               r_1(2,i)=ya(iA(k))
 131               r_1(3,i)=za(iA(k))
 132               r_2(1,i)=xb(iB(k))
 133               r_2(2,i)=yb(iB(k))
 134               r_2(3,i)=zb(iB(k))
 135               LL=LL+1
 136               ka=ka+1
 137               k_ali(ka)=k
 138            enddo
 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
 141               armsd=dsqrt(rms/LL)
 142               Rcomm=armsd
 143            endif
 144            do j=1,nseqA
 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)
 148            enddo
 149            d=d0_search-1
 150            call score_fun       !init, get scores, n_cut+i_ali(i) for iteration
 151            if(score_max.lt.score)then
 152               score_max=score
 153               ka0=ka
 154               do i=1,ka0
 155                  k_ali0(i)=k_ali(i)
 156               enddo
 157            endif
 158 ***   iteration for extending ---------------------------------->
 159            d=d0_search+1
 160            do 301 it=1,n_it
 161               LL=0
 162               ka=0
 163               do i=1,n_cut
 164                  m=i_ali(i)     ![1,n_ali]
 165                  r_1(1,i)=xa(iA(m))
 166                  r_1(2,i)=ya(iA(m))
 167                  r_1(3,i)=za(iA(m))
 168                  r_2(1,i)=xb(iB(m))
 169                  r_2(2,i)=yb(iB(m))
 170                  r_2(3,i)=zb(iB(m))
 171                  ka=ka+1
 172                  k_ali(ka)=m
 173                  LL=LL+1
 174               enddo
 175               call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2
 176               do j=1,nseqA
 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)
 180               enddo
 181               call score_fun    !get scores, n_cut+i_ali(i) for iteration
 182               if(score_max.lt.score)then
 183                  score_max=score
 184                  ka0=ka
 185                  do i=1,ka
 186                     k_ali0(i)=k_ali(i)
 187                  enddo
 188               endif
 189               if(it.eq.n_it)goto 302
 190               if(n_cut.eq.ka)then
 191                  neq=0
 192                  do i=1,n_cut
 193                     if(i_ali(i).eq.k_ali(i))neq=neq+1
 194                  enddo
 195                  if(n_cut.eq.neq)goto 302
 196               endif
 197  301       continue             !for iteration
 198  302       continue
 199  300    continue                !for shift
 200  333  continue                  !for initial length, L_ali/M
 201
 202 ******** return the final rotation ****************
 203       LL=0
 204       do i=1,ka0
 205          m=k_ali0(i)            !record of the best alignment
 206          r_1(1,i)=xa(iA(m))
 207          r_1(2,i)=ya(iA(m))
 208          r_1(3,i)=za(iA(m))
 209          r_2(1,i)=xb(iB(m))
 210          r_2(2,i)=yb(iB(m))
 211          r_2(3,i)=zb(iB(m))
 212          LL=LL+1
 213       enddo
 214       call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2
 215       do j=1,nseqA
 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)
 219       enddo
 220       TM=score_max
 221
 222 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 223       return
 224       END
 225
 226 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 227 c     1, collect those residues with dis<d;
 228 c     2, calculate score_GDT, score_maxsub, score_TM
 229 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 230       subroutine score_fun
 231       include 'DIMENSIONS'
 232       PARAMETER(nmax=maxres)
 233
 234       common/stru/xa(nmax),ya(nmax),za(nmax),xb(nmax),yb(nmax),zb(nmax)
 235       common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
 236       common/para/d,d0
 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
 239       common/scores/score
 240       double precision score,score_max
 241
 242       d_tmp=d
 243  21   n_cut=0                   !number of residue-pairs dis<d, for iteration
 244       score_sum=0               !TMscore
 245       do k=1,n_ali
 246          i=iA(k)                ![1,nseqA] reoder number of structureA
 247          j=iB(k)                ![1,nseqB]
 248          dis=sqrt((xa(i)-xb(j))**2+(ya(i)-yb(j))**2+(za(i)-zb(j))**2)
 249          if(dis.lt.d_tmp)then
 250             n_cut=n_cut+1
 251             i_ali(n_cut)=k      ![1,n_ali], mark the residue-pairs in dis<d
 252          endif
 253          score_sum=score_sum+1/(1+(dis/d0)**2)
 254       enddo
 255       if(n_cut.lt.3.and.n_ali.gt.3)then
 256          d_tmp=d_tmp+.5
 257          goto 21
 258       endif
 259       score=score_sum/float(nseqB) !TM-score
 260
 261       return
 262       end
 263
 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,*)
 281       integer n, mode
 282
 283       double precision rms, u(3,3), t(3)
 284       integer ier
 285
 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
 292
 293       double precision sqrt3, tol, zero
 294
 295       data sqrt3 / 1.73205080756888d+00 /
 296       data tol / 1.0d-2 /
 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 /
 300
 301       wc = zero
 302       rms = zero
 303       e0 = zero
 304
 305       do i=1, 3
 306          xc(i) = zero
 307          yc(i) = zero
 308          t(i) = zero
 309          do j=1, 3
 310             r(i,j) = zero
 311             u(i,j) = zero
 312             a(i,j) = zero
 313             if( i .eq. j ) then
 314                u(i,j) = 1.0
 315                a(i,j) = 1.0
 316             end if
 317          end do
 318       end do
 319
 320       ier = -1
 321       if( n .lt. 1 ) return
 322       ier = -2
 323       do m=1, n
 324          if( w(m) .lt. 0.0 ) return
 325          wc = wc + w(m)
 326          do i=1, 3
 327             xc(i) = xc(i) + w(m)*x(i,m)
 328             yc(i) = yc(i) + w(m)*y(i,m)
 329          end do
 330       end do
 331       if( wc .le. zero ) return
 332       do i=1, 3
 333          xc(i) = xc(i) / wc
 334          yc(i) = yc(i) / wc
 335       end do
 336
 337       do m=1, n
 338          do i=1, 3
 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) )
 341             do j=1, 3
 342                r(i,j) = r(i,j) + d*( x(j,m) - xc(j) )
 343             end do
 344          end do
 345       end do
 346
 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)) )
 350
 351       sigma = det
 352
 353       m = 0;
 354       do j=1, 3
 355          do i=1, j
 356             m = m+1
 357             rr(m) = r(1,i)*r(1,j) + r(2,i)*r(2,j) + r(3,i)*r(3,j)
 358          end do
 359       end do
 360
 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
 364       det = det*det
 365
 366       do i=1, 3
 367          e(i) = spur
 368       end do
 369       if( spur .le. zero ) goto 40
 370       d = spur*spur
 371       h = d - cof
 372       g = (spur*cof - det)/2.0 - spur*h
 373       if( h .le. zero ) then
 374          if( mode .eq. 0 ) then
 375             goto 50
 376          else
 377             goto 30
 378          end if
 379       end if
 380       sqrth = dsqrt(h)
 381       d = h*h*h - g*g
 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
 389
 390       if( mode .eq. 0 ) then
 391          goto 50
 392       end if
 393
 394       do l=1, 3, 2
 395          d = e(l)
 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)
 402
 403          if( dabs(ss(1)) .ge. dabs(ss(3)) ) then
 404             j=1
 405             if( dabs(ss(1)) .lt. dabs(ss(6)) ) j = 3
 406          else if( dabs(ss(3)) .ge. dabs(ss(6)) ) then
 407             j = 2
 408          else
 409             j = 3
 410          end if
 411
 412          d = zero
 413          j = 3 * (j - 1)
 414
 415          do i=1, 3
 416             k = ip(i+j)
 417             a(i,l) = ss(k)
 418             d = d + ss(k)*ss(k)
 419          end do
 420          if( d .gt. zero ) d = 1.0 / dsqrt(d)
 421          do i=1, 3
 422             a(i,l) = a(i,l) * d
 423          end do
 424       end do
 425
 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
 428          m1 = 3
 429          m = 1
 430       else
 431          m1 = 1
 432          m = 3
 433       endif
 434
 435       p = zero
 436       do i=1, 3
 437          a(i,m1) = a(i,m1) - d*a(i,m)
 438          p = p + a(i,m1)**2
 439       end do
 440       if( p .le. tol ) then
 441          p = 1.0
 442          do i=1, 3
 443             if (p .lt. dabs(a(i,m))) cycle
 444             p = dabs( a(i,m) )
 445             j = i
 446          end do
 447          k = ip2312(j)
 448          l = ip2312(j+1)
 449          p = dsqrt( a(k,m)**2 + a(l,m)**2 )
 450          if( p .le. tol ) goto 40
 451          a(j,m1) = zero
 452          a(k,m1) = -a(l,m)/p
 453          a(l,m1) =  a(k,m)/p
 454       else
 455          p = 1.0 / dsqrt(p)
 456          do i=1, 3
 457             a(i,m1) = a(i,m1)*p
 458          end do
 459       end if
 460
 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)
 464
 465  30   do l=1, 2
 466          d = zero
 467          do i=1, 3
 468             b(i,l) = r(i,1)*a(1,l) + r(i,2)*a(2,l) + r(i,3)*a(3,l)
 469             d = d + b(i,l)**2
 470          end do
 471          if( d .gt. zero ) d = 1.0 / dsqrt(d)
 472          do i=1, 3
 473             b(i,l) = b(i,l)*d
 474          end do
 475       end do
 476       d = b(1,1)*b(1,2) + b(2,1)*b(2,2) + b(3,1)*b(3,2)
 477       p = zero
 478
 479       do i=1, 3
 480          b(i,2) = b(i,2) - d*b(i,1)
 481          p = p + b(i,2)**2
 482       end do
 483       if( p .le. tol ) then
 484          p = 1.0
 485          do i=1, 3
 486             if( p .lt. dabs(b(i,1)) ) cycle
 487             p = dabs( b(i,1) )
 488             j = i
 489          end do
 490          k = ip2312(j)
 491          l = ip2312(j+1)
 492          p = dsqrt( b(k,1)**2 + b(l,1)**2 )
 493          if( p .le. tol ) goto 40
 494          b(j,2) = zero
 495          b(k,2) = -b(l,1)/p
 496          b(l,2) =  b(k,1)/p
 497       else
 498          p = 1.0 / dsqrt(p)
 499          do i=1, 3
 500             b(i,2) = b(i,2)*p
 501          end do
 502       end if
 503
 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)
 507
 508       do i=1, 3
 509          do j=1, 3
 510             u(i,j) = b(i,1)*a(j,1) + b(i,2)*a(j,2) + b(i,3)*a(j,3)
 511          end do
 512       end do
 513
 514  40   do i=1, 3
 515          t(i) = ((yc(i) - u(i,1)*xc(1)) - u(i,2)*xc(2)) - u(i,3)*xc(3)
 516       end do
 517  50   do i=1, 3
 518          if( e(i) .lt. zero ) e(i) = zero
 519          e(i) = dsqrt( e(i) )
 520       end do
 521
 522       ier = 0
 523       if( e(2) .le. (e(1) * 1.0d-05) ) ier = -1
 524
 525       d = e(3)
 526       if( sigma .lt. 0.0 ) then
 527          d = - d
 528          if( (e(2) - e(3)) .le. (e(1) * 1.0d-05) ) ier = -1
 529       end if
 530       d = (d + e(2)) + e(1)
 531
 532       rms = (e0 - d) - d
 533       if( rms .lt. 0.0 ) rms = 0.0
 534
 535       return
 536       end