File size: 16,060 Bytes
0dab632
 
 
 
 
 
 
 
 
 
 
016c3b8
0dab632
016c3b8
0dab632
 
016c3b8
 
 
 
 
 
0dab632
 
016c3b8
 
 
0dab632
 
016c3b8
 
0dab632
016c3b8
 
 
 
0dab632
016c3b8
0dab632
 
016c3b8
0dab632
016c3b8
0dab632
 
016c3b8
 
0dab632
016c3b8
 
0dab632
 
016c3b8
 
 
 
 
 
 
0dab632
 
 
 
 
 
 
 
 
 
 
016c3b8
 
 
 
 
 
6756e43
016c3b8
 
6756e43
016c3b8
0dab632
 
 
 
 
 
 
6756e43
016c3b8
0dab632
6756e43
0dab632
 
 
 
016c3b8
0dab632
6756e43
0dab632
 
016c3b8
 
 
 
 
0dab632
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
6756e43
0dab632
 
 
 
 
 
 
 
 
 
 
 
 
6756e43
0dab632
 
 
 
 
 
 
 
 
 
 
6756e43
0dab632
 
6756e43
0dab632
 
6756e43
0dab632
 
 
6756e43
0dab632
 
 
 
 
 
 
 
 
 
 
016c3b8
0dab632
 
6756e43
0dab632
 
 
016c3b8
 
 
 
 
 
 
 
 
 
 
 
 
 
6756e43
016c3b8
 
 
 
0dab632
 
 
 
6756e43
0dab632
 
6756e43
0dab632
 
016c3b8
0dab632
6756e43
0dab632
 
 
 
 
 
 
 
 
 
016c3b8
6756e43
016c3b8
0dab632
 
 
 
 
 
 
 
 
 
 
 
 
016c3b8
 
 
 
 
0dab632
 
 
 
 
 
 
 
 
 
 
 
016c3b8
0dab632
 
 
6756e43
016c3b8
 
0dab632
 
 
 
6756e43
0dab632
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
6756e43
0dab632
 
016c3b8
 
 
 
 
 
 
0dab632
016c3b8
6756e43
016c3b8
 
 
 
 
6756e43
016c3b8
 
 
 
6756e43
 
016c3b8
6756e43
016c3b8
6756e43
016c3b8
 
0dab632
6756e43
0dab632
016c3b8
6756e43
0dab632
 
 
 
 
 
 
 
 
 
6756e43
0dab632
 
 
 
 
 
 
 
6756e43
0dab632
 
016c3b8
 
 
 
0dab632
6756e43
0dab632
016c3b8
6756e43
016c3b8
6756e43
 
016c3b8
6756e43
016c3b8
6756e43
0dab632
016c3b8
0dab632
6756e43
0dab632
016c3b8
6756e43
0dab632
 
 
 
6756e43
0dab632
 
 
6756e43
0dab632
 
 
 
 
 
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
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
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
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
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
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396


from dataclasses import dataclass
import numpy as np
import scipy.linalg as la
from scipy.signal import find_peaks
from math import ceil




def thin_peaks(peak_list, dmin=10, voxel_size=(1,1,1), return_larger_peaks=False):
    """
    Remove peaks within a specified distance of each other, retaining the peak with the highest intensity.

    Args:
    - peak_list (list of PeakData): Each element contains:
        - pos (list of float): 3D coordinates of the peak.
        - intensity (float): The intensity value of the peak.
        - key (tuple): A unique identifier or index for the peak (#trace, #peak)
    - dmin (float, optional): Minimum distance between peaks. peaks closer than this threshold will be thinned. Defaults to 10.
    - return_larger_peaks (bool, optional): Indicate larger peak for each thinned peak

    Returns:
    - list of tuples: A list containing keys of the removed peaks.
    if return_larger_peaks
    - list of tuples: A list containing the keys of the larger peak causing the peak to be removed

    Notes:
    - The function uses the L2 norm (Euclidean distance) to compute the distance between peaks.
    - When two peaks are within `dmin` distance, the peak with the lower intensity is removed.
    """
    removed_peaks = []
    removed_larger_peaks = []
    for i in range(len(peak_list)):
        if peak_list[i].key in removed_peaks:
            continue
        for j in range(len(peak_list)):
            if i==j:
                continue
            if peak_list[j].key in removed_peaks:
                continue
            d = (np.array(peak_list[i].pos) - np.array(peak_list[j].pos))*np.array(voxel_size)
            d = la.norm(d)
            if d<dmin:
                hi = peak_list[i].intensity
                hj = peak_list[j].intensity
                if hi<hj:
                    removed_peaks.append(peak_list[i].key)
                    removed_larger_peaks.append(peak_list[j].key)
                    break
                else:
                    removed_peaks.append(peak_list[j].key)
                    removed_larger_peaks.append(peak_list[i].key)

    if return_larger_peaks:
        return removed_peaks, removed_larger_peaks
    else:
        return removed_peaks


