PROGRAM gethel c find helices in a protein implicit none integer*4 nat,iag,iac,iah,iap,iax,ia,ix,iargc real*8,allocatable::r(:) integer*4,allocatable::ity(:),ipep(:,:),bt(:,:),nt(:),us(:), 1ipeh(:,:),ipec(:,:),ipex(:,:),ippp(:,:) character*80 fn WRITE(6,64000) 64000 FORMAT(' GETHEL: FINDs HELICES IN A PROTEIN',/) c if(iargc().ne.1)then write(6,*)'usage: gethel ' stop endif call getarg(1,fn) OPEN(20,FILE=fn,STATUS='OLD') read(20,*) read(20,*)nat allocate(r(3*nat),ity(nat),us(nat),bt(nat,7),nt(nat), 1ipep(nat,8),ipeh(nat,8),ipec(nat,8),ipex(nat,8), 1ippp(nat,8)) do 41 ia=1,nat read(20,*)ity(ia),(r(3*(ia-1)+ix),ix=1,3),(bt(ia,ix),ix=1,7) nt(ia)=7 us(ia)=1 do 411 ix=1,7 if(bt(ia,ix).eq.0)then nt(ia)=ix-1 goto 41 endif 411 continue 41 continue close(20) c identify amide groups call getamides(nat,ity,us,nt,bt,r,iag,iac,iah,iap,iax, 1ipep,ipeh,ipec,ipex,ippp) c find secondary structure call dohel(iag,ipep,nat,r,ippp,iap,iac,ipec) end subroutine dohel(iag,ipep,nat,r,ippp,iap,iac,ipec) implicit none integer*4 ia,iag,nat,ipep(nat,8),ica,icn,ib,ifr,ifa,igc,ift,ix, 1infa,iap,ippp(nat,8),iagg,iac,ipec(nat,8),i0,i1,i2,i3,i4,i5,i6, 1i,in,io,ih,nh,ic,ih0,iam real*8 r(*),fi,psi,w,ta,fi0(5),psi0(5),w0(5),di,dist,dd, 1xo,yo,zo,xn,yn,zn c standard angle values for c 1 alpha helix c 2 3_10 helix c 3 beta sheet c 4 PPII c 5 PPI data fi0/-57.0d0,-60.0d0,-140.0d0,-75.0d0,-75.0d0/ data psi0/-47.0d0,-30.0d0, 130.0d0,150.0d0,150.0d0/ data w0/180.0d0,180.0d0, 180.0d0,180.0d0, 0.0d0/ character*5 styp(6) data styp/'alpha','3_10 ','beta ',' PII ',' PI ',' XXX '/ integer*4,allocatable::iact(:),iant(:),iaff(:),nati(:), 1itpp(:,:),ipt(:),ist(:),nath(:),lat(:,:) parameter (ih0=10) character*10 HC(ih0) data HC/' H-TX ',' H-TY ',' H-TZ ',' H-RX ', 1 ' H-RY ',' H-RZ ',' H-stretch',' H-breath.', 1 ' H-def. ',' H-unwind.'/ c c iagg - number of all recognized amide groups: c ipt - type of amide group c ipt: c number of ordinary amide groups: iag 1 c number of terminal amide groups: iah 2 c number of proline amide groups: iap 3 c number of strange amide groups: iax 4 c number of cis amide groups: iac c allocate(iact(nat),iant(nat),iaff(nat),nati(nat), 1itpp(nat,8),ipt(nat),ist(nat),nath(nat),lat(nat,nat)) c helix coordinates: c put trans cis and proline amide groups together iagg=0 do 101 ia=1,iag iagg=iagg+1 do 102 ix=1,8 102 itpp(iagg,ix)=ipep(ia,ix) 101 ipt(iagg)=1 do 111 ia=1,iac iagg=iagg+1 do 112 ix=1,8 112 itpp(iagg,ix)=ipec(ia,ix) 111 ipt(iagg)=1 do 105 ia=1,iap iagg=iagg+1 do 106 ix=1,8 106 itpp(iagg,ix)=ippp(ia,ix) 105 ipt(iagg)=3 c find if amide groups continues write(6,6001) 6001 format(' aminogroup type N-end C-end (to am.gr. #)') do 1 ia=1,iagg c c-alpha and C on nitrogen: ica=itpp(ia,6) icn=itpp(ia,5) c c C(5) - N(3)H(4) - [c(2)=O(1)] - c(6) - c icn ica c (iant) (iact) c continuation on N and C ends: c iant(ia)=0 iact(ia)=0 iaff(ia)=0 do 107 ib=1,iagg if(ib.ne.ia)then if(icn.eq.itpp(ib,6))iant(ia)=icn if(ica.eq.itpp(ib,5))then iact(ia)=ica iaff(ia)=ib endif endif 107 continue 1 write(6,6002)ia,ipt(ia),iant(ia),iact(ia),iaff(ia) 6002 format(4i12,'(',i4,')') open(33,file='HELIX.LOG') nh=0 c ifr ... number of the chain: ifr=0 do 2 ia=1,iagg if(iant(ia).eq.0)then ifr=ifr+1 write(6,*)ia,'start' c N-terminal is free, initiate chain c c igc,ia current amide group to add: igc=ia c ifa ... number of amadie groups already in the chain: 77 ifa=0 c ift ... number of atoms in the chain: ift=0 3333 ifa=ifa+1 c nati atom list in the fragment: c add C-(N-H)-C=O- c 5 3 4 2 1 do 3 ix=1,5 ift=ift+1 3 nati(ift)=itpp(igc,ix) c number of the following aminoacid: infa=iaff(igc) write(6,*)'following',infa,iact(igc),ifa if(ifa.gt.1)then i1=nati(ift-5) i2=nati(ift-7) i3=nati(ift-8) i4=nati(ift ) i5=nati(ift-2) i6=nati(ift-3) w =ta(i1,i2,i3,i4,r) psi=ta(i2,i3,i4,i5,r) fi =ta(i3,i4,i5,i6,r) c determine sec structure between am groups ia ifa, save in ist: i0=6 di=1.0d9 do 4 i=1,5 dist=dsqrt(dd(fi,fi0(i))+dd(psi,psi0(i))+dd(w,w0(i))) if(dist.lt.di.and.dist.lt.40.0d0)then i0=i di=dist endif 4 continue write(6,600)fi,psi,w,styp(i0) 600 format(' fi psi w:',3f10.1,1x,a5) ist(ifa)=i0 endif c look for helical loops: io=nati(ift-4) xo=r(1+3*(io-1)) yo=r(2+3*(io-1)) zo=r(3+3*(io-1)) if(ifa.gt.2)then in=nati(ift-12) xn=r(1+3*(in-1)) yn=r(2+3*(in-1)) zn=r(3+3*(in-1)) dist=dsqrt((xo-xn)**2+(yo-yn)**2+(zo-zn)**2) if(dist.lt.3.3d0)then write(6,*)'3-10 loop' c change types back: ist(ifa )=2 ist(ifa-1)=2 ist(ifa-2)=2 endif endif if(ifa.gt.3)then in=nati(ift-17) xn=r(1+3*(in-1)) yn=r(2+3*(in-1)) zn=r(3+3*(in-1)) dist=dsqrt((xo-xn)**2+(yo-yn)**2+(zo-zn)**2) if(dist.lt.3.3d0)then write(6,*)'alpha helical loop' c change types back: ist(ifa )=1 ist(ifa-1)=1 ist(ifa-2)=1 ist(ifa-3)=1 endif endif c capture end of helix HHXX : if(ifa.gt.4)then if(( ist(ifa ).ne.1.and.ist(ifa ).ne.2 .and. 1 ist(ifa-1).ne.1.and.ist(ifa-1).ne.2).and. 2 ((ist(ifa-2).eq.1.or. ist(ifa-2).eq.2).and. 3 (ist(ifa-3).eq.1.or. ist(ifa-3).eq.2)))then nh=nh+1 nath(nh)=0 do 8 iam=1,ifa if(ist(iam).eq.1.or.ist(iam).eq.2)then lat(nh,nath(nh)+1)=nati(1+5*(iam-1)) lat(nh,nath(nh)+2)=nati(2+5*(iam-1)) lat(nh,nath(nh)+3)=nati(3+5*(iam-1)) lat(nh,nath(nh)+4)=nati(4+5*(iam-1)) lat(nh,nath(nh)+5)=nati(5+5*(iam-1)) nath(nh)=nath(nh)+5 endif 8 continue write(6,331)nh,nath(nh) 331 format(' Helix number',i6,',',i6,' atoms:') do 5 ix=1,nath(nh) write(6,332)lat(nh,ix) 332 format(i6,$) 5 if(mod(ix,15).eq.0)write(6,*) write(6,*) c start a new fragment if(iact(igc).ne.0.and.infa.ne.0)then igc=infa ifr=ifr+1 goto 77 endif endif endif if(iact(igc).eq.0.or.infa.eq.0)then c terminate chain building goto 222 else igc=infa ifr=ifr+1 goto 3333 endif c iact(igc)=0 222 if(ifa.gt.3)then write(6,*)' chain ',ifr,',',ifa,' amides,',ift,' atoms:' do 7 ix=1,ift write(6,332)nati(ix) 7 if(mod(ix,15).eq.0)write(6,*) write(6,*) endif endif c iant(igc)=0 2 continue c write to file open(9,file='HEL.TXT') do 6 ih=1,nh write(9,80)ih 80 format(i4) c helical/cigar FILE.UMA format do 6 ic=1,ih0 WRITE(9,7778)13+ic,nath(ih),nath(ih),HC(ic),ic 7778 FORMAT(3I6,4(5x,'0'),4x,'1.',6('0'),A10,3x,'0 ',3x,'0, coord',I6) c write the atom list only for the first coordinate: 6 if(ic.eq.1)write(9,7772)(lat(ih,ix),ix=1,nath(ih)) 7772 format(12i6) close(9) write(6,*)'HEL.TXT written' return end function ta(l1,l2,l3,l4,r) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4(I-N) DIMENSION XX(3),YY(3),A(3),P(3),E0(3),B0(3), 1AP(3),r(*) PI=3.141592653589D0 do 23 ix=1,3 A(ix) =r(ix+3*(l3-1))-r(ix+3*(l4-1)) P(ix) =r(ix+3*(l1-1))-r(ix+3*(l2-1)) 23 E0(ix)=r(ix+3*(l2-1))-r(ix+3*(l3-1)) call vp(B0,E0,A) SSB0B0=sp(B0,B0) SSE0E0=sp(E0,E0) DO 8 KK=1,3 B0(KK)=B0(KK)/DSQRT(SSB0B0) 8 E0(KK)=E0(KK)/DSQRT(SSE0E0) SSAE0=sp(A,E0) SSPE0=sp(P,E0) DO 9 KK=1,3 XX(KK)=A(KK)-E0(KK)*SSAE0 9 YY(KK)=P(KK)-E0(KK)*SSPE0 SSXXXX=sp(XX,XX) SSYYYY=sp(YY,YY) SSXXYY=sp(XX,YY) call vp(AP,XX,YY) SSE0AP=sp(AP,E0) SIGN=1.0D0 IF (dabs(SSE0AP).gt.1.0D-9)SIGN=-SSE0AP/ABS(SSE0AP) CA=0.0D0 IF (dabs(SSXXXX).gt.1.0D-9.AND.dabs(SSYYYY).gt.1.0D-9) 1CA=-SSXXYY/DSQRT(SSXXXX)/DSQRT(SSYYYY) ANGLE=0.0D0 IF (CA.LE.-1.0D0) ANGLE=PI IF (ABS(CA).LT.1.0D0) ANGLE=ACOS(CA) ANGLE=180.0D0* ANGLE/PI ANGLE=SIGN*ANGLE if(ANGLE.lt.-180.0d0)ANGLE= 360.0d0+ANGLE if(ANGLE.gt.+180.0d0)ANGLE=-360.0d0+ANGLE ta=ANGLE return end subroutine getamides(nat,ity,us,nt,bt,r,iag,iac,iah,iap,iax, 1ipep,ipeh,ipec,ipex,ippp) implicit none integer*4 nat,ity(*),us(*),nt(*),bt(nat,7),ix, 1iag,iac,iah,iap,iax,ia,ia2,ia3,ia4,ia5,ia6,ia7,ia8,ia9, 1ix2,ix3,ix4,ix5,ix6,ix7,ix8,ix9,ipep(nat,8),ipeh(nat,8), 1ipec(nat,8),ipex(nat,8),ippp(nat,8) logical lcc,lcon real*8 r(*),xi,yi,zi,xj,yj,zj,d2 c number of ordinary trans amide groups: iag=0 c number of ordinary cis amide groups: iac=0 c number of terminal amide groups: iah=0 c number of proline amide groups: iap=0 c number of strange amide groups: iax=0 do 9 ia=1,nat if(ity(ia).eq.8.and.us(ia).eq.1.and.nt(ia).eq.1)then do 124 ix2=1,nt(ia) ia2=bt(ia,ix2) if(ity(ia2).eq.6.and.us(ia2).eq.1)then do 125 ix3=1,nt(ia2) ia3=bt(ia2,ix3) if(ity(ia3).eq.7.and.us(ia3).eq.1)then c ia ia2 ia3 c O = C - N c c carbon on C: lcc=.false. do 135 ix6=1,nt(ia2) ia6=bt(ia2,ix6) if(ity(ia6).eq.6)then lcc=.true. goto 1003 endif 135 continue 1003 continue c do 126 ix4=1,nt(ia3) ia4=bt(ia3,ix4) if(ity(ia4).eq.1)then c ia ia2 ia3 ia4 c O = C - N - H xi=r(3*(ia -1)+1) yi=r(3*(ia -1)+2) zi=r(3*(ia -1)+3) xj=r(3*(ia4-1)+1) yj=r(3*(ia4-1)+2) zj=r(3*(ia4-1)+3) d2=(xi-xj)**2+(yi-yj)**2+(zi-zj)**2 if(d2.gt.7.0d0)then c trans do 127 ix5=1,nt(ia3) ia5=bt(ia3,ix5) if(ity(ia5).eq.6.and.ia5.ne.ia2)then c O = C - N(-C) - H if(lcc)then c O = C(-C) - N(-C) - H call swallow(nat,us,iag,ipep, 1 ia,ia2,ia3,ia4,ia5,ia6,ia6,ia6) goto 9 endif endif if(ity(ia5).eq.1.and.ia5.ne.ia4)then c O = C - NH2 if(lcc)then c O = C(-C) - NH2 call swallow(nat,us,iah,ipeh, 1 ia,ia2,ia3,ia4,ia5,ia6,ia6,ia6) goto 9 endif endif 127 continue else c cis do 138 ix5=1,nt(ia3) ia5=bt(ia3,ix5) if(ity(ia5).eq.6.and.ia5.ne.ia2)then c O = C - N(-C) - H if(lcc)then c O = C(-C) - N(-C) - H call swallow(nat,us,iac,ipec, 1 ia,ia2,ia3,ia4,ia5,ia6,ia6,ia6) goto 9 endif endif 138 continue endif endif 126 continue do 129 ix4=1,nt(ia3) ia4=bt(ia3,ix4) if(ity(ia4).eq.6.and.ia4.ne.ia2)then c ia ia2 ia3 ia4 c O = C - N - C do 131 ix5=1,nt(ia3) ia5=bt(ia3,ix5) if(ity(ia5).eq.6.and.ia5.ne.ia2.and.ia5.ne.ia4)then c O = C - N(-C) - C c ensure that ia5 is the main chain continuation carbon: lcon=.false. do 136 ix7=1,nt(ia5) ia7=bt(ia5,ix7) if(ity(ia7).eq.6)then do 137 ix8=1,nt(ia7) ia8=bt(ia7,ix8) 137 if(ity(ia8).eq.8)lcon=.true. endif 136 continue if(.not.lcon)then ix7=ia4 ia4=ia5 ia5=ix7 endif if(lcc)then c O = C(-C) - N(-C) - C do 132 ix7=1,nt(ia5) ia7=bt(ia5,ix7) if(ity(ia7).eq.6)then do 133 ix8=1,nt(ia7) ia8=bt(ia7,ix8) if(ity(ia8).eq.6)then do 134 ix9=1,nt(ia8) ia9=bt(ia8,ix9) if(ia9.eq.ia4)then c proline amide call swallow(nat,us,iap,ippp, 1 ia,ia2,ia3,ia4,ia5,ia6,ia7,ia8) goto 9 endif 134 continue endif 133 continue endif 132 continue endif endif 131 continue endif 129 continue c unknown O = C - N group call swallow(nat,us,iax,ipex,ia,ia2,ia3,0,0,0,0,0) endif 125 continue endif 124 continue endif 9 continue if(iag.gt.0)then open(11,file='AMIDE.LST') write( 6,*)iag,' trans O = C(-C) - N(-C) - H amide groups:' write(11,*)iag,' trans O = C(-C) - N(-C) - H amide groups:' do 4 ia=1,iag write(11,6009)ia,(ipep(ia,ix),ix=1,8) 4 write(6,6009)ia,(ipep(ia,ix),ix=1,6) 6009 format(i4,4x,11i5) if(iac.eq.0)close(11) endif if(iac.gt.0)then open(11,file='AMIDE.LST') write( 6,*)iac,' cis O = C(-C) - N(-C) - H amide groups:' write(11,*)iac,' cis O = C(-C) - N(-C) - H amide groups:' do 71 ia=1,iac write(11,6009)ia,(ipec(ia,ix),ix=1,8) 71 write(6,6009)ia,(ipec(ia,ix),ix=1,6) close(11) endif if(iap.gt.0)then open(12,file='PROLINE.LST') write(6,*)iap,' proline O = C(-C) - N (amide) groups:' write(12,*)iap,' proline O = C(-C) - N (amide) groups:' do 11 ia=1,iap write(12,6009)ia,(ippp(ia,ix),ix=1,8) 11 write(6,6009)ia,(ippp(ia,ix),ix=1,3) close(12) endif if(iah.gt.0)then open(13,file='TERMINAL.LST') write(6,*)iah,' terminal O = C(-C) - NH2 amide groups:' write(13,*)iah,' terminal O = C(-C) - NH2 amide groups:' do 12 ia=1,iah write(13,6009)ia,(ipeh(ia,ix),ix=1,8) 12 write(6,6009)ia,(ipeh(ia,ix),ix=1,3) close(13) endif if(iax.gt.0)then open(14,file='UNKNOWN.LST') write(6,*)iax,' unknown O = C - N (amide) groups:' write(14,*)iax,' unknown O = C - N (amide) groups:' do 10 ia=1,iax write(14,6009)ia,(ipex(ia,ix),ix=1,8) 10 write(6,6009)ia,(ipex(ia,ix),ix=1,3) close(14) endif return end subroutine swallow(n,u,i,j,i1,i2,i3,i4,i5,i6,i7,i8) implicit none integer n,i,j(n,8),i1,i2,i3,i4,i5,i6,i7,i8,u(*),ii i=i+1 j(i,1)=i1 j(i,2)=i2 j(i,3)=i3 j(i,4)=i4 j(i,5)=i5 j(i,6)=i6 j(i,7)=i7 j(i,8)=i8 do 1 ii=1,8 1 u(j(i,ii))=0 return end function sp(a,b) IMPLICIT none real*8 a(*),b(*),sp sp=a(1)*b(1)+a(2)*b(2)+a(3)*b(3) return end subroutine vp(A,X,Y) implicit none real*8 A(*),X(*),Y(*) A(1)=X(2)*Y(3)-X(3)*Y(2) A(2)=X(3)*Y(1)-X(1)*Y(3) A(3)=X(1)*Y(2)-X(2)*Y(1) return end function dd(a,b) implicit none real*8 a,b,dd,d d=a-b if(dabs(d-360.0d0).lt.dabs(d))d=d-360.0d0 if(dabs(d+360.0d0).lt.dabs(d))d=d+360.0d0 dd=d**2 return end