X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=binaries%2Fsrc%2Ftcoffee%2Ft_coffee_source%2FTMalign.f;fp=binaries%2Fsrc%2Ftcoffee%2Ft_coffee_source%2FTMalign.f;h=0000000000000000000000000000000000000000;hb=8140043fefd2101326b3e651e3b7e2d4ac806b99;hp=5920db11a00f26591360e58932ca0164b0aec48c;hpb=2eb05b9319960f2346045e5e2119d112c5909490;p=jabaws.git diff --git a/binaries/src/tcoffee/t_coffee_source/TMalign.f b/binaries/src/tcoffee/t_coffee_source/TMalign.f deleted file mode 100644 index 5920db1..0000000 --- a/binaries/src/tcoffee/t_coffee_source/TMalign.f +++ /dev/null @@ -1,2554 +0,0 @@ -************************************************************************* -* This program is to identify the best alignment of two protein -* structures to give the best TM-score. By default, TM-score is -* normalized by the second protein. The program can be freely -* copied or modified or redistributed. -* (For comments, please email to: yzhang@ku.edu) -* -* Reference: -* Yang Zhang, Jeffrey Skolnick, Nucl. Acid Res. 2005 33: 2303-9 -* -************************ updating history ******************************* -* 2005/06/01: A small bug of two-point superposition was fixed. -* 2005/10/19: the program was reformed so that the alignment -* results are not dependent on the specific compilers. -* 2006/06/20: select 'A' if there is altLoc when reading PDB file. -* 2007/02/27: rotation matrix from Chain-1 to Chain-2 is added. -* 2007/04/18: add options with TM-score normalized by average -* length, shorter length, or longer length of two -* structures. -* 2007/05/23: add additional output file 'TM.sup_all' for showing -* all atoms while 'TM.sup' is only for aligned atoms -* 2007/09/19: add a new feature alignment to deal with the problem -* of aligning fractional structures (e.g. protein -* interfaces). -************************************************************************* - - 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/d0/d0,anseq - common/d0min/d0_min - common/d00/d00,d002 - - character*100 fnam,pdb(100),outname - 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) - character aseq1(nmax2),aseq2(nmax2),aseq3(nmax2) - - dimension xx(nmax),yy(nmax),zz(nmax) - dimension m1(nmax),m2(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/n1n2/n1(nmax),n2(nmax) - common/d8/d8 - common/initial4/mm1(nmax),mm2(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 - - 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'/ - - call getarg(1,fnam) - if(fnam.eq.' '.or.fnam.eq.'?'.or.fnam.eq.'-h')then - write(*,*) - write(*,*)'Brief instruction for running TM-align program:' - write(*,*)'(For detail: Zhang & Skolnick, Nucl. Acid Res.', - & '2005 33, 2303)' - write(*,*) - write(*,*)'1. Align ''structure.pdb'' to ''target.pdb''' - write(*,*)' (By default, TM-score is normalized by the ', - & 'length of ''target.pdb'')' - write(*,*)' >TMalign structure.pdb target.pdb' - write(*,*) - write(*,*)'2. Run TM-align and output the superposition ', - & 'to ''TM.sup'' and ''TM.sup_all'':' - write(*,*)' >TMalign structure.pdb target.pdb -o TM.sup' - write(*,*)' To view the superimposed structures of the', - & ' aligned regions by rasmol:' - write(*,*)' >rasmol -script TM.sup)' - write(*,*)' To view the superimposed structures of all', - & ' regions by rasmol:' - write(*,*)' >rasmol -script TM.sup_all)' - write(*,*) - write(*,*)'3. If you want TM-score normalized by ', - & 'an assigned length, e.g. 100 aa:' - write(*,*)' >TMalign structure.pdb target.pdb -L 100' - write(*,*)' If you want TM-score normalized by the ', - & 'average length of two structures:' - write(*,*)' >TMalign structure.pdb target.pdb -a' - write(*,*)' If you want TM-score normalized by the ', - & 'shorter length of two structures:' - write(*,*)' >TMalign structure.pdb target.pdb -b' - write(*,*)' If you want TM-score normalized by the ', - & 'longer length of two structures:' - write(*,*)' >TMalign structure.pdb target.pdb -c' - write(*,*) -c write(*,*)'5. If you want to set a minimum cutoff for d0, ', -c & 'e.g. d0>3.0' -c write(*,*)' (By default d0>0.5, this option need ', -c & 'be considered only when L<35 aa)' -c write(*,*)' >TMalign structure.pdb target.pdb -dmin 3.0' -c write(*,*) - write(*,*)'(All above options does 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 - else - j=j+1 - pdb(j)=fnam - endif - if(i.lt.narg)goto 115 - -ccccccccc read data from first CA file: - open(unit=10,file=pdb(1),status='old') - i=0 - 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(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))seq1(i)=slc(j) - enddo - ss1(i)=aanam - if(i.ge.nmax)goto 1010 - endif - endif - endif - 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 - 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(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))seq2(i)=slc(j) - enddo - ss2(i)=aanam - if(i.ge.nmax)goto 1011 - endif - endif - endif - enddo - 1011 continue - close(10) - nseq2=i -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 - -***** do alignment ************************** - CALL super_align !to find invmap(j) - -************************************************************ -*** resuperpose to find residues of dis - n_al=0 - do j=1,nseq2 - if(invmap0(j).gt.0)then - i=invmap0(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 - 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 - -*** 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 - seq_id=float(n_eq)/(n_al+0.00000001) - n8_al=j - 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 - TM8=TM8*n8_al/anseq !TM-score after cutoff - -********* for output summary ****************************** - write(*,*) - write(*,*)'*****************************************************', - & '*********************' - write(*,*)'* TM-align ', - & ' *' - write(*,*)'* A protein structural alignment algorithm based on T', - & 'M-score *' - write(*,*)'* Reference: Y. Zhang and J. Skolnick, Nucl. Acids Re', - & 's. 2005 33, 2302-9 *' - write(*,*)'* Comments on the program, please email to: yzhang@ku', - & '.edu *' - write(*,*)'*****************************************************', - & '*********************' - write(*,*) - write(*,101)pdb(1),nseq1 - 101 format('Chain 1:',A10,' Size=',I4) - write(*,102)pdb(2),nseq2,int(anseq) - 102 format('Chain 2:',A10,' Size=',I4, - & ' (TM-score is normalized by ',I4,')') - write(*,*) - 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(*,*) - -********* extract rotation matrix ------------> - L=0 - do i=1,n8_al - k=m1(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)=xtm1(i) - r_2(2,L)=ytm1(i) - r_2(3,L)=ztm1(i) - enddo - if(L.gt.3)then - call u3b(w,r_1,r_2,L,1,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 -c do i=1,nseq1 -c ax=t(1)+u(1,1)*xa(1,i,0)+u(1,2)*xa(2,i,0)+u(1,3)*xa(3,i,0) -c ay=t(2)+u(2,1)*xa(1,i,0)+u(2,2)*xa(2,i,0)+u(2,3)*xa(3,i,0) -c az=t(3)+u(3,1)*xa(1,i,0)+u(3,2)*xa(2,i,0)+u(3,3)*xa(3,i,0) -c 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<2000' - write(7,900)'wireframe .45' - write(7,900)'select none' - write(7,900)'select atomno>2000' - write(7,900)'wireframe .20' - 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.5)then - write(7,901)m1(i) - write(7,900)'color red' - write(7,901)2000+m2(i) - write(7,900)'color red' - 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: - do i=1,n8_al - write(7,1237)m1(i),ss1(m1(i)),mm1(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)2000+m2(i),ss2(m2(i)),mm2(m2(i)), - $ xtm2(i),ytm2(i),ztm2(i) - enddo - write(7,1238) - do i=2,n8_al - write(7,1239)2000+m2(i-1),2000+m2(i) - 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<2000' - write(8,900)'wireframe .45' - write(8,900)'select none' - write(8,900)'select atomno>2000' - 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.5)then - write(8,901)m1(i) - write(8,900)'color red' - write(8,901)2000+m2(i) - write(8,900)'color red' - 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: - 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 - write(8,1239)i-1,i - enddo -*** chain2: - do i=1,nseq2 - write(8,1237)2000+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 - write(8,1239)2000+i-1,2000+i - enddo - close(8) - endif -*^^^^^^^^^^^^^^^^^^ output finished ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -************ output aligned sequences ************************** - ii=0 - i1_old=1 - i2_old=1 - do i=1,n8_al - do j=i1_old,m1(i)-1 - ii=ii+1 - aseq1(ii)=seq1(j) - aseq2(ii)='-' - aseq3(ii)=' ' - enddo - do j=i2_old,m2(i)-1 - ii=ii+1 - aseq1(ii)='-' - aseq2(ii)=seq2(j) - aseq3(ii)=' ' - enddo - ii=ii+1 - aseq1(ii)=seq1(m1(i)) - aseq2(ii)=seq2(m2(i)) - dis2=sqrt((xtm1(i)-xtm2(i))**2+ - & (ytm1(i)-ytm2(i))**2+(ztm1(i)-ztm2(i))**2) - if(dis2.le.5)then - aseq3(ii)=':' - else - aseq3(ii)='.' - endif - i1_old=m1(i)+1 - i2_old=m2(i)+1 - enddo - do i=i1_old,nseq1 - ii=ii+1 - aseq1(ii)=seq1(i) - aseq2(ii)='-' - aseq3(ii)=' ' - enddo - do i=i2_old,nseq2 - ii=ii+1 - aseq1(ii)='-' - aseq2(ii)=seq2(i) - aseq3(ii)=' ' - enddo - 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(*,*) - -c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - 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) - dimension gapp(100) - - TMmax=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) -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) - if(TM.gt.TMmax)then - TMmax=TM - do j=1,nseq2 - invmap0(j)=invmap(j) - enddo - endif - -***************************************************************** -* 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) -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 - -*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) - if(TM.gt.TMmax)then - TMmax=TM - do j=1,nseq2 - invmap0(j)=invmap(j) - enddo - endif - -***************************************************************** -* 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) -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 - -*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 - 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 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) -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 - -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) - - aL=min(nseq1,nseq2) - idel=aL/2.5 !minimum size of considered fragment - if(idel.le.5)idel=5 - n1=-nseq2+idel - n2=nseq1-idel - 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)=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 - - 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 - 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 - if(dis.gt.dcu_min)then - contin=.true. - endif - endif - elseif(mm(k,i).eq.(mm(k,i-1)+1))then - if(dis.lt.dcu)then - if(dis.gt.dcu_min)then - contin=.true. - endif - 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 - -************************************************************** -* 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 - 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 - -**************************************************************** -* 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 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) - - 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 - 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,1,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) - - 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,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,1,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,1,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 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_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