Commit ad9e6189 authored by cmaffeo2's avatar cmaffeo2
Browse files

Updated code that assigns num_nt to bead types to be more maintainable and fail-safe

parent 7bd92d97
......@@ -1600,42 +1600,48 @@ class SegmentModel(ArbdModel):
if self.DEBUG: print("Assigning bead types")
beadtype_s = dict()
def _assign_bead_type(bead, num_nt, decimals):
num_nt0 = bead.num_nt
bead.num_nt = np.around( np.float32(num_nt), decimals=decimals )
key = (bead.type_.name[0].upper(), bead.num_nt)
if key in beadtype_s:
bead.type_ = beadtype_s[key]
else:
t = deepcopy(bead.type_)
t.__dict__["nts"] = bead.num_nt*2 if t.name[0].upper() in ("D","O") else bead.num_nt
t.name = t.name + "%03d" % (t.nts*10**decimals)
print( "{} --> {} ({})".format(num_nt0, bead.num_nt, t.name) )
beadtype_s[key] = bead.type_ = t
# (cluster_size[c-1])
import scipy.cluster.hierarchy as hcluster
beads = [b for s in segments for b in s if b.type_.name[0].upper() in ("D","O")]
data = np.array([b.num_nt for b in beads])[:,np.newaxis]
clusters = hcluster.fclusterdata(data, float(max_basepairs_per_bead)/80, criterion="distance")
cluster_size = [np.mean(data[clusters == i]) for i in np.unique(clusters)]
order = int(2-np.log10(2*max_basepairs_per_bead)//1)
try:
clusters = hcluster.fclusterdata(data, float(max_basepairs_per_bead)/80, criterion="distance")
cluster_size = [np.mean(data[clusters == i]) for i in np.unique(clusters)]
except:
## TODO debug
clusters = data
cluster_size = np.arange(len(data))+1
for b,c in zip(beads,clusters):
num_nt0 = b.num_nt
b.num_nt = np.around( np.float32(cluster_size[c-1]), decimals=order )
key = (b.type_.name[0].upper(), b.num_nt)
if key in beadtype_s:
b.type_ = beadtype_s[key]
else:
t = deepcopy(b.type_)
t.__dict__["nts"] = b.num_nt*2
t.name = t.name + "%03d" % (t.nts*10**order)
print( "{} --> {} ({})".format(num_nt0, b.num_nt, t.name) )
beadtype_s[key] = b.type_ = t
_assign_bead_type(b, cluster_size[c-1], decimals=order)
beads = [b for s in segments for b in s if b.type_.name[0].upper() == "S"]
beads = [b for s in segments for b in s if b.type_.name[0].upper() in ("S")]
data = np.array([b.num_nt for b in beads])[:,np.newaxis]
clusters = hcluster.fclusterdata(data, float(max_nucleotides_per_bead)/40, criterion="distance")
cluster_size = [np.mean(data[clusters == i]) for i in np.unique(clusters)]
order = int(2-np.log10(max_nucleotides_per_bead)//1)
try:
clusters = hcluster.fclusterdata(data, float(max_basepairs_per_bead)/40, criterion="distance")
cluster_size = [np.mean(data[clusters == i]) for i in np.unique(clusters)]
except:
## TODO debug
clusters = data
cluster_size = np.arange(len(data))+1
for b,c in zip(beads,clusters):
num_nt0 = b.num_nt
b.num_nt = np.around( np.float32(cluster_size[c-1]), decimals=order )
key = (b.type_.name[0].upper(), b.num_nt)
if key in beadtype_s:
b.type_ = beadtype_s[key]
else:
t = deepcopy(b.type_)
t.__dict__["nts"] = b.num_nt
t.name = t.name + "%03d" % (t.nts*10**order)
print( "{} --> {} ({})".format(num_nt0, b.num_nt, t.name) )
beadtype_s[key] = b.type_ = t
_assign_bead_type(b, cluster_size[c-1], decimals=order)
""" Update bead indices """
self._countParticleTypes() # probably not needed here
......
Supports Markdown
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