<filename>diana/toolbox/wrappers.py
#######################################################################
# Recipies / wrapper functions using toolbox methods for disease,
# drug and network analysis
# e.g. 10/2015
#######################################################################
import diana.toolbox.network_utilities as network_utilities # For DIANA
import diana.toolbox.stat_utilities as stat_utilities # For DIANA
import diana.toolbox.dict_utilities as dict_utilities # For DIANA
import diana.toolbox.text_utilities as text_utilities # For DIANA
import diana.toolbox.TsvReader # For DIANA
import diana.toolbox.functional_enrichment # For DIANA
import diana.toolbox.parse_umls # For DIANA
import diana.toolbox.parse_msigdb # For DIANA
import diana.toolbox.parse_uniprot # For DIANA
import diana.toolbox.parse_ncbi # For DIANA
import diana.toolbox.parse_do # For DIANA
import diana.toolbox.parse_medic # For DIANA
import diana.toolbox.parse_disgenet # For DIANA
import diana.toolbox.parse_drugbank # For DIANA
import diana.toolbox.parse_medi # For DIANA
import diana.toolbox.parse_hetionet # For DIANA
import csv, numpy, os, cPickle
import random
try:
    from toolbox.external.diamond import diamond
except:
    pass # For DIANA
    #print("DIAMOnD not found and thus will not be available!")

##### Id mapping related #####

def get_mapping(file_name, from_column = None, to_column = None, delim=None, one_to_one=True):
    """
    Assumes header, maps first column to the second
    """
    key_to_value = TsvReader.get_from_to_mapping(file_name, from_column = from_column, to_column = to_column, delim=delim, inner_delim = None, filter_column = None, exclude_value = None, include_value = None, one_to_one=one_to_one) 
    return key_to_value

def convert_to_geneid(file_name, id_type, id_mapping_file):
    """
    Expects a file where each line is a gene name / uniprot id
    id_type: symbol | uniprot
    """
    genes = [ line.strip("\n") for line in open(file_name) ]
    if id_type == "symbol":
        geneid_to_name, name_to_geneid = get_geneid_symbol_mapping(id_mapping_file)
    elif id_type == "uniprot":
        name_to_geneid = get_uniprot_to_id(id_mapping_file, uniprot_ids=genes, only_min=True, key_function=int)
    else:
        raise ValueError("Uknown id type: %s" % id_type)
    geneids = set([ name_to_geneid[gene] for gene in genes if gene in name_to_geneid ])
    genes_non = set([ gene for gene in genes if gene not in name_to_geneid ])
    print("Not found genes: {}".format(genes_non))
    return geneids


def get_uniprot_to_id(uniprot_file, uniprot_ids=None, only_min=True, key_function=int):
    """
    uniprot_file = %(data_dir)s/uniprot/idmapping.tab or idmapping.tab.symbol or idmapping.tab.mouse
    Can be used to convert to geneids (key_function=int) as well as symbols (key_function=len) depending on the input file
    """
    uniprot_to_gene = parse_uniprot.get_uniprot_to_geneid(uniprot_file, uniprot_ids, only_min, key_function)
    return uniprot_to_gene


def get_geneid_symbol_mapping(mapping_file):
    """
    id_mapping_file = %(data_dir)s/ncbi/geneid_to_symbol.txt
    """
    geneid_to_name, name_to_geneid = parse_ncbi.get_geneid_symbol_mapping(mapping_file)
    return geneid_to_name, name_to_geneid


def get_mesh_id_mapping(desc_file, rel_file, dump_file = None):
    """ 
    Get all concept id - mesh id mapping (also gets entry names in addition to main header)
    """
    #dump_file = CONFIG.get("umls_dir") + "/mapping.pcl"
    mesh_id_to_name, concept_id_to_mesh_id, mesh_id_to_name_with_synonyms = parse_umls.get_mesh_id_mapping(desc_file, rel_file, dump_file = dump_file)
    return mesh_id_to_name, concept_id_to_mesh_id, mesh_id_to_name_with_synonyms


def get_mesh_disease_ontology(desc_file, rel_file, dump_file = None):
    #dump_file = CONFIG.get("umls_dir") + "/ontology.pcl"
    g = parse_umls.get_mesh_disease_ontology(desc_file, rel_file, dump_file = dump_file)
    return g


def get_medic_mesh_id_mapping(medic_file):
    name_to_id, id_to_mesh_ids = get_medic_mesh_id_mapping(medic_file)
    return name_to_id, id_to_mesh_ids 


def get_mesh_disease_category_mapping(desc_file, rel_file, dump_file = None):
    #dump_file = CONFIG.get("umls_dir") + "/ontology.pcl"
    mesh_id_to_top_ids = parse_umls.get_mesh_id_to_disease_category(desc_file, rel_file, dump_file)
    mesh_id_to_name, concept_id_to_mesh_id, mesh_id_to_name_with_synonyms = get_mesh_id_mapping(desc_file, rel_file, dump_file)
    mesh_name_to_parents = {}
    for child, parents in mesh_id_to_top_ids.iteritems():
        name_child = mesh_id_to_name[concept_id_to_mesh_id[child]]
        values = []
        for parent in parents:
            name_parent = mesh_id_to_name[concept_id_to_mesh_id[parent]]
            values.append(name_parent.lower())
        values.sort()
        mesh_name_to_parents[name_child.lower()] = values
    #print mesh_name_to_parents["asthma"], mesh_name_to_parents["psoriasis"]
    return mesh_name_to_parents


def get_icd_to_mesh_ids(disease_ontology_file, id_type="ICD9CM"):
    icd_to_mesh_ids = parse_do.get_icd_to_mesh_ids(disease_ontology_file, id_type)
    return icd_to_mesh_ids


def get_do_to_mesh_ids(disease_ontology_file):
    name_to_do_id, do_to_mesh_ids, mesh_id_to_type_to_ids = parse_do.get_do_mesh_id_mapping(disease_ontology_file)
    return do_to_mesh_ids


def get_homology_mapping(homologene_file, tax_id="10090", from_tax_id="9606", symbol_type="geneid"):
    """
    symbol_type: geneid | symbol
    """
    geneid_to_geneid, group_to_taxid_to_geneid = parse_ncbi.get_homology_mapping(homologene_file, tax_id, from_tax_id=from_tax_id, symbol_type=symbol_type)
    return geneid_to_geneid 


##### Network related #####

def get_network(network_file, only_lcc):
    network = network_utilities.create_network_from_sif_file(network_file, use_edge_data = False, delim = None, include_unconnected=True)
    #print len(network.nodes()), len(network.edges())
    if only_lcc and not network_file.endswith(".lcc"):
        print("Shrinking network to its LCC: {} nodes, {} edges".format(len(network.nodes()), len(network.edges())))
        components = network_utilities.get_connected_components(network, False)
        network = network_utilities.get_subgraph(network, components[0])
        print("Final shape: {} nodes, {} edges".format(len(network.nodes()), len(network.edges())))
        #print len(network.nodes()), len(network.edges())
        network_lcc_file = network_file + ".lcc"
        if not os.path.exists(network_lcc_file ):
            f = open(network_lcc_file, 'w')
            for u,v in network.edges():
                f.write("%s 1 %s\n" % (u, v))
            f.close()
    return network


def create_functional_network(links_file, mapping_file, cutoff = 900): 
    #string_dir = CONFIG.get("string_dir") + "/"
    #links_file = string_dir + CONFIG.get("string_links_file")
    #mapping_file = string_dir + CONFIG.get("string_mapping_file")
    output_file = CONFIG.get("network_file")
    parse_string.get_interactions(links_file, mapping_file, output_file, cutoff) #, include_score=True) 
    #network = get_network()
    #print len(network.nodes()), len(network.edges())
    return


def calculate_lcc_significance(network, nodes, nodes_random=None, bins=None, n_random=1000, min_bin_size=100, seed=452456):
    # Degree matching problematic for small bin sizes
    #if bins is None and nodes_random is None:
    #        bins = network_utilities.get_degree_binning(network, min_bin_size) 
    random.seed(seed)
    if nodes_random is None:
        network_nodes = list(network.nodes())
        #nodes_random = get_random_nodes(nodes, network, bins = bins, n_random = n_random, min_bin_size = min_bin_size, seed = seed)
        nodes_random = []
        for i in xrange(n_random):
            random.shuffle(network_nodes)
            nodes_random.append(network_nodes[:len(nodes)])
    network_sub = network.subgraph(nodes)
    component_nodes = network_utilities.get_connected_components(network_sub, False)
    #print component_nodes 
    d = len(component_nodes[0])
    values = numpy.empty(len(nodes_random)) 
    for i, nodes in enumerate(nodes_random):
        network_sub = network.subgraph(nodes)
        component_nodes = network_utilities.get_connected_components(network_sub, False)[0]
        values[i] = len(component_nodes)
    m, s = numpy.mean(values), numpy.std(values)
    if s == 0:
        z = 0.0
    else:
        z = (d - m) / s
    return d, z, (m, s), values


##### Gene expression related #####

def get_expression_info(gexp_file, process=None, delim=',', quote='"', R_header=False, dump_file=None):
    """
    To get gene expression info
    process: a set(["log2", "z", "abs"]) or None
    """
    if dump_file is not None and os.path.exists(dump_file):
        gexp, gene_to_idx, cell_line_to_idx = cPickle.load(open(dump_file))
        return gexp, gene_to_idx, cell_line_to_idx
    #gene_to_values = {}
    f = open(gexp_file)
    reader = csv.reader(f, delimiter=delim, quotechar=quote)
    header = reader.next()
    #print len(header), header
    if R_header == False:
        header = header[1:]
    cell_line_to_idx = dict([ (cell_line, i) for i, cell_line in enumerate(header) ])
    gene_to_idx = {}
    values_arr = []
    for i, row in enumerate(reader):
        gene = row[0]
        values = map(float, row[1:])
        #gene_to_values[gene] = values
        gene_to_idx[gene] = i
        values_arr.append(values)
    f.close()
    gexp = numpy.array(values_arr)
    if process is not None:
        if "log2" in process:
            gexp = numpy.log2(gexp)
        if "z" in process:
            gexp = (gexp - gexp.mean(axis=1)[:, numpy.newaxis]) / gexp.std(axis=1, ddof=1)[:, numpy.newaxis]
        if "abs" in process:
            gexp = numpy.abs(gexp)
        #if "na.rm" in process:
        #    idx = numpy.where(numpy.isnan(a)) # need to remove rows with NAs
        #print gexp.shape, gexp_norm.shape
        #print gexp[0,0], gexp_norm[0,0]
        #return gene_to_values, cell_line_to_idx
    if dump_file is not None:
        values = gexp, gene_to_idx, cell_line_to_idx
        cPickle.dump(values, open(dump_file, 'w')) 
    return gexp, gene_to_idx, cell_line_to_idx


