add

mrdna/RELEASE-VERSION

-1.0a.dev71
+1.0a.dev74

mrdna/readers/.ipynb_checkpoints/__init__-checkpoint.py

+## TODO: make module this package conform to a single style for input/output
+
+def read_cadnano(json_file, sequence=None, fill_sequence='T', **model_parameters):
+    from .cadnano_segments import read_json_file
+    from .cadnano_segments import read_model as model_from_cadnano_json
+
+    data = read_json_file(json_file)
+    return model_from_cadnano_json(data, sequence, fill_sequence, **model_parameters)
+
+def read_vhelix(maya_file, **model_parameters):
+    from .polygon_mesh import parse_maya_file, convert_maya_bases_to_segment_model
+
+    maya_bases = parse_maya_file(maya_file)
+    model = convert_maya_bases_to_segment_model( maya_bases, **model_parameters )
+    return model
+
+def read_list(infile,**model_parameters):
+    import numpy as np
+    try:
+        data = np.loadtxt(infile)
+    except:
+        data = np.loadtxt(infile,skiprows=1) # strip header
+    coords = data[:,:3]*10
+    bp = data[:,3]
+    stack = data[:,4]
+    three_prime = data[:,5]
+
+    from .segmentmodel_from_lists import model_from_basepair_stack_3prime
+    return model_from_basepair_stack_3prime(coords, bp, stack, three_prime)
+
+
+def read_atomic_pdb(pdb_file, **model_parameters):
+    from .segmentmodel_from_pdb import SegmentModelFromPdb
+    return SegmentModelFromPdb(pdb_file)
+
+def read_oxdna(coordinate_file, topology_file, idealized_coordinate_file=None, **model_parameters):
+    """ Idealized coordinate file: coordinates for detecting stacks and base pairs, defaults to coordinate_file """ 
+    from .segmentmodel_from_oxdna import mrdna_model_from_oxdna
+    return mrdna_model_from_oxdna(coordinate_file, topology_file, idealized_coordinate_file, **model_parameters)

mrdna/readers/.ipynb_checkpoints/segmentmodel_from_lists-checkpoint.py

