def have_identical_haplotypes(v1, v2, ref):
    """Check if two variant produce the same haplotype / variant sequence.

    - v1 and v2: PyVCF variants to compare
    - ref: PySAM FastaFile
    """

    assert (v1.is_indel or v1.is_snp) and (v2.is_indel or v2.is_snp)

    if v1.CHROM != v2.CHROM:
        return False

    if v1.is_snp and v2.is_snp:
        assert v1.REF.upper() == v2.REF.upper()
        return str(v1.ALT[0]).upper() == str(v2.ALT[0]).upper()
    if v1.is_snp or v2.is_snp:
        # one snp one indel: can't produce identical results
        return False

    assert v1.is_indel and v2.is_indel
    # only on allele per variant allowed
    assert len(v1.ALT) == 1 and len(v2.ALT) == 1, (
    +        "Can't handle multi-allelic entries")

    # get the sequence context whic fully overlaps both variants.
    # note: pyvcf is one-based, but start and end are zero-based half-open
    start = min([v1.POS, v2.POS])-1
    end = max([v1.POS + max([len(v1.REF), len(v1.ALT[0])]),
               v2.POS + max([len(v2.REF), len(v2.ALT[0])])
               ])
    chrom = v1.CHROM# made sure before they are identical before
    seq = list(ref.fetch(chrom, start, end).upper())

    if len(seq) != end-start:
        # FIXME how to handle?
        sys.stderr.write("WARN: Couldn't fetch full sequence window. Skipping"
                     " allele-aware comparison, otherwise indices would"
                     " be off\n")
        raise NotImplementedError

    v1_offset = v1.POS-1-start
    v2_offset = v2.POS-1-start
    # lower() in replacement for debugging purposes only
    v1_seq = seq[:v1_offset] + list(str(v1.ALT[0]).lower()) + seq[v1_offset+len(v1.REF):]
    v2_seq = seq[:v2_offset] + list(str(v2.ALT[0]).lower()) + seq[v2_offset+len(v2.REF):]
    if False:
        print "reference sequence context\t%s" % (''.join(seq))
        print "v1 (offset %d) %s\t%s" % (v1_offset, v1, ''.join(v1_seq))
        print "v2 (offset %d) %s\t%s" % (v2_offset, v2, ''.join(v2_seq))
        print

    try:
        assert seq[v1_offset] == v1.REF[0].upper()
        assert seq[v2_offset] == v2.REF[0].upper()
        assert len(v1_seq) == len(seq) - len(v1.REF) + len(v1.ALT[0])
        assert len(v2_seq) == len(seq) - len(v2.REF) + len(v2.ALT[0])
    except __HOLE__:
        #import pdb; pdb.set_trace()
        raise

    #if ''.join(v1_seq).upper() == ''.join(v2_seq).upper():
    #    print ''.join(v1_seq).upper()
    return ''.join(v1_seq).upper() == ''.join(v2_seq).upper()
dataset/ETHPy150Open adamewing/bamsurgeon/scripts/evaluator.py/have_identical_haplotypes