<jupyter_start><jupyter_text>Distances<jupyter_code># Preparation
PPI_nx = nx.read_edgelist(f"{PPI_DIRECTORY}/BioGRID_sc.txt")
GDV = graco.orbits(PPI_nx)
GCV = graco.coefficients(GDV)<jupyter_output><empty_output><jupyter_text>TVD  Classic GCV<jupyter_code>for tvd in ['0', '1', '2', '3']:
    t1 = time.time()
    D = cdist(np.array(GCV['-1'][tvd]), np.array(GCV['-1'][tvd]), graco.functions.tvd)
    t2 = time.time()
    print(f'{tvd}: {t2-t1:.2f}sec')  
    np.savetxt(f"{MATRIX_DIRECTORY}/sc_BioGRID_tvd{tvd}.txt", D, 
               fmt='%.7f', header=' '.join(PPI_nx), comments='')<jupyter_output>0: 157.81sec
1: 157.68sec
2: 161.01sec
3: 161.96sec<jupyter_text>Higher order GCV<jupyter_code>t1 = time.time()

for (a,b) in [(0,0), (1,1), (3,3), (1,2), (2,1)]:
    D = cdist(np.array(GCV[str(a)][str(b)]), np.array(GCV[str(a)][str(b)]), graco.functions.tvd)
    t2 = time.time()
    print(f'{a}-{b}:{t2-t1:.2f}sec')  
    np.savetxt(f"{MATRIX_DIRECTORY}/sc_BioGRID_{a}tvd{b}.txt", D, 
               fmt='%.7f', header=' '.join(PPI_nx), comments='')<jupyter_output>0-0:158.06sec
1-1:321.87sec
3-3:486.69sec
1-2:650.43sec
2-1:814.10sec<jupyter_text>Hellinger<jupyter_code>_SQRT2 = np.sqrt(2) 

def hellinger(p, q):
    return np.sqrt(np.sum((np.sqrt(p) - np.sqrt(q)) ** 2)) / _SQRT2<jupyter_output><empty_output><jupyter_text>Classical GCV<jupyter_code>for tvd in ['0', '1', '2', '3']:
    t1 = time.time()
    D = cdist(np.array(GCV['-1'][tvd]), np.array(GCV['-1'][tvd]), hellinger)
    t2 = time.time()
    print(f'{tvd}: {t2-t1:.2f}sec')  
    np.savetxt(f"{MATRIX_DIRECTORY}/sc_BioGRID_hell{tvd}.txt", D, 
               fmt='%.7f', header=' '.join(PPI_nx), comments='')<jupyter_output>0: 233.81sec
1: 231.59sec
2: 231.01sec
3: 232.45sec<jupyter_text>Higher order GCV<jupyter_code>t1 = time.time()
for (a,b) in [(0,0), (1,1), (3,3), (1,2), (2,1)]:
    D = cdist(np.array(GCV[str(a)][str(b)]), np.array(GCV[str(a)][str(b)]), hellinger)
    t2 = time.time()
    print(f'{a}-{b}: {t2-t1:.2f}sec')  
    np.savetxt(f"{MATRIX_DIRECTORY}/sc_BioGRID_{a}hell{b}.txt", D, 
               fmt='%.7f', header=' '.join(PPI_nx), comments='')<jupyter_output>0-0: 230.43sec
1-1: 473.35sec
3-3: 718.71sec
1-2: 956.85sec
2-1: 1200.53sec<jupyter_text>Mahalanobis<jupyter_code>gcv = GCV.fillna(GCV.mean())

t1 = time.time()
D = cdist(np.array(gcv), np.array(gcv), 'mahalanobis')
t2 = time.time()
print(f'{t2-t1:.2f}sec')  
np.savetxt(f"{MATRIX_DIRECTORY}/sc_BioGRID_GCV_malahanobis_mean.txt", D, 
           fmt='%.7f', header=' '.join(PPI_nx), comments='')
gcv = GCV.T.dropna().T

t1 = time.time()
D = cdist(np.array(gcv), np.array(gcv), 'mahalanobis')
t2 = time.time()
print(f'{t2-t1:.2f}sec')  
np.savetxt(f"{MATRIX_DIRECTORY}/sc_BioGRID_malahanobisNa.txt", D, 
           fmt='%.7f', header=' '.join(PPI_nx), comments='')<jupyter_output>1.84sec<jupyter_text>Clustering<jupyter_code>METHOD = "kmedoids"

