Difference between revisions of "Recycling Corner/DNA Generator"

From Jmol
Jump to navigation Jump to search
(Fix RNA group names DU instead of U etc.)
(Avoid "axis," a newly reserved word)
 
(13 intermediate revisions by the same user not shown)
Line 1: Line 1:
 
'''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.
 
'''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 [[User:Remig/plico|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.
 
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.
 
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 '3' then the string is considered as 3' to 5'.
 
If prepended with 'R' then RNA is generated instead of DNA.
 
If prepended with 'R' then RNA is generated instead of DNA.
Line 9: Line 12:
 
Multiple prepends are allowed (though 'M' would be inconsistent with 'R' or 'S').  
 
Multiple prepends are allowed (though 'M' would be inconsistent with 'R' or 'S').  
  
If the 3d 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 preceding single character.
+
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.''
  
A polynucleotide may be added onto by subsequent runs if the previous helices are not moved away. Note that a single chain helix could then be added to a double chain or RNA to DNA or whatever. Have fun...
+
''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 IUPAC/IUBMB 1 letter code is used: A=Adenine C=Cytosine G=Guanine T=Thymine U=Uracil
Line 19: Line 26:
 
The top level function '''plicoGenHelix''' accepts a string as a parameter.
 
The top level function '''plicoGenHelix''' accepts a string as a parameter.
  
Polymeraze is a member of the Plico suite of protein folding tools. It may be installed and accessed as a macro with the file:
+
Polymeraze is a member of the Plico suite of protein folding tools described [[User:Remig/plico|here]]. It may be installed and accessed as a macro with the file:
 
<pre>Title=PLICO Generate Polynucleotide
 
<pre>Title=PLICO Generate Polynucleotide
Script=script <path to your script folder>/polymeraze.spt;plicoGenNT</pre>
+
Script=script <path to your script directory>/polymeraze.spt;plico_gen_nt</pre>
saved as plicoGenNT.macro in your .jmol/macros folder as described in [[Macro]].
+
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.
 
<pre>#  POLYMERAZE - Jmol script by Ron Mignery
 
<pre>#  POLYMERAZE - Jmol script by Ron Mignery
#  v1.4 beta    2/18/2014 - fix RNA residue mislabel
+
#  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
 
#  and generates a corresponding double helix one nt at a time from the
 
#  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.
 
#  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.
 
#  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
 
#  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 '3' then the string is considered as 3' to 5'
#  If prepended with 'R' then RNA is generated instead of DNA
+
#  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 'S' then a single strand helix is produced
 
#  If prepended with 'M' then a mixed helix is produced where the first
 
#  If prepended with 'M' then a mixed helix is produced where the first
 
#  strand is DNA and the second RNA - multiple prepends are allowed
 
#  strand is DNA and the second RNA - multiple prepends are allowed
#  though 'M' is inconsistent with 'R' or 'S'  
+
#  though 'M' is inconsistent with 'R' or 'S'
 
#
 
#
#  If the 3d character is ':' then the two chains are labeled by the
+
#  If the 3d (4th if 1st is '+') character is ':' then the two chains are labeled
#  two preceding characters instead of the default 'A' and 'B'
+
by the two preceding characters instead of the default 'A' and 'B'
 
#  Likewise if the 2d character is ':' then the presumably single chain is
 
#  Likewise if the 2d character is ':' then the presumably single chain is
 
#  labeled by the single preceding character
 
#  labeled by the single preceding character
Line 48: Line 61:
 
#  A=Adenine C=Cytosine G=Guanine T=Thymine U=Uracil
 
#  A=Adenine C=Cytosine G=Guanine T=Thymine U=Uracil
  
 +
kPolymeraze = 1
 
# The following constant values determine the pitch of the helices
 
# The following constant values determine the pitch of the helices
var kC5O5PO3 = -27.0
+
kC5O5PO3 = -27.0
var kO5PO3C3 = -117.8
+
kO5PO3C3 = -117.8
var kPO3C3C4 = -171.9
+
kPO3C3C4 = -171.9
var gCHAIN1 = 'A'    # The chain id
+
kO3C3C4C5 = 121
var gCHAIN2 = 'B'    # The complementary chain id
+
kC3C4C5O5 = 54
var gA = ""
+
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
 
# Lookup 3 letter code from 1 letter code
var kNt3from1 = {"A":" DA", "C":" DC", "G":" DG", "T":" DT", "U":" DU"}
+
kNt3from1 = {"A":" DA", "C":" DC", "G":" DG", "T":" DT", "U":" DU", "D":" DD", "X":" DX"}
var kNtComp = {"A":"T", "C":"G", "G":"C", "T":"A", "U":"A"}
+
kNtComp = {"A":"T", "C":"G", "G":"C", "T":"A", "U":"A", "D":"G", "X":"C"}
  
 
# Generate PDB atom record
 
# Generate PDB atom record
 
# Writes gNa or gNb
 
# Writes gNa or gNb
function genAtom(atomname, group, resno, xyz, comp) {
+
function gen_atom_nt(atomname, group, resno, xyz, comp) {
 
     # Fixed column format:
 
     # Fixed column format:
 
     #ATOM    500  O4'  DA B  29      -3.745  7.211  45.474
 
     #ATOM    500  O4'  DA B  29      -3.745  7.211  45.474
Line 69: Line 91:
 
     }
 
     }
 
     var a =  format("ATOM  %5d %4s %3s ", (comp ? gNb : gNa), atomname, group )
 
     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("%s%4d    %8.3f", (comp ? gChain2 : gChain1), resno, xyz[1] )
 
     a +=  format("%8.3f%8.3f\n", xyz[2], xyz[3] )
 
     a +=  format("%8.3f%8.3f\n", xyz[2], xyz[3] )
 
     if (comp) gNb++; else gNa++
 
     if (comp) gNb++; else gNa++
       
+
 
 
     return a
 
     return a
 
};
 
};
  
 
# Generate a PDB nucleotide record set
 
