evaluate.py 26.3 KB
Newer Older
1
2
3
4
5
6
7
#!/usr/bin/env python3


import sys
import argparse
from termcolor import colored
import networkx as nx
8
sys.setrecursionlimit(10000)
9
10
11


def parse_args():
12
    parser = argparse.ArgumentParser(description='Process a d2 graph (complete graph or path) to evaluate its quality.')
13
    parser.add_argument('filename', type=str,
Yoann Dufresne's avatar
Yoann Dufresne committed
14
                        help='The file to evalute')
Yoann Dufresne's avatar
Yoann Dufresne committed
15
    parser.add_argument('--type', '-t', choices=["d2", "path", "bg", "d2-2annotate", "dgraphs"], default="path", required=True,
16
                        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
17
18
    parser.add_argument('--light-print', '-l', action='store_true',
                        help='Print only wrong nodes and paths')
19
    parser.add_argument('--max_gap', '-g', type=int, default=0, help="Allow to jump over max_gap nodes during the increasing path search")
Yoann Dufresne's avatar
Yoann Dufresne committed
20
    parser.add_argument('--barcode_graph', '-b', help="Path to the barcode graph corresponding to the d2_graph to analyse.")
21
    parser.add_argument('--optimization_file', '-o',
Yoann Dufresne's avatar
Yoann Dufresne committed
22
                        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.")
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()

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

Yoann Dufresne's avatar
Yoann Dufresne committed
46
47
48
49
50
51
52
53
54
55
56
57
58
59
def transform_bg(graph):
    idx = 0
    node_names = {}
    nx.set_node_attributes(graph, 0, 'score')
    nx.set_node_attributes(graph, "", 'barcode_edges')
    for node in graph.nodes():
        graph.nodes[node]['udg'] = f"[{node}][][]"
        node_names[node] = str(idx)
        idx += 1

    graph = nx.relabel_nodes(graph, node_names)

    return graph

60
61


62

Yoann Dufresne's avatar
Yoann Dufresne committed
63
64
# ------------- Path Graph -------------

Yoann Dufresne's avatar
Yoann Dufresne committed
65
66
67
68
def mols_from_node(node_name):
    return [int(idx) for idx in node_name.split(":")[1].split(".")[0].split("_")]


69
70
71
72
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
73
74
        :param graph: The networkx graph representing the deconvolved graph
        :return: A tuple containing two dictionaries. The first one with theoretical frequencies of each node, the second one with observed frequencies.
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
    """
    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
93
def parse_path_graph_frequencies(graph, barcode_graph):
94
95
96
97
    """ 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
98
        :param barcode_graph: The barcode graph
99
100
        :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
101
    # Compute origin nodes formatted as `{idx}:{mol1_id}_{mol2_id}_...`
Yoann Dufresne's avatar
Yoann Dufresne committed
102
    observed_frequencies = {}
Yoann Dufresne's avatar
Yoann Dufresne committed
103
    real_frequencies = {}
104
    origin_node_names = []
Yoann Dufresne's avatar
Yoann Dufresne committed
105
106
    node_per_barcode = {}

Yoann Dufresne's avatar
Yoann Dufresne committed
107
    for node, data in graph.nodes(data=True):
108
109
        parsed = parse_dg_name(graph, node)
        origin_name = parsed[0][1]
Yoann Dufresne's avatar
Yoann Dufresne committed
110
111

        if origin_name not in node_per_barcode:
Yoann Dufresne's avatar
Yoann Dufresne committed
112
113
            node_per_barcode[origin_name] = []
        node_per_barcode[origin_name].append(node)
114
115

        # Count frequency
Yoann Dufresne's avatar
Yoann Dufresne committed
116
        if origin_name not in observed_frequencies:
Yoann Dufresne's avatar
Yoann Dufresne committed
117
            observed_frequencies[origin_name] = 0
118
            origin_node_names.append(origin_name)
Yoann Dufresne's avatar
Yoann Dufresne committed
119
        observed_frequencies[origin_name] += 1
120

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

Yoann Dufresne's avatar
Yoann Dufresne committed
124
    return real_frequencies, observed_frequencies, node_per_barcode
125
126


Yoann Dufresne's avatar
Yoann Dufresne committed
127
def parse_graph_path(graph):
128
129
130
    """ This function aims to look for direct molecule neighbors.
        If a node has more than 2 direct neighbors, it's not rightly split
    """
Yoann Dufresne's avatar
Yoann Dufresne committed
131
132
133
134
135
    neighborhood = {}

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

Yoann Dufresne's avatar
Yoann Dufresne committed
137
138
139
140
141
142
143
144
        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)
145

Yoann Dufresne's avatar
Yoann Dufresne committed
146
    return neighborhood
147
148


Yoann Dufresne's avatar
Yoann Dufresne committed
149
150
def print_path_summary(frequencies, light_print=False, file_pointer=sys.stdout):
    if file_pointer is None:
Yoann Dufresne's avatar
Yoann Dufresne committed
151
152
        return

Yoann Dufresne's avatar
Yoann Dufresne committed
153
    print("--- Nodes analysis ---", file=file_pointer)
Yoann Dufresne's avatar
Yoann Dufresne committed
154
155
156
    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
157
158
159
160
161
        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
162
        if light_print and obs == the:
Yoann Dufresne's avatar
Yoann Dufresne committed
163
164
165
166
167
            continue

        print(result, file=file_pointer)


168
169
    print("--- Global summary ---", file=file_pointer)

Yoann Dufresne's avatar
Yoann Dufresne committed
170
171
    # --- Frequency usage ---
    # Tags
Yoann Dufresne's avatar
Yoann Dufresne committed
172
173
174
    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
175
    # molecules
Yoann Dufresne's avatar
Yoann Dufresne committed
176
177
178
    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
179
180
181
    # Wrong splits
    over_split = 0
    under_split = 0
Yoann Dufresne's avatar
Yoann Dufresne committed
182
183
184
185
186
    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
187
    print(f"Under/Over splitting: {under_split} - {over_split}")
188
189


190
191
192
193
def print_dgraphs_summary(frequencies, light_print=False):
    pass


Yoann Dufresne's avatar
Yoann Dufresne committed
194
195
# ------------- D2 Graph -------------

196
197
def str_to_udg_lists(s):
    udg = s.replace("]", "").replace(' [', '[')
198
    return udg.split('[')[1:]
199

Rayan  CHIKHI's avatar
Rayan CHIKHI committed
200
201
202
203
# speeds up networkx access to attributes
cached_udg_attr = None
cached_score_attr = None

204
def parse_dg_name(gr, name):
Rayan  CHIKHI's avatar
Rayan CHIKHI committed
205
206
207
208
    global cached_udg_attr, cached_score_attr 
    if cached_udg_attr is None:
        cached_udg_attr = nx.get_node_attributes(gr, 'udg')
    udg = cached_udg_attr[name]
Rayan Chikhi's avatar
Rayan Chikhi committed
209
210
211
212
    res = str_to_udg_lists(udg)
    if len(res) != 3:
        print("parsing problem:",res)
    central, h1, h2 = res
Yoann Dufresne's avatar
Yoann Dufresne committed
213
    
214
    idx = name
Rayan  CHIKHI's avatar
Rayan CHIKHI committed
215
216
217
    if cached_score_attr is None:
        cached_score_attr = nx.get_node_attributes(gr, 'score')
    score = cached_score_attr[name]
Yoann Dufresne's avatar
Yoann Dufresne committed
218
219
220
221
222
223
224
225

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

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


226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
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


Yoann Dufresne's avatar
Yoann Dufresne committed
246
def print_d2_summary(connected_components, longest_path, coverage_vars=(0, 0), light_print=False):
Yoann Dufresne's avatar
Yoann Dufresne committed
247
248
249
    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
250
251
    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
252
253
254
    print(f"The 5 largest components: {[len(x) for x in connected_components][:5]}")

    print("--- Largest component analysis ---")
255
256
257
258
259
260
261
262
263
    # 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))
Yoann Dufresne's avatar
Yoann Dufresne committed
264
265
266

    print(f"Number of usable coverage variables: {len(coverage_vars[1])}")
    print(f"Coverage: {len(coverage_vars[0])}/{len(coverage_vars[1])}")
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
    if not light_print:
        print(f"Missing coverage variables:\n{coverage_vars[1]-coverage_vars[0]}")

def _get_distant_neighbors(graph, node, dist):
    neighbors = set()

    to_compute = [node]
    for _ in range(dist):
        next_compute = []
        for node in to_compute:
            for neighbor in graph[node]:
                if neighbor not in neighbors:
                    neighbors.add(neighbor)
                    next_compute.append(neighbor)
        to_compute = next_compute
Yoann Dufresne's avatar
Yoann Dufresne committed
282

283
    return neighbors
Yoann Dufresne's avatar
Yoann Dufresne committed
284

285
def compute_next_nodes(d2_component, max_jumps=0):
286
287
288
289
290
    # First parse dg names
    dg_names = {}
    for node in d2_component.nodes():
        dg_names[node] = parse_dg_name(d2_component,node)

Yoann Dufresne's avatar
Yoann Dufresne committed
291
292
293
294
    next_nodes = {}

    for node in d2_component.nodes():
        # Parse the current node name
295
        head, h1, h2 = dg_names[node]
Yoann Dufresne's avatar
Yoann Dufresne committed
296
297
298
299
        next_nodes[node] = {}

        # Get the current molecule idxs
        molecule_idxs = mols_from_node(head[1])
Rayan  CHIKHI's avatar
Rayan CHIKHI committed
300
301
        #print("node",node,"dg name",dg_names[node],"mol idxs",molecule_idxs)

Yoann Dufresne's avatar
Yoann Dufresne committed
302
303
        for mol_idx in molecule_idxs:
            nexts = []
304
305
            # for neighbor in d2_component[node]:
            for neighbor in _get_distant_neighbors(d2_component, node, max_jumps+1):
Rayan  CHIKHI's avatar
Rayan CHIKHI committed
306
                # nei_head: central node of the neighbor of 'node'
307
                nei_head, _, _ = dg_names[neighbor]
Yoann Dufresne's avatar
Yoann Dufresne committed
308
                nei_mols = mols_from_node(nei_head[1])
Rayan  CHIKHI's avatar
Rayan CHIKHI committed
309
                # only consider neighbor molecules that are strictly bigger than the current molecule idx considered (from 'node')
Yoann Dufresne's avatar
Yoann Dufresne committed
310
311
312
313
314
                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)
