ChromAlignNet/utilsAlignment.py at master · mili7522/ChromAlignNet · 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
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
import numpy as np
import matplotlib.pyplot as plt
import scipy.spatial
import scipy.cluster


### Group and cluster
def getDistances(prediction):
    """
    Turns a vector of probabilities into a vector of distances

    Arguments:
        prediction -- Numpy array giving a column vector of probabilities (from main output of ChromAlignNet)

    Returns:
        distances -- Numpy array giving a column vector of distances
    """
    distances = 1 / prediction
    return distances


def getDistanceMatrix(comparisons, number_of_peaks, prediction, clip = 10, info_df = None):
    """
    Produces a matrix of distances between all peaks

    Arguments:
        comparisons -- Numpy array with two columns - x1 and x2 - containing the IDs of the two peaks being compared
        number_of_peaks -- Int: Total number of peaks under consideration
        prediction -- Numpy array giving a column vector of probabilities (from main output of ChromAlignNet)
        clip -- Int or Float: Maximum value of the distance matrix
        info_df -- DataFrame containing information about each peak, in particular the peak times and file number

    Returns:
        distance_matrix -- 2D numpy array of distances between all peaks
    """
    distances = getDistances(prediction)

    distance_matrix = np.empty((number_of_peaks, number_of_peaks))
    distance_matrix.fill(clip)  # Clip value

    probability_matrix = np.zeros((number_of_peaks, number_of_peaks))

    for i, (x1, x2) in enumerate(comparisons):
        if info_df is not None and info_df.loc[x1, 'File'] == info_df.loc[x2, 'File']:
            val = min(distances[i] * 2, clip)
        elif info_df is not None and np.abs(info_df.loc[x1, 'peakMaxTime'] - info_df.loc[x2, 'peakMaxTime']) > 0.5:
            val = min(distances[i] * 2, clip)
        else:
            val = min(distances[i], clip)
        distance_matrix[x1, x2] = distance_matrix[x2, x1] = val
        probability_matrix[x1, x2] = probability_matrix[x2, x1] = prediction[i]

    for i in range(number_of_peaks):
        distance_matrix[i,i] = 0
        probability_matrix[i,i] = 1

    return distance_matrix, probability_matrix


def assignGroups(distance_matrix, threshold = 2, plot_dendrogram = False):
    """
    Assigns a group number to peaks based on what should be aligned together

    Arguments:
        distance_matrix -- 2D numpy array of distances between all peaks
        threshold -- Int or Float: Value to cut the dendrogram to obtain the clusters
        plot_dendrogram -- Boolean: To create a new figure for a plot of the dendrogram

    Returns:
        groups -- Dictionary with group ID as keys and sets of peak IDs for each group
    """
    sqform = scipy.spatial.distance.squareform(distance_matrix)
    mergings = scipy.cluster.hierarchy.linkage(sqform, method = 'average')
    if plot_dendrogram:
        plt.figure()
        scipy.cluster.hierarchy.dendrogram(mergings, leaf_font_size = 3, color_threshold = threshold)
    labels = scipy.cluster.hierarchy.fcluster(mergings, threshold, criterion = 'distance')

    groups = {}
    for i in range(max(labels)):
        groups[i] = set(np.where(labels == i + 1)[0])  # labels start at 1

    return groups

def postprocessGroups(groups, info_df):
    """
    Adjusts the groups assigned to each peak

    Arguments:
        groups -- Dictionary with group ID as keys and sets of peak IDs for each group
        info_df -- DataFrame containing information about each peak, in particular the peak times and file number

    Returns:
        new_groups -- Dictionary with new group ID as keys and sets of peak IDs for each group
    """
    max_group = len(groups) - 1  # max(groups.keys())
    new_groups = {}
    for i, group in groups.items():
        group_df = info_df.loc[group].copy()
        group_df.sort_values(by = 'peakMaxTime', axis = 0, inplace = True)
        files = group_df['File']
        files_count = dict()
        new_groups[i] = set()
        max_group_increment = 0
        for peak, file in files.iteritems():
            if file in files_count:
                if max_group + files_count[file] not in new_groups:
                    new_groups[max_group + files_count[file]] = set()
                new_groups[max_group + files_count[file]].add(peak)
                if files_count[file] > max_group_increment:
                    max_group_increment += 1
                files_count[file] += 1
            else:
                files_count[file] = 1
                new_groups[i].add(peak)

        max_group += max_group_increment

    return new_groups


def alignTimes(groups, info_df, peak_intensity, align_to):
    """
    Creates a column in the info_df Dataframe with the aligned times of all peaks.
    Aligns the peaks in each group according to the average time of their peakMaxTime, weighted by the peak intensity

    Arguments:
        groups -- Dictionary with group ID as keys and sets of peak IDs for each group
        info_df -- DataFrame containing information about each peak, in particular the peak times
        peak_intensity -- pandas Series of peak intensities
        align_to -- String: Giving the name of the new column in the info_df Dataframe which will contain the aligned times
    """
    info_df[align_to] = info_df['peakMaxTime']
    for group in groups.values():
        times = info_df.loc[group, 'peakMaxTime']
        peak_values = peak_intensity.loc[group]
        average_time = np.average(times, weights = peak_values)
#        average_time = np.mean(times)
        info_df.loc[group, align_to] = average_time