User:Remig/plico/plicoCommonNT
< User:Remig | plico
Jump to navigation
Jump to search
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.2 beta 6/9/2014 -add stuff
#
# 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, aC1.chain]
}
# end of plicoNTcommon.spt
