357 lines
17 KiB
Python
Executable File
357 lines
17 KiB
Python
Executable File
#!/usr/bin/env python
|
|
|
|
import sys, itertools, FlatFileTable
|
|
from enum import Enum
|
|
|
|
def subset_list_by_indices(indices, list):
|
|
subset = []
|
|
for index in indices:
|
|
subset.append(list[index])
|
|
return subset
|
|
|
|
def chunk_generator(record_gen, key_fields):
|
|
"""Input:
|
|
record_gen: generator that produces dictionaries (records in database speak)
|
|
key_fields: keys in each dictionary used to determine chunk membership
|
|
Output:
|
|
locus_chunk: list of consecutive lines that have the same key_fields"""
|
|
|
|
locus_chunk = []
|
|
last_key = ""
|
|
first_record = True
|
|
for record in record_gen:
|
|
key = [record[f] for f in key_fields]
|
|
if key == last_key or first_record:
|
|
locus_chunk.append(record)
|
|
first_record = False
|
|
else:
|
|
if locus_chunk != []:
|
|
yield locus_chunk
|
|
locus_chunk = [record]
|
|
last_key = key
|
|
yield locus_chunk
|
|
|
|
|
|
class call_stats:
|
|
def __init__(self, call_type, coverage): #, acc_conf_calls, conf_call_rate, cum_corr_calls, cum_calls, coverage):
|
|
self.call_type = call_type
|
|
self.coverage = coverage
|
|
self.calls = 0
|
|
self.conf_ref_calls = 0
|
|
self.conf_het_calls = 0
|
|
self.conf_hom_calls = 0
|
|
self.conf_var_calls = 0
|
|
self.conf_genotype_calls = 0
|
|
self.conf_refvar_calls = 0
|
|
self.correct_genotype = 0
|
|
self.correct_refvar = 0
|
|
|
|
|
|
def add_stat(self, calls, conf_ref_calls, conf_het_calls, conf_hom_calls, conf_var_calls, conf_genotype_calls, conf_refvar_calls, correct_genotype, correct_refvar):
|
|
#print self, calls, conf_ref_calls, conf_het_calls, conf_hom_calls, conf_var_calls, conf_genotype_calls, conf_refvar_calls, correct_genotype, correct_refvar
|
|
self.calls += calls
|
|
self.conf_ref_calls += conf_ref_calls
|
|
self.conf_het_calls += conf_het_calls
|
|
self.conf_hom_calls += conf_hom_calls
|
|
self.conf_var_calls += conf_var_calls
|
|
self.conf_genotype_calls += conf_genotype_calls
|
|
self.conf_refvar_calls += conf_refvar_calls
|
|
self.correct_genotype += correct_genotype
|
|
self.correct_refvar += correct_refvar
|
|
|
|
@staticmethod
|
|
def calc_discovery_stats(chunk):
|
|
conf_calls = 0
|
|
correct_genotype = 0
|
|
calls = 0
|
|
for record in chunk:
|
|
calls += 1
|
|
if float(record["BtrLod"]) >= 5 or call_type.from_record_2_state(record["ReferenceBase"][0],record["BestGenotype"]) == call_type.call_types_2_state.Reference and float(record["BtnbLod"]) >= 5:
|
|
conf_calls += 1
|
|
if call_type.discovery_call_correct(record):
|
|
correct_genotype += 1
|
|
|
|
return correct_genotype, conf_calls, calls
|
|
|
|
@staticmethod
|
|
def calc_genotyping_stats(chunk):
|
|
conf_calls = 0
|
|
correct_genotype = 0
|
|
calls = 0
|
|
for record in chunk:
|
|
calls += 1
|
|
if float(record["BtnbLod"]) >= 5:
|
|
conf_calls += 1
|
|
if call_type.genotyping_call_correct(record):
|
|
correct_genotype += 1
|
|
|
|
return correct_genotype, conf_calls, calls
|
|
|
|
@staticmethod
|
|
def stats_header():
|
|
return "TrueGenotype,Coverage,AccuracyConfidentGenotypingCalls,ConfidentGenotypingCallRate,AccuracyConfidentDiscoveryCalls,ConfidentDiscoveryCallRate,Calls,ConfRefCalls,ConfHetCalls,ConfHomCalls,ConfGenotypeCalls,CorrectGenotypes,ConfVarCalls,ConfDiscoveryCalls,CorrectDiscovery"
|
|
|
|
def __str__(self):
|
|
return ",".join(map(str, (self.calls, self.conf_ref_calls, self.conf_het_calls, self.conf_hom_calls, self.conf_genotype_calls, self.correct_genotype, self.conf_var_calls, self.conf_refvar_calls, self.correct_refvar)))
|
|
|
|
def stats_str(self):
|
|
return "%s,%d,%.5f,%.5f,%.5f,%.5f,%s" % (self.call_type, self.coverage, float(self.correct_genotype) / max(self.conf_genotype_calls,1), float(self.conf_genotype_calls) / max(self.calls,1), float(self.correct_refvar) / max(self.conf_refvar_calls,1), float(self.conf_refvar_calls) / max(self.calls,1), self.__str__())
|
|
|
|
class call_type:
|
|
"""Class that returns an Enum with the type of call provided by a record"""
|
|
call_types_3_state = Enum("HomozygousSNP","HeterozygousSNP","HomozygousReference")
|
|
call_types_3_state_short = Enum("Hom","Het","Ref")
|
|
call_types_2_state = Enum("Variant","Reference")
|
|
|
|
@staticmethod
|
|
def from_record_3_state(ref, genotype): # record):
|
|
"""Given reference base as string, determine whether called genotype is homref, het, homvar"""
|
|
#ref = record["ReferenceBase"][0]
|
|
#genotype = record["HapmapChipGenotype"]
|
|
return call_type.call_types_3_state[genotype.count(ref)]
|
|
|
|
@staticmethod
|
|
def from_record_2_state(ref, genotype):
|
|
"""Given reference base as string, determine whether called genotype is ref or var"""
|
|
#ref = record["ReferenceBase"][0]
|
|
#genotype = record["HapmapChipGenotype"]
|
|
return call_type.call_types_2_state[0] if genotype.count(ref) < 2 else call_type.call_types_2_state[1]
|
|
|
|
@staticmethod
|
|
def genotyping_call_correct(record):
|
|
return record["HapmapChipGenotype"] == record["BestGenotype"]
|
|
|
|
|
|
@staticmethod
|
|
def discovery_call_correct(record):
|
|
return call_type.from_record_2_state(record["ReferenceBase"][0], record["HapmapChipGenotype"]) == call_type.from_record_2_state(record["ReferenceBase"][0], record["BestGenotype"])
|
|
|
|
def print_record_debug(rec):
|
|
print "TEST Genotyping:", call_type.genotyping_call_correct(rec)
|
|
print "TEST Discovery:", call_type.discovery_call_correct(rec)
|
|
print "DIFF:", call_type.genotyping_call_correct(rec) != call_type.discovery_call_correct(rec)
|
|
print "\n".join([" %20s => '%s'" % (k,v) for k,v in sorted(rec.items())])
|
|
print call_type.from_record_2_state(rec["ReferenceBase"][0],rec["HapmapChipGenotype"])
|
|
print call_type.from_record_2_state(rec["ReferenceBase"][0],rec["BestGenotype"])
|
|
print
|
|
|
|
def aggregate_stats(filename, max_loci, table_filename, debug):
|
|
aggregate = dict()
|
|
fout = open(table_filename,"w")
|
|
fout.write(call_stats.stats_header()+"\n")
|
|
|
|
locus_gen = chunk_generator(FlatFileTable.record_generator(filename, None), ("Sequence","Position"))
|
|
for index, locus_chunk in enumerate(locus_gen):
|
|
if index >= max_loci:
|
|
break
|
|
if (index % 1000) == 0:
|
|
sys.stderr.write( str(index)+" loci processed, at: "+locus_chunk[0]["Sequence"]+":"+locus_chunk[0]["Position"]+"\n")
|
|
|
|
covs = dict()
|
|
coverage_chunk_gen = chunk_generator(locus_chunk, ("DownsampledCoverage", "Sequence", "Position"))
|
|
for cov_chunk in coverage_chunk_gen:
|
|
first_record = cov_chunk[0]
|
|
#stat = call_stats.calc_stats(cov_chunk)
|
|
|
|
for record in cov_chunk:
|
|
hapmap_genotyping_call_type = call_type.from_record_3_state(record["ReferenceBase"][0],record["HapmapChipGenotype"])
|
|
key = hapmap_genotyping_call_type, int(record["DownsampledCoverage"])
|
|
#key = call_type.from_record_3_state(record)#, int(record["DownsampledCoverage"])
|
|
record["DownsampledCoverage"] = int(record["DownsampledCoverage"])
|
|
record["HapmapChipCallType"] = hapmap_genotyping_call_type
|
|
value = record
|
|
|
|
correct_genotype, conf_genotype_calls, genotype_calls = call_stats.calc_genotyping_stats([record])
|
|
correct_refvar, conf_refvar_calls, refvar_calls = call_stats.calc_discovery_stats([record])
|
|
assert(genotype_calls == refvar_calls)
|
|
|
|
conf_ref_calls = 0
|
|
conf_het_calls = 0
|
|
conf_hom_calls = 0
|
|
best_genotyping_call_type = call_type.from_record_3_state(record["ReferenceBase"][0],record["BestGenotype"])
|
|
if conf_genotype_calls:
|
|
if best_genotyping_call_type.index == 0: conf_hom_calls = 1
|
|
if best_genotyping_call_type.index == 1: conf_het_calls = 1
|
|
if best_genotyping_call_type.index == 2: conf_ref_calls = 1
|
|
|
|
conf_var_calls = 0
|
|
if conf_refvar_calls:
|
|
this_variant_call_type = call_type.from_record_2_state(record["ReferenceBase"][0],record["BestGenotype"])
|
|
conf_var_calls = 1 if this_variant_call_type.index == 0 else 0
|
|
|
|
aggregate.setdefault(key, call_stats(*key))
|
|
#print ",".join(map(str,(genotype_calls, conf_ref_calls, conf_het_calls, conf_hom_calls, conf_var_calls, conf_genotype_calls, conf_refvar_calls, correct_genotype, correct_refvar)))
|
|
aggregate[key].add_stat(genotype_calls, conf_ref_calls, conf_het_calls, conf_hom_calls, conf_var_calls, conf_genotype_calls, conf_refvar_calls, correct_genotype, correct_refvar)
|
|
|
|
if debug:# and conf_refvar_calls:
|
|
print "KEYS:",key
|
|
print_record_debug(record)
|
|
|
|
|
|
for key, records in sorted(aggregate.items()):
|
|
fout.write(records.stats_str()+"\n")
|
|
|
|
fout.close()
|
|
#return aggregate
|
|
|
|
class weighted_avg:
|
|
|
|
def __init__(self):
|
|
self.sum = 0.0
|
|
self.count = 0
|
|
|
|
def add(self, value, counts):
|
|
self.sum += value*counts
|
|
self.count += counts
|
|
|
|
def return_avg(self):
|
|
return float(self.sum) / max(self.count,1)
|
|
|
|
def stats_from_hist(options, depth_hist_filename, stats_filename, variant_eval_dir, depth_multiplier=1.0):
|
|
|
|
#hist_zero = {"CallType" : ,"Coverage","AccuracyConfidentCalls","ConfidentCallRate","CumCorrectCalls","CumCalls","CumgCorrectFraction"}
|
|
#prob_genotype = [1e-5, 1e-3, 1-1e-3]
|
|
#prob_genotype = [0.37, 0.62, .0]
|
|
#prob_genotype = [0.203, 0.304, .491]
|
|
#prob_genotype = [0.216, 0.302, .481]
|
|
prob_genotype = [0.205, 0.306, 0.491] # Based on CEU NA12878 actual hapmap chip calls
|
|
#prob_genotype = [0.213, 0.313, 0.474] # Based on YRB NA19240 actual hapmap chip calls
|
|
theta = 1.0/1850 # expected heterozygosity
|
|
prob_genotype = [0.5 * theta, 1.0 * theta, 1 - 1.5*theta] # Based on CEU NA12878 actual hapmap chip calls
|
|
|
|
|
|
hist = []
|
|
hist_gen = FlatFileTable.record_generator(depth_hist_filename, sep=" ", skip_n_lines=3)
|
|
for index, record in enumerate(hist_gen):
|
|
assert int(record["depth"]) == index
|
|
hist.append(int(record["count"]))
|
|
|
|
# If upsampling is not done in the CoverageEval GATK module, the number of observations of reads
|
|
# with high depth of coverage can be very low giving
|
|
|
|
stats_dict = dict()
|
|
stats_gen = FlatFileTable.record_generator(stats_filename, sep=",")
|
|
for record in stats_gen:
|
|
key1 = int(record["Coverage"])
|
|
key2 = record["TrueGenotype"]
|
|
#print key2
|
|
stats_dict.setdefault(key1, dict()) # create a nested dict if it doesn't exist
|
|
stats_dict[key1][key2] = record # create an entry for these keys
|
|
|
|
#highest_homref_calls = 0
|
|
#if record["TrueGenotype"] == "HomozygousReference":
|
|
# calls = int(record["Calls"])
|
|
# if calls > highest_homref_calls:
|
|
# highest_homref_calls = calls
|
|
|
|
acc = dict()
|
|
call_rate = dict()
|
|
conf_calls = dict()
|
|
|
|
start = 1
|
|
end = 1000
|
|
max_usable_depth = 40 # Depth of coverage beyond which stats are not sampled enough and we take the stat at this depth instead
|
|
for depth, depth_count in enumerate(hist[start:end+1],start): # For Cd = depth count
|
|
#print "DEPTH: "+str(depth)
|
|
try:
|
|
depth = max(int(float(depth*depth_multiplier)),1)
|
|
if depth > max_usable_depth: # Ensure that high entries with bad stats use a good stat from a depth that we now is well sampled
|
|
depth = max_usable_depth
|
|
depth_entries = stats_dict[depth]
|
|
except KeyError:
|
|
print "Stopped on depth",depth
|
|
break
|
|
if True:
|
|
for true_genotype, stat in depth_entries.items(): # For t (SNP type) = true_genotype
|
|
#print "TRUE_GENOTYPE: "+str(true_genotype)
|
|
for genotype in call_type.call_types_3_state:
|
|
conf_calls.setdefault(genotype, 0.0)
|
|
prob_conf_x_call = float(stat["Conf"+str(call_type.call_types_3_state_short[genotype.index])+"Calls"])/float(stat["Calls"])
|
|
conf_calls[genotype] += depth_count * prob_conf_x_call * prob_genotype[genotype.index]
|
|
#if genotype.index == 1:
|
|
# print "%.5f " % prob_conf_x_call, depth, depth_count, conf_calls[genotype], int(stat["Conf"+str(call_type.call_types_3_state_short[genotype.index])+"Calls"]), int(stat["Calls"])
|
|
|
|
acc.setdefault(true_genotype,weighted_avg())
|
|
call_rate.setdefault(true_genotype,weighted_avg())
|
|
acc[true_genotype].add(float(stat["AccuracyConfidentGenotypingCalls"]), depth_count)
|
|
call_rate[true_genotype].add(float(stat["ConfidentGenotypingCallRate"]), depth_count)
|
|
|
|
import numpy
|
|
for genotype in call_type.call_types_3_state:
|
|
print "%19s accuracy : %.3f" % (str(genotype), acc[str(genotype)].return_avg())
|
|
print "%19s call rate: %.3f" % (str(genotype), call_rate[str(genotype)].return_avg())
|
|
|
|
print "\nExpected performance given perfect accuracy and call rate:"
|
|
print "%19s %7s %7s %7s" % ("", "Actual", "Perfect", "Diff.")
|
|
total_hist_sites = numpy.sum(hist)
|
|
total_predicted = 0
|
|
for genotype in call_type.call_types_3_state:
|
|
predicted = conf_calls[genotype]
|
|
total_predicted += predicted
|
|
perfect = prob_genotype[genotype.index]*total_hist_sites
|
|
diff = perfect - predicted
|
|
print "%19s calls: %7.0f %7.0f %7.0f" % (genotype, predicted, perfect, diff)
|
|
#repl_string += "%s %.0f\n" % (genotype, predicted)
|
|
print " Total calls: %7d" % total_predicted
|
|
|
|
#stats_gen = FlatFileTable.record_generator(stats_filename, sep=",")
|
|
#for chunk in chunk_generator(stats_gen, key_fields=("True_Genotype")):
|
|
|
|
print "\nCoverage histogram mean: %.2f" % numpy.average(range(len(hist)), weights=hist)
|
|
|
|
# STEM AND LEAF PLOT
|
|
|
|
# If VariantEval directory given, compare with those results
|
|
if options.variant_eval_file != None:
|
|
vareval_file = open(options.variant_eval_file)
|
|
num_hets = None; num_homs = None
|
|
for line in vareval_file:
|
|
if "UNKNOWN_CALLED_VAR_HET_NO_SITES" in line: num_hets = line.rstrip().split()[2]
|
|
if "UNKNOWN_CALLED_VAR_HOM_NO_SITES" in line: num_homs = line.rstrip().split()[2]
|
|
|
|
if options.oneline_stats != None:
|
|
oneline_stats = open(options.oneline_stats, "w")
|
|
pred_hom = conf_calls[call_type.call_types_3_state[0]]
|
|
pred_het = conf_calls[call_type.call_types_3_state[1]]
|
|
print >>oneline_stats, "%s %.0f %s %.0f %s" % (options.oneline_stats, pred_hom, num_homs, pred_het, num_hets)
|
|
|
|
def usage(parser):
|
|
parser.print_usage()
|
|
sys.exit()
|
|
|
|
if __name__ == "__main__":
|
|
from optparse import OptionParser
|
|
|
|
parser = OptionParser()
|
|
parser.add_option("-c", "--genotype_call_file", help="GELI file to use in generating coverage stat table", dest="genotype_filename", )
|
|
parser.add_option("-s", "--stats_file", help="file to containing empirical genotyping stats", dest="stats_filename", default=None)
|
|
parser.add_option("-g", "--histogram_file", help="file containing counts of each depth of coverage", dest="hist_filename")
|
|
parser.add_option("-m", "--max_loci", help="maximum number of loci to parse (for debugging)", default=sys.maxint, dest="max_loci", type="int")
|
|
parser.add_option("-e", "--evaluate", help="evaluate genotypes; requires a stats file and a histogram file", default=False, dest="evaluate_genotypes", action="store_true")
|
|
parser.add_option("-d", "--debug", help="provide debugging output", default=False, dest="debug", action="store_true")
|
|
parser.add_option("-p", "--depth_multiplier", help="multiply all depths in histogram by this value; for \"downsampling\" depth", default=1.0, dest="depth_multiplier", type=float)
|
|
parser.add_option("-v", "--variant_eval_file", help="file with output of VariantEval genotype concordance to compare this prediction to", default=None, dest="variant_eval_file")
|
|
parser.add_option("-o", "--oneline_stats_file", help="output single, tabular line of stats to this file", default=None, dest="oneline_stats")
|
|
|
|
(options, args) = parser.parse_args()
|
|
|
|
if not options.evaluate_genotypes:
|
|
print "Creating performance tables from genotypes file"
|
|
#if options.stats_filename == None:
|
|
# sys.exit("Must provide -s stats fliname option")
|
|
if options.genotype_filename == None:
|
|
sys.exit("Must provide -c genotype call filename option")
|
|
stats_filename = options.stats_filename if options.stats_filename != None else options.genotype_filename+".stats"
|
|
aggregate_stats(options.genotype_filename, options.max_loci, options.stats_filename, options.debug)
|
|
else:
|
|
print "Evaluating genotypes"
|
|
print "Depth multiplier:",options.depth_multiplier
|
|
if options.hist_filename == None:
|
|
sys.exit("Must provide -g histogram filename option")
|
|
if options.stats_filename == None:
|
|
sys.exit("Must provide -s stats fliname option")
|
|
stats_from_hist(options, options.hist_filename, options.stats_filename, options.variant_eval_file, options.depth_multiplier)
|
|
|
|
|
|
|