annotate utilities/vcf_compare_OLD.py @ 2:8a739c944dbf draft

planemo upload commit e96b43f96afce6a7b7dfd4499933aad7d05c955e-dirty

1 #!/usr/bin/env python
2 # encoding: utf-8
3 
4 """ **************************************************
5 
6 vcf_compare.py
7 
8 - compare vcf file produced by workflow to golden vcf produced by simulator
9 
10 Written by:		Zach Stephens
11 Date:			January 20, 2015
12 Contact:		[e-mail]
13 
14 ************************************************** """
15 
16 import sys
17 import os
18 import copy
19 import time
20 import bisect
21 import re
22 import numpy as np
23 import optparse
24 
25 
26 EV_BPRANGE = 50		# how far to either side of a particular variant location do we want to check for equivalents?
27 
28 DEFAULT_QUAL = -666	# if we can't find a qual score, use this instead so we know it's missing
29 
30 MAX_VAL = 9999999999999	# an unreasonably large value that no reference fasta could concievably be longer than
31 
32 DESC   = """%prog: vcf comparison script."""
33 VERS   = 0.1
34 
35 PARSER = optparse.OptionParser('python %prog [options] -r <ref.fa> -g <golden.vcf> -w <workflow.vcf>',description=DESC,version="%prog v"+str(VERS))
36 
37 PARSER.add_option('-r', help='* Reference Fasta', dest='REFF', action='store', metavar='<ref.fa>')
38 PARSER.add_option('-g', help='* Golden VCF',      dest='GVCF', action='store', metavar='<golden.vcf>')
39 PARSER.add_option('-w', help='* Workflow VCF',    dest='WVCF', action='store', metavar='<workflow.vcf>')
40 PARSER.add_option('-o', help='* Output Prefix',   dest='OUTF', action='store', metavar='<prefix>')
41 PARSER.add_option('-m', help='Mappability Track', dest='MTRK', action='store', metavar='<track.bed>')
42 PARSER.add_option('-M', help='Maptrack Min Len',  dest='MTMM', action='store', metavar='<int>')
43 PARSER.add_option('-t', help='Targetted Regions', dest='TREG', action='store', metavar='<regions.bed>')
44 PARSER.add_option('-T', help='Min Region Len',    dest='MTRL', action='store', metavar='<int>')
45 PARSER.add_option('-c', help='Coverage Filter Threshold [%default]',       dest='DP_THRESH', default=15, action='store', metavar='<int>')
46 PARSER.add_option('-a', help='Allele Freq Filter Threshold [%default]',    dest='AF_THRESH', default=0.3, action='store', metavar='<float>')
47 
48 PARSER.add_option('--vcf-out',   help="Output Match/FN/FP variants [%default]",       dest='VCF_OUT', default=False, action='store_true')
49 PARSER.add_option('--no-plot',   help="No plotting [%default]",                       dest='NO_PLOT', default=False, action='store_true')
50 PARSER.add_option('--incl-homs', help="Include homozygous ref calls [%default]",      dest='INCL_H',  default=False, action='store_true')
51 PARSER.add_option('--incl-fail', help="Include calls that failed filters [%default]", dest='INCL_F',  default=False, action='store_true')
52 PARSER.add_option('--fast',      help="No equivalent variant detection [%default]",   dest='FAST',    default=False, action='store_true')
53 
54 (OPTS,ARGS) = PARSER.parse_args()
55 
56 REFERENCE    = OPTS.REFF
57 GOLDEN_VCF   = OPTS.GVCF
58 WORKFLOW_VCF = OPTS.WVCF
59 OUT_PREFIX   = OPTS.OUTF
60 MAPTRACK     = OPTS.MTRK
61 MIN_READLEN  = OPTS.MTMM
62 BEDFILE      = OPTS.TREG
63 DP_THRESH    = int(OPTS.DP_THRESH)
64 AF_THRESH    = float(OPTS.AF_THRESH)
65 
66 VCF_OUT      = OPTS.VCF_OUT
67 NO_PLOT      = OPTS.NO_PLOT
68 INCLUDE_HOMS = OPTS.INCL_H
69 INCLUDE_FAIL = OPTS.INCL_F
70 FAST         = OPTS.FAST
71 
72 if len(sys.argv[1:]) == 0:
73 	PARSER.print_help()
74 	exit(1)
75 
76 if OPTS.MTRL != None:
77 	MINREGIONLEN = int(OPTS.MTRL)
78 else:
79 	MINREGIONLEN = None
80 
81 if MIN_READLEN == None:
82 	MIN_READLEN = 0
83 else:
84 	MIN_READLEN = int(MIN_READLEN)
85 
86 if REFERENCE == None:
87 	print 'Error: No reference provided.'
88 	exit(1)
89 if GOLDEN_VCF == None:
90 	print 'Error: No golden VCF provided.'
91 	exit(1)
92 if WORKFLOW_VCF == None:
93 	print 'Error: No workflow VCF provided.'
94 	exit(1)
95 if OUT_PREFIX == None:
96 	print 'Error: No output prefix provided.'
97 	exit(1)
98 if (BEDFILE != None and MINREGIONLEN == None) or (BEDFILE == None and MINREGIONLEN != None):
99 	print 'Error: Both -t and -T must be specified'
00 	exit(1)
01 
02 if NO_PLOT == False:
03 	import matplotlib
04 	matplotlib.use('Agg')
05 	import matplotlib.pyplot as mpl
06 	from matplotlib_venn import venn2, venn3
07 	import warnings
08 	warnings.filterwarnings("ignore", category=UserWarning, module='matplotlib_venn')
09 
10 AF_STEPS = 20
11 AF_KEYS  = np.linspace(0.0,1.0,AF_STEPS+1)
12 
13 def quantize_AF(af):
14 	if af >= 1.0:
15 		return AF_STEPS
16 	elif af <= 0.0:
17 		return 0
18 	else:
19 		return int(af*AF_STEPS)
20 
21 VCF_HEADER = '##fileformat=VCFv4.1\n##reference='+REFERENCE+'##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">\n##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">\n'
22 
23 DP_TOKENS = ['DP','DPU','DPI']	# in the order that we'll look for them
24 
25 def parseLine(splt,colDict,colSamp):
26 
27 	#	check if we want to proceed..
28 	ra = splt[colDict['REF']]
29 	aa = splt[colDict['ALT']]
30 	if not(INCLUDE_HOMS) and (aa == '.' or aa == '' or aa == ra):
31 		return None
32 	if not(INCLUDE_FAIL) and (splt[colDict['FILTER']] != 'PASS' and splt[colDict['FILTER']] != '.'):
33 		return None
34 
35 	#	default vals
36 	cov  = None
37 	qual = DEFAULT_QUAL
38 	alt_alleles = []
39 	alt_freqs   = [None]
40 
41 	#	any alt alleles?
42 	alt_split = aa.split(',')
43 	if len(alt_split) > 1:
44 		alt_alleles = alt_split
45 
46 	#	cov
47 	for dp_tok in DP_TOKENS:
48 		#	check INFO for DP first
49 		if 'INFO' in colDict and dp_tok+'=' in splt[colDict['INFO']]:
50 			cov = int(re.findall(re.escape(dp_tok)+r"=[0-9]+",splt[colDict['INFO']])[0][3:])
51 		#	check FORMAT/SAMPLE for DP second:
52 		elif 'FORMAT' in colDict and len(colSamp):
53 			format = splt[colDict['FORMAT']]+':'
54 			if ':'+dp_tok+':' in format:
55 				dpInd = format.split(':').index(dp_tok)
56 				cov   = int(splt[colSamp[0]].split(':')[dpInd])
57 		if cov != None:
58 			break
59 
60 	#	check INFO for AF first
61 	af = None
62 	if 'INFO' in colDict and ';AF=' in ';'+splt[colDict['INFO']]:
63 		info = splt[colDict['INFO']]+';'
64 		af  = re.findall(r"AF=.*?(?=;)",info)[0][3:]
65 	#	check FORMAT/SAMPLE for AF second:
66 	elif 'FORMAT' in colDict and len(colSamp):
67 		format = splt[colDict['FORMAT']]+':'
68 		if ':AF:' in format:
69 			afInd = splt[colDict['FORMAT']].split(':').index('AF')
70 			af    = splt[colSamp[0]].split(':')[afInd]
71 
72 	if af != None:
73 		af_splt = af.split(',')
74 		while(len(af_splt) < len(alt_alleles)):	# are we lacking enough AF values for some reason?
75 			af_splt.append(af_splt[-1])			# phone it in.
76 		if len(af_splt) != 0 and af_splt[0] != '.' and af_splt[0] != '':		# missing data, yay
77 			alt_freqs = [float(n) for n in af_splt]
78 	else:
79 		alt_freqs = [None]*max([len(alt_alleles),1])
80 
81 	#	get QUAL if it's interesting
82 	if 'QUAL' in colDict and splt[colDict['QUAL']] != '.':
83 		qual = float(splt[colDict['QUAL']])
84 
85 	return (cov, qual, alt_alleles, alt_freqs)
86 
87 
88 def parseVCF(VCF_FILENAME,refName,targRegionsFl,outFile,outBool):
89 	v_Hashed      = {}
90 	v_posHash     = {}
91 	v_Alts        = {}
92 	v_Cov         = {}
93 	v_AF          = {}
94 	v_Qual        = {}
95 	v_TargLen     = {}
96 	nBelowMinRLen = 0
97 	line_unique   = 0	# number of lines in vcf file containing unique variant
98 	hash_coll     = 0	# number of times we saw a hash collision ("per line" so non-unique alt alleles don't get counted multiple times)
99 	var_filtered  = 0	# number of variants excluded due to filters (e.g. hom-refs, qual)
00 	var_merged    = 0	# number of variants we merged into another due to having the same position specified
01 	colDict       = {}
02 	colSamp       = []
03 	for line in open(VCF_FILENAME,'r'):
04 		if line[0] != '#':
05 			if len(colDict) == 0:
06 				print '\n\nError: VCF has no header?\n'+VCF_FILENAME+'\n\n'
07 				exit(1)
08 			splt = line[:-1].split('\t')
09 			if splt[0] == refName:
10 
11 				var  = (int(splt[1]),splt[3],splt[4])
12 				targInd = bisect.bisect(targRegionsFl,var[0])
13 
14 				if targInd%2 == 1:
15 					targLen = targRegionsFl[targInd]-targRegionsFl[targInd-1]
16 					if (BEDFILE != None and targLen >= MINREGIONLEN) or BEDFILE == None:
17 						
18 						pl_out = parseLine(splt,colDict,colSamp)
19 						if pl_out == None:
20 							var_filtered += 1
21 							continue
22 						(cov, qual, aa, af) = pl_out
23 
24 						if var not in v_Hashed:
25 
26 							vpos = var[0]
27 							if vpos in v_posHash:
28 								if len(aa) == 0:
29 									aa = [var[2]]
30 								aa.extend([n[2] for n in v_Hashed.keys() if n[0] == vpos])
31 								var_merged += 1
32 							v_posHash[vpos] = 1
33 							
34 							if len(aa):
35 								allVars = [(var[0],var[1],n) for n in aa]
36 								for i in xrange(len(allVars)):
37 									v_Hashed[allVars[i]] = 1
38 									#if allVars[i] not in v_Alts:
39 									#	v_Alts[allVars[i]] = []
40 									#v_Alts[allVars[i]].extend(allVars)
41 									v_Alts[allVars[i]] = allVars
42 							else:
43 								v_Hashed[var] = 1
44 
45 							if cov != None:
46 								v_Cov[var] = cov
47 							v_AF[var]      = af[0]		# only use first AF, even if multiple. fix this later?
48 							v_Qual[var]    = qual
49 							v_TargLen[var] = targLen
50 							line_unique += 1
51 
52 						else:
53 							hash_coll += 1
54 
55 					else:
56 						nBelowMinRLen += 1
57 		else:
58 			if line[1] != '#':
59 				cols = line[1:-1].split('\t')
60 				for i in xrange(len(cols)):
61 					if 'FORMAT' in colDict:
62 						colSamp.append(i)
63 					colDict[cols[i]] = i
64 				if VCF_OUT and outBool:
65 					outBool = False
66 					outFile.write(line)
67 
68 	return (v_Hashed, v_Alts, v_Cov, v_AF, v_Qual, v_TargLen, nBelowMinRLen, line_unique, var_filtered, var_merged, hash_coll)
69 
70 
71 def condenseByPos(listIn):
72 	varListOfInterest = [n for n in listIn]
73 	indCount = {}
74 	for n in varListOfInterest:
75 		c = n[0]
76 		if c not in indCount:
77 			indCount[c] = 0
78 		indCount[c] += 1
79 	#nonUniqueDict = {n:[] for n in sorted(indCount.keys()) if indCount[n] > 1}		# the python 2.7 way
80 	nonUniqueDict = {}
81 	for n in sorted(indCount.keys()):
82 		if indCount[n] > 1:
83 			nonUniqueDict[n] = []
84 	delList = []
85 	for i in xrange(len(varListOfInterest)):
86 		if varListOfInterest[i][0] in nonUniqueDict:
87 			nonUniqueDict[varListOfInterest[i][0]].append(varListOfInterest[i])
88 			delList.append(i)
89 	delList = sorted(delList,reverse=True)
90 	for di in delList:
91 		del varListOfInterest[di]
92 	for v in nonUniqueDict.values():
93 		var = (v[0][0],v[0][1],','.join([n[2] for n in v[::-1]]))
94 		varListOfInterest.append(var)
95 	return varListOfInterest
96 
97 
98 def main():
99 
00 	ref = []
01 	f = open(REFERENCE,'r')
02 	nLines = 0
03 	prevR = None
04 	prevP = None
05 	ref_inds = []
06 	sys.stdout.write('\nindexing reference fasta... ')
07 	sys.stdout.flush()
08 	tt = time.time()
09 	while 1:
10 		nLines += 1
11 		data = f.readline()
12 		if not data:
13 			ref_inds.append( (prevR, prevP, f.tell()-len(data)) )
14 			break
15 		if data[0] == '>':
16 			if prevP != None:
17 				ref_inds.append( (prevR, prevP, f.tell()-len(data)) )
18 			prevP = f.tell()
19 			prevR = data[1:-1]
20 	print '{0:.3f} (sec)'.format(time.time()-tt)
21 	#ref_inds = [('chrM', 6, 16909), ('chr1', 16915, 254252549), ('chr2', 254252555, 502315916), ('chr3', 502315922, 704298801), ('chr4', 704298807, 899276169), ('chr5', 899276175, 1083809741), ('chr6', 1083809747, 1258347116), ('chr7', 1258347122, 1420668559), ('chr8', 1420668565, 1569959868), ('chr9', 1569959874, 1713997574), ('chr10', 1713997581, 1852243023), ('chr11', 1852243030, 1989949677), ('chr12', 1989949684, 2126478617), ('chr13', 2126478624, 2243951900), ('chr14', 2243951907, 2353448438), ('chr15', 2353448445, 2458030465), ('chr16', 2458030472, 2550192321), ('chr17', 2550192328, 2633011443), ('chr18', 2633011450, 2712650243), ('chr19', 2712650250, 2772961813), ('chr20', 2772961820, 2837247851), ('chr21', 2837247858, 2886340351), ('chr22', 2886340358, 2938671016), ('chrX', 2938671022, 3097046994), ('chrY', 3097047000, 3157608038)]
22 
23 	ztV = 0	# total golden variants
24 	ztW = 0	# total workflow variants
25 	znP = 0	# total perfect matches
26 	zfP = 0	# total false positives
27 	znF = 0	# total false negatives
28 	znE = 0	# total equivalent variants detected
29 	zgF = 0	# total golden variants that were filtered and excluded
30 	zgR = 0	# total golden variants that were excluded for being redundant
31 	zgM = 0	# total golden variants that were merged into a single position
32 	zwF = 0	# total workflow variants that were filtered and excluded
33 	zwR = 0	# total workflow variants that were excluded for being redundant
34 	zwM = 0	# total workflow variants that were merged into a single position
35 	if BEDFILE != None:
36 		zbM = 0
37 
38 	mappability_vs_FN = {0:0, 1:0}	# [0] = # of FNs that were in mappable regions, [1] = # of FNs that were in unmappable regions
39 	coverage_vs_FN    = {}			# [C] = # of FNs that were covered by C reads
40 	alleleBal_vs_FN   = {}			# [AF] = # of FNs that were heterozygous genotypes with allele freq AF (quantized to multiples of 1/AF_STEPS)
41 	for n in AF_KEYS:
42 		alleleBal_vs_FN[n] = 0
43 
44 	#
45 	#	read in mappability track
46 	#
47 	mappability_tracks = {}		# indexed by chr string (e.g. 'chr1'), has boolean array
48 	prevRef            = ''
49 	relevantRegions    = []
50 	if MAPTRACK != None:
51 		mtf = open(MAPTRACK,'r')
52 		for line in mtf:
53 			splt = line.strip().split('\t')
54 			if prevRef != '' and splt[0] != prevRef:
55 				# do stuff
56 				if len(relevantRegions):
57 					myTrack = [0]*(relevantRegions[-1][1]+100)
58 					for r in relevantRegions:
59 						for ri in xrange(r[0],r[1]):
60 							myTrack[ri] = 1
61 					mappability_tracks[prevRef] = [n for n in myTrack]
62 				#
63 				relevantRegions = []
64 			if int(splt[3]) >= MIN_READLEN:
65 				relevantRegions.append((int(splt[1]),int(splt[2])))
66 			prevRef = splt[0]
67 		mtf.close()
68 		# do stuff
69 		if len(relevantRegions):
70 			myTrack = [0]*(relevantRegions[-1][1]+100)
71 			for r in relevantRegions:
72 				for ri in xrange(r[0],r[1]):
73 					myTrack[ri] = 1
74 			mappability_tracks[prevRef] = [n for n in myTrack]
75 		#
76 
77 	#
78 	#	init vcf output, if desired
79 	#
80 	vcfo2 = None
81 	vcfo3 = None
82 	global vcfo2_firstTime
83 	global vcfo3_firstTime
84 	vcfo2_firstTime = False
85 	vcfo3_firstTime = False
86 	if VCF_OUT:
87 		vcfo2 = open(OUT_PREFIX+'_FN.vcf','w')
88 		vcfo3 = open(OUT_PREFIX+'_FP.vcf','w')
89 		vcfo2_firstTime = True
90 		vcfo3_firstTime = True
91 
92 	#
93 	#	data for plotting FN analysis
94 	#
95 	set1 = []
96 	set2 = []
97 	set3 = []
98 	varAdj = 0
99 
00 	#
01 	#
02 	#	For each sequence in reference fasta...
03 	#
04 	#
05 	for n_RI in ref_inds:
06 
07 		refName = n_RI[0]
08 		if FAST == False:
09 			f.seek(n_RI[1])
10 			print 'reading '+refName+'...',
11 			myDat  = f.read(n_RI[2]-n_RI[1]).split('\n')
12 			myLen  = sum([len(m) for m in myDat])
13 			if sys.version_info >= (2,7):
14 				print '{:,} bp'.format(myLen)
15 			else:
16 				print '{0:} bp'.format(myLen)
17 			inWidth = len(myDat[0])
18 			if len(myDat[-1]) == 0:	# if last line is empty, remove it.
19 				del myDat[-1]
20 			if inWidth*(len(myDat)-1)+len(myDat[-1]) != myLen:
21 				print 'fasta column-width not consistent.'
22 				print myLen, inWidth*(len(myDat)-1)+len(myDat[-1])
23 				for i in xrange(len(myDat)):
24 					if len(myDat[i]) != inWidth:
25 						print i, len(myDat[i]), inWidth
26 				exit(1)
27 
28 			myDat = bytearray(''.join(myDat)).upper()
29 			myLen = len(myDat)
30 
31 		#
32 		#	Parse relevant targeted regions
33 		#
34 		targRegionsFl = []
35 		if BEDFILE != None:
36 			bedfile = open(BEDFILE,'r')
37 			for line in bedfile:
38 				splt = line.split('\t')
39 				if splt[0] == refName:
40 					targRegionsFl.extend((int(splt[1]),int(splt[2])))
41 			bedfile.close()
42 		else:
43 			targRegionsFl = [-1,MAX_VAL+1]
44 
45 		#
46 		#	Parse vcf files
47 		#
48 		sys.stdout.write('comparing variation in '+refName+'... ')
49 		sys.stdout.flush()
50 		tt = time.time()
51 
52 		(correctHashed, correctAlts, correctCov, correctAF, correctQual, correctTargLen, correctBelowMinRLen, correctUnique, gFiltered, gMerged, gRedundant)        = parseVCF(GOLDEN_VCF, refName, targRegionsFl, vcfo2, vcfo2_firstTime)
53 		(workflowHashed, workflowAlts, workflowCov, workflowAF, workflowQual, workflowTarLen, workflowBelowMinRLen, workflowUnique, wFiltered, wMerged, wRedundant) = parseVCF(WORKFLOW_VCF, refName, targRegionsFl, vcfo3, vcfo3_firstTime)
54 		zgF += gFiltered
55 		zgR += gRedundant
56 		zgM += gMerged
57 		zwF += wFiltered
58 		zwR += wRedundant
59 		zwM += wMerged
60 
61 		#
62 		#	Deduce which variants are FP / FN
63 		#
64 		solvedInds = {}
65 		for var in correctHashed.keys():
66 			if var in workflowHashed or var[0] in solvedInds:
67 				correctHashed[var]  = 2
68 				workflowHashed[var] = 2
69 				solvedInds[var[0]] = True
70 		for var in correctHashed.keys()+workflowHashed.keys():
71 			if var[0] in solvedInds:
72 				correctHashed[var]  = 2
73 				workflowHashed[var] = 2
74 		nPerfect = len(solvedInds)
75 		
76 		# correctHashed[var] = 1: were not found
77 		#                    = 2: should be discluded because we were found
78 		#                    = 3: should be discluded because an alt was found
79 		notFound   = [n for n in sorted(correctHashed.keys()) if correctHashed[n] == 1]
80 		FPvariants = [n for n in sorted(workflowHashed.keys()) if workflowHashed[n] == 1]
81 
82 		#
83 		#	condense all variants who have alternate alleles and were *not* found to have perfect matches
84 		#	into a single variant again. These will not be included in the candidates for equivalency checking. Sorry!
85 		#
86 		notFound   = condenseByPos(notFound)
87 		FPvariants = condenseByPos(FPvariants)
88 
89 		#
90 		#	tally up some values, if there are no golden variants lets save some CPU cycles and move to the next ref
91 		#
92 		totalGoldenVariants   = nPerfect + len(notFound)
93 		totalWorkflowVariants = nPerfect + len(FPvariants)
94 		if totalGoldenVariants == 0:
95 			zfP += len(FPvariants)
96 			ztW += totalWorkflowVariants
97 			print '{0:.3f} (sec)'.format(time.time()-tt)
98 			continue
99 
00 		#
01 		#	let's check for equivalent variants
02 		#
03 		if FAST == False:
04 			delList_i = []
05 			delList_j = []
06 			regionsToCheck = []
07 			for i in xrange(len(FPvariants)):
08 				pos = FPvariants[i][0]
09 				regionsToCheck.append((max([pos-EV_BPRANGE-1,0]),min([pos+EV_BPRANGE,len(myDat)-1])))
10 
11 			for n in regionsToCheck:
12 				refSection = myDat[n[0]:n[1]]
13 
14 				fpWithin = []
15 				for i in xrange(len(FPvariants)):
16 					m  = FPvariants[i]
17 					if (m[0] > n[0] and m[0] < n[1]):
18 						fpWithin.append((m,i))
19 				fpWithin = sorted(fpWithin)
20 				adj = 0
21 				altSection = copy.deepcopy(refSection)
22 				for (m,i) in fpWithin:
23 					lr = len(m[1])
24 					la = len(m[2])
25 					dpos = m[0]-n[0]+adj
26 					altSection = altSection[:dpos-1] + m[2] + altSection[dpos-1+lr:]
27 					adj += la-lr
28 
29 				nfWithin = []
30 				for j in xrange(len(notFound)):
31 					m = notFound[j]
32 					if (m[0] > n[0] and m[0] < n[1]):
33 						nfWithin.append((m,j))
34 				nfWithin = sorted(nfWithin)
35 				adj = 0
36 				altSection2 = copy.deepcopy(refSection)
37 				for (m,j) in nfWithin:
38 					lr = len(m[1])
39 					la = len(m[2])
40 					dpos = m[0]-n[0]+adj
41 					altSection2 = altSection2[:dpos-1] + m[2] + altSection2[dpos-1+lr:]
42 					adj += la-lr
43 
44 				if altSection == altSection2:
45 					for (m,i) in fpWithin:
46 						if i not in delList_i:
47 							delList_i.append(i)
48 					for (m,j) in nfWithin:
49 						if j not in delList_j:
50 							delList_j.append(j)
51 
52 			nEquiv = 0
53 			for i in sorted(list(set(delList_i)),reverse=True):
54 				del FPvariants[i]
55 			for j in sorted(list(set(delList_j)),reverse=True):
56 				del notFound[j]
57 				nEquiv += 1
58 			nPerfect += nEquiv
59 
60 		#
61 		#	Tally up errors and whatnot
62 		#
63 		ztV += totalGoldenVariants
64 		ztW += totalWorkflowVariants
65 		znP += nPerfect
66 		zfP += len(FPvariants)
67 		znF += len(notFound)
68 		if FAST == False:
69 			znE += nEquiv
70 		if BEDFILE != None:
71 			zbM += correctBelowMinRLen
72 
73 		#
74 		#	try to identify a reason for FN variants:
75 		#
76 
77 		venn_data = [[0,0,0] for n in notFound]		# [i] = (unmappable, low cov, low het)
78 		for i in xrange(len(notFound)):
79 			var = notFound[i]
80 
81 			noReason = True
82 
83 			#	mappability?
84 			if MAPTRACK != None:
85 				if refName in mappability_tracks and var[0] < len(mappability_tracks[refName]):
86 					if mappability_tracks[refName][var[0]]:
87 						mappability_vs_FN[1] += 1
88 						venn_data[i][0] = 1
89 						noReason = False
90 					else:
91 						mappability_vs_FN[0] += 1
92 
93 			#	coverage?
94 			if var in correctCov:
95 				c = correctCov[var]
96 				if c != None:
97 					if c not in coverage_vs_FN:
98 						coverage_vs_FN[c] = 0
99 					coverage_vs_FN[c] += 1
00 					if c < DP_THRESH:
01 						venn_data[i][1] = 1
02 						noReason = False
03 
04 			#	heterozygous genotype messing things up?
05 			#if var in correctAF:
06 			#	a = correctAF[var]
07 			#	if a != None:
08 			#		a = AF_KEYS[quantize_AF(a)]
09 			#		if a not in alleleBal_vs_FN:
10 			#			alleleBal_vs_FN[a] = 0
11 			#		alleleBal_vs_FN[a] += 1
12 			#		if a < AF_THRESH:
13 			#			venn_data[i][2] = 1
14 
15 			#	no reason?
16 			if noReason:
17 				venn_data[i][2] += 1
18 
19 		for i in xrange(len(notFound)):
20 			if venn_data[i][0]: set1.append(i+varAdj)
21 			if venn_data[i][1]: set2.append(i+varAdj)
22 			if venn_data[i][2]: set3.append(i+varAdj)
23 		varAdj += len(notFound)
24 
25 		#
26 		#	if desired, write out vcf files.
27 		#
28 		notFound   = sorted(notFound)
29 		FPvariants = sorted(FPvariants)
30 		if VCF_OUT:
31 			for line in open(GOLDEN_VCF,'r'):
32 				if line[0] != '#':
33 					splt = line.split('\t')
34 					if splt[0] == refName:
35 						var  = (int(splt[1]),splt[3],splt[4])
36 						if var in notFound:
37 							vcfo2.write(line)
38 			for line in open(WORKFLOW_VCF,'r'):
39 				if line[0] != '#':
40 					splt = line.split('\t')
41 					if splt[0] == refName:
42 						var  = (int(splt[1]),splt[3],splt[4])
43 						if var in FPvariants:
44 							vcfo3.write(line)
45 
46 		print '{0:.3f} (sec)'.format(time.time()-tt)
47 
48 	#
49 	#	close vcf output
50 	#
51 	print ''
52 	if VCF_OUT:
53 		print OUT_PREFIX+'_FN.vcf'
54 		print OUT_PREFIX+'_FP.vcf'
55 		vcfo2.close()
56 		vcfo3.close()
57 
58 	#
59 	#	plot some FN stuff
60 	#
61 	if NO_PLOT == False:
62 		nDetected = len(set(set1+set2+set3))
63 		set1 = set(set1)
64 		set2 = set(set2)
65 		set3 = set(set3)
66 
67 		if len(set1): s1 = 'Unmappable'
68 		else: s1 = ''
69 		if len(set2): s2 = 'DP < '+str(DP_THRESH)
70 		else: s2 = ''
71 		#if len(set3): s3 = 'AF < '+str(AF_THRESH)
72 		if len(set3): s3 = 'Unknown'
73 		else: s3 = ''
74 
75 		mpl.figure(0)
76 		tstr1 = 'False Negative Variants (Missed Detections)'
77 		#tstr2 = str(nDetected)+' / '+str(znF)+' FN variants categorized'
78 		tstr2 = ''
79 		if MAPTRACK != None:
80 			v = venn3([set1, set2, set3], (s1, s2, s3))
81 		else:
82 			v = venn2([set2, set3], (s2, s3))
83 		mpl.figtext(0.5,0.95,tstr1,fontdict={'size':14,'weight':'bold'},horizontalalignment='center')
84 		mpl.figtext(0.5,0.03,tstr2,fontdict={'size':14,'weight':'bold'},horizontalalignment='center')
85 
86 		ouf = OUT_PREFIX+'_FNvenn.pdf'
87 		print ouf
88 		mpl.savefig(ouf)
89 
90 	#
91 	#	spit out results to console
92 	#
93 	print '\n**********************************\n'
94 	if BEDFILE != None:
95 		print 'ONLY CONSIDERING VARIANTS FOUND WITHIN TARGETED REGIONS\n\n'
96 	print 'Total Golden Variants:  ',ztV,'\t[',zgF,'filtered,',zgM,'merged,',zgR,'redundant ]'
97 	print 'Total Workflow Variants:',ztW,'\t[',zwF,'filtered,',zwM,'merged,',zwR,'redundant ]'
98 	print ''
99 	if ztV > 0 and ztW > 0:
00 		print 'Perfect Matches:',znP,'({0:.2f}%)'.format(100.*float(znP)/ztV)
01 		print 'FN variants:    ',znF,'({0:.2f}%)'.format(100.*float(znF)/ztV)
02 		print 'FP variants:    ',zfP#,'({0:.2f}%)'.format(100.*float(zfP)/ztW)
03 	if FAST == False:
04 		print '\nNumber of equivalent variants denoted differently between the two vcfs:',znE
05 	if BEDFILE != None:
06 		print '\nNumber of golden variants located in targeted regions that were too small to be sampled from:',zbM
07 	if FAST:
08 		print "\nWarning! Running with '--fast' means that identical variants denoted differently between the two vcfs will not be detected! The values above may be lower than the true accuracy."
09 	#if NO_PLOT:
10 	if True:
11 		print '\n#unmappable:  ',len(set1)
12 		print '#low_coverage:',len(set2)
13 		print '#unknown:     ',len(set3)
14 	print '\n**********************************\n'
15 
16 
17 
18 
19 
20 if __name__ == '__main__':
21 	main()