def buildAlleles(sequence, variants, reference_start=0, phased=True):

    def _delete(allele, del_start, del_end, variant, sequence, startoffset, endoffset, feature_start, feature_end):
        '''little helper: update ``allele`` with a deletion ``del_start:del_end``.
        '''

        # truncate variant according to the feature
        variant = variant[startoffset:len(variant) - endoffset]

        n = variant.count("-")
        if n:
            if variant.startswith("-"):
                del_start += n
                variant = variant[n:]
            else:
                del_end -= n
                variant = variant[:-n]

        # due to gaps, the variant is not actually within the feauture
        if del_start >= del_end:
            return

        refseq = sequence[del_start:del_end].upper()

        assert refseq == variant, \
            'reference base mismatch at deletion: expected %s %s %s, got %s[%i:%i] at feature=%i-%i, variant=%i-%i, relative=%i-%i, del=%i-%i, action=%s' % \
            (sequence[del_start - 10:del_start],
             refseq,
             sequence[del_end:del_end + 10],
             variant, startoffset, len(variant) - endoffset,
             feature_start, feature_end,
             var_start, var_end,
             rel_start, rel_end,
             del_start, del_end,
             action)

        l = del_end - del_start

        # assert len("".join(allele[del_start:del_end])) == l, \
        #     "deletion conflicts with other indels: " \
        #     "got %s[%i:%i] (ref=%s, allele=%s) at feature=%i-%i, variant=%i-%i, relative=%i-%i, del=%i-%i, action=%s" % \
        #     (variant, startoffset, len(variant)-endoffset,
        #      refseq, str(allele[del_start:del_end]),
        #      feature_start, feature_end,
        #      var_start, var_end,
        #      rel_start, rel_end,
        #      del_start, del_end,
        #      action)

        allele[del_start:del_end] = [""] * l

    allele1 = list(sequence.lower())
    allele2 = list(sequence.lower())

    if reference_start is None:
        feature_start = 0
    else:
        feature_start = reference_start

    feature_end = feature_start + len(sequence)

    # main loop: insert variants into allele sequences
    for var_start, var_end, reference, action, has_wildtype, variantseqs in variants:

        # skip variants that are out-of-range
        if var_end <= feature_start or var_start >= feature_end:
            continue

        is_homozygous = len(variantseqs) == 1 and not has_wildtype

        rel_start, rel_end = var_start - feature_start, var_end - feature_start
        startoffset = max(0, feature_start - var_start)
        endoffset = max(0, var_end - feature_end)
        pruned_start, pruned_end = max(
            0, rel_start), min(len(sequence), rel_end)

        if action == "=":

            if E.global_options.loglevel >= 10:
                E.debug("adding SNP at postition %i: reference=%s variants=%s" %
                        (var_start, reference, variantseqs))

            if allele1[rel_start] == "" or allele2[rel_start] == "":
                # these can be cases, where a base is deleted in one allele,
                # but recorded as a homozygous substitution in another allele.
                E.warn("substitution conflicts with a deletion - ignored: %s" %
                       str((var_start, var_end, reference, action, has_wildtype, variantseqs)))
                continue

            assert rel_start >= 0
            assert sequence[rel_start].upper() == reference, \
                'reference base mismatch: expected %s %s %s, got %s at feature=%i-%i, variant=%i-%i, relative=%i-%i, pruned=%i-%i, action=%s' % \
                (sequence[rel_start - 10:rel_start],
                 sequence[rel_start].upper(),
                 sequence[rel_start + 1:rel_start + 10],
                 reference,
                 feature_start, feature_end,
                 var_start, var_end,
                 rel_start, rel_end,
                 pruned_start, pruned_end,
                 action)

            if phased:
                allele1[rel_start] = variantseqs[0] + allele1[rel_start][1:]
                allele2[rel_start] = variantseqs[1] + allele2[rel_start][1:]
            elif is_homozygous:
                allele1[rel_start] = variantseqs[0] + allele1[rel_start][1:]
                allele2[rel_start] = variantseqs[0] + allele2[rel_start][1:]
            else:
                if has_wildtype:
                    if reference == variantseqs[0]:
                        allele2[rel_start] = variantseqs[
                            1] + allele2[rel_start][1:]
                    else:
                        allele2[rel_start] = variantseqs[
                            0] + allele2[rel_start][1:]
                else:
                    allele1[rel_start] = variantseqs[
                        0] + allele1[rel_start][1:]
                    allele2[rel_start] = variantseqs[
                        1] + allele2[rel_start][1:]

        elif action == "-":
            if phased:
                _delete(allele1, pruned_start, pruned_end, variantseqs[0],
                        sequence, startoffset, endoffset, feature_start, feature_end)
                _delete(allele2, pruned_start, pruned_end, variantseqs[1],
                        sequence, startoffset, endoffset, feature_start, feature_end)
            elif is_homozygous:
                _delete(allele1, pruned_start, pruned_end, variantseqs[0],
                        sequence, startoffset, endoffset, feature_start, feature_end)
                _delete(allele2, pruned_start, pruned_end, variantseqs[0],
                        sequence, startoffset, endoffset, feature_start, feature_end)
            else:
                if has_wildtype:
                    _delete(allele2, pruned_start, pruned_end, variantseqs[0],
                            sequence, startoffset, endoffset, feature_start, feature_end)
                else:
                    _delete(allele1, pruned_start, pruned_end, variantseqs[0],
                            sequence, startoffset, endoffset, feature_start, feature_end)
                    _delete(allele2, pruned_start, pruned_end, variantseqs[1],
                            sequence, startoffset, endoffset, feature_start, feature_end)

        elif action == "+":
            # ignore insertions at position -1
            if rel_start < 0:
                continue

            if phased:
                allele1[rel_start] += variantseqs[0].upper()
                allele2[rel_start] += variantseqs[1].upper()
            elif is_homozygous:
                allele1[rel_start] += variantseqs[0].upper()
                allele2[rel_start] += variantseqs[0].upper()
            else:
                if has_wildtype:
                    allele2[rel_start] += variantseqs[0].upper()
                else:
                    allele1[rel_start] += variantseqs[0].upper()
                    allele2[rel_start] += variantseqs[1].upper()

        elif action == ">":
            # indel
            if rel_start >= 0:
                allele1[rel_start] += variantseqs[0].upper()
            _delete(allele2, pruned_start, pruned_end, variantseqs[1],
                    sequence, startoffset, endoffset, feature_start, feature_end)

        elif action == "<":
            # delin
            if rel_start >= 0:
                allele2[rel_start] += variantseqs[1].upper()
            _delete(allele1, pruned_start, pruned_end, variantseqs[0],
                    sequence, startoffset, endoffset, feature_start, feature_end)

    assert len(sequence) == len(allele1)
    assert len(sequence) == len(allele2)

    return (allele1, allele2)
Returns object definitions for object name(s)
dataset/ETHPy150Open CGATOxford/cgat/CGAT/Variants.py/buildAlleles--crossover with--dataset/ETHPy150Open ooici/pyon/pyon/container/shell_api.py/obj_defs