make a rbh network for all pair-wise genome comparisons
    - filter network based on normalized bit score (default, automatic) or specified e-value / bit score
    - evaluate the scores for each genome pair compared
    - build second network that is filtered so that clusters have only one ORF per genome
    - evaluate the scores in this network
        - consider including fbhs under threshold if they don't violate cluster
    - compare 'connected-ness' for each genome

def network(thresholds, fastas, algorithm = 'usearch'):
    """
    make a rbh network for all pair-wise genome comparisons
    - filter network based on normalized bit score (default, automatic) or specified e-value / bit score
    - evaluate the scores for each genome pair compared
    - build second network that is filtered so that clusters have only one ORF per genome
    - evaluate the scores in this network
        - consider including fbhs under threshold if they don't violate cluster
    - compare 'connected-ness' for each genome
    """
    id2desc = get_descriptions(fastas)
            # get [fasta, description, length] for ORF id
    id2desc = self_compare(fastas, id2desc, algorithm)
            # get best possible bit score for each ORF 
            # (comparing with itself) [fasta, description, length, bestbit]
    hits = compare_genomes(fastas, id2desc, algorithm)
            # pair wise genome comparisons {genome: {id: [match_type = 'rbh' or 'fbh', scores]}}
    calc_thresholds(hits, file_name = 'fbh.scores.summary.txt')
    rbh_network(id2desc, hits, file_name = 'fbh.network.edges.txt')
    hits, rbh = find_rbh(hits, id2desc)
            # remove hits that are not reciprocal best blast hits
    thresholds = calc_thresholds(rbh, 'rbh.scores.summary.txt', thresholds)
            # print rbh score summary to rbh_score_summary.txt and
            # calculate normalized bit score cutoff for each pair of
            # genomes, if desired
    g = rbh_network(id2desc, rbh, file_name = 'rbh.network.edges.txt')
    filtered_g, filtered_rbh = rbh_network(id2desc, rbh, 'rbh.filtered.network.edges.txt', thresholds)
    calc_thresholds(filtered_rbh, file_name = 'rbh.filtered.scores.summary.txt')
    print_summary(filtered_g, fastas, id2desc, file_name = 'rbh.filtered.network.nodes.txt')
    print_network_matrix(filtered_g, fastas, id2desc, file_name = 'rbh.filtered.network.matrix.txt')
    print_genome_matrix(filtered_rbh, fastas, id2desc, file_name = 'rbh.filtered.network.genome_matrix.txt')
    split_g = split_network(filtered_g, id2desc, file_name = 'rbh.filtered.split.network.edges.txt')
    print_summary(split_g, fastas, id2desc, file_name = 'rbh.filtered.split.network.nodes.txt')
    print_network_matrix(split_g, fastas, id2desc, file_name = 'rbh.filtered.split.network.matrix.txt')
    return split_g