[med-svn] [tm-align] 02/06: Imported Upstream version 20160521+dfsg
Andreas Tille
tille at debian.org
Mon May 30 07:47:44 UTC 2016
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository tm-align.
commit 6ad8311e1b263831038caaa52f03850ec245acd9
Author: Andreas Tille <tille at debian.org>
Date: Mon May 30 09:27:08 2016 +0200
Imported Upstream version 20160521+dfsg
---
TMalign.f | 48 ++++-
TMscore.f | 713 +++++++++++++++++++++++++++++++++++++++++++++++---------------
2 files changed, 586 insertions(+), 175 deletions(-)
diff --git a/TMalign.f b/TMalign.f
index 65d3f7f..c186d95 100755
--- a/TMalign.f
+++ b/TMalign.f
@@ -72,7 +72,8 @@
* indicators and residue insertions.
* 2013/05/11: Fixed a bug in array overflow.
* 2014/06/01: Added 'TM.sup_all_atm_lig' to display ligand structures
-* 2015/09/14: optimize I/O and increased speed by ~100%
+* 2015/09/14: optimized I/O which increased speed by ~100%
+* 2016/05/21: fixed a bug on conformation output
**************************************************************************
program TMalign
@@ -198,7 +199,7 @@ ccc
goto 9999
endif
- version='20150914'
+ version='20160521'
if(fnam.eq.'-v')then
write(*,*)'TM-align Version ',version
goto 9999
@@ -256,6 +257,11 @@ ccc
pdb(j)=fnam
endif
if(i.lt.narg)goto 115
+
+ irmx=0 !no conformation output
+ if(m_matrix.eq.1.or.m_out.eq.1)then !we need to extract rotation matrix
+ irmx=1
+ endif
ccccccccc read data from first CA file:
if(m_out.eq.-1)then !no output so no need to read all atoms (time-saving)
@@ -592,7 +598,7 @@ c write(*,*)'---------',rms,n8_al,RMSD,u(1,1),t(1)
if(d0.lt.d0_min)d0=d0_min
d0_input=d0
call TMscore(d0_input,n8_al,xtm1,ytm1,ztm1,n8_al,
- & xtm2,ytm2,ztm2,TM8,Rcomm,Lcomm) !normal TMscore
+ & xtm2,ytm2,ztm2,TM8,Rcomm,Lcomm,0) !normal TMscore
TM1=TM8*n8_al/anseq
* Based on Chain_2===>
anseq=nseq2
@@ -604,7 +610,7 @@ c write(*,*)'---------',rms,n8_al,RMSD,u(1,1),t(1)
if(d0.lt.d0_min)d0=d0_min
d0_input=d0
call TMscore(d0_input,n8_al,xtm1,ytm1,ztm1,n8_al,
- & xtm2,ytm2,ztm2,TM8,Rcomm,Lcomm) !normal TMscore
+ & xtm2,ytm2,ztm2,TM8,Rcomm,Lcomm,irmx) !normal TMscore
TM2=TM8*n8_al/anseq
* Based on Average length===>
if(m_ave.eq.1)then
@@ -617,7 +623,7 @@ c write(*,*)'---------',rms,n8_al,RMSD,u(1,1),t(1)
if(d0.lt.d0_min)d0=d0_min
d0_input=d0
call TMscore(d0_input,n8_al,xtm1,ytm1,ztm1,n8_al,
- & xtm2,ytm2,ztm2,TM8,Rcomm,Lcomm) !normal TMscore
+ & xtm2,ytm2,ztm2,TM8,Rcomm,Lcomm,irmx) !normal TMscore
TM12=TM8*n8_al/anseq
endif
* Based on assigned length===>
@@ -631,7 +637,7 @@ c write(*,*)'---------',rms,n8_al,RMSD,u(1,1),t(1)
if(d0.lt.d0_min)d0=d0_min
d0_input=d0
call TMscore(d0_input,n8_al,xtm1,ytm1,ztm1,n8_al,
- & xtm2,ytm2,ztm2,TM8,Rcomm,Lcomm) !normal TMscore
+ & xtm2,ytm2,ztm2,TM8,Rcomm,Lcomm,irmx) !normal TMscore
TML=TM8*n8_al/anseq
endif
* Based on user-specified d0===>
@@ -640,7 +646,7 @@ c write(*,*)'---------',rms,n8_al,RMSD,u(1,1),t(1)
d0_out=d0_fix
d0_input=d0
call TMscore(d0_input,n8_al,xtm1,ytm1,ztm1,n8_al,
- & xtm2,ytm2,ztm2,TM8,Rcomm,Lcomm) !normal TMscore
+ & xtm2,ytm2,ztm2,TM8,Rcomm,Lcomm,irmx) !normal TMscore
TMfix=TM8*n8_al/nseq2
endif
@@ -1086,6 +1092,9 @@ ccc stick to the initial alignment -------------->
if(L2.lt.L)L=L2
do 6661 i=1,L
if(sequence(1)(i:i).ne.'-')i1=i1+1
+ if(i1.gt.nseq1)then
+ goto 6661
+ endif
if(sequence(2)(i:i).ne.'-')then
i2=i2+1
if(i2.gt.nseq2)then
@@ -2183,7 +2192,7 @@ cccc collect aligned residues from invmap ------->
c write(*,*)'---------',rms,n_al,RMSD,u(1,1),t(1)
call TMscore(d0_input,n_al,xtm1,ytm1,ztm1,n_al,
- & xtm2,ytm2,ztm2,TM8,Rcomm,Lcomm) !normal TMscore
+ & xtm2,ytm2,ztm2,TM8,Rcomm,Lcomm,0) !normal TMscore
TM2=TM8*n_al/anseq
c^^^^^^^^^^ TM-score calculation is done ^^^^^^^^^^^^^^^^^^^^^^^^
@@ -2661,7 +2670,7 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*************************************************************************
*** normal TM-score:
subroutine TMscore(dx,L1,x1,y1,z1,L2,x2,y2,z2,
- & TM,Rcomm,Lcomm)
+ & TM,Rcomm,Lcomm,irmx)
PARAMETER(nmax=5000)
common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax)
common/nres/nseqA,nseqB
@@ -2811,6 +2820,27 @@ c d0=1.24*(nseqB-15)**(1.0/3.0)-1.8
300 continue !for shift
333 continue !for initial length, L_ali/M
+******** return the final rotation ****************
+ if(irmx.eq.1)then !we need coordinates for output structure
+ 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
+ endif
+
TM=score_max
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
diff --git a/TMscore.f b/TMscore.f
index 030c28f..883eeae 100644
--- a/TMscore.f
+++ b/TMscore.f
@@ -25,10 +25,26 @@
* 2010/08/02: A new RMSD matrix was used and obsolete statement removed.
* 2011/01/03: The length of pdb file names were extended to 500.
* 2011/01/30: An open source license is attached to the program.
+* 2012/05/07: Improved RMSD calculation subroutine which speeds up
+* TM-score program by 30%.
+* 2012/06/05: Added option '-l L' which calculates TM-score (and maxsub
+* and GDT scores) normalized by a specific length 'L'.
+* 2012/12/17: Added 'TM.sup_atm' to superpose full-atom structures.
+* The former superposition is for CA-trace only.
+* 2013/05/08: Update TM-score so that it can read all alternate location
+* indicators and residue insertions.
+* 2013/05/11: Fix a bug in array overflow.
+* 2016/03/23: Extended the program to allow calculating TM-score for for
+* complex structure comparisons, where multiple-chains are
+* merged into a single chain. Chain ID is now included in
+* the output files.
*************************************************************************
+c 1 2 3 4 5 6 7 !
+c 3456789012345678901234567890123456789012345678901234567890123456789012345678
+
program TMscore
- PARAMETER(nmax=3000)
+ PARAMETER(nmax=5000)
common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax)
common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
@@ -40,15 +56,28 @@
character*500 fnam,pdb(100),outname
character*3 aa(-1:20),seqA(nmax),seqB(nmax)
character*500 s,du
- character seq1A(nmax),seq1B(nmax),ali(nmax)
+ character*1 chA(nmax),chB(nmax)
+ character*1 chain1(90000),chain2(90000)
+ character seq1A(nmax),seq1B(nmax),ali(nmax),du1
character sequenceA(nmax),sequenceB(nmax),sequenceM(nmax)
-
+ character ins1(nmax),ins2(nmax),ains1(90000),ains2(90000)
+
dimension L_ini(100),iq(nmax)
common/scores/score,score_maxsub,score_fix,score10
common/GDT/n_GDT05,n_GDT1,n_GDT2,n_GDT4,n_GDT8
double precision score,score_max,score_fix,score_fix_max
double precision score_maxsub,score10
dimension xa(nmax),ya(nmax),za(nmax)
+
+ character*10 aa1,ra1,aa2,ra2
+ dimension ia1(90000),aa1(90000),ra1(90000),ir1(90000)
+ dimension xa1(90000),ya1(90000),za1(90000)
+ dimension ia2(90000),aa2(90000),ra2(90000),ir2(90000)
+ dimension xa2(90000),ya2(90000),za2(90000)
+
+ dimension nres1(nmax,32:122),nres2(nmax,32:122) !number of atoms
+ character*500 atom1(nmax,30),atom2(nmax,30) !atom name
+c here atom1(i,j) should be atom2(i,j,k) but will beyond dimension
ccc RMSD:
double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)
@@ -76,19 +105,29 @@ ccc
write(*,*)'(For detail: Zhang & Skolnick, Proteins, 2004',
& ' 57:702-10)'
write(*,*)
- write(*,*)'1. Run TM-score to compare ''model'' and ''native',
- & ':'
+ write(*,*)'1. Run TM-score to compare ''model'' and ',
+ & '''native'':'
write(*,*)' >TMscore model native'
write(*,*)
- write(*,*)'2. Run TM-score with an assigned d0, e.g. 5',
- & ' Angstroms:'
+ write(*,*)'2. Run TM-score to compare two complex structures',
+ & ' with multiple chains'
+ write(*,*)' >TMscore -c model native'
+ write(*,*)
+ write(*,*)'2. TM-score normalized with an assigned scale d0',
+ & ' e.g. 5 A:'
write(*,*)' >TMscore model native -d 5'
write(*,*)
- write(*,*)'3. Run TM-score with superposition output, e.g. ',
+ write(*,*)'3. TM-score normalized by a specific length, ',
+ & 'e.g. 120 AA:'
+ write(*,*)' >TMscore model native -l 120'
+ write(*,*)
+ write(*,*)'4. TM-score with superposition output, e.g. ',
& '''TM.sup'':'
write(*,*)' >TMscore model native -o TM.sup'
- write(*,*)' To view the superimposed structures by rasmol:'
+ write(*,*)' To view the superimposed CA-traces by rasmol:'
write(*,*)' >rasmol -script TM.sup'
+ write(*,*)' To view superimposed atomic models by rasmol:'
+ write(*,*)' >rasmol -script TM.sup_atm'
write(*,*)
goto 9999
endif
@@ -96,9 +135,11 @@ ccc
******* options ----------->
m_out=-1
m_fix=-1
+ m_len=-1
narg=iargc()
i=0
j=0
+ m_complex=0
115 continue
i=i+1
call getarg(i,fnam)
@@ -111,102 +152,244 @@ ccc
i=i+1
call getarg(i,fnam)
read(fnam,*)d0_fix
+ elseif(fnam.eq.'-c')then
+ m_complex=1
+ elseif(fnam.eq.'-l')then
+ m_len=1
+ i=i+1
+ call getarg(i,fnam)
+ read(fnam,*)l0_fix
else
j=j+1
pdb(j)=fnam
endif
if(i.lt.narg)goto 115
-
+
ccccccccc read data from first CA file:
+
open(unit=10,file=pdb(1),status='old')
i=0
+ na1=0
+ ntt=0
101 read(10,104,end=102) s
- if(s(1:3).eq.'TER') goto 102
+ if(m_complex.eq.0)then
+ if(s(1:3).eq.'TER') goto 102
+ endif
if(s(1:4).eq.'ATOM')then
- if(s(13:16).eq.'CA '.or.s(13:16).eq.' CA '.or.s(13:16).
- & eq.' CA')then
- if(s(17:17).eq.' '.or.s(17:17).eq.'A')then
- i=i+1
- read(s,103)du,seqA(i),du,nresA(i),du,xa(i),ya(i),za(i)
- do j=-1,20
- if(seqA(i).eq.aa(j))then
- seq1A(i)=slc(j)
- goto 21
+***
+ if(s(27:27).eq.' ')then
+ du1=' '
+ else
+ read(s(27:27),*)du1 !Residue_Insertion_tag
+ endif
+ mk=1
+ if(s(17:17).ne.' ')then !ignor alternative residues/atoms
+ read(s(13:16),*)du
+ read(s(23:26),*)i8
+ do i1=1,nres1(i8,ichar(du1))
+ if(du.eq.atom1(i8,i1))then !atom is a ALTloc
+ mk=-1
endif
enddo
- seq1A(i)=slc(-1)
- 21 continue
endif
+ if(mk.eq.1)then
+ read(s(13:16),*)du
+ read(s(23:26),*)i8
+ if(nres1(i8,ichar(du1)).lt.30)then
+ nres1(i8,ichar(du1))=nres1(i8,ichar(du1))+1
+ atom1(i8,nres1(i8,ichar(du1)))=du
+ endif
+***
+ ntt=ntt+1
+ if(ntt.ge.90000)goto 102
+ na1=na1+1
+ read(s,8999)du,ia1(na1),du,aa1(na1),du,ra1(na1),du,
+ & chain1(na1),ir1(na1),du,xa1(na1),ya1(na1),za1(na1)
+ ains1(na1)=du1
+ if(s(13:16).eq.'CA '.or.s(13:16).eq.' CA '.or.s(13:16).
+ & eq.' CA')then
+ i=i+1
+ read(s,103)du,seqA(i),du,chA(i),nresA(i),du,
+ & xa(i),ya(i),za(i)
+ ins1(i)=du1
+ do j=-1,20
+ if(seqA(i).eq.aa(j))then
+ seq1A(i)=slc(j)
+ goto 21
+ endif
+ enddo
+ seq1A(i)=slc(-1)
+ 21 continue
+ if(i.ge.nmax)goto 102
+ endif
endif
endif
goto 101
102 continue
- 103 format(A17,A3,A2,i4,A4,3F8.3)
+ 103 format(A17,A3,A1,A1,i4,A4,3F8.3)
104 format(A100)
+ 8999 format(a6,I5,a1,A4,a1,A3,a1,a1,I4,a4,3F8.3)
close(10)
nseqA=i
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-
+
ccccccccc read data from first CA file:
open(unit=10,file=pdb(2),status='old')
i=0
+ na2=0
+ ntt=0
201 read(10,204,end=202) s
- if(s(1:3).eq.'TER') goto 202
+ if(m_complex.eq.0)then
+ if(s(1:3).eq.'TER') goto 202
+ endif
if(s(1:4).eq.'ATOM')then
- if(s(13:16).eq.'CA '.or.s(13:16).eq.' CA '.or.s(13:16).
- & eq.' CA')then
- if(s(17:17).eq.' '.or.s(17:17).eq.'A')then
- i=i+1
- read(s,203)du,seqB(i),du,nresB(i),du,xb(i),yb(i),zb(i)
- do j=-1,20
- if(seqB(i).eq.aa(j))then
- seq1B(i)=slc(j)
- goto 22
+***
+ if(s(27:27).eq.' ')then
+ du1=' '
+ else
+ read(s(27:27),*)du1
+ endif
+ mk=1
+ if(s(17:17).ne.' ')then
+ read(s(13:16),*)du
+ read(s(23:26),*)i8
+ do i1=1,nres2(i8,ichar(du1))
+ if(du.eq.atom2(i8,i1))then
+ mk=-1
endif
enddo
- seq1B(i)=slc(-1)
- 22 continue
endif
+ if(mk.eq.1)then
+ read(s(13:16),*)du
+ read(s(23:26),*)i8
+ if(nres2(i8,ichar(du1)).lt.30)then
+ nres2(i8,ichar(du1))=nres2(i8,ichar(du1))+1
+ atom2(i8,nres2(i8,ichar(du1)))=du
+ endif
+***
+ ntt=ntt+1
+ if(ntt.ge.90000)goto 202
+ na2=na2+1
+ read(s,8999)du,ia2(na2),du,aa2(na2),du,ra2(na2),du,
+ & chain2(na2),ir2(na2),du,xa2(na2),ya2(na2),za2(na2)
+ ains2(na2)=du1
+ if(s(13:16).eq.'CA '.or.s(13:16).eq.' CA '.or.s(13:16).
+ & eq.' CA')then
+ i=i+1
+ read(s,203)du,seqB(i),du,chB(i),nresB(i),du,
+ & xb(i),yb(i),zb(i)
+ ins2(i)=du1
+ do j=-1,20
+ if(seqB(i).eq.aa(j))then
+ seq1B(i)=slc(j)
+ goto 22
+ endif
+ enddo
+ seq1B(i)=slc(-1)
+ 22 continue
+ if(i.ge.nmax)goto 202
+ endif
endif
endif
goto 201
202 continue
- 203 format(A17,A3,A2,i4,A4,3F8.3)
+ 203 format(A17,A3,A1,A1,i4,A4,3F8.3)
204 format(A100)
close(10)
nseqB=i
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-
+
******************************************************************
* 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
+ if(m_complex.eq.0)then !monomer
+ k=0
+ do i=1,nseqA
+ do j=1,nseqB
+ if(nresA(i).eq.nresB(j))then
+ if(ins1(i).eq.ins2(j))then
+ k=k+1
+ iA(k)=i
+ iB(k)=j
+ goto 205
+ endif
+ endif
+ enddo
+ 205 continue
enddo
- 205 continue
- enddo
+ else !complex
+ k=0
+ do i=1,nseqA
+ do j=1,nseqB
+ if(nresA(i).eq.nresB(j).and.chA(i).eq.chB(j))then
+ if(ins1(i).eq.ins2(j))then !Residue_insertion_code
+ k=k+1
+ iA(k)=i
+ iB(k)=j
+ goto 206
+ endif
+ endif
+ enddo
+ 206 continue
+ enddo
+ endif
n_ali=k !number of aligned residues
if(n_ali.lt.1)then
- write(*,*)'There is no common residues in the input structures'
- goto 9999
+ write(*,*)'There is no common residues in the input structures'
+ goto 9999
+ endif
+
+******************************************************************
+* check the residue serial numbers ------------->
+ if(m_complex.eq.0)then
+ nA_repeat=0
+ do i=1,nseqA
+ do j=i+1,nseqA
+ if(nresA(i).eq.nresA(j))then
+ if(ins1(i).eq.ins1(j))then
+ nA_repeat=nA_repeat+1
+ endif
+ endif
+ enddo
+ enddo
+ if(nA_repeat.gt.0)then
+ write(*,380)nA_repeat
+ endif
+ 380 format('Warning: TMscore calculation can be wrong because ',
+ & 'there are ',I3,' residues with the serial number same ',
+ & 'as others in first Chain. Please modify the PDB file ',
+ & 'and rerun the program!!')
+ nB_repeat=0
+ do i=1,nseqB
+ do j=i+1,nseqB
+ if(nresB(i).eq.nresB(j))then
+ if(ins2(i).eq.ins2(j))then
+ nB_repeat=nB_repeat+1
+ endif
+ endif
+ enddo
+ enddo
+ if(nB_repeat.gt.0)then
+ write(*,381)nB_repeat
+ endif
+ 381 format('Warning: TMscore calculation can be wrong because ',
+ & 'there are ',I3,' residues with the serial number same ',
+ & 'as others in second Chain. Please modify the PDB file ',
+ & 'and rerun the program!!')
endif
************/////
* parameters:
*****************
*** d0------------->
- if(nseqB.gt.21)then
+ if(nseqB.gt.15)then
d0=1.24*(nseqB-15)**(1.0/3.0)-1.8
else
d0=0.5
endif
+ if(m_len.eq.1)then
+ d0=1.24*(l0_fix-15)**(1.0/3.0)-1.8
+ endif
if(d0.lt.0.5)d0=0.5
if(m_fix.eq.1)d0=d0_fix
*** d0_search ----->
@@ -262,10 +445,12 @@ c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
k_ali(ka)=k
LL=LL+1
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
+ call u3b(w,r_1,r_2,LL,2,rms,u,t,ier) !0:rmsd; 1:u,t; 2:rmsd,u,t
armsd=dsqrt(rms/LL)
rmsd_ali=armsd
+ else
+ call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2
endif
do j=1,nseqA
xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
@@ -340,7 +525,12 @@ c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
302 continue
300 continue !for shift
333 continue !for initial length, L_ali/M
-
+
+ ratio=1
+ if(m_len.gt.0)then
+ ratio=float(nseqB)/float(l0_fix)
+ endif
+
******************************************************************
* Output
******************************************************************
@@ -367,7 +557,14 @@ c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
write(*,*)
write(*,501)pdb(1),nseqA
501 format('Structure1: ',A10,' Length= ',I4)
- write(*,502)pdb(2),nseqB
+ if(m_len.eq.1)then
+ write(*,411)pdb(2),nseqB
+ write(*,412)l0_fix
+ else
+ write(*,502)pdb(2),nseqB
+ endif
+ 411 format('Structure2: ',A10,' Length= ',I4)
+ 412 format('TM-score is notmalized by ',I4)
502 format('Structure2: ',A10,' Length= ',I4,
& ' (by which all scores are normalized)')
write(*,503)n_ali
@@ -375,20 +572,25 @@ c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
write(*,513)rmsd_ali
513 format('RMSD of the common residues= ',F8.3)
write(*,*)
- write(*,504)score_max,d0,score10_max
- 504 format('TM-score = ',f6.4,' (d0=',f5.2,',',' TM10= ',f6.4,')')
- write(*,505)score_maxsub_max
+ if(m_len.eq.1)then
+ score_max=score_max*float(nseqB)/float(l0_fix)
+ endif
+ write(*,504)score_max,d0
+ 504 format('TM-score = ',f6.4,' (d0=',f5.2,')')
+ write(*,505)score_maxsub_max*ratio
505 format('MaxSub-score= ',f6.4,' (d0= 3.50)')
score_GDT=(n_GDT1_max+n_GDT2_max+n_GDT4_max+n_GDT8_max)
& /float(4*nseqB)
- write(*,506)score_GDT,n_GDT1_max/float(nseqB),n_GDT2_max
- & /float(nseqB),n_GDT4_max/float(nseqB),n_GDT8_max/float(nseqB)
+ write(*,506)score_GDT*ratio,n_GDT1_max/float(nseqB)*ratio,
+ & n_GDT2_max/float(nseqB)*ratio,n_GDT4_max/float(nseqB)*ratio,
+ & n_GDT8_max/float(nseqB)*ratio
506 format('GDT-TS-score= ',f6.4,' %(d<1)=',f6.4,' %(d<2)=',f6.4,
$ ' %(d<4)=',f6.4,' %(d<8)=',f6.4)
score_GDT_HA=(n_GDT05_max+n_GDT1_max+n_GDT2_max+n_GDT4_max)
& /float(4*nseqB)
- write(*,507)score_GDT_HA,n_GDT05_max/float(nseqB),n_GDT1_max
- & /float(nseqB),n_GDT2_max/float(nseqB),n_GDT4_max/float(nseqB)
+ write(*,507)score_GDT_HA*ratio,n_GDT05_max/float(nseqB)*ratio,
+ & n_GDT1_max/float(nseqB)*ratio,n_GDT2_max/float(nseqB)*ratio,
+ & n_GDT4_max/float(nseqB)*ratio
507 format('GDT-HA-score= ',f6.4,' %(d<0.5)=',f6.4,' %(d<1)=',f6.4,
$ ' %(d<2)=',f6.4,' %(d<4)=',f6.4)
write(*,*)
@@ -432,37 +634,25 @@ c enddo
********* rmsd in superposed regions --------------->
d=d_output !for output
call score_fun() !give i_ali(i), score_max=score now
- LL=0
- do i=1,n_cut
- m=i_ali(i) ![1,nseqA]
- 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,0,rms,u,t,ier)
- armsd=dsqrt(rms/LL)
- rmsd=armsd
*** output rotated chain1 + chain2----->
if(m_out.ne.1)goto 999
+ccc output CA-trace superposition------>
OPEN(unit=7,file=outname,status='unknown') !pdb1.aln + pdb2.aln
900 format(A)
- 901 format('select ',I4)
+ 901 format('select atomno= ',I5)
+ 902 format('select atomno= ',I5)
write(7,900)'load inline'
- write(7,900)'select atomno<1000'
-c write(7,900)'color [255,20,147]'
+ write(7,900)'select atomno<50001'
write(7,900)'wireframe .45'
write(7,900)'select none'
- write(7,900)'select atomno>1000'
-c write(7,900)'color [100,149,237]'
+ write(7,900)'select atomno>50000'
write(7,900)'wireframe .15'
- write(7,900)'color white'
+ write(7,900)'color white' !all above atoms are white
do i=1,n_cut
- write(7,901)nresA(iA(i_ali(i)))
+ write(7,901)iA(i_ali(i))
+ write(7,900)'color red'
+ write(7,902)50000+iB(i_ali(i))
write(7,900)'color red'
enddo
write(7,900)'select all'
@@ -472,22 +662,54 @@ c write(7,900)'color [100,149,237]'
write(7,515)score_max,d0
515 format('REMARK TM-score=',f6.4,' (d0=',f5.2,')')
do i=1,nseqA
- write(7,1237)nresA(i),seqA(i),nresA(i),xt(i),yt(i),zt(i)
+ write(7,1236)i,seqA(i),chA(i),nresA(i),ins1(i),
+ & xt(i),yt(i),zt(i)
enddo
write(7,1238)
do i=2,nseqA
- write(7,1239)nresA(i-1),nresA(i)
+ write(7,1239)i-1,i
enddo
do i=1,nseqB
- write(7,1237)2000+nresB(i),seqB(i),nresB(i),xb(i),yb(i),zb(i)
+ write(7,1236)50000+i,seqB(i),chB(i),nresB(i),ins2(i),
+ & xb(i),yb(i),zb(i)
enddo
write(7,1238)
do i=2,nseqB
- write(7,1239)2000+nresB(i-1),2000+nresB(i)
+ write(7,1239)50000+i-1,50000+i
enddo
- 1237 format('ATOM ',i5,' CA ',A3,I6,4X,3F8.3)
+ close(7)
+ 1236 format('ATOM ',i5,' CA ',A3,' ',A1,I4,A1,3X,3F8.3)
1238 format('TER')
1239 format('CONECT',I5,I5)
+ccc output all-atom superposition------>
+ outname=trim(outname)//'_atm'
+ OPEN(unit=7,file=outname,status='unknown') !pdb1.aln + pdb2.aln
+ write(7,900)'load inline'
+ write(7,900)'select atomno<50001'
+ write(7,900)'color blue'
+ write(7,900)'select atomno>50000'
+ write(7,900)'color red'
+ write(7,900)'select all'
+ write(7,900)'cartoon'
+ write(7,900)'exit'
+ write(7,514)rmsd_ali
+ write(7,515)score_max,d0
+*** chain1:
+ do i=1,na1
+ ax=t(1)+u(1,1)*xa1(i)+u(1,2)*ya1(i)+u(1,3)*za1(i)
+ ay=t(2)+u(2,1)*xa1(i)+u(2,2)*ya1(i)+u(2,3)*za1(i)
+ az=t(3)+u(3,1)*xa1(i)+u(3,2)*ya1(i)+u(3,3)*za1(i)
+ write(7,8888)i,aa1(i),ra1(i),chain1(i),ir1(i),ains1(i),ax,ay,az
+ enddo
+ write(7,1238) !TER
+*** chain2:
+ do i=1,na2
+ write(7,8888)50000+i,aa2(i),ra2(i),chain2(i),ir2(i),ains2(i),
+ & xa2(i),ya2(i),za2(i)
+ enddo
+ write(7,1238) !TER
+ close(7)
+ 8888 format('ATOM ',I5,1x,A4,1x,A3,' ',A1,I4,A1,3x,3F8.3)
999 continue
*** record aligned residues by i=[1,nseqA], for sequenceM()------------>
@@ -498,61 +720,165 @@ c write(7,900)'color [100,149,237]'
j=iA(i_ali(i)) ![1,nseqA]
k=iB(i_ali(i)) ![1,nseqB]
dis=sqrt((xt(j)-xb(k))**2+(yt(j)-yb(k))**2+(zt(j)-zb(k))**2)
+c write(*,*)i,j,k,dis,d_output,'1--'
if(dis.lt.d_output)then
iq(j)=1
endif
enddo
*******************************************************************
*** output aligned sequences
- k=0
- i=1
- j=1
- 800 continue
- if(i.gt.nseqA.and.j.gt.nseqB)goto 802
- if(i.gt.nseqA.and.j.le.nseqB)then
- k=k+1
- sequenceA(k)='-'
- sequenceB(k)=seq1B(j)
- sequenceM(k)=' '
- j=j+1
- goto 800
- endif
- if(i.le.nseqA.and.j.gt.nseqB)then
- k=k+1
- sequenceA(k)=seq1A(i)
- sequenceB(k)='-'
- sequenceM(k)=' '
- i=i+1
- goto 800
- endif
- if(nresA(i).eq.nresB(j))then
- k=k+1
- sequenceA(k)=seq1A(i)
- sequenceB(k)=seq1B(j)
- if(iq(i).eq.1)then
- sequenceM(k)=':'
- else
+ if(m_complex.eq.0)then
+ k=0
+ i=1
+ j=1
+ 800 continue
+ if(i.gt.nseqA.and.j.gt.nseqB)goto 802
+ if(i.gt.nseqA.and.j.le.nseqB)then
+ k=k+1
+ sequenceA(k)='-'
+ sequenceB(k)=seq1B(j)
sequenceM(k)=' '
+ j=j+1
+ goto 800
endif
- i=i+1
- j=j+1
- goto 800
- elseif(nresA(i).lt.nresB(j))then
- k=k+1
- sequenceA(k)=seq1A(i)
- sequenceB(k)='-'
- sequenceM(k)=' '
- i=i+1
- goto 800
- elseif(nresB(j).lt.nresA(i))then
- k=k+1
- sequenceA(k)='-'
- sequenceB(k)=seq1B(j)
- sequenceM(k)=' '
- j=j+1
- goto 800
+ if(i.le.nseqA.and.j.gt.nseqB)then
+ k=k+1
+ sequenceA(k)=seq1A(i)
+ sequenceB(k)='-'
+ sequenceM(k)=' '
+ i=i+1
+ goto 800
+ endif
+ if(nresA(i).eq.nresB(j))then
+ k=k+1
+ sequenceA(k)=seq1A(i)
+ sequenceB(k)=seq1B(j)
+ if(iq(i).eq.1)then
+ sequenceM(k)=':'
+ else
+ sequenceM(k)=' '
+ endif
+ i=i+1
+ j=j+1
+ goto 800
+ elseif(nresA(i).lt.nresB(j))then
+ k=k+1
+ sequenceA(k)=seq1A(i)
+ sequenceB(k)='-'
+ sequenceM(k)=' '
+ i=i+1
+ goto 800
+ elseif(nresB(j).lt.nresA(i))then
+ k=k+1
+ sequenceA(k)='-'
+ sequenceB(k)=seq1B(j)
+ sequenceM(k)=' '
+ j=j+1
+ goto 800
+ endif
+ 802 continue
+ else
+ k=0
+ i=1 !
+ j=1
+ 803 continue
+ if(i.gt.nseqA.and.j.gt.nseqB)goto 804
+ if(i.gt.nseqA.and.j.le.nseqB)then
+ k=k+1
+ sequenceA(k)='-'
+ sequenceB(k)=seq1B(j)
+ sequenceM(k)=' '
+ j=j+1
+ goto 803
+ endif
+ if(i.le.nseqA.and.j.gt.nseqB)then
+ k=k+1
+ sequenceA(k)=seq1A(i)
+ sequenceB(k)='-'
+ sequenceM(k)=' '
+ i=i+1
+ goto 803
+ endif
+ if(chA(i).ne.chB(j))then
+ kk=0 !check if chA(i) match later chains
+ do i2=j,nseqB
+ if(chA(i).eq.chB(i2))then
+ kk=kk+1
+ endif
+ enddo
+ if(kk.eq.0)then !chA(i) is not matched
+ k=k+1
+ sequenceA(k)=seq1A(i)
+ sequenceB(k)='-'
+ sequenceM(k)=' '
+ i=i+1
+ goto 803
+ endif
+
+ kk=0 !check if chB(j) match later chains
+ do i1=i,nseqA
+ if(chA(i1).eq.chB(j))then
+ kk=kk+1
+ endif
+ enddo
+ if(kk.eq.0)then !chB(j) is not matched
+ k=k+1
+ sequenceA(k)='-'
+ sequenceB(k)=seq1B(j)
+ sequenceM(k)=' '
+ j=j+1
+ goto 803
+ endif
+ write(*,*)'Chains in complexes have crossed order, ',
+ & 'therefore, no alignment is output'
+ stop
+ else
+ if(nresA(i).eq.nresB(j))then
+ k=k+1
+ sequenceA(k)=seq1A(i)
+ sequenceB(k)=seq1B(j)
+ if(iq(i).eq.1)then
+ sequenceM(k)=':'
+ else
+ sequenceM(k)=' '
+ endif
+ i=i+1
+ j=j+1
+ goto 803
+ elseif(nresA(i).lt.nresB(j))then
+ k=k+1
+ sequenceA(k)=seq1A(i)
+ sequenceB(k)='-'
+ sequenceM(k)=' '
+ i=i+1
+ goto 803
+ elseif(nresB(j).lt.nresA(i))then
+ k=k+1
+ sequenceA(k)='-'
+ sequenceB(k)=seq1B(j)
+ sequenceM(k)=' '
+ j=j+1
+ goto 803
+ endif
+ endif
+ 804 continue
endif
- 802 continue
+
+ccc RMSD (d<5.0)-------->
+ LL=0
+ do i=1,n_cut
+ m=i_ali(i) ![1,nseqA]
+ 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,0,rms,u,t,ier)
+ armsd=dsqrt(rms/LL)
+ rmsd=armsd
write(*,600)d_output,n_cut,rmsd
600 format('Superposition in the TM-score: Length(d<',f3.1,
@@ -576,7 +902,7 @@ c 1, collect those residues with dis<d;
c 2, calculate score_GDT, score_maxsub, score_TM
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine score_fun
- PARAMETER(nmax=3000)
+ PARAMETER(nmax=5000)
common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax)
common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
@@ -646,13 +972,14 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
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 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,short)
+c 1:calculate u,t only (given,medium)
+c 2:calculate rms,u,t (given,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)
@@ -673,6 +1000,9 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
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 c1x, c1y, c1z, c2x, c2y, c2z
+ double precision s1x, s1y, s1z, s2x, s2y, s2z
+ double precision sxx, sxy, sxz, syx, syy, syz, szx, szy, szz
double precision sqrt3, tol, zero
@@ -682,9 +1012,24 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
data ip / 1, 2, 4, 2, 3, 5, 4, 5, 6 /
data ip2312 / 2, 3, 1, 2 /
- wc = zero
+ wc = zero
rms = zero
- e0 = zero
+ e0 = zero
+ s1x = zero
+ s1y = zero
+ s1z = zero
+ s2x = zero
+ s2y = zero
+ s2z = zero
+ sxx = zero
+ sxy = zero
+ sxz = zero
+ syx = zero
+ syy = zero
+ syz = zero
+ szx = zero
+ szy = zero
+ szz = zero
do i=1, 3
xc(i) = zero
@@ -704,29 +1049,61 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
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
+ c1x=x(1, m)
+ c1y=x(2, m)
+ c1z=x(3, m)
+
+ c2x=y(1, m)
+ c2y=y(2, m)
+ c2z=y(3, m)
+
+ s1x = s1x + c1x
+ s1y = s1y + c1y;
+ s1z = s1z + c1z;
+
+ s2x = s2x + c2x;
+ s2y = s2y + c2y;
+ s2z = s2z + c2z;
+
+ sxx = sxx + c1x*c2x;
+ sxy = sxy + c1x*c2y;
+ sxz = sxz + c1x*c2z;
+
+ syx = syx + c1y*c2x;
+ syy = syy + c1y*c2y;
+ syz = syz + c1y*c2z;
+
+ szx = szx + c1z*c2x;
+ szy = szy + c1z*c2y;
+ szz = szz + c1z*c2z;
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
+ xc(1) = s1x/n;
+ xc(2) = s1y/n;
+ xc(3) = s1z/n;
+
+ yc(1) = s2x/n;
+ yc(2) = s2y/n;
+ yc(3) = s2z/n;
+ if(mode.eq.2.or.mode.eq.0) then ! need rmsd
+ do m=1, n
+ do i=1, 3
+ e0 = e0+ (x(i, m)-xc(i))**2 + (y(i, m)-yc(i))**2
+ end do
end do
- end do
+ endif
+
+ r(1, 1) = sxx-s1x*s2x/n;
+ r(2, 1) = sxy-s1x*s2y/n;
+ r(3, 1) = sxz-s1x*s2z/n;
+ r(1, 2) = syx-s1y*s2x/n;
+ r(2, 2) = syy-s1y*s2y/n;
+ r(3, 2) = syz-s1y*s2z/n;
+ r(1, 3) = szx-s1z*s2x/n;
+ r(2, 3) = szy-s1z*s2y/n;
+ r(3, 3) = szz-s1z*s2z/n;
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)) )
@@ -770,7 +1147,7 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
e(1) = (spur + cth) + cth
e(2) = (spur - cth) + sth
e(3) = (spur - cth) - sth
-
+
if( mode .eq. 0 ) then
goto 50
end if
@@ -913,8 +1290,12 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
end if
d = (d + e(2)) + e(1)
- rms = (e0 - d) - d
- if( rms .lt. 0.0 ) rms = 0.0
+ if(mode .eq. 2.or.mode.eq.0) then ! need rmsd
+ rms = (e0 - d) - d
+ if( rms .lt. 0.0 ) rms = 0.0
+ endif
return
end
+
+
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/tm-align.git
More information about the debian-med-commit
mailing list