def get_de_genes(file_name, cutoff_adj = 0.05, cutoff_logfc=0.585, n_top=None, id_type = "GeneID"):
    """
    For parsing DE file generated using R PEPPER package
    """
    fields_to_include = [id_type, "P.Value", "logFC", "adj.P.Val"]
    parser = TsvReader.TsvReader(file_name, delim="\t", inner_delim=None)
    header_to_idx, id_to_values = parser.read(fields_to_include, keys_to_include=None, merge_inner_values=False)
    if "" in id_to_values:
        del id_to_values[""]
    #print len(id_to_values)
    #gene = "10458"
    #if gene in id_to_values:
    #    print id_to_values[gene] 
    genes = set()
    genes_all = set()
    genes_up, genes_down = set(), set()
    values_gene = []
    for gene, values in id_to_values.iteritems():
        include = False
        positive = False
        for val in values:
            pval = val[header_to_idx["adj.p.val"]] # "p.value"]]
            if pval == "NA":
                continue
            fc = float(val[header_to_idx["logfc"]])
            if float(pval) <= cutoff_adj:
                if abs(fc) >= cutoff_logfc: 
                    include = True
                if fc >= 0: 
                    positive = True
                if n_top is not None:
                    values_gene.append((abs(fc), gene))
        for word in gene.split("///"):
            word = word.strip()
            if word == "---":
                continue
            if include:
                genes.add(word)
                if positive:
                    genes_up.add(word)
                else:
                    genes_down.add(word)
            else:
                genes_all.add(word)
    if n_top is not None:
        values_gene.sort()
        genes = set([ word.strip() for fc, gene in values_gene[-n_top:] for word in gene.split("///") ])
    return genes, genes_all, genes_up, genes_down


def get_z_genes(file_name, cutoff_z = 2):
    """
    For parsing DE-Z file generated using R PEPPER package
    """
    gexp, gene_to_idx, cell_line_to_idx = get_expression_info(file_name, process=None, delim='\t', R_header=True) #, quote='"', dump_file=None)
    genes = gene_to_idx.items()
    genes.sort(key=lambda x: x[1])
    genes = numpy.array(zip(*genes)[0])
    sample_to_genes = {}
    for cell_line, idx in cell_line_to_idx.iteritems():
        indices = numpy.abs(gexp[:,idx]) > cutoff_z 
        sample_to_genes[cell_line] = genes[indices]
        #if cell_line in ["GSM734834", "GSM734833"]: 
        #    print cell_line, len(genes[indices]), genes[indices] 
    return sample_to_genes


def get_sample_mapping(file_name, labels_case, labels_control=None):
    f = open(file_name)
    labels_case = set(labels_case)
    if labels_control is not None:
        labels_control = set(labels_control)
    samples_case = []
    samples_control = []
    for line in f:
        sample, label = line.strip("\n").split("\t")
        label = label.strip()
        if label in labels_case:
            samples_case.append(sample)
        else:
            if labels_control is None or label in labels_control:
                samples_control.append(sample)
    return samples_case, samples_control


##### Disease, pathway, comorbidity, symptom info related #####

def get_pathway_info(pathway_file, prefix=None, nodes=None, max_pathway_size=None, inner_delim=None):
    """
    Assumes a tab separated file containing pathway name, link, geneids
    nodes to filter geneids that are not in the network
    prefix: kegg | reactome | biocarta
    inner_delim: None for tab separated geneids, " " for space separated geneids
    """
    pathway_to_geneids, geneid_to_pathways = parse_msigdb.get_msigdb_info(pathway_file, prefix, inner_delim=inner_delim)
    if nodes is not None or max_pathway_size is not None:
        pathway_to_geneids_mod = {}
        for pathway, geneids in pathway_to_geneids.iteritems():
            if max_pathway_size is not None:
                if len(geneids) > max_pathway_size:
                    continue
            if nodes is not None:
                geneids &= nodes
                if len(geneids) == 0:
                    continue
            pathway_to_geneids_mod[pathway] = geneids
        pathway_to_geneids = pathway_to_geneids_mod
    return pathway_to_geneids


def get_diseasome_genes(diseasome_file, nodes=None, network=None):
    """
    If nodes is not None, keep only nodes in the network
    If network is not None, keep only LCC
    """
    disease_to_genes = {}
    disease_to_category = {}
    for line in open(diseasome_file):
        words = line.strip("\n").split("\t")
        disease = words[1].strip('"')
        category = words[0]
        genes = set(words[2:])
        if nodes is not None:
            genes &= nodes
            if len(genes) == 0:
                continue
        if network is not None:
            network_sub = network.subgraph(genes)
            genes = network_utilities.get_connected_components(network_sub, False)[0]
        disease_to_genes[disease] = genes
        disease_to_category[disease] = category
    return disease_to_genes, disease_to_category


def get_disgenet_genes(file_name):
    disease_to_genes, disease_to_sources, cui_to_disease = parse_disgenet.get_disgenet_genes(file_name)
    return disease_to_genes, disease_to_sources, cui_to_disease 


def get_comorbidity_info(comorbidity_file, disease_ontology_file, mesh_dump, correlation_type="RR", only_significant=False):
    """
    Parse HuDiNe data from AllNet3 and map ICD9 to MeSH using DO
    correlation_type: phi (pearson correlation) | RR (favors rare disease pairs)
    """
    icd_to_mesh_ids = get_icd_to_mesh_ids(disease_ontology_file, id_type="ICD9CM")
    #print len(icd_to_mesh_ids), icd_to_mesh_ids.items()[:5]
    #print [ icd_to_mesh_ids[val] for val in ("289", "502", "578", "579", "542", "543")]
    mesh_id_to_name, concept_id_to_mesh_id, mesh_id_to_name_with_synonyms = get_mesh_id_mapping(None, None, dump_file = mesh_dump)
    #print mesh_id_to_name.items()[:5]
    f = open(comorbidity_file)
    header_to_idx = dict((word, i) for i, word in enumerate(f.readline().strip().split("\t")))
    disease_to_disease_comorbidity = {}
    for line in f:
        words = line.strip().split("\t")
        icd1, icd2 = words[:2]
        if only_significant and significance == "0":
            continue
        if icd1 not in icd_to_mesh_ids or icd2 not in icd_to_mesh_ids:
            #print "Not in DO mapping:", icd1, icd2
            continue
        val = float(words[header_to_idx[correlation_type]]) # idx:5
        for mesh1 in icd_to_mesh_ids[icd1]:
            if mesh1 not in mesh_id_to_name:
                #print "Not in name mapping:", mesh1
                continue
            disease1 = mesh_id_to_name[mesh1].lower()
            for mesh2 in icd_to_mesh_ids[icd2]:
                if mesh2 not in mesh_id_to_name:
                    #print "Not in name mapping:", mesh2
                    continue
                disease2 = mesh_id_to_name[mesh2].lower()
                disease1_mod, disease2_mod = sorted((disease1, disease2))
                d = disease_to_disease_comorbidity.setdefault(disease1, {})
                if disease2 in d:
                    if d[disease2] > val: # skip if the existing comorbidity value is higher
                        continue
                d[disease2] = val 
                d = disease_to_disease_comorbidity.setdefault(disease2, {})
                d[disease1] = val
                #print icd1, mesh1, disease1, icd2, mesh2, disease2, val 
    #print len(disease_to_disease_comorbidity), disease_to_disease_comorbidity.values()[0].items()[:5]
    return disease_to_disease_comorbidity 
    # Parse HuDiNe data from potentially buggy comorbidity_new.tsv
    #comorbidity_file = CONFIG.get("comorbidity_file")
    f = open(comorbidity_file)
    header_to_idx = dict((word, i) for i, word in enumerate(f.readline().strip().split("\t")))
    disease_to_disease_comorbidity = {}
    for line in f:
        words = line.strip().split("\t")
        disease1, disease2 = words[:2]
        significance = words[header_to_idx["sign_"+correlation_type]]
        if only_significant and significance == "0":
            continue
        val = float(words[header_to_idx[correlation_type]])
        disease_to_disease_comorbidity.setdefault(disease1, {})[disease2] = (val, significance)
        disease_to_disease_comorbidity.setdefault(disease2, {})[disease1] = (val, significance)
    f.close()
    return disease_to_disease_comorbidity


def get_symptom_info(symptom_file, tfidf_cutoff=None):
    """
    Parse Zhou et al supplementary s4. A cutoff of 3.5 is likely to filter spurious associations.
    """
    disease_to_symptoms = {}
    symptom_to_diseases = {}
    disease_to_symptom_to_score = {}
    #symptom_file = CONFIG.get("symptom_file")
    f = open(symptom_file)
    f.readline()
    for line in f:
        words = line.strip("\n").split("\t")
        symptom, disease, n, score = words
        symptom = symptom.lower()
        disease = disease.lower()
        if tfidf_cutoff is not None and not float(score) >= tfidf_cutoff:
            continue
        disease_to_symptoms.setdefault(disease, set()).add(symptom)
        symptom_to_diseases.setdefault(symptom, set()).add(disease)
        d = disease_to_symptom_to_score.setdefault(disease, {})
        d[symptom] = float(score)
    return disease_to_symptoms, symptom_to_diseases, disease_to_symptom_to_score


##### Drug related info #####

def get_drugbank(drugbank_file):
    dump_file = drugbank_file + ".pcl"
    if os.path.exists(dump_file):
        parser = cPickle.load(open(dump_file))
    else:
        parser = parse_drugbank.DrugBankXMLParser(drugbank_file)
        parser.parse()
        cPickle.dump(parser, open(dump_file, 'w'))
    return parser


def get_medi_indications(medi_file, drugbank_file, mesh_dump, disease_ontology_file, only_hps=True):
    dump_file = medi_file + ".pcl"
    if os.path.exists(dump_file):
        drug_to_diseases = cPickle.load(open(dump_file))
        return drug_to_diseases 
    parser = get_drugbank(drugbank_file)
    name_to_drug, synonym_to_drug = parser.get_synonyms(selected_drugs=None, only_synonyms=False)
    name_to_icd_and_confidences = parse_medi.get_medi_mapping(medi_file)
    mesh_id_to_name, concept_id_to_mesh_id, mesh_id_to_name_with_synonyms = get_mesh_id_mapping(None, None, dump_file = mesh_dump)
    icd_to_mesh_ids = get_icd_to_mesh_ids(disease_ontology_file, id_type="ICD9CM")
    drug_to_indications = parse_medi.get_drug_disease_mapping(name_to_icd_and_confidences, name_to_drug, synonym_to_drug, icd_to_mesh_ids, mesh_id_to_name, dump_file = None) 
    drug_to_diseases = {}
    for drug, values in drug_to_indications.iteritems():
        for phenotype, dui, val in values:
            if only_hps and val <= 0.5:
                continue
            drug_to_diseases.setdefault(drug, set()).add(phenotype)
    # Drug name to disease name textual mapping
    #name_to_indication_and_confidences = parse_medi.get_medi_mapping(medi_file, textual_indication=True)
    #drug_to_diseases = {}
    #for name, values in name_to_indication_and_confidences.iteritems():
    #        for indication, confidence in values:
    #            if only_hps and confidence <= 0.5:
    #                continue
    #            drug_to_diseases.setdefault(name, set()).add(indication.lower())
    # Get disease to name mapping
    #phenotype_to_mesh_id = dict((name, mesh_id) for mesh_id, name in mesh_id_to_name.iteritems())
    #disease_to_drugs = parse_medi.get_disease_specific_drugs(drug_to_diseases, phenotype_to_mesh_id)
    cPickle.dump(drug_to_diseases, open(dump_file, 'w'))
    return drug_to_diseases


