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

From Jmol
Jump to navigation Jump to search
(Script for Jmol application)
m (RIBOZOME - a Polypeptide Alpha Helix Generator script)
 
(14 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]]
+
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 plicoGenPP asks user for a peptide sequence in a pop-up.
+
The top-level function plico_gen_pp asks user for a peptide sequence in a pop-up.
  
The top-level function plicoGenAlpha accepts a sequence string as a parameter.
+
The top-level function plico_gen_alpha accepts a sequence string as a parameter.
  
Ribozome is a member of the Plico suite of protein folding tools. It may be installed and accessed as a macro with the file:
+
''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
 
<pre>Title=PLICO Generate Polypeptide
Script=script <path to your script folder>/ribozome.spt;plicogenpp
+
Script=script <path to your script folder>/ribozome.spt;plico_gen_pp
</pre> saved as ribozome.macro in your .jmol/macros folder as described in[[Macro]].
+
</pre> 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.
 
<pre>#  RIBOZOME - Jmol script by Ron Mignery
 
<pre>#  RIBOZOME - Jmol script by Ron Mignery
#  v1.7 beta    1/23/2014 - add plico recording
+
#  v1.19 beta    4/12/2016 -axis is now a reserved word
 
#
 
#
 
#  RIBOZOME takes a string message encoding an amino acid (aa) sequence
 
#  RIBOZOME takes a string message encoding an amino acid (aa) sequence
Line 23: 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 35: 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
var kPHI = -57   # Dihedral angle of N-CA bond (nominally -57)
+
kPHI = -64   # Dihedral angle of N-CA bond (nominally -64)
var kPSI = -47   # Dihedral angle of CA-C bond (nominally -47)
+
kPSI = -39   # Dihedral angle of CA-C bond (nominally -39)
var kOMEGA = 180    # Dihedral angle of the peptide bond (nominally 180)
+
kOMEGA = 180    # Dihedral angle of the peptide bond (nominally 180)
var kPEPTIDE_ANGLE = 110   # C-N-CA angle (nominally 110)
+
kPEPTIDE_ANGLE = 120   # C-N-CA angle (nominally 120)
var kPRO_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
var gCHAIN = 'A'    # The chain id
+
gCHAIN = 'A'    # The chain id
var gA = ""
+
gA = ""
 +
gPdbAddHydrogens = false
 +
gAppendNew = true
 +
gNewFrame = false
 +
gSeq = ""
  
 
# Lookup 3 letter code from 1 letter code
 
# Lookup 3 letter code from 1 letter code
var 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++
Line 64: Line 77:
  
 
# 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 72: 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 82: Line 95:
 
     var H1 = [ 4.090, 6.173, -0.166 ]
 
     var H1 = [ 4.090, 6.173, -0.166 ]
 
     var H2 = [ 5.565, 4.707, -1.229 ]
 
     var H2 = [ 5.565, 4.707, -1.229 ]
   
+
 
 
     var D1dn = [ 2.955, 2.229, -1.250 ]
 
     var D1dn = [ 2.955, 2.229, -1.250 ]
 
     var D2dn = [ 2.859, 0.552, 0.102 ]
 
     var D2dn = [ 2.859, 0.552, 0.102 ]
   
+
 
 
     var E1eq = [ 3.821, 3.528, -0.382 ]
 
     var E1eq = [ 3.821, 3.528, -0.382 ]
 
     var E2eq = [ 3.337, 2.634, -2.293 ]
 
     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 99: 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 109: 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
 
     var a3 = g3from1[aa]
 
     var a3 = g3from1[aa]
Line 118: Line 131:
 
     }
 
     }
 
     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 131: 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, E2eq)    # GLN with Es switched
+
         gA += gen_atom(" OE1", a3, i, E2eq)    # GLN with Es switched
         gA += genAtom(" NE2", a3, i, E1eq)
+
         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, D1dn)
+
         gA += gen_atom(" OD1", a3, i, D1dn)
         gA += genAtom(" OD2", a3, i, D2dn)
+
         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, E1eq)
+
         gA += gen_atom(" OE1", a3, i, E1eq)
         gA += genAtom(" OE2", a3, i, E2eq)
+
         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, D1dn)
+
         gA += gen_atom(" OD1", a3, i, D1dn)
         gA += genAtom(" ND2", a3, i, D2dn)
+
         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, E1eq)
+
         gA += gen_atom(" OE1", a3, i, E1eq)
         gA += genAtom(" NE2", a3, i, E2eq)
+
         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, D2dn)    # ASN with Ds switched
+
         gA += gen_atom(" OD1", a3, i, D2dn)    # ASN with Ds switched
         gA += genAtom(" ND2", a3, i, D1dn)
+
         gA += gen_atom(" ND2", a3, i, D1dn)
 
         break;
 
         break;
 
     default :
 
     default :
Line 261: Line 274:
  
 
# Generate an alpha helix
 
