Complex_predictor/communities.py at master · PasqM/Complex_predictor · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
import pandas as pd
from  tqdm import tqdm
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.metrics import roc_curve, auc
import sqlite3
import gzip
import numpy as np
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster
import networkx as nx
import community as community_louvain
import pickle
from sklearn.metrics import precision_recall_curve
from sklearn.metrics import average_precision_score
from sklearn.metrics import confusion_matrix
from cdlib import algorithms, evaluation

class Communities:
    def __init__(self, communities):
        self.communities = communities

with open("proteome_list.txt") as f:
    proteome = f.read().splitlines()
proteome = set(proteome)

def fix(c):
    if 'Q9BXH1' in c:
        c.remove('Q9BXH1')
        c.add('Q96PG8')
    if 'Q13765' in c:
        c.remove('Q13765')
        c.add('E9PAV3')
    if 'P58400' in c:
        c.remove('P58400')
        c.add('Q9ULB1')
    if 'O96033' in c:
        c.remove('O96033')
        c.add('O96007')
    if 'P0CG05' in c:
        c.remove('P0CG05')
        c.add('P0DOY2')
    return c

complexes = pd.read_csv("/projects/bioinformatics/DB/ComplexPortal/09_2024/9606.tsv", sep="\t")
complexes['Identifiers (and stoichiometry) of molecules in complex'] = complexes['Identifiers (and stoichiometry) of molecules in complex'].str.replace(',', '|')
complexes['Expanded participant list'] = complexes['Expanded participant list'].str.replace(',', '|')
complexes['complex_ids'] = complexes['Identifiers (and stoichiometry) of molecules in complex'].str.split('|')
complexes['extended_complex_ids'] = complexes['Expanded participant list'].str.split('|')
complexes['complex_ids'] = complexes['complex_ids'].apply(lambda x: set([i.split('(')[0] for i in x]))
complexes['extended_complex_ids'] = complexes['extended_complex_ids'].apply(lambda x: set([i.split('(')[0] for i in x]))
complexes['complex_ids'] = complexes['complex_ids'].apply(lambda x: set([i.split('-')[0] for i in x]))
complexes['extended_complex_ids'] = complexes['extended_complex_ids'].apply(lambda x: set([i.split('-')[0] for i in x]))
complexes['complex_ids'] = complexes['complex_ids'].apply(lambda x: set([i.strip('[').strip(']') for i in x]))
complexes['extended_complex_ids'] = complexes['extended_complex_ids'].apply(lambda x: set([i.strip('[').strip(']') for i in x]))
complexes['extended_complex_ids'] = complexes['extended_complex_ids'].apply(lambda x: fix(x).intersection(proteome))
#complexes = complexes[complexes['extended_complex_ids'].apply(lambda x: len(x) > 1)]
complexes['extended_complex_ids'] = complexes['extended_complex_ids'].apply(lambda x: tuple(sorted(list(x))))

# Merge complexes with the same extended_complex_ids

complexes = complexes.drop_duplicates(subset='extended_complex_ids', keep='first')
complexes['extended_complex_ids'] = complexes['extended_complex_ids'].apply(lambda x: set(x))

all_complexes = set()
for complex in complexes['extended_complex_ids']:
    all_complexes = all_complexes.union(complex)

df = pd.read_feather('predictions_new.feather')
G = nx.from_pandas_edgelist(df, 'prot1', 'prot2', create_using=nx.Graph(), edge_attr='prediction')
print("Loaded predictions...")

true_comms = []
pairs = set()
seen = set()
unseen = set()
for cmpx in complexes['extended_complex_ids']:
    c = list(cmpx)
    true_comms.append([])
    for i, prot1 in enumerate(c):
        if prot1 not in proteome:
            unseen.add(prot1)
            continue
        true_comms[-1].append(prot1)
        for prot2 in c[i+1:]:
            if prot2 not in proteome:
                unseen.add(prot2)
                continue
            pairs.add((min(prot1, prot2), max(prot1, prot2)))
            seen.add(prot1)
            seen.add(prot2)

for i in range(1,100):
    coms = algorithms.louvain(G, weight='prediction', resolution=i)#.ipca(G, weights='label', t_in=i/100)
    print(i, coms.communities, flush=True)
    true_comms = []
    valid_proteins = set(G.nodes())
    for j in all_complexes.difference(set([j for c in coms.communities for j in c])):
        coms.communities.append([j])
    for c in complexes['extended_complex_ids']:
        true_comms.append([])
        for prot in c.intersection(valid_proteins):
            true_comms[-1].append(prot)
    true_comms = Communities(true_comms)
    for j in range(len(coms.communities)):
        coms.communities[j] = set(coms.communities[j])
        coms.communities[j] = coms.communities[j].intersection(all_complexes)
        coms.communities[j] = list(coms.communities[j])
    print(i, len(coms.communities), len(true_comms.communities), evaluation.geometric_accuracy(coms,true_comms), flush=True)