Difference between revisions of "Recycling Corner/Alpha Helix Generator"

From Jmol
Jump to navigation Jump to search
(moving long sections to subpages)
 
m (RIBOZOME - a Polypeptide Alpha Helix Generator script)
 
(22 intermediate revisions by 2 users not shown)
Line 1: Line 1:
== RIBOZOME - an Alpha Helix Generator script ==
+
== RIBOZOME - a Polypeptide Alpha Helix Generator script ==
 
Creates an alpha helix from a user input string.
 
Creates an alpha helix from a user input string.
 +
 +
Jmol script by Ron Mignery with help from [[User:AngelHerraez|Angel Herráez]]
  
 
__TOC__
 
__TOC__
=== Script for Jmol application ===
+
=== Script for the Jmol application ===
 +
 
 +
The top-level function plico_gen_pp asks user for a peptide sequence in a pop-up.
 +
 
 +
The top-level function plico_gen_alpha accepts a sequence string as a parameter.
 +
 
 +
''With this script you can generate a polypeptide alpha helix chain from a sequence string in 1-char amino-acid coding.  You can optionally give it a chain label other than the default :A or start at a residue number other than the default 1. You can also add to an existing chain, add new chains to an existing model or now create the chain in a new frame by prepending a '+' to the sequence string.''
 +
 
 +
''Chains start from the origin and extend along the Z-axis as they are built.  If a chain with the same chain label as the new chain is at the origin, the new sequence is added to the old.  If a chain with a different chain label is at the origin, all existing chains are shifted 20 angstroms to the right along the X axis until the origin is clear.''
 +
 
 +
Ribozome is a member of the Plico suite of protein folding tools described  [[User:Remig/plico|here]]. It may be installed and accessed as a macro with the file:
 +
<pre>Title=PLICO Generate Polypeptide
 +
Script=script <path to your script folder>/ribozome.spt;plico_gen_pp
 +
</pre> saved as ribozome.macro in your .jmol/macros directory as described in [[Macro]].
  
<pre>#  RIBOZOME - Jmol script by Ron Mignery with help from [[User:AngelHerraez|Angel Herráez]]
+
Copy and paste the following into a text editor and save in your scripts directory as ribozome.spt.
#  v1.2 beta    10/23/2013
+
<pre>#  RIBOZOME - Jmol script by Ron Mignery
#
+
#  v1.19 beta    4/12/2016 -axis is now a reserved word
#  New in v1.2 beta:
 
#      No longer zaps existing atoms
 
#      Now adds on to existing helix on subsequent runs
 
#      Fixes proline cross-link bug with XXXXPXP
 
#      Localizes variables and globalizes constants
 
 
#
 
#
 
#  RIBOZOME takes a string message encoding an amino acid (aa) sequence
 
#  RIBOZOME takes a string message encoding an amino acid (aa) sequence
Line 18: Line 28:
 
#  N terminus to the C terminus rotating the emerging helix as it goes.
 
#  N terminus to the C terminus rotating the emerging helix as it goes.
 
#
 
#
#  The message is a string entered by the user at a prompt.
+
#  The message is a string entered by the user at a prompt. <+A:><n...>[A-Z]...
 
#  It may be typed in or pasted in and be of any length
 
#  It may be typed in or pasted in and be of any length
#  If the message is prepended with <C>: (where C is any single letter)
+
#  If the message is prepended with "+" then the chain is created in a new
 +
#  frame. Otherwise it is added to the current frame
 +
#  If the message is then prepended with <C>: (where C is any single letter)
 
#  then the chain is so labeled and separated from existing chains
 
#  then the chain is so labeled and separated from existing chains
#  if different from the first chain.  
+
#  if different from the first chain.
 +
#  If the message is then prepended with digits (-0-9), that number becomes
 +
#  the first resno (even if negative).  Otherwise the first is 1.
 
#
 
#
 
#  The IUPAC/IUBMB 1 letter code is used:
 
#  The IUPAC/IUBMB 1 letter code is used:
Line 30: Line 44:
 
#  Q=GLutamiNe R=ARGinine S=SERine T=THReonine U=SElenoCysteine
 
#  Q=GLutamiNe R=ARGinine S=SERine T=THReonine U=SElenoCysteine
 
#  V=VALine W=TRyPtophan X=UNKnown Y=TYRosine Z=ASpar?X**
 
#  V=VALine W=TRyPtophan X=UNKnown Y=TYRosine Z=ASpar?X**
#    *Either GLU or GLN: drawn as GLN with chi3 flipped  
+
#    *Either GLU or GLN: drawn as GLN with chi3 flipped
 
#    **Either ASP or ASN: drawn as ASN with chi3 flipped
 
#    **Either ASP or ASN: drawn as ASN with chi3 flipped
 
#  ***Not supported: drawn as ALA
 
#  ***Not supported: drawn as ALA
  
 
# The following constant values determine the pitch of the alpha helix
 
# The following constant values determine the pitch of the alpha helix
gPHI = -57   # Dihedral angle of N-CA bond (nominally -57)
+
kPHI = -64   # Dihedral angle of N-CA bond (nominally -64)
gPSI = -47   # Dihedral angle of CA-C bond (nominally -47)
+
kPSI = -39   # Dihedral angle of CA-C bond (nominally -39)
gOMEGA = 180    # Dihedral angle of the peptide bond (nominally 180)
+
kOMEGA = 180    # Dihedral angle of the peptide bond (nominally 180)
gPEPTIDE_ANGLE = 110   # C-N-CA angle (nominally 110)
+
kPEPTIDE_ANGLE = 120   # C-N-CA angle (nominally 120)
gPRO_BUMP = -10 # Psi angle change increment to N-3psi when N is Pro
+
kPRO_BUMP = -10 # Psi angle change increment to N-3psi when N is Pro
 
gCHAIN = 'A'    # The chain id
 
gCHAIN = 'A'    # The chain id
 +
gA = ""
 +
gPdbAddHydrogens = false
 +
gAppendNew = true
 +
gNewFrame = false
 +
gSeq = ""
  
 
# Lookup 3 letter code from 1 letter code
 
# Lookup 3 letter code from 1 letter code
 
g3from1 = {"A":"ALA", "B":"GLX","C":"CYS", "D":"ASP","E":"GLU", "F":"PHE",
 
g3from1 = {"A":"ALA", "B":"GLX","C":"CYS", "D":"ASP","E":"GLU", "F":"PHE",
 
     "G":"GLY", "H":"HIS","I":"ILE", "K":"LYS","L":"LEU", "M":"MET",
 
     "G":"GLY", "H":"HIS","I":"ILE", "K":"LYS","L":"LEU", "M":"MET",
     "N":"GLN", "O":"PYL","P":"PRO", "Q":"GLN","R":"ARG", "S":"SER",
+
     "N":"ASN", "O":"PYL","P":"PRO", "Q":"GLN","R":"ARG", "S":"SER",
 
     "T":"THR", "U":"SEC","V":"VAL", "W":"TRP","X":"UNK", "Y":"TYR", "Z":"ASX"}
 
     "T":"THR", "U":"SEC","V":"VAL", "W":"TRP","X":"UNK", "Y":"TYR", "Z":"ASX"}
  
 
# Generate PDB atom record
 
# Generate PDB atom record
function genAtom(e, aa, i, xyz) {
+
function gen_atom(e, aa, i, xyz) {
     gA =  format("ATOM  %5d %4s %3s ", gN, e, aa )        
+
     gA =  format("ATOM  %5d %4s %3s ", gN, e, aa )
     gA +=  format("%s%4d    %8.3f", gCHAIN, i, xyz[1] )        
+
     gA +=  format("%s%4d    %8.3f", gCHAIN, i, xyz[1] )
 
     gA +=  format("%8.3f%8.3f\n", xyz[2], xyz[3] )
 
     gA +=  format("%8.3f%8.3f\n", xyz[2], xyz[3] )
     gn++
+
     gN++
 
     return gA
 
     return gA
 
};
 
};
  
 
# Generate a PDB amino acid record set
 
# Generate a PDB amino acid record set
function genAA(i, aa) {    # Writes globals gA and gN
+
function gen_aa(i, aa) {    # Writes globals gA and gN
  
 
     # From constructed AAs
 
     # From constructed AAs
Line 66: Line 85:
 
     var O  = [ -1.320, 1.693, 2.62 ]
 
     var O  = [ -1.320, 1.693, 2.62 ]
 
     var CB = [ 1.062, 2.1950, 0.230 ]
 
     var CB = [ 1.062, 2.1950, 0.230 ]
   
+
 
 
     var G1 = [ 2.359, 1.607, -0.344]
 
     var G1 = [ 2.359, 1.607, -0.344]
 
     var G2 = [ 1.363, 3.336, 1.157 ]
 
     var G2 = [ 1.363, 3.336, 1.157 ]
Line 73: Line 92:
 
     var E1 = [ 3.645, 3.749, -0.167 ]
 
     var E1 = [ 3.645, 3.749, -0.167 ]
 
     var E2 = [ 2.491, 3.234, -2.249 ]
 
     var E2 = [ 2.491, 3.234, -2.249 ]
     var Z  = [ 4.470, 4.717, -0.885 ]
+
     var Z  = [ 4.474, 4.857, -0.565 ]
     var H1 = [ 4.450, 6.006, -0.220 ]
+
     var H1 = [ 4.090, 6.173, -0.166 ]
     var H2 = [5.833, 4.228, -0.984 ]
+
     var H2 = [ 5.565, 4.707, -1.229 ]
      
+
 
 +
    var D1dn = [ 2.955, 2.229, -1.250 ]
 +
    var D2dn = [ 2.859, 0.552, 0.102 ]
 +
 
 +
     var E1eq = [ 3.821, 3.528, -0.382 ]
 +
    var E2eq = [ 3.337, 2.634, -2.293 ]
 +
 
 
     var Gp = [ 2.008, 1.24, -0.46 ]
 
     var Gp = [ 2.008, 1.24, -0.46 ]
     var Dp = [1.022, 0.213, -1.031 ]
+
     var Dp = [ 1.022, 0.213, -1.031 ]
   
+
 
 
     var Gfy  = [ 2.368, 1.471, -0.0152 ]
 
     var Gfy  = [ 2.368, 1.471, -0.0152 ]
 
     var D1fy = [ 3.346, 1.524, 0.921 ]
 
     var D1fy = [ 3.346, 1.524, 0.921 ]
Line 87: Line 112:
 
     var Zfy  = [ 4.588, -0.285, -0.168 ]
 
     var Zfy  = [ 4.588, -0.285, -0.168 ]
 
     var Hfy = [ 5.738, -1.245, -0.233 ]
 
     var Hfy = [ 5.738, -1.245, -0.233 ]
   
+
 
 
     var Ghw  = [ 2.406, 1.626, -0.134 ]
 
     var Ghw  = [ 2.406, 1.626, -0.134 ]
 
     var D1hw = [3.498, 1.936, 0.675]
 
     var D1hw = [3.498, 1.936, 0.675]
Line 97: Line 122:
 
     var Z3hw = [ 5.014, 2.550, 2.503 ]
 
     var Z3hw = [ 5.014, 2.550, 2.503 ]
 
     var H2hw = [ 6.153, 1.846, 1.844 ]
 
     var H2hw = [ 6.153, 1.846, 1.844 ]
   
+
 
 
     #N1 = [ 2.069, -2.122, -0.554]
 
     #N1 = [ 2.069, -2.122, -0.554]
   
+
 
 
     # Build PDB atom records common to all AAs
 
     # Build PDB atom records common to all AAs
     a3 = g3from1[aa]
+
     var a3 = g3from1[aa]
 
     if (a3 = "") {
 
     if (a3 = "") {
 
         a3 = "UNK"
 
         a3 = "UNK"
 
     }
 
     }
 
     print format("+ %s%d/%d", a3, i, gSeq.count + gResno)
 
     print format("+ %s%d/%d", a3, i, gSeq.count + gResno)
     gA = genAtom(" N  ", a3, i, N0)
+
     gA = gen_atom(" N  ", a3, i, N0)
     gA += genAtom(" CA ", a3, i, CA)
+
     gA += gen_atom(" CA ", a3, i, CA)
     gA += genAtom(" C  ", a3, i, C)
+
     gA += gen_atom(" C  ", a3, i, C)
     gA += genAtom(" O  ", a3, i, O)
+
     gA += gen_atom(" O  ", a3, i, O)
 
     if ((aa != 'G') && (aa != 'X')) {
 
     if ((aa != 'G') && (aa != 'X')) {
         gA += genAtom(" CB ", a3, i, CB)
+
         gA += gen_atom(" CB ", a3, i, CB)
 
     }
 
     }
  
Line 119: Line 144:
 
         break;
 
         break;
 
     case 'B' :
 
     case 'B' :
         gA += genAtom(" CG ", a3, i, G1)
+
         gA += gen_atom(" CG ", a3, i, G1)
         gA += genAtom(" CD ", a3, i, D1)
+
         gA += gen_atom(" CD ", a3, i, D1)
         gA += genAtom(" OE1", a3, i, E2)    # GLN with Es switched
+
         gA += gen_atom(" OE1", a3, i, E2eq)    # GLN with Es switched
         gA += genAtom(" NE2", a3, i, E1)
+
         gA += gen_atom(" NE2", a3, i, E1eq)
 
         break;
 
         break;
 
     case 'C' :
 
     case 'C' :
         gA += genAtom(" SG ", a3, i, G2)
+
         gA += gen_atom(" SG ", a3, i, G2)
 
         break;
 
         break;
 
     case 'D' :
 
     case 'D' :
         gA += genAtom(" CG ", a3, i, G1)
+
         gA += gen_atom(" CG ", a3, i, G1)
         gA += genAtom(" OD1", a3, i, D1)
+
         gA += gen_atom(" OD1", a3, i, D1dn)
         gA += genAtom(" OD2", a3, i, D2)
+
         gA += gen_atom(" OD2", a3, i, D2dn)
 
         break;
 
         break;
 
     case 'E' :
 
     case 'E' :
         gA += genAtom(" CG ", a3, i, G1)
+
         gA += gen_atom(" CG ", a3, i, G1)
         gA += genAtom(" CD ", a3, i, D1)
+
         gA += gen_atom(" CD ", a3, i, D1)
         gA += genAtom(" OE1", a3, i, E1)
+
         gA += gen_atom(" OE1", a3, i, E1eq)
         gA += genAtom(" OE2", a3, i, E2)
+
         gA += gen_atom(" OE2", a3, i, E2eq)
 
         break;
 
         break;
 
     case 'F' :
 
     case 'F' :
         gA += genAtom(" CG ", a3, i, Gfy)
+
         gA += gen_atom(" CG ", a3, i, Gfy)
         gA += genAtom(" CD1", a3, i, D1fy)
+
         gA += gen_atom(" CD1", a3, i, D1fy)
         gA += genAtom(" CD2", a3, i, D2fy)
+
         gA += gen_atom(" CD2", a3, i, D2fy)
         gA += genAtom(" CE1", a3, i, E1fy)
+
         gA += gen_atom(" CE1", a3, i, E1fy)
         gA += genAtom(" CE2", a3, i, E2fy)
+
         gA += gen_atom(" CE2", a3, i, E2fy)
         gA += genAtom(" CZ ", a3, i, Zfy)
+
         gA += gen_atom(" CZ ", a3, i, Zfy)
 
         break;
 
         break;
 
     case 'G' :
 
     case 'G' :
 
         break;
 
         break;
 
     case 'H' :
 
     case 'H' :
         gA += genAtom(" CG ", a3, i, Ghw)
+
         gA += gen_atom(" CG ", a3, i, Ghw)
         gA += genAtom(" ND1", a3, i, D1hw)
+
         gA += gen_atom(" ND1", a3, i, D1hw)
         gA += genAtom(" CD2", a3, i, D2hw)
+
         gA += gen_atom(" CD2", a3, i, D2hw)
         gA += genAtom(" CE1", a3, i, E2hw)
+
         gA += gen_atom(" CE1", a3, i, E2hw)
         gA += genAtom(" NE2", a3, i, E1hw)
+
         gA += gen_atom(" NE2", a3, i, E1hw)
 
         break;
 
         break;
 
     case 'I' :
 
     case 'I' :
         gA += genAtom(" CG1", a3, i, G1)
+
         gA += gen_atom(" CG1", a3, i, G1)
         gA += genAtom(" CG2", a3, i, G2)
+
         gA += gen_atom(" CG2", a3, i, G2)
         gA += genAtom(" CD1", a3, i, D1)
+
         gA += gen_atom(" CD1", a3, i, D1)
 
         break;
 
         break;
 
     case 'K' :
 
     case 'K' :
         gA += genAtom(" CG ", a3, i, G1)
+
         gA += gen_atom(" CG ", a3, i, G1)
         gA += genAtom(" CD ", a3, i, D1)
+
         gA += gen_atom(" CD ", a3, i, D1)
         gA += genAtom(" CE ", a3, i, E1)
+
         gA += gen_atom(" CE ", a3, i, E1)
         gA += genAtom(" NZ ", a3, i, Z)
+
         gA += gen_atom(" NZ ", a3, i, Z)
 
         break;
 
         break;
 
     case 'L' :
 
     case 'L' :
         gA += genAtom(" CG1", a3, i, G1)
+
         gA += gen_atom(" CG ", a3, i, G1)
         gA += genAtom(" CD1", a3, i, D1)
+
         gA += gen_atom(" CD1", a3, i, D1)
         gA += genAtom(" CD2", a3, i, D2)
+
         gA += gen_atom(" CD2", a3, i, D2)
 
         break;
 
         break;
 
     case 'M' :
 
     case 'M' :
         gA += genAtom(" CG ", a3, i, G1)
+
         gA += gen_atom(" CG ", a3, i, G1)
         gA += genAtom(" SD ", a3, i, D1)
+
         gA += gen_atom(" SD ", a3, i, D1)
         gA += genAtom(" CE ", a3, i, E1)
+
         gA += gen_atom(" CE ", a3, i, E1)
 
         break;
 
         break;
 
     case 'N' :
 
     case 'N' :
         gA += genAtom(" CG ", a3, i, G1)
+
         gA += gen_atom(" CG ", a3, i, G1)
         gA += genAtom(" OD1", a3, i, D1)
+
         gA += gen_atom(" OD1", a3, i, D1dn)
         gA += genAtom(" ND2", a3, i, D2)
+
         gA += gen_atom(" ND2", a3, i, D2dn)
 
         break;
 
         break;
 
     case 'P' :
 
     case 'P' :
         gA += genAtom(" CG ", a3, i, GP)
+
         gA += gen_atom(" CG ", a3, i, GP)
         gA += genAtom(" CD ", a3, i, DP)
+
         gA += gen_atom(" CD ", a3, i, DP)
 
         break;
 
         break;
 
     case 'Q' :
 
     case 'Q' :
         gA += genAtom(" CG ", a3, i, G1)
+
         gA += gen_atom(" CG ", a3, i, G1)
         gA += genAtom(" CD ", a3, i, D1)
+
         gA += gen_atom(" CD ", a3, i, D1)
         gA += genAtom(" OE1", a3, i, E1)
+
         gA += gen_atom(" OE1", a3, i, E1eq)
         gA += genAtom(" NE2", a3, i, E2)
+
         gA += gen_atom(" NE2", a3, i, E2eq)
 
         break;
 
         break;
 
     case 'R' :
 
     case 'R' :
         gA += genAtom(" CG ", a3, i, G1)
+
         gA += gen_atom(" CG ", a3, i, G1)
         gA += genAtom(" CD ", a3, i, D1)
+
         gA += gen_atom(" CD ", a3, i, D1)
         gA += genAtom(" NE ", a3, i, E1)
+
         gA += gen_atom(" NE ", a3, i, E1)
         gA += genAtom(" CZ ", a3, i, Z)
+
         gA += gen_atom(" CZ ", a3, i, Z)
         gA += genAtom(" NH1", a3, i, H1)
+
         gA += gen_atom(" NH1", a3, i, H1)
         gA += genAtom(" NH2", a3, i, H2)
+
         gA += gen_atom(" NH2", a3, i, H2)
 
         break;
 
         break;
 
     case 'S' :
 
     case 'S' :
         gA += genAtom(" OG ", a3, i, G1)
+
         gA += gen_atom(" OG ", a3, i, G1)
 
         break;
 
         break;
 
     case 'T' :
 
     case 'T' :
         gA += genAtom(" OG1", a3, i, G1)
+
         gA += gen_atom(" OG1", a3, i, G1)
         gA += genAtom(" CG2", a3, i, G2)
+
         gA += gen_atom(" CG2", a3, i, G2)
 
         break;
 
         break;
 
     case 'U' :
 
     case 'U' :
         gA += genAtom("SeG ", a3, i, G1)
+
         gA += gen_atom("SeG ", a3, i, G1)
 
         break;
 
         break;
 
     case 'V' :
 
     case 'V' :
         gA += genAtom(" CG1", a3, i, G1)
+
         gA += gen_atom(" CG1", a3, i, G1)
         gA += genAtom(" CG2", a3, i, G2)
+
         gA += gen_atom(" CG2", a3, i, G2)
 
         break;
 
         break;
 
     case 'W' :
 
     case 'W' :
         gA += genAtom(" CG ", a3, i, Ghw)
+
         gA += gen_atom(" CG ", a3, i, Ghw)
         gA += genAtom(" CD1", a3, i, D1hw)
+
         gA += gen_atom(" CD2", a3, i, D1hw)
         gA += genAtom(" CD2", a3, i, D2hw)
+
         gA += gen_atom(" CD1", a3, i, D2hw)
         gA += genAtom(" CE2", a3, i, E2hw)
+
         gA += gen_atom(" CE2", a3, i, E2hw)
         gA += genAtom(" NE1", a3, i, E1hw)
+
         gA += gen_atom(" NE1", a3, i, E1hw)
         gA += genAtom(" CE3", a3, i, E3hw)
+
         gA += gen_atom(" CE3", a3, i, E3hw)
         gA += genAtom(" CZ2", a3, i, Z2hw)
+
         gA += gen_atom(" CZ2", a3, i, Z2hw)
         gA += genAtom(" CZ3", a3, i, Z3hw)
+
         gA += gen_atom(" CZ3", a3, i, Z3hw)
         gA += genAtom(" CH2", a3, i, H2hw)
+
         gA += gen_atom(" CH2", a3, i, H2hw)
 
         break;
 
         break;
 
     case 'X' :
 
     case 'X' :
         gA += genAtom(" Xx ", a3, i, CB)
+
         gA += gen_atom(" Xx ", a3, i, CB)
 
         break;
 
         break;
 
     case 'Y' :
 
     case 'Y' :
         gA += genAtom(" CG ", a3, i, Gfy)
+
         gA += gen_atom(" CG ", a3, i, Gfy)
         gA += genAtom(" CD1", a3, i, D1fy)
+
         gA += gen_atom(" CD1", a3, i, D1fy)
         gA += genAtom(" CD2", a3, i, D2fy)
+
         gA += gen_atom(" CD2", a3, i, D2fy)
         gA += genAtom(" CE1", a3, i, E1fy)
+
         gA += gen_atom(" CE1", a3, i, E1fy)
         gA += genAtom(" CE2", a3, i, E2fy)
+
         gA += gen_atom(" CE2", a3, i, E2fy)
         gA += genAtom(" CZ ", a3, i, Zfy)
+
         gA += gen_atom(" CZ ", a3, i, Zfy)
         gA += genAtom(" OH ", a3, i, Hfy)
+
         gA += gen_atom(" OH ", a3, i, Hfy)
 
         break;
 
         break;
 
     case 'Z' :
 
     case 'Z' :
         gA += genAtom(" CG ", a3, i, G1)
+
         gA += gen_atom(" CG ", a3, i, G1)
         gA += genAtom(" OD1", a3, i, D2)    # ASN with Ds switched
+
         gA += gen_atom(" OD1", a3, i, D2dn)    # ASN with Ds switched
         gA += genAtom(" ND2", a3, i, D1)
+
         gA += gen_atom(" ND2", a3, i, D1dn)
 
         break;
 
         break;
 
     default :
 
     default :
Line 249: Line 274:
  
 
# Generate an alpha helix
 
# Generate an alpha helix
function genAlpha(gSeq) {
+
function plico_gen_alpha(seq) {
 
+
    try {
     gN = all.count    + 1 # global new N atom index
+
        gPdbAddHydrogens = pdbAddHydrogens
 
+
        set pdbAddHydrogens false
    # Find last linkable N if any
+
        if (gPlicoRecord != "") {
    gResno = 0    # global pre-existing AA count
+
            var g = format("show file \"%s\"", gPlicoRecord)
    var pn = 1    # previous gN
+
            var ls = script(g)
     for (var i = all.count-1; i > 0; i--) {
+
            if (ls.find("FileNotFoundException")) {
 +
                ls = ""
 +
            }
 +
            ls += format("plico_gen_alpha(\"%s\");", seq)
 +
            write var ls @gPlicoRecord
 +
        }
 +
      
 +
        seq = seq%9999%0
 +
        gNewFrame = false
 +
        if (seq[1] == '+') {
 +
            gNewFrame = true
 +
            seq = seq[2][9999]
 +
        }
 +
        if (seq[2] == ':') {
 +
            gCHAIN = seq[1]
 +
            seq = seq[3][9999]
 +
        }
 +
        var resStart = ""
 +
        while ("_0123456789-".find(seq[1]) >  0) {
 +
            resStart += ""+seq[1]
 +
            seq = seq[2][9999]
 +
        }
 +
        gSeq = seq
 +
   
 +
        var c = gCHAIN
 +
        if (resStart.size > 0) {
 +
            gResno = -1 + resStart    # global pre-existing AA count
 +
        }
 +
        else {
 +
            gResno = 0    # global pre-existing AA count
 +
        }
 +
        var pn = 1    # previous gN
 +
        gAppendNew = appendNew
 +
        gN = 1
 +
        if (gNewFrame) {
 +
            appendNew = true
 +
        }
 +
        else {
 +
            appendNew = false
 +
   
 +
            # If not new
 +
            if (m > 0) {
 +
   
 +
                # Get the largest atomno in frame
 +
                gN = {thisModel}.atomno.max
 +
   
 +
                # While there is an atom at the origin
 +
                while (true) {
 +
                    var ai = {within(0.1, {0 0 0}) and thisModel}
 +
                    if (ai) {
 +
   
 +
                        # If on the same chain
 +
                        if (ai[1].chain = c) {
 +
   
 +
                            # Delete OXT
 +
                            delete {(atomName="OXT") and thisModel
 +
                                and (chain=c)}
 +
                            gResno = ai.resno
 +
                            pn = ai.atomno
 +
                            break
 +
                        }
 +
                        # Else move all from the new chain rightward on X axis
 +
                        else {
 +
                            select {thisModel and (chain!=c)}
 +
                            translateselected {20, 0, 0 }
 +
                            gN++
 +
                        }
 +
                    }
 +
                    else {
 +
                        break
 +
                    }
 +
                } # endwhile
 +
            }
 +
        }
 +
      
 +
        # For each aa
 +
        var nn = gN    # new N
 +
        for (var i = 1; i <= gSeq.count; i++) {
 +
   
 +
            if ((m > 0) and not appendNew) {
 
      
 
      
        # If found
+
                # Move polypeptide C to bond distance from new AA N
        if (distance({atomno=i}, {0,0,0})  < 0.1) {
+
                select {thisModel and (chain=c)}
            pn = i
+
                 fix none
           
 
            # If new chain, separate from existing chain
 
            if ({atomno=i}.chain != gCHAIN) {
 
                 select all
 
 
                 translateselected {2.069, -2.122, -0.554 } #N1
 
                 translateselected {2.069, -2.122, -0.554 } #N1
 
             }
 
             }
             else {
+
   
                 gResno = {atomno=i}.resno
+
            # Gen AA
 +
            gA = "data \"append aa\"\n" # global PDB atom record
 +
            gA += gen_aa(i + gResno, gSeq[i]) # gN is updated in subroutine
 +
            gA += "end \"append aa\""
 +
            script inline @{gA}
 +
   
 +
            var f = (_frameID/1000000)
 +
            var m = (_frameID%1000000)
 +
            appendNew = false
 +
   
 +
            # If PRO ahead
 +
            var pb = 0
 +
            if ((gSeq.count - i) >= 2) {
 +
                if (gSeq[i + 2] == 'P') {
 +
                    pb = kPRO_BUMP
 +
                }
 +
            }
 +
   
 +
            # If not first AA
 +
             if (nn > 1) {
 +
   
 +
                # Gen axis on new N perpendicular to the plane
 +
                # containing atoms nn, pn+2 and nn+1
 +
                var aNn = {(atomno=nn) and thisModel and (chain=c)}
 +
                 var aNn1 = {(atomno=@{nn+1}) and thisModel and (chain=c)}
 +
                var aPn1 = {(atomno=@{pn+1}) and thisModel and (chain=c)}
 +
                var aPn2 = {(atomno=@{pn+2}) and thisModel and (chain=c)}
 +
                var aPn3 = {(atomno=@{pn+3}) and thisModel and (chain=c)}
 +
                var v1=aPn2.xyz - aNn.xyz
 +
                var v2=aNn1.xyz - aNn.xyz
 +
                var caxis = cross(v1, v2)
 +
   
 +
                # Center on atom previous C
 +
                caxis += aPn2.xyz
 +
   
 +
                # Rotate the polypeptide N on the new AA C to tetrahedral (nominally 110)
 +
                select {(atomno < nn) and thisModel and (chain=c)}
 +
                fix {(atomno >= nn) and thisModel and (chain=c)}
 +
                rotateselected @caxis @aNn @{kPEPTIDE_ANGLE - 65.5}
 +
   
 +
                # Make omega dihedral = kOMEGA (nominally 180)
 +
                rotateselected @aPn2 @aNn @{kOMEGA - 157.8}
 +
   
 +
                # Make the new phi dihedral = kPHI (nominally -60)
 +
                rotateselected @aNn @aNn1 @{-kPHI - 116.5}
 +
   
 +
                # Make the old psi dihedral = kPSI (nominally -50)
 +
                select {(atomno < nn) and thisModel and (chain=c)
 +
                    and not aPn2 and not aPn3}
 +
                rotateselected @aPn1 @aPn2 @{-kPSI - 60.6 + pb}
 +
   
 +
                # Make the peptide bond
 +
                connect @aNn @aPn2
 +
                refresh
 
             }
 
             }
            break;
+
      
        }
+
            # Step new and previous N
     }
+
            pn = nn
 
+
            nn = gN
    # For each aa
+
      
    set appendnew false
+
         } # endfor i
    var nn = gN   # new N
+
   
     for (var i = 1; i <= gSeq.count; i++) {
+
         # Add terminal O
 
 
        # Move polypeptide C to bond distance from new AA N
 
         select all
 
        fix none
 
        translateselected {2.069, -2.122, -0.554 } #N1
 
 
 
         # Gen AA
 
 
         gA = "data \"append aa\"\n"    # global PDB atom record
 
         gA = "data \"append aa\"\n"    # global PDB atom record
         gA += genAA(i + gResno, gSeq[i]);    # gN is updated in subroutine
+
         gA += gen_atom(" OXT", g3from1[gSeq[gResno+gSeq.count]],
 +
            gResno+gSeq.count, [ -2.142, 2.057, 0.574])
 
         gA += "end \"append aa\""
 
         gA += "end \"append aa\""
 
         script inline @{gA}
 
         script inline @{gA}
 
+
        refresh
         # If PRO ahead
+
   
         var pb = 0
+
         connect
        if ((gSeq.count - i) >= 2) {
+
         var xx = {(element="Xx") and thisModel and (chain=c)}
            if (gSeq[i + 2] == 'P') {
+
         for (var i = 1; i <= xx.size; i++) {
                pb = gPRO_BUMP
+
              connect 1.8 {(atomindex=@{xx[i].atomIndex}) and thisModel and (chain=c)}
            }
 
         }
 
 
 
        # If not first AA
 
        if (nn > 1) {
 
 
 
            # Gen axis on new N perpendicular to the plane
 
            # containing atoms nn, pn+2 and nn+1
 
            var v1={atomno = @{pn+2}}.xyz - {atomno = nn}.xyz
 
            var v2={atomno = @{nn+1}}.xyz - {atomno = nn}.xyz
 
            var axis = cross(v1, v2)
 
           
 
            # Center on atom previous C
 
            axis += {atomno = @{pn+2}}.xyz
 
 
 
            # Rotate the polypeptide N on the new AA C to tetrahedral (nominally 110)
 
            select atomno < nn
 
            fix atomno >= nn
 
            rotateselected @axis {atomno = nn} @{gPEPTIDE_ANGLE - 69.3}
 
           
 
            # Make omega dihedral = gOMEGA (nominally 180)
 
            rotateselected {atomno=@{pn+2}} {atomno=nn} @{gOMEGA - 147.4}
 
 
 
            # Make the new phi dihedral = gPHI (nominally -57)
 
            rotateselected {atomno = nn} {atomno = @{nn+1}} @{gPHI - 8.7}
 
 
 
            # Make the old psi dihedral = gPSI (nominally -47)
 
            select atomno < @{pn+2} && atomno != @{pn+3}
 
            rotateselected {atomno=@{pn+1}} {atomno=@{pn+2}} @{gPSI + 33.4 + pb}
 
 
         }
 
         }
          
+
   
         # Step new and previous N
+
         # Clean up
         pn = nn
+
         connect ([UNK].CA) ([UNK].Xx and within(group, _1) and thisModel and (chain=c))
         nn = gN
+
        select (thisModel)
 
+
        fix none
        # Make the peptide bond
+
         print format("%d atoms generated", gN-1)
         connect
+
        appendNew = gAppendNew
 +
         pdbAddHydrogens = gPdbAddHydrogens
 +
    }
 +
    catch {
 +
         print "error"
 
     }
 
     }
   
 
    # Clean up
 
    connect ([UNK].CA) ([UNK].Xx and within(group, _1))
 
    select all
 
    fix none
 
    print format("%d atoms generated", gN)
 
 
}
 
}
  
echo Generating Alpha Helix
+
function plico_gen_pp {
 +
    print "Generating Alpha Helix..."
  
# Get the sequence from the user
+
    # Get the sequence from the user
gSeq = prompt("Enter AA sequence (1 char coded)", "")%9999%0
+
    gSeq = prompt("Enter AA sequence (<+A:><n...>[A-Z]...)", "")%9999%0
if (gSeq.count > 0) {
+
    if ((gSeq != "NULL") and (gSeq.count > 0)) {
    if (gSeq[2] == ':') {
+
        print format ("Sequence=%s  phi=%d  psi=%d", gSeq, kPHI, kPSI)
        gCHAIN = gSeq[1]
+
         plico_gen_alpha(gSeq)
         gSeq[1] = ' '
 
        gSeq[2] = ' '
 
        gSeq = gSeq%0
 
 
     }
 
     }
    print format ("Sequence=%s  phi=%d  psi=%d", gSeq, gPHI, gPSI)
 
    print format ("chain=%s peptide angle=%d  pro bump=%d", gCHAIN, gPEPTIDE_ANGLE, gPRO_BUMP)
 
    genAlpha(gSeq)
 
 
}
 
}
</pre>
+
# end of ribozome.spt</pre>
  
 
=== In a webpage ===
 
=== In a webpage ===
Using a J(S)mol object.
+
''(Updated to RIBOZOME v1.16 script, now works both in Jmol-Java and JSmol-HTML5)''
 +
 
 +
Sequence input is managed from the webpage. In addition, the helix parameters (which are fixed in the original script) may be changed in the page by a user.
 +
The model is generated and displayed in a Jmol object within the page (JmolApplet or JSmol HTML5 object).
  
(coming soon...)
+
See it [http://biomodel.uah.es/Jmol/helix_builder/ in action]

Latest revision as of 19:34, 9 October 2019

RIBOZOME - a Polypeptide Alpha Helix Generator script

Creates an alpha helix from a user input string.

Jmol script by Ron Mignery with help from Angel Herráez

Script for the Jmol application

The top-level function plico_gen_pp asks user for a peptide sequence in a pop-up.

The top-level function plico_gen_alpha accepts a sequence string as a parameter.

With this script you can generate a polypeptide alpha helix chain from a sequence string in 1-char amino-acid coding. You can optionally give it a chain label other than the default :A or start at a residue number other than the default 1. You can also add to an existing chain, add new chains to an existing model or now create the chain in a new frame by prepending a '+' to the sequence string.

Chains start from the origin and extend along the Z-axis as they are built. If a chain with the same chain label as the new chain is at the origin, the new sequence is added to the old. If a chain with a different chain label is at the origin, all existing chains are shifted 20 angstroms to the right along the X axis until the origin is clear.

Ribozome is a member of the Plico suite of protein folding tools described here. It may be installed and accessed as a macro with the file:

Title=PLICO Generate Polypeptide
Script=script <path to your script folder>/ribozome.spt;plico_gen_pp

saved as ribozome.macro in your .jmol/macros directory as described in Macro.

Copy and paste the following into a text editor and save in your scripts directory as ribozome.spt.

#   RIBOZOME - Jmol script by Ron Mignery
#   v1.19 beta    4/12/2016 -axis is now a reserved word
#
#   RIBOZOME takes a string message encoding an amino acid (aa) sequence
#   and generates a corresponding alpha helix one aa at a time from the
#   N terminus to the C terminus rotating the emerging helix as it goes.
#
#   The message is a string entered by the user at a prompt. <+A:><n...>[A-Z]...
#   It may be typed in or pasted in and be of any length
#   If the message is prepended with "+" then the chain is created in a new
#   frame. Otherwise it is added to the current frame
#   If the message is then prepended with <C>: (where C is any single letter)
#   then the chain is so labeled and separated from existing chains
#   if different from the first chain.
#   If the message is then prepended with digits (-0-9), that number becomes
#   the first resno (even if negative).  Otherwise the first is 1.
#
#   The IUPAC/IUBMB 1 letter code is used:
#   A=ALAnine B=GLutam?X* C=CYSteine D=ASPartate E=GLUtamate
#   F=PHEnylalanine G=GLYcine H=HIStidine I=IsoLEucine K=LYSine
#   L=LEUcine M=METhionine N=ASparagiNe O=PYrroLysine*** P=PROline
#   Q=GLutamiNe R=ARGinine S=SERine T=THReonine U=SElenoCysteine
#   V=VALine W=TRyPtophan X=UNKnown Y=TYRosine Z=ASpar?X**
#     *Either GLU or GLN: drawn as GLN with chi3 flipped
#    **Either ASP or ASN: drawn as ASN with chi3 flipped
#   ***Not supported: drawn as ALA

# The following constant values determine the pitch of the alpha helix
kPHI = -64    # Dihedral angle of N-CA bond (nominally -64)
kPSI = -39    # Dihedral angle of CA-C bond (nominally -39)
kOMEGA = 180    # Dihedral angle of the peptide bond (nominally 180)
kPEPTIDE_ANGLE = 120    # C-N-CA angle (nominally 120)
kPRO_BUMP = -10 # Psi angle change increment to N-3psi when N is Pro
gCHAIN = 'A'    # The chain id
gA = ""
gPdbAddHydrogens = false
gAppendNew = true
gNewFrame = false
gSeq = ""

# Lookup 3 letter code from 1 letter code
g3from1 = {"A":"ALA", "B":"GLX","C":"CYS", "D":"ASP","E":"GLU", "F":"PHE",
    "G":"GLY", "H":"HIS","I":"ILE", "K":"LYS","L":"LEU", "M":"MET",
    "N":"ASN", "O":"PYL","P":"PRO", "Q":"GLN","R":"ARG", "S":"SER",
    "T":"THR", "U":"SEC","V":"VAL", "W":"TRP","X":"UNK", "Y":"TYR", "Z":"ASX"}

# Generate PDB atom record
function gen_atom(e, aa, i, xyz) {
    gA =  format("ATOM  %5d %4s %3s ", gN, e, aa )
    gA +=  format("%s%4d    %8.3f", gCHAIN, i, xyz[1] )
    gA +=  format("%8.3f%8.3f\n", xyz[2], xyz[3] )
    gN++
    return gA
};

# Generate a PDB amino acid record set
function gen_aa(i, aa) {    # Writes globals gA and gN

    # From constructed AAs
    var N0 = [0.0, 0.0, 0.0]
    var CA = [ 0.200, 1.174, 0.911 ]
    var C  = [ -1.110, 1.668, 1.425 ]
    var O  = [ -1.320, 1.693, 2.62 ]
    var CB = [ 1.062, 2.1950, 0.230 ]

    var G1 = [ 2.359, 1.607, -0.344]
    var G2 = [ 1.363, 3.336, 1.157 ]
    var D1 = [ 3.222, 2.656, -1.048 ]
    var D2 = [ 3.143, 0.904, 0.725 ]
    var E1 = [ 3.645, 3.749, -0.167 ]
    var E2 = [ 2.491, 3.234, -2.249 ]
    var Z  = [ 4.474, 4.857, -0.565 ]
    var H1 = [ 4.090, 6.173, -0.166 ]
    var H2 = [ 5.565, 4.707, -1.229 ]

    var D1dn = [ 2.955, 2.229, -1.250 ]
    var D2dn = [ 2.859, 0.552, 0.102 ]

    var E1eq = [ 3.821, 3.528, -0.382 ]
    var E2eq = [ 3.337, 2.634, -2.293 ]

    var Gp = [ 2.008, 1.24, -0.46 ]
    var Dp = [ 1.022, 0.213, -1.031 ]

    var Gfy  = [ 2.368, 1.471, -0.0152 ]
    var D1fy = [ 3.346, 1.524, 0.921 ]
    var D2fy = [ 2.493, 0.516, -1.151 ]
    var E1fy = [ 4.513, 0.615, 0.8244 ]
    var E2fy = [ 3.528, -0.336, -1.206 ]
    var Zfy  = [ 4.588, -0.285, -0.168 ]
    var Hfy = [ 5.738, -1.245, -0.233 ]

    var Ghw  = [ 2.406, 1.626, -0.134 ]
    var D1hw = [3.498, 1.936, 0.675]
    var D2hw = [ 2.713, 0.901, -1.211 ]
    var E1hw = [ 4.160, 0.518, -1.178 ]
    var E2hw = [ 4.622, 1.160, 0.0816 ]
    var E3hw = [ 3.789, 2.523, 1.944 ]
    var Z2hw = [ 5.973, 1.177, 0.689 ]
    var Z3hw = [ 5.014, 2.550, 2.503 ]
    var H2hw = [ 6.153, 1.846, 1.844 ]

    #N1 = [ 2.069, -2.122, -0.554]

    # Build PDB atom records common to all AAs
    var a3 = g3from1[aa]
    if (a3 = "") {
        a3 = "UNK"
    }
    print format("+ %s%d/%d", a3, i, gSeq.count + gResno)
    gA = gen_atom(" N  ", a3, i, N0)
    gA += gen_atom(" CA ", a3, i, CA)
    gA += gen_atom(" C  ", a3, i, C)
    gA += gen_atom(" O  ", a3, i, O)
    if ((aa != 'G') && (aa != 'X')) {
        gA += gen_atom(" CB ", a3, i, CB)
    }

    # Now add AA specific atom records
    switch (aa) {
    case 'A' :
        break;
    case 'B' :
        gA += gen_atom(" CG ", a3, i, G1)
        gA += gen_atom(" CD ", a3, i, D1)
        gA += gen_atom(" OE1", a3, i, E2eq)    # GLN with Es switched
        gA += gen_atom(" NE2", a3, i, E1eq)
        break;
    case 'C' :
        gA += gen_atom(" SG ", a3, i, G2)
        break;
    case 'D' :
        gA += gen_atom(" CG ", a3, i, G1)
        gA += gen_atom(" OD1", a3, i, D1dn)
        gA += gen_atom(" OD2", a3, i, D2dn)
        break;
    case 'E' :
        gA += gen_atom(" CG ", a3, i, G1)
        gA += gen_atom(" CD ", a3, i, D1)
        gA += gen_atom(" OE1", a3, i, E1eq)
        gA += gen_atom(" OE2", a3, i, E2eq)
        break;
    case 'F' :
        gA += gen_atom(" CG ", a3, i, Gfy)
        gA += gen_atom(" CD1", a3, i, D1fy)
        gA += gen_atom(" CD2", a3, i, D2fy)
        gA += gen_atom(" CE1", a3, i, E1fy)
        gA += gen_atom(" CE2", a3, i, E2fy)
        gA += gen_atom(" CZ ", a3, i, Zfy)
        break;
    case 'G' :
        break;
    case 'H' :
        gA += gen_atom(" CG ", a3, i, Ghw)
        gA += gen_atom(" ND1", a3, i, D1hw)
        gA += gen_atom(" CD2", a3, i, D2hw)
        gA += gen_atom(" CE1", a3, i, E2hw)
        gA += gen_atom(" NE2", a3, i, E1hw)
        break;
    case 'I' :
        gA += gen_atom(" CG1", a3, i, G1)
        gA += gen_atom(" CG2", a3, i, G2)
        gA += gen_atom(" CD1", a3, i, D1)
        break;
    case 'K' :
        gA += gen_atom(" CG ", a3, i, G1)
        gA += gen_atom(" CD ", a3, i, D1)
        gA += gen_atom(" CE ", a3, i, E1)
        gA += gen_atom(" NZ ", a3, i, Z)
        break;
    case 'L' :
        gA += gen_atom(" CG ", a3, i, G1)
        gA += gen_atom(" CD1", a3, i, D1)
        gA += gen_atom(" CD2", a3, i, D2)
        break;
    case 'M' :
        gA += gen_atom(" CG ", a3, i, G1)
        gA += gen_atom(" SD ", a3, i, D1)
        gA += gen_atom(" CE ", a3, i, E1)
        break;
    case 'N' :
        gA += gen_atom(" CG ", a3, i, G1)
        gA += gen_atom(" OD1", a3, i, D1dn)
        gA += gen_atom(" ND2", a3, i, D2dn)
        break;
    case 'P' :
        gA += gen_atom(" CG ", a3, i, GP)
        gA += gen_atom(" CD ", a3, i, DP)
        break;
    case 'Q' :
        gA += gen_atom(" CG ", a3, i, G1)
        gA += gen_atom(" CD ", a3, i, D1)
        gA += gen_atom(" OE1", a3, i, E1eq)
        gA += gen_atom(" NE2", a3, i, E2eq)
        break;
    case 'R' :
        gA += gen_atom(" CG ", a3, i, G1)
        gA += gen_atom(" CD ", a3, i, D1)
        gA += gen_atom(" NE ", a3, i, E1)
        gA += gen_atom(" CZ ", a3, i, Z)
        gA += gen_atom(" NH1", a3, i, H1)
        gA += gen_atom(" NH2", a3, i, H2)
        break;
    case 'S' :
        gA += gen_atom(" OG ", a3, i, G1)
        break;
    case 'T' :
        gA += gen_atom(" OG1", a3, i, G1)
        gA += gen_atom(" CG2", a3, i, G2)
        break;
    case 'U' :
        gA += gen_atom("SeG ", a3, i, G1)
        break;
    case 'V' :
        gA += gen_atom(" CG1", a3, i, G1)
        gA += gen_atom(" CG2", a3, i, G2)
        break;
    case 'W' :
        gA += gen_atom(" CG ", a3, i, Ghw)
        gA += gen_atom(" CD2", a3, i, D1hw)
        gA += gen_atom(" CD1", a3, i, D2hw)
        gA += gen_atom(" CE2", a3, i, E2hw)
        gA += gen_atom(" NE1", a3, i, E1hw)
        gA += gen_atom(" CE3", a3, i, E3hw)
        gA += gen_atom(" CZ2", a3, i, Z2hw)
        gA += gen_atom(" CZ3", a3, i, Z3hw)
        gA += gen_atom(" CH2", a3, i, H2hw)
        break;
    case 'X' :
        gA += gen_atom(" Xx ", a3, i, CB)
        break;
    case 'Y' :
        gA += gen_atom(" CG ", a3, i, Gfy)
        gA += gen_atom(" CD1", a3, i, D1fy)
        gA += gen_atom(" CD2", a3, i, D2fy)
        gA += gen_atom(" CE1", a3, i, E1fy)
        gA += gen_atom(" CE2", a3, i, E2fy)
        gA += gen_atom(" CZ ", a3, i, Zfy)
        gA += gen_atom(" OH ", a3, i, Hfy)
        break;
    case 'Z' :
        gA += gen_atom(" CG ", a3, i, G1)
        gA += gen_atom(" OD1", a3, i, D2dn)    # ASN with Ds switched
        gA += gen_atom(" ND2", a3, i, D1dn)
        break;
    default :
        break;
    }

    return gA
};

# Generate an alpha helix
function plico_gen_alpha(seq) {
    try {
        gPdbAddHydrogens = pdbAddHydrogens
        set pdbAddHydrogens false
        if (gPlicoRecord != "") {
            var g = format("show file \"%s\"", gPlicoRecord)
            var ls = script(g)
            if (ls.find("FileNotFoundException")) {
                ls = ""
            }
            ls += format("plico_gen_alpha(\"%s\");", seq)
            write var ls @gPlicoRecord
        }
    
        seq = seq%9999%0
        gNewFrame = false
        if (seq[1] == '+') {
            gNewFrame = true
            seq = seq[2][9999]
        }
        if (seq[2] == ':') {
            gCHAIN = seq[1]
            seq = seq[3][9999]
        }
        var resStart = ""
        while ("_0123456789-".find(seq[1]) >  0) {
            resStart += ""+seq[1]
            seq = seq[2][9999]
        }
        gSeq = seq
    
        var c = gCHAIN
        if (resStart.size > 0) {
            gResno = -1 + resStart    # global pre-existing AA count
        }
        else {
            gResno = 0    # global pre-existing AA count
        }
        var pn = 1    # previous gN
        gAppendNew = appendNew
        gN = 1
        if (gNewFrame) {
            appendNew = true
        }
        else {
            appendNew = false
    
            # If not new
            if (m > 0) {
    
                # Get the largest atomno in frame
                gN = {thisModel}.atomno.max
    
                # While there is an atom at the origin
                while (true) {
                    var ai = {within(0.1, {0 0 0}) and thisModel}
                    if (ai) {
    
                        # If on the same chain
                        if (ai[1].chain = c) {
    
                            # Delete OXT
                            delete {(atomName="OXT") and thisModel
                                and (chain=c)}
                            gResno = ai.resno
                            pn = ai.atomno
                            break
                        }
                        # Else move all from the new chain rightward on X axis
                        else {
                            select {thisModel and (chain!=c)}
                            translateselected {20, 0, 0 }
                            gN++
                        }
                    }
                    else {
                        break
                    }
                } # endwhile
            }
        }
    
        # For each aa
        var nn = gN    # new N
        for (var i = 1; i <= gSeq.count; i++) {
    
            if ((m > 0) and not appendNew) {
    
                # Move polypeptide C to bond distance from new AA N
                select {thisModel and (chain=c)}
                fix none
                translateselected {2.069, -2.122, -0.554 } #N1
            }
    
            # Gen AA
            gA = "data \"append aa\"\n" # global PDB atom record
            gA += gen_aa(i + gResno, gSeq[i]) # gN is updated in subroutine
            gA += "end \"append aa\""
            script inline @{gA}
    
            var f = (_frameID/1000000)
            var m = (_frameID%1000000)
            appendNew = false
    
            # If PRO ahead
            var pb = 0
            if ((gSeq.count - i) >= 2) {
                if (gSeq[i + 2] == 'P') {
                    pb = kPRO_BUMP
                }
            }
    
            # If not first AA
            if (nn > 1) {
    
                # Gen axis on new N perpendicular to the plane
                # containing atoms nn, pn+2 and nn+1
                var aNn = {(atomno=nn) and thisModel and (chain=c)}
                var aNn1 = {(atomno=@{nn+1}) and thisModel and (chain=c)}
                var aPn1 = {(atomno=@{pn+1}) and thisModel and (chain=c)}
                var aPn2 = {(atomno=@{pn+2}) and thisModel and (chain=c)}
                var aPn3 = {(atomno=@{pn+3}) and thisModel and (chain=c)}
                var v1=aPn2.xyz - aNn.xyz
                var v2=aNn1.xyz - aNn.xyz
                var caxis = cross(v1, v2)
    
                # Center on atom previous C
                caxis += aPn2.xyz
    
                # Rotate the polypeptide N on the new AA C to tetrahedral (nominally 110)
                select {(atomno < nn) and thisModel and (chain=c)}
                fix {(atomno >= nn) and thisModel and (chain=c)}
                rotateselected @caxis @aNn @{kPEPTIDE_ANGLE - 65.5}
    
                # Make omega dihedral = kOMEGA (nominally 180)
                rotateselected @aPn2 @aNn @{kOMEGA - 157.8}
    
                # Make the new phi dihedral = kPHI (nominally -60)
                rotateselected @aNn @aNn1 @{-kPHI - 116.5}
    
                # Make the old psi dihedral = kPSI (nominally -50)
                select {(atomno < nn) and thisModel and (chain=c)
                    and not aPn2 and not aPn3}
                rotateselected @aPn1 @aPn2 @{-kPSI - 60.6 + pb}
    
                # Make the peptide bond
                connect @aNn @aPn2
                refresh
            }
    
            # Step new and previous N
            pn = nn
            nn = gN
    
        } # endfor i
    
        # Add terminal O
        gA = "data \"append aa\"\n"    # global PDB atom record
        gA += gen_atom(" OXT", g3from1[gSeq[gResno+gSeq.count]],
            gResno+gSeq.count, [ -2.142, 2.057, 0.574])
        gA += "end \"append aa\""
        script inline @{gA}
        refresh
    
        connect
        var xx = {(element="Xx") and thisModel and (chain=c)}
        for (var i = 1; i <= xx.size; i++) {
              connect 1.8 {(atomindex=@{xx[i].atomIndex}) and thisModel and (chain=c)}
        }
    
        # Clean up
        connect ([UNK].CA) ([UNK].Xx and within(group, _1) and thisModel and (chain=c))
        select (thisModel)
        fix none
        print format("%d atoms generated", gN-1)
        appendNew = gAppendNew
        pdbAddHydrogens = gPdbAddHydrogens
    }
    catch {
        print "error"
    }
}

function plico_gen_pp {
    print "Generating Alpha Helix..."

    # Get the sequence from the user
    gSeq = prompt("Enter AA sequence (<+A:><n...>[A-Z]...)", "")%9999%0
    if ((gSeq != "NULL") and (gSeq.count > 0)) {
        print format ("Sequence=%s  phi=%d  psi=%d", gSeq, kPHI, kPSI)
        plico_gen_alpha(gSeq)
    }
}
# end of ribozome.spt

In a webpage

(Updated to RIBOZOME v1.16 script, now works both in Jmol-Java and JSmol-HTML5)

Sequence input is managed from the webpage. In addition, the helix parameters (which are fixed in the original script) may be changed in the page by a user. The model is generated and displayed in a Jmol object within the page (JmolApplet or JSmol HTML5 object).

See it in action

Contributors

Remig, AngelHerraez