[med-svn] [Git][med-team/tm-align][master] 7 commits: New upstream version 20190708+dfsg
Steffen Möller
gitlab at salsa.debian.org
Mon Jul 29 22:13:26 BST 2019
Steffen Möller pushed to branch master at Debian Med / tm-align
Commits:
d848d030 by Steffen Moeller at 2019-07-29T21:00:44Z
New upstream version 20190708+dfsg
- - - - -
4b941a23 by Steffen Moeller at 2019-07-29T21:00:44Z
Update upstream source from tag 'upstream/20190708+dfsg'
Update to upstream version '20190708+dfsg'
with Debian dir 8aa336566445c0987325a88e6966a8052e2cd611
- - - - -
1fc394be by Steffen Moeller at 2019-07-29T21:00:44Z
New upstream version
- - - - -
31cd74bb by Steffen Moeller at 2019-07-29T21:00:45Z
debhelper-compat 12
- - - - -
f2bbc56b by Steffen Moeller at 2019-07-29T21:00:47Z
Standards-Version: 4.4.0
- - - - -
02d3e5b9 by Steffen Moeller at 2019-07-29T21:01:06Z
Upload to unstable
- - - - -
ee122aa0 by Steffen Moeller at 2019-07-29T21:09:18Z
Have a bit of verbosity while building
- - - - -
6 changed files:
- TMalign.f
- TMscore.f
- debian/changelog
- − debian/compat
- debian/control
- debian/rules
Changes:
=====================================
TMalign.f
=====================================
@@ -72,11 +72,16 @@
* 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: optimized I/O which increased speed by ~100%
-* 2016/05/21: fixed a bug on conformation output
+* 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
+* 2019/07/08: Enable TM-align to support both PDB and mmCIF formats, and
+* fixed a bug on file output
**************************************************************************
+c 1 2 3 4 5 6 7 !
+c 3456789012345678901234567890123456789012345678901234567890123456789012345678
+
program TMalign
PARAMETER(nmax=5000) !maximum length of the sequence
PARAMETER(nmax2=10000) !for alignment output
@@ -92,27 +97,37 @@
common/alignment/m_alignment,sequence(10),TM_ali,L_ali,rmsd_ali
common/alignment1/m_alignment_stick
character*10000 sequence
- character ins1(nmax),ins2(nmax),ains1(90000),ains2(90000)
-
- common/sequence/seq1(0:nmax),seq2(0:nmax)
- character seq1,seq2,du1
+ character seq1(0:nmax),seq2(0:nmax),du1,du3,du4*2
character*500 fnam,pdb(100),outname,falign,fmatrix
- character*3 aa(-1:20),aanam,ss1(nmax),ss2(nmax)
+ character*3 aa(-1:20),aanam,ss(2,nmax)
character*500 s,du,dum1,dum2
character*504 outnameall_tmp,outnameall
character aseq1(nmax2),aseq2(nmax2),aseq3(nmax2)
character*8 version
character*5000 s1 !maximum length of protein is 5000
+
+ccc mmCIF
+ character*500 ctmp(1000),ch(nmax),ent(nmax),mn(nmax)
+ character*500 ch_t,ent_t,mn_t
+ character*20 Aatom(2,90000)
+ character Agroup(2,90000)*6
+ character Ares(2,90000)*3,Aalt(2,90000)
+ character Ains(2,9000),Ach(2,90000),Aent(2,90000)
+ character Cins(2,nmax),Cch(2)
+
+ integer Aatomi(2,90000),Aresi(2,90000),n_cut(2)
+ dimension iform(10),xx(2,90000),yy(2,90000),zz(2,90000),nL(2)
+ccc
- dimension m1(nmax),m2(nmax)
+ dimension m1(nmax),m2(nmax),m12(2,nmax)
dimension xtm1(nmax),ytm1(nmax),ztm1(nmax)
dimension xtm2(nmax),ytm2(nmax),ztm2(nmax)
common/init/invmap_i(nmax)
common/TM/TM,TMmax
common/d8/d8
- common/initial4/mm1(nmax),mm2(nmax)
+ common/initial4/mm(2,nmax)
character*10 aa1,ra1,aa2,ra2,du2
dimension ia1(90000),aa1(90000),ra1(90000),ir1(90000)
@@ -123,10 +138,9 @@
dimension ma1(nmax),ma2(nmax)
dimension nc1(nmax),nc2(nmax)
- 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
-
+ dimension nres1(2,nmax,32:122),nres2(2,nmax,32:122,32:122) !number of atoms
+ character*5 atom1(50) !atom name
+
ccc RMSD:
double precision r_1(3,nmax),r_2(3,nmax),w(nmax)
double precision u(3,3),t(3),rms !armsd is real
@@ -200,7 +214,7 @@ ccc
goto 9999
endif
- version='20170708'
+ version='20190708'
if(fnam.eq.'-v')then
write(*,*)'TM-align Version ',version
goto 9999
@@ -264,218 +278,186 @@ ccc
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)
- 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
+**********decide file format (PDB or mmCIF) ------------------->
+ do j=1,2
+ iform(j)=0 !format of pdb(j)
+ open(unit=10,file=pdb(j),status='old')
+ do while(iform(j).eq.0)
+ read(10,*)s
+ if(s(1:6).eq.'HEADER'.or.s(1:4).eq.'ATOM'.or.
+ & s(1:6).eq.'REMARK')then
+ iform(j)=1 !PDB format
+ elseif(s(1:4).eq.'data'.or.s(1:1).eq.'#'.or.
+ & s(1:5).eq.'loop_')then
+ iform(j)=2 !mmCIF format
endif
enddo
- 1013 continue
+ if(iform(j).eq.0)then
+ write(*,*)'error: file must in PDB or mmCIF format!'
+ endif
close(10)
- nseq1=i
-
- open(unit=10,file=pdb(2),status='old')
+ enddo
+*******^^^^ format is decided ^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+cccccccccRead data from CA file ---------------->
+c we only need to read following (keep first chain, keep only one altLoc):
+c xa(i,j,k)----(x,y,z)
+c mm(2,nmax)---residue order number, for gapless threading
+c ss(2,nmax)---residue name ('GLY') for seq_ID calculation and output
+c seq1(0,nmax),seq2(0,nmax)----single characters ('G', 'N')
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+ do 1001 ic=1,2 !ic=1,2 for file1 and file2
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
+ open(unit=10,file=pdb(ic),status='old')
+ if(iform(ic).eq.1)then !file in PDB format------->
+ do while (.true.) !start do-while ---------->
+ read(10,'(A500)',end=1013) s
+ if(i.gt.0)then
+ if(s(1:3).eq.'TER'.or.s(1:3).eq.'MOD'.
+ & or.s(1:3).eq.'END')then
+ goto 1013 !only read the first chain
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
- ntt=0
- do while (.true.)
- read(10,9001,end=1010) s
- if(i.gt.0.and.s(1:3).eq.'TER')goto 1010
- if(s(1:3).eq.'ATO')then
+ if(s(1:4).eq.'ATOM')then !read 'ATOM' =======>
+ read(s(13:16),*)du4 !will remove space before 'CA'
+ if(du4.eq.'CA')then !read 'CA' ---------->
+ du1=s(27:27) !residue insertion tag
+*******remove repeated altLoc atoms -------->
+ mk=1
+ if(s(17:17).ne.' ')then !with Alternate atom
+ read(s(23:26),*)i8 !res ID
+ if(nres1(ic,i8,ichar(du1)).ge.1)then !since CA, cannot be >1
+ mk=-1 !this residue was already read, skip it
+ endif
+ endif
+c^^^^^^^^^altLoc checked (mk.ne.1 will be skipped) ^^^^^^^^^
+ if(mk.eq.1)then !mk=1 ------------>
+ i=i+1
+ read(s,'(a6,I5,a6,A3,A2,i4,A1,a3,3F8.3)')
+ & du,itmp,du,aanam,du,mm(ic,i),
+ & Cins(ic,i),du,
+ & xa(1,i,ic-1),xa(2,i,ic-1),xa(3,i,ic-1)
***
- du1=s(27:27) !read insertion tag
- mk=1
- 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
+ i8=mm(ic,i)
+ nres1(ic,i8,ichar(du1))=
+ & nres1(ic,i8,ichar(du1))+1 !nres1 only for check altLoc
***
- 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
+ ss(ic,i)=aanam !residue name, 'GLY', for seq_ID
+ if(i.ge.nmax)goto 1013
+ endif !<-----mk=1
+ endif !<---- read 'CA'
+ endif !<==== read 'ATOM'
+ enddo !<----end do-while
+ else !<----mmCIF format,if(iform(1).eq.2)
+ in=0 !number of entries
+ do while (.true.) !start do-while ---------->
+ read(10,'(A500)',end=1013) s
+ if(i.gt.0)then !skip unuseful read to save time
+ if(s(1:1).eq.'#'.or.s(1:5).eq.'loop_'.or.
+ & s(1:6).eq.'HETATM')goto 1013
endif
- 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,9000)du,ma1(i),du,aanam,du,mm1(i),du,
- $ xa(1,i,0),xa(2,i,0),xa(3,i,0)
- ins1(i)=du1
- do j=-1,20
- if(aanam.eq.aa(j))then
- seq1(i)=slc(j)
- goto 21
- endif
- enddo
- seq1(i)=slc(-1)
- 21 continue
- ss1(i)=aanam
- if(i.ge.nmax)goto 1010
+ if(s(1:11).eq.'_atom_site.')then
+ in=in+1
+ read(s,*)du
+ if(du.eq.'_atom_site.label_atom_id')i_atom=in !CA,O, no need
+ if(du.eq.'_atom_site.label_alt_id')i_alt=in !'.',A,B
+ if(du.eq.'_atom_site.label_comp_id')i_res=in !GLY,LEU
+ if(du.eq.'_atom_site.label_asym_id')i_ch=in !A,B
+ if(du.eq.'_atom_site.label_entity_id')i_ent=in !1,2,a,b
+ if(du.eq.'_atom_site.label_seq_id')i_resi=in !1,2,3
+ if(du.eq.'_atom_site.pdbx_PDB_ins_code')i_ins=in !A,?
+ if(du.eq.'_atom_site.Cartn_x')i_x=in !x, 1.234
+ if(du.eq.'_atom_site.Cartn_y')i_y=in !y, 1.234
+ if(du.eq.'_atom_site.Cartn_z')i_z=in !z, 1.234
+ if(du.eq.'_atom_site.pdbx_PDB_model_num')i_mn=in !model number, 1,2,3
endif
- endif
- endif
- enddo
- 1010 continue
- 8999 format(a6,I5,a1,A4,a1,A3,a2,I4,a4,3F8.3)
- 9000 format(a6,I5,a6,A3,A2,i4,A4,3F8.3)
- 9001 format(A100)
- close(10)
- nseq1=i
-c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-
-ccccccccc read data from the second structure file:
- open(unit=10,file=pdb(2),status='old')
- i=0
- na2=0
- ntt=0
- do while (.true.)
- read(10,9001,end=1011) s
- if(i.gt.0.and.s(1:3).eq.'TER')goto 1011
- if(s(1:3).eq.'ATO')then
+ if(s(1:4).eq.'ATOM')then !read 'ATOM' =======>
+ read(s,*)(ctmp(j),j=1,in) !no space before characters
+ if(ctmp(i_atom).eq.'CA')then !read 'CA' ---------->
+ if(i.gt.0)then
+ if(i_mn.gt.1)then !sometimes it may have no i_mn
+ read(ctmp(i_mn),*)mn_t
+ if(mn_t.ne.mn(i)) goto 1013 !only read first model
+ endif
+ read(ctmp(i_ch),*)ch_t
+ if(ch_t.ne.ch(i)) goto 1013 !only read first chain
+ read(ctmp(i_ent),*)ent_t
+ if(ent_t.ne.ent(i)) goto 1013 !only read first entity
+ endif
+
+*******remove repeated altLoc atoms (this does not care inserted residues) -------->
+ mk=1
+ du1=ctmp(i_ins) !read insertion tag
+ if(ctmp(i_alt).ne.'.')then !with Alternate atom
+ read(ctmp(i_resi),*)i8 !res ID
+ if(nres1(ic,i8,ichar(du1)).ge.1)then !since CA, cannot be >1
+ mk=-1 !this residue was already read, skip it
+ endif
+ endif
+c^^^^^^^^^altLoc checked (mk.ne.1 will be skipped) ^^^^^^^^^
+ if(mk.eq.1)then !mk=1 ------------>
+ i=i+1
+ read(ctmp(i_ch),*)ch(i) !for check other chain
+ read(ctmp(i_ent),*)ent(i) !for check other entity
+ read(ctmp(i_mn),*)mn(i) !for check model number
+***
+ read(ctmp(i_x),*)xa(1,i,ic-1)
+ read(ctmp(i_y),*)xa(2,i,ic-1)
+ read(ctmp(i_z),*)xa(3,i,ic-1)
+ read(ctmp(i_resi),*)mm(ic,i) !residue order, 4,5,6
+ read(ctmp(i_ins),*)Cins(ic,i) !residue insertion, A, ?
+ if(Cins(ic,i).eq.'?')Cins(ic,i)=' '
+ ss(ic,i)=ctmp(i_res) !residue name, 'GLY', for seq_ID
+***
+ i8=mm(ic,i)
+ nres1(ic,i8,ichar(du1))=
+ & nres1(ic,i8,ichar(du1))+1 !nres1 only for check altLoc
***
- du1=s(27:27) !read insertion tag
- mk=1
- if(s(17:17).ne.' ')then
- du2=s(13:16)
- read(s(23:26),*)i8
- do i1=1,nres2(i8,ichar(du1))
- if(du2.eq.atom2(i8,i1))then
- mk=-1
+ if(i.ge.nmax)goto 1013
+ endif !<-----mk=1
+ endif !<-----read 'CA'
+ endif !<==== read 'ATOM'
+ enddo !<----end do-while
+ endif !if(iform(ic).eq.2)
+ 1013 continue
+ close(10)
+
+c-------convert 'GLY' to 'G' ------------>
+ if(ic.eq.1)then
+ nseq1=i
+ do i=1,nseq1
+ do j=-1,20
+ if(ss(1,i).eq.aa(j))then
+ seq1(i)=slc(j)
+ goto 121
endif
enddo
- endif
- if(mk.eq.1)then
-***
- 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
- 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)
- ins2(i)=du1
- do j=-1,20
- if(aanam.eq.aa(j))then
- seq2(i)=slc(j)
- goto 22
- endif
- enddo
- seq2(i)=slc(-1)
- 22 continue
- ss2(i)=aanam
- if(i.ge.nmax)goto 1011
- endif
- endif
+ seq1(i)=slc(-1)
+ 121 continue
+ enddo
+ else
+ nseq2=i
+ do i=1,nseq2
+ do j=-1,20
+ if(ss(2,i).eq.aa(j))then
+ seq2(i)=slc(j)
+ goto 122
+ endif
+ enddo
+ seq2(i)=slc(-1)
+ 122 continue
+ enddo
endif
- enddo
- 1011 continue
- close(10)
- nseq2=i
-c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
- 1017 continue
+ 1001 continue !ic=1,2 for file1 and file2
+
+c^^^^^^^^^^^^^read input files completed ^^^^^^^^^^^^^^^^^^^^^^
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')
n_p=0
do while (.true.)
- read(10,9002,end=1012)s1
-c write(*,*)'s1=',trim(s1)
+ read(10,'(A5000)',end=1012)s1
if(s1(1:1).eq.">")then
n_p=n_p+1
sequence(n_p)=''
@@ -484,14 +466,11 @@ c write(*,*)'s1=',trim(s1)
if(n_p.gt.0)then
sequence(n_p)=sequence(n_p)(1:len_trim(sequence(n_p)))
& //s1(1:len_trim(s1))
-c write(*,*)n_p,trim(sequence(n_p))
endif
endif
enddo
1012 continue
close(10)
-c write(*,*)trim(sequence(1))
-c write(*,*)trim(sequence(2))
if(n_p.lt.2)then
write(*,*)'ERROR: FASTA format is wrong, two proteins',
& ' should be included'
@@ -503,7 +482,6 @@ c write(*,*)trim(sequence(2))
& ' run it anyway'
endif
endif
- 9002 format(A5000)
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
* Scale of TM-score in search based on length of smaller protein --------->
@@ -517,7 +495,6 @@ c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
endif
d0_min=d0+0.8 !best for search, this should be changed when calculate real TM-score
if(d0.lt.d0_min)d0=d0_min !min d0 in search=0.968, min d0 in output=0.5
-c write(*,*)'d0 in search=',d0
d00=d0 !for quickly calculate TM-score in searching
if(d00.gt.8)d00=8
if(d00.lt.4.5)d00=4.5
@@ -570,7 +547,7 @@ c write(*,*)'d0 in search=',d0
m1(j)=m1(i) !record alignment
m2(j)=m2(i)
- if(ss1(m1(i)).eq.ss2(m2(i)))then
+ if(ss(1,m1(i)).eq.ss(2,m2(i)))then
n_eq=n_eq+1
endif
endif
@@ -579,7 +556,6 @@ c write(*,*)'d0 in search=',d0
seq_id=float(n_eq)/(n8_al+0.00000001)
call u3b(w,r_1,r_2,n8_al,0,rms,u,t,ier)
rmsd=dsqrt(rms/n8_al)
-c write(*,*)'---------',rms,n8_al,RMSD,u(1,1),t(1)
*^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
*^^^^^^^ alignment is done, all cutoffs were based on shorter chain^^^^^^^^
@@ -651,7 +627,9 @@ c write(*,*)'---------',rms,n8_al,RMSD,u(1,1),t(1)
TMfix=TM8*n8_al/nseq2
endif
-*********for output summary ******************************
+*****************************************
+* screen output of TM-align results *
+*****************************************
write(*,*)
write(*,*)'*****************************************************',
& '*********************'
@@ -755,244 +733,9 @@ c write(*,*)'---------',rms,n8_al,RMSD,u(1,1),t(1)
write(1,*)' enddo'
write(1,*)
close(1)
- 209 format(I2,f18.10,f15.10,f15.10,f15.10)
endif
+ 209 format(I2,f18.10,f15.10,f15.10,f15.10)
-********* output superposition ******************************
- if(m_out.eq.1)then
-c 11, output superimpostion in the aligned regions ------------->
- 1236 format('ATOM ',i5,' CA ',A3,' A',I4,A1,3X,3F8.3)
- 1237 format('ATOM ',i5,' CA ',A3,' B',I4,A1,3X,3F8.3)
- 1238 format('TER')
- 1239 format('CONECT',I5,I5)
- 900 format(A)
- 902 format('select ',I4,':A,',I4,':B')
- 903 format('REMARK TM-align Version ',A8,'')
- 104 format('REMARK Chain 1:',A10,' Size=',I4)
- 105 format('REMARK Chain 2:',A10,' Size=',I4,
- & ' (TM-score is normalized by ',I4,', d0=',f6.2,')')
- 106 format('REMARK Aligned length=',I4,', RMSD=',f6.2,
- & ', TM-score=',f7.5,', ID=',f5.3)
- OPEN(unit=7,file=outname,status='unknown') !TM.sup=pdb1.aln + pdb2.aln
-*** script:
- write(7,900)'load inline'
- write(7,900)'select *A'
- write(7,900)'wireframe .45'
- write(7,900)'select *B'
- write(7,900)'wireframe .20'
- write(7,900)'select all'
- write(7,900)'color white'
- do i=1,n8_al
- dis2=sqrt((xtm1(i)-xtm2(i))**2+
- & (ytm1(i)-ytm2(i))**2+(ztm1(i)-ztm2(i))**2)
- if(dis2.le.d0_out)then
- write(7,902)mm1(m1(i)),mm2(m2(i))
- write(7,900)'color red'
- endif
- enddo
- write(7,900)'select all'
- write(7,900)'exit'
- write(7,903)version
- write(7,104)pdb(1),nseq1
- write(7,105)pdb(2),nseq2,int(anseq),d0
- write(7,106)n8_al,rmsd,TM2,seq_id
-*** chain1:
- do i=1,n8_al
- write(7,1236)m1(i),ss1(m1(i)),mm1(m1(i)),ins1(m1(i)),
- & xtm1(i),ytm1(i),ztm1(i)
- enddo
- write(7,1238) !TER
- do i=2,n8_al
- write(7,1239)m1(i-1),m1(i) !connect atoms
- enddo
-*** chain2:
- do i=1,n8_al
- write(7,1237)5000+m2(i),ss2(m2(i)),mm2(m2(i)),ins2(m2(i)),
- $ xtm2(i),ytm2(i),ztm2(i)
- enddo
- write(7,1238)
- do i=2,n8_al
- write(7,1239)5000+m2(i-1),5000+m2(i)
- enddo
- close(7)
-c 22, output CA-trace of whole chain in 'TM.sup_all' -------->
- outnameall=outname(1:len_trim(outname))//'_all'
- OPEN(unit=8,file=outnameall,status='unknown') !pdb1.aln + pdb2.aln
-*** script:
- write(8,900)'load inline'
- write(8,900)'select *A'
- write(8,900)'wireframe .45'
- write(8,900)'select none'
- write(8,900)'select *B'
- write(8,900)'wireframe .20'
- write(8,900)'color white'
- do i=1,n8_al
- dis2=sqrt((xtm1(i)-xtm2(i))**2+
- & (ytm1(i)-ytm2(i))**2+(ztm1(i)-ztm2(i))**2)
- if(dis2.le.d0_out)then
- write(8,902)mm1(m1(i)),mm2(m2(i)) !select residue
- write(8,900)'color red'
- endif
- enddo
- write(8,900)'select all'
- write(8,900)'exit'
- write(8,903)version
- write(8,104)pdb(1),nseq1
- write(8,105)pdb(2),nseq2,int(anseq),d0
- write(8,106)n8_al,rmsd,TM2,seq_id
-*** chain1:
- do i=1,nseq1
- ax=t(1)+u(1,1)*xa(1,i,0)+u(1,2)*xa(2,i,0)+u(1,3)*xa(3,i,0)
- ay=t(2)+u(2,1)*xa(1,i,0)+u(2,2)*xa(2,i,0)+u(2,3)*xa(3,i,0)
- az=t(3)+u(3,1)*xa(1,i,0)+u(3,2)*xa(2,i,0)+u(3,3)*xa(3,i,0)
- write(8,1236)i,ss1(i),mm1(i),ins1(i),ax,ay,az
- enddo
- write(8,1238) !TER
- do i=2,nseq1
- write(8,1239)i-1,i !CONECT atom numbers
- enddo
-*** chain2:
- do i=1,nseq2
- write(8,1237)5000+i,ss2(i),mm2(i),ins2(i),
- $ xa(1,i,1),xa(2,i,1),xa(3,i,1)
- enddo
- write(8,1238)
- do i=2,nseq2
- write(8,1239)5000+i-1,5000+i
- enddo
- close(8)
-c 33, output full-atomic structure of whole chain in 'TM.sup_atm' -------->
- outnameall=outname(1:len_trim(outname))//'_atm'
- OPEN(unit=8,file=outnameall,status='unknown') !pdb1.aln + pdb2.aln
-*** script:
- write(8,900)'load inline'
- write(8,900)'select *A'
- write(8,900)'color blue'
- write(8,900)'select *B'
- write(8,900)'color red'
- write(8,900)'select all'
- write(8,900)'cartoon'
- write(8,900)'exit'
- write(8,903)version
- write(8,104)pdb(1),nseq1
- write(8,105)pdb(2),nseq2,int(anseq),d0
- write(8,106)n8_al,rmsd,TM2,seq_id
-*** chain1:
- do i=1,na1
- do j=1,n8_al
- if(ir1(i).eq.mm1(m1(j)))then !aligned residues
- if(ains1(i).eq.ins1(m1(j)))then
- 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(8,8888)ia1(i),aa1(i),ra1(i),ir1(i),ains1(i),
- & ax,ay,az
- endif
- endif
- enddo
- enddo
- write(8,1238) !TER
-*** chain2:
- do i=1,na2
- do j=1,n8_al
- if(ir2(i).eq.mm2(m2(j)))then !aligned residues
- if(ains2(i).eq.ins2(m2(j)))then
- write(8,8889)ia2(i),aa2(i),ra2(i),ir2(i),ains2(i),
- & xa2(i),ya2(i),za2(i)
- endif
- endif
- enddo
- enddo
- write(8,1238) !TER
- close(8)
-c 44, output full-atomic structure of whole chain in 'TM.sup_all_atm' -------->
- outnameall=outname(1:len_trim(outname))//'_all_atm'
- OPEN(unit=8,file=outnameall,status='unknown') !pdb1.aln + pdb2.aln
-*** script:
- write(8,900)'load inline'
- write(8,900)'select *A'
- write(8,900)'color blue'
- write(8,900)'select *B'
- write(8,900)'color red'
- write(8,900)'select all'
- write(8,900)'cartoon'
- write(8,900)'exit'
- write(8,903)version
- write(8,104)pdb(1),nseq1
- write(8,105)pdb(2),nseq2,int(anseq),d0
- write(8,106)n8_al,rmsd,TM2,seq_id
-*** 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(8,8888)ia1(i),aa1(i),ra1(i),ir1(i),ains1(i),
- & ax,ay,az
- enddo
- write(8,1238) !TER
-*** chain2:
- do i=1,na2
- write(8,8889)ia2(i),aa2(i),ra2(i),ir2(i),ains2(i),
- & xa2(i),ya2(i),za2(i)
- enddo
- write(8,1238) !TER
- close(8)
-c 55, output full-atomic structure of whole chain in 'TM.sup_all_atm_lig' ---->
- outnameall=outname(1:len_trim(outname))//'_all_atm_lig'
- OPEN(unit=8,file=outnameall,status='unknown') !pdb1.aln + pdb2.aln
-*** script:
- write(8,900)'load inline'
- write(8,900)'select all'
- write(8,900)'cartoon'
- write(8,900)'select *A'
- write(8,900)'color blue'
- write(8,900)'select *B'
- write(8,900)'color red'
- write(8,900)'select ligand'
- write(8,900)'wireframe 0.25'
- write(8,900)'select solvent'
- write(8,900)'spacefill 0.25'
- write(8,900)'select all'
- write(8,900)'exit'
- write(8,903)version
- write(8,104)pdb(1),nseq1
- write(8,105)pdb(2),nseq2,int(anseq),d0
- write(8,106)n8_al,rmsd,TM2,seq_id
-*** chain1:
- open(unit=10,file=pdb(1),status='old')
- do while (.true.)
- read(10,9001,end=1015) s
- if(s(1:6).eq.'ATOM '.or.s(1:6).eq.'HETATM')then
- read(s,8900)du,xxx,yyy,zzz
- xxn=t(1)+u(1,1)*xxx+u(1,2)*yyy+u(1,3)*zzz
- yyn=t(2)+u(2,1)*xxx+u(2,2)*yyy+u(2,3)*zzz
- zzn=t(3)+u(3,1)*xxx+u(3,2)*yyy+u(3,3)*zzz
- write(8,8901)du(1:21),'A',du(23:30),xxn,yyn,zzn
- endif
- enddo
- 1015 continue
- close(10)
- write(8,1238) !TER
-*** chain2:
- open(unit=10,file=pdb(2),status='old')
- do while (.true.)
- read(10,9001,end=1016) s
- if(s(1:6).eq.'ATOM '.or.s(1:6).eq.'HETATM')then
- read(s,8900)du,xxx,yyy,zzz
- write(8,8901)du(1:21),'B',du(23:30),xxx,yyy,zzz
- endif
- enddo
- 1016 continue
- close(10)
- write(8,1238) !TER
- close(8)
- 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)
- 8900 format(a30,3F8.3)
- 8901 format(a21,a1,a8,3F8.3)
- endif
-*^^^^^^^^^^^^^^^^^^ output finished ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-
************ output aligned sequences **************************
ii=0
i1_old=1
@@ -1044,6 +787,369 @@ c 55, output full-atomic structure of whole chain in 'TM.sup_all_atm_lig' --
10 format(10000A1)
write(*,*)
+c^^^^^^^^^^screen output is done ^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+ if(m_out.ne.1)then
+ stop
+ endif
+
+*******************************************************
+* output alignment structure for Rasmal review *
+*******************************************************
+
+c*****************************************************************
+c************* Step-1: read all-atom structures ******************
+c*****************************************************************
+ do 1002 ic=1,2 !ic=1,2 for file1 and file2
+ nL(ic)=0 !number of lines in PDB file
+ n_cut(ic)=0 !end of first chain, decide to show aaa_all_atm_lig
+ open(unit=10,file=pdb(ic),status='old')
+ if(iform(ic).eq.1)then !file in PDB format %%%%%%%%%%------->
+ do while (.true.)
+ read(10,'(A500)',end=1015) s
+ if(n_cut(ic).eq.0.and.nL(ic).gt.0)then !decide n_cut
+ if(s(1:3).eq.'TER'.or.s(1:3).eq.'MOD'.
+ & or.s(1:3).eq.'END')then
+ n_cut(ic)=nL(ic)
+ endif
+ endif
+
+*** skip model_number >=2 ----------->
+ if(s(1:5).eq.'MODEL')then
+ read(s,*)du,model_num
+ if(model_num.gt.1)then
+ do while(.true.)
+ read(10,'(A500)',end=1015) s
+ if(s(1:6).eq.'ENDMDL')goto 1014
+ enddo
+ endif
+ endif
+*^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+ if(s(1:6).eq.'ATOM '.or.s(1:6).eq.'HETATM')then !read ATOM/HETATM ====>
+*******remove repeated altLoc atoms -------->
+ mk=1
+ du1=s(27:27) !read insertion tag
+ du3=s(22:22) !chain ID
+ if(s(17:17).ne.' ')then !with Alternate atom
+ du2=s(13:16) !atom ID
+c read(s(13:16),*)du2
+ read(s(23:26),*)i8 !res ID
+ do i1=1,nres2(ic,i8,ichar(du1),ichar(du3)) !#of atoms for res_insert
+ if(du2.eq.atom1(i1))then !such atom was already read
+ mk=-1
+ endif
+ enddo
+ endif
+c^^^^^^^^^altLoc checked (mk.ne.1) will be skipped) ^^^^^^^^^
+
+ if(mk.eq.1)then !mk=1--------->
+ nL(ic)=nL(ic)+1
+ read(s,8999)Agroup(ic,nL(ic)),Aatomi(ic,nL(ic)),
+ & du,Aatom(ic,nL(ic)),Aalt(ic,nL(ic)),
+ & Ares(ic,nL(ic)),du,Ach(ic,nL(ic)),
+ & Aresi(ic,nL(ic)),Ains(ic,nL(ic)),du,
+ & xx(ic,nL(ic)),yy(ic,nL(ic)),zz(ic,nL(ic))
+c read(Aatom(ic,nL(ic)),*)Aatom(ic,nL(ic)) !remove space before ' CA'
+***
+ i8=Aresi(ic,nL(ic))
+ if(nres2(ic,i8,ichar(du1),ichar(du3)).lt.50)then
+ nres2(ic,i8,ichar(du1),ichar(du3))=
+ & nres2(ic,i8,ichar(du1),ichar(du3))+1 !#of atoms for res_ins
+ atom1(nres2(ic,i8,ichar(du1),ichar(du3)))=
+ & Aatom(ic,nL(ic)) !atom ID
+ endif
+***
+ if(nL(ic).ge.90000)goto 1015
+ endif !<------mk=1
+ endif !<======read ATOM/HETATM
+ 1014 continue !skip model_num >1
+ enddo !do while
+ else !<----mmCIF format,if(iform(1).eq.2) %%%%%%%%%%%%%%%%
+ in=0 !number of entries
+ do while (.true.) !start do-while ---------->
+ read(10,'(A500)',end=1015) s
+ if(nL(ic).gt.0)then !skip unuseful read to save time
+ if(s(1:1).eq.'#'.or.s(1:5).eq.'loop_')goto 1015
+ endif
+ if(s(1:11).eq.'_atom_site.')then
+ in=in+1
+ read(s,*)du
+ if(du.eq.'_atom_site.group_PDB')i_group=in !'ATOM','HETATM'
+ if(du.eq.'_atom_site.id')i_atomi=in !1,2,3,4
+ if(du.eq.'_atom_site.type_symbol')i_type=in !N,C,O
+ if(du.eq.'_atom_site.label_atom_id')i_atom=in !CA,O
+ if(du.eq.'_atom_site.label_alt_id')i_alt=in !'.',A,B
+
+ if(du.eq.'_atom_site.label_comp_id')i_res=in !GLY,LEU
+ if(du.eq.'_atom_site.label_asym_id')i_ch=in !A,B
+ if(du.eq.'_atom_site.label_entity_id')i_ent=in !1,2,a,b
+ if(du.eq.'_atom_site.label_seq_id')i_resi=in !1,2,3
+ if(du.eq.'_atom_site.pdbx_PDB_ins_code')i_ins=in !A,?
+
+ if(du.eq.'_atom_site.Cartn_x')i_x=in !x, 1.234
+ if(du.eq.'_atom_site.Cartn_y')i_y=in !y, 1.234
+ if(du.eq.'_atom_site.Cartn_z')i_z=in !z, 1.234
+ if(du.eq.'_atom_site.pdbx_PDB_model_num')i_mn=in !model number, 1,2,3
+ endif
+
+ if(s(1:4).eq.'ATOM'.or.s(1:6).eq.'HETATM')then !read 'ATOM' =======>
+ read(s,*)(ctmp(j),j=1,in)
+*** skip model_num >1 --------------------------------------->
+ if(i_mn>1)then
+ read(ctmp(i_mn),*)model_num
+ if(model_num.gt.1)goto 1016 !skip model_num >1
+ endif
+*^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+ if(n_cut(ic).eq.0.and.nL(ic).gt.0)then !decide n_cut
+ read(ctmp(i_ch),*)ch_t
+ read(ctmp(i_ent),*)ent_t
+ if(ch_t.ne.Ach(ic,nL(ic)).or.ent_t.ne.
+ & Aent(ic,nL(ic)))then
+ n_cut(ic)=nL(ic)
+ endif
+ endif
+
+*******remove repeated altLoc atoms (this does not care inserted residues) -------->
+ mk=1
+ du1=ctmp(i_ins) !read insertion tag
+ du3=ctmp(i_ch) !chain ID
+ if(ctmp(i_alt).ne.'.')then !with Alternate atom
+ du2=ctmp(i_atom) !atom ID
+ read(ctmp(i_resi),*)i8 !res ID
+ do i1=1,nres2(ic,i8,ichar(du1),ichar(du3)) !#of atoms for res_insert
+ if(du2.eq.atom1(i1))then !such atom was already read
+ mk=-1
+ endif
+ enddo
+ endif
+c^^^^^^^^^altLoc checked (mk.ne.1) will be skipped) ^^^^^^^^^
+
+ if(mk.eq.1)then !mk=1 -------------->
+ nL(ic)=nL(ic)+1
+ read(ctmp(i_group),*)Agroup(ic,nL(ic))
+ read(ctmp(i_atomi),*)Aatomi(ic,nL(ic))
+ read(ctmp(i_atom),*)Aatom(ic,nL(ic))
+********add space before Aatom(ic,nL(ic)) for output ------>
+ if(Aatom(ic,nL(ic))(4:4).eq.' ')then
+ Aatom(ic,nL(ic))=' '//Aatom(ic,nL(ic))
+ endif
+c read(ctmp(i_alt),*)Aalt(ic,nL(ic)) ! not used, because we alway use 1 atom for altLoc
+c if(Aalt(ic,nL(ic)).eq.'.')Aalt(ic,nL(ic))=' '
+
+ read(ctmp(i_res),*)Ares(ic,nL(ic))
+ read(ctmp(i_ch),*)Ach(ic,nL(ic)) !for check other chain
+ read(ctmp(i_ent),*)Aent(ic,nL(ic)) !for check other entity
+ if(ctmp(i_resi)(1:1).eq.'.')ctmp(i_resi)='0' !for HETATM
+ read(ctmp(i_resi),*)Aresi(ic,nL(ic))
+ read(ctmp(i_ins),*)Ains(ic,nL(ic))
+ if(Ains(ic,nL(ic)).eq.'?')Ains(ic,nL(ic))=' '
+***
+ read(ctmp(i_x),*)xx(ic,nL(ic))
+ read(ctmp(i_y),*)yy(ic,nL(ic))
+ read(ctmp(i_z),*)zz(ic,nL(ic))
+***
+ i8=Aresi(ic,nL(ic))
+ if(nres2(ic,i8,ichar(du1),ichar(du3)).lt.50)then
+ nres2(ic,i8,ichar(du1),ichar(du3))=
+ & nres2(ic,i8,ichar(du1),ichar(du3))+1 !#of atoms for res_ins
+ atom1(nres2(ic,i8,ichar(du1),ichar(du3)))=
+ & Aatom(ic,nL(ic)) !atom ID
+ endif
+***
+ if(nL(ic).ge.90000)goto 1015
+ endif !<---- mk=1
+ 1016 continue !skip model_num >1
+ endif !<==== read 'ATOM'
+ enddo !<----end do-while
+ endif !<----mmCIF format,if(iform(1).eq.2) %%%%%%%%%
+ 1015 continue
+ close(10)
+ 1002 continue !ic=1,2 for file1 and file2
+
+ 8999 format(a6,I5,a1,A4,a1,A3,a1,A1,I4,A1,a3,3F8.3)
+
+********** mark aligned residues --------------->
+ do i=1,n8_al
+ m12(1,i)=m1(i)
+ m12(2,i)=m2(i)
+ enddo
+ Cch(1)='A'
+ Cch(2)='B'
+c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+c*********************************************************************
+c************* Step-2: output structures for Rasmol ******************
+c*********************************************************************
+ do 1019 ip=1,5
+ if(ip.eq.1) outnameall=outname(1:len_trim(outname)) !'aaa', aligned CA
+ if(ip.eq.2) outnameall=outname(1:len_trim(outname))//'_all' !'aaa_all', all CA
+ if(ip.eq.3) outnameall=outname(1:len_trim(outname))//'_atm' !'aaa_atm', all aligned atoms
+ if(ip.eq.4) outnameall=outname(1:len_trim(outname))
+ & //'_all_atm' !'aaa_all_atm', all atoms
+ if(ip.eq.5) outnameall=outname(1:len_trim(outname))
+ & //'_all_atm_lig' !'aaa_all_atm_lig', all atom + ligand
+
+ OPEN(unit=10,file=outnameall,status='unknown')
+*** script:
+ if(ip.eq.1.or.ip.eq.2)then !'aaa','aaa_all'
+ write(10,'(A)')'load inline'
+ write(10,'(A)')'select *A'
+ write(10,'(A)')'wireframe .45'
+ write(10,'(A)')'select *B'
+ write(10,'(A)')'wireframe .20'
+ write(10,'(A)')'select all'
+ write(10,'(A)')'color white'
+ do i=1,n8_al
+ dis2=sqrt((xtm1(i)-xtm2(i))**2+
+ & (ytm1(i)-ytm2(i))**2+(ztm1(i)-ztm2(i))**2)
+ if(dis2.le.d0_out)then
+ write(10,902)mm(1,m1(i)),mm(2,m2(i)) !select residue
+ write(10,'(A)')'color red'
+ endif
+ enddo
+ write(10,'(A)')'select all'
+ elseif(ip.eq.3.or.ip.eq.4)then !'aaa_atm', 'aaa_all_atm'
+ write(10,'(A)')'load inline'
+ write(10,'(A)')'select *A'
+ write(10,'(A)')'color blue'
+ write(10,'(A)')'select *B'
+ write(10,'(A)')'color red'
+ write(10,'(A)')'select all'
+ write(10,'(A)')'cartoon'
+ else !'aaa_all_atm_lig
+ write(10,'(A)')'load inline'
+ write(10,'(A)')'select all'
+ write(10,'(A)')'cartoon'
+ write(10,'(A)')'select *A'
+ write(10,'(A)')'color blue'
+ write(10,'(A)')'select *B'
+ write(10,'(A)')'color red'
+ write(10,'(A)')'select ligand'
+ write(10,'(A)')'wireframe 0.25'
+ write(10,'(A)')'select solvent'
+ write(10,'(A)')'spacefill 0.25'
+ write(10,'(A)')'select all'
+ endif
+ write(10,'(A)')'exit'
+ write(10,903)version
+ write(10,104)pdb(1),nseq1
+ write(10,105)pdb(2),nseq2,int(anseq),d0
+ write(10,106)n8_al,rmsd,TM2,seq_id
+
+c------------output coordinates ----------->
+ do ic=1,2
+ if(ic.eq.1)then !na=1-5000
+ na=0 !for CONECT
+ else !na=5001-10000
+ na=5000 !for CONECT
+ endif
+c write(*,*)'ip=',ip,'ic=',ic,nL(ic)
+
+ if(n_cut(ic).eq.0)then
+ n_cut(ic)=nL(ic) !there is no cut
+ endif
+ do i=1,nL(ic)
+c---------- decide if we should output this line ------------>
+ mk=0
+ if(ip.eq.1)then !'aaa', aligned CA
+ if(i.le.n_cut(ic))then
+ if(Aatom(ic,i).eq.' CA')then
+ do j=1,n8_al
+ if(Aresi(ic,i).eq.mm(ic,m12(ic,j)))then
+ if(Ains(ic,i).eq.Cins(ic,m12(ic,j)))then
+ mk=1
+ goto 120 !each line output once
+ endif
+ endif
+ enddo
+ endif
+ endif
+ elseif(ip.eq.2)then !'aaa_all', all CA
+ if(i.le.n_cut(ic))then
+ if(Aatom(ic,i).eq.' CA')then
+ mk=1
+ endif
+ endif
+ elseif(ip.eq.3)then !'aaa_atm', all aligned atoms
+ if(i.le.n_cut(ic))then
+ do j=1,n8_al
+ if(Aresi(ic,i).eq.mm(ic,m12(ic,j)))then
+ if(Ains(ic,i).eq.Cins(ic,m12(ic,j)))then
+ mk=1
+ goto 120 !each line output once
+ endif
+ endif
+ enddo
+ endif
+ elseif(ip.eq.4)then !'aaa_all_atm', all atoms
+ if(i.le.n_cut(ic))then
+ mk=1
+ endif
+ elseif(ip.eq.5)then !'aaa_all_atm_lig', all atoms and ligands
+ mk=1
+ endif
+ 120 continue
+
+c^^^^^^^^^^
+c write(*,*)ip,ic,i,'mk=',mk,n_cut(ic)
+ if(mk.eq.1)then !mk=1----------->
+ if(ic.eq.1)then
+ ax=t(1)+u(1,1)*xx(ic,i)+u(1,2)*yy(ic,i)+
+ & u(1,3)*zz(ic,i)
+ ay=t(2)+u(2,1)*xx(ic,i)+u(2,2)*yy(ic,i)+
+ & u(2,3)*zz(ic,i)
+ az=t(3)+u(3,1)*xx(ic,i)+u(3,2)*yy(ic,i)+
+ & u(3,3)*zz(ic,i)
+ else
+ ax=xx(ic,i)
+ ay=yy(ic,i)
+ az=zz(ic,i)
+ endif
+
+c--------------output (x,y,z) --------------------------------->
+ na=na+1 !number of atoms for CONECT
+ if(ip.eq.1.or.ip.eq.2)then !need CONECT
+ write(10,1236)Agroup(ic,i),na,'',
+ & Aatom(ic,i),'',Ares(ic,i),'',Cch(ic),
+ & Aresi(ic,i),Ains(ic,i),'',ax,ay,az
+ else !all atoms
+ write(10,1236)Agroup(ic,i),Aatomi(ic,i),'',
+ & Aatom(ic,i),'',Ares(ic,i),'',Cch(ic),
+ & Aresi(ic,i),Ains(ic,i),'',ax,ay,az
+ endif
+ endif !<----mk=1
+ enddo !do i=1,nL(ic)
+
+ write(10,'(A)')'TER' !TER
+ if(ip.eq.1.or.ip.eq.2)then
+ if(ic.eq.1)then
+ do i=2,na
+ write(10,'(A6,I5,I5)')'CONECT',i-1,i !CONECT atom numbers
+ enddo
+ else
+ do i=5002,na
+ write(10,'(A6,I5,I5)')'CONECT',i-1,i !CONECT atom numbers
+ enddo
+ endif
+ endif
+ enddo !do ic=1,2
+ close(10)
+ 1019 continue !ip=1,5
+
+ 902 format('select ',I4,':A,',I4,':B')
+ 903 format('REMARK TM-align Version ',A8,'')
+ 104 format('REMARK Chain 1:',A10,' Size=',I4)
+ 105 format('REMARK Chain 2:',A10,' Size=',I4,
+ & ' (TM-score is normalized by ',I4,', d0=',f6.2,')')
+ 106 format('REMARK Aligned length=',I4,', RMSD=',f6.2,
+ & ', TM-score=',f7.5,', ID=',f5.3)
+ 1236 format(A6,I5,a1,A4,a1,A3,a1,A1,I4,A1,a3,3F8.3)
+
+*^^^^^^^^^^^^^^^^^^ output finished ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
9999 END
@@ -1573,23 +1679,16 @@ c 1->coil, 2->helix, 3->turn, 4->strand
common/zscore/zrms,n_al,rmsd_al
common/TM/TM,TMmax
common/init/invmap_i(nmax)
- common/initial4/mm1(nmax),mm2(nmax)
+ common/initial4/mm(2,nmax)
logical contin
dimension ifr2(2,nmax,nmax),Lfr2(2,nmax),i_fr2(2)
dimension ifr(nmax)
- dimension mm(2,nmax)
fra_min=4 !>=4,minimum fragment for search
fra_min1=fra_min-1 !cutoff for shift, save time
dcu0=4.25
- do i=1,nseq1
- mm(1,i)=mm1(i)
- enddo
- do i=1,nseq2
- mm(2,i)=mm2(i)
- enddo
GL_max=0
c do k=1,2 !k=1, fragment from protein1; k=2, fragment from protein2
@@ -1766,7 +1865,6 @@ c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
common/d0/d0,anseq
common/d0min/d0_min
common/init/invmap_i(nmax)
- common/initial4/mm1(nmax),mm2(nmax)
common/sec/isec(nmax),jsec(nmax)
double precision r_1(3,nmax),r_2(3,nmax),w(nmax)
double precision u(3,3),t(3),rms
@@ -2904,7 +3002,7 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
* 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'
+* For example, '1se9A.pdb' and '2edpA.pdb' (with score(i,j)=0 or 1 in SS)
********************************************************************
SUBROUTINE DP(NSEQ1,NSEQ2)
PARAMETER(nmax=5000)
=====================================
TMscore.f
=====================================
@@ -38,13 +38,16 @@
* complex structure comparisons, where multiple-chains are
* merged into a single chain. Chain ID is now included in
* the output files.
+* 2019/07/08: Enabled TM-score to support both PDB and mmCIF formats,
+* and updated structure reading which makes program faster.
*************************************************************************
c 1 2 3 4 5 6 7 !
c 3456789012345678901234567890123456789012345678901234567890123456789012345678
-
+
program TMscore
- PARAMETER(nmax=5000)
+ PARAMETER(nmax=5000) !maximum number of residues
+ PARAMETER(namax=50000) !maximum of atoms
common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax)
common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
@@ -54,13 +57,13 @@ c 3456789012345678901234567890123456789012345678901234567890123456789012345678
dimension k_ali(nmax),k_ali0(nmax)
character*500 fnam,pdb(100),outname
- character*3 aa(-1:20),seqA(nmax),seqB(nmax)
+ character*3 aa(-1:20),seqA(nmax),seqB(nmax),aanam,ss(2,nmax)
character*500 s,du
- character*1 chA(nmax),chB(nmax)
- character*1 chain1(90000),chain2(90000)
- character seq1A(nmax),seq1B(nmax),ali(nmax),du1
+ character*1 chA(nmax),chB(nmax),ch(nmax)
+ character*1 chain1(namax),chain2(namax)
+ character seq1A(nmax),seq1B(nmax),ali(nmax),du1,du3,du4*2
character sequenceA(nmax),sequenceB(nmax),sequenceM(nmax)
- character ins1(nmax),ins2(nmax),ains1(90000),ains2(90000)
+ character ins1(nmax),ins2(nmax),ains1(namax),ains2(namax)
dimension L_ini(100),iq(nmax)
common/scores/score,score_maxsub,score_fix,score10
@@ -68,22 +71,35 @@ c 3456789012345678901234567890123456789012345678901234567890123456789012345678
double precision score,score_max,score_fix,score_fix_max
double precision score_maxsub,score10
dimension xa(nmax),ya(nmax),za(nmax)
+ dimension x2(2,nmax),y2(2,nmax),z2(2,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)
+ character*10 aa1,ra1,aa2,ra2,du2
+ integer mm(2,nmax)
- 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
+ dimension nres1(2,nmax,32:122),nres2(2,nmax,32:122,32:122) !number of atoms
+ character*5 atom1(50) !atom name
+ dimension m12(2,nmax)
+
+ccc mmCIF
+ character*500 ctmp(1000),ent(nmax),mn(nmax)
+ character*500 ch_t,ent_t,mn_t
+ character*20 Aatom(2,namax)
+ character Agroup(2,namax)*6
+ character Ares(2,namax)*3,Aalt(2,namax)
+ character Ains(2,9000),Ach(2,namax),Aent(2,namax)
+ character Cins(2,nmax)
+
+ integer Aatomi(2,namax),Aresi(2,namax)
+ dimension iform(10),xx(2,namax),yy(2,namax),zz(2,namax),nL(2)
+ccc
ccc RMSD:
double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)
double precision u(3,3),t(3),rms,drms !armsd is real
data w /nmax*1.0/
ccc
+
+ double precision u2(3,3),t2(3) !for output
data aa/ 'BCK','GLY','ALA','SER','CYS',
& 'VAL','THR','ILE','PRO','MET',
@@ -111,17 +127,19 @@ ccc
write(*,*)
write(*,*)'2. Run TM-score to compare two complex structures',
& ' with multiple chains'
+ write(*,*)' (Compare all chains with the same chain',
+ & ' identifier)'
write(*,*)' >TMscore -c model native'
write(*,*)
- write(*,*)'2. TM-score normalized with an assigned scale d0',
+ write(*,*)'3. TM-score normalized with an assigned scale d0',
& ' e.g. 5 A:'
write(*,*)' >TMscore model native -d 5'
write(*,*)
- write(*,*)'3. TM-score normalized by a specific length, ',
+ write(*,*)'4. 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(*,*)'5. 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:'
@@ -164,141 +182,226 @@ ccc
pdb(j)=fnam
endif
if(i.lt.narg)goto 115
+
+**********decide file format (PDB or mmCIF) ------------------->
+ do j=1,2
+ iform(j)=0 !format of pdb(j)
+ open(unit=10,file=pdb(j),status='old')
+ do while(iform(j).eq.0)
+ read(10,*)s
+ if(s(1:6).eq.'HEADER'.or.s(1:4).eq.'ATOM'.or.
+ & s(1:6).eq.'REMARK')then
+ iform(j)=1 !PDB format
+ elseif(s(1:4).eq.'data'.or.s(1:1).eq.'#'.or.
+ & s(1:5).eq.'loop_')then
+ iform(j)=2 !mmCIF format
+ endif
+ enddo
+ if(iform(j).eq.0)then
+ write(*,*)'error: file must in PDB or mmCIF format!'
+ endif
+ close(10)
+ enddo
+*******^^^^ format is decided ^^^^^^^^^^^^^^^^^^^^^^^^^^
-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(m_complex.eq.0)then
- if(s(1:3).eq.'TER') goto 102
- endif
- if(s(1:4).eq.'ATOM')then
+c write(*,*)'m_complex=',m_complex
+
+cccccccccRead data from CA file ---------------->
+c we only need to read following (keep first chain, keep only one altLoc):
+c x,y,z2(ic,i)----(x,y,z)
+c mm(2,nmax)---residue order number, for gapless threading
+c ss(2,nmax)---residue name ('GLY') for seq_ID calculation and output
+c Cins(ic,i)---insertion
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+ do 1001 ic=1,2 !ic=1,2 for file1 and file2
+ i=0
+ open(unit=10,file=pdb(ic),status='old')
+ if(iform(ic).eq.1)then !file in PDB format------->
+ do while (.true.) !start do-while ---------->
+ read(10,'(A500)',end=1013) s
+ if(i.gt.0)then
+ if(m_complex.eq.0)then
+ if(s(1:3).eq.'TER'.or.s(1:3).eq.'MOD'.
+ & or.s(1:3).eq.'END')then
+ goto 1013 !only read the first chain
+ endif
+ else
+*** skip model_number >=2 ----------->
+ if(s(1:5).eq.'MODEL')then
+ read(s,*)du,model_num
+ if(model_num.gt.1)then
+ do while(.true.)
+ read(10,'(A500)',end=1013) s
+ if(s(1:6).eq.'ENDMDL')goto 1014
+ enddo
+ endif
+ endif
+*^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+ endif
+ endif
+ if(s(1:4).eq.'ATOM')then !read 'ATOM' =======>
+ read(s(13:16),*)du4 !will remove space before 'CA'
+ if(du4.eq.'CA')then !read 'CA' ---------->
+ du1=s(27:27) !residue insertion tag
+*******remove repeated altLoc atoms (this does not care inserted residues) -------->
+ mk=1
+ if(s(17:17).ne.' ')then !with Alternate atom
+ read(s(23:26),*)i8 !res ID
+ if(nres1(ic,i8,ichar(du1)).ge.1)then !since CA, cannot be >1
+ mk=-1 !this residue was already read, skip it
+ endif
+ endif
+c^^^^^^^^^altLoc checked (mk.ne.1 will be skipped) ^^^^^^^^^
+ if(mk.eq.1)then !mk=1 ------------>
+ i=i+1
+ read(s,'(a6,I5,a6,A3,a1,A1,I4,A1,a3,3F8.3)') !variable include space, different from read(*,*)
+ & du,itmp,du,aanam,du,ch(i),mm(ic,i),
+ & Cins(ic,i),du,x2(ic,i),y2(ic,i),z2(ic,i)
***
- 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
+ i8=mm(ic,i)
+ nres1(ic,i8,ichar(du1))=
+ & nres1(ic,i8,ichar(du1))+1 !nres1 only for check altLoc
+***
+ ss(ic,i)=aanam !residue name, 'GLY', for seq_ID
+ if(i.ge.nmax)goto 1013
+ endif !<-----mk=1
+ endif !<---- read 'CA'
+ endif !<==== read 'ATOM'
+ 1014 continue
+ enddo !<----end do-while
+ else !<----mmCIF format,if(iform(1).eq.2)
+ in=0 !number of entries
+ do while (.true.) !start do-while ---------->
+ read(10,'(A500)',end=1013) s
+ if(i.gt.0)then !skip unuseful read to save time, suppose '#' appear before complex completed
+ if(s(1:1).eq.'#'.or.s(1:5).eq.'loop_'.or.
+ & s(1:6).eq.'HETATM')goto 1013
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
+ if(s(1:11).eq.'_atom_site.')then
+ in=in+1
+ read(s,*)du
+ if(du.eq.'_atom_site.label_atom_id')i_atom=in !CA,O, no need
+ if(du.eq.'_atom_site.label_alt_id')i_alt=in !'.',A,B
+ if(du.eq.'_atom_site.label_comp_id')i_res=in !GLY,LEU
+ if(du.eq.'_atom_site.label_asym_id')i_ch=in !A,B
+ if(du.eq.'_atom_site.label_entity_id')i_ent=in !1,2,a,b
+ if(du.eq.'_atom_site.label_seq_id')i_resi=in !1,2,3
+ if(du.eq.'_atom_site.pdbx_PDB_ins_code')i_ins=in !A,?
+ if(du.eq.'_atom_site.Cartn_x')i_x=in !x, 1.234
+ if(du.eq.'_atom_site.Cartn_y')i_y=in !y, 1.234
+ if(du.eq.'_atom_site.Cartn_z')i_z=in !z, 1.234
+ if(du.eq.'_atom_site.pdbx_PDB_model_num')i_mn=in !model number, 1,2,3
+ endif
+ if(s(1:4).eq.'ATOM')then !read 'ATOM' =======>
+ read(s,*)(ctmp(j),j=1,in) !no space before characters
+ if(ctmp(i_atom).eq.'CA')then !read 'CA' ---------->
+ if(i.gt.0)then
+ if(m_complex.eq.0)then
+ if(i_mn.gt.1)then !sometimes it may have no i_mn
+ read(ctmp(i_mn),*)mn_t
+ if(mn_t.ne.mn(i)) goto 1013 !only read first model
+ endif
+ read(ctmp(i_ch),*)ch_t
+ if(ch_t.ne.ch(i)) goto 1013 !only read first chain
+ read(ctmp(i_ent),*)ent_t
+ if(ent_t.ne.ent(i)) goto 1013 !only read first entity
+ else
+ if(i_mn.gt.1)then !sometimes it may have no i_mn
+ read(ctmp(i_mn),*)mn_t
+ if(mn_t.ne.mn(i)) goto 1016 !only read first model
+ endif
+ endif
+ endif
+
+*******remove repeated altLoc atoms (this does not care inserted residues) -------->
+ mk=1
+ du1=ctmp(i_ins) !read insertion tag
+
+ if(ctmp(i_alt).ne.'.')then !with Alternate atom
+ read(ctmp(i_resi),*)i8 !res ID
+ if(nres1(ic,i8,ichar(du1)).ge.1)then !since CA, cannot be >1
+ mk=-1 !this residue was already read, skip it
+ endif
+ endif
+c^^^^^^^^^altLoc checked (mk.ne.1 will be skipped) ^^^^^^^^^
+ if(mk.eq.1)then !mk=1 ------------>
+ i=i+1
+ read(ctmp(i_ch),*)ch(i) !for check other chain
+ read(ctmp(i_ent),*)ent(i) !for check other entity
+ read(ctmp(i_mn),*)mn(i) !for check model number
***
- 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
+ read(ctmp(i_x),*)x2(ic,i)
+ read(ctmp(i_y),*)y2(ic,i)
+ read(ctmp(i_z),*)z2(ic,i)
+ read(ctmp(i_resi),*)mm(ic,i) !residue order, 4,5,6
+ read(ctmp(i_ins),*)Cins(ic,i) !residue insertion, A, ?
+ if(Cins(ic,i).eq.'?')Cins(ic,i)=' '
+ ss(ic,i)=ctmp(i_res) !residue name, 'GLY', for seq_ID
+***
+ i8=mm(ic,i)
+ nres1(ic,i8,ichar(du1))=
+ & nres1(ic,i8,ichar(du1))+1 !nres1 only for check altLoc
+***
+ if(i.ge.nmax)goto 1013
+ endif !<-----mk=1
+ endif !<-----read 'CA'
+ endif !<==== read 'ATOM'
+ 1016 continue !skip model_num >1
+ enddo !<----end do-while
+ endif !if(iform(ic).eq.2)
+ 1013 continue
+ close(10)
+
+c-------convert 'GLY' to 'G', etc ------------>
+ if(ic.eq.1)then
+ nseqA=i
+ do i=1,nseqA
+ chA(i)=ch(i)
+ nresA(i)=mm(ic,i)
+ seqA(i)=ss(ic,i) !GLY
+ ins1(i)=Cins(1,i)
+ xa(i)=x2(1,i)
+ ya(i)=y2(1,i)
+ za(i)=z2(1,i)
do j=-1,20
- if(seqA(i).eq.aa(j))then
- seq1A(i)=slc(j)
- goto 21
+ if(ss(1,i).eq.aa(j))then
+ seq1A(i)=slc(j) !G
+ goto 121
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,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(m_complex.eq.0)then
- if(s(1:3).eq.'TER') goto 202
- endif
- 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
- endif
+ 121 continue
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 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
+ else
+ nseqB=i
+ do i=1,nseqB
+ chB(i)=ch(i)
+ nresB(i)=mm(ic,i)
+ seqB(i)=ss(ic,i) !'GLY'
+ ins2(i)=Cins(2,i)
+ xb(i)=x2(2,i)
+ yb(i)=y2(2,i)
+ zb(i)=z2(2,i)
do j=-1,20
- if(seqB(i).eq.aa(j))then
- seq1B(i)=slc(j)
- goto 22
+ if(ss(2,i).eq.aa(j))then
+ seq1B(i)=slc(j) !'G'
+ goto 122
endif
enddo
seq1B(i)=slc(-1)
- 22 continue
- if(i.ge.nmax)goto 202
- endif
+ 122 continue
+ enddo
endif
- endif
- goto 201
- 202 continue
- 203 format(A17,A3,A1,A1,i4,A4,3F8.3)
- 204 format(A100)
- close(10)
- nseqB=i
-c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-
+ 1001 continue !ic=1,2 for file1 and file2
+c^^^^^^^^^^^^^read input files completed ^^^^^^^^^^^^^^^^^^^^^^
+
+c do i=1,nseqA
+c write(*,*)i,xa(i),seqA(i)
+c enddo
+c do i=1,nseqB
+c write(*,*)i,xb(i),seqB(i)
+c enddo
+
******************************************************************
* pickup the aligned residues:
******************************************************************
@@ -339,6 +442,10 @@ c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
goto 9999
endif
+c do i=1,n_ali
+c write(*,*)i,iA(i),iB(i)
+c enddo
+
******************************************************************
* check the residue serial numbers ------------->
if(m_complex.eq.0)then
@@ -635,82 +742,16 @@ c enddo
d=d_output !for output
call score_fun() !give i_ali(i), score_max=score now
-*** 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 atomno= ',I5)
- 902 format('select atomno= ',I5)
- write(7,900)'load inline'
- write(7,900)'select atomno<50001'
- write(7,900)'wireframe .45'
- write(7,900)'select none'
- write(7,900)'select atomno>50000'
- write(7,900)'wireframe .15'
- write(7,900)'color white' !all above atoms are white
- do i=1,n_cut
- 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'
- write(7,900)'exit'
- write(7,514)rmsd_ali
- 514 format('REMARK RMSD of the common residues=',F8.3)
- write(7,515)score_max,d0
- 515 format('REMARK TM-score=',f6.4,' (d0=',f5.2,')')
- do i=1,nseqA
- 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)i-1,i
- enddo
- do i=1,nseqB
- 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)50000+i-1,50000+i
- enddo
- 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
+******* for output ---------->
+ if(m_out.eq.1)then
+ do i=1,3
+ t2(i)=t(i)
+ do j=1,3
+ u2(i,j)=u(i,j)
+ enddo
+ enddo
+ endif
+*******************************
*** record aligned residues by i=[1,nseqA], for sequenceM()------------>
do i=1,nseqA
@@ -894,7 +935,274 @@ ccc RMSD (d<5.0)-------->
602 format(2000I1)
write(*,*)
-c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+c^^^^^^^^^^screen output is done ^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+ if(m_out.ne.1)then
+ stop
+ endif
+
+*******************************************************
+* output alignment structure for Rasmal review *
+*******************************************************
+
+c*****************************************************************
+c************* Step-1: read all-atom structures ******************
+c*****************************************************************
+ do 1002 ic=1,2 !ic=1,2 for file1 and file2
+ nL(ic)=0 !number of lines in PDB file
+ open(unit=10,file=pdb(ic),status='old')
+ if(iform(ic).eq.1)then !file in PDB format %%%%%%%%%%------->
+ do while (.true.)
+ read(10,'(A500)',end=1015) s
+ if(nL(ic).gt.0)then !decide n_cut
+ if(m_complex.eq.0)then
+ if(s(1:3).eq.'TER'.or.s(1:3).eq.'MOD'.
+ & or.s(1:3).eq.'END')then
+ goto 1015 !no ligand allowed
+ endif
+ else
+*** skip model_number >=2 ----------->
+ if(s(1:5).eq.'MODEL')then
+ read(s,*)du,model_num
+ if(model_num.gt.1)then
+ do while(.true.)
+ read(10,'(A500)',end=1015) s
+ if(s(1:6).eq.'ENDMDL')goto 1017
+ enddo
+ endif
+ endif
+*^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+ endif
+ endif
+ if(s(1:6).eq.'ATOM ')then !read ATOM ====>
+*******remove repeated altLoc atoms -------->
+ mk=1
+ du1=s(27:27) !read insertion tag
+ du3=s(22:22) !chain ID
+ 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,nres2(ic,i8,ichar(du1),ichar(du3)) !#of atoms for res_insert
+ if(du2.eq.atom1(i1))then !such atom was already read
+ mk=-1
+ endif
+ enddo
+ endif
+c^^^^^^^^^altLoc checked (mk.ne.1) will be skipped) ^^^^^^^^^
+
+ if(mk.eq.1)then !mk=1--------->
+ nL(ic)=nL(ic)+1
+ read(s,8999)Agroup(ic,nL(ic)),Aatomi(ic,nL(ic)),
+ & du,Aatom(ic,nL(ic)),Aalt(ic,nL(ic)),
+ & Ares(ic,nL(ic)),du,Ach(ic,nL(ic)),
+ & Aresi(ic,nL(ic)),Ains(ic,nL(ic)),du,
+ & xx(ic,nL(ic)),yy(ic,nL(ic)),zz(ic,nL(ic))
+***
+ i8=Aresi(ic,nL(ic))
+ if(nres2(ic,i8,ichar(du1),ichar(du3)).lt.50)then
+ nres2(ic,i8,ichar(du1),ichar(du3))=
+ & nres2(ic,i8,ichar(du1),ichar(du3))+1 !#of atoms for res_ins
+ atom1(nres2(ic,i8,ichar(du1),ichar(du3)))=
+ & Aatom(ic,nL(ic)) !atom ID
+ endif
+***
+ if(nL(ic).ge.namax)goto 1015
+ endif !<------mk=1
+ endif !<======read ATOM
+ 1017 continue
+ enddo !<--- do while
+ else !<----mmCIF format,if(iform(1).eq.2) %%%%%%%%%%%%%%%%
+ in=0 !number of entries
+ do while (.true.) !start do-while ---------->
+ read(10,'(A500)',end=1015) s
+ if(nL(ic).gt.0)then !skip unuseful read to save time
+ if(s(1:1).eq.'#'.or.s(1:5).eq.'loop_')goto 1015
+ endif
+ if(s(1:11).eq.'_atom_site.')then
+ in=in+1
+ read(s,*)du
+ if(du.eq.'_atom_site.group_PDB')i_group=in !'ATOM','HETATM'
+ if(du.eq.'_atom_site.id')i_atomi=in !1,2,3,4
+ if(du.eq.'_atom_site.type_symbol')i_type=in !N,C,O
+ if(du.eq.'_atom_site.label_atom_id')i_atom=in !CA,O
+ if(du.eq.'_atom_site.label_alt_id')i_alt=in !'.',A,B
+
+ if(du.eq.'_atom_site.label_comp_id')i_res=in !GLY,LEU
+ if(du.eq.'_atom_site.label_asym_id')i_ch=in !A,B
+ if(du.eq.'_atom_site.label_entity_id')i_ent=in !1,2,a,b
+ if(du.eq.'_atom_site.label_seq_id')i_resi=in !1,2,3
+ if(du.eq.'_atom_site.pdbx_PDB_ins_code')i_ins=in !A,?
+
+ if(du.eq.'_atom_site.Cartn_x')i_x=in !x, 1.234
+ if(du.eq.'_atom_site.Cartn_y')i_y=in !y, 1.234
+ if(du.eq.'_atom_site.Cartn_z')i_z=in !z, 1.234
+ if(du.eq.'_atom_site.pdbx_PDB_model_num')i_mn=in !model number, 1,2,3
+ endif
+ if(s(1:4).eq.'ATOM')then !read 'ATOM' =======>
+ read(s,*)(ctmp(j),j=1,in)
+ if(nL(ic).gt.0)then
+ if(m_complex.eq.0)then !monomer
+ if(i_mn.gt.1)then !sometimes it may have no i_mn
+ read(ctmp(i_mn),*)mn_t
+ if(mn_t.ne.mn(i)) goto 1015 !only read first model
+ endif
+ read(ctmp(i_ch),*)ch_t
+ if(ch_t.ne.Ach(ic,nL(ic))) goto 1015
+ read(ctmp(i_ent),*)ent_t
+ if(ent_t.ne.Aent(ic,nL(ic))) goto 1015
+ else
+ if(i_mn.gt.1)then !sometimes it may have no i_mn
+ read(ctmp(i_mn),*)mn_t
+ if(mn_t.ne.mn(i)) goto 1018 !only read first model
+ endif
+ endif
+ endif
+
+*******remove repeated altLoc atoms -------->
+ mk=1
+ du1=ctmp(i_ins) !read insertion tag
+ du3=ctmp(i_ch) !chain ID
+ if(ctmp(i_alt).ne.'.')then !with Alternate atom
+ du2=ctmp(i_atom) !atom ID
+ read(ctmp(i_resi),*)i8 !res ID
+ do i1=1,nres2(ic,i8,ichar(du1),ichar(du3)) !#of atoms for res_insert
+ if(du2.eq.atom1(i1))then !such atom was already read
+ mk=-1
+ endif
+ enddo
+ endif
+c^^^^^^^^^altLoc checked (mk.ne.1) will be skipped) ^^^^^^^^^
+
+ if(mk.eq.1)then !mk=1 -------------->
+ nL(ic)=nL(ic)+1
+ read(ctmp(i_group),*)Agroup(ic,nL(ic))
+ read(ctmp(i_atomi),*)Aatomi(ic,nL(ic))
+ read(ctmp(i_atom),*)Aatom(ic,nL(ic))
+********add space before Aatom(ic,nL(ic)) for output ------>
+ if(Aatom(ic,nL(ic))(4:4).eq.' ')then
+ Aatom(ic,nL(ic))=' '//Aatom(ic,nL(ic))
+ endif
+c read(ctmp(i_alt),*)Aalt(ic,nL(ic)) ! not used, because we alway use 1 atom for altLoc
+c if(Aalt(ic,nL(ic)).eq.'.')Aalt(ic,nL(ic))=' '
+
+ read(ctmp(i_res),*)Ares(ic,nL(ic))
+ read(ctmp(i_ch),*)Ach(ic,nL(ic)) !for check other chain
+ read(ctmp(i_ent),*)Aent(ic,nL(ic)) !for check other entity
+ read(ctmp(i_mn),*)mn(i) !for check model number
+ read(ctmp(i_resi),*)Aresi(ic,nL(ic))
+ read(ctmp(i_ins),*)Ains(ic,nL(ic))
+ if(Ains(ic,nL(ic)).eq.'?')Ains(ic,nL(ic))=' '
+***
+ read(ctmp(i_x),*)xx(ic,nL(ic))
+ read(ctmp(i_y),*)yy(ic,nL(ic))
+ read(ctmp(i_z),*)zz(ic,nL(ic))
+
+***
+ i8=Aresi(ic,nL(ic))
+ if(nres2(ic,i8,ichar(du1),ichar(du3)).lt.50)then
+ nres2(ic,i8,ichar(du1),ichar(du3))=
+ & nres2(ic,i8,ichar(du1),ichar(du3))+1 !#of atoms for res_ins
+ atom1(nres2(ic,i8,ichar(du1),ichar(du3)))=
+ & Aatom(ic,nL(ic)) !atom ID
+ endif
+***
+ if(nL(ic).ge.namax)goto 1015
+ endif !<---- mk=1
+ endif !<==== read 'ATOM'
+ 1018 continue
+ enddo !<----end do-while
+ endif !<----mmCIF format,if(iform(1).eq.2) %%%%%%%%%
+ 1015 continue
+ close(10)
+ 1002 continue !ic=1,2 for file1 and file2
+
+ 8999 format(a6,I5,a1,A4,a1,A3,a1,A1,I4,A1,a3,3F8.3)
+
+c**************************************************************
+c************* Step-2: output 'aaa' (CA only) ***************
+c**************************************************************
+ OPEN(unit=7,file=outname,status='unknown') !pdb1.aln + pdb2.aln
+ 900 format(A)
+ 901 format('select atomno= ',I5)
+ 902 format('select atomno= ',I5)
+ write(7,900)'load inline'
+ write(7,900)'select atomno<50001'
+ write(7,900)'wireframe .45'
+ write(7,900)'select atomno>50000'
+ write(7,900)'wireframe .15'
+ write(7,900)'select all'
+ write(7,900)'color white' !all above atoms are white
+ do i=1,n_cut
+ 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'
+ write(7,900)'exit'
+ write(7,514)rmsd_ali
+ 514 format('REMARK RMSD of the common residues=',F8.3)
+ write(7,515)score_max,d0
+ 515 format('REMARK TM-score=',f6.4,' (d0=',f5.2,')')
+ do i=1,nseqA
+ 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)i-1,i
+ enddo
+ do i=1,nseqB
+ 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)50000+i-1,50000+i
+ enddo
+ close(7)
+ 1236 format('ATOM ',i5,' CA ',A3,' ',A1,I4,A1,3X,3F8.3)
+ 1238 format('TER')
+ 1239 format('CONECT',I5,I5)
+
+
+c**************************************************************
+c************* Step-3: output 'aaa_atm' (all-atom) ***********
+c**************************************************************
+c outname=trim(outname)//'_atm'
+ outname=outname(1:len_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,nL(1)
+ ax=t2(1)+u2(1,1)*xx(1,i)+u2(1,2)*yy(1,i)+u2(1,3)*zz(1,i)
+ ay=t2(2)+u2(2,1)*xx(1,i)+u2(2,2)*yy(1,i)+u2(2,3)*zz(1,i)
+ az=t2(3)+u2(3,1)*xx(1,i)+u2(3,2)*yy(1,i)+u2(3,3)*zz(1,i)
+
+ write(7,8888)i,Aatom(1,i),Ares(1,i),Ach(1,i),
+ & Aresi(1,i),Ains(1,i),ax,ay,az
+ enddo
+ write(7,1238) !TER
+*** chain2:
+ do i=1,nL(2)
+ write(7,8888)50000+i,Aatom(2,i),Ares(2,i),Ach(2,i),
+ & Aresi(2,i),Ains(2,i),xx(2,i),yy(2,i),zz(2,i)
+ enddo
+ write(7,1238) !TER
+ close(7)
+ 8888 format('ATOM ',I5,1x,A4,1x,A3,' ',A1,I4,A1,3x,3F8.3)
+
+*^^^^^^^^^^^^^^^^^^ output finished ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
9999 END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
=====================================
debian/changelog
=====================================
@@ -1,3 +1,12 @@
+tm-align (20190708+dfsg-1) unstable; urgency=medium
+
+ * New upstream version
+ * Packaged with routine-update script
+ * debhelper-compat 12
+ * Standards-Version: 4.4.0
+
+ -- Steffen Moeller <moeller at debian.org> Mon, 29 Jul 2019 23:00:47 +0200
+
tm-align (20170708+dfsg-2) unstable; urgency=medium
* Point Vcs fields to salsa.debian.org
=====================================
debian/compat deleted
=====================================
@@ -1 +0,0 @@
-11
=====================================
debian/control
=====================================
@@ -5,9 +5,9 @@ Uploaders: Steffen Moeller <moeller at debian.org>,
Andreas Tille <tille at debian.org>
Section: science
Priority: optional
-Build-Depends: debhelper (>= 11~),
+Build-Depends: debhelper-compat (= 12),
gfortran
-Standards-Version: 4.2.1
+Standards-Version: 4.4.0
Vcs-Browser: https://salsa.debian.org/med-team/tm-align
Vcs-Git: https://salsa.debian.org/med-team/tm-align.git
Homepage: http://zhanglab.ccmb.med.umich.edu/TM-align/
=====================================
debian/rules
=====================================
@@ -1,9 +1,5 @@
#!/usr/bin/make -f
-# debian/rules for tm-align
-
-# Uncomment this to turn on verbose mode.
-#export DH_VERBOSE=1
-
+export DH_VERBOSE=1
export DEB_BUILD_MAINT_OPTIONS = hardening=+all
%:
@@ -11,7 +7,7 @@ export DEB_BUILD_MAINT_OPTIONS = hardening=+all
override_dh_auto_build:
for src in `ls *.f` ; do \
- gfortran -O3 -ffast-math $(CPPFLAGS) $(LDFLAGS) -lm -o `basename $${src} .f` $${src} ; \
+ gfortran -O3 -g -ffast-math $(CPPFLAGS) $(LDFLAGS) -lm -o `basename $${src} .f` $${src} ; \
done
override_dh_auto_clean:
View it on GitLab: https://salsa.debian.org/med-team/tm-align/compare/8e8034d1eabe3ad26b4b061880bae25b48e36d1e...ee122aa039791b74992629dbc3c56415bd766564
--
View it on GitLab: https://salsa.debian.org/med-team/tm-align/compare/8e8034d1eabe3ad26b4b061880bae25b48e36d1e...ee122aa039791b74992629dbc3c56415bd766564
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/debian-med-commit/attachments/20190729/e4272a7b/attachment-0001.html>
More information about the debian-med-commit
mailing list