Skip to content
Snippets Groups Projects
Commit 0d1ed546 authored by cmaffeo2's avatar cmaffeo2
Browse files

Updated polygon_mesh reader to consider helix direction and to add 3/5-prime ends

parent aea81edd
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dnarbd/readers/polygon_mesh.py
+ 67
− 34
View file @ 0d1ed546
......@@ -209,7 +209,7 @@ def parse_maya_file(maya_file):
if conn is not None:
new_connections.append(conn)
connections = new_connections
## Assign base connectivity
for c in connections:
h1 = c.helix1
......@@ -274,21 +274,61 @@ def convert_maya_to_segments(helices):
# for b in upaired:
# print(hn, b.name, b.end3, b.end5)
z = [b.position[2] for b in paired_bases]
bot = paired_bases[ np.argmin(z) ]
top = paired_bases[ np.argmax(z) ]
if bot.name[0] == "f": bot = bot.basepair
if top.name[0] == "b": top = top.basepair
## Find helical axis
# if len(paired_bases) < 6:
# raise Exception("Helices shorter than 3 bp are not supported")
# coords = np.array([b.position for b in paired_bases])
# u,s,vh = np.linalg.svd(coords)
# axis = vh[0] / np.linalg.norm(vh[0])
fwd_bases = list(filter(lambda b:b.name[0] == 'f', paired_bases))
# rank = [np.dot(b.position,axis) for b in fwd_bases]
rank = [b.position[2] for b in fwd_bases]
order = np.argsort(rank)
ordered_fwd_bases = [fwd_bases[o] for o in order]
next_fwd_is_3end = [b1.end3==b2 for b1,b2 in
zip(ordered_fwd_bases[:-1],ordered_fwd_bases[1:])]
## Check if axis points the wrong way
assert( np.mean(next_fwd_is_3end) > 0.5 )
# if np.mean( next_fwd_is_3end ) < 0.5:
# # axis = axis * -1
# rank = rank[order[-1]]-rank
# order = reversed(order)
# ordered_fwd_bases = reversed(ordered_fwd_bases)
top = ordered_fwd_bases[-1]
bot = ordered_fwd_bases[0].basepair
assert( bot.name[0] == 'b' )
assert( top.name[0] == 'f' )
botPos = 0.5 * (bot.get_position() + bot.basepair.get_position())
topPos = 0.5 * (top.get_position() + top.basepair.get_position())
segments[hn] = DoubleStrandedSegment(hn, num_bp = len(z)//2,
segments[hn] = DoubleStrandedSegment(hn, num_bp = len(fwd_bases),
start_position = botPos,
end_position = topPos)
## Add 3prime/5prime ends
nt = 0
def add_primes(b,nt,fwd):
if b.end3 is None:
if not fwd: print("3prime at {}[{},{}]".format(hn,nt,fwd))
segments[hn].add_3prime(nt,fwd)
if b.end5 is None:
if not fwd: print("5prime at {}[{},{}]".format(hn,nt,fwd))
segments[hn].add_5prime(nt,fwd)
for b in ordered_fwd_bases:
add_primes(b,nt,True)
add_primes(b.basepair,nt,False)
nt+=1
segmentsBot[hn] = bot
segmentsTop[hn] = top
assert( bot.end3 not in paired_bases )
assert( top.end3 not in paired_bases )
""" Sanity check """
def end_crosses(end,direction):
......@@ -322,20 +362,28 @@ def convert_maya_to_segments(helices):
# print(bot.parent.name,bot.name, bot.end3.parent.name, bot.end3.name)
## TODO: clean this up
def get_end_positions(ssSeg, is_fwd):
start_idx, end_idx = (0,-1)
if not is_fwd:
start_idx, end_idx = (-1,0)
r0 = ssSeg[0].get_position()
if ssSeg[0].end5 is not None:
r0 = 0.75*r0 + 0.25*ssSeg[0].end5.get_position()
r1 = ssSeg[-1].get_position()
if ssSeg[-1].end3 is not None:
r1 = 0.75*r1 + 0.25*ssSeg[-1].end3.get_position()
return r0,r1
if top not in connectionToSingleStrand:
ssSeg = get_end3_ssDNA(top)
if len(ssSeg) > 0:
r0 = ssSeg[0].get_position()
if ssSeg[0].end5 is not None:
r0 = 0.75*r0 + 0.25*ssSeg[0].end5.get_position()
r1 = ssSeg[-1].get_position()
if ssSeg[-1].end3 is not None:
r1 = 0.75*r1 + 0.25*ssSeg[-1].end3.get_position()
num_nt = len(ssSeg)
if num_nt > 0:
r0,r1 = get_end_positions(ssSeg,is_fwd=True)
seg2 = SingleStrandedSegment("S{:03d}".format(len(ss_segments)),
start_position = r0,
end_position = r1,
num_nt = len(ssSeg))
num_nt = num_nt)
seg.connect_end3(seg2)
connectionToSingleStrand[top] = seg2
other = ssSeg[-1].end3
......@@ -360,12 +408,7 @@ def convert_maya_to_segments(helices):
if top.basepair not in connectionToSingleStrand:
ssSeg = get_end5_ssDNA(top.basepair)[::-1]
if len(ssSeg) > 0:
r0 = ssSeg[-1].get_position()
if ssSeg[-1].end5 is not None:
r0 = 0.75*r0 + 0.25*ssSeg[-1].end5.get_position()
r1 = ssSeg[0].get_position()
if ssSeg[0].end3 is not None:
r1 = 0.75*r1 + 0.25*ssSeg[0].end3.get_position()
r0,r1 = get_end_positions(ssSeg,is_fwd=False)
seg2 = SingleStrandedSegment("S{:03d}".format(len(ss_segments)),
start_position = r0,
end_position = r1,
......@@ -394,12 +437,7 @@ def convert_maya_to_segments(helices):
if bot not in connectionToSingleStrand:
ssSeg = get_end3_ssDNA(bot)
if len(ssSeg) > 0:
r0 = ssSeg[0].get_position()
if ssSeg[0].end5 is not None:
r0 = 0.75*r0 + 0.25*ssSeg[0].end5.get_position()
r1 = ssSeg[-1].get_position()
if ssSeg[-1].end3 is not None:
r1 = 0.75*r1 + 0.25*ssSeg[-1].end3.get_position()
r0,r1 = get_end_positions(ssSeg,is_fwd=True)
seg2 = SingleStrandedSegment("S{:03d}".format(len(ss_segments)),
start_position = r0,
end_position = r1,
......@@ -424,12 +462,7 @@ def convert_maya_to_segments(helices):
if bot.basepair not in connectionToSingleStrand:
ssSeg = get_end5_ssDNA(bot.basepair)[::-1]
if len(ssSeg) > 0:
r0 = ssSeg[-1].get_position()
if ssSeg[-1].end5 is not None:
r0 = 0.75*r0 + 0.25*ssSeg[-1].end5.get_position()
r1 = ssSeg[0].get_position()
if ssSeg[0].end3 is not None:
r1 = 0.75*r1 + 0.25*ssSeg[0].end3.get_position()
r0,r1 = get_end_positions(ssSeg,is_fwd=False)
seg2 = SingleStrandedSegment("S{:03d}".format(len(ss_segments)),
start_position = r0,
end_position = r1,
......
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