OSCR

Brain PDGFRβ+ Cells Exhibit Diverse Reactive Phenotypes After Stroke Without Requiring KLF4

Code ↔ Paper

21 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 21 matches
  1. [1] § Materials and Methods › Topological Data Analysis (TDA) ↔ Widefield_20x_ROIs_Pdgfrb_TDA.ipynb, lines 314 to 395 · score 0.80 · PointClouds, persistence diagrams, persistent homology, topological features, Ripser, dimension · Report a mismatch
  2. [2] § Results › Morphologically Distinct PDGFRβ + Cells Are Present at the Injured Tissue ↔ FlowCytometry_Pdgfrb_Processing.qmd, lines 477 to 562 · score 0.79 · Live Dead, geom_density2d, geom_stats, geom_gate, FACS_Pdgfrb, dead cells · Report a mismatch
  3. [3] § Results › Morphologically Distinct PDGFRβ + Cells Are Present at the Injured Tissue ↔ FlowCytometry_Pdgfrb_Processing.qmd, lines 476 to 561 · score 0.79 · Live Dead, geom_density2d, geom_stats, geom_gate, FACS_Pdgfrb, dead cells · Report a mismatch
  4. [4] § Materials and Methods › Fluorescence‐Activated Cell Sorting (FACS) ↔ FlowCytometry_Pdgfrb_Processing.qmd, lines 48 to 92 · score 0.75 · flowAI, flowCore, compensation controls, OSF repository, ggcyto, Bioconductor · Report a mismatch
  5. [5] § Materials and Methods › Fluorescence‐Activated Cell Sorting (FACS) ↔ FlowCytometry_Pdgfrb_Processing.qmd, lines 48 to 92 · score 0.75 · flowAI, flowCore, compensation controls, OSF repository, ggcyto, Bioconductor · Report a mismatch
  6. [6] § Materials and Methods › Data Analysis ↔ Widefield_5x_Ipsilateral_Gfap_Convex.qmd, lines 234 to 342 · score 0.73 · adapt_delta, max_treedepth, parameter space, Bayesian modeling, brms package, Loo · Report a mismatch
  7. [7] § Materials and Methods › Data Analysis ↔ Widefield_5x_Ipsilateral_Gfap_Convex.qmd, lines 234 to 342 · score 0.73 · adapt_delta, max_treedepth, parameter space, Bayesian modeling, brms package, Loo · Report a mismatch
  8. [8] § Results › Spatial and Topological Organization of Reactive PDGFRβ + Cells Upon Injury ↔ Widefield_20x_ROIs_Pdgfrb_TDA.ipynb, lines 205 to 253 · score 0.69 · animated vietoris rips, vietoris rips complexes, pdgfrb tda, widefield 20x, point clouds, filtration · Report a mismatch
  9. [9] § Results › KLF4 Induction in Reactive PDGFRβ + Cells Is Minimal Upon Ischemic Stroke ↔ Widefield_10x_Ipsilateral_KO_Klf4_Exp.qmd, lines 654 to 688 · score 0.66 · sigma Genotype, spread_draws, Probability density, stat_halfeye, KLF4 expression, KLF4 KO · Report a mismatch
  10. [10] § Results › KLF4 Induction in Reactive PDGFRβ + Cells Is Minimal Upon Ischemic Stroke ↔ Widefield_10x_Ipsilateral_KO_Klf4_Exp.qmd, lines 653 to 686 · score 0.66 · sigma Genotype, spread_draws, Probability density, stat_halfeye, KLF4 expression, KLF4 KO · Report a mismatch
  11. [11] § Results › Spatial and Topological Organization of Reactive PDGFRβ + Cells Upon Injury ↔ Widefield_20x_ROIs_Pdgfrb_TDA.ipynb, lines 700 to 791 · score 0.61 · dimensional homology, median distances, bottleneck distances, peri, CI, striatum · Report a mismatch
  12. [12] § Results › Reactive PDGFRβ + Cells Comprise Morphologically Diverse Populations ↔ Confocal_40x_ROIs_Pdgfrb_Morpho_Analysis.qmd, lines 317 to 349 · score 0.59 · Cartesian plane, opposite sides, negatively correlated, Cos2, Biplot, PCA · Report a mismatch
  13. [13] § Results › Reactive PDGFRβ + Cells Comprise Morphologically Diverse Populations ↔ Confocal_40x_ROIs_Pdgfrb_Morpho_Analysis.qmd, lines 317 to 349 · score 0.59 · Cartesian plane, opposite sides, negatively correlated, Cos2, Biplot, PCA · Report a mismatch
  14. [14] § Materials and Methods › Data Analysis ↔ Widefield_5x_Ipsilateral_Gfap_Convex.qmd, lines 234 to 342 · score 0.59 · Markov chain Monte, model comparison, brms package, Carlo, MCMC, statistical modeling · Report a mismatch
  15. [15] § Materials and Methods › Data Analysis ↔ Widefield_5x_Ipsilateral_Gfap_Convex.qmd, lines 234 to 342 · score 0.59 · Markov chain Monte, model comparison, brms package, Carlo, MCMC, statistical modeling · Report a mismatch
  16. [16] § Materials and Methods › Point Pattern Analysis (PPA) ↔ Widefield_10x_Ipsilateral_Ki67_PPA.qmd, lines 206 to 276 · score 0.57 · PointPatterns, spatstat package, hyperframes, ppp, PPA, rds · Report a mismatch
  17. [17] § Materials and Methods › Point Pattern Analysis (PPA) ↔ Widefield_10x_Ipsilateral_KO_Klf4_Exp.qmd, lines 198 to 267 · score 0.56 · PointPatterns, spatstat package, hyperframes, ppp, PPA, rds · Report a mismatch
  18. [18] § Results › Spatial and Topological Organization of Reactive PDGFRβ + Cells Upon Injury ↔ Widefield_20x_ROIs_Pdgfrb_Haralick.qmd, lines 651 to 729 · score 0.55 · local variation, inverse, Haralick, IDM, entropy, Ctx · Report a mismatch
  19. [19] § Results › Spatial and Topological Organization of Reactive PDGFRβ + Cells Upon Injury ↔ Widefield_20x_ROIs_Pdgfrb_Haralick.qmd, lines 651 to 729 · score 0.55 · local variation, inverse, Haralick, IDM, entropy, Ctx · Report a mismatch
  20. [20] § Materials and Methods ↔ Widefield_10x_ROIs_CD31-Pdgfrb_Coloc.qmd, lines 195 to 236 · score 0.53 · middle cerebral artery, sham animals, post ischemic, tamoxifen, occlusion, MCAo · Report a mismatch
  21. [21] § Results › Spatial and Topological Organization of Reactive PDGFRβ + Cells Upon Injury ↔ Widefield_10x_ROIs_KO_CD31-ColIV_Coloc.qmd, lines 414 to 493 · score 0.51 · geom_point, add_predicted_draws, stat_halfeye, CD31, sigma, KLF4 · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Jupyter notebook · 1,065 lines · 45 KB · no license · 3 matches

  1. # %% [markdown]
  2. # # Topological Data Analysis (TDA) of PDGFR-β+ cells
  4. # %% [markdown]
  5. # ## Preview
  6. #
  7. # This notebook reports the topological data analysis (TDA) for PDGFR-β cells in specific brain regions (injured striatum and cortex, and perilesional tissue). Please note that this notebook runs Python code.
  8. #
  9. # **Parent dataset:** PDGFR-β(Td-tomato) ROIs imaged at 20x in different Z-planes. Images were z-projected for analysis. Samples are grouped at 0 (Sham), 3, 7, 14, and 30 days post-ischemia (DPI). The raw images and pre-processing scripts (if applicable) are available at the Zenodo repository (10.5281/zenodo.10553084) under the name `Widefield_20x_ROIs_Pdgfrb.zip`.
  10. #
  11. # **Working dataset**: The `Raw_Data/Widefield_20x_ROIs_Pdgfrb/CellProfiler/Cells.csv`data frame containing the individual xy coordinates of PDGFR-β+ cells generated with CellProfiler [@stirling2021].
  12. #
  13. # **Note:** This notebook runs using python code.
  15. # %% [markdown]
  16. # ## Load required libraries
  17. #
  18. # We import the required Python libraries.
  20. # %%
  21. import pandas as pd
  22. import numpy as np
  23. import seaborn as sns
  24. import os
  25. import csv
  27. import matplotlib.pyplot as plt
  28. from matplotlib.collections import LineCollection
  29. from matplotlib.animation import FuncAnimation
  30. from matplotlib.patches import Circle
  32. import gudhi as gd
  33. from gudhi import RipsComplex, plot_persistence_diagram
  34. from gudhi import bottleneck_distance
  35. from gudhi.wasserstein import wasserstein_distance
  37. import itertools
  38. from ripser import ripser
  39. from persim import plot_diagrams
  40. from random import choices
  41. from scipy.spatial import ConvexHull
  43. # %% [markdown]
  44. # ## Load and handle the dataset
  45. #
  46. # We load the `Raw_Data/Widefield_20x_ROIs_Pdgfrb/CellProfiler/Cells.csv` dataset containing the cell coordinates. We handle the data to obtain and organize the relevant metadata from the file name.
  48. # %%
  49. # Load the dataset
  50. file_path = 'Data_Raw/Widefield_20x_ROIs_Pdgfrb/CellProfiler/Cells.csv'
  51. data = pd.read_csv(file_path)
  53. # Keep only the specified columns
  54. data = data[['Location_Center_X', 'Location_Center_Y', 'FileName_Raw']]
  56. # Define a function to parse the FileName_Raw column
  57. def parse_filename(filename):
  58. parts = filename.split('_')
  59. # Extract the needed parts based on their position
  60. MouseID = parts[0]
  61. DPI = parts[1]
  62. Condition = parts[2]
  63. Lesion = parts[3]
  64. Region = parts[4]
  65. return MouseID, DPI, Condition, Lesion, Region
  67. # Apply the function to the DataFrame
  68. data[['MouseID', 'DPI', 'Condition', 'Lesion', 'Region']] = data.apply(lambda row: parse_filename(row['FileName_Raw']), axis=1, result_type='expand')
  70. # Drop the FileName_Raw column
  71. data.drop('FileName_Raw', axis=1, inplace=True)
  73. # Reorder columns to move Location_Center_X and Location_Center_Y to the right
  74. column_order = ['MouseID', 'DPI', 'Condition', 'Lesion', 'Region', 'Location_Center_X', 'Location_Center_Y']
  75. data = data[column_order]
  77. # Save the updated dataset to a new file
  78. output_file_path = 'Data_Processed/Widefield_20x_ROIs_Pdgfrb/Widefield_20x_ROIs_Pdgfrb_Coords.csv'
  79. data.to_csv(output_file_path, index=False)
  81. print("Dataset updated and saved to:", output_file_path)
  84. # %% [markdown]
  85. # ## Calculate Vietoris-Rips Complex
  86. #
  87. # In this section, we generate Vietoris-Rips complex for elucidating the topological features embedded within high-dimensional datasets. We construct the complexes from a finite set of points in a metric space to facilitate the analysis of their spatial relationships through a simplicial complex. The construction follows this scheme:
  88. #
  89. # Let $S$ be a finite set of points in a metric space. The Vietoris-Rips complex, denoted as $ VR_{\epsilon}(S) $, is defined for a chosen distance threshold $ \epsilon $. The construction process is as follows:
  90. #
  91. # 1. **Vertices**: Each point in $ S $ is considered a vertex, or 0-simplex.
  92. #
  93. # 2. **Edges**: An edge, or 1-simplex, is formed between any pair of vertices $ (u, v) $ if the distance between them, denoted $ d(u, v) $, satisfies the condition $ d(u, v) \leq \epsilon $.
  94. #
  95. # 3. **Higher-Dimensional Simplices**: A $ k $-simplex is formed by a set of $ k+1 $ vertices if every pair of vertices within this set is connected by an edge. Mathematically, for a subset $ \{v_0, v_1, \ldots, v_k\} \subseteq S $, a $ k $-simplex exists if:
  96. #
  97. # $$
  98. # \forall i, j \in \{0, 1, \ldots, k\},\, d(v_i, v_j) \leq \epsilon
  99. # $$
  100. #
  101. # 4. **Complex Assembly**: The Vietoris-Rips complex $ VR_{\epsilon}(S) $ comprises all simplices, from 0-simplices to the highest dimension formed under the given threshold $ \epsilon $.
  102. #
  103. # We load the cell coordinates to build the Vietoris-Rips complexes. We set `100` as the maximum filtration value to facilitate visualization.
  105. # %%
  106. # Load the organized dataset
  107. data_path = 'Data_Processed/Widefield_20x_ROIs_Pdgfrb/Widefield_20x_ROIs_Pdgfrb_Coords.csv' # Adjust this path as necessary
  108. data = pd.read_csv(data_path)
  110. # Group the dataset by MouseID and Region, then extract the point clouds
  111. point_clouds_by_mouse_region = data.groupby(['MouseID', 'DPI', 'Region']).apply(lambda df: df[['Location_Center_X', 'Location_Center_Y']].values).to_dict()
  114. def calculate_vietoris_rips_complex(point_cloud, max_dim, max_radius):
  115. """
  116. Calculate the Vietoris-Rips complex for a given point cloud.
  118. Parameters:
  119. - point_cloud: numpy array of shape (n_points, n_dimensions) representing the point cloud.
  120. - max_dim: int, the maximal dimension of the simplices in the Vietoris-Rips complex.
  121. - max_radius: float, the maximal radius for simplices in the Vietoris-Rips complex.
  123. Returns:
  124. - simplex_tree: a GUDHI SimplexTree object representing the Vietoris-Rips complex.
  125. """
  126. rips_complex = gd.RipsComplex(points=point_cloud, max_edge_length=max_radius)
  127. simplex_tree = rips_complex.create_simplex_tree(max_dimension=max_dim)
  128. return simplex_tree
  130. # Parameters for Vietoris-Rips complex calculation
  131. max_dim = 3
  132. max_radius = 100
  134. # Calculate Vietoris-Rips complexes for each region of each animal
  135. vietoris_rips_complexes = {}
  136. for (mouse_id, dpi, region), point_cloud in point_clouds_by_mouse_region.items():
  137. key = (mouse_id, dpi, region) # Use a tuple of MouseID and Region as the key
  138. vietoris_rips_complexes[key] = calculate_vietoris_rips_complex(point_cloud, max_dim, max_radius)
  140. # %% [markdown]
  141. # ### Visualize and save Vietoris-Rips Complex
  142. #
  143. # Once calculated, we save the images in a designated folder:
  145. # %%
  146. def plot_vietoris_rips_complex(point_cloud, simplex_tree, filename, save_dir, x_range=None, y_range=None):
  147. """
  148. Plot and save the Vietoris-Rips complex for a given 2D point cloud, with specified x and y axis ranges.
  150. Parameters:
  151. - point_cloud: numpy array of shape (n_points, 2) representing the 2D point cloud.
  152. - simplex_tree: a gudhi object representing the Vietoris-Rips complex.
  153. - filename: str, the name of the file to save the image.
  154. - save_dir: str, directory to save the image.
  155. - x_range: tuple, the range for the x-axis (optional).
  156. - y_range: tuple, the range for the y-axis (optional).
  157. """
  158. fig, ax = plt.subplots(figsize=(6, 4))
  160. # Plot points
  161. ax.scatter(point_cloud[:, 0], point_cloud[:, 1], c='black', s=20, alpha=0.8)
  163. # Plot edges
  164. edges = []
  165. for simplex in simplex_tree.get_skeleton(2): # Get up to 2-simplices (edges)
  166. if len(simplex[0]) == 2: # Edges
  167. p1, p2 = point_cloud[simplex[0][0]], point_cloud[simplex[0][1]]
  168. edges.append([p1, p2])
  169. lc = LineCollection(edges, colors='green', linewidths=1, alpha=0.7)
  170. ax.add_collection(lc)
  172. # Adjust plot settings
  173. ax.set_xticks([])
  174. ax.set_yticks([])
  175. ax.set_xlabel('')
  176. ax.set_ylabel('')
  178. if x_range is not None:
  179. ax.set_xlim(x_range)
  180. if y_range is not None:
  181. ax.set_ylim(y_range)
  183. plt.tight_layout()
  185. if not os.path.exists(save_dir):
  186. os.makedirs(save_dir)
  187. plt.savefig(os.path.join(save_dir, f"{filename}.png"))
  188. plt.close()
  191. # Directory to save the images
  192. output_dir = "D:/Research/Stroke_PDGFR-B_Reactivity/Images/Widefield_20x-ROIs_Pdgfrb_TDA/VietorisRips_Static" # Update this path
  193. x_range = (0, 2000) # Example x-axis range
  194. y_range = (0, 1500) # Example y-axis range
  195. aspect_ratio = 'auto' # Adjust this based on your data or preferences
  198. # Assuming vietoris_rips_complexes and point_clouds_by_mouse have been computed accordingly
  199. for (mouse_id, dpi, region), simplex_tree in vietoris_rips_complexes.items():
  200. filename = f"{mouse_id}_{dpi}_{region}_VietorisRips"
  201. point_cloud = point_clouds_by_mouse_region[(mouse_id, dpi, region)]
  202. plot_vietoris_rips_complex(point_cloud, simplex_tree, filename, output_dir, x_range, y_range)
  205. # %% [markdown]
  206. # ### Visualize and save Vietoris-Rips Complex (animated)
  207. #
  208. # We created animated Vietoris-rips for illustration in the scientific article. We save animations with a maximum of `300` filtration value.
  210. # %%
  211. def animate_vietoris_rips_mp4(point_cloud, simplex_tree, filename, save_dir, max_radius, step_size=1):
  212. """
  213. Create an .mp4 animation illustrating the growth of the Vietoris-Rips complex for 2D data.
  214. """
  215. fig, ax = plt.subplots(figsize=(7, 7))
  216. ax.scatter(point_cloud[:, 0], point_cloud[:, 1], c='black', s=20, alpha=0.8)
  218. # Function to update the plot for each frame
  219. def update(threshold):
  220. ax.cla() # Clear current plot details
  221. ax.scatter(point_cloud[:, 0], point_cloud[:, 1], c='black', s=10, alpha=0.7)
  222. edges = []
  223. for simplex in simplex_tree.get_skeleton(2):
  224. if simplex[1] <= threshold:
  225. if len(simplex[0]) == 2: # Edges
  226. p1, p2 = point_cloud[simplex[0][0]], point_cloud[simplex[0][1]]
  227. edges.append([p1, p2])
  228. lc = LineCollection(edges, colors='green', alpha=0.7)
  229. ax.add_collection(lc)
  231. ax.set_title(f"Threshold: {threshold:.2f}", fontsize=20)
  232. ax.set_xlabel('X', fontsize=20, labelpad=15)
  233. ax.set_ylabel('Y', fontsize=20, labelpad=15)
  234. ax.tick_params(axis='both', which='major', labelsize=20)
  235. ax.axis('equal')
  237. ani = FuncAnimation(fig, update, frames=np.arange(0, max_radius, step_size), blit=False)
  238. ani_file = os.path.join(save_dir, f"{filename}_Animation.mp4")
  239. ani.save(ani_file, writer='ffmpeg', fps=2, extra_args=['-vcodec', 'libx264'])
  241. plt.close()
  242. return ani_file # Return the file path of the saved animation
  244. # Example usage
  245. output_dir = "D:/Research/Stroke_PDGFR-B_Reactivity/Images/Widefield_20x-ROIs_Pdgfrb_TDA/VietorisRips_Animated"
  246. max_radius = 300 # Adjust as necessary
  247. step_size = 20 # Adjust as necessary
  249. for (mouse_id, dpi, region), simplex_tree in vietoris_rips_complexes.items():
  250. filename = f"{mouse_id}_{dpi}_{region}_VietorisRips"
  251. point_cloud = point_clouds_by_mouse_region[(mouse_id, dpi, region)]
  252. ani_file = animate_vietoris_rips_mp4(point_cloud, simplex_tree, filename, output_dir, max_radius, step_size)
  253. print(f"Saved animation: {ani_file}")
  255. # %% [markdown]
  256. # ## Generate point clouds
  257. #
  258. # We generate point clouds and store them as individual files to use along the TDA. First, we load again the dataset in case is not present in the environment.
  260. # %%
  261. # Load your dataset
  262. data_path = 'Data_Processed/Widefield_20x_ROIs_Pdgfrb/Widefield_20x_ROIs_Pdgfrb_Coords.csv'
  263. data = pd.read_csv(data_path)
  265. # Group data by MouseID, Region, and DPI
  266. point_clouds_by_group = data.groupby(['MouseID', 'DPI', 'Region',]).apply(
  267. lambda df: df[['Location_Center_X', 'Location_Center_Y']].values
  268. ).to_dict()
  270. # %% [markdown]
  271. # Next, we verify the structire of the data to make sure is grouped as intented.
  273. # %%
  274. if not point_clouds_by_group:
  275. print("The point_clouds_by_group dictionary is empty.")
  276. else:
  277. print("Point clouds loaded:", len(point_clouds_by_group))
  279. # Inspect the first few items to ensure they're structured correctly
  280. for key, value in list(point_clouds_by_group.items())[:5]:
  281. print(key, value.shape)
  283. # %% [markdown]
  284. # Then, we save point clouds as .npy files (Numpy array) in a designated folder
  286. # %%
  287. def save_point_clouds_as_npy(point_clouds_by_group, data_directory):
  288. """
  289. Saves point clouds to .npy files in the specified directory.
  291. Parameters:
  292. point_clouds_by_group (dict): Dictionary with keys as identifiers
  293. (e.g., ('MouseID', 'Region', 'DPI')) and values as numpy arrays representing point clouds.
  294. data_directory (str): Path to the directory where .npy files will be saved.
  295. """
  296. # Ensure the directory exists
  297. if not os.path.exists(data_directory):
  298. os.makedirs(data_directory)
  300. # Iterate through the point clouds and save each as an .npy file
  301. for identifier, point_cloud in point_clouds_by_group.items():
  302. # Create a filename from the identifier
  303. filename = "_".join(identifier) + ".npy"
  304. filepath = os.path.join(data_directory, filename)
  305. np.save(filepath, point_cloud)
  307. print("Point clouds saved as .npy files.")
  309. data_directory = r"PointClouds/Widefield_20x_ROIs_Pdgfrb/"
  311. # Call the function with your data and desired save location
  312. save_point_clouds_as_npy(point_clouds_by_group, data_directory)
  314. # %% [markdown]
  315. # ## Computation of persistent homology
  317. # %% [markdown]
  318. # In this section, we analize the generated point clouds using persisten homology. We use the same data the `process_and_visualize_persistence` function to process the same dataset containing xy coordinates of PDGFR-B+ cells. We compute persistent diagrams using the `ripster` libary. These diagrams capture the birth and death of topological features (like loops or voids) in the data as one varies a filtering parameter. The persistence diagrams are then visualized using the `persim` library.
  319. #
  320. # This analysis provide insights into the topological structure of the data, indicating features that persist across scales. In TDA, features that persist longer (i.e., have a wider gap between birth and death) are deemed more significant as they represent genuine structures rather than noise. Then, we save each persistence diagram as a .png image.
  322. # %% [markdown]
  323. # ### Compute persistence diagrams
  325. # %% [markdown]
  326. # We compute persistence diagrams and save then as .npy file for further analysis and as .png images for visualization
  328. # %%
  329. # Disable LaTeX usage in matplotlib
  330. plt.rcParams['text.usetex'] = False
  332. def process_npy_for_persistence_diagrams(data_directory, save_directory, save_diag_directory, img_size=(6, 6), x_range=None, y_range=None):
  333. """
  334. Batch process .npy files (containing point clouds) to compute, visualize, and save persistence diagrams.
  335. """
  336. # Ensure save directories exist
  337. if not os.path.exists(save_directory):
  338. os.makedirs(save_directory)
  339. if not os.path.exists(save_diag_directory):
  340. os.makedirs(save_diag_directory)
  342. # List all .npy files in the data directory
  343. npy_files = [f for f in os.listdir(data_directory) if f.endswith('.npy')]
  345. for file in npy_files:
  346. # Load the point cloud from the .npy file
  347. points_array = np.load(os.path.join(data_directory, file))
  349. # Compute the persistence diagrams using Ripser
  350. diagrams = ripser(points_array, maxdim=1)['dgms']
  352. # Save each dimension's persistence diagram as a separate .npy file
  353. for dim, diag in enumerate(diagrams):
  354. diag_filename = os.path.splitext(file)[0] + f"_persistence_diag_dim{dim}.npy"
  355. np.save(os.path.join(save_diag_directory, diag_filename), diag)
  357. # Visualize and save the persistence diagrams
  358. fig, ax = plt.subplots(figsize=img_size)
  359. plot_diagrams(diagrams, show=False, ax=ax)
  360. ax.set_title("")
  361. ax.set_xlabel('Birth', fontsize=25)
  362. ax.set_ylabel('Death', fontsize=25)
  363. ax.tick_params(axis='both', which='major', labelsize=25)
  365. # Optionally adjust the legend
  366. legend = ax.legend()
  367. if legend:
  368. for label in legend.get_texts():
  369. label.set_fontsize(20)
  371. plt.tight_layout()
  373. if x_range is not None:
  374. ax.set_xlim(x_range)
  375. if y_range is not None:
  376. ax.set_ylim(y_range)
  378. # Save the figure
  379. filename = os.path.splitext(file)[0] + "_persistence_diagram.png"
  380. fig.savefig(os.path.join(save_directory, filename))
  381. plt.close()
  383. print(f"Persistence diagrams saved to {save_directory}.")
  384. print(f"Persistence diagrams data saved to {save_diag_directory}.")
  386. # Paths to your directories
  387. data_directory = r"PointClouds/Widefield_20x_ROIs_Pdgfrb/"
  388. save_directory = r"D:/Research/Stroke_PDGFR-B_Reactivity/Images/Widefield_20x-ROIs_Pdgfrb_TDA/PersistenceDiagrams"
  389. save_diag_directory = r"Data_Raw/Widefield_20x_ROIs_Pdgfrb/PersistenceDiagrams/"
  391. #x_range = (0, 500)
  392. #y_range = (0, 500)
  394. # Call the function to process, visualize, and save persistence diagrams
  395. process_npy_for_persistence_diagrams(data_directory, save_directory, save_diag_directory)
  398. # %% [markdown]
  399. # ## Calculate Betti curves for the persistent diagrams
  401. # %% [markdown]
  402. # Now, we use Betti curves to encapsulate the evolution of topological features across a spectrum of different scales. We can visualize how connected components ($\beta_0$) and loops ($\beta_1$) change as a function of a scaling parameter,
  403. # The Betti curve for each dimension is a function mapping the scaling parameter to the corresponding Betti number (not shown here).
  404. #
  405. # Here, we compute Betti curves in different chunks for 0D-Homology and 1D-homology. To eliminte bias for the number of connected components, we normalize 0D-hoology Betti curves bubt not 1D-homology (loops). We calculate the Betti curvas by region and plot different lines for each DPI.
  407. # %% [markdown]
  408. # ### Betti curves for 0D homology
  410. # %%
  411. def load_persistence_diagrams(dim0_path):
  412. files = [f for f in os.listdir(dim0_path) if 'dim0' in f]
  413. data = {}
  414. for file in files:
  415. parts = file.split('_')
  416. day, region = parts[1], parts[2]
  417. diagram = np.load(os.path.join(dim0_path, file))
  418. key = (day, region)
  419. data[key] = diagram
  420. return data
  422. def compute_betti_0_curve(persistence_diagram, max_filtration_value, n_steps):
  423. betti_0_curve = []
  424. filtration_values = np.linspace(0, max_filtration_value, n_steps)
  425. for value in filtration_values:
  426. betti_0 = sum(1 for point in persistence_diagram if point[0] <= value and (point[1] > value or point[1] == float('inf')))
  427. betti_0_curve.append(betti_0)
  429. # Normalize the curve to start at 1
  430. if betti_0_curve[0] != 0: # Ensure there is no division by zero
  431. initial_value = betti_0_curve[0]
  432. normalized_betti_0_curve = [b / initial_value for b in betti_0_curve]
  433. else:
  434. # Handle the case where the initial value is 0
  435. normalized_betti_0_curve = betti_0_curve
  437. return normalized_betti_0_curve
  439. dim0_path = "Data_Raw/Widefield_20x_ROIs_Pdgfrb/PersistenceDiagrams/"
  440. dim0_data = load_persistence_diagrams(dim0_path)
  442. n_steps = 300
  443. max_filtration_value = 300
  444. betti_0_curves = {key: compute_betti_0_curve(diagram, max_filtration_value, n_steps) for key, diagram in dim0_data.items()}
  446. # Aggregating Betti curves by region
  447. aggregated_curves_by_region_0 = {}
  448. for (day, region), curve in betti_0_curves.items():
  449. if region not in aggregated_curves_by_region_0:
  450. aggregated_curves_by_region_0[region] = {}
  451. aggregated_curves_by_region_0[region][day] = curve
  453. def plot_aggregated_betti_curves_0(aggregated_curves_by_region_0, save_directory):
  454. if not os.path.exists(save_directory):
  455. os.makedirs(save_directory)
  457. days_order = ['0D', '3D', '7D', '14D', '30D']
  458. day_colors = {
  459. '0D': 'red',
  460. '3D': 'blue',
  461. '7D': 'green',
  462. '14D': 'purple',
  463. '30D': 'orange'
  464. }
  466. for region, curves in aggregated_curves_by_region_0.items():
  467. plt.figure(figsize=(5, 5))
  468. for day in days_order:
  469. if day in curves:
  470. plt.plot(curves[day], label=day, color=day_colors[day])
  472. plt.xlabel("Filtration Value", fontsize=22)
  473. plt.ylabel("Normalized Betti-0", fontsize=22)
  474. plt.xticks(fontsize=22, rotation=45)
  475. plt.yticks(fontsize=22)
  476. leg = plt.legend(title="Day", fontsize=16)
  477. plt.setp(leg.get_title(), fontsize=16)
  478. plt.tight_layout()
  480. plot_path = os.path.join(save_directory, f"Betti0_Curves_{region}.png")
  481. plt.savefig(plot_path)
  482. print(f"Saved Betti-0 curve plot for {region} at {plot_path}")
  484. plt.show()
  485. plt.close()
  487. save_directory = "D:/Research/Stroke_PDGFR-B_Reactivity/Pdgfrb_Reactivity_DataAnalysis/Stroke_Pdgfrb_Reactivity/Plots/Widefield_20x_ROIs_pdgfrb_TDA"
  488. plot_aggregated_betti_curves_0(aggregated_curves_by_region_0, save_directory)
  490. # %% [markdown]
  491. # ### Betti curves for 1D homology
  492. #
  493. # Now, we calculate the Betti curves for 1D homoogy.
  495. # %%
  496. def load_persistence_diagrams(dim1_path):
  497. files = [f for f in os.listdir(dim1_path) if 'dim1' in f]
  498. data = {}
  499. for file in files:
  500. parts = file.split('_')
  501. day, region = parts[1], parts[2]
  502. diagram = np.load(os.path.join(dim1_path, file))
  503. key = (day, region)
  504. data[key] = diagram
  505. return data
  507. dim1_path = "Data_Raw/Widefield_20x_ROIs_Pdgfrb/PersistenceDiagrams/"
  508. dim1_data = load_persistence_diagrams(dim1_path)
  511. def compute_betti_1_curve(persistence_diagram, max_filtration_value, n_steps):
  512. betti_1_curve = []
  513. filtration_values = np.linspace(0, max_filtration_value, n_steps)
  514. for value in filtration_values:
  515. betti_1 = sum(1 for point in persistence_diagram if point[0] <= value and (point[1] > value or point[1] == float('inf')))
  516. betti_1_curve.append(betti_1)
  517. return betti_1_curve
  519. n_steps = 100
  520. max_filtration_value = 700
  521. betti_1_curves = {key: compute_betti_1_curve(diagram, max_filtration_value, n_steps) for key, diagram in dim1_data.items()}
  523. # Aggregating Betti curves by region
  524. aggregated_curves_by_region_1 = {}
  525. for (day, region), curve in betti_1_curves.items():
  526. if region not in aggregated_curves_by_region_1:
  527. aggregated_curves_by_region_1[region] = {}
  528. aggregated_curves_by_region_1[region][day] = curve
  531. def plot_aggregated_betti_curves_1(aggregated_curves_by_region_1, save_directory):
  532. if not os.path.exists(save_directory):
  533. os.makedirs(save_directory)
  535. days_order = ['0D', '3D', '7D', '14D', '30D']
  536. day_colors = {
  537. '0D': 'red',
  538. '3D': 'blue',
  539. '7D': 'green',
  540. '14D': 'purple',
  541. '30D': 'orange'
  542. }
  544. for region, curves in aggregated_curves_by_region_1.items():
  545. plt.figure(figsize=(5, 5)) # Adjust figure size as needed
  546. for day in days_order:
  547. if day in curves:
  548. plt.plot(curves[day], label=day, color=day_colors[day])
  549. plt.xlabel("Filtration Value", fontsize=22)
  550. plt.ylabel("Betti-1", fontsize=22)
  551. plt.xticks(fontsize=22) # Increase tick label size
  552. plt.yticks(fontsize=22) # Increase tick label size
  553. leg = plt.legend(title="Day", fontsize=16)
  554. plt.setp(leg.get_title(), fontsize=16) # Increase legend title size
  555. plt.tight_layout()
  557. # Save the plot before calling plt.show()
  558. plot_path = os.path.join(save_directory, f"Betti1_Curves_{region}.png")
  559. plt.savefig(plot_path)
  560. print(f"Saved Betti-1 curve plot for {region} at {plot_path}")
  562. # Show the plot inline
  563. plt.show()
  564. plt.close()
  566. # Specify the directory where you want to save the plots
  567. save_directory = "Plots/Widefield_20x_ROIs_pdgfrb_TDA/"
  569. # Call the function with your aggregated Betti curves and the save directory
  570. plot_aggregated_betti_curves_1(aggregated_curves_by_region_1, save_directory)
  573. # %% [markdown]
  574. # ## Calculate Bottleneck distances
  575. #
  576. # Now, we use the bottleneck distance to compare the persistent diagrams. Given two persistence diagrams D1 and D2, the bottleneck distance measures the minimum amount one needs to move points in D1 to match points in D2 and vice versa, such that the worst-case (or maximum) movement is minimized. Mathematically, the bottleneck distance is defined as:
  577. #
  578. # $$
  579. # \text{bottleneck}(D_1, D_2) = \inf_{\phi: D_1 \to D_2} \sup_{x \in D_1} ||x - \phi(x)||_{\infty}
  580. # $$
  581. #
  582. # Performing bootstrapping, we estimate the distribution of the bottleneck distances by resampling with replacement from the data. Then, we repeat the sampling 1000 times to compute 95% confidence intervals.
  583. #
  584. # In our analysis, we load organize persistence diagrams corresponding to various days, regions, and proteins from designated file paths. The diagrams are gouped based on the mentioned factors and perform specific pairwise day comparisons. We use bootstrap iterations to obtain bottleneck distances between the persistence diagrams of two distinct days, representing the topological disparity between them. These results are stored in a CSV file and the visualized as a heat map.
  586. # %% [markdown]
  587. # ### Load persistance diagrams
  588. #
  589. # We load the persistent diagrams by **0D** or **1D** homology Interchangeably. From each diagram, we obtain the grouping metadata and perform the bootstraping procedure to calculate Bottlenecks distances. Please note this procedure is expensive and make take hours to complete.
  590. #
  591. # **Note:** To shift focus to dimension 1 persistence diagrams, one can replace 'dim0.npy' with 'dim1.npy' in the code.
  593. # %%
  594. ### To calculate dimension 1, change the dim0 to dim1
  596. # Define the path to persistence diagrams
  597. path = "Data_Raw/Widefield_20x_ROIs_Pdgfrb/PersistenceDiagrams/"
  599. # Load .npy files for dimension 1 diagrams
  600. files = [f for f in os.listdir(path) if f.endswith('dim1.npy')]
  602. # Organize data into a nested dictionary {day: {region: [diagrams, ...], ...}, ...}
  603. data = {}
  604. for file in files:
  605. # Extract information from the filename
  606. parts = file.split("_")
  607. day, region = parts[1], parts[2]
  608. # Load the persistence diagram
  609. diagram = np.load(os.path.join(path, file))
  610. # Initialize day if not already in data
  611. if day not in data:
  612. data[day] = {}
  613. # Append the diagram to the correct day and region
  614. data[day].setdefault(region, []).append(diagram)
  616. # %% [markdown]
  617. # Next, we verifying the loaded data structure
  619. # %%
  620. for day in data:
  621. print(f"Day: {day}")
  622. for region in data[day]:
  623. num_diagrams = len(data[day][region])
  624. print(f" Region: {region}, Number of diagrams: {num_diagrams}")
  625. # Optionally, print out some details of the first diagram as an example
  626. if num_diagrams > 0:
  627. example_diagram = data[day][region][0]
  628. print(f" Example diagram size: {example_diagram.shape}")
  629. print(f" Example diagram content (first few pairs): {example_diagram[:5]}")
  630. else:
  631. print(" No diagrams for this region and day.")
  632. print("\n")
  634. # %% [markdown]
  635. # ### Computing distances by boostraping and saving metrics
  636. #
  637. # Now we perform the analysis between persistence diagrams for specified pairwise day comparisons.First, we define the specific pairwise day comparisons of interest (0D-3D, 3D-7D, 7D-14D, 14D-30D) and extract unique days, regions, and proteins from the data. For each unique combination of Region and Protein, we extract persistence diagrams associated with each day. The region Peri (Perilesional) does not contain information from 0D. Therefore, we skip to the next iteration.
  639. # %%
  640. # Define the specified pairwise day comparisons
  641. day_pairs = [("0D", "3D"), ("3D", "7D"), ("7D", "14D"), ("14D", "30D")]
  643. # Extract unique days and regions from the data
  644. unique_days = set(data.keys())
  645. unique_regions = set(region for day in data for region in data[day])
  647. results = []
  649. # Define number of bootstrap replications
  650. num_replications = 1000
  652. # For each unique combination of Day and Region
  653. for region in unique_regions:
  654. # Perform the specified pairwise comparison of diagrams from different days
  655. for day1, day2 in day_pairs:
  656. bootstrap_distances = []
  658. for _ in range(num_replications):
  659. # Randomly sample diagrams from day1 and day2 with replacement
  660. bootstrap_sample1 = choices(data[day1].get(region, []), k=len(data[day1].get(region, [])))
  661. bootstrap_sample2 = choices(data[day2].get(region, []), k=len(data[day2].get(region, [])))
  663. # Compute pairwise bottleneck distances for the bootstrap samples
  664. distances = [bottleneck_distance(diag1, diag2) for diag1 in bootstrap_sample1 for diag2 in bootstrap_sample2 if diag1.size and diag2.size]
  666. # Only proceed if there are distances to calculate
  667. if distances:
  668. # Compute the median distance for this bootstrap replication
  669. bootstrap_distances.append(np.median(distances))
  671. # Skip if no distances were calculated
  672. if not bootstrap_distances:
  673. continue
  675. # Compute metrics
  676. median_distance = np.median(bootstrap_distances)
  677. confidence_lower = np.percentile(bootstrap_distances, 2.5)
  678. confidence_upper = np.percentile(bootstrap_distances, 97.5)
  680. # Store the results
  681. results.append((day1, day2, region, median_distance, confidence_lower, confidence_upper))
  683. # Define the results path and save results to CSV
  684. results_path = "Data_Processed//Widefield_20x_ROIs_Pdgfrb/"
  686. # Change the name for 1D
  687. csv_file = os.path.join(results_path, "Widefield_20x_ROIs_Pdgfrb_Bottleneck_0d.csv")
  689. with open(csv_file, 'w', newline='') as csvfile:
  690. csvwriter = csv.writer(csvfile)
  691. # Write header
  692. csvwriter.writerow(['Day1', 'Day2', 'Region', 'Median Bottleneck Distance', '95% CI Lower', '95% CI Upper'])
  693. # Write the results
  694. for row in results:
  695. csvwriter.writerow(row)
  697. print(f"Bootstrap results saved to {csv_file}.")
  700. # %% [markdown]
  701. # ### Visualization of bottleneck distances by heatmap
  703. # %% [markdown]
  704. # We vislauzie the bottleneck distances by heatmap using seaborn:
  706. # %%
  707. # Load CSV data for 0-dimension and 1-dimension homology
  708. csv_path_dim0 = os.path.join("Data_Processed/Widefield_20x_ROIs_Pdgfrb/Widefield_20x_ROIs_Pdgfrb_Bottleneck_0d.csv")
  709. csv_path_dim1 = os.path.join("Data_Processed/Widefield_20x_ROIs_Pdgfrb/Widefield_20x_ROIs_Pdgfrb_Bottleneck_1d.csv")
  710. data_dim0 = pd.read_csv(csv_path_dim0)
  711. data_dim1 = pd.read_csv(csv_path_dim1)
  713. # Assuming the structure of your CSV files matches the expected columns
  715. # Define label mappings and structure
  716. regions = ["Peri", "Str", "Ctx"]
  717. region_labels = ["Perilesion", "Striatum", "Cortex"]
  718. days = ["0D-3D", "3D-7D", "7D-14D", "14D-30D"]
  719. dimensions = ["dim0", "dim1"] # Representing 0-dimension and 1-dimension homology
  721. # Initialize an empty dictionary to store heatmap data for each dimension
  722. heatmap_data = {dim: pd.DataFrame(index=regions, columns=days) for dim in dimensions}
  723. heatmap_annotations = {dim: pd.DataFrame(index=regions, columns=days) for dim in dimensions}
  725. # Function to populate matrices with data from CSV
  726. def populate_heatmap(data, dim):
  727. for _, row in data.iterrows():
  728. day_pair, region, median_distance, lower, upper = row['Day1']+'-'+row['Day2'], row['Region'], round(row['Median Bottleneck Distance']), round(row['95% CI Lower']), round(row['95% CI Upper'])
  729. if region in regions:
  730. heatmap_data[dim].loc[region, day_pair] = median_distance
  731. # Adjust annotation format here to round and remove decimals
  732. annotation = f"{median_distance}\n({lower} - {upper})"
  733. heatmap_annotations[dim].loc[region, day_pair] = annotation
  735. populate_heatmap(data_dim0, "dim0")
  736. populate_heatmap(data_dim1, "dim1")
  738. # Convert the data to numeric and handle NaN values
  739. for dim in dimensions:
  740. heatmap_data[dim] = heatmap_data[dim].astype(float)
  741. heatmap_data[dim].fillna(-1, inplace=True) # Fill NaN values with -1 for missing data
  743. # Visualization
  744. fig, axes = plt.subplots(1, 2, figsize=(16, 7), sharey=True)
  745. for ax, dim in zip(axes, dimensions):
  746. sns.heatmap(heatmap_data[dim], annot=heatmap_annotations[dim], cmap='viridis', ax=ax,
  747. cbar_kws={'label': 'Bottleneck Distance', 'shrink': 0.5},
  748. mask=(heatmap_data[dim] == -1), fmt="", annot_kws={"size": 16})
  749. ax.set_title(f"Dimension: {dim[-1]}", fontsize=16)
  750. ax.set_xlabel("Day Pair", fontsize=16)
  751. ax.set_ylabel("Region", fontsize=16)
  752. ax.set_yticklabels(region_labels, fontsize=16, rotation=0)
  753. ax.tick_params(axis='x', labelsize=16)
  755. # Adjust colorbar label size
  756. cbar = ax.collections[0].colorbar
  757. cbar.ax.tick_params(labelsize=16)
  758. cbar.ax.set_ylabel('Bottleneck Distance', fontsize=16)
  760. plt.tight_layout()
  761. plt.show()
  763. # Update the output directory path for saving individual images
  764. output_dir = "Plots/Widefield_20x_ROIs_pdgfrb_TDA"
  766. # Ensure the output directory exists
  767. if not os.path.exists(output_dir):
  768. os.makedirs(output_dir)
  770. # Save individual heatmap images for each dimension
  771. for dim in dimensions:
  772. fig, ax = plt.subplots(figsize=(7, 5))
  773. sns.heatmap(heatmap_data[dim], annot=heatmap_annotations[dim], cmap='viridis', ax=ax,
  774. cbar_kws={'label': 'Bottleneck Distance'},
  775. mask=(heatmap_data[dim] == -1), fmt="s")
  776. ax.set_title(f"Bottleneck Distances (Dimension: {dim[-1]})", fontsize=16)
  777. ax.set_xlabel("Day Pair", fontsize=16)
  778. ax.set_ylabel("Region", fontsize=16)
  779. ax.set_yticklabels(region_labels, fontsize=16)
  780. ax.tick_params(axis='x', labelsize=16)
  782. # Adjust colorbar label size
  783. cbar = ax.collections[0].colorbar
  784. cbar.ax.tick_params(labelsize=16)
  785. cbar.ax.set_ylabel('Bottleneck Distance', fontsize=16)
  787. plt.tight_layout()
  789. # Save the figure
  790. fig.savefig(os.path.join(output_dir, f"Bottleneck_Heatmap_{dim}.png"))
  791. plt.close(fig)
  794. # %% [markdown]
  795. # In the graphs, we observed the **median** Bottleneck distances and 95% onfidence intervals for day pairs.
  797. # %% [markdown]
  798. # ## Calculate Wasserstain distances
  800. # %% [markdown]
  801. # We calculate as well Wasserstein distances, also known as the Earth Mover's or Kantorovich-Rubinstein distance. This metric comprises all possible pairs of points between the two diagrams and weighs them by their distance. This represents the minimum "effort" or "cost" to transform one persistence diagram into another.
  803. # %% [markdown]
  804. # ### Load persistance diagrams
  806. # %% [markdown]
  807. # **Noe** Load the data from the Bottleneck distances section. Please make sure to change dim0 or dim1 as needed.
  809. # %% [markdown]
  810. # ### Compute Wasserstein distances by boostrapping and saving metrics
  812. # %%
  813. # Define the specified pairwise day comparisons
  814. day_pairs = [("0D", "3D"), ("3D", "7D"), ("7D", "14D"), ("14D", "30D")]
  816. # Extract unique days and regions from the data
  817. unique_days = set(data.keys())
  818. unique_regions = set(region for day in data for region in data[day])
  820. results = []
  822. # Define number of bootstrap replications
  823. num_replications = 1000
  825. # For each unique combination of Day and Region
  826. for region in unique_regions:
  827. # Perform the specified pairwise comparison of diagrams from different days
  828. for day1, day2 in day_pairs:
  829. if day1 in data and day2 in data and region in data[day1] and region in data[day2]:
  830. bootstrap_distances = []
  832. for _ in range(num_replications):
  833. # Randomly sample diagrams from day1 and day2 with replacement
  834. bootstrap_sample1 = choices(data[day1].get(region, []), k=len(data[day1].get(region, [])))
  835. bootstrap_sample2 = choices(data[day2].get(region, []), k=len(data[day2].get(region, [])))
  837. # Compute pairwise Wasserstein distances for the bootstrap samples
  838. distances = [wasserstein_distance(diag1, diag2, order=2) for diag1 in bootstrap_sample1 for diag2 in bootstrap_sample2 if diag1.size and diag2.size]
  840. # Only proceed if there are distances to calculate
  841. if distances:
  842. # Compute the median distance for this bootstrap replication
  843. bootstrap_distances.append(np.median(distances))
  845. # Skip if no distances were calculated
  846. if not bootstrap_distances:
  847. continue
  849. # Compute metrics
  850. median_distance = np.median(bootstrap_distances)
  851. confidence_lower = np.percentile(bootstrap_distances, 2.5)
  852. confidence_upper = np.percentile(bootstrap_distances, 97.5)
  854. # Store the results
  855. results.append((day1, day2, region, median_distance, confidence_lower, confidence_upper))
  857. # Define the results path and save results to CSV
  858. results_path = "D:/Research/Stroke_PDGFR-B_Reactivity/Pdgfrb_Reactivity_DataAnalysis/Stroke_Pdgfrb_Reactivity/Data_Raw/Widefield_20x_ROIs_Pdgfrb/"
  859. csv_file = os.path.join(results_path, "Widefield_20x_ROIs_Pdgfrb_Wasserstein_1d.csv")
  861. with open(csv_file, 'w', newline='') as csvfile:
  862. csvwriter = csv.writer(csvfile)
  863. # Write header
  864. csvwriter.writerow(['Day1', 'Day2', 'Region', 'Median Wasserstein Distance', '95% CI Lower', '95% CI Upper'])
  865. # Write the results
  866. for row in results:
  867. csvwriter.writerow(row)
  869. print(f"Wasserstein distance bootstrap results saved to {csv_file}.")
  872. # %% [markdown]
  873. # ### Visualization of Wasserstein distances by heatmap
  875. # %% [markdown]
  876. # We plot the calculated Wasserstein distances and save corresponding plots in the designated folder.
  878. # %%
  879. # Load CSV data for 0-dimension and 1-dimension homology
  880. csv_path_dim0 = os.path.join("Data_Processed/Widefield_20x_ROIs_Pdgfrb/Widefield_20x_ROIs_Pdgfrb_Wasserstein_0d.csv")
  881. csv_path_dim1 = os.path.join("Data_Processed/Widefield_20x_ROIs_Pdgfrb/Widefield_20x_ROIs_Pdgfrb_Wasserstein_1d.csv")
  882. data_dim0 = pd.read_csv(csv_path_dim0)
  883. data_dim1 = pd.read_csv(csv_path_dim1)
  885. # Assuming the structure of your CSV files matches the expected columns
  887. # Define label mappings and structure
  888. regions = ["Peri", "Str", "Ctx"]
  889. region_labels = ["Perilesion", "Striatum", "Cortex"]
  890. days = ["0D-3D", "3D-7D", "7D-14D", "14D-30D"]
  891. dimensions = ["dim0", "dim1"] # Representing 0-dimension and 1-dimension homology
  893. # Initialize an empty dictionary to store heatmap data for each dimension
  894. heatmap_data = {dim: pd.DataFrame(index=regions, columns=days) for dim in dimensions}
  895. heatmap_annotations = {dim: pd.DataFrame(index=regions, columns=days) for dim in dimensions}
  897. # Function to populate matrices with data from CSV
  898. def populate_heatmap(data, dim):
  899. for _, row in data.iterrows():
  900. day_pair, region, median_distance, lower, upper = row['Day1']+'-'+row['Day2'], row['Region'], round(row['Median Wasserstein Distance']), round(row['95% CI Lower']), round(row['95% CI Upper'])
  901. if region in regions:
  902. heatmap_data[dim].loc[region, day_pair] = median_distance
  903. # Adjust annotation format here to round and remove decimals
  904. annotation = f"{median_distance}\n({lower} - {upper})"
  905. heatmap_annotations[dim].loc[region, day_pair] = annotation
  907. populate_heatmap(data_dim0, "dim0")
  908. populate_heatmap(data_dim1, "dim1")
  910. # Convert the data to numeric and handle NaN values
  911. for dim in dimensions:
  912. heatmap_data[dim] = heatmap_data[dim].astype(float)
  913. heatmap_data[dim].fillna(-1, inplace=True) # Fill NaN values with -1 for missing data
  915. # Visualization
  916. fig, axes = plt.subplots(1, 2, figsize=(16, 7), sharey=True)
  917. for ax, dim in zip(axes, dimensions):
  918. sns.heatmap(heatmap_data[dim], annot=heatmap_annotations[dim], cmap='viridis', ax=ax,
  919. cbar_kws={'label': 'Wasserstein Distance', 'shrink': 0.5},
  920. mask=(heatmap_data[dim] == -1), fmt="", annot_kws={"size": 16})
  921. ax.set_title(f"Dimension: {dim[-1]}", fontsize=16)
  922. ax.set_xlabel("Day Pair", fontsize=16)
  923. ax.set_ylabel("Region", fontsize=16)
  924. ax.set_yticklabels(region_labels, fontsize=16, rotation=0)
  925. ax.tick_params(axis='x', labelsize=16)
  927. # Adjust colorbar label size
  928. cbar = ax.collections[0].colorbar
  929. cbar.ax.tick_params(labelsize=16)
  930. cbar.ax.set_ylabel('Wasserstein Distance', fontsize=16)
  932. plt.tight_layout()
  933. plt.show()
  935. # Update the output directory path for saving individual images
  936. output_dir = "Plots/Widefield_20x_ROIs_pdgfrb_TDA"
  938. # Ensure the output directory exists
  939. if not os.path.exists(output_dir):
  940. os.makedirs(output_dir)
  942. # Save individual heatmap images for each dimension
  943. for dim in dimensions:
  944. fig, ax = plt.subplots(figsize=(7, 5))
  945. sns.heatmap(heatmap_data[dim], annot=heatmap_annotations[dim], cmap='viridis', ax=ax,
  946. cbar_kws={'label': 'Wasserstein Distance'},
  947. mask=(heatmap_data[dim] == -1), fmt="s")
  948. ax.set_title(f"Wasserstein Distances (Dimension: {dim[-1]})", fontsize=16)
  949. ax.set_xlabel("Day Pair", fontsize=16)
  950. ax.set_ylabel("Region", fontsize=16)
  951. ax.set_yticklabels(region_labels, fontsize=16)
  952. ax.tick_params(axis='x', labelsize=16)
  954. # Adjust colorbar label size
  955. cbar = ax.collections[0].colorbar
  956. cbar.ax.tick_params(labelsize=16)
  957. cbar.ax.set_ylabel('Wasserstein Distance', fontsize=16)
  959. plt.tight_layout()
  961. # Save the figure
  962. fig.savefig(os.path.join(output_dir, f"Wasserstein_Heatmap_{dim}.png"))
  963. plt.close(fig)
  966. # %% [markdown]
  967. # Analogous to Bottleneck distances, here we appreacite the estimates and its uncertainty for each day pair.
  969. # %% [markdown]
  970. # ## Addon: Generate another set of images for illustration
  972. # %% [markdown]
  973. # Here, we create an additional visualization of Vietoris-rips complexes plotting the filtration circles at different values. We plot representative images and then store a full set of images in a designated location.
  975. # %%
  976. def load_point_cloud(file_path):
  977. return np.load(file_path, allow_pickle=True)
  979. def get_edges_at_filtration(points, filtration_value):
  980. rips_complex = RipsComplex(points=points, max_edge_length=filtration_value)
  981. simplex_tree = rips_complex.create_simplex_tree(max_dimension=1)
  983. edges = []
  984. for simplex, f_value in simplex_tree.get_filtration():
  985. if len(simplex) == 2 and f_value <= filtration_value:
  986. edges.append(simplex)
  988. return edges
  990. def plot_point_cloud_with_structures(points, edges, ax, filtration_value):
  991. ax.scatter(points[:, 0], points[:, 1], color='green', s=10, zorder=1)
  992. for edge in edges:
  993. point1 = points[edge[0]]
  994. point2 = points[edge[1]]
  995. ax.plot([point1[0], point2[0]], [point1[1], point2[1]], 'black', alpha=0.9, zorder=1.5)
  996. for point in points:
  997. circle = Circle(point, radius=filtration_value, color='skyblue', fill=True, alpha=0.2)
  998. ax.add_patch(circle)
  999. ax.set_xticks([])
  1000. ax.set_yticks([])
  1001. ax.set_xlabel('')
  1002. ax.set_ylabel('')
  1004. file_paths = ['PointClouds/Widefield_20x_ROIs_Pdgfrb/Td047_30D_Peri.npy',
  1005. 'PointClouds/Widefield_20x_ROIs_Pdgfrb/Td047_30D_Str.npy',
  1006. 'PointClouds/Widefield_20x_ROIs_Pdgfrb/Td047_30D_Ctx.npy']
  1007. filtration_values = [0, 50, 100, 200]
  1009. # Ensure each point cloud gets its subplot for each filtration value
  1010. fig, axs = plt.subplots(len(file_paths), len(filtration_values), figsize=(15, 5 * len(file_paths)))
  1012. for file_index, file_path in enumerate(file_paths):
  1013. points = load_point_cloud(file_path)
  1014. for ax_index, filtration_value in enumerate(filtration_values):
  1015. ax = axs[file_index, ax_index] if len(file_paths) > 1 else axs[ax_index]
  1016. edges = get_edges_at_filtration(points, filtration_value)
  1017. plot_point_cloud_with_structures(points, edges, ax, filtration_value)
  1019. plt.tight_layout()
  1020. plt.show()
  1023. # %% [markdown]
  1024. # Finally, we save a series f images in a destinated folder for illustration in the scientific article.
  1026. # %%
  1027. import os
  1029. # Load all .npy files in a folder
  1030. def load_point_clouds_from_folder(folder_path):
  1031. point_cloud_files = [os.path.join(folder_path, f) for f in os.listdir(folder_path) if f.endswith('.npy')]
  1032. return point_cloud_files
  1034. # Save the plot for a given point cloud and filtration values
  1035. def save_plots_for_point_cloud(file_path, filtration_values, save_folder):
  1036. points = load_point_cloud(file_path)
  1037. for filtration_value in filtration_values:
  1038. edges = get_edges_at_filtration(points, filtration_value)
  1040. fig, ax = plt.subplots(figsize=(5, 5)) # Adjust size as needed
  1041. plot_point_cloud_with_structures(points, edges, ax, filtration_value)
  1043. # Save plot to the designated folder
  1044. base_filename = os.path.splitext(os.path.basename(file_path))[0]
  1045. save_path = os.path.join(save_folder, f'{base_filename}_filtration_{filtration_value}.png')
  1046. plt.savefig(save_path)
  1047. plt.close(fig) # Close the figure after saving to free up memory
  1049. # Specify folder paths
  1050. point_clouds_folder = 'PointClouds/Widefield_20x_ROIs_Pdgfrb/'
  1051. save_folder = 'D:/Research/Stroke_PDGFR-B_Reactivity/Images/Widefield_20x-ROIs_Pdgfrb_TDA/FiltrationValues'
  1052. filtration_values = [0, 50, 100, 200]
  1054. # Ensure the save folder exists
  1055. if not os.path.exists(save_folder):
  1056. os.makedirs(save_folder)
  1058. # Load and process each point cloud file
  1059. point_cloud_files = load_point_clouds_from_folder(point_clouds_folder)
  1060. for file_path in point_cloud_files:
  1061. save_plots_for_point_cloud(file_path, filtration_values, save_folder)
  1064. # %%

Widefield_20x_ROIs_Pdgfrb_TDA.ipynb at commit ae087d5, no license · at the source

Overview

Authors: Daniel Manrique‐Castano1,2, Ayman ElAli1,2
ORCID iDs: Daniel Manrique‐Castano, Ayman ElAli
  1. Neuroscience Axis, Research Center of CHU de Québec‐Université Laval Quebec City Quebec Canada
  2. Department of Psychiatry and Neuroscience, Faculty of Medicine Université Laval Quebec City Quebec Canada
Journal: Glia, volume 75, issue 1, article e70232
Dates: received 2 April 2026; accepted 12 September 2026; published online 31 October 2026; in print January 2027
Type: Research article · Language: English
License: CC BY
Identifiers: DOI · PMCID PMC13723038
Status: code verified (how the code is found and verified)
Categories: stroke (population), cellular / molecular (subfield)
Methods: Smoothing, state filtering, decompositions, Machine learning, Statistics, fMRI & imaging
Keywords: ischemic stroke, KLF4, PDGFRβ + cells, point pattern analysis (PPA), tissue scar, topological data analysis (TDA)
Funding: Gouvernement du Canada | Canadian Institutes of Health Research (IRSC) (169062, 186148); Heart and Stroke Foundation of Canada (Heart and Stroke Foundation) (G‐23‐0035047); Fonds de Recherche du Québec - Santé (318466)
Citations: not cited yet (Europe PMC); 79 references in the paper
Research resources: PDGFRα immunolabeling RRID:AB_2, RRID:AB_2227288, RRID:AB_777165, RRID:AB_839504, some GFAP‐stained cytoskeletons RRID:AB_86543, homozygous for B6.129S6‐Klf4tm1khk/Mmmh RRID:MMRRC_029877

Abstract

The pathobiology of ischemic stroke is associated with the activation of brain platelet‐derived growth factor receptor (PDGFR)β + perivascular cells. While Krüppel‐like factor (KLF)4 is a well‐established mediator of peripheral PDGFRβ + cell activation, its role in brain PDGFRβ + cells remains elusive. Herein, we aimed to characterize the spatiotemporal reactivity of brain PDGFRβ + cells, while assessing KLF4 contribution upon experimental ischemic stroke using transgenic mice that enable tracking PDGFRβ + cells or conditionally depleting KLF4 in PDGFRβ + cells. Moving beyond the qualitative descriptors, we sought to develop point pattern analysis (PPA) and topological data analysis (TDA) workflows to generate robust and reusable mathematical models to assess spatial cell responses and tissue scar topological changes. First, we report the emergence of phenotypically distinct reactive parenchymal PDGFRβ + cells, which cannot be causally and exclusively linked, as it is currently proposed, to local cell proliferation or adjacent vascular detachment. Furthermore, we show that KLF4 is barely expressed in brain PDGFRβ + cells under normal conditions and that it is only marginally induced in reactive PDGFRβ + cells at the injured brain. Notably, the conditional attenuation of KLF4 expression specifically in reactive PDGFRβ + cells does not alter spatial cellular distribution nor scar organization. This suggests that, in contrast to peripheral PDGFRβ + cells upon ischemic insults, KLF4 is not a main regulator of the reactivity of brain cells in response to ischemic stroke. Ultimately, our findings supported by open‐source spatiotemporal models indicate that the reactivity of brain PDGFRβ + perivascular cells is highly heterogeneous and fundamentally distinct from peripheral PDGFRβ + cells.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 21 matches between paragraphs and lines of code.

