User talk:Remig
Jump to navigation
Jump to search
I use Jmol to study protein folding. The script here before that generates polypeptide helices is now available in the Recycling Corner: Recycling_Corner/Alpha_Helix_Generator so I have removed it from here. In its place is a similar script I wrote that accepts a nucleotide sequence (1 letter encoding: "TACGAAC...GCT" for example) and generates a DNA or RNA single or double helix using the Model Kit:
# POLYMERAZE - Jmol script by Ron Mignery
# v1.0 beta 10/31/2013
#
# POLYMERAZE takes a string message encoding a nucleotide (nt) sequence
# and generates a corresponding double helix one nt at a time from the
# 5' terminus to the 3' 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 prepended with '3' then the string is considered as 3' to 5'
# If prepended with 'R' then RNA is generated instead of DNA
# If prepended with 'S' then a single strand helix is produced
#
# The IUPAC/IUBMB 1 letter code is used:
# A=Adenine C=Cytosine G=Guanine T=Thymine U=Uracil
# The following constant values determine the pitch of the helices
gO4C4C5O5 = -82.3 - 10 # Values from 30-31:B of 3BSE
gC4C5O5P = 172.5 + 0 # with the indicated tweaks to make
gC5O5PO3 = -44.0 + 9 # it work (sort of - yes, I know
gO5PO3C3 = -109.8 + 0 # the B chain P is a bit screwy)
gPO3C3C4 = -172.9 - 1 #
gCHAIN1 = 'A' # The chain id
gCHAIN2 = 'B' # The complementary chain id
# Lookup 3 letter code from 1 letter code
g3from1 = {"A":" DA", "C":" DC","G":" DG", "T":" DT", "U":" DU"}
gComp = {"A":"T", "C":"G","G":"C", "T":"A", "U":"A"}
# Generate PDB atom record
# Writes gNa or gNb
function genAtom(e, aa, i, xyz, comp) {
# Fixed column format:
#ATOM 500 O4' DA B 29 -3.745 7.211 45.474
var a = format("ATOM %5d %4s %3s ", (comp ? gNb : gNa), e, aa )
a += format("%s%4d %8.3f", (comp ? gCHAIN2 : gCHAIN1), i, xyz[1] )
a += format("%8.3f%8.3f\n", xyz[2], xyz[3] )
if (comp) gNb++; else gNa++
return a
};
# Generate a PDB nucleotide record set
# Calls genAtom that writes gNa or gNb
function genNT(i, nt, rna, comp) {
# From constructed nucleotides
var P0 = (comp ? [16.753, 4.551, 8.935] : [0.000, 0.000, 0.000])
var OP1= (comp ? [17.916, 5.432, 9.209] : [-0.973, 0.363, -1.067])
var OP2= (comp ? [16.030, 4.050, 10.108] : [0.297, -1.428, 0.272])
var O5p= (comp ? [16.110, 5.278, 7.954] : [1.351, 0.795, -0.286])
var C5p= (comp ? [16.536, 5.450, 6.610] : [1.345, 2.211, -0.125])
var C4p= (comp ? [15.487, 6.211, 5.838] : [2.739, 2.779, -0.195])
var O4p= (comp ? [14.427, 5.302, 5.454] : [3.469, 2.629, 1.048])
var C3p= (comp ? [14.830, 7.309, 6.676] : [3.624, 2.197, -1.288])
var O3p= (comp ? [14.657, 8.490, 5.878] : [4.212, 3.282, -1.983])
var C2p= (comp ? [13.509, 6.688, 7.119] : [4.692, 1.432, -0.523])
var O2p= (comp ? [12.530, 7.743, 7.153] : [5.994, 1.461, -1.226])
var C1p= (comp ? [13.175, 5.739, 5.976] : [4.797, 2.255, 0.746])
var N1ct=(comp ? [12.404, 4.556, 6.401] : [5.353, 1.551, 1.907])
var C2ct=(comp ? [11.423, 4.111, 5.551] : [6.351, 2.183, 2.606])
var O2ct=(comp ? [11.137, 4.685, 4.516] : [6.789, 3.274, 2.289])
var N3ct=(comp ? [10.784, 2.966, 5.956] : [6.823, 1.489, 3.692])
var C4ct=(comp ? [11.025, 2.241, 7.106] : [6.404, 0.251, 4.135])
var NO4ct=(comp ?[10.382, 1.213, 7.329] : [6.908, -0.242, 5.146])
var C5ct=(comp ? [12.061, 2.780, 7.971] : [5.360, -0.370, 3.337])
var C6ct=(comp ? [12.692, 3.894, 7.575] : [4.892, 0.307, 2.279])
var nC7ct=(comp ?[12.421, 2.090, 9.243] : [4.862, -1.727, 3.721])
var N9ag=(comp ? [12.426, 4.545, 6.367] : [5.333, 1.584, 1.923])
var C8ag=(comp ? [12.706, 3.662, 7.382] : [4.870, 0.450, 2.545])
var N7ag=(comp ? [11.856, 2.668, 7.467] : [5.595, 0.076, 3.570])
var C5ag=(comp ? [10.952, 2.911, 6.441] : [6.608, 1.025, 3.629])
var C6ag=(comp ? [9.818, 2.215, 5.997] : [7.694, 1.194, 4.500])
var NO6ag=(comp ?[9.379, 1.083, 6.552] : [7.955, 0.379, 5.525])
var N1ag=(comp ? [9.137, 2.727, 4.945] : [8.517, 2.246, 4.283])
var C2ag=(comp ? [9.573, 3.862, 4.390] : [8.259, 3.061, 3.256])
var nN2ag=(comp ?[8.847, 4.313, 3.345] : [9.119, 4.090, 3.100])
var N3ag=(comp ? [10.630, 4.605, 4.715] : [7.265, 3.009, 2.370])
var C4ag=(comp ? [11.285, 4.068, 5.759] : [6.465, 1.956, 2.616])
# Build PDB atom records common to all NTs
n3 = g3from1[nt]
if (n3 = "") {
n3 = " D?"
}
print format("+ %s%d/%d", n3, i, gSeq.count + gResno)
var a = genAtom(" P ", n3, i, P0, comp)
a += genAtom(" OP1", n3, i, OP1, comp)
a += genAtom(" OP2", n3, i, OP2, comp)
a += genAtom(" O5'", n3, i, O5p, comp)
a += genAtom(" C5'", n3, i, C5p, comp)
a += genAtom(" C4'", n3, i, C4p, comp)
a += genAtom(" O4'", n3, i, O4p, comp)
a += genAtom(" C3'", n3, i, C3p, comp)
a += genAtom(" O3'", n3, i, O3p, comp)
a += genAtom(" C2'", n3, i, C2p, comp)
a += genAtom(" C1'", n3, i, C1p, comp)
if (rna) {
a += genAtom(" O2'", n3, i, O2p, comp)
}
# Now add NT specific atom records
switch (nt) {
case 'A' :
a += genAtom(" N9 ", n3, i, N9ag, comp)
a += genAtom(" C8 ", n3, i, C8ag, comp)
a += genAtom(" N7 ", n3, i, N7ag, comp)
a += genAtom(" C5 ", n3, i, C5ag, comp)
a += genAtom(" C6 ", n3, i, C6ag, comp)
a += genAtom(" N6 ", n3, i, NO6ag, comp)
a += genAtom(" N1 ", n3, i, N1ag, comp)
a += genAtom(" C2 ", n3, i, C2ag, comp)
a += genAtom(" N3 ", n3, i, N3ag, comp)
a += genAtom(" C4 ", n3, i, C4ag, comp)
break;
case 'C' :
a += genAtom(" N1 ", n3, i, N1ct, comp)
a += genAtom(" C2 ", n3, i, C2ct, comp)
a += genAtom(" O2 ", n3, i, O2ct, comp)
a += genAtom(" N3 ", n3, i, N3ct, comp)
a += genAtom(" C4 ", n3, i, C4ct, comp)
a += genAtom(" N4 ", n3, i, NO4ct, comp)
a += genAtom(" C5 ", n3, i, C5ct, comp)
a += genAtom(" C6 ", n3, i, C6ct, comp)
break;
case 'G' :
a += genAtom(" N9 ", n3, i, N9ag, comp)
a += genAtom(" C8 ", n3, i, C8ag, comp)
a += genAtom(" N7 ", n3, i, N7ag, comp)
a += genAtom(" C5 ", n3, i, C5ag, comp)
a += genAtom(" C6 ", n3, i, C6ag, comp)
a += genAtom(" O6 ", n3, i, NO6ag, comp)
a += genAtom(" N1 ", n3, i, N1ag, comp)
a += genAtom(" C2 ", n3, i, C2ag, comp)
a += genAtom(" N2 ", n3, i, nN2ag, comp)
a += genAtom(" N3 ", n3, i, N3ag, comp)
a += genAtom(" C4 ", n3, i, C4ag, comp)
break;
case 'T' :
a += genAtom(" N1 ", n3, i, N1ct, comp)
a += genAtom(" C2 ", n3, i, C2ct, comp)
a += genAtom(" O2 ", n3, i, O2ct, comp)
a += genAtom(" N3 ", n3, i, N3ct, comp)
a += genAtom(" C4 ", n3, i, C4ct, comp)
a += genAtom(" O4 ", n3, i, NO4ct, comp)
a += genAtom(" C5 ", n3, i, C5ct, comp)
a += genAtom(" C6 ", n3, i, C6ct, comp)
a += genAtom(" C7 ", n3, i, nC7ct, comp)
break;
case 'U' :
a += genAtom(" N1 ", n3, i, N1ct, comp)
a += genAtom(" C2 ", n3, i, C2ct, comp)
a += genAtom(" O2 ", n3, i, O2ct, comp)
a += genAtom(" N3 ", n3, i, N3ct, comp)
a += genAtom(" C4 ", n3, i, C4ct, comp)
a += genAtom(" O4 ", n3, i, NO4ct, comp)
a += genAtom(" C5 ", n3, i, C5ct, comp)
a += genAtom(" C6 ", n3, i, C6ct, comp)
break;
default :
break;
}
return a
};
# Rotate a1 on a2 in the plane of a1, a2 and a3 to the given angle
# a1 and all connected except by a2 must be selected
function setAngle (a1, a2, a3, toangle) {
var v1={atomno = a1}.xyz - {atomno = a2}.xyz
var v2={atomno = a3}.xyz - {atomno = a2}.xyz
var axis = cross(v1, v2) + {atomno = a2}.xyz
var curangle = angle({atomno=a1}, {atomno=a2}, {atomno=a3})
rotateselected @axis {atomno = a2} @{curangle-toangle}
}
# Set the dihedral to the given angle
# a1 (or a4) and all connected except by a2 (or a3) must be selected
# If selected < unselected ==> a2 < a3 and vice versa
function setDihedral (a1, a2, a3, a4, toangle) {
var curangle = angle({atomno=a1}, {atomno=a2}, {atomno=a3}, {atomno=a4})
rotateselected {atomno=a2} {atomno=a3} @{toangle-curangle}
}
function countAtoms(seq, rna) {
var ntc = {"A":21, "C":20, "G":22, "T":20, "U":19}
var cnt = 0
for (var i = 1; i <= seq.count; i++) {
cnt += (ntc[seq[i]] + (rna ? 1 : 0))
}
return cnt
}
# Generate a helix
function genHelixStrand(gSeq, reverse, rna, double) {
var cha = ":A"
var chb = ":B"
var seq = ""
if (reverse) {
for (var i = gSeq.count; i > 0; i--) {
seq += gSeq[i]
}
}
else {
seq = gSeq
}
gNa = all.count + 1 # global new N atom index
gNb = (double ? (gNa + countAtoms(seq) + ((gResno > 0) ? 0 : 1)) : 0)
print format("gNa=%d gNb = %d", gNa, gNb)
# Find last linkable P if any
gResno = 0 # global pre-existing NT count
var pna = 1 # previous gN
for (var i = all.count-1; i > 0; i--) {
# If found
if (distance({atomno=i}, {0,0,0}) < 0.1) {
if ({atomno=i}.chain == chain[2]) {
pna = i
}
gResno = {atomno=i}.resno
break;
}
}
# For each nt
set appendnew false
var nna = gNa # new P
var nnb = gNb # new P
for (var i = 1; i <= seq.count; i++) {
if (seq[i] == "") {
continue
}
gA = "data \"append nt\"\n" # global PDB atom record
# Move polynucleotide O3p to bond distance from new nt P
var pO3 = {-0.521, 0.638, 1.234}
select all
if ((i + gResno) > 1) {
var nO3 = {atomno=@{pna+8}}.xyz
var xyz = @{pO3 - nO3}
translateselected @xyz
}
# Else 1st 5' nt so add OP3
else {
var O3n = {pO3n}.xyz
gA += genAtom(" OP3", g3from1[seq[i]], i + gResno, O3n, FALSE)
nna++
}
# Gen NT ==================================================
gA += genNT(i + gResno, seq[i], rna, FALSE); # gNa updated
if (double) {
gA += genNT(i + gResno + seq.count, gComp[seq[i]], rna, TRUE); # gNb updated
}
gA += "end \"append nt\""
script inline @{gA} # <== new atoms added here
# First shape up the comp side
if (double) {
select (@chb and (atomno < @{nnb + 6}) && (atomno >= nnb))
setDihedral(nnb+6, nnb+5, nnb+4, nnb+3, gO4C4C5O5)
select selected and (atomno != @{nnb + 5})
setDihedral(nnb+5, nnb+4, nnb+3, nnb, gC4C5O5P)
}
# Adjust link dihedrals of the new
select (@cha and (atomno > @{nna+4}) or (@chb and (atomno >= nnb)))
setDihedral(nna+3, nna+4, nna+5, nna+6, gO4C4C5O5)
setDihedral(nna, nna+3, nna+4, nna+5, gC4C5O5P)
# If any older
if (i > 1) {
# Now move the old
select (@cha and (atomno < nna) or (@chb and (atomno < nnb)))
setAngle(nna, pna+8, pna+7, 120.0)
select (@cha and (atomno < @{nna+3}) or (@chb and (atomno < nnb)))
setDihedral(nna+4, nna+3, nna, pna+8, gC5O5PO3)
select (@cha and (atomno < nna) or (@chb and (atomno < nnb)))
setDihedral(nna+3, nna, pna+8, pna+7, gO5PO3C3)
setDihedral(nna, pna+8, pna+7, pna+5, gPO3C3C4)
}
# Step new and previous N
pna = nna; pnb = nnb
nna = gNa; nnb = gNb
}
# Make the nucleotide bonds
connect
# Clean up
select all
print format("%d atoms generated for chain %s", gNa+gNb,
(comp ? gCHAIN2 : gCHAIN1))
}
# Generate a helix or two
function genHelix(gSeq) {
var single = FALSE
var reverse = FALSE
var rna = FALSE
var done = FALSE
if (gSeq[2] == ':') {
gCHAIN1 = gSeq[1]
gSeq[1] = ' '; gSeq[2] = ' '
gSeq = gSeq%0
}
else if (gSeq[3] == ':') {
gCHAIN1 = gSeq[1]
gCHAIN2 = gSeq[2]
gSeq[1] = ' '; gSeq[2] = ' '; gSeq[3] = ' '
gSeq = gSeq%0
}
while (done == FALSE) {
done = TRUE;
if (gSeq[1] == 'S') {
single = TRUE;
done = FALSE;
}
else if (gSeq[1] == '3') {
reverse = TRUE;
done = FALSE;
}
else if (gSeq[1] == 'R') {
rna = TRUE;
done = FALSE;
}
if (done == FALSE) {
gSeq[1] = ' '
gSeq = gSeq%0
}
}
print format ("Sequence=%s single=%s reverse=%s", gSeq, single, reverse)
print format ("rna=%s", rna)
# Gen first strand
genHelixStrand(gSeq, reverse, rna, single ? FALSE : TRUE)
}
# ==============================================
echo Generating Alpha Helix
# Get the sequence from the user
gSeq = prompt("Enter NT sequence (ACGTU)", "")%9999%0
if (gSeq.count > 0) {
genHelix(gSeq)
}
I thought adapting the Ribozome script to polynucleotides would be easy since there are only four types rather than 20. I was not easy and this script has its problems. With five rotors in a row for each nucleotide, finding correct values for the bond angles was beyond me. If anyone can tweak them better, please do and let me know. Regardless I plan to develop some more utility scripts to make that process easier and hope to update this script in the hopefully not too distant future.
