if thresholds are not specififed, calculate based on the distribution of normalized bit scores

def calc_thresholds(rbh, file_name, thresholds = [False, False, False, False], stdevs = 2):
    """
    if thresholds are not specififed, calculate based on the distribution of normalized bit scores
    """
    calc_threshold = thresholds[-1]
    norm_threshold = {}
    for pair in itertools.permutations([i for i in rbh], 2):
        if pair[0] not in norm_threshold:
            norm_threshold[pair[0]] = {}
        norm_threshold[pair[0]][pair[1]] = {}
    out = open(file_name, 'w')
    print('#### summary of rbh comparisons\n', file=out)
    comparisons = []
    for genome in rbh:
        for compare in rbh[genome]:
            pair = ''.join(sorted([genome, compare]))
            if pair in comparisons:
                continue
            comparisons.append(pair)
            scores = {'percent identity': [], 'e-value': [], 'bit score': [], 'normalized bit score': [], 'alignment length fraction': []}
            print('### blast between %s and %s\n' % (genome, compare), file=out)
            for id in rbh[genome][compare]:
                pident, length_fraction, e, bit, norm_bit = rbh[genome][compare][id][3:]
                scores['percent identity'].append(pident)
                scores['alignment length fraction'].append(length_fraction)
                scores['e-value'].append(e)
                scores['bit score'].append(bit)
                scores['normalized bit score'].append(norm_bit)
            if calc_threshold is True:
                norms = scores['normalized bit score']
                average = numpy.average(norms) 
                std = numpy.std(norms)
                normal_thresh = average - (std * stdevs)
                print('## average normalized bit score: %s' % average, file=out)
                print('## standard deviation of normalized bit scores: %s' % std, file=out)
                print('## normalized bit score threshold set to: %s\n' % (normal_thresh), file=out)
                norm_threshold[genome][compare], norm_threshold[compare][genome] = normal_thresh, normal_thresh
            for score in scores:
                print('## %s' % (score), file=out)
                if len(scores[score]) > 0:
                    print('## average: %s' % numpy.average(scores[score]), file=out)
#                    hist = histogram(scores[score], [])
#                    for line in hist:
#                        print >> out, line
                print('', file=out)
    out.close()
    if calc_threshold is True:
        return thresholds[0:-1] + [norm_threshold]
    else:
        return thresholds