User:Remig/plico/plicoCommonNT

This script contains routines used by some other scripts of the Plico suite involved with polynucleotide manipulation. It must be located in the same directory as any script that uses these routines.

Copy and paste the following into a text editor and save in your scripts folder as plicoNTcommon.spt.

#   plicoNTcommon - Jmol script by Ron Mignery
#   v1.3 beta    6/11/2014 -bug fix
#
#   Routines and values common to Plico suite scripts that work with nucleotides
#   Must be present in the same directory as other Plico scripts that use it
kNTcommon = 3
kC5O5PO3B = -71.0
kO5PO3C3B = -106.0
kPO3C3C4B = -160.67
kO3C3C4C5B = 125.44
kC3C4C5O5B = 55.65
kC4C5O5PB = 169.0

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

kPuB = 59.0
kPyB = 61.0

kC5O5PO3A = -59.3
kO5PO3C3A = -63.1
kPO3C3C4A = -157.4
kO3C3C4C5A = 75.5
kC3C4C5O5A = 49.55
kC4C5O5PA = 169.2

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

kPuA = 13.5
kPyA = 16.5

gChain1 = "A"
gChain2 = ""

function set_p_angle_res(cres, ang, iChain) {
print format("setPangleaRes(%d, %5.2f, %s)", cres, ang, iChain)
    var pres = cres-1
    var stator = {(resno=cres) and (atomName="C1\'") and (chain=iChain)}
    var pivot = {(resno=cres) and (atomName="P") and (chain=iChain)}
    var rotor = {(resno=pres) and (atomName="C1\'") and (chain=iChain)}
    select {(resno < cres) and (chain=iChain)}
    set_angle_atoms(stator, pivot, rotor, ang)
}
                
function fix_p_res(cres, iChain) {
    var pres = cres-1
    var aP = {(resno=cres) and (atomName="P") and (chain=iChain)}
    var aO5 = {(resno=cres) and (atomName="O5\'") and (chain=iChain)}
    var aC5 = {(resno=cres) and (atomName="C5\'") and (chain=iChain)}
    var aC4 = {(resno=cres) and (atomName="C4\'") and (chain=iChain)}
    var aOP1 = {(resno=cres) and (atomName="OP1") and (chain=iChain)}
    var aOP2 = {(resno=cres) and (atomName="OP2") and (chain=iChain)}
    var aO3p = {(resno=pres) and (atomName="O3\'") and (chain=iChain)}
    var aC3p = {(resno=pres) and (atomName="C3\'") and (chain=iChain)}
    
    # If collision
    if (distance(aC3p, aC5) <= kCtolerance) {
        # Push away
        select {(resno <= @{aC5.resno}) and (chain=iChain)}
        set_distance_atoms(aC3p, aC5, kCtolerance)
    }
    
    select aO5
    var dist = distance(aO3p, aO5)
    var widen = (dist < 2.85)
    var dir = (widen ? -1 : 1)
    var first = TRUE
    while (abs(dist-2.85) > kDtolerance) {
        rotateSelected @aC4 @aC5 @dir
        var newdist = distance(aO3p, aO5)
        if (widen ? (newdist < dist) : (newdist > dist)) {
            if (first) {
                dir = -dir
                rotateSelected @aC4 @aC5 @dir
            }
            else {
                break
            }
        }
        dist=newdist
        first = FALSE
    }
    select aP
    set_distance_atoms(aO5, aP, 1.73)
    set_angle_atoms(aC5, aO5, aP, 110.1)
    set_dihedral_atoms(aC4, aC5, aO5, aP, 150.3)
    aOP2.xyz = get_tet_idx(aO3p.atomIndex, aP.atomIndex, aO5.atomIndex, 1.73)
    aOP1.xyz = get_tet_idx(aO5.atomIndex, aP.atomIndex, aO3p.atomIndex, 1.73)
    #minimize select {connected(aP) or aP}
}
                