@dataclass
class CellData(object):
    """Represents data related to a single cell.

    Attributes:
        pathdata_list (list): A list of PathData objects representing the various paths associated with the cell.
    """
    pathdata_list: list

@dataclass
class RemovedPeakData(object):
    """Represents data related to a removed peak

    Attributes:
        idx (int): Index of peak along path
        screening_peak (tuple): (path_idx, position along path) for screening peak
    """
    idx: int 
    screening_peak: tuple 

@dataclass
class PathData(object):
    """Represents data related to a specific path in the cell.

    This dataclass encapsulates information about the peaks, 
    the defining points, the fluorescence values, and the path length of a specific path.

    Attributes: peaks (list): List of peaks in the path (indicies of positions in points, o_intensity).
        removed_peaks (list): List of peaks in the path which have been removed because of a nearby larger peak
        points (list): List of points defining the path.
        o_intensity (list): List of (unnormalized) fluorescence intensity values along the path
        SC_length (float): Length of the path.

    """
    peaks: list
    removed_peaks: list
    points: list
    o_intensity: list
    SC_length: float

@dataclass
class PeakData(object):
    pos: tuple
    intensity: float
    key: tuple


def find_peaks2(v, distance=5,  prominence=0.5):
    """
    Find peaks in a 1D array with extended boundary handling.

    The function pads the input array at both ends to handle boundary peaks. It then identifies peaks in the extended array
    and maps them back to the original input array.

    Args:
    - v (numpy.ndarray): 1D input array in which to find peaks.
    - distance (int, optional): Minimum number of array elements that separate two peaks. Defaults to 5.
    - prominence (float, optional): Minimum prominence required for a peak to be identified. Defaults to 0.5.

    Returns:
    - list of int: List containing the indices of the identified peaks in the original input array.
    - dict: Information about the properties of the identified peaks (as returned by scipy.signal.find_peaks).

    """
    pad = int(ceil(distance))+1
    v_ext = np.concatenate([np.ones((pad,), dtype=v.dtype)*np.min(v), v, np.ones((pad,), dtype=v.dtype)*np.min(v)])

    assert(len(v_ext) == len(v)+2*pad)
    peaks, _ = find_peaks(v_ext, distance=distance, prominence=prominence)
    peaks = peaks - pad
    n_peaks = []
    for i in peaks:
        if 0<=i<len(v):
            n_peaks.append(i)
        else:
            raise Exception
    return n_peaks, _
          

def process_cell_traces(all_paths, path_lengths, measured_trace_fluorescence, dmin=10):
    """
    Process traces of cells to extract peak information and organize the data.

    The function normalizes fluorescence data, finds peaks, refines peak information, 
    removes unwanted peaks that might be due to close proximity of bright peaks from 
    other paths, and organizes all the information into a structured data format.

    Args:
        all_paths (list of list of tuples): A list containing paths, where each path is 
                                            represented as a list of 3D coordinate tuples.
        path_lengths (list of float): List of path lengths corresponding to the provided paths.
        measured_trace_fluorescence (list of list of float): A list containing fluorescence 
                                                            data corresponding to each path point.
        dmin (float): Distance below which brighter peaks screen less bright ones.

    Returns:
        CellData: An object containing organized peak and path data for a given cell.

    Note:
        - The function assumes that each path and its corresponding length and fluorescence data 
          are positioned at the same index in their respective lists.
    """
    
    cell_peaks = []

    for points, o_intensity in zip(all_paths, measured_trace_fluorescence):
                     
        # For peak determination normalize each trace to have mean zero and s.d. 1
        intensity_normalized = (o_intensity - np.mean(o_intensity))/np.std(o_intensity)
        
        # Find peaks - these will be further refined later
        p,_ = find_peaks2(intensity_normalized, distance=5,  prominence=0.5*np.std(intensity_normalized))
        peaks = np.array(p, dtype=np.int32)

        # Store peak data - using original values, not normalized ones
        peak_mean_heights = [ o_intensity[u] for u in peaks ]
        peak_points = [ points[u] for u in peaks ]
        
        cell_peaks.append((peaks, peak_points, peak_mean_heights))
        
    # Eliminate peaks which have another larger peak nearby (in 3D space, on any chromosome).
    # This aims to remove small peaks in the mean intensity generated when an  SC passes close
    # to a bright peak on another SC - this is nearby in space, but brighter.

    to_thin = []
    for k in range(len(cell_peaks)):
        for u in range(len(cell_peaks[k][0])):
            to_thin.append(PeakData(pos=cell_peaks[k][1][u], intensity=cell_peaks[k][2][u], key=(k, u)))
    
    # Exclude any peak with a nearby brighter peak (on any SC)
    removed_peaks, removed_larger_peaks = thin_peaks(to_thin, return_larger_peaks=True, dmin=dmin)

    # Clean up and remove these peaks
    new_cell_peaks = []
    removed_cell_peaks = []
    removed_cell_peaks_larger = []
    for path_idx in range(len(cell_peaks)):
        path_retained_peaks = []
        path_removed_peaks = []
        path_peaks = cell_peaks[path_idx][0]

        for peak_idx in range(len(path_peaks)):
            if (path_idx, peak_idx) not in removed_peaks:
                path_retained_peaks.append(path_peaks[peak_idx])
            else:
                # What's the larger point?
                idx = removed_peaks.index((path_idx, peak_idx))
                larger_path, larger_idx = removed_larger_peaks[idx]        
                path_removed_peaks.append(RemovedPeakData(idx=path_peaks[peak_idx], screening_peak=(larger_path, cell_peaks[larger_path][0][larger_idx])))
                ###
                 
        new_cell_peaks.append(path_retained_peaks)
        removed_cell_peaks.append(path_removed_peaks)
        
    cell_peaks = new_cell_peaks
    pd_list = []
    
    # Save peak positions, absolute intensity intensities, and length for each SC
    for k in range(len(all_paths)):
        
        points, o_intensity = all_paths[k], measured_trace_fluorescence[k]

        peaks = cell_peaks[k]
        removed_peaks = removed_cell_peaks[k]
        
        pd = PathData(peaks=peaks, removed_peaks=removed_peaks, points=points, o_intensity=o_intensity, SC_length=path_lengths[k])
        pd_list.append(pd)

    cd = CellData(pathdata_list=pd_list)

    return cd


alpha_max = 0.4


# Criterion used for identifying peak as a focus - normalized (with mean and s.d.)
# intensity levels being above 0.4 time maximum peak level
def focus_criterion(pos, v, alpha=alpha_max):
    """
    Identify and return positions where values in the array `v` exceed a certain threshold.

    The threshold is computed as `alpha` times the maximum value in `v`.

    Args:
    - pos (numpy.ndarray): Array of positions.
    - v (numpy.ndarray): 1D array of values, e.g., intensities.
    - alpha (float, optional): A scaling factor for the threshold. Defaults to `alpha_max`.

    Returns:
    - numpy.ndarray: Array of positions where corresponding values in `v` exceed the threshold.
    """
    if len(v):
        idx = (v>=alpha*np.max(v))
        return np.array(pos[idx])
    else:
        return np.array([], dtype=np.int32)

def analyse_celldata(cell_data, config):
    """
    Analyse the provided cell data to extract focus-related information.

    Args:
        cd (CellData): An instance of the CellData class containing path data information.
        config (dictionary): Configuration dictionary containing 'peak_threshold' and 'threshold_type'
                             'peak_threshold' (float) - threshold for calling peaks as foci
                             'threshold_type' (str) = 'per-trace', 'per-foci'

    Returns:
        tuple: A tuple containing:
            - foci_rel_intensity (list): List of relative intensities for the detected foci.
            - foci_pos (list): List of absolute positions of the detected foci.
            - foci_pos_index (list): List of indices of the detected foci.
            - screened_foci_data (list): List of RemovedPeakData indicating positions of removed peaks and the index of the larger peak
            - trace_median_intensities (list): Per-trace median intensity
            - trace_thresholds (list): Per-trace absolute threshold for calling peaks as foci
    """
    foci_abs_intensity = []
    foci_pos = []
    foci_pos_index = []
    screened_foci_data = []
    trace_median_intensities = []
    trace_thresholds = []
    
    peak_threshold = config['peak_threshold']

    threshold_type = config['threshold_type']

    if threshold_type == 'per-trace':
        """
        Call extracted peaks as foci if intensity - trace_mean > peak_threshold * (trace_max_foci_intensity - trace_mean)
        """
        
        for path_data in cell_data.pathdata_list:
            peaks = np.array(path_data.peaks, dtype=np.int32)

            # Normalize extracted fluorescent intensities by subtracting mean (and dividing
            # by standard deviation - note that the latter should have no effect on the results).
            h = np.array(path_data.o_intensity)
            h = h - np.mean(h)
            h = h/np.std(h)
            # Extract foci according to criterion
            foci_idx = focus_criterion(peaks, h[peaks], peak_threshold) 
            
            #
            removed_peaks = path_data.removed_peaks
            removed_peaks_idx = np.array([u.idx for u in removed_peaks], dtype=np.int32)

            
            if len(peaks):
                trace_thresholds.append((1-peak_threshold)*np.mean(path_data.o_intensity) + peak_threshold*np.max(np.array(path_data.o_intensity)[peaks]))
            else:
                trace_thresholds.append(None)

            if len(removed_peaks):
                if len(peaks):
                    threshold = (1-peak_threshold)*np.mean(path_data.o_intensity) + peak_threshold*np.max(np.array(path_data.o_intensity)[peaks])
                else:
                    threshold = float('-inf')

                
                removed_peak_heights = np.array(path_data.o_intensity)[removed_peaks_idx]
                screened_foci_idx = np.where(removed_peak_heights>threshold)[0]
                
                screened_foci_data.append([removed_peaks[i] for i in screened_foci_idx])
            else:
                screened_foci_data.append([])
                
            pos_abs = (foci_idx/len(path_data.points))*path_data.SC_length
            foci_pos.append(pos_abs)
            foci_abs_intensity.append(np.array(path_data.o_intensity)[foci_idx])
            
            foci_pos_index.append(foci_idx)
            trace_median_intensities.append(np.median(path_data.o_intensity))
            
    elif threshold_type == 'per-cell':
        """
        Call extracted peaks as foci if intensity - trace_mean > peak_threshold * max(intensity - trace_mean)
        """
        max_cell_intensity = float("-inf")
        for path_data in cell_data.pathdata_list:

            # Normalize extracted fluorescent intensities by subtracting mean (and dividing
            # by standard deviation - note that the latter should have no effect on the results).
            h = np.array(path_data.o_intensity)
            h = h - np.mean(h)
            max_cell_intensity = max(max_cell_intensity, np.max(h))

        for path_data in cell_data.pathdata_list:
            peaks = np.array(path_data.peaks, dtype=np.int32)

            # Normalize extracted fluorescent intensities by subtracting mean (and dividing
            # by standard deviation - note that the latter should have no effect on the results).
            h = np.array(path_data.o_intensity)
            h = h - np.mean(h)

            foci_idx = peaks[h[peaks]>peak_threshold*max_cell_intensity]

            removed_peaks = path_data.removed_peaks
            removed_peaks_idx = np.array([u.idx for u in removed_peaks], dtype=np.int32)

            trace_thresholds.append(np.mean(path_data.o_intensity) + peak_threshold*max_cell_intensity)

            if len(removed_peaks):
                threshold = np.mean(path_data.o_intensity) + peak_threshold*max_cell_intensity

                removed_peak_heights = np.array(path_data.o_intensity)[removed_peaks_idx]
                screened_foci_idx = np.where(removed_peak_heights>threshold)[0]
                
                screened_foci_data.append([removed_peaks[i] for i in screened_foci_idx])
            else:
                screened_foci_data.append([])

            pos_abs = (foci_idx/len(path_data.points))*path_data.SC_length
            foci_pos.append(pos_abs)
            foci_abs_intensity.append(np.array(path_data.o_intensity)[foci_idx])
            
            foci_pos_index.append(foci_idx)
            trace_median_intensities.append(np.median(path_data.o_intensity))            
            
    else:
        raise NotImplementedError
    
    return foci_abs_intensity, foci_pos, foci_pos_index, screened_foci_data, trace_median_intensities, trace_thresholds

def analyse_traces(all_paths, path_lengths, measured_trace_fluorescence, config):
    
    cd = process_cell_traces(all_paths, path_lengths, measured_trace_fluorescence, dmin=config['screening_distance'])

    return analyse_celldata(cd, config)