Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
help="Create edge within certain distance [default: %default]")
p.set_verbose(help="Print verbose reports to stdout")
opts, args = p.parse_args(args)
if len(args) != 2:
sys.exit(not p.print_help())
blastfile, subjectfasta = args
clique = opts.clique
maxdist = opts.maxdist
sort([blastfile, "--query"])
blast = BlastSlow(blastfile, sorted=True)
g = BiGraph()
for query, blines in groupby(blast, key=lambda x: x.query):
blines = list(blines)
iterator = combinations(blines, 2) if clique else pairwise(blines)
for a, b in iterator:
asub, bsub = a.subject, b.subject
if asub == bsub:
continue
arange = (a.query, a.qstart, a.qstop, "+")
brange = (b.query, b.qstart, b.qstop, "+")
dist, oo = range_distance(arange, brange, distmode="ee")
if dist > maxdist:
continue
atag = ">" if a.orientation == "+" else "<"
btag = ">" if b.orientation == "+" else "<"
g.add_edge(asub, bsub, atag, btag)
graph_to_agp(g, blastfile, subjectfasta, verbose=opts.verbose)
def happy_edges(row, prefix=None):
"""
Convert a row in HAPPY file and yield edges.
"""
trans = maketrans("[](){}", " ")
row = row.strip().strip("+")
row = row.translate(trans)
scfs = [x.strip("+") for x in row.split(":")]
for a, b in pairwise(scfs):
oa = '<' if a.strip()[0] == '-' else '>'
ob = '<' if b.strip()[0] == '-' else '>'
is_uncertain = a[-1] == ' ' or b[0] == ' '
a = a.strip().strip('-')
b = b.strip().strip('-')
if prefix:
a = prefix + a
b = prefix + b
yield (a, b, oa, ob), is_uncertain
frags = base + ".frags.fasta"
pairs = base + ".pairs.fasta"
if fastafile.endswith(".gz"):
frags += ".gz"
pairs += ".gz"
fragsfw = must_open(frags, "w")
pairsfw = must_open(pairs, "w")
N = opts.rclip
strip_name = lambda x: x[:-N] if N else str
skipflag = False # controls the iterator skip
fastaiter = SeqIO.parse(fastafile, "fasta")
for a, b in pairwise(fastaiter):
aid, bid = [strip_name(x) for x in (a.id, b.id)]
if skipflag:
skipflag = False
continue
if aid == bid:
SeqIO.write([a, b], pairsfw, "fasta")
skipflag = True
else:
SeqIO.write([a], fragsfw, "fasta")
# don't forget the last one, when b is None
if not skipflag:
SeqIO.write([a], fragsfw, "fasta")
def path_to_agp(g, path, object, sizes, status):
lines = []
for (a, ao), (b, bo) in pairwise(path):
ao = get_orientation(ao, status)
e = g.get_edge(a.v, b.v)
cline = AGPLine.cline(object, a.v, sizes, ao)
gline = AGPLine.gline(object, e.length)
lines.append(cline)
lines.append(gline)
# Do not forget the last one
z, zo = path[-1]
zo = get_orientation(zo, status)
cline = AGPLine.cline(object, z.v, sizes, zo)
lines.append(cline)
return lines
sol = model.cbGetSolution([model._vars[i] for i in range(nedges)])
selected = [all_edges[i] for i, x in enumerate(sol) if x > .5]
# find the shortest cycle in the selected edge list
current_nedges = 0
for salesman in salesmen:
incoming, outgoing, nodes = node_to_edge(edges)
edge_store = dict(((a, b), i + current_nedges) \
for i, (a, b, w) in salesman)
salesmen_selected = [(idx[a], idx[b]) for a, b, w in selected \
if a in nodes and b in nodes]
tour = subtour(salesmen_selected)
if len(tour) == len(nodes):
return
# add a subtour elimination constraint
c = tour
incident = [edge_store[a, b] for a, b in pairwise(c + [c[0]])]
model.cbLazy(quicksum(model._vars[x] for x in incident) <= len(tour) - 1)
current_nedges += len(salesman)
def subtourelim(model, where):
if where != GRB.callback.MIPSOL:
return
selected = []
# make a list of edges selected in the solution
sol = model.cbGetSolution([model._vars[i] for i in range(nedges)])
selected = [edges[i] for i, x in enumerate(sol) if x > .5]
selected = [(idx[a], idx[b]) for a, b, w in selected]
# find the shortest cycle in the selected edge list
tour = subtour(selected)
if len(tour) == n:
return
# add a subtour elimination constraint
c = tour
incident = [edge_store[a, b] for a, b in pairwise(c + [c[0]])]
model.cbLazy(quicksum(model._vars[x]
for x in incident) <= len(tour) - 1)
start, end = b.sstart, b.sstop
cstart, cend = min(size, clip), max(0, size - clip)
if start > cstart and end < cend:
continue
blasts.append(b)
key = lambda x: x.query
blasts.sort(key=key)
g = BiGraph()
for query, bb in groupby(blasts, key=key):
bb = sorted(bb, key=lambda x: x.qstart)
nsubjects = len(set(x.subject for x in bb))
if nsubjects == 1:
continue
print("\n".join(str(x) for x in bb))
for a, b in pairwise(bb):
astart, astop = a.qstart, a.qstop
bstart, bstop = b.qstart, b.qstop
if a.subject == b.subject:
continue
arange = astart, astop
brange = bstart, bstop
ov = range_intersect(arange, brange)
alen = astop - astart + 1
blen = bstop - bstart + 1
if ov:
ostart, ostop = ov
ov = ostop - ostart + 1
print(ov, alen, blen)
if ov and (ov > alen / 2 or ov > blen / 2):
r2 = Rectangle((c, yb - ytip), d - c, 2 * ytip, fc='r', lw=0, zorder=2)
r3 = Rectangle((a, ya - ytip), b - a, 2 * ytip, fill=False, zorder=3)
r4 = Rectangle((c, yb - ytip), d - c, 2 * ytip, fill=False, zorder=3)
r5 = Polygon(((a, ya - ytip), (c, yb + ytip),
(d, yb + ytip), (b, ya - ytip)),
fc='r', alpha=.2)
rr = (r1, r2, r3, r4, r5)
if i == 2:
rr = rr[:-1]
for r in rr:
root.add_patch(r)
# Gap pairs
hspa, hspb = zip(*myhsps)
gapa, gapb = [], []
for (a, b), (c, d) in pairwise(hspa):
gapa.append((b + 1, c - 1))
for (a, b), (c, d) in pairwise(hspb):
gapb.append((b + 1, c - 1))
gaps = zip(gapa, gapb)
tpos = titlepos[-1]
yy = tpos - .05
for i, ((a, b), (c, d)) in enumerate(gaps):
i += 1
a, b, c, d = [m(x) for x in (a, b, c, d)]
xx = (a + b + c + d) / 4
TextCircle(root, xx, yy, str(i))
# Bites
ystart = .24
ytip = .05
def graph(self):
from jcvi.algorithms.graph import BiGraph
g = BiGraph()
for ob, lines in self.iter_object():
components = [x for x in lines if not x.is_gap]
gaps = [x for x in lines if x.is_gap]
for i, (a, b) in enumerate(pairwise(components)):
g.add_edge(a.component_id, b.component_id,
a.orientation, b.orientation,
length=gaps[i].gap_length)
if len(components) == 1: # Singleton object
a = components[0]
g.add_node(a.component_id)
return g