+# -*- coding: utf-8 -*-
+
+import pdb
+import numpy as np
+import os,sys
+import scipy
+
+from mrdna import logger, devlogger
+from ..segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment
+from ..arbdmodel.coords import quaternion_from_matrix, rotationAboutAxis, quaternion_slerp
+from .. import get_resource_path
+
+ref_stack_position = np.array((-2.41851735, -0.259761333, 3.39999978))
+
+def _three_prime_list_to_five_prime(three_prime):
+    five_prime = -np.ones(three_prime.shape, dtype=int)
+    has_three_prime = np.where(three_prime >= 0)[0]
+    five_prime[three_prime[has_three_prime]] = has_three_prime
+    return five_prime  
+
+def _primes_list_to_strands(three_prime, five_prime):
+    five_prime_ends = np.where(five_prime < 0)[0]
+    strands = []
+    strand_is_circular = []
+    
+    idx_to_strand = -np.ones(three_prime.shape, dtype=int)
+
+    def build_strand(nt_idx, conditional):
+        strand = [nt_idx]
+        idx_to_strand[nt_idx] = len(strands)
+        while conditional(nt_idx):
+            nt_idx = three_prime[nt_idx]
+            strand.append(nt_idx)
+            idx_to_strand[nt_idx] = len(strands)
+        strands.append( np.array(strand, dtype=int) )
+
+    for nt_idx in five_prime_ends:
+        build_strand(nt_idx,
+                     lambda nt: three_prime[nt] >= 0)
+        strand_is_circular.append(False)
+
+    while True:
+        ## print("WARNING: working on circular strand {}".format(len(strands)))
+        ids = np.where(idx_to_strand < 0)[0]
+        if len(ids) == 0: break
+        build_strand(ids[0],
+                     lambda nt: three_prime[nt] >= 0 and \
+                     idx_to_strand[three_prime[nt]] < 0)
+        strand_is_circular.append(True)
+
+    return strands, strand_is_circular
+
+def find_stacks(centers, transforms):
+
+    ## Find orientation and center of each nucleotide
+    expected_stack_positions = []
+    for R,c in zip(transforms,centers):
+        expected_stack_positions.append( c + ref_stack_position.dot(R) )
+
+    expected_stack_positions = np.array(expected_stack_positions, dtype=np.float32)
+
+    dists = scipy.spatial.distance_matrix(expected_stack_positions, centers)
+    dists = dists + 5*np.eye(len(dists))
+    idx1, idx2 = np.where(dists < 3.5)
+
+    ## Convert distances to stacks
+    stacks_above = -np.ones(len(centers), dtype=int)
+    _z = np.array((0,0,1))
+    for i in np.unique(idx1):
+        js = idx2[ idx1 == i ]
+        with np.errstate(divide='ignore',invalid='ignore'):
+            angles = [np.arccos( transforms[j].T.dot( transforms[i].dot(_z) ).dot( _z ) ) for j in js]
+        angles = np.array( angles )
+        tmp = np.argmin(dists[i][js] + 1.0*angles)
+        j = js[tmp]
+        stacks_above[i] = j
+
+    return stacks_above
+
+def basepairs_and_stacks_to_helixmap(basepairs,stacks_above):
+
+    helixmap = -np.ones(basepairs.shape, dtype=int)
+    helixrank = -np.ones(basepairs.shape)
+    is_fwd = np.ones(basepairs.shape, dtype=int)
+    
+    ## Remove stacks with nts lacking a basepairs
+    nobp = np.where(basepairs < 0)[0]
+    stacks_above[nobp] = -1
+    stacks_with_nobp = np.in1d(stacks_above, nobp)
+    stacks_above[stacks_with_nobp] = -1
+
+    end_ids = np.where( (stacks_above < 0)*(basepairs >= 0) )[0]
+
+    hid = 0
+    for end in end_ids:
+        if helixmap[end] >= 0:
+            continue
+        rank = 0
+        nt = basepairs[end]
+        bp = basepairs[nt]
+        assert( bp == end )
+        if helixmap[nt] >= 0 or helixmap[bp] >= 0:
+            logger.warning(f'Ill-formed helix: problematic basepair or stacking data near nucleotide {nt} or {bp}... skipping')
+            continue
+        # assert(helixmap[nt] == -1)
+        # assert(helixmap[bp] == -1)
+        helixmap[nt] = helixmap[bp] = hid
+        helixrank[nt] = helixrank[bp] = rank
+        is_fwd[bp] = 0
+        rank +=1
+
+        _tmp = [(nt,bp)]
+        
+        while stacks_above[nt] >= 0:
+            nt = stacks_above[nt]
+            if basepairs[nt] < 0: break
+            bp = basepairs[nt]
+            if helixmap[nt] >= 0 or helixmap[bp] >= 0:
+                logger.warning(f'Ill-formed helix: problematic basepair or stacking data near nucleotide {nt} or {bp}... skipping')
+                break
+            helixmap[nt] = helixmap[bp] = hid
+            helixrank[nt] = helixrank[bp] = rank
+            is_fwd[bp] = 0
+            _tmp.append((nt,bp))
+            rank +=1
+
+        hid += 1
+
+    ## Create "helix" for each circular segment
+    intrahelical = []
+    processed = set()
+    unclaimed_bases = np.where( (basepairs >= 0)*(helixmap == -1) )[0]
+    for nt0 in unclaimed_bases:
+        if nt0 in processed: continue
+
+        nt = nt0
+        all_nts = [nt]
+
+        rank = 0
+        nt = nt0
+        bp = basepairs[nt]
+        if helixmap[nt] >= 0 or helixmap[bp] >= 0:
+            logger.warning(f'Ill-formed cylic helix: problematic basepair or stacking data near nucleotide {nt} or {bp}... skipping')
+            continue
+        helixmap[nt] = helixmap[bp] = hid
+        helixrank[nt] = helixrank[bp] = rank
+        is_fwd[bp] = 0
+        rank +=1
+        processed.add(nt)
+        processed.add(bp)
+
+        counter = 0
+        while stacks_above[nt] >= 0:
+            lastnt = nt
+            nt = stacks_above[nt]
+            bp = basepairs[nt]
+            if nt == nt0 or nt == basepairs[nt0]:
+                intrahelical.append((lastnt,nt0))
+                break
+                
+            assert( bp >= 0 )
+            if helixmap[nt] >= 0 or helixmap[bp] >= 0:
+                logger.warning(f'Ill-formed cyclic helix: problematic basepair or stacking data near nucleotide {nt} or {bp}... skipping')
+                break
+            
+            helixmap[nt] = helixmap[bp] = hid
+            helixrank[nt] = helixrank[bp] = rank
+            is_fwd[bp] = 0
+            processed.add(nt)
+            processed.add(bp)
+            rank +=1
+        hid += 1
+
+    return helixmap, helixrank, is_fwd, intrahelical
+
+
+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]
+        pos = np.mean( [coordinate[r,:] for r in ids ], axis=0 )
+        coords = [pos,pos]
+        contours = [0,1]
+        if orientation is not None:
+            ids = np.where((hmap == hid) * fwd)[0]
+            assert( len(ids) == 1 )
+            q = quaternion_from_matrix( orientation[ids[0]] )
+            quats = [q, q]
+            coords[-1] = pos + orientation[ids[0]].dot(np.array((0,0,1)))
+
+    else:
+        coords,contours,quats = [[],[],[]]
+        last_q = None
+        for rank in range(int(maxrank)+1):
+            ids = np.where((hmap == hid) * (hrank == rank))[0]
+                
+            coords.append(np.mean( [coordinate[r,:] for r in ids ], axis=0 ))
+            contours.append( float(rank+0.5)/(maxrank+1) )
+            if orientation is not None:
+                ids = np.where((hmap == hid) * (hrank == rank) * fwd)[0]
+                assert(len(ids) == 1)
+                q = quaternion_from_matrix( orientation[ids[0]] )
+
+                if last_q is not None and last_q.dot(q) < 0:
+                    q = -q
+
+                ## Average quaterion with reverse direction
+                bp = basepair[ids[0]]
+                if bp >= 0:
+                    bp_o = orientation[bp].dot(rotationAboutAxis(np.array((1,0,0)),180))
+                    q2 = quaternion_from_matrix( bp_o )
+                    if q.dot(q2) < 0:
+                        q2 = -q2
+
+                    ## probably good enough, but slerp is better: q = (q + q2)*0.5
+                    q = quaternion_slerp(q,q2,0.5)
+
+                quats.append(q)
+                last_q = q
+
+    coords = np.array(coords)
+    seg.set_splines(contours,coords)
+    if orientation is not None:
+        quats = np.array(quats)
+        seg.set_orientation_splines(contours,quats)
+
+    seg.start_position = coords[0,:]
+    seg.end_position = coords[-1,:]
+
+
+def model_from_basepair_stack_3prime(coordinate, basepair, stack, three_prime,
+                                     sequence=None, orientation=None,
+                                     max_basepairs_per_bead = 5,
+                                     max_nucleotides_per_bead = 5,
+                                     local_twist = False,
+                                     dimensions=(5000,5000,5000),
+                                     **model_parameters):
+    """ 
+    Creates a SegmentModel object from lists of each nucleotide's
+    basepair, its stack (on 3' side) and its 3'-connected nucleotide
+
+    The first argument should be an N-by-3 numpy array containing the
+    coordinate of each nucleotide, where N is the number of
+    nucleotides. The following three arguments should be integer lists
+    where the i-th element corresponds to the i-th nucleotide; the
+    list element should the integer index of the corresponding
+    basepaired / stacked / phosphodiester-bonded nucleotide. If there
+    is no such nucleotide, the value should be -1.
+
+    Args:
+        basepair:  List of each nucleotide's basepair's index
+        stack:  List containing index of the nucleotide stacked on the 3' of each nucleotide
+        three_prime:  List of each nucleotide's the 3' end of each nucleotide
+
+    Returns:
+        SegmentModel
+    """
+
+    """ Validate Input """
+    inputs = (basepair,three_prime)
+    try:
+        basepair,three_prime = [np.array(a,dtype=int) for a in inputs]
+    except:
+        raise TypeError("One or more of the input lists could not be converted into a numpy array")
+    inputs = (basepair,three_prime)
+    coordinate = np.array(coordinate)
+
+    if np.any( [len(a.shape) > 1 for a in inputs] ):
+        raise ValueError("One or more of the input lists has the wrong dimensionality")
+
+    if len(coordinate.shape) != 2:
+        raise ValueError("Coordinate array has the wrong dimensionality")
+
+    inputs = (coordinate,basepair,three_prime)
+    if not np.all(np.diff([len(a) for a in inputs]) == 0):
+        raise ValueError("Inputs are not the same length")
+        
+    num_nt = len(basepair)
+    if sequence is not None and len(sequence) != num_nt:
+        raise ValueError("The 'sequence' parameter is the wrong length {} != {}".format(len(sequence),num_nt))
+
+    if orientation is not None:
+        orientation = np.array(orientation)
+        if len(orientation.shape) != 3:
+            raise ValueError("The 'orientation' array has the wrong dimensionality (should be Nx3x3)")
+        if orientation.shape != (num_nt,3,3):
+            raise ValueError("The 'orientation' array is not properly formatted")
+
+    if stack is None:
+        if orientation is not None:
+            stack = find_stacks(coordinate, orientation)
+        else:
+            ## Guess stacking based on 3' connectivity
+            stack = np.array(three_prime,dtype=int) # Assume nts on 3' ends are stacked
+            _stack_below = _three_prime_list_to_five_prime(stack)
+            _has_bp = (basepair >= 0)
+            _nostack = np.where( (stack == -1)*_has_bp )[0]
+            _has_stack_below = _stack_below[basepair[_nostack]] >= 0
+            _nostack2 = _nostack[_has_stack_below]
+            stack[_nostack2] = basepair[_stack_below[basepair[_nostack2]]]
+
+    else:
+        try:
+            stack = np.array(stack,dtype=int)
+        except:
+            raise TypeError("The 'stack' array could not be converted into a numpy integer array")
+
+        if len(stack.shape) != 1:
+            raise ValueError("The 'stack' array has the wrong dimensionality")
+
+        if len(stack) != num_nt:
+            raise ValueError("The length of the 'stack' array does not match other inputs")
+
+    bps = basepair              # alias
+
+    """ Fix stacks: require that the stack of a bp of a base's stack is its bp """
+    _has_bp =  (bps >= 0)
+    _has_stack = (stack >= 0)
+    _stack_has_basepair = (bps[stack] >= 0) * _has_stack
+    stack = np.where( (stack[bps[stack]] == bps) * _has_bp * _has_stack * _has_bp,
+                      stack, -np.ones(len(stack),dtype=int) )
+
+    five_prime = _three_prime_list_to_five_prime(three_prime)
+
+    """ Build map of dsDNA helices and strands """
+    hmap,hrank,fwd,intrahelical = basepairs_and_stacks_to_helixmap(bps,stack)
+    double_stranded_helices = np.unique(hmap[hmap >= 0])    
+    strands, strand_is_circular = _primes_list_to_strands(three_prime, five_prime)
+
+    """ Add ssDNA to hmap """
+    if len(double_stranded_helices) > 0:
+        hid = double_stranded_helices[-1]+1
+    else:
+        hid = 0
+    ss_residues = hmap < 0
+    #
+    if np.any(bps[ss_residues] != -1):
+        logger.warning(f'{np.sum(bps[ss_residues] != -1)} ssDNA nucleotides appear to have basepairs... ignoring')
+    
+    for s,c in zip(strands, strand_is_circular):
+        strand_segment_ends = [i for i in np.where( np.diff(hmap[s]) != 0 )[0]] + [len(s)-1]
+        seg_start = 0
+        for i in strand_segment_ends:
+            if hmap[s[i]] < 0:
+                ## Found single-stranded segment
+                ids = s[seg_start:i+1]
+                assert( np.all(hmap[ids] == -1) )
+                hmap[ids] = hid
+                hrank[ids] = np.arange(i+1-seg_start)
+                hid+=1
+            seg_start = i+1
+
+    if len(double_stranded_helices) > 0:
+        single_stranded_helices = np.arange(double_stranded_helices[-1]+1,hid)
+    else:
+        single_stranded_helices = np.arange(hid)
+
+    ## Create double-stranded segments
+    doubleSegments = []
+    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, basepair, orientation)
+
+        assert(hid == len(doubleSegments))
+        doubleSegments.append(seg)
+
+    ## Create single-stranded segments
+    singleSegments = []
+    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, basepair, orientation)
+
+        assert(hid == len(doubleSegments) + len(singleSegments))
+        singleSegments.append(seg)
+
+    ## Find crossovers and 5prime/3prime ends
+    crossovers,prime5,prime3 = [[],[],[]]
+    for s,c in zip(strands,strand_is_circular):
+        tmp = np.where(np.diff(hmap[s]) != 0)[0]
+        for i in tmp:
+            crossovers.append( (s[i],s[i+1]) )
+        if c:
+            if hmap[s[-1]] != hmap[s[0]]:
+                crossovers.append( (s[-1],s[0]) )
+        else:
+            prime5.append(s[0])
+            prime3.append(s[-1])
+
+    ## Add connections
+    allSegments = doubleSegments+singleSegments
+
+    for r1,r2 in crossovers:
+        seg1,seg2 = [allSegments[hmap[i]] for i in (r1,r2)]
+        nt1,nt2 = [hrank[i] for i in (r1,r2)]
+        f1,f2 = [fwd[i] for i in (r1,r2)]
+
+        ## Handle connections at the ends
+        is_terminal1 = (nt1,f1) in ((0,0),(seg1.num_nt-1,1))
+        is_terminal2 = (nt2,f2) in ((0,1),(seg2.num_nt-1,0))
+
+        print(seg1,seg2, r1, r2, is_terminal1, is_terminal2)
+        if is_terminal1 or is_terminal2:
+            """ Ensure that we don't have three-way dsDNA junctions """
+            if is_terminal1 and (bps[r1] >= 0) and (five_prime[bps[r1]] >= 0) and (three_prime[r1] >= 0):
+                if (bps[five_prime[bps[r1]]] >= 0) and (bps[three_prime[r1]] >= 0):
+                    # is_terminal1 = (three_prime[r1] == bps[five_prime[bps[r1]]])
+                    is_terminal1 = hmap[five_prime[bps[r1]]] == hmap[three_prime[r1]]
+            if is_terminal2 and (bps[r2] >= 0) and (three_prime[bps[r2]] >= 0) and (five_prime[r2] >= 0):
+                if (bps[three_prime[bps[r2]]] >= 0) and (bps[five_prime[r2]] >= 0):
+                    # is_terminal2 = (five_prime[r2] == bps[three_prime[bps[r2]]])
+                    is_terminal2 = hmap[three_prime[bps[r2]]] == hmap[five_prime[r2]]
+                
+            """ Place connection """
+            if is_terminal1 and is_terminal2:
+                end1 = seg1.end3 if f1 else seg1.start3
+                end2 = seg2.start5 if f2 else seg2.end5
+                seg1._connect_ends( end1, end2, type_='intrahelical')
+            else:
+                seg1.add_crossover(nt1,seg2,nt2,[f1,f2],type_="terminal_crossover")
+        else:
+            seg1.add_crossover(nt1,seg2,nt2,[f1,f2])
+
+    ## Add 5prime/3prime ends
+    for r in prime5:
+        seg = allSegments[hmap[r]]
+        seg.add_5prime(hrank[r],fwd[r])
+    for r in prime3:
+        seg = allSegments[hmap[r]]
+        seg.add_3prime(hrank[r],fwd[r])
+
+    ## Add intrahelical connections to circular helical sections
+    for nt0,nt1 in intrahelical:
+        seg = allSegments[hmap[nt0]]
+        assert( seg is allSegments[hmap[nt1]] )
+        if three_prime[nt0] >= 0:
+            if hmap[nt0] == hmap[three_prime[nt0]]:
+                seg.connect_end3(seg.start5)
+
+        bp0,bp1 = [bps[nt] for nt in (nt0,nt1)]
+        if three_prime[bp1] >= 0:
+            if hmap[bp1] == hmap[three_prime[bp1]]:
+                seg.connect_start3(seg.end5)
+
+    ## Assign sequence
+    if sequence is not None:
+        for hid in range(len(allSegments)):
+            resids = np.where( (hmap==hid)*(fwd==1) )[0]
+            s = allSegments[hid]
+            s.sequence = [sequence[r] for r in sorted(resids,key=lambda x: hrank[x])]
+
+
+    ## Build model
+    model = SegmentModel( allSegments,
+                          max_basepairs_per_bead = max_basepairs_per_bead,
+                          max_nucleotides_per_bead = max_nucleotides_per_bead,
+                          local_twist = local_twist,
+                          dimensions = dimensions,
+                          **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_five_prime = five_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()
+
+    return model
+
+def UNIT_circular():
+    """ Make a circular DNA strand, with dsDNA in the middle """
+    coordinate = [(0,0,3.4*i) for i in range(7)]
+    stack = -1*np.ones(len(coordinate))
+    three_prime = [ 1, 2, 4,-1, 6, 3, 0]
+    basepair    = [-1,-1, 3, 2, 5, 4,-1]
+    stack       = [-1,-1, 4,-1,-1, 3,-1]
+    for i in [3,5]:
+        coordinate[i] = (1,0,3.4*i)
+
+    model = model_from_basepair_stack_3prime(coordinate, basepair, stack, three_prime,
+                                             max_basepairs_per_bead=1,
+                                             max_nucleotides_per_bead=1,
+                                             local_twist=True)
+
+    model.simulate(output_name='circular', directory='unittest')

mrdna/readers/.ipynb_checkpoints/segmentmodel_from_oxdna_pinyi-checkpoint.py

+from mrdna import logger, devlogger
+from .segmentmodel_from_lists import model_from_basepair_stack_3prime
+from ..arbdmodel.coords import rotationAboutAxis
+
+from oxlibs import *
+import numpy as np
+from scipy.spatial import distance_matrix
+
+_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, virt2nuc=None, **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)
+    
+    if idealized_coordinate_file is not None:
+        import pickle
+        df = pickle.load(virt2nuc)
+        vh_vb2nuc_rev, vhelix_pattern=df
+        
+
+
+    ## 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)])
+    
+    def _get_bp(sequence=None):
+        dists = distance_matrix(r,basepair_pos) + np.eye(len(r))*1000
+        dists = 0.5*(dists + dists.T)
+        
+        bp = np.array([np.argmin(da) for da in dists])
+
+        for i,j in enumerate(bp):
+            if j == -1: continue
+            # devlogger.info(f'bp {i} {j} {dists[i,j]}')
+            if dists[i,j] > 2:
+                bp[i] = -1
+            elif bp[j] != i:
+                bpj = bp[j]
+                logger.warning( " ".join([str(_x) for _x in ["Bad pair", i, j, bp[i], bp[j], dists[i,j], dists[j,i], dists[bpj,j], dists[j,bpj]]]) )
+
+        for i,j in enumerate(bp):
+            if j == -1: continue
+            if bp[j] != i:
+                bpj = bp[j]
+                logger.warning( " ".join([str(_x) for _x in ["Bad pair2", i, j, bp[i], bp[j], dists[i,j], dists[j,i], dists[bpj,j], dists[j,bpj]]]) )
+                raise Exception
+
+        if sequence is not None:
+            seq = sequence
+            bp_seq = sequence[bp]
+            bp_seq[bp==-1] = 'X'
+            bad_bps = np.where( (bp >= 0) &
+                                (((seq == 'C') & (bp_seq != 'G')) |
+                                 ((seq == 'G') & (bp_seq != 'C')) |
+                                 ((seq == 'T') & (bp_seq != 'A')) | 
+                                 ((seq == 'U') & (bp_seq != 'A')) |
+                                 ((seq == 'A') & ((bp_seq != 'T') | (bp_seq != 'U')))
+                                 ) )[0]
+            bp[bp[bad_bps]] = -1
+            bp[bad_bps] = -1
+            
+        return bp
+
+    def _get_stack():
+        dists = distance_matrix( r + 3.5*normal_dir + 2.1*perp_dir -1*base_dir, r ) + np.eye(len(r))*1000
+        stack = np.array([np.argmin(da) for da in dists])
+        for i,j in enumerate(stack):
+            if dists[i,j] > 8:
+                stack[i] = -1
+            elif i < 10:
+                ## development info
+                # devlogger.info([i,j,dists[i,j]])
+                # dr = r[j] - (r[i] - normal_dir[i]*3.4 + perp_dir[i]*1 + base_dir[i]*1)
+                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()
+      
+    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
+
+    def _debug_write_bonds():
+        from ..arbdmodel import ParticleType, PointParticle, ArbdModel, Group
+        bond = tuple()
+        b_t = ParticleType('BASE')
+        p_t = ParticleType('PHOS')
+
+        parts = []
+        for i,(r0,r_bp,three_prime0,bp0,stack0,seq0) in enumerate(zip(r,basepair_pos, three_prime, bp, stack, seq)):
+            p = PointParticle(p_t, name='PHOS', position = r0, resid=i)
+            b = PointParticle(b_t, name=seq0, position = 0.5*(r0+r_bp), resid=i)
+            parts.extend((p,b))
+
+        model = ArbdModel(parts)
+        model.writePdb('test.pdb')
+
+        for i,(r0,r_bp,three_prime0,bp0,stack0) in enumerate(zip(r,basepair_pos, three_prime, bp, stack)):
+            model.add_bond(parts[2*i],parts[2*i+1],bond)
+            j = three_prime0
+            if j >= 0:
+                model.add_bond(parts[2*i],parts[2*j],bond)
+            j = bp0
+            if j >= 0:
+                model.add_bond(parts[2*i+1],parts[2*j+1],bond)
+        model.writePsf('test.psf')
+
+        model.bonds = []
+        for i,(r0,r_bp,three_prime0,bp0,stack0) in enumerate(zip(r,basepair_pos, three_prime, bp, stack)):
+            j = stack0
+            if j >= 0:
+                model.add_bond(parts[2*i],parts[2*j],bond)
+        model.writePsf('test.stack.psf')
+    ## _debug_write_bonds()
+
+    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 )
+
+    """
+    model.DEBUG = True
+    model.generate_bead_model(1,1,False,True,one_bead_per_monomer=True)
+    for seg in model.segments:
+        for bead in seg.beads:
+            bead.position = bead.position + np.random.standard_normal(3)
+    
+    simulate( model, output_name='test', directory='test4' )
+    """
+    return model
+
+if __name__ == "__main__":
+    mrdna_model_from_oxdna("0-from-collab/nanopore.oxdna","0-from-collab/nanopore.top")
+    # mrdna_model_from_oxdna("2-oxdna.manual/output/from_mrdna-oxdna-min.last.conf","0-from-collab/nanopore.top")

mrdna/readers/segmentmodel_from_oxdna_pinyi.py

@@ -2,6 +2,7 @@ from mrdna import logger, devlogger
 from .segmentmodel_from_lists import model_from_basepair_stack_3prime
 from ..arbdmodel.coords import rotationAboutAxis
 
+from oxlibs import *
 import numpy as np
 from scipy.spatial import distance_matrix
 
@@ -10,11 +11,8 @@ _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, **model_parameters):
+def mrdna_model_from_oxdna(coordinate_file, topology_file, ,idealized_coordinate_file=None,virt2nuc=None, **model_parameters):
     """ Construct an mrdna model from oxDNA coordinate and topology files """
-    if idealized_coordinate_file is None:
-        idealized_coordinate_file = coordinate_file
-
     top_data = np.loadtxt(topology_file, skiprows=1,
                           unpack=True,
                           dtype=np.dtype('i4,U1,i4,i4')
@@ -141,6 +139,16 @@ def mrdna_model_from_oxdna(coordinate_file, topology_file, idealized_coordinate_
     
     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)