if not os.path.exists(f"{CLUSTERS_DIRECTORY}/{METHOD}"):
    os.makedirs(f"{CLUSTERS_DIRECTORY}/{METHOD}")<jupyter_output><empty_output><jupyter_text>Mahalanobis<jupyter_code># Individual
MATRIX_NAME = "sc_BioGRID_GCV_malahanobis_mean"
D_df = pd.read_csv(f"{MATRIX_DIRECTORY}/{MATRIX_NAME}.txt", delimiter=' ')

D_df.index = D_df.columns
d = ~D_df.isna().all()
D_df = D_df[d].T[d]

int2gene = dict(enumerate(D_df.columns))

t1 = time.time()
for n_clusters in range(2,100):
    initial_medoids = range(n_clusters)
    kmedoids_instance = kmedoids(np.array(D_df), initial_medoids, data_type='distance_matrix')
    kmedoids_instance.process()

    with open(f"{CLUSTERS_DIRECTORY}/{METHOD}/{MATRIX_NAME}_{n_clusters}.txt", 'w') as f:
        for cluster in kmedoids_instance.get_clusters():
            f.write(' '.join(map(int2gene.get,cluster)) + '\n')
    t2 = time.time()
    print(f'{n_clusters}: {t2-t1:.2f}sec', end='\r')
# Individual
MATRIX_NAME = "sc_BioGRID_malahanobisNa"
D_df = pd.read_csv(f"{MATRIX_DIRECTORY}/{MATRIX_NAME}.txt", delimiter=' ')

D_df.index = D_df.columns
d = ~D_df.isna().all()
D_df = D_df[d].T[d]

int2gene = dict(enumerate(D_df.columns))

t1 = time.time()
for n_clusters in range(2,100):
    initial_medoids = range(n_clusters)
    kmedoids_instance = kmedoids(np.array(D_df), initial_medoids, data_type='distance_matrix')
    kmedoids_instance.process()

    with open(f"{CLUSTERS_DIRECTORY}/{METHOD}/{MATRIX_NAME}_{n_clusters}.txt", 'w') as f:
        for cluster in kmedoids_instance.get_clusters():
            f.write(' '.join(map(int2gene.get,cluster)) + '\n')
    t2 = time.time()
    print(f'{n_clusters}: {t2-t1:.2f}sec', end='\r')<jupyter_output>99: 602.76sec<jupyter_text>TVD<jupyter_code># Individual
MATRIX_NAME = "sc_BioGRID_all_tvd"
D0_df = pd.read_csv(f"{MATRIX_DIRECTORY}/sc_BioGRID_tvd0.txt", delimiter=' ')
D1_df = pd.read_csv(f"{MATRIX_DIRECTORY}/sc_BioGRID_tvd1.txt", delimiter=' ')
D2_df = pd.read_csv(f"{MATRIX_DIRECTORY}/sc_BioGRID_tvd2.txt", delimiter=' ')
D3_df = pd.read_csv(f"{MATRIX_DIRECTORY}/sc_BioGRID_tvd3.txt", delimiter=' ')

D00_df = pd.read_csv(f"{MATRIX_DIRECTORY}/sc_BioGRID_0tvd0.txt", delimiter=' ')
D11_df = pd.read_csv(f"{MATRIX_DIRECTORY}/sc_BioGRID_1tvd1.txt", delimiter=' ')
D33_df = pd.read_csv(f"{MATRIX_DIRECTORY}/sc_BioGRID_3tvd3.txt", delimiter=' ')
D12_df = pd.read_csv(f"{MATRIX_DIRECTORY}/sc_BioGRID_1tvd2.txt", delimiter=' ')
D21_df = pd.read_csv(f"{MATRIX_DIRECTORY}/sc_BioGRID_2tvd1.txt", delimiter=' ')

D = np.nanmean([D0_df.values, D1_df.values, D2_df.values, D3_df.values,
                D00_df.values, D11_df.values, D33_df.values, D12_df.values, D21_df.values], axis=0)
D_df = pd.DataFrame(D, index=D0_df.columns, columns=D0_df.columns)
D_df.index = D_df.columns
d = ~D_df.isna().all()
D_df = D_df[d].T[d]

int2gene = dict(enumerate(D_df.columns))

