diff -Nru tm-align-20140601+dfsg/debian/changelog tm-align-20150914+dfsg/debian/changelog --- tm-align-20140601+dfsg/debian/changelog 2014-08-07 14:24:05.000000000 +0000 +++ tm-align-20150914+dfsg/debian/changelog 2015-09-23 08:41:37.000000000 +0000 @@ -1,3 +1,12 @@ +tm-align (20150914+dfsg-1) unstable; urgency=medium + + * New homepage + * New upstream version + * cme fix dpkg-control + * dversionmangle in d/watch + + -- Andreas Tille Wed, 23 Sep 2015 10:01:59 +0200 + tm-align (20140601+dfsg-1) unstable; urgency=medium * New upstream version diff -Nru tm-align-20140601+dfsg/debian/control tm-align-20150914+dfsg/debian/control --- tm-align-20140601+dfsg/debian/control 2014-02-13 16:04:40.000000000 +0000 +++ tm-align-20150914+dfsg/debian/control 2015-09-23 08:02:10.000000000 +0000 @@ -7,10 +7,10 @@ Priority: optional Build-Depends: debhelper (>= 9), gfortran -Standards-Version: 3.9.5 +Standards-Version: 3.9.6 Vcs-Browser: http://anonscm.debian.org/viewvc/debian-med/trunk/packages/zhanglab/tm-align/trunk/ Vcs-Svn: svn://anonscm.debian.org/debian-med/trunk/packages/zhanglab/tm-align/trunk/ -Homepage: http://zhang.bioinformatics.ku.edu/TM-align/ +Homepage: http://zhanglab.ccmb.med.umich.edu/TM-align/ Package: tm-align Architecture: any diff -Nru tm-align-20140601+dfsg/debian/watch tm-align-20150914+dfsg/debian/watch --- tm-align-20140601+dfsg/debian/watch 2014-08-07 14:22:11.000000000 +0000 +++ tm-align-20150914+dfsg/debian/watch 2015-09-23 08:41:19.000000000 +0000 @@ -1,3 +1,3 @@ version=3 -opts=uversionmangle=s/$/+dfsg/ \ -http://zhanglab.ccmb.med.umich.edu/TM-align/TMtools([\d]*)\.tar\.gz +opts="repacksuffix=+dfsg,dversionmangle=s/\+dfsg//g" \ + http://zhanglab.ccmb.med.umich.edu/TM-align/TMtools([\d]*)\.tar\.gz diff -Nru tm-align-20140601+dfsg/TMalign.f tm-align-20150914+dfsg/TMalign.f --- tm-align-20140601+dfsg/TMalign.f 2014-06-07 14:27:52.000000000 +0000 +++ tm-align-20150914+dfsg/TMalign.f 2015-09-15 01:51:46.000000000 +0000 @@ -72,6 +72,7 @@ * 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% ************************************************************************** program TMalign @@ -111,7 +112,7 @@ common/d8/d8 common/initial4/mm1(nmax),mm2(nmax) - character*10 aa1,ra1,aa2,ra2 + character*10 aa1,ra1,aa2,ra2,du2 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) @@ -141,7 +142,7 @@ & 'D','N','L','K','E', & 'Q','R','H','F','Y', & 'W','C'/ - + call getarg(1,fnam) if(fnam.eq.' '.or.fnam.eq.'?'.or.fnam.eq.'-h')then write(*,*) @@ -197,7 +198,7 @@ goto 9999 endif - version='20140601' + version='20150914' if(fnam.eq.'-v')then write(*,*)'TM-align Version ',version goto 9999 @@ -257,6 +258,89 @@ if(i.lt.narg)goto 115 ccccccccc read data from first CA file: + if(m_out.eq.-1)then !no output so no need to read all atoms (time-saving) + open(unit=10,file=pdb(1),status='old') + i=0 + do while (.true.) + read(10,9001,end=1013) s + if(i.gt.0.and.s(1:3).eq.'TER')goto 1013 + if(s(1:3).eq.'ATO')then + if(s(13:16).eq.'CA '.or.s(13:16).eq.' CA '.or + & .s(13:16).eq.' CA')then + du1=s(27:27) !read insertion tag + mk=1 + if(s(17:17).ne.' ')then !with Alternate atom + read(s(23:26),*)i8 !res ID + if(nres1(i8,ichar(du1)).ge.1)then + mk=-1 !this residue was already read + endif + endif + if(mk.eq.1)then + i=i+1 + read(s,9000)du,ma1(i),du,aanam,du,mm1(i),du, + $ xa(1,i,0),xa(2,i,0),xa(3,i,0) + nres1(mm1(i),ichar(du1))=nres1(mm1(i),ichar(du1))+1 + do j=-1,20 + if(aanam.eq.aa(j))then + seq1(i)=slc(j) + goto 121 + endif + enddo + seq1(i)=slc(-1) + 121 continue + ss1(i)=aanam + if(i.ge.nmax)goto 1013 + endif + endif + endif + enddo + 1013 continue + close(10) + nseq1=i + + open(unit=10,file=pdb(2),status='old') + i=0 + do while (.true.) + read(10,9001,end=1014) s + if(i.gt.0.and.s(1:3).eq.'TER')goto 1014 + if(s(1:3).eq.'ATO')then + if(s(13:16).eq.'CA '.or.s(13:16).eq.' CA '.or + & .s(13:16).eq.' CA')then + du1=s(27:27) !read insertion tag + mk=1 + if(s(17:17).ne.' ')then !with Alternate atom + read(s(23:26),*)i8 !res ID + if(nres2(i8,ichar(du1)).ge.1)then + mk=-1 !this residue was already read + endif + endif + if(mk.eq.1)then + i=i+1 + read(s,9000)du,ma2(i),du,aanam,du,mm2(i),du, + $ xa(1,i,1),xa(2,i,1),xa(3,i,1) + nres2(mm2(i),ichar(du1))=nres2(mm2(i),ichar(du1))+1 + do j=-1,20 + if(aanam.eq.aa(j))then + seq2(i)=slc(j) + goto 122 + endif + enddo + seq2(i)=slc(-1) + 122 continue + ss2(i)=aanam + if(i.ge.nmax)goto 1014 + endif + endif + endif + enddo + 1014 continue + close(10) + nseq2=i + goto 1017 + endif +c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +cccccccc read full-atom structure that is needed only for output ----> open(unit=10,file=pdb(1),status='old') i=0 na1=0 @@ -266,34 +350,29 @@ if(i.gt.0.and.s(1:3).eq.'TER')goto 1010 if(s(1:3).eq.'ATO')then *** - if(s(27:27).eq.' ')then - du1=' ' - else - read(s(27:27),*)du1 - endif + du1=s(27:27) !read insertion tag mk=1 - if(s(17:17).ne.' ')then - read(s(13:16),*)du - read(s(23:26),*)i8 - do i1=1,nres1(i8,ichar(du1)) - if(du.eq.atom1(i8,i1))then + if(s(17:17).ne.' ')then !with Alternate atom + du2=s(13:16) !atom ID + read(s(23:26),*)i8 !res ID + do i1=1,nres1(i8,ichar(du1)) !#of atoms for res_insert + if(du2.eq.atom1(i8,i1))then !such atom was already read mk=-1 endif enddo 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 1010 na1=na1+1 read(s,8999)du,ia1(na1),du,aa1(na1),du,ra1(na1),du, & ir1(na1),du,xa1(na1),ya1(na1),za1(na1) + i8=ir1(na1) + if(nres1(i8,ichar(du1)).lt.30)then + nres1(i8,ichar(du1))=nres1(i8,ichar(du1))+1 !#of atoms for res_ins + atom1(i8,nres1(i8,ichar(du1)))=aa1(na1) !atom ID + endif ains1(na1)=du1 if(s(13:16).eq.'CA '.or.s(13:16).eq.' CA '.or & .s(13:16).eq.' CA')then @@ -323,7 +402,7 @@ nseq1=i c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -ccccccccc read data from the second CA file: +ccccccccc read data from the second structure file: open(unit=10,file=pdb(2),status='old') i=0 na2=0 @@ -333,34 +412,29 @@ if(i.gt.0.and.s(1:3).eq.'TER')goto 1011 if(s(1:3).eq.'ATO')then *** - if(s(27:27).eq.' ')then - du1=' ' - else - read(s(27:27),*)du1 - endif + du1=s(27:27) !read insertion tag mk=1 if(s(17:17).ne.' ')then - read(s(13:16),*)du + du2=s(13:16) read(s(23:26),*)i8 do i1=1,nres2(i8,ichar(du1)) - if(du.eq.atom2(i8,i1))then + if(du2.eq.atom2(i8,i1))then mk=-1 endif enddo 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 1011 na2=na2+1 read(s,8999)du,ia2(na2),du,aa2(na2),du,ra2(na2),du, & ir2(na2),du,xa2(na2),ya2(na2),za2(na2) + i8=ir2(na2) + if(nres2(i8,ichar(du1)).lt.30)then + nres2(i8,ichar(du1))=nres2(i8,ichar(du1))+1 + atom2(i8,nres2(i8,ichar(du1)))=aa2(na2) + endif ains2(na2)=du1 if(s(13:16).eq.'CA '.or.s(13:16).eq.' CA '.or. & s(13:16).eq.' CA')then @@ -386,7 +460,8 @@ close(10) nseq2=i c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - + 1017 continue + ccccccccc read initial alignment file from 'alignment.txt': if(m_alignment.eq.1.or.m_alignment_stick.eq.1)then open(unit=10,file=falign,status='old') @@ -2538,22 +2613,25 @@ common/d8/d8 d_tmp=d + d_tmp2=d*d + d82=d8*d8 21 n_cut=0 !number of residue-pairs disTMscore model native' write(*,*) - write(*,*)'2. TM-score normalized with an assigned scale d0', - & ' e.g. 5 A:' + write(*,*)'2. Run TM-score with an assigned d0, e.g. 5', + & ' Angstroms:' write(*,*)' >TMscore model native -d 5' write(*,*) - 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. ', + write(*,*)'3. Run TM-score with superposition output, e.g. ', & '''TM.sup'':' write(*,*)' >TMscore model native -o TM.sup' - write(*,*)' To view the superimposed CA-traces by rasmol:' + write(*,*)' To view the superimposed structures 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 @@ -122,7 +96,6 @@ ******* options -----------> m_out=-1 m_fix=-1 - m_len=-1 narg=iargc() i=0 j=0 @@ -138,11 +111,6 @@ i=i+1 call getarg(i,fnam) read(fnam,*)d0_fix - 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 @@ -152,120 +120,53 @@ 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(s(1:4).eq.'ATOM')then -*** - 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,nres1(i8,ichar(du1)) - if(du.eq.atom1(i8,i1))then !atom is a ALTloc - mk=-1 + 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 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, - & 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,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) 104 format(A100) - 8999 format(a6,I5,a1,A4,a1,A3,a2,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(s(1:4).eq.'ATOM')then -*** - 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 + 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 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, - & 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,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 @@ -275,7 +176,7 @@ close(10) nseqB=i c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - + ****************************************************************** * pickup the aligned residues: ****************************************************************** @@ -283,74 +184,29 @@ 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 + k=k+1 + iA(k)=i + iB(k)=j + goto 205 endif enddo 205 continue enddo 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 -------------> - 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!!') - - - ************///// * parameters: ***************** *** d0-------------> - if(nseqB.gt.15)then + if(nseqB.gt.21)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 -----> @@ -406,12 +262,10 @@ 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) @@ -486,12 +340,7 @@ 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 ****************************************************************** @@ -518,14 +367,7 @@ write(*,*) write(*,501)pdb(1),nseqA 501 format('Structure1: ',A10,' Length= ',I4) - 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) + write(*,502)pdb(2),nseqB 502 format('Structure2: ',A10,' Length= ',I4, & ' (by which all scores are normalized)') write(*,503)n_ali @@ -533,25 +375,20 @@ write(*,513)rmsd_ali 513 format('RMSD of the common residues= ',F8.3) write(*,*) - 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 + write(*,504)score_max,d0,score10_max + 504 format('TM-score = ',f6.4,' (d0=',f5.2,',',' TM10= ',f6.4,')') + write(*,505)score_maxsub_max 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*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 + write(*,506)score_GDT,n_GDT1_max/float(nseqB),n_GDT2_max + & /float(nseqB),n_GDT4_max/float(nseqB),n_GDT8_max/float(nseqB) 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*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 + 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) 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(*,*) @@ -595,10 +432,23 @@ ********* 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) @@ -622,56 +472,22 @@ write(7,515)score_max,d0 515 format('REMARK TM-score=',f6.4,' (d0=',f5.2,')') do i=1,nseqA - write(7,1236)nresA(i),seqA(i),nresA(i),ins1(i), - & xt(i),yt(i),zt(i) + write(7,1237)nresA(i),seqA(i),nresA(i),xt(i),yt(i),zt(i) enddo write(7,1238) do i=2,nseqA write(7,1239)nresA(i-1),nresA(i) enddo do i=1,nseqB - write(7,1237)2000+nresB(i),seqB(i),nresB(i), - & ins2(i),xb(i),yb(i),zb(i) + write(7,1237)2000+nresB(i),seqB(i),nresB(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) enddo - close(7) - 1236 format('ATOM ',i5,' CA ',A3,' A',I4,A1,3X,3F8.3) - 1237 format('ATOM ',i5,' CA ',A3,' B',I4,A1,3X,3F8.3) + 1237 format('ATOM ',i5,' CA ',A3,I6,4X,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 *A' - write(7,900)'color blue' - write(7,900)'select *B' - 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)ia1(i),aa1(i),ra1(i),ir1(i),ains1(i),ax,ay,az - enddo - write(7,1238) !TER -*** chain2: - do i=1,na2 - write(7,8889)ia2(i),aa2(i),ra2(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,' A',I4,A1,3x,3F8.3) - 8889 format('ATOM ',I5,1x,A4,1x,A3,' B',I4,A1,3x,3F8.3) 999 continue *** record aligned residues by i=[1,nseqA], for sequenceM()------------> @@ -738,22 +554,6 @@ 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, $ ')=',i3,' RMSD=',f6.2) @@ -776,7 +576,7 @@ c 2, calculate score_GDT, score_maxsub, score_TM ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine score_fun - PARAMETER(nmax=5000) + PARAMETER(nmax=3000) common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax) common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB @@ -846,14 +646,13 @@ 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,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 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 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) @@ -874,9 +673,6 @@ 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 @@ -886,24 +682,9 @@ 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 - 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 + e0 = zero do i=1, 3 xc(i) = zero @@ -923,61 +704,29 @@ ier = -1 if( n .lt. 1 ) return ier = -2 - do m=1, n - 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; + 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 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 + 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 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; + end do 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)) ) @@ -1021,7 +770,7 @@ e(1) = (spur + cth) + cth e(2) = (spur - cth) + sth e(3) = (spur - cth) - sth - + if( mode .eq. 0 ) then goto 50 end if @@ -1164,12 +913,8 @@ end if d = (d + e(2)) + e(1) - if(mode .eq. 2.or.mode.eq.0) then ! need rmsd - rms = (e0 - d) - d - if( rms .lt. 0.0 ) rms = 0.0 - endif + rms = (e0 - d) - d + if( rms .lt. 0.0 ) rms = 0.0 return end - -