New upstream version 20170708+dfsg (8d066904) · Commits · Debian Med / tm-align

TMalign.f

+115 −106

Original line number	Diff line number	Diff line
		@@ -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)