Difference between revisions of "User talk:Remig"

From Jmol
Jump to navigation Jump to search
(Alpha helix generater)
 
(Ribozome - an Alpha Helix Generator)
Line 1: Line 1:
 
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:
 
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:
−
<pre>function get3from1(c) {
+
<pre># RIBOZOME - Jmol script by Ron Mignery
 +
# v1.0 beta 10/19/2013
 +
#
 +
# 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
 +
#
 +
# The IUPAC/IUBMB 1 letter code is used:
 +
# A=ALAnine B=GLUtamine?* 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=ASparagiNe?**
 +
#   *Either GLU or GLN: drawn as GLN with chi3 flipped
 +
#   **Either ASP or ASN: drawn as ASN with chi3 flipped
 +
# ***Not supported
 +
# ****Unknown aa will be drawn as ALA
 +
#
 +
 
 +
# The following constant values determine the pitch of the alpha helix
 +
var PHI = -50 # Dihedral angle of N-CA bond (nominally -50)
 +
var PSI = -60 # Dihedral angle of CA-C bond (nominally -60)
 +
var OMEGA = 180 # Dihedral angle of the peptide bond (nominally 180)
 +
var PEPTIDE_ANGLE = 110 # C-N-CA angle (nominally 110)
 +
var PRO_BUMP = -15 # Psi angle change to N-3psi when N is Pro
 +
 
 +
# Lookup 3 letter code from 1 letter code
 +
function get3from1(c) {
 
ret = ""
 
ret = ""
 
switch (c) {
 
switch (c) {
 
case 'A':
 
case 'A':
 +
case 'X':
 
ret = "ALA";
 
ret = "ALA";
 
break;
 
break;
Line 37: Line 68:
 
break;
 
break;
 
case 'N':
 
case 'N':
 +
case 'Z':
 
ret = "ASN";
 
ret = "ASN";
 
break;
 
break;
Line 42: Line 74:
 
ret = "PRO";
 
ret = "PRO";
 
break;
 
break;
 +
case 'B':
 
case 'Q':
 
case 'Q':
 
ret = "GLN";
 
ret = "GLN";
Line 70: Line 103:
 
};
 
};
  
 +
# Generate PDM atom record
 
function genAtom(n, e, aa, i, xyz) {
 
function genAtom(n, e, aa, i, xyz) {
 
a =  format("ATOM  %5d  %3s %3s A", n, e, aa )           
 
a =  format("ATOM  %5d  %3s %3s A", n, e, aa )           
Line 77: Line 111:
 
};
 
};
  
 +
# Generate an amino acid record set
 
function genAA(i, aa, x) {
 
function genAA(i, aa, x) {
 
n = x
 
n = x
  
 +
# From constructed AAs
 
N0 = [0.0, 0.0, 0.0]
 
N0 = [0.0, 0.0, 0.0]
 
CA = [ 0.200, 1.174, 0.911 ]
 
CA = [ 0.200, 1.174, 0.911 ]
−
C = [ -1.129, 1.783, 1.241 ]
+
C = [ -1.110, 1.668, 1.425 ] #[ -1.129, 1.783, 1.241 ]
−
O = [ -1.241, 1.967, 2.726 ]
+
O = [ -1.320, 1.693, 2.62 ] #[ -1.241, 1.967, 2.726 ]
 
CB = [ 1.062, 2.1950, 0.230 ]
 
CB = [ 1.062, 2.1950, 0.230 ]
−
G1 = [ 2.396, 1.588, -0.091 ]
+
−
G2 = [ 0.680, 3.652, 0.423]
+
G1 = [ 2.359, 1.607, -0.344] #2.396, 1.588, -0.091 ]
−
Gfy = [ 2.368, 1.471, -0.0152 ]
+
G2 = [ 1.363, 3.336, 1.157 ] #0.680, 3.652, 0.423]
−
Ghw = [ 2.406, 1.626, -0.134 ]
+
D1 = [ 3.222, 2.656, -1.048 ] #[ 3.225, 2.340, -1.096]
−
D1 = [ 3.225, 2.340, -1.096]
+
D2 = [ 3.143, 0.904, 0.725 ] #[ 3.189, 1.093, 1.087]
−
D1hw = [3.498, 1.936, 0.675]
+
E1 = [ 3.645, 3.749, -0.167 ] #[ 3.652, 3.503, -0.111 ]
 +
E2 = [ 2.491, 3.234, -2.249 ] #[ 4.342, 1.591, -1.456 ]
 +
Z  = [ 4.470, 4.717, -0.885 ] #[ 4.115, 3.339, 1.403 ]
 +
H1 = [ 4.450, 6.006, -0.220 ] #[4.087, 4.572, 2.139]
 +
H2 = [5.833, 4.228, -0.984 ] #[5.469, 2.866, 1.296]
 +
 +
Gp = [ 2.008, 1.24, -0.46 ]
 +
Dp = [1.022, 0.213, -1.031 ]
 +
 +
Gfy  = [ 2.368, 1.471, -0.0152 ]
 
D1fy = [ 3.346, 1.524, 0.921 ]
 
D1fy = [ 3.346, 1.524, 0.921 ]
−
D2 = [ 3.189, 1.093, 1.087]
 
−
D2hw = [ 2.713, 0.901, -1.211 ]
 
 
D2fy = [ 2.493, 0.516, -1.151 ]
 
D2fy = [ 2.493, 0.516, -1.151 ]
 
E1fy = [ 4.513, 0.615, 0.8244 ]
 
E1fy = [ 4.513, 0.615, 0.8244 ]
 +
E2fy = [ 3.528, -0.336, -1.206 ]
 +
Zfy  = [ 4.588, -0.285, -0.168 ]
 +
Hfy = [ 5.738, -1.245, -0.233 ]
 +
 +
Ghw  = [ 2.406, 1.626, -0.134 ]
 +
D1hw = [3.498, 1.936, 0.675]
 +
D2hw = [ 2.713, 0.901, -1.211 ]
 
E1hw = [ 4.160, 0.518, -1.178 ]
 
E1hw = [ 4.160, 0.518, -1.178 ]
−
E2fy = [ 3.528, -0.336, -1.206 ]
 
 
E2hw = [ 4.622, 1.160, 0.0816 ]
 
E2hw = [ 4.622, 1.160, 0.0816 ]
 
E3hw = [ 3.789, 2.523, 1.944 ]
 
E3hw = [ 3.789, 2.523, 1.944 ]
Line 103: Line 152:
 
Z3hw = [ 5.014, 2.550, 2.503 ]
 
Z3hw = [ 5.014, 2.550, 2.503 ]
 
H2hw = [ 6.153, 1.846, 1.844 ]
 
H2hw = [ 6.153, 1.846, 1.844 ]
−
Zfy = [ 4.588, -0.285, -0.168 ]
 
−
Hfy = [ 5.738, -1.245, -0.233 ]
 
−
N1 = [ -1.965, 2.307, 0.206 ]
 
−
Gp = [ 2.008, 1.24, -0.46 ]
 
−
Dp = [1.022, 0.213, -1.031 ]
 
−
E1 = [ 3.652, 3.503, -0.111 ]
 
−
E2 = [ 4.342, 1.591, -1.456 ]
 
−
Z = [ 4.115, 3.339, 1.403 ]
 
−
H1 = [4.087, 4.572, 2.139]
 
−
H2 = [5.469, 2.866, 1.296]
 
 
 
 +
N1 = [ 2.069, -2.122, -0.554] #[ -1.965, 2.307, 0.206 ]
 +
 +
# Build PDB atom records common to all AAs
 
a3 = get3from1(seq[i])
 
a3 = get3from1(seq[i])
 +
if (a3 == "") {
 +
a3 = "UNK"
 +
}
 +
print format("+ %s%d", a3, i)
 
a = genAtom(n++, "N  ", a3, i, N0)
 
a = genAtom(n++, "N  ", a3, i, N0)
 
a += genAtom(n++, "CA ", a3, i, CA)
 
a += genAtom(n++, "CA ", a3, i, CA)
Line 123: Line 169:
 
}
 
}
  
 +
# Now add AA specific atom records
 
switch (aa) {
 
switch (aa) {
 
case 'A' :
 
case 'A' :
 +
case 'X' :
 +
break;
 +
case 'B' :
 +
a += genAtom(n++, "CG ", a3, i, G1)
 +
a += genAtom(n++, "CD ", a3, i, D1)
 +
a += genAtom(n++, "OE1", a3, i, E2) # GLN with Es switched
 +
a += genAtom(n++, "NE2", a3, i, E1)
 
break;
 
break;
 
case 'C' :
 
case 'C' :
−
a += genAtom(n++, "SG ", a3, i, G1)
+
a += genAtom(n++, "SG ", a3, i, G2)
 
break;
 
break;
 
case 'D' :
 
case 'D' :
Line 133: Line 187:
 
a += genAtom(n++, "OD1", a3, i, D1)
 
a += genAtom(n++, "OD1", a3, i, D1)
 
a += genAtom(n++, "OD2", a3, i, D2)
 
a += genAtom(n++, "OD2", a3, i, D2)
−
break;
+
break;
 
case 'E' :
 
case 'E' :
 
a += genAtom(n++, "CG ", a3, i, G1)
 
a += genAtom(n++, "CG ", a3, i, G1)
Line 234: Line 288:
 
a += genAtom(n++, "CZ ", a3, i, Zfy)
 
a += genAtom(n++, "CZ ", a3, i, Zfy)
 
a += genAtom(n++, "OH ", a3, i, Hfy)
 
a += genAtom(n++, "OH ", a3, i, Hfy)
 +
break;
 +
case 'Z' :
 +
a += genAtom(n++, "CG ", a3, i, G1)
 +
a += genAtom(n++, "OD1", a3, i, D2) # ASN with Ds switched
 +
a += genAtom(n++, "ND2", a3, i, D1)
 
break;
 
break;
 
default :
 
default :
Line 240: Line 299:
  
 
return a
 
return a
−
}
+
};
 +
 
 +
function rotateNward (a1, a2, angle) {
 +
select atomno<a1
 +
fix atomno>=a2
 +
rotateselected {atomno=a1} {atomno=a2} @angle
 +
#print format("a1=%d a2=%d angle=%d", a1, a2, angle) #DEBUG
 +
 
 +
};
  
 
# GenAlph
 
# GenAlph
−
function genAlpha(seq, pa) {
+
function genAlpha(seq, PHI, PSI, OMEGA, PEPTIDE_ANGLE, PRO_BUMP) {
  
 
# For each aa
 
# For each aa
Line 250: Line 317:
 
pn = 0
 
pn = 0
 
pc= 0 # previous C
 
pc= 0 # previous C
−
for (var i = 1; i <= seq.count; i++) {
+
ca1 = 0; ca2 = 0; ca3 = 0
 +
for (var i = 1; i <= seq.count; i++) {
  
 +
# Step previous N
 
pn = n
 
pn = n
−
m = pn + 1
 
  
−
# Anticipate PRO
+
# Move polypeptide C to bond distance from new AA N
−
dpsi = 0
 
−
dphi = 0
 
−
if ((seq.count - i) >= 4) {
 
−
if (seq[i + 4] == 'P') {
 
−
dpsi = =9
 
−
}
 
−
}
 
−
if ((seq.count - i) >= 3) {
 
−
if (seq[i + 3] == 'P') {
 
−
dphi = 10
 
−
}
 
−
}
 
−
if ((seq.count - i) >= 2) {
 
−
if (seq[i + 2] == 'P') {
 
−
dpsi = -11
 
−
#dphi = -0
 
−
}
 
−
}
 
−
if ((seq.count - i) >= 1) {
 
−
if (seq[i + 1] == 'P') {
 
−
dpsi = 25
 
−
#dphi = -81
 
−
}
 
−
}
 
−
 
−
# Move polypeptide C to bond distance to next AA N
 
 
select all
 
select all
 
fix none
 
fix none
−
translateselected { 1.955, -1.993, 0.000 }
+
translateselected {2.069, -2.122, -0.554 } #N1
  
 
# Gen AA
 
# Gen AA
Line 291: Line 333:
 
a += "end \"append aa\""
 
a += "end \"append aa\""
 
script inline @{a}
 
script inline @{a}
 +
 +
# If PRO
 +
# Adjust i-3psi as rigid PRO D4 bumps Oi-4 as i-3psi
 +
# is the most torqued psi and psis are the easiest to twist
 +
if (seq[i] == 'P') {
 +
rotateNward (ca3, ca3+1, PRO_BUMP)
 +
}
  
 
# If not first AA
 
# If not first AA
 
if (pc > 0) {
 
if (pc > 0) {
 +
 +
# Gen axis on previous n perpendicular to the plane
 +
# containing atoms pc, pn and pn+1
 
v1={atomno=pc}.xyz - {atomno=pn}.xyz
 
v1={atomno=pc}.xyz - {atomno=pn}.xyz
−
v2={atomno=m}.xyz - {atomno=pn}.xyz
+
v2={atomno=@{pn+1}}.xyz - {atomno=pn}.xyz
−
# Gen normalized axis perpendicular to the plane contaoning atoms pc, pn and m
+
axis = cross(v1, v2)
−
axis = cross(v1, v2) - {atomno=pn}.xyz
+
−
 
+
# Center on atom previous n
−
# Center on atom pn
 
 
axis += {atomno=pn}.xyz
 
axis += {atomno=pn}.xyz
  
−
# Rotate the polypeptide on the new AA
+
# Rotate the polypeptide N on the new AA C to the
 +
# desired angle (nominally 110)
 
select atomno<pn
 
select atomno<pn
 
fix atomno>=pn
 
fix atomno>=pn
 +
rotateselected @axis {atomno=pn} @{PEPTIDE_ANGLE - 77}
  
−
# pull peptide  angle trigonal
+
# Make omega dihedral OMEGA (nominally 180)
−
rotateselected @axis {atomno=pn} @{pa-62}
+
rotateselected {atomno=pc} {atomno=pn} @{OMEGA - 123}
  
−
# Make peptide bond dihedral 180
+
# Make the phi PHI (nominally -50)
−
rotateselected {atomno=pc} {atomno=pn} 52
+
rotateselected {atomno=pn} {atomno=@{pn+1}} @{PHI - 30}
  
−
# Make the phi -60
+
# Make the psi PSI (nominally -60)
−
rotateselected {atomno=pn} {atomno=m} @{dphi - 68}
 
−
 
 
−
# Make the psi -45
 
−
pca = pc - 1
 
 
fix atomno>pca
 
fix atomno>pca
 
select atomno<pn and (atomno != @{pc+1})
 
select atomno<pn and (atomno != @{pc+1})
−
rotateselected {atomno=pc} {atomno=pca} @{dpsi}
+
rotateselected {atomno=pc} {atomno=@{pc-1}} @{PSI + 60}
  
 
# If aromatic go trans on chi 1
 
# If aromatic go trans on chi 1
 
select atomno>@{pn+4} and atomno<n
 
select atomno>@{pn+4} and atomno<n
 
if ((seq[i] == 'H') || (seq[i] == 'W') || (seq[i] == 'F') || (seq[i] == 'Y')) {
 
if ((seq[i] == 'H') || (seq[i] == 'W') || (seq[i] == 'F') || (seq[i] == 'Y')) {
−
rotateselected {atomno=@{pn+1}} {atomno=@{pn+4}} -120
+
#rotateselected {atomno=@{pn+1}} {atomno=@{pn+4}} -60
 
}
 
}
  
 +
# Save last three CAs for proline bumps
 +
ca3 = ca2; ca2 = ca1; ca1 = pn + 1
 
}
 
}
 +
 +
# Step previous C
 
pc =pn + 2
 
pc =pn + 2
  
 +
# Make the peptide bond
 
connect
 
connect
 
}
 
}
 +
 +
# Clean up
 
select all
 
select all
 
fix none
 
fix none
 
}
 
}
−
 
  
 
echo Generating Alpha Helix
 
echo Generating Alpha Helix
−
zap
 
  
 
# Get the sequence from the user
 
# Get the sequence from the user
−
seq = prompt("Enter AA sequence (1 char coded)", "")%9999
+
seq = prompt("*** Any existing will be cleared ***\nEnter AA sequence (1 char coded)", "")%9999
−
pa = prompt("Enter the peptide angle", "120")
+
if (seq.count > 0) {
−
print format ("seq=%s at %d deg", seq, pa)
+
zap # disable to keep existing structures
−
genAlpha(seq , pa)
+
print format ("seq=%s  phi=%d  psi=%d", seq, PHI, PSI)
−
 
+
print format ("phi=%d  psi=%d  peptide angle=%d  pro bump=%d", PEPTIDE_ANGLE, PRO_BUMP)
−
 
+
genAlpha(seq, PHI,PSI, OMEGA, PEPTIDE_ANGLE, PRO_BUMP) # defined at top of scripts
−
 
+
}
 
</pre>
 
</pre>

Revision as of 03:07, 20 October 2013

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 
#	v1.0 beta	10/19/2013
#
#	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
#
#	The IUPAC/IUBMB 1 letter code is used:
#	A=ALAnine B=GLUtamine?* 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=ASparagiNe?**
#	   *Either GLU or GLN: drawn as GLN with chi3 flipped 
#	  **Either ASP or ASN: drawn as ASN with chi3 flipped
#	 ***Not supported
#	****Unknown aa will be drawn as ALA 
#

# The following constant values determine the pitch of the alpha helix
var PHI = -50	# Dihedral angle of N-CA bond (nominally -50)
var PSI = -60	# Dihedral angle of CA-C bond (nominally -60)
var OMEGA = 180	# Dihedral angle of the peptide bond (nominally 180)
var PEPTIDE_ANGLE = 110	# C-N-CA angle (nominally 110)
var PRO_BUMP = -15 # Psi angle change to N-3psi when N is Pro

# Lookup 3 letter code from 1 letter code
function get3from1(c) {
	ret = ""
	switch (c) {
	case 'A':
	case 'X':
		ret = "ALA";
		break;
	case 'C':
		ret = "CYS";
		break;
	case 'D':
		ret = "ASP";
		break;
	case 'E':
		ret = "GLU";
		break;
	case 'F':
		ret = "PHE";
		break;
	case 'G':
		ret = "GLY";
		break;
	case 'H':
		ret = "HIS";
		break;
	case 'I':
		ret = "ILE";
		break;
	case 'K':
		ret = "LYS";
		break;
	case 'L':
		ret = "LEU";
		break;
	case 'M':
		ret = "MET";
		break;
	case 'N':
	case 'Z':
		ret = "ASN";
		break;
	case 'P':
		ret = "PRO";
		break;
	case 'B':
	case 'Q':
		ret = "GLN";
		break;
	case 'R':
		ret = "ARG";
		break;
	case 'S':
		ret = "SER";
		break;
	case 'T':
		ret = "THR";
		break;
	case 'U':
		ret = "SEC";
		break;
	case 'V':
		ret = "VAL";
		break;
	case 'W':
		ret = "TRP";
		break;
	case 'Y':
		ret = "TYR";
		break;
	}
	return ret
};

# Generate PDM atom record
function genAtom(n, e, aa, i, xyz) {
		a =  format("ATOM  %5d  %3s %3s A", n, e, aa )          
		a +=  format("%4d    %8.3f", i, xyz[1] )          
		a +=  format("%8.3f%8.3f\n", xyz[2], xyz[3] )
		return a
};

# Generate an amino acid record set
function genAA(i, aa, x) {
	n = x

	# From constructed AAs
	N0 = [0.0, 0.0, 0.0]
	CA = [ 0.200, 1.174, 0.911 ]
	C  = [ -1.110, 1.668, 1.425 ] #[ -1.129, 1.783, 1.241 ]
	O  = [ -1.320, 1.693, 2.62 ] #[ -1.241, 1.967, 2.726 ]
	CB = [ 1.062, 2.1950, 0.230 ]
	
	G1 = [ 2.359, 1.607, -0.344] #2.396, 1.588, -0.091 ]
	G2 = [ 1.363, 3.336, 1.157 ] #0.680, 3.652, 0.423]
	D1 = [ 3.222, 2.656, -1.048 ] #[ 3.225, 2.340, -1.096]
	D2 = [ 3.143, 0.904, 0.725 ] #[ 3.189, 1.093, 1.087]
	E1 = [ 3.645, 3.749, -0.167 ] #[ 3.652, 3.503, -0.111 ]
	E2 = [ 2.491, 3.234, -2.249 ] #[ 4.342, 1.591, -1.456 ]
	Z  = [ 4.470, 4.717, -0.885 ] #[ 4.115, 3.339, 1.403 ]
	H1 = [ 4.450, 6.006, -0.220 ] #[4.087, 4.572, 2.139]
	H2 = [5.833, 4.228, -0.984 ] #[5.469, 2.866, 1.296]
	
	Gp = [ 2.008, 1.24, -0.46 ]
	Dp = [1.022, 0.213, -1.031 ]
	
	Gfy  = [ 2.368, 1.471, -0.0152 ]
	D1fy = [ 3.346, 1.524, 0.921 ]
	D2fy = [ 2.493, 0.516, -1.151 ]
	E1fy = [ 4.513, 0.615, 0.8244 ]
	E2fy = [ 3.528, -0.336, -1.206 ]
	Zfy  = [ 4.588, -0.285, -0.168 ]
	Hfy = [ 5.738, -1.245, -0.233 ]
	
	Ghw  = [ 2.406, 1.626, -0.134 ]
	D1hw = [3.498, 1.936, 0.675]
	D2hw = [ 2.713, 0.901, -1.211 ]
	E1hw = [ 4.160, 0.518, -1.178 ]
	E2hw = [ 4.622, 1.160, 0.0816 ]
	E3hw = [ 3.789, 2.523, 1.944 ]
	Z2hw = [ 5.973, 1.177, 0.689 ]
	Z3hw = [ 5.014, 2.550, 2.503 ]
	H2hw = [ 6.153, 1.846, 1.844 ]
	
	N1 = [ 2.069, -2.122, -0.554] #[ -1.965, 2.307, 0.206 ]
	
	# Build PDB atom records common to all AAs
	a3 = get3from1(seq[i])
	if (a3 == "") {
		a3 = "UNK"
	}
	print format("+ %s%d", a3, i)
	a = genAtom(n++, "N  ", a3, i, N0)
	a += genAtom(n++, "CA ", a3, i, CA)
	a += genAtom(n++, "C  ", a3, i, C)
	a += genAtom(n++, "O  ", a3, i, O)
	if (seq[i] != 'G') {
		a += genAtom(n++, "CB ", a3, i, CB)
	}

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

	return a
};

function rotateNward (a1, a2, angle) {
	select atomno<a1
	fix atomno>=a2
	rotateselected {atomno=a1} {atomno=a2} @angle
#print format("a1=%d a2=%d angle=%d", a1, a2, angle)	#DEBUG

};

# GenAlph
function genAlpha(seq, PHI, PSI, OMEGA, PEPTIDE_ANGLE, PRO_BUMP) {

	# For each aa
	set appendnew false
	n = 1
	pn = 0
	pc= 0 # previous C
	ca1 = 0; ca2 = 0; ca3 = 0
		for (var i = 1; i <= seq.count; i++) {

		# Step previous N
		pn = n

		# Move polypeptide C to bond distance from new AA N
		select all
		fix none
		translateselected {2.069, -2.122, -0.554 } #N1

		# Gen AA
		a = "data \"append aa\"\n"
		a += genAA(i, seq[i], n);
		a += "end \"append aa\""
		script inline @{a}

		# If PRO
		# Adjust i-3psi as rigid PRO D4 bumps Oi-4 as i-3psi
		# is the most torqued psi and psis are the easiest to twist
		if (seq[i] == 'P') {
			rotateNward (ca3, ca3+1, PRO_BUMP)
		}

		# If not first AA
		if (pc > 0) {

			# Gen axis on previous n perpendicular to the plane
			# containing atoms pc, pn and pn+1
			v1={atomno=pc}.xyz - {atomno=pn}.xyz
			v2={atomno=@{pn+1}}.xyz - {atomno=pn}.xyz
			axis = cross(v1, v2)
			
			# Center on atom previous n
			axis += {atomno=pn}.xyz

			# Rotate the polypeptide N on the new AA C to the
			# desired angle (nominally 110)
			select atomno<pn
			fix atomno>=pn
			rotateselected @axis {atomno=pn} @{PEPTIDE_ANGLE - 77}

			# Make omega dihedral OMEGA (nominally 180)
			rotateselected {atomno=pc} {atomno=pn} @{OMEGA - 123}

			# Make the phi PHI (nominally -50)
			rotateselected {atomno=pn} {atomno=@{pn+1}} @{PHI - 30}

			# Make the psi PSI (nominally -60)
			fix atomno>pca
			select atomno<pn and (atomno != @{pc+1})
			rotateselected {atomno=pc} {atomno=@{pc-1}} @{PSI + 60}

			# If aromatic go trans on chi 1
			select atomno>@{pn+4} and atomno<n
			if ((seq[i] == 'H') || (seq[i] == 'W') || (seq[i] == 'F') || (seq[i] == 'Y')) {
				#rotateselected {atomno=@{pn+1}} {atomno=@{pn+4}} -60
			}

			# Save last three CAs for proline bumps
			ca3 = ca2; ca2 = ca1; ca1 = pn + 1
		}
		
		# Step previous C
		pc =pn + 2

		# Make the peptide bond
		connect
	}
	
	# Clean up
	select all
	fix none
}

echo Generating Alpha Helix

# Get the sequence from the user
seq = prompt("*** Any existing will be cleared ***\nEnter AA sequence (1 char coded)", "")%9999
if (seq.count > 0) {
	zap # disable to keep existing structures
	print format ("seq=%s  phi=%d  psi=%d", seq, PHI, PSI)
	print format ("phi=%d  psi=%d  peptide angle=%d  pro bump=%d", PEPTIDE_ANGLE, PRO_BUMP)
	genAlpha(seq, PHI,PSI, OMEGA, PEPTIDE_ANGLE, PRO_BUMP)	# defined at top of scripts
}

Contributors

Remig

Retrieved from "https://wiki.jmol.org/index.php?title=User_talk:Remig&oldid=8127"