evaluate.py 18.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
#!/usr/bin/env python3


import sys
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
13
                        help='The file to evalute')
14
    parser.add_argument('--type', '-t', choices=["d2", "path", "d2-2annotate", "dgraphs"], default="path", required=True,
15
                        help="Define the data type to evaluate. Must be 'd2' or 'path' or 'd2-2annotate' (Rayan's hack).")
Yoann Dufresne's avatar
Yoann Dufresne committed
16
17
    parser.add_argument('--light-print', '-l', action='store_true',
                        help='Print only wrong nodes and paths')
Yoann Dufresne's avatar
Yoann Dufresne committed
18
    parser.add_argument('--barcode_graph', '-b', help="Path to the barcode graph corresponding to the d2_graph to analyse.")
19
    parser.add_argument('--optimization_file', '-o',
Yoann Dufresne's avatar
Yoann Dufresne committed
20
                        help="If the main file is a d2, a file formatted for optimization can be set. This file will be used to compute the coverage of the longest path on the barcode graph.")
21
22
23
24
25
26
27
28
29
30
31
32
33
34

    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()

35
36
37
38
39
40
41
42
43
44
45
def save_graph(g, filename):
    if filename.endswith('.graphml'):
        return nx.write_graphml(g,filename)
    elif filename.endswith('.gexf'):
        return nx.write_gexf(g,filename)
    else:
        print("Wrong file format. Require graphml or gefx format", file=sys.stderr)
        exit()



46

Yoann Dufresne's avatar
Yoann Dufresne committed
47
48
# ------------- Path Graph -------------

Yoann Dufresne's avatar
Yoann Dufresne committed
49
50
51
52
def mols_from_node(node_name):
    return [int(idx) for idx in node_name.split(":")[1].split(".")[0].split("_")]


53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
def parse_udg_qualities(graph):
    """ Compute the quality for the best udgs present in the graph.
        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
        :return: A tuple containing two dictionaries. The first one with theoritical frequencies of each node, the second one with observed frequencies.
    """
    dg_per_node = {}

    for node, data in graph.nodes(data=True):
        str_udg = data["udg"]
        central, h1, h2 = str_to_udg_lists(str_udg)

        if central not in dg_per_node:
            dg_per_node[central] = []
        dg_per_node[central].append(data["udg"])

    for node in dg_per_node:
        print(node, dg_per_node[node])
        print(len(dg_per_node))

    return dg_per_node


Yoann Dufresne's avatar
Yoann Dufresne committed
77
def parse_path_graph_frequencies(graph, barcode_graph):
78
79
80
81
82
83
84
85
    """ 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 representing 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.
        :return: A tuple containing two dictionaries. The first one with theoretical frequencies of each node, the second one with observed frequencies.
    """
Yoann Dufresne's avatar
Yoann Dufresne committed
86
    # Compute origin nodes formatted as `{idx}:{mol1_id}_{mol2_id}_...`
Yoann Dufresne's avatar
Yoann Dufresne committed
87
    observed_frequencies = {}
Yoann Dufresne's avatar
Yoann Dufresne committed
88
    real_frequencies = {}
89
    origin_node_names = []
Yoann Dufresne's avatar
Yoann Dufresne committed
90
91
    node_per_barcode = {}

Yoann Dufresne's avatar
Yoann Dufresne committed
92
    for node, data in graph.nodes(data=True):
93
94
        parsed = parse_dg_name(graph, node)
        origin_name = parsed[0][1]
Yoann Dufresne's avatar
Yoann Dufresne committed
95
96

        if origin_name not in node_per_barcode:
Yoann Dufresne's avatar
Yoann Dufresne committed
97
98
            node_per_barcode[origin_name] = []
        node_per_barcode[origin_name].append(node)
99
100

        # Count frequency
Yoann Dufresne's avatar
Yoann Dufresne committed
101
        if origin_name not in observed_frequencies:
Yoann Dufresne's avatar
Yoann Dufresne committed
102
            observed_frequencies[origin_name] = 0
