Skip to content
Snippets Groups Projects
Commit 12943abc authored by cmaffeo2's avatar cmaffeo2
Browse files

Add oxDNA reader

parent a18783e2
No related branches found
No related tags found
No related merge requests found
from mrdna import logger, devlogger
from .segmentmodel_from_lists import model_from_basepair_stack_3prime
from ..arbdmodel.coords import rotationAboutAxis
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, **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(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)])
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)}')
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")
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment