1 *************************************************************************
\r
2 * This program is to identify the best alignment of two protein
\r
3 * structures to give the best TM-score. By default, TM-score is
\r
4 * normalized by the second protein. The program can be freely
\r
5 * copied or modified or redistributed.
\r
6 * (For comments, please email to: yzhang@ku.edu)
\r
9 * Yang Zhang, Jeffrey Skolnick, Nucl. Acid Res. 2005 33: 2303-9
\r
11 ************************ updating history *******************************
\r
12 * 2005/06/01: A small bug of two-point superposition was fixed.
\r
13 * 2005/10/19: the program was reformed so that the alignment
\r
14 * results are not dependent on the specific compilers.
\r
15 * 2006/06/20: select 'A' if there is altLoc when reading PDB file.
\r
16 * 2007/02/27: rotation matrix from Chain-1 to Chain-2 is added.
\r
17 * 2007/04/18: add options with TM-score normalized by average
\r
18 * length, shorter length, or longer length of two
\r
20 * 2007/05/23: add additional output file 'TM.sup_all' for showing
\r
21 * all atoms while 'TM.sup' is only for aligned atoms
\r
22 * 2007/09/19: add a new feature alignment to deal with the problem
\r
23 * of aligning fractional structures (e.g. protein
\r
25 *************************************************************************
\r
28 PARAMETER(nmax=5000)
\r
29 PARAMETER(nmax2=10000)
\r
31 COMMON/BACKBONE/XA(3,nmax,0:1)
\r
32 common/dpc/score(nmax,nmax),gap_open,invmap(nmax)
\r
33 common/alignrst/invmap0(nmax)
\r
34 common/length/nseq1,nseq2
\r
39 character*100 fnam,pdb(100),outname
\r
40 character*3 aa(-1:20),aanam,ss1(nmax),ss2(nmax)
\r
42 character*200 outnameall_tmp,outnameall
\r
43 character seq1(0:nmax),seq2(0:nmax)
\r
44 character aseq1(nmax2),aseq2(nmax2),aseq3(nmax2)
\r
46 dimension xx(nmax),yy(nmax),zz(nmax)
\r
47 dimension m1(nmax),m2(nmax)
\r
48 dimension xtm1(nmax),ytm1(nmax),ztm1(nmax)
\r
49 dimension xtm2(nmax),ytm2(nmax),ztm2(nmax)
\r
50 common/init/invmap_i(nmax)
\r
53 common/n1n2/n1(nmax),n2(nmax)
\r
55 common/initial4/mm1(nmax),mm2(nmax)
\r
58 double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)
\r
59 double precision u(3,3),t(3),rms,drms !armsd is real
\r
63 data aa/ 'BCK','GLY','ALA','SER','CYS','VAL','THR','ILE',
\r
64 & 'PRO','MET','ASP','ASN','LEU',
\r
65 & 'LYS','GLU','GLN','ARG',
\r
66 & 'HIS','PHE','TYR','TRP','CYX'/
\r
67 character*1 slc(-1:20)
\r
68 data slc/'X','G','A','S','C','V','T','I',
\r
69 & 'P','M','D','N','L','K','E','Q','R',
\r
70 & 'H','F','Y','W','C'/
\r
73 if(fnam.eq.' '.or.fnam.eq.'?'.or.fnam.eq.'-h')then
\r
75 write(*,*)'Brief instruction for running TM-align program:'
\r
76 write(*,*)'(For detail: Zhang & Skolnick, Nucl. Acid Res.',
\r
79 write(*,*)'1. Align ''structure.pdb'' to ''target.pdb'''
\r
80 write(*,*)' (By default, TM-score is normalized by the ',
\r
81 & 'length of ''target.pdb'')'
\r
82 write(*,*)' >TMalign structure.pdb target.pdb'
\r
84 write(*,*)'2. Run TM-align and output the superposition ',
\r
85 & 'to ''TM.sup'' and ''TM.sup_all'':'
\r
86 write(*,*)' >TMalign structure.pdb target.pdb -o TM.sup'
\r
87 write(*,*)' To view the superimposed structures of the',
\r
88 & ' aligned regions by rasmol:'
\r
89 write(*,*)' >rasmol -script TM.sup)'
\r
90 write(*,*)' To view the superimposed structures of all',
\r
91 & ' regions by rasmol:'
\r
92 write(*,*)' >rasmol -script TM.sup_all)'
\r
94 write(*,*)'3. If you want TM-score normalized by ',
\r
95 & 'an assigned length, e.g. 100 aa:'
\r
96 write(*,*)' >TMalign structure.pdb target.pdb -L 100'
\r
97 write(*,*)' If you want TM-score normalized by the ',
\r
98 & 'average length of two structures:'
\r
99 write(*,*)' >TMalign structure.pdb target.pdb -a'
\r
100 write(*,*)' If you want TM-score normalized by the ',
\r
101 & 'shorter length of two structures:'
\r
102 write(*,*)' >TMalign structure.pdb target.pdb -b'
\r
103 write(*,*)' If you want TM-score normalized by the ',
\r
104 & 'longer length of two structures:'
\r
105 write(*,*)' >TMalign structure.pdb target.pdb -c'
\r
107 c write(*,*)'5. If you want to set a minimum cutoff for d0, ',
\r
109 c write(*,*)' (By default d0>0.5, this option need ',
\r
110 c & 'be considered only when L<35 aa)'
\r
111 c write(*,*)' >TMalign structure.pdb target.pdb -dmin 3.0'
\r
113 write(*,*)'(All above options does not change the ',
\r
114 & 'final structure alignment result)'
\r
119 ******* options ----------->
\r
120 m_out=-1 !decided output
\r
121 m_fix=-1 !fixed length-scale only for output
\r
122 m_ave=-1 !using average length
\r
123 m_d0_min=-1 !diminum d0 for search
\r
124 m_d0=-1 !given d0 for both search and output
\r
130 call getarg(i,fnam)
\r
131 if(fnam.eq.'-o')then
\r
134 call getarg(i,outname)
\r
135 elseif(fnam.eq.'-L')then !change both L_all and d0
\r
138 call getarg(i,fnam)
\r
140 elseif(fnam.eq.'-dmin')then
\r
143 call getarg(i,fnam)
\r
144 read(fnam,*)d0_min_input
\r
145 elseif(fnam.eq.'-d0')then
\r
148 call getarg(i,fnam)
\r
150 elseif(fnam.eq.'-a')then ! this will change the superposed output but not the alignment
\r
153 elseif(fnam.eq.'-b')then
\r
156 elseif(fnam.eq.'-c')then
\r
163 if(i.lt.narg)goto 115
\r
165 ccccccccc read data from first CA file:
\r
166 open(unit=10,file=pdb(1),status='old')
\r
169 read(10,9001,end=1010) s
\r
170 if(i.gt.0.and.s(1:3).eq.'TER')goto 1010
\r
171 if(s(1:3).eq.'ATO')then
\r
172 if(s(13:16).eq.'CA '.or.s(13:16).eq.' CA '.or
\r
173 & .s(13:16).eq.' CA')then
\r
174 if(s(17:17).eq.' '.or.s(17:17).eq.'A')then
\r
176 read(s,9000)du,aanam,du,mm1(i),du,
\r
177 $ xa(1,i,0),xa(2,i,0),xa(3,i,0)
\r
179 if(aanam.eq.aa(j))seq1(i)=slc(j)
\r
182 if(i.ge.nmax)goto 1010
\r
188 9000 format(A17,A3,A2,i4,A4,3F8.3)
\r
192 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\r
194 ccccccccc read data from the second CA file:
\r
195 open(unit=10,file=pdb(2),status='old')
\r
198 read(10,9001,end=1011) s
\r
199 if(i.gt.0.and.s(1:3).eq.'TER')goto 1011
\r
200 if(s(1:3).eq.'ATO')then
\r
201 if(s(13:16).eq.'CA '.or.s(13:16).eq.' CA '.or.
\r
202 & s(13:16).eq.' CA')then
\r
203 if(s(17:17).eq.' '.or.s(17:17).eq.'A')then
\r
205 read(s,9000)du,aanam,du,mm2(i),du,
\r
206 $ xa(1,i,1),xa(2,i,1),xa(3,i,1)
\r
208 if(aanam.eq.aa(j))seq2(i)=slc(j)
\r
211 if(i.ge.nmax)goto 1011
\r
219 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\r
221 *!!! Scale of TM-score in search is based on the smaller protein --------->
\r
223 if(m_d0_min.eq.1)then
\r
224 d0_min=d0_min_input !for search
\r
226 anseq_min=min(nseq1,nseq2)
\r
227 anseq=anseq_min !length for defining TMscore in search
\r
228 d8=1.5*anseq_min**0.3+3.5 !remove pairs with dis>d8 during search & final
\r
229 if(anseq.gt.15)then
\r
230 d0=1.24*(anseq-15)**(1.0/3.0)-1.8 !scale for defining TM-score
\r
234 if(d0.lt.d0_min)d0=d0_min
\r
235 if(m_d0.eq.1)d0=d0_fix
\r
236 d00=d0 !for quickly calculate TM-score in searching
\r
238 if(d00.lt.4.5)d00=4.5
\r
240 nseq=max(nseq1,nseq2)
\r
246 ***** do alignment **************************
\r
247 CALL super_align !to find invmap(j)
\r
249 ************************************************************
\r
250 *** resuperpose to find residues of dis<d8 ------------------------>
\r
253 if(invmap0(j).gt.0)then
\r
256 xtm1(n_al)=xa(1,i,0)
\r
257 ytm1(n_al)=xa(2,i,0)
\r
258 ztm1(n_al)=xa(3,i,0)
\r
259 xtm2(n_al)=xa(1,j,1)
\r
260 ytm2(n_al)=xa(2,j,1)
\r
261 ztm2(n_al)=xa(3,j,1)
\r
262 m1(n_al)=i !for recording residue order
\r
267 call TMscore8(d0_input,n_al,xtm1,ytm1,ztm1,n1,n_al,
\r
268 & xtm2,ytm2,ztm2,n2,TM,Rcomm,Lcomm) !TM-score with dis<d8 only
\r
270 *!!! Output TM-score is based on the second protein------------------>
\r
271 d0_min=0.5 !for output
\r
272 anseq=nseq2 !length for defining final TMscore
\r
273 if(m_ave.eq.1)anseq=(nseq1+nseq2)/2.0 !<L>
\r
274 if(m_ave.eq.2)anseq=min(nseq1,nseq2)
\r
275 if(m_ave.eq.3)anseq=max(nseq1,nseq2)
\r
276 if(anseq.lt.anseq_min)anseq=anseq_min
\r
277 if(m_fix.eq.1)anseq=L_fix !input length
\r
278 if(anseq.gt.15)then
\r
279 d0=1.24*(anseq-15)**(1.0/3.0)-1.8 !scale for defining TM-score
\r
283 if(d0.lt.d0_min)d0=d0_min
\r
284 if(m_d0.eq.1)d0=d0_fix
\r
286 *** remove dis>d8 in normal TM-score calculation for final report----->
\r
290 dis2=sqrt((xtm1(i)-xtm2(i))**2+(ytm1(i)-ytm2(i))**2+
\r
291 & (ztm1(i)-ztm2(i))**2)
\r
302 if(ss1(m1(i)).eq.ss2(m2(i)))then
\r
307 seq_id=float(n_eq)/(n_al+0.00000001)
\r
310 call TMscore(d0_input,n8_al,xtm1,ytm1,ztm1,n1,n8_al,
\r
311 & xtm2,ytm2,ztm2,n2,TM8,Rcomm,Lcomm) !normal TMscore
\r
313 TM8=TM8*n8_al/anseq !TM-score after cutoff
\r
315 ********* for output summary ******************************
\r
317 write(*,*)'*****************************************************',
\r
318 & '*********************'
\r
319 write(*,*)'* TM-align ',
\r
321 write(*,*)'* A protein structural alignment algorithm based on T',
\r
323 write(*,*)'* Reference: Y. Zhang and J. Skolnick, Nucl. Acids Re',
\r
324 & 's. 2005 33, 2302-9 *'
\r
325 write(*,*)'* Comments on the program, please email to: yzhang@ku',
\r
327 write(*,*)'*****************************************************',
\r
328 & '*********************'
\r
330 write(*,101)pdb(1),nseq1
\r
331 101 format('Chain 1:',A10,' Size=',I4)
\r
332 write(*,102)pdb(2),nseq2,int(anseq)
\r
333 102 format('Chain 2:',A10,' Size=',I4,
\r
334 & ' (TM-score is normalized by ',I4,')')
\r
336 write(*,103)n8_al,rmsd8_al,TM8,seq_id
\r
337 103 format('Aligned length=',I4,', RMSD=',f6.2,
\r
338 & ', TM-score=',f7.5,', ID=',f5.3)
\r
341 ********* extract rotation matrix ------------>
\r
354 call u3b(w,r_1,r_2,L,1,rms,u,t,ier) !u rotate r_1 to r_2
\r
356 write(*,*)'-------- rotation matrix to rotate Chain-1 to ',
\r
358 write(*,*)'i t(i) u(i,1) u(i,2) ',
\r
361 write(*,204)i,t(i),u(i,1),u(i,2),u(i,3)
\r
364 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)
\r
365 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)
\r
366 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)
\r
370 204 format(I2,f18.10,f15.10,f15.10,f15.10)
\r
372 ********* for output superposition ******************************
\r
374 1237 format('ATOM ',i5,' CA ',A3,I6,4X,3F8.3)
\r
376 1239 format('CONECT',I5,I5)
\r
378 901 format('select atomno=',I4)
\r
379 104 format('REMARK Chain 1:',A10,' Size=',I4)
\r
380 105 format('REMARK Chain 2:',A10,' Size=',I4,
\r
381 & ' (TM-score is normalized by ',I4,')')
\r
382 106 format('REMARK Aligned length=',I4,', RMSD=',f6.2,
\r
383 & ', TM-score=',f7.5,', ID=',f5.3)
\r
384 OPEN(unit=7,file=outname,status='unknown') !pdb1.aln + pdb2.aln
\r
386 write(7,900)'load inline'
\r
387 write(7,900)'select atomno<2000'
\r
388 write(7,900)'wireframe .45'
\r
389 write(7,900)'select none'
\r
390 write(7,900)'select atomno>2000'
\r
391 write(7,900)'wireframe .20'
\r
392 write(7,900)'color white'
\r
394 dis2=sqrt((xtm1(i)-xtm2(i))**2+
\r
395 & (ytm1(i)-ytm2(i))**2+(ztm1(i)-ztm2(i))**2)
\r
398 write(7,900)'color red'
\r
399 write(7,901)2000+m2(i)
\r
400 write(7,900)'color red'
\r
403 write(7,900)'select all'
\r
405 write(7,104)pdb(1),nseq1
\r
406 write(7,105)pdb(2),nseq2,int(anseq)
\r
407 write(7,106)n8_al,rmsd8_al,TM8,seq_id
\r
410 write(7,1237)m1(i),ss1(m1(i)),mm1(m1(i)),
\r
411 & xtm1(i),ytm1(i),ztm1(i)
\r
415 write(7,1239)m1(i-1),m1(i) !connect atoms
\r
419 write(7,1237)2000+m2(i),ss2(m2(i)),mm2(m2(i)),
\r
420 $ xtm2(i),ytm2(i),ztm2(i)
\r
424 write(7,1239)2000+m2(i-1),2000+m2(i)
\r
429 outnameall_tmp=outname//'_all'
\r
432 if(outnameall_tmp(i:i).ne.' ')then
\r
434 outnameall(k:k)=outnameall_tmp(i:i)
\r
437 OPEN(unit=8,file=outnameall,status='unknown') !pdb1.aln + pdb2.aln
\r
439 write(8,900)'load inline'
\r
440 write(8,900)'select atomno<2000'
\r
441 write(8,900)'wireframe .45'
\r
442 write(8,900)'select none'
\r
443 write(8,900)'select atomno>2000'
\r
444 write(8,900)'wireframe .20'
\r
445 write(8,900)'color white'
\r
447 dis2=sqrt((xtm1(i)-xtm2(i))**2+
\r
448 & (ytm1(i)-ytm2(i))**2+(ztm1(i)-ztm2(i))**2)
\r
451 write(8,900)'color red'
\r
452 write(8,901)2000+m2(i)
\r
453 write(8,900)'color red'
\r
456 write(8,900)'select all'
\r
458 write(8,104)pdb(1),nseq1
\r
459 write(8,105)pdb(2),nseq2,int(anseq)
\r
460 write(8,106)n8_al,rmsd8_al,TM8,seq_id
\r
463 ax=t(1)+u(1,1)*xa(1,i,0)+u(1,2)*xa(2,i,0)+u(1,3)*xa(3,i,0)
\r
464 ay=t(2)+u(2,1)*xa(1,i,0)+u(2,2)*xa(2,i,0)+u(2,3)*xa(3,i,0)
\r
465 az=t(3)+u(3,1)*xa(1,i,0)+u(3,2)*xa(2,i,0)+u(3,3)*xa(3,i,0)
\r
466 write(8,1237)i,ss1(i),mm1(i),ax,ay,az
\r
474 write(8,1237)2000+i,ss2(i),mm2(i),
\r
475 $ xa(1,i,1),xa(2,i,1),xa(3,i,1)
\r
479 write(8,1239)2000+i-1,2000+i
\r
483 *^^^^^^^^^^^^^^^^^^ output finished ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\r
485 ************ output aligned sequences **************************
\r
490 do j=i1_old,m1(i)-1
\r
496 do j=i2_old,m2(i)-1
\r
503 aseq1(ii)=seq1(m1(i))
\r
504 aseq2(ii)=seq2(m2(i))
\r
505 dis2=sqrt((xtm1(i)-xtm2(i))**2+
\r
506 & (ytm1(i)-ytm2(i))**2+(ztm1(i)-ztm2(i))**2)
\r
528 50 format('(":" denotes the residue pairs of distance < 5.0 ',
\r
530 write(*,10)(aseq1(i),i=1,ii)
\r
531 write(*,10)(aseq3(i),i=1,ii)
\r
532 write(*,10)(aseq2(i),i=1,ii)
\r
536 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\r
539 ***********************************************************************
\r
540 ***********************************************************************
\r
541 * Structure superposition
\r
542 ***********************************************************************
\r
543 ***********************************************************************
\r
544 ***********************************************************************
\r
545 SUBROUTINE super_align
\r
546 PARAMETER(nmax=5000)
\r
547 COMMON/BACKBONE/XA(3,nmax,0:1)
\r
548 common/length/nseq1,nseq2
\r
549 common/dpc/score(nmax,nmax),gap_open,invmap(nmax)
\r
550 common/alignrst/invmap0(nmax)
\r
551 common/zscore/zrms,n_al,rmsd_al
\r
553 common/init/invmap_i(nmax)
\r
554 dimension gapp(100)
\r
563 c gapp(i)=-(n_gapp-i)
\r
566 *11111111111111111111111111111111111111111111111111111111
\r
567 * get initial alignment from gapless threading
\r
568 **********************************************************
\r
569 call get_initial !gapless threading
\r
571 invmap(i)=invmap_i(i) !with highest zcore
\r
573 call get_score !TM, matrix score(i,j)
\r
574 if(TM.gt.TMmax)then
\r
577 invmap0(j)=invmap(j)
\r
581 *****************************************************************
\r
582 * initerative alignment, for different gap_open:
\r
583 *****************************************************************
\r
584 DO 111 i_gapp=1,n_gapp !different gap panalties
\r
585 GAP_OPEN=gapp(i_gapp) !gap panalty
\r
586 do 222 id=1,30 !maximum interation is 200
\r
587 call DP(NSEQ1,NSEQ2) !produce alignment invmap(j)
\r
588 * Input: score(i,j), and gap_open
\r
589 * Output: invmap(j)
\r
591 call get_score !calculate TM-score, score(i,j)
\r
592 c record the best alignment in whole search ---------->
\r
593 if(TM.gt.TMmax)then
\r
596 invmap0(j)=invmap(j)
\r
600 diff=abs(TM-TM_old)
\r
601 if(diff.lt.0.000001)goto 33
\r
608 *222222222222222222222222222222222222222222222222222222222
\r
609 * get initial alignment from secondary structure alignment
\r
610 **********************************************************
\r
611 call get_initial2 !DP for secondary structure
\r
613 invmap(i)=invmap_i(i) !with highest zcore
\r
615 call get_score !TM, score(i,j)
\r
616 if(TM.gt.TMmax)then
\r
619 invmap0(j)=invmap(j)
\r
623 *****************************************************************
\r
624 * initerative alignment, for different gap_open:
\r
625 *****************************************************************
\r
626 DO 1111 i_gapp=1,n_gapp !different gap panalties
\r
627 GAP_OPEN=gapp(i_gapp) !gap panalty
\r
628 do 2222 id=1,30 !maximum interation is 200
\r
629 call DP(NSEQ1,NSEQ2) !produce alignment invmap(j)
\r
630 * Input: score(i,j), and gap_open
\r
631 * Output: invmap(j)
\r
633 call get_score !calculate TM-score, score(i,j)
\r
634 c write(*,21)gap_open,rmsd_al,n_al,TM
\r
635 c record the best alignment in whole search ---------->
\r
636 if(TM.gt.TMmax)then
\r
639 invmap0(j)=invmap(j)
\r
643 diff=abs(TM-TM_old)
\r
644 if(diff.lt.0.000001)goto 333
\r
651 *333333333333333333333333333333333333333333333333333333333333
\r
652 * get initial alignment from invmap0+SS
\r
653 *************************************************************
\r
654 call get_initial3 !invmap0+SS
\r
656 invmap(i)=invmap_i(i) !with highest zcore
\r
658 call get_score !TM, score(i,j)
\r
659 if(TM.gt.TMmax)then
\r
662 invmap0(j)=invmap(j)
\r
666 *****************************************************************
\r
667 * initerative alignment, for different gap_open:
\r
668 *****************************************************************
\r
669 DO 1110 i_gapp=1,n_gapp !different gap panalties
\r
670 GAP_OPEN=gapp(i_gapp) !gap panalty
\r
671 do 2220 id=1,30 !maximum interation is 200
\r
672 call DP(NSEQ1,NSEQ2) !produce alignment invmap(j)
\r
673 * Input: score(i,j), and gap_open
\r
674 * Output: invmap(j)
\r
676 call get_score !calculate TM-score, score(i,j)
\r
677 c write(*,21)gap_open,rmsd_al,n_al,TM
\r
678 c record the best alignment in whole search ---------->
\r
679 if(TM.gt.TMmax)then
\r
682 invmap0(j)=invmap(j)
\r
686 diff=abs(TM-TM_old)
\r
687 if(diff.lt.0.000001)goto 330
\r
694 *444444444444444444444444444444444444444444444444444444444
\r
695 * get initial alignment of pieces from gapless threading
\r
696 **********************************************************
\r
697 call get_initial4 !gapless threading
\r
699 invmap(i)=invmap_i(i) !with highest zcore
\r
701 call get_score !TM, matrix score(i,j)
\r
702 if(TM.gt.TMmax)then
\r
705 invmap0(j)=invmap(j)
\r
709 *****************************************************************
\r
710 * initerative alignment, for different gap_open:
\r
711 *****************************************************************
\r
712 DO 44 i_gapp=2,n_gapp !different gap panalties
\r
713 GAP_OPEN=gapp(i_gapp) !gap panalty
\r
714 do 444 id=1,2 !maximum interation is 200
\r
715 call DP(NSEQ1,NSEQ2) !produce alignment invmap(j)
\r
716 * Input: score(i,j), and gap_open
\r
717 * Output: invmap(j)
\r
719 call get_score !calculate TM-score, score(i,j)
\r
720 c record the best alignment in whole search ---------->
\r
721 if(TM.gt.TMmax)then
\r
724 invmap0(j)=invmap(j)
\r
730 c^^^^^^^^^^^^^^^ best alignment invmap0(j) found ^^^^^^^^^^^^^^^^^^
\r
734 **************************************************************
\r
735 * get initial alignment invmap0(i) from gapless threading
\r
736 **************************************************************
\r
737 subroutine get_initial
\r
738 PARAMETER(nmax=5000)
\r
739 COMMON/BACKBONE/XA(3,nmax,0:1)
\r
740 common/length/nseq1,nseq2
\r
741 common/dpc/score(nmax,nmax),gap_open,invmap(nmax)
\r
742 common/alignrst/invmap0(nmax)
\r
743 common/zscore/zrms,n_al,rmsd_al
\r
745 common/init/invmap_i(nmax)
\r
747 aL=min(nseq1,nseq2)
\r
748 idel=aL/2.5 !minimum size of considered fragment
\r
749 if(idel.le.5)idel=5
\r
757 if(i.ge.1.and.i.le.nseq1)then
\r
766 if(GL.gt.GL_max)then
\r
769 invmap_i(i)=invmap(i)
\r
778 **************************************************************
\r
779 * get initial alignment invmap0(i) from secondary structure
\r
780 **************************************************************
\r
781 subroutine get_initial2
\r
782 PARAMETER(nmax=5000)
\r
783 COMMON/BACKBONE/XA(3,nmax,0:1)
\r
784 common/length/nseq1,nseq2
\r
785 common/dpc/score(nmax,nmax),gap_open,invmap(nmax)
\r
786 common/alignrst/invmap0(nmax)
\r
787 common/zscore/zrms,n_al,rmsd_al
\r
789 common/sec/isec(nmax),jsec(nmax)
\r
790 common/init/invmap_i(nmax)
\r
792 ********** assign secondary structures ***************
\r
793 c 1->coil, 2->helix, 3->turn, 4->strand
\r
801 if(j1.ge.1.and.j5.le.nseq1)then
\r
802 dis13=diszy(0,j1,j3)
\r
803 dis14=diszy(0,j1,j4)
\r
804 dis15=diszy(0,j1,j5)
\r
805 dis24=diszy(0,j2,j4)
\r
806 dis25=diszy(0,j2,j5)
\r
807 dis35=diszy(0,j3,j5)
\r
808 isec(i)=make_sec(dis13,dis14,dis15,dis24,dis25,dis35)
\r
818 if(j1.ge.1.and.j5.le.nseq2)then
\r
819 dis13=diszy(1,j1,j3)
\r
820 dis14=diszy(1,j1,j4)
\r
821 dis15=diszy(1,j1,j5)
\r
822 dis24=diszy(1,j2,j4)
\r
823 dis25=diszy(1,j2,j5)
\r
824 dis35=diszy(1,j3,j5)
\r
825 jsec(i)=make_sec(dis13,dis14,dis15,dis24,dis25,dis35)
\r
828 call smooth !smooth the assignment
\r
830 ********** score matrix **************************
\r
833 if(isec(i).eq.jsec(j))then
\r
841 ********** find initial alignment: invmap(j) ************
\r
842 gap_open=-1.0 !should be -1
\r
843 call DP(NSEQ1,NSEQ2) !produce alignment invmap(j)
\r
845 invmap_i(i)=invmap(i)
\r
848 *^^^^^^^^^^^^ initial alignment done ^^^^^^^^^^^^^^^^^^^^^^
\r
852 **************************************************************
\r
853 * get initial alignment invmap0(i) from secondary structure
\r
854 * and previous alignments
\r
855 **************************************************************
\r
856 subroutine get_initial3
\r
857 PARAMETER(nmax=5000)
\r
858 COMMON/BACKBONE/XA(3,nmax,0:1)
\r
859 common/length/nseq1,nseq2
\r
860 common/dpc/score(nmax,nmax),gap_open,invmap(nmax)
\r
861 common/alignrst/invmap0(nmax)
\r
862 common/zscore/zrms,n_al,rmsd_al
\r
864 common/sec/isec(nmax),jsec(nmax)
\r
865 common/init/invmap_i(nmax)
\r
867 ********** score matrix **************************
\r
869 invmap(i)=invmap0(i)
\r
871 call get_score1 !get score(i,j) using RMSD martix
\r
874 if(isec(i).eq.jsec(j))then
\r
875 score(i,j)=0.5+score(i,j)
\r
877 score(i,j)=score(i,j)
\r
882 ********** find initial alignment: invmap(j) ************
\r
883 gap_open=-1.0 !should be -1
\r
884 call DP(NSEQ1,NSEQ2) !produce alignment invmap(j)
\r
886 invmap_i(i)=invmap(i)
\r
889 *^^^^^^^^^^^^ initial alignment done ^^^^^^^^^^^^^^^^^^^^^^
\r
893 **************************************************************
\r
894 * get initial alignment invmap0(i) from fragment gapless threading
\r
895 **************************************************************
\r
896 subroutine get_initial4
\r
897 PARAMETER(nmax=5000)
\r
898 COMMON/BACKBONE/XA(3,nmax,0:1)
\r
899 common/length/nseq1,nseq2
\r
900 common/dpc/score(nmax,nmax),gap_open,invmap(nmax)
\r
901 common/alignrst/invmap0(nmax)
\r
902 common/zscore/zrms,n_al,rmsd_al
\r
904 common/init/invmap_i(nmax)
\r
905 common/initial4/mm1(nmax),mm2(nmax)
\r
908 dimension ifr2(2,nmax,nmax),Lfr2(2,nmax),Lfr_max2(2),i_fr2(2)
\r
909 dimension ifr(nmax)
\r
910 dimension mm(2,nmax)
\r
912 fra_min=4 !>=4,minimum fragment for search
\r
913 fra_min1=fra_min-1 !cutoff for shift, save time
\r
917 ccc Find the smallest continuous fragments -------->
\r
928 r_min=nseq1/3.0 !minimum fragment, in case too small protein
\r
931 r_min=nseq2/3.0 !minimum fragment, in case too small protein
\r
933 if(r_min.gt.fra_min)r_min=fra_min
\r
934 20 nfr=1 !number of fragments
\r
935 j=1 !number of residue at nf-fragment
\r
936 ifr2(k,nfr,j)=1 !what residue
\r
937 Lfr2(k,nfr)=j !length of the fragment
\r
939 dis=diszy(k-1,i-1,i)
\r
941 if(dcu.gt.dcu0)then
\r
943 if(dis.gt.dcu_min)then
\r
947 elseif(mm(k,i).eq.(mm(k,i-1)+1))then
\r
949 if(dis.gt.dcu_min)then
\r
966 i_fr2(k)=1 !ID of the maximum piece
\r
968 if(Lfr_max.lt.Lfr2(k,i))then
\r
973 if(Lfr_max.lt.r_min)then
\r
978 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\r
980 ccc select what piece will be used (this may araise ansysmetry, but
\r
981 ccc only when L1=L2 and Lfr1=Lfr2 and L1 ne Lfr1
\r
982 ccc if L1=Lfr1 and L2=Lfr2 (normal proteins), it will be the same as initial1
\r
984 if(Lfr2(1,i_fr2(1)).lt.Lfr2(2,i_fr2(2)))then
\r
986 elseif(Lfr2(1,i_fr2(1)).gt.Lfr2(2,i_fr2(2)))then
\r
989 if(nseq1.le.nseq2)then
\r
996 L_fr=Lfr2(mark,i_fr2(mark))
\r
998 ifr(i)=ifr2(mark,i_fr2(mark),i)
\r
1001 if(mark.eq.1)then !non-redundant to get_initial1
\r
1006 if(L_fr.eq.nseq0)then
\r
1007 n1=int(nseq0*0.1) !0
\r
1008 n2=int(nseq0*0.89) !2
\r
1017 ccc get initial ------------->
\r
1018 if(mark.eq.1)then !nseq1 as the smallest one
\r
1020 aL=min(nseq1_,nseq2)
\r
1021 idel=aL/2.5 !minimum size of considered fragment
\r
1022 if(idel.le.fra_min1)idel=fra_min1
\r
1023 n1=-nseq2+idel !shift1
\r
1024 n2=nseq1_-idel !shift2
\r
1030 if(i.ge.1.and.i.le.nseq1_)then
\r
1039 if(GL.gt.GL_max)then
\r
1042 invmap_i(i)=invmap(i)
\r
1047 else !@@@@@@@@@@@@@@@@@@@@
\r
1049 aL=min(nseq1,nseq2_)
\r
1050 idel=aL/2.5 !minimum size of considered fragment
\r
1051 if(idel.le.fra_min1)idel=fra_min1
\r
1062 if(i.ge.1.and.i.le.nseq1)then
\r
1069 if(GL.gt.GL_max)then
\r
1072 invmap_i(i)=invmap(i)
\r
1082 **************************************************************
\r
1083 * smooth the secondary structure assignment
\r
1084 **************************************************************
\r
1086 PARAMETER(nmax=5000)
\r
1087 common/sec/isec(nmax),jsec(nmax)
\r
1088 common/length/nseq1,nseq2
\r
1090 *** smooth single -------------->
\r
1091 *** --x-- => -----
\r
1093 if(isec(i).eq.2.or.isec(i).eq.4)then
\r
1095 if(isec(i-2).ne.j)then
\r
1096 if(isec(i-1).ne.j)then
\r
1097 if(isec(i+1).ne.j)then
\r
1098 if(isec(i+1).ne.j)then
\r
1107 if(jsec(i).eq.2.or.jsec(i).eq.4)then
\r
1109 if(jsec(i-2).ne.j)then
\r
1110 if(jsec(i-1).ne.j)then
\r
1111 if(jsec(i+1).ne.j)then
\r
1112 if(jsec(i+1).ne.j)then
\r
1121 *** smooth double -------------->
\r
1122 *** --xx-- => ------
\r
1124 if(isec(i).ne.2)then
\r
1125 if(isec(i+1).ne.2)then
\r
1126 if(isec(i+2).eq.2)then
\r
1127 if(isec(i+3).eq.2)then
\r
1128 if(isec(i+4).ne.2)then
\r
1129 if(isec(i+5).ne.2)then
\r
1139 if(isec(i).ne.4)then
\r
1140 if(isec(i+1).ne.4)then
\r
1141 if(isec(i+2).eq.4)then
\r
1142 if(isec(i+3).eq.4)then
\r
1143 if(isec(i+4).ne.4)then
\r
1144 if(isec(i+5).ne.4)then
\r
1155 if(jsec(i).ne.2)then
\r
1156 if(jsec(i+1).ne.2)then
\r
1157 if(jsec(i+2).eq.2)then
\r
1158 if(jsec(i+3).eq.2)then
\r
1159 if(jsec(i+4).ne.2)then
\r
1160 if(jsec(i+5).ne.2)then
\r
1170 if(jsec(i).ne.4)then
\r
1171 if(jsec(i+1).ne.4)then
\r
1172 if(jsec(i+2).eq.4)then
\r
1173 if(jsec(i+3).eq.4)then
\r
1174 if(jsec(i+4).ne.4)then
\r
1175 if(jsec(i+5).ne.4)then
\r
1186 *** connect -------------->
\r
1189 if(isec(i).eq.2)then
\r
1190 if(isec(i+1).ne.2)then
\r
1191 if(isec(i+2).eq.2)then
\r
1197 if(isec(i).eq.4)then
\r
1198 if(isec(i+1).ne.4)then
\r
1199 if(isec(i+2).eq.4)then
\r
1206 if(jsec(i).eq.2)then
\r
1207 if(jsec(i+1).ne.2)then
\r
1208 if(jsec(i+2).eq.2)then
\r
1214 if(jsec(i).eq.4)then
\r
1215 if(jsec(i+1).ne.4)then
\r
1216 if(jsec(i+2).eq.4)then
\r
1226 *************************************************************
\r
1227 * assign secondary structure:
\r
1228 *************************************************************
\r
1229 function diszy(i,i1,i2)
\r
1230 PARAMETER(nmax=5000)
\r
1231 COMMON/BACKBONE/XA(3,nmax,0:1)
\r
1232 diszy=sqrt((xa(1,i1,i)-xa(1,i2,i))**2
\r
1233 & +(xa(2,i1,i)-xa(2,i2,i))**2
\r
1234 & +(xa(3,i1,i)-xa(3,i2,i))**2)
\r
1238 *************************************************************
\r
1239 * assign secondary structure:
\r
1240 *************************************************************
\r
1241 function make_sec(dis13,dis14,dis15,dis24,dis25,dis35)
\r
1244 if(abs(dis15-6.37).lt.delta)then
\r
1245 if(abs(dis14-5.18).lt.delta)then
\r
1246 if(abs(dis25-5.18).lt.delta)then
\r
1247 if(abs(dis13-5.45).lt.delta)then
\r
1248 if(abs(dis24-5.45).lt.delta)then
\r
1249 if(abs(dis35-5.45).lt.delta)then
\r
1259 if(abs(dis15-13).lt.delta)then
\r
1260 if(abs(dis14-10.4).lt.delta)then
\r
1261 if(abs(dis25-10.4).lt.delta)then
\r
1262 if(abs(dis13-6.1).lt.delta)then
\r
1263 if(abs(dis24-6.1).lt.delta)then
\r
1264 if(abs(dis35-6.1).lt.delta)then
\r
1265 make_sec=4 !strand
\r
1273 if(dis15.lt.8)then
\r
1280 ****************************************************************
\r
1281 * quickly calculate TM-score with given invmap(i) in 3 iterations
\r
1282 ****************************************************************
\r
1283 subroutine get_GL(GL)
\r
1284 PARAMETER(nmax=5000)
\r
1285 common/length/nseq1,nseq2
\r
1286 COMMON/BACKBONE/XA(3,nmax,0:1)
\r
1287 common/dpc/score(nmax,nmax),gap_open,invmap(nmax)
\r
1288 common/zscore/zrms,n_al,rmsd_al
\r
1289 common/d0/d0,anseq
\r
1290 dimension xtm1(nmax),ytm1(nmax),ztm1(nmax)
\r
1291 dimension xtm2(nmax),ytm2(nmax),ztm2(nmax)
\r
1292 common/TM/TM,TMmax
\r
1293 common/n1n2/n1(nmax),n2(nmax)
\r
1294 common/d00/d00,d002
\r
1296 dimension xo1(nmax),yo1(nmax),zo1(nmax)
\r
1297 dimension xo2(nmax),yo2(nmax),zo2(nmax)
\r
1298 dimension dis2(nmax)
\r
1301 double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)
\r
1302 double precision u(3,3),t(3),rms,drms !armsd is real
\r
1306 c calculate RMSD between aligned structures and rotate the structures -->
\r
1309 i=invmap(j) !j aligned to i
\r
1312 r_1(1,n_al)=xa(1,i,0)
\r
1313 r_1(2,n_al)=xa(2,i,0)
\r
1314 r_1(3,n_al)=xa(3,i,0)
\r
1315 r_2(1,n_al)=xa(1,j,1)
\r
1316 r_2(2,n_al)=xa(2,j,1)
\r
1317 r_2(3,n_al)=xa(3,j,1)
\r
1318 xo1(n_al)=xa(1,i,0)
\r
1319 yo1(n_al)=xa(2,i,0)
\r
1320 zo1(n_al)=xa(3,i,0)
\r
1321 xo2(n_al)=xa(1,j,1)
\r
1322 yo2(n_al)=xa(2,j,1)
\r
1323 zo2(n_al)=xa(3,j,1)
\r
1326 call u3b(w,r_1,r_2,n_al,1,rms,u,t,ier) !u rotate r_1 to r_2
\r
1329 xx=t(1)+u(1,1)*xo1(i)+u(1,2)*yo1(i)+u(1,3)*zo1(i)
\r
1330 yy=t(2)+u(2,1)*xo1(i)+u(2,2)*yo1(i)+u(2,3)*zo1(i)
\r
1331 zz=t(3)+u(3,1)*xo1(i)+u(3,2)*yo1(i)+u(3,3)*zo1(i)
\r
1332 dis2(i)=(xx-xo2(i))**2+(yy-yo2(i))**2+(zz-zo2(i))**2
\r
1333 GL=GL+1/(1+dis2(i)/(d0**2))
\r
1335 ccc for next iteration------------->
\r
1339 if(dis2(i).le.d002t)then
\r
1349 if(j.lt.3.and.n_al.gt.3)then
\r
1354 call u3b(w,r_1,r_2,L,1,rms,u,t,ier) !u rotate r_1 to r_2
\r
1357 xx=t(1)+u(1,1)*xo1(i)+u(1,2)*yo1(i)+u(1,3)*zo1(i)
\r
1358 yy=t(2)+u(2,1)*xo1(i)+u(2,2)*yo1(i)+u(2,3)*zo1(i)
\r
1359 zz=t(3)+u(3,1)*xo1(i)+u(3,2)*yo1(i)+u(3,3)*zo1(i)
\r
1360 dis2(i)=(xx-xo2(i))**2+(yy-yo2(i))**2+(zz-zo2(i))**2
\r
1361 G2=G2+1/(1+dis2(i)/(d0**2))
\r
1363 ccc for next iteration------------->
\r
1367 if(dis2(i).le.d002t)then
\r
1377 if(j.lt.3.and.n_al.gt.3)then
\r
1382 call u3b(w,r_1,r_2,L,1,rms,u,t,ier) !u rotate r_1 to r_2
\r
1385 xx=t(1)+u(1,1)*xo1(i)+u(1,2)*yo1(i)+u(1,3)*zo1(i)
\r
1386 yy=t(2)+u(2,1)*xo1(i)+u(2,2)*yo1(i)+u(2,3)*zo1(i)
\r
1387 zz=t(3)+u(3,1)*xo1(i)+u(3,2)*yo1(i)+u(3,3)*zo1(i)
\r
1388 dis2(i)=(xx-xo2(i))**2+(yy-yo2(i))**2+(zz-zo2(i))**2
\r
1389 G3=G3+1/(1+dis2(i)/(d0**2))
\r
1394 c^^^^^^^^^^^^^^^^ GL done ^^^^^^^^^^^^^^^^^^^^^^^^^^^
\r
1398 ****************************************************************
\r
1399 * with invmap(i) calculate TM-score and martix score(i,j) for rotation
\r
1400 ****************************************************************
\r
1401 subroutine get_score
\r
1402 PARAMETER(nmax=5000)
\r
1403 common/length/nseq1,nseq2
\r
1404 COMMON/BACKBONE/XA(3,nmax,0:1)
\r
1405 common/dpc/score(nmax,nmax),gap_open,invmap(nmax)
\r
1406 common/zscore/zrms,n_al,rmsd_al
\r
1407 common/d0/d0,anseq
\r
1408 dimension xtm1(nmax),ytm1(nmax),ztm1(nmax)
\r
1409 dimension xtm2(nmax),ytm2(nmax),ztm2(nmax)
\r
1410 common/TM/TM,TMmax
\r
1411 common/n1n2/n1(nmax),n2(nmax)
\r
1414 double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)
\r
1415 double precision u(3,3),t(3),rms,drms !armsd is real
\r
1419 c calculate RMSD between aligned structures and rotate the structures -->
\r
1422 i=invmap(j) !j aligned to i
\r
1426 xtm1(n_al)=xa(1,i,0) !for TM-score
\r
1427 ytm1(n_al)=xa(2,i,0)
\r
1428 ztm1(n_al)=xa(3,i,0)
\r
1429 xtm2(n_al)=xa(1,j,1)
\r
1430 ytm2(n_al)=xa(2,j,1)
\r
1431 ztm2(n_al)=xa(3,j,1)
\r
1432 ccc for rotation matrix:
\r
1433 r_1(1,n_al)=xa(1,i,0)
\r
1434 r_1(2,n_al)=xa(2,i,0)
\r
1435 r_1(3,n_al)=xa(3,i,0)
\r
1438 *** calculate TM-score for the given alignment----------->
\r
1440 call TMscore8_search(d0_input,n_al,xtm1,ytm1,ztm1,n1,
\r
1441 & n_al,xtm2,ytm2,ztm2,n2,TM,Rcomm,Lcomm) !simplified search engine
\r
1442 TM=TM*n_al/anseq !TM-score
\r
1443 *** calculate score matrix score(i,j)------------------>
\r
1449 call u3b(w,r_1,r_2,n_al,1,rms,u,t,ier) !u rotate r_1 to r_2
\r
1451 xx=t(1)+u(1,1)*xa(1,i,0)+u(1,2)*xa(2,i,0)+u(1,3)*xa(3,i,0)
\r
1452 yy=t(2)+u(2,1)*xa(1,i,0)+u(2,2)*xa(2,i,0)+u(2,3)*xa(3,i,0)
\r
1453 zz=t(3)+u(3,1)*xa(1,i,0)+u(3,2)*xa(2,i,0)+u(3,3)*xa(3,i,0)
\r
1455 dd=(xx-xa(1,j,1))**2+(yy-xa(2,j,1))**2+(zz-xa(3,j,1))**2
\r
1456 score(i,j)=1/(1+dd/d0**2)
\r
1460 c^^^^^^^^^^^^^^^^ score(i,j) done ^^^^^^^^^^^^^^^^^^^^^^^^^^^
\r
1464 ****************************************************************
\r
1465 * with invmap(i) calculate score(i,j) using RMSD rotation
\r
1466 ****************************************************************
\r
1467 subroutine get_score1
\r
1468 PARAMETER(nmax=5000)
\r
1469 common/length/nseq1,nseq2
\r
1470 COMMON/BACKBONE/XA(3,nmax,0:1)
\r
1471 common/dpc/score(nmax,nmax),gap_open,invmap(nmax)
\r
1472 common/zscore/zrms,n_al,rmsd_al
\r
1473 common/d0/d0,anseq
\r
1474 common/d0min/d0_min
\r
1475 dimension xtm1(nmax),ytm1(nmax),ztm1(nmax)
\r
1476 dimension xtm2(nmax),ytm2(nmax),ztm2(nmax)
\r
1477 common/TM/TM,TMmax
\r
1478 common/n1n2/n1(nmax),n2(nmax)
\r
1481 double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)
\r
1482 double precision u(3,3),t(3),rms,drms !armsd is real
\r
1486 c calculate RMSD between aligned structures and rotate the structures -->
\r
1489 i=invmap(j) !j aligned to i
\r
1492 ccc for rotation matrix:
\r
1493 r_1(1,n_al)=xa(1,i,0)
\r
1494 r_1(2,n_al)=xa(2,i,0)
\r
1495 r_1(3,n_al)=xa(3,i,0)
\r
1496 r_2(1,n_al)=xa(1,j,1)
\r
1497 r_2(2,n_al)=xa(2,j,1)
\r
1498 r_2(3,n_al)=xa(3,j,1)
\r
1501 *** calculate score matrix score(i,j)------------------>
\r
1502 call u3b(w,r_1,r_2,n_al,1,rms,u,t,ier) !u rotate r_1 to r_2
\r
1504 if(d01.lt.d0_min)d01=d0_min
\r
1507 xx=t(1)+u(1,1)*xa(1,i,0)+u(1,2)*xa(2,i,0)+u(1,3)*xa(3,i,0)
\r
1508 yy=t(2)+u(2,1)*xa(1,i,0)+u(2,2)*xa(2,i,0)+u(2,3)*xa(3,i,0)
\r
1509 zz=t(3)+u(3,1)*xa(1,i,0)+u(3,2)*xa(2,i,0)+u(3,3)*xa(3,i,0)
\r
1511 dd=(xx-xa(1,j,1))**2+(yy-xa(2,j,1))**2+(zz-xa(3,j,1))**2
\r
1512 score(i,j)=1/(1+dd/d02)
\r
1516 c^^^^^^^^^^^^^^^^ score(i,j) done ^^^^^^^^^^^^^^^^^^^^^^^^^^^
\r
1521 *************************************************************************
\r
1522 *************************************************************************
\r
1523 * This is a subroutine to compare two structures and find the
\r
1524 * superposition that has the maximum TM-score.
\r
1526 * L1--Length of the first structure
\r
1527 * (x1(i),y1(i),z1(i))--coordinates of i'th residue at the first structure
\r
1528 * n1(i)--Residue sequence number of i'th residue at the first structure
\r
1529 * L2--Length of the second structure
\r
1530 * (x2(i),y2(i),z2(i))--coordinates of i'th residue at the second structure
\r
1531 * n2(i)--Residue sequence number of i'th residue at the second structure
\r
1532 * TM--TM-score of the comparison
\r
1533 * Rcomm--RMSD of two structures in the common aligned residues
\r
1534 * Lcomm--Length of the common aligned regions
\r
1537 * 1, Always put native as the second structure, by which TM-score
\r
1539 * 2, The returned (x1(i),y1(i),z1(i)) are the rotated structure after
\r
1540 * TM-score superposition.
\r
1541 *************************************************************************
\r
1542 *************************************************************************
\r
1543 *** dis<8, simplified search engine
\r
1544 subroutine TMscore8_search(dx,L1,x1,y1,z1,n1,L2,x2,y2,z2,n2,
\r
1546 PARAMETER(nmax=5000)
\r
1547 common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax)
\r
1548 common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
\r
1550 common/d0min/d0_min
\r
1551 common/align/n_ali,iA(nmax),iB(nmax)
\r
1552 common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score
\r
1553 dimension k_ali(nmax),k_ali0(nmax)
\r
1554 dimension L_ini(100),iq(nmax)
\r
1555 common/scores/score
\r
1556 double precision score,score_max
\r
1557 dimension xa(nmax),ya(nmax),za(nmax)
\r
1558 dimension iL0(nmax)
\r
1560 dimension x1(nmax),y1(nmax),z1(nmax),n1(nmax)
\r
1561 dimension x2(nmax),y2(nmax),z2(nmax),n2(nmax)
\r
1564 double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)
\r
1565 double precision u(3,3),t(3),rms,drms !armsd is real
\r
1569 ********* convert input data ***************************
\r
1570 * because L1=L2 in this special case---------->
\r
1585 n_ali=L1 !number of aligned residues
\r
1591 *** d0------------->
\r
1593 if(d0.lt.d0_min)d0=d0_min
\r
1594 *** d0_search ----->
\r
1596 if(d0_search.gt.8)d0_search=8
\r
1597 if(d0_search.lt.4.5)d0_search=4.5
\r
1598 *** iterative parameters ----->
\r
1599 n_it=20 !maximum number of iterations
\r
1600 d_output=5 !for output alignment
\r
1601 n_init_max=6 !maximum number of L_init
\r
1604 if(n_ali.lt.4)L_ini_min=n_ali
\r
1605 do i=1,n_init_max-1
\r
1607 L_ini(n_init)=n_ali/2**(n_init-1)
\r
1608 if(L_ini(n_init).le.L_ini_min)then
\r
1609 L_ini(n_init)=L_ini_min
\r
1614 L_ini(n_init)=L_ini_min
\r
1617 ******************************************************************
\r
1618 * find the maximum score starting from local structures superposition
\r
1619 ******************************************************************
\r
1620 score_max=-1 !TM-score
\r
1621 do 333 i_init=1,n_init
\r
1622 L_init=L_ini(i_init)
\r
1623 iL_max=n_ali-L_init+1
\r
1625 do i=1,iL_max,40 !this is the simplification!
\r
1629 if(iL0(k).lt.iL_max)then
\r
1634 do 300 i_shift=1,n_shift
\r
1639 k=iL+i-1 ![1,n_ali] common aligned
\r
1640 r_1(1,i)=xa(iA(k))
\r
1641 r_1(2,i)=ya(iA(k))
\r
1642 r_1(3,i)=za(iA(k))
\r
1643 r_2(1,i)=xb(iB(k))
\r
1644 r_2(2,i)=yb(iB(k))
\r
1645 r_2(3,i)=zb(iB(k))
\r
1650 call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2
\r
1651 if(i_init.eq.1)then !global superposition
\r
1652 armsd=dsqrt(rms/LL)
\r
1656 xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
\r
1657 yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
\r
1658 zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
\r
1661 call score_fun8 !init, get scores, n_cut+i_ali(i) for iteration
\r
1662 if(score_max.lt.score)then
\r
1666 k_ali0(i)=k_ali(i)
\r
1669 *** iteration for extending ---------------------------------->
\r
1675 m=i_ali(i) ![1,n_ali]
\r
1676 r_1(1,i)=xa(iA(m))
\r
1677 r_1(2,i)=ya(iA(m))
\r
1678 r_1(3,i)=za(iA(m))
\r
1679 r_2(1,i)=xb(iB(m))
\r
1680 r_2(2,i)=yb(iB(m))
\r
1681 r_2(3,i)=zb(iB(m))
\r
1686 call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2
\r
1688 xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
\r
1689 yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
\r
1690 zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
\r
1692 call score_fun8 !get scores, n_cut+i_ali(i) for iteration
\r
1693 if(score_max.lt.score)then
\r
1697 k_ali0(i)=k_ali(i)
\r
1700 if(it.eq.n_it)goto 302
\r
1701 if(n_cut.eq.ka)then
\r
1704 if(i_ali(i).eq.k_ali(i))neq=neq+1
\r
1706 if(n_cut.eq.neq)goto 302
\r
1708 301 continue !for iteration
\r
1710 300 continue !for shift
\r
1711 333 continue !for initial length, L_ali/M
\r
1713 ******** return the final rotation ****************
\r
1716 m=k_ali0(i) !record of the best alignment
\r
1717 r_1(1,i)=xa(iA(m))
\r
1718 r_1(2,i)=ya(iA(m))
\r
1719 r_1(3,i)=za(iA(m))
\r
1720 r_2(1,i)=xb(iB(m))
\r
1721 r_2(2,i)=yb(iB(m))
\r
1722 r_2(3,i)=zb(iB(m))
\r
1725 call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2
\r
1727 x1(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
\r
1728 y1(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
\r
1729 z1(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
\r
1733 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\r
1737 *************************************************************************
\r
1738 *************************************************************************
\r
1739 * This is a subroutine to compare two structures and find the
\r
1740 * superposition that has the maximum TM-score.
\r
1742 * L1--Length of the first structure
\r
1743 * (x1(i),y1(i),z1(i))--coordinates of i'th residue at the first structure
\r
1744 * n1(i)--Residue sequence number of i'th residue at the first structure
\r
1745 * L2--Length of the second structure
\r
1746 * (x2(i),y2(i),z2(i))--coordinates of i'th residue at the second structure
\r
1747 * n2(i)--Residue sequence number of i'th residue at the second structure
\r
1748 * TM--TM-score of the comparison
\r
1749 * Rcomm--RMSD of two structures in the common aligned residues
\r
1750 * Lcomm--Length of the common aligned regions
\r
1753 * 1, Always put native as the second structure, by which TM-score
\r
1755 * 2, The returned (x1(i),y1(i),z1(i)) are the rotated structure after
\r
1756 * TM-score superposition.
\r
1757 *************************************************************************
\r
1758 *************************************************************************
\r
1759 *** dis<8, but same search engine
\r
1760 subroutine TMscore8(dx,L1,x1,y1,z1,n1,L2,x2,y2,z2,n2,
\r
1762 PARAMETER(nmax=5000)
\r
1763 common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax)
\r
1764 common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
\r
1766 common/d0min/d0_min
\r
1767 common/align/n_ali,iA(nmax),iB(nmax)
\r
1768 common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score
\r
1769 dimension k_ali(nmax),k_ali0(nmax)
\r
1770 dimension L_ini(100),iq(nmax)
\r
1771 common/scores/score
\r
1772 double precision score,score_max
\r
1773 dimension xa(nmax),ya(nmax),za(nmax)
\r
1775 dimension x1(nmax),y1(nmax),z1(nmax),n1(nmax)
\r
1776 dimension x2(nmax),y2(nmax),z2(nmax),n2(nmax)
\r
1779 double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)
\r
1780 double precision u(3,3),t(3),rms,drms !armsd is real
\r
1784 ********* convert input data ***************************
\r
1785 * because L1=L2 in this special case---------->
\r
1800 n_ali=L1 !number of aligned residues
\r
1806 *** d0------------->
\r
1808 if(d0.lt.d0_min)d0=d0_min
\r
1809 *** d0_search ----->
\r
1811 if(d0_search.gt.8)d0_search=8
\r
1812 if(d0_search.lt.4.5)d0_search=4.5
\r
1813 *** iterative parameters ----->
\r
1814 n_it=20 !maximum number of iterations
\r
1815 d_output=5 !for output alignment
\r
1816 n_init_max=6 !maximum number of L_init
\r
1819 if(n_ali.lt.4)L_ini_min=n_ali
\r
1820 do i=1,n_init_max-1
\r
1822 L_ini(n_init)=n_ali/2**(n_init-1)
\r
1823 if(L_ini(n_init).le.L_ini_min)then
\r
1824 L_ini(n_init)=L_ini_min
\r
1829 L_ini(n_init)=L_ini_min
\r
1832 ******************************************************************
\r
1833 * find the maximum score starting from local structures superposition
\r
1834 ******************************************************************
\r
1835 score_max=-1 !TM-score
\r
1836 do 333 i_init=1,n_init
\r
1837 L_init=L_ini(i_init)
\r
1838 iL_max=n_ali-L_init+1
\r
1839 do 300 iL=1,iL_max !on aligned residues, [1,nseqA]
\r
1843 k=iL+i-1 ![1,n_ali] common aligned
\r
1844 r_1(1,i)=xa(iA(k))
\r
1845 r_1(2,i)=ya(iA(k))
\r
1846 r_1(3,i)=za(iA(k))
\r
1847 r_2(1,i)=xb(iB(k))
\r
1848 r_2(2,i)=yb(iB(k))
\r
1849 r_2(3,i)=zb(iB(k))
\r
1854 call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2
\r
1855 if(i_init.eq.1)then !global superposition
\r
1856 armsd=dsqrt(rms/LL)
\r
1860 xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
\r
1861 yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
\r
1862 zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
\r
1865 call score_fun8 !init, get scores, n_cut+i_ali(i) for iteration
\r
1866 if(score_max.lt.score)then
\r
1870 k_ali0(i)=k_ali(i)
\r
1873 *** iteration for extending ---------------------------------->
\r
1879 m=i_ali(i) ![1,n_ali]
\r
1880 r_1(1,i)=xa(iA(m))
\r
1881 r_1(2,i)=ya(iA(m))
\r
1882 r_1(3,i)=za(iA(m))
\r
1883 r_2(1,i)=xb(iB(m))
\r
1884 r_2(2,i)=yb(iB(m))
\r
1885 r_2(3,i)=zb(iB(m))
\r
1890 call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2
\r
1892 xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
\r
1893 yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
\r
1894 zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
\r
1896 call score_fun8 !get scores, n_cut+i_ali(i) for iteration
\r
1897 if(score_max.lt.score)then
\r
1901 k_ali0(i)=k_ali(i)
\r
1904 if(it.eq.n_it)goto 302
\r
1905 if(n_cut.eq.ka)then
\r
1908 if(i_ali(i).eq.k_ali(i))neq=neq+1
\r
1910 if(n_cut.eq.neq)goto 302
\r
1912 301 continue !for iteration
\r
1914 300 continue !for shift
\r
1915 333 continue !for initial length, L_ali/M
\r
1917 ******** return the final rotation ****************
\r
1920 m=k_ali0(i) !record of the best alignment
\r
1921 r_1(1,i)=xa(iA(m))
\r
1922 r_1(2,i)=ya(iA(m))
\r
1923 r_1(3,i)=za(iA(m))
\r
1924 r_2(1,i)=xb(iB(m))
\r
1925 r_2(2,i)=yb(iB(m))
\r
1926 r_2(3,i)=zb(iB(m))
\r
1929 call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2
\r
1931 x1(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
\r
1932 y1(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
\r
1933 z1(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
\r
1937 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\r
1941 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
\r
1942 c 1, collect those residues with dis<d;
\r
1943 c 2, calculate score_GDT, score_maxsub, score_TM
\r
1944 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
\r
1945 subroutine score_fun8
\r
1946 PARAMETER(nmax=5000)
\r
1948 common/stru/xa(nmax),ya(nmax),za(nmax),xb(nmax),yb(nmax),zb(nmax)
\r
1949 common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
\r
1951 common/align/n_ali,iA(nmax),iB(nmax)
\r
1952 common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score
\r
1953 common/scores/score
\r
1954 double precision score,score_max
\r
1958 21 n_cut=0 !number of residue-pairs dis<d, for iteration
\r
1959 score_sum=0 !TMscore
\r
1961 i=iA(k) ![1,nseqA] reoder number of structureA
\r
1962 j=iB(k) ![1,nseqB]
\r
1963 dis=sqrt((xa(i)-xb(j))**2+(ya(i)-yb(j))**2+(za(i)-zb(j))**2)
\r
1964 if(dis.lt.d_tmp)then
\r
1966 i_ali(n_cut)=k ![1,n_ali], mark the residue-pairs in dis<d
\r
1969 score_sum=score_sum+1/(1+(dis/d0)**2)
\r
1972 if(n_cut.lt.3.and.n_ali.gt.3)then
\r
1976 score=score_sum/float(nseqB) !TM-score
\r
1981 *************************************************************************
\r
1982 *************************************************************************
\r
1983 * This is a subroutine to compare two structures and find the
\r
1984 * superposition that has the maximum TM-score.
\r
1986 * L1--Length of the first structure
\r
1987 * (x1(i),y1(i),z1(i))--coordinates of i'th residue at the first structure
\r
1988 * n1(i)--Residue sequence number of i'th residue at the first structure
\r
1989 * L2--Length of the second structure
\r
1990 * (x2(i),y2(i),z2(i))--coordinates of i'th residue at the second structure
\r
1991 * n2(i)--Residue sequence number of i'th residue at the second structure
\r
1992 * TM--TM-score of the comparison
\r
1993 * Rcomm--RMSD of two structures in the common aligned residues
\r
1994 * Lcomm--Length of the common aligned regions
\r
1997 * 1, Always put native as the second structure, by which TM-score
\r
1999 * 2, The returned (x1(i),y1(i),z1(i)) are the rotated structure after
\r
2000 * TM-score superposition.
\r
2001 *************************************************************************
\r
2002 *************************************************************************
\r
2003 *** normal TM-score:
\r
2004 subroutine TMscore(dx,L1,x1,y1,z1,n1,L2,x2,y2,z2,n2,
\r
2006 PARAMETER(nmax=5000)
\r
2007 common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax)
\r
2008 common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
\r
2010 common/d0min/d0_min
\r
2011 common/align/n_ali,iA(nmax),iB(nmax)
\r
2012 common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score
\r
2013 dimension k_ali(nmax),k_ali0(nmax)
\r
2014 dimension L_ini(100),iq(nmax)
\r
2015 common/scores/score
\r
2016 double precision score,score_max
\r
2017 dimension xa(nmax),ya(nmax),za(nmax)
\r
2019 dimension x1(nmax),y1(nmax),z1(nmax),n1(nmax)
\r
2020 dimension x2(nmax),y2(nmax),z2(nmax),n2(nmax)
\r
2023 double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)
\r
2024 double precision u(3,3),t(3),rms,drms !armsd is real
\r
2028 ********* convert input data ***************************
\r
2029 * because L1=L2 in this special case---------->
\r
2044 n_ali=L1 !number of aligned residues
\r
2050 *** d0------------->
\r
2051 c d0=1.24*(nseqB-15)**(1.0/3.0)-1.8
\r
2053 if(d0.lt.d0_min)d0=d0_min
\r
2054 *** d0_search ----->
\r
2056 if(d0_search.gt.8)d0_search=8
\r
2057 if(d0_search.lt.4.5)d0_search=4.5
\r
2058 *** iterative parameters ----->
\r
2059 n_it=20 !maximum number of iterations
\r
2060 d_output=5 !for output alignment
\r
2061 n_init_max=6 !maximum number of L_init
\r
2064 if(n_ali.lt.4)L_ini_min=n_ali
\r
2065 do i=1,n_init_max-1
\r
2067 L_ini(n_init)=n_ali/2**(n_init-1)
\r
2068 if(L_ini(n_init).le.L_ini_min)then
\r
2069 L_ini(n_init)=L_ini_min
\r
2074 L_ini(n_init)=L_ini_min
\r
2077 ******************************************************************
\r
2078 * find the maximum score starting from local structures superposition
\r
2079 ******************************************************************
\r
2080 score_max=-1 !TM-score
\r
2081 do 333 i_init=1,n_init
\r
2082 L_init=L_ini(i_init)
\r
2083 iL_max=n_ali-L_init+1
\r
2084 do 300 iL=1,iL_max !on aligned residues, [1,nseqA]
\r
2088 k=iL+i-1 ![1,n_ali] common aligned
\r
2089 r_1(1,i)=xa(iA(k))
\r
2090 r_1(2,i)=ya(iA(k))
\r
2091 r_1(3,i)=za(iA(k))
\r
2092 r_2(1,i)=xb(iB(k))
\r
2093 r_2(2,i)=yb(iB(k))
\r
2094 r_2(3,i)=zb(iB(k))
\r
2099 call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2
\r
2100 if(i_init.eq.1)then !global superposition
\r
2101 armsd=dsqrt(rms/LL)
\r
2105 xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
\r
2106 yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
\r
2107 zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
\r
2110 call score_fun !init, get scores, n_cut+i_ali(i) for iteration
\r
2111 if(score_max.lt.score)then
\r
2115 k_ali0(i)=k_ali(i)
\r
2118 *** iteration for extending ---------------------------------->
\r
2124 m=i_ali(i) ![1,n_ali]
\r
2125 r_1(1,i)=xa(iA(m))
\r
2126 r_1(2,i)=ya(iA(m))
\r
2127 r_1(3,i)=za(iA(m))
\r
2128 r_2(1,i)=xb(iB(m))
\r
2129 r_2(2,i)=yb(iB(m))
\r
2130 r_2(3,i)=zb(iB(m))
\r
2135 call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2
\r
2137 xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
\r
2138 yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
\r
2139 zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
\r
2141 call score_fun !get scores, n_cut+i_ali(i) for iteration
\r
2142 if(score_max.lt.score)then
\r
2146 k_ali0(i)=k_ali(i)
\r
2149 if(it.eq.n_it)goto 302
\r
2150 if(n_cut.eq.ka)then
\r
2153 if(i_ali(i).eq.k_ali(i))neq=neq+1
\r
2155 if(n_cut.eq.neq)goto 302
\r
2157 301 continue !for iteration
\r
2159 300 continue !for shift
\r
2160 333 continue !for initial length, L_ali/M
\r
2162 ******** return the final rotation ****************
\r
2165 m=k_ali0(i) !record of the best alignment
\r
2166 r_1(1,i)=xa(iA(m))
\r
2167 r_1(2,i)=ya(iA(m))
\r
2168 r_1(3,i)=za(iA(m))
\r
2169 r_2(1,i)=xb(iB(m))
\r
2170 r_2(2,i)=yb(iB(m))
\r
2171 r_2(3,i)=zb(iB(m))
\r
2174 call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2
\r
2176 x1(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
\r
2177 y1(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
\r
2178 z1(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
\r
2182 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
\r
2186 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
\r
2187 c 1, collect those residues with dis<d;
\r
2188 c 2, calculate score_GDT, score_maxsub, score_TM
\r
2189 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
\r
2190 subroutine score_fun
\r
2191 PARAMETER(nmax=5000)
\r
2193 common/stru/xa(nmax),ya(nmax),za(nmax),xb(nmax),yb(nmax),zb(nmax)
\r
2194 common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
\r
2196 common/align/n_ali,iA(nmax),iB(nmax)
\r
2197 common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score
\r
2198 common/scores/score
\r
2199 double precision score
\r
2202 21 n_cut=0 !number of residue-pairs dis<d, for iteration
\r
2203 score_sum=0 !TMscore
\r
2205 i=iA(k) ![1,nseqA] reoder number of structureA
\r
2206 j=iB(k) ![1,nseqB]
\r
2207 dis=sqrt((xa(i)-xb(j))**2+(ya(i)-yb(j))**2+(za(i)-zb(j))**2)
\r
2208 if(dis.lt.d_tmp)then
\r
2210 i_ali(n_cut)=k ![1,n_ali], mark the residue-pairs in dis<d
\r
2212 score_sum=score_sum+1/(1+(dis/d0)**2)
\r
2214 if(n_cut.lt.3.and.n_ali.gt.3)then
\r
2218 score=score_sum/float(nseqB) !TM-score
\r
2223 ********************************************************************
\r
2224 * Dynamic programming for alignment.
\r
2225 * Input: score(i,j), and gap_open
\r
2226 * Output: invmap(j)
\r
2228 * Please note this subroutine is not a correct implementation of
\r
2229 * the N-W dynamic programming because the score tracks back only
\r
2230 * one layer of the matrix. This code was exploited in TM-align
\r
2231 * because it is about 1.5 times faster than a complete N-W code
\r
2232 * and does not influence much the final structure alignment result.
\r
2233 ********************************************************************
\r
2234 SUBROUTINE DP(NSEQ1,NSEQ2)
\r
2235 PARAMETER(nmax=5000)
\r
2237 common/dpc/score(nmax,nmax),gap_open,invmap(nmax)
\r
2238 dimension DIR(0:nmax,0:nmax),VAL(0:nmax,0:nmax)
\r
2241 *** initialize the matrix:
\r
2253 *** decide matrix and path:
\r
2256 D=VAL(i-1,j-1)+SCORE(i,j)
\r
2258 if(DIR(i-1,j))H=H+GAP_OPEN
\r
2260 if(DIR(i,j-1))V=V+GAP_OPEN
\r
2262 IF((D.GE.H).AND.(D.GE.V)) THEN
\r
2276 *** extract the alignment:
\r
2279 DO WHILE((i.GT.0).AND.(j.GT.0))
\r
2286 if(DIR(i-1,j))H=H+GAP_OPEN
\r
2288 if(DIR(i,j-1))V=V+GAP_OPEN
\r
2297 c^^^^^^^^^^^^^^^Dynamical programming done ^^^^^^^^^^^^^^^^^^^
\r
2301 cccccccccccccccc Calculate sum of (r_d-r_m)^2 cccccccccccccccccccccccccc
\r
2302 c w - w(m) is weight for atom pair c m (given)
\r
2303 c x - x(i,m) are coordinates of atom c m in set x (given)
\r
2304 c y - y(i,m) are coordinates of atom c m in set y (given)
\r
2305 c n - n is number of atom pairs (given)
\r
2306 c mode - 0:calculate rms only (given)
\r
2307 c 1:calculate rms,u,t (takes longer)
\r
2308 c rms - sum of w*(ux+t-y)**2 over all atom pairs (result)
\r
2309 c u - u(i,j) is rotation matrix for best superposition (result)
\r
2310 c t - t(i) is translation vector for best superposition (result)
\r
2311 c ier - 0: a unique optimal superposition has been determined(result)
\r
2312 c -1: superposition is not unique but optimal
\r
2313 c -2: no result obtained because of negative weights w
\r
2314 c or all weights equal to zero.
\r
2315 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
\r
2316 subroutine u3b(w, x, y, n, mode, rms, u, t, ier)
\r
2317 integer ip(9), ip2312(4), i, j, k, l, m1, m, ier, n, mode
\r
2318 double precision w(1), x(3, 1), y(3, 1), u(3, 3), t(3), rms, sigma
\r
2319 double precision r(3, 3), xc(3), yc(3), wc, a(3, 3), b(3, 3), e0,
\r
2320 &e(3), e1, e2, e3, d, spur, det, cof, h, g, cth, sth, sqrth, p, tol
\r
2321 &, rr(6), rr1, rr2, rr3, rr4, rr5, rr6, ss(6), ss1, ss2, ss3, ss4,
\r
2322 &ss5, ss6, zero, one, two, three, sqrt3
\r
2323 equivalence (rr(1), rr1), (rr(2), rr2), (rr(3), rr3), (rr(4), rr4)
\r
2324 &, (rr(5), rr5), (rr(6), rr6), (ss(1), ss1), (ss(2), ss2), (ss(3),
\r
2325 &ss3), (ss(4), ss4), (ss(5), ss5), (ss(6), ss6), (e(1), e1), (e(2)
\r
2326 &, e2), (e(3), e3)
\r
2327 data sqrt3 / 1.73205080756888d+00 /
\r
2328 data tol / 1.0d-2 /
\r
2329 data zero / 0.0d+00 /
\r
2330 data one / 1.0d+00 /
\r
2331 data two / 2.0d+00 /
\r
2332 data three / 3.0d+00 /
\r
2333 data ip / 1, 2, 4, 2, 3, 5, 4, 5, 6 /
\r
2334 data ip2312 / 2, 3, 1, 2 /
\r
2345 if (i .eq. j) d = one
\r
2350 c**** DETERMINE CENTROIDS OF BOTH VECTOR SETS X AND Y
\r
2352 if (n .lt. 1) return
\r
2356 if (w(m) .lt. 0.0) return
\r
2359 xc(i) = xc(i) + (w(m) * x(i,m))
\r
2360 2 yc(i) = yc(i) + (w(m) * y(i,m))
\r
2361 if (wc .le. zero) return
\r
2363 xc(i) = xc(i) / wc
\r
2364 c**** DETERMINE CORRELATION MATRIX R BETWEEN VECTOR SETS Y AND X
\r
2366 3 yc(i) = yc(i) / wc
\r
2370 e0 = e0 + (w(m) * (((x(i,m) - xc(i)) ** 2) + ((y(i,m) - yc(i)) **
\r
2373 d = w(m) * (y(i,m) - yc(i))
\r
2375 c**** CALCULATE DETERMINANT OF R(I,J)
\r
2377 4 r(i,j) = r(i,j) + (d * (x(j,m) - xc(j)))
\r
2379 det = ((r(1,1) * ((r(2,2) * r(3,3)) - (r(2,3) * r(3,2)))) - (r(1,2
\r
2380 &) * ((r(2,1) * r(3,3)) - (r(2,3) * r(3,1))))) + (r(1,3) * ((r(2,1)
\r
2381 & * r(3,2)) - (r(2,2) * r(3,1))))
\r
2382 c**** FORM UPPER TRIANGLE OF TRANSPOSED(R)*R
\r
2390 c***************** EIGENVALUES *****************************************
\r
2391 c**** FORM CHARACTERISTIC CUBIC X**3-3*SPUR*X**2+3*COF*X-DET=0
\r
2393 5 rr(m) = ((r(1,i) * r(1,j)) + (r(2,i) * r(2,j))) + (r(3,i) * r(3,j)
\r
2396 spur = ((rr1 + rr3) + rr6) / three
\r
2397 cof = ((((((rr3 * rr6) - (rr5 * rr5)) + (rr1 * rr6)) - (rr4 * rr4)
\r
2398 &) + (rr1 * rr3)) - (rr2 * rr2)) / three
\r
2403 c**** REDUCE CUBIC TO STANDARD FORM Y**3-3HY+2G=0 BY PUTTING X=Y+SPUR
\r
2405 if (spur .le. zero) goto 40
\r
2409 c**** SOLVE CUBIC. ROOTS ARE E1,E2,E3 IN DECREASING ORDER
\r
2411 g = (((spur * cof) - det) / two) - (spur * h)
\r
2413 if (h .le. zero) goto 8
\r
2415 d = ((h * h) * h) - (g * g)
\r
2416 if (d .lt. zero) d = zero
\r
2417 d = datan2(dsqrt(d),- g) / three
\r
2418 cth = sqrth * dcos(d)
\r
2419 sth = (sqrth * sqrt3) * dsin(d)
\r
2420 e1 = (spur + cth) + cth
\r
2421 e2 = (spur - cth) + sth
\r
2422 e3 = (spur - cth) - sth
\r
2423 c.....HANDLE SPECIAL CASE OF 3 IDENTICAL ROOTS
\r
2425 if (mode) 10, 50, 10
\r
2426 c**************** EIGENVECTORS *****************************************
\r
2428 8 if (mode) 30, 50, 30
\r
2430 10 do 15 l = 1, 3, 2
\r
2432 ss1 = ((d - rr3) * (d - rr6)) - (rr5 * rr5)
\r
2433 ss2 = ((d - rr6) * rr2) + (rr4 * rr5)
\r
2434 ss3 = ((d - rr1) * (d - rr6)) - (rr4 * rr4)
\r
2435 ss4 = ((d - rr3) * rr4) + (rr2 * rr5)
\r
2436 ss5 = ((d - rr1) * rr5) + (rr2 * rr4)
\r
2437 ss6 = ((d - rr1) * (d - rr3)) - (rr2 * rr2)
\r
2439 if (dabs(ss1) .ge. dabs(ss3)) goto 12
\r
2441 if (dabs(ss3) .ge. dabs(ss6)) goto 13
\r
2444 12 if (dabs(ss1) .lt. dabs(ss6)) goto 11
\r
2450 14 d = d + (ss(k) * ss(k))
\r
2451 if (d .gt. zero) d = one / dsqrt(d)
\r
2453 15 a(i,l) = a(i,l) * d
\r
2454 d = ((a(1,1) * a(1,3)) + (a(2,1) * a(2,3))) + (a(3,1) * a(3,3))
\r
2457 if ((e1 - e2) .gt. (e2 - e3)) goto 16
\r
2462 a(i,m1) = a(i,m1) - (d * a(i,m))
\r
2463 17 p = p + (a(i,m1) ** 2)
\r
2464 if (p .le. tol) goto 19
\r
2465 p = one / dsqrt(p)
\r
2467 18 a(i,m1) = a(i,m1) * p
\r
2471 if (p .lt. dabs(a(i,m))) goto 20
\r
2477 p = dsqrt((a(k,m) ** 2) + (a(l,m) ** 2))
\r
2478 if (p .le. tol) goto 40
\r
2480 a(k,m1) = - (a(l,m) / p)
\r
2481 a(l,m1) = a(k,m) / p
\r
2482 21 a(1,2) = (a(2,3) * a(3,1)) - (a(2,1) * a(3,3))
\r
2483 a(2,2) = (a(3,3) * a(1,1)) - (a(3,1) * a(1,3))
\r
2484 c****************** ROTATION MATRIX ************************************
\r
2486 a(3,2) = (a(1,3) * a(2,1)) - (a(1,1) * a(2,3))
\r
2491 b(i,l) = ((r(i,1) * a(1,l)) + (r(i,2) * a(2,l))) + (r(i,3) * a(3,l
\r
2494 31 d = d + (b(i,l) ** 2)
\r
2495 if (d .gt. zero) d = one / dsqrt(d)
\r
2497 32 b(i,l) = b(i,l) * d
\r
2498 d = ((b(1,1) * b(1,2)) + (b(2,1) * b(2,2))) + (b(3,1) * b(3,2))
\r
2501 b(i,2) = b(i,2) - (d * b(i,1))
\r
2502 33 p = p + (b(i,2) ** 2)
\r
2503 if (p .le. tol) goto 35
\r
2504 p = one / dsqrt(p)
\r
2506 34 b(i,2) = b(i,2) * p
\r
2510 if (p .lt. dabs(b(i,1))) goto 36
\r
2516 p = dsqrt((b(k,1) ** 2) + (b(l,1) ** 2))
\r
2517 if (p .le. tol) goto 40
\r
2519 b(k,2) = - (b(l,1) / p)
\r
2520 b(l,2) = b(k,1) / p
\r
2521 37 b(1,3) = (b(2,1) * b(3,2)) - (b(2,2) * b(3,1))
\r
2522 b(2,3) = (b(3,1) * b(1,2)) - (b(3,2) * b(1,1))
\r
2523 b(3,3) = (b(1,1) * b(2,2)) - (b(1,2) * b(2,1))
\r
2526 c****************** TRANSLATION VECTOR *********************************
\r
2528 39 u(i,j) = ((b(i,1) * a(j,1)) + (b(i,2) * a(j,2))) + (b(i,3) * a(j,3
\r
2531 c********************** RMS ERROR **************************************
\r
2533 41 t(i) = ((yc(i) - (u(i,1) * xc(1))) - (u(i,2) * xc(2))) - (u(i,3)
\r
2536 if (e(i) .lt. zero) e(i) = zero
\r
2537 51 e(i) = dsqrt(e(i))
\r
2539 if (e2 .le. (e1 * 1.0d-05)) ier = -1
\r
2541 if (sigma .ge. 0.0) goto 52
\r
2543 if ((e2 - e3) .le. (e1 * 1.0d-05)) ier = -1
\r
2544 52 d = (d + e2) + e1
\r
2545 rms = (e0 - d) - d
\r
2546 if (rms .lt. 0.0) rms = 0.0
\r
2548 c.....END U3B...........................................................
\r
2549 c----------------------------------------------------------
\r
2551 c----------------------------------------------------------
\r