JFoz/test_path_analysis
0
1 2import lxml.etree as ET3import gzip4import tifffile5import matplotlib.pyplot as plt6import numpy as np7from PIL import Image, ImageDraw8import pandas as pd9from itertools import cycle10from .data_preprocess import analyse_traces11import math12import scipy.linalg as la13 14 15def get_paths_from_traces_file(traces_file):16 """17 Parses the specified traces file and extracts paths and their lengths.18 19 Args:20 traces_file (str): Path to the XML traces file.21 22 Returns:23 tuple: A tuple containing a list of paths (each path is a list of tuples representing points)24 and a list of corresponding path lengths.25 """26 tree = ET.parse(traces_file)27 root = tree.getroot()28 all_paths = []29 path_lengths = []30 for path in root.findall('path'):31 length=path.get('reallength')32 path_points = []33 for point in path:34 path_points.append((int(point.get('x')), int(point.get('y')), int(point.get('z'))))35 all_paths.append(path_points)36 path_lengths.append(float(length))37 return all_paths, path_lengths38 39 40 41def calculate_path_length_partials(point_list, voxel_size=(1,1,1)):42 """43 Calculate the partial path length of a series of points.44 45 Args:46 point_list (list of tuple): List of points, each represented as a tuple of coordinates (x, y, z).47 voxel_size (tuple, optional): Size of the voxel in each dimension (x, y, z). Defaults to (1, 1, 1).48 49 Returns:50 numpy.ndarray: Array of cumulative partial path lengths at each point.51 """52 # Simple calculation53 section_lengths = [0.0]54 s = np.array(voxel_size)55 for i in range(len(point_list)-1):56 # Euclidean distance between successive points57 section_lengths.append(la.norm(s * (np.array(point_list[i+1]) - np.array(point_list[i]))))58 return np.cumsum(section_lengths)59 60 61def visualise_ordering(points_list, dim, wr=5, wc=5):62 """63 Visualize the ordering of points in an image.64 65 Args:66 points_list (list): List of points to be visualized.67 dim (tuple): Dimensions of the image (rows, columns, channels).68 wr (int, optional): Width of the region to visualize around the point in the row direction. Defaults to 5.69 wc (int, optional): Width of the region to visualize around the point in the column direction. Defaults to 5.70 71 Returns:72 np.array: An image array with visualized points.73 """74 # Visualizes the ordering of the points in the list on a blank image.75 rdim, cdim, _ = dim76 vis = np.zeros((rdim, cdim, 3), dtype=np.uint8)77 78 def get_col(i):79 r = int(255 * i/len(points_list))80 g = 255 - r81 return r, g, 082 83 for n, p in enumerate(points_list):84 c, r, _ = map(int, p)85 vis[max(0,r-wr):min(rdim,r+wr+1),max(0,c-wc):min(cdim,c+wc+1)] = get_col(n)86 87 return vis88 89# A color map for paths90col_map = [(255,0,0), (0,255,0), (0,0,255), (255,255,0), (255,0,255), (0,255,255),91 (255,127,0), (255, 0, 127), (127, 255, 0), (0, 255, 127), (127,0,255), (0,127,255)]92 93def draw_paths(all_paths, foci_stack, foci_index=None, r=3, screened_foci_data=None):94 """95 Draws paths on the provided image stack and overlays markers for the foci96 97 Args:98 all_paths (list): List of paths where each path is a list of points.99 foci_stack (np.array): 3D numpy array representing the image stack.100 foci_index (list, optional): List of list of focus indices (along each path). Defaults to None.101 r (int, optional): Radius for the ellipse or line drawing around the focus. Defaults to 3.102 screened_foci_data (list, optional): List of RemovedPeakData for screened foci103 Returns:104 PIL.Image.Image: An image with the drawn paths.105 """106 im = np.max(foci_stack, axis=0)107 im = (im/np.max(im)*255).astype(np.uint8)108 im = np.dstack((im,)*3)109 im = Image.fromarray(im) 110 draw = ImageDraw.Draw(im)111 for i, (p, col) in enumerate(zip(all_paths, cycle(col_map))):112 draw.line([(u[0], u[1]) for u in p], fill=col)113 draw.text((p[0][0], p[0][1]), str(i+1), fill=col)114 115 if screened_foci_data is not None:116 for i, removed_peaks in enumerate(screened_foci_data):117 for p in removed_peaks:118 u = all_paths[i][p.idx]119 v = all_paths[p.screening_peak[0]][p.screening_peak[1]]120 draw.line((int(u[0]), int(u[1]), int(v[0]), int(v[1])), fill=(127,127,127), width=2)121 122 if foci_index is not None:123 for i, (idx, p, col) in enumerate(zip(foci_index, all_paths, cycle(col_map))):124 if len(idx):125 for j in idx:126 draw.line((int(p[j][0]-r), int(p[j][1]), int(p[j][0]+r), int(p[j][1])), fill=col, width=2)127 draw.line((int(p[j][0]), int(p[j][1]-r), int(p[j][0]), int(p[j][1]+r)), fill=col, width=2)128 return im129 130 131def measure_from_mask(mask, measure_stack):132 """133 Compute the sum of measure_stack values where the mask is equal to 1.134 135 Args:136 mask (numpy.ndarray): Binary mask where the measurement should be applied.137 measure_stack (numpy.ndarray): Stack of measurements.138 139 Returns:140 measure_stack.dtype: Sum of measure_stack values where the mask is 1.141 """142 return np.sum(mask * measure_stack)143 144# Max of measure_stack over region where mask==1145def max_from_mask(mask, measure_stack):146 """147 Compute the maximum of measure_stack values where the mask is equal to 1.148 149 Args:150 mask (numpy.ndarray): Binary mask where the measurement should be applied.151 measure_stack (numpy.ndarray): Stack of measurements.152 153 Returns:154 measure_stack.dtype: Maximum value of measure_stack where the mask is 1.155 """156 return np.max(mask * measure_stack)157 158def make_mask_s(p, melem, measure_stack):159 """160 Translate a mask to point p, ensuring correct treatment near the edges of the measure_stack.161 162 Args:163 p (tuple): Target point (r, c, z).164 melem (numpy.ndarray): Structuring element for the mask.165 measure_stack (numpy.ndarray): Stack of measurements.166 167 Returns:168 tuple: A tuple containing the translated mask and a section of the measure_stack.169 """170 171 172 # 173 174 R = [u//2 for u in melem.shape]175 176 r, c, z = p177 178 mask = np.zeros(melem.shape)179 180 m_data = np.zeros(melem.shape)181 s = measure_stack.shape182 o_1, o_2, o_3 = max(R[0]-r, 0), max(R[1]-c, 0), max(R[2]-z,0)183 e_1, e_2, e_3 = min(R[0]-r+s[0], 2*R[0]+1), min(R[1]-c+s[1], 2*R[1]+1), min(R[2]-z+s[2], 2*R[2]+1)184 m_data[o_1:e_1,o_2:e_2,o_3:e_3] = measure_stack[max(r-R[0],0):min(r+R[0]+1,s[0]),max(c-R[1],0):min(c+R[1]+1,s[1]),max(z-R[2],0):min(z+R[2]+1, s[2])]185 mask[o_1:e_1,o_2:e_2,o_3:e_3] = melem[o_1:e_1,o_2:e_2,o_3:e_3] 186 187 188 return mask, m_data189 190 191def measure_at_point(p, melem, measure_stack, op='mean'):192 """193 Measure the mean or max value of measure_stack around a specific point using a structuring element.194 195 Args:196 p (tuple): Target point (r, c, z).197 melem (numpy.ndarray): Structuring element for the mask.198 measure_stack (numpy.ndarray): Stack of measurements.199 op (str, optional): Operation to be applied; either 'mean' or 'max'. Default is 'mean'.200 201 Returns:202 float: Measured value based on the specified operation.203 """204 205 p = map(int, p)206 if op=='mean':207 mask, m_data = make_mask_s(p, melem, measure_stack)208 melem_size = np.sum(mask)209 return float(measure_from_mask(mask, m_data) / melem_size)210 else:211 mask, m_data = make_mask_s(p, melem, measure_stack)212 return float(max_from_mask(mask, m_data))213 214# Generate spherical region215def make_sphere(R=5, z_scale_ratio=2.3):216 """217 Generate a binary representation of a sphere in 3D space.218 219 Args:220 R (int, optional): Radius of the sphere. Default is 5. Centred on the centre of the middle voxel.221 Includes all voxels whose centre is precisely R from the middle voxel.222 z_scale_ratio (float, optional): Scaling factor for the z-axis. Default is 2.3.223 224 Returns:225 numpy.ndarray: Binary representation of the sphere.226 """227 R_z = int(math.ceil(R/z_scale_ratio))228 x, y, z = np.ogrid[-R:R+1, -R:R+1, -R_z:R_z+1]229 sphere = x**2 + y**2 + (z_scale_ratio * z)**2 <= R**2230 return sphere231 232# Measure the values of measure_stack at each of the points of points_list in turn.233# Measurement is the mean / max (specified by op) on the spherical region about each point234def measure_all_with_sphere(points_list, measure_stack, op='mean', R=5, z_scale_ratio=2.3):235 """236 Measure the values of measure_stack at each point in a list using a spherical region.237 238 Args:239 points_list (list): List of points (r, c, z) to be measured.240 measure_stack (numpy.ndarray): Stack of measurements.241 op (str, optional): Operation to be applied; either 'mean' or 'max'. Default is 'mean'.242 R (int, optional): Radius of the sphere. Default is 5.243 z_scale_ratio (float, optional): Scaling factor for the z-axis. Default is 2.3.244 245 Returns:246 list: List of measured values for each point.247 """248 melem = make_sphere(R, z_scale_ratio)249 measure_func = lambda p: measure_at_point(p, melem, measure_stack, op)250 return list(map(measure_func, points_list))251 252 253# Measure fluorescence levels along ordered skeleton254def measure_chrom2(path, intensity, config):255 """256 Measure fluorescence levels along an ordered skeleton.257 258 Args:259 path (list): List of ordered path points (r, c, z).260 intensity (numpy.ndarray): 3D fluorescence data.261 config (dict): Configuration dictionary containing 'z_res', 'xy_res', and 'sphere_radius' values.262 263 Returns:264 tuple: A tuple containing the visualization, mean measurements, and max measurements along the path.265 """266 # Calculate size of spheroid used for measurement267 scale_ratio = config['z_res']/config['xy_res']268 sphere_xy_radius = int(math.ceil(config['sphere_radius']/config['xy_res']))269 270 vis = visualise_ordering(path, dim=intensity.shape, wr=sphere_xy_radius, wc=sphere_xy_radius)271 272 measurements = measure_all_with_sphere(path, intensity, op='mean', R=sphere_xy_radius, z_scale_ratio=scale_ratio)273 measurements_max = measure_all_with_sphere(path, intensity, op='max', R=sphere_xy_radius, z_scale_ratio=scale_ratio)274 275 276 return vis, measurements, measurements_max277 278def extract_peaks(cell_id, all_paths, path_lengths, measured_traces, config):279 """280 Extract peak information from given traces and compile them into a DataFrame.281 282 Args:283 - cell_id (int or str): Identifier for the cell being analyzed.284 - all_paths (list of lists): Contains ordered path points for multiple paths.285 - path_lengths (list of floats): List containing lengths of each path in all_paths.286 - measured_traces (list of lists): Contains fluorescence measurement values along the paths.287 - config (dict): Configuration dictionary containing:288 - 'peak_threshold': Threshold value to determine a peak in the trace.289 - 'sphere_radius': Radius of the sphere used in fluorescence measurement.290 291 Returns:292 - pd.DataFrame: DataFrame containing peak information for each path.293 - list of lists: Absolute intensities of the detected foci.294 - list of lists: Index positions of the detected foci.295 - list of lists: Absolute focus intensity threshold for each trace.296 - list of numpy.ndarray: For each trace, distances of each point from start of trace in microns297 """298 299 n_paths = len(all_paths)300 301 data = []302 foci_absolute_intensity, foci_position, foci_position_index, screened_foci_data, trace_median_intensities, trace_thresholds = analyse_traces(all_paths, path_lengths, measured_traces, config)303 304 # Normalize foci intensities (for quantification) using trace medians as estimates of background305 foci_intensities = []306 for path_foci_abs_int, tmi in zip(foci_absolute_intensity, trace_median_intensities):307 foci_intensities.extend(list(path_foci_abs_int - tmi))308 309 # Divide all foci intensities by the mean within the cell310 mean_intensity = np.mean(foci_intensities)311 trace_positions = []312 313 for i in range(n_paths):314 315 # Calculate real (Euclidean) distance of each point along the traced path316 pl = calculate_path_length_partials(all_paths[i], (config['xy_res'], config['xy_res'], config['z_res']))317 318 319 path_data = { 'Cell_ID':cell_id,320 'Trace': i+1,321 'SNT_trace_length(um)': path_lengths[i],322 'Measured_trace_length(um)': pl[-1],323 'Trace_median_intensity': trace_median_intensities[i],324 'Detection_sphere_radius(um)': config['sphere_radius'],325 'Screening_distance(voxels)': config['screening_distance'],326 'Foci_ID_threshold': config['peak_threshold'],327 'Trace_foci_number': len(foci_position_index[i]) }328 for j, (idx, u,v) in enumerate(zip(foci_position_index[i], foci_position[i], foci_absolute_intensity[i])):329 if config['use_corrected_positions']:330 # Use the calculated position along the traced path331 path_data[f'Foci_{j+1}_position(um)'] = pl[idx]332 else:333 # Use the measured trace length (from SNT), and assume all steps of path are approximately the same length334 path_data[f'Foci_{j+1}_position(um)'] = u335 # The original measured intensity (mean in spheroid around detected peak)336 path_data[f'Foci_{j+1}_absolute_intensity'] = v337 # Measure relative intensity by removing per-trace background and dividing by cell total338 path_data[f'Foci_{j+1}_relative_intensity'] = (v - trace_median_intensities[i])/mean_intensity339 data.append(path_data)340 trace_positions.append(pl)341 return pd.DataFrame(data), foci_absolute_intensity, foci_position_index, screened_foci_data, trace_thresholds, trace_positions342 343 344def analyse_paths(cell_id,345 foci_file,346 traces_file,347 config348 ):349 """350 Analyzes paths for the given cell ID using provided foci and trace files.351 352 Args:353 cell_id (int/str): Identifier for the cell.354 foci_file (str): Path to the foci image file.355 traces_file (str): Path to the XML traces file.356 config (dict): Configuration dictionary containing necessary parameters such as resolutions and thresholds.357 358 Returns:359 tuple: A tuple containing an overlay image of the traces, visualization images for each trace,360 a figure with plotted measurements, and a dataframe with extracted peaks.361 """362 363 364 # Read stack365 366 foci_stack = tifffile.imread(foci_file)367 368 # If 2D add additional (z) dimension369 if foci_stack.ndim==2:370 foci_stack = foci_stack[None,:,:]371 372 all_paths, path_lengths = get_paths_from_traces_file(traces_file)373 374 all_trace_vis = [] # Per-path visualizations375 all_m = [] # Per-path measured intensities376 for p in all_paths:377 # Measure intensity along path - transpose the stack ZYX -> XYZ378 vis, m, _ = measure_chrom2(p,foci_stack.transpose(2,1,0), config)379 all_trace_vis.append(vis)380 all_m.append(m)381 382 383 # Extract all data from paths and traces384 extracted_peaks, foci_absolute_intensity, foci_pos_index, screened_foci_data, trace_thresholds, trace_positions = extract_peaks(cell_id, all_paths, path_lengths, all_m, config)385 386 # Plot per-path measured intensities and indicate foci387 n_cols = 2388 n_rows = (len(all_paths)+n_cols-1)//n_cols389 fig, ax = plt.subplots(n_rows,n_cols, figsize=(5*n_cols, 3*n_rows))390 ax = ax.flatten()391 392 for i, m in enumerate(all_m):393 ax[i].set_title(f'Trace {i+1}')394 ax[i].plot(trace_positions[i], m)395 if len(foci_pos_index[i]):396 # Plot detected foci397 ax[i].plot(trace_positions[i][foci_pos_index[i]], np.array(m)[foci_pos_index[i]], 'rx')398 399 if len(screened_foci_data[i]):400 # Indicate screened foci by gray circles on plots401 screened_foci_pos_index = [u.idx for u in screened_foci_data[i]]402 ax[i].plot(trace_positions[i][screened_foci_pos_index], np.array(m)[screened_foci_pos_index], color=(0.5,0.5,0.5), marker='o', linestyle='None')403 404 # Show per-trace intensity thresholds with red dotted lines405 if trace_thresholds[i] is not None:406 ax[i].axhline(trace_thresholds[i], c='r', ls=':')407 ax[i].set_xlabel('Distance from start (um)')408 ax[i].set_ylabel('Intensity')409 # Hide excess plots410 for i in range(len(all_m), n_cols*n_rows):411 ax[i].axis('off')412 413 plt.tight_layout()414 trace_overlay = draw_paths(all_paths, foci_stack, foci_index=foci_pos_index, screened_foci_data=screened_foci_data)415 416 return trace_overlay, all_trace_vis, fig, extracted_peaks417