User:Remig/plico/toABnt

From Jmol
< User:Remig‎ | plico
Revision as of 16:09, 3 April 2014 by Remig (talk | contribs) (Do not var globals)
Jump to navigation Jump to search

ToABnt allows the user to convert a polynucleotide chain from A-form to B-form and vice versa. When you ALT-CLICK on any atom, you are prompted to convert its chain from its current form to either the A-form or the B-form. The two forms differ primarily in the pucker conformation of their sugar moieties. The A-form has its sugars in the C3'-endo conformation whereas the B-form has it in the C2'-endo conformation. The A-form is characteristic of RNA while the B-form is characteristic of DNA.

This routine is called automatically by the Plico routine Polymeraze when a polynucleotide is generated.

ToABnt 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 toAB NT
Script=script <path to your script directory>/toABnt.spt;plicotoABnt

saved as toABnt.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 toABnt.spt.

#   toABnt - Jmol script by Ron Mignery
#   v1.1 beta    4/3/2014 for Jmol 14 -do not var globals for Jmol 14.0.13+
#
#   Convert a polynucleotide chain from A-form to B-form or vice versa
#
kC5O5PO3B = -71.0
kO5PO3C3B = -106.50
kPO3C3C4B = -160.67
kO3C3C4C5B = 125.44
kC3C4C5O5B = 55.65
kC4C5O5PB = 169.10

kO4C4C3C2B = 15.92
kC4O4C1C2B = -41.7
kC4O4C1NxB = -159.03
kC5C4O4C1B = 146.31
kC3C1C2O2B = -167.9

kPuB = 59.0
kPyB = 61.0

kC5O5PO3A = -59.3
kO5PO3C3A = -56.9
kPO3C3C4A = -157.4
kO3C3C4C5A = 75.5
kC3C4C5O5A = 49.55
kC4C5O5PA = 170.8

kO4C4C3C2A = -35.55
kC4O4C1C2A = 3.8
kC4O4C1NxA = -117.4
kC5C4O4C1A = 144.85
kC3C1C2O2A = 116.3

kPuA = 15
kPyA = 15

gToab = FALSE
gToabAlt = FALSE
gToA = FALSE

function getInterbaseRotors(stator, mover) {
    var rotors = array()
    var sRes = stator.resno
    var mRes = mover.resno
    var iChain = mover.chain
    var mC4 = {(resno=mRes) and (chain=iChain) and (atomName="C4\'")}
    var mC3 = {(resno=mRes) and (chain=iChain) and (atomName="C3\'")}
    var mO3 = {(resno=mRes) and (chain=iChain) and (atomName="O3\'")}
    var sP = {(resno=sRes) and (chain=iChain) and (atomName="P")}
    var sO5 = {(resno=sRes) and (chain=iChain) and (atomName="O5\'")}
    var sC5 = {(resno=sRes) and (chain=iChain) and (atomName="C5\'")}
    var sC4 = {(resno=sRes) and (chain=iChain) and (atomName="C4\'")}
    var sC3 = {(resno=sRes) and (chain=iChain) and (atomName="C3\'")}

    rotors += [mC4.atomIndex, mC3.atomIndex, mO3.atomIndex, sP.atomIndex]
    rotors += [mC3.atomIndex, mO3.atomIndex, sP.atomIndex, sO5.atomIndex]
    rotors += [mO3.atomIndex, sP.atomIndex, sO5.atomIndex, sC5.atomIndex]
    rotors += [sP.atomIndex, sO5.atomIndex, sC5.atomIndex, sC4.atomIndex]
    rotors = [sO5.atomIndex, sC5.atomIndex, sC4.atomIndex, sC3.atomIndex]
    return rotors
}

function setResidueDistance(stator, mover, dist) {
    var selsave = {selected}
    var cp = mover.xyz
    select mover
    setDistanceAtoms(stator, mover, dist)
    var pt = mover.xyz
    mover.xyz = cp
    select selsave
    var rotors = getInterbaseRotors(stator, mover)
    toabTrackIdx(mover.atomIndex, pt, rotors)
    toabTrackIdx(mover.atomIndex, pt, rotors)
    toabTrackIdx(mover.atomIndex, pt, rotors)
}