t1 = time.time()
for n_clusters in range(2,100):
    initial_medoids = range(n_clusters)
    kmedoids_instance = kmedoids(np.array(D_df), initial_medoids, data_type='distance_matrix')
    kmedoids_instance.process()

    with open(f"{CLUSTERS_DIRECTORY}/{METHOD}/{MATRIX_NAME}_{n_clusters}.txt", 'w') as f:
        for cluster in kmedoids_instance.get_clusters():
            f.write(' '.join(map(int2gene.get,cluster)) + '\n')
    t2 = time.time()
    print(f'{n_clusters}: {t2-t1:.2f}sec', end='\r')
# Individual
MATRIX_NAME = "sc_BioGRID_0tvd0"
D_df = pd.read_csv(f"{MATRIX_DIRECTORY}/{MATRIX_NAME}.txt", delimiter=' ')

D_df.index = D_df.columns
d = ~D_df.isna().all()
D_df = D_df[d].T[d]

int2gene = dict(enumerate(D_df.columns))

t1 = time.time()
for n_clusters in range(2,100):
    initial_medoids = range(n_clusters)
    kmedoids_instance = kmedoids(np.array(D_df), initial_medoids, data_type='distance_matrix')
    kmedoids_instance.process()

    with open(f"{CLUSTERS_DIRECTORY}/{METHOD}/{MATRIX_NAME}_{n_clusters}.txt", 'w') as f:
        for cluster in kmedoids_instance.get_clusters():
            f.write(' '.join(map(int2gene.get,cluster)) + '\n')
    t2 = time.time()
    print(f'{n_clusters}: {t2-t1:.2f}sec', end='\r')
# Individual
MATRIX_NAME = "sc_BioGRID_1tvd1"
D_df = pd.read_csv(f"{MATRIX_DIRECTORY}/{MATRIX_NAME}.txt", delimiter=' ')

D_df.index = D_df.columns
d = ~D_df.isna().all()
D_df = D_df[d].T[d]

int2gene = dict(enumerate(D_df.columns))

t1 = time.time()
for n_clusters in range(2,100):
    initial_medoids = range(n_clusters)
    kmedoids_instance = kmedoids(np.array(D_df), initial_medoids, data_type='distance_matrix')
    kmedoids_instance.process()

    with open(f"{CLUSTERS_DIRECTORY}/{METHOD}/{MATRIX_NAME}_{n_clusters}.txt", 'w') as f:
        for cluster in kmedoids_instance.get_clusters():
            f.write(' '.join(map(int2gene.get,cluster)) + '\n')
    t2 = time.time()
    print(f'{n_clusters}: {t2-t1:.2f}sec', end='\r')
# Individual
MATRIX_NAME = "sc_BioGRID_3tvd3"
D_df = pd.read_csv(f"{MATRIX_DIRECTORY}/{MATRIX_NAME}.txt", delimiter=' ')

D_df.index = D_df.columns
d = ~D_df.isna().all()
D_df = D_df[d].T[d]

int2gene = dict(enumerate(D_df.columns))

t1 = time.time()
for n_clusters in range(2,100):
    initial_medoids = range(n_clusters)
    kmedoids_instance = kmedoids(np.array(D_df), initial_medoids, data_type='distance_matrix')
    kmedoids_instance.process()

    with open(f"{CLUSTERS_DIRECTORY}/{METHOD}/{MATRIX_NAME}_{n_clusters}.txt", 'w') as f:
        for cluster in kmedoids_instance.get_clusters():
            f.write(' '.join(map(int2gene.get,cluster)) + '\n')
    t2 = time.time()
    print(f'{n_clusters}: {t2-t1:.2f}sec', end='\r')
# Individual
MATRIX_NAME = "sc_BioGRID_1tvd2"
D_df = pd.read_csv(f"{MATRIX_DIRECTORY}/{MATRIX_NAME}.txt", delimiter=' ')

D_df.index = D_df.columns
d = ~D_df.isna().all()
D_df = D_df[d].T[d]

int2gene = dict(enumerate(D_df.columns))

t1 = time.time()
for n_clusters in range(2,100):
    initial_medoids = range(n_clusters)
    kmedoids_instance = kmedoids(np.array(D_df), initial_medoids, data_type='distance_matrix')
    kmedoids_instance.process()

    with open(f"{CLUSTERS_DIRECTORY}/{METHOD}/{MATRIX_NAME}_{n_clusters}.txt", 'w') as f:
        for cluster in kmedoids_instance.get_clusters():
            f.write(' '.join(map(int2gene.get,cluster)) + '\n')
    t2 = time.time()
    print(f'{n_clusters}: {t2-t1:.2f}sec', end='\r')