def get_hetionet_indications(hetionet_file, mesh_dump, disease_ontology_file):
    dump_file = hetionet_file + ".pcl"
    if os.path.exists(dump_file):
        drug_to_diseases = cPickle.load(open(dump_file))
        return drug_to_diseases 
    drug_to_do_ids = parse_hetionet.get_hetionet_mapping(hetionet_file, metaedge="CtD")
    do_to_mesh_ids = get_do_to_mesh_ids(disease_ontology_file)
    mesh_id_to_name, concept_id_to_mesh_id, mesh_id_to_name_with_synonyms = get_mesh_id_mapping(None, None, dump_file = mesh_dump)
    drug_to_indications = parse_hetionet.get_drug_disease_mapping(drug_to_do_ids, do_to_mesh_ids, mesh_id_to_name, dump_file = None)
    drug_to_diseases = {}
    for drug, values in drug_to_indications.iteritems():
        for phenotype, dui, val in values:
            drug_to_diseases.setdefault(drug, set()).add(phenotype)
    cPickle.dump(drug_to_diseases, open(dump_file, 'w'))
    return drug_to_diseases 


##### Statistics related #####

def overlap_significance(geneids1, geneids2, nodes, method="hyper"):
    """
    method: hyper(geometric) | fishers (two-sided version of hypergeometric) | jaccard | jaccard_max | overlap
    """
    n1, n2 = len(geneids1), len(geneids2)
    n = len(geneids1 & geneids2)
    N = len(nodes)
    if method == "hyper":
        val = stat_utilities.hypergeometric_test_numeric(n, n1, N, n2)
    elif method == "fishers":
        oddsratio, val = stat_utilities.fisher_exact(n, n1 - n, n2 -n, N - n1 - n2 + n, alternative="two-sided")
    elif method == "jaccard":
        val = stat_utilities.jaccard(geneids1, geneids2)
    elif method == "jaccard_max":
        val = stat_utilities.jaccard_max(geneids1, geneids2)
    elif method == "overlap":
        val = n
    else:
        raise ValueError("Uknown method: %s" % method)
    return n, n1, n2, val


##### Proximity related #####

def calculate_proximity(network, nodes_from, nodes_to, nodes_from_random=None, nodes_to_random=None, bins=None, n_random=1000, min_bin_size=100, seed=452456, lengths=None):
    """
    Calculate proximity from nodes_from to nodes_to
    If degree binning or random nodes are not given, they are generated
    lengths: precalculated shortest path length dictionary
    """
    #distance = "closest"
    #lengths = network_utilities.get_shortest_path_lengths(network, "../data/toy.sif.pcl")
    #d = network_utilities.get_separation(network, lengths, nodes_from, nodes_to, distance, parameters = {})
    nodes_network = set(network.nodes())
    nodes_from = set(nodes_from) & nodes_network 
    nodes_to = set(nodes_to) & nodes_network
    if len(nodes_from) == 0 or len(nodes_to) == 0:
        return None # At least one of the node group not in network
    d = calculate_closest_distance(network, nodes_from, nodes_to, lengths)
    if bins is None and (nodes_from_random is None or nodes_to_random is None):
        bins = network_utilities.get_degree_binning(network, min_bin_size, lengths) # if lengths is given, it will only use those nodes
    if nodes_from_random is None:
        nodes_from_random = get_random_nodes(nodes_from, network, bins = bins, n_random = n_random, min_bin_size = min_bin_size, seed = seed)
    if nodes_to_random is None:
        nodes_to_random = get_random_nodes(nodes_to, network, bins = bins, n_random = n_random, min_bin_size = min_bin_size, seed = seed)
    random_values_list = zip(nodes_from_random, nodes_to_random)
    values = numpy.empty(len(nodes_from_random)) #n_random
    for i, values_random in enumerate(random_values_list):
        nodes_from, nodes_to = values_random
        #values[i] = network_utilities.get_separation(network, lengths, nodes_from, nodes_to, distance, parameters = {})
        values[i] = calculate_closest_distance(network, nodes_from, nodes_to, lengths)
    #pval = float(sum(values <= d)) / len(values) # needs high number of n_random
    m, s = numpy.mean(values), numpy.std(values)
    if s == 0:
        z = 0.0
    else:
        z = (d - m) / s
    return d, z, (m, s) #(z, pval)


def calculate_proximity_multiple(network, from_file=None, to_file=None, n_random=1000, min_bin_size=100, seed=452456, lengths=None, out_file="output.txt"):
    """
    Run proximity on each entries of from and to files in a pairwise manner
    output is saved in out_file (e.g., output.txt)
    """
    nodes = set(network.nodes())
    drug_to_targets, drug_to_category = get_diseasome_genes(from_file, nodes = nodes)
    #drug_to_targets = dict((drug, nodes & targets) for drug, targets in drug_to_targets.iteritems())
    disease_to_genes, disease_to_category = get_diseasome_genes(to_file, nodes = nodes)
    # Calculate proximity values
    print(len(drug_to_targets), len(disease_to_genes))
    # Get degree binning 
    bins = network_utilities.get_degree_binning(network, min_bin_size)
    f = open(out_file, 'w')
    f.write("source\ttarget\tn.source\tn.target\td\tz\n")
    for drug, nodes_from in drug_to_targets.iteritems():
        values = []
        for disease, nodes_to in disease_to_genes.iteritems():
            print(drug, disease)
            d, z, (m, s) = calculate_proximity(network, nodes_from, nodes_to, nodes_from_random=None, nodes_to_random=None, bins=bins, n_random=n_random, min_bin_size=min_bin_size, seed=seed, lengths=lengths)
            values.append((drug, disease, z, len(nodes_from), len(nodes_to), d, m, s))
            #f.write("%s\t%s\t%f\t%f\t%f\t%f\n" % (drug, disease, z, d, m, s))
        values.sort(key=lambda x: x[2])
        for drug, disease, z, k, l, d, m, s in values:
            #f.write("%s\t%s\t%f\t%d\t%d\t%f\t%f\t%f\n" % (drug, disease, z, k, l, d, m, s))
            f.write("%s\t%s\t%d\t%d\t%f\t%f\n" % (drug, disease, k, l, d, z))
    f.close()
    return 


def calculate_closest_distance(network, nodes_from, nodes_to, lengths=None):
    values_outer = []
    if lengths is None:
        for node_from in nodes_from:
            values = []
            for node_to in nodes_to:
                val = network_utilities.get_shortest_path_length_between(network, node_from, node_to)
                values.append(val)
            d = min(values)
            #print d,
            values_outer.append(d)
    else:
        for node_from in nodes_from:
            values = []
            vals = lengths[node_from]
            for node_to in nodes_to:
                val = vals[node_to]
                values.append(val)
            d = min(values)
            values_outer.append(d)
    d = numpy.mean(values_outer)
    #print d
    return d


def get_random_nodes(nodes, network, bins=None, n_random=1000, min_bin_size=100, degree_aware=True, seed=None):
    if bins is None:
        # Get degree bins of the network
        bins = network_utilities.get_degree_binning(network, min_bin_size) 
    nodes_random = network_utilities.pick_random_nodes_matching_selected(network, bins, nodes, n_random, degree_aware, seed=seed) 
    return nodes_random


### Separation related

def calculate_separation_proximity(network, nodes_from, nodes_to, nodes_from_random=None, nodes_to_random=None, bins=None, n_random=1000, min_bin_size=100, seed=452456, lengths=None):
    """
    Calculate proximity from nodes_from to nodes_to
    If degree binning or random nodes are not given, they are generated
    lengths: precalculated shortest path length dictionary
    """
    nodes_network = set(network.nodes())
    if len(set(nodes_from) & nodes_network) == 0 or len(set(nodes_to) & nodes_network) == 0:
        return None # At least one of the node group not in network
    d = get_separation(network, nodes_from, nodes_to, lengths)
    if bins is None and (nodes_from_random is None or nodes_to_random is None):
        bins = network_utilities.get_degree_binning(network, min_bin_size, lengths) # if lengths is given, it will only use those nodes
    if nodes_from_random is None:
        nodes_from_random = get_random_nodes(nodes_from, network, bins = bins, n_random = n_random, min_bin_size = min_bin_size, seed = seed)
    if nodes_to_random is None:
        nodes_to_random = get_random_nodes(nodes_to, network, bins = bins, n_random = n_random, min_bin_size = min_bin_size, seed = seed)
    random_values_list = zip(nodes_from_random, nodes_to_random)
    values = numpy.empty(len(nodes_from_random)) #n_random
    for i, values_random in enumerate(random_values_list):
        nodes_from, nodes_to = values_random
        values[i] = get_separation(network, nodes_from, nodes_to, lengths)
    m, s = numpy.mean(values), numpy.std(values)
    if s == 0:
        z = 0.0
    else:
        z = (d - m) / s
    return d, z, (m, s) #(z, pval)


def get_separation(network, nodes_from, nodes_to, lengths=None):
    dAA = numpy.mean(get_separation_within_set(network, nodes_from, lengths))
    dBB = numpy.mean(get_separation_within_set(network, nodes_to, lengths))
    dAB = numpy.mean(get_separation_between_sets(network, nodes_from, nodes_to, lengths))
    d = dAB - (dAA + dBB) / 2.0
    return d


def get_separation_between_sets(network, nodes_from, nodes_to, lengths=None):
    """
    Calculate dAB in separation metric proposed by Menche et al. 2015
    """
    values = []
    target_to_values = {}
    source_to_values = {}
    for source_id in nodes_from:
        for target_id in nodes_to:
            if lengths is not None:
                d = lengths[source_id][target_id] 
            else:
                d = network_utilities.get_shortest_path_length_between(network, source_id, target_id)
            source_to_values.setdefault(source_id, []).append(d)
            target_to_values.setdefault(target_id, []).append(d)
    # Distances to closest node in nodes_to (B) from nodes_from (A)
    for source_id in nodes_from:
        inner_values = source_to_values[source_id]
        values.append(numpy.min(inner_values))
    # Distances to closest node in nodes_from (A) from nodes_to (B)
    for target_id in nodes_to:
        inner_values = target_to_values[target_id]
        values.append(numpy.min(inner_values))
    return values


def get_separation_within_set(network, nodes_from, lengths=None):
    """
    Calculate dAA or dBB in separation metric proposed by Menche et al. 2015
    """
    if len(nodes_from) == 1:
        return [ 0 ]
    values = []
    # Distance to closest node within the set (A or B)
    for source_id in nodes_from:
        inner_values = []
        for target_id in nodes_from:
            if source_id == target_id:
                continue
            if lengths is not None:
                d = lengths[source_id][target_id] 
            else:
                d = network_utilities.get_shortest_path_length_between(network, source_id, target_id)
            inner_values.append(d)
        values.append(numpy.min(inner_values))
    return values



### GUILD related ###

def create_node_file(node_to_score, nodes, node_file, background_score = 0.01):
    """
    Simplified method for creating guild node score files
    """
    f = open(node_file, 'w')
    for node in nodes:
        if node in node_to_score:
            score = node_to_score[node]
        else:
            score = background_score
        f.write("%s %f\n" % (node, score))
    f.close()
    return


