cadnano,virt2nuc

mrdna/readers/segmentmodel_from_cadnano.py

@@ -173,7 +173,7 @@ def gen_prop_table(part):
     nt_prop=pd.DataFrame(id_series)
     nt_prop.reset_index(names=list(range(len(nt_prop.index))),inplace=True)
     nt_prop["seq"]=-1
-    ind_tuple=[(nt_prop["vh"][i],nt_prop["zid"][i],nt_prop["fwd"][i]) for i in nt_prop.index]
+    ind_tuple=list(zip(nt_prop["vh"],nt_prop["zid"],nt_prop["fwd"]))
     stacks=[]
     for i in list(nt_prop["stack_tuple"]):
         if i ==-1:
@@ -189,19 +189,22 @@ def gen_prop_table(part):
             tprime.append(ind_tuple.index(i))
     nt_prop["threeprime"]=tprime
     vhzid=list(zip(nt_prop["vh"],nt_prop["zid"]))
-    nt_prop["bp"]=-1
     nt_prop["orientation"]=[get_helix_angle(part, helix_id, int(float(indices))) for helix_id,indices in vhzid]
     nt_prop=nt_prop.fillna(-1)
-    for i in range(int(len(nt_prop.index)/2)):
+    counter=-1
+    bp=-np.ones(len(nt_prop.index),dtype=int)
+    bp_map=dict(zip(ind_tuple,nt_prop.index))
+    for i,j,k in ind_tuple:
+        counter+=1
         try:
-            bp1,bp2=(i,1+i+vhzid[i+1:].index(vhzid[i]))
-            nt_prop["bp"][bp1]=bp2
-            nt_prop["bp"][bp2]=bp1
+            bp[counter]=bp_map[(i,j,not(k))]
         except:
             pass
+    nt_prop["bp"]=bp
+
     return nt_prop
 
-def mrdna_model_from_cadnano(json_file,seq=None,**model_parameters):
+def mrdna_model_from_cadnano(json_file,seq=None,return_nt=True,**model_parameters):
     part=read_json_file(json_file)
     
     nt_prop=gen_prop_table(part)
@@ -219,4 +222,7 @@ def mrdna_model_from_cadnano(json_file,seq=None,**model_parameters):
     stack=np.array((list(nt_prop["stack"])))
     orientation=np.array(list(nt_prop["orientation"]))
     model = model_from_basepair_stack_3prime( r, bp, stack, three_prime, seq, orientation, **model_parameters )
-    return model
+    if return_nt==True:
+        return model,nt_prop
+    else:
+        return model

mrdna/readers/segmentmodel_from_lists.py

@@ -4,6 +4,7 @@ import pdb
 import numpy as np
 import os,sys
 import scipy
+import pandas as pd
 
 from mrdna import logger, devlogger
 from mrdna.segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment
@@ -233,7 +234,7 @@ def model_from_basepair_stack_3prime(coordinate, basepair, stack, three_prime,
                                      max_basepairs_per_bead = 5,
                                      max_nucleotides_per_bead = 5,
                                      local_twist = False,
-                                     dimensions=(5000,5000,5000),
+                                     dimensions=(5000,5000,5000),return_prop=False
                                      **model_parameters):
     """ 
     Creates a SegmentModel object from lists of each nucleotide's
@@ -474,8 +475,11 @@ def model_from_basepair_stack_3prime(coordinate, basepair, stack, three_prime,
     if sequence is None:
         for s in model.segments:
             s.randomize_unset_sequence()
-
-    return model
+    if return_prop is True:
+        nt_prop=pd.DataFrame({"r":coordinate,"bp":basepair,"stack":stack,"fwd":fwd,"threeprime":three_prime,"seq":sequence,"orientation":orientation})
+        return model,nt_prop
+    else:
+        return model
 
 def UNIT_circular():
     """ Make a circular DNA strand, with dsDNA in the middle """

mrdna/readers/segmentmodel_from_oxdna_pinyi.py

 from mrdna import logger, devlogger
 from .segmentmodel_from_lists import model_from_basepair_stack_3prime
 from ..arbdmodel.coords import rotationAboutAxis
+import pickle
+import pandas as pd
 
 from oxlibs import *
 import numpy as np
 from scipy.spatial import distance_matrix
+pd.options.mode.chained_assignment = None  # default='warn'
 
 _seq_to_int_dict = dict(A=0,T=1,C=2,G=3)
 _seq_to_int_dict = {k:str(v) for k,v in _seq_to_int_dict.items()}
 
 _yrot = rotationAboutAxis(axis=(0,1,0), angle=180).dot(rotationAboutAxis(axis=(0,0,1),angle=-40))
 
-def mrdna_model_from_oxdna(coordinate_file, topology_file, ,idealized_coordinate_file=None,virt2nuc=None, **model_parameters):
+def mrdna_model_from_oxdna(coordinate_file, topology_file,virt2nuc=None,get_nt_prop=False, **model_parameters):
     """ Construct an mrdna model from oxDNA coordinate and topology files """
     top_data = np.loadtxt(topology_file, skiprows=1,
                           unpack=True,
                           dtype=np.dtype('i4,U1,i4,i4')
                           )
-    conf_data = np.loadtxt(idealized_coordinate_file, skiprows=3)
-
-    ## Reverse direction so indices run 5'-to-3'
-    top_data = [a[::-1] for a in top_data]
-    conf_data = conf_data[::-1,:]
-
-    r = conf_data[:,:3] * 8.518
-    base_dir = conf_data[:,3:6]
-    # basepair_pos = r + base_dir*6.0
-    basepair_pos = r + base_dir*10.0
-    normal_dir = -conf_data[:,6:9]
-    perp_dir = np.cross(base_dir, normal_dir)
-    orientation = np.array([np.array(o).T.dot(_yrot) for o in zip(perp_dir,-base_dir,-normal_dir)])
-    
+    conf_data = np.loadtxt(coordinate_file, skiprows=3)
     def _get_bp(sequence=None):
         dists = distance_matrix(r,basepair_pos) + np.eye(len(r))*1000
         dists = 0.5*(dists + dists.T)
@@ -82,15 +72,84 @@ def mrdna_model_from_oxdna(coordinate_file, topology_file, ,idealized_coordinate
                 dr = r[j] - r[i]
                 # devlogger.info([normal_dir[i].dot(dr), perp_dir[i].dot(dr), base_dir[i].dot(dr)])
         return np.array(stack)
-
-    seq = top_data[1]
-    bp = _get_bp(seq)
-    stack = _get_stack()
+    
+    def _find_vh_vb_table(s,is_scaf):
+        L=[]
+        for i in list(s.keys()):
+            vh,zid=i
+            strand,indices=s[i]
+            if len(indices)==0:
+                continue
+            else:
+                if len(indices)==1:
+                    zids=[str(zid)]
+                else:
+                    zids=[str(zid)+"."+str(j) for j in range(len(indices))]
+                for index,z in zip(indices,zids):
+                    L.append(pd.Series({"index":index,"vh":vh,"zid":z,"strand":strand,"is_scaf":bool(is_scaf)}))
+        return L
+    def get_virt2nuc(virt2nuc,top_data):
+        vh_vb,pattern=pd.read_pickle(virt2nuc)
+        L1=_find_vh_vb_table(vh_vb._scaf,1)
+        L2=_find_vh_vb_table(vh_vb._stap,0)
+        nt_prop=pd.DataFrame(L1+L2)
+        nt_prop.set_index("index",inplace=True)
+        nt_prop.sort_index(inplace=True)
+        nt_prop["threeprime"]=top_data[2]
+        nt_prop["seq"]=top_data[1]
+        nt_prop["stack"]=top_data[2]
+        for i in nt_prop.index:
+            if nt_prop.loc[i]["threeprime"] in nt_prop.index:
+                if nt_prop.loc[nt_prop.loc[i]["threeprime"]]["vh"]!=nt_prop.loc[i]["vh"]:
+                    nt_prop["stack"][i]=-1
+        bp_map=dict(zip(zip(nt_prop["vh"],nt_prop["zid"],nt_prop["is_scaf"]),nt_prop.index))
+        bp=-np.ones(len(nt_prop.index),dtype=int)
+        counter=0
+        for i,j,k in zip(nt_prop["vh"],nt_prop["zid"],nt_prop["is_scaf"]):
+            try:
+                bp[counter]=bp_map[(i,j,not(k))]
+            except:
+                pass
+            counter+=1
+        nt_prop["bp"]=bp
+        return nt_prop
+    try:
+        nt_prop=get_virt2nuc(virt2nuc,top_data)
+        r=conf_data[:,:3] * 8.518
+        base_dir = conf_data[:,3:6]
+        # basepair_pos = r + base_dir*6.0
+        basepair_pos = r + base_dir*10.0
+        normal_dir = -conf_data[:,6:9]
+        perp_dir = np.cross(base_dir, normal_dir)
+        orientation = np.array([np.array(o).T.dot(_yrot) for o in zip(perp_dir,-base_dir,-normal_dir)])
+        seq=nt_prop["seq"]
+        bp=nt_prop["bp"]
+        stack=nt_prop["stack"]
+        three_prime=nt_prop["threeprime"]
+        nt_prop["r"]=r
+        nt_prop["orientation"]=orientation
+
+    except:
+        ## Reverse direction so indices run 5'-to-3'
+        top_data = [a[::-1] for a in top_data]
+        conf_data = conf_data[::-1,:]
+    
+        r = conf_data[:,:3] * 8.518
+        base_dir = conf_data[:,3:6]
+        # basepair_pos = r + base_dir*6.0
+        basepair_pos = r + base_dir*10.0
+        normal_dir = -conf_data[:,6:9]
+        perp_dir = np.cross(base_dir, normal_dir)
+        orientation = np.array([np.array(o).T.dot(_yrot) for o in zip(perp_dir,-base_dir,-normal_dir)])
+        seq = top_data[1]
+        bp = _get_bp(seq)
+        stack = _get_stack()
       
-    three_prime = len(r) - top_data[2] -1
-    five_prime = len(r) - top_data[3] -1
-    three_prime[three_prime >= len(r)] = -1
-    five_prime[five_prime >= len(r)] = -1
+        three_prime = len(r) - top_data[2] -1
+        five_prime = len(r) - top_data[3] -1
+        three_prime[three_prime >= len(r)] = -1
+        five_prime[five_prime >= len(r)] = -1
+        nt_prop=pd.DataFrame({"r":r,"bp":bp,"stack":stack,"threeprime":three_prime, "seq":seq,"orientation":orientation})
 
     def _debug_write_bonds():
         from ..arbdmodel import ParticleType, PointParticle, ArbdModel, Group
@@ -127,28 +186,9 @@ def mrdna_model_from_oxdna(coordinate_file, topology_file, ,idealized_coordinate
 
     logger.info(f'mrdna_model_from_oxdna: num_bp, num_ss_nt, num_stacked: {np.sum(bp>=0)//2} {np.sum(bp<0)} {np.sum(stack >= 0)}')
 
-    if coordinate_file != idealized_coordinate_file:
-        conf_data = conf_data[::-1,:]
 
-        r = conf_data[:,:3] * 8.518
-        base_dir = conf_data[:,3:6]
-        basepair_pos = r + base_dir*10.0
-        normal_dir = -conf_data[:,6:9]
-        perp_dir = np.cross(base_dir, normal_dir)
-        orientation = np.array([np.array(o).T.dot(_yrot) for o in zip(perp_dir,-base_dir,-normal_dir)])
-    
     model = model_from_basepair_stack_3prime( r, bp, stack, three_prime, seq, orientation, **model_parameters )
-
-    def find_cadnano_strands(virt2nuc):
-        import pickle
-        if virt2nuc is not None:
-            df = pickle.load(virt2nuc)
-            vh_vb2nuc_rev, vhelix_pattern=df #vhelix_pattern is the position on yzplane
-            
-
-        else:
-            print("virt2nuc not found")
-        
+     
     """
     model.DEBUG = True
     model.generate_bead_model(1,1,False,True,one_bead_per_monomer=True)
@@ -158,7 +198,10 @@ def mrdna_model_from_oxdna(coordinate_file, topology_file, ,idealized_coordinate
     
     simulate( model, output_name='test', directory='test4' )
     """
-    return model
+    if get_nt_prop is True:
+        return nt_prop,model
+    else:
+        return model
 
 if __name__ == "__main__":
     mrdna_model_from_oxdna("0-from-collab/nanopore.oxdna","0-from-collab/nanopore.top")