User talk:Remig
Jump to navigation
Jump to search
I use Jmol to study protein folding. Here is a script I wrote that accepts an amino acid sequence (1 letter encoding: "AAC...FYW" for example) and generates an alpha helix using the Model Kit:
# RIBOZOME - Jmol script by Ron Mignery with help from Dr. Angel Herráez
# v1.2 beta 10/23/2013
#
# 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
# 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
gPHI = -57 # Dihedral angle of N-CA bond (nominally -57)
gPSI = -47 # Dihedral angle of CA-C bond (nominally -47)
gOMEGA = 180 # Dihedral angle of the peptide bond (nominally 180)
gPEPTIDE_ANGLE = 110 # C-N-CA angle (nominally 110)
gPRO_BUMP = -10 # Psi angle change increment to N-3psi when N is Pro
gCHAIN = 'A' # The chain id
# 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":"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"}
# 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.470, 4.717, -0.885 ]
var H1 = [ 4.450, 6.006, -0.220 ]
var H2 = [5.833, 4.228, -0.984 ]
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
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, E2) # GLN with Es switched
gA += genAtom(" NE2", a3, i, E1)
break;
case 'C' :
gA += genAtom(" SG ", a3, i, G2)
break;
case 'D' :
gA += genAtom(" CG ", a3, i, G1)
gA += genAtom(" OD1", a3, i, D1)
gA += genAtom(" OD2", a3, i, D2)
break;
case 'E' :
gA += genAtom(" CG ", a3, i, G1)
gA += genAtom(" CD ", a3, i, D1)
gA += genAtom(" OE1", a3, i, E1)
gA += genAtom(" OE2", a3, i, E2)
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(" CG1", 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, D1)
gA += genAtom(" ND2", a3, i, D2)
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, E1)
gA += genAtom(" NE2", a3, i, E2)
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(" CD1", a3, i, D1hw)
gA += genAtom(" CD2", 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, D2) # ASN with Ds switched
gA += genAtom(" ND2", a3, i, D1)
break;
default :
break;
}
return gA
};
# Generate an alpha helix
function genAlpha(gSeq) {
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
# 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
var pb = 0
if ((gSeq.count - i) >= 2) {
if (gSeq[i + 2] == 'P') {
pb = gPRO_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} @{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
pn = nn
nn = gN
# Make the peptide bond
connect
}
# 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
# Get the sequence from the user
gSeq = prompt("Enter AA sequence (1 char coded)", "")%9999%0
if (gSeq.count > 0) {
if (gSeq[2] == ':') {
gCHAIN = gSeq[1]
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)
}