Rayan  CHIKHI's avatar
Rayan CHIKHI committed
315
                    # append to the neighbors of (node,mol_idx) that 'neighbor' if it contains a molecule that's bigger than mol_idx
Yoann Dufresne's avatar
Yoann Dufresne committed
316
317
318
319
                    nexts.append((next_nei_mol, neighbor))

            nexts.sort(key=lambda x: x[0])
            next_nodes[node][mol_idx] = nexts
Yoann Dufresne's avatar
Yoann Dufresne committed
320
            # print("next nodes of node",node,"mol idx",mol_idx,":",next_nodes)
Yoann Dufresne's avatar
Yoann Dufresne committed
321
322
323
324

    return next_nodes


325
326
def compute_longest_increasing_paths(d2_component, max_gap=0):
    next_nodes = compute_next_nodes(d2_component, max_jumps=max_gap)
Yoann Dufresne's avatar
Yoann Dufresne committed
327
328
329
330
331
332

    # 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]:
333
            recursive_longest_path(start_node, mol_idx, next_nodes, longest_paths)
Yoann Dufresne's avatar
Yoann Dufresne committed
334

Rayan  CHIKHI's avatar
Rayan CHIKHI committed
335
336
    # Get the longest path size,
    # across all barcode graph nodes and all molecules in these barcodes