# Individual
MATRIX_NAME = "sc_BioGRID_2tvd1"
D_df = pd.read_csv(f"{MATRIX_DIRECTORY}/{MATRIX_NAME}.txt", delimiter=' ')

D_df.index = D_df.columns
d = ~D_df.isna().all()
D_df = D_df[d].T[d]

int2gene = dict(enumerate(D_df.columns))

t1 = time.time()
for n_clusters in range(2,100):
    initial_medoids = range(n_clusters)
    kmedoids_instance = kmedoids(np.array(D_df), initial_medoids, data_type='distance_matrix')
    kmedoids_instance.process()

    with open(f"{CLUSTERS_DIRECTORY}/{METHOD}/{MATRIX_NAME}_{n_clusters}.txt", 'w') as f:
        for cluster in kmedoids_instance.get_clusters():
            f.write(' '.join(map(int2gene.get,cluster)) + '\n')
    t2 = time.time()
    print(f'{n_clusters}: {t2-t1:.2f}sec', end='\r')<jupyter_output>99: 468.56sec<jupyter_text>Hellinger<jupyter_code># Individual
MATRIX_NAME = "sc_BioGRID_all_hell"
D0_df = pd.read_csv(f"{MATRIX_DIRECTORY}/sc_BioGRID_hell0.txt", delimiter=' ')
D1_df = pd.read_csv(f"{MATRIX_DIRECTORY}/sc_BioGRID_hell1.txt", delimiter=' ')
D2_df = pd.read_csv(f"{MATRIX_DIRECTORY}/sc_BioGRID_hell2.txt", delimiter=' ')
D3_df = pd.read_csv(f"{MATRIX_DIRECTORY}/sc_BioGRID_hell3.txt", delimiter=' ')

D00_df = pd.read_csv(f"{MATRIX_DIRECTORY}/sc_BioGRID_0hell0.txt", delimiter=' ')
D11_df = pd.read_csv(f"{MATRIX_DIRECTORY}/sc_BioGRID_1hell1.txt", delimiter=' ')
D33_df = pd.read_csv(f"{MATRIX_DIRECTORY}/sc_BioGRID_3hell3.txt", delimiter=' ')
D12_df = pd.read_csv(f"{MATRIX_DIRECTORY}/sc_BioGRID_1hell2.txt", delimiter=' ')
D21_df = pd.read_csv(f"{MATRIX_DIRECTORY}/sc_BioGRID_2hell1.txt", delimiter=' ')

D = np.nanmean([D0_df.values, D1_df.values, D2_df.values, D3_df.values,
                D00_df.values, D11_df.values, D33_df.values, D12_df.values, D21_df.values], axis=0)

D_df = pd.DataFrame(D, index=D0_df.columns, columns=D0_df.columns)
D_df.index = D_df.columns
d = ~D_df.isna().all()
D_df = D_df[d].T[d]

int2gene = dict(enumerate(D_df.columns))

t1 = time.time()
for n_clusters in range(2,100):
    initial_medoids = range(n_clusters)
    kmedoids_instance = kmedoids(np.array(D_df), initial_medoids, data_type='distance_matrix')
    kmedoids_instance.process()

    with open(f"{CLUSTERS_DIRECTORY}/{METHOD}/{MATRIX_NAME}_{n_clusters}.txt", 'w') as f:
        for cluster in kmedoids_instance.get_clusters():
            f.write(' '.join(map(int2gene.get,cluster)) + '\n')
    t2 = time.time()
    print(f'{n_clusters}: {t2-t1:.2f}sec', end='\r')
# Individual
MATRIX_NAME = "sc_BioGRID_0hell0"
D_df = pd.read_csv(f"{MATRIX_DIRECTORY}/{MATRIX_NAME}.txt", delimiter=' ')

D_df.index = D_df.columns
d = ~D_df.isna().all()
D_df = D_df[d].T[d]

int2gene = dict(enumerate(D_df.columns))

