TMalign.f

@@ -8,7 +8,7 @@
 *     please contact email: zhng@umich.edu.
 *     
 *     Reference to cite:
-*     Yang Zhang, Jeffrey Skolnick, Nucl. Acid Res. 2005 33: 2303-9
+*     Yang Zhang, Jeffrey Skolnick, Nucl. Acid Res. 2005 33: 2302-9
 *     
 *     Permission to use, copy, modify, and distribute this program for 
 *     any purpose, with or without fee, is hereby granted, provided that
@@ -74,6 +74,7 @@
 *     2014/06/01: Added 'TM.sup_all_atm_lig' to display ligand structures
 *     2015/09/14: optimized I/O which increased speed by ~100%
 *     2016/05/21: fixed a bug on conformation output
+*     2017/07/08: Added one iteration in initial4 to avoid asymmetric alignment
 **************************************************************************
       
       program TMalign
@@ -149,7 +150,7 @@ ccc
          write(*,*)
          write(*,*)'Brief instruction for running TM-align program:'
          write(*,*)'(For detail: Zhang & Skolnick, Nucl. Acid. Res.',
-     &     ' 33: 2303, 2005)'
+     &     ' 33: 2302-9, 2005)'
          write(*,*)
          write(*,*)'1. Align ''chain_1.pdb'' and ''chain_2.pdb'':'
          write(*,*)'   >TMalign chain_1.pdb chain_2.pdb'
@@ -199,7 +200,7 @@ ccc
          goto 9999
       endif
       
-      version='20160521'
+      version='20170708'
       if(fnam.eq.'-v')then
          write(*,*)'TM-align Version ',version
          goto 9999
@@ -1118,7 +1119,7 @@ ccc
       endif
       
 *11111111111111111111111111111111111111111111111111111111
-*     get initial alignment from gapless threading
+*     get initial alignment from global gapless threading
 **********************************************************
       call get_initial1          !gapless threading
       do i=1,nseq2
@@ -1293,7 +1294,7 @@ c     record the best alignment in whole search ---------->
  3333 continue
 
 *444444444444444444444444444444444444444444444444444444444
-*     get initial alignment of pieces from gapless threading
+*     initial alignment from gapless threading on largest continous fragments
 **********************************************************
       call get_initial4         !gapless threading
       do i=1,nseq2
@@ -1583,15 +1584,18 @@ c     1->coil, 2->helix, 3->turn, 4->strand
       fra_min1=fra_min-1        !cutoff for shift, save time
       dcu0=4.25
       
-ccc   Find the smallest continuous fragments -------->
       do i=1,nseq1
          mm(1,i)=mm1(i)
       enddo
       do i=1,nseq2
          mm(2,i)=mm2(i)
       enddo
-      do k=1,2
-         dcu=dcu0
+      GL_max=0
+      
+c      do k=1,2                  !k=1, fragment from protein1; k=2, fragment from protein2
+      do k=2,1,-1               !k=1, fragment from protein1; k=2, fragment from protein2
+ccc   Find the smallest continuous fragments on protein-k -------->
+         dcu=dcu0               !breaking bond-length
          if(k.eq.1)then
             nseq0=nseq1
             r_min=nseq1/3.0     !minimum fragment, in case too small protein
@@ -1600,12 +1604,14 @@ ccc   Find the smallest continuous fragments -------->
             r_min=nseq2/3.0     !minimum fragment, in case too small protein
          endif
          if(r_min.gt.fra_min)r_min=fra_min
+         
  20      nfr=1                  !number of fragments
-         j=1                    !number of residue at nf-fragment
-         ifr2(k,nfr,j)=1        !what residue
+         j=1                    !number of residues at nfr-fragment
+         ifr2(k,nfr,j)=1        !residue ID of nfr-fragment
          Lfr2(k,nfr)=j          !length of the fragment
+         
          do i=2,nseq0
-            dis=diszy(k-1,i-1,i)
+            dis=diszy(k-1,i-1,i) !str,res,res
             contin=.false.
             if(dcu.gt.dcu0)then
                if(dis.lt.dcu)then
@@ -1639,38 +1645,17 @@ ccc   Find the smallest continuous fragments -------->
             dcu=dcu+0.01
             goto 20
          endif
-      enddo
-c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
          
