From 671fe4b1b714aefa4c066f6f1ee76cd5bee2ac5c Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Fri, 3 Feb 2023 11:56:40 -0600
Subject: [PATCH] Improve pdb and list readers
---
mrdna/RELEASE-VERSION | 2 +-
mrdna/readers/segmentmodel_from_lists.py | 32 ++++++++++++++++++++----
mrdna/readers/segmentmodel_from_pdb.py | 13 ++++------
3 files changed, 33 insertions(+), 14 deletions(-)
diff --git a/mrdna/RELEASE-VERSION b/mrdna/RELEASE-VERSION
index 5f778f8..6644d6c 100644
--- a/mrdna/RELEASE-VERSION
+++ b/mrdna/RELEASE-VERSION
@@ -1 +1 @@
-1.0a.dev55
+1.0a.dev57
diff --git a/mrdna/readers/segmentmodel_from_lists.py b/mrdna/readers/segmentmodel_from_lists.py
index c1ce49b..25dbafb 100644
--- a/mrdna/readers/segmentmodel_from_lists.py
+++ b/mrdna/readers/segmentmodel_from_lists.py
@@ -6,7 +6,7 @@ import os,sys
import scipy
from ..segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment
-from ..arbdmodel.coords import quaternion_from_matrix
+from ..arbdmodel.coords import quaternion_from_matrix, rotationAboutAxis
from .. import get_resource_path
ref_stack_position = np.array((-2.41851735, -0.259761333, 3.39999978))
@@ -162,7 +162,7 @@ def basepairs_and_stacks_to_helixmap(basepairs,stacks_above):
return helixmap, helixrank, is_fwd, intrahelical
-def set_splines(seg, coordinate, hid, hmap, hrank, fwd, orientation=None):
+def set_splines(seg, coordinate, hid, hmap, hrank, fwd, basepair, orientation=None):
maxrank = np.max( hrank[hmap==hid] )
if maxrank == 0:
ids = np.where((hmap == hid))[0]
@@ -189,8 +189,20 @@ def set_splines(seg, coordinate, hid, hmap, hrank, fwd, orientation=None):
q = quaternion_from_matrix( orientation[ids[0]] )
if last_q is not None and last_q[1:].dot(q[1:]) < 0:
q = -1*q
+
+ ## Average quaterion with reverse direction
+ bp = basepair[ids[0]]
+ if bp >= 0:
+ q2 = quaternion_from_matrix( orientation[bp].dot(rotationAboutAxis(np.array((1,0,0)),180)) )
+ if last_q is None:
+ if q[1:].dot(q2[1:]) < 0:
+ q2 = -1*q2
+ else:
+ if last_q[1:].dot(q2[1:]) < 0:
+ q2 = -1*q2
+ q = (q + q2)*0.5
quats.append(q)
- ## TODO: average quaterion with reverse direction
+ last_q = q
coords = np.array(coords)
seg.set_splines(contours,coords)
@@ -326,7 +338,7 @@ def model_from_basepair_stack_3prime(coordinate, basepair, stack, three_prime,
for hid in double_stranded_helices:
seg = DoubleStrandedSegment(name=str(hid),
num_bp = np.sum(hmap==hid)//2)
- set_splines(seg, coordinate, hid, hmap, hrank, fwd, orientation)
+ set_splines(seg, coordinate, hid, hmap, hrank, fwd, basepair, orientation)
assert(hid == len(doubleSegments))
doubleSegments.append(seg)
@@ -336,7 +348,7 @@ def model_from_basepair_stack_3prime(coordinate, basepair, stack, three_prime,
for hid in single_stranded_helices:
seg = SingleStrandedSegment(name=str(hid),
num_nt = np.sum(hmap==hid))
- set_splines(seg, coordinate, hid, hmap, hrank, fwd, orientation)
+ set_splines(seg, coordinate, hid, hmap, hrank, fwd, basepair, orientation)
assert(hid == len(doubleSegments) + len(singleSegments))
singleSegments.append(seg)
@@ -422,6 +434,16 @@ def model_from_basepair_stack_3prime(coordinate, basepair, stack, three_prime,
**model_parameters )
+ model._reader_list_coordinates = coordinate
+ model._reader_list_basepair = basepair
+ model._reader_list_stack = stack
+ model._reader_list_three_prime = three_prime
+ model._reader_list_sequence = sequence
+ model._reader_list_orientation = orientation
+ model._reader_list_hmap = hmap
+ model._reader_list_fwd = fwd
+ model._reader_list_hrank = hrank
+
if sequence is None:
for s in model.segments:
s.randomize_unset_sequence()
diff --git a/mrdna/readers/segmentmodel_from_pdb.py b/mrdna/readers/segmentmodel_from_pdb.py
index 6460027..5f8e67b 100644
--- a/mrdna/readers/segmentmodel_from_pdb.py
+++ b/mrdna/readers/segmentmodel_from_pdb.py
@@ -38,8 +38,7 @@ for key,val in refUnis.items():
ref_res = val.residues[0]
names = ref_res.atoms.select_atoms( stack_text ).atoms.names
ref_name_ids_inv[key] = np.argsort(np.argsort(names))
-
-
+ val.atoms.positions = val.atoms.positions.dot(CanonicalNucleotideFactory.DefaultOrientation) # return to default orientation
## TODO: possibly improve ref_stack_position
# avg_stack_position = [ru.residues[0].atoms.select_atoms(stack_text).center_of_mass()
@@ -297,7 +296,10 @@ def SegmentModelFromPdb(*args,**kwargs):
three_prime = -np.ones(bps.shape, dtype=np.int)
for s in u.segments:
r = s.atoms.select_atoms("name C1'").resindices
+ if u.segments[0] == s: print(r)
+ # three_prime[r[:-1]] = r[1:]
three_prime[r[:-1]] = r[1:]
+ if u.segments[0] == s: print(three_prime[r])
""" ## Remove connections between resids that are very far away
dist = centers[r[1:]]-centers[r[:-1]]
@@ -321,17 +323,12 @@ def SegmentModelFromPdb(*args,**kwargs):
seq = [resname_to_key[rn]
for rn in u.select_atoms("name C1'").resnames]
- ## Get orientations
- orientations = []
- for R in transforms:
- orientations.append( R.T.dot(CanonicalNucleotideFactory.DefaultOrientation) )
-
return model_from_basepair_stack_3prime(centers,
bps,
stacks,
three_prime,
seq,
- np.array(orientations)
+ np.array([t.T for t in transforms])
)
if __name__ == "__main__":
--
GitLab