def run_guild(phenotype, node_to_score, network_nodes, network_file, output_dir, executable_path = None, background_score = 0.01, qname=None, method='s'): 
    # Create node file
    node_file = "%s%s.node" % (output_dir, phenotype)
    create_node_file(node_to_score, network_nodes, node_file, background_score)
    output_file = "%s%s.n%s" % (output_dir, phenotype, method)
    # Get and run the GUILD command
    #print strftime("%H:%M:%S - %d %b %Y") #, score_command
    if method == 's': 
        n_repetition = 3 
        n_iteration = 2 
        score_command = ' -s s -n "%s" -e "%s" -o "%s" -r %d -i %d' % (node_file, network_file, output_file, n_repetition, n_iteration)
    elif method == 'd':
        score_command = ' -s d -n "%s" -e "%s" -o "%s"' % (node_file, network_file, output_file)
    elif method == 'r':
        n_iteration = 50 
        score_command = ' -s r -n "%s" -e "%s" -o "%s" -i %d' % (node_file, network_file, output_file, n_iteration)
    elif method == 'p':
        score_command = ' "%s" "%s" "%s" 1' % (node_file, network_file, output_file)
    elif method == 'w':
        score_command = ' "%s" "%s" "%s"' % (node_file, network_file, output_file)
    else:
        raise NotImplementedError("method %s" % method) 
    if qname is None:
        if executable_path is None:
            if method in ["s", "r", "d"]:
                executable_path = "guild" # assuming accessible guild executable
            else:
                executable_path = "netwalk.sh" # assuming R and netwalk.sh is accessible  
        score_command = executable_path + score_command
        print(score_command)
        os.system(score_command)
    else:
        #os.system("qsub -cwd -o out -e err -q %s -N %s -b y %s" % (qname, scoring_type, score_command))
        #print "qsub -cwd -o out -e err -q %s -N guild_%s -b y %s" % (qname, drug, score_command)
        print("{}".format(score_command.replace('"', '')))
    return score_command


def guildify_multiple(network_file, to_file, output_dir, from_file=None, out_file="guild.txt", method="s", executable_path=None):
    """
    to_file: seeds
    If from_file is not None, returns a dictionary containing average z scores of targets to source, otherwise returns empty dictionary
    method: d | s | r | w | p
    (netshort | netscore | page rank | random walk | propagation)
    """
    if from_file is not None and os.path.exists(out_file):
        target_to_source_score = dict(line.strip("\n").split() for line in open(out_file).readlines())
        return target_to_source_score 
    target_to_source_score = {}
    network = get_network(network_file, only_lcc = True) # using LCC 
    if network_file.endswith(".lcc"):
        network_lcc_file = network_file
    else:
        network_lcc_file = network_file + ".lcc"
    nodes = set(network.nodes())
    disease_to_genes, disease_to_category = get_diseasome_genes(to_file, nodes = nodes)
    if not os.path.exists(output_dir):
        print("Creating output directory: {}".format(output_dir))
        os.makedirs(output_dir)
    # Generate background file (for P-value calculation)
    if not os.path.exists(output_dir + "/background.node"):
        node_to_degree = dict(network.degree())
        n = max(map(len, disease_to_genes.values()))
        values = node_to_degree.items()
        values.sort(key=lambda x: -x[1])
        #k = 1.0 * max(node_to_degree.values())
        values = set(zip(*values[:n])[0])
        f = open(output_dir + "/background.node", 'w')
        for node, degree in node_to_degree.iteritems():
            #score = degree/k
            if node in values: score = 1
            else: score = 0.01
            f.write("%s %f\n" % (node, score))
        f.close()
    if from_file is not None:
        drug_to_targets, drug_to_category = get_diseasome_genes(from_file, nodes = nodes)
        f = open(out_file, 'w')
        f.write("source\ttarget\tscore\n")
    for target, geneids in disease_to_genes.iteritems():
        #print target, len(geneids)
        target_mod = text_utilities.convert_to_R_string(target)
        target_to_score = dict((gene, 1.0) for gene in geneids)
        node_file = output_dir + "%s.n%s" % (target_mod, method) 
        if os.path.exists(node_file):
            print("Skipping existing: {}".format(node_file))
            continue
        run_guild(target_mod, target_to_score, nodes, network_lcc_file, output_dir, executable_path, background_score = 0.01, qname = "print", method = method) #!
        node_to_score = dict(line.strip("\n").split() for line in open(node_file).readlines())
        if from_file is not None:
            values = map(float, numpy.array(node_to_score.values()))
            m = numpy.mean(values)
            s = numpy.std(values)
            for source, geneids in drug_to_targets.iteritems():
                score = -numpy.mean([(float(node_to_score[gene]) - m) / s for gene in geneids])
                f.write("%s\t%s\t%f\n" % (source, target, score))
                d = target_to_source_score.setdefault(target, {})
                d[source] = score
    if from_file is not None:
        f.close()
    return target_to_source_score


def get_scores(score_file):
    """
        Parses scores from a scoring file created by GUILD (node <whitespace> score), returns a dictionary where the values are floats.
    """
    nodes, dummy, node_to_score, dummy = network_utilities.get_nodes_and_edges_from_sif_file(file_name = score_file, store_edge_type = False, delim=None, data_to_float=True)
    return node_to_score



### DIAMOnD related ###

def get_diamond_genes(network_file, seeds, file_name, only_lcc=True):
    network = get_network(network_file, only_lcc=only_lcc) 
    nodes = set(network.nodes())
    seeds = set(seeds) & nodes
    #print len(seeds)
    n_iteration = 500
    if not os.path.exists(file_name):
        diamond.DIAMOnD(network, seeds, n_iteration, alpha = 1, outfile = file_name)
    f = open(file_name)
    f.readline()
    genes = []
    for line in f:
        rank, geneid = line.strip("\n").split()
        genes.append(geneid)
    f.close()
    if not os.path.exists(file_name + ".coverage"):
        f_out = open(file_name + ".coverage", 'w')
        n = float(len(seeds))
        component = network.subgraph(seeds)
        #component = max(networkx.connected_components(component), key=len)
        components = max(network_utilities.get_connected_components(network, False), key=len)
        f_out.write("%s %f\n" % ("0", len(component & seeds)/n))
        for i, gene in enumerate(genes):
            rank = i + 1
            component = network.subgraph(genes[:rank] + list(seeds))
            #component = max(networkx.connected_components(component), key=len)
            components = max(network_utilities.get_connected_components(network, False), key=len)
            f_out.write("%s %f\n" % (rank, len(component & seeds)/n))
        f_out.close()
    return genes, nodes


### Functional enrichment related ###

def check_functional_enrichment(id_list, background_id_weights = None, id_type = "genesymbol", species = "Homo sapiens", mode="unordered", evidences = None, out_file_name = None, tex_format = False):
    """
    id_type = "geneid" # "uniprotacession" # "genesymbol"
    evidences = ['EXP', 'IDA', 'IEP', 'IGI', 'IMP', 'ISA', 'ISM', 'ISO', 'ISS', 'IGC'] # 'IPI'
    evidences = None corresponds to ['EXP', 'IC', 'IDA', 'IEA', 'IEP', 'IGC', 'IGI', 'IMP', 'IPI', 'ISA', 'ISM', 'ISO', 'ISS', 'NAS', 'RCA', 'TAS']
    for custom associations: association = [["GO:0006509", "351"], ["GO:0048167", "348", "5663", "5664", "23621"], ["GO:0097458", "1005", "1006", "1007"], ["GO:0048487", "1", "2", "351"], ["GO:0048488", map(str, range(1000,2000))]]
    for backgroud id weights, such as occurrence frequency: (gene, weight) 
    """
    if out_file_name is not None:
        f_output = open(out_file_name, 'w').write
    else:
        from sys import stdout 
        f_output = stdout.write
    return functional_enrichment.check_functional_enrichment(id_list, background_id_weights, id_type, f_output, species = species, mode = mode, tex_format = tex_format, support = evidences, associations = None) 




<filename>diana/toolbox/wrappers.py
#######################################################################
# Recipies / wrapper functions using toolbox methods for disease,
# drug and network analysis
# e.g. 10/2015
#######################################################################
import diana.toolbox.network_utilities as network_utilities # For DIANA
import diana.toolbox.stat_utilities as stat_utilities # For DIANA
import diana.toolbox.dict_utilities as dict_utilities # For DIANA
import diana.toolbox.text_utilities as text_utilities # For DIANA
import diana.toolbox.TsvReader # For DIANA
import diana.toolbox.functional_enrichment # For DIANA
import diana.toolbox.parse_umls # For DIANA
import diana.toolbox.parse_msigdb # For DIANA
import diana.toolbox.parse_uniprot # For DIANA
import diana.toolbox.parse_ncbi # For DIANA
import diana.toolbox.parse_do # For DIANA
import diana.toolbox.parse_medic # For DIANA
import diana.toolbox.parse_disgenet # For DIANA
import diana.toolbox.parse_drugbank # For DIANA
import diana.toolbox.parse_medi # For DIANA
import diana.toolbox.parse_hetionet # For DIANA
import csv, numpy, os, cPickle
import random
try:
    from toolbox.external.diamond import diamond
except:
    pass # For DIANA
    #print("DIAMOnD not found and thus will not be available!")

##### Id mapping related #####

def get_mapping(file_name, from_column = None, to_column = None, delim=None, one_to_one=True):
    """
    Assumes header, maps first column to the second
    """
    key_to_value = TsvReader.get_from_to_mapping(file_name, from_column = from_column, to_column = to_column, delim=delim, inner_delim = None, filter_column = None, exclude_value = None, include_value = None, one_to_one=one_to_one) 
    return key_to_value

def convert_to_geneid(file_name, id_type, id_mapping_file):
    """
    Expects a file where each line is a gene name / uniprot id
    id_type: symbol | uniprot
    """
    genes = [ line.strip("\n") for line in open(file_name) ]
    if id_type == "symbol":
        geneid_to_name, name_to_geneid = get_geneid_symbol_mapping(id_mapping_file)
    elif id_type == "uniprot":
        name_to_geneid = get_uniprot_to_id(id_mapping_file, uniprot_ids=genes, only_min=True, key_function=int)
    else:
        raise ValueError("Uknown id type: %s" % id_type)
    geneids = set([ name_to_geneid[gene] for gene in genes if gene in name_to_geneid ])
    genes_non = set([ gene for gene in genes if gene not in name_to_geneid ])
    print("Not found genes: {}".format(genes_non))
    return geneids


def get_uniprot_to_id(uniprot_file, uniprot_ids=None, only_min=True, key_function=int):
    """
    uniprot_file = %(data_dir)s/uniprot/idmapping.tab or idmapping.tab.symbol or idmapping.tab.mouse
    Can be used to convert to geneids (key_function=int) as well as symbols (key_function=len) depending on the input file
    """
    uniprot_to_gene = parse_uniprot.get_uniprot_to_geneid(uniprot_file, uniprot_ids, only_min, key_function)
    return uniprot_to_gene


def get_geneid_symbol_mapping(mapping_file):
    """
    id_mapping_file = %(data_dir)s/ncbi/geneid_to_symbol.txt
    """
    geneid_to_name, name_to_geneid = parse_ncbi.get_geneid_symbol_mapping(mapping_file)
    return geneid_to_name, name_to_geneid


def get_mesh_id_mapping(desc_file, rel_file, dump_file = None):
    """ 
    Get all concept id - mesh id mapping (also gets entry names in addition to main header)
    """
    #dump_file = CONFIG.get("umls_dir") + "/mapping.pcl"
    mesh_id_to_name, concept_id_to_mesh_id, mesh_id_to_name_with_synonyms = parse_umls.get_mesh_id_mapping(desc_file, rel_file, dump_file = dump_file)
    return mesh_id_to_name, concept_id_to_mesh_id, mesh_id_to_name_with_synonyms


def get_mesh_disease_ontology(desc_file, rel_file, dump_file = None):
    #dump_file = CONFIG.get("umls_dir") + "/ontology.pcl"
    g = parse_umls.get_mesh_disease_ontology(desc_file, rel_file, dump_file = dump_file)
    return g


def get_medic_mesh_id_mapping(medic_file):
    name_to_id, id_to_mesh_ids = get_medic_mesh_id_mapping(medic_file)
    return name_to_id, id_to_mesh_ids 


