Skip to content
Snippets Groups Projects
Commit f72d8f48 authored by pinyili2's avatar pinyili2
Browse files

remove redundant files

parent ab0e3b72
No related branches found
No related tags found
1 merge request!1Pinyili2
Hide whitespace changes
Inline Side-by-side
mrdna/readers/test/segmentmodel_from_cadnano_old.py deleted 100644 → 0
+ 0
− 207
View file @ ab0e3b72
# -*- coding: utf-8 -*-
import pdb
import numpy as np
import os,sys
from glob import glob
import re
import pandas as pd
pd.options.mode.chained_assignment = None # default='warn'
from ..segmentmodel_from_lists import model_from_basepair_stack_3prime
from ...arbdmodel.coords import readArbdCoords, readAvgArbdCoords, rotationAboutAxis
from ...segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment
from ...model.dna_sequence import m13 as m13seq
import json
import re
## TODO: separate SegmentModel from ArbdModel so multiple parts can be combined
## TODO: catch circular strands in "get_5prime" cadnano calls
## TODO: handle special motifs
## - doubly-nicked helices
## - helices that should be stacked across an empty region (crossovers from and end in the helix to another end in the helix)
## - circular constructs
def read_json_file(filename):
try:
with open(filename) as ch:
json_data = json.load(ch)
except:
with open(filename) as ch:
content = ""
for l in ch:
l = re.sub(r"'", r'"', l)
# https://stackoverflow.com/questions/4033633/handling-lazy-json-in-python-expecting-property-name
# l = re.sub(r"{\s*(\w)", r'{"\1', l)
# l = re.sub(r",\s*(\w)", r',"\1', l)
# l = re.sub(r"(\w):", r'\1":', l)
content += l+"\n"
json_data = json.loads(content)
return json_data
def get_lattice(part):
lattice_type = None
_gt = part.getGridType()
try:
lattice_type = _gt.name.lower()
except:
if _gt == 1:
lattice_type = 'square'
elif _gt == 2:
lattice_type = 'honeycomb'
else:
print("WARNING: unable to determine cadnano part lattice type")
return lattice_type
def cadnano_parts(json_data):
import cadnano
from cadnano.document import Document
try:
doc = Document()
cadnano.fileio.v3decode.decode(doc, json_data)
decoder = 3
except:
doc = Document()
cadnano.fileio.v2decode.decode(doc, json_data)
decoder = 2
parts = [p for p in doc.getParts()]
if len(parts) != 1:
raise Exception("Only documents containing a single cadnano part are implemented at this time.")
part = parts[0]
if decoder == 2:
""" It seems cadnano2.5 (as of ce6ff019) does not set the EulerZ for square lattice structures correctly, doing so here """
l = get_lattice(part)
if l == 'square':
for id_num in part.getIdNums():
if part.vh_properties.loc[id_num,'eulerZ'] == 0:
part.vh_properties.loc[id_num,'eulerZ'] = 360*(6/10.5)
return part,decoder
def get_helix_angle(part, helix_id, indices):
""" Get "start_orientation" for helix """
# import ipdb
# ipdb.set_trace()
""" FROM CADNANO2.5
+ angle is CCW
- angle is CW
Right handed DNA rotates clockwise from 5' to 3'
we use the convention the 5' end starts at 0 degrees
and it's pair is minor_groove_angle degrees away
direction, hence the minus signs. eulerZ
"""
hp, bpr, tpr, eulerZ, mgroove = part.vh_properties.loc[helix_id,
['helical_pitch',
'bases_per_repeat',
'turns_per_repeat',
'eulerZ',
'minor_groove_angle']]
twist_per_base = tpr*360./bpr
# angle = eulerZ - twist_per_base*indices + 0.5*mgroove + 180
angle = eulerZ + twist_per_base*indices - 0.5*mgroove
return rotationAboutAxis(np.array((0,0,1)),angle)
def nttype_v2(scafs):
def judge(i):
if i ==[-1,-1,-1,-1]:
return 0
else: return 1
n=np.array([judge(i) for i in scafs])
return n
def gen_prop_table_cad2(json_data,part):
df=pd.DataFrame(json_data["vstrands"])
vslist=df.set_index("num")
props = part.getModelProperties().copy()
try:
if props.get('point_type') == PointType.ARBITRARY:
# TODO add code to encode Parts with ARBITRARY point configurations
raise NotImplementedError("Not implemented")
except:
try:
vh_props, origins = part.helixPropertiesAndOrigins()
except:
origins = {hid:part.getVirtualHelixOrigin(hid)[:2] for hid in part.getidNums()}
scaf_id=[nttype_v2(vslist['scaf'][i]) for i in vslist.index]
stap_id=[nttype_v2(vslist['stap'][i]) for i in vslist.index]
cad_bps=part.getIndices(0)
vslist["scafnt"]=np.sum(np.array(scaf_id),axis=1)
vslist["stapnt"]=np.sum(np.array(stap_id),axis=1)
totnt=np.sum(vslist["scafnt"])+np.sum(vslist["stapnt"])
is_scaf=np.zeros(totnt,dtype=bool)
is_scaf[0:np.sum(vslist["scafnt"])]=1
nt_prop=pd.DataFrame(index=range(totnt),columns=["vh","zid","is_scaf","r","bp","stack","threeprime","seq","orientation"])
nt_prop["is_scaf"]=is_scaf
tot_id=scaf_id+stap_id
vhi,zidi=np.where(np.array(scaf_id)==1)
vhj,zidj=np.where(np.array(stap_id)==1)
vhi=vslist.index[vhi]
vhj=vslist.index[vhj]
nt_prop["vh"]=list(vhi)+list(vhj)
nt_prop["zid"]=list(zidi)+list(zidj)
vhzid=list(zip(nt_prop["vh"],nt_prop["zid"]))
nt_prop["r"]=[part.getCoordinate(helix_id, indices) for helix_id, indices in zip(nt_prop["vh"],nt_prop["zid"])]
nt_prop["orientation"]=[get_helix_angle(part, helix_id, indices) for helix_id,indices in zip(nt_prop["vh"],nt_prop["zid"])]
nt_prop=nt_prop.fillna(-1)
for i in range(int(len(vhzid)/2)):
try:
bp1,bp2=(i,1+i+vhzid[i+1:].index(vhzid[i]))
nt_prop["bp"][bp1]=bp2
nt_prop["bp"][bp2]=bp1
except:
pass
tprime_list=-np.ones(len(nt_prop.index),dtype=int)
index2=list(zip(vhzid,nt_prop["is_scaf"]))
for i in range(len(nt_prop.index)):
((m,n),p)=list(zip(vhzid,nt_prop["is_scaf"]))[i]
if p==True:
k,l=(vslist["scaf"][m])[n][2:]
if k!=-1 and l!=-1:
n=index2.index(((k,l),True))
tprime_list[i]=int(n)
else:
k,l=(vslist["stap"][m])[n][2:]
if k!=-1 and l!=-1:
n=index2.index(((k,l),False))
tprime_list[i]=int(n)
nt_prop["threeprime"]=tprime_list
#(n,)=np.where(nt_prop["threeprime"]==-1)
#stackid=nt_prop["bp"][[list(nt_prop["threeprime"]).index(i) for i in n]]
index3=dict(zip(nt_prop.index,zip(nt_prop["vh"],nt_prop["is_scaf"])))
index3[-1]=(-1,-1)
nt_prop["stack"]=np.where([index3[i]==index3[nt_prop["threeprime"][i]] for i in (nt_prop.index)],nt_prop["threeprime"],-1)
return nt_prop
def mrdna_model_from_cadnano(json_file,seq=None,**model_parameters):
json_data=read_json_file(json_file)
part,decoder=cadnano_parts(json_data)
if decoder==2:
nt_prop=gen_prop_table_cad2(json_data,part)
elif decoder==3:
nt_prop=gen_prop_table_cad3(json_data,part)
if seq is None:
if nt_prop["seq"][0]==-1:
seq=None
else:
seq=nt_prop["seq"]
r=np.array(list(nt_prop['r']))
bp=np.array(list(nt_prop['bp']))
three_prime=np.array((list(nt_prop["threeprime"])))
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=None, **model_parameters )
return model
mrdna/readers/test/segmentmodel_from_oxdna_pinyi.py deleted 100644 → 0
+ 0
− 205
View file @ ab0e3b72
from mrdna import logger, devlogger
from .segmentmodel_from_lists import model_from_basepair_stack_3prime
from ..arbdmodel.coords import rotationAboutAxis
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,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(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)
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)
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
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
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' )
"""
model._dataframe=nt_prop
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% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment