source/cluster/wham/src-M-homology/TMscore.F

   1 *************************************************************************
   2 *     This program is to compare two protein structures and identify the
   3 *     best superposition that has the highest TM-score. Input structures
   4 *     must be in the PDB format. By default, TM-score is normalized by
   5 *     the second protein. Users can obtain a brief instruction by simply
   6 *     running the program without arguments. For comments/suggestions,
   7 *     please contact email: zhng@umich.edu.
   8 *
   9 *     Reference:
  10 *     Yang Zhang, Jeffrey Skolnick, Proteins, 2004 57:702-10.
  11 *
  12 *     Permission to use, copy, modify, and distribute this program for
  13 *     any purpose, with or without fee, is hereby granted, provided that
  14 *     the notices on the head, the reference information, and this
  15 *     copyright notice appear in all copies or substantial portions of
  16 *     the Software. It is provided "as is" without express or implied
  17 *     warranty.
  18 ******************* Updating history ************************************
  19 *     2005/10/19: the program was reformed so that the score values.
  20 *                 are not dependent on the specific compilers.
  21 *     2006/06/20: selected 'A' if there is altLoc when reading PDB file.
  22 *     2007/02/05: fixed a bug with length<15 in TMscore_32.
  23 *     2007/02/27: rotation matrix from Chain-1 to Chain-2 was added.
  24 *     2007/12/06: GDT-HA score was added, fixed a bug for reading PDB.
  25 *     2010/08/02: A new RMSD matrix was used and obsolete statement removed.
  26 *     2011/01/03: The length of pdb file names were extended to 500.
  27 *     2011/01/30: An open source license is attached to the program.
  28 *     2012/05/07: Improved RMSD calculation subroutine which speeds up
  29 *                 TM-score program by 30%.
  30 *     2012/06/05: Added option '-l L' which calculates TM-score (and maxsub
  31 *                 and GDT scores) normalized by a specific length 'L'.
  32 *************************************************************************
  33
  34 c      program TMscore
  35       subroutine TMscore_sub(rmsd,gdt_ts,gdt_ha,tmscore,cfname,lprint)
  36       include 'DIMENSIONS'
  37       PARAMETER(nmax=5000)
  38       include 'COMMON.CONTROL'
  39       include 'COMMON.IOUNITS'
  40       include 'COMMON.CHAIN'
  41       include 'COMMON.INTERACT'
  42
  43       real*8 rmsd,gdt_ts,gdt_ha,tmscore
  44       common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax)
  45       common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
  46       common/para/d,d0,d0_fix
  47       common/align/n_ali,iA(nmax),iB(nmax)
  48       common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score
  49       dimension k_ali(nmax),k_ali0(nmax)
  50
  51       character*500 fnam,pdb(100)!,outname
  52       character*80 cfname
  53       character*3 aa(-1:20),seqA(nmax),seqB(nmax)
  54       character*500 s,du
  55       character seq1A(nmax),seq1B(nmax),ali(nmax)
  56       character sequenceA(nmax),sequenceB(nmax),sequenceM(nmax)
  57
  58       dimension L_ini(100),iq(nmax)
  59       common/scores/score,score_maxsub,score_fix,score10
  60       common/GDT/n_GDT05,n_GDT1,n_GDT2,n_GDT4,n_GDT8
  61       double precision score,score_max,score_fix,score_fix_max
  62       double precision score_maxsub,score10
  63       dimension xa(nmax),ya(nmax),za(nmax)
  64
  65 ccc   RMSD:
  66       double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)
  67       double precision u(3,3),tt(3),rms,drms !armsd is real
  68       data w /nmax*1.0/
  69
  70       logical lprint
  71 ccc
  72
  73       data aa/ 'BCK','GLY','ALA','SER','CYS',
  74      &     'VAL','THR','ILE','PRO','MET',
  75      &     'ASP','ASN','LEU','LYS','GLU',
  76      &     'GLN','ARG','HIS','PHE','TYR',
  77      &     'TRP','CYX'/
  78       character*1 slc(-1:20)
  79       data slc/'X','G','A','S','C',
  80      &     'V','T','I','P','M',
  81      &     'D','N','L','K','E',
  82      &     'Q','R','H','F','Y',
  83      &     'W','C'/
  84
  85 *****instructions ----------------->
  86 c      call getarg(1,fnam)
  87 c      if(fnam.eq.' '.or.fnam.eq.'?'.or.fnam.eq.'-h')then
  88 c         write(*,*)
  89 c         write(*,*)'Brief instruction for running TM-score program:'
  90 c         write(*,*)'(For detail: Zhang & Skolnick,  Proteins, 2004',
  91 c     &        ' 57:702-10)'
  92 c         write(*,*)
  93 c         write(*,*)'1. Run TM-score to compare ''model'' and ',
  94 c     &        '''native'':'
  95 c         write(*,*)'   >TMscore model native'
  96 c         write(*,*)
  97 c         write(*,*)'2. TM-score normalized with an assigned scale d0',
  98 c     &        ' e.g. 5 A:'
  99 c         write(*,*)'   >TMscore model native -d 5'
 100 c         write(*,*)
 101 c         write(*,*)'3. TM-score normalized by a specific length, ',
 102 c     &        'e.g. 120 AA:'
 103 c         write(*,*)'   >TMscore model native -l 120'
 104 c         write(*,*)
 105 c         write(*,*)'4. TM-score with superposition output, e.g. ',
 106 c     &        '''TM.sup'':'
 107 c         write(*,*)'   >TMscore model native -o TM.sup'
 108 c         write(*,*)'   To view the superimposed structures by rasmol:'
 109 c         write(*,*)'      >rasmol -script TM.sup'
 110 c         write(*,*)
 111 c         goto 9999
 112 c      endif
 113
 114       pdb(1)=cfname
 115       pdb(2)=pdbfile
 116 ******* options ----------->
 117       m_out=-1
 118       m_fix=-1
 119       m_len=-1
 120 c      narg=iargc()
 121 c      i=0
 122 c      j=0
 123 c 115  continue
 124 c      i=i+1
 125 c      call getarg(i,fnam)
 126 c      if(fnam.eq.'-o')then
 127 c         m_out=1
 128 c         i=i+1
 129 c         call getarg(i,outname)
 130 c      elseif(fnam.eq.'-d')then
 131 c         m_fix=1
 132 c         i=i+1
 133 c         call getarg(i,fnam)
 134 c         read(fnam,*)d0_fix
 135 c      elseif(fnam.eq.'-l')then
 136 c         m_len=1
 137 c         i=i+1
 138 c         call getarg(i,fnam)
 139 c         read(fnam,*)l0_fix
 140 c      else
 141 c         j=j+1
 142 c         pdb(j)=fnam
 143 c      endif
 144 c      if(i.lt.narg)goto 115
 145 c
 146 ccccccccc read data from first CA file:
 147 c      open(unit=10,file=pdb(1),status='old')
 148 c      i=0
 149 c 101  read(10,104,end=102) s
 150 c      if(s(1:3).eq.'TER') goto 102
 151 c      if(s(1:4).eq.'ATOM')then
 152 c         if(s(13:16).eq.'CA  '.or.s(13:16).eq.' CA '.or.s(13:16).
 153 c     &        eq.'  CA')then
 154 c         if(s(17:17).eq.' '.or.s(17:17).eq.'A')then
 155 c            i=i+1
 156 c            read(s,103)du,seqA(i),du,nresA(i),du,xa(i),ya(i),za(i)
 157 c            do j=-1,20
 158 c               if(seqA(i).eq.aa(j))then
 159 c                  seq1A(i)=slc(j)
 160 c                  goto 21
 161 c               endif
 162 c            enddo
 163 c            seq1A(i)=slc(-1)
 164 c 21         continue
 165 c         endif
 166 c         endif
 167 c      endif
 168 c      goto 101
 169 c 102  continue
 170 c 103  format(A17,A3,A2,i4,A4,3F8.3)
 171 c 104  format(A100)
 172 c      close(10)
 173 c      nseqA=i
 174 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 175 c
 176 ccccccccc read data from first CA file:
 177 c      open(unit=10,file=pdb(2),status='old')
 178 c      i=0
 179 c 201  read(10,204,end=202) s
 180 c      if(s(1:3).eq.'TER') goto 202
 181 c      if(s(1:4).eq.'ATOM')then
 182 c         if(s(13:16).eq.'CA  '.or.s(13:16).eq.' CA '.or.s(13:16).
 183 c     &        eq.'  CA')then
 184 c         if(s(17:17).eq.' '.or.s(17:17).eq.'A')then
 185 c            i=i+1
 186 c            read(s,203)du,seqB(i),du,nresB(i),du,xb(i),yb(i),zb(i)
 187 c            do j=-1,20
 188 c               if(seqB(i).eq.aa(j))then
 189 c                  seq1B(i)=slc(j)
 190 c                  goto 22
 191 c               endif
 192 c            enddo
 193 c            seq1B(i)=slc(-1)
 194 c 22         continue
 195 c         endif
 196 c         endif
 197 c      endif
 198 c      goto 201
 199 c 202  continue
 200 c 203  format(A17,A3,A2,i4,A4,3F8.3)
 201 c 204  format(A100)
 202 c      close(10)
 203 c      nseqB=i
 204 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 205
 206 ******************************************************************
 207 *     pickup the aligned residues:
 208 ******************************************************************
 209 c      k=0
 210 c      do i=1,nseqA
 211 c         do j=1,nseqB
 212 c            if(nresA(i).eq.nresB(j))then
 213 c               k=k+1
 214 c               iA(k)=i
 215 c               iB(k)=j
 216 c               goto 205
 217 c            endif
 218 c         enddo
 219 c 205     continue
 220 c      enddo
 221 c      n_ali=k                   !number of aligned residues
 222 c      if(n_ali.lt.1)then
 223 c        write(*,*)'There is no common residues in the input structures'
 224 c        goto 9999
 225 c      endif
 226 c
 227 ************/////
 228 *     parameters:
 229 *****************
 230       noverlap=nres
 231       if (nres.gt.nsup+nnt-1) noverlap=nsup+nnt-1
 232       nnsup=0
 233       do i=1,noverlap
 234         if (itype(i).ne.ntyp1) then
 235           nnsup=nnsup+1
 236           iA(nnsup)=nnsup
 237           iB(nnsup)=nnsup
 238         endif
 239       enddo
 240       nseqA=nnsup
 241       nseqB=nnsup
 242       n_ali=nnsup
 243 ***   d0------------->
 244       if(nseqB.gt.15)then
 245          d0=1.24*(nseqB-15)**(1.0/3.0)-1.8
 246       else
 247          d0=0.5
 248       endif
 249       if(m_len.eq.1)then
 250          d0=1.24*(l0_fix-15)**(1.0/3.0)-1.8
 251       endif
 252       if(d0.lt.0.5)d0=0.5
 253       if(m_fix.eq.1)d0=d0_fix
 254 ***   d0_search ----->
 255       d0_search=d0
 256       if(d0_search.gt.8)d0_search=8
 257       if(d0_search.lt.4.5)d0_search=4.5
 258 ***   iterative parameters ----->
 259       n_it=20                   !maximum number of iterations
 260       d_output=5                !for output alignment
 261       if(m_fix.eq.1)d_output=d0_fix
 262       n_init_max=6              !maximum number of L_init
 263       n_init=0
 264       L_ini_min=4
 265       if(n_ali.lt.4)L_ini_min=n_ali
 266       do i=1,n_init_max-1
 267          n_init=n_init+1
 268          L_ini(n_init)=n_ali/2**(n_init-1)
 269          if(L_ini(n_init).le.L_ini_min)then
 270             L_ini(n_init)=L_ini_min
 271             goto 402
 272          endif
 273       enddo
 274       n_init=n_init+1
 275       L_ini(n_init)=L_ini_min
 276  402  continue
 277
 278 ******************************************************************
 279 *     find the maximum score starting from local structures superposition
 280 ******************************************************************
 281       score_max=-1              !TM-score
 282       score_maxsub_max=-1       !MaxSub-score
 283       score10_max=-1            !TM-score10
 284       n_GDT05_max=-1            !number of residues<0.5
 285       n_GDT1_max=-1             !number of residues<1
 286       n_GDT2_max=-1             !number of residues<2
 287       n_GDT4_max=-1             !number of residues<4
 288       n_GDT8_max=-1             !number of residues<8
 289
 290 #ifdef DEBUG
 291       write (iout,*) "cref and ccref"
 292 #endif
 293       noverlap=nres
 294       if (nres.gt.nsup+nnt-1) noverlap=nsup+nnt-1
 295       nnsup=0
 296       do i=1,noverlap
 297         if (itype(i).ne.ntyp1) then
 298           nnsup=nnsup+1
 299           xa(nnsup)=c(1,i)
 300           ya(nnsup)=c(2,i)
 301           za(nnsup)=c(3,i)
 302           xb(nnsup)=cref_pdb(1,i,1)
 303           yb(nnsup)=cref_pdb(2,i,1)
 304           zb(nnsup)=cref_pdb(3,i,1)
 305 c          do j=1,3
 306 c            cc(j,nnsup)=c(j,i)
 307 c            ccref(j,nnsup)=cref_pdb(j,i,1)
 308 c          enddo
 309 #ifdef DEBUG
 310           write (iout,'(i5,3f10.5,5x,3f10.5)') nnsup,
 311      &     xa(nnsup),ya(nnsup),za(nnsup),xb(nnsup),yb(nnsup),zb(nnsup)
 312 #endif
 313         endif
 314       enddo
 315
 316       do 333 i_init=1,n_init
 317         L_init=L_ini(i_init)
 318         iL_max=n_ali-L_init+1
 319         do 300 iL=1,iL_max      !on aligned residues, [1,nseqA]
 320            LL=0
 321            ka=0
 322            do i=1,L_init
 323               k=iL+i-1          ![1,n_ali] common aligned
 324               r_1(1,i)=xa(iA(k))
 325               r_1(2,i)=ya(iA(k))
 326               r_1(3,i)=za(iA(k))
 327               r_2(1,i)=xb(iB(k))
 328               r_2(2,i)=yb(iB(k))
 329               r_2(3,i)=zb(iB(k))
 330               ka=ka+1
 331               k_ali(ka)=k
 332               LL=LL+1
 333            enddo
 334            if(i_init.eq.1)then  !global superposition
 335               call u3b(w,r_1,r_2,LL,2,rms,u,tt,ier) !0:rmsd; 1:u,t; 2:rmsd,u,t
 336               armsd=dsqrt(rms/LL)
 337               rmsd_ali=armsd
 338            else
 339               call u3b(w,r_1,r_2,LL,1,rms,u,tt,ier) !u rotate r_1 to r_2
 340            endif
 341            do j=1,nseqA
 342               xt(j)=tt(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
 343               yt(j)=tt(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
 344               zt(j)=tt(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
 345            enddo
 346            d=d0_search-1
 347            call score_fun       !init, get scores, n_cut+i_ali(i) for iteration
 348            if(score_max.lt.score)then
 349               score_max=score
 350               ka0=ka
 351               do i=1,ka0
 352                  k_ali0(i)=k_ali(i)
 353               enddo
 354            endif
 355            if(score10_max.lt.score10)score10_max=score10
 356            if(score_maxsub_max.lt.score_maxsub)score_maxsub_max=
 357      &          score_maxsub
 358            if(n_GDT05_max.lt.n_GDT05)n_GDT05_max=n_GDT05
 359            if(n_GDT1_max.lt.n_GDT1)n_GDT1_max=n_GDT1
 360            if(n_GDT2_max.lt.n_GDT2)n_GDT2_max=n_GDT2
 361            if(n_GDT4_max.lt.n_GDT4)n_GDT4_max=n_GDT4
 362            if(n_GDT8_max.lt.n_GDT8)n_GDT8_max=n_GDT8
 363 ***   iteration for extending ---------------------------------->
 364            d=d0_search+1
 365            do 301 it=1,n_it
 366               LL=0
 367               ka=0
 368               do i=1,n_cut
 369                  m=i_ali(i)     ![1,n_ali]
 370                  r_1(1,i)=xa(iA(m))
 371                  r_1(2,i)=ya(iA(m))
 372                  r_1(3,i)=za(iA(m))
 373                  r_2(1,i)=xb(iB(m))
 374                  r_2(2,i)=yb(iB(m))
 375                  r_2(3,i)=zb(iB(m))
 376                  ka=ka+1
 377                  k_ali(ka)=m
 378                  LL=LL+1
 379               enddo
 380               call u3b(w,r_1,r_2,LL,1,rms,u,tt,ier) !u rotate r_1 to r_2
 381               do j=1,nseqA
 382                  xt(j)=tt(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
 383                  yt(j)=tt(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
 384                  zt(j)=tt(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
 385               enddo
 386               call score_fun    !get scores, n_cut+i_ali(i) for iteration
 387               if(score_max.lt.score)then
 388                  score_max=score
 389                  ka0=ka
 390                  do i=1,ka
 391                     k_ali0(i)=k_ali(i)
 392                  enddo
 393               endif
 394               if(score10_max.lt.score10)score10_max=score10
 395               if(score_maxsub_max.lt.score_maxsub)score_maxsub_max
 396      &             =score_maxsub
 397               if(n_GDT05_max.lt.n_GDT05)n_GDT05_max=n_GDT05
 398               if(n_GDT1_max.lt.n_GDT1)n_GDT1_max=n_GDT1
 399               if(n_GDT2_max.lt.n_GDT2)n_GDT2_max=n_GDT2
 400               if(n_GDT4_max.lt.n_GDT4)n_GDT4_max=n_GDT4
 401               if(n_GDT8_max.lt.n_GDT8)n_GDT8_max=n_GDT8
 402               if(it.eq.n_it)goto 302
 403               if(n_cut.eq.ka)then
 404                  neq=0
 405                  do i=1,n_cut
 406                     if(i_ali(i).eq.k_ali(i))neq=neq+1
 407                  enddo
 408                  if(n_cut.eq.neq)goto 302
 409               endif
 410  301       continue             !for iteration
 411  302       continue
 412  300    continue                !for shift
 413  333  continue                  !for initial length, L_ali/M
 414
 415       ratio=1
 416       if(m_len.gt.0)then
 417          ratio=float(nseqB)/float(l0_fix)
 418       endif
 419       if(m_len.eq.1)then
 420          score_max=score_max*float(nseqB)/float(l0_fix)
 421       endif
 422       score_GDT=(n_GDT1_max+n_GDT2_max+n_GDT4_max+n_GDT8_max)
 423      &     /float(4*nseqB)
 424       score_GDT_HA=(n_GDT05_max+n_GDT1_max+n_GDT2_max+n_GDT4_max)
 425      &     /float(4*nseqB)
 426       tmscore=score_max
 427       gdt_ts=score_GDT*ratio
 428       gdt_ha=score_GDT_HA*ratio
 429       rmsd=rmsd_ali
 430 ******************************************************************
 431 *     Output
 432 ******************************************************************
 433 ***   output TM-scale ---------------------------->
 434
 435       if (lprint) then
 436
 437       write(iout,*)
 438       write(iout,*)'**************************************************',
 439      &     '***************************'
 440       write(iout,*)'*                              TM-SCORE           ',
 441      &     '                          *'
 442       write(iout,*)'* A scoring function to assess the similarity of p',
 443      &     'rotein structures         *'
 444       write(iout,*)'* Based on statistics:                            ',
 445      &     '                          *'
 446       write(iout,*)'*       0.0 < TM-score < 0.17, random structural s',
 447      &     'imilarity                 *'
 448       write(iout,*)'*       0.5 < TM-score < 1.00, in about the same f',
 449      &     'old                    *'
 450       write(iout,*)'* Reference: Yang Zhang and Jeffrey Skolnick, ',
 451      &     'Proteins 2004 57: 702-710     *'
 452       write(iout,*)'* For comments, please email to: zhng@umich.edu   ',
 453      &     '                          *'
 454       write(iout,*)'**************************************************',
 455      &     '***************************'
 456       write(iout,*)
 457       write(iout,501)pdb(1),nseqA
 458  501  format('Structure1: ',A10,'  Length= ',I4)
 459       if(m_len.eq.1)then
 460          write(iout,411)pdb(2),nseqB
 461          write(iout,412)l0_fix
 462       else
 463          write(iout,502)pdb(2),nseqB
 464       endif
 465  411  format('Structure2: ',A10,'  Length= ',I4)
 466  412  format('TM-score is notmalized by ',I4)
 467  502  format('Structure2: ',A10,'  Length= ',I4,
 468      &     ' (by which all scores are normalized)')
 469       write(iout,503)n_ali
 470  503  format('Number of residues in common= ',I4)
 471       write(iout,513)rmsd_ali
 472  513  format('RMSD of  the common residues= ',F8.3)
 473       write(iout,*)
 474       write(iout,504)score_max,d0
 475  504  format('TM-score    = ',f6.4,'  (d0=',f5.2,')')
 476       write(iout,505)score_maxsub_max*ratio
 477  505  format('MaxSub-score= ',f6.4,'  (d0= 3.50)')
 478       write(iout,506)score_GDT*ratio,n_GDT1_max/float(nseqB)*ratio,
 479      &     n_GDT2_max/float(nseqB)*ratio,n_GDT4_max/float(nseqB)*ratio,
 480      &     n_GDT8_max/float(nseqB)*ratio
 481  506  format('GDT-TS-score= ',f6.4,' %(d<1)=',f6.4,' %(d<2)=',f6.4,
 482      $     ' %(d<4)=',f6.4,' %(d<8)=',f6.4)
 483       write(iout,507)score_GDT_HA*ratio,n_GDT05_max/float(nseqB)*ratio,
 484      &     n_GDT1_max/float(nseqB)*ratio,n_GDT2_max/float(nseqB)*ratio,
 485      &     n_GDT4_max/float(nseqB)*ratio
 486  507  format('GDT-HA-score= ',f6.4,' %(d<0.5)=',f6.4,' %(d<1)=',f6.4,
 487      $     ' %(d<2)=',f6.4,' %(d<4)=',f6.4)
 488       write(iout,*)
 489
 490       endif
 491
 492       return
 493       end
 494 c------------------------------------------------------------------------
 495 ***   recall and output the superposition of maxiumum TM-score:
 496 c      LL=0
 497 c      do i=1,ka0
 498 c         m=k_ali0(i)            !record of the best alignment
 499 c         r_1(1,i)=xa(iA(m))
 500 c         r_1(2,i)=ya(iA(m))
 501 c         r_1(3,i)=za(iA(m))
 502 c         r_2(1,i)=xb(iB(m))
 503 c         r_2(2,i)=yb(iB(m))
 504 c         r_2(3,i)=zb(iB(m))
 505 c         LL=LL+1
 506 c      enddo
 507 c      call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2
 508 c      do j=1,nseqA
 509 c         xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
 510 c         yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
 511 c         zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
 512 c      enddo
 513 c
 514 c********* extract rotation matrix ------------>
 515 c      write(*,*)'-------- rotation matrix to rotate Chain-1 to ',
 516 c     &     'Chain-2 ------'
 517 c      write(*,*)'i          t(i)         u(i,1)         u(i,2) ',
 518 c     &     '        u(i,3)'
 519 c      do i=1,3
 520 c         write(*,304)i,t(i),u(i,1),u(i,2),u(i,3)
 521 c      enddo
 522 cc      do j=1,nseqA
 523 cc         xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
 524 cc         yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
 525 cc         zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
 526 cc         write(*,*)j,xt(j),yt(j),zt(j)
 527 cc      enddo
 528 c      write(*,*)
 529 c 304  format(I2,f18.10,f15.10,f15.10,f15.10)
 530 c
 531 c********* rmsd in superposed regions --------------->
 532 c      d=d_output                !for output
 533 c      call score_fun()          !give i_ali(i), score_max=score now
 534 c      LL=0
 535 c      do i=1,n_cut
 536 c         m=i_ali(i)             ![1,nseqA]
 537 c         r_1(1,i)=xa(iA(m))
 538 c         r_1(2,i)=ya(iA(m))
 539 c         r_1(3,i)=za(iA(m))
 540 c         r_2(1,i)=xb(iB(m))
 541 c         r_2(2,i)=yb(iB(m))
 542 c         r_2(3,i)=zb(iB(m))
 543 c         LL=LL+1
 544 c      enddo
 545 c      call u3b(w,r_1,r_2,LL,0,rms,u,t,ier)
 546 c      armsd=dsqrt(rms/LL)
 547 c      rmsd=armsd
 548 c
 549 c***   output rotated chain1 + chain2----->
 550 c      if(m_out.ne.1)goto 999
 551 c      OPEN(unit=7,file=outname,status='unknown') !pdb1.aln + pdb2.aln
 552 c 900  format(A)
 553 c 901  format('select ',I4)
 554 c      write(7,900)'load inline'
 555 c      write(7,900)'select atomno<1000'
 556 cc      write(7,900)'color [255,20,147]'
 557 c      write(7,900)'wireframe .45'
 558 c      write(7,900)'select none'
 559 c      write(7,900)'select atomno>1000'
 560 cc      write(7,900)'color [100,149,237]'
 561 c      write(7,900)'wireframe .15'
 562 c      write(7,900)'color white'
 563 c      do i=1,n_cut
 564 c         write(7,901)nresA(iA(i_ali(i)))
 565 c         write(7,900)'color red'
 566 c      enddo
 567 c      write(7,900)'select all'
 568 c      write(7,900)'exit'
 569 c      write(7,514)rmsd_ali
 570 c 514  format('REMARK  RMSD of the common residues=',F8.3)
 571 c      write(7,515)score_max,d0
 572 c 515  format('REMARK  TM-score=',f6.4,' (d0=',f5.2,')')
 573 c      do i=1,nseqA
 574 c         write(7,1237)nresA(i),seqA(i),nresA(i),xt(i),yt(i),zt(i)
 575 c      enddo
 576 c      write(7,1238)
 577 c      do i=2,nseqA
 578 c         write(7,1239)nresA(i-1),nresA(i)
 579 c      enddo
 580 c      do i=1,nseqB
 581 c         write(7,1237)2000+nresB(i),seqB(i),nresB(i),xb(i),yb(i),zb(i)
 582 c      enddo
 583 c      write(7,1238)
 584 c      do i=2,nseqB
 585 c         write(7,1239)2000+nresB(i-1),2000+nresB(i)
 586 c      enddo
 587 c 1237 format('ATOM  ',i5,'  CA  ',A3,I6,4X,3F8.3)
 588 c 1238 format('TER')
 589 c 1239 format('CONECT',I5,I5)
 590 c 999  continue
 591 c
 592 c***   record aligned residues by i=[1,nseqA], for sequenceM()------------>
 593 c      do i=1,nseqA
 594 c         iq(i)=0
 595 c      enddo
 596 c      do i=1,n_cut
 597 c         j=iA(i_ali(i))         ![1,nseqA]
 598 c         k=iB(i_ali(i))         ![1,nseqB]
 599 c         dis=sqrt((xt(j)-xb(k))**2+(yt(j)-yb(k))**2+(zt(j)-zb(k))**2)
 600 c         if(dis.lt.d_output)then
 601 c            iq(j)=1
 602 c         endif
 603 c      enddo
 604 c*******************************************************************
 605 c***   output aligned sequences
 606 c      k=0
 607 c      i=1
 608 c      j=1
 609 c 800  continue
 610 c      if(i.gt.nseqA.and.j.gt.nseqB)goto 802
 611 c      if(i.gt.nseqA.and.j.le.nseqB)then
 612 c         k=k+1
 613 c         sequenceA(k)='-'
 614 c         sequenceB(k)=seq1B(j)
 615 c         sequenceM(k)=' '
 616 c         j=j+1
 617 c         goto 800
 618 c      endif
 619 c      if(i.le.nseqA.and.j.gt.nseqB)then
 620 c         k=k+1
 621 c         sequenceA(k)=seq1A(i)
 622 c         sequenceB(k)='-'
 623 c         sequenceM(k)=' '
 624 c         i=i+1
 625 c         goto 800
 626 c      endif
 627 c      if(nresA(i).eq.nresB(j))then
 628 c         k=k+1
 629 c         sequenceA(k)=seq1A(i)
 630 c         sequenceB(k)=seq1B(j)
 631 c         if(iq(i).eq.1)then
 632 c            sequenceM(k)=':'
 633 c         else
 634 c            sequenceM(k)=' '
 635 c         endif
 636 c         i=i+1
 637 c         j=j+1
 638 c         goto 800
 639 c      elseif(nresA(i).lt.nresB(j))then
 640 c         k=k+1
 641 c         sequenceA(k)=seq1A(i)
 642 c         sequenceB(k)='-'
 643 c         sequenceM(k)=' '
 644 c         i=i+1
 645 c         goto 800
 646 c      elseif(nresB(j).lt.nresA(i))then
 647 c         k=k+1
 648 c         sequenceA(k)='-'
 649 c         sequenceB(k)=seq1B(j)
 650 c         sequenceM(k)=' '
 651 c         j=j+1
 652 c         goto 800
 653 c      endif
 654 c 802  continue
 655 c
 656 c      write(*,600)d_output,n_cut,rmsd
 657 c 600  format('Superposition in the TM-score: Length(d<',f3.1,
 658 c     $     ')=',i3,'  RMSD=',f6.2)
 659 c      write(*,603)d_output
 660 c 603  format('(":" denotes the residue pairs of distance < ',f3.1,
 661 c     &     ' Angstrom)')
 662 c      write(*,601)(sequenceA(i),i=1,k)
 663 c      write(*,601)(sequenceM(i),i=1,k)
 664 c      write(*,601)(sequenceB(i),i=1,k)
 665 c      write(*,602)(mod(i,10),i=1,k)
 666 c 601  format(2000A1)
 667 c 602  format(2000I1)
 668 c      write(*,*)
 669 c
 670 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 671 c 9999 END
 672
 673 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 674 c     1, collect those residues with dis<d;
 675 c     2, calculate score_GDT, score_maxsub, score_TM
 676 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 677       subroutine score_fun
 678       PARAMETER(nmax=5000)
 679
 680       common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax)
 681       common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
 682       common/para/d,d0,d0_fix
 683       common/align/n_ali,iA(nmax),iB(nmax)
 684       common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score
 685       common/scores/score,score_maxsub,score_fix,score10
 686       common/GDT/n_GDT05,n_GDT1,n_GDT2,n_GDT4,n_GDT8
 687       double precision score,score_max,score_fix,score_fix_max
 688       double precision score_maxsub,score10
 689
 690       d_tmp=d
 691  21   n_cut=0                   !number of residue-pairs dis<d, for iteration
 692       n_GDT05=0                 !for GDT-score, # of dis<0.5
 693       n_GDT1=0                  !for GDT-score, # of dis<1
 694       n_GDT2=0                  !for GDT-score, # of dis<2
 695       n_GDT4=0                  !for GDT-score, # of dis<4
 696       n_GDT8=0                  !for GDT-score, # of dis<8
 697       score_maxsub_sum=0        !Maxsub-score
 698       score_sum=0               !TMscore
 699       score_sum10=0             !TMscore10
 700       do k=1,n_ali
 701          i=iA(k)                ![1,nseqA] reoder number of structureA
 702          j=iB(k)                ![1,nseqB]
 703          dis=sqrt((xt(i)-xb(j))**2+(yt(i)-yb(j))**2+(zt(i)-zb(j))**2)
 704 ***   for iteration:
 705          if(dis.lt.d_tmp)then
 706             n_cut=n_cut+1
 707             i_ali(n_cut)=k      ![1,n_ali], mark the residue-pairs in dis<d
 708          endif
 709 ***   for GDT-score:
 710          if(dis.le.8)then
 711             n_GDT8=n_GDT8+1
 712             if(dis.le.4)then
 713                n_GDT4=n_GDT4+1
 714                if(dis.le.2)then
 715                   n_GDT2=n_GDT2+1
 716                   if(dis.le.1)then
 717                      n_GDT1=n_GDT1+1
 718                      if(dis.le.0.5)then
 719                         n_GDT05=n_GDT05+1
 720                      endif
 721                   endif
 722                endif
 723             endif
 724          endif
 725 ***   for MAXsub-score:
 726          if(dis.lt.3.5)then
 727             score_maxsub_sum=score_maxsub_sum+1/(1+(dis/3.5)**2)
 728          endif
 729 ***   for TM-score:
 730          score_sum=score_sum+1/(1+(dis/d0)**2)
 731 ***   for TM-score10:
 732          if(dis.lt.10)then
 733             score_sum10=score_sum10+1/(1+(dis/d0)**2)
 734          endif
 735       enddo
 736       if(n_cut.lt.3.and.n_ali.gt.3)then
 737          d_tmp=d_tmp+.5
 738          goto 21
 739       endif
 740       score_maxsub=score_maxsub_sum/float(nseqB) !MAXsub-score
 741       score=score_sum/float(nseqB) !TM-score
 742       score10=score_sum10/float(nseqB) !TM-score10
 743
 744       return
 745       end
 746
 747 cccccccccccccccc Calculate sum of (r_d-r_m)^2 cccccccccccccccccccccccccc
 748 c  w    - w(m) is weight for atom pair  c m                    (given)
 749 c  x    - x(i,m) are coordinates of atom c m in set x          (given)
 750 c  y    - y(i,m) are coordinates of atom c m in set y          (given)
 751 c  n    - n is number of atom pairs                            (given)
 752 c  mode  - 0:calculate rms     only                            (given,short)
 753 c          1:calculate     u,t only                            (given,medium)
 754 c          2:calculate rms,u,t                                 (given,longer)
 755 c  rms   - sum of w*(ux+t-y)**2 over all atom pairs            (result)
 756 c  u    - u(i,j) is   rotation  matrix for best superposition  (result)
 757 c  t    - t(i)   is translation vector for best superposition  (result)
 758 c  ier  - 0: a unique optimal superposition has been determined(result)
 759 c       -1: superposition is not unique but optimal
 760 c       -2: no result obtained because of negative weights w
 761 c           or all weights equal to zero.
 762 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 763       subroutine u3b(w, x, y, n, mode, rms, u, t, ier)
 764       double precision w(*), x(3,*), y(3,*)
 765       integer n, mode
 766
 767       double precision rms, u(3,3), t(3)
 768       integer ier
 769
 770       integer i, j, k, l, m1, m
 771       integer ip(9), ip2312(4)
 772       double precision r(3,3), xc(3), yc(3), wc
 773       double precision a(3,3), b(3,3), e(3), rr(6), ss(6)
 774       double precision e0, d, spur, det, cof, h, g
 775       double precision cth, sth, sqrth, p, sigma
 776       double precision c1x, c1y, c1z, c2x, c2y, c2z
 777       double precision s1x, s1y, s1z, s2x, s2y, s2z
 778       double precision sxx, sxy, sxz, syx, syy, syz, szx, szy, szz
 779
 780       double precision sqrt3, tol, zero
 781
 782       data sqrt3 / 1.73205080756888d+00 /
 783       data tol / 1.0d-2 /
 784       data zero / 0.0d+00 /
 785       data ip / 1, 2, 4, 2, 3, 5, 4, 5, 6 /
 786       data ip2312 / 2, 3, 1, 2 /
 787
 788       wc  = zero
 789       rms = zero
 790       e0  = zero
 791       s1x = zero
 792       s1y = zero
 793       s1z = zero
 794       s2x = zero
 795       s2y = zero
 796       s2z = zero
 797       sxx = zero
 798       sxy = zero
 799       sxz = zero
 800       syx = zero
 801       syy = zero
 802       syz = zero
 803       szx = zero
 804       szy = zero
 805       szz = zero
 806
 807       do i=1, 3
 808          xc(i) = zero
 809          yc(i) = zero
 810          t(i) = zero
 811          do j=1, 3
 812             r(i,j) = zero
 813             u(i,j) = zero
 814             a(i,j) = zero
 815             if( i .eq. j ) then
 816                u(i,j) = 1.0
 817                a(i,j) = 1.0
 818             end if
 819          end do
 820       end do
 821
 822       ier = -1
 823       if( n .lt. 1 ) return
 824       ier = -2
 825
 826       do m=1, n
 827          c1x=x(1, m)
 828          c1y=x(2, m)
 829          c1z=x(3, m)
 830
 831          c2x=y(1, m)
 832          c2y=y(2, m)
 833          c2z=y(3, m)
 834
 835          s1x = s1x + c1x
 836          s1y = s1y + c1y;
 837          s1z = s1z + c1z;
 838
 839          s2x = s2x + c2x;
 840          s2y = s2y + c2y;
 841          s2z = s2z + c2z;
 842
 843          sxx = sxx + c1x*c2x;
 844          sxy = sxy + c1x*c2y;
 845          sxz = sxz + c1x*c2z;
 846
 847          syx = syx + c1y*c2x;
 848          syy = syy + c1y*c2y;
 849          syz = syz + c1y*c2z;
 850
 851          szx = szx + c1z*c2x;
 852          szy = szy + c1z*c2y;
 853          szz = szz + c1z*c2z;
 854       end do
 855
 856       xc(1) = s1x/n;
 857       xc(2) = s1y/n;
 858       xc(3) = s1z/n;
 859
 860       yc(1) = s2x/n;
 861       yc(2) = s2y/n;
 862       yc(3) = s2z/n;
 863       if(mode.eq.2.or.mode.eq.0) then ! need rmsd
 864          do m=1, n
 865             do i=1, 3
 866                e0 = e0+ (x(i, m)-xc(i))**2 + (y(i, m)-yc(i))**2
 867             end do
 868          end do
 869       endif
 870
 871       r(1, 1) = sxx-s1x*s2x/n;
 872       r(2, 1) = sxy-s1x*s2y/n;
 873       r(3, 1) = sxz-s1x*s2z/n;
 874       r(1, 2) = syx-s1y*s2x/n;
 875       r(2, 2) = syy-s1y*s2y/n;
 876       r(3, 2) = syz-s1y*s2z/n;
 877       r(1, 3) = szx-s1z*s2x/n;
 878       r(2, 3) = szy-s1z*s2y/n;
 879       r(3, 3) = szz-s1z*s2z/n;
 880
 881       det = r(1,1) * ( (r(2,2)*r(3,3)) - (r(2,3)*r(3,2)) )
 882      &     - r(1,2) * ( (r(2,1)*r(3,3)) - (r(2,3)*r(3,1)) )
 883      &     + r(1,3) * ( (r(2,1)*r(3,2)) - (r(2,2)*r(3,1)) )
 884
 885       sigma = det
 886
 887       m = 0
 888       do j=1, 3
 889          do i=1, j
 890             m = m+1
 891             rr(m) = r(1,i)*r(1,j) + r(2,i)*r(2,j) + r(3,i)*r(3,j)
 892          end do
 893       end do
 894
 895       spur = (rr(1)+rr(3)+rr(6)) / 3.0
 896       cof = (((((rr(3)*rr(6) - rr(5)*rr(5)) + rr(1)*rr(6))
 897      &     - rr(4)*rr(4)) + rr(1)*rr(3)) - rr(2)*rr(2)) / 3.0
 898       det = det*det
 899
 900       do i=1, 3
 901          e(i) = spur
 902       end do
 903       if( spur .le. zero ) goto 40
 904       d = spur*spur
 905       h = d - cof
 906       g = (spur*cof - det)/2.0 - spur*h
 907       if( h .le. zero ) then
 908          if( mode .eq. 0 ) then
 909             goto 50
 910          else
 911             goto 30
 912          end if
 913       end if
 914       sqrth = dsqrt(h)
 915       d = h*h*h - g*g
 916       if( d .lt. zero ) d = zero
 917       d = datan2( dsqrt(d), -g ) / 3.0
 918       cth = sqrth * dcos(d)
 919       sth = sqrth*sqrt3*dsin(d)
 920       e(1) = (spur + cth) + cth
 921       e(2) = (spur - cth) + sth
 922       e(3) = (spur - cth) - sth
 923
 924       if( mode .eq. 0 ) then
 925          goto 50
 926       end if
 927
 928       do l=1, 3, 2
 929          d = e(l)
 930          ss(1) = (d-rr(3)) * (d-rr(6))  - rr(5)*rr(5)
 931          ss(2) = (d-rr(6)) * rr(2)      + rr(4)*rr(5)
 932          ss(3) = (d-rr(1)) * (d-rr(6))  - rr(4)*rr(4)
 933          ss(4) = (d-rr(3)) * rr(4)      + rr(2)*rr(5)
 934          ss(5) = (d-rr(1)) * rr(5)      + rr(2)*rr(4)
 935          ss(6) = (d-rr(1)) * (d-rr(3))  - rr(2)*rr(2)
 936
 937          if( dabs(ss(1)) .ge. dabs(ss(3)) ) then
 938             j=1
 939             if( dabs(ss(1)) .lt. dabs(ss(6)) ) j = 3
 940          else if( dabs(ss(3)) .ge. dabs(ss(6)) ) then
 941             j = 2
 942          else
 943             j = 3
 944          end if
 945
 946          d = zero
 947          j = 3 * (j - 1)
 948
 949          do i=1, 3
 950             k = ip(i+j)
 951             a(i,l) = ss(k)
 952             d = d + ss(k)*ss(k)
 953          end do
 954          if( d .gt. zero ) d = 1.0 / dsqrt(d)
 955          do i=1, 3
 956             a(i,l) = a(i,l) * d
 957          end do
 958       end do
 959
 960       d = a(1,1)*a(1,3) + a(2,1)*a(2,3) + a(3,1)*a(3,3)
 961       if ((e(1) - e(2)) .gt. (e(2) - e(3))) then
 962          m1 = 3
 963          m = 1
 964       else
 965          m1 = 1
 966          m = 3
 967       endif
 968
 969       p = zero
 970       do i=1, 3
 971          a(i,m1) = a(i,m1) - d*a(i,m)
 972          p = p + a(i,m1)**2
 973       end do
 974       if( p .le. tol ) then
 975          p = 1.0
 976          do 21 i=1, 3
 977             if (p .lt. dabs(a(i,m))) goto 21
 978             p = dabs( a(i,m) )
 979             j = i
 980  21      continue
 981          k = ip2312(j)
 982          l = ip2312(j+1)
 983          p = dsqrt( a(k,m)**2 + a(l,m)**2 )
 984          if( p .le. tol ) goto 40
 985          a(j,m1) = zero
 986          a(k,m1) = -a(l,m)/p
 987          a(l,m1) =  a(k,m)/p
 988       else
 989          p = 1.0 / dsqrt(p)
 990          do i=1, 3
 991             a(i,m1) = a(i,m1)*p
 992          end do
 993       end if
 994
 995       a(1,2) = a(2,3)*a(3,1) - a(2,1)*a(3,3)
 996       a(2,2) = a(3,3)*a(1,1) - a(3,1)*a(1,3)
 997       a(3,2) = a(1,3)*a(2,1) - a(1,1)*a(2,3)
 998
 999  30   do l=1, 2
1000          d = zero
1001          do i=1, 3
1002             b(i,l) = r(i,1)*a(1,l) + r(i,2)*a(2,l) + r(i,3)*a(3,l)
1003             d = d + b(i,l)**2
1004          end do
1005          if( d .gt. zero ) d = 1.0 / dsqrt(d)
1006          do i=1, 3
1007             b(i,l) = b(i,l)*d
1008          end do
1009       end do
1010       d = b(1,1)*b(1,2) + b(2,1)*b(2,2) + b(3,1)*b(3,2)
1011       p = zero
1012
1013       do i=1, 3
1014          b(i,2) = b(i,2) - d*b(i,1)
1015          p = p + b(i,2)**2
1016       end do
1017       if( p .le. tol ) then
1018          p = 1.0
1019          do 22 i=1, 3
1020             if(p.lt.dabs(b(i,1)))goto 22
1021             p = dabs( b(i,1) )
1022             j = i
1023  22      continue
1024          k = ip2312(j)
1025          l = ip2312(j+1)
1026          p = dsqrt( b(k,1)**2 + b(l,1)**2 )
1027          if( p .le. tol ) goto 40
1028          b(j,2) = zero
1029          b(k,2) = -b(l,1)/p
1030          b(l,2) =  b(k,1)/p
1031       else
1032          p = 1.0 / dsqrt(p)
1033          do i=1, 3
1034             b(i,2) = b(i,2)*p
1035          end do
1036       end if
1037
1038       b(1,3) = b(2,1)*b(3,2) - b(2,2)*b(3,1)
1039       b(2,3) = b(3,1)*b(1,2) - b(3,2)*b(1,1)
1040       b(3,3) = b(1,1)*b(2,2) - b(1,2)*b(2,1)
1041
1042       do i=1, 3
1043          do j=1, 3
1044             u(i,j) = b(i,1)*a(j,1) + b(i,2)*a(j,2) + b(i,3)*a(j,3)
1045          end do
1046       end do
1047
1048  40   do i=1, 3
1049          t(i) = ((yc(i) - u(i,1)*xc(1)) - u(i,2)*xc(2)) - u(i,3)*xc(3)
1050       end do
1051  50   do i=1, 3
1052          if( e(i) .lt. zero ) e(i) = zero
1053          e(i) = dsqrt( e(i) )
1054       end do
1055
1056       ier = 0
1057       if( e(2) .le. (e(1) * 1.0d-05) ) ier = -1
1058
1059       d = e(3)
1060       if( sigma .lt. 0.0 ) then
1061          d = - d
1062          if( (e(2) - e(3)) .le. (e(1) * 1.0d-05) ) ier = -1
1063       end if
1064       d = (d + e(2)) + e(1)
1065
1066       if(mode .eq. 2.or.mode.eq.0) then ! need rmsd
1067          rms = (e0 - d) - d
1068          if( rms .lt. 0.0 ) rms = 0.0
1069       endif
1070
1071       return
1072       end
1073
1074