# Generate an alpha helix
function plicoGenAlpha(gSeq) {
+
function plico_gen_alpha(seq) {
     set pdbAddHydrogens FALSE   
+
     try {
    if (gPlicoRecord != "") {
+
        gPdbAddHydrogens = pdbAddHydrogens
        var g = format("show file \"%s\"", gPlicoRecord)
+
        set pdbAddHydrogens false
        var ls = script(g)
+
        if (gPlicoRecord != "") {
        if (ls.find("FileNotFoundException")) {
+
            var g = format("show file \"%s\"", gPlicoRecord)
             ls = ""
+
            var ls = script(g)
 +
            if (ls.find("FileNotFoundException")) {
 +
                ls = ""
 +
            }
 +
             ls += format("plico_gen_alpha(\"%s\");", seq)
 +
            write var ls @gPlicoRecord
 
         }
 
         }
        ls += format("plicoGenAlpha(\"%s\");", gSeq)
 
        write var ls @gPlicoRecord
 
    }
 
 
      
 
      
    gSeq = gSeq%9999%0
+
        seq = seq%9999%0
    if (gSeq[2] == ':') {
+
        gNewFrame = false
        gCHAIN = gSeq[1]
+
        if (seq[1] == '+') {
         gSeq[1] = ' '
+
            gNewFrame = true
        gSeq[2] = ' '
+
            seq = seq[2][9999]
         gSeq = gSeq%0
+
        }
     }
+
        if (seq[2] == ':') {
 
+
            gCHAIN = seq[1]
    gN = all.count + 1 # global new N atom index
+
            seq = seq[3][9999]
 
+
        }
    # Find last linkable N if any
+
        var resStart = ""
    gResno = 0    # global pre-existing AA count
+
         while ("_0123456789-".find(seq[1]) >  0) {
    var pn = 1    # previous gN
+
            resStart += ""+seq[1]
     for (var i = all.count-1; i > 0; i--) {
+
            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
           
 
            # Delete O
 
            delete {atomIndex=@{all.atomIndex.max}}
 
           
 
            # 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
 +
                }
 
             }
 
             }
             break;
+
   
        }
+
             # If not first AA
     }
+
            if (nn > 1) {
 
+
      
    # For each aa
+
                # Gen axis on new N perpendicular to the plane
    set appendnew false
+
                # containing atoms nn, pn+2 and nn+1
    var nn = gN    # new N
+
                var aNn = {(atomno=nn) and thisModel and (chain=c)}
    for (var i = 1; i <= gSeq.count; i++) {
+
                var aNn1 = {(atomno=@{nn+1}) and thisModel and (chain=c)}
 
+
                var aPn1 = {(atomno=@{pn+1}) and thisModel and (chain=c)}
        # Move polypeptide C to bond distance from new AA N
+
                var aPn2 = {(atomno=@{pn+2}) and thisModel and (chain=c)}
        select all
+
                var aPn3 = {(atomno=@{pn+3}) and thisModel and (chain=c)}
        fix none
+
                var v1=aPn2.xyz - aNn.xyz
        translateselected {2.069, -2.122, -0.554 } #N1
+
                var v2=aNn1.xyz - aNn.xyz
 
+
                var caxis = cross(v1, v2)
         # Gen AA
+
   
 +
                # 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 = "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 = kPRO_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} @{kPEPTIDE_ANGLE - 69.3}
 
           
 
            # Make omega dihedral = kOMEGA (nominally 180)
 
            rotateselected {atomno=@{pn+2}} {atomno=nn} @{kOMEGA - 147.4}
 
 
 
            # Make the new phi dihedral = kPHI (nominally -57)
 
            rotateselected {atomno = nn} {atomno = @{nn+1}} @{kPHI - 8.7}
 
 
 
            # Make the old psi dihedral = kPSI (nominally -47)
 
            select atomno < nn && atomno != @{pn+2} && atomno != @{pn+3}
 
            rotateselected {atomno=@{pn+1}} {atomno=@{pn+2}} @{kPSI + 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"
 
     }
 
     }
   
 
    # Add terminal O
 
    gA = "data \"append aa\"\n"    # global PDB atom record
 
    gA += genAtom(" OXT", g3from1[gSeq[gResno+gSeq.count]],
 
        gResno+gSeq.count, [ -2.142, 2.057, 0.574])
 
    gA += "end \"append aa\""
 
    script inline @{gA}
 
    connect
 
   
 
    # Clean up
 
    connect ([UNK].CA) ([UNK].Xx and within(group, _1))
 
    select all
 
    fix none
 
    print format("%d atoms generated", gN-1)
 
 
}
 
}
  
 +
function plico_gen_pp {
 +
    print "Generating Alpha Helix..."
  
function plicoGenPP {
 
    echo Generating Alpha Helix
 
    kPHI = -57
 
    kPSI = -47
 
    kOMEGA = 180
 
    kPEPTIDE_ANGLE = 110
 
    kPRO_BUMP = -10
 
    gCHAIN = 'A'
 
    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":"GLN", "O":"PYL","P":"PRO", "Q":"GLN","R":"ARG", "S":"SER",
 
        "T":"THR", "U":"SEC","V":"VAL", "W":"TRP","X":"UNK", "Y":"TYR", "Z":"ASX"}
 
   
 
 
     # 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 != "NULL") and (gSeq.count > 0)) {
 
     if ((gSeq != "NULL") and (gSeq.count > 0)) {
 
         print format ("Sequence=%s  phi=%d  psi=%d", gSeq, kPHI, kPSI)
 
         print format ("Sequence=%s  phi=%d  psi=%d", gSeq, kPHI, kPSI)
         plicoGenAlpha(gSeq)
+
         plico_gen_alpha(gSeq)
 
     }
 
     }
 
}
 
}
# end of ribozome.spt
+
# end of ribozome.spt</pre>
</pre>
 
  
 
=== In a webpage ===
 
=== In a webpage ===
The main script is a small adaptation of the one above.
+
''(Updated to RIBOZOME v1.16 script, now works both in Jmol-Java and JSmol-HTML5)''
Sequence input has been moved to 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 in the page (JmolApplet or JSmol HTML5 object).
 
  
The JSmol variant is not working yet. We need debugging the script, something in it fails to run properly in JSmol and no good helix is formed.
+
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 [http://biomodel.uah.es/Jmol/helix_builder/ in action]
 
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