def evaluate(submission, truth, vtype='SNV', reffa=None, ignorechroms=None, ignorepass=False, 
             fp_vcf=None, fn_vcf=None, tp_vcf=None,
             debug=False):
    ''' return stats on sensitivity, specificity, balanced accuracy '''

    assert vtype in ('SNV', 'SV', 'INDEL')
    subvcfh = vcf.Reader(filename=submission)
    truvcfh = vcf.Reader(filename=truth)

    fpvcfh = fnvcfh = tpvcfh = None
    if fp_vcf:
        fpvcfh = vcf.Writer(open(fp_vcf, 'w'), template=subvcfh)
    if fn_vcf:
        fnvcfh = vcf.Writer(open(fn_vcf, 'w'), template=subvcfh)
    if tp_vcf:
        tpvcfh = vcf.Writer(open(tp_vcf, 'w'), template=subvcfh)

    reffa_fh = None
    if reffa:
        reffa_fh  = pysam.Fastafile(reffa)
        if debug:
            print "DEBUG: Using haplotype aware indel comparison"
    
    tpcount = 0
    fpcount = 0
    subrecs = 0
    trurecs = 0

    truchroms = {}
    fns = OrderedDict()

    ''' count records in truth vcf, track contigs/chromosomes '''
    for trurec in truvcfh:
        if relevant(trurec, vtype, ignorechroms):
            trurecs += 1
            truchroms[trurec.CHROM] = True
            fns[str(trurec)] = trurec
            
    used_truth = {} # keep track of 'truth' sites used, they should only be usable once

    ''' parse submission vcf, compare to truth '''
    for subrec in subvcfh:
        if passfilter(subrec, disabled=ignorepass):
            if subrec.is_snp and vtype == 'SNV':
                if not svmask(subrec, truvcfh, truchroms):
                    subrecs += 1
            if subrec.is_sv and vtype == 'SV':
                subrecs += 1
            if subrec.is_indel and vtype == 'INDEL':
                subrecs += 1

        matched = False

        startpos, endpos = subrec.start, subrec.end

        if vtype == 'SV' and subrec.is_sv:
            startpos, endpos = expand_sv_ends(subrec)
        try:
            if relevant(subrec, vtype, ignorechroms) and passfilter(subrec, disabled=ignorepass) and subrec.CHROM in truchroms:
                for trurec in truvcfh.fetch(subrec.CHROM, startpos, end=endpos):
                    if match(subrec, trurec, vtype=vtype) and str(trurec) not in used_truth:
                        matched = True
                if not matched and subrec.is_indel and reffa_fh:# try haplotype aware comparison
                    window = 100
                    for (trurec, _) in get_close_matches(subrec, truvcfh, window, indels_only=True):
                        if str(trurec) in used_truth:
                            continue
                        if have_identical_haplotypes(subrec, trurec, reffa_fh):
                            matched = True
                            if debug:
                                print "DEBUG: Rescuing %s which has same haplotype as %s" % (subrec, trurec)
                            break
            if matched:
                used_truth[str(trurec)] = True
                    
        except __HOLE__ as e:
            sys.stderr.write("Warning: " + str(e) + "\n")

        if matched:
            tpcount += 1
            if tpvcfh:
                tpvcfh.write_record(subrec)
            if fns.has_key(str(trurec)):
                del fns[str(trurec)]            
        else:
            if relevant(subrec, vtype, ignorechroms) and passfilter(subrec, disabled=ignorepass) and not svmask(subrec, truvcfh, truchroms): 
                fpcount += 1 # FP counting method needs to change for real tumors
                if fpvcfh:
                    fpvcfh.write_record(subrec)


    if fnvcfh:
        for fn in fns.values():
            fnvcfh.write_record(fn)


    print "tpcount, fpcount, subrecs, trurecs:"
    print tpcount, fpcount, subrecs, trurecs

    recall    = float(tpcount) / float(trurecs)
    if tpcount+fpcount > 0:
        precision = float(tpcount) / float(tpcount + fpcount)
    else:
        precision = 0.0
    #fdr       = 1.0 - float(fpcount) / float(subrecs)
    f1score   = 0.0 if tpcount == 0 else 2.0*(precision*recall)/(precision+recall)

    for fh in [fpvcfh, fnvcfh, tpvcfh]:
        if fh:
            fh.close()
            
    return precision, recall, f1score
dataset/ETHPy150Open adamewing/bamsurgeon/scripts/evaluator.py/evaluate