Difference between revisions of "Recycling Corner/Alpha Helix Generator"
(Add Plico recording support) |
m (→Ribozome.spt, a script for the Jmol application) |
||
Line 17: | Line 17: | ||
<pre># RIBOZOME - Jmol script by Ron Mignery | <pre># RIBOZOME - Jmol script by Ron Mignery | ||
− | # v1. | + | # v1.8 beta 2/7/2014 - fix ASN mislabel |
# | # | ||
# RIBOZOME takes a string message encoding an amino acid (aa) sequence | # RIBOZOME takes a string message encoding an amino acid (aa) sequence | ||
Line 51: | Line 51: | ||
var g3from1 = {"A":"ALA", "B":"GLX","C":"CYS", "D":"ASP","E":"GLU", "F":"PHE", | var 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":" | + | "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"} | ||
Line 179: | Line 179: | ||
break; | break; | ||
case 'L' : | case 'L' : | ||
− | gA += genAtom(" | + | gA += genAtom(" CG ", a3, i, G1) |
gA += genAtom(" CD1", a3, i, D1) | gA += genAtom(" CD1", a3, i, D1) | ||
gA += genAtom(" CD2", a3, i, D2) | gA += genAtom(" CD2", a3, i, D2) | ||
Line 227: | Line 227: | ||
case 'W' : | case 'W' : | ||
gA += genAtom(" CG ", a3, i, Ghw) | gA += genAtom(" CG ", a3, i, Ghw) | ||
− | gA += genAtom(" | + | gA += genAtom(" CD2", a3, i, D1hw) |
− | gA += genAtom(" | + | gA += genAtom(" CD1", a3, i, D2hw) |
gA += genAtom(" CE2", a3, i, E2hw) | gA += genAtom(" CE2", a3, i, E2hw) | ||
gA += genAtom(" NE1", a3, i, E1hw) | gA += genAtom(" NE1", a3, i, E1hw) | ||
Line 374: | Line 374: | ||
script inline @{gA} | script inline @{gA} | ||
connect | connect | ||
+ | var xx = {element="Xx"} | ||
+ | for (var i = 1; i <= xx.size; i++) { | ||
+ | connect 1.8 {atomindex=@{xx[i].atomIndex}} | ||
+ | } | ||
# Clean up | # Clean up | ||
Line 393: | Line 397: | ||
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":" | + | "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"} | ||
Revision as of 17:45, 8 February 2014
RIBOZOME - an Alpha Helix Generator script
Creates an alpha helix from a user input string.
Jmol script by Ron Mignery with help from Angel Herráez
Contents
Script for Jmol application
The top-level function plicoGenPP asks user for a peptide sequence in a pop-up.
The top-level function plicoGenAlpha accepts a sequence string as a parameter.
Ribozome is a member of the Plico suite of protein folding tools described in User:Remig/plico. 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;plicogenpp
saved as ribozome.macro in your .jmol/macros folder as described in Macro.
# RIBOZOME - Jmol script by Ron Mignery # v1.8 beta 2/7/2014 - fix ASN mislabel # # 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. # 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) # then the chain is so labeled and separated from existing chains # if different from the first chain. # # 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 var kPHI = -57 # Dihedral angle of N-CA bond (nominally -57) var kPSI = -47 # Dihedral angle of CA-C bond (nominally -47) var kOMEGA = 180 # Dihedral angle of the peptide bond (nominally 180) var kPEPTIDE_ANGLE = 110 # C-N-CA angle (nominally 110) var kPRO_BUMP = -10 # Psi angle change increment to N-3psi when N is Pro var gCHAIN = 'A' # The chain id var gA = "" # Lookup 3 letter code from 1 letter code var 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 genAtom(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 genAA(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 = genAtom(" N ", a3, i, N0) gA += genAtom(" CA ", a3, i, CA) gA += genAtom(" C ", a3, i, C) gA += genAtom(" O ", a3, i, O) if ((aa != 'G') && (aa != 'X')) { gA += genAtom(" CB ", a3, i, CB) } # Now add AA specific atom records switch (aa) { case 'A' : break; case 'B' : gA += genAtom(" CG ", a3, i, G1) gA += genAtom(" CD ", a3, i, D1) gA += genAtom(" OE1", a3, i, E2eq) # GLN with Es switched gA += genAtom(" NE2", a3, i, E1eq) break; case 'C' : gA += genAtom(" SG ", a3, i, G2) break; case 'D' : gA += genAtom(" CG ", a3, i, G1) gA += genAtom(" OD1", a3, i, D1dn) gA += genAtom(" OD2", a3, i, D2dn) break; case 'E' : gA += genAtom(" CG ", a3, i, G1) gA += genAtom(" CD ", a3, i, D1) gA += genAtom(" OE1", a3, i, E1eq) gA += genAtom(" OE2", a3, i, E2eq) break; case 'F' : gA += genAtom(" CG ", a3, i, Gfy) gA += genAtom(" CD1", a3, i, D1fy) gA += genAtom(" CD2", a3, i, D2fy) gA += genAtom(" CE1", a3, i, E1fy) gA += genAtom(" CE2", a3, i, E2fy) gA += genAtom(" CZ ", a3, i, Zfy) break; case 'G' : break; case 'H' : gA += genAtom(" CG ", a3, i, Ghw) gA += genAtom(" ND1", a3, i, D1hw) gA += genAtom(" CD2", a3, i, D2hw) gA += genAtom(" CE1", a3, i, E2hw) gA += genAtom(" NE2", a3, i, E1hw) break; case 'I' : gA += genAtom(" CG1", a3, i, G1) gA += genAtom(" CG2", a3, i, G2) gA += genAtom(" CD1", a3, i, D1) break; case 'K' : gA += genAtom(" CG ", a3, i, G1) gA += genAtom(" CD ", a3, i, D1) gA += genAtom(" CE ", a3, i, E1) gA += genAtom(" NZ ", a3, i, Z) break; case 'L' : gA += genAtom(" CG ", a3, i, G1) gA += genAtom(" CD1", a3, i, D1) gA += genAtom(" CD2", a3, i, D2) break; case 'M' : gA += genAtom(" CG ", a3, i, G1) gA += genAtom(" SD ", a3, i, D1) gA += genAtom(" CE ", a3, i, E1) break; case 'N' : gA += genAtom(" CG ", a3, i, G1) gA += genAtom(" OD1", a3, i, D1dn) gA += genAtom(" ND2", a3, i, D2dn) break; case 'P' : gA += genAtom(" CG ", a3, i, GP) gA += genAtom(" CD ", a3, i, DP) break; case 'Q' : gA += genAtom(" CG ", a3, i, G1) gA += genAtom(" CD ", a3, i, D1) gA += genAtom(" OE1", a3, i, E1eq) gA += genAtom(" NE2", a3, i, E2eq) break; case 'R' : gA += genAtom(" CG ", a3, i, G1) gA += genAtom(" CD ", a3, i, D1) gA += genAtom(" NE ", a3, i, E1) gA += genAtom(" CZ ", a3, i, Z) gA += genAtom(" NH1", a3, i, H1) gA += genAtom(" NH2", a3, i, H2) break; case 'S' : gA += genAtom(" OG ", a3, i, G1) break; case 'T' : gA += genAtom(" OG1", a3, i, G1) gA += genAtom(" CG2", a3, i, G2) break; case 'U' : gA += genAtom("SeG ", a3, i, G1) break; case 'V' : gA += genAtom(" CG1", a3, i, G1) gA += genAtom(" CG2", a3, i, G2) break; case 'W' : gA += genAtom(" CG ", a3, i, Ghw) gA += genAtom(" CD2", a3, i, D1hw) gA += genAtom(" CD1", a3, i, D2hw) gA += genAtom(" CE2", a3, i, E2hw) gA += genAtom(" NE1", a3, i, E1hw) gA += genAtom(" CE3", a3, i, E3hw) gA += genAtom(" CZ2", a3, i, Z2hw) gA += genAtom(" CZ3", a3, i, Z3hw) gA += genAtom(" CH2", a3, i, H2hw) break; case 'X' : gA += genAtom(" Xx ", a3, i, CB) break; case 'Y' : gA += genAtom(" CG ", a3, i, Gfy) gA += genAtom(" CD1", a3, i, D1fy) gA += genAtom(" CD2", a3, i, D2fy) gA += genAtom(" CE1", a3, i, E1fy) gA += genAtom(" CE2", a3, i, E2fy) gA += genAtom(" CZ ", a3, i, Zfy) gA += genAtom(" OH ", a3, i, Hfy) break; case 'Z' : gA += genAtom(" CG ", a3, i, G1) gA += genAtom(" OD1", a3, i, D2dn) # ASN with Ds switched gA += genAtom(" ND2", a3, i, D1dn) break; default : break; } return gA }; # Generate an alpha helix function plicoGenAlpha(gSeq) { set pdbAddHydrogens FALSE if (gPlicoRecord != "") { var g = format("show file \"%s\"", gPlicoRecord) var ls = script(g) if (ls.find("FileNotFoundException")) { ls = "" } ls += format("plicoGenAlpha(\"%s\");", gSeq) write var ls @gPlicoRecord } gSeq = gSeq%9999%0 if (gSeq[2] == ':') { gCHAIN = gSeq[1] gSeq[1] = ' ' gSeq[2] = ' ' gSeq = gSeq%0 } gN = all.count + 1 # global new N atom index # Find last linkable N if any gResno = 0 # global pre-existing AA count var pn = 1 # previous gN for (var i = all.count-1; i > 0; i--) { # If found if (distance({atomno=i}, {0,0,0}) < 0.1) { pn = i # 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 } else { gResno = {atomno=i}.resno } break; } } # For each aa set appendnew false var nn = gN # new N for (var i = 1; i <= gSeq.count; i++) { # 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 += genAA(i + gResno, gSeq[i]); # gN is updated in subroutine gA += "end \"append aa\"" script inline @{gA} # 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 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 pn = nn nn = gN # Make the peptide bond connect } # 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 var xx = {element="Xx"} for (var i = 1; i <= xx.size; i++) { connect 1.8 {atomindex=@{xx[i].atomIndex}} } # Clean up connect ([UNK].CA) ([UNK].Xx and within(group, _1)) select all fix none print format("%d atoms generated", gN-1) } 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":"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"} # Get the sequence from the user gSeq = prompt("Enter AA sequence (1 char coded)", "")%9999%0 if ((gSeq != "NULL") and (gSeq.count > 0)) { print format ("Sequence=%s phi=%d psi=%d", gSeq, kPHI, kPSI) plicoGenAlpha(gSeq) } } # end of ribozome.spt
In a webpage
The main script is a small adaptation of the one above. 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.
See it in action