Difference between revisions of "User:Remig"

From Jmol
Jump to navigation Jump to search
(DNA helix generator script)
(Add links)
 
(4 intermediate revisions by the same user not shown)
Line 1: Line 1:
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.
+
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 for the general amusement. Some are now available here [[User:Remig/plico]].
 
 
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:
 
 
 
<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.
 

Latest revision as of 21:52, 4 February 2014

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 for the general amusement. Some are now available here User:Remig/plico.

Contributors

Remig