Difference between revisions of "User:Remig/plico/plicoCommonNT"
< User:Remig | plico
Jump to navigation
Jump to search
(Add support for fold_trna.spt) |
(Support fold_trna.spt) |
||
| Line 3: | Line 3: | ||
Copy and paste the following into a text editor and save in your scripts folder as plicoNTcommon.spt. | Copy and paste the following into a text editor and save in your scripts folder as plicoNTcommon.spt. | ||
<pre># plicoNTcommon - Jmol script by Ron Mignery | <pre># plicoNTcommon - Jmol script by Ron Mignery | ||
| − | # v1.8 beta 7/13/2015 - | + | # v1.8 beta 7/13/2015 -support fold_trna.spt |
# | # | ||
# Routines and values common to Plico suite scripts that work with nucleotides | # Routines and values common to Plico suite scripts that work with nucleotides | ||
| Line 85: | Line 85: | ||
aOP2.xyz = get_tet_idx(aO3p.atomIndex, aP.atomIndex, aO5.atomIndex, 1.73) | 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) | aOP1.xyz = get_tet_idx(aO5.atomIndex, aP.atomIndex, aO3p.atomIndex, 1.73) | ||
| − | + | plico_minimize( {(connected(aP) or aP) and not aO3p}) | |
| − | |||
| − | |||
| − | |||
select selsave | select selsave | ||
} | } | ||
| Line 177: | Line 174: | ||
aOP1.xyz = get_tet_idx(aO5.atomIndex, aP.atomIndex, aO3p.atomIndex, 1.73) | aOP1.xyz = get_tet_idx(aO5.atomIndex, aP.atomIndex, aO3p.atomIndex, 1.73) | ||
if (force) { | if (force) { | ||
| − | + | plico_minimize ({(connected(aP) or aP) and not aO3p}) | |
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
} | } | ||
} # endfor i | } # endfor i | ||
| Line 198: | Line 189: | ||
function fix_all_nt_collisions( iChain) { | function fix_all_nt_collisions( iChain) { | ||
var selsave = {selected} | var selsave = {selected} | ||
| − | |||
| − | |||
| − | |||
| − | |||
chset = count_collisions(true) | chset = count_collisions(true) | ||
for (var i = 1; i <= chset.size; i++) { | for (var i = 1; i <= chset.size; i++) { | ||
| Line 214: | Line 201: | ||
if ({c and cset and not base}.size > 0) { | if ({c and cset and not base}.size > 0) { | ||
select {@{rset} and base} | select {@{rset} and base} | ||
| − | + | plico_minimize( {selected}) | |
} | } | ||
else if (c.atomname[1][2] == "OP") { | else if (c.atomname[1][2] == "OP") { | ||
| Line 220: | Line 207: | ||
} | } | ||
else { | else { | ||
| − | + | plico_minimize(rset) | |
| − | |||
} | } | ||
} | } | ||
measure off | measure off | ||
| − | |||
| − | |||
select selsave | select selsave | ||
} | } | ||
| Line 489: | Line 473: | ||
select {((resno=res5) or (resno=@{res5+1})) | select {((resno=res5) or (resno=@{res5+1})) | ||
and (chain=iChain) and (file=f) and (model=m)} | and (chain=iChain) and (file=f) and (model=m)} | ||
| − | + | plico_minimize( {selected}) | |
fix_p_res(res5+1, iChain, true) | fix_p_res(res5+1, iChain, true) | ||
} | } | ||
| Line 522: | Line 506: | ||
select {((resno=res5) or (resno=@{res5+1})) | select {((resno=res5) or (resno=@{res5+1})) | ||
and (chain=iChain) and (file=f) and (model=m)} | and (chain=iChain) and (file=f) and (model=m)} | ||
| − | + | plico_minimize( {selected}) | |
fix_p_res(res5+1, iChain, true) | fix_p_res(res5+1, iChain, true) | ||
fix_p_res(res5, iChain, true) | fix_p_res(res5, iChain, true) | ||
| Line 895: | Line 879: | ||
function adjust_nts(res5, res3, iChain, toab, a, s) { | function adjust_nts(res5, res3, iChain, toab, a, s) { | ||
| − | |||
| − | |||
# Collect any pairing | # Collect any pairing | ||
| Line 934: | Line 916: | ||
} | } | ||
} | } | ||
| − | |||
} | } | ||
| Line 957: | Line 938: | ||
function repair_p_res(res, iChain) { | function repair_p_res(res, iChain) { | ||
var aP = get_atom_rcn( res, iChain, "P") | var aP = get_atom_rcn( res, iChain, "P") | ||
| − | + | plico_minimize( {connected(aP) or aP}) | |
| − | |||
| − | |||
| − | |||
} | } | ||
| Line 1,002: | Line 980: | ||
var pyang = abs(angle(aC6or4, near[i], aN1or3)) | var pyang = abs(angle(aC6or4, near[i], aN1or3)) | ||
var dihedral = abs(angle(aC2p, near[i], aN1or3, aC2)) | var dihedral = abs(angle(aC2p, near[i], aN1or3, aC2)) | ||
| − | |||
if ((puang > 110) and (pyang > 110) and (dihedral < 40)) { | if ((puang > 110) and (pyang > 110) and (dihedral < 40)) { | ||
return [near[i].resno, near[i].chain, dist, puang] | return [near[i].resno, near[i].chain, dist, puang] | ||
| Line 1,493: | Line 1,470: | ||
connected({resno=r}) and {(chain=iChain) and (file=f) and (model=m)}) | connected({resno=r}) and {(chain=iChain) and (file=f) and (model=m)}) | ||
if (cset.size > 0) { | if (cset.size > 0) { | ||
| − | |||
| − | |||
| − | |||
| − | |||
for (var i = 0; i < cset.size; i++) { | for (var i = 0; i < cset.size; i++) { | ||
| − | + | plico_minimize( {(resno=r) or (resno=@{cset[i].resno})}) | |
} | } | ||
| − | |||
| − | |||
} | } | ||
| Line 1,509: | Line 1,480: | ||
and {(chain=iChain) and (file=f) and (model=m)}) | and {(chain=iChain) and (file=f) and (model=m)}) | ||
if (cset.size > 0) { | if (cset.size > 0) { | ||
| − | |||
| − | |||
| − | |||
| − | |||
for (var i = 0; i < cset.size; i++) { | for (var i = 0; i < cset.size; i++) { | ||
| − | + | plico_minimize( {(resno=r) or (resno=@{cset[i].resno})}) | |
} | } | ||
| − | |||
| − | |||
} | } | ||
} | } | ||
| Line 1,579: | Line 1,544: | ||
function replane_bases(iChain) { | function replane_bases(iChain) { | ||
| − | |||
| − | |||
for (var i = {chain=iChain}.resno.min; i <= {chain=iChain}.resno.max; i++) { | for (var i = {chain=iChain}.resno.min; i <= {chain=iChain}.resno.max; i++) { | ||
| − | + | plico_minimize( {resno=i and base}) | |
} | } | ||
| − | |||
} | } | ||
Revision as of 18:27, 24 July 2015
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.8 beta 7/13/2015 -support fold_trna.spt
#
# 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 them
kNTcommon = 5
kC5O5PO3B = -71.0
kO5PO3C3B = -107.0
kPO3C3C4B = -161.5
kO3C3C4C5B = 140.0
kC3C4C5O5B = 55.65
kC4C5O5PB = 169.0
kO4C4C3C2B = 15.92
kC4O4C1C2B = -19.9 #-41.7 1bna minimized
kC2O4C1NxB = -122.6 #-159.0 1bna minimized
kC5C4O4C1B = 122.2 #146.3 1bna minimized
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
kC2O4C1NxA = -131.0
kC5C4O4C1A = 144.85
kC3C1C2O2A = 116.3
kPuA = 13.5
kPyA = 16.5
gChain1 = "A"
gChain2 = ""
# Select before calling
function force_p_res(cres, iChain) {
var pres = cres-1
var aP = get_atom_rcn( cres, iChain, "P")
var aO5 = get_atom_rcn( cres, iChain, "O5\'")
var aC5 = get_atom_rcn( cres, iChain, "C5\'")
var aC4 = get_atom_rcn( cres, iChain, "C4\'")
var aOP1 = get_atom_rcn( cres, iChain, "OP1")
var aOP2 = get_atom_rcn( cres, iChain, "OP2")
var aO3p = get_atom_rcn( pres, iChain, "O3\'")
var aC3p = get_atom_rcn( pres, iChain, "C3\'")
if (aO3p.size > 0) {
var selsave = {selected}
set_distance_atoms(aP3p, aC5, 3.1)
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)
plico_minimize( {(connected(aP) or aP) and not aO3p})
select selsave
}
}
function fix_p_res(cres, iChain, force) {
var f = (_frameID/1000000)
var m = (_frameID%1000000)
var pres = cres-1
var aP = get_atom_rcn( cres, iChain, "P")
var aO5 = get_atom_rcn( cres, iChain, "O5\'")
var aC5 = get_atom_rcn( cres, iChain, "C5\'")
var aC4 = get_atom_rcn( cres, iChain, "C4\'")
var aC1 = get_atom_rcn( cres, iChain, "C1\'")
var aOP1 = get_atom_rcn( cres, iChain, "OP1")
var aOP2 = get_atom_rcn( cres, iChain, "OP2")
var aO3p = get_atom_rcn( pres, iChain, "O3\'")
var aC3p = get_atom_rcn( pres, iChain, "C3\'")
if ((aO3p.size > 0)and (aC4.size > 0)) {
var selsave = {selected}
select {aP or aOP1 or aOP2}
set_angle_atoms(aC3p, aO3p, aP, 109.0)
# Twice through
for (var i = 0; i < 2; i++) {
# If collision
if (force and distance(aC3p, aC5) <= kCtolerance) {
# Push away
select {(resno <= @{aC5.resno}) and (chain=iChain)
and (file=f) and (model=m)}
set_distance_atoms(aC3p, aC5, kCtolerance)
}
select aO5
set_dihedral_atoms(aC1, aC4, aC5, aO5, -40)
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
}
if (force and (abs(distance(aO3p, aC5)-4.11) < kDtolerance)) {
# Push away
select {(resno <= @{aC5.resno}) and (chain=iChain)
and (file=f) and (model=m)}
set_distance_atoms(aO3p, aC5, 4.11)
}
select aP
set_distance_atoms(aO5, aP, 1.73)
set_angle_atoms(aC5, aO5, aP, 110.1)
dist = distance(aO3p, aP)
widen = (dist < 1.73)
dir = (widen ? -1 : 1)
first = true
while (abs(dist-1.73) > kDtolerance) {
rotateSelected @aC5 @aO5 @dir
var newdist = distance(aO3p, aP)
if (widen ? (newdist < dist) : (newdist > dist)) {
if (first) {
dir = -dir
rotateSelected @aC5 @aO5 @dir
}
else {
break
}
}
dist=newdist
first = false
}
#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)
if (force) {
plico_minimize ({(connected(aP) or aP) and not aO3p})
}
} # endfor i
select selsave
}
}
function fix_p_res_range(res5, res3, iChain, force) {
for (var i = res5; i <= res3; i++) {
fix_p_res(i, iChain, force)
}
}
function fix_all_nt_collisions( iChain) {
var selsave = {selected}
chset = count_collisions(true)
for (var i = 1; i <= chset.size; i++) {
c = chset[i]
cset = (within(kCtolerance, c) and not c
and not connected(chset[i]))
rset = [{resno=@{c.resno}}]
for (var j = 1; j <= cset.size; j++) {
rset += {resno=@{cset[j].resno}}
}
if ({c and cset and not base}.size > 0) {
select {@{rset} and base}
plico_minimize( {selected})
}
else if (c.atomname[1][2] == "OP") {
fix_p_res(chset[i].resno, iChain, true)
}
else {
plico_minimize(rset)
}
}
measure off
select selsave
}
# The following functions position one nt relative to another:
# Common positioning functions:
function get_rotors_res(res) {
var f = (_frameID/1000000)
var m = (_frameID%1000000)
var rotors = array()
var sRes = res
var mRes = sRes-1
var iChain = {(resno=res) and (atomName="P")
and (file=f) and (model=m)}.chain
var mC4 = get_atom_rcn( mRes, iChain, "C4\'")
var mC3 = get_atom_rcn( mRes, iChain, "C3\'")
var mO3 = get_atom_rcn( mRes, iChain, "O3\'")
var sP = get_atom_rcn( sRes, iChain, "P" )
var sO5 = get_atom_rcn( sRes, iChain, "O5\'")
var sC5 = get_atom_rcn( sRes, iChain, "C5\'")
var sC4 = get_atom_rcn( sRes, iChain, "C4\'")
var sC3 = get_atom_rcn( sRes, iChain, "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_nt_chi_rotor_res(res, iChain) {
var rotors = array()
var aO4 = get_atom_rcn( res, iChain, "O4\'")
var aC1 = get_atom_rcn( res, iChain, "C1\'")
var isR = ((aC1 and {purine}).size > 0)
var N1or9 = (isR ? "N9" : "N1")
var C6or8 = (isR ? "C8" : "C6")
var aN = get_atom_rcn(res, iChain, N1or9)
var aC = get_atom_rcn(res, iChain, C6or8)
rotors = [aO4.atomIndex, aC1.atomIndex, aN.atomIndex, aC.atomIndex]
return rotors
}
function get_nt_ab_rotor_res(res, iChain) {
var rotors = array()
var aC5 = get_atom_rcn(res, iChain, "C5\'")
var aC4 = get_atom_rcn(res, iChain, "C4\'")
var aC3 = get_atom_rcn(res, iChain, "C3\'")
var aO3 = get_atom_rcn(res, iChain, "O3\'")
rotors = [aO3.atomIndex, aC3.atomIndex, aC4.atomIndex, aC5.atomIndex]
return rotors
}
function gen_nt_rotors(res5, res3, iChain) {
var rotors = array()
for (var i = res5+1; i <= res3; i++) {
rotors += get_rotors_res(i, iChain)
}
return rotors
}
# ri=moved rj=fixed
function position_nt_by_vs(ri, rj, ares, vs, iChain, jChain) {
var f = (_frameID/1000000)
var m = (_frameID%1000000)
if (ri > rj) {
var as = gen_as(ri, rj, iChain, jChain)
select {(resno < ares) and (resno >= ri)
and (chain=iChain) and (file=f) and (model=m)}
}
else {
var as = gen_as(rj, ri, jChain, iChain)
select {(resno > ares) and (resno <= ri)
and (chain=iChain) and (file=f) and (model=m)}
}
move_it(as, vs)
}
# Moved object must be selected, fixed object not
# as[6] = fixed[1-3] moved[4-6] chis [7-8]
# 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])
# If chis
if (vs.size > 6) {
var sel = {selected}
select {(resno=@{as[4].resno}) and base}
set_dihedral_atoms(as[6], as[5], as[4], as[8], vs[8])
}
}
# ri=moved rj=fixed
function gen_as(ri, rj, iChain, jChain) {
var as = array()
as[1] = get_atom_rcn(rj, jChain, "C4\'")
as[2] = get_atom_rcn(rj, jChain, "C1\'")
as[3] = connected(as[2]) and {element="N"}
as[5] = get_atom_rcn(ri, iChain, "C1\'")
as[6] = get_atom_rcn(ri, iChain, "C4\'")
as[4] = connected(as[5]) and {element="N"}
as[7] = get_atom_rcn(rj, jChain,
(((as[1] and {purine}).size > 0) ? "C8" : "C6"))
as[8] = get_atom_rcn(ri, iChain,
(((as[6] and {purine}).size > 0) ? "C8" : "C6"))
return as
}
# Specific positioning functions:
# Pair res5 on res3 moving res5 <= res3
function pair_it_res(res5, res3, ares, iChain, jChain) {
var f = (_frameID/1000000)
var m = (_frameID%1000000)
var as = gen_as(res5, res3, iChain, jChain)
var isA = is_form_a(res5, iChain)
var vs = array()
vs[1] = 8.83 # distance res5 N9or1 and res3 N9or1
vs[2] = 126.12 # angle res5 N9or1 and res3 N9or1 C1
vs[3] = (isA ? 160.0 : -134.97) # dihedral res5 N9or1 and res3 N9or1 C1 C4
vs[4] = 125.32 # angle res5 C1 N9or1 and res3 N9or1
vs[5] = (isA ? 160.0 : -141.46) # dihedral res5 C4 N9or1 C1 and res3 N9or1
vs[6] = (isA ? -5.0 : -17.87) # dihedral res5 N9or1 C1 and res3 N9or1 C1
vs[7] = (isA ? -20.0 : 38) # dihedral chi res3
vs[8] = (isA ? -20.0 : 38) # dihedral chi res5
if (ares < 0) {
select ((resno=res5) and (chain=iChain) and (file=f) and (model=m))
}
else if (ares > 0) {
select ((resno <= ares) and (chain=iChain) and (file=f) and (model=m))
}
move_it(as, vs)
fix_p_res(res5, iChain, true)
fix_p_res(res5+1, iChain, true)
}
# Flatstack res5 on res3 moving just res5
function single_flatstack_res5_on_res3(res5, res3, iChain, jChain) {
var f = (_frameID/1000000)
var m = (_frameID%1000000)
var as = gen_as(res5, res3, iChain, jChain)
var vs = array()
vs[1] = 7.00 # distance res5 N9or1 and res3 N9or1
vs[2] = 89.1 # angle res5 N9or1 and res3 N9or1 C1
vs[3] = -49.9 # dihedral res5 N9or1 and res3 N9or1 C1 C4
vs[4] = 83.4 # angle res5 C1 N9or1 and res3 N9or1
vs[5] = 125.7 # dihedral res5 C4 N9or1 C1 and res3 N9or1
vs[6] = 5.8 # dihedral res5 N9or1 C1 and res3 N9or1 C1
select {(resno=res5) and (chain=iChain) and (file=f) and (model=m)}
move_it(as, vs)
#force_p_res(res3, jChain)
#move_it(as, vs)
#fix_p_res(res3, jChain, true)
#force_p_res(res3, jChain)
}
# Outstack res5 on res3 moving just res5
function single_outstack_res5_on_res3(res5, res3, iChain, jChain) {
var f = (_frameID/1000000)
var m = (_frameID%1000000)
var as = gen_as(res5, res3, iChain, jChain)
var vs = array()
vs[1] = 8.23 # distance res5 N9or1 and res3 N9or1
vs[2] = 32.4 # angle res5 N9or1 and res3 N9or1 C1
vs[3] = -26.8 # dihedral res5 N9or1 and res3 N9or1 C1 C4
vs[4] = 99.6 # angle res5 C1 N9or1 and res3 N9or1
vs[5] = 57.4 # dihedral res5 C4 N9or1 C1 and res3 N9or1
vs[6] = 179.1 # dihedral res5 N9or1 C1 and res3 N9or1 C1
select {(resno=res5) and (chain=iChain) and (file=f) and (model=m)}
move_it(as, vs)
force_p_res(res3, jChain)
move_it(as, vs)
##fix_p_res(res3, jChain, true)
}
# Flatstack res5 on res3 moving just res5
function single_flatstack_res5_on_res3(res3, res5, iChain, jChain) {
var f = (_frameID/1000000)
var m = (_frameID%1000000)
var as = gen_as(res5, res3, iChain, jChain)
vs = array()
vs[1] = 6.00 # distance res5 N9or1 and res3 N9or1
vs[2] = 90#75.1 # angle res5 N9or1 and res3 N9or1 C1
vs[3] = 90#135.3 # dihedral res5 N9or1 and res3 N9or1 C1 C4
vs[4] = 90#89.9 # angle res5 C1 N9or1 and res3 N9or1
vs[5] = -90#-47.3 # dihedral res5 C4 N9or1 C1 and res3 N9or1
vs[6] = 0#1.7 # dihedral res5 N9or1 C1 and res3 N9or1 C1
select {(resno=res5) and (chain=iChain) and (file=f) and (model=m)}
move_it(as, vs)
force_p_res(res5, jChain)
}
# Outstack res3 on res5 moving just res5
function single_outstack_res3_on_res5(res5, res3, iChain, jChain) {
var f = (_frameID/1000000)
var m = (_frameID%1000000)
var as = gen_as(res5, res3, iChain, jChain)
var vs = array()
vs[1] = 8.9 # distance res5 N9or1 and res3 N9or1
vs[2] = 65.3 # angle res5 N9or1 and res3 N9or1 C1
vs[3] = 55.7 # dihedral res5 N9or1 and res3 N9or1 C1 C4
vs[4] = 61.2 # angle res5 C1 N9or1 and res3 N9or1
vs[5] = -41.2 # dihedral res5 C4 N9or1 C1 and res3 N9or1
vs[6] = -138.4 # dihedral res5 N9or1 C1 and res3 N9or1 C1
select {(resno=res5) and (chain=iChain) and (file=f) and (model=m)}
move_it(as, vs)
force_p_res(res5, jChain)
move_it(as, vs)
##fix_p_res(res5, jChain, true)
}
function make_major_groove_triplex(res5, res3, iChain, jChain) {
var f = (_frameID/1000000)
var m = (_frameID%1000000)
var as = gen_as(res5, res3, iChain, jChain)
var vs = array()
vs[1] = 8.11 # distance res5 N9or1 and res3 N9or1
vs[2] = 166.2 # angle res5 N9or1 and res3 N9or1 C1
vs[3] = 0.3 # dihedral res5 N9or1 and res3 N9or1 C1 C4
vs[4] = 162.8 # angle res5 C1 N9or1 and res3 N9or1
vs[5] = 114.6 # dihedral res5 C4 N9or1 C1 and res3 N9or1
vs[6] = 138.1 # dihedral res5 N9or1 C1 and res3 N9or1 C1
# Move the nt into final position
select {(resno <= res5) and (chain=iChain) and (file=f) and (model=m)}
move_it(as, vs)
# Rotate ribose to hbond O2' to res5+1 N7
select {(resno = res5) and not base
and (chain=iChain) and (file=f) and (model=m)}
var aC1 = get_atom_rcn( res5, iChain, "C1\'")
var aN9 = get_atom_rcn( res5, iChain, "N9")
rotate selected @aC1 @aN9 -40.0
# Fix up
select {((resno=res5) or (resno=@{res5+1}))
and (chain=iChain) and (file=f) and (model=m)}
plico_minimize( {selected})
fix_p_res(res5+1, iChain, true)
}
# Pair U res5 on A res3 Hoogsteen N3-N7, N6-O2, O4-O1p moving res5 => res3
function make_hoogsteen_pair_yr(res5, res3, iChain, jChain) {
var f = (_frameID/1000000)
var m = (_frameID%1000000)
var as = gen_as(res5, res3, iChain, jChain)
var vs = array()
var cp = as[6].xyz
vs[1] = 7.05 # distance res5 N9or1 and res3 N9or1
vs[2] = 150.7 # angle res5 N9or1 and res3 N9or1 C1
vs[3] = -33.1 # dihedral res5 N9or1 and res3 N9or1 C1 C4
vs[4] = 143.0 # angle res5 C1 N9or1 and res3 N9or1
vs[5] = -173.9 # dihedral res5 C4 N9or1 C1 and res3 N9or1
vs[6] = -179.8 # dihedral res5 N9or1 C1 and res3 N9or1 C1
# Move the nt into final position
select {(resno <= res5) and (chain=iChain) and (file=f) and (model=m)}
move_it(as, vs)
# Rotate 5 end out of the way
select {(resno < res5) and (chain=iChain) and (file=f) and (model=m)}
var aC4 = get_atom_rcn( res5, iChain, "C4\'")
var aC5 = get_atom_rcn( res5, iChain, "C5\'")
rotate selected @aC4 @aC5 160.0
# Fix up
select {((resno=res5) or (resno=@{res5+1}))
and (chain=iChain) and (file=f) and (model=m)}
plico_minimize( {selected})
fix_p_res(res5+1, iChain, true)
fix_p_res(res5, iChain, true)
}
function level_base(rMove, rFixed, iChain, jChain) {
var selsave = {selected}
var f = (_frameID/1000000)
var m = (_frameID%1000000)
var mC1 = get_atom_rcn(rMove, iChain, "C1\'")
var mIsR = ({mC1 and purine}.size > 0)
var m9or1 = (mIsR ? "N9" : "N1")
var m6or8 = (mIsR ? "C8" : "C6")
var m4or2 = (mIsR ? "C4" : "C2")
var mN = get_atom_rcn(rMove, iChain, m9or1)
var mC6or8 = get_atom_rcn(rMove, iChain, m6or8)
var mC4or2 = get_atom_rcn(rMove, iChain, m4or2)
var fC1 = get_atom_rcn(rFixed, jChain, "C1\'")
var fIsR = ({fC1 and purine}.size > 0)
var f9or1 = (fIsR ? "N9" : "N1")
var f6or8 = (fIsR ? "C8" : "C6")
var f4or2 = (fIsR ? "C4" : "C2")
var fN = get_atom_rcn(rFixed, jChain, f9or1)
var fC6or8 = get_atom_rcn(rFixed, jChain, f6or8)
var fC4or2 = get_atom_rcn(rFixed, jChain, f4or2)
var dist = abs(distance(fc4or2, mc4or2) - distance(fc6or8, mc6or8))
var newdist = dist
var dir = 0.1
select {(resno=rMove) and (chain=iChain) and base and (file=f) and (model=m)}
while(newdist > 0.01) {
if (newdist > dist) {
if (dir == 0.1) {
dir = -0.1
}
else {
rotateSelected @mC1 @mN @{-dir}
break
}
}
dist = newdist
rotateSelected @mC1 @mN @{dir}
newdist = abs(distance(fc4or2, mc4or2) - distance(fc6or8, mc6or8))
}
select selsave
}
# Stack res rMove on res rFixed
function base_stack_res( rMove, rFixed, iChain, jChain, sep , ang, single) {
var isA = is_form_a(rMove, iChain)
var is3on5 = (rMove > rFixed)
var j = rFixed
var i = rMove
var as = array()
var vs = array()
as = gen_as(rMove, rFixed, iChain, jChain)
var f = as[1].file
var m = as[1].model
if (single) {
select {(resno = i) and (chain=iChain) and (file=f) and (model=m)}
}
else {
if (is3on5) {
select {(resno >= i) and (chain=iChain) and (file=f) and (model=m)}
}
else {
select {(resno <= i) and (chain=iChain) and (file=f) and (model=m)}
}
}
# Set distance of fres N1or9 from mres N1or9 (1tna=4.2)
vs[1] = sep
# Set angle fres C1' N1or9 and mres N1or9 (A=6tna B=1ana)
vs[2] = (isA ? (is3on5 ? 92.8 : 113.9) : (is3on5 ? 83.3 : 115.23))#78.3 : 110.23
# Set dihedral fres C4' C1' N1or9 and mres N1or9 (A=6tna B=1ana)
vs[3] = (isA ? (is3on5 ? 110.0 : -71.2) : (is3on5 ? 165.92 : -28.31))
# Set angle fres N1or9 and mres N1or9 C1' (A=6tna B=1ana)
vs[4] = (isA ? (is3on5 ? 113.9 : 92.8) : (is3on5 ? 115.23 : 83.3))
# Set dihedral fres N1or9 and mres N1or9 C1' C4' (A=6tna B=1ana)
vs[5] = (isA ? (is3on5 ? -71.2 : 110.0) : (is3on5 ? -28.31 : 165.92))
# Set dihedral of fres C5 N1or9 and mres N1or9 C5 (1tna=20)
vs[6] = ang
move_it(as, vs)
#select {(resno=rMove) or (resno=rFixed)}
#minimize {selected}
#force_p_res(i, iChain)
##fix_p_res(i+1, iChain, true)
}
# Rotate rotor set to move target atom to its proper place
# ares is the 5ward res limit exclusive of the mobile
function move_atom_nt(targetIdx, targetPt, ares, rotors, force) {
var set3
var f = (_frameID/1000000)
var m = (_frameID%1000000)
var pt = targetPt
var targetNo = {atomIndex=targetIdx}.atomno
var targetRes = {atomIndex=targetIdx}.resno
var iChain = {atomIndex=targetIdx}.chain
gOK = false
var dist = distance(pt, {atomIndex=targetIdx}.xyz)
# If target is a C1' atom, collect its base
var tBase = ({})
var i1 = 0
var i2 = 0
var i3 = 0
var i4 = 0
if ({atomIndex=targetIdx}.atomName == "C1\'") {
tBase = {(resno = targetRes) and base}
}
# 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) {
i2 = rotors[ri+1]
i3 = rotors[ri+2]
i4 = rotors[ri+3]
if ((i2 != targetIdx) and (i3 != targetIdx) and (i4 != targetIdx)) {
if ({blocked and {atomIndex=i2}}.count == 0) {
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
}
i1 = rotors[imax+0]
i2 = rotors[imax+1]
i3 = rotors[imax+2]
i4 = rotors[imax+3]
# Get dihedral of rotor with target point
var dt = (angle({atomIndex=targetIdx}, {atomIndex=i2},
{atomIndex=i3}, pt)/(rotors.size/20))
# Select and rotate
if (ares > targetRes) {
select_3ward_atom({atomIndex=i3}, ares, iChain)
res3 = {atomIndex=i4}.resno
}
else {
select_5ward_atom({atomIndex=i3}, ares, iChain)
res3 = {atomIndex=i1}.resno
}
select remove tbase
#***************************************************
rotateSelected {atomIndex=i3} {atomIndex=i2} @{-dt}
# If collisions
set3 = (within(kCtolerance, {selected}) and not {selected}
and not connected({selected}))
if ((force == false) and (set3.size > 0)) {
# Binary undo until fixed
while ((abs(dt) > kDtolerance)
and ((set5 and within(kCtolerance, set3)).size > 0)) {
dt /= 2.0
rotateSelected {atomIndex=i3} {atomIndex=i2} @{dt}
}
while ((abs(dt) > kDtolerance)
and ((set5 and within(kCtolerance, set3)).size > 0)) {
dt /= 2.0
rotateSelected {atomIndex=i3} {atomIndex=i2} @{-dt}
}
rotateSelected {atomIndex=i3} {atomIndex=i2} @{dt}
}
# If close enough, stop
dist = distance(pt, {atomIndex=targetIdx})
if (dist < kDtolerance) {
gOK = true
gTargetPt = pt
break
}
# Block rotor
blocked |= {atomIndex=i2}
} # endfor num rotors passes
if (gOK) {
break
}
} # endfor 20 passes
return set3
}
# Counter rotate rotor set to move target atom to its proper place
function move_atom_by_cr_nt(targetIdx, targetPt, ares, iRotors) {
var f = (_frameID/1000000)
var m = (_frameID%1000000)
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)
# If target is a C1' atom, collect its base
var tBase = ({})
var i1 = 0
var i2 = 0
var i3 = 0
var i4 = 0
if ({atomIndex=targetIdx}.atomName == "C1\'") {
tBase = {(resno = targetRes) and base}
}
# For all C4'-C5' axes
for (var ri = 1; ri < rotors.size; ri += 4) {
if ({atomIndex=@{rotors[ri]}}.atomName == "C4\'") {
# While distance lessens
var dist = distance(pt, {atomIndex=targetIdx})
var first = true
var dt = 5.0
while (dist > kDtolerance) {
# Counter rotate C4'-C5' and O5'-P axes
var i1 = rotors[ri+8]
var i2 = rotors[ri+9]
var i3 = rotors[ri+10]
var i4 = rotors[ri+11]
var x2 = rotors[ri+17]
var x3 = rotors[ri+18]
var x4 = rotors[ri+19]
var res3 = 0
# Select and rotate
if (ares > targetRes) {
select_3ward_atom({atomIndex=i3}, ares, iChain)
res3 = {atomIndex=i4}.resno
}
else {
select_5ward_atom({atomIndex=i3}, ares, iChain)
res3 = {atomIndex=i1}.resno
}
select remove tbase
rotateSelected {atomIndex=i3} {atomIndex=i2} @{dt}
# Select and counter rotate
if (ares > targetRes) {
select_3ward_atom({atomIndex=x3}, ares, iChain)
}
else {
select_5ward_atom({atomIndex=x3}, ares, iChain)
}
select remove tbase
rotateSelected {atomIndex=x3} {atomIndex=x2} @{-dt}
# If first and worse, reverse
var newdist = distance(pt, {atomIndex=targetIdx})
if (newdist > dist) {
if (first) {
dt = -dt
}
else {
break
}
}
first = false
dist = newdist
/*# If collisions
var res5 = res3-1
var set3 = {(resno=res3) and (atomName!="P") and (atomName!="OP1")
and (file=f) and (model=m)}
var set5 = {(resno=res5) and (atomName!="P") and (atomName!="OP1")
and (file=f) and (model=m)}
if ((set5 and within(kCtolerance, set3)).size > 0) {
# Binary undo until fixed
while ((abs(dt) > kDtolerance)
and ((set5 and within(kCtolerance, set3)).size > 0)) {
dt /= 2.0
rotateSelected {atomIndex=i3} {atomIndex=i2} @{dt}
}
while ((abs(dt) > kDtolerance)
and ((set5 and within(kCtolerance, set3)).size > 0)) {
dt /= 2.0
rotateSelected {atomIndex=i3} {atomIndex=i2} @{-dt}
}
rotateSelected {atomIndex=i3} {atomIndex=i2} @{dt}
}*/
} # endwhile
}
} # endfor rotors
}
# If ares < 0 then adjust iRes only
function to_ab_nt_res(res, ares, iChain, toA) {
var selsave = {selected}
var f = (_frameID/1000000)
var m = (_frameID%1000000)
var aO3 = get_atom_rcn( res, iChain, "O3\'")
var aC3 = get_atom_rcn( res, iChain, "C3\'")
var aC4 = get_atom_rcn( res, iChain, "C4\'")
var aC5 = get_atom_rcn( res, iChain, "C5\'")
var aC1 = get_atom_rcn( res, iChain, "C1\'")
var aC2 = get_atom_rcn( res, iChain, "C2\'")
var aO2 = get_atom_rcn( res, iChain, "O2\'")
var aO4 = get_atom_rcn( res, iChain, "O4\'")
if (ares < 0) {
select ((resno=res) and (chain=iChain) and (file=f) and (model=m)
and not aO3 and not aC3 and not aC4)
}
else {
select ((resno >= ares) and (resno <= res) and (chain=iChain)
and (file=f) and (model=m) and not aO3 and not aC3 and not aC4)
}
set_dihedral_atoms(aO3, aC3, aC4, aC5, (toA ? kO3C3C4C5A : kO3C3C4C5B))
# Set chi
var aNx = -1
var aCx = -1
var ang = 0.0
select {(resno=res) and (chain=iChain) and (file=f) and (model=m) and base}
if ((aC1 and {purine}).size > 0) {
aNx = get_atom_rcn( res, iChain, "N9")
aCx = get_atom_rcn( res, iChain, "C8")
ang = (toA ? kPuA : kPuB)
}
else {
aNx = get_atom_rcn(res, iChain, "N1")
aCx = get_atom_rcn(res, iChain, "C6")
ang = (toA ? kPyA : kPyB)
}
set_dihedral_atoms(aO4, aC1, aNx, aCx, ang)
# Set pucker 3' endo or 2' endo
var pSet = {aC1 or aC2 or aO2}
select pSet or {(resno=res) and (chain=iChain)
and (file=f) and (model=m) and base}
set_dihedral_atoms(aC5, aC4, aO4, aC1, (toA ? kC5C4O4C1A : kC5C4O4C1B))
set_dihedral_atoms(aC4, aO4, aC1, aC2, (toA ? kC4O4C1C2A : kC4O4C1C2B))
set_dihedral_atoms(aC2, aO4, aC1, aNx, (toA ? kC2O4C1NxA : kC2O4C1NxB))
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)
select selsave
}
function adjust_nts(res5, res3, iChain, toab, a, s) {
# Collect any pairing
var w = array()
for (var i = res5; i <= res3; i++) {
w = w + [who_pairs(i, iChain)]
}
# Twist and turn
for (var i = res3; i >= res5; i--) {
var j = i-res5+1
if (toab.size > 0) {
to_ab_nt_res(i, -1, iChain, (toab == "A"))
if ((w[j])[1] >= 0) {
to_ab_nt_res((w[j])[1], -1, (w[j])[2], (toab == "A"))
}
}
}
for (var i = res5; i < res3; i++) {
base_stack_res(i, i+1, iChain, iChain, s, a, false)
}
# Restore pairings
for (var i = res3; i >= res5; i--) {
var j = i-res5+1
if ((w[j])[1] >= 0) {
pair_it_res((w[j])[1], i, -1, (w[j])[2], iChain)
}
}
# Clean up
for (var i = res3; i >= res5; i--) {
var j = i-res5+1
fix_p_res(i, iChain, true)
if ((w[j])[1] >= 0) {
fix_p_res((w[j])[1], (w[j])[2], true)
}
}
}
#########################################################
### STAND ALONE GENERAL PURPOSE FUNCTIONS ###
#########################################################
function is_form_a( iResno, iChain) {
var aO4 = get_atom_rcn( iResno, iChain, "O4\'")
var aC1 = get_atom_rcn( iResno, iChain, "C1\'")
var aC2 = get_atom_rcn( iResno, iChain, "C2\'")
var aC3 = get_atom_rcn( iResno, iChain, "C3\'")
return (angle(aO4, aC1, aC2, aC3) < 0.0)
}
function is_r_res( iResno, iChain) {
var f = (_frameID/1000000)
var m = (_frameID%1000000)
return ({(resno=iResno) and (chain=iChain) and (file=f)
and (model=m) and purine}.size > 0)
}
function repair_p_res(res, iChain) {
var aP = get_atom_rcn( res, iChain, "P")
plico_minimize( {connected(aP) or aP})
}
function who_pairs(iRes, iChain) {
var aC4or6 = get_atom_rcn( iRes, iChain, "C4")
var aN1or3 = get_atom_rcn( iRes, iChain, "N1")
if ({aN1or3 and purine}.size = 0) {
aC4or6 = get_atom_rcn( iRes, iChain, "C6")
aN1or3 = get_atom_rcn( iRes, iChain, "N3")
}
if (aN1or3.size > 0) {
var near = within(3.2, aN1or3) and {resno!=iRes} and {element="N"}
for (var i = 1; i <= near.size; i++) {
var dist = distance(near[i], aN1or3)
var ang = abs(angle(near[i], aN1or3, aC4or6))
if (ang > 150) {
return [near[i].resno, near[i].chain, dist, ang]
}
}
}
return [-1, aC4or6.chain, -1, -1]
}
function who_almost_pairs(iRes, iChain) {
var aC2 = get_atom_rcn( iRes, iChain, "C2")
var aC4or6 = get_atom_rcn( iRes, iChain, "C4")
var aN1or3 = get_atom_rcn( iRes, iChain, "N1")
var pname = "C6"
if ({aN1or3 and purine}.size = 0) {
aC4or6 = get_atom_rcn( iRes, iChain, "C6")
aN1or3 = get_atom_rcn( iRes, iChain, "N3")
pname = "C4"
}
if (aN1or3.size > 0) {
var near = within(3.4, aN1or3) and {resno!=iRes} and {element="N"}
for (var i = 1; i <= near.size; i++) {
var aC6or4 = get_atom_rcn(near[i].resno, near[i].chain, pname)
var aC2p = get_atom_rcn(near[i].resno, near[i].chain, "C2")
var dist = distance(near[i], aN1or3)
var puang = abs(angle(near[i], aN1or3, aC4or6))
var pyang = abs(angle(aC6or4, near[i], aN1or3))
var dihedral = abs(angle(aC2p, near[i], aN1or3, aC2))
if ((puang > 110) and (pyang > 110) and (dihedral < 40)) {
return [near[i].resno, near[i].chain, dist, puang]
}
}
}
return [-1, aC4or6.chain, -1, -1]
}
function who_stacks(iRes, iChain) {
var ret = array()
var aNear = ((within(4.0, {(resno=iRes) and base}) )
and {base} and {not resno=iRes})
var done = array()
for (var i = 1; i <= aNear.size; i++) {
var jRes = aNear[i].resno
if (not done.find(jRes)) {
var jChain = aNear[i].chain
var as = gen_as(iRes, jRes, iChain, jChain)
var d = distance({(resno=iRes) and base}, {(resno=jRes) and base})
var a1 = angle(as[2], as[3], as[4])
var a2 = angle(as[5], as[4], as[3])
var dh = angle(as[5], as[4], as[3], as[2])
var bset = ((connected(as[3]) and not as[2])
or (connected(as[4]) and not as[5]))
var a3 = angle(bset[1], bset[2], bset[3])
var a4 = angle(bset[2], bset[3], bset[4])
var isStacked = true
# Bases are parallel as sin(a1) = sin(a2) and sin(a3) = sin(a4)
if (abs(sin(a1)-sin(a2)) > 20) {
isStacked = false
}
if (abs(sin(a3)-sin(a4)) > 20) {
isStacked = false
}
# Bases are stacked as d*sin(a1) < 6.0 and d3 = 0.0
if (d*sin(a1) > 6.2) {
isStacked = false
}
if (abs(dh) > 30) {
#isStacked = false
}
if (isStacked) {
ret += aNear[i].resno
}
done += jRes
}
}
return ret
}
function match_nt(mask, nt) {
var ret = false
switch (mask) {
case "A":
case "U":
case "C":
case "G":
ret = (mask = nt)
break
case "N":
ret = true
break
case "M":
ret = ((nt=="A") or (nt=="C"))
break
case "Y":
ret = ((nt=="U") or (nt=="C"))
break
case "R":
ret = ((nt=="A") or (nt=="G"))
break
}
return ret
}
# Calls function match_nt above
function select_seqs(seq, r5, r3, iChain, f, m) {
select none
for (var i = r5; i < r3; i++) {
var j = 0
for (; j < seq.size; j++) {
if ((gNTres[i+j])[2] >= 0) {
break
}
if ((i+j) >= r3) {
break
}
if (not match_nt(seq[j+1], gSeq[i+j])) { # CALL
break
}
}
if (j == seq.size) {
print format("%s at %d (%s-%s-%s)%s", seq, i,
gSeq[i-1],
gSeq[i][i+j-1],
gSeq[i+j], aster)
var rset = {(resno=i) and (chain=iChain) and (file=f) and (model=m)}
rset.selected = true
}
}
}
function find_tetras(seq, r5, r3, based) {
select none
for (var i = r5; i < r3; i++) {
var j = 0
for (; j < seq.size; j++) {
if ((i+j) >= r3) {
break
}
if (not match_nt(seq[j+1], gSeq[i+j])) { # CALL
break
}
}
if (j == seq.size) {
if (based) {
var ends = gSeq[i-1] + gSeq[i+j]
var bads = ["AA","GG","AG","GA"]
if (bads.count(ends) > 0) {
continue
}
}
print format("%s at %d (%s-%s-%s)%s", seq, i,
gSeq[i-1],
gSeq[i][i+j-1],
gSeq[i+j], aster)
}
}
}
# Calls is_form_a
function select_b_form_nts(iChain) {
var f = (_frameID/1000000)
var m = (_frameID%1000000)
select none
for (var i = {(chain=iChain) and (file=f) and (model=m)}.resno.min;
i <= {(chain=iChain) and (file=f) and (model=m)}.resno.max; i++) {
if (not is_form_a(i, iChain)) { # <== external call
print format("Res %d is form B", i)
var rset = {(resno=i) and (chain=iChain) and (file=f) and (model=m)}
rset.selected = true
}
}
}
function select_3ward_atom(ar3, ares, iChain) {
var f = (_frameID/1000000)
var m = (_frameID%1000000)
var i = ar3.resno
var aP = get_atom_rcn( i, iChain, "P")
switch(ar3.atomName) {
case "O3\'" :
select {(resno>i) and (resno<ares) and (chain=iChain)
and (file=f) and (model=m)}
break
case "P" :
select {(resno>=i) and (resno<ares) and (chain=iChain)
and (file=f) and (model=m)}
break
case "O5\'" :
case "C5\'" :
case "C4\'" :
select {(resno>=i) and (resno<ares) and (chain=iChain)
and (file=f) and (model=m) and not (connected(aP) or aP)}
break
case "C3\'" :
var aO3 = get_atom_rcn( i, iChain, "O3\'")
select {((resno>i) and (resno<ares) and (chain=iChain)
and (file=f) and (model=m)) or aO3}
break
}
}
function select_5ward_atom(ar5, ares, iChain) {
var f = (_frameID/1000000)
var m = (_frameID%1000000)
var i = ar5.resno
var aP = get_atom_rcn( i, iChain, "P")
switch(ar5.atomName) {
case "O3\'" :
select {(resno<=i) and (resno>ares) and (chain=iChain)
and (file=f) and (model=m)}
break
case "P" :
case "O5\'" :
case "C5\'" :
select {((resno<i) and (resno>ares) and (chain=iChain)
and (file=f) and (model=m)) or (connected(aP) or aP)}
break
case "C4\'" :
var aC5 = get_atom_rcn( i, iChain, "C5\'")
select {((resno<i) and (resno>ares) and (chain=iChain)
and (file=f) and (model=m)) or (connected(aP) or aP or aC5)}
break
}
}
function plot_ab_chi( iChain) {
var f = (_frameID/1000000)
var m = (_frameID%1000000)
select none
for (var i = {(chain=iChain) and (file=f) and (model=m)}.resno.min;
i <= {(chain=iChain) and (file=f) and (model=m)}.resno.max; i++) {
var aO4 = get_atom_rcn(i, iChain, "O4\'")
var aC1 = get_atom_rcn(i, iChain, "C1\'")
var isR = ({aC1 and purine}.size > 0)
var a1or9 = (isR ? "N9" : "N1")
var a6or8 = (isR ? "C8" : "C6")
var aN = get_atom_rcn(i, iChain, a1or9)
var aC = get_atom_rcn(i, iChain, a6or8)
var aO3 = get_atom_rcn(i, iChain, "O3\'")
var aC3 = get_atom_rcn(i, iChain, "C3\'")
var aC4 = get_atom_rcn(i, iChain, "C4\'")
var aC5 = get_atom_rcn(i, iChain, "C5\'")
var chi = angle(aO4, aC1, aN, aC)
aN.vx = chi
var aorb = angle(aO3, aC3, aC4, aC5)
aN.vy = aorb
select ADD aN
}
plot properties vx vy resno
set echo top left
echo "vx = base chi angle vy = a ==> b form"
}
function move_p_to_close_o3( aP, aO3, ares, force) {
var cp = aP.xyz
var aC3 = {connected(aO3) and (atomname="C3\'")}
select aP
set_distance_atoms(aO3, aP, 1.74)
set_angle_atoms(aC3, aO3, aP, 109.0)
var pt = aP.xyz
aP.xyz = cp
var rotors = gen_nt_rotors(aP.resno, ares, aP.chain)
move_atom_nt(aP.atomIndex, pt, aP.resno-1, rotors, force)
##fix_p_res(aP.resno, aP.chain, true)
}
function move_o3_to_close_p( aO3, aP, ares, force) {
var cp = aO3.xyz
var aO5 = {connected(aP) and (atomname="O5\'")}
select aO3
set_distance_atoms(aP, aO3, 1.74)
set_angle_atoms(aO5, aP, aO3, 109.0)
var pt = aO3.xyz
aO3.xyz = cp
var rotors = gen_nt_rotors(ares, aO3.resno, aO3.chain)
move_atom_nt(aO3.atomIndex, pt, aO3.resno+1, rotors, force)
##fix_p_res(aP.resno, aP.chain, true)
}
# Select mobile before calling
function pivot_to_close_atoms( aMov, aStat, aPivot, dist) {
var axis = cross(aStat.xyz, aMov.xyz) - aPivot.xyz
var dir = 1
var d = distance(aMov, aStat)
if (d > dist) {
while (d > dist) {
rotateSelected @aPivot @axis @dir
var nd = distance(aMov, aStat)
if (nd > d) {
rotateSelected @aPivot @axis @{-dir}
if (dir == 1) {
rotateSelected @aPivot @axis @{-dir}
dir = -1
continue
}
else {
break
}
}
d = nd
}
}
else {
while (d < dist) {
rotateSelected @aPivot @axis @dir
var nd = distance(aMov, aStat)
if (nd < d) {
rotateSelected @aPivot @axis @{-dir}
if (dir == 1) {
rotateSelected @aPivot @axis @{-dir}
dir = -1
continue
}
else {
break
}
}
d = nd
}
}
}
function print_adjacent_vs( iChain) {
var f = (_frameID/1000000)
var m = (_frameID%1000000)
var tNN = 0
var taA = 0
var tbA = 0
var taD = 0
var tbD = 0
var tabD = 0
var tc = 0.0
var rmin = {(chain=iChain) and (file=f) and (model=m)}.resno.min
var rmax = {(chain=iChain) and (file=f) and (model=m)}.resno.max
for (var i = rmin; i < rmax; i++) {
var aC4 = get_atom_rcn(i, iChain, "C4\'")
var aC1 = get_atom_rcn(i, iChain, "C1\'")
var isR = ({aC1 and purine}.size > 0)
var a1or9 = (isR ? "N9" : "N1")
var aN = get_atom_rcn(i, iChain, a1or9)
var aC4p = get_atom_rcn(i+1, iChain, "C4\'")
var aC1p = get_atom_rcn(i+1, iChain, "C1\'")
isR = ({aC1p and purine}.size > 0)
a1or9 = (isR ? "N9" : "N1")
var aNp = get_atom_rcn(i+1, iChain, a1or9)
if ((aC4.size > 0) and (aC4p.size > 0)) {
var NN = distance(aN, aNp)
tNN += NN
var aA = angle(aC1, aN, aNp)
taA += aA
var bA = angle(aC1p, aNp, aN)
tbA += bA
var aD = angle(aC4, aC1, aN, aNp)
if ((aD < 0) and (taD > 0)) {
aD += 360
}
taD += aD
var bD = angle(aC4p, aC1p, aNp, aN)
if ((bD < 0) and (tbD > 0)) {
bD += 360
}
tbD += bD
var abD = angle(aC1, aN, aNp, aC1p)
tabD += abD
tc++
print format("%6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %d%s %s",
NN, aA, aD, bA, bD, abD, tc, aC1.group, aC1p.group)
}
}
print format("v1=%6.2f v2=%6.2f v3=%6.2f v4=%6.2f v5=%6.2f v6=%6.2f",
tNN/tc, taA/tc, taD/tc, tbA/tc, tbD/tc, tabD/tc)
}
function measure_adjacent_vs(r5, iChain) {
var aC4 = get_atom_rcn(r5, iChain, "C4\'")
var aC1 = get_atom_rcn(r5, iChain, "C1\'")
var isR = ({aC1 and purine}.size > 0)
var aN1or9 = get_atom_rcn(r5, iChain, (isR ? "N9" : "N1"))
var aC8or6 = get_atom_rcn(r5, iChain, (isR ? "C8" : "C6"))
var aC4p = get_atom_rcn(r5+1, iChain, "C4\'")
var aC1p = get_atom_rcn(r5+1, iChain, "C1\'")
isR = ({aC1p and purine}.size > 0)
var aN1or9p = get_atom_rcn(r5+1, iChain, (isR ? "N9" : "N1"))
var aC8or6p = get_atom_rcn(r5+1, iChain, (isR ? "C8" : "C6"))
measure @aN1or9 @aN1or9p
measure @aC1 @aN1or9 @aN1or9p
measure @aC4 @aC1 @aN1or9 @aN1or9p
measure @aC1p @aN1or9p @aN1or9
measure @aC4p @aC1p @aN1or9p @aN1or9
measure @aC1p @aN1or9p @aN1or9 @aC1
measure @aC4 @aC1 @aN1or9 @aC8or6
measure @aC4p @aC1p @aN1or9p @aC8or6p
}
function print_pair_vs( iChain) {
var f = (_frameID/1000000)
var m = (_frameID%1000000)
var tNN = 0.0
var taA = 0.0
var tbA = 0.0
var taD = 0.0
var tbD = 0.0
var tabD = 0.0
var tc = 0.0
var taChi = 0.0
var tbChi = 0.0
var rmin = {(chain=iChain) and (file=f) and (model=m)}.resno.min
var rmax = {(chain=iChain) and (file=f) and (model=m)}.resno.max
for (var i = rmin; i < rmax; i++) {
var aC4 = get_atom_rcn(i, iChain, "C4\'")
if (aC4.size > 0) {
var w = who_pairs(i, iChain) # CALL
var aC1 = get_atom_rcn(i, iChain, "C1\'")
var isR = ({aC1 and purine}.size > 0)
var a1or9 = (isR ? "N9" : "N1")
var aN = get_atom_rcn(i, iChain, a1or9)
var a6or8 = (isR ? "C8" : "C6")
var aC = get_atom_rcn(i, iChain, a6or8)
var aC4p = get_atom_rcn(w[1], w[2], "C4\'")
var aC1p = get_atom_rcn(w[1], w[2], "C1\'")
a1or9 = (isR ? "N1" : "N9") # rev
var aNp = get_atom_rcn(w[1], w[2], a1or9)
a6or8 = (isR ? "C6" : "C8") # rev
var aCp = get_atom_rcn(w[1], w[2], a6or8)
var NN = distance(aN, aNp)
tNN += NN
var aA = angle(aC1, aN, aNp)
taA += aA
var bA = angle(aC1p, aNp, aN)
tbA += bA
var aD = angle(aC4, aC1, aN, aNp)
if ((aD < 0) and (taD > 0)) {
aD += 360
}
taD += aD
var bD = angle(aC4p, aC1p, aNp, aN)
if ((bD < 0) and (tbD > 0)) {
bD += 360
}
tbD += bD
var abD = angle(aC1, aN, aNp, aC1p)
tabD += abD
var aChi = angle(aC4, aC1, aN, aC)
taChi += aChi
var bChi = angle(aC4p, aC1p, aNp, aCp)
tbChi += bChi
tc++
print format(
"%6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %d%s %d%s",
NN, aA, aD, bA, bD, abD, aChi, bChi, tc, aC1.group, w[1],
aC1p.group)
}
}
print format(
"v1=%6.2f v2=%6.2f v3=%6.2f v4=%6.2f v5=%6.2f v6=%6.2f v7=%6.2f v8=%6.2f",
tNN/tc, taA/tc, taD/tc, tbA/tc, tbD/tc, tabD/tc, taChi/tc, tbChi/tc)
}
function measure_pair_vs(i, j, iChain, jChain) {
var aC4 = get_atom_rcn(i, iChain, "C4\'")
var aC4p = get_atom_rcn(j, jChain, "C4\'")
if ((aC4.size > 0) and (aC4p.size > 0)) {
var aC1 = get_atom_rcn(i, iChain, "C1\'")
var isR = ({aC1 and purine}.size > 0)
var a1or9 = (isR ? "N9" : "N1")
var aN = get_atom_rcn(i, iChain, a1or9)
var a6or8 = (isR ? "C8" : "C6")
var aC = get_atom_rcn(i, iChain, a6or8)
var aC1p = get_atom_rcn(j, jChain, "C1\'")
isR = ({aC1p and purine}.size > 0)
a1or9 = (isR ? "N9" : "N1")
var aNp = get_atom_rcn(j, jChain, a1or9)
a6or8 = (isR ? "C8" : "C6")
var aCp = get_atom_rcn(j, jChain, a6or8)
measure @aN @aNp
measure @aC1p @aNp @aN
measure @aC4p @aC1p @aNp @aN
measure @aC1 @aN @aNp
measure @aC4 @aC1 @aN @aNp
measure @aC1p @aNp @aN @aC1
measure @aC4p @aC1p @aNp @aCp
measure @aC4 @aC1 @aN @aC
}
else {
print "No pair found"
}
}
function minimize_for_collision( r, iChain) {
var f = (_frameID/1000000)
var m = (_frameID%1000000)
fix_p_res( r, iChain, true)
# External
var cset = (within(kCtolerance, {resno=r}) and not {resno=r} and not
connected({resno=r}) and {(chain=iChain) and (file=f) and (model=m)})
if (cset.size > 0) {
for (var i = 0; i < cset.size; i++) {
plico_minimize( {(resno=r) or (resno=@{cset[i].resno})})
}
}
# Internal
cset = (within(kCtolerance, {(resno=r) and (atomName="OP?")}) and not
{(atomName="OP?") or (atomName="P")}
and {(chain=iChain) and (file=f) and (model=m)})
if (cset.size > 0) {
for (var i = 0; i < cset.size; i++) {
plico_minimize( {(resno=r) or (resno=@{cset[i].resno})})
}
}
}
function eval_pairing( seq, r5, r3, len) {
var val = 0
for (var i = 0; i < len; i++) {
var c5 = seq[r5+i]
var c3 = seq[r3-i]
if ((c5 == "A") and (c3 == "U")) {
val += 2
}
else if ((c5 == "U") and (c3 == "A")) {
val += 2
}
else if ((c5 == "G") and (c3 == "C")) {
val += 3
}
else if ((c5 == "C") and (c3 == "G")) {
val += 3
}
else if ((c5 == "G") and (c3 == "U")) {
val += 2
}
else if ((c5 == "U") and (c3 == "G")) {
val += 1
}
else if ((c5 == "A") and (c3 == "A")) {
val -= 3
}
else if ((c5 == "G") and (c3 == "G")) {
val -= 3
}
}
return val
}
function update_atomnos(iChain) {
print "Update atomnos..."
var f = {chain=iChain}.file
var m = {chain=iChain}.model
var minres = {chain=iChain}.resno.min
var maxres = {chain=iChain}.resno.max
var b = {(chain=iChain) and (file=f) and (model=m)}.atomno.min
for (var r = {chain=iChain}.resno.min; r <= {chain=iChain}.resno.max; r++) {
var rset = {(resno=r) and (chain=iChain) and (file=f) and (model=m)}
for (var n = rset.atomno.min; n <= rset.atomno.max; n++) {
var a = {(resno=r) and (chain=iChain) and (file=f)
and (model=m) and (atomno=n)}
a.atomno = -b
b++
}
}
for (var i = -b; i < 0; i++) {
var a = {(chain=iChain) and (file=f) and (model=m)
and (atomno=i)}
a.atomno *= -1
}
print "Updated"
}
function replane_bases(iChain) {
for (var i = {chain=iChain}.resno.min; i <= {chain=iChain}.resno.max; i++) {
plico_minimize( {resno=i and base})
}
}
function measure_p_dihedrals(i, iChain) {
var j = i-1
var as = array()
as += get_atom_rcn( j, iChain, "C4\'")
as += get_atom_rcn( j, iChain, "C3\'")
as += get_atom_rcn( j, iChain, "O3\'")
as += get_atom_rcn( i, iChain, "P")
as += get_atom_rcn( i, iChain, "O5\'")
as += get_atom_rcn( i, iChain, "C5\'")
as += get_atom_rcn( i, iChain, "C4\'")
as += get_atom_rcn( i, iChain, "C3\'")
for (var k = 1; k <= (as.size-3); k++) {
measure @{as[k+0]} @{as[k+1]} @{as[k+2]} @{as[k+3]}
}
}
function pair_stem( r5, r3, ar5, iChain) {
# Pair entire stem
var c = r3 - r5 + 1
for (var k = (c\2)-1 ; k >= 0; k--) {
pair_it_res(r5+k, r3-k, ar5, iChain, iChain) # CALL
}
}
# From 2LU0
function make_uncg_loop(res5, ares5, iChain) {
var f = (_frameID/1000000)
var m = (_frameID%1000000)
var va = array()
va[1] = [4.46, 100.5, -83.4, 114.4, 99.8, 39.2, -22.2, -11.8]
va[2] = [7.89, 93.6, 50.7, 29.7, 75.3, 163.6, 14.2 -22.2]
va[3] = [7.57, 90.9, -11.9, 59.5, 75.6, 74.1, 33.5, 14.2]
va[4] = [6.09, 124.9, -5.8, 46.4, -72.6, -27.1, -138.3, 33.5]
va[5] = [5.1, 64.0, -102.5, 98.2, 87.0, 58.6, -9.2, -138.3]
for (var i = 1; i <= 5; i++) {
var as = gen_as(res5+i-2, res5+i-1, iChain, iChain)
var vs = array()
vs[1] = (va[i])[1] # distance res5 N9or1 and res3 N9or1
vs[2] = (va[i])[2] # angle res5 N9or1 and res3 N9or1 C1
vs[3] = (va[i])[3] # dihedral res5 N9or1 and res3 N9or1 C1 C4
vs[4] = (va[i])[4] # angle res5 C1 N9or1 and res3 N9or1
vs[5] = (va[i])[5] # dihedral res5 C4 N9or1 C1 and res3 N9or1
vs[6] = (va[i])[6] # dihedral res5 N9or1 C1 and res3 N9or1 C1
vs[7] = (va[i])[7] # dihedral chi res3
vs[8] = (va[i])[8] # dihedral chi res5
select {(resno < @{res5+i-1}) and (resno > ares5)
and (chain=iChain) and (file=f) and (model=m)}
move_it(as, vs)
fix_p_res(res5+i-1, iChain, true)
}
}
# From 2LU0
function make_gnra_loop(res5, ares5, iChain) {
var f = (_frameID/1000000)
var m = (_frameID%1000000)
var va = array()
va[1] = [4.37, 93.9, -84.9, 117.9, 99.1, 42.0, -19.4, 3.2]
va[2] = [7.94, 77.5, 47.5, 58.9, 87.1, 179.5, 4.5 -19.4]
va[3] = [4.9, 102.9, -98.6, 79.2, 112.4, 61.2, -29.6, 4.5]
va[4] = [4.97, 124.5, -70.0, 99.8, 93.2, 34.0, 2.5, -29.6]
va[5] = [4.31, 76.0, -98.4, 106.6, 97.2, 45.6, -15.5, 2.5]
for (var i = 1; i <= 5; i++) {
var as = gen_as(res5+i-2, res5+i-1, iChain, iChain)
var vs = array()
vs[1] = (va[i])[1] # distance res5 N9or1 and res3 N9or1
vs[2] = (va[i])[2] # angle res5 N9or1 and res3 N9or1 C1
vs[3] = (va[i])[3] # dihedral res5 N9or1 and res3 N9or1 C1 C4
vs[4] = (va[i])[4] # angle res5 C1 N9or1 and res3 N9or1
vs[5] = (va[i])[5] # dihedral res5 C4 N9or1 C1 and res3 N9or1
vs[6] = (va[i])[6] # dihedral res5 N9or1 C1 and res3 N9or1 C1
vs[7] = (va[i])[7] # dihedral chi res3
vs[8] = (va[i])[8] # dihedral chi res5
select {(resno < @{res5+i-1}) and (resno > ares5)
and (chain=iChain) and (file=f) and (model=m)}
move_it(as, vs)
fix_p_res(res5+i-1, iChain, true)
}
}
# end of plicoNTcommon.spt