-ccc   select what piece will be used (this may araise ansysmetry, but
-ccc   only when L1=L2 and Lfr1=Lfr2 and L1 ne Lfr1
-ccc   if L1=Lfr1 and L2=Lfr2 (normal proteins), it will be the same as initial1
-      mark=1
-      if(Lfr2(1,i_fr2(1)).lt.Lfr2(2,i_fr2(2)))then
-         mark=1
-      elseif(Lfr2(1,i_fr2(1)).gt.Lfr2(2,i_fr2(2)))then
-         mark=2
-      else                      !Lfr1=Lfr2
-         if(nseq1.le.nseq2)then
-            mark=1
-         else
-            mark=2
-         endif
-      endif
-ccc   
-      L_fr=Lfr2(mark,i_fr2(mark))
+         L_fr=Lfr2(k,i_fr2(k))  !length of the maximum fragment
          do i=1,L_fr
-         ifr(i)=ifr2(mark,i_fr2(mark),i)
+            ifr(i)=ifr2(k,i_fr2(k),i)
          enddo
-ccc   
-      if(mark.eq.1)then         !non-redundant to get_initial1
-         nseq0=nseq1
-      else
-         nseq0=nseq2
-      endif
-      if(L_fr.eq.nseq0)then
-         n1=int(nseq0*0.1)      !0
-         n2=int(nseq0*0.89)     !2
+         
+ccc   find the best initial alignment------------>
+         if(k.eq.1)then         !using fragment from protein-1
+            if(L_fr.eq.nseq1)then !to make it different from initial1
+               n1=int(nseq1*0.1) !0
+               n2=int(nseq1*0.89) !2
                j=0
                do i=n1,n2
                   j=j+1
@@ -1679,15 +1664,13 @@ ccc
                L_fr=j
             endif
             
-ccc   get initial ------------->
-      if(mark.eq.1)then    !nseq1 as the smallest one
             nseq1_=L_fr
             aL=min(nseq1_,nseq2)
             idel=aL/2.5         !minimum size of considered fragment
             if(idel.le.fra_min1)idel=fra_min1
             n1=-nseq2+idel      !shift1
             n2=nseq1_-idel      !shift2
-         GL_max=0
+c            write(*,*)idel,aL,n1,n2,'aaa---'
             do ishift=n1,n2
                L=0
                do j=1,nseq2
@@ -1709,14 +1692,31 @@ ccc   get initial ------------->
                   endif
                endif
             enddo
+            
+c            write(*,*)'GL_max=',GL_max
+c            do i=1,nseq2
+c               write(*,*)i,invmap_i(i),'111'
+c            enddo
+            
          else
+            if(L_fr.eq.nseq2)then !to make it different from initial1
+               n1=int(nseq2*0.1) !0
+               n2=int(nseq2*0.89) !2
+               j=0
+               do i=n1,n2
+                  j=j+1
+                  ifr(j)=ifr(n1+j)
+               enddo
+               L_fr=j
+            endif
+            
             nseq2_=L_fr
             aL=min(nseq1,nseq2_)
             idel=aL/2.5         !minimum size of considered fragment
             if(idel.le.fra_min1)idel=fra_min1
             n1=-nseq2_+idel
             n2=nseq1-idel
-         GL_max=0
+c            write(*,*)idel,aL,n1,n2,'bbb----'
             do ishift=n1,n2
                L=0
                do j=1,nseq2
@@ -1739,7 +1739,14 @@ ccc   get initial ------------->
                   endif
                endif
             enddo
+
+c            write(*,*)'GL_max=',GL_max
+c            do i=1,nseq2
+c               write(*,*)i,invmap_i(i),'1112222'
+c            enddo
          endif
+      enddo
+c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
       
       return
       end
@@ -2896,6 +2903,8 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 *     one layer of the matrix. This code was exploited in TM-align 
 *     because it is about 1.5 times faster than a complete N-W code
 *     and does not influence much the final structure alignment result.
+*     In 1/1000 case, it may result in asymmetry, i.e. A_to_B!=B_to_A
+*     For example, '1se9A.pdb' and '2edpA.pdb'
 ********************************************************************
       SUBROUTINE DP(NSEQ1,NSEQ2)
       PARAMETER(nmax=5000)

debian/changelog

+tm-align (20170708+dfsg-1) unstable; urgency=medium
+
+  * New upstream version
+  * Standards-Version: 4.1.3
+  * debhelper 11
+  * d/watch: version=4 (despite upstream page supports https, watch file
+    fails if using https there)
+
+ -- Andreas Tille <tille@debian.org>  Sun, 18 Feb 2018 17:32:48 +0100
+
 tm-align (20160521+dfsg-2) unstable; urgency=medium
 
   * Team upload.

debian/compat

-9
+11

debian/control

@@ -5,9 +5,9 @@ Uploaders: Steffen Moeller <moeller@debian.org>,
            Andreas Tille <tille@debian.org>
 Section: science
 Priority: optional
-Build-Depends: debhelper (>= 9),
+Build-Depends: debhelper (>= 11~),
                gfortran
-Standards-Version: 3.9.8
+Standards-Version: 4.1.3
 Vcs-Browser: https://anonscm.debian.org/cgit/debian-med/tm-align.git
 Vcs-Git: https://anonscm.debian.org/git/debian-med/tm-align.git
 Homepage: http://zhanglab.ccmb.med.umich.edu/TM-align/

debian/watch

-version=3
+version=4
 opts="repacksuffix=+dfsg,dversionmangle=s/\+dfsg//g,repack,compression=xz" \
   http://zhanglab.ccmb.med.umich.edu/TM-align/TMtools([\d]*)\.tar\.gz