Yoann Dufresne's avatar
Yoann Dufresne committed
337
338
339
340
341
342
343
344
345
346
347
348
349
350
    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


Yoann Dufresne's avatar
Yoann Dufresne committed
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
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
def compute_shortest_edit_path(path):
    min_mol = float("inf")
    max_mol = 0
    molecules = {}

    path_extremes = []
    node_names = {}

    # Parse molecules
    for node_name, data in path.nodes(data=True):
        udg = data["udg"]
        central_node = udg.split(']')[0][1:]

        node_names[node_name] = central_node
        if len(list(path[node_name])) == 1:
            path_extremes.append(node_name)

        mol_names = central_node.split(":")[1].split('_')
        for mol_name in mol_names:
            mol_idx = int(mol_name)
            if mol_idx < min_mol: min_mol = mol_idx
            if mol_idx >= max_mol: max_mol = mol_idx
            molecules[mol_idx] = central_node

    # create barcode path from molecules
    molecule_order = []
    for idx in range(min_mol, max_mol+1):
        if idx in molecules:
            molecule_order.append(molecules[idx])

    # create barcode path from d2_path
    first = path_extremes[0]
    last = path_extremes[1]
    d2_path_order = [first, list(path[first])[0]]
    while d2_path_order[-1] != last:
        neighbors = list(path[d2_path_order[-1]])

        if neighbors[0] == d2_path_order[-2]:
            d2_path_order.append(neighbors[1])
        else:
            d2_path_order.append(neighbors[0])
    d2_path_order = [node_names[x] for x in d2_path_order]

    edit_path = edit_distance(molecule_order, d2_path_order)
    filtered_edit_path = [x[0] for x in edit_path if x[0] == x[1]]
    reverse_edit_path = edit_distance(molecule_order, d2_path_order[::-1])
Yoann Dufresne's avatar
Yoann Dufresne committed
397
    reverse_filtered_edit_path = [x[0] for x in reverse_edit_path if x[0] == x[1]]
Yoann Dufresne's avatar
Yoann Dufresne committed
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453

    if len(filtered_edit_path) > len(reverse_filtered_edit_path):
        return filtered_edit_path
    else:
        return reverse_filtered_edit_path



def edit_distance(array_vertical, array_horizontal):
    dynamic = [[float("inf") for column in range(len(array_horizontal)+1)] for row in range(len(array_vertical)+1)]
    # Fill init line and column
    for i in range(len(array_horizontal)+1):
        dynamic[0][i] = i
    for i in range(len(array_vertical)+1):
        dynamic[i][0] = i

    # Fill the array
    for row in range(1, len(array_vertical)+1):
        for column in range(1, len(array_horizontal)+1):
            if array_vertical[row-1] == array_horizontal[column-1]:
                dynamic[row][column] = min(dynamic[row-1][column-1], dynamic[row-1][column]+1, dynamic[row][column-1]+1)
            else:
                dynamic[row][column] = min(dynamic[row-1][column-1]+1, dynamic[row-1][column]+1, dynamic[row][column-1]+1)

    # Compute alignment
    row = len(array_vertical)
    column = len(array_horizontal)
    path = [(array_vertical[row-1], array_horizontal[column-1])]

    while row != 0 and column != 0:
        score = dynamic[row][column]
        # Match
        if array_vertical[row-1] == array_horizontal[column-1] and dynamic[row-1][column-1] == dynamic[row][column]:
            row -= 1
            column -= 1
        elif array_vertical[row-1] != array_horizontal[column-1] and dynamic[row-1][column-1] == dynamic[row][column] - 1:
            row -= 1
            column -= 1
        elif dynamic[row-1][column] == dynamic[row][column] - 1:
            row -= 1
        elif dynamic[row][column-1] == dynamic[row][column] - 1:
            column -= 1
        else:
            print("Huge problem in edit distance !", file=sys.stderr)
            return []

        path.append((array_vertical[row-1], array_horizontal[column-1]))

    for finalize_row in range(row-1, -1, -1):
        path.append((array_vertical[finalize_row], array_horizontal[column]))
    for finalize_column in range(column-1, -1, -1):
        path.append((array_vertical[row], array_horizontal[finalize_column]))

    return path[::-1]


Yoann Dufresne's avatar
Yoann Dufresne committed
454
def backtrack_longest_path(node, molecule, longest_paths, path=[]):
455
    if node is None:
Yoann Dufresne's avatar
Yoann Dufresne committed
456
457
        return path

458
    path.append((molecule, node))
Yoann Dufresne's avatar
Yoann Dufresne committed
459
460
461
462
463
464
465
    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]:
Rayan  CHIKHI's avatar
Rayan CHIKHI committed
466
        #print("getting cached result for node",current_node,"mol",current_molecule,longest_paths[current_node][current_molecule])
Yoann Dufresne's avatar
Yoann Dufresne committed
467
468
469
470
        return longest_paths[current_node][current_molecule]

    longest = 0
    longest_next = None
471
    min_mol_idx = float('inf') 
Yoann Dufresne's avatar
Yoann Dufresne committed
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
    # 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]


492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
# def longest_common_subsequence(barcode_true_path, barcoded_graph):
#     """ Assume that the two graphs have an attribute barcode for each node and a unique node name"""
#     path_nodes = []
#     path_nodes_barcodes = []
#     for node, data in barcode_true_path.nodes(data=True):
#         path_nodes.append(node)
#         path_nodes_barcodes.append(data["barcode"])
#     path_nodes_to_idx = {n: idx for idx, n in enumerate(path_nodes)}
#
#     graph_nodes = []
#     graph_nodes_barcodes = []
#     for node, data in barcoded_graph.nodes(data=True):
#         graph_nodes.append(node)
#         graph_nodes_barcodes.append(data["barcode"])
#     graph_nodes_to_idx = {n: idx for idx, n in enumerate(graph_nodes)}
#
#     dynamic_array = [[0 for _ in range(len(graph_nodes)+1)] for _ in range(len(path_nodes)+1)]
#     for row in range(1, len(path_nodes)+1):
#         path_node = path_nodes[row-1]
#         path_barcode = path_nodes_barcodes[row-1]
#
#         for column in range(1, len(graph_nodes)):
#             graph_node = graph_nodes[column-1]
#             graph_barcode = graph_nodes_barcodes[column-1]
#
#             prev_scores = [dynamic_array[row-1][column]]
#             for neighbor_node in barcoded_graph[graph_node]:
#                 neighbor_idx = graph_nodes_to_idx[neighbor_node]
#                 prev_scores.append(dynamic_array[row-1][neighbor_idx])
#
#             match_point = 1 if path_barcode == graph_barcode else 0
#             dynamic_array[row][column] = max(prev_scores) + match_point


Yoann Dufresne's avatar
Yoann Dufresne committed
526
527
528
529
def compute_covered_variables(graph, path):
    path_nodes = set()
    for mol, node_name in path:
        path_nodes.add(node_name)
530

Yoann Dufresne's avatar
Yoann Dufresne committed
531
532
533
534
535
536
537
    used_vars = set()
    total_vars = set()
    for node, data in graph.nodes(data=True):
        vars = data["barcode_edges"].split(" ")
        total_vars.update(vars)
        if node in path_nodes:
            used_vars.update(vars)
538

Yoann Dufresne's avatar
Yoann Dufresne committed
539
    return used_vars, total_vars
540

541
# returns True iff there exist x in mol1 such that there exists y in mol2 and |x-y| <= some_value
542
543
544
545
546
547
548
549
550
551
552
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
553
        head, c1, c2 = parse_dg_name(d2_component,node)
554
555

        # Construct the molecule(s) that this udg really 'reflects'
556
557
558
559
560
        # i.e. the udg has a central node and two cliques
        # 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
        # (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
561
562
563
564
        udg_molecules = set()

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

        # Examine molecule idx's of cliques to see which are close to the central node
568
        # rationale: c1/c2 contain nearby molecule id's
569
570
571
572
573
        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])
574
                    nei_mols = [x for x in nei_mols if x > mol_idx]  # fixme: also look at the <= molecules for more robustness
575
576
577
578
579
580
581
582
                    
                    # 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)])
583
            if quality > 0.6:  # eyeballed but still arbitrary
584
                udg_molecules.add(mol_idx)
585
            # print("mol",mol_idx,molecule_idxs,"quality",quality,"nexts",nexts)
586
587
588
589
590
591
       
        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
592
        head, c1, c2 = parse_dg_name(d2_component,n1)
593
594
        node_udg_molecules = udg_molecules_dict[head[0]]
        
595
        n_head, n_c1, n_c2 = parse_dg_name(d2_component,n2)
596
597
598
599
600
601
602
        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
603
        # print("edge",node_udg_molecules,neighbor_udg_molecules,color)
604
605
606
607
    
    # also, annotate nodes by their putative molecule found
    for n, data in d2_component.nodes(data=True):
        # Parse the current node name
608
        head, c1, c2 = parse_dg_name(d2_component,n)
609
610
611
        node_udg_molecules = udg_molecules_dict[head[0]]
        data['udg_molecule']= '_'.join(list(map(str,node_udg_molecules)))
        
612
613
614
615
616
617
618
619
    # 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
620
            head, c1, c2 = parse_dg_name(d2_component,n)
621
622
623
            if "_" in data['udg_molecule'] or data['udg_molecule'] == '':
                if "_" in data['udg_molecule']:
                    m1, m2 = list(map(int,data['udg_molecule'].split("_")))
624
                    if abs(m2-m1) < 30: continue  # don't remove that kind of nodes
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
                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
641

642
643
644
645
def main():
    args = parse_args()
    graph = load_graph(args.filename)

Yoann Dufresne's avatar
Yoann Dufresne committed
646
    if args.type == "path":
Yoann Dufresne's avatar
Yoann Dufresne committed
647
648
649
650
        if args.barcode_graph is None:
            print("--barcode_graph is required for path analysis", file=sys.stderr)
            exit(0)

Yoann Dufresne's avatar
Yoann Dufresne committed
651
        barcode_graph = load_graph(args.barcode_graph)
Yoann Dufresne's avatar
Yoann Dufresne committed
652
653
654
        # if len(list(nx.connected_components(graph))) != 1:
        #     print([len(x) for x in list(nx.connected_components(graph))])
        #     exit("when running evaluate.py --type path, the graph should have a single connected component (it's supposed to be a path, after all)")
Rayan  CHIKHI's avatar
Rayan CHIKHI committed
655
656

        # compute LIS over the path
Yoann Dufresne's avatar
Yoann Dufresne committed
657
658
        # longest_path = compute_longest_increasing_paths(graph)
        longest_path = compute_shortest_edit_path(graph)
Rayan  CHIKHI's avatar
Rayan CHIKHI committed
659
660
        print("--- Largest component analysis ---")
        print(f"Size of the longest path: {len(longest_path)}")
Yoann Dufresne's avatar
Yoann Dufresne committed
661
662
        if not args.light_print:
            print(longest_path)
Yoann Dufresne's avatar
Yoann Dufresne committed
663

Rayan  CHIKHI's avatar
Rayan CHIKHI committed
664
665
666
        # get over/under counted molecules
        print("--- Under/over molecule counts ---")
        frequencies = parse_path_graph_frequencies(graph, barcode_graph)
Yoann Dufresne's avatar
Yoann Dufresne committed
667
        print_path_summary(frequencies, light_print=args.light_print)
Yoann Dufresne's avatar
Yoann Dufresne committed
668
        print(f"Size of the longest path: {len(longest_path)}")
669
670
671
    elif args.type == "dgraphs":
        udg_per_node = parse_udg_qualities(graph)
        # print(udg_per_node)
Yoann Dufresne's avatar
Yoann Dufresne committed
672
673
674
675
676
    elif args.type == "bg" or args.type == "d2":

        if args.type == "bg":
            graph = transform_bg(graph)

Yoann Dufresne's avatar
Yoann Dufresne committed
677
678
679
680
        components = list(nx.connected_components(graph))
        components.sort(key=lambda x: -len(x))

        component = graph.subgraph(components[0])
681
        longest_path = compute_longest_increasing_paths(component, max_gap=args.max_gap)
Yoann Dufresne's avatar
Yoann Dufresne committed
682
683
        vars = compute_covered_variables(graph, longest_path)
        print_d2_summary(components, longest_path, coverage_vars=vars, light_print=args.light_print)
684
685
686
687
688
689
690
691
692
    
    # 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])
 
693
        component = verify_graph_edges(component)
694

Yoann Dufresne's avatar
Yoann Dufresne committed
695
696
697
698
        extension = args.filename.split('.')[-1]
        base_filename = '.'.join(args.filename.split('.')[:-1])
        save_graph(component, base_filename+".verified."+extension)

699
700
701

if __name__ == "__main__":
    main()