evaluate.py 13.4 KB
Newer Older
1
2
3
4
#!/usr/bin/env python3


import sys
5
import csv
6
7
8
9
10
11
12
13
14
import argparse
from termcolor import colored

import networkx as nx


def parse_args():
    parser = argparse.ArgumentParser(description='Process some integers.')
    parser.add_argument('filename', type=str,
Yoann Dufresne's avatar
Yoann Dufresne committed
15
16
17
                        help='The file to evalute')
    parser.add_argument('--type', '-t', choices=["d2", "path"], default="path", required=True,
                        help="Define the data type to evaluate. Must be 'd2' or 'path'.")
Yoann Dufresne's avatar
Yoann Dufresne committed
18
19
    parser.add_argument('--light-print', '-l', action='store_true',
                        help='Print only wrong nodes and paths')
20
21
    parser.add_argument('--optimization_file', '-o',
                        help="If the main file is a d2, a file formated for optimization can be set. This file will be used to compute the coverage of the longest path on the barcode graph.")
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

    args = parser.parse_args()
    return args


def load_graph(filename):
    if filename.endswith('.graphml'):
        return nx.read_graphml(filename)
    elif filename.endswith('.gexf'):
        return nx.read_gexf(filename)
    else:
        print("Wrong file format. Require graphml or gefx format", file=sys.stderr)
        exit()


Yoann Dufresne's avatar
Yoann Dufresne committed
37
38
# ------------- Path Graph -------------

Yoann Dufresne's avatar
Yoann Dufresne committed
39
40
41
42
def mols_from_node(node_name):
    return [int(idx) for idx in node_name.split(":")[1].split(".")[0].split("_")]


43
44
45
46
47
48
""" Compute appearance frequencies from node names.
    All the node names must be under the format :
    {idx}:{mol1_id}_{mol2_id}_...{molx_id}.other_things_here
    @param graph The networkx graph representinf the deconvolved graph
    @param only_wong If True, don't print correct nodes
    @param file_pointer Where to print the output. If set to stdout, then pretty print. If set to None, don't print anything.
Yoann Dufresne's avatar
Yoann Dufresne committed
49
    @return A tuple containing two dictionaries. The first one with theoritical frequencies of each node, the second one with observed frequencies.
50
"""
Yoann Dufresne's avatar
Yoann Dufresne committed
51
def parse_path_graph_frequencies(graph):
52
    # Compute origin nodes formated as `{idx}:{mol1_id}_{mol2_id}_...`
Yoann Dufresne's avatar
Yoann Dufresne committed
53
    observed_frequencies = {}
54
    origin_node_names = []
Yoann Dufresne's avatar
Yoann Dufresne committed
55
    node_per_barcode = {}
56
    for node in graph.nodes():
Yoann Dufresne's avatar
Yoann Dufresne committed
57
58
59
60
61
        origin_name = node.split(".")[0]

        if not origin_name in node_per_barcode:
            node_per_barcode[origin_name] = []
        node_per_barcode[origin_name].append(node)
62
63

        # Count frequency
Yoann Dufresne's avatar
Yoann Dufresne committed
64
65
        if not origin_name in observed_frequencies:
            observed_frequencies[origin_name] = 0
66
            origin_node_names.append(origin_name)
Yoann Dufresne's avatar
Yoann Dufresne committed
67
        observed_frequencies[origin_name] += 1
68
69
70
71
72
73
74
75
76
77
    
    # Compute wanted frequencies
    theoritical_frequencies = {}
    for node_name in origin_node_names:
        _, composition = node_name.split(':')

        mol_ids = composition.split('_')
        # The node should be splited into the number of molecules inside itself
        theoritical_frequencies[node_name] = len(mol_ids)

Yoann Dufresne's avatar
Yoann Dufresne committed
78
    return theoritical_frequencies, observed_frequencies, node_per_barcode
79
80