Zenodo 13948420

License: CC-BY-4.0
State: the link answers, verified on 1 November 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: “Data Availability Statement”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: tidyverse (27 files), ggplot2 (26 files), brms (23 files), easystats (23 files), patchwork (23 files), emmeans (15 files), broom (5 files), NumPy (4 files), pandas (4 files), Matplotlib (3 files), reshape2 (3 files), SciPy (3 files), seaborn (2 files), data.table (1 file), scikit-image (1 file), tifffile (1 file)
Availability: 1 check, the latest on 1 November 2026: the link answers (HTTP 200)
  • 1 November 2026: the link answers (HTTP 200)
35 files
At the source: doi.org/10.5281/zenodo.13948420

elalilab/stroke_pdgfr-b_reactivity

License: none: the authors keep all their rights
State: the link answers, verified on 1 November 2026
Evidence: files inventoried
Commit: ae087d5288616f7e66b467932c2f2627aa7fb910, 26 March 2026
Languages: Quarto (33), Jupyter (2)
Size: 41 files, 35 scripts
Software Heritage: not archived
Found in: the Zenodo archive record
Holds: README, 35 notebooks
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: tidyverse (28 files), ggplot2 (27 files), brms (24 files), easystats (24 files), patchwork (24 files), emmeans (15 files), broom (6 files), NumPy (4 files), pandas (4 files), reshape2 (4 files), Matplotlib (3 files), SciPy (3 files), seaborn (2 files), data.table (1 file), scikit-image (1 file), tifffile (1 file)
Availability: 1 check, the latest on 1 November 2026: the link answers
  • 1 November 2026: the link answers