def get_mesh_disease_category_mapping(desc_file, rel_file, dump_file = None):
    #dump_file = CONFIG.get("umls_dir") + "/ontology.pcl"
    mesh_id_to_top_ids = parse_umls.get_mesh_id_to_disease_category(desc_file, rel_file, dump_file)
    mesh_id_to_name, concept_id_to_mesh_id, mesh_id_to_name_with_synonyms = get_mesh_id_mapping(desc_file, rel_file, dump_file)
    mesh_name_to_parents = {}
    for child, parents in mesh_id_to_top_ids.iteritems():
        name_child = mesh_id_to_name[concept_id_to_mesh_id[child]]
        values = []
        for parent in parents:
            name_parent = mesh_id_to_name[concept_id_to_mesh_id[parent]]
            values.append(name_parent.lower())
        values.sort()
        mesh_name_to_parents[name_child.lower()] = values
    #print mesh_name_to_parents["asthma"], mesh_name_to_parents["psoriasis"]
    return mesh_name_to_parents


def get_icd_to_mesh_ids(disease_ontology_file, id_type="ICD9CM"):
    icd_to_mesh_ids = parse_do.get_icd_to_mesh_ids(disease_ontology_file, id_type)
    return icd_to_mesh_ids


def get_do_to_mesh_ids(disease_ontology_file):
    name_to_do_id, do_to_mesh_ids, mesh_id_to_type_to_ids = parse_do.get_do_mesh_id_mapping(disease_ontology_file)
    return do_to_mesh_ids


def get_homology_mapping(homologene_file, tax_id="10090", from_tax_id="9606", symbol_type="geneid"):
    """
    symbol_type: geneid | symbol
    """
    geneid_to_geneid, group_to_taxid_to_geneid = parse_ncbi.get_homology_mapping(homologene_file, tax_id, from_tax_id=from_tax_id, symbol_type=symbol_type)
    return geneid_to_geneid 


##### Network related #####

def get_network(network_file, only_lcc):
    network = network_utilities.create_network_from_sif_file(network_file, use_edge_data = False, delim = None, include_unconnected=True)
    #print len(network.nodes()), len(network.edges())
    if only_lcc and not network_file.endswith(".lcc"):
        print("Shrinking network to its LCC: {} nodes, {} edges".format(len(network.nodes()), len(network.edges())))
        components = network_utilities.get_connected_components(network, False)
        network = network_utilities.get_subgraph(network, components[0])
        print("Final shape: {} nodes, {} edges".format(len(network.nodes()), len(network.edges())))
        #print len(network.nodes()), len(network.edges())
        network_lcc_file = network_file + ".lcc"
        if not os.path.exists(network_lcc_file ):
            f = open(network_lcc_file, 'w')
            for u,v in network.edges():
                f.write("%s 1 %s\n" % (u, v))
            f.close()
    return network


def create_functional_network(links_file, mapping_file, cutoff = 900): 
    #string_dir = CONFIG.get("string_dir") + "/"
    #links_file = string_dir + CONFIG.get("string_links_file")
    #mapping_file = string_dir + CONFIG.get("string_mapping_file")
    output_file = CONFIG.get("network_file")
    parse_string.get_interactions(links_file, mapping_file, output_file, cutoff) #, include_score=True) 
    #network = get_network()
    #print len(network.nodes()), len(network.edges())
    return


def calculate_lcc_significance(network, nodes, nodes_random=None, bins=None, n_random=1000, min_bin_size=100, seed=452456):
    # Degree matching problematic for small bin sizes
    #if bins is None and nodes_random is None:
    #        bins = network_utilities.get_degree_binning(network, min_bin_size) 
    random.seed(seed)
    if nodes_random is None:
        network_nodes = list(network.nodes())
        #nodes_random = get_random_nodes(nodes, network, bins = bins, n_random = n_random, min_bin_size = min_bin_size, seed = seed)
        nodes_random = []
        for i in xrange(n_random):
            random.shuffle(network_nodes)
            nodes_random.append(network_nodes[:len(nodes)])
    network_sub = network.subgraph(nodes)
    component_nodes = network_utilities.get_connected_components(network_sub, False)
    #print component_nodes 
    d = len(component_nodes[0])
    values = numpy.empty(len(nodes_random)) 
    for i, nodes in enumerate(nodes_random):
        network_sub = network.subgraph(nodes)
        component_nodes = network_utilities.get_connected_components(network_sub, False)[0]
        values[i] = len(component_nodes)
    m, s = numpy.mean(values), numpy.std(values)
    if s == 0:
        z = 0.0
    else:
        z = (d - m) / s
    return d, z, (m, s), values


##### Gene expression related #####

def get_expression_info(gexp_file, process=None, delim=',', quote='"', R_header=False, dump_file=None):
    """
    To get gene expression info
    process: a set(["log2", "z", "abs"]) or None
    """
    if dump_file is not None and os.path.exists(dump_file):
        gexp, gene_to_idx, cell_line_to_idx = cPickle.load(open(dump_file))
        return gexp, gene_to_idx, cell_line_to_idx
    #gene_to_values = {}
    f = open(gexp_file)
    reader = csv.reader(f, delimiter=delim, quotechar=quote)
    header = reader.next()
    #print len(header), header
    if R_header == False:
        header = header[1:]
    cell_line_to_idx = dict([ (cell_line, i) for i, cell_line in enumerate(header) ])
    gene_to_idx = {}
    values_arr = []
    for i, row in enumerate(reader):
        gene = row[0]
        values = map(float, row[1:])
        #gene_to_values[gene] = values
        gene_to_idx[gene] = i
        values_arr.append(values)
    f.close()
    gexp = numpy.array(values_arr)
    if process is not None:
        if "log2" in process:
            gexp = numpy.log2(gexp)
        if "z" in process:
            gexp = (gexp - gexp.mean(axis=1)[:, numpy.newaxis]) / gexp.std(axis=1, ddof=1)[:, numpy.newaxis]
        if "abs" in process:
            gexp = numpy.abs(gexp)
        #if "na.rm" in process:
        #    idx = numpy.where(numpy.isnan(a)) # need to remove rows with NAs
        #print gexp.shape, gexp_norm.shape
        #print gexp[0,0], gexp_norm[0,0]
        #return gene_to_values, cell_line_to_idx
    if dump_file is not None:
        values = gexp, gene_to_idx, cell_line_to_idx
        cPickle.dump(values, open(dump_file, 'w')) 
    return gexp, gene_to_idx, cell_line_to_idx


def get_de_genes(file_name, cutoff_adj = 0.05, cutoff_logfc=0.585, n_top=None, id_type = "GeneID"):
    """
    For parsing DE file generated using R PEPPER package
    """
    fields_to_include = [id_type, "P.Value", "logFC", "adj.P.Val"]
    parser = TsvReader.TsvReader(file_name, delim="\t", inner_delim=None)
    header_to_idx, id_to_values = parser.read(fields_to_include, keys_to_include=None, merge_inner_values=False)
    if "" in id_to_values:
        del id_to_values[""]
    #print len(id_to_values)
    #gene = "10458"
    #if gene in id_to_values:
    #    print id_to_values[gene] 
    genes = set()
    genes_all = set()
    genes_up, genes_down = set(), set()
    values_gene = []
    for gene, values in id_to_values.iteritems():
        include = False
        positive = False
        for val in values:
            pval = val[header_to_idx["adj.p.val"]] # "p.value"]]
            if pval == "NA":
                continue
            fc = float(val[header_to_idx["logfc"]])
            if float(pval) <= cutoff_adj:
                if abs(fc) >= cutoff_logfc: 
                    include = True
                if fc >= 0: 
                    positive = True
                if n_top is not None:
                    values_gene.append((abs(fc), gene))
        for word in gene.split("///"):
            word = word.strip()
            if word == "---":
                continue
            if include:
                genes.add(word)
                if positive:
                    genes_up.add(word)
                else:
                    genes_down.add(word)
            else:
                genes_all.add(word)
    if n_top is not None:
        values_gene.sort()
        genes = set([ word.strip() for fc, gene in values_gene[-n_top:] for word in gene.split("///") ])
    return genes, genes_all, genes_up, genes_down


def get_z_genes(file_name, cutoff_z = 2):
    """
    For parsing DE-Z file generated using R PEPPER package
    """
    gexp, gene_to_idx, cell_line_to_idx = get_expression_info(file_name, process=None, delim='\t', R_header=True) #, quote='"', dump_file=None)
    genes = gene_to_idx.items()
    genes.sort(key=lambda x: x[1])
    genes = numpy.array(zip(*genes)[0])
    sample_to_genes = {}
    for cell_line, idx in cell_line_to_idx.iteritems():
        indices = numpy.abs(gexp[:,idx]) > cutoff_z 
        sample_to_genes[cell_line] = genes[indices]
        #if cell_line in ["GSM734834", "GSM734833"]: 
        #    print cell_line, len(genes[indices]), genes[indices] 
    return sample_to_genes


def get_sample_mapping(file_name, labels_case, labels_control=None):
    f = open(file_name)
    labels_case = set(labels_case)
    if labels_control is not None:
        labels_control = set(labels_control)
    samples_case = []
    samples_control = []
    for line in f:
        sample, label = line.strip("\n").split("\t")
        label = label.strip()
        if label in labels_case:
            samples_case.append(sample)
        else:
            if labels_control is None or label in labels_control:
                samples_control.append(sample)
    return samples_case, samples_control


##### Disease, pathway, comorbidity, symptom info related #####

def get_pathway_info(pathway_file, prefix=None, nodes=None, max_pathway_size=None, inner_delim=None):
    """
    Assumes a tab separated file containing pathway name, link, geneids
    nodes to filter geneids that are not in the network
    prefix: kegg | reactome | biocarta
    inner_delim: None for tab separated geneids, " " for space separated geneids
    """
    pathway_to_geneids, geneid_to_pathways = parse_msigdb.get_msigdb_info(pathway_file, prefix, inner_delim=inner_delim)
    if nodes is not None or max_pathway_size is not None:
        pathway_to_geneids_mod = {}
        for pathway, geneids in pathway_to_geneids.iteritems():
            if max_pathway_size is not None:
                if len(geneids) > max_pathway_size:
                    continue
            if nodes is not None:
                geneids &= nodes
                if len(geneids) == 0:
                    continue
            pathway_to_geneids_mod[pathway] = geneids
        pathway_to_geneids = pathway_to_geneids_mod
    return pathway_to_geneids


def get_diseasome_genes(diseasome_file, nodes=None, network=None):
    """
    If nodes is not None, keep only nodes in the network
    If network is not None, keep only LCC
    """
    disease_to_genes = {}
    disease_to_category = {}
    for line in open(diseasome_file):
        words = line.strip("\n").split("\t")
        disease = words[1].strip('"')
        category = words[0]
        genes = set(words[2:])
        if nodes is not None:
            genes &= nodes
            if len(genes) == 0:
                continue
        if network is not None:
            network_sub = network.subgraph(genes)
            genes = network_utilities.get_connected_components(network_sub, False)[0]
        disease_to_genes[disease] = genes
        disease_to_category[disease] = category
    return disease_to_genes, disease_to_category


def get_disgenet_genes(file_name):
    disease_to_genes, disease_to_sources, cui_to_disease = parse_disgenet.get_disgenet_genes(file_name)
    return disease_to_genes, disease_to_sources, cui_to_disease 


