Difference between revisions of "Recycling Corner/DNA Generator"
(O2' now numbers before C1') |
(Avoid "axis," a newly reserved word) |
||
(One intermediate revision by the same user not shown) | |||
Line 33: | Line 33: | ||
Copy and paste the following into a text editor and save in your scripts directory as polymeraze.spt. | Copy and paste the following into a text editor and save in your scripts directory as polymeraze.spt. | ||
<pre># POLYMERAZE - Jmol script by Ron Mignery | <pre># POLYMERAZE - Jmol script by Ron Mignery | ||
− | # v1. | + | # v1.16 beta 4/12/2016 -axis is now a reserved word |
# | # | ||
# POLYMERAZE takes a string message encoding a nucleotide (nt) sequence | # POLYMERAZE takes a string message encoding a nucleotide (nt) sequence | ||
Line 236: | Line 236: | ||
# a1 and all connected except by a2 must be selected | # a1 and all connected except by a2 must be selected | ||
function set_angle_no (a1no, a2no, a3no, toangle) { | function set_angle_no (a1no, a2no, a3no, toangle) { | ||
− | |||
− | |||
var c = gChain1 | var c = gChain1 | ||
− | var a1 = {(atomno=a1no) and (chain=c) and | + | var a1 = {(atomno=a1no) and (chain=c) and thisModel} |
− | var a2 = {(atomno=a2no) and (chain=c) and | + | var a2 = {(atomno=a2no) and (chain=c) and thisModel} |
− | var a3 = {(atomno=a3no) and (chain=c) and | + | var a3 = {(atomno=a3no) and (chain=c) and thisModel} |
var v1 = (a1.xyz - a2.xyz) | var v1 = (a1.xyz - a2.xyz) | ||
var v2 = (a3.xyz - a2.xyz) | var v2 = (a3.xyz - a2.xyz) | ||
− | var | + | var caxis = cross(v1, v2) + a2.xyz |
var curangle = angle(a1, a2, a3) | var curangle = angle(a1, a2, a3) | ||
− | rotateselected @ | + | rotateselected @caxis @a2 @{curangle-toangle} |
} | } | ||
Line 253: | Line 251: | ||
# If selected < unselected ==> a2 < a3 and vice versa | # If selected < unselected ==> a2 < a3 and vice versa | ||
function set_dihedral_no (a1no, a2no, a3no, a4no, toangle) { | function set_dihedral_no (a1no, a2no, a3no, a4no, toangle) { | ||
− | |||
− | |||
var c = gChain1 | var c = gChain1 | ||
− | var a1 = {(atomno=a1no) and (chain=c) and | + | var a1 = {(atomno=a1no) and (chain=c) and thisModel} |
− | var a2 = {(atomno=a2no) and (chain=c) and | + | var a2 = {(atomno=a2no) and (chain=c) and thisModel} |
− | var a3 = {(atomno=a3no) and (chain=c) and | + | var a3 = {(atomno=a3no) and (chain=c) and thisModel} |
− | var a4 = {(atomno=a4no) and (chain=c) and | + | var a4 = {(atomno=a4no) and (chain=c) and thisModel} |
var curangle = angle(a1, a2, a3, a4) | var curangle = angle(a1, a2, a3, a4) | ||
rotateselected @a2 @a3 @{toangle-curangle} | rotateselected @a2 @a3 @{toangle-curangle} | ||
Line 282: | Line 278: | ||
if (reverse) { | if (reverse) { | ||
for (var i = gSeq.count; i > 0; i--) { | for (var i = gSeq.count; i > 0; i--) { | ||
− | seq += gSeq[i]%9999%0 | + | if ("AUCGT".find(gSeq[i]) > 0) { |
+ | seq += gSeq[i]%9999%0 | ||
+ | } | ||
+ | else { | ||
+ | print format("Unknown nucleotide %s!", seq[i]) | ||
+ | } | ||
} | } | ||
} | } | ||
else { | else { | ||
− | seq = gSeq%9999%0 | + | for (var i = 1; i <= gSeq.count; i++) { |
+ | if ("AUCGT".find(gSeq[i]) > 0) { | ||
+ | seq += gSeq[i]%9999%0 | ||
+ | } | ||
+ | else { | ||
+ | print format("Unknown nucleotide %s!", seq[i]) | ||
+ | } | ||
+ | } | ||
} | } | ||
Line 299: | Line 307: | ||
# global new P atom index for chain A | # global new P atom index for chain A | ||
− | gNa = ((m > 0) ? { | + | gNa = ((m > 0) ? {thisModel}.atomno.max + 1 : 1) |
gNb = 0 | gNb = 0 | ||
if (double) { | if (double) { | ||
Line 317: | Line 325: | ||
appendNew = false | appendNew = false | ||
− | # | + | # If there is an atom at the origin |
while (m > 0) { | while (m > 0) { | ||
− | var ai = {within(0.1, {0 0 0}) and | + | var ai = {within(0.1, {0 0 0}) and thisModel} |
− | if (ai | + | if (ai) { |
# If on the same chain | # If on the same chain | ||
Line 329: | Line 337: | ||
echo "Adding to existing strand..." | echo "Adding to existing strand..." | ||
pNa = ai.atomno | pNa = ai.atomno | ||
− | aResno = | + | aResno = get_resno_max(gChain1) + 1 |
− | + | gNa = {(chain=gChain1) and thisModel}.atomno.max + 1 | |
− | gNa = {(chain=gChain1) and | ||
− | |||
gNb += gNa | gNb += gNa | ||
Line 340: | Line 346: | ||
gA = "data \"append nt\"\n" # global PDB atom record | gA = "data \"append nt\"\n" # global PDB atom record | ||
for (var j = 1; j <= all.atomno.max; j++) { | for (var j = 1; j <= all.atomno.max; j++) { | ||
− | var aj = {(atomno=j) and | + | var aj = {(atomno=j) and thisModel |
and (chain == gChain2)} | and (chain == gChain2)} | ||
− | if (aj | + | if (aj) { |
gA += gen_atom_nt(aj.atomName, aj.group, | gA += gen_atom_nt(aj.atomName, aj.group, | ||
(aj.resno+seq.count+cSeq.count), | (aj.resno+seq.count+cSeq.count), | ||
Line 356: | Line 362: | ||
# Else move all from the new chain rightward on X axis | # Else move all from the new chain rightward on X axis | ||
else { | else { | ||
− | select { | + | select {thisModel and (chain!=gChain1) |
and (chain!=gChain2)} | and (chain!=gChain2)} | ||
translateselected {0, -20, -20 } | translateselected {0, -20, -20 } | ||
Line 373: | Line 379: | ||
# For each NT | # For each NT | ||
for (var i = 1; i <= seq.count; i++) { | for (var i = 1; i <= seq.count; i++) { | ||
− | |||
− | |||
− | |||
− | |||
# Move polynucleotide O3p to bond distance 1.59 from new nt P | # Move polynucleotide O3p to bond distance 1.59 from new nt P | ||
var pO3 = { -0.759, 0.925, 1.048} | var pO3 = { -0.759, 0.925, 1.048} | ||
if (double) { | if (double) { | ||
− | select ((@cha or @chb) and | + | select ((@cha or @chb) and thisModel) |
} | } | ||
else { | else { | ||
− | select (@cha and | + | select (@cha and thisModel) |
} | } | ||
if ((i + aResno) > 2) { | if ((i + aResno) > 2) { | ||
var nO3 = {@cha and (atomno=@{pNa+8}) | var nO3 = {@cha and (atomno=@{pNa+8}) | ||
− | and | + | and thisModel}.xyz |
var xyz = @{pO3 - nO3} | var xyz = @{pO3 - nO3} | ||
− | translateselected @xyz | + | if (xyz) { |
+ | translateselected @xyz | ||
+ | } | ||
} | } | ||
Line 407: | Line 411: | ||
script inline @{gA} # <== new atoms added here | script inline @{gA} # <== new atoms added here | ||
− | |||
− | |||
appendNew = false | appendNew = false | ||
# Flip comp to comp strand | # Flip comp to comp strand | ||
if (double) { | if (double) { | ||
+ | f = (_frameID/1000000) | ||
+ | m = (_frameID%1000000) | ||
select @{("" + bResno + chb + "/" + f + "." + m)} | select @{("" + bResno + chb + "/" + f + "." + m)} | ||
var v1={8.238, 2.809, 6.004} | var v1={8.238, 2.809, 6.004} | ||
Line 424: | Line 428: | ||
# Set the angles between the new NT and the old NTs | # Set the angles between the new NT and the old NTs | ||
select (((@cha and (atomno < nNa)) or (@chb and (resno != bResno))) | select (((@cha and (atomno < nNa)) or (@chb and (resno != bResno))) | ||
− | and | + | and thisModel) |
set_angle_no(nNa, pNa+8, pNa+7, 120.0) | set_angle_no(nNa, pNa+8, pNa+7, 120.0) | ||
select (((@cha and (atomno < @{nNa+3})) | select (((@cha and (atomno < @{nNa+3})) | ||
− | or (@chb and (resno != bResno))) and | + | or (@chb and (resno != bResno))) and thisModel) |
set_dihedral_no(nNa+4, nNa+3, nNa, pNa+8, kC5O5PO3) | set_dihedral_no(nNa+4, nNa+3, nNa, pNa+8, kC5O5PO3) | ||
select (((@cha and (atomno < nNa)) or (@chb and (resno != bResno))) | select (((@cha and (atomno < nNa)) or (@chb and (resno != bResno))) | ||
− | and | + | and thisModel) |
set_dihedral_no(nNa+3, nNa, pNa+8, pNa+7, kO5PO3C3) | set_dihedral_no(nNa+3, nNa, pNa+8, pNa+7, kO5PO3C3) | ||
Line 449: | Line 453: | ||
appendNew = gAppendNew | appendNew = gAppendNew | ||
echo | echo | ||
− | select ( | + | select (thisModel) |
refresh | refresh | ||
Latest revision as of 17:02, 12 April 2016
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 resulting polynucleotide defaults to an open B-form. If the script "to_ab_nt" described here is available, it prompts to converts to a regular B-form if DNA or to a regular A-form if RNA or mixed.
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 first character is '+' then the polynucleotide is created in a new frame. 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. If prepended with 'M' then a mixed helix is produced where the first strand is DNA and the second RNA. Multiple prepends are allowed (though 'M' would be inconsistent with 'R' or 'S').
If the 3d character is ':' (4th if the first is '+') then the two chains are labeled by the two preceding characters instead of the default 'A' and 'B'. Likewise if the 2d character is ':' then the presumably single chain is labeled by the preceding single character.
With this script you can generate a polynucleotide helix chain from a sequence string in 1-char NT coding. You can optionally give it chain labels other than the default :A :B. You can also add to an existing chain(s), add new chains to an existing model or now create the chain(s) in a new frame by prepending a '+' to the sequence string.
Chains start from the origin and extend along the -X+Y+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 left along the YZ axis until the origin is clear.
Note that a single chain helix could then be added to a double chain or RNA to DNA or whatever. Have fun...
The IUPAC/IUBMB 1 letter code is used: A=Adenine C=Cytosine G=Guanine T=Thymine U=Uracil
The top level function plicoGenNt prompts the user for input.
The top level function plicoGenHelix accepts a string as a parameter.
Polymeraze 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 Polynucleotide Script=script <path to your script directory>/polymeraze.spt;plico_gen_nt
saved as plicoGenNT.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 polymeraze.spt.
# POLYMERAZE - Jmol script by Ron Mignery # v1.16 beta 4/12/2016 -axis is now a reserved word # # 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 resulting polynucleotide defaults to an open B-form. If double-stranded # and the script "to_ab_nt" is available, it prompts to converts to a regular # B-form if DNA or to a regular A-form if RNA or mixed. # # 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 1st character is '+' then the chain is created in a new frame. # If prepended with '3' then the string is considered as 3' to 5' # If prepended with 'R' then RNA is generated instead of DNA (Ts convert to Us) # If prepended with 'S' then a single strand helix is produced # If prepended with 'M' then a mixed helix is produced where the first # strand is DNA and the second RNA - multiple prepends are allowed # though 'M' is inconsistent with 'R' or 'S' # # If the 3d (4th if 1st is '+') character is ':' then the two chains are labeled # by the two preceding characters instead of the default 'A' and 'B' # Likewise if the 2d character is ':' then the presumably single chain is # labeled by the single preceding character # # The IUPAC/IUBMB 1 letter code is used: # A=Adenine C=Cytosine G=Guanine T=Thymine U=Uracil kPolymeraze = 1 # The following constant values determine the pitch of the helices kC5O5PO3 = -27.0 kO5PO3C3 = -117.8 kPO3C3C4 = -171.9 kO3C3C4C5 = 121 kC3C4C5O5 = 54 kC4C5O5P = 164 kPu = 65 kPy = 52 gChain1 = 'A' # The default chain id gChain2 = 'B' # The default complementary chain id gA = "" gSeq = "" gAppendNew = false gNewFrame = false # Lookup 3 letter code from 1 letter code kNt3from1 = {"A":" DA", "C":" DC", "G":" DG", "T":" DT", "U":" DU", "D":" DD", "X":" DX"} kNtComp = {"A":"T", "C":"G", "G":"C", "T":"A", "U":"A", "D":"G", "X":"C"} # Generate PDB atom record # Writes gNa or gNb function gen_atom_nt(atomname, group, resno, xyz, comp) { # Fixed column format: #ATOM 500 O4' DA B 29 -3.745 7.211 45.474 while (atomname.size < 3) { atomname += " "; } var a = format("ATOM %5d %4s %3s ", (comp ? gNb : gNa), atomname, group ) a += format("%s%4d %8.3f", (comp ? gChain2 : gChain1), resno, 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 gen_atom_nt that writes gNa or gNb function gen_nt(i, nt, rna, comp) { # From constructed nucleotides var P0 = [0.000, 0.000, 0.000] var OP1= [-0.973,0.363,-1.067] var OP2= [0.297,-1.428, 0.272] var O5p= [1.351, 0.795,-0.286] var C5p= [1.345, 2.211,-0.125] var C4p= [2.732, 2.786,-0.255] var O4p= [3.413, 2.900, 1.019] var C3p= [3.670, 2.020,-1.178] var O3p= [4.269, 2.960,-2.051] var C2p= [4.717, 1.445,-0.238] var O2p= [6.046, 1.365,-0.884] var C1p= [4.758, 2.505, 0.846] var N1ct= [5.277, 2.056, 2.143] var C2ct= [6.236, 2.836, 2.740] var O2ct= [6.670, 3.853, 2.230] var N3ct= [6.674, 2.381, 3.958] var C4ct= [6.256, 1.245, 4.622] var NO4ct=[6.726, 0.972, 5.728] var C5ct= [5.255, 0.455, 3.924] var C6ct= [4.820, 0.900, 2.737] var nC7ct=[4.762,-0.811, 4.551] var N9ag= [5.256, 2.091, 2.152] var C8ag= [4.867, 1.016, 2.913] var N7ag= [5.532, 0.894, 4.035] var C5ag= [6.425, 1.959, 4.013] var C6ag= [7.401, 2.391, 4.922] var NO6ag=[7.656, 1.780, 6.081] var N1ag= [8.118, 3.493, 4.599] var C2ag= [7.865, 4.104, 3.438] var nN2ag=[8.616, 5.197, 3.181] var N3ag= [6.968, 3.796, 2.503] var C4ag= [6.271, 2.701, 2.856] # Build PDB atom records common to all NTs var n3 = kNt3from1[nt] if (n3 = "") { n3 = " D?" } if (rna) { if (n3 == " DD") { n3 = " D" } else { n3 = n3.replace('D', ' ') } } var a = gen_atom_nt(" P ", n3, i, P0, comp) a += gen_atom_nt(" OP1", n3, i, OP1, comp) a += gen_atom_nt(" OP2", n3, i, OP2, comp) a += gen_atom_nt(" O5'", n3, i, O5p, comp) a += gen_atom_nt(" C5'", n3, i, C5p, comp) a += gen_atom_nt(" C4'", n3, i, C4p, comp) a += gen_atom_nt(" O4'", n3, i, O4p, comp) a += gen_atom_nt(" C3'", n3, i, C3p, comp) a += gen_atom_nt(" O3'", n3, i, O3p, comp) a += gen_atom_nt(" C2'", n3, i, C2p, comp) if (rna) { a += gen_atom_nt(" O2'", n3, i, O2p, comp) } a += gen_atom_nt(" C1'", n3, i, C1p, comp) # Now add NT specific atom records switch (nt) { case 'A' : a += gen_atom_nt(" N9 ", n3, i, N9ag, comp) a += gen_atom_nt(" C8 ", n3, i, C8ag, comp) a += gen_atom_nt(" N7 ", n3, i, N7ag, comp) a += gen_atom_nt(" C5 ", n3, i, C5ag, comp) a += gen_atom_nt(" C6 ", n3, i, C6ag, comp) a += gen_atom_nt(" N6 ", n3, i, NO6ag, comp) a += gen_atom_nt(" N1 ", n3, i, N1ag, comp) a += gen_atom_nt(" C2 ", n3, i, C2ag, comp) a += gen_atom_nt(" N3 ", n3, i, N3ag, comp) a += gen_atom_nt(" C4 ", n3, i, C4ag, comp) break; case 'C' : a += gen_atom_nt(" N1 ", n3, i, N1ct, comp) a += gen_atom_nt(" C2 ", n3, i, C2ct, comp) a += gen_atom_nt(" O2 ", n3, i, O2ct, comp) a += gen_atom_nt(" N3 ", n3, i, N3ct, comp) a += gen_atom_nt(" C4 ", n3, i, C4ct, comp) a += gen_atom_nt(" N4 ", n3, i, NO4ct, comp) a += gen_atom_nt(" C5 ", n3, i, C5ct, comp) a += gen_atom_nt(" C6 ", n3, i, C6ct, comp) break; case 'X' : case 'G' : a += gen_atom_nt(" N9 ", n3, i, N9ag, comp) a += gen_atom_nt(" C8 ", n3, i, C8ag, comp) a += gen_atom_nt(" N7 ", n3, i, N7ag, comp) a += gen_atom_nt(" C5 ", n3, i, C5ag, comp) a += gen_atom_nt(" C6 ", n3, i, C6ag, comp) a += gen_atom_nt(" O6 ", n3, i, NO6ag, comp) a += gen_atom_nt(" N1 ", n3, i, N1ag, comp) a += gen_atom_nt(" C2 ", n3, i, C2ag, comp) a += gen_atom_nt(" N2 ", n3, i, nN2ag, comp) a += gen_atom_nt(" N3 ", n3, i, N3ag, comp) a += gen_atom_nt(" C4 ", n3, i, C4ag, comp) break; case 'T' : a += gen_atom_nt(" N1 ", n3, i, N1ct, comp) a += gen_atom_nt(" C2 ", n3, i, C2ct, comp) a += gen_atom_nt(" O2 ", n3, i, O2ct, comp) a += gen_atom_nt(" N3 ", n3, i, N3ct, comp) a += gen_atom_nt(" C4 ", n3, i, C4ct, comp) a += gen_atom_nt(" O4 ", n3, i, NO4ct, comp) a += gen_atom_nt(" C5 ", n3, i, C5ct, comp) a += gen_atom_nt(" C6 ", n3, i, C6ct, comp) a += gen_atom_nt(" C7 ", n3, i, nC7ct, comp) break; case 'D' : case 'U' : a += gen_atom_nt(" N1 ", n3, i, N1ct, comp) a += gen_atom_nt(" C2 ", n3, i, C2ct, comp) a += gen_atom_nt(" O2 ", n3, i, O2ct, comp) a += gen_atom_nt(" N3 ", n3, i, N3ct, comp) a += gen_atom_nt(" C4 ", n3, i, C4ct, comp) a += gen_atom_nt(" O4 ", n3, i, NO4ct, comp) a += gen_atom_nt(" C5 ", n3, i, C5ct, comp) a += gen_atom_nt(" 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 set_angle_no (a1no, a2no, a3no, toangle) { var c = gChain1 var a1 = {(atomno=a1no) and (chain=c) and thisModel} var a2 = {(atomno=a2no) and (chain=c) and thisModel} var a3 = {(atomno=a3no) and (chain=c) and thisModel} var v1 = (a1.xyz - a2.xyz) var v2 = (a3.xyz - a2.xyz) var caxis = cross(v1, v2) + a2.xyz var curangle = angle(a1, a2, a3) rotateselected @caxis @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 set_dihedral_no (a1no, a2no, a3no, a4no, toangle) { var c = gChain1 var a1 = {(atomno=a1no) and (chain=c) and thisModel} var a2 = {(atomno=a2no) and (chain=c) and thisModel} var a3 = {(atomno=a3no) and (chain=c) and thisModel} var a4 = {(atomno=a4no) and (chain=c) and thisModel} var curangle = angle(a1, a2, a3, a4) rotateselected @a2 @a3 @{toangle-curangle} } function count_atoms(seq, rna, start, finish) { var ntc = {"A":21, "C":20, "G":22, "T":20, "U":19} var cnt = 0 for (var i = start; i <= finish; i++) { cnt += (ntc[seq[i]] + (rna ? 1 : 0)) } return cnt } # Generate a helix function gen_helix_strand(reverse, drm, double) { var f = (_frameID/1000000) var m = (_frameID%1000000) var cha = ":" + gChain1 var chb = ":" + gChain2 var seq = "" if (reverse) { for (var i = gSeq.count; i > 0; i--) { if ("AUCGT".find(gSeq[i]) > 0) { seq += gSeq[i]%9999%0 } else { print format("Unknown nucleotide %s!", seq[i]) } } } else { for (var i = 1; i <= gSeq.count; i++) { if ("AUCGT".find(gSeq[i]) > 0) { seq += gSeq[i]%9999%0 } else { print format("Unknown nucleotide %s!", seq[i]) } } } var cSeq = "" if (double) { for (var i = seq.count; i > 0; i--) { cSeq += ((seq[i] == 'A') and (drm > 0)) ? "U" : kNtComp[seq[i]] } } var aAtomCount = count_atoms(seq, (drm == 1), 1, seq.count) var bAtomCount = count_atoms(cSeq, (drm > 0), 1, cSeq.count) # global new P atom index for chain A gNa = ((m > 0) ? {thisModel}.atomno.max + 1 : 1) gNb = 0 if (double) { gNb = (aAtomCount + bAtomCount - count_atoms(cSeq, (drm>0), cSeq.count, cSeq.count)) # last P in cSeq } var aResno = 1 var pNa = 1 # previous gNa gAppendNew = appendNew if (gNewFrame) { appendNew = true gNa = 1 } else { appendNew = false # If there is an atom at the origin while (m > 0) { var ai = {within(0.1, {0 0 0}) and thisModel} if (ai) { # If on the same chain if (ai[1].chain = gChain1) { # Add to existing strand echo "Adding to existing strand..." pNa = ai.atomno aResno = get_resno_max(gChain1) + 1 gNa = {(chain=gChain1) and thisModel}.atomno.max + 1 gNb += gNa # Bump up all B chain atomno and resno savNb = gNb gNb = aAtomCount + bAtomCount + gNa gA = "data \"append nt\"\n" # global PDB atom record for (var j = 1; j <= all.atomno.max; j++) { var aj = {(atomno=j) and thisModel and (chain == gChain2)} if (aj) { gA += gen_atom_nt(aj.atomName, aj.group, (aj.resno+seq.count+cSeq.count), array(aj.x, aj.y, aj.z), true) } } gA += "end \"append nt\"" delete @chb script inline @{gA} # <== new atoms added here gNb = savNb break; } # Else move all from the new chain rightward on X axis else { select {thisModel and (chain!=gChain1) and (chain!=gChain2)} translateselected {0, -20, -20 } } } else { break } } # endwhile } var bResno = aResno + seq.count + cSeq.count - 1 var nNa = gNa # new P var nNb = 0#bBase # new comp P # For each NT for (var i = 1; i <= seq.count; i++) { # Move polynucleotide O3p to bond distance 1.59 from new nt P var pO3 = { -0.759, 0.925, 1.048} if (double) { select ((@cha or @chb) and thisModel) } else { select (@cha and thisModel) } if ((i + aResno) > 2) { var nO3 = {@cha and (atomno=@{pNa+8}) and thisModel}.xyz var xyz = @{pO3 - nO3} if (xyz) { translateselected @xyz } } # Gen NT ================================================== gA = "data \"append nt\"\n" # global PDB atom record gA += gen_nt(aResno, seq[i], (drm == 1), false); # gNa updated if (double) { nNb = gNb var nti = cSeq.count-i+1 gA += gen_nt(bResno, cSeq[nti], (drm > 0), true); # gNb++ if (i > 0) { gNb -= count_atoms(cSeq, (drm>0), nti-1, nti) } } gA += "end \"append nt\"" script inline @{gA} # <== new atoms added here appendNew = false # Flip comp to comp strand if (double) { f = (_frameID/1000000) m = (_frameID%1000000) select @{("" + bResno + chb + "/" + f + "." + m)} var v1={8.238, 2.809, 6.004} var v2={8.461, 4.646, 4.125} rotateSelected @v2 @v1 180.0 } # If any older NTs if ((i + aResno) > 2) { # Set the angles between the new NT and the old NTs select (((@cha and (atomno < nNa)) or (@chb and (resno != bResno))) and thisModel) set_angle_no(nNa, pNa+8, pNa+7, 120.0) select (((@cha and (atomno < @{nNa+3})) or (@chb and (resno != bResno))) and thisModel) set_dihedral_no(nNa+4, nNa+3, nNa, pNa+8, kC5O5PO3) select (((@cha and (atomno < nNa)) or (@chb and (resno != bResno))) and thisModel) set_dihedral_no(nNa+3, nNa, pNa+8, pNa+7, kO5PO3C3) set_dihedral_no(nNa, pNa+8, pNa+7, pNa+5, kPO3C3C4) } # Step new and previous N aResno++; bResno-- pNa = nNa nNa = gNa; nNb = gNb } # Make the nucleotide bonds connect # Clean up appendNew = gAppendNew echo select (thisModel) refresh # Convert to A-form if RNA or mixed else B-form try { script $SCRIPT_PATH$toabNT.spt var s = format("Convert to %s-form?", ((drm > 0) ? "A" : "tight B")) var p = prompt(s, "Yes|No", true) if (p = "Yes") { to_ab_nt_auto(gChain1, (drm > 0)) } } catch { } } # Generate a helix or two function plico_gen_helix(seq) { if (gPlicoRecord != "") { var g = format("show file \"%s\"", gPlicoRecord) var ls = script(g) if (ls.find("FileNotFoundException")) { ls = "" } ls += format("plico_gen_helix(\"%s\");", gSeq) write var ls @gPlicoRecord } var single = false var reverse = false var drm = 0 var done = false gSeq = seq%9999%0 print format ("Seq=%s", gSeq) gNewFrame = false if (gSeq[1] == '+') { gNewFrame = true gSeq = gSeq[2][9999] } if (gSeq[2] == ':') { gChain1 = gSeq[1] gSeq = gSeq[3][9999] } else if (gSeq[3] == ':') { gChain1 = gSeq[1] gChain2 = gSeq[2] gSeq = gSeq[4][9999] } 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') { drm = 1; done = false; } else if (gSeq[1] == 'M') { drm = 2; done = false; } if (done == false) { gSeq = gSeq[2][9999] } } if (drm = 1) { gSeq = gSeq.replace('T', 'U') } else { gSeq = gSeq.replace('U', 'T') } # Generate gen_helix_strand(reverse, drm, single ? false : true) } function plico_gen_nt { echo Generating Nucleotide Helix # Get the sequence from the user var seq = prompt("Enter NT sequence (<+AB:3RSM>ACGTU...)", "")%9999%0 if ((seq != "NULL") and (seq.count > 0)) { plico_gen_helix(seq) } } # end of polymeraze.spt