Yoann Dufresne's avatar
Yoann Dufresne committed
81
82
83
84
85
86
87
88
89
""" This function aims to look for direct molecule neighbors.
    If a node has more than 2 direct neighbors, it's not rightly splitted
"""
def parse_graph_path(graph):
    neighborhood = {}

    for node in graph.nodes():
        molecules = mols_from_node(node)
        neighbors = list(graph.neighbors(node))
90

Yoann Dufresne's avatar
Yoann Dufresne committed
91
92
93
94
95
96
97
98
        neighborhood[node] = []
        for mol in molecules:
            for nei in neighbors:
                nei_mols = mols_from_node(nei)
                if mol-1 in nei_mols:
                    neighborhood[node].append(nei)
                if mol+1 in nei_mols:
                    neighborhood[node].append(nei)
99

Yoann Dufresne's avatar
Yoann Dufresne committed
100
    return neighborhood
101
102


Yoann Dufresne's avatar
Yoann Dufresne committed
103
def print_path_summary(frequencies, neighborhood, light_print=False, file_pointer=sys.stdout):
Yoann Dufresne's avatar
Yoann Dufresne committed
104
105
106
    if file_pointer == None:
        return

Yoann Dufresne's avatar
Yoann Dufresne committed
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
    print("--- Nodes analysis ---", file=file_pointer)
    theoritical_frequencies, observed_frequencies,  node_per_barcode = frequencies
    for key in theoritical_frequencies:
        obs, the = observed_frequencies[key], theoritical_frequencies[key]
        result = f"{key}: {obs}/{the}"

        if file_pointer == sys.stdout:
            result = colored(result, 'green' if obs==the else 'red')

        # Compute neighborhood correctness
        neighborhood_ok = True
        for node in node_per_barcode[key]:
            if len(neighborhood[node]) != 2:
                neighborhood_ok = False

Yoann Dufresne's avatar
Yoann Dufresne committed
122
        if light_print and obs==the and neighborhood_ok:
Yoann Dufresne's avatar
Yoann Dufresne committed
123
124
125
126
127
128
129
130
131
132
133
134
            continue

        print(result, file=file_pointer)
        for node in node_per_barcode[key]:
            text = f"\t{node}\t{' '.join(neighborhood[node])}"

            if file_pointer == sys.stdout:
                text = colored(text, 'green' if len(neighborhood[node]) == 2 else 'yellow')

            print(text, file=file_pointer)


135
136
    print("--- Global summary ---", file=file_pointer)

Yoann Dufresne's avatar
Yoann Dufresne committed
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
    # --- Frequency usage ---
    # Tags
    distinct_theoritical_nodes = len(frequencies[0])
    distinct_observed_nodes = len(frequencies[1])
    print(f"Distinct barcodes: {distinct_observed_nodes}/{distinct_theoritical_nodes}", file=file_pointer)
    # molecules
    cumulative_theoritical_nodes = sum(frequencies[0].values())
    cumulative_observed_nodes = sum(frequencies[1].values())
    print(f"Molecules: {cumulative_observed_nodes}/{cumulative_theoritical_nodes}", file=file_pointer)
    # Wrong splits
    over_split = 0
    under_split = 0
    for barcode in frequencies[0]:
        observed = frequencies[1][barcode]
        theoritic = frequencies[0][barcode]
        over_split += max(observed-theoritic, 0)
        under_split += max(theoritic-observed, 0)
    print(f"Under/Over splitting: {under_split} - {over_split}")
155
156


Yoann Dufresne's avatar
Yoann Dufresne committed
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
# ------------- D2 Graph -------------

def parse_dg_name(name):
    name = name.replace("]", "").replace(' [', '[')
    header, h1, h2 = name.split('[')
    
    # Parse header
    header = header.split(" ")
    idx = central = score = -1
    if len(header) == 3:
        idx, central, score = header
    else:
        central, score = header

    # Parse hands
    h1 = h1.split(', ')
    h2 = h2.split(', ')

    return (idx, central, score), h1, h2