def get_comorbidity_info(comorbidity_file, disease_ontology_file, mesh_dump, correlation_type="RR", only_significant=False):
    """
    Parse HuDiNe data from AllNet3 and map ICD9 to MeSH using DO
    correlation_type: phi (pearson correlation) | RR (favors rare disease pairs)
    """
    icd_to_mesh_ids = get_icd_to_mesh_ids(disease_ontology_file, id_type="ICD9CM")
    #print len(icd_to_mesh_ids), icd_to_mesh_ids.items()[:5]
    #print [ icd_to_mesh_ids[val] for val in ("289", "502", "578", "579", "542", "543")]
    mesh_id_to_name, concept_id_to_mesh_id, mesh_id_to_name_with_synonyms = get_mesh_id_mapping(None, None, dump_file = mesh_dump)
    #print mesh_id_to_name.items()[:5]
    f = open(comorbidity_file)
    header_to_idx = dict((word, i) for i, word in enumerate(f.readline().strip().split("\t")))
    disease_to_disease_comorbidity = {}
    for line in f:
        words = line.strip().split("\t")
        icd1, icd2 = words[:2]
        if only_significant and significance == "0":
            continue
        if icd1 not in icd_to_mesh_ids or icd2 not in icd_to_mesh_ids:
            #print "Not in DO mapping:", icd1, icd2
            continue
        val = float(words[header_to_idx[correlation_type]]) # idx:5
        for mesh1 in icd_to_mesh_ids[icd1]:
            if mesh1 not in mesh_id_to_name:
                #print "Not in name mapping:", mesh1
                continue
            disease1 = mesh_id_to_name[mesh1].lower()
            for mesh2 in icd_to_mesh_ids[icd2]:
                if mesh2 not in mesh_id_to_name:
                    #print "Not in name mapping:", mesh2
                    continue
                disease2 = mesh_id_to_name[mesh2].lower()
                disease1_mod, disease2_mod = sorted((disease1, disease2))
                d = disease_to_disease_comorbidity.setdefault(disease1, {})
                if disease2 in d:
                    if d[disease2] > val: # skip if the existing comorbidity value is higher
                        continue
                d[disease2] = val 
                d = disease_to_disease_comorbidity.setdefault(disease2, {})
                d[disease1] = val
                #print icd1, mesh1, disease1, icd2, mesh2, disease2, val 
    #print len(disease_to_disease_comorbidity), disease_to_disease_comorbidity.values()[0].items()[:5]
    return disease_to_disease_comorbidity 
    # Parse HuDiNe data from potentially buggy comorbidity_new.tsv
    #comorbidity_file = CONFIG.get("comorbidity_file")
    f = open(comorbidity_file)
    header_to_idx = dict((word, i) for i, word in enumerate(f.readline().strip().split("\t")))
    disease_to_disease_comorbidity = {}
    for line in f:
        words = line.strip().split("\t")
        disease1, disease2 = words[:2]
        significance = words[header_to_idx["sign_"+correlation_type]]
        if only_significant and significance == "0":
            continue
        val = float(words[header_to_idx[correlation_type]])
        disease_to_disease_comorbidity.setdefault(disease1, {})[disease2] = (val, significance)
        disease_to_disease_comorbidity.setdefault(disease2, {})[disease1] = (val, significance)
    f.close()
    return disease_to_disease_comorbidity


def get_symptom_info(symptom_file, tfidf_cutoff=None):
    """
    Parse Zhou et al supplementary s4. A cutoff of 3.5 is likely to filter spurious associations.
    """
    disease_to_symptoms = {}
    symptom_to_diseases = {}
    disease_to_symptom_to_score = {}
    #symptom_file = CONFIG.get("symptom_file")
    f = open(symptom_file)
    f.readline()
    for line in f:
        words = line.strip("\n").split("\t")
        symptom, disease, n, score = words
        symptom = symptom.lower()
        disease = disease.lower()
        if tfidf_cutoff is not None and not float(score) >= tfidf_cutoff:
            continue
        disease_to_symptoms.setdefault(disease, set()).add(symptom)
        symptom_to_diseases.setdefault(symptom, set()).add(disease)
        d = disease_to_symptom_to_score.setdefault(disease, {})
        d[symptom] = float(score)
    return disease_to_symptoms, symptom_to_diseases, disease_to_symptom_to_score


##### Drug related info #####

def get_drugbank(drugbank_file):
    dump_file = drugbank_file + ".pcl"
    if os.path.exists(dump_file):
        parser = cPickle.load(open(dump_file))
    else:
        parser = parse_drugbank.DrugBankXMLParser(drugbank_file)
        parser.parse()
        cPickle.dump(parser, open(dump_file, 'w'))
    return parser


def get_medi_indications(medi_file, drugbank_file, mesh_dump, disease_ontology_file, only_hps=True):
    dump_file = medi_file + ".pcl"
    if os.path.exists(dump_file):
        drug_to_diseases = cPickle.load(open(dump_file))
        return drug_to_diseases 
    parser = get_drugbank(drugbank_file)
    name_to_drug, synonym_to_drug = parser.get_synonyms(selected_drugs=None, only_synonyms=False)
    name_to_icd_and_confidences = parse_medi.get_medi_mapping(medi_file)
    mesh_id_to_name, concept_id_to_mesh_id, mesh_id_to_name_with_synonyms = get_mesh_id_mapping(None, None, dump_file = mesh_dump)
    icd_to_mesh_ids = get_icd_to_mesh_ids(disease_ontology_file, id_type="ICD9CM")
    drug_to_indications = parse_medi.get_drug_disease_mapping(name_to_icd_and_confidences, name_to_drug, synonym_to_drug, icd_to_mesh_ids, mesh_id_to_name, dump_file = None) 
    drug_to_diseases = {}
    for drug, values in drug_to_indications.iteritems():
        for phenotype, dui, val in values:
            if only_hps and val <= 0.5:
                continue
            drug_to_diseases.setdefault(drug, set()).add(phenotype)
    # Drug name to disease name textual mapping
    #name_to_indication_and_confidences = parse_medi.get_medi_mapping(medi_file, textual_indication=True)
    #drug_to_diseases = {}
    #for name, values in name_to_indication_and_confidences.iteritems():
    #        for indication, confidence in values:
    #            if only_hps and confidence <= 0.5:
    #                continue
    #            drug_to_diseases.setdefault(name, set()).add(indication.lower())
    # Get disease to name mapping
    #phenotype_to_mesh_id = dict((name, mesh_id) for mesh_id, name in mesh_id_to_name.iteritems())
    #disease_to_drugs = parse_medi.get_disease_specific_drugs(drug_to_diseases, phenotype_to_mesh_id)
    cPickle.dump(drug_to_diseases, open(dump_file, 'w'))
    return drug_to_diseases


def get_hetionet_indications(hetionet_file, mesh_dump, disease_ontology_file):
    dump_file = hetionet_file + ".pcl"
    if os.path.exists(dump_file):
        drug_to_diseases = cPickle.load(open(dump_file))
        return drug_to_diseases 
    drug_to_do_ids = parse_hetionet.get_hetionet_mapping(hetionet_file, metaedge="CtD")
    do_to_mesh_ids = get_do_to_mesh_ids(disease_ontology_file)
    mesh_id_to_name, concept_id_to_mesh_id, mesh_id_to_name_with_synonyms = get_mesh_id_mapping(None, None, dump_file = mesh_dump)
    drug_to_indications = parse_hetionet.get_drug_disease_mapping(drug_to_do_ids, do_to_mesh_ids, mesh_id_to_name, dump_file = None)
    drug_to_diseases = {}
    for drug, values in drug_to_indications.iteritems():
        for phenotype, dui, val in values:
            drug_to_diseases.setdefault(drug, set()).add(phenotype)
    cPickle.dump(drug_to_diseases, open(dump_file, 'w'))
    return drug_to_diseases 


##### Statistics related #####

def overlap_significance(geneids1, geneids2, nodes, method="hyper"):
    """
    method: hyper(geometric) | fishers (two-sided version of hypergeometric) | jaccard | jaccard_max | overlap
    """
    n1, n2 = len(geneids1), len(geneids2)
    n = len(geneids1 & geneids2)
    N = len(nodes)
    if method == "hyper":
        val = stat_utilities.hypergeometric_test_numeric(n, n1, N, n2)
    elif method == "fishers":
        oddsratio, val = stat_utilities.fisher_exact(n, n1 - n, n2 -n, N - n1 - n2 + n, alternative="two-sided")
    elif method == "jaccard":
        val = stat_utilities.jaccard(geneids1, geneids2)
    elif method == "jaccard_max":
        val = stat_utilities.jaccard_max(geneids1, geneids2)
    elif method == "overlap":
        val = n
    else:
        raise ValueError("Uknown method: %s" % method)
    return n, n1, n2, val


##### Proximity related #####

def calculate_proximity(network, nodes_from, nodes_to, nodes_from_random=None, nodes_to_random=None, bins=None, n_random=1000, min_bin_size=100, seed=452456, lengths=None):
    """
    Calculate proximity from nodes_from to nodes_to
    If degree binning or random nodes are not given, they are generated
    lengths: precalculated shortest path length dictionary
    """
    #distance = "closest"
    #lengths = network_utilities.get_shortest_path_lengths(network, "../data/toy.sif.pcl")
    #d = network_utilities.get_separation(network, lengths, nodes_from, nodes_to, distance, parameters = {})
    nodes_network = set(network.nodes())
    nodes_from = set(nodes_from) & nodes_network 
    nodes_to = set(nodes_to) & nodes_network
    if len(nodes_from) == 0 or len(nodes_to) == 0:
        return None # At least one of the node group not in network
    d = calculate_closest_distance(network, nodes_from, nodes_to, lengths)
    if bins is None and (nodes_from_random is None or nodes_to_random is None):
        bins = network_utilities.get_degree_binning(network, min_bin_size, lengths) # if lengths is given, it will only use those nodes
    if nodes_from_random is None:
        nodes_from_random = get_random_nodes(nodes_from, network, bins = bins, n_random = n_random, min_bin_size = min_bin_size, seed = seed)
    if nodes_to_random is None:
        nodes_to_random = get_random_nodes(nodes_to, network, bins = bins, n_random = n_random, min_bin_size = min_bin_size, seed = seed)
    random_values_list = zip(nodes_from_random, nodes_to_random)
    values = numpy.empty(len(nodes_from_random)) #n_random
    for i, values_random in enumerate(random_values_list):
        nodes_from, nodes_to = values_random
        #values[i] = network_utilities.get_separation(network, lengths, nodes_from, nodes_to, distance, parameters = {})
        values[i] = calculate_closest_distance(network, nodes_from, nodes_to, lengths)
    #pval = float(sum(values <= d)) / len(values) # needs high number of n_random
    m, s = numpy.mean(values), numpy.std(values)
    if s == 0:
        z = 0.0
    else:
        z = (d - m) / s
    return d, z, (m, s) #(z, pval)


def calculate_proximity_multiple(network, from_file=None, to_file=None, n_random=1000, min_bin_size=100, seed=452456, lengths=None, out_file="output.txt"):
    """
    Run proximity on each entries of from and to files in a pairwise manner
    output is saved in out_file (e.g., output.txt)
    """
    nodes = set(network.nodes())
    drug_to_targets, drug_to_category = get_diseasome_genes(from_file, nodes = nodes)
    #drug_to_targets = dict((drug, nodes & targets) for drug, targets in drug_to_targets.iteritems())
    disease_to_genes, disease_to_category = get_diseasome_genes(to_file, nodes = nodes)
    # Calculate proximity values
    print(len(drug_to_targets), len(disease_to_genes))
    # Get degree binning 
    bins = network_utilities.get_degree_binning(network, min_bin_size)
    f = open(out_file, 'w')
    f.write("source\ttarget\tn.source\tn.target\td\tz\n")
    for drug, nodes_from in drug_to_targets.iteritems():
        values = []
        for disease, nodes_to in disease_to_genes.iteritems():
            print(drug, disease)
            d, z, (m, s) = calculate_proximity(network, nodes_from, nodes_to, nodes_from_random=None, nodes_to_random=None, bins=bins, n_random=n_random, min_bin_size=min_bin_size, seed=seed, lengths=lengths)
            values.append((drug, disease, z, len(nodes_from), len(nodes_to), d, m, s))
            #f.write("%s\t%s\t%f\t%f\t%f\t%f\n" % (drug, disease, z, d, m, s))
        values.sort(key=lambda x: x[2])
        for drug, disease, z, k, l, d, m, s in values:
            #f.write("%s\t%s\t%f\t%d\t%d\t%f\t%f\t%f\n" % (drug, disease, z, k, l, d, m, s))
            f.write("%s\t%s\t%d\t%d\t%f\t%f\n" % (drug, disease, k, l, d, z))
    f.close()
    return 


def calculate_closest_distance(network, nodes_from, nodes_to, lengths=None):
    values_outer = []
    if lengths is None:
        for node_from in nodes_from:
            values = []
            for node_to in nodes_to:
                val = network_utilities.get_shortest_path_length_between(network, node_from, node_to)
                values.append(val)
            d = min(values)
            #print d,
            values_outer.append(d)
    else:
        for node_from in nodes_from:
            values = []
            vals = lengths[node_from]
            for node_to in nodes_to:
                val = vals[node_to]
                values.append(val)
            d = min(values)
            values_outer.append(d)
    d = numpy.mean(values_outer)
    #print d
    return d


def get_random_nodes(nodes, network, bins=None, n_random=1000, min_bin_size=100, degree_aware=True, seed=None):
    if bins is None:
        # Get degree bins of the network
        bins = network_utilities.get_degree_binning(network, min_bin_size) 
    nodes_random = network_utilities.pick_random_nodes_matching_selected(network, bins, nodes, n_random, degree_aware, seed=seed) 
    return nodes_random


### Separation related

def calculate_separation_proximity(network, nodes_from, nodes_to, nodes_from_random=None, nodes_to_random=None, bins=None, n_random=1000, min_bin_size=100, seed=452456, lengths=None):
    """
    Calculate proximity from nodes_from to nodes_to
    If degree binning or random nodes are not given, they are generated
    lengths: precalculated shortest path length dictionary
    """
    nodes_network = set(network.nodes())
    if len(set(nodes_from) & nodes_network) == 0 or len(set(nodes_to) & nodes_network) == 0:
        return None # At least one of the node group not in network
    d = get_separation(network, nodes_from, nodes_to, lengths)
    if bins is None and (nodes_from_random is None or nodes_to_random is None):
        bins = network_utilities.get_degree_binning(network, min_bin_size, lengths) # if lengths is given, it will only use those nodes
    if nodes_from_random is None:
        nodes_from_random = get_random_nodes(nodes_from, network, bins = bins, n_random = n_random, min_bin_size = min_bin_size, seed = seed)
    if nodes_to_random is None:
        nodes_to_random = get_random_nodes(nodes_to, network, bins = bins, n_random = n_random, min_bin_size = min_bin_size, seed = seed)
    random_values_list = zip(nodes_from_random, nodes_to_random)
    values = numpy.empty(len(nodes_from_random)) #n_random
    for i, values_random in enumerate(random_values_list):
        nodes_from, nodes_to = values_random
        values[i] = get_separation(network, nodes_from, nodes_to, lengths)
    m, s = numpy.mean(values), numpy.std(values)
    if s == 0:
        z = 0.0
    else:
        z = (d - m) / s
    return d, z, (m, s) #(z, pval)


def get_separation(network, nodes_from, nodes_to, lengths=None):
    dAA = numpy.mean(get_separation_within_set(network, nodes_from, lengths))
    dBB = numpy.mean(get_separation_within_set(network, nodes_to, lengths))
    dAB = numpy.mean(get_separation_between_sets(network, nodes_from, nodes_to, lengths))
    d = dAB - (dAA + dBB) / 2.0
    return d


def get_separation_between_sets(network, nodes_from, nodes_to, lengths=None):
    """
    Calculate dAB in separation metric proposed by Menche et al. 2015
    """
    values = []
    target_to_values = {}
    source_to_values = {}
    for source_id in nodes_from:
        for target_id in nodes_to:
            if lengths is not None:
                d = lengths[source_id][target_id] 
            else:
                d = network_utilities.get_shortest_path_length_between(network, source_id, target_id)
            source_to_values.setdefault(source_id, []).append(d)
            target_to_values.setdefault(target_id, []).append(d)
    # Distances to closest node in nodes_to (B) from nodes_from (A)
    for source_id in nodes_from:
        inner_values = source_to_values[source_id]
        values.append(numpy.min(inner_values))
    # Distances to closest node in nodes_from (A) from nodes_to (B)
    for target_id in nodes_to:
        inner_values = target_to_values[target_id]
        values.append(numpy.min(inner_values))
    return values


def get_separation_within_set(network, nodes_from, lengths=None):
    """
    Calculate dAA or dBB in separation metric proposed by Menche et al. 2015
    """
    if len(nodes_from) == 1:
        return [ 0 ]
    values = []
    # Distance to closest node within the set (A or B)
    for source_id in nodes_from:
        inner_values = []
        for target_id in nodes_from:
            if source_id == target_id:
                continue
            if lengths is not None:
                d = lengths[source_id][target_id] 
            else:
                d = network_utilities.get_shortest_path_length_between(network, source_id, target_id)
            inner_values.append(d)
        values.append(numpy.min(inner_values))
    return values



### GUILD related ###

def create_node_file(node_to_score, nodes, node_file, background_score = 0.01):
    """
    Simplified method for creating guild node score files
    """
    f = open(node_file, 'w')
    for node in nodes:
        if node in node_to_score:
            score = node_to_score[node]
        else:
            score = background_score
        f.write("%s %f\n" % (node, score))
    f.close()
    return


def run_guild(phenotype, node_to_score, network_nodes, network_file, output_dir, executable_path = None, background_score = 0.01, qname=None, method='s'): 
    # Create node file
    node_file = "%s%s.node" % (output_dir, phenotype)
    create_node_file(node_to_score, network_nodes, node_file, background_score)
    output_file = "%s%s.n%s" % (output_dir, phenotype, method)
    # Get and run the GUILD command
    #print strftime("%H:%M:%S - %d %b %Y") #, score_command
    if method == 's': 
        n_repetition = 3 
        n_iteration = 2 
        score_command = ' -s s -n "%s" -e "%s" -o "%s" -r %d -i %d' % (node_file, network_file, output_file, n_repetition, n_iteration)
    elif method == 'd':
        score_command = ' -s d -n "%s" -e "%s" -o "%s"' % (node_file, network_file, output_file)
    elif method == 'r':
        n_iteration = 50 
        score_command = ' -s r -n "%s" -e "%s" -o "%s" -i %d' % (node_file, network_file, output_file, n_iteration)
    elif method == 'p':
        score_command = ' "%s" "%s" "%s" 1' % (node_file, network_file, output_file)
    elif method == 'w':
        score_command = ' "%s" "%s" "%s"' % (node_file, network_file, output_file)
    else:
        raise NotImplementedError("method %s" % method) 
    if qname is None:
        if executable_path is None:
            if method in ["s", "r", "d"]:
                executable_path = "guild" # assuming accessible guild executable
            else:
                executable_path = "netwalk.sh" # assuming R and netwalk.sh is accessible  
        score_command = executable_path + score_command
        print(score_command)
        os.system(score_command)
    else:
        #os.system("qsub -cwd -o out -e err -q %s -N %s -b y %s" % (qname, scoring_type, score_command))
        #print "qsub -cwd -o out -e err -q %s -N guild_%s -b y %s" % (qname, drug, score_command)
        print("{}".format(score_command.replace('"', '')))
    return score_command


def guildify_multiple(network_file, to_file, output_dir, from_file=None, out_file="guild.txt", method="s", executable_path=None):
    """
    to_file: seeds
    If from_file is not None, returns a dictionary containing average z scores of targets to source, otherwise returns empty dictionary
    method: d | s | r | w | p
    (netshort | netscore | page rank | random walk | propagation)
    """
    if from_file is not None and os.path.exists(out_file):
        target_to_source_score = dict(line.strip("\n").split() for line in open(out_file).readlines())
        return target_to_source_score 
    target_to_source_score = {}
    network = get_network(network_file, only_lcc = True) # using LCC 
    if network_file.endswith(".lcc"):
        network_lcc_file = network_file
    else:
        network_lcc_file = network_file + ".lcc"
    nodes = set(network.nodes())
    disease_to_genes, disease_to_category = get_diseasome_genes(to_file, nodes = nodes)
    if not os.path.exists(output_dir):
        print("Creating output directory: {}".format(output_dir))
        os.makedirs(output_dir)
    # Generate background file (for P-value calculation)
    if not os.path.exists(output_dir + "/background.node"):
        node_to_degree = dict(network.degree())
        n = max(map(len, disease_to_genes.values()))
        values = node_to_degree.items()
        values.sort(key=lambda x: -x[1])
        #k = 1.0 * max(node_to_degree.values())
        values = set(zip(*values[:n])[0])
        f = open(output_dir + "/background.node", 'w')
        for node, degree in node_to_degree.iteritems():
            #score = degree/k
            if node in values: score = 1
            else: score = 0.01
            f.write("%s %f\n" % (node, score))
        f.close()
    if from_file is not None:
        drug_to_targets, drug_to_category = get_diseasome_genes(from_file, nodes = nodes)
        f = open(out_file, 'w')
        f.write("source\ttarget\tscore\n")
    for target, geneids in disease_to_genes.iteritems():
        #print target, len(geneids)
        target_mod = text_utilities.convert_to_R_string(target)
        target_to_score = dict((gene, 1.0) for gene in geneids)
        node_file = output_dir + "%s.n%s" % (target_mod, method) 
        if os.path.exists(node_file):
            print("Skipping existing: {}".format(node_file))
            continue
        run_guild(target_mod, target_to_score, nodes, network_lcc_file, output_dir, executable_path, background_score = 0.01, qname = "print", method = method) #!
        node_to_score = dict(line.strip("\n").split() for line in open(node_file).readlines())
        if from_file is not None:
            values = map(float, numpy.array(node_to_score.values()))
            m = numpy.mean(values)
            s = numpy.std(values)
            for source, geneids in drug_to_targets.iteritems():
                score = -numpy.mean([(float(node_to_score[gene]) - m) / s for gene in geneids])
                f.write("%s\t%s\t%f\n" % (source, target, score))
                d = target_to_source_score.setdefault(target, {})
                d[source] = score
    if from_file is not None:
        f.close()
    return target_to_source_score


def get_scores(score_file):
    """
        Parses scores from a scoring file created by GUILD (node <whitespace> score), returns a dictionary where the values are floats.
    """
    nodes, dummy, node_to_score, dummy = network_utilities.get_nodes_and_edges_from_sif_file(file_name = score_file, store_edge_type = False, delim=None, data_to_float=True)
    return node_to_score