t1 = time.time()
for n_clusters in range(2,100):
    initial_medoids = range(n_clusters)
    kmedoids_instance = kmedoids(np.array(D_df), initial_medoids, data_type='distance_matrix')
    kmedoids_instance.process()

    with open(f"{CLUSTERS_DIRECTORY}/{METHOD}/{MATRIX_NAME}_{n_clusters}.txt", 'w') as f:
        for cluster in kmedoids_instance.get_clusters():
            f.write(' '.join(map(int2gene.get,cluster)) + '\n')
    t2 = time.time()
    print(f'{n_clusters}: {t2-t1:.2f}sec', end='\r')
# Individual
MATRIX_NAME = "sc_BioGRID_1hell1"
D_df = pd.read_csv(f"{MATRIX_DIRECTORY}/{MATRIX_NAME}.txt", delimiter=' ')

D_df.index = D_df.columns
d = ~D_df.isna().all()
D_df = D_df[d].T[d]

int2gene = dict(enumerate(D_df.columns))

t1 = time.time()
for n_clusters in range(2,100):
    initial_medoids = range(n_clusters)
    kmedoids_instance = kmedoids(np.array(D_df), initial_medoids, data_type='distance_matrix')
    kmedoids_instance.process()

    with open(f"{CLUSTERS_DIRECTORY}/{METHOD}/{MATRIX_NAME}_{n_clusters}.txt", 'w') as f:
        for cluster in kmedoids_instance.get_clusters():
            f.write(' '.join(map(int2gene.get,cluster)) + '\n')
    t2 = time.time()
    print(f'{n_clusters}: {t2-t1:.2f}sec', end='\r')
# Individual
MATRIX_NAME = "sc_BioGRID_3hell3"
D_df = pd.read_csv(f"{MATRIX_DIRECTORY}/{MATRIX_NAME}.txt", delimiter=' ')

D_df.index = D_df.columns
d = ~D_df.isna().all()
D_df = D_df[d].T[d]

int2gene = dict(enumerate(D_df.columns))

t1 = time.time()
for n_clusters in range(2,100):
    initial_medoids = range(n_clusters)
    kmedoids_instance = kmedoids(np.array(D_df), initial_medoids, data_type='distance_matrix')
    kmedoids_instance.process()

    with open(f"{CLUSTERS_DIRECTORY}/{METHOD}/{MATRIX_NAME}_{n_clusters}.txt", 'w') as f:
        for cluster in kmedoids_instance.get_clusters():
            f.write(' '.join(map(int2gene.get,cluster)) + '\n')
    t2 = time.time()
    print(f'{n_clusters}: {t2-t1:.2f}sec', end='\r')
# Individual
MATRIX_NAME = "sc_BioGRID_1hell2"
D_df = pd.read_csv(f"{MATRIX_DIRECTORY}/{MATRIX_NAME}.txt", delimiter=' ')

D_df.index = D_df.columns
d = ~D_df.isna().all()
D_df = D_df[d].T[d]

int2gene = dict(enumerate(D_df.columns))

t1 = time.time()
for n_clusters in range(2,100):
    initial_medoids = range(n_clusters)
    kmedoids_instance = kmedoids(np.array(D_df), initial_medoids, data_type='distance_matrix')
    kmedoids_instance.process()

    with open(f"{CLUSTERS_DIRECTORY}/{METHOD}/{MATRIX_NAME}_{n_clusters}.txt", 'w') as f:
        for cluster in kmedoids_instance.get_clusters():
            f.write(' '.join(map(int2gene.get,cluster)) + '\n')
    t2 = time.time()
    print(f'{n_clusters}: {t2-t1:.2f}sec', end='\r')
# Individual
MATRIX_NAME = "sc_BioGRID_2hell1"
D_df = pd.read_csv(f"{MATRIX_DIRECTORY}/{MATRIX_NAME}.txt", delimiter=' ')

D_df.index = D_df.columns
d = ~D_df.isna().all()
D_df = D_df[d].T[d]

int2gene = dict(enumerate(D_df.columns))

t1 = time.time()
for n_clusters in range(2,100):
    initial_medoids = range(n_clusters)
    kmedoids_instance = kmedoids(np.array(D_df), initial_medoids, data_type='distance_matrix')
    kmedoids_instance.process()

    with open(f"{CLUSTERS_DIRECTORY}/{METHOD}/{MATRIX_NAME}_{n_clusters}.txt", 'w') as f:
        for cluster in kmedoids_instance.get_clusters():
            f.write(' '.join(map(int2gene.get,cluster)) + '\n')
    t2 = time.time()
    print(f'{n_clusters}: {t2-t1:.2f}sec', end='\r')<jupyter_output>99: 487.40sec<jupyter_text>Enrichement  Annotation<jupyter_code>distance = "all_hell"

METHOD = "kmedoids"
MATRIX_NAME = f"sc_BioGRID_{distance}"

with open(f"{CLUSTERS_DIRECTORY}/{METHOD}/{MATRIX_NAME}_2.txt", 'r') as f:
            clusters = list(map(str.split, f))
population = {gene for cluster in clusters for gene in cluster}<jupyter_output><empty_output><jupyter_text>Preparation<jupyter_code># Load PPI
PPI = nx.read_edgelist(f"{PPI_DIRECTORY}/BioGRID_sc.txt")
PPI_induced = PPI.subgraph(population)

# Load obo file
GO_FILENAME = "go-basic.obo"
GO_FILEPATH = f"{RAW_DATA_DIRECTORY}/{GO_FILENAME}"
go_dag = obo_parser.GODag(GO_FILEPATH)

all_BP_annotations_df = pd.read_csv(f"{ANNOTATIONS_DIRECTORY}/SGD_BP_sc.csv")
all_MF_annotations_df = pd.read_csv(f"{ANNOTATIONS_DIRECTORY}/SGD_MF_sc.csv")
all_CC_annotations_df = pd.read_csv(f"{ANNOTATIONS_DIRECTORY}/SGD_CC_sc.csv")
PPI_BP_annotations_df = all_BP_annotations_df[all_BP_annotations_df.Systematic_ID.isin(population)]
PPI_MF_annotations_df = all_MF_annotations_df[all_MF_annotations_df.Systematic_ID.isin(population)]
PPI_CC_annotations_df = all_CC_annotations_df[all_CC_annotations_df.Systematic_ID.isin(population)]<jupyter_output><empty_output><jupyter_text>Definition<jupyter_code># Chose namespace
annotation_df = PPI_BP_annotations_df

GO_population = {go_id for go_id in set(annotation_df.GO_ID) 
                           if (5 <= len(annotation_df[annotation_df.GO_ID == go_id]) <= 500 and
                               go_dag[go_id].level > -1)}

annotation_df = annotation_df[annotation_df.GO_ID.isin(GO_population)]
# Conversion dictionaries
int2GO = dict(enumerate(GO_population))
GO2int = dict(zip(int2GO.values(), int2GO.keys()))

GO2genes = {go_id:set(annotation_df.Systematic_ID[annotation_df.GO_ID == go_id])
                for go_id in GO_population}

gene2GOs = {gene :set(annotation_df.GO_ID[annotation_df.Systematic_ID == gene]) 
                for gene in PPI}<jupyter_output><empty_output><jupyter_text>Here we GO<jupyter_code>def gene_enriched_in_cluster(gene, cluster, enrichment):
    return bool(gene2GOs[gene] & set(GO_index[enrichment[cluster]]))

def get_enrichment_df(alpha, p_values):
    m = p_values.size
    c = np.log(m) + np.euler_gamma + 1/(2*m)
    sorted_p_values = np.sort(p_values.values.flatten())
    for k,P_k in enumerate(sorted_p_values,1):
        if P_k > k/(m*c) * alpha:
            break
    threshold = sorted_p_values[k-2]
    return p_values_df < threshold
# List of success states
list_of_success_states = list(GO2genes.values())

# This will be our K, see below. Reshped to fit the shape of k 'array_of_observed_successes'
array_of_total_successes = np.array(list(map(len,list_of_success_states))).reshape(-1,1)
MIN_CLUSTERS = 2
MAX_CLUSTERS = 100
confidences = {'hc':0.01, 
               'mc':0.05, 
               'lc':0.1}

cluster_coverages = defaultdict(dict)
GO_coverages = defaultdict(dict)
gene_coverages = defaultdict(dict)
for distance in [
                 'GCV_malahanobis_mean'
                ]:
    print(distance)
    MATRIX_NAME = f"sc_BioGRID_{distance}"
    
    
    for confidence in confidences:
        cluster_coverages[distance][confidence] = np.zeros(MAX_CLUSTERS-MIN_CLUSTERS)
        GO_coverages[distance][confidence]      = np.zeros(MAX_CLUSTERS-MIN_CLUSTERS)
        gene_coverages[distance][confidence]    = np.zeros(MAX_CLUSTERS-MIN_CLUSTERS)

    t1 = time.time()
    for i, n_clusters in enumerate(range(2, MAX_CLUSTERS)):

        with open(f"{CLUSTERS_DIRECTORY}/{METHOD}/{MATRIX_NAME}_{n_clusters}.txt", 'r') as f:
                list_of_experiments = [set(line.split()) for line in f] 

        clusters = dict(enumerate(list_of_experiments))

        # For each GO term and cluster we get an experiment 
        array_of_observed_successes = np.array([[len(draws & success_states) for draws in list_of_experiments]
                                            for success_states in list_of_success_states])



        K = array_of_total_successes             # defined in section 'Preparation'
        n = list(map(len, list_of_experiments))  # cluster lengths
        k = array_of_observed_successes          # number of annotated genes found in cluster
        N = sum(n)                               # PPI size, i.e. number of all genes that appear in a cluster

        # scipy has a really messed up nomeclature... 
        p_values_array = 1-hypergeom.cdf(k=k-1, M=N, N=n, n=K)
        p_values_df    = pd.DataFrame(p_values_array, index=GO_population)
        GO_index = p_values_df.index

        m = p_values_array.size
        enrichment_dfs = {confidence:get_enrichment_df(alpha,p_values_df) 
                                  for confidence, alpha in confidences.items()}

        for confidence in confidences:
            cluster_coverages[distance][confidence][i] = sum(enrichment_dfs[confidence].any())      /n_clusters
            GO_coverages[distance][confidence][i] = sum(enrichment_dfs[confidence].any(axis=1))/len(GO_population)
            gene_coverages[distance][confidence][i] = sum(1 for cluster in clusters for gene in clusters[cluster] 
                                            if gene_enriched_in_cluster(gene, cluster, enrichment_dfs[confidence]))/N

        t2 = time.time()
        print(f'{n_clusters}: {t2-t1:.2f}sec', end='\r')
    print()<jupyter_output>GCV_malahanobis_mean
99: 160.66sec<jupyter_text>Plots<jupyter_code>plot_distances = ["malahanobis0", "all_hell", "GDV_mahalanobis"]

name2string = {"tvd0": "TVD_{0}",
               "tvd1": "TVD_{1}",
               "tvd2": "TVD_{2}",
               "tvd3": "TVD_{3}",
               "2tvd0": "2TVD_{0}",
               "tvd0123": "TVD_{0123}",
               "all_tvd": "all \;TVDs",
               
               "hell0": "hellinger_{0}",
               "hell1": "hellinger_{1}",
               "hell2": "hellinger_{2}",
               "hell3": "hellinger_{3}",
               "2hell0": "2hellinger_{0}",
               "hell0123": "hellinger_{0123}",
               "all_hell": "all \;hellingers",
              
               "tijana": "Tijana",
               "GDV_mahalanobis":"GDV_{Mahalanobis}",
               "GCV_malahanobis_mean":"GCV_{Mahalanobis}-mean",
               "malahanobis0":"GCV_{Mahalanobis}-0",
               "malahanobisNa":"GCV_{Mahalanobis}"}
#Cluster coverage
figname = 'all'
fig, ax = plt.subplots(figsize=(12,9))
fig.patch.set_alpha(0)
fig.subplots_adjust(hspace = 0.4)

Blues = iter(sns.color_palette("Blues",6)[::-1])
Reds  = iter(sns.color_palette("Reds", 6)[::-1])

for distance in plot_distances:
    ax.plot(range(2,MAX_CLUSTERS), 100*cluster_coverages[distance]['mc'], 
            label=f'${name2string[distance]}$', 
            linewidth=2.5,
            alpha=0.75);
    
    ax.fill_between(range(2,MAX_CLUSTERS), 
                     100*cluster_coverages[distance]['hc'], 
                     100*cluster_coverages[distance]['lc'], 
                     alpha=0.1,);
    ax.set_title('Clusters enriched', fontsize=28)
     
    ax.patch.set_alpha(0)

    ax.set_xlabel('# clusters', fontsize=24)
    ax.set_ylabel('% enriched', fontsize=24)
    ax.tick_params(axis='both', which='major', labelsize=24)

    ax.spines['left'].set_linewidth(2.5)
    ax.spines['left'].set_color('black')

    ax.spines['bottom'].set_linewidth(2.5)
    ax.spines['bottom'].set_color('black')
    
    ax.legend(fontsize=18, shadow=True, facecolor=[0.95, 0.95, 0.95, 0]);
