edison/ordering.py at main · Noble-Lab/edison · 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
def get_edges(position_table: dict, use_orientation: bool = False) -> list:
    """
    Get the edges present in an assembly based on the position table of the
    adjacency graph.

    Args:
        position_table: The position table from the create_table() step.
        use_orientation: Boolean, we can reuse this function between ordering
            and orientation.

    Returns:
        A list of tuples which describe all contig adjacencies in an assembly or
        reference.
    """
    edges = []
    for index in position_table:
        edge = set()
        for contig, orientation in position_table[index]:
            if use_orientation:
                contig = (contig, orientation)
            edge.add(contig)

        # Skip edges with telomeres.
        if len(edge) < 2:
            continue

        edges.append(edge)

    return edges


def count(first: tuple, second: tuple, orientation: bool = False) -> float:
    """
    Compute the length weighted intersection of two edge lists.

    Args:
        first: The reference genome as (contig_table, position_table).
        second: The assembly as (contig_table, position, table).

    Returns:
        A fraction which represents how many contigs were correctly ordered or
        oriented.
    """
    first_contigs, first_positions = first
    second_contigs, second_positions = second

    first_connections = get_edges(first_positions, orientation)
    second_connections = get_edges(second_positions, orientation)

    weighted_first = 0
    weighted_second = 0
    for edge in first_connections:

        # Add the edge weight.
        current_weight = 0
        for contig in edge:
            if orientation:
                contig = contig[0]
            current_weight += first_contigs[contig]["length"]
        weighted_first += current_weight

        # If the edge is in the other assembly, count it and weight it.
        if edge in second_connections:
            weighted_second += current_weight

    if weighted_first == 0:
        return 0

    return weighted_second / weighted_first