178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
def path_to_jumps(path):
    chuncks = []
    prev_start = -1000
    current_molecule = -1000

    for mol, node in path:
        # If there is a gap
        if mol > current_molecule + 1:
            chuncks.append((prev_start, current_molecule))
            prev_start = mol

        current_molecule = mol

    # Add the last piece
    chuncks.append((prev_start, current_molecule))

    del chuncks[0]
    return chuncks


def print_d2_summary(connected_components, longest_path, covered_vars={}, light_print=False):
Yoann Dufresne's avatar
Yoann Dufresne committed
199
200
201
202
203
204
    print("--- Global summary ---")
    print(f"Number of connected components: {len(connected_components)}")
    print(f"Total number of nodes: {sum([len(x) for x in connected_components])}")
    print(f"The 5 largest components: {[len(x) for x in connected_components][:5]}")

    print("--- Largest component analysis ---")
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
    # Get the list of node idx
    path_dg_idx = [int(x[1].split(" ")[0]) for x in longest_path]
    # print("\n".join(longest_path))
    if not light_print:
        print("Longest path for increasing molecule number:")
        print(path_dg_idx)
    print(f"Size of the longest path: {len(longest_path)}")
    print("Jumps in central nodes:")
    print(path_to_jumps(longest_path))
    print(f"Number of optimization variable coverage: {len(covered_vars)}")
    nb_true = 0
    falses = []
    for idx, val in covered_vars.items():
        if val:
            nb_true += 1
        else:
            falses.append(idx)
    print(f"Coverage: {nb_true}/{len(covered_vars)}")
    print(f"Uncovered_values:\n{falses}")
Yoann Dufresne's avatar
Yoann Dufresne committed
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
    


# def component_to_nearest_neighbor_graph(component):
#     nng = nx.Graph()
#     nng.add_nodes_from(component.nodes())

#     for edge in component.edges():
#         node1, node2 = edge
#         node1 = parse_dg_name(node1)
#         node2 = parse_dg_name(node2)

#         central1 = mols_from_node(node1[0][1])
#         central2 = frozenset(mols_from_node(node2[0][1]))

#         for mol1 in central1:
#             if mol1-1 in central2 or mol1+1 in central2:
#                 nng.add_edge(edge[0], edge[1])

#     componnents = list(nx.connected_components(nng))
#     print([len(x) for x in componnents])
#     componnents.sort(key=lambda x: -len(x))
#     componnents = [nng.subgraph(x) for x in componnents]
#     nx.write_gexf(componnents[0], "data/d2_reducted.gexf")

#     return nng, componnents


def compute_next_nodes(d2_component):
    next_nodes = {}

    for node in d2_component.nodes():
        # Parse the current node name
        head, h1, h2 = parse_dg_name(node)
        next_nodes[node] = {}
        neighbors = list(d2_component.neighbors(node))

        # Get the current molecule idxs
        molecule_idxs = mols_from_node(head[1])
        for mol_idx in molecule_idxs:
            nexts = []
            for neighbor in neighbors:
                nei_head, _, _ = parse_dg_name(neighbor)
                nei_mols = mols_from_node(nei_head[1])
                nei_mols = [x for x in nei_mols if x > mol_idx]
                
                # If there are molecule next
                if len(nei_mols) > 0:
                    next_nei_mol = min(nei_mols)
                    nexts.append((next_nei_mol, neighbor))

            nexts.sort(key=lambda x: x[0])
            next_nodes[node][mol_idx] = nexts
            # print(next_nodes)

    return next_nodes


