JFoz/test_path_analysis
0
1 2 3from dataclasses import dataclass4import numpy as np5import scipy.linalg as la6from scipy.signal import find_peaks7from math import ceil8 9 10 11 12def thin_peaks(peak_list, dmin=10, voxel_size=(1,1,1), return_larger_peaks=False):13 """14 Remove peaks within a specified distance of each other, retaining the peak with the highest intensity.15 16 Args:17 - peak_list (list of PeakData): Each element contains:18 - pos (list of float): 3D coordinates of the peak.19 - intensity (float): The intensity value of the peak.20 - key (tuple): A unique identifier or index for the peak (#trace, #peak)21 - dmin (float, optional): Minimum distance between peaks. peaks closer than this threshold will be thinned. Defaults to 10.22 - return_larger_peaks (bool, optional): Indicate larger peak for each thinned peak23 24 Returns:25 - list of tuples: A list containing keys of the removed peaks.26 if return_larger_peaks27 - list of tuples: A list containing the keys of the larger peak causing the peak to be removed28 29 Notes:30 - The function uses the L2 norm (Euclidean distance) to compute the distance between peaks.31 - When two peaks are within `dmin` distance, the peak with the lower intensity is removed.32 """33 removed_peaks = []34 removed_larger_peaks = []35 for i in range(len(peak_list)):36 if peak_list[i].key in removed_peaks:37 continue38 for j in range(len(peak_list)):39 if i==j:40 continue41 if peak_list[j].key in removed_peaks:42 continue43 d = (np.array(peak_list[i].pos) - np.array(peak_list[j].pos))*np.array(voxel_size)44 d = la.norm(d)45 if d<dmin:46 hi = peak_list[i].intensity47 hj = peak_list[j].intensity48 if hi<hj:49 removed_peaks.append(peak_list[i].key)50 removed_larger_peaks.append(peak_list[j].key)51 break52 else:53 removed_peaks.append(peak_list[j].key)54 removed_larger_peaks.append(peak_list[i].key)55 56 if return_larger_peaks:57 return removed_peaks, removed_larger_peaks58 else:59 return removed_peaks60 61 62@dataclass63class CellData(object):64 """Represents data related to a single cell.65 66 Attributes:67 pathdata_list (list): A list of PathData objects representing the various paths associated with the cell.68 """69 pathdata_list: list70 71@dataclass72class RemovedPeakData(object):73 """Represents data related to a removed peak74 75 Attributes:76 idx (int): Index of peak along path77 screening_peak (tuple): (path_idx, position along path) for screening peak78 """79 idx: int 80 screening_peak: tuple 81 82@dataclass83class PathData(object):84 """Represents data related to a specific path in the cell.85 86 This dataclass encapsulates information about the peaks, 87 the defining points, the fluorescence values, and the path length of a specific path.88 89 Attributes: peaks (list): List of peaks in the path (indicies of positions in points, o_intensity).90 removed_peaks (list): List of peaks in the path which have been removed because of a nearby larger peak91 points (list): List of points defining the path.92 o_intensity (list): List of (unnormalized) fluorescence intensity values along the path93 SC_length (float): Length of the path.94 95 """96 peaks: list97 removed_peaks: list98 points: list99 o_intensity: list100 SC_length: float101 102@dataclass103class PeakData(object):104 pos: tuple105 intensity: float106 key: tuple107 108 109def find_peaks2(v, distance=5, prominence=0.5):110 """111 Find peaks in a 1D array with extended boundary handling.112 113 The function pads the input array at both ends to handle boundary peaks. It then identifies peaks in the extended array114 and maps them back to the original input array.115 116 Args:117 - v (numpy.ndarray): 1D input array in which to find peaks.118 - distance (int, optional): Minimum number of array elements that separate two peaks. Defaults to 5.119 - prominence (float, optional): Minimum prominence required for a peak to be identified. Defaults to 0.5.120 121 Returns:122 - list of int: List containing the indices of the identified peaks in the original input array.123 - dict: Information about the properties of the identified peaks (as returned by scipy.signal.find_peaks).124 125 """126 pad = int(ceil(distance))+1127 v_ext = np.concatenate([np.ones((pad,), dtype=v.dtype)*np.min(v), v, np.ones((pad,), dtype=v.dtype)*np.min(v)])128 129 assert(len(v_ext) == len(v)+2*pad)130 peaks, _ = find_peaks(v_ext, distance=distance, prominence=prominence)131 peaks = peaks - pad132 n_peaks = []133 for i in peaks:134 if 0<=i<len(v):135 n_peaks.append(i)136 else:137 raise Exception138 return n_peaks, _139 140 141def process_cell_traces(all_paths, path_lengths, measured_trace_fluorescence, dmin=10):142 """143 Process traces of cells to extract peak information and organize the data.144 145 The function normalizes fluorescence data, finds peaks, refines peak information, 146 removes unwanted peaks that might be due to close proximity of bright peaks from 147 other paths, and organizes all the information into a structured data format.148 149 Args:150 all_paths (list of list of tuples): A list containing paths, where each path is 151 represented as a list of 3D coordinate tuples.152 path_lengths (list of float): List of path lengths corresponding to the provided paths.153 measured_trace_fluorescence (list of list of float): A list containing fluorescence 154 data corresponding to each path point.155 dmin (float): Distance below which brighter peaks screen less bright ones.156 157 Returns:158 CellData: An object containing organized peak and path data for a given cell.159 160 Note:161 - The function assumes that each path and its corresponding length and fluorescence data 162 are positioned at the same index in their respective lists.163 """164 165 cell_peaks = []166 167 for points, o_intensity in zip(all_paths, measured_trace_fluorescence):168 169 # For peak determination normalize each trace to have mean zero and s.d. 1170 intensity_normalized = (o_intensity - np.mean(o_intensity))/np.std(o_intensity)171 172 # Find peaks - these will be further refined later173 p,_ = find_peaks2(intensity_normalized, distance=5, prominence=0.5*np.std(intensity_normalized))174 peaks = np.array(p, dtype=np.int32)175 176 # Store peak data - using original values, not normalized ones177 peak_mean_heights = [ o_intensity[u] for u in peaks ]178 peak_points = [ points[u] for u in peaks ]179 180 cell_peaks.append((peaks, peak_points, peak_mean_heights))181 182 # Eliminate peaks which have another larger peak nearby (in 3D space, on any chromosome).183 # This aims to remove small peaks in the mean intensity generated when an SC passes close184 # to a bright peak on another SC - this is nearby in space, but brighter.185 186 to_thin = []187 for k in range(len(cell_peaks)):188 for u in range(len(cell_peaks[k][0])):189 to_thin.append(PeakData(pos=cell_peaks[k][1][u], intensity=cell_peaks[k][2][u], key=(k, u)))190 191 # Exclude any peak with a nearby brighter peak (on any SC)192 removed_peaks, removed_larger_peaks = thin_peaks(to_thin, return_larger_peaks=True, dmin=dmin)193 194 # Clean up and remove these peaks195 new_cell_peaks = []196 removed_cell_peaks = []197 removed_cell_peaks_larger = []198 for path_idx in range(len(cell_peaks)):199 path_retained_peaks = []200 path_removed_peaks = []201 path_peaks = cell_peaks[path_idx][0]202 203 for peak_idx in range(len(path_peaks)):204 if (path_idx, peak_idx) not in removed_peaks:205 path_retained_peaks.append(path_peaks[peak_idx])206 else:207 # What's the larger point?208 idx = removed_peaks.index((path_idx, peak_idx))209 larger_path, larger_idx = removed_larger_peaks[idx] 210 path_removed_peaks.append(RemovedPeakData(idx=path_peaks[peak_idx], screening_peak=(larger_path, cell_peaks[larger_path][0][larger_idx])))211 ###212 213 new_cell_peaks.append(path_retained_peaks)214 removed_cell_peaks.append(path_removed_peaks)215 216 cell_peaks = new_cell_peaks217 pd_list = []218 219 # Save peak positions, absolute intensity intensities, and length for each SC220 for k in range(len(all_paths)):221 222 points, o_intensity = all_paths[k], measured_trace_fluorescence[k]223 224 peaks = cell_peaks[k]225 removed_peaks = removed_cell_peaks[k]226 227 pd = PathData(peaks=peaks, removed_peaks=removed_peaks, points=points, o_intensity=o_intensity, SC_length=path_lengths[k])228 pd_list.append(pd)229 230 cd = CellData(pathdata_list=pd_list)231 232 return cd233 234 235alpha_max = 0.4236 237 238# Criterion used for identifying peak as a focus - normalized (with mean and s.d.)239# intensity levels being above 0.4 time maximum peak level240def focus_criterion(pos, v, alpha=alpha_max):241 """242 Identify and return positions where values in the array `v` exceed a certain threshold.243 244 The threshold is computed as `alpha` times the maximum value in `v`.245 246 Args:247 - pos (numpy.ndarray): Array of positions.248 - v (numpy.ndarray): 1D array of values, e.g., intensities.249 - alpha (float, optional): A scaling factor for the threshold. Defaults to `alpha_max`.250 251 Returns:252 - numpy.ndarray: Array of positions where corresponding values in `v` exceed the threshold.253 """254 if len(v):255 idx = (v>=alpha*np.max(v))256 return np.array(pos[idx])257 else:258 return np.array([], dtype=np.int32)259 260def analyse_celldata(cell_data, config):261 """262 Analyse the provided cell data to extract focus-related information.263 264 Args:265 cd (CellData): An instance of the CellData class containing path data information.266 config (dictionary): Configuration dictionary containing 'peak_threshold' and 'threshold_type'267 'peak_threshold' (float) - threshold for calling peaks as foci268 'threshold_type' (str) = 'per-trace', 'per-foci'269 270 Returns:271 tuple: A tuple containing:272 - foci_rel_intensity (list): List of relative intensities for the detected foci.273 - foci_pos (list): List of absolute positions of the detected foci.274 - foci_pos_index (list): List of indices of the detected foci.275 - screened_foci_data (list): List of RemovedPeakData indicating positions of removed peaks and the index of the larger peak276 - trace_median_intensities (list): Per-trace median intensity277 - trace_thresholds (list): Per-trace absolute threshold for calling peaks as foci278 """279 foci_abs_intensity = []280 foci_pos = []281 foci_pos_index = []282 screened_foci_data = []283 trace_median_intensities = []284 trace_thresholds = []285 286 peak_threshold = config['peak_threshold']287 288 threshold_type = config['threshold_type']289 290 if threshold_type == 'per-trace':291 """292 Call extracted peaks as foci if intensity - trace_mean > peak_threshold * (trace_max_foci_intensity - trace_mean)293 """294 295 for path_data in cell_data.pathdata_list:296 peaks = np.array(path_data.peaks, dtype=np.int32)297 298 # Normalize extracted fluorescent intensities by subtracting mean (and dividing299 # by standard deviation - note that the latter should have no effect on the results).300 h = np.array(path_data.o_intensity)301 h = h - np.mean(h)302 h = h/np.std(h)303 # Extract foci according to criterion304 foci_idx = focus_criterion(peaks, h[peaks], peak_threshold) 305 306 #307 removed_peaks = path_data.removed_peaks308 removed_peaks_idx = np.array([u.idx for u in removed_peaks], dtype=np.int32)309 310 311 if len(peaks):312 trace_thresholds.append((1-peak_threshold)*np.mean(path_data.o_intensity) + peak_threshold*np.max(np.array(path_data.o_intensity)[peaks]))313 else:314 trace_thresholds.append(None)315 316 if len(removed_peaks):317 if len(peaks):318 threshold = (1-peak_threshold)*np.mean(path_data.o_intensity) + peak_threshold*np.max(np.array(path_data.o_intensity)[peaks])319 else:320 threshold = float('-inf')321 322 323 removed_peak_heights = np.array(path_data.o_intensity)[removed_peaks_idx]324 screened_foci_idx = np.where(removed_peak_heights>threshold)[0]325 326 screened_foci_data.append([removed_peaks[i] for i in screened_foci_idx])327 else:328 screened_foci_data.append([])329 330 pos_abs = (foci_idx/len(path_data.points))*path_data.SC_length331 foci_pos.append(pos_abs)332 foci_abs_intensity.append(np.array(path_data.o_intensity)[foci_idx])333 334 foci_pos_index.append(foci_idx)335 trace_median_intensities.append(np.median(path_data.o_intensity))336 337 elif threshold_type == 'per-cell':338 """339 Call extracted peaks as foci if intensity - trace_mean > peak_threshold * max(intensity - trace_mean)340 """341 max_cell_intensity = float("-inf")342 for path_data in cell_data.pathdata_list:343 344 # Normalize extracted fluorescent intensities by subtracting mean (and dividing345 # by standard deviation - note that the latter should have no effect on the results).346 h = np.array(path_data.o_intensity)347 h = h - np.mean(h)348 max_cell_intensity = max(max_cell_intensity, np.max(h))349 350 for path_data in cell_data.pathdata_list:351 peaks = np.array(path_data.peaks, dtype=np.int32)352 353 # Normalize extracted fluorescent intensities by subtracting mean (and dividing354 # by standard deviation - note that the latter should have no effect on the results).355 h = np.array(path_data.o_intensity)356 h = h - np.mean(h)357 358 foci_idx = peaks[h[peaks]>peak_threshold*max_cell_intensity]359 360 removed_peaks = path_data.removed_peaks361 removed_peaks_idx = np.array([u.idx for u in removed_peaks], dtype=np.int32)362 363 trace_thresholds.append(np.mean(path_data.o_intensity) + peak_threshold*max_cell_intensity)364 365 if len(removed_peaks):366 threshold = np.mean(path_data.o_intensity) + peak_threshold*max_cell_intensity367 368 removed_peak_heights = np.array(path_data.o_intensity)[removed_peaks_idx]369 screened_foci_idx = np.where(removed_peak_heights>threshold)[0]370 371 screened_foci_data.append([removed_peaks[i] for i in screened_foci_idx])372 else:373 screened_foci_data.append([])374 375 pos_abs = (foci_idx/len(path_data.points))*path_data.SC_length376 foci_pos.append(pos_abs)377 foci_abs_intensity.append(np.array(path_data.o_intensity)[foci_idx])378 379 foci_pos_index.append(foci_idx)380 trace_median_intensities.append(np.median(path_data.o_intensity)) 381 382 else:383 raise NotImplementedError384 385 return foci_abs_intensity, foci_pos, foci_pos_index, screened_foci_data, trace_median_intensities, trace_thresholds386 387def analyse_traces(all_paths, path_lengths, measured_trace_fluorescence, config):388 389 cd = process_cell_traces(all_paths, path_lengths, measured_trace_fluorescence, dmin=config['screening_distance'])390 391 return analyse_celldata(cd, config)392 393 394 395 396