How to use the jcvi.utils.iter.pairwise function in jcvi

To help you get started, we’ve selected a few jcvi examples, based on popular ways it is used in public projects.

Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.

github tanghaibao / jcvi / jcvi / assembly / syntenypath.py View on Github external
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)
github tanghaibao / jcvi / jcvi / assembly / syntenypath.py View on Github external
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
github tanghaibao / jcvi / jcvi / formats / fasta.py View on Github external
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")
github tanghaibao / jcvi / jcvi / assembly / sspace.py View on Github external
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
github tanghaibao / jcvi / algorithms / synctsp.py View on Github external
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)
github tanghaibao / jcvi / jcvi / algorithms / lpsolve.py View on Github external
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)
github tanghaibao / jcvi / jcvi / assembly / syntenypath.py View on Github external
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):
github tanghaibao / jcvi / jcvi / projects / bites.py View on Github external
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
github tanghaibao / jcvi / jcvi / formats / agp.py View on Github external
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