******************************************************************************* * This program is to identify the best alignment of two protein complex * structures and output the alignment with best TM-score. By default, * TM-score is normalized by the second complex. For comments/suggestions, * please contact email: srayanta@umich.edu or zhng@umich.edu) * * Reference to cite: * Srayanta Mukherjee, Yang Zhang. Nucleic Acids Research 2009; 37:e83 * * Permission to use, copy, modify, and distribute this program for * any purpose, with or without fee, is hereby granted, provided that * the notices on the head, the reference information, and this * copyright notice appear in all copies or substantial portions of * the Software. It is provided "as is" without express or implied * warranty. *********************** updating history ************************************** * 2009/01/27: Extension from dimer to multimer align * 2009/09/30: A bug in the output display was fixed * 2013/08/15: Fixed a number of compiler bugs ******************************************************************************* program compares PARAMETER(nmax=5000) PARAMETER(nmax2=10000) COMMON/BACKBONE/XA(3,nmax,0:1) common/dpc/score(nmax,nmax),gap_open,invmap(nmax) common/alignrst/invmap0(nmax) common/length/nseq1,nseq2 common/length1/nseq11,nseq22,nseq11_p,nseq22_p common/d0/d0,anseq common/d0min/d0_min common/d00/d00,d002 common/Tmsin/TMsin(nmax,nmax) character*100 fnam,pdb(100),outname,pdb1,pdb2 character*3 aa(-1:20),aanam,ss1(nmax),ss2(nmax) character*100 s,du character*200 outnameall_tmp,outnameall character seq1(0:nmax),seq2(0:nmax),sequ character aseq1(nmax2),aseq2(nmax2),aseq3(nmax2) character*1 temp character*2 tmp dimension xx(nmax),yy(nmax),zz(nmax) dimension m1(nmax),m2(nmax) dimension mf1(nmax),mf2(nmax) dimension xtm1(nmax),ytm1(nmax),ztm1(nmax) dimension xtm2(nmax),ytm2(nmax),ztm2(nmax) dimension xtmf1(nmax),ytmf1(nmax),ztmf1(nmax) dimension xtmf2(nmax),ytmf2(nmax),ztmf2(nmax) common/init/invmap_i(nmax) common/cutoff/int_cut common/TM/TM,TMmax,TM8temp common/n1n2/n1(nmax),n2(nmax) common/d8/d8 common/initial4/mm1(nmax),mm2(nmax) common/initial5/inter1(nmax),inter2(nmax),inter3(nmax) common/initial5/inter4(nmax) common /lens/ len_first1,len_first2,len_second1,len_second2 common/index/ind 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 data aa/ 'BCK','GLY','ALA','SER','CYS', & 'VAL','THR','ILE','PRO','MET', & 'ASP','ASN','LEU','LYS','GLU', & 'GLN','ARG','HIS','PHE','TYR', & 'TRP','CYX'/ character*1 slc(-1:20) data slc/'X','G','A','S','C', & 'V','T','I','P','M', & 'D','N','L','K','E', & 'Q','R','H','F','Y', & 'W','C'/ common/choice/ichachoice,ichain,ii integer a(nmax2) !common npm common/arr/ib(1000000,50),ic(1000000,50),id(1000000,50) common/arr1/ia(1000000,50) common/chain/ibig,ismall COMMON/BACKBONE/iXA(3,nmax,0:1) COMMON/BACKBONE/iXB(3,nmax,0:1) common/length/inseq1,inseq2 common/initial4/imm1(nmax),imm2(nmax),imm3(nmax),imm4(nmax) common/initial6/ifcounter1,ifcounter2 common/intscore/n_int character*100 ipdb*(100) character*100 is,idu character*3 iaa(-1:20),iaanam,iss1(nmax),iss2(nmax),iss3(nmax) character*3 iss4(nmax) character iseq1(0:nmax),iseq2(0:nmax),iseq3(0:nmax),iseq4(0:nmax) dimension ixx(nmax),iyy(nmax),izz(nmax) dimension im1(nmax),im2(nmax),im3(nmax),im4(nmax) common/flag/intflag,domflag common/initial4/mempo(nmax) common/chains/ichainterm1(0:nmax),ichainterm2(0:nmax) common/nchains/itercount,itercount1,icombo1,icombo2 Logical mtc integer ii,jj,chk real*8 dist,xaa,yaa,zaa,cov ccc data iaa/ 'BCK','GLY','ALA','SER','CYS', & 'VAL','THR','ILE','PRO','MET', & 'ASP','ASN','LEU','LYS','GLU', & 'GLN','ARG','HIS','PHE','TYR', & 'TRP','CYX'/ character*1 islc(-1:20) data islc/'X','G','A','S','C', & 'V','T','I','P','M', & '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(*,*) write(*,*)'Brief instruction for running MM-align program:' write(*,*)'(For detail: Mukherjee & Zhang, Nucl. Acid. Res.', & ' 37:e83, 2009)' write(*,*) write(*,*)'Please make sure that the two chains of both '// & 'complexes are seperated by a "TER".' write(*,*) write(*,*)'1. Align ''complex1.pdb'' to ''complex2.pdb''' write(*,*)' (By default, TM-score is normalized by the ', & 'length of ''complex2.pdb'')' write(*,*)' >MMalign complex1.pdb complex2.pdb' write(*,*) write(*,*)'2. Run MM-align and output the superposition ', & 'to ''MM.sup'' and ''MM.sup_all'':' write(*,*)' >MMalign complex1.pdb complex2.pdb -o MM.sup' write(*,*)' To view the superimposed structures of the', & ' aligned regions by rasmol:' write(*,*)' >rasmol -script MM.sup' write(*,*)' To view the superimposed structures of all', & ' regions by rasmol:' write(*,*)' >rasmol -script MM.sup_all' write(*,*) write(*,*)'3. If you want TM-score normalized by ', & 'an assigned length, e.g. 100 aa:' write(*,*)' >MMalign complex1.pdb complex2.pdb -L 100' write(*,*)' If you want TM-score normalized by the ', & 'average length of two structures:' write(*,*)' >MMalign complex1.pdb complex2.pdb -a' write(*,*)' If you want TM-score normalized by the ', & 'shorter length of two structures:' write(*,*)' >MMalign complex1.pdb complex2.pdb -b' write(*,*)' If you want TM-score normalized by the ', & 'longer length of two structures:' write(*,*)' >MMalign complex1.pdb complex2.pdb -c' write(*,*) write(*,*)'4. If you want MM-align to enforce interface'// & ' alignment defined by a cutoff e.g. 10A :' write(*,*)' >MMalign complex1.pdb complex2.pdb -I 10' write(*,*)' (Distance cutoff is used to determine interface', & ' resdiues. Default 8A.)' write(*,*) write(*,*)'5. If you want MM-align to allow cross'// & ' alignment e.g: Domain swap, Gene Fusion :' write(*,*)' >MMalign complex1.pdb complex2.pdb -DG ' write(*,*) write(*,*)'(First three options,do not change the', & ' final structure alignment result)' write(*,*) goto 9999 endif ******* options -----------> m_out=-1 !decided output m_fix=-1 !fixed length-scale only for output m_ave=-1 !using average length m_d0_min=-1 !diminum d0 for search m_d0=-1 !given d0 for both search and output narg=iargc() i=0 j=0 115 continue i=i+1 call getarg(i,fnam) if(fnam.eq.'-o')then m_out=1 i=i+1 call getarg(i,outname) elseif(fnam.eq.'-L')then !change both L_all and d0 m_fix=1 i=i+1 call getarg(i,fnam) read(fnam,*)L_fix elseif(fnam.eq.'-dmin')then m_d0_min=1 i=i+1 call getarg(i,fnam) read(fnam,*)d0_min_input elseif(fnam.eq.'-d0')then m_d0=1 i=i+1 call getarg(i,fnam) read(fnam,*)d0_fix elseif(fnam.eq.'-a')then ! this will change the superposed output but not the alignment m_ave=1 i=i+1 elseif(fnam.eq.'-b')then m_ave=2 i=i+1 elseif(fnam.eq.'-c')then m_ave=3 i=i+1 elseif(fnam.eq.'-I')then intflag=1 i=i+1 c m_out=1 call getarg(i,fnam) if(fnam.eq. ' ')then int_cut=8 else write(*,*)fnam read(fnam,*)int_cut endif elseif(fnam.eq.'-DG')then dom_flag=1 else j=j+1 pdb(j)=fnam endif if(i.lt.narg)goto 115 ccccc Searching for interface residues ************* if (intflag.eq.1)then pdb1=pdb(1) pdb2=pdb(2) call interface(pdb1,pdb2,int_cut) endif ccccc Determine no of chains ccccccccc ichainterm1(0)=1 ichainterm2(0)=1 open(unit=10,file=pdb(1),status='old') i=0 itercount=0 do while (.true.) read (10,9001,end=1234) s if (s(1:3).eq.'TER')then itercount=itercount+1 ichainterm1(itercount)=i+1 !-ichainterm1(itercount-1) endif 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 if(s(17:17).eq.' '.or.s(17:17).eq.'A')then i=i+1 endif endif endif enddo 1234 continue open(unit=10,file=pdb(2),status='old') i=0 itercount1=0 do while (.true.) read (10,9001,end=1235) s if (s(1:3).eq.'TER')then itercount1=itercount1+1 ichainterm2(itercount1)=i+1 endif 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 if(s(17:17).eq.' '.or.s(17:17).eq.'A')then i=i+1 endif endif endif enddo 1235 continue ccccccccc read data from first CA file: temp="" open(unit=10,file=pdb(1),status='old') i=0 l1=0 flag=0 terflag=0 counter=1 do while (.true.) read(10,9001,end=1010) s if(i.gt.0.and.s(1:3).eq.'TER'.and.terflag.eq.itercount)then goto 1010 endif if(i.gt.0.and.s(1:3).eq.'TER'.and.terflag.lt.itercount)then terflag=terflag+1 goto 1017 endif 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 if(s(17:17).eq.' '.or.s(17:17).eq.'A')then i=i+1 read(s,9000)du,aanam,du,mm1(i),du, $ xa(1,i,0),xa(2,i,0),xa(3,i,0) 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 endif if(s(10:11) .eq. " 2" .and. flag .eq. 0)then temp=s(22:22) flag=1 endif endif endif 1017 continue enddo 1010 continue 9000 format(A17,A3,A2,i4,A4,3F8.3) 9001 format(A100) close(10) nseq1=i c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ccccccccc read data from the second CA file: open(unit=10,file=pdb(2),status='old') i=0 l2=0 flag1=0 terflag1=0 do while (.true.) read(10,9001,end=1011) s if(i.gt.0.and.s(1:3).eq.'TER'.and.terflag1.eq.itercount1)then goto 1011 endif if(i.gt.0.and.s(1:3).eq.'TER'.and.terflag1.lt.itercount1)then terflag1=terflag1+1 goto 1018 endif 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 if(s(17:17).eq.' '.or.s(17:17).eq.'A')then i=i+1 read(s,9000)du,aanam,du,mm2(i),du, $ xa(1,i,1),xa(2,i,1),xa(3,i,1) 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 if(s(10:11) .eq. " 2" .and. flag1 .eq. 0)then temp=s(22:22) flag1=1 endif endif endif 1018 continue enddo 1011 continue close(10) nseq2=i ccccccError message in case there is only 1 chain if (terflag1.gt.50.and.terflag.gt.50)then write(*,*) write(*,*)'* Error. Too many chains in Complex files.'// & ' Please split chains and try again *' write(*,*) goto 9999 endif if (terflag.gt.50)then write(*,*) write(*,*)'* Error: Too many chains in complex 1.'// & ' Please split chains and try again *' write(*,*) goto 9999 endif if (terflag1.gt.50)then write(*,*) write(*,*)'* Error: Too many chains in complex 2.'// & ' Please split chains and try again *' write(*,*) goto 9999 endif c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ *!!! Scale of TM-score in search is based on the smaller protein ---------> d0_min=0.5 if(m_d0_min.eq.1)then d0_min=d0_min_input !for search endif anseq_min=min(nseq1,nseq2) anseq=anseq_min !length for defining TMscore in search d8=1.5*anseq_min**0.3+3.5 !remove pairs with dis>d8 during search & final if(anseq.gt.15)then d0=1.24*(anseq-15)**(1.0/3.0)-1.8 !scale for defining TM-score else d0=d0_min endif if(d0.lt.d0_min)d0=d0_min if(m_d0.eq.1)d0=d0_fix d00=d0 !for quickly calculate TM-score in searching if(d00.gt.8)d00=8 if(d00.lt.4.5)d00=4.5 d002=d00**2 nseq=max(nseq1,nseq2) do i=1,nseq n1(i)=i n2(i)=i enddo ******** Quickly calculate TM-score FOR ALL CHAINS INDIVIDUALLY ********* ichainterm1(0)=1 ichainterm2(0)=1 do ii=1,itercount do jj=1,itercount1 nseq11=ichainterm1(ii)-ichainterm1(ii-1) nseq22=ichainterm2(jj)-ichainterm2(jj-1) nseq11_p=ichainterm1(ii-1) nseq22_p=ichainterm2(jj-1) !print *,nseq11,nseq22 **** RUN INITIAL1 TO GET TM-SCORE d01=d0+1.5 if(d01.lt.d0_min)d01=d0_min d02=d01*d01 aL=min(nseq11,nseq22) idel=aL/2.5 !minimum size of considered fragment if(idel.le.5)idel=5 n11=-nseq22+idel n22=nseq11-idel GL_max=0 !print *,nseq22_p,nseq22+nseq22_p do ishift=n11,n22 L=0 K=0 nseq=nseq22+nseq22_p do j=nseq22_p,nseq k=k+1 i=k+ishift if(i.ge.1.and.i.le.nseq11)then L=L+1 invmap(j)=i+nseq11_p else invmap(j)=-1 endif enddo !print *,L if(L.ge.idel)then call get_GL11(GL) !print *,GL,GL_max if(GL.gt.GL_max)then GL_max=GL do i=nseq22_p,nseq22+nseq22_p invmap_i(i)=invmap(i) enddo endif endif enddo GL_max=GL_max/nseq22 TMsin(ii,jj)=GL_max !stop enddo !stop enddo !stop ********* iterating over all chain combos ***************** if(itercount.le.itercount1)then ibig = itercount1 ismall = itercount else ibig = itercount ismall = itercount1 endif ************************************************************* mtc=.false. icomb1=0 18 call nexper(ismall,a,mtc,even) if(a(1).ne.0)then icomb1=icomb1+1 do i=1,ismall ia(icomb1,i)=a(i) enddo endif if(icomb1.eq.2)then c write(*,*) " Too many Possible Combinations of Complex 1"// c & " Only First 1 million are tried out" goto 19 endif if(mtc) goto 18 19 continue cccc the smaller of the two proteins will be kept fixed cccc and various chain combinations of the bigger to be tried *** Generating combinations of the bigger chain: nPm **** *** where n = no of chains in bigger proteins **** *** m = no of chains in smaller proteins **** *********************************************************** mtc=.false. icomb=0 17 call nexper(ibig,a,mtc,even) if(a(1).ne.0)then icomb=icomb+1 do i=1,ibig ib(icomb,i)=a(i) enddo endif if(icomb.eq.1000)then cc write(*,*) " Too many Possible Combinations of Complex 2"// c & " Only First 1 million are tried out" c icomb=2 goto 16 endif if(mtc) goto 17 16 continue icnt=0 do i=1,icomb do j=1,ismall ic(i,j)=ib(i,j) enddo do ii=0,i-1 do jj=1,ismall if(ic(i,jj).eq.ic(ii,jj))then icnt=icnt+1 endif enddo if(icnt.eq.ismall)then do jj=1,ismall ic(i,jj)=0 enddo icnt=0 goto 117 endif icnt=0 enddo 117 continue icnt=0 enddo ind=0 do i=1,icomb if(ic(i,1).ne.0)then ind=ind+1 do j=1,ismall id(ind,j)=ic(i,j) enddo endif enddo *** calculating npm ******* !icomb=icomb+1 fact1=1 fact2=1 fact3=1 fnpm=1 do k=2,ibig fact1=fact1*k enddo do k=2,ismall fact2=fact2*k enddo icdiff=ibig-ismall if(icdiff.le.1)then fact3=1 else do k=2,icdiff fact3=fact3*k enddo endif if(ibig.eq.ismall)then fnpm=fact1 ! print *,"here" else ! print *,"no here" fnpm=real(fact1)/real(fact3) endif tempTM=0 n_al=0 itemp=0 rmsdtp=0 seqtp=0 !print *,ibig,ismall,fact1,fact3,fnpm !stop if(fnpm.gt.100)then fnpm=100 endif iss=0 ibb=0 sTM=0 sTMmax=0 iteration=0 icombo1=1 do icombo2=1,fnpm c write(*,*)icombo2 ichain=0 sTM=0 do i=1,ismall iss=ia(icombo1,i) ibb=id(icombo2,i) sTM=TMsin(iss,ibb)+sTM enddo !print *,sTM,sTMmax if(sTM.lt.sTMmax*0.95)goto 9998 if(sTM.ge.sTMmax)then sTMmax=sTM endif iteration=iteration + 1 if(fnpm.gt.5.and.iteration.eq.5)goto 9997 do j=1,nseq1 xtm1(j)=0 ytm1(j)=0 ztm1(j)=0 enddo do j=1,nseq2 xtm2(j)=0 ytm2(j)=0 ztm2(j)=0 enddo ***** do alignment ************************** CALL super_align !to find invmap(j) !enddo ************************************************************ *** resuperpose to find residues of dis in_align=0 n_int=0 n_int1=0 n_int2=0 !itemp=0 ix=0 iy=0 cov=0 n_al=0 do j=1,nseq2 if(invmap0(j).gt.0)then i=invmap0(j) c print *,i,j n_al=n_al+1 xtm1(n_al)=xa(1,i,0) ytm1(n_al)=xa(2,i,0) ztm1(n_al)=xa(3,i,0) xtm2(n_al)=xa(1,j,1) ytm2(n_al)=xa(2,j,1) ztm2(n_al)=xa(3,j,1) m1(n_al)=i !for recording residue order m2(n_al)=j endif enddo !stop d0_input=d0 call TMscore8(d0_input,n_al,xtm1,ytm1,ztm1,n1,n_al, & xtm2,ytm2,ztm2,n2,TM,Rcomm,Lcomm) !TM-score with dis d0_min=0.5 !for output anseq=nseq2 !length for defining final TMscore if(m_ave.eq.1)anseq=(nseq1+nseq2)/2.0 ! if(m_ave.eq.2)anseq=min(nseq1,nseq2) if(m_ave.eq.3)anseq=max(nseq1,nseq2) if(anseq.lt.anseq_min)anseq=anseq_min if(m_fix.eq.1)anseq=L_fix !input length if(anseq.gt.15)then d0=1.24*(anseq-15)**(1.0/3.0)-1.8 !scale for defining TM-score else d0=d0_min endif if(d0.lt.d0_min)d0=d0_min if(m_d0.eq.1)d0=d0_fix *** calculating interface aligned residues for final report----> if(intflag.eq.1)then iy=0 ix=0 i=0 icount1=0 icount2=0 dis2=0 icount1=ifcounter1*2 icount2=ifcounter2*2 do j=0,len_first2 if(invmap0(j).gt.0)then i=invmap0(j) if(inter1(i).eq.1 .and. inter3(j).eq.1)then dis2=sqrt((xa(1,i,0)-xa(1,j,1))**2+(xa(1,i,0)- & xa(2,j,1))**2+(xa(3,i,0)-xa(3,j,0))**2) iy=iy+1 if(dis2.le.d8)then ix=ix+1 endif endif endif enddo dis2=0 i=0 do j=len_first2,nseq2 if(invmap0(j).gt.0)then i=invmap0(j) if(inter2(i).eq.1 .and. inter4(j).eq.1)then dis2=sqrt((xa(1,i,0)-xa(1,j,1))**2+(xa(1,i,0)- & xa(2,j,1))**2+(xa(3,i,0)-xa(3,j,0))**2) iy=iy+1 if(dis2.le.d8)then ix=ix+1 endif endif endif enddo n_int=iy !n_int1+n_int2 cov=real(iy)/real(icount1) n_int5=ix cov5=real(ix)/real(icount1) endif *** remove dis>d8 in normal TM-score calculation for final report-----> j=0 n_eq=0 do i=1,n_al dis2=sqrt((xtm1(i)-xtm2(i))**2+(ytm1(i)-ytm2(i))**2+ & (ztm1(i)-ztm2(i))**2) if(dis2.le.d8)then j=j+1 xtm1(j)=xtm1(i) ytm1(j)=ytm1(i) ztm1(j)=ztm1(i) xtm2(j)=xtm2(i) ytm2(j)=ytm2(i) ztm2(j)=ztm2(i) m1(j)=m1(i) m2(j)=m2(i) if(ss1(m1(i)).eq.ss2(m2(i)))then n_eq=n_eq+1 endif endif enddo n8_al=j seq_id=float(n_eq)/(n_al+0.00000001) d0_input=d0 call TMscore(d0_input,n8_al,xtm1,ytm1,ztm1,n1,n8_al, & xtm2,ytm2,ztm2,n2,TM8,Rcomm,Lcomm) !normal TMscore rmsd8_al=Rcomm TM8temp=TM8*n8_al/anseq if(TM8temp.gt.tempTM)then tempTM=TM8temp itemp=n8_al rmsdtp=rmsd8_al seqtp=seq_id do i=1,itemp xtmf1(i)=xtm1(i) ytmf1(i)=ytm1(i) ztmf1(i)=ztm1(i) xtmf2(i)=xtm2(i) ytmf2(i)=ytm2(i) ztmf2(i)=ztm2(i) mf1(i)=m1(i) mf2(i)=m2(i) enddo endif 9998 continue enddo 9997 continue !enddo TM8=tempTM rmsd8_al=rmsdtp n8_al=itemp seq_id=seqtp ********* for output summary ****************************** write(*,*) write(*,*)'*****************************************************', & '*********************' write(*,*)'* MM-align ', & ' *' write(*,*)'* Aligning protein complex structure by Dynamic', & ' Programming *' write(*,*)'* Comments on the program, please email to: z', & 'hng@umich.edu *' write(*,*)'*****************************************************', & '*********************' write(*,*) write(*,101)pdb(1),nseq1 101 format('Protein 1:',A10,' Size=',I4) write(*,102)pdb(2),nseq2,int(anseq) 102 format('Protein 2:',A10,' Size=',I4, & ' (TM-score is normalized by ',I4,')') write(*,*) if (intflag.eq.1)then write(*,111)du,ifcounter1,ifcounter1 111 format('No. Of Interface Residue:'A4,'Protein 1,Chain 1 :',I4 & ' | Protein 1,Chain 2 :',I4) write(*,112)du,ifcounter2,ifcounter2 112 format(A29,'Protein 2,Chain 1 :',I4,' | Protein 2,Chain 2 :',I4) write(*,*) write(*,113)n8_al,rmsd8_al,TM8,seq_id 113 format('Aligned length=',I4, & ', RMSD=',f6.2,', TM-score=',f7.5,', ID=',f5.3) write(*,*) else write(*,103)n8_al,rmsd8_al,TM8,seq_id 103 format('Aligned length=',I4, & ', RMSD=',f6.2,', TM-score=',f7.5,', ID=',f5.3) write(*,*) endif ********* extract rotation matrix ------------> L=0 do i=1,n8_al k=mf1(i) L=L+1 r_1(1,L)=xa(1,k,0) r_1(2,L)=xa(2,k,0) r_1(3,L)=xa(3,k,0) r_2(1,L)=xtmf1(i) r_2(2,L)=ytmf1(i) r_2(3,L)=ztmf1(i) enddo if(L.gt.3)then call u3b(w,r_1,r_2,L,2,rms,u,t,ier) !u rotate r_1 to r_2 armsd=dsqrt(rms/L) write(*,*)'-------- rotation matrix to rotate Chain-1 to ', & 'Chain-2 ------' write(*,*)'i t(i) u(i,1) u(i,2) ', & ' u(i,3)' do i=1,3 write(*,204)i,t(i),u(i,1),u(i,2),u(i,3) enddo write(*,*) endif 204 format(I2,f18.10,f15.10,f15.10,f15.10) ********* for output superposition ****************************** if(m_out.eq.1)then 1237 format('ATOM ',i5,' CA ',A3,I6,4X,3F8.3) 1238 format('TER') 1239 format('CONECT',I5,I5) 900 format(A) 901 format('select atomno=',I4) 104 format('REMARK Chain 1:',A10,' Size=',I4) 105 format('REMARK Chain 2:',A10,' Size=',I4, & ' (TM-score is normalized by ',I4,')') 106 format('REMARK Aligned length=',I4,', RMSD=',f6.2, & ', TM-score=',f7.5,', ID=',f5.3) OPEN(unit=7,file=outname,status='unknown') !pdb1.aln + pdb2.aln *** script: write(7,900)'load inline' write(7,900)'select atomno<4000' write(7,900)'wireframe .45' write(7,900)'select none' write(7,900)'select atomno>4000' write(7,900)'wireframe .15' write(7,900)'color white' write(7,900)'background white' ikk1=0 ikk2=0 il1=0 il2=0 do i=1,n8_al if(mf1(i).ge.ichainterm1(ikk1))then ikk1=ikk1+1 endif il1=mod(ikk1,2) if(mf2(i).ge.ichainterm2(ikk2))then ikk2=ikk2+1 endif il2=mod(ikk2,2) dis2=sqrt((xtmf1(i)-xtmf2(i))**2+ & (ytmf1(i)-ytmf2(i))**2+(ztmf1(i)-ztmf2(i))**2) if(dis2.le.5)then if(mf1(i).ge.ichainterm1(ikk1-1).and.mf1(i).lt. & ichainterm1(ikk1).and.il1.ne.0)then write(7,901)mf1(i) write(7,900)'color red' else write(7,901)mf1(i) write(7,900)'color blue' endif if(mf2(i).ge.ichainterm2(ikk2-1).and.mf2(i).lt. & ichainterm2(ikk2).and.il2.ne.0)then write(7,901)4000+mf2(i) write(7,900)'color cyan' else write(7,901)4000+mf2(i) write(7,900)'color magenta' endif endif enddo write(7,900)'select all' write(7,900)'exit' write(7,104)pdb(1),nseq1 write(7,105)pdb(2),nseq2,int(anseq) write(7,106)n8_al,rmsd8_al,TM8,seq_id *** chain1: ik1=1 ik2=1 do i=1,n8_al write(7,1237)mf1(i),ss1(mf1(i)),mm1(mf1(i)), & xtmf1(i),ytmf1(i),ztmf1(i) enddo write(7,1238) !TER do i=2,n8_al if(mf1(i-1).le.ichainterm1(ik1).and. & mf1(i).gt.ichainterm1(ik1))then ik1=ik1+1 else write(7,1239)mf1(i-1),mf1(i) !connect atoms endif 1414 continue enddo *** chain2: do i=1,n8_al write(7,1237)4000+mf2(i),ss2(mf2(i)),mm2(mf2(i)), $ xtmf2(i),ytmf2(i),ztmf2(i) enddo write(7,1238) do i=2,n8_al if(mf2(i-1).le.ichainterm2(ik2).and. & mf2(i).gt.ichainterm2(ik2))then ik2=ik2+1 else write(7,1239)4000+mf2(i-1),4000+mf2(i) endif 1515 continue enddo close(7) ccc k=0 outnameall_tmp=outname//'_all' outnameall='' do i=1,200 if(outnameall_tmp(i:i).ne.' ')then k=k+1 outnameall(k:k)=outnameall_tmp(i:i) endif enddo OPEN(unit=8,file=outnameall,status='unknown') !pdb1.aln + pdb2.aln *** script: write(8,900)'load inline' write(8,900)'select atomno<4000' write(8,900)'wireframe .45' write(8,900)'select none' write(8,900)'select atomno>4000' write(8,900)'wireframe .15' write(8,900)'background white' ikk1=0 ikk2=0 il1=0 il2=0 do i=1,nseq1 if(i.ge.ichainterm1(ikk1))then ikk1=ikk1+1 endif il1=mod(ikk1,2) if(i.ge.ichainterm1(ikk1-1).and.i.lt. & ichainterm1(ikk1).and.il1.ne.0)then write(8,901)i write(8,900)'color red' else write(8,901)i write(8,900)'color blue' endif enddo do j=1,nseq2 if(j.ge.ichainterm2(ikk2))then ikk2=ikk2+1 endif il2=mod(ikk2,2) if(j.ge.ichainterm2(ikk2-1).and.j.lt. & ichainterm2(ikk2).and.il2.ne.0)then write(8,901)4000+j write(8,900)'color cyan' else write(8,901)4000+j write(8,900)'color magenta' endif enddo write(8,900)'select all' write(8,900)'exit' write(8,104)pdb(1),nseq1 write(8,105)pdb(2),nseq2,int(anseq) write(8,106)n8_al,rmsd8_al,TM8,seq_id *** chain1: ik1=1 ik2=1 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,1237)i,ss1(i),mm1(i),ax,ay,az enddo write(8,1238) !TER do i=2,nseq1 ij=i-1 if(ij.eq.ichainterm1(ik1)-1)then !print *,i,ichainterm1(ik1) ik1=ik1+1 else write(8,1239)i-1,i endif 1212 continue enddo *** chain2: do i=1,nseq2 write(8,1237)4000+i,ss2(i),mm2(i), & xa(1,i,1),xa(2,i,1),xa(3,i,1) enddo write(8,1238) do i=2,nseq2 ij=i-1 if(ij.eq.ichainterm2(ik2)-1)then !print *,i,ichainterm2(ik2) ik2=ik2+1 else write(8,1239)4000+i-1,4000+i endif 1313 continue enddo close(8) endif *^^^^^^^^^^^^^^^^^^ output finished ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ************ output aligned sequences ************************** ii=0 i1_old=1 i2_old=1 itimes=0 iflag1=0 iflag2=0 iflag3=0 iflag4=0 iflag5=0 iflag6=0 l=0 !k1=0 !k2=0 ll=0 l1=0 l2=0 ik1=1 ik2=1 ll1=0 ll2=0 do i=1,n8_al if(i1_old-1.lt.ichainterm1(ik1).and. & mf1(i)-1.ge.ichainterm1(ik1)-1)then itill1=ichainterm1(ik1)-1 iflag1=1 else itill1=mf1(i)-1 iflag1=0 endif if(i2_old-1.lt.ichainterm2(ik2).and. & mf2(i)-1.ge.ichainterm2(ik2)-1)then itill2=ichainterm2(ik2)-1 iflag2=1 else itill2=mf2(i)-1 iflag2=0 endif do j=i1_old,itill1 !mf1(i)-1 if(j.ge.ichainterm1(ik1))then ik1=ik1+1 goto 7777 endif ii=ii+1 l1=mod(ik1,2) if(j.ge.ichainterm1(ik1-1).and.j.lt.ichainterm1 & (ik1).and. l1.ne.0)then isequence1=ichar(seq1(j)) isequence1=isequence1+32 aseq1(ii)=char(isequence1) else aseq1(ii)=seq1(j) endif aseq2(ii)='-' aseq3(ii)=' ' enddo 7777 continue do j=i2_old,itill2 !mf2(i)-1 if(j.ge.ichainterm2(ik2))then ik2=ik2+1 !print *,"going out" !goto 7774 endif ii=ii+1 aseq1(ii)='-' l2=mod(ik2,2) if(j.ge.ichainterm2(ik2-1).and.j.lt.ichainterm2(ik2).and. & l2.ne.0)then isequence=ichar(seq2(j)) isequence=isequence+32 aseq2(ii)=char(isequence) else aseq2(ii)=seq2(j) endif aseq3(ii)=' ' enddo 7774 continue 7776 continue if(iflag1.eq.1)then do j=itill1+1,mf1(i)-1 iflag3=0 if(j.ge.ichainterm1(ik1))then ik1=ik1+1 itill1=j-1 iflag3=1 goto 7773 endif ii=ii+1 l1=mod(ik1,2) if(j.ge.ichainterm1(ik1-1).and.j.lt.ichainterm1 & (ik1).and. l1.ne.0)then isequence1=ichar(seq1(j)) isequence1=isequence1+32 aseq1(ii)=char(isequence1) else aseq1(ii)=seq1(j) endif aseq2(ii)='-' aseq3(ii)=' ' enddo !iflag1=0 endif iflag1=0 7773 continue 7778 continue if(iflag2.eq.1)then do j=itill2+1,mf2(i)-1 iflag4=0 if(j.ge.ichainterm2(ik2))then ik2=ik2+1 itill2=j-1 iflag4=1 goto 7775 endif ii=ii+1 aseq1(ii)='-' l2=mod(ik2,2) if(j.ge.ichainterm2(ik2-1).and.j.lt.ichainterm2(ik2).and. & l2.ne.0)then isequence=ichar(seq2(j)) isequence=isequence+32 aseq2(ii)=char(isequence) else aseq2(ii)=seq2(j) endif aseq3(ii)=' ' enddo !iflag2=0 endif iflag2=0 7775 continue if(iflag3.eq.1.and.iflag1.eq.1.and.itill1.lt.mf1(i)-1)goto 7776 if(iflag4.eq.1.and.iflag2.eq.1.and.itill2.lt.mf2(i)-1)goto 7778 ii=ii+1 if(mf1(i).ge.ichainterm1(ik1))then ik1=ik1+1 endif l=mod(ik1,2) if(mf2(i).ge.ichainterm2(ik2))then ik2=ik2+1 endif ll=mod(ik2,2) k=0 if(mf1(i).ge.ichainterm1(ik1-1).and.mf1(i).lt.ichainterm1(ik1) & .and.l.ne.0)then isequence1=ichar(seq1(mf1(i))) isequence1=isequence1+32 aseq1(ii)=char(isequence1) else aseq1(ii)=seq1(mf1(i)) endif if(mf2(i).ge.ichainterm2(ik2-1).and.mf2(i).lt.ichainterm2(ik2) & .and.ll.ne.0)then isequence=ichar(seq2(mf2(i))) isequence=isequence+32 aseq2(ii)=char(isequence) else aseq2(ii)=seq2(mf2(i)) endif dis2=sqrt((xtmf1(i)-xtmf2(i))**2+ & (ytmf1(i)-ytmf2(i))**2+(ztmf1(i)-ztmf2(i))**2) if(dis2.le.5)then aseq3(ii)=':' else aseq3(ii)='.' endif i1_old=mf1(i)+1 i2_old=mf2(i)+1 enddo ********************************** if(ichainterm1(ik1).lt.nseq1)then itill1=ichainterm1(ik1) iflag1=1 else itill1=nseq1 iflag1=0 endif if(ichainterm2(ik2).lt.nseq2)then itill2=ichainterm2(ik2) iflag2=1 else itill2=nseq2 iflag2=0 endif do i=i1_old,itill1 if(i.ge.ichainterm1(ik1))then ik1=ik1+1 goto 6669 endif ii=ii+1 ll1=mod(ik1,2) !print *,i,ii,seq1(i),ik1 if(i.ge.ichainterm1(ik1-1).and.i.lt.ichainterm1(ik1) & .and.ll1.ne.0)then isequence1=ichar(seq1(i)) isequence1=isequence1+32 aseq1(ii)=char(isequence1) aseq2(ii)='-' aseq3(ii)=' ' !print *,"HERE",aseq1(ii) else isequence1=ichar(seq1(i)) aseq1(ii)=char(isequence1) aseq2(ii)='-' aseq3(ii)=' ' !print *,aseq1(ii),"HERE",i,ichainterm1(ik1),ll1 endif enddo 6669 continue do i=i2_old,itill2 if(i.ge.ichainterm2(ik2))then ik2=ik2+1 goto 6668 endif ii=ii+1 ll2=mod(ik2,2) if(i.ge.ichainterm2(ik2-1).and.i.lt.ichainterm2(ik2) & .and.ll2.ne.0)then aseq1(ii)='-' isequence=ichar(seq2(i)) isequence=isequence+32 aseq2(ii)=char(isequence) aseq3(ii)=' ' else aseq1(ii)='-' isequence=ichar(seq2(i)) aseq2(ii)=char(isequence) aseq3(ii)=' ' endif enddo 6668 continue 6667 continue if(iflag1.eq.1)then do i=itill1,nseq1 iflag5=0 if(i.ge.ichainterm1(ik1))then ik1=ik1+1 itill1=i iflag5=1 goto 6664 endif ii=ii+1 ll1=mod(ik1,2) !print *,i,ii,seq1(i),ik1 if(i.ge.ichainterm1(ik1-1).and.i.lt.ichainterm1(ik1) & .and.ll1.ne.0)then isequence1=ichar(seq1(i)) isequence1=isequence1+32 aseq1(ii)=char(isequence1) aseq2(ii)='-' aseq3(ii)=' ' else isequence1=ichar(seq1(i)) aseq1(ii)=char(isequence1) !print *, "No HERe" ,aseq1(ii),i,ichainterm1(ik1-1),ll1 aseq2(ii)='-' aseq3(ii)=' ' endif enddo endif iflag1=0 6664 continue 6665 continue if(iflag2.eq.1)then do i=itill2,nseq2 iflag6=0 if(i.ge.ichainterm2(ik2))then ik2=ik2+1 itill2=i iflag6=1 goto 6666 endif ii=ii+1 ll2=mod(ik2,2) if(i.ge.ichainterm2(ik2-1).and.i.lt.ichainterm2(ik2) & .and.ll2.ne.0)then aseq1(ii)='-' isequence=ichar(seq2(i)) isequence=isequence+32 aseq2(ii)=char(isequence) aseq3(ii)=' ' else aseq1(ii)='-' isequence=ichar(seq2(i)) aseq2(ii)=char(isequence) aseq3(ii)=' ' endif enddo endif iflag2=0 6666 continue if(iflag5.eq.1.and.itill1.lt.nseq1-1.and.iflag1.eq.1)goto 6667 if(iflag6.eq.1.and.itill2.lt.nseq2-1.and.iflag2.eq.1)goto 6665 write(*,50) 50 format('(":" denotes the residue pairs of distance < 5.0 ', & 'Angstrom)') write(*,10)(aseq1(i),i=1,ii) write(*,10)(aseq3(i),i=1,ii) write(*,10)(aseq2(i),i=1,ii) 10 format(10000A1) write(*,*) if(intflag.eq.1)then write(*,2222) cov 2222 format('Total Aligned Interface Coverage= ',f0.3) write(*,*) write(*,*)'(Odd no chains of both proteins are in Lower '// & 'case and Even no chains of both proteins are in '// & 'UPPER case)' else write(*,*)'(Odd no chains of both proteins are in lower '// & 'case and Even no chains of both proteins are in '// & 'UPPER case)' endif write(*,*) 9999 END *********************************************************************** *********************************************************************** * Structure superposition *********************************************************************** *********************************************************************** *********************************************************************** SUBROUTINE super_align PARAMETER(nmax=5000) COMMON/BACKBONE/XA(3,nmax,0:1) common/length/nseq1,nseq2 common/dpc/score(nmax,nmax),gap_open,invmap(nmax) common/alignrst/invmap0(nmax) common/zscore/zrms,n_al,rmsd_al common/TM/TM,TMmax common/init/invmap_i(nmax) common/init1/invmap_ii(nmax) dimension gapp(100) TMmax=0 TMtemp=0 n_gapp=2 gapp(1)=-0.6 gapp(2)=0 c n_gapp=11 c do i=1,n_gapp c gapp(i)=-(n_gapp-i) c enddo *11111111111111111111111111111111111111111111111111111111 * get initial alignment from gapless threading ********************************************************** call get_initial !gapless threading do i=1,nseq2 invmap(i)=invmap_i(i) !with highest zcore enddo call get_score !TM, matrix score(i,j) if(TM.gt.TMmax)then TMmax=TM do j=1,nseq2 invmap0(j)=invmap(j) enddo endif ***************************************************************** * initerative alignment, for different gap_open: ***************************************************************** DO 111 i_gapp=1,n_gapp !different gap panalties GAP_OPEN=gapp(i_gapp) !gap panalty do 222 id=1,30 !maximum interation is 200 call DP(NSEQ1,NSEQ2) !produce alignment invmap(j) * Input: score(i,j), and gap_open * Output: invmap(j) call get_score !calculate TM-score, score(i,j) !print *,"In1a",TM,TMmax c record the best alignment in whole search ----------> if(TM.gt.TMmax)then TMmax=TM do j=1,nseq2 invmap0(j)=invmap(j) enddo endif if(id.gt.1)then diff=abs(TM-TM_old) if(diff.lt.0.000001)goto 33 endif TM_old=TM 222 continue 33 continue 111 continue *222222222222222222222222222222222222222222222222222222222 * get initial alignment from secondary structure alignment ********************************************************** call get_initial2 !DP for secondary structure do i=1,nseq2 invmap(i)=invmap_i(i) !with highest zcore enddo call get_score !TM, score(i,j) !print *,"Ini2",TM,TMmax if(TM.gt.TMmax)then TMmax=TM do j=1,nseq2 invmap0(j)=invmap(j) enddo endif if(TM.gt.TMmax*0.2)then ***************************************************************** * initerative alignment, for different gap_open: ***************************************************************** DO 1111 i_gapp=1,n_gapp !different gap panalties GAP_OPEN=gapp(i_gapp) !gap panalty do 2222 id=1,30 !maximum interation is 200 call DP(NSEQ1,NSEQ2) !produce alignment invmap(j) * Input: score(i,j), and gap_open * Output: invmap(j) call get_score !calculate TM-score, score(i,j) ! print *,"Ini2a",TM,TMmax c write(*,21)gap_open,rmsd_al,n_al,TM c record the best alignment in whole search ----------> if(TM.gt.TMmax)then TMmax=TM do j=1,nseq2 invmap0(j)=invmap(j) enddo endif if(id.gt.1)then diff=abs(TM-TM_old) if(diff.lt.0.000001)goto 333 endif TM_old=TM 2222 continue 333 continue 1111 continue endif ******************************************************************* * get initial alignment of local structure superposition ******************************************************************* call get_initial5 do i=1,nseq2 invmap(i)=invmap_i(i) !with highest zcore enddo call get_score !TM, matrix score(i,j) if(TM.gt.TMmax)then TMmax=TM do j=1,nseq2 invmap0(j)=invmap(j) enddo endif if(TM.gt.TMmax*0.2)then ***************************************************************** * initerative alignment, for different gap_open: ***************************************************************** DO 88 i_gapp=2,n_gapp !different gap panalties GAP_OPEN=gapp(i_gapp) !gap panalty do 888 id=1,2 !maximum interation is 200 call DP(NSEQ1,NSEQ2) !produce alignment invmap(j) * Input: score(i,j), and gap_open * Output: invmap(j) call get_score !calculate TM-score, score(i,j) c print *,"Ini5a",TM,TMmax c record the best alignment in whole search ----------> if(TM.gt.TMmax)then TMmax=TM do j=1,nseq2 invmap0(j)=invmap(j) enddo endif if(id.gt.1)then diff=abs(TM-TM_old) if(diff.lt.0.000001)goto 880 endif TM_old=TM 888 continue 880 continue 88 continue endif *333333333333333333333333333333333333333333333333333333333333 * get initial alignment from invmap0+SS ************************************************************* call get_initial3 !invmap0+SS do i=1,nseq2 invmap(i)=invmap_i(i) !with highest zcore enddo call get_score !TM, score(i,j) ! print *,"Ini3",TM,TMmax if(TM.gt.TMmax)then TMmax=TM do j=1,nseq2 invmap0(j)=invmap(j) enddo endif if(TM.gt.TMmax*0.2)then ***************************************************************** * initerative alignment, for different gap_open: ***************************************************************** DO 1110 i_gapp=1,n_gapp !different gap panalties GAP_OPEN=gapp(i_gapp) !gap panalty do 2220 id=1,30 !maximum interation is 200 call DP(NSEQ1,NSEQ2) !produce alignment invmap(j) * Input: score(i,j), and gap_open * Output: invmap(j) call get_score !calculate TM-score, score(i,j) !print *,"Ini3a",TM,TMmax c write(*,21)gap_open,rmsd_al,n_al,TM c record the best alignment in whole search ----------> if(TM.gt.TMmax)then TMmax=TM do j=1,nseq2 invmap0(j)=invmap(j) enddo endif if(id.gt.1)then diff=abs(TM-TM_old) if(diff.lt.0.000001)goto 330 endif TM_old=TM 2220 continue 330 continue 1110 continue endif *444444444444444444444444444444444444444444444444444444444 * get initial alignment of pieces from gapless threading ********************************************************** call get_initial4 !gapless threading do i=1,nseq2 invmap(i)=invmap_i(i) !with highest zcore enddo do i=1,nseq2 do j=1,nseq1 score(i,j)=0 enddo enddo call get_score !TM, matrix score(i,j) !print *,"Ini4",TM,TMmax if(TM.gt.TMmax)then TMmax=TM do j=1,nseq2 invmap0(j)=invmap(j) enddo endif if(TM.gt.TMmax*0.2)then ***************************************************************** * initerative alignment, for different gap_open: ***************************************************************** DO 44 i_gapp=2,n_gapp !different gap panalties GAP_OPEN=gapp(i_gapp) !gap panalty do 444 id=1,2 !maximum interation is 200 call DP(NSEQ1,NSEQ2) !produce alignment invmap(j) * Input: score(i,j), and gap_open * Output: invmap(j) call get_score !calculate TM-score, score(i,j) ! print *,"Ini4a",TM,TMmax c record the best alignment in whole search ----------> if(TM.gt.TMmax)then TMmax=TM do j=1,nseq2 invmap0(j)=invmap(j) enddo endif 444 continue 44 continue endif c^^^^^^^^^^^^^^^ best alignment invmap0(j) found ^^^^^^^^^^^^^^^^^^ RETURN END ************************************************************** * get initial alignment invmap0(i) from gapless threading ************************************************************** subroutine get_initial PARAMETER(nmax=5000) COMMON/BACKBONE/XA(3,nmax,0:1) common/length/nseq1,nseq2 common/dpc/score(nmax,nmax),gap_open,invmap(nmax) common/alignrst/invmap0(nmax) common/zscore/zrms,n_al,rmsd_al common/TM/TM,TMmax common/init/invmap_i(nmax) common/init1/invmap_ii(nmax) common/d0/d0,anseq common/d0min/d0_min common/speed/align_len1 common/nchains/itercount,itercount1,icombo1,icombo2 common/chain/ibig,ismall common/chains/ichainterm1(0:nmax),ichainterm2(0:nmax) real n_cut integer n_all 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 data w /nmax*1.0/ dimension xtmm1(nmax),ytmm1(nmax),ztmm1(nmax) dimension xtmm2(nmax),ytmm2(nmax),ztmm2(nmax) c d01=d0+1.5 c if(d01.lt.d0_min)d01=d0_min c d02=d01*d01 c n_jump=5 c n_cut=0.75 c aL=min(nseq1,nseq2) c idel=aL/2.3 !minimum size of considered fragment c if(idel.le.5)idel=5 c n1=-nseq2+idel c n2=nseq1-idel c GL_max=0 c GL_maxA=0 c GL_max_cut=0.95 c do ishift=n1,n2,n_jump c L=0 c do j=1,nseq2 c i=j+ishift c if(i.ge.1.and.i.le.nseq1)then c L=L+1 c invmap(j)=i c xtmm1(L)=xa(1,i,0) c ytmm1(L)=xa(1,i,0) c ztmm1(L)=xa(1,i,0) c xtmm2(L)=xa(1,j,1) c ytmm2(L)=xa(1,j,1) c ztmm2(L)=xa(1,j,1) c else c invmap(j)=-1 c endif c enddo do ii=1,itercount do jj=1,itercount1 nseq11=ichainterm1(ii)-ichainterm1(ii-1) nseq22=ichainterm2(jj)-ichainterm2(jj-1) nseq11_p=ichainterm1(ii-1) nseq22_p=ichainterm2(jj-1) !print *,nseq11,nseq22 **** RUN INITIAL1 TO GET TM-SCORE d01=d0+1.5 if(d01.lt.d0_min)d01=d0_min d02=d01*d01 aL=min(nseq11,nseq22) idel=aL/2.5 !minimum size of considered fragment if(idel.le.5)idel=5 n11=-nseq22+idel n22=nseq11-idel GL_max=0 !print *,nseq22_p,nseq22+nseq22_p do ishift=n11,n22 L=0 K=0 nseq=nseq22+nseq22_p do j=nseq22_p,nseq k=k+1 i=k+ishift if(i.ge.1.and.i.le.nseq11)then L=L+1 invmap(j)=i+nseq11_p xtmm1(L)=xa(1,i,0) ytmm1(L)=xa(1,i,0) ztmm1(L)=xa(1,i,0) xtmm2(L)=xa(1,j,1) ytmm2(L)=xa(1,j,1) ztmm2(L)=xa(1,j,1) else invmap(j)=-1 endif enddo if(L.ge.idel)then call get_GL(GL) if(GL.gt.GL_max)then GL_max=GL do i=1,nseq2 invmap_i(i)=invmap(i) enddo endif endif enddo enddo enddo return end ************************************************************** * get initial alignment invmap0(i) from secondary structure ************************************************************** subroutine get_initial2 PARAMETER(nmax=5000) COMMON/BACKBONE/XA(3,nmax,0:1) common/length/nseq1,nseq2 common/dpc/score(nmax,nmax),gap_open,invmap(nmax) common/alignrst/invmap0(nmax) common/zscore/zrms,n_al,rmsd_al common/TM/TM,TMmax common/sec/isec(nmax),jsec(nmax) common/init/invmap_i(nmax) ********** assign secondary structures *************** c 1->coil, 2->helix, 3->turn, 4->strand do i=1,nseq1 isec(i)=1 j1=i-2 j2=i-1 j3=i j4=i+1 j5=i+2 if(j1.ge.1.and.j5.le.nseq1)then dis13=diszy(0,j1,j3) dis14=diszy(0,j1,j4) dis15=diszy(0,j1,j5) dis24=diszy(0,j2,j4) dis25=diszy(0,j2,j5) dis35=diszy(0,j3,j5) isec(i)=make_sec(dis13,dis14,dis15,dis24,dis25,dis35) endif enddo do i=1,nseq2 jsec(i)=1 j1=i-2 j2=i-1 j3=i j4=i+1 j5=i+2 if(j1.ge.1.and.j5.le.nseq2)then dis13=diszy(1,j1,j3) dis14=diszy(1,j1,j4) dis15=diszy(1,j1,j5) dis24=diszy(1,j2,j4) dis25=diszy(1,j2,j5) dis35=diszy(1,j3,j5) jsec(i)=make_sec(dis13,dis14,dis15,dis24,dis25,dis35) endif enddo call smooth !smooth the assignment ********** score matrix ************************** do i=1,nseq1 do j=1,nseq2 if(isec(i).eq.jsec(j))then score(i,j)=1 else score(i,j)=0 endif enddo enddo ********** find initial alignment: invmap(j) ************ gap_open=-1.0 !should be -1 call DP(NSEQ1,NSEQ2) !produce alignment invmap(j) do i=1,nseq2 invmap_i(i)=invmap(i) enddo *^^^^^^^^^^^^ initial alignment done ^^^^^^^^^^^^^^^^^^^^^^ return end ************************************************************** * get initial alignment invmap0(i) from secondary structure * and previous alignments ************************************************************** subroutine get_initial3 PARAMETER(nmax=5000) COMMON/BACKBONE/XA(3,nmax,0:1) common/length/nseq1,nseq2 common/dpc/score(nmax,nmax),gap_open,invmap(nmax) common/alignrst/invmap0(nmax) common/zscore/zrms,n_al,rmsd_al common/TM/TM,TMmax common/sec/isec(nmax),jsec(nmax) common/init/invmap_i(nmax) ********** score matrix ************************** do i=1,nseq2 invmap(i)=invmap0(i) enddo call get_score1 !get score(i,j) using RMSD martix do i=1,nseq1 do j=1,nseq2 if(isec(i).eq.jsec(j))then score(i,j)=0.5+score(i,j) else score(i,j)=score(i,j) endif enddo enddo ********** find initial alignment: invmap(j) ************ gap_open=-1.0 !should be -1 call DP(NSEQ1,NSEQ2) !produce alignment invmap(j) do i=1,nseq2 invmap_i(i)=invmap(i) enddo *^^^^^^^^^^^^ initial alignment done ^^^^^^^^^^^^^^^^^^^^^^ return end ************************* ************************************* * get initial alignment invmap0(i) from fragment gapless threading ************************************************************** subroutine get_initial4 PARAMETER(nmax=5000) COMMON/BACKBONE/XA(3,nmax,0:1) common/length/nseq1,nseq2 common/dpc/score(nmax,nmax),gap_open,invmap(nmax) common/alignrst/invmap0(nmax) common/zscore/zrms,n_al,rmsd_al common/TM/TM,TMmax common/init/invmap_i(nmax) common/initial4/mm1(nmax),mm2(nmax) logical contin dimension ifr2(2,nmax,nmax),Lfr2(2,nmax),Lfr_max2(2),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=3.85 c dcu_min=3.65 ccc Find the smallest continuous fragments --------> do i=1,nseq1 mm(1,i)=mm1(i) enddo do i=1,nseq2 mm(2,i)=mm2(i) enddo do k=1,2 dcu=dcu0 if(k.eq.1)then nseq0=nseq1 r_min=nseq1/3.0 !minimum fragment, in case too small protein else nseq0=nseq2 r_min=nseq2/3.0 !minimum fragment, in case too small protein endif if(r_min.gt.fra_min)r_min=fra_min 20 nfr=1 !number of fragments j=1 !number of residue at nf-fragment ifr2(k,nfr,j)=1 !what residue Lfr2(k,nfr)=j !length of the fragment do i=2,nseq0 dis=diszy(k-1,i-1,i) contin=.false. if(dcu.gt.dcu0)then if(dis.lt.dcu)then contin=.true. endif elseif(mm(k,i).eq.(mm(k,i-1)+1))then if(dis.lt.dcu)then contin=.true. endif endif if(contin)then j=j+1 ifr2(k,nfr,j)=i Lfr2(k,nfr)=j else nfr=nfr+1 j=1 ifr2(k,nfr,j)=i Lfr2(k,nfr)=j endif enddo Lfr_max=0 i_fr2(k)=1 !ID of the maximum piece do i=1,nfr if(Lfr_max.lt.Lfr2(k,i))then Lfr_max=Lfr2(k,i) i_fr2(k)=i endif enddo if(Lfr_max.lt.r_min)then dcu=1.1*dcu goto 20 endif enddo c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ccc select what piece will be used (this may araise ansysmetry, but ccc only when L1=L2 and Lfr1=Lfr2 and L1 ne Lfr1 ccc if L1=Lfr1 and L2=Lfr2 (normal proteins), it will be the same as initial1 mark=1 if(Lfr2(1,i_fr2(1)).lt.Lfr2(2,i_fr2(2)))then mark=1 elseif(Lfr2(1,i_fr2(1)).gt.Lfr2(2,i_fr2(2)))then mark=2 else !Lfr1=Lfr2 if(nseq1.le.nseq2)then mark=1 else mark=2 endif endif ccc L_fr=Lfr2(mark,i_fr2(mark)) do i=1,L_fr ifr(i)=ifr2(mark,i_fr2(mark),i) enddo ccc if(mark.eq.1)then !non-redundant to get_initial1 nseq0=nseq1 else nseq0=nseq2 endif if(L_fr.eq.nseq0)then n1=int(nseq0*0.1) !0 n2=int(nseq0*0.89) !2 j=0 do i=n1,n2 j=j+1 ifr(j)=ifr(n1+j) enddo L_fr=j endif ccc get initial -------------> if(mark.eq.1)then !nseq1 as the smallest one nseq1_=L_fr aL=min(nseq1_,nseq2) idel=aL/2.5 !minimum size of considered fragment if(idel.le.fra_min1)idel=fra_min1 n1=-nseq2+idel !shift1 n2=nseq1_-idel !shift2 GL_max=0 do ishift=n1,n2 L=0 do j=1,nseq2 i=j+ishift if(i.ge.1.and.i.le.nseq1_)then L=L+1 invmap(j)=ifr(i) else invmap(j)=-1 endif enddo if(L.ge.idel)then call get_GL(GL) if(GL.gt.GL_max)then GL_max=GL do i=1,nseq2 invmap_i(i)=invmap(i) enddo endif endif enddo else !@@@@@@@@@@@@@@@@@@@@ nseq2_=L_fr aL=min(nseq1,nseq2_) idel=aL/2.5 !minimum size of considered fragment if(idel.le.fra_min1)idel=fra_min1 n1=-nseq2_+idel n2=nseq1-idel GL_max=0 do ishift=n1,n2 L=0 do j=1,nseq2 invmap(j)=-1 enddo do j=1,nseq2_ i=j+ishift if(i.ge.1.and.i.le.nseq1)then L=L+1 invmap(ifr(j))=i endif enddo if(L.ge.idel)then call get_GL(GL) if(GL.gt.GL_max)then GL_max=GL do i=1,nseq2 invmap_i(i)=invmap(i) enddo endif endif enddo endif return end ************************************************************** * fifth initial alignement. Using local structure super * * position. * ************************************************************** subroutine get_initial5 PARAMETER (nmax=5000) COMMON/BACKBONE/XA(3,nmax,0:1) common/length/nseq1,nseq2 common/dpc/score(nmax,nmax),gap_open,invmap(nmax) common/alignrst/invmap0(nmax) common/zscore/zrms,n_al,rmsd_al common/TM/TM,TMmax common/d0/d0,anseq common/d0min/d0_min common/init/invmap_i(nmax) common/initial4/mm1(nmax),mm2(nmax) common/tinit/invmap_ii(nmax) common/local/lflag common/sec/isec(nmax),jsec(nmax) logical contin integer n_all,secmatch,coilcnt real n_max_cut,n_all_max,n_cutt 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 data w /nmax*1.0/ common/inv/invmap_a(nmax) common/speed/align_len1 ************* points needed to tuned ****************** c 1, reduce number of iteration in get_GL from 3 to 2 or 1? c this should reduce the time for get_initial1 as well. c 2, put get_DP into get_initial1 as well? Then you need to c reduced the number of gapless threading by jump. c 3, in this subroutine, n_frag, n_jump, GLmaxA_cut, d02, n_frag2 c need to be tuned c 4, there are serveral bugs in your codes such as the edge-condition c Fixing these bugs will also improve the results. c Also, always put the common calculations outside of the circles ************************************************************ ***** setting parameters ************************************ n_frag=5 c n_frag=15 n_jump=n_frag c n_jump=n_frag/1.5 c n_jump=1 n_all=0 n_all_max=0 n_max_cut=0.7 n_cutt=0 d01=d0+1.5 if(d01.lt.d0_min)d01=d0_min d02=d01*d01 GLmaxB=0 GLmaxA=0 GLmaxA_cut=0.12 n_frag1=n_frag-1 n_frag2=n_frag/2 m1=nseq1-n_frag+1 m2=nseq2-n_frag+1 do ii=1,m1,n_jump do jj=1,m2,n_jump secmatch=0 coilcnt=0 do k=1,n_frag iii=ii+k-1 jjj=jj+k-1 if(isec(iii).eq.jsec(jjj))then secmatch=secmatch+1 endif if(isec(iii).eq.1)then !.or.isec(iii).ne.4)then coilcnt=coilcnt+1 endif enddo if(secmatch.gt.n_frag2.and.coilcnt.lt.n_frag2)then do j=1,nseq2 invmap(j)=-1 enddo do k=0,n_frag1 i=ii+k j=jj+k invmap(j)=i enddo call get_GL(GL) if(GL.gt.GLmaxA*GLmaxA_cut)then !GLmaxB=GL do k=1,n_frag iii=ii+k-1 jjj=jj+k-1 r_1(1,k)=xa(1,iii,0) r_1(2,k)=xa(2,iii,0) r_1(3,k)=xa(3,iii,0) r_2(1,k)=xa(1,jjj,1) r_2(2,k)=xa(2,jjj,1) r_2(3,k)=xa(3,jjj,1) enddo *********superpose the two structures and rotate it ***************** call u3b(w,r_1,r_2,n_frag,1,rms,u,t,ier) !u rotate r_1 to r_2 do i=1,nseq1 xx=t(1)+u(1,1)*xa(1,i,0)+u(1,2)*xa(2,i,0)+u(1,3)*xa(3,i,0) yy=t(2)+u(2,1)*xa(1,i,0)+u(2,2)*xa(2,i,0)+u(2,3)*xa(3,i,0) zz=t(3)+u(3,1)*xa(1,i,0)+u(3,2)*xa(2,i,0)+u(3,3)*xa(3,i,0) do j=1,nseq2 dd=(xx-xa(1,j,1))**2+(yy-xa(2,j,1))**2+(zz-xa(3,j,1))**2 score(i,j)=1/(1+dd/d02) enddo enddo *********extract alignement with score(i,j) ***************** align_len1=0 call DP(NSEQ1,NSEQ2) !print *,align_len1 !n_all=0 !do j=1,nseq2 ! if(invmap(j).gt.0)then ! n_all=n_all+1 ! endif !enddo n_cutt=n_all_max*n_max_cut !print *,n_all,n_cutt,n_all_max,n_max_cut if(align_len1.gt.n_cutt)then if(align_len1.gt.n_all_max)then n_all_max=align_len1 endif call get_GL(GL) if(GL.gt.GLmaxA)then GLmaxA=GL do j=1,nseq2 invmap_i(j)=invmap(j) enddo endif endif endif endif enddo enddo return end ************************************************************** * smooth the secondary structure assignment ************************************************************** subroutine smooth PARAMETER(nmax=5000) common/sec/isec(nmax),jsec(nmax) common/length/nseq1,nseq2 *** smooth single --------------> *** --x-- => ----- do i=1,nseq1 if(isec(i).eq.2.or.isec(i).eq.4)then j=isec(i) if(isec(i-2).ne.j)then if(isec(i-1).ne.j)then if(isec(i+1).ne.j)then if(isec(i+1).ne.j)then isec(i)=1 endif endif endif endif endif enddo do i=1,nseq2 if(jsec(i).eq.2.or.jsec(i).eq.4)then j=jsec(i) if(jsec(i-2).ne.j)then if(jsec(i-1).ne.j)then if(jsec(i+1).ne.j)then if(jsec(i+1).ne.j)then jsec(i)=1 endif endif endif endif endif enddo *** smooth double --------------> *** --xx-- => ------ do i=1,nseq1 if(isec(i).ne.2)then if(isec(i+1).ne.2)then if(isec(i+2).eq.2)then if(isec(i+3).eq.2)then if(isec(i+4).ne.2)then if(isec(i+5).ne.2)then isec(i+2)=1 isec(i+3)=1 endif endif endif endif endif endif if(isec(i).ne.4)then if(isec(i+1).ne.4)then if(isec(i+2).eq.4)then if(isec(i+3).eq.4)then if(isec(i+4).ne.4)then if(isec(i+5).ne.4)then isec(i+2)=1 isec(i+3)=1 endif endif endif endif endif endif enddo do i=1,nseq2 if(jsec(i).ne.2)then if(jsec(i+1).ne.2)then if(jsec(i+2).eq.2)then if(jsec(i+3).eq.2)then if(jsec(i+4).ne.2)then if(jsec(i+5).ne.2)then jsec(i+2)=1 jsec(i+3)=1 endif endif endif endif endif endif if(jsec(i).ne.4)then if(jsec(i+1).ne.4)then if(jsec(i+2).eq.4)then if(jsec(i+3).eq.4)then if(jsec(i+4).ne.4)then if(jsec(i+5).ne.4)then jsec(i+2)=1 jsec(i+3)=1 endif endif endif endif endif endif enddo *** connect --------------> *** x-x => xxx do i=1,nseq1 if(isec(i).eq.2)then if(isec(i+1).ne.2)then if(isec(i+2).eq.2)then isec(i+1)=2 endif endif endif if(isec(i).eq.4)then if(isec(i+1).ne.4)then if(isec(i+2).eq.4)then isec(i+1)=4 endif endif endif enddo do i=1,nseq2 if(jsec(i).eq.2)then if(jsec(i+1).ne.2)then if(jsec(i+2).eq.2)then jsec(i+1)=2 endif endif endif if(jsec(i).eq.4)then if(jsec(i+1).ne.4)then if(jsec(i+2).eq.4)then jsec(i+1)=4 endif endif endif enddo return end ************************************************************* * assign secondary structure: ************************************************************* function diszy(i,i1,i2) PARAMETER(nmax=5000) COMMON/BACKBONE/XA(3,nmax,0:1) diszy=sqrt((xa(1,i1,i)-xa(1,i2,i))**2 & +(xa(2,i1,i)-xa(2,i2,i))**2 & +(xa(3,i1,i)-xa(3,i2,i))**2) return end ************************************************************* * assign secondary structure: ************************************************************* function make_sec(dis13,dis14,dis15,dis24,dis25,dis35) make_sec=1 delta=2.1 if(abs(dis15-6.37).lt.delta)then if(abs(dis14-5.18).lt.delta)then if(abs(dis25-5.18).lt.delta)then if(abs(dis13-5.45).lt.delta)then if(abs(dis24-5.45).lt.delta)then if(abs(dis35-5.45).lt.delta)then make_sec=2 !helix return endif endif endif endif endif endif delta=1.42 if(abs(dis15-13).lt.delta)then if(abs(dis14-10.4).lt.delta)then if(abs(dis25-10.4).lt.delta)then if(abs(dis13-6.1).lt.delta)then if(abs(dis24-6.1).lt.delta)then if(abs(dis35-6.1).lt.delta)then make_sec=4 !strand return endif endif endif endif endif endif if(dis15.lt.8)then make_sec=3 endif return end **************************************************************** * quickly calculate TM-score with given invmap(i) in 3 iterations **************************************************************** subroutine get_GL(GL) PARAMETER(nmax=5000) common/length/nseq1,nseq2 COMMON/BACKBONE/XA(3,nmax,0:1) common/dpc/score(nmax,nmax),gap_open,invmap(nmax) common/zscore/zrms,n_al,rmsd_al common/d0/d0,anseq dimension xtm1(nmax),ytm1(nmax),ztm1(nmax) dimension xtm2(nmax),ytm2(nmax),ztm2(nmax) common/TM/TM,TMmax common/n1n2/n1(nmax),n2(nmax) common/d00/d00,d002 dimension xo1(nmax),yo1(nmax),zo1(nmax) dimension xo2(nmax),yo2(nmax),zo2(nmax) dimension dis2(nmax) 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 c calculate RMSD between aligned structures and rotate the structures --> n_al=0 do j=1,NSEQ2 i=invmap(j) !j aligned to i !print *,j,i if(i.gt.0)then n_al=n_al+1 r_1(1,n_al)=xa(1,i,0) r_1(2,n_al)=xa(2,i,0) r_1(3,n_al)=xa(3,i,0) r_2(1,n_al)=xa(1,j,1) r_2(2,n_al)=xa(2,j,1) r_2(3,n_al)=xa(3,j,1) xo1(n_al)=xa(1,i,0) yo1(n_al)=xa(2,i,0) zo1(n_al)=xa(3,i,0) xo2(n_al)=xa(1,j,1) yo2(n_al)=xa(2,j,1) zo2(n_al)=xa(3,j,1) endif enddo call u3b(w,r_1,r_2,n_al,1,rms,u,t,ier) !u rotate r_1 to r_2 GL=0 do i=1,n_al xx=t(1)+u(1,1)*xo1(i)+u(1,2)*yo1(i)+u(1,3)*zo1(i) yy=t(2)+u(2,1)*xo1(i)+u(2,2)*yo1(i)+u(2,3)*zo1(i) zz=t(3)+u(3,1)*xo1(i)+u(3,2)*yo1(i)+u(3,3)*zo1(i) dis2(i)=(xx-xo2(i))**2+(yy-yo2(i))**2+(zz-zo2(i))**2 GL=GL+1/(1+dis2(i)/(d0**2)) enddo ccc for next iteration-------------> d002t=d002 21 j=0 do i=1,n_al if(dis2(i).le.d002t)then j=j+1 r_1(1,j)=xo1(i) r_1(2,j)=yo1(i) r_1(3,j)=zo1(i) r_2(1,j)=xo2(i) r_2(2,j)=yo2(i) r_2(3,j)=zo2(i) endif enddo if(j.lt.3.and.n_al.gt.3)then d002t=d002t+.5 goto 21 endif L=j call u3b(w,r_1,r_2,L,1,rms,u,t,ier) !u rotate r_1 to r_2 G2=0 do i=1,n_al xx=t(1)+u(1,1)*xo1(i)+u(1,2)*yo1(i)+u(1,3)*zo1(i) yy=t(2)+u(2,1)*xo1(i)+u(2,2)*yo1(i)+u(2,3)*zo1(i) zz=t(3)+u(3,1)*xo1(i)+u(3,2)*yo1(i)+u(3,3)*zo1(i) dis2(i)=(xx-xo2(i))**2+(yy-yo2(i))**2+(zz-zo2(i))**2 G2=G2+1/(1+dis2(i)/(d0**2)) enddo ccc for next iteration-------------> d002t=d002+1 22 j=0 do i=1,n_al if(dis2(i).le.d002t)then j=j+1 r_1(1,j)=xo1(i) r_1(2,j)=yo1(i) r_1(3,j)=zo1(i) r_2(1,j)=xo2(i) r_2(2,j)=yo2(i) r_2(3,j)=zo2(i) endif enddo if(j.lt.3.and.n_al.gt.3)then d002t=d002t+.5 goto 22 endif L=j call u3b(w,r_1,r_2,L,1,rms,u,t,ier) !u rotate r_1 to r_2 G3=0 do i=1,n_al xx=t(1)+u(1,1)*xo1(i)+u(1,2)*yo1(i)+u(1,3)*zo1(i) yy=t(2)+u(2,1)*xo1(i)+u(2,2)*yo1(i)+u(2,3)*zo1(i) zz=t(3)+u(3,1)*xo1(i)+u(3,2)*yo1(i)+u(3,3)*zo1(i) dis2(i)=(xx-xo2(i))**2+(yy-yo2(i))**2+(zz-zo2(i))**2 G3=G3+1/(1+dis2(i)/(d0**2)) enddo if(G2.gt.GL)GL=G2 if(G3.gt.GL)GL=G3 c^^^^^^^^^^^^^^^^ GL done ^^^^^^^^^^^^^^^^^^^^^^^^^^^ return end **************************************************************** * quickly calculate TM-score with given invmap(i) in 3 iterations **************************************************************** subroutine get_GL11(GL) PARAMETER(nmax=5000) common/length1/nseq11,nseq22,nseq11_p,nseq22_p COMMON/BACKBONE/XA(3,nmax,0:1) common/dpc/score(nmax,nmax),gap_open,invmap(nmax) common/zscore/zrms,n_al,rmsd_al common/d0/d0,anseq dimension xtm1(nmax),ytm1(nmax),ztm1(nmax) dimension xtm2(nmax),ytm2(nmax),ztm2(nmax) common/TM/TM,TMmax common/n1n2/n1(nmax),n2(nmax) common/d00/d00,d002 dimension xo1(nmax),yo1(nmax),zo1(nmax) dimension xo2(nmax),yo2(nmax),zo2(nmax) dimension dis2(nmax) 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 c calculate RMSD between aligned structures and rotate the structures --> n_al=0 do j=nseq22_p,NSEQ22+nseq22_p i=invmap(j) !j aligned to i if(i.gt.0)then n_al=n_al+1 r_1(1,n_al)=xa(1,i,0) r_1(2,n_al)=xa(2,i,0) r_1(3,n_al)=xa(3,i,0) r_2(1,n_al)=xa(1,j,1) r_2(2,n_al)=xa(2,j,1) r_2(3,n_al)=xa(3,j,1) xo1(n_al)=xa(1,i,0) yo1(n_al)=xa(2,i,0) zo1(n_al)=xa(3,i,0) xo2(n_al)=xa(1,j,1) yo2(n_al)=xa(2,j,1) zo2(n_al)=xa(3,j,1) endif enddo !print *,n_al call u3b(w,r_1,r_2,n_al,1,rms,u,t,ier) !u rotate r_1 to r_2 GL=0 do i=1,n_al xx=t(1)+u(1,1)*xo1(i)+u(1,2)*yo1(i)+u(1,3)*zo1(i) yy=t(2)+u(2,1)*xo1(i)+u(2,2)*yo1(i)+u(2,3)*zo1(i) zz=t(3)+u(3,1)*xo1(i)+u(3,2)*yo1(i)+u(3,3)*zo1(i) dis2(i)=(xx-xo2(i))**2+(yy-yo2(i))**2+(zz-zo2(i))**2 GL=GL+1/(1+dis2(i)/(d0**2)) enddo ccc for next iteration-------------> d002t=d002 21 j=0 do i=1,n_al if(dis2(i).le.d002t)then j=j+1 r_1(1,j)=xo1(i) r_1(2,j)=yo1(i) r_1(3,j)=zo1(i) r_2(1,j)=xo2(i) r_2(2,j)=yo2(i) r_2(3,j)=zo2(i) endif enddo if(j.lt.3.and.n_al.gt.3)then d002t=d002t+.5 goto 21 endif L=j call u3b(w,r_1,r_2,L,1,rms,u,t,ier) !u rotate r_1 to r_2 G2=0 do i=1,n_al xx=t(1)+u(1,1)*xo1(i)+u(1,2)*yo1(i)+u(1,3)*zo1(i) yy=t(2)+u(2,1)*xo1(i)+u(2,2)*yo1(i)+u(2,3)*zo1(i) zz=t(3)+u(3,1)*xo1(i)+u(3,2)*yo1(i)+u(3,3)*zo1(i) dis2(i)=(xx-xo2(i))**2+(yy-yo2(i))**2+(zz-zo2(i))**2 G2=G2+1/(1+dis2(i)/(d0**2)) enddo ccc for next iteration-------------> d002t=d002+1 22 j=0 do i=1,n_al if(dis2(i).le.d002t)then j=j+1 r_1(1,j)=xo1(i) r_1(2,j)=yo1(i) r_1(3,j)=zo1(i) r_2(1,j)=xo2(i) r_2(2,j)=yo2(i) r_2(3,j)=zo2(i) endif enddo if(j.lt.3.and.n_al.gt.3)then d002t=d002t+.5 goto 22 endif L=j call u3b(w,r_1,r_2,L,1,rms,u,t,ier) !u rotate r_1 to r_2 G3=0 do i=1,n_al xx=t(1)+u(1,1)*xo1(i)+u(1,2)*yo1(i)+u(1,3)*zo1(i) yy=t(2)+u(2,1)*xo1(i)+u(2,2)*yo1(i)+u(2,3)*zo1(i) zz=t(3)+u(3,1)*xo1(i)+u(3,2)*yo1(i)+u(3,3)*zo1(i) dis2(i)=(xx-xo2(i))**2+(yy-yo2(i))**2+(zz-zo2(i))**2 G3=G3+1/(1+dis2(i)/(d0**2)) enddo if(G2.gt.GL)GL=G2 if(G3.gt.GL)GL=G3 !print *,n_al c^^^^^^^^^^^^^^^^ GL done ^^^^^^^^^^^^^^^^^^^^^^^^^^^ return end **************************************************************** * with invmap(i) calculate TM-score and martix score(i,j) for rotation **************************************************************** subroutine get_score PARAMETER(nmax=5000) common/length/nseq1,nseq2 COMMON/BACKBONE/XA(3,nmax,0:1) common/dpc/score(nmax,nmax),gap_open,invmap(nmax) common/zscore/zrms,n_al,rmsd_al common/d0/d0,anseq dimension xtm1(nmax),ytm1(nmax),ztm1(nmax) dimension xtm2(nmax),ytm2(nmax),ztm2(nmax) common/TM/TM,TMmax common/n1n2/n1(nmax),n2(nmax) 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 c calculate RMSD between aligned structures and rotate the structures --> n_al=0 do j=1,NSEQ2 i=invmap(j) !j aligned to i if(i.gt.0)then n_al=n_al+1 ccc for TM-score: xtm1(n_al)=xa(1,i,0) !for TM-score ytm1(n_al)=xa(2,i,0) ztm1(n_al)=xa(3,i,0) xtm2(n_al)=xa(1,j,1) ytm2(n_al)=xa(2,j,1) ztm2(n_al)=xa(3,j,1) ccc for rotation matrix: r_1(1,n_al)=xa(1,i,0) r_1(2,n_al)=xa(2,i,0) r_1(3,n_al)=xa(3,i,0) endif enddo *** calculate TM-score for the given alignment-----------> d0_input=d0 call TMscore8_search(d0_input,n_al,xtm1,ytm1,ztm1,n1, & n_al,xtm2,ytm2,ztm2,n2,TM,Rcomm,Lcomm) !simplified search engine TM=TM*n_al/anseq !TM-score *** calculate score matrix score(i,j)------------------> do i=1,n_al r_2(1,i)=xtm1(i) r_2(2,i)=ytm1(i) r_2(3,i)=ztm1(i) enddo call u3b(w,r_1,r_2,n_al,1,rms,u,t,ier) !u rotate r_1 to r_2 do i=1,nseq1 xx=t(1)+u(1,1)*xa(1,i,0)+u(1,2)*xa(2,i,0)+u(1,3)*xa(3,i,0) yy=t(2)+u(2,1)*xa(1,i,0)+u(2,2)*xa(2,i,0)+u(2,3)*xa(3,i,0) zz=t(3)+u(3,1)*xa(1,i,0)+u(3,2)*xa(2,i,0)+u(3,3)*xa(3,i,0) do j=1,nseq2 dd=(xx-xa(1,j,1))**2+(yy-xa(2,j,1))**2+(zz-xa(3,j,1))**2 score(i,j)=1/(1+dd/d0**2) enddo enddo c^^^^^^^^^^^^^^^^ score(i,j) done ^^^^^^^^^^^^^^^^^^^^^^^^^^^ return end **************************************************************** * with invmap(i) calculate TM-score and martix score(i,j) for rotation **************************************************************** subroutine get_score11 PARAMETER(nmax=5000) common/length1/nseq11,nseq22 COMMON/BACKBONE/XA(3,nmax,0:1) common/dpc/score(nmax,nmax),gap_open,invmap(nmax) common/zscore/zrms,n_al,rmsd_al common/d0/d0,anseq dimension xtm1(nmax),ytm1(nmax),ztm1(nmax) dimension xtm2(nmax),ytm2(nmax),ztm2(nmax) common/TM/TM,TMmax common/n1n2/n1(nmax),n2(nmax) 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 c calculate RMSD between aligned structures and rotate the structures --> n_al=0 do j=1,NSEQ22 i=invmap(j) !j aligned to i if(i.gt.0)then n_al=n_al+1 ccc for TM-score: xtm1(n_al)=xa(1,i,0) !for TM-score ytm1(n_al)=xa(2,i,0) ztm1(n_al)=xa(3,i,0) xtm2(n_al)=xa(1,j,1) ytm2(n_al)=xa(2,j,1) ztm2(n_al)=xa(3,j,1) ccc for rotation matrix: r_1(1,n_al)=xa(1,i,0) r_1(2,n_al)=xa(2,i,0) r_1(3,n_al)=xa(3,i,0) endif enddo *** calculate TM-score for the given alignment-----------> d0_input=d0 call TMscore8_search(d0_input,n_al,xtm1,ytm1,ztm1,n1, & n_al,xtm2,ytm2,ztm2,n2,TM,Rcomm,Lcomm) !simplified search engine TM=TM*n_al/nseq22 !TM-score *** calculate score matrix score(i,j)------------------> do i=1,n_al r_2(1,i)=xtm1(i) r_2(2,i)=ytm1(i) r_2(3,i)=ztm1(i) enddo call u3b(w,r_1,r_2,n_al,1,rms,u,t,ier) !u rotate r_1 to r_2 do i=1,nseq11 xx=t(1)+u(1,1)*xa(1,i,0)+u(1,2)*xa(2,i,0)+u(1,3)*xa(3,i,0) yy=t(2)+u(2,1)*xa(1,i,0)+u(2,2)*xa(2,i,0)+u(2,3)*xa(3,i,0) zz=t(3)+u(3,1)*xa(1,i,0)+u(3,2)*xa(2,i,0)+u(3,3)*xa(3,i,0) do j=1,nseq22 dd=(xx-xa(1,j,1))**2+(yy-xa(2,j,1))**2+(zz-xa(3,j,1))**2 score(i,j)=1/(1+dd/d0**2) enddo enddo c^^^^^^^^^^^^^^^^ score(i,j) done ^^^^^^^^^^^^^^^^^^^^^^^^^^^ return end **************************************************************** * with invmap(i) calculate score(i,j) using RMSD rotation **************************************************************** subroutine get_score1 PARAMETER(nmax=5000) common/length/nseq1,nseq2 COMMON/BACKBONE/XA(3,nmax,0:1) common/dpc/score(nmax,nmax),gap_open,invmap(nmax) common/zscore/zrms,n_al,rmsd_al common/d0/d0,anseq common/d0min/d0_min dimension xtm1(nmax),ytm1(nmax),ztm1(nmax) dimension xtm2(nmax),ytm2(nmax),ztm2(nmax) common/TM/TM,TMmax common/n1n2/n1(nmax),n2(nmax) 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 c calculate RMSD between aligned structures and rotate the structures --> n_al=0 do j=1,NSEQ2 i=invmap(j) !j aligned to i if(i.gt.0)then n_al=n_al+1 ccc for rotation matrix: r_1(1,n_al)=xa(1,i,0) r_1(2,n_al)=xa(2,i,0) r_1(3,n_al)=xa(3,i,0) r_2(1,n_al)=xa(1,j,1) r_2(2,n_al)=xa(2,j,1) r_2(3,n_al)=xa(3,j,1) endif enddo *** calculate score matrix score(i,j)------------------> call u3b(w,r_1,r_2,n_al,1,rms,u,t,ier) !u rotate r_1 to r_2 d01=d0+1.5 if(d01.lt.d0_min)d01=d0_min d02=d01*d01 do i=1,nseq1 xx=t(1)+u(1,1)*xa(1,i,0)+u(1,2)*xa(2,i,0)+u(1,3)*xa(3,i,0) yy=t(2)+u(2,1)*xa(1,i,0)+u(2,2)*xa(2,i,0)+u(2,3)*xa(3,i,0) zz=t(3)+u(3,1)*xa(1,i,0)+u(3,2)*xa(2,i,0)+u(3,3)*xa(3,i,0) do j=1,nseq2 dd=(xx-xa(1,j,1))**2+(yy-xa(2,j,1))**2+(zz-xa(3,j,1))**2 score(i,j)=1/(1+dd/d02) enddo enddo c^^^^^^^^^^^^^^^^ score(i,j) done ^^^^^^^^^^^^^^^^^^^^^^^^^^^ return end ************************************************************************* ************************************************************************* * This is a subroutine to compare two structures and find the * superposition that has the maximum TM-score. * * L1--Length of the first structure * (x1(i),y1(i),z1(i))--coordinates of i'th residue at the first structure * n1(i)--Residue sequence number of i'th residue at the first structure * L2--Length of the second structure * (x2(i),y2(i),z2(i))--coordinates of i'th residue at the second structure * n2(i)--Residue sequence number of i'th residue at the second structure * TM--TM-score of the comparison * Rcomm--RMSD of two structures in the common aligned residues * Lcomm--Length of the common aligned regions * * Note: * 1, Always put native as the second structure, by which TM-score * is normalized. * 2, The returned (x1(i),y1(i),z1(i)) are the rotated structure after * TM-score superposition. ************************************************************************* ************************************************************************* *** dis<8, simplified search engine subroutine TMscore8_search(dx,L1,x1,y1,z1,n1,L2,x2,y2,z2,n2, & TM,Rcomm,Lcomm) 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 common/para/d,d0 common/d0min/d0_min common/align/n_ali,iA(nmax),iB(nmax) common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score dimension k_ali(nmax),k_ali0(nmax) dimension L_ini(100),iq(nmax) common/scores/score double precision score,score_max dimension xa(nmax),ya(nmax),za(nmax) dimension iL0(nmax) common/flag/intflag dimension x1(nmax),y1(nmax),z1(nmax),n1(nmax) dimension x2(nmax),y2(nmax),z2(nmax),n2(nmax) ccc RMSD:arrays 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 ********* convert input data *************************** * because L1=L2 in this special case----------> nseqA=L1 nseqB=L2 do i=1,nseqA xa(i)=x1(i) ya(i)=y1(i) za(i)=z1(i) nresA(i)=n1(i) xb(i)=x2(i) yb(i)=y2(i) zb(i)=z2(i) nresB(i)=n2(i) iA(i)=i iB(i)=i enddo n_ali=L1 !number of aligned residues Lcomm=L1 ************///// * parameters: ***************** *** d0-------------> d0=dx if(d0.lt.d0_min)d0=d0_min *** d0_search -----> d0_search=d0 if(d0_search.gt.8)d0_search=8 if(d0_search.lt.4.5)d0_search=4.5 *** iterative parameters -----> n_it=20 !maximum number of iterations d_output=5 !for output alignment n_init_max=6 !maximum number of L_init n_init=0 L_ini_min=4 if(n_ali.lt.4)L_ini_min=n_ali do i=1,n_init_max-1 n_init=n_init+1 L_ini(n_init)=n_ali/2**(n_init-1) if(L_ini(n_init).le.L_ini_min)then L_ini(n_init)=L_ini_min goto 402 endif enddo n_init=n_init+1 L_ini(n_init)=L_ini_min 402 continue ****************************************************************** * find the maximum scarrays ore starting from local structures superposition ****************************************************************** score_max=-1 !TM-score do 333 i_init=1,n_init L_init=L_ini(i_init) iL_max=n_ali-L_init+1 k=0 do i=1,iL_max,40 !this is the simplification! k=k+1 iL0(k)=i enddo if(iL0(k).lt.iL_max)then k=k+1 iL0(k)=iL_max endif n_shift=k do 300 i_shift=1,n_shift iL=iL0(i_shift) LL=0 ka=0 do i=1,L_init k=iL+i-1 ![1,n_ali] common aligned r_1(1,i)=xa(iA(k)) r_1(2,i)=ya(iA(k)) r_1(3,i)=za(iA(k)) r_2(1,i)=xb(iB(k)) r_2(2,i)=yb(iB(k)) r_2(3,i)=zb(iB(k)) LL=LL+1 ka=ka+1 k_ali(ka)=k enddo call u3b(w,r_1,r_2,LL,2,rms,u,t,ier) !u rotate r_1 to r_2 if(i_init.eq.1)then !global superposition armsd=dsqrt(rms/LL) Rcomm=armsd endif do j=1,nseqA xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j) yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j) zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j) enddo d=d0_search-1 call score_fun8 !init, get scores, n_cut+i_ali(i) for iteration if(score_max.lt.score)then score_max=score ka0=ka do i=1,ka0 k_ali0(i)=k_ali(i) enddo endif *** iteration for extending ----------------------------------> d=d0_search+1 do 301 it=1,n_it LL=0 ka=0 do i=1,n_cut m=i_ali(i) ![1,n_ali] 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)) ka=ka+1 k_ali(ka)=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 xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j) yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j) zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j) enddo call score_fun8 !get scores, n_cut+i_ali(i) for iteration if(score_max.lt.score)then score_max=score ka0=ka do i=1,ka k_ali0(i)=k_ali(i) enddo endif if(it.eq.n_it)goto 302 if(n_cut.eq.ka)then neq=0 do i=1,n_cut if(i_ali(i).eq.k_ali(i))neq=neq+1 enddo if(n_cut.eq.neq)goto 302 endif 301 continue !for iteration 302 continue 300 continue !for shift 333 continue !for initial length, L_ali/M ******** return the final rotation **************** 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 TM=score_max c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ return END ************************************************************************* ************************************************************************* * This is a subroutine to compare two structures and find the * superposition that has the maximum TM-score. * * L1--Length of the first structure * (x1(i),y1(i),z1(i))--coordinates of i'th residue at the first structure * n1(i)--Residue sequence number of i'th residue at the first structure * L2--Length of the second structure * (x2(i),y2(i),z2(i))--coordinates of i'th residue at the second structure * n2(i)--Residue sequence number of i'th residue at the second structure * TM--TM-score of the comparison * Rcomm--RMSD of two structures in the common aligned residues * Lcomm--Length of the common aligned regions * * Note: * 1, Always put native as the second structure, by which TM-score * is normalized. * 2, The returned (x1(i),y1(i),z1(i)) are the rotated structure after * TM-score superposition. ************************************************************************* ************************************************************************* *** dis<8, but same search engine subroutine TMscore8(dx,L1,x1,y1,z1,n1,L2,x2,y2,z2,n2, & TM,Rcomm,Lcomm) 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 common/para/d,d0 common/d0min/d0_min common/align/n_ali,iA(nmax),iB(nmax) common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score dimension k_ali(nmax),k_ali0(nmax) dimension L_ini(100),iq(nmax) common/scores/score double precision score,score_max dimension xa(nmax),ya(nmax),za(nmax) common/flag/intflag dimension x1(nmax),y1(nmax),z1(nmax),n1(nmax) dimension x2(nmax),y2(nmax),z2(nmax),n2(nmax) 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 ********* convert input data *************************** * because L1=L2 in this special case----------> nseqA=L1 nseqB=L2 do i=1,nseqA xa(i)=x1(i) ya(i)=y1(i) za(i)=z1(i) nresA(i)=n1(i) xb(i)=x2(i) yb(i)=y2(i) zb(i)=z2(i) nresB(i)=n2(i) iA(i)=i iB(i)=i enddo n_ali=L1 !number of aligned residues Lcomm=L1 ************///// * parameters: ***************** *** d0-------------> d0=dx if(d0.lt.d0_min)d0=d0_min *** d0_search -----> d0_search=d0 if(d0_search.gt.8)d0_search=8 if(d0_search.lt.4.5)d0_search=4.5 *** iterative parameters -----> n_it=20 !maximum number of iterations d_output=5 !for output alignment n_init_max=6 !maximum number of L_init n_init=0 L_ini_min=4 if(n_ali.lt.4)L_ini_min=n_ali do i=1,n_init_max-1 n_init=n_init+1 L_ini(n_init)=n_ali/2**(n_init-1) if(L_ini(n_init).le.L_ini_min)then L_ini(n_init)=L_ini_min goto 402 endif enddo n_init=n_init+1 L_ini(n_init)=L_ini_min 402 continue ****************************************************************** * find the maximum score starting from local structures superposition ****************************************************************** score_max=-1 !TM-score do 333 i_init=1,n_init L_init=L_ini(i_init) iL_max=n_ali-L_init+1 do 300 iL=1,iL_max !on aligned residues, [1,nsene.jqA] LL=0 ka=0 do i=1,L_init k=iL+i-1 ![1,n_ali] common aligned r_1(1,i)=xa(iA(k)) r_1(2,i)=ya(iA(k)) r_1(3,i)=za(iA(k)) r_2(1,i)=xb(iB(k)) r_2(2,i)=yb(iB(k)) r_2(3,i)=zb(iB(k)) LL=LL+1 ka=ka+1 k_ali(ka)=k enddo call u3b(w,r_1,r_2,LL,2,rms,u,t,ier) !u rotate r_1 to r_2 if(i_init.eq.1)then !global superposition armsd=dsqrt(rms/LL) Rcomm=armsd endif do j=1,nseqA xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j) yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j) zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j) enddo d=d0_search-1 call score_fun8 !init, get scores, n_cut+i_ali(i) for iteration if(score_max.lt.score)then score_max=score ka0=ka do i=1,ka0 k_ali0(i)=k_ali(i) enddo endif *** iteration for extending ----------------------------------> d=d0_search+1 do 301 it=1,n_it LL=0 ka=0 do i=1,n_cut m=i_ali(i) ![1,n_ali] 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)) ka=ka+1 k_ali(ka)=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 xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j) yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j) zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j) enddo call score_fun8 !get scores, n_cut+i_ali(i) for iteration if(score_max.lt.score)then score_max=score ka0=ka do i=1,ka k_ali0(i)=k_ali(i) enddo endif if(it.eq.n_it)goto 302 if(n_cut.eq.ka)then neq=0 do i=1,n_cut if(i_ali(i).eq.k_ali(i))neq=neq+1 enddo if(n_cut.eq.neq)goto 302 endif 301 continue !for iteration 302 continue 300 continue !for shift 333 continue !for initial length, L_ali/M ******** return the final rotation **************** 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 TM=score_max c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ return END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c 1, collect those residues with dis nseqA=L1 nseqB=L2 do i=1,nseqA xa(i)=x1(i) ya(i)=y1(i) za(i)=z1(i) nresA(i)=n1(i) xb(i)=x2(i) yb(i)=y2(i) zb(i)=z2(i) nresB(i)=n2(i) iA(i)=i iB(i)=i enddo n_ali=L1 !number of aligned residues Lcomm=L1 ************///// * parameters: ***************** *** d0-------------> c d0=1.24*(nseqB-15)**(1.0/3.0)-1.8 d0=dx if(d0.lt.d0_min)d0=d0_min *** d0_search -----> d0_search=d0 if(d0_search.gt.8)d0_search=8 if(d0_search.lt.4.5)d0_search=4.5 *** iterative parameters -----> n_it=20 !maximum number of iterations d_output=5 !for output alignment n_init_max=6 !maximum number of L_init n_init=0 L_ini_min=4 if(n_ali.lt.4)L_ini_min=n_ali do i=1,n_init_max-1 n_init=n_init+1 L_ini(n_init)=n_ali/2**(n_init-1) if(L_ini(n_init).le.L_ini_min)then L_ini(n_init)=L_ini_min goto 402 endif enddo n_init=n_init+1 L_ini(n_init)=L_ini_min 402 continue ****************************************************************** * find the maximum score starting from local structures superposition ****************************************************************** score_max=-1 !TM-score do 333 i_init=1,n_init L_init=L_ini(i_init) iL_max=n_ali-L_init+1 do 300 iL=1,iL_max !on aligned residues, [1,nseqA] LL=0 ka=0 do i=1,L_init k=iL+i-1 ![1,n_ali] common aligned r_1(1,i)=xa(iA(k)) r_1(2,i)=ya(iA(k)) r_1(3,i)=za(iA(k)) r_2(1,i)=xb(iB(k)) r_2(2,i)=yb(iB(k)) r_2(3,i)=zb(iB(k)) LL=LL+1 ka=ka+1 k_ali(ka)=k enddo call u3b(w,r_1,r_2,LL,2,rms,u,t,ier) !u rotate r_1 to r_2 if(i_init.eq.1)then !global superposition armsd=dsqrt(rms/LL) Rcomm=armsd endif do j=1,nseqA xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j) yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j) zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j) enddo d=d0_search-1 call score_fun !init, get scores, n_cut+i_ali(i) for iter if(score_max.lt.score)then score_max=score ka0=ka do i=1,ka0 k_ali0(i)=k_ali(i) enddo endif *** iteration for extending ----------------------------------> d=d0_search+1 do 301 it=1,n_it LL=0 ka=0 do i=1,n_cut m=i_ali(i) ![1,n_ali] 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)) ka=ka+1 k_ali(ka)=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 xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j) yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j) zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j) enddo call score_fun !get scores, n_cut+i_ali(i) for iteration if(score_max.lt.score)then score_max=score ka0=ka do i=1,ka k_ali0(i)=k_ali(i) enddo endif if(it.eq.n_it)goto 302 if(n_cut.eq.ka)then neq=0 do i=1,n_cut if(i_ali(i).eq.k_ali(i))neq=neq+1 enddo if(n_cut.eq.neq)goto 302 endif 301 continue !for iteration 302 continue 300 continue !for shift 333 continue !for initial length, L_ali/M ******** return the final rotation **************** 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 TM=score_max c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ return END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c 1, collect those residues with dis