### DIAMOnD related ###

def get_diamond_genes(network_file, seeds, file_name, only_lcc=True):
    network = get_network(network_file, only_lcc=only_lcc) 
    nodes = set(network.nodes())
    seeds = set(seeds) & nodes
    #print len(seeds)
    n_iteration = 500
    if not os.path.exists(file_name):
        diamond.DIAMOnD(network, seeds, n_iteration, alpha = 1, outfile = file_name)
    f = open(file_name)
    f.readline()
    genes = []
    for line in f:
        rank, geneid = line.strip("\n").split()
        genes.append(geneid)
    f.close()
    if not os.path.exists(file_name + ".coverage"):
        f_out = open(file_name + ".coverage", 'w')
        n = float(len(seeds))
        component = network.subgraph(seeds)
        #component = max(networkx.connected_components(component), key=len)
        components = max(network_utilities.get_connected_components(network, False), key=len)
        f_out.write("%s %f\n" % ("0", len(component & seeds)/n))
        for i, gene in enumerate(genes):
            rank = i + 1
            component = network.subgraph(genes[:rank] + list(seeds))
            #component = max(networkx.connected_components(component), key=len)
            components = max(network_utilities.get_connected_components(network, False), key=len)
            f_out.write("%s %f\n" % (rank, len(component & seeds)/n))
        f_out.close()
    return genes, nodes


### Functional enrichment related ###

def check_functional_enrichment(id_list, background_id_weights = None, id_type = "genesymbol", species = "Homo sapiens", mode="unordered", evidences = None, out_file_name = None, tex_format = False):
    """
    id_type = "geneid" # "uniprotacession" # "genesymbol"
    evidences = ['EXP', 'IDA', 'IEP', 'IGI', 'IMP', 'ISA', 'ISM', 'ISO', 'ISS', 'IGC'] # 'IPI'
    evidences = None corresponds to ['EXP', 'IC', 'IDA', 'IEA', 'IEP', 'IGC', 'IGI', 'IMP', 'IPI', 'ISA', 'ISM', 'ISO', 'ISS', 'NAS', 'RCA', 'TAS']
    for custom associations: association = [["GO:0006509", "351"], ["GO:0048167", "348", "5663", "5664", "23621"], ["GO:0097458", "1005", "1006", "1007"], ["GO:0048487", "1", "2", "351"], ["GO:0048488", map(str, range(1000,2000))]]
    for backgroud id weights, such as occurrence frequency: (gene, weight) 
    """
    if out_file_name is not None:
        f_output = open(out_file_name, 'w').write
    else:
        from sys import stdout 
        f_output = stdout.write
    return functional_enrichment.check_functional_enrichment(id_list, background_id_weights, id_type, f_output, species = species, mode = mode, tex_format = tex_format, support = evidences, associations = None) 




####################################################################### # Recipies / wrapper functions using toolbox methods for disease, # drug and network analysis # e.g. 10/2015 ####################################################################### # For DIANA # For DIANA # For DIANA # For DIANA # For DIANA # For DIANA # For DIANA # For DIANA # For DIANA # For DIANA # For DIANA # For DIANA # For DIANA # For DIANA # For DIANA # For DIANA # For DIANA #print("DIAMOnD not found and thus will not be available!") ##### Id mapping related ##### Assumes header, maps first column to the second Expects a file where each line is a gene name / uniprot id     id_type: symbol | uniprot uniprot_file = %(data_dir)s/uniprot/idmapping.tab or idmapping.tab.symbol or idmapping.tab.mouse     Can be used to convert to geneids (key_function=int) as well as symbols (key_function=len) depending on the input file id_mapping_file = %(data_dir)s/ncbi/geneid_to_symbol.txt Get all concept id - mesh id mapping (also gets entry names in addition to main header) #dump_file = CONFIG.get("umls_dir") + "/mapping.pcl" #dump_file = CONFIG.get("umls_dir") + "/ontology.pcl" #dump_file = CONFIG.get("umls_dir") + "/ontology.pcl" #print mesh_name_to_parents["asthma"], mesh_name_to_parents["psoriasis"] symbol_type: geneid | symbol ##### Network related ##### #print len(network.nodes()), len(network.edges()) #print len(network.nodes()), len(network.edges()) #string_dir = CONFIG.get("string_dir") + "/" #links_file = string_dir + CONFIG.get("string_links_file") #mapping_file = string_dir + CONFIG.get("string_mapping_file") #, include_score=True) #network = get_network() #print len(network.nodes()), len(network.edges()) # Degree matching problematic for small bin sizes #if bins is None and nodes_random is None: #        bins = network_utilities.get_degree_binning(network, min_bin_size) #nodes_random = get_random_nodes(nodes, network, bins = bins, n_random = n_random, min_bin_size = min_bin_size, seed = seed) #print component_nodes ##### Gene expression related ##### To get gene expression info     process: a set(["log2", "z", "abs"]) or None #gene_to_values = {} #print len(header), header #gene_to_values[gene] = values #if "na.rm" in process: #    idx = numpy.where(numpy.isnan(a)) # need to remove rows with NAs #print gexp.shape, gexp_norm.shape #print gexp[0,0], gexp_norm[0,0] #return gene_to_values, cell_line_to_idx For parsing DE file generated using R PEPPER package #print len(id_to_values) #gene = "10458" #if gene in id_to_values: #    print id_to_values[gene] # "p.value"]] For parsing DE-Z file generated using R PEPPER package #, quote='"', dump_file=None) #if cell_line in ["GSM734834", "GSM734833"]: #    print cell_line, len(genes[indices]), genes[indices] ##### Disease, pathway, comorbidity, symptom info related ##### Assumes a tab separated file containing pathway name, link, geneids     nodes to filter geneids that are not in the network     prefix: kegg | reactome | biocarta     inner_delim: None for tab separated geneids, " " for space separated geneids If nodes is not None, keep only nodes in the network     If network is not None, keep only LCC Parse HuDiNe data from AllNet3 and map ICD9 to MeSH using DO     correlation_type: phi (pearson correlation) | RR (favors rare disease pairs) #print len(icd_to_mesh_ids), icd_to_mesh_ids.items()[:5] #print [ icd_to_mesh_ids[val] for val in ("289", "502", "578", "579", "542", "543")] #print mesh_id_to_name.items()[:5] #print "Not in DO mapping:", icd1, icd2 # idx:5 #print "Not in name mapping:", mesh1 #print "Not in name mapping:", mesh2 # skip if the existing comorbidity value is higher #print icd1, mesh1, disease1, icd2, mesh2, disease2, val #print len(disease_to_disease_comorbidity), disease_to_disease_comorbidity.values()[0].items()[:5] # Parse HuDiNe data from potentially buggy comorbidity_new.tsv #comorbidity_file = CONFIG.get("comorbidity_file") Parse Zhou et al supplementary s4. A cutoff of 3.5 is likely to filter spurious associations. #symptom_file = CONFIG.get("symptom_file") ##### Drug related info ##### # Drug name to disease name textual mapping #name_to_indication_and_confidences = parse_medi.get_medi_mapping(medi_file, textual_indication=True) #drug_to_diseases = {} #for name, values in name_to_indication_and_confidences.iteritems(): #        for indication, confidence in values: #            if only_hps and confidence <= 0.5: #                continue #            drug_to_diseases.setdefault(name, set()).add(indication.lower()) # Get disease to name mapping #phenotype_to_mesh_id = dict((name, mesh_id) for mesh_id, name in mesh_id_to_name.iteritems()) #disease_to_drugs = parse_medi.get_disease_specific_drugs(drug_to_diseases, phenotype_to_mesh_id) ##### Statistics related ##### method: hyper(geometric) | fishers (two-sided version of hypergeometric) | jaccard | jaccard_max | overlap ##### Proximity related ##### Calculate proximity from nodes_from to nodes_to     If degree binning or random nodes are not given, they are generated     lengths: precalculated shortest path length dictionary #distance = "closest" #lengths = network_utilities.get_shortest_path_lengths(network, "../data/toy.sif.pcl") #d = network_utilities.get_separation(network, lengths, nodes_from, nodes_to, distance, parameters = {}) # At least one of the node group not in network # if lengths is given, it will only use those nodes #n_random #values[i] = network_utilities.get_separation(network, lengths, nodes_from, nodes_to, distance, parameters = {}) #pval = float(sum(values <= d)) / len(values) # needs high number of n_random #(z, pval) Run proximity on each entries of from and to files in a pairwise manner     output is saved in out_file (e.g., output.txt) #drug_to_targets = dict((drug, nodes & targets) for drug, targets in drug_to_targets.iteritems()) # Calculate proximity values # Get degree binning #f.write("%s\t%s\t%f\t%f\t%f\t%f\n" % (drug, disease, z, d, m, s)) #f.write("%s\t%s\t%f\t%d\t%d\t%f\t%f\t%f\n" % (drug, disease, z, k, l, d, m, s)) #print d, #print d # Get degree bins of the network ### Separation related Calculate proximity from nodes_from to nodes_to     If degree binning or random nodes are not given, they are generated     lengths: precalculated shortest path length dictionary # At least one of the node group not in network # if lengths is given, it will only use those nodes #n_random #(z, pval) Calculate dAB in separation metric proposed by Menche et al. 2015 # Distances to closest node in nodes_to (B) from nodes_from (A) # Distances to closest node in nodes_from (A) from nodes_to (B) Calculate dAA or dBB in separation metric proposed by Menche et al. 2015 # Distance to closest node within the set (A or B) ### GUILD related ### Simplified method for creating guild node score files # Create node file # Get and run the GUILD command #print strftime("%H:%M:%S - %d %b %Y") #, score_command # assuming accessible guild executable # assuming R and netwalk.sh is accessible #os.system("qsub -cwd -o out -e err -q %s -N %s -b y %s" % (qname, scoring_type, score_command)) #print "qsub -cwd -o out -e err -q %s -N guild_%s -b y %s" % (qname, drug, score_command) to_file: seeds     If from_file is not None, returns a dictionary containing average z scores of targets to source, otherwise returns empty dictionary     method: d | s | r | w | p     (netshort | netscore | page rank | random walk | propagation) # using LCC # Generate background file (for P-value calculation) #k = 1.0 * max(node_to_degree.values()) #score = degree/k #print target, len(geneids) #! Parses scores from a scoring file created by GUILD (node <whitespace> score), returns a dictionary where the values are floats. ### DIAMOnD related ### #print len(seeds) #component = max(networkx.connected_components(component), key=len) #component = max(networkx.connected_components(component), key=len) ### Functional enrichment related ### id_type = "geneid" # "uniprotacession" # "genesymbol"     evidences = ['EXP', 'IDA', 'IEP', 'IGI', 'IMP', 'ISA', 'ISM', 'ISO', 'ISS', 'IGC'] # 'IPI'     evidences = None corresponds to ['EXP', 'IC', 'IDA', 'IEA', 'IEP', 'IGC', 'IGI', 'IMP', 'IPI', 'ISA', 'ISM', 'ISO', 'ISS', 'NAS', 'RCA', 'TAS']     for custom associations: association = [["GO:0006509", "351"], ["GO:0048167", "348", "5663", "5664", "23621"], ["GO:0097458", "1005", "1006", "1007"], ["GO:0048487", "1", "2", "351"], ["GO:0048488", map(str, range(1000,2000))]]     for backgroud id weights, such as occurrence frequency: (gene, weight)