fig.savefig(f"{DATA_DIRECTORY}/plots/{figname}_cluster.png")
#Cluster coverage

fig, ax = plt.subplots(figsize=(12,9))
fig.patch.set_alpha(0)
fig.subplots_adjust(hspace = 0.4)

for distance in plot_distances:
    
    ax.plot(range(2,MAX_CLUSTERS), 100*GO_coverages[distance]['mc'], 
             label=f'${name2string[distance]}$', 
             linewidth=2.5);
    ax.fill_between(range(2,MAX_CLUSTERS), 
                     100*GO_coverages[distance]['hc'], 
                     100*GO_coverages[distance]['lc'], 
                     alpha=0.1);
    ax.set_title('GO-terms enriched', fontsize=28)
    
     
    ax.patch.set_alpha(0)

    ax.set_xlabel('# clusters', fontsize=24)
    ax.set_ylabel('% enriched', fontsize=24)
    ax.tick_params(axis='both', which='major', labelsize=24)

    ax.spines['left'].set_linewidth(2.5)
    ax.spines['left'].set_color('black')

    ax.spines['bottom'].set_linewidth(2.5)
    ax.spines['bottom'].set_color('black')
    
    ax.legend(fontsize=18, shadow=True, facecolor=[0.95, 0.95, 0.95, 0]);
fig.savefig(f"{DATA_DIRECTORY}/plots/{figname}_go.png")
#Cluster coverage

fig, ax = plt.subplots(figsize=(12,9))
fig.patch.set_alpha(0)
fig.subplots_adjust(hspace = 0.4)

for distance in plot_distances:
    
    ax.plot(range(2,MAX_CLUSTERS), 100*gene_coverages[distance]['mc'], 
             label=f'${name2string[distance]}$', 
             linewidth=2.5);
    ax.fill_between(range(2,MAX_CLUSTERS), 
                     100*gene_coverages[distance]['hc'], 
                     100*gene_coverages[distance]['lc'], 
                     alpha=0.1);
    ax.set_title('gene-terms enriched', fontsize=28)
    
     
    ax.patch.set_alpha(0)

    ax.set_xlabel('# clusters', fontsize=24)
    ax.set_ylabel('% enriched', fontsize=24)
    ax.tick_params(axis='both', which='major', labelsize=24)

    ax.spines['left'].set_linewidth(2.5)
    ax.spines['left'].set_color('black')

    ax.spines['bottom'].set_linewidth(2.5)
    ax.spines['bottom'].set_color('black')
    
    ax.legend(fontsize=18, shadow=True, facecolor=[0.95, 0.95, 0.95, 0]);
fig.savefig(f"{DATA_DIRECTORY}/plots/{figname}_gene.png")<jupyter_output><empty_output><jupyter_text>Test<jupyter_code>p_values_df
len(population)
with open(f"{CLUSTERS_DIRECTORY}/{METHOD}/sc_BioGRID_all_hell_2.txt", 'r') as f:
    list_of_experiments = [set(line.split()) for line in f] 

clusters = dict(enumerate(list_of_experiments))

# For each GO term and cluster we get an experiment 
array_of_observed_successes = np.array([[len(draws & success_states) for draws in list_of_experiments]
                                    for success_states in list_of_success_states])



K = array_of_total_successes             # defined in section 'Preparation'
n = list(map(len, list_of_experiments))  # cluster lengths
k = array_of_observed_successes          # number of annotated genes found in cluster
N = sum(n)                               # PPI size, i.e. number of all genes that appear in a cluster

# scipy has a really messed up nomeclature... 
p_values_array = 1-hypergeom.cdf(k=k-1, M=N, N=n, n=K)
p_values_df    = pd.DataFrame(p_values_array, index=GO_population)
GO_index = p_values_df.index

m = p_values_array.size
enrichment_dfs = {confidence:get_enrichment_df(alpha,p_values_df) 
                          for confidence, alpha in confidences.items()}
N<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1615184647144277, 'fraction_numerical': 0.027639727807936387, 'mean_word_length': 4.135257165292501, 'pattern_counts': {'":': 19, '<': 63, '<?xml version=': 0, '>': 62, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '491357', 'license': 'Apache-2.0', 'n_tokens_mistral': 11005, 'n_tokens_neox': 10480, 'n_words': 1464, 'path': 'individual PPIs.ipynb', 'repo_name': 'kiri93/graco', 'size': '26158'}