36 files
At the source: github.com/elalilab/Stroke_PDGFR-B_Reactivity/tree/v1.0.0

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 69 scripts, each with its path and the digest of its content;
  • 21 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data Availability Statement

All the data associated with the publication are made fully available via https://zenodo.org/records/10553084, which includes original microscopy images stored in Zenodo (10.5281/zenodo.10553084 (https://doi.org/10.5281/zenodo.10553084)) as well as access to complementary data stored in Open Science Framework (OSF) repository (10.17605/OSF.IO/74MQN (https://doi.org/10.17605/OSF.IO/74MQN)), and reproducible workflows for image processing and analysis are available on GitHub (10.5281/zenodo.13948420 (https://doi.org/10.5281/zenodo.13948420)).

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 1 November 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 2 authors, 6 keywords, 3 funders, 70 references, 6 RRIDs.

Cite

This paper

Manrique‐Castano, D., & ElAli, A. (2026). Brain PDGFRβ+ Cells Exhibit Diverse Reactive Phenotypes After Stroke Without Requiring KLF4. Glia, 75(1), e70232.

BibTeX

@article{manriquecastano2026brain,
author = {Manrique‐Castano, Daniel and ElAli, Ayman},
title = {{Brain PDGFRβ+ Cells Exhibit Diverse Reactive Phenotypes After Stroke Without Requiring KLF4}},
journal = {Glia},
year = {2026},
month = oct,
volume = {75},
number = {1},
pages = {e70232},
publisher = {Wiley},
issn = {0894-1491},
pmcid = {PMC13723038}
}