# Rotate rotor set to move target atom to its proper place
function toabTrackIdx(targetIdx, targetPt, iRotors) {
    gOK = FALSE
    var pt = targetPt
    var dist = distance(pt, {atomIndex=targetIdx}.xyz)
    var baseSet = {selected}
    var rotors = iRotors

    # For a number of passes
    for (var pass1 = 0; pass1 < 20; pass1++) {
        var blocked = ({})
        for (var pass2 = 0; pass2 < (rotors.size/4); pass2++) {

            var v1 = {atomIndex=targetIdx}.xyz - pt

            # Find the most orthgonal unused rotor
            var imax = 0
            var smax = 0.5
            for (var ri = 1; ri < rotors.size; ri += 4) {
                var i2 = rotors[ri+1]
                var i3 = rotors[ri+2]
                var i4 = rotors[ri+3]
                if ((i2 != targetIdx) and (i3 != targetIdx) and (i4 != targetIdx)) {
                    if ({blocked and {atomIndex=i2}}.count == 0) {
                        var v2 = {atomIndex=i3}.xyz - {atomIndex=i2}.xyz

                        var s = sin(abs(angle(v1, {0 0 0}, v2)))
                        if (s > smax) {
                            smax = s
                            imax = ri
                        }
                    }
                }
            }

            # If no more rotors, break to next full try
            if (imax == 0) {
               break
            }
            var i1 = rotors[imax+0]
            var i2 = rotors[imax+1]
            var i3 = rotors[imax+2]
            var i4 = rotors[imax+3]

            # Get dihedral of rotor with target point
            var dt = angle({atomIndex=targetIdx}, {atomIndex=i2}, {atomIndex=i3}, pt)

            # Add rotor stators to select
            select baseSet
            for (var ri = 1; ri <= imax; ri += 4) {
                {atomIndex=@{rotors[ri]}}.selected = TRUE
                {atomIndex=@{rotors[ri+1]}}.selected = TRUE
            }

            # Rotate to minimize vector ====================
            rotateSelected {atomIndex=i2} {atomIndex=i3} @dt

            # If close enough, stop
            if (distance(pt, {atomIndex=targetIdx}) < (kDtolerance/4)) {
                gOK = TRUE
                gTargetPt = pt
                break
            }

            # Block rotor
            blocked |= {atomIndex=i2}

        }   # endfor num rotors passes

        if (gOK) {
            break
        }
    }   # endfor 20 passes
}

function toabNtAuto(iChain, toA) {

    # Load common functions if not already
    if (kCommon < 1) {
        script $SCRIPT_PATH$plicoCommon.spt
        if (kCommon < 1) {
            prompt ("A newer version of plicoCommon.SPT is required")
            quit
        }
    }

    gChain1 = iChain
    gToA = toA
    toabNt(TRUE)
}