# Generate a PDB nucleotide record set
# Calls genAtom that writes gNa or gNb
+
# Calls gen_atom_nt that writes gNa or gNb
function genNT(i, nt, rna, comp) {
+
function gen_nt(i, nt, rna, comp) {
  
 
     # From constructed nucleotides
 
     # From constructed nucleotides
Line 93: Line 115:
 
     var O2p=  [6.046, 1.365,-0.884]
 
     var O2p=  [6.046, 1.365,-0.884]
 
     var C1p=  [4.758, 2.505, 0.846]
 
     var C1p=  [4.758, 2.505, 0.846]
   
+
 
 
     var N1ct= [5.277, 2.056, 2.143]
 
     var N1ct= [5.277, 2.056, 2.143]
 
     var C2ct= [6.236, 2.836, 2.740]
 
     var C2ct= [6.236, 2.836, 2.740]
Line 99: Line 121:
 
     var N3ct= [6.674, 2.381, 3.958]
 
     var N3ct= [6.674, 2.381, 3.958]
 
     var C4ct= [6.256, 1.245, 4.622]
 
     var C4ct= [6.256, 1.245, 4.622]
     var NO4ct [6.726, 0.972, 5.728]
+
     var NO4ct=[6.726, 0.972, 5.728]
 
     var C5ct= [5.255, 0.455, 3.924]
 
     var C5ct= [5.255, 0.455, 3.924]
 
     var C6ct= [4.820, 0.900, 2.737]
 
     var C6ct= [4.820, 0.900, 2.737]
     var nC7ct [4.762,-0.811, 4.551]
+
     var nC7ct=[4.762,-0.811, 4.551]
   
+
 
 
     var N9ag= [5.256, 2.091, 2.152]
 
     var N9ag= [5.256, 2.091, 2.152]
 
     var C8ag= [4.867, 1.016, 2.913]
 
     var C8ag= [4.867, 1.016, 2.913]
Line 112: Line 134:
 
     var N1ag= [8.118, 3.493, 4.599]
 
     var N1ag= [8.118, 3.493, 4.599]
 
     var C2ag= [7.865, 4.104, 3.438]
 
     var C2ag= [7.865, 4.104, 3.438]
     var nN2ag [8.616, 5.197, 3.181]
+
     var nN2ag=[8.616, 5.197, 3.181]
 
     var N3ag= [6.968, 3.796, 2.503]
 
     var N3ag= [6.968, 3.796, 2.503]
 
     var C4ag= [6.271, 2.701, 2.856]
 
     var C4ag= [6.271, 2.701, 2.856]
   
+
 
 
     # Build PDB atom records common to all NTs
 
     # Build PDB atom records common to all NTs
 
     var n3 = kNt3from1[nt]
 
     var n3 = kNt3from1[nt]
Line 122: Line 144:
 
     }
 
     }
 
     if (rna) {
 
     if (rna) {
         n3 = n3.replace("D", " ")
+
         if (n3 == " DD") {
 +
            n3 = "  D"
 +
        }
 +
        else {
 +
            n3 = n3.replace('D', ' ')
 +
        }
 
     }
 
     }
     var a = genAtom(" P  ", n3, i, P0, comp)
+
     var a = gen_atom_nt(" P  ", n3, i, P0, comp)
     a += genAtom(" OP1", n3, i, OP1, comp)
+
     a += gen_atom_nt(" OP1", n3, i, OP1, comp)
     a += genAtom(" OP2", n3, i, OP2, comp)
+
     a += gen_atom_nt(" OP2", n3, i, OP2, comp)
     a += genAtom(" O5'", n3, i, O5p, comp)
+
     a += gen_atom_nt(" O5'", n3, i, O5p, comp)
     a += genAtom(" C5'", n3, i, C5p, comp)
+
     a += gen_atom_nt(" C5'", n3, i, C5p, comp)
     a += genAtom(" C4'", n3, i, C4p, comp)
+
     a += gen_atom_nt(" C4'", n3, i, C4p, comp)
     a += genAtom(" O4'", n3, i, O4p, comp)
+
     a += gen_atom_nt(" O4'", n3, i, O4p, comp)
     a += genAtom(" C3'", n3, i, C3p, comp)
+
     a += gen_atom_nt(" C3'", n3, i, C3p, comp)
     a += genAtom(" O3'", n3, i, O3p, comp)
+
     a += gen_atom_nt(" O3'", n3, i, O3p, comp)
     a += genAtom(" C2'", n3, i, C2p, comp)
