Doing calculations

" Copyright 2003-2006 Alexander V. Diemand
This file is part of MolTalk.
MolTalk is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
MolTalk is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with MolTalk; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
"

[ |

" calculate Phi, Psi, and Chi angles of every residue in chain"

" plotting with gnuplot:
* extract the last part of the output and put it in file test.phipsi
* start gnuplot and enter: plot [-180:180] [-180:180] 'test.phipsi'
"

main

| arg strxcode chaincode |

   strxcode := nil.
   chaincode := 0.
   1 to: (ARGS count) - 1 do: [ :argc |     " ARGS is a NSArray: 9.3     "
      arg := (ARGS objectAtIndex:(argc intValue)).     " it holds objects of type NSString: 9.1     "
           " check for input file name "
      (arg = '-strx' ) ifTrue: [ targc := argc + 1. strxcode := (ARGS objectAtIndex:(targc)). ].
      (arg = '-chain' ) ifTrue: [ targc := argc + 1. chaincode := (ARGS objectAtIndex:(targc)) characterAtIndex:0. ].
   ].
   ((strxcode isNil) or: (chaincode < 10)) ifTrue: [ self usage. ]     " check for acceptable parameters "
   ifFalse: [ self continue:strxcode with:chaincode. ].     " only continue if we know the structure. "
   ^self

!
usage

   Transcript showLine: 'Calculate Phi/Psi angles around backbone of all residues in structural chain.' .
   Transcript showLine: 'Usage: MolTalk scripts/calc_phi_psi.st -strx [strxcode] -chain [chaincode]' .
   Transcript showLine: '' .

^self

!
continue:strxcode with:chaincode

| strx chain |

   Transcript showLine: ( 'got:' , strxcode).     " Transcript 9.7     "
   strx := MTStructureFactory newStructureFromPDBDirectory: strxcode.     " MTStructureFactory 3.1     "
   (strx notNil) ifTrue: [     " MTStructure 3.2     "
      Transcript showLine: 'have structure' .
      chain := strx getChain: chaincode.     " MTChain 3.4     "
      (chain notNil) ifTrue: [
         Transcript showLine: 'have chain' .

self calculatePhiPsiForChain: chain. " do the calculation for all residues in chain "

]

ifFalse: [ Transcript showLine: 'no such structure found.' ].

]

ifFalse: [ Transcript showLine: 'no such chain found.' ].

^self

!
calculatePhiPsiForChain:chain

| resphi respsi resp dist Phi Psi CAatom Catom C2atom Natom N2atom |

   resp := nil.
   resphi := NSMutableDictionary new.     " 9.5     "
   respsi := NSMutableDictionary new.
   (chain allResidues allObjects) do: [ :residue |     " MTChain: 3.4     MTResidue: 3.6     "
      (resp notNil) ifTrue: [
         atm1 := residue getCA.     " MTAtom: 3.8     "
         atm2 := resp getCA.
         dist := atm1 distanceTo: atm2.
         Transcript showLine:( ' CA-CA distance: ' ,dist stringValue).
         (dist < 4.0) ifTrue: [
                 "calculate Psi (around Ca-C)"
            CAatm := resp getCA.
            Natm := resp getAtomWithName: 'N' .
            Catm := resp getAtomWithName: 'C' .
            N2atm := residue getAtomWithName: 'N' .
            (((CAatm notNil) and: (Catm notNil)) and: ((N2atm notNil) and: (Natm notNil))) ifTrue: [
               Psi := self calculateDihedralAngleBetween:Natm and:CAatm to:Catm and:N2atm.
               Transcript showLine: ('Residue: ',resp description, ' Psi=',Psi stringValue).
               respsi setObject:Psi forKey:(resp description).
            ].
                 " calculate Phi (around Ca-N) "
            CAatm := residue getCA.
            Catm := residue getAtomWithName: 'C' .
            Natm := residue getAtomWithName: 'N' .
            C2atm := resp getAtomWithName: 'C' .
            (((CAatm notNil) and: (Catm notNil)) and: ((C2atm notNil) and: (Natm notNil))) ifTrue: [
               Phi := self calculateDihedralAngleBetween:C2atm and:Natm to:CAatm and:Catm.
               Transcript showLine: ( 'Residue: ' ,residue description, ' Phi=' ,Phi stringValue).
               resphi setObject:Phi forKey:(residue description).
            ].
         ]
         ifFalse: [ Transcript showLine: 'CA-CA distance > 4.0 A, do not assume direct bonding between consecutive amino acids in sequence.' ].
      ]     " if we know the previous residue "
      ifFalse: [ Transcript showLine: 'previous residue not found.' ].
      resp := residue.     " remember this residue in next iteration as previous residue "
   ].     " for all residues "
   Transcript showLine: 'done.' .
   Transcript showLine: '* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *' .
   (resphi allKeys) do: [:residue |
      Phi := resphi objectForKey: residue.
      Psi := respsi objectForKey: residue.
      (Psi notNil) ifTrue: [     " output only if Phi and Psi known "
         Transcript showLine: (Phi stringValue, ' ' , Psi stringValue).
      ].
   ].
   ^self

!
calculateDihedralAngleBetween:P1 and:P2 to:P3 and:P4

" This function calculates the dihedral angle between two vectors
(each given by its starting and ending point)
"

| angle V21 V23 V32 V34 N1 N2 Sign Spat |

   V21 := P1 differenceTo: P2.     " points define vectors; 4.1     "
   V23 := P3 differenceTo: P2.
   V32 := P2 differenceTo: P3.
   V34 := P4 differenceTo: P3.
   N1 := V21 vectorProductBy: V23.
   N2 := V32 vectorProductBy: V34.
   Spat := V23 mixedProductBy: N1 and: N2.
   (Spat >= 0.0) ifTrue: [ Sign := -1.0 ] ifFalse: [ Sign := 1.0 ].
   angle := (N1 angleBetween: N2) * Sign.
   ^angle

]