function toabNt(auto) {
    select all
    g5cargoIdx = {(atomno=@{{chain=gChain1}.atomno.min})
        and (chain=gChain1)}.atomIndex
    g3cargoIdx = {(atomno=@{{chain=gChain1}.atomno.max})
        and (chain=gChain1)}.atomIndex

    # If new drag
    if (gNewDrag) {
        gNewDrag = FALSE
        save state gState
    }

    # Move the cargo
    select {gCargoSet}
    gChain1 = {atomIndex=g3cargoIdx}.chain
    gChain2 = ""

    # Move between A-form and B-form
    var res5 = {atomIndex=g5cargoIdx}.resno
    var res3 = {atomIndex=g3cargoIdx}.resno
    var pSet = ({})
    var bSet = ({})
    var k = -1
    var ang = 0
    var pang = 0
    var isP = FALSE
    for (var i = res5; i <= res3; i++) {
        var j = i-1
        var aO3 = {(resno=i) and (chain=gChain1) and (atomName="O3\'")}
        var aC3 = {(resno=i) and (chain=gChain1) and (atomName="C3\'")}
        var aC4 = {(resno=i) and (chain=gChain1) and (atomName="C4\'")}
        var aC5 = {(resno=i) and (chain=gChain1) and (atomName="C5\'")}
        var aO5 = {(resno=i) and (chain=gChain1) and (atomName="O5\'")}
        var aP =  {(resno=i) and (chain=gChain1) and (atomName="P")}
        if (i == res5) {
            if (auto == FALSE) {
                isForm = "A"
                if (angle(aO3, aC3, aC4, aC5) > ((kO3C3C4C5A + kO3C3C4C5B)/2)) {
                    isForm = "B"
                }
                p = prompt(format("Convert chain %s from %s-form to:", gChain1, isForm),
                    "A-form|B-form", TRUE)
                if (p == "A-form") {
                    gToA = TRUE
                }
                else if (p == "B-form") {
                    gToA = FALSE
                }
                else {
                    color {all} @gScheme
                    break
                }
                print format("Converting to %s", p)
            }

            # If paired
            lcAtoms = (within(3.0, FALSE,
                {(resno=i) and (chain=gChain1) and base})
                and not {chain=gChain1})
            isP = (lcAtoms.size > 0)
            gChain2 = (isP ? lcAtoms[1].chain : "")
            k = (isP ? lcAtoms[1].resno : -1)
        }
        pSet = {(resno=k) and (chain=gChain2)}
        var bO3 = {(resno=k) and (chain=gChain2) and (atomName="O3\'")}
        var bC3 = {(resno=k) and (chain=gChain2) and (atomName="C3\'")}
        var bC4 = {(resno=k) and (chain=gChain2) and (atomName="C4\'")}
        var bC5 = {(resno=k) and (chain=gChain2) and (atomName="C5\'")}
        var bO5 = {(resno=k) and (chain=gChain2) and (atomName="O5\'")}
        var bP =  {(resno=k) and (chain=gChain2) and (atomName="P")}

        # If not 5' terminus
        var aC4p = ({})
        var aC3p = ({})
        var aO3p = ({})
        var bO5p = ({})
        var bC4p = ({})
        var bC3p = ({})
        var bO3p = ({})
        var bPp = ({})
        if (i > res5) {
            aC4p = {(resno=j) and (chain=gChain1) and (atomName="C4\'")}
            aC3p = {(resno=j) and (chain=gChain1) and (atomName="C3\'")}
            aO3p = {(resno=j) and (chain=gChain1) and (atomName="O3\'")}

            bO5p = {(resno=@{k+1}) and (chain=gChain2) and (atomName="O5\'")}
            bC4p = {(resno=@{k+1}) and (chain=gChain2) and (atomName="C4\'")}
            bC3p = {(resno=@{k+1}) and (chain=gChain2) and (atomName="C3\'")}
            bO3p = {(resno=@{k+1}) and (chain=gChain2) and (atomName="O3\'")}
            bPp = {(resno=@{k+1}) and (chain=gChain2) and (atomName="P")}

            select bSet or {(atomno < @{aP.atomno}) and
                (atomno >= @{{atomIndex=g5cargoIdx}.atomno}) and (chain=gChain1)}
            setDihedralAtoms(aP, aO3p, aC3p, aC4p, (gToA ? kPO3C3C4A : kPO3C3C4B))

            setDihedralAtoms(aO5, aP, aO3p, aC3p, (gToA ? kO5PO3C3A : kO5PO3C3B))

            select bSet or {(atomno < @{aO5.atomno}) and
                (atomno >= @{{atomIndex=g5cargoIdx}.atomno}) and (chain=gChain1)}
            setDihedralAtoms(aC5, aO5, aP, aO3p, (gToA ? kC5O5PO3A : kC5O5PO3B))
        }
        select bSet or {(atomno < @{aO5.atomno}) and
            (atomno >= @{{atomIndex=g5cargoIdx}.atomno}) and (chain=gChain1)}
        setDihedralAtoms(aC4, aC5, aO5, aP, (gToA ? kC4C5O5PA : kC4C5O5PB))
        select (bP or (connected(bP) and {resno=k}))
        setDihedralAtoms(bC4, bC5, bO5, bP, (gToA ? kC4C5O5PA : kC4C5O5PB))

        select bSet or {(atomno < @{aC5.atomno}) and
            (atomno >= @{{atomIndex=g5cargoIdx}.atomno}) and (chain=gChain1)}
        setDihedralAtoms(aC3, aC4, aC5, aO5, (gToA ? kC3C4C5O5A : kC3C4C5O5B))
        select (bP or (connected(bP) and {resno=k}))
        setDihedralAtoms(bC3, bC4, bC5, bO5, (gToA ? kC3C4C5O5A : kC3C4C5O5B))

        select (({(resno >= k) and (chain=gChain2)}
            or {(resno <= i) and (chain=gChain1)}) and not aO3)
        setDihedralAtoms(aO3, aC3, aC4, aC5, (gToA ? kO3C3C4C5A : kO3C3C4C5B))
        select bO3
        setDihedralAtoms(bC5, bC4, bC3, bO3, (gToA ? kO3C3C4C5A : kO3C3C4C5B))


        var aC1 =  {(resno=i) and (chain=gChain1) and (atomName="C1\'")}
        var aC2 =  {(resno=i) and (chain=gChain1) and (atomName="C2\'")}
        var aO2 =  {(resno=i) and (chain=gChain1) and (atomName="O2\'")}
        var aO4 =  {(resno=i) and (chain=gChain1) and (atomName="O4\'")}

        var bC1 =  {(resno=k) and (chain=gChain2) and (atomName="C1\'")}
        var bC2 =  {(resno=k) and (chain=gChain2) and (atomName="C2\'")}
        var bO2 =  {(resno=k) and (chain=gChain2) and (atomName="O2\'")}
        var bO4 =  {(resno=k) and (chain=gChain2) and (atomName="O4\'")}

        # Set chi
        var aNx = -1
        var aCx = -1
        var bNx = -1
        var bCx = -1
        if ((aC1 and {purine}).size > 0) {
            aNx =  {(resno=i) and (chain=gChain1) and (atomName="N9")}
            aCx =  {(resno=i) and (chain=gChain1) and (atomName="C8")}
            bNx =  {(resno=k) and (chain=gChain2) and (atomName="N1")}
            bCx =  {(resno=k) and (chain=gChain2) and (atomName="C6")}
            ang = (gToA ? kPuA : kPuB)
            pang = (gToA ? kPyA : kPyB)
        }
        else {
            aNx =  {(resno=i) and (chain=gChain1) and (atomName="N1")}
            aCx =  {(resno=i) and (chain=gChain1) and (atomName="C6")}
            bNx =  {(resno=k) and (chain=gChain2) and (atomName="N9")}
            bCx =  {(resno=k) and (chain=gChain2) and (atomName="C8")}
            ang = (gToA ? kPyA : kPyB)
            pang = (gToA ? kPuA : kPuB)
        }
        select pSet or {(resno=i) and (chain=gChain1) and base}
        setDihedralAtoms(aO4, aC1, aNx, aCx, ang)
        select {pSet and not base}
        setDihedralAtoms(bCx, bNx, bC1, bO4, pang)

        # Set pucker 3' endo or 2' endo
        select pSet or {(resno=i) and (chain=gChain1) and base}
        setDihedralAtoms(aC4, aO4, aC1, aNx, (gToA ? kC4O4C1NxA : kC4O4C1NxB))
        select {pSet and not base}
        setDihedralAtoms(bNx, bC1, bO4, bC4, (gToA ? kC4O4C1NxA : kC4O4C1NxB))

        select (pSet or
            {((resno=i) and (chain=gChain1) and base) or aC1})
        setDihedralAtoms(aC5, aC4, aO4, aC1, (gToA ? kC5C4O4C1A : kC5C4O4C1B))
        select {pSet and not base and not bC1}
        setDihedralAtoms(bC1, bO4, bC4, bC5, (gToA ? kC5C4O4C1A : kC5C4O4C1B))

        select aC2 or aO2
        setDihedralAtoms(aC4, aO4, aC1, aC2, (gToA ? kC4O4C1C2A : kC4O4C1C2B))
        if (aO2.size > 0) {
            ang = (gToA ? kC3C1C2O2A : kC3C1C2O2B)
            setDihedralAtoms(aC3, aC1, aC2, aO2, (gToA ? kC3C1C2O2A : kC3C1C2O2B))
        }
        setDistanceAtoms(aC3, aC2, 1.52)
        setDistanceAtoms(aC1, aC2, 1.52)
        select bC2 or bO2
        setDihedralAtoms(bC4, bO4, bC1, bC2, (gToA ? kC4O4C1C2A : kC4O4C1C2B))
        if (bO2.size > 0) {
            ang = (gToA ? kC3C1C2O2A : kC3C1C2O2B)
            setDihedralAtoms(bC3, bC1, bC2, bO2, (gToA ? kC3C1C2O2A : kC3C1C2O2B))
        }
        setDistanceAtoms(bC3, bC2, 1.52)
        setDistanceAtoms(bC1, bC2, 1.52)

        # compress
        if (i > res5) {
            select ({resno<i} or {resno>k}) and not aP
            #setResidueDistance(aC4, aC4p, 4.0)
        }

        # If paired, make minor adjustments
        if (isP and (i > res5)) {

            var cp = bPp.xyz
            select bPp
            setDistanceAtoms(bO3, bPp, 1.59)
            setAngleAtoms(bC3, bO3, bPp, 118.2)
            setDihedralAtoms(bC4, bC3, bO3, bPp, (gToA ? kPO3C3C4A : kPO3C3C4B))
            var pt = bPp.xyz
            bPp.xyz = cp
            rotors = [aC4p.atomIndex, aC3p.atomIndex, aO3p.atomIndex, aP.atomIndex]
            rotors += [aC3p.atomIndex, aO3p.atomIndex, aP.atomIndex, aO5.atomIndex]
            rotors += [aO3p.atomIndex, aP.atomIndex, aO5.atomIndex, aC5.atomIndex]
            rotors += [aP.atomIndex, aO5.atomIndex, aC5.atomIndex, aC4.atomIndex]
            rotors += [aO5.atomIndex, aC5.atomIndex, aC4.atomIndex, aC3.atomIndex]
            select ({((resno < i) and (chain=gChain1))} or
                {(resno > k) and (chain=gChain2)})
            toabTrackIdx(bPp.atomIndex, pt, rotors)
            var bOP1 =  {(resno=@{bPp.resno}) and (chain=gChain2) and (atomName="OP1")}
            var bOP2 =  {(resno=@{bPp.resno}) and (chain=gChain2) and (atomName="OP2")}
            bOP1.xyz = getTetIdx(bO3.atomIndex, bPp.atomIndex, bO5p.atomIndex, 1.5)
            bOP2.xyz = getTetIdx(bO5p.atomIndex, bPp.atomIndex, bO3.atomIndex, 1.5)
            var aOP1 =  {(resno=@{aP.resno}) and (chain=gChain1) and (atomName="OP1")}
            var aOP2 =  {(resno=@{aP.resno}) and (chain=gChain1) and (atomName="OP2")}
            aOP1.xyz = getTetIdx(aO3p.atomIndex, aP.atomIndex, aO5.atomIndex, 1.5)
            aOP2.xyz = getTetIdx(aO5.atomIndex, aP.atomIndex, aO3p.atomIndex, 1.5)
        }

        bSet = bset or pSet
        if (k > 0) {
            k--
        }

    } # endfor

}

# Bound to ALT-LEFT-CLICK by plicoToabNt
function toabChainMB() {
    color {all} @gScheme
    gChain1 = {atomIndex=_atomPicked}.chain
    color {chain=gChain1} @gAltScheme
    refresh
    toabNt(FALSE)
}

# Top level of ToABnt
function plicoToabNT() {

    # Load common functions if not already
    if (kCommon < 1) {
        script $SCRIPT_PATH$plicoCommon.spt
        if (kCommon < 1) {
            prompt ("A newer version of plicoCommon.SPT is required")
            quit
        }
    }

    gPlico = "TO A-FORM/B-FORM"
    plicoPrelim()

    gEcho = ("________A <==> B_______|ALT-CLICK=mark chain|DOUBLE-CLICK=exit")
    echo @gEcho
    gChain = ""
    unbind

    bind "ALT-LEFT-CLICK" "_pickAtom";
    bind "ALT-LEFT-CLICK" "+:toabChainMB";
    bind "DOUBLE" "plicoExit";
}
# End of TOABNT.SPT

Contributors

Remig