RIS

TY - JOUR
AU - Manrique‐Castano, Daniel
AU - ElAli, Ayman
TI - Brain PDGFRβ+ Cells Exhibit Diverse Reactive Phenotypes After Stroke Without Requiring KLF4
T2 - Glia
J2 - Glia
PY - 2026
DA - 2026/10/31
VL - 75
IS - 1
SP - e70232
SN - 0894-1491
PB - Wiley
LA - en
ER -

CSL-JSON

{
"id": "pmcid:PMC13723038",
"type": "article-journal",
"title": "Brain PDGFRβ+ Cells Exhibit Diverse Reactive Phenotypes After Stroke Without Requiring KLF4",
"container-title": "Glia",
"author": [
{
"family": "Manrique‐Castano",
"given": "Daniel"
},
{
"family": "ElAli",
"given": "Ayman"
}
],
"container-title-short": "Glia",
"volume": "75",
"issue": "1",
"page": "e70232",
"PMCID": "PMC13723038",
"ISSN": "0894-1491",
"publisher": "Wiley",
"language": "en",
"issued": {
"date-parts": [
[
2026,
10,
31
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1007/s00018-024-05244-w [code]
Endothelial PDGF-D contributes to neurovascular protection after ischemic stroke by rescuing pericyte functions.
Journal: Cellular and molecular life sciences : CMLS
In common: broom, patchwork, ggplot2, 1 other tool, stroke, cellular / molecular, 8 references, author Ayman ElAli
[2] doi:10.1016/j.xpro.2024.102989 [code]
Unbiased quantification of the spatial distribution of murine cells using point pattern analysis.
Journal: STAR protocols
In common: brms, easystats, emmeans, 5 other tools, cellular / molecular, 4 references, author Ayman ElAli
[3] doi:10.1016/j.xcrm.2026.102766 [code]
A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: tifffile, patchwork, scikit-image, 9 other tools, cellular / molecular, 2 references
[4] doi:10.1016/j.nicl.2025.103881 [code]
Brain age gap reduction following exercise mirrors clinical improvements in schizophrenia spectrum disorders.
Journal: NeuroImage. Clinical
In common: brms, easystats, emmeans, 9 other tools, 1 reference
[5] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: brms, broom, patchwork, 10 other tools, cellular / molecular
[6] doi:10.1016/j.celrep.2026.117505 [code]
Impaired spatial coding and neuronal hyperactivity in the medial entorhinal cortex of aged APP knock-in mice.
Journal: Cell reports
In common: brms, broom, easystats, 8 other tools, 1 reference
[7] doi:10.1523/jneurosci.1285-25.2025 [code]
Hidden State Inference or Continuous Belief Updating during a Dynamic Visuomotor Skill.
Journal: The Journal of neuroscience : the official journal of the Society for Neuroscience
In common: brms, easystats, emmeans, 7 other tools, 2 references
[8] doi:10.1038/s41467-026-73865-9 [code]
Histamine shapes the neurocomputational dynamics of human learning.
Journal: Nature communications
In common: brms, broom, easystats, 9 other tools
[9] doi:10.1186/s11689-025-09651-3 [code]
The association between infant EEG aperiodic exponent and the trajectory of restricted and repetitive behaviors for toddlers with and without autism.
Journal: Journal of neurodevelopmental disorders
In common: brms, broom, easystats, 9 other tools
[10] doi:10.1038/s41586-022-05279-8 [code]
Inferring and perturbing cell fate regulomes in human brain organoids.
Journal: Nature
In common: brms, easystats, patchwork, 8 other tools, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Add the badge to the code's README

OSCR: the badge of this paper

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown
<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/repo/zenodo:13948420.svg)](https://openscicode.org/badge/repo/zenodo:13948420)
<!-- oscr-badge:end -->
reStructuredText
.. oscr-badge:start

.. image:: https://openscicode.org/badge/repo/zenodo:13948420.svg
   :target: https://openscicode.org/badge/repo/zenodo:13948420
   :alt: OSCR

.. oscr-badge:end
HTML
<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/repo/zenodo:13948420"><img src="https://openscicode.org/badge/repo/zenodo:13948420.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

On GitHub, open the README, choose Edit (the pencil), paste the snippet under the title, then "Commit changes".

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

Versions of its code: the releases of its code tied to a version of this paper (the preprint, the accepted manuscript, the version of record…), each with the tracing map versioned with it.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.