103
            origin_node_names.append(origin_name)
Yoann Dufresne's avatar
Yoann Dufresne committed
104
        observed_frequencies[origin_name] += 1
105

Yoann Dufresne's avatar
Yoann Dufresne committed
106
107
    # Theoretical frequencies
    real_frequencies = {node_id: len(node_id.split(":")[1].split("_")) for node_id in barcode_graph.nodes()}
108

Yoann Dufresne's avatar
Yoann Dufresne committed
109
    return real_frequencies, observed_frequencies, node_per_barcode
110
111


Yoann Dufresne's avatar
Yoann Dufresne committed
112
113
114
115
116
117
118
119
120
""" 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))
121

Yoann Dufresne's avatar
Yoann Dufresne committed
122
123
124
125
126
127
128
129
        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)
130

Yoann Dufresne's avatar
Yoann Dufresne committed
131
    return neighborhood
132
133


Yoann Dufresne's avatar
Yoann Dufresne committed
134
135
def print_path_summary(frequencies, light_print=False, file_pointer=sys.stdout):
    if file_pointer is None:
Yoann Dufresne's avatar
Yoann Dufresne committed
136
137
        return

Yoann Dufresne's avatar
Yoann Dufresne committed
138
    print("--- Nodes analysis ---", file=file_pointer)
Yoann Dufresne's avatar
Yoann Dufresne committed
139
140
141
    theoretical_frequencies, observed_frequencies, node_per_barcode = frequencies
    for key in theoretical_frequencies:
        obs, the = observed_frequencies[key] if key in observed_frequencies else 0, theoretical_frequencies[key]
Yoann Dufresne's avatar
Yoann Dufresne committed
142
143
144
145
146
        result = f"{key}: {obs}/{the}"

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

Yoann Dufresne's avatar
Yoann Dufresne committed
147
        if light_print and obs == the:
Yoann Dufresne's avatar
Yoann Dufresne committed
148
149
150
151
152
            continue

        print(result, file=file_pointer)


153
154
    print("--- Global summary ---", file=file_pointer)

Yoann Dufresne's avatar
Yoann Dufresne committed
155
156
    # --- Frequency usage ---
    # Tags
Yoann Dufresne's avatar
Yoann Dufresne committed
157
158
159
    distinct_theoretical_nodes = len(theoretical_frequencies)
    distinct_observed_nodes = len(observed_frequencies)
    print(f"Distinct barcodes: {distinct_observed_nodes}/{distinct_theoretical_nodes}", file=file_pointer)
Yoann Dufresne's avatar
Yoann Dufresne committed
160
    # molecules
Yoann Dufresne's avatar
Yoann Dufresne committed
161
162
163
    cumulative_theoretical_nodes = sum(theoretical_frequencies.values())
    cumulative_observed_nodes = sum(observed_frequencies.values())
    print(f"Molecules: {cumulative_observed_nodes}/{cumulative_theoretical_nodes}", file=file_pointer)
Yoann Dufresne's avatar
Yoann Dufresne committed
164
165
166
    # Wrong splits
    over_split = 0
    under_split = 0
Yoann Dufresne's avatar
Yoann Dufresne committed
167
168
169
170
171
    for barcode in theoretical_frequencies:
        observed = observed_frequencies[barcode] if barcode in observed_frequencies else 0
        theoretic = theoretical_frequencies[barcode]
        over_split += max(observed-theoretic, 0)
        under_split += max(theoretic-observed, 0)
Yoann Dufresne's avatar
Yoann Dufresne committed
172
    print(f"Under/Over splitting: {under_split} - {over_split}")
173
174


175
176
177
178
def print_dgraphs_summary(frequencies, light_print=False):
    pass


Yoann Dufresne's avatar
Yoann Dufresne committed
179
180
# ------------- D2 Graph -------------

181
182
183
184
def str_to_udg_lists(s):
    udg = s.replace("]", "").replace(' [', '[')
    return udg.split('[')

185
186
def parse_dg_name(gr, name):
    udg = nx.get_node_attributes(gr, 'udg')[name]
187
    central, h1, h2 = str_to_udg_lists(udg)
Yoann Dufresne's avatar
Yoann Dufresne committed
188
    
189
190
    idx = name
    score = nx.get_node_attributes(gr, 'score')[name]
Yoann Dufresne's avatar
Yoann Dufresne committed
191
192
193
194
195
196
197
198

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

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


199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
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
220
221
222
    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])}")
Yoann Dufresne's avatar
Yoann Dufresne committed
223
224
    no_singleton = [x for x in connected_components if len(x) > 1]
    print(f"There are {len(no_singleton)} connex components with at least 2 nodes")
Yoann Dufresne's avatar
Yoann Dufresne committed
225
226
227
    print(f"The 5 largest components: {[len(x) for x in connected_components][:5]}")

    print("--- Largest component analysis ---")
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
    # 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
247
248
249
250
251
252
253


def compute_next_nodes(d2_component):
    next_nodes = {}

    for node in d2_component.nodes():
        # Parse the current node name
254
        head, h1, h2 = parse_dg_name(d2_component,node)
Yoann Dufresne's avatar
Yoann Dufresne committed
255
256
257
258
259
260
261
262
        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:
263
                nei_head, _, _ = parse_dg_name(d2_component,neighbor)
Yoann Dufresne's avatar
Yoann Dufresne committed
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
                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

308
    path.append((molecule, node))
Yoann Dufresne's avatar
Yoann Dufresne committed
309
310
311
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
    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]


341
342
343
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
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

376
377
378
379
380
381
382
383
384
385
386
387
# returns True iff there exist x in mol1 such that there exists y in mol2 and |x-y| <= some_value
def nearby_udg_molecules(mols1, mols2):
    for x in mols1:
        for y in mols2:
            if abs(x-y) <= 5:
                return True
    return False

def verify_graph_edges(d2_component):
    udg_molecules_dict=dict()
    for node in d2_component.nodes():
        # Parse the current node name
388
        head, c1, c2 = parse_dg_name(d2_component,node)
389
390
391

        # Construct the molecule(s) that this udg really 'reflects'
        # i.e. the udg has a central node and two cliques
rchikhi's avatar
rchikhi committed
392
393
        # that central node is the result of merging of several molecules
        # ideally, only one of those molecules is connected to the molecules of the cliques
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
        # (there could be more than one though; in that case the udg is 'ambiguous')
        # udg_molecules aims to reflect the molecule(s) underlying this udg 
        udg_molecules = set()

        # Get the current molecule idxs of central node
        molecule_idxs = mols_from_node(head[1])
        #print("mol idxs", molecule_idxs)

        # Examine molecule idx's of cliques to see which are close to the central node
        # rationale: c1/c2 contain nearby molecule id's
        for mol_idx in molecule_idxs:
            nexts = []
            for c in [c1,c2]:
                for c_node in c:
                    nei_mols = mols_from_node(c_node.split()[0])
                    nei_mols = [x for x in nei_mols if x > mol_idx] # fixme: also look at the <= molecules for more robustness
                    
                    # If there are molecule next
                    if len(nei_mols) > 0:
                        next_nei_mol = min(nei_mols)
                        nexts.append((next_nei_mol))

            nexts.sort(key=lambda x: x)
            quality = sum([1.0/x if mol_idx+x in nexts else 0 for x in range(1,6)]) / sum([1.0/x for x in range(1,6)])
rchikhi's avatar
rchikhi committed
418
            if quality > 0.6: # eyeballed but still arbitrary
419
420
421
422
423
424
425
426
                udg_molecules.add(mol_idx)
            #print("mol",mol_idx,molecule_idxs,"quality",quality,"nexts",nexts)
       
        udg_molecules_dict[head[0]]=udg_molecules

    # Then: annotate edges as to whether they're real (their udg_molecule(s) are nearby) or fake
    for n1, n2 , data in d2_component.edges(data=True):
        # Parse the current node name
427
        head, c1, c2 = parse_dg_name(d2_component,n1)
428
429
        node_udg_molecules = udg_molecules_dict[head[0]]
        
430
        n_head, n_c1, n_c2 = parse_dg_name(d2_component,n2)
431
432
433
434
435
436
437
        neighbor_udg_molecules = udg_molecules_dict[n_head[0]]
       
        if nearby_udg_molecules(node_udg_molecules, neighbor_udg_molecules):
            color = 'green'
        else:
            color = 'red'
        data['color'] = color
438
        #print("edge",node_udg_molecules,neighbor_udg_molecules,color)
439
440
441
442
    
    # also, annotate nodes by their putative molecule found
    for n, data in d2_component.nodes(data=True):
        # Parse the current node name
443
        head, c1, c2 = parse_dg_name(d2_component,n)
444
445
446
        node_udg_molecules = udg_molecules_dict[head[0]]
        data['udg_molecule']= '_'.join(list(map(str,node_udg_molecules)))
        
447
448
449
450
451
452
453
454
    # aggressive: delete nodes which have either no found udg_molecule, or two udg_molecules
    # turns out it's not a good strategy as the nodes with two udg_molecules are important to connect portions of graph
    # but what if we magically keep those where the two adjacent molecules are close together
    if True:
        d2_component = d2_component.copy()
        nodes_to_remove = []
        for n, data in d2_component.nodes(data=True):
            # Parse the current node name
455
            head, c1, c2 = parse_dg_name(d2_component,n)
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
            if "_" in data['udg_molecule'] or data['udg_molecule'] == '':
                if "_" in data['udg_molecule']:
                    m1, m2 = list(map(int,data['udg_molecule'].split("_")))
                    if abs(m2-m1) < 30: continue # don't remove that kind of nodes
                nodes_to_remove += [n]
        d2_component.remove_nodes_from(nodes_to_remove)
        print("removed",len(nodes_to_remove),"bad nodes")

    # aggressive: delete red edges
    if True:
        d2_component = d2_component.copy()
        edges_to_remove = []
        for n1, n2, data in d2_component.edges(data=True):
            if data['color'] == 'red':
                edges_to_remove += [(n1,n2)]
        d2_component.remove_edges_from(edges_to_remove)
        print("removed",len(edges_to_remove),"bad edges")


    return d2_component
476

477
478
479
480
def main():
    args = parse_args()
    graph = load_graph(args.filename)

Yoann Dufresne's avatar
Yoann Dufresne committed
481
    if args.type == "path":
Yoann Dufresne's avatar
Yoann Dufresne committed
482
483
        barcode_graph = load_graph(args.barcode_graph)
        frequencies = parse_path_graph_frequencies(graph, barcode_graph)
Yoann Dufresne's avatar
Yoann Dufresne committed
484

Yoann Dufresne's avatar
Yoann Dufresne committed
485
        print_path_summary(frequencies, light_print=args.light_print)
486
487
488
    elif args.type == "dgraphs":
        udg_per_node = parse_udg_qualities(graph)
        # print(udg_per_node)
Yoann Dufresne's avatar
Yoann Dufresne committed
489
490
491
492
493
494
    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)
495
496
497
498
        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)
499
500
501
502
503
504
505
506
507
    
    # added by Rayan
    # example:
    # python evaluate.py --type d2-2annotate ~/Dropbox/cnrs/projects/10x-barcodes/yoann_to_cedric_ilp/d2_graph.gexf
    elif args.type == "d2-2annotate":
        components = list(nx.connected_components(graph))
        components.sort(key=lambda x: -len(x))
        component = graph.subgraph(components[0])
 
508
        component = verify_graph_edges(component)
509

Yoann Dufresne's avatar
Yoann Dufresne committed
510
511
512
513
        extension = args.filename.split('.')[-1]
        base_filename = '.'.join(args.filename.split('.')[:-1])
        save_graph(component, base_filename+".verified."+extension)

514
515
516

if __name__ == "__main__":
    main()