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