User:Remig

My interest is in protein folding for which I have written software tools in Java over the years. I had always displayed in JMol and only now realize I could have written in JMol script directly and have folded and viewed in the same application. I am now doing so and will present scripts of possible interest to others in the discussion tab for the general amusement.

The script shown 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:

1. POLYMERAZE - Jmol script by Ron Mignery
2. v1.0 beta 10/31/2013
4. POLYMERAZE takes a string message encoding a nucleotide (nt) sequence
5. and generates a corresponding double helix one nt at a time from the
6. 5' terminus to the 3' terminus rotating the emerging helix as it goes.
8. The message is a string entered by the user at a prompt.
9. It may be typed in or pasted in and be of any length
10. If prepended with '3' then the string is considered as 3' to 5'
11. If prepended with 'R' then RNA is generated instead of DNA
12. If prepended with 'S' then a single strand helix is produced
14. Multiple runs append to the previous helix if any unless it is moved
16. The IUPAC/IUBMB 1 letter code is used:
17. A=Adenine C=Cytosine G=Guanine T=Thymine U=Uracil

1. 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

1. 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"}

1. Generate PDB atom record
2. 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

};

1. Generate a PDB nucleotide record set
2. 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

};

1. Rotate a1 on a2 in the plane of a1, a2 and a3 to the given angle
2. 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}

}

1. Set the dihedral to the given angle
2. a1 (or a4) and all connected except by a2 (or a3) must be selected
3. 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

}

1. 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)

   # 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 == cha[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))

}

1. 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)

}

1. ==============================================

echo Generating Alpha Helix

1. 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.

Contributors

Remig