+
     a += gen_atom_nt(" C2'", n3, i, C2p, comp)
    a += genAtom(" C1'", n3, i, C1p, comp)
 
 
     if (rna) {
 
     if (rna) {
         a += genAtom(" O2'", n3, i, O2p, comp)
+
         a += gen_atom_nt(" O2'", n3, i, O2p, comp)
 
     }
 
     }
 +
    a += gen_atom_nt(" C1'", n3, i, C1p, comp)
  
 
     # Now add NT specific atom records
 
     # Now add NT specific atom records
 
     switch (nt) {
 
     switch (nt) {
 
     case 'A' :
 
     case 'A' :
         a += genAtom(" N9 ", n3, i, N9ag, comp)
+
         a += gen_atom_nt(" N9 ", n3, i, N9ag, comp)
         a += genAtom(" C8 ", n3, i, C8ag, comp)
+
         a += gen_atom_nt(" C8 ", n3, i, C8ag, comp)
         a += genAtom(" N7 ", n3, i, N7ag, comp)
+
         a += gen_atom_nt(" N7 ", n3, i, N7ag, comp)
         a += genAtom(" C5 ", n3, i, C5ag, comp)
+
         a += gen_atom_nt(" C5 ", n3, i, C5ag, comp)
         a += genAtom(" C6 ", n3, i, C6ag, comp)
+
         a += gen_atom_nt(" C6 ", n3, i, C6ag, comp)
         a += genAtom(" N6 ", n3, i, NO6ag, comp)
+
         a += gen_atom_nt(" N6 ", n3, i, NO6ag, comp)
         a += genAtom(" N1 ", n3, i, N1ag, comp)
+
         a += gen_atom_nt(" N1 ", n3, i, N1ag, comp)
         a += genAtom(" C2 ", n3, i, C2ag, comp)
+
         a += gen_atom_nt(" C2 ", n3, i, C2ag, comp)
         a += genAtom(" N3 ", n3, i, N3ag, comp)
+
         a += gen_atom_nt(" N3 ", n3, i, N3ag, comp)
         a += genAtom(" C4 ", n3, i, C4ag, comp)
+
         a += gen_atom_nt(" C4 ", n3, i, C4ag, comp)
 
         break;
 
         break;
 
     case 'C' :
 
     case 'C' :
         a += genAtom(" N1 ", n3, i, N1ct, comp)
+
         a += gen_atom_nt(" N1 ", n3, i, N1ct, comp)
         a += genAtom(" C2 ", n3, i, C2ct, comp)
+
         a += gen_atom_nt(" C2 ", n3, i, C2ct, comp)
         a += genAtom(" O2 ", n3, i, O2ct, comp)
+
         a += gen_atom_nt(" O2 ", n3, i, O2ct, comp)
         a += genAtom(" N3 ", n3, i, N3ct, comp)
+
         a += gen_atom_nt(" N3 ", n3, i, N3ct, comp)
         a += genAtom(" C4 ", n3, i, C4ct, comp)
+
         a += gen_atom_nt(" C4 ", n3, i, C4ct, comp)
         a += genAtom(" N4 ", n3, i, NO4ct, comp)
+
         a += gen_atom_nt(" N4 ", n3, i, NO4ct, comp)
         a += genAtom(" C5 ", n3, i, C5ct, comp)
+
         a += gen_atom_nt(" C5 ", n3, i, C5ct, comp)
         a += genAtom(" C6 ", n3, i, C6ct, comp)
+
         a += gen_atom_nt(" C6 ", n3, i, C6ct, comp)
 
         break;
 
         break;
 +
    case 'X' :
 
     case 'G' :
 
     case 'G' :
         a += genAtom(" N9 ", n3, i, N9ag, comp)
+
         a += gen_atom_nt(" N9 ", n3, i, N9ag, comp)
         a += genAtom(" C8 ", n3, i, C8ag, comp)
+
         a += gen_atom_nt(" C8 ", n3, i, C8ag, comp)
         a += genAtom(" N7 ", n3, i, N7ag, comp)
+
         a += gen_atom_nt(" N7 ", n3, i, N7ag, comp)
         a += genAtom(" C5 ", n3, i, C5ag, comp)
+
         a += gen_atom_nt(" C5 ", n3, i, C5ag, comp)
         a += genAtom(" C6 ", n3, i, C6ag, comp)
+
         a += gen_atom_nt(" C6 ", n3, i, C6ag, comp)
         a += genAtom(" O6 ", n3, i, NO6ag, comp)
+
         a += gen_atom_nt(" O6 ", n3, i, NO6ag, comp)
         a += genAtom(" N1 ", n3, i, N1ag, comp)
+
         a += gen_atom_nt(" N1 ", n3, i, N1ag, comp)
         a += genAtom(" C2 ", n3, i, C2ag, comp)
+
         a += gen_atom_nt(" C2 ", n3, i, C2ag, comp)
         a += genAtom(" N2 ", n3, i, nN2ag, comp)
+
         a += gen_atom_nt(" N2 ", n3, i, nN2ag, comp)
         a += genAtom(" N3 ", n3, i, N3ag, comp)
+
         a += gen_atom_nt(" N3 ", n3, i, N3ag, comp)
         a += genAtom(" C4 ", n3, i, C4ag, comp)
+
         a += gen_atom_nt(" C4 ", n3, i, C4ag, comp)
 
         break;
 
         break;
 
     case 'T' :
 
     case 'T' :
         a += genAtom(" N1 ", n3, i, N1ct, comp)
+
         a += gen_atom_nt(" N1 ", n3, i, N1ct, comp)
         a += genAtom(" C2 ", n3, i, C2ct, comp)
+
         a += gen_atom_nt(" C2 ", n3, i, C2ct, comp)
         a += genAtom(" O2 ", n3, i, O2ct, comp)
+
         a += gen_atom_nt(" O2 ", n3, i, O2ct, comp)
         a += genAtom(" N3 ", n3, i, N3ct, comp)
+
         a += gen_atom_nt(" N3 ", n3, i, N3ct, comp)
         a += genAtom(" C4 ", n3, i, C4ct, comp)
+
         a += gen_atom_nt(" C4 ", n3, i, C4ct, comp)
         a += genAtom(" O4 ", n3, i, NO4ct, comp)
+
         a += gen_atom_nt(" O4 ", n3, i, NO4ct, comp)
         a += genAtom(" C5 ", n3, i, C5ct, comp)
+
         a += gen_atom_nt(" C5 ", n3, i, C5ct, comp)
         a += genAtom(" C6 ", n3, i, C6ct, comp)
+
         a += gen_atom_nt(" C6 ", n3, i, C6ct, comp)
         a += genAtom(" C7 ", n3, i, nC7ct, comp)
+
         a += gen_atom_nt(" C7 ", n3, i, nC7ct, comp)
 
         break;
 
         break;
 +
    case 'D' :
 
     case 'U' :
 
     case 'U' :
         a += genAtom(" N1 ", n3, i, N1ct, comp)
+
         a += gen_atom_nt(" N1 ", n3, i, N1ct, comp)
         a += genAtom(" C2 ", n3, i, C2ct, comp)
+
         a += gen_atom_nt(" C2 ", n3, i, C2ct, comp)
         a += genAtom(" O2 ", n3, i, O2ct, comp)
+
         a += gen_atom_nt(" O2 ", n3, i, O2ct, comp)
         a += genAtom(" N3 ", n3, i, N3ct, comp)
+
         a += gen_atom_nt(" N3 ", n3, i, N3ct, comp)
         a += genAtom(" C4 ", n3, i, C4ct, comp)
+
         a += gen_atom_nt(" C4 ", n3, i, C4ct, comp)
         a += genAtom(" O4 ", n3, i, NO4ct, comp)
+
         a += gen_atom_nt(" O4 ", n3, i, NO4ct, comp)
         a += genAtom(" C5 ", n3, i, C5ct, comp)
+
         a += gen_atom_nt(" C5 ", n3, i, C5ct, comp)
         a += genAtom(" C6 ", n3, i, C6ct, comp)
+
         a += gen_atom_nt(" C6 ", n3, i, C6ct, comp)
 
         break;
 
         break;
 
     default :
 
     default :
Line 205: Line 234:
  
 
# Rotate a1 on a2 in the plane of a1, a2 and a3 to the given angle
 
# 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  
+
# a1 and all connected except by a2 must be selected
function setAngle (a1, a2, a3, toangle) {
+
function set_angle_no (a1no, a2no, a3no, toangle) {
     var v1={atomno = a1}.xyz - {atomno = a2}.xyz
+
     var c = gChain1
     var v2={atomno = a3}.xyz - {atomno = a2}.xyz
+
    var a1 = {(atomno=a1no) and (chain=c) and thisModel}
     var axis = cross(v1, v2) + {atomno = a2}.xyz
+
    var a2 =  {(atomno=a2no) and (chain=c) and thisModel}
     var curangle =  angle({atomno=a1}, {atomno=a2}, {atomno=a3})
+
     var a3 = {(atomno=a3no) and (chain=c) and thisModel}
     rotateselected @axis {atomno = a2} @{curangle-toangle}
+
    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
 
# Set the dihedral to the given angle
# a1 (or a4) and all connected except by a2 (or a3) must be selected  
+
# a1 (or a4) and all connected except by a2 (or a3) must be selected
 
# If selected < unselected ==> a2 < a3 and vice versa
 
# If selected < unselected ==> a2 < a3 and vice versa
function setDihedral (a1, a2, a3, a4, toangle) {
+
function set_dihedral_no (a1no, a2no, a3no, a4no, toangle) {
     var curangle angle({atomno=a1}, {atomno=a2}, {atomno=a3}, {atomno=a4})
+
     var c = gChain1
     rotateselected {atomno=a2} {atomno=a3} @{toangle-curangle}
+
    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 countAtoms(seq, rna, start, finish) {
+
function count_atoms(seq, rna, start, finish) {
 
     var ntc = {"A":21, "C":20, "G":22, "T":20, "U":19}
 
     var ntc = {"A":21, "C":20, "G":22, "T":20, "U":19}
 
     var cnt = 0
 
     var cnt = 0
Line 232: Line 270:
  
 
# Generate a helix
 
# Generate a helix
function genHelixStrand(gSeq, reverse, drm, double) {
+
function gen_helix_strand(reverse, drm, double) {
     var cha = ":" + gCHAIN1
+
    var f = (_frameID/1000000)
     var chb = ":" + gCHAIN2
+
    var m = (_frameID%1000000)
 +
     var cha = ":" + gChain1
 +
     var chb = ":" + gChain2
 
     var seq = ""
 
     var seq = ""
 
     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 248: Line 300:
 
     if (double) {
 
     if (double) {
 
         for (var i = seq.count; i > 0; i--) {
 
         for (var i = seq.count; i > 0; i--) {
             cSeq += ((seq[i] == 'A') and (dr > 0)) ? "U" : kNtComp[seq[i]]
+
             cSeq += ((seq[i] == 'A') and (drm > 0)) ? "U" : kNtComp[seq[i]]
 
         }
 
         }
 
     }
 
     }
     var aAtomCount = countAtoms(seq, (drm == 1), 1, seq.count)
+
     var aAtomCount = count_atoms(seq, (drm == 1), 1, seq.count)
     var bAtomCount = countAtoms(cSeq, (drm > 0), 1, cSeq.count)
+
     var bAtomCount = count_atoms(cSeq, (drm > 0), 1, cSeq.count)
   
+
 
     gNa = 1 # global new P atom index for chain A
+
     # global new P atom index for chain A
 +
    gNa = ((m > 0) ? {thisModel}.atomno.max + 1 : 1)
 
     gNb = 0
 
     gNb = 0
 
     if (double) {
 
     if (double) {
 
         gNb = (aAtomCount + bAtomCount
 
         gNb = (aAtomCount + bAtomCount
         - countAtoms(cSeq, (drm>0), cSeq.count, cSeq.count)) # last P in cSeq
+
         - count_atoms(cSeq, (drm>0), cSeq.count, cSeq.count)) # last P in cSeq
 
     }
 
     }
    #var bBase = (all.count + countAtoms(seq, (drm > 0), seq.count) + 1)
 
  
    # Find last linkable P if any
 
 
     var aResno = 1
 
     var aResno = 1
 
     var pNa = 1    # previous gNa
 
     var pNa = 1    # previous gNa
    for (var i = all.count; i  > 0; i--) {
 
  
         # If A strand found at {0,0,0}
+
    gAppendNew = appendNew
         if (distance({atomno=i}, {0,0,0}) < 0.1) {
+
    if (gNewFrame) {
            if ({atomno=i}.chain == gCHAIN1) {
+
        appendNew = true
           
+
        gNa = 1
                # Add to existing strand
+
    }
                echo "Adding to existing strand..."
+
    else {
                pNa = i
+
        appendNew = false
                aResno = {chain=gCHAIN1}.resno.max + 1
+
 
                gNa = {chain=gCHAIN1}.atomno.max + 1
+
         # If there is an atom at the origin
                gNb += gNa
+
         while (m > 0) {
               
+
 
                # Bump up all B chain atomno and resno
+
            var ai = {within(0.1, {0 0 0}) and thisModel}
                # KLUDGE to work-around of Jmol's lack of resno rewrite
+
            if (ai) {
                savNb = gNb
+
 
                gNb = aAtomCount + bAtomCount + gNa
+
                # If on the same chain
                gA = "data \"append nt\"\n"    # global PDB atom record
+
                if (ai[1].chain = gChain1) {
                for (j = 1; j <= all.atomno.max; j++) {
+
 
                    if ({atomno=j}.chain == gCHAIN2) {
+
                    # Add to existing strand
                        gA += genAtom({atomno=j}.atomName, {atomno=j}.group,
+
                    echo "Adding to existing strand..."
                            ({atomno=j}.resno+seq.count+cSeq.count),
+
                    pNa = ai.atomno
                            array({atomno=j}.x, {atomno=j}.y, {atomno=j}.z), true)
+
                    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 }
 
                 }
 
                 }
                gA += "end \"append nt\""
 
                delete @chb
 
                script inline @{gA} # <== new atoms added here
 
                gNb = savNb
 
                break;
 
 
             }
 
             }
         }
+
            else {
 +
                break
 +
            }
 +
         } # endwhile
 
     }
 
     }
 +
    var bResno = aResno + seq.count + cSeq.count - 1
  
    var bResno = aResno + seq.count + cSeq.count - 1
 
   
 
 
     var nNa = gNa    # new P
 
     var nNa = gNa    # new P
 
     var nNb = 0#bBase  # new comp P
 
     var nNb = 0#bBase  # new comp P
  
 
     # For each NT
 
     # For each NT
    set appendnew false
 
 
     for (var i = 1; i <= seq.count; i++) {
 
     for (var i = 1; i <= seq.count; i++) {
          
+
 
         if (seq[i] == "") {
+
         # Move polynucleotide O3p to bond distance 1.59 from new nt P
             continue
+
        var pO3 = {  -0.759, 0.925, 1.048}
 +
         if (double) {
 +
            select ((@cha or @chb) and thisModel)
 +
        }
 +
        else {
 +
             select (@cha and thisModel)
 
         }
 
         }
         
 
        # Move polynucleotide O3p to bond distance from new nt P
 
        var pO3 = {-0.521, 0.638, 1.234}
 
        select all
 
 
         if ((i + aResno) > 2) {
 
         if ((i + aResno) > 2) {
             var nO3 =  {atomno=@{pNa+8}}.xyz
+
             var nO3 =  {@cha and (atomno=@{pNa+8})
 +
                and thisModel}.xyz
 
             var xyz = @{pO3 - nO3}
 
             var xyz = @{pO3 - nO3}
             translateselected @xyz
+
             if (xyz) {
 +
                translateselected @xyz
 +
            }
 
         }
 
         }
       
+
 
 
         # Gen NT ==================================================
 
         # Gen NT ==================================================
 
         gA = "data \"append nt\"\n"    # global PDB atom record
 
         gA = "data \"append nt\"\n"    # global PDB atom record
         gA += genNT(aResno, seq[i], (drm == 1), FALSE); # gNa updated
+
         gA += gen_nt(aResno, seq[i], (drm == 1), false); # gNa updated
 
         if (double) {
 
         if (double) {
 
             nNb = gNb
 
             nNb = gNb
 
             var nti = cSeq.count-i+1
 
             var nti = cSeq.count-i+1
             gA += genNT(bResno, cSeq[nti], (drm > 0), TRUE); # gNb++
+
             gA += gen_nt(bResno, cSeq[nti], (drm > 0), true); # gNb++
 
             if (i > 0) {
 
             if (i > 0) {
                 gNb -= countAtoms(cSeq, (drm>0), nti-1, nti)
+
                 gNb -= count_atoms(cSeq, (drm>0), nti-1, nti)
 
             }
 
             }
 
         }
 
         }
 
         gA += "end \"append nt\""
 
         gA += "end \"append nt\""
 
         script inline @{gA} # <== new atoms added here
 
         script inline @{gA} # <== new atoms added here
          
+
 
 +
         appendNew = false
 +
 
 
         # Flip comp to comp strand
 
         # Flip comp to comp strand
 
         if (double) {
 
         if (double) {
             select @{"" + bResno + chb}
+
            f = (_frameID/1000000)
 +
            m = (_frameID%1000000)
 +
             select @{("" + bResno + chb + "/" + f + "." + m)}
 
             var v1={8.238, 2.809, 6.004}
 
             var v1={8.238, 2.809, 6.004}
 
             var v2={8.461, 4.646, 4.125}
 
             var v2={8.461, 4.646, 4.125}
 
             rotateSelected @v2 @v1 180.0
 
             rotateSelected @v2 @v1 180.0
 
         }
 
         }
       
+
 
 
         # If any older NTs
 
         # If any older NTs
 
         if ((i + aResno) > 2) {
 
         if ((i + aResno) > 2) {
 +
 
             # 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)))
             setAngle(nNa, pNa+8, pNa+7, 120.0)
+
                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+3}) or (@chb and (resno != bResno)))
+
             select (((@cha and (atomno < nNa)) or (@chb and (resno != bResno)))
            setDihedral(nNa+4, nNa+3, nNa, pNa+8, kC5O5PO3)
+
                and thisModel)
           
+
             set_dihedral_no(nNa+3, nNa, pNa+8, pNa+7, kO5PO3C3)
            select (@cha and (atomno < nNa) or (@chb and (resno != bResno)))
+
 
             setDihedral(nNa+3, nNa, pNa+8, pNa+7, kO5PO3C3)
+
             set_dihedral_no(nNa, pNa+8, pNa+7, pNa+5, kPO3C3C4)
           
 
             setDihedral(nNa, pNa+8, pNa+7, pNa+5, kPO3C3C4)
 
 
         }
 
         }
       
+
 
 
         # Step new and previous N
 
         # Step new and previous N
 
         aResno++; bResno--
 
         aResno++; bResno--
Line 363: Line 446:
 
         nNa = gNa; nNb = gNb
 
         nNa = gNa; nNb = gNb
 
     }
 
     }
   
+
 
 
     # Make the nucleotide bonds
 
     # Make the nucleotide bonds
 
     connect
 
     connect
   
+
 
 
     # Clean up
 
     # Clean up
     select all
+
    appendNew = gAppendNew
     print format("%d atoms generated for chain %s", gNa+gNb,
+
    echo
         (comp ? gCHAIN2 : gCHAIN1))
+
     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
 
# Generate a helix or two
function plicoGenHelix(gSeq) {
+
function plico_gen_helix(seq) {
  
 
     if (gPlicoRecord != "") {
 
     if (gPlicoRecord != "") {
Line 382: Line 478:
 
             ls = ""
 
             ls = ""
 
         }
 
         }
         ls += format("plicoGenHelix(\"%s\");", gSeq)
+
         ls += format("plico_gen_helix(\"%s\");", gSeq)
 
         write var ls @gPlicoRecord
 
         write var ls @gPlicoRecord
 
     }
 
     }
   
+
 
     var single = FALSE
+
     var single = false
     var reverse = FALSE
+
     var reverse = false
 
     var drm = 0
 
     var drm = 0
     var done = FALSE
+
     var done = false
     gSeq = gSeq%9999%0
+
     gSeq = seq%9999%0
 
     print format ("Seq=%s", gSeq)
 
     print format ("Seq=%s", gSeq)
      
+
 
 +
     gNewFrame = false
 +
    if (gSeq[1] == '+') {
 +
        gNewFrame = true
 +
        gSeq = gSeq[2][9999]
 +
    }
 
     if (gSeq[2] == ':') {
 
     if (gSeq[2] == ':') {
         gCHAIN1 = gSeq[1]
+
         gChain1 = gSeq[1]
         gSeq[1] = ' '; gSeq[2] = ' '
+
         gSeq = gSeq[3][9999]
        gSeq = gSeq%0
 
 
     }
 
     }
 
     else if (gSeq[3] == ':') {
 
     else if (gSeq[3] == ':') {
         gCHAIN1 = gSeq[1]
+
         gChain1 = gSeq[1]
         gCHAIN2 = gSeq[2]
+
         gChain2 = gSeq[2]
         gSeq[1] = ' '; gSeq[2] = ' '; gSeq[3] = ' '
+
         gSeq = gSeq[4][9999]
        gSeq = gSeq%0
 
 
     }
 
     }
     while (done == FALSE) {
+
     while (done == false) {
         done = TRUE;
+
         done = true;
 
         if (gSeq[1] == 'S') {
 
         if (gSeq[1] == 'S') {
             single = TRUE;
+
             single = true;
             done = FALSE;
+
             done = false;
 
         }
 
         }
 
         else if (gSeq[1] == '3') {
 
         else if (gSeq[1] == '3') {
             reverse = TRUE;
+
             reverse = true;
             done = FALSE;
+
             done = false;
 
         }
 
         }
 
         else if (gSeq[1] == 'R') {
 
         else if (gSeq[1] == 'R') {
 
             drm = 1;
 
             drm = 1;
             done = FALSE;
+
             done = false;
 
         }
 
         }
 
         else if (gSeq[1] == 'M') {
 
         else if (gSeq[1] == 'M') {
 
             drm = 2;
 
             drm = 2;
             done = FALSE;
+
             done = false;
 
         }
 
         }
         if (done == FALSE) {
+
         if (done == false) {
             gSeq[1] = ' '
+
             gSeq = gSeq[2][9999]
            gSeq = gSeq%0
 
 
         }
 
         }
 
     }
 
     }
      
+
 
 +
     if (drm = 1) {
 +
        gSeq = gSeq.replace('T', 'U')
 +
    }
 +
    else {
 +
        gSeq = gSeq.replace('U', 'T')
 +
    }
 +
 
 
     # Generate
 
     # Generate
     genHelixStrand(gSeq, reverse, drm, single ? FALSE : TRUE)
+
     gen_helix_strand(reverse, drm, single ? false : true)
  
 
}
 
}
  
function plicoGenNT {
+
function plico_gen_nt {
     echo Generating Nuleotide Helix
+
     echo Generating Nucleotide Helix
   
+
 
 
     # Get the sequence from the user
 
     # Get the sequence from the user
     gSeq = prompt("Enter NT sequence (<3RSM>ACGTU)", "")%9999%0
+
     var seq = prompt("Enter NT sequence (<+AB:3RSM>ACGTU...)", "")%9999%0
     if ((gSeq != "NULL") and (gSeq.count > 0)) {
+
     if ((seq != "NULL") and (seq.count > 0)) {
         plicoGenHelix(gSeq)
+
         plico_gen_helix(seq)
 
     }
 
     }
 
}
 
}
# end of polymeraze.spt
+
# end of polymeraze.spt</pre>
</pre>
 

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

Contributors

Remig