def compute_longest_increasing_paths(d2_component):
    next_nodes = compute_next_nodes(d2_component)

    # Compute the longest path for each node
    longest_paths = {}
    for idx, start_node in enumerate(next_nodes):
        # print(f"{idx}/{len(next_nodes)}")
        for mol_idx in next_nodes[start_node]:
            recursive_longest_path(start_node, mol_idx , next_nodes, longest_paths)

    # Get the longest path size
    max_len, node_val, mol_idx = 0, None, -1
    for node in longest_paths:
        for mol in longest_paths[node]:
            length, _, _ = longest_paths[node][mol]
            if max_len < length:
                max_len = length
                node_val = node
                mol_idx = mol

    # Backtrack the longest path
    path = backtrack_longest_path(node_val, mol_idx, longest_paths)
    return path


def backtrack_longest_path(node, molecule, longest_paths, path=[]):
    if node == None:
        return path

311
    path.append((molecule, node))
Yoann Dufresne's avatar
Yoann Dufresne committed
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
    length, next_node, next_mol = longest_paths[node][molecule]
    return backtrack_longest_path(next_node, next_mol, longest_paths, path)


def recursive_longest_path(current_node, current_molecule, next_nodes, longest_paths):
    # Dynamic programming
    if current_node in longest_paths and current_molecule in longest_paths[current_node]:
        return longest_paths[current_node][current_molecule]

    longest = 0
    longest_next = None
    min_mol_idx = current_molecule + 10000
    # Recursively compute the longest path
    for mol_idx, node in next_nodes[current_node][current_molecule]:
        size, _, _ = recursive_longest_path(node, mol_idx, next_nodes, longest_paths)
        if size + 1 > longest:
            longest = size + 1
            longest_next = node
            min_mol_idx = mol_idx
        # If there is an alternative path with shorter distance
        elif size + 1 == longest and mol_idx < min_mol_idx:
            longest = size + 1
            longest_next = node
            min_mol_idx = mol_idx
    
    # Save the result
    if not current_node in longest_paths:
        longest_paths[current_node] = {}
    longest_paths[current_node][current_molecule] = (longest, longest_next, min_mol_idx)
    return longest_paths[current_node][current_molecule]


344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
def compute_covered_variables(optimization_file, path):
    vars = None
    var_assignments = {}

    # Read optimization variables
    with open(optimization_file) as of:
        header = of.readline()
        header = [int(x) for x in header.split(" ")]
        nb_nodes = header[0]
        nb_vars = header[1]
        vars = {x:False for x in range(nb_vars)}
        # nb_true = 0
        # for x in vars.values():
        #     if x: nb_true += 1
        # print(nb_true)
        # exit()

        for idx, line in enumerate(of):
            # Stop at the end of nodes
            if idx >= nb_nodes:
                break

            parsed = [int(x) for x in line.split(' ')]
            var_assignments[parsed[0]] = parsed[1:]

    print(var_assignments[0])

    # Read the path to cover the variables
    for node in path:
        node_idx = int(node[1].split(" ")[0])
        for var_idx in var_assignments[node_idx]:
            vars[var_idx] = True

    return vars


380
381
382
383
def main():
    args = parse_args()
    graph = load_graph(args.filename)

Yoann Dufresne's avatar
Yoann Dufresne committed
384
385
386
387
388
389
390
391
392
393
394
    if args.type == "path":
        frequencies = parse_path_graph_frequencies(graph)
        neighborhood = parse_graph_path(graph)

        print_path_summary(frequencies, neighborhood, light_print=args.light_print)
    elif args.type == "d2":
        components = list(nx.connected_components(graph))
        components.sort(key=lambda x: -len(x))

        component = graph.subgraph(components[0])
        longest_path = compute_longest_increasing_paths(component)
395
396
397
398
        covered_vars = {}
        if args.optimization_file and len(args.optimization_file) > 0:
            covered_vars = compute_covered_variables(args.optimization_file, longest_path)
        print_d2_summary(components, longest_path, covered_vars=covered_vars, light_print=args.light_print)
399
400
401
402


if __name__ == "__main__":
    main()