function get_interbase_rotors(aP) {
    var rotors = array()
    var sRes = aP.resno
    var mRes = sRes-1
    var iChain = aP.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 get_chi_rotor_res(res, iChain) {
    var rotors = array()
    var aO4 = {(resno=res) and (chain=iChain) and (atomName="O4\'")}
    var aC1 = {(resno=res) and (chain=iChain) and (atomName="C1\'")}
    var isR  = ((aC1 and {purine}).size > 0)
    var N1or9 = (isR ? "N9" : "N1")
    var C6or8 = (isR ? "C8" : "C6")
    
    var aN = {(resno=res) and (chain=iChain) and (atomName=N1or9)}
    var aC = {(resno=res) and (chain=iChain) and (atomName=C6or8)}

    rotors = [aO4.atomIndex, aC1.atomIndex, aN.atomIndex, aC.atomIndex]
    return rotors
}

function gen_nt_rotors(res5, res3, iChain) {
print format("gen_nt_rotors(res5=%d, res3=%d, iChain=%s)", res5, res3, iChain)#DEBUG
    var rotors = array()
    for (var i = res5+1; i <= res3; i++) {
        var aP = {(resno=i) and (chain=iChain) and (atomName="P")}
        rotors += get_interbase_rotors(aP)        
    }
    return rotors
}

function set_res_distance(stator, mover, dist, rotors) {
print format("set_res_distance(stator=%s, mover=%s, dist=%5.2f, rotors=%s)",
stator, mover, dist, rotors)#DEBUG
    var selsave = {selected}
    var cp = mover.xyz
    select mover
    set_distance_atoms(stator, mover, dist)
    var pt = mover.xyz
    mover.xyz = cp
    select selsave
    toab_track_idx(mover.atomIndex, pt, rotors)
    toab_track_idx(mover.atomIndex, pt, rotors)
    toab_track_idx(mover.atomIndex, pt, rotors)
}

function set_distance_nt_atoms( static, mobile, desired) {
print format("set_distance_nt_atoms( static=%s, mobile=%s, desired=%5.2f)",
static, mobile, desired)#DEBUG
    var rotors = gen_nt_rotors(mobile.resno, static.resno, static.chain)
    set_res_distance(static, mobile, desired, rotors) 
}

function set_distance_nt_idx( staticIdx, mobileIdx, desired) {
    set_distance_nt_atoms({atomIndex=staticIdx}, {atomIndex=mobileIdx}, desired)
}

# Moved object must be selected, fixed object not
# as[6] = fixed[1-3] moved[4-6]
# vs[6] = [distance(as[3-4]), angle(as[2-4]),
#  dihedral(as[1-4]), angle(as[5-3], dihedral(as[6-3],
#  dihedral(as[2-5]
function move_it(as, vs) {
    
    # Distance, angle, dihedral positions atom[4] to a point
    set_distance_atoms(as[3], as[4], vs[1])
    set_angle_atoms(as[2], as[3], as[4], vs[2])
    set_dihedral_atoms(as[1], as[2], as[3], as[4], vs[3])

    # Angle and dihedral orients atom[4]'s object
    set_angle_atoms(as[3], as[4], as[5], vs[4])
    set_dihedral_atoms(as[3], as[4], as[5], as[6], vs[5])

    # Dihedral sets TBD
    set_dihedral_atoms(as[2], as[3], as[4], as[5], vs[6])
}

# Pair res i on res j moving res <= i
function pair_it_res(i, j, iChain, jChain) {

    var as = array()
    var vs = array()
    as[1] = {(resno=j) and (atomName="C4\'") and (chain=jChain)}
    as[2] = {(resno=j) and (atomName="C1\'") and (chain=jChain)}
    as[3] = connected(as[2]) and {element="N"}
    as[5] = {(resno=i) and (atomName="C1\'") and (chain=iChain)}
    as[6] = {(resno=i) and (atomName="C4\'") and (chain=iChain)}
    as[4] = connected(as[5]) and {element="N"}
    select {(resno <= i) and (chain=iChain)}

    # Set distance of iN from jN (1ana=9.00)
    vs[1] = 9.00

    # Set angle of iN from jN and jC1 (1ana=124.6)
    vs[2] = 124.6

    # Set dihedral of iN from jN and jC1 and jC4 (1ana=160.0)
    vs[3] = 160.0

    # Set angle of iC1 from iN nad jN (1ana=124.6)
    vs[4] = 124.6

    # Set dihedral of iC4 from iN and iC1 and jN (1ana=160.0)
    vs[5] = 160.0

    # Set dihedral of iN from iC1 and jN and jC1 (1ana=-5.0)
    vs[6] = -5.0

    move_it(as, vs)
    fix_p_res(i, iChain)
}

# Pair A res i on A res j Hogsteen (N6-N7)2 moving res <= i
function pair_it_h_aa(i, j, iChain, jChain) {
    var as = array()
    var vs = array()
    as[1] = {(resno=j) and (atomName="N6") and (chain=jChain)}
    as[2] = {(resno=j) and (atomName="C1\'") and (chain=jChain)}
    as[3] = {(resno=j) and (atomName="N7") and (chain=jChain)}
    as[4] = {(resno=i) and (atomName="N6") and (chain=iChain)}
    as[5] = {(resno=i) and (atomName="C1\'") and (chain=iChain)}
    as[6] = {(resno=i) and (atomName="N7") and (chain=iChain)}
    select {(resno <= i) and (chain=iChain)}

    # Set distance of iN6 from jN7 (1tna=2.92)
    vs[1] = 2.92

    # Set angle of iN6 from jN7 and jC1 (1tna=123.1)
    vs[2] = 123.1

    # Set dihedral of iN6 from jN7 and jC1 and jN6 (1tna= 154.9)
    vs[3] =  154.9

    # Set angle of iN7 from iN6 and jC1 (1ana=98.2)
    vs[4] =  98.2

    # Set dihedral of iC4 from iN and iC1 and jN (1tna=18.2)
    vs[5] = 18.2
    
    # Set dihedral of iN7 from iN6 and jN7 and jC1 (1tna=177.6)
    vs[6] = 177.6

    move_it(as, vs)
    fix_p_res(i, iChain)
}

# Stack res rMove on res rFixed in A form
function base_stack_a( rMove, rFixed, iChain, jChain, sep , ang) {
    var j = rFixed
    var i = rMove
    var as = array()
    var vs = array()
    as[1] = {(resno=j) and (atomName="O3\'") and (chain=jChain)}
    as[2] = {(resno=j) and (atomName="C5\'") and (chain=jChain)}
    var Njx = (((as[1] and {purine}).size > 0) ? "N9" : "N1")
    as[3] = {(resno=j) and (atomName=Njx) and (chain=jChain)}
    as[5] = {(resno=i) and (atomName="C5\'") and (chain=iChain)}
    as[6] = {(resno=i) and (atomName="O3\'") and (chain=iChain)}
    var Nix = (((as[5] and {purine}).size > 0) ? "N9" : "N1")
    as[4] = {(resno=i) and (atomName=Nix) and (chain=iChain)}
    select {(resno <= i) and (chain=iChain)}

    # Set distance of iNx from jNx (1tna=4.2)
    vs[1] = sep

    # Set angle Njx Nix C5i (1ana=85.7)
    vs[2] =  85.7

    # Set dihedral Njx Nix C5i O3i (1tna=179.9)
    vs[3] = 179.9
    
    # Set angle C5j Njx Nxif (1tna=112.9)
    vs[4] = 112.9

    # Set dihedral O5j C5j Njx Nix (1tna= -20)
    vs[5] =  -20

    # Set dihedral of C5j Njx Nix C5i (1tna=23.7)
    vs[6] = ang

    move_it(as, vs)
    fix_p_res(i, iChain)
}

function rotate_selected_cd_atoms(a1, a2, ang) {
print format("rotate_selected_cd_atoms(a1=%s, a2=%s, ang=%5.2f", a1, a2, ang)
    #rotateSelected a1 a2 @ang#TBD
    var rang = abs(ang)
    var iang = 1.0
    var dir = ((ang < 0) ? -iang : iang)
    while (rang > 0) {
        rotateSelected @a1 @a2 @dir
        if (is_collision_in_select()) {
        #ca = count_collision_in_select(TRUE)
        #if (ca.size > 0) {
            rotateSelected @a1 @a2 @{-dir}
            break
        }
        rang -= iang
    }
}

# Rotate rotor set to move target atom to its proper place
function toab_track_idx(targetIdx, targetPt, iRotors) {
    var pt = targetPt
    var rotors = iRotors
    var targetNo = {atomIndex=targetIdx}.atomno
    var targetRes = {atomIndex=targetIdx}.resno
    var iChain = {atomIndex=targetIdx}.chain
    gOK = FALSE
    var dist = distance(pt, {atomIndex=targetIdx}.xyz)

    # For idx 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)

            # Rotate to minimize vector ====================
            #select {resno <= targetRes} or connected(i2)
            select (((atomno >= targetNo) and (chain = iChain) and
                (atomno <= @{{atomIndex=i3}.atomno}))
                or connected({atomIndex=i2}))            
            #rotateSelected {atomIndex=i3} {atomIndex=i2} @dt
            rotate_selected_cd_atoms({atomIndex=i3}, {atomIndex=i2}, -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 is_form_a( iResno, iChain) {
    var aO4 = {(resno=iResno) and (chain=iChain) and (atomName="O4\'")}
    var aC1 = {(resno=iResno) and (chain=iChain) and (atomName="C1\'")}
    var aC2 = {(resno=iResno) and (chain=iChain) and (atomName="C2\'")}
    var aC3 = {(resno=iResno) and (chain=iChain) and (atomName="C3\'")}
    return (angle(aO4, aC1, aC2, aC3) < 0.0)
}

function repair_p_res(res, iChain) {
    var aP = {((resno=res) and (chain=iChain) and (atomName="P"))}
    minimize select {connected(aP) or aP}
}

function pivot_180(res5, res3, iChain) {
    var aO5 = {(resno=res5) and (atomName="O5\'") and (chain=iChain)}
    var bO3 = {(resno=@{res3-1}) and (atomName="O3\'") and (chain=iChain)}
    var aP = {((resno=res5) and (chain=iChain) and (atomName="P"))}
    select {(resno>=res5) and (resno<res3) and not (connected(aP) or aP)}
    fix not selected
    rotateSelected @bO3 @aO5 180.0
    fix none
    fix_p_res(res5, iChain)
}

function res_to_ab(iRes, iChain, toA) {
    var i = iRes
    var aO3 = {(resno=i) and (chain=iChain) and (atomName="O3\'")}
    var aC3 = {(resno=i) and (chain=iChain) and (atomName="C3\'")}
    var aC4 = {(resno=i) and (chain=iChain) and (atomName="C4\'")}
    var aC5 = {(resno=i) and (chain=iChain) and (atomName="C5\'")}

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

    select {resno <= i} and not aO3 and not aC3
    set_dihedral_atoms(aO3, aC3, aC4, aC5, (toA ? kO3C3C4C5A : kO3C3C4C5B))
    
    # Set chi
    var aNx = -1
    var aCx = -1
    select {(resno=i) and base}
    if ((aC1 and {purine}).size > 0) {
        aNx =  {(resno=i) and (chain=iChain) and (atomName="N9")}
        aCx =  {(resno=i) and (chain=iChain) and (atomName="C8")}
        ang = (toA ? kPuA : kPuB)
        pang = (toA ? kPyA : kPyB)
    }
    else {
        aNx =  {(resno=i) and (chain=iChain) and (atomName="N1")}
        aCx =  {(resno=i) and (chain=iChain) and (atomName="C6")}
        ang = (toA ? kPyA : kPyB)
        pang = (toA ? kPuA : kPuB)
    }
    set_dihedral_atoms(aO4, aC1, aNx, aCx, ang)
    
    # Set pucker 3' endo or 2' endo
    pSet = {aC1 or aC2 or aO2}
    select pSet or {(resno=i) and (chain=iChain) and base}
    set_dihedral_atoms(aC5, aC4, aO4, aC1, (toA ? kC5C4O4C1A : kC5C4O4C1B))
    set_dihedral_atoms(aC4, aO4, aC1, aNx, (toA ? kC4O4C1NxA : kC4O4C1NxB))
    set_dihedral_atoms(aC4, aO4, aC1, aC2, (toA ? kC4O4C1C2A : kC4O4C1C2B))
    if (aO2.size > 0) {
        select aO2 or aC2
        ang = (toA ? kC3C1C2O2A : kC3C1C2O2B)
        set_dihedral_atoms(aC3, aC1, aC2, aO2, (toA ? kC3C1C2O2A : kC3C1C2O2B))
    }
    set_distance_atoms(aC3, aC2, 1.52)
    set_distance_atoms(aC1, aC2, 1.52)
}


function who_pairs(iRes, iChain) {
    var aC4or6 =  {(resno=iRes) and (chain=iChain) and (atomName="C4")}
    var aN1or3 =  {(resno=iRes) and (chain=iChain) and (atomName="N1")}
    if ({aN1or3 and purine}.size = 0) {
        aC4or6 =  {(resno=iRes) and (chain=iChain) and (atomName="C6")}
        aN1or3 =  {(resno=iRes) and (chain=iChain) and (atomName="N3")}
    }
    var near = within(3.1, aN1or3) and {resno!=iRes} and {element="N"}
    for (var i = 1; i <= near.size; i++) {
        if (angle(near[i], aN1or3, aC4or6) > 150) {
            return [near[i].resno, near[i].chain]
        } 
    }
    return [-1, aC4or6.chain]
}
  
# end of plicoNTcommon.spt

Contributors

Remig