A method for identifying cross-modal features from spatially resolved datasets

Abstract

A method for identifying cross-modal features from two or more spatially resolved datasets is disclosed, the method comprising: (a) registering the two or more spatially resolved datasets to generate an aligned feature image comprising the two or more spatially resolved datasets that are spatially aligned; and (b) extracting the cross-modal features from the aligned feature image. TIFF2024537615000147.tif157166

Description

[Technical field]

[0001] FIELD OF THEINVENTION The present application relates to methods and systems for identifying diagnostic, prognostic, or theranostic information from one or more correlates identified from aligned, spatially resolved data sets. [Background technology]

[0002] background The development of spatially resolved detection modalities has revolutionized diagnosis, prognosis, and theranostics, however, their multimodal application potential remains largely unrealized because each modality is typically analyzed independently of the others.

[0003] New methods utilizing multiple spatially resolved detection modalities for identifying multimodal diagnostics, prognostics, and theranostics are needed. Summary of the Invention

[0004] In one aspect, the present invention provides a method of providing a diagnosis, prognosis, or theranostics for a disease state from three or more image modalities obtained from a biopsy sample from a subject, the method comprising comparing a plurality of cross-modal features to identify a correlation between at least one cross-modal feature parameter and the disease state to identify a diagnosis, prognosis, or theranostics, wherein the plurality of cross-modal features comprises: (a) registering the three or more spatially resolved data sets to generate an aligned feature image, the aligned feature image comprising the three or more spatially resolved data sets that are spatially aligned; and (b) Extracting cross-modal features from the aligned feature images. wherein each cross-modal feature comprises a cross-modal feature parameter, and the three or more spatially resolved data sets are output by a corresponding imaging modality selected from the group consisting of three or more imaging modalities.

[0005] In some embodiments, at least one of the three or more spatially resolved data sets comprises data on the abundance and spatial distribution of cells. In some embodiments, at least one of the three or more spatially resolved data sets comprises data on the abundance and spatial distribution of tissue structures. In some embodiments, at least one of the three or more spatially resolved data sets comprises data on the abundance and spatial distribution of one or more molecular analytes. In some embodiments, the one or more molecular analytes are selected from the group consisting of cells. In some embodiments, the one or more molecular analytes are selected from the group consisting of proteins, antibodies, nucleic acids, lipids, metabolites, carbohydrates, and therapeutic compounds.

[0006] In some embodiments, the biopsy sample is from a subject having or suspected of having a disease whose condition is to be determined. In some embodiments, the disease is type 2 diabetes. In some embodiments, the diagnosis is for diabetic foot ulcers. In some embodiments, the one or more molecular analytes include median distances between distinct immune cell populations. In some embodiments, the one or more molecular analytes include median distances between distinct immune cell populations and tissue structures or diseased cells (e.g., cancer cells). In some embodiments, the one or more molecular analytes include median distances of NK cells from suppressive macrophages, abundance of mature B cells compared to adjacent healthy tissue, and levels of mass spectrometry analytes corresponding to complement proteins, lipoproteins, and metabolites associated with bacteria compared to naturally healing wounds. In some embodiments, the disease is cancer. In some embodiments, the cancer is prostate cancer, lung cancer, kidney cancer, ovarian cancer, or mesothelioma. In some embodiments, the one or more molecular analytes include proteins and analytes associated with immune activity or genomic instability.

[0007] In some embodiments, the method is multiplexed. In some embodiments, the method allows for matching at least 10 molecular analytes. In some embodiments, the method allows for matching at least 20 molecular analytes.

[0008] In one aspect, the present invention provides a method for identifying cross-modal features from two or more spatially resolved datasets, the method comprising: (a) registering the two or more spatially resolved datasets to generate an aligned feature image comprising the two or more spatially resolved datasets that are spatially aligned; and (b) extracting cross-modal features from the aligned feature image.

[0009] In some embodiments, step (a) includes dimensionality reduction for each of the two or more datasets. In some embodiments, the dimensionality reduction is performed by Uniform Manifold Approximation and Projection (UMAP), Isometric Mapping (Isomap), t-SNE, potential of heat diffusion for affinity-based transition embedding (PHATE), Principal Component Analysis (PCA), Diffusion Map, or Non-negative Matrix Factorization (NMF). In some embodiments, the dimensionality reduction is performed by Uniform Manifold Approximation and Projection (UMAP). In some embodiments, step (a) includes optimizing a global spatial alignment in the aligned feature images. In some embodiments, step (a) includes optimizing a local alignment in the aligned feature images.

[0010] In some embodiments, the method further comprises clustering the two or more spatially resolved data sets to complement the data sets with a similarity matrix representing inter-data point similarities. In some embodiments, the clustering comprises extracting a high-dimensional graph from the aligned feature images. In some embodiments, the clustering is performed by the Leiden algorithm, the Louvain algorithm, random walk graph partitioning, spectral clustering, or affinity propagation. In some embodiments, the method comprises predicting cluster assignments to unseen data. In some embodiments, the method comprises modeling cluster-cluster spatial interactions. In some embodiments, the method comprises intensity-based analysis. In some embodiments, the method comprises analysis of cell type abundance or heterogeneity of a given region in the data. In some embodiments, the method comprises analysis of spatial interactions between objects. In some embodiments, the method comprises analysis of type-specific neighborhood interactions. In some embodiments, the method comprises analysis of higher order spatial interactions. In some embodiments, the method comprises analysis of spatial niche predictions.

[0011] In some embodiments, the method further comprises classifying the data, hi some embodiments, the classifying is performed by a hard classifier, a soft classifier, or a fuzzy classifier.

[0012] In some embodiments, the method further comprises defining one or more spatially resolved objects in the aligned feature images. In some embodiments, the method further comprises analyzing the spatially resolved objects. In some embodiments, analyzing the spatially resolved objects comprises segmentation. In some embodiments, the method further comprises inputting one or more landmarks into the aligned feature images.

[0013] In some embodiments, step (b) comprises permutation testing for enrichment or depletion of cross-modal features. In some embodiments, the permutation testing generates a list of p-values ​​and / or identities of enriched or depleted factors. In some embodiments, the permutation testing is performed by a mean permutation test.

[0014] In some embodiments, step (b) includes a multi-domain transformation. In some embodiments, the multi-domain transformation generates a trained model or a predicted output based on the cross-modal features. In some embodiments, the multi-domain transformation is performed by a generative adversarial network or an adversarial autoencoder.

[0015] In some embodiments, at least one of the two or more spatially resolved data sets is an image from immunohistochemistry, imaging mass cytometry, multiplexed ion beam imaging, mass spectrometry imaging, cell staining, RNA-ISH, spatial transcriptome analysis, or co-detection by index imaging. In some embodiments, at least one of the spatially resolved measurement modalities is immunofluorescence imaging. In some embodiments, at least one of the spatially resolved measurement modalities is imaging mass cytometry. In some embodiments, at least one of the spatially resolved measurement modalities is multiplexed ion beam imaging. In some embodiments, at least one of the spatially resolved measurement modalities is mass spectrometry imaging, which is MALDI imaging, DESI imaging, or SIMS imaging. In some embodiments, at least one of the spatially resolved measurement modalities is cell staining, which is H&E, toluidine blue, or fluorescent staining. In some embodiments, at least one of the spatially resolved measurement modalities is RNA-ISH, which is RNAScope. In some embodiments, at least one of the spatially resolved measurement modalities is spatial transcriptome analysis. In some embodiments, at least one of the spatially resolved measurement modalities is co-detection with index imaging.

[0016] In another aspect, the present invention provides a method of identifying a diagnosis, prognosis, or theranostics for a disease state from two or more imaging modalities, the method comprising comparing a plurality of cross-modal features to identify a correlation between at least one cross-modal feature parameter and the disease state to identify a diagnosis, prognosis, or theranostics, wherein the plurality of cross-modal features are identified by the methods described herein, each cross-modal feature comprises a cross-modal feature parameter, and the two or more spatially resolved data sets are output by a corresponding imaging modality selected from the group consisting of two or more imaging modalities.

[0017] In some embodiments, the cross-modal feature parameter is a molecular signature, a single molecule marker, or the abundance of a marker.In some embodiments, the diagnosis, prognosis, or theranostics is individualized to the individual that is the source of two or more spatially resolved data sets.In some embodiments, the diagnosis, prognosis, or theranostics is a population-level diagnosis, prognosis, or theranostics.

[0018] In yet another aspect, the present invention provides a method for identifying trends in a parameter of interest within a plurality of aligned feature images identified by the methods described herein, the method comprising the steps of identifying a parameter of interest in the plurality of aligned feature images, and comparing the parameter of interest between the plurality of aligned feature images to identify trends.

[0019] In still yet another aspect, the present invention provides a computer readable storage medium having stored thereon a computer program for identifying cross-modal features from two or more spatially resolved data sets, the computer program comprising a routine set of instructions for causing a computer to perform steps of the methods described herein.

[0020] In a further aspect, the present invention provides a computer readable storage medium having stored thereon a computer program for identifying a diagnosis, prognosis, or theranostics of a disease state from two or more imaging modalities, the computer program comprising a routine set of instructions for causing a computer to perform the steps of the methods described herein.

[0021] In yet a further aspect, the present invention provides a computer readable storage medium having stored thereon a computer program for identifying trends in a parameter of interest within a plurality of aligned feature images identified by the methods described herein, the computer program comprising a routine set of instructions for causing a computer to perform steps of the methods described herein.

[0022] In an even further aspect, the present invention provides a method of identifying a vaccine, the method comprising the steps of: (a) providing a first dataset of cytometric markers for a disease-naive population; (b) providing a second dataset of cytometric markers for a population suffering from the disease; (c) identifying from the first and second datasets one or more markers that correlate with a clinical or phenotypic measure of the disease; and (d) identifying as the vaccine a composition that (1) is capable of inducing one or more markers that directly correlate with a positive clinical or phenotypic measure of the disease; or (2) identifying as the vaccine a composition that is capable of suppressing one or more markers that directly correlate with a negative clinical or phenotypic measure of the disease. [Brief description of the drawings]

[0023] [Figure 1] FIG. 1 is a diagrammatic representation showing the process of imaging diabetic foot ulcer (DFU) biopsy tissue with multiple modalities, e.g., H&E staining, mass spectrometry imaging (MSI) and imaging mass cytometry (IMC), followed by processing and analysis of the multimodal image dataset using an integrated analysis pipeline. [Diagram 2]FIG. 2A is a high-resolution scanned image showing a DFU biopsy tissue section on a microscopy glass slide. FIG. 2B is a schematic showing a DFU biopsy tissue section on a glass slide before treatment with a spray matrix solution (optimized for each analyte type) containing 40% 2,5-dihydroxybenzoic acid (DHB) in 50:50 v / v acetonitrile:0.1% TFA in water. FIG. 2C is a schematic showing a DFU biopsy tissue section on a glass slide after treatment with a spray matrix solution (optimized for each analyte type) containing 40% 2,5-dihydroxybenzoic acid (DHB) in 50:50 v / v acetonitrile:0.1% TFA in water. FIG. 2D is a graph showing the mass-to-charge average spectrum of a region of DFU tissue obtained after laser desorption, ionization and characterization using mass spectrometry. [Diagram 3] Schematic showing the process underlying imaging of DFU biopsy tissue or cell lines with IMC. Following sample pretreatment, staining with metal-labeled antibodies is performed. Laser ablation of the sample produces aerosolized droplets that are directly transported to the inductively coupled plasma torch of the instrument, producing atomized and ionized sample components. Filtering of unwanted components is performed in a quadrupole ion deflector, where low mass ions and photons are filtered out. High mass ions, mainly representing metal ions associated with labeled antibodies, are further pushed to a time-of-flight (TOF) detector, from which the flight time of each ion is recorded based on the mass-to-charge ratio of each ion, thereby identifying and quantifying the metals present in the sample. Each isotopically labeled sample component is then represented by an isotopic intensity profile, where each peak represents the abundance of each isotope in the sample. Multidimensional analysis is then performed to visualize the data. [Figure 4] 1 is a flow chart summarizing the steps involved in acquiring a multimodal image dataset and extracting molecular signatures from the multimodal dataset. [Diagram 5]Figures 5A-5F are a series of graphs showing the estimation of the intrinsic dimensionality of the MSI dataset using the dimensionality reduction methods t-SNE, Uniform Manifold Approximation and Projection (UMAP), potential of heat diffusion for affinity-based transition embedding (PHATE), Isometry (Isomap), Nonnegative Matrix Factorization (NMF), and Principal Component Analysis (PCA). The convergence to the embedding error value indicated that increasing the dimensionality of the resulting embedding no longer improved the algorithm's ability to capture the complexity of the data. Nonlinear dimensionality reduction methods, e.g., t-SNE, UMAP, PHATE, and Isomap, converged to a much lower intrinsic dimensionality than those of linear methods, e.g., NMF and PCA, indicating that much fewer dimensions are needed to accurately represent the dataset. [Figure 6] Figures 6A and 6B are graphs showing the computational run time for each algorithm across embedding dimensions from 1 to 10. For each method, we plot the mean and standard deviation (n=5) across each number of dimensions. The results show that the nonlinear methods t-SNE and Isomap require longer run times than the nonlinear methods PHATE and UMAP. [Figure 7A] Figure 1 shows a graph comparing the mutual information captured by each dimensionality reduction method tested between a 3D embedding of a grayscale version of MSI data and the corresponding H&E stained tissue section. Mutual information is defined as greater than or equal to zero, with negative values ​​being consistent with a minimum of the cost function in the registration process. Results show that Isomap and UMAP consistently share more information with the H&E images than the other methods tested. [Figure 7B]1 is a scheme illustrating the key technical steps of the analysis described herein. We evaluated the ability of each tested dimensionality reduction method to recover data connectivity (manifold structure) using both the full dataset (noisy) or the denoised dataset (peak-picked). We calculated the denoised manifold-preserving (DeMaP) metric

[18] between the Euclidean distance of the resulting embeddings corresponding to the non-peak-picked data and the geodesic distance in the ambient space (not dimensionally reduced after peak picking) of the corresponding peak-picked data. [Figure 7C] Figure 7B shows the mean and standard deviation DeMaP metric (Spearman's rho correlation coefficient) for all dimensionality reduction methods tested (n=5). The figure shows the correlation results described in Figure 7B. The nonlinear methods Isomap, PHATE and UMAP all consistently preserve the manifold structure without pre-filtering of the data, with consistent correlations above 0.85 across 2-10 dimensions. [Figure 8] FIG. 1 is a schematic flow chart showing the steps from mass spectrometry data and image reconstruction to dimensionality reduction using UMAP and data visualization through a pixelated embedded representation of the mass spectrometry data. [Figure 9] Figure 1 illustrates the mapping of a 3D embedding of MSI data after dimensionality reduction with UMAP onto the original DFU tissue section, where each of the three UMAP dimensions is colored either red (U1), green (U2) or blue (U3). The merged image (RGB image) contains the overlap of all three pseudocolor images. Conversion of the RGB image to grayscale is achieved by adding the pixel intensities for each of the three pseudocolor channels as shown in the equation. For visualization purposes, a weighting factor can be added to each channel (x1, x2, x3) to adjust the signal contribution per channel. A representative grayscale image for the dataset is shown in the pseudocolor image. [Figure 10] 1 is a series of grayscale images of a DFU biopsy tissue sample showing a comparison of various linear and nonlinear dimensionality reduction methods. [Figure 11] A set of images of DFU biopsy tissue acquired by brightfield microscopy (H&E), MSI, and IMC. The spatial resolution of the three imaging modalities is displayed to convey the difference in imaging resolution between brightfield microscopy, MSI, and IMC images. [Figure 12] 1 is a flow chart of representative grayscale DFU biopsy tissue images illustrating the process of image registration between imaging modalities. [Figure 13] 1 is a flowchart illustrating a process for aligning multimodal images with a local region of interest (ROI) approach. [Figure 14] 1 is a flow chart of representative grayscale DFU biopsy tissue images showing the process of fine-tuning local-scale registration. Regions of interest in the toluidine blue image corresponding to each MSI image were selected for local-scale registration. [Figure 15] A series of MSI images (A–C and A''–C'') and IMC images (A'–C' and A''''–C'''') showing three different regions of interest (ROIs) in a DFU biopsy tissue section. Single-cell coordinates on each ROI were identified by segmentation using IMC parameters, and cell types (cell types 1–12) were defined using subsequent clustering analysis on the IMC profiles of the extracted single-cell measurements. These single-cell coordinates were used to extract the corresponding MSI data. Panels A, B, and C show the spatial distribution of the MSI parameters identified through permutation testing. Panels A', B', and C' show the spatial distribution of the IMC markers of interest before single-cell segmentation. Panels A'', B'', and C'' show the overlap of panels A+A', B+B', and C+C'. Panels A'''', B'''', and C'' show the single-cell masks (ROIs defined by single-cell pixel coordinates) identified by segmentation. The color scheme represents the cell types identified by clustering single-cell measurements with respect to the IMC parameters. [Figure 16] Figure 1 illustrates an exemplary workflow of image modality integration (boxes labeled (C)) and modeling of complex tissue states with MIAAIM. Inputs and outputs (boxes labeled (A)) are connected to key modules (shaded boxes) or exploration analysis modules (dashed arrows) through MIAAIM's Nextflow implementation (solid arrows). Algorithms specific to MIAAIM (boxes labeled (D)) are detailed in the corresponding figures (black bold text). Included are methods built in for application to single channel image data types and external software tools interfacing with MIAAIM (white boxes). [Figure 17-1] 17A and 17B illustrate HDIprep compression and HDIreg manifold alignment, respectively. HDIprep compression process may include (i) high-dimensional modality, (ii) subsampling, (iii) data manifold. Edge bundle connectivity of the manifold is shown on two axes of the resulting steady state embedding (*fractal-like structures may not reflect biologically relevant features). (iv) High connectivity landmarks identified by spectral clustering. (v) Landmarks are embedded in broad dimensions and exponential regression identifies steady state dimensions. The pixel locations are used to reconstruct the compressed image. HDIreg manifold alignment may include: (i) A spatial transformation is optimized to align the video to the fixed image. α-MI is calculated using KNN graph lengths between the resampled points (yellow). The edge length distribution panel shows the Shannon MI between the distribution of edge lengths in the graph at the resampling positions before and after alignment (α-MI converges to Shannon MI as α → 1). The MI value indicates the increase in the amount of information shared between the images after alignment. The KNN graph join shows the agreement between the modalities. (ii) The optimized transformation aligns the images. The results of the transformation from H&E images (green) to IMC (red) are shown. [Figure 17-2]FIG. 17C demonstrates an exemplary alignment: (i) full tissue MSI to H&E registration produces T0; (ii) H&E is transformed to the IMC full tissue reference to produce T1; (iii) ROI coordinates are extracted to the underlying MSI and IMC data in IMC reference space; (iv) H&E ROI is corrected in the IMC domain by transformation to produce T2; final alignment applies modality specific transformation; results are shown for IMC ROI. [Figure 18-1]Figures 18A-18J provide an overview of the performance of dimensionality reduction algorithms for summarizing mass spectrometry imaging data of diabetic foot ulcers. Figure 18A: Three mass spectrometry peaks highlighting tissue morphology were manually selected (top) and used to create an RGB image representation of the MSI data, which was converted to a grayscale image. The MSI grayscale image was then registered to its corresponding grayscale-converted hematoxylin-eosin (H&E)-stained section. The deformation field described by the determinant of its spatial Jacobian matrix (center) was saved for downstream use as a control registration. A 3D Euclidean embedding of the MSI data was then created using random initialization of each dimensionality reduction algorithm (bottom). These embeddings were then used to create RGB images following the procedure described above. A spatial transformation created by registering the manually identified peaks with the H&E image was then applied to the dimensionality-reduced grayscale image, and each was aligned to the grayscale H&E image. Figure 18B: Mutual information between each aligned grayscale embedding image (n=5 per method) and the grayscale H&E image was calculated using Parzen window histogram density estimation with a histogram bin width of 64. The plots are oriented so that the results match the concept of a "cost function" in an optimization situation, where the objective is to minimize the cost. Thus, more negative values ​​represent higher mutual information. UMAP consistently captures the multimodal information content for the H&E data. Figure 18C: Optimization of image registration between the grayscale version of the manually identified mass spectrometry peaks and the grayscale H&E image using mutual information as a cost function (Figure 18A, top) and external validation using Dice scores on seven manually annotated regions. The registration parameters used for the final registration used in Figure 18A are shown in dashed lines. Registration was performed by first aligning the images with a multiresolution affine registration (left).The converted greyscale version of the manually identified mass spectrometry peaks were then registered to the greyscale H&E image using nonlinear multiresolution registration. Figure 18D: Average neighborhood entropy (n=5) of each pixel calculated within a 10 pixel disk across dimensionality reduction algorithms. The results show that UMAP has the ability to enhance structures in tissue sections. Figure 18E: Manual annotation of greyscale H&E image used to verify registration quality with the control deformation field used for mutual information calculation in Figure 18B. Figure 18F: Cutout region with the same spatial coordinates as Figure 18E of the manually annotated region used to calculate Dice score in Figure 18C. The results show good spatial overlap across the different annotations. Figure 18G: Radar plot showing a comparison of the performance of dimensionality reduction algorithms - linear, nonlinear, local and global data structure preserving (t-SNE, UMAP, PHATE, Isomap, NMF, PCA) across a wide range of data representations. Shown are the average values ​​(n=5) of algorithm execution time (top, log-transformed), estimated steady-state manifold embedding dimension (right), noise tolerance (bottom), and multimodal mutual information with DFU MSI data (left). All plots are oriented such that larger values ​​represent better algorithm performance. The results show that UMAP has the ability to efficiently capture data complexity with a small number of degrees of freedom while balancing noise tolerance with the multimodal information content contained in histology images. Figure 18H: Intrinsic dimensionality of MSI data estimated by each dimensionality reduction method. Embedding errors (y-axis) are not comparable between plots. We plot the mean and standard deviation (n=5) embedding errors across embedding dimensions from 1 to 10. Convergence to the y-axis indicates that the ability of the algorithm to capture data complexity no longer improves with increasing dimensionality of the resulting embedding.Results show that the intrinsic dimensionality estimated by nonlinear methods (t-SNE, UMAP, PHATE, Isomap) is much less than that of linear methods (NMF, PCA), implying that fewer dimensions are needed to accurately represent the dataset. Figure 18I: Denoised manifold preserving (DeMaP) metric between the Euclidean distance of the resulting embedding corresponding to non-peak picked data and the geodesic distance in ambient space (no dimensionality reduction after peak picking) of the corresponding peak picked data. Results showing the mean and standard deviation DeMaP metric (Spearman's rho correlation coefficient) for all dimensionality reduction methods tested (n=5). The nonlinear methods Isomap, PHATE and UMAP all consistently preserve the manifold structure without pre-filtering of the data, with consistent correlations above 0.85 across 2-10 dimensions. Figure 18J: Computational run time per algorithm across embedding dimensions from 1 to 10. Mean and standard deviation (n=5) are plotted across each number of dimensions for each method. The nonlinear methods t-SNE and Isomap require longer run times than the nonlinear methods PHATE and UMAP. The linear methods require minimal run times; however, these methods cannot succinctly capture data complexity. [Figure 18-2] See description of Figure 18-1. [Figure 18-3] See description of Figure 18-1. [Figure 18-4] See description of Figure 18-1. [Figure 19-1]Figures 19A-19H provide an overview of the performance of dimensionality reduction algorithms for summarizing mass spectrometry imaging data of prostate cancer. Figure 19A: Same as Figure 18A, except for prostate cancer tissue biopsies. Figure 18B: Same as Figure 18B, except for prostate cancer tissue biopsies. Figure 19C: Optimization of image registration between the grayscale version of the manually identified mass spectrometry peaks and the grayscale H&E image using mutual information as a cost function (Figure 19A, top). The registration parameters used for the final registration used in Figure 19A are shown in dashed lines. Registration was performed by first aligning the images with a multi-resolution affine registration (left). The transformed grayscale version of the manually identified mass spectrometry peaks were then registered to the grayscale H&E image using a nonlinear multi-resolution registration. Figure 19D: Same as Figure 18D, except for prostate cancer tissue biopsies. Figure 19E: Same as Figure 18G, except for prostate cancer tissue biopsies. Figure 19F: Same as Figure 18H, except for prostate cancer tissue biopsies. Figure 19G: Same as Figure 18I, but for prostate cancer tissue biopsies. Nonlinear methods Isomap, PHATE, and UMAP all consistently preserve the manifold structure without pre-filtering of the data, with consistent correlations above 0.75 across 2-10 dimensions. Figure 19H: Results showing computational run times for each algorithm across embedding dimensions of 1-10. Means and standard deviations (n=5) are plotted across each number of dimensions for each method. Results show that nonlinear methods t-SNE, PHATE, and Isomap require longer run times than UMAP. Linear methods require minimal run times; however, these methods cannot succinctly capture data complexity and are not noise tolerant. [Figure 19-2] See description of Figure 19-1. [Figure 19-3] See description of Figure 19-1. [Figure 19-4] See description of Figure 19-1. [Figure 20-1]Figures 20A-20H provide an overview of the performance of dimensionality reduction algorithms for summarizing mass spectrometry imaging data of tonsils. Figure 20A: Same as Figure 18A except for tonsil tissue biopsies. Figure 20B: Same as Figure 18B except for tonsil tissue biopsies. Isomap and NMF consistently capture multimodal information content for H&E data. Figure 20C: Same as Figure 19C except for tonsil tissue biopsies. Figure 20D: Same as Figure 18D except for tonsil tissue biopsies. Figure 20E: Same as Figure 18G except for tonsil tissue biopsies. Figure 30F: Same as Figure 18H except for tonsil tissue biopsies. Figure 20G: Same as Figure 18I except for tonsil tissue biopsies. Figure 20H: Same as Figure 18J except for tonsil tissue biopsies. [Figure 20-2] See description of Figure 20-1. [Figure 20-3] See description of Figure 20-1. [Figure 20-4] See description of Figure 20-1. [Figure 21A] Figures 21A and 21B demonstrate that spectral-centered landmarks reproduce the steady-state manifold embedding dimension across tissue types and imaging techniques. Figure 21A: Sum-squared error of exponential regressions fit to steady-state embedding dimension selection from spectral landmarks compared to the full mass spectrometry imaging dataset across tissue types. The difference between exponential regressions fit to the cross entropy of landmark-centered embeddings and the full dataset embedding approaches zero as the number of landmarks increases. The dashed line indicates MIAAIM's default selection of 3,000 landmarks for computing the steady-state manifold embedding dimension. Figure 21B: Same as Figure 21A except pixels were subsampled in the imaging mass cytometry region of interest. [Figure 21B] See legend to Figure 21A. [Figure 22A]Figures 22A and 22B demonstrate that UMAP embedding of spatially subsampled imaging mass cytometry data with off-sample projection reproduces full data embedding (Figure 22B) while reducing run times in a diabetic foot ulcer sample (Figure 22A). [Figure 22B] See legend to Figure 22A. [Figure 23A] Figures 23A and 23B demonstrate that UMAP embedding of spatially subsampled imaging mass cytometry data with out-of-sample projection reproduces full data embedding (Figure 23B) while reducing run time in a prostate cancer sample (Figure 23A). [Figure 23B] See legend to Figure 23A. [Figure 24A] Figures 24A and 24B demonstrate that UMAP embedding of spatially subsampled imaging mass cytometry data with off-sample projection reproduces full data embedding (Figure 24B) while reducing run times in a tonsil sample (Figure 24A). [Figure 24B] See legend to Figure 24A. [Diagram 25] Figures 25A and 25B show the scaling of MIAAIM image compression to wide-field and high-resolution multiple image datasets by incorporating parametric UMAP. Figure 25A: Multiplexed CyCIF images of lung adenocarcinoma metastasis to lymph node (n = 100 million pixels, 0.65 μm / pixel resolution, 44 channels, 27 antibodies) and corresponding steady-state UMAP embedding and spatial reconstruction (showing 3 UMAP channels of 4-channel steady-state embedding). Parametric UMAP compresses millions of pixels and preserves tissue structure across multiple length scales. Figure 25B: Same as Figure 25A except for tonsil CyCIF data (n = 256 million pixels, 0.65 μm / pixel resolution). [Figure 26-1]Figures 26A-26I show that microenvironment correlation network analysis (MCNA) links protein expression to molecular distribution in the DFU niche. Figure 26A: MCNA UMAP of m / z peaks grouped into modules. Figure 26B: Exponentially weighted moving average of normalized ion intensity for top 5 positive and negative protein correlates. Color scheme indicates module assignment. Heatmap (right) shows Spearman's rho. Figure 26C: Exponentially weighted moving average of normalized mean ion intensity per module ordered away from the center of the DFU wound. Figure 26D: Raw IMC nuclei (Ir) and CD3 staining (left) in ROI (scale bar = 80 μm). Mask showing CD3 expression (center left). Aligned MSI showing one of the top CD3 correlates (center right). Overlap of CD3 expression with top molecular correlates (right). Figure 26E: Same as Figure 26D but with different ROI. Figure 26F: Unsupervised phenotyping. Shaded boxes indicate CD3+ populations. Heatmaps show normalized protein expression. Figure 26G: MCNA UMAP colored to reflect correlation of ions with Ki-67 in CD3+ and CD3- populations. Color scheme indicates Spearman's rho and point size indicates negative log-transformed Benjamini-Hochberg corrected P-value of correlation. Figure 26H: Tornado plot showing top 5 Ki-67 and CD3+ differential negativity and positivity correlations compared to CD3- cell population. X-axis shows CD3+ specific Ki-67 values. Color scheme of each bar indicates change in correlation from CD3- to CD3+ population. Figure 26I: Box plot showing ion intensity of top ions (top, positive; bottom, negative) differentially correlated with CD3+ specific Ki-67 expression across ROIs on DFU. Tissue maps of the top differentially associated CD3+ Ki-67 correlates (top, positive; bottom, negative) show ROIs on tissues containing CD3+ cells boxed (white). [Figure 26-2] See description of Figure 26-1. [Figure 26-3] See description of Figure 26-1. [Figure 26-4] See description of Figure 26-1. [Figure 27-1] Figures 27A-27H illustrate cobordism projection and domain transfer by (i-)PatchMAP. Figure 27A: Schematic showing PatchMAP stitching between bounded manifolds (reference and query data) to form cobordism (grey), information transfer across cobordism geodesics (top) and cobordism projection visualization (bottom). Figure 27B: Bounded manifold stitching simulation. PatchMAP projection (manually drawn dashed lines indicate stitching) and UMAP projection of integrated data are shown with NN values ​​that maximized SC for each method. Figure 27C: Data transfer from MSI to IMC by i-PatchMAP. Line plot shows Spearman's rho between predicted and true spatial autocorrelation values. Figure 27D: MSI to IMC data transfer benchmark. Figure 27E: CBMC multimodal CITE-seq data transfer benchmark. Figure 27F: PatchMAP of DFU single cells (blue) and DFU (red), tonsil (green) and prostate (orange) pixels based on MSI profile. Individual plots show IMC display for DFU single cells (right). Figure 27G: MSI to IMC data transfer from DFU single cells to whole tissue. Figure 27H: MSI to IMC data transfer from DFU single cells to tonsil tissue. [Figure 27-2] See description of Figure 27-1. [Figure 27-3] See description of Figure 27-1. [Figure 27-4] See description of Figure 27-1. [Figure 28]Figures 28A and 28B show that PatchMAP preserves bounded manifold structure while accurately embedding the relations between bounded manifolds in the cobordism. Figure 28A: PatchMAP embedding of the MNIST digit dataset (n=70,000) randomly split into two equally sized bounded manifolds. For smaller nearest neighbor values, PatchMAP is similar to the UMAP embedding, since open neighborhoods are preserved after intersecting pairwise nearest neighbor queries. Under these conditions, the intersection operation is similar to the fuzzy union performed by UMAP. For larger nearest neighbor values, PatchMAP captures the manifold relations in the cobordism while preserving the bounded manifold structure. Here, PatchMAP aligns the bounded manifolds along their principal axes to produce an approximate mirror image. This results in an equivalent split of the data in half, which is captured in the cobordism geodesic distance. Figure 28B: Validation against the full MNIST digit dataset of Figure 27B. Each digit in the dataset is considered to be a bounded manifold. When the nearest neighbor value is smaller, it is similar to the UMAP embedding, and when the nearest neighbor value is larger, PatchMAP is able to accurately model the cobordism geodesic distance. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0024] Detailed Description Overall, the present invention provides methods and computer-readable storage media for processing two or more spatially resolved data sets to identify cross-modal features, to identify diagnostics, prognoses, or theranostics for a disease state, or to identify trends in parameters of interest.

[0025] The term "theranostics" as used herein refers to diagnostic therapeutic methods. For example, theranostic approaches can be used for personalized therapy.

[0026] The method is designed as a general framework to collate spatially resolved datasets of widely diverse origins (e.g., laboratory samples, various imaging modalities, geographic information system data) together with other aligned data to identify cross-modal features that can be used as high-value or actionable indicators (e.g., biomarkers or prognostic images) consisting of one or more parameters that are uniquely revealed through the creation and analysis of multidimensional maps.

[0027] Advantageously, the methods described herein can be multiplexed to allow for the simultaneous interrogation of multiple (eg, at least 10 or at least 20) molecular analytes.

[0028] The method of the invention may be a method of providing a diagnosis, prognosis, or theranostics for a condition from three or more image modalities obtained from a biopsy sample from a subject, the method comprising comparing a plurality of cross-modal features to identify a correlation between at least one cross-modal feature parameter and the condition to identify a diagnosis, prognosis, or theranostics, wherein the plurality of cross-modal features comprises: (a) registering the three or more spatially resolved data sets to generate an aligned feature image, the aligned feature image comprising the three or more spatially resolved data sets that are spatially aligned; and (b) Extracting cross-modal features from the aligned feature images. wherein each cross-modal feature comprises a cross-modal feature parameter, and the three or more spatially resolved data sets are output by a corresponding imaging modality selected from the group consisting of three or more imaging modalities.

[0029] A method of the present invention may be a method of identifying cross-modal features from two or more spatially resolved datasets by the steps of: (a) registering the two or more spatially resolved datasets to generate an aligned feature image comprising the two or more spatially resolved datasets that are spatially aligned; and (b) extracting cross-modal features from the aligned feature image.

[0030] The method of the present invention can be a method of identifying a diagnosis, prognosis, or theranostics for a disease state from two or more imaging modalities. The method includes comparing a plurality of cross-modal features to identify a correlation between at least one cross-modal feature parameter and the disease state to identify a diagnosis, prognosis, or theranostics. The plurality of cross-modal features can be identified as described herein. In the method described herein, each cross-modal feature includes a cross-modal feature parameter. The two or more spatially resolved data sets are outputs by corresponding imaging modalities selected from the group consisting of two or more imaging modalities described herein.

[0031] A method of the present invention may be a method of identifying a trend of a parameter of interest within a plurality of aligned feature images identified by the methods described herein, the method comprising identifying a parameter of interest in the plurality of aligned feature images, and comparing the parameter of interest between the plurality of aligned feature images to identify a trend.

[0032] Figure 4 summarizes the required and optional steps for identifying cross-modal features. Step 1 is the spatial alignment of all modalities of interest. Steps 2-4 can be performed in parallel and are complementary approaches used to identify trends in expression / abundance of parameters of interest for modeling and predicting biological processes at multiple scales: cellular niches (fine-scale local context), local tissue heterogeneity (local population context), whole-tissue heterogeneity and trending features (global context), and disease / tissue states (combination of local and global tissue contexts).

[0033] With respect to data derived from biological samples relevant to biomedical and research applications, the method is expected to have broad applicability to data from a wide variety of tissue-based data acquisition techniques, including but not limited to RNAscope [1], multiplex ion beam imaging (MIBI) [2], cyclic immunofluorescence (CyCIF) [3], tissue-CyCIF [4], spatial transcriptome analysis [5], mass spectrometry imaging [6], codetection by index imaging (CODEX) [7], and imaging mass cytometry (IMC) [8].

[0034] The present invention also provides a computer readable storage medium. A computer program for identifying cross-modal features from two or more spatially resolved data sets as described herein can be stored in a computer readable storage medium, the computer program comprising a routine set of instructions for making a computer perform the steps of a method for identifying cross-modal features from two or more spatially resolved data sets. A computer program for identifying a diagnosis, prognosis, or theranostics of a disease state from two or more imaging modalities can be stored in a computer readable storage medium, the computer program comprising a routine set of instructions for making a computer perform the steps of the corresponding method described herein. A computer program for identifying a trend of a parameter of interest in a plurality of aligned feature images identified by the corresponding method described herein can be stored in a computer readable storage medium, the computer program comprising a routine set of instructions for making a computer perform the steps of the corresponding method described herein.

[0035] All computer-readable storage media described herein exclude any transitory media (e.g., volatile memory, carrier waves, e.g., data signals embodied in carrier waves in a network, e.g., the Internet). Examples of computer-readable storage media include non-volatile memory media, such as magnetic storage devices (e.g., traditional "hard drives", RAID arrays, floppy disks), optical storage devices (e.g., compact disks (CDs) or digital video disks (DVDs)), or integrated circuit devices, such as solid-state drives (SSDs) or USB flash drives.

[0036] Registration of spatially resolved data sets Integration of spatially resolved datasets (e.g., high-parameter spatially resolved datasets from various imaging modalities) poses challenges due to different spatial resolutions between modalities, possible presence of spatial transformations and misalignments, technical variations within modalities, and the uncertain existence of statistical relationships between different modalities, given the goal of discovering new relationships. Thus, the systems, methods, and computer-readable storage media disclosed herein provide a general approach for accurately integrating datasets from a wide variety of imaging modalities.

[0037] The method is demonstrated on an exemplary dataset designed for the integration of imaging mass cytometry (IMC), mass spectrometry imaging (MSI), and hematoxylin and eosin (H&E) datasets.

[0038] Image registration is often perceived as a fitting problem where a quality function is iteratively optimized by applying transformations to one or more images to spatially align them. In practice, image registration frameworks typically consist of successive pairwise or group-wise registration to a chosen reference image; the latter has been proposed as a way to register multiple images in a single optimization step and eliminate biases incurred by the choice of reference image and therefore reference modality [9,10]. Recently, both these frameworks have been extended to learning-based registration, which can handle large datasets using spatial transformer networks [11,12,13,14]. In our search for a suitable registration pipeline, we recognize the potential of group-wise registration schemes and learning-based models, respectively, especially in situations where tissue morphology varies significantly between adjacent sections (as with glandular prostate tissue) or where there is a large amount of data.

[0039] The methods disclosed herein are centered around a sequential pairwise registration scheme that can be guided and optimized at each step. Thus, the methods disclosed herein provide a platform for not only one-time image registration, but also for multi-sample registration in data sets across acquisition techniques and tissue types.

[0040] Image Registration Dimensionality reduction High-parameter datasets, often confounded by technical variations and noise, pose challenges regarding their analysis and mutual integration. The spatial integration of each modality currently requires that a representative image be presented that allows statistical correspondence with the other modalities in an image registration scheme. For the dataset under consideration, manual identification of such images quickly becomes difficult due to the number of parameters acquired and the complex relationships between these parameters.

[0041] The method of the present invention includes registering two or more spatially resolved datasets to generate a feature image comprising two or more spatially resolved datasets that are spatially aligned. Automatic definition of image features can be achieved using techniques that embed data into a space with a metric adapted to construct an entropy spanning graph. Such techniques include dimensionality reduction and compression techniques that embed high-dimensional data points (e.g., pixels) into Euclidean space. Non-limiting examples of dimensionality reduction techniques include Uniform Manifold Approximation and Projection (UMAP)

[15] , Isometric Mapping (Isomap)

[16] , t-SNE

[17] , potential of heat diffusion for affinity-based transition embedding (PHATE)

[18] , Principal Component Analysis (PCA)

[19] , Diffusion Map

[20] , and Nonnegative Matrix Factorization (NMF)

[21] , which are used to condense the dimensionality of data into a concise representation of the full set.

[0042] Uniform Manifold Approximation and Projection (UMAP) is a machine learning technique for dimensionality reduction. UMAP consists of a theoretical framework based on Riemannian geometry and algebraic topology. It results in a practical and scalable algorithm that can be applied to real-world data. The UMAP algorithm competes with t-SNE in terms of visualization quality and in some cases better preserves global data structures with superior run-time performance. Furthermore, UMAP has no computational constraints on the embedding dimension, which makes UMAP promising as a versatile machine learning dimensionality reduction technique.

[0043] Isomap is a nonlinear dimensionality reduction method used to compute a quasi-isometric low-dimensional embedding of a set of high-dimensional data points. The method allows to estimate the intrinsic geometry of the data manifold based on an approximation of each data point's neighbors on the manifold.

[0044] t-SNE is a machine learning algorithm for nonlinear dimensionality reduction that allows us to represent high-dimensional data in a low-dimensional space of 2D or 3D for better visualization. Specifically, it models each high-dimensional object by 2D or 3D points such that, with high probability, similar objects are modeled by nearby points and dissimilar objects are modeled by distant points.

[0045] PHATE (Potential of heat diffusion for affinity-based transition embedding) is an unsupervised low-dimensional embedding of high-dimensional data.

[0046] Principal component analysis (PCA) is a technique for the dimensionality reduction of large data sets by creating new uncorrelated variables that successively maximize the variance.

[0047] Diffusion mapping is a dimensionality reduction or feature extraction method that computes an embedding of a dataset into a (often low-dimensional) Euclidean space, whose coordinates can be computed from the eigenvectors and eigenvalues ​​of a diffusion operator on the data. The Euclidean distance between points in the embedded space is equivalent to the diffusion distance between the probability distributions centered on those points. Diffusion mapping is a nonlinear dimensionality reduction method that focuses on discovering the underlying manifold from which the data was sampled.

[0048] Nonnegative matrix factorization (NMF) is a dimensionality reduction method that decomposes a nonnegative matrix into the product of two nonnegative matrices.

[0049] This dimensionality reduction process is often data-dependent, and an adequate representation of the dataset requires observation of the performance of the chosen algorithm. In the exemplary dataset, our chosen method for dimensionality reduction is the Uniform Manifold Approximation and Projection (UMAP) algorithm

[17] . Our results (Figures 5, 6, 7A, 7B, and 7C) show that this algorithm, a manifold-based nonlinear technique, provides the best MSI data representation among the methods considered for the multimodal comparison with H&E, based on standard image registration and computational complexity tests, robustness to noise, and ability to capture information in the low-dimensional embedding. The dimensionality reduction process listed above can be applied to all datasets under consideration, but allows for manual curation of representative features of the modalities and is considered a "guided" dimensionality reduction.

[0050] To represent a compressed high-dimensional dataset as an image with a foreground and background, each pixel in the compressed high-dimensional image is considered as an n-dimensional vector, and the corresponding image is pixelated by referencing the spatial location of each pixel in the original dataset. This process results in an image with a number of channels equal to the embedding dimension. Dimensionality reduction algorithms typically compress the data into an n-dimensional Euclidean vector space, where n is the chosen embedding dimension. By definition, this space contains zero vectors, so pixels / data points are not guaranteed to be distinguishable from the image background (typically taken to be zero-valued). To avoid this, each channel is linearly rescaled to the range zero to one following the process of

[23] , allowing for the distinction between foreground (spatial locations containing acquired data) and background (spatial locations with no information).

[0051] Entering Landmarks The image registration step may include, for example, user-directed input of landmarks. User-directed input of landmarks is not a required step to complete image registration. Instead, this step may be included for quality improvement of results, for example, in cases where unsupervised automated image registration does not produce optimal results (e.g., different adjacent tissue sections, histological artifacts, etc.). In such cases, the methods described herein may include providing one or more user-defined landmarks. The user-defined landmarks may be input prior to optimization of registration parameters.

[0052] In certain preferred embodiments, user input is incorporated after dimensionality reduction. Alternatively, user input can be incorporated before dimensionality reduction by using the spatial coordinates of the raw data. In practice, user-defined landmarks can be stored within image visualization software (e.g., Image J, available at imagej.nih.gov).

[0053] Optimizing registration parameters Once features have been selected for modality registration by dimensionality reduction, the parameters of the alignment process can be optimized semi-automatically by hyperparameter grid search and also, for example, by manual validation. The computation for the registration procedure in the current implementation (separated from the dimensionality reduction step) can be performed, for example, in the open source Elastix software

[22] , which introduces a modular design to our framework. Thus, the pipeline can incorporate multiple registration parameters, cost functions (dissimilarity measures optimized during registration), and deformation models (transformations applied to pixels to align spatial locations from multiple images), allowing the alignment of images with any number of dimensions (from dimensionality reduction), the incorporation of manual landmark setting (for difficult registration problems), as well as the configuration of multiple transformations that allow fine-tuning and registration of datasets acquired with more than two imaging modalities (e.g., MSI, IMC, IHC, H&E, etc.).

[0054] Global spatial alignment optimization The image registration process may include optimizing a global spatial alignment of registration parameters. The optimization of the global spatial alignment may be performed on two or more data sets after reducing their dimensionality.

[0055] Using hyperparameter grid search, registration parameters can be optimized to ensure proper alignment of each modality at the whole-tissue scale, for example, for coarse-grained analysis (e.g., tissue-wide gradient calculation of markers of interest, marker / cell heterogeneity across tissue, identification of regions of interest (ROI) for further examination, etc.). In some embodiments, spatial alignment of the datasets can be performed in a propagative manner by registering full tissue sections (e.g., MSI, H&E, and toluidine blue stained images) for each dataset. The spatial coordinates for the ROI (e.g., IMC ROI taken from the toluidine blue stained image) can then be used to correct any local deformations that require further adjustment for fine-grained analysis (Figures 14 and 15).

[0056] In the exemplary dataset described herein, the spatial resolution of each modality was as follows: MSI approx. 50 μm, H&E approx. 0.2 μm, and IMC approx. 1 μm.

[0057] The methods described herein may preserve the spatial coordinates of high-dimensional high-resolution structure and tissue morphology. Thus, in some methods described herein, the higher-resolution ROI may remain unchanged at each step of the registration scheme (e.g., the exemplary registration scheme described herein). Such higher-resolution ROI may, for example, serve as the final reference image to which all other images are aligned. MSI data has been shown to reflect tissue morphology present in traditional tissue stains

[24] . Given this correspondence, coupled with the ability of tissue (H&E) stains to capture the spatial organization of cells, we choose to look at H&E images as a medium between the MSI and IMC datasets and as a linchpin for spatially aligning all modalities. Due to limitations in computational resources, a resolution of approximately 1.2 μm per pixel is used for H&E images in the registration process.

[0058] However, the use of a hierarchical multi-resolution registration scheme similar to the one implemented for our dataset also has the potential to register datasets of any resolution.

[0059] Optimizing local alignment for fine-grained spatial overlap The methods described herein may include a secondary fine-tuning of image alignment for smaller sized ROIs, which may be performed, for example, after all modalities have been aligned at the tissue level (global registration).

[0060] In the exemplary dataset described herein, the lack of morphological information currently available at the whole-tissue scale on the IMC images after acquisition, a result of the destructive nature of the IMC technique, necessitates this additional step of correcting for local deformations occurring within each ROI. To that end, single-cell multiplex imaging techniques that allow acquisition of full-tissue data, such as tissue-based cyclic immunofluorescence (t-CyCIF) [4] and coincidence detection by indexing (CODEX) [7], provide both large-scale coarse analysis of the heterogeneity of the specimen and local analysis of the ROI; however, the dilution of single-cell relationships resulting from that tissue-wide heterogeneity, combined with potential encounters with artifacts on the edges of the full-tissue specimen, often requires finer analysis of regions of interest (ROIs) within the full tissue. As a result, in full-tissue specimens, slides are often scanned for coarse morphological characteristics using low-magnification fields before obtaining finer analysis at the cellular level at higher magnifications.

[0061] In this light, our iterative full-tissue-to-ROI approach is generalizable to any multi-imaging technique, whether whole tissue or with predefined ROIs, as seen in our exemplary dataset. Our propagation registration pipeline allows for the correction of local deformations smaller than the grid spacing used in our hierarchical B-spline transformation model at the whole-tissue scale. It is well known that the number of degrees of freedom, and therefore the computational complexity and flexibility, of a deformation model increases with the resolution of the uniform control point grid spacing

[25] . The control point grid spacing of a nonlinear deformation model represents the spacing between the nodes that fix the deformation surface in the transformed image. When used with a multiresolution registration approach, the uniform control point spacing for the nonlinear deformation often scales with the image resolution. Thus, the coarse nonlinear deformation is corrected before the finer high-resolution registration at the local scale. While our pyramidal approach to full-tissue registration attempts to mitigate misalignment due to very fine or coarse grid spacing, we ultimately choose to ensure fine structural registration of each ROI by reducing the sampling space of the cost function from a global tissue-wide cost to one centered on each ROI after registering the full tissue.

[0062] The final registration proceeds by composing the resulting transformations in a propagation scheme after the steps of dimensionality reduction, global spatial alignment optimization, and local alignment optimization. The original data corresponding to each modality is then spatially aligned with all the others by applying the respective transformation arrays to each of the channels.

[0063] Manifold-based data clustering / annotation Once all modalities are spatially aligned through dimensionality reduction, analysis can proceed at the pixel level or at the level of spatially resolved objects (see Analysis of Predefined Spatially Resolved Objects). At the pixel level, data from each modality is aligned, but parsing from the amount of data present at the individual pixel level can be challenging (similar issues faced when selecting feature images for registration). Clustering is a method where similar data points (e.g., pixels, cells, etc.) are grouped together to reduce data complexity and preserve the overall data structure. Through this approach, individual pixels of an image can be grouped together to cluster homogenous regions of tissue, providing a more interpretable discretized version of the full image, reducing the complexity of the analysis from millions of individual pixels to a defined number of clusters (e.g., tens to hundreds). When used in conjunction with a heatmap or another form of data visualization, an overview of each cluster, or tissue region, can be visualized in a single image, aiding in the rapid interpretation of the profile of each region.

[0064] For the exemplary datasets described herein (FIGS. 7B and 7C), the UMAP algorithm proved robust to noisy (variable) features, and its computational efficiency allowed iterative partitioning of the data in a reasonable period of time. As a result of UMAP's robustness to noise and its ability to capture complexity, we found it to be most suitable for constructing mathematical representations of very high dimensional data, such as data obtained from MSI or similar methods, where hundreds to thousands of channels per image are available.

[0065] The dimensionality reduction part of the UMAP algorithm works by maximizing the information content contained in a low-dimensional graph representation of a dataset compared to its high-dimensional counterpart

[15] . In certain preferred embodiments, the dimensionality reduction optimization scheme can generalize over the high-dimensional graph itself. As a result, we extract the high-dimensional graph (simplicial set) and use it as input for community detection (clustering) methods (e.g., Leiden algorithm

[28] , Leuven algorithm

[29] , random walk graph partitioning

[34] , spectral clustering

[35] , affinity propagation

[36] , etc.) as opposed to clustering into the embedded data space itself as in

[30] . This graph-based approach can be applied to any algorithm that builds pairwise similarity matrices (e.g., UMAP

[15] , Isomap

[16] , PHATE

[18] , etc.). The method described herein performs clustering of the high-dimensional graph before the actual data dimensionality reduction (embedding) so that clusters are formed based on configurations that represent global manifold structures. The exemplary clustering approach used herein preserves global features of the data

[32] , in contrast to embeddings generated by local dimensionality reduction using methods such as t-SNE or UMAP (preferably t-SNE)

[18] . Compared to clustering approaches on graph structures obtained from reduced data spaces such as

[31] , the approach adopted on our exemplary dataset reduces the burden of identifying principal components from the raw data prior to clustering, but we found it to be susceptible to noise when using large or noisy datasets (e.g., the full MSI dataset from the image registration section above).

[0066] Regardless of the choice of clustering algorithm, the simplified representation of the data through the process allows the subsequent performance of numerous analyses ranging from predicting cluster assignments for unseen data, directly modeling cluster-cluster spatial interactions, to performing traditional intensity-based analyses independent of the spatial context. The choice of analysis depends on whether one is interested in features outside the scope of the study and / or task-spatial context at hand (cell type abundance in the data, heterogeneity of a given region, etc.) or whether one is focused on spatial interactions between objects (e.g. type-specific neighborhood interactions

[26] , higher-order spatial interactions - an extension of first-order interactions [7], predicting spatial niches

[27] ). The resulting analyses and predictions can then be used as diagnostic and prognostic features for disease, as well as indicators of biological processes of interest for purely scientific reasons.

[0067] Clustering allows data to be collated in an unsupervised manner. However, it is just as easy to manually annotate pixels on an image to identify a set of features that correspond to annotations of interest. In our UMAP embedded representation of an exemplary dataset from diabetic foot ulcer biopsy tissue, for example, two opposing extremes of tissue health can be easily identified. These tissue states can be combined after labeling to provide the same analysis listed above. In either case, the annotations and cluster identities serve as a set of discretized labels that can be further analyzed.

[0068] classification Classification algorithms can then be run after the image clustering or manual annotation portion to extend the cluster assignments to unseen data. These algorithms assign or predict data to groups based on their values ​​for the parameters used to build the classifier. "Hard" classifiers are algorithms that create defined margins between labels in a dataset; in contrast, "soft" classifiers form "fuzzy" boundaries between categories in a dataset that represent the conditional probability of class assignment based on parameter values ​​for the given data.

[0069] When using soft classifiers (e.g., conditional probabilities generated by random forests, neural nets with sigmoid final activation functions, etc.), further creation of a probability map of, for example, diseased / healthy tissue regions-diagnosis can be derived. This concept of probability map is best exemplified by the pixel classification workflow in the image analysis software Ilastik

[38] . After classification with a random forest classifier, relevant features can then be extracted and used to make predictions for better understanding. For example, the MSI parameters that had the greatest impact on the cluster conditional probabilities in our random forest classification were used to identify distinctive features between tissue regions.

[0070] On the other hand, hard classifiers allow for a clear assignment of classes to the data, and thus help impose constraints on when a clear category assignment (decision) is required. In our exemplary dataset, the MSI dataset was clustered at the pixel level using the UMAP-based method described above, and a random forest classifier was used to extend the cluster assignment to new pixels by assigning them to the maximum probability cluster (hard classification). This prescription was dictated by computational constraints and efficiency, as well as the ability to identify nonlinear decision boundaries generated by our manifold clustering scheme that are robust to parameter selection

[37] .

[0071] Analysis of predefined spatially resolved objects (cells, tissue structures, etc.) In tissue specimens, there are many objects of interest that are cells or other morphological features (e.g., blood vessels, nerves, extracellular matrix; or entire structures such as hair follicles or tumors). The spatial coordinates of these objects are important to identify in order to subsequently understand the imaging dataset at a level higher than the pixel level. In the exemplary datasets described herein, the IMC modality contains data at single-cell resolution, and the goal of the analysis is to connect this single-cell information to parameters of other modalities. In single-cell multiplex imaging analysis, computer vision and / or machine learning techniques can be applied to locate the coordinates of cells on the image and use those coordinates to extract aggregate pixel-level data, which can then be used to analyze the data at the single-cell level rather than the pixel level. This process is called "segmentation," and a wide variety of single-cell segmentation software and pipelines are available, such as Ilastik

[38] , Watershed Segmentation

[39] , UNet

[40] , and DeepCell

[41] . However, this segmentation process can be applied to any object of interest, and the coordinates obtained from the process can be used to aggregate data for the application of any of the above analyses (e.g., clustering, spatial analysis, etc.). Importantly, for our application, this segmentation allows for the aggregation of pixel-level data for every single cell, allowing for the clustering of cells regardless of spatial location. This process allows for the formation of cellular identities based on traditional surface or activation marker staining with the IMC modality alone. A similar approach is applicable to any object, under conditions where analysis and aggregation of pixel-level data is warranted.

[0072] Multimodal Data Feature Extraction and Analysis The methods described herein may include, for example, comparing data from different modalities for spatially resolved objects by using their spatial coordinates. The process of image registration spatially aligns all imaging modalities, so that an object can be defined in any one of the modalities used and still accurately maintain relevant features across all modalities. In the example described herein (FIG. 15), the IMC dataset was used to identify single cell coordinates, which were then used to extract features for single cells from both the aligned MSI pixel level data and the IMC pixel level data itself. The data was then clustered based on single cell measurements from the IMC modality alone and the MSI modality alone. The clustering of IMC single cell measurements can be used to determine cell types. This ability to integrate multiple imaging modalities allowed permutation tests to be performed for the enrichment or depletion of certain features in the MSI modality as a function of the corresponding cell types defined in the IMC dataset. Alternatively, the methods described herein can identify which IMC features are depleted or enriched based on cell types defined by MSI modalities. This type of cross-modal analysis can be extended to any number of parameters and any number of modalities. Permutation testing evaluates the randomized average value of each parameter for its observations regardless of modality, allowing for one-to-all comparisons, where the measurements evaluated are aggregated for a single modality by label. Also, as opposed to testing for enrichment or depletion using a statistical significance cutoff, one can ask how parameters from other modalities affect or correlate with the values ​​obtained in the current modality of interest. To address this question, correlation analysis can be performed across modalities, and models can be created that take multiple modalities into account. To this end, the tools listed above, such as random forest classifiers, can be used for the task of predictive modeling of objects based on multimodal portraits.Subsequent partitions of the classifier weights can then be extracted, as described above, to understand the relative influence of each parameter in each modality on the predictive task at hand.

[0073] Advantageously, the integration of these spatially resolved imaging data sets provides flexibility in analysis: analysis pipelines may be drawn from and used for many of the imaging modalities enumerated independently. Considering cross-modal analysis in this light reveals opportunities to validate exciting new multi-modal analysis techniques, in addition to proving their usefulness with new insights.

[0074] Additional potential applications Pixel-level calculations and analysis Most of the above analyses focus on either identifying spatially resolved objects or clustering pixel-level data to discretize the data set for review. Alternatively, if one wishes to analyze the registered images when they are at the pixel level, the trends of the parameters of interest can be collected across tissues or tissue regions of interest. For example, the gradient of the parameters of interest across the image can be visualized by calculating parameter density estimates. The resulting smooth representation of pixel-level data resembles a continuous gradient and can be visualized as a contour map or heat map. In our exemplary data set, we generated this visualization by calculating smooth versions of the markers of interest in the IMC data relative to each other, which shows that the overall trends of these parameters are relative to each other. This analysis is not limited to a single modality. As a result of the registration process and spatial alignment across modalities, gradients across modalities can also be calculated. These continuous representations, when formally implemented in spatial gradient models such as in

[49] , can also be used to provide numerical solutions to the effects of attractions and repulsions that intra- or inter-modality parameters have on each other. Combined with time-dependent analyses, these numerical solutions and equations offer the potential to develop tissue-level cross-modal simulation models. For example, given the data acquisition sensitivity required to identify single molecules with high confidence in MSI, our datasets can also be combined with cross-modal gradient relationships and known attractions and repulsions between single molecules to simulate biological processes in tissues.

[0075] Following the above discussion, spatial regression models are commonly used in geographic systems analysis [42,43] and can also be used to analyze relationships in multimodal biological tissue data at the pixel level as well as for spatially resolved objects. The utility of pixel-oriented analysis is best demonstrated in

[33] , where spatial variance component analysis is used to infer the contribution and effect of parameters calculated at the pixel level on cells (spatially resolved objects).

[0076] Multi-Domain Transformation Recently, there has been a continuous advancement in computer vision and artificial intelligence algorithms in both classification tasks and generative modeling. Of note are models capable of learning and generating the underlying distributions that create / represent the dataset at hand, such as generative adversarial networks

[44] and adversarial autoencoders

[45] . These models have the ability to predict and transfer knowledge gathered from one image / modality to the other. This concept of image-to-image, in our case, domain-to-domain translation, is best demonstrated by cycle-consistent generative adversarial networks

[46] . From this point, any modality can be translated from one to the other, provided that there exists a relationship between each other for training. This process, which we consider label-free, is an extension of the application of generative modeling in biological image prediction, such as that exemplified in [47,48], where we generate IMC images from a trained generative model on another modality.

[0077] The following examples are intended to illustrate the invention. They are not intended to limit the invention in any way. EXAMPLES

[0078] Example 1. Multimodal imaging and analysis of diabetic foot ulcer tissue Diabetic foot ulcer (DFU) biopsies, including the completely excised ulcer and surrounding healthy margin tissue, were performed, followed by tissue processing in preparation for multimodal imaging. Serial sections of DFU biopsies were imaged by matrix-assisted laser desorption / ionization (MALDI) mass spectrometry imaging (MSI), imaging mass cytometry (IMC), and light microscopy. Following multimodal imaging, the acquired high-dimensional data were processed using an integrated analysis pipeline to characterize molecular signatures (Figures 1 and 4). Each slice of the DFU biopsy was stained using hematoxylin and eosin (H&E) and imaged using bright-field microscopy scanning. To prepare DFU biopsy slices for MSI (Figure 2A), the slices were sprayed with matrix solution (optimized for each type of analyte). In this example, a matrix containing 40% 2,5-dihydroxybenzoic acid (DHB) in 50:50 v / v acetonitrile:0.1% TFA aqueous solution was used to preferentially image small molecules and lipids (Figures 2B and 2C). Imaging was performed using a Bruker Rapiflex™ MALDI-TOF mass spectrometry imaging system in positive ion mode, 10 kHz, 86% laser and 50 μm raster to obtain mass / charge (m / z) ratio spectra with peaks representing the molecular composition of the DFU biopsy slices (Figure 2D). Imaging mass cytometry was performed on regions of interest within the H&E stained and MSI imaged DFU biopsy slices. Following tissue or cell culture pretreatment, the samples were stained with metal-labeled antibodies (Figure 3). The labeled molecular markers in the samples were then ablated using an ultraviolet laser connected to the mass cytometer system (Figure 3). In the mass cytometer, the cells of the sample were vaporized, atomized, ionized, and filtered through a quadrupole ion filter. The isotope intensities were profiled using time-of-flight (TOF) mass spectrometry, and the atomic composition of each labeled marker in the sample was reconstructed and analyzed based on the isotope intensity profile (Figure 3).

[0079] Example 2. Multimodal and high-dimensional data processing and analysis Multimodal imaging data acquired using any combination of modalities including MSI, IMC, immunohistochemistry (IHC), H&E staining, etc. were processed using an integrated analysis pipeline (Figure 4). The analysis pipeline was designed as a generalizable framework to collate spatially resolved datasets of widely diverse origins (e.g., laboratory samples, various imaging modalities, geographic information system data) together with other aligned data to identify high-value or actionable indicators (e.g., biomarkers, or prognostic features) consisting of one or more parameters that are uniquely revealed through the creation and analysis of a multidimensional map. To create such a multidimensional map, a series of steps were employed to process the multimodal imaging data. First, spatial alignment of all modalities was performed in a process referred to as image registration (Figure 4). Steps 2-4 (2) image segmentation, (3) manifold-based clustering and annotation at pixel level, and (4) multimodal data feature extraction and analysis were performed in parallel and were complementary approaches used to identify trends in expression or abundance of parameters of interest for modeling and predicting biological processes at multiple scales: cellular niches (fine-scale local context), local tissue heterogeneity (local population context), whole-tissue heterogeneity and trending features (global context), and disease / tissue states (combination of local and global tissue contexts).

[0080] Example 3. Comparison of execution time and data dimension estimation using multiple dimensionality reduction methods Characterization of the runtime of multiple dimensionality reduction methods was performed on multimodal and high-dimensional imaging MSI datasets to develop a rapid and accurate method to (1) identify non-healing diabetic foot ulcers (DFUs) at presentation and (2) evaluate the efficacy of debridement procedures in DFU wound healing. The dimensionality reduction techniques uniform manifold approximation and projection (UMAP), isometric mapping (Isomap), t-SNE, potential of heat diffusion for affinity-based transition embedding (PHATE), principal component analysis (PCA), and non-negative matrix factorization (NMF) were used to condense the dimensionality of the MSI datasets (Figure 5). The intrinsic dimensionality of the MSI data was estimated by each dimensionality reduction method (Figure 5). The embedding error as the mean and standard deviation (n=5) was plotted as a function of 1-10 dimensions for all methods. Convergence to the embedding error values ​​indicated that increasing the dimensionality of the resulting embedding no longer improved the algorithm's ability to capture the complexity of the data. We observed that nonlinear dimensionality reduction methods, e.g., t-SNE, UMAP, PHATE and Isomap, converged to a much lower intrinsic dimension than linear methods, e.g., NMF and PCA, indicating that much fewer dimensions are needed to accurately represent the dataset. The computational run time for each algorithm was measured and plotted as the average run time and standard deviation across each number of dimensions (Figure 6). The nonlinear methods, t-SNE and Isomap, required longer run times than the nonlinear methods PHATE and UMAP. The linear methods required the least run time, but also failed to capture data complexity succinctly. The results showed that the UMAP algorithm, a manifold-based nonlinear technique, provided the best MSI data representation compared to other methods based on standard image registration and computational complexity tests, robustness to noise, and ability to capture information in low-dimensional embedding.

[0081] Example 4. Comparison of the mutual information captured by each dimension reduction method tested The mutual information between the 3D embedding of the grayscale version of the MSI data and the corresponding H&E stained tissue sections was characterized for nonlinear, e.g., t-SNE, UMAP, PHATE, and Isomap, as well as linear, e.g., NMF and PCA, dimensionality reduction methods (Figure 7). A standard of image registration for grayscale multimodal image alignment using mutual information as a cost function was implemented. Images resulting from each dimensionality reduction method were treated with an equivalent deformation field to facilitate spatial alignment with the same H&E image. The mutual information between the H&E grayscale image and each 3D embedding was then calculated. Mutual information is defined as greater than or equal to zero, where negative values ​​are consistent with the cost function being minimized in the registration process. The results showed that Isomap and UMAP consistently shared more information with the H&E grayscale image than the other tested methods. (Figures 7A, 7B, and 7C).

[0082] Example 5. Dimensionality reduction process pipeline Dimensionality reduction using UMAP was performed on the DFU biopsy MSI dataset (Figures 8-9). Each UMAP dimension in the 3D embedding was pseudocolored, for example, red for dimension U1, green for dimension U2, and blue for dimension U3 (Figure 9). Overlap of the three channels resulted in a composite grayscale image that was used for further analysis, including registration and feature extraction methods. Figure 8 illustrates this process, where the raw MSI m / z data (left panel) is subjected to dimensionality reduction using UMAP in this example to three dimensions (middle panel). The embedding dimensions can be assigned arbitrary colors to better visualize the projection of the data along the third dimension. After UMAP 3D embedding, each pixel of the dataset, now colored by a UMAP dimension, can be mapped back to its original location on the DFU image (right panel). This allows the visualization of arbitrary structures in the high dimensional dataset with respect to the collected tissue section.

[0083] Example 6. Relative evaluation of the robustness to noise of selected dimensionality reduction methods Linear dimensionality reduction methods, e.g., NMF and PCA, suffer from the problem of overestimating the intrinsic dimensionality of the data and are sensitive to noisy channels. We performed dimensionality reduction for linear and nonlinear methods and visualized the first two dimensions of the 4-dimensional embedding for each method (Figure 10). Linear methods require a larger number of features to capture the complexity of the dataset, and often the captured features are confounded by noise, with some features being dedicated only to representing noise. To further evaluate the robustness to noise of nonlinear, e.g., t-SNE, UMAP, PHATE and Isomap, and linear, e.g., NMF and PCA, dimensionality reduction methods, we characterized the manifold structure of full mass spectrometry imaging (MSI) data (noisy) and denoised MSI data (peak picked) using the denoised manifold preserving (DeMaP) metric. We calculated the DeMaP metric between the Euclidean distance of the resulting embedding corresponding to the noisy MSI data and the geodesic distance of the corresponding peak picked data. The mean and standard deviation DeMaP metrics for all dimensionality reduction methods tested were plotted across dimensions 1–10 (Figure 7C).

[0084] Example 7. Multi-scale image registration pipeline A multi-scale iterative registration approach was implemented to first spatially align the multimodal image dataset at the whole tissue level (referred to as global registration), followed by higher resolution registration at subset regions of interest (ROIs) (referred to as local registration). The spatial resolution of the imaging modalities varies widely between them. For example, MSI resolution is about 50 μm, H&E and toluidine blue resolution is about 0.2 μm, and IMC resolution is about 1.0 μm (FIG. 11). To preserve the spatial coordinates of high-dimensional high-resolution structures and tissue morphology during multimodal image registration, we maintained a high-resolution image that served as a reference image that remained unchanged at each step of the registration scheme, to which all other images were aligned.

[0085] Global grayscale image registration for DFU biopsy tissue imaged with MSI, H&E staining and toluidine blue staining was performed in a multi-step process using the Elastix registration toolkit (Figure 12). First, the MSI image was processed for dimensionality reduction using UMAP. The resulting MSI image after UMAP dimensionality reduction (referred to as MSI0) was registered to its corresponding H&E0 image to generate a transformed MSI1 image (Figures 12 and 13). This transformation (T1) warps the MSI image while keeping the H&E image fixed. This results in a transformed MSI image (MSI1), which is aligned to the H&E image. In parallel, the H&E0 image was registered to its corresponding toluidine blue0 image, a separate adjacent tissue section of the same DFU biopsy, which was used for IMC imaging. Toluidine blue0 contained the spatial coordinates of the IMC region of interest that served as reference coordinates for the subsequent local transformation of the images. This transformation (T2) warps the H&E image while keeping the toluidine blue image fixed. Finally, we apply the transformation T2 to the already transformed MSI1, resulting in an MSI image (MSI2) that is registered to toluidine blue 0. This process is summarized in two equations: MSI-fis the final transformed MSI image used for downstream analysis, T1 is the registration transformation from the MSI image to the H&E image, and T2 is the registration transformation from the H&E image to the toluidine blue (IMC) image. MSI-f = T2(T1); and T H&E-f is the final transformed H&E image used for downstream analysis, and T2 is the registration transformation from the H&E image to the toluidine blue (IMC) image, as described above. H&E-f = T2.

[0086] After spatially aligning the images from all modalities at a global level, we incorporated a secondary fine-tuning step of image alignment for smaller sized ROIs (Figure 13). As a result of the destructive nature of IMC imaging, there is a need to add spatial information about the imaged sample using a reference image of the same sample collected before IMC. The reference image is obtained from a toluidine blue stained image, providing the ability to correct for local deformations occurring in the tissue sample within each ROI. Fine-tuning of the registration at a local scale was performed by selecting a region of interest in the toluidine blue image corresponding to each MSI image. The overall registration for a single ROI is performed by appropriate (modality-dependent) successive transformations, first at a global level, followed by local transformations (Figure 14).

[0087] Example 8. Feature extraction and analysis of multimodal data Spatially aligned images from the multimodal dataset were analyzed to identify objects in a process called segmentation. Once the spatially resolved objects were identified, we set out to compare data from different modalities for that object by using their spatial coordinates. We compared features from registered images containing data from the IMC dataset (used to identify single cell coordinates due to its relatively high spatial resolution) and the MSI dataset (images A-C and A''-C'' in FIG. 15). We then clustered the data based on single cell measurements in the IMC modality alone and in the MSI modality alone. Clustering of IMC single cell measurements was used to determine cell type (images A'-C' and A''''-C'''' in FIG. 15). The ability to integrate multiple imaging modalities allowed us to perform permutation tests for enrichment or depletion of certain features in the MSI modality as a function of the corresponding cell type defined in the IMC dataset.

[0088] Example 9. Multi-omics Image Alignment and Analysis by Information Manifold (MIAAIM) MIAAIM is a sequential workflow that aims to provide a comprehensive portrait of tissue status. It includes four processing stages: (i) image preprocessing with the High-Dimensional Image Creation (HDIprep) workflow, (ii) image registration with the High-Dimensional Image Registration (HDIreg) workflow, (iii) tissue state transition modeling with Cobordism Approximation and Projection (PatchMAP), and (iv) cross-modality information transfer with i-PatchMAP (Figure 16). Image integration in MIAAIM starts with two or more assembled images (Level 2 data) or spatially resolved raster datasets (assembled images, Figure 16). The size and standardized format of the assembled images vary depending on the technique. For example, cycle fluorescence-based methods (e.g., CODEX, CyCIF) assemble BioFormats / OME-compatible 20-60-plex whole-tissue mosaic images after illumination correction (e.g., BaSiC) and tile stitching (e.g., ASHLAR); other methods acquire 20-100-plex data directly in regions of interest (ROIs) (e.g., MIBI, IMC). Further methods quantify thousands of parameters at rasterized locations on the full tissue or ROIs and are not stored in a BioFormats / OME-compatible format. For example, the imzML format, which relies on the mzML format used by the Human Proteome Organization, often stores MSI data.

[0089] Regardless of technique, assembled images contain a large number of non-uniformly distributed parameters, which makes comprehensive manually guided image alignment impossible. In addition, high-dimensional imaging generates a large feature space that challenges methods typically used in unsupervised settings. The HDIprep workflow in MIAAIM creates compressed images that preserve multiple salient features, enabling inter-technique statistical comparisons while minimizing computational complexity (HDIprep, Figure 16). For images acquired from tissue staining, HDIprep offers parallelized smoothing and morphological operations that can be applied consecutively for preprocessing. Image registration with HDIreg generates transformations to combine modalities in the same spatial domain (HDIreg, Figure 16). HDIreg uses Elastix, a parallelized image registration library for computing transformations, and is optimized to transform large multi-channel images with minimal memory usage while supporting tissue staining. HDIreg performs automatic image resizing, padding, and border cropping before applying image transformations.

[0090] The aligned data are well suited for established single-cell and spatial neighborhood analyses - they can be segmented to capture multimodal single-cell measurements (level 3 and 4 data) such as average protein expression or spatial features of cells, or analyzed at the pixel level. However, a common goal in pathology is to utilize synthetic tissue portraits to map the transition from health to disease. Similarities between systems - tissue state levels can be visualized with the PatchMAP workflow (PatchMAP, Figure 16). PatchMAP models tissue states as smooth manifolds that are stitched together to form a higher-order manifold called cobordism. It is a nested model that captures nonlinear intrasystem states and intersystem continuity as a result. This paradigm can be applied as a tissue-based atlas mapping tool to transfer information between modalities with i-PatchMAP (i-PatchMAP, Figure 16).

[0091] MIAAIM's workflow is nonparametric, using probability distributions supported by manifolds rather than training data models. MIAAIM is therefore technology-independent and generalizes to multiple imaging systems (Table 1). However, nonparametric image registration is often an iterative parameter tuning process rather than a "black box" solution. This poses significant challenges to reproducible data integration across institutions and computing architectures. Therefore, we have Docker containerized MIAAIM's data integration workflow and developed a Nextflow implementation to document human-in-the-loop processing and remove language-specific dependencies by following FAIR (Findable, Accessible, Interoperable, and Reusable) data stewardship principles.

[0092] [Table 1]

[0093] High-dimensional image compression with HDIprep. To compress high-parameter images, HDIprep performs dimensionality reduction on pixels using uniform manifold approximation and projection (UMAP) (Figure 17A). We performed a rigorous comparison using novel imaging datasets of diverse tissue biopsies with high complexity of cellular states including human DFU, prostate cancer and healthy tonsils acquired using MSI, IMC and hematoxylin-eosin (H&E). Based on dimensionality reduction benchmarks, UMAP consistently outperformed competing linear, nonlinear, global and local information-preserving algorithms in terms of its robustness to noise and ability to efficiently preserve data complexity while capturing morphological structure (Figures 18A-18J, 19A-19H, and 20A-20H).

[0094] HDIprep preserves global data complexity with the fewest degree of freedom requirements by finding a steady-state manifold embedding. To identify the steady-state dimensionality, we compute the information captured by the UMAP pixel embedding over a wide range of embedding dimensions (cross-entropy, Definition 1, Methods) and select the first dimension where the observed cross-entropy approaches the asymptote of the exponential regression fit. The computation of the steady-state embedding scales quadratically with the number of pixels, therefore HDIprep embeds the spectral landmarks in a pixel manifold that represents their global structure (Figures 21A and 21B).

[0095] Pixel-level dimensionality reduction is computationally expensive for large images, i.e., at high resolution (e.g., 1 μm / pixel). To reduce compression time while preserving quality, we developed a subsampling scheme that embeds a pixel subset of the spatial representation before spectral landmark selection and projects out-of-sample pixels into the embedding (Figures 22A, 22B, 23A, 23B, 24A, and 24B). HDIprep also combines all optimizations with a recent neural network UMAP implementation to scale to whole-tissue images. We demonstrate its effectiveness on publicly available 44-channel CyCIF images containing ~100 and ~256 million pixels (Figure 25). Thus, HDIprep presents an objective pixel-level compression method that is applicable to multiple modalities (Algorithm 1, Methods).

[0096] High-dimensional image registration (HDIreg). MIAAIM bridges the HDIprep and HDIreg workflows with a manifold alignment scheme parameterized by spatial transformations. We developed a theory for computing manifold α-entropy using entropy graphs for UMAP embeddings and applied it to image registration using entropy graph-based Renyi α-mutual information (α-MI) (HDIreg, Methods). HDIreg generates a transformation that maximizes the image-to-image (manifold-to-manifold) α-MI (Fig. 17B). This image similarity measure generalizes to Euclidean embeddings of any dimension by considering the distribution of k-nearest neighbor (KNN) graph lengths of the compressed pixels, rather than directly comparing the pixels themselves. Combining HDIprep compression with KNN α-MI extends intensity-based registration to complex images, even without the corresponding counterstaining between techniques.

[0097] Proof of principle 1: MIAAIM yields information on cell phenotype, molecular ion distribution and tissue state at all scales. To highlight the utility of high-dimensional image integration, we applied the HDIPrep and HDIreg workflows to MALDI-TOF MSI, H&E and IMC data from a DFU tissue biopsy containing a wide range of tissue states from the necrotic core of an ulcer to healthy margins. Image acquisition was performed on a 1.2 cm CT scanner for H&E and MSI data. 2 Molecular imaging with MSI enabled untargeted mapping of lipids and small metabolites in the 400-1000 m / z range across the entire specimen at a resolution of 50 μm / pixel. Tissue morphology was captured by H&E at 0.2 μm / pixel, and 27-plex IMC data were acquired from seven ROIs on adjacent sections at a resolution of 1 μm / pixel.

[0098] Cross-modality alignment was performed in a global-to-local fashion (FIG. 17C). We utilized HDIPrep compression for the high-parameter data and HDIreg manifold alignment for the registration of the compressed images. Destructive nature of IMC acquisition in small ROIs2 Due to this, we first aligned full tissue data from MSI, H&E (downsampled to approximately 3.5 μm / pixel) and IMC reference images. Local deformations that were not captured at the whole tissue scale within each ROI were corrected using manual landmark guidance. Serial sectioning deformations were taken into account by nonlinear transformations. Prior to nonlinear correction, registration was initialized by performing an affine transformation on the coarse alignment. Resolution differences were taken into account by a multiresolution smoothing scheme. The final alignment proceeded by configuring both modality and ROI-specific transformations.

[0099] After segmentation, quantification with the image processing software MCMICRO, and antibody staining quality control, the registered images yielded the following information for 7,114 cells: (i) the average expression of 14 proteins, including markers of lymphocytes, macrophages, fibroblasts, keratinocytes, and endothelial cells, as well as extracellular matrix proteins such as collagen and smooth muscle actin; (ii) morphological features, e.g., cell eccentricity, stereotacticity, extent, and area, spatial positioning of each cell center; and (iii) the distribution of 9,753 m / z MSI peaks across the complete tissue. We also quantified the distance from each MSI pixel and IMC ROI to the ulcer center, identified by manual inspection of H&E. Through the integration of these modalities, MIAAIM provided cross-modal information that could not be gathered by any single imaging system alone, such as single-cell protein expression and molecular abundance profiling of the microenvironment.

[0100] Proof of principle 2: Identification of molecular microenvironment niches correlated with cells and pathology through multiple-omics networking. We confirmed the presence of cross-modal associations from proof of principle 1 by performing microenvironment correlation network analysis (MCNA) on the registered IMC and MSI data (Figures 26A-26I). We performed community detection (i.e., clustering) on ​​MSI analytes (m / z peaks) based on single-cell protein measurements and their correlation with defined microenvironment correlation network modules (MCNMs; different colors in Figure 26A). Examination of MCNMs with top correlations with IMC-identified protein levels revealed that collections of molecules, rather than individual peaks, were associated with cellular protein expression (Figure 26B). MCNM clustered moderately positively correlated cellular markers indicative of inflammation and cell death (CD68, activated caspase-3) on a separate axis, and moderately positively correlated markers of immune regulation (CD163, CD4, FoxP3) and vasculature (CD31) on a separate axis. Some proteins, such as CD14 (a myeloid cell marker) and the cell proliferation marker Ki-67, did not correlate strongly with any m / z peak in all cells.

[0101] To gain insight into the association of molecular distribution with tissue health, we plotted the ion intensity distribution of MCNMs against their proximity to the ulcer epicenter (Figure 26C). This analysis revealed a shift in the molecular profile approximately 6 mm from the ulcer epicenter as the tissue condition progressed from healthy to damaged. We validated our observations and the performance of HDIreg to align microscale structures by visualizing the distribution of top correlated ions within the cellular microenvironment (Figures 26D and 26E).

[0102] The advantage of our analysis is the potential to identify molecular variations in one modality (here, MSI) that correlate with cell states identified using a different modality (here, IMC). We investigated whether m / z peaks are differentially associated with cell proliferation (Ki-67 marker in IMC). We performed unsupervised clustering on segmented cell-level expression patterns of IMC to identify cell phenotypes (Figure 26F), and performed differential correlation network analysis between phenotypes within well-separated CD3+ clusters (likely to identify infiltrating T cells at the wound site) and CD3- cell populations (Figure 26G). Interestingly, we found that correlations with Ki-67 expression shifted near significance (2σ) between CD3- and CD3+ populations for multiple m / z peaks (Fisher transformed one-sided z-statistics; Bonferroni corrected P-values) (Figure 26H).

[0103] We then took advantage of the spatial context preserved by MIAAIM to observe that the ion intensity of m / z peaks positively correlated with Ki-67 in CD3+ cells increased with distance from the wound, while molecules negatively correlated with Ki-67 specific to CD3+ cells showed the opposite trend (Figure 26I). This suggests that CD3+ T cell proliferation occurs primarily near the healthy margins of DFUs, confirming that molecular correlates of T-cell proliferation can be identified through this unbiased analysis. Collectively, these results provide insight into the molecular microenvironments associated with distinct functional and metabolic states of specific cell subtypes, as well as how these microenvironments are distributed across spatial contexts in a gradient from injured to healthy tissue.

[0104] Mapping tissue state transitions via cobordism approximation and projection (PatchMAP). To model transitions between tissue states, such as healthy or damaged, we generalize manifold learning and dimensionality reduction to higher-dimensional situations by developing a new algorithm called PatchMAP that integrates mutual neighborhood calculations with UMAP (Figure 27A and Algorithm 2, Methods). We hypothesize that the phase space of system-level transitions is nonlinear and can be parameterized consistently by manifold learning. PatchMAP therefore represents disjoints of manifolds (i.e., system states) as boundaries (i.e., state transitions) of a higher-dimensional manifold called a cobordism. Overlapping patches are connected by pairwise dependent neighborhood queries that represent geodesics of the cobordism between bounded manifolds, and stitched with t-norm, making their metrics compatible. Interpretation and implementation of PatchMAP embedding is similar to existing dimensionality reduction algorithms - similar data within or between bounded manifolds are located close to each other, while disparate data are further apart. PatchMAP incorporates both the topological structure of the bounded manifold and the continuity across the bounded manifold to generate a cobordism.

[0105] Currently, there is no method to form cobordism - the closest method to achieve this is the dataset integration algorithm from the single cell biology community. Therefore, to benchmark PatchMAP's manifold stitching, we compared PatchMAP with the data integration methods BBKNN, Seurat v3, and Scanorama in stitching simulations using the "digit" machine learning method development dataset (Figure 27B). We split the data by labels into bounded manifolds, and then applied each method to restitch the full dataset. In this task, full stitching resulted in a perfect separation of the projected bounded manifolds, which we quantified using the Silhouette Coefficient (SC). For visualization control, UMAP was used after data integration to provide similar embeddings to PatchMAP for all benchmark methods.

[0106] PatchMAP was robust to overlapping bounded manifolds and outperformed data integration methods at higher neighbor (NN) counts. All other methods mixed bounded manifolds incorrectly when there was no overlap, which is expected given that their assumptions were violated by the lack of manifold connectivity. On the other hand, PatchMAP stitching uses fuzzy intersection sets, which strongly weights accurate connections while pruning incorrectly connected data across manifolds. We also verified whether PatchMAP preserves bounded manifold structure while embedding higher-order structure between similar bounded manifolds (Figures 28A and 28B). At low NN values ​​and when bounded manifolds are similar, PatchMAP is similar to UMAP projection (Figures 28A and 28B). At higher NN values, manifold annotations are strongly weighted, which results in less mixing and better manifold separation.

[0107] Information transfer between imaging techniques and tissues (i-PatchMAP). We hypothesized that information transfer between biological states should account for successive transitions similarly and should be robust to the strength of manifold connections (including lack thereof). Therefore, the i-PatchMAP workflow uses PatchMAP as a method of paired domain transfer and quality control visualization to propagate information between different samples (information transfer, Fig. 27A). To do this, i-PatchMAP first normalizes the connections between bounded manifolds of the "reference" and "query" data to define local one-step Markov chain transition probabilities (transition probabilities, Fig. 27A), and then linearly interpolates measurements from the reference data to the query data (information propagation, Fig. 27A). Quality control of i-PatchMAP can be performed by visualizing the connections between bounded manifolds in the PatchMAP embedding (visualization of manifold connections, Fig. 27A).

[0108] To benchmark i-PatchMAP, we compare it with another non-parametric domain transfer tool, Seurat v3, and a modified version of UMAP (UMAP+) that incorporates transition probability-based interpolation similar to i-PatchMAP, on data from proof-of-principle 1 and a publicly available cord blood mononuclear cell (CBMC) CITE-seq dataset. 11We compared against i-PatchMAP. UMAP+ utilizes a directed NN graph from query data to reference data for data interpolation, rather than PatchMAP's metric-fit stitching. It therefore acts as a control for PatchMAP. We constructed 23 evaluation cases by tiling ROIs from proof-of-principle 1 and performed leave-one-out cross-validation with single-cell MSI profiles to predict IMC protein expression. We assessed accuracy using Spearman correlation between predicted spatial autocorrelation (Moran's I) and true autocorrelation per parameter, accounting for differences in resolution between imaging modalities. i-PatchMAP outperformed the tested methods in terms of its ability to transfer IMC measurements based on MSI profiles to the query data (Figure 27B), but all implemented methods consistently underperformed for parameters that do not have original spatial autocorrelation within the tile (TGF-β, FoxP3, CD163). For the CITE-seq dataset, we generated 15 evaluation cases and predicted the abundance of antibody-derived tags (ADTs) using single-cell RNA profiles. We quantified the performance using Pearson's correlation between true and predicted ADT values ​​(Figure 27C) and found that our i-PatchMAP implementation performed better or slightly better than other tested methods for all parameters.

[0109] Proof of principle 3: i-PatchMAP transfers multiplexed protein distributions between tissues based on molecular microenvironment profiles. To assess whether i-PatchMAP can be used to transfer molecular signature information between imaging modalities and even between different tissue samples, we used single-cell IMC / MSI protein measurements (see Proof of principle 1) to extrapolate IMC information based on MSI profiles to full DFU samples, as well as to separate prostate tumor and tonsil specimens. PatchMAP embedding of single cells and individual pixels into tissue-wide DFU ROIs based on MSI parameters revealed that the single-cell molecular microenvironment in the DFU ROI provided a good representation of the entire DFU molecular profile (Figure 27F). Thus, we used i-PatchMAP to transfer DFU single-cell protein measurements to full DFU tissues based on molecular similarity. i-PatchMAP predicted that the wound zone of DFU tissues would show high expression levels for CD68 (a marker of proinflammatory macrophages) and activated caspase-3 (a marker of apoptotic cell death). On the one hand, healthy margins of DFU biopsies were predicted to contain higher levels of CD4 (indicative of infiltrating T cells), and the cell proliferation marker Ki-67. Interestingly, PatchMAP visualization revealed that molecular microenvironments corresponding to certain single-cell measures (e.g., CD4) in DFUs had strong connections with MSI pixels in tonsil tissue (Figure 27F). In tonsil tissue, a lymphocyte-rich tissue, i-PatchMAP predictions for CD4 matched well with lymphocyte structures, and areas devoid of cellular content were accurately predicted not to contain CD4. On the other hand, there was no strong connection between the molecular profiles of prostate cancer samples and DFU biopsies. Thus, in the current dataset, the strong cellular and molecular connections between samples are likely supported by the common presence of certain immune cell populations. Indeed, IMC examination of the prostate biopsies used here showed poor immune cell infiltration.

[0110] method MIAAIM implementation. MIAAIM workflows are implemented in Python and connected via the Nextflow pipeline language, allowing automated results caching, dynamic processing resumption after changes in workflow parameters, and streamlining parallelized processing of multiple images. MIAAIM is also available as a Python package. Each data integration workflow is containerized to enable a reproducible environment and to eliminate any language-specific dependencies. MIAAIM output interfaces with many existing image analysis software tools (see Supplementary Note 1, Combining MIAAIM with existing bioimaging software). MIAAIM therefore complements existing tools rather than replacing them.

[0111] High-dimensional image compression and preprocessing (HDIprep). HDIprep is implemented by specifying successive processing steps. Options include image compression for high-parameter data, as well as filtering and morphological operations for single-channel images. Processed images were exported as 32-bit NIfTI-1 images using Python's NiBabel library. NIfTI-1 was chosen as the default file format for many of MIAAIM's operations due to its compatibility with Elastix, ImageJ for visualization, and its memory mapping capabilities in Python.

[0112] To compress high-parameter images, HDIprep identifies a steady-state embedding dimension for pixel-level data. Compression starts with arbitrary spatially guided subsampling to reduce the dataset size. We then implement UMAP to build a graph representing the data manifold and its underlying topological structure (FuzzySimplicialSet, Algorithm 1). UMAP aims to optimize the embedding of high-dimensional fuzzy simplicial sets (i.e., weighted undirected graphs) such that the fuzzy set cross-entropy between the embedded simplicial set and its high-dimensional counterpart is minimized, where the fuzzy set cross-entropy is defined as 35 .

[0113] Definition 1. Given a reference set A and membership functions u:A→[0,1] and v:A→[0,1], the fuzzy set cross entropy C of (A,u) and (A,v) is Defined as TIFF2024537615000003.tif11128.

[0114] Fuzzy set cross-entropy is a global measure of agreement between simplicial sets, aggregated among members of a reference set A (here, graph edges). Computation of its exact value scales quadratically with the number of data points, limiting its use to large data sets. Therefore, the current implementation of UMAP cannot compute the exact cross-entropy during its optimization of the low-dimensional embedding. Instead, it relies on stochastic edge sampling and negative sampling to reduce the running time for large data sets. 35 Together, to identify the steady-state embedding dimension, we compute patches on the data manifold that represent its global structure, and we use these to compute accurate cross-entropy after projecting these patches across a wide range of dimensions using UMAP. The result is a global estimate of the dimensionality required to accurately capture the manifold complexity.

[0115] To identify patches of global representation on the data manifold, we subjected the fuzzy simplicial sets to a variant of spectral clustering. We repeatedly project the spectral centers into a Euclidean space of increasing dimension using UMAP, calculate the fuzzy set cross entropy in each case, and then minimax normalize the obtained values. To identify the steady-state embedding dimension, we fit a least-squares exponential regression to the normalized cross entropy according to the dimension, and then simulate samples along the regression line to identify the first dimension that falls within the 95% confidence interval of the exponential asymptote. The subsampled data is embedded into the steady-state dimension, and the out-of-sample pixels are projected into this embedding using a native neighborhood-based method of UMAP (transform () function). Finally, all pixels are mapped back to their original spatial coordinates, and a compressed image is constructed whose number of channels is equal to the steady-state embedding dimension. These steps are summarized in the following pseudocode:

[0116] Algorithm 1: Image compression Input: multi-channel image (X), SVD dimension (b), k-means clusters (k), embedding dimension (n) Output: Compressed image (I) Function Compression TIFF2024537615000004.tif161141

[0117] Subsampling of image data. Subsampling is performed at the pixel level and is optional for image compression. Implemented options include a uniformly spaced grid in the (x,y) plane, random coordinate selection, and random selection starting from a uniformly spaced grid ("pseudo-random"). HDIprep also supports the specification of a mask for sampling regions, which can be useful for extremely large datasets.

[0118] By default, images with fewer than 50,000 pixels are not subsampled, images with 50,000-100,000 pixels are subsampled using 55% pseudorandom sampling starting from a uniformly spaced grid of 2x2 pixels, images with 100,000-150,000 pixels are subsampled using 15% pseudorandom sampling starting from a 3x3 pixel grid, and images with more than 150,000 pixels are subsampled on a 3x3 pixel grid. These default values ​​are based on empirical studies (Figures 22A, 22B, 23A, 23B, 24A, and 24B).

[0119] No subsampling was used for the presented MSI data. The subsampling ratios used for the presented IMC data were determined case by case from empirical studies and aligned with those used in the spectral landmark sampling experiments. Subsampling on a uniformly spaced grid of 10 × 10 pixels was used for CyCIF data compression.

[0120] Fuzzy simplicial set generation. To construct a pixel-level data manifold, we represent each pixel as a d-dimensional vector, where d is the number of channels in a given high-parameter image (i.e., discarding spatial information). We then implement the UMAP algorithm and extract a fuzzy simplicial set that represents the resulting manifold structure of these d-dimensional points. For all presented results, we used default UMAP parameters to generate this manifold: 15 nearest neighbors and a Euclidean metric.

[0121] Manifold landmark selection via spectral clustering. Spectral landmarks are identified using a variant of spectral clustering. We use randomized singular value decomposition (SVD) followed by mini-batch k-means to scale spectral clustering to large datasets following a procedure implemented in the PHATE (potential of heat diffusion for affinity-based transition embedding) algorithm. Given a symmetric adjacency matrix A that represents pairwise similarities between nodes (here, pixels) originating from TIFF2024537615000005.tif4128, we first compute the eigenvectors corresponding to the largest eigenvalue k of A. We then perform mini-batch k-means on the nodes of A using these eigenvectors k as features. Spectral landmarks are then defined as the d-dimensional centers of the resulting clusters.

[0122] By default, the input data is reduced to 100 components using randomized SVD and then partitioned into 3,000 clusters using mini-batch k-means. These default parameter values ​​are based on empirical studies (Figures 21A and 21B). Because steady-state embeddings of MSI and IMC data are only available after experimental testing, no landmark selection was used to process or determine the optimal embedding dimension for these datasets. Instead, full or subsampled datasets were used. All other steady-state embeddings for image data were compressed using the default parameters listed above.

[0123] Steady-state UMAP embedding dimension. By default, HDIprep embeds the spectral landmarks in a 1- to 10-dimensional Euclidean space to determine the steady-state embedding dimension. Exponential regression on the spectral landmark fuzzy set cross entropy is performed using a built-in function from the Scipy Python library. These default parameters were used for all presented data.

[0124] Tissue image preprocessing. HDIprep processing options for hematoxylin and eosin (H&E) stained tissue and other low channel tissue stains include image filters (e.g., median), thresholding (e.g., manually set or automatic), and sequential morphological operations (e.g., thresholding, opening and closing). Presented H&E and toluidine blue stained images were processed using a median filter to remove salt-and-pepper noise, followed by Otsu thresholding to create a binary mask representing the foreground. Successive morphological operations were then applied to the mask, including morphological opening to remove small connected foreground components, morphological closing to fill small holes in the foreground, and filling to close large holes in the foreground.

[0125] Image compression with UMAP parameterized by a neuronal network. We implemented parametric UMAP with default parameters and neural architecture in the TensorFlow backend. The default architecture consisted of a fully connected neuronal network with 3 layers and 100 neurons. Training was performed using gradient descent with a batch size of 1,000 edges and the Adam optimizer with a learning rate of 0.001.

[0126] High-Dimensional Image Registration (HDIreg). HDIreg is a containerized workflow that implements the open-source Elastix software together with a bespoke Python module to automate image resizing, padding and cropping that are often applied prior to registration. HDIreg incorporates several different registration parameters, cost functions and deformation models, and additionally allows manual definition of point correspondences for challenging problems, as well as configuration of transformations for fine-tuning (see Supplementary Material 2, notes on predicted performance of the HDIreg workflow).

[0127] High-parameter images are registered using a manifold alignment scheme parameterized by a spatial transformation that aims to maximize image similarity. Formally, we consider registration as a solution to the optimization problem 40 Think of it as such.

[0128] Domain Omega F A fixed d-dimensional image I F : TIFF2024537615000006.tif3128 and domain Ω m A moving q-dimensional image I M : Given TIFF2024537615000007.tif3128, we can calculate The aim is to optimize TIFF2024537615000008.tif7128, where T μ : TIFF2024537615000009.tif4128 is the parameter where S is the smooth transformation defined by the vectors in TIFF2024537615000010.tif4128, TIFF2024537615000011.tif4128 and I F is the similarity measure that is maximized when the

[0129] Differential Geometry and Manifold Learning: The MIAAIM manifold alignment scheme uses the entropy graph-based Renyi α-Mutual Information (α-MI) as the similarity measure S in Equation 1, which is extended to manifold representations of images embedded in Euclidean spaces (i.e., compressed images) with potentially different dimensions. This measure is justified in the HDIreg manifold alignment scheme through the notion of intrinsic manifold information (i.e., entropy). In the following, we introduce fundamental differential geometric concepts that can extend the existing foundations of intrinsic manifold entropy estimation to the UMAP algorithm.

[0130] Definition 2: Let X and Y be topological spaces. A function f:X→Y is each point TIFF2024537615000012.tif3128 and For each open neighborhood N in TIFF2024537615000013.tif4128, set f -1 (N) A function ∫:X→Y is a homeomorphism if it is continuous in one-to-one and upward correspondence and has a continuous inverse. If a homeomorphism exists between spaces X and Y, then they are called homeomorphism spaces.

[0131] Definition 3. A manifold M of dimension n (i.e., an n-manifold) is a second countable Hausdorff space such that each point in the n-dimensional Euclidean space TIFF2024537615000015.tif3128 has an open neighborhood of homeomorphism. Any open set For TIFF2024537615000016.tif3128, we can define the chart (φ,U), where: TIFF2024537615000017.tif4128 is a homeomorphism. (φ,U) can be said to act as a local coordinate system for M, and we If TIFF2024537615000018.tif3128 is nonempty, we define the transition between two charts (φ,U) and (ω,V) as It can be defined as TIFF2024537615000019.tif4128.

[0132] Definition 4. A smooth manifold is a manifold in which there exists a smooth transition map between each chart of M. The Riemannian metric g is a function of each point In TIFF2024537615000020.tif4128, the dot product g between the vectors tangent to M at y y ( . , . ) is a mapping in which the tangent vector of y is defined as T y We denote this as M. A Riemannian manifold, written (M,g), is a smooth manifold M with a Riemannian metric g. Given a Riemannian manifold, Riemann volume elements provide a means of incorporating volume functions in local coordinates. Given (M,g), we define a volume element ω in units of metric g, and the local coordinates of a point χ=χ1,…,χ n of TIFF2024537615000021.tif4128, where g(χ)>0 and ∧ indicate the wedge product. The volume of M under this volume formalism is Vol(M)=∫ M is given by ω.

[0133] Definition 5. An immersion of a smooth n-manifold M into N is a differentiable map Ψ:M→N such that dΨ p :T p M→T Ψ(p) N is all points is injective in TIFF2024537615000022.tif4128. Therefore, Ψ is an immersion if its derivative is injective everywhere.

[0134] Definition 6. An embedding between smooth manifolds M and N is a smooth function f:M→N, where f is an immersion and its continuous functions are embeddings of a topological space (i.e., injective homeomorphisms). A closed embedding between M and N is TIFF2024537615000023.tif4128 is the closed embedding.

[0135] (M, g) is the ambient Let be a compact n-dimensional Riemannian manifold embedded in TIFF2024537615000024.tif3128, Let TIFF2024537615000025.tif4128 be a set of independent identically distributed random vectors with values ​​drawn from the distribution supported by M. TIFF2024537615000026.tif4128 We define the neighborhood of the elements of TIFF2024537615000027.tif4128 as an open neighborhood. The goal of manifold learning is to find that the measure of distortion D is TIFF2024537615000028.tif4128 and The goal of the manifold learning problem is therefore to approximate the embedding f such that it is minimal between TIFF2024537615000029.tif5128 and TIFF2024537615000029.tif5128. It can be exported as TIFF2024537615000030.tif8128, where: TIFF2024537615000031.tif3128 is represents a family of possible measurable functions that take on the vector Open neighborhood for TIFF2024537615000033.tif4128 TIFF2024537615000034.tif4128 is often defined as the geodesic distance (or a probabilistic encoding thereof) approximated with a positive definite kernel, allowing the computation of the inner product in a Riemannian framework (as opposed to a pseudo-Riemannian framework, which is required to be non-positive definite). The distortion measures vary from algorithm to algorithm (see Supplementary Note 3, e.g. HDIPreprep dimensionality reduction validation). Our interest in this discussion is the measures induced by the embedded geodesics through the volume elements of the gluing of these coordinates. These provide the necessary components to quantify the intrinsic Renyi α-entropy of the embedded data manifold.

[0136] Entropy graph estimator. The Lebesgue density f, and Identical random vectors X1,…,X with values ​​in a compact subset of TIFF2024537615000035.tif3128 n Given, the external Renyi α-entropy of f is given by TIFF2024537615000036.tif10128, Where: The file is TIFF2024537615000037.tif4128.

[0137] Definition 7 (from Costa and Hero) 38 ). TIFF2024537615000038.tif4128 Let t2024537615000039.tif3128 be an equally distributed random vector with values ​​in a compact subset of t2024537615000039.tif3128, then under the Euclidean metric The nearest neighbor of TIFF2024537615000040.tif4128 is Given by TIFF2024537615000041.tif7128.

[0138] The k-nearest neighbor (KNN) graph is Represents the edges between TIFF2024537615000042.tif4128 and its k-nearest neighbors. Let TIFF2024537615000043.tif4128 be the set of k-nearest neighbors of TIFF2024537615000044.tif3128. Then, the total edge length of the KNN graph with respect to TIFF2024537615000045.tif4128 is given by TIFF2024537615000046.tif14128, where γ > 0 is the power weighting constant.

[0139] In fact, the extrinsic Renyi α-entropy of f can be appropriately approximated using a type of graph known as a continuous pseudo-additive graph that includes the k-nearest neighbor (KNN) Euclidean graph, and as the number of feature vectors increases, its edge length asymptotically converges to the Renyi α-entropy of the feature distribution. This property leads to the convergence of the KNN Euclidean edge length to the extrinsic Renyi α-entropy of the set of random vectors having values in a compact subset of TIFF2024537615000047.tif3128. This is a direct and inevitable consequence of the Beardwood-Halton-Hammersley theorem described below. 50 Let (M,g) be a compact Riemannian m-manifold embedded in TIFF2024537615000048.tif3128. Let TIFF2024537615000049.tif4128 be

[0140] Beardwood-Halton-Hammersley (BHH) theorem. Let (M,g) be an ambient compact Riemannian m-manifold embedded in TIFF2024537615000048.tif3128. Let TIFF2024537615000049.tif4128 be a set of i.i.d. random vectors having values in a compact subset of TIFF2024537615000050.tif3128, and consider the Lebesgue density to be f. Assume d ≥ 2, 1 ≤ γ < d, Define TIFF2024537615000051.tif6128. Then, with probability 1, TIFF2024537615000052.tif becomes 14128. The value that determines the right side of the limit line in Equation 6 is the external Renyi α - entropy given by Equation 4. The identically - distributed random vectors are ambient When restricted to a compact and smooth m - manifold M in TIFF2024537615000053.tif at 3128, the BHH theorem is generalized to the intrinsic Renyi α - entropy of the multivariate density f with respect to M defined by TIFF2024537615000054.tif at 10128 the estimation of TIFF2024537615000055.tif at 4128 becomes possible by incorporating the measure μ naturally induced by the Riemannian metric via the Riemannian volume element g This is formalized by the following given by Costa and Hero.

[0141] Theorem 1 (Costa and Hero): Let (M,g) be a compact Riemannian m - manifold embedded in ambient TIFF2024537615000056.tif at 3128. Let TIFF2024537615000057.tif at 4128 be an identically - distributed random vector on M with a bounded density f relative to the differential volume element μ induced by the metric g g Assume m≧2, 1≦γ<m, and Define TIFF2024537615000058.tif at 6128. Then, with probability 1, TIFF2024537615000059.tif at 22128 holds. Where β m,γ,k is a constant independent of f and (M,g). Similarly, the expected value TIFF2024537615000060.tif at 5128 converges to the same limit line.

[0142] The quantity that determines the limit when d'=m is the intrinsic Renyi alpha entropy of f given by Equation 7. Theorem 1 has been used together with the manifold learning algorithms Isomap and modified C-Isomap to estimate the intrinsic dimension of the embedded manifold. 39 In contrast to these results, which use all pairwise geodesic approximations to each point of the dataset to estimate the α-entropy, we aim to provide a similar formulation that exploits the local information contained in the data manifold, following the results of our dimensionality reduction benchmarks showing that algorithms that preserve local information are well suited to the task of high-dimensional image data compression (Figures 18A-18J, 19A-19H, and 20A-20H). Indeed, the information density of the volume of a contiguous region of the model family (i.e., the output embedding space or the set family of input points) is recognized in the definition of the information geometry of statistical manifold learning.

[0143] Entropy Graph Estimator of Local Information of Embedded Manifolds: In the following, we exploit two concepts to show that the intrinsic information of multivariate probability distributions supported by embedded manifolds in Euclidean space with the UMAP algorithm can be approximated using the BHH theorem: (i.) the compactness of the constructed manifold and (ii.) the preservation of Riemannian volume elements. We address (i.) with a straightforward proof, and we address (ii.) by providing a motivating example of the preservation of volume elements with UMAP.

[0144] Definition 8. A topological space X is compact if every open covering A of X contains a finite set of subsets that also cover X. An open covering means that the elements of A are open and that the union of the elements of A is equal to X: TIFF2024537615000061.tif3128.

[0145] Proposition 1. Let n>d and M be the ambient Consider a compact manifold of dimension γ (γ≦d) embedded in TIFF2024537615000062.tif3128. As a result, the projection f: The image f(M) of M below in TIFF2024537615000063.tif3128 is compact. Proof. Let (M,g) be the ambient Let TIFF2024537615000064.tif3128 be a compact Riemannian manifold (e.g., a manifold constructed by UMAP) with metric g, and let f be a function of M to Let f be a projection up to TIFF2024537615000065.tif3128. Since f is a projection, it is continuous and maps compact sets to compact sets.

[0146] Proposition 1 states that the d-dimensional Euclidean projection of a compact Riemannian manifold is It shows that if we map a value into a compact subset of TIFF2024537615000066.tif3128, this is a sufficient condition in the BHH theorem. The UMAP algorithm considers fuzzy simplicial sets, i.e., manifolds, constructed from finite extended pseudometric spaces (finite fuzzy realizing functors, see Definition 7). Finite means that these extended pseudometric spaces are constructed from a finite set of points. Given this finiteness condition, the compactness of a UMAP manifold follows naturally from Definition 8 - given an open covering on a manifold, we can find a finite subcover.

[0147] Therefore, the UMAP projection is compact according to Proposition 1. To extend the BHH theorem to the computation of the intrinsic α-entropy of a UMAP embedding as in Equation 7, we must show that the volume elements induced via the embedding are well approximated. Note that these results apply to any dimensionality reduction algorithm that can clearly preserve distances within an open neighborhood when embedding a compact manifold into Euclidean space. In what follows, we do not provide a proof that UMAP preserves distances within an open neighborhood around a point, although this would be the ideal scenario. Rather, we assume that this ideal scenario exists, and we explain how to find the optimal dimension to project the data to in order to satisfy this assumption.

[0148] In contrast to global data-preserving algorithms such as Isomap, which compute all pairwise geodesic distances or their approximations using a landmark-based approach, UMAP approximates the geodesic distances with an open neighborhood close to each point (see Lemma 2 below). Let Y be an identically distributed random vector with Lebesgue density μ, whose values ​​are constrained to lie on a compact, uniformly distributed Riemannian manifold M. 1,…, Y n Given, the sample Y drawn from μ i and Y j The geodesic curve between is probabilistically encoded using UMAP and a scaled exponential distribution: TIFF2024537615000067.tif22128, where ρ i is the vector Y i is the distance from to its nearest neighbor, and σ i is an adaptively selected normalization factor. Using the terminology of Equation 2, the embedding objective in UMAP is given by minimizing the fuzzy simplicial set cross entropy, which represents the distortion D (Definition 1). Formally, for a sample Y i and Y j A probability distribution P that encodes the geodesic between ij and sample f (Y i ) and f (Y j) a probability distribution Q that encodes the distance between ij Given TIFF2024537615000068.tif11128, where Q ij teeth, The embedded vector f (Y i ) and f (Y j ), where a, b are user-defined parameters for controlling the embedding spread.

[0149] Minimizing Equation 11 is generally not a convex optimization problem. The family from Equation 2 The optimization over TIFF2024537615000070.tif3128 is restricted to a subset rather than the complete family and therefore represents a local optimum at best. We include a larger family of measurable functions to more accurately approximate an optimal embedding of geodesic distances within an open neighborhood of a vector through identification of a steady-state embedding dimension, as outlined in the HDIprep workflow of a "quasi-global" optimization procedure.

[0150] Using the notation of Equation 2, the embedding The volume element induced by TIFF2024537615000071.tif4128 is similarly σ=g(χ)dχ which is TIFF2024537615000072.tif4128 By local coordinate χ 1,…, χ n minimizes the distortion between the volume elements of M given by Local coordinates by TIFF2024537615000073.tif4128 Given for TIFF2024537615000074.tif4128 The global optimal solution to Equation 7 (i.e., P ij =Q ij ) the volume elements are preserved: TIFF2024537615000076.tif5128. Moser proves that there exists a volume that preserves diffeomorphisms for the compact manifold under consideration. 53 .

[0151] To identify a manifold embedding that minimizes the distortion of the volume elements, we scale the increase in dimensionality with an exponential increase in the potential locations of points (i.e., copies of the real line). We observed a natural increase in the dimensionality of real-valued data by modeling it as a function of the size of the manifold. The relationship between the radius of an open neighborhood of a manifold in a learning environment and the volume and density was shown by Narayan et al. 52 is formalized in applications where density is preserved at a given embedding dimension by varying the open neighborhood radius; however, we can easily extend this to volume-preserving scenarios fitted under fixed radii to infer dimensions that satisfy the geodesic distance-preserving assumption. Consider the following example.

[0152] TIFF2024537615000078.tif4128, ambient Let M be the vector of the manifold embedded in TIFF2024537615000079.tif3128, and Y i The k-nearest neighbors of d Sphere In TIFF2024537615000080.tif4128, proportional volume Let us assume that the image is uniformly distributed in TIFF2024537615000081.tif4128. i The open neighborhood of TIFF2024537615000082.tif4128 with radius r mm-ball of a manifold N of TIFF2024537615000083.tif4128 and at the same time, the uniform distribution and the induced Riemann volume element V m We assume that we preserve the structure including 52 According to the above, the inventors have Using TIFF2024537615000084.tif4128, the local radius r of the embedding space m and Power law relationship between the original radius of TIFF2024537615000085.tif4128 It can be inferred that there is a TIFF2024537615000086.tif4128.

[0153] To extend this example, r m andr d (We consider that the radius is fixed in the original space, and the radius in the embedding space is Q in Eq. 11.) ij (We can assume that the ambient metrics in the embedding space and the original space are the same, and that these are the same as those in the native UMAP method.) TIFF2024537615000087.tif4128 and Suppose we want to produce geodesic distances in TIFF2024537615000088.tif4128 (Lemma 2 below). The ambient metric and radius are preserved, TIFF2024537615000089.tif4128, so this is TIFF2024537615000090.tif4128 and Geodesic distance δ between points in TIFF2024537615000091.tif4128 m and δ d Also, there is a power law relationship TIFF2024537615000092.tif4128. If we further assume TIFF2024537615000093.tif4128 (i.e., an ideal scenario in which the geodesic distances between points in an open neighborhood are preserved), we can solve the relationship between the geodesic curve and the dimension m. Specifically, Consider TIFF2024537615000094.tif4128. TIFF2024537615000095.tif4128 to δ m We obtain TIFF2024537615000096.tif37128, which means that the geodesics in an open neighborhood of a point in the original space are exponentially related to the dimension m.

[0154] Using the power-law relationship between the volumes of open neighborhoods in UMAP manifolds and their embedded counterparts, we can attempt to identify the dimension m such that the geodesics in the open neighborhoods are preserved by exponential regression. Note that geodesic distance preservation is stronger than volume preservation; however, volume preservation is also implied. As a result, the steady-state manifold embedding provides the Euclidean dimension for approximating the manifold geodesics of vectors sampled from M, and thus the volume elements of M. The KNN graph function computed in the steady-state embedding space provides the necessary machinery to compute the intrinsic α-entropy of the data manifold embedded in MIAAIM by applying the BHH theorem using the induced measure in the gluing of all coordinates as in Theorem 1.

[0155] However, note that given the assumptions introduced in our example, there is no guarantee that the distance outside the open neighborhood of a point is accurately modeled in the embedding space using UMAP. Therefore, applying Theorem 1 together with UMAP does not yield the same results as Costa and Hero, by replacing the KNN graph length with that obtained using the length function of a geodesic minimum spanning tree (GMST), another type of entropy graph. 39One should not expect to be able to reproduce the intrinsic entropy originally reported by. Our main contribution here is to combine the KNN intrinsic entropy estimator with a dimensionality reduction algorithm that preserves local information. We wish to extend these results to the setting where two such manifolds are compared, as a basis for our image registration application. An entropy graph-based estimator of α-MI in the image registration setting is described by 40 .

[0156] z(χ i )=[z i (χ i ), … ,z d (χ i )] at point χ i Let Z be a d-dimensional vector that encodes the features of f (χ)={z f (χ1), … ,z f (χ N )} is the feature set of the fixed image, and Z m (T μ (χ))={z m (T μ (χ1)), … ,z m (T μ (χ N ))} to T μ Let the feature set of the transformed video at (χ) points be TIFF2024537615000097.tif6128,χ i The result is a concatenation of the feature vectors of the fixed and transformed video in TIFF2024537615000098.tif20128 is a graph-based estimator for α-MI, where γ=d(1-α), 0<α<1, and the three graphs TIFF2024537615000099.tif52128 is the Euclidean graph function (length) of feature vector z over the k considered nearest neighbors to its pth nearest neighbor.

[0157] Renyi α-MI provides a quantitative measure of relatedness between the intrinsic structure of multiple manifold embeddings constructed with the UMAP algorithm. The Renyi α-MI measure is extended to feature spaces of arbitrary dimensionality, and is then used by MIAAIM in conjunction with its image compression methods to quantify the similarity between steady-state embeddings of image pixels in potentially different dimensions.

[0158] Proof-of-concept study. Since data acquisition was clear of tissue context with regions of interest on the IMC full-tissue reference image, we first aligned full-tissue sections and then used the coordinates of the IMC regions to extract data from all modalities for fine-tuning. A custom Python script was used to propagate the alignment across imaging modalities. Manual landmark correspondence was used when unsupervised alignment proved suboptimal. We accounted for alignment errors around the IMC region following full-tissue registration by padding regions with extra pixels before cropping. All registrations involving MSI or IMC data were performed using KNNα-MI with α=0.99 and 15 nearest neighbors as in Equation 12. All registrations aligning low channel slides (IMC reference toluidine blue image and H&E) were performed using histogram-based MI after grayscale conversion for rapid processing.

[0159] For full tissue images, we implemented a two-step registration process by first aligning the images using an affine model for the vector of parameters μ (Equation 1), followed by aligning the images with a nonlinear model parameterized by a B-spline. A hierarchical Gaussian smoothing pyramid was used to account for the resolution differences between image modalities, and stochastic gradient descent with random coordinate sampling was used for optimization. We further optimized the final control point grid spacing and number of hierarchical levels for the B-spline model to add to the pyramid smoothing to H&E alignment of each MSI data set separately (FIGS. 18A-18J, 19A-19H, and 20A-20H). A final control point spacing of 300 pixels for the nonlinear B-spline registration of MSI data to the corresponding H&E data struck a balance between correct alignment and unrealistic distortion, and then we identified spatial Jacobian matrices whose values ​​substantially deviated from 1 by eye and inspection. H&E and IMC reference tissue registration utilized a final grid spacing of 5 pixels. A similar optimization of the number of pyramid levels was performed on these data. All data that underwent image registration were exported and saved as 32-bit NIfTI-1 images. IMC data was not converted and was kept in 16-bit OME-TIF(F) format.

[0160] Cobordism Approximation and Projection (PatchMAP). PatchMAP is an algorithm that constructs smooth manifolds by gluing Riemannian manifolds to their boundaries, and projects higher-dimensional manifolds into lower-dimensional spaces for visualization. The higher-dimensional manifolds produced by PatchMAP can be understood as cobordisms, which are described by the following set of definitions:

[0161] Definition 9. Family of sets Given TIFF2024537615000100.tif4128 and a set of indices I, The disjoint union represented by TIFF2024537615000101.tif4128 is iFor an injective function φ i :S i →It is the set obtained by adding S. Disjoint union corresponds to the coproduct of sets.

[0162] Definition 10. Two closed n-manifolds M and N are given by A manifold is cobordant if its disjoint union, represented by TIFF2024537615000102.tif4128, is a boundary of some manifold W. We call a manifold W a cobordist. A boundary of an n-manifold is a boundary of the upper half-plane of the manifold W. We denote the boundary of W as ∂W.

[0163] PatchMAP approaches cobordism learning in a semi-supervised fashion where the data is assumed to follow a nonlinear cobordistic structure, and our task is to glue a lower-dimensional manifold to the boundary of a higher-dimensional manifold to generate a cobordism. Here, we want the coordinate transformation across the cobordism to have its own geometry independent of the metric of the bounded manifold. Indeed, this property allows us to explore data within the bounded manifold without relying on the specific geometry of the cobordism. Finally, cobordism geodesics are a fundamental building block for downstream applications such as the i-PatchMAP workflow. Furthermore, we want the cobordism to highlight bounded manifolds whose points overlap with high confidence - such overlaps can give rise to interesting nonlinearities in higher-dimensional spaces. A natural way to satisfy both of these conditions is to use the fuzzy set theoretical foundations of the UMAP algorithm.

[0164] The main goal of PatchMAP is then to identify smooth manifolds whose boundaries are disjoint unions of lower dimensional smooth manifolds and have a metric that is independent of the metric of each bounded manifold we choose to represent. We address this with a two-step algorithm that decouples the computation of bounded manifolds from cobordism. First, we compute the bounded manifolds by applying the UMAP algorithm to each dataset with a user-provided metric. In effect, this step results in a symmetric weighted graph that represents the geodesics within each bounded manifold. Our task is to find a finite set of n bounded manifolds F={(M1,g1) … (M n ,g n )} to the manifold (M F , g), where the metric g is a disjoint union. The inventors of the present invention have F Approximate the geodesic of For each element of TIFF2024537615000105.tif4128, the inner product g p Tangent space T with p M F We wish to approximate

[0165] Lemma 2 (McInnes and Healy 35 ). (M,g) is the ambient Let the Riemannian manifold in TIFF2024537615000106.tif3128 be Let TIFF2024537615000107.tif3128 be the point. g is an open neighborhood If p is locally constant around p in TIFF2024537615000108.tif3128, so that g is a constant diagonal matrix in the ambient coordinates, then Sphere centered on p with TIFF2024537615000109.tif9128 In TIFF2024537615000110.tif3128, The geodesic distance from p to any point in TIFF2024537615000111.tif3128 is TIFF2024537615000112.tif6128, where r is the ambient space is the radius of the sphere in TIFF2024537615000113.tif3128, TIFF2024537615000114.tif4128 is a metric relative to ambient space.

[0166] Assuming we can compare the data points across the bounded manifolds with the appropriate metric, we can use Lemma 2 for its disjoint union to compute the geodesics between points on the projection of each bounded manifold under a user-provided ambient metric. i and M. j Given, we have TIFF2024537615000115.tif4128 and The point is TIFF2024537615000116.tif4128 Calculate the pairwise geodesic distance between TIFF2024537615000117.tif4128, We can obtain the components needed to construct a metric for TIFF2024537615000118.tif4128. By extension, we can obtain the disjoint To TIFF2024537615000119.tif4128, A bounded manifold M with TIFF2024537615000120.tif4128 i and M. jThe concatenation of all pairwise distance calculations between and provides the ingredients for constructing geodesics on the complete cobordism over all pairwise combinations of bounded manifolds. However, as a result of using Lemma 2 to approximate the projection of manifold geodesics, we have an oriented, non-conforming view of geodesics to and from bounded manifolds over the entire cobordism. We can interpret these oriented geodesics as being defined on an oriented cobordism. We aim to encode the oriented geodesics in oriented cobordisms and the bounded manifold geodesics in a single data representation.

[0167] Our goal then is to construct an undirected cobordism, where the directed geodesics of the directed cobordism are decomposed into a single symmetric matrix representation. To this end, we can convert each of the above-mentioned extended pseudometric spaces into a fuzzy simplicial set using a fuzzy single-set functor (see Definition 9), which captures both the topological representation of the directed cobordism and the underlying metric information. The incompatibility in the directed cobordism geodesics can be resolved using our chosen norm. A natural choice for our chosen fuzzy set representation is the t-norm (also known as the fuzzy intersection). If we interpret the fuzzy simplicial set representation of the directed cobordism probabilistically, the intersection corresponds to a joint distribution of the directed cobordism metric space, which highlights directed cobordism geodesics occurring in both directions with high probability.

[0168] The final step is to integrate the bounded manifold geodesic with the symmetric cobordism geodesic obtained from the fuzzy intersection. We can do this by employing a fuzzy union (probabilistic t-conorm) over the genus of the extended pseudometric space, as in the original UMAP implementation. The result is a cobordism containing its own geometry captured in the cobordism geodesic, in addition to the individual bounded manifolds containing their own geometry.

[0169] Optimization of the low-dimensional representation of the cobordism can be achieved in a number of ways - we choose to optimize the embedding using fuzzy set cross-entropy (Definition 1) as in the original UMAP implementation for consistency. Note that since our algorithm produces symmetric matrices, PatchMAP can also be applied iteratively to construct "nested" cobordisms of hierarchical dimensions.

[0170] PatchMAP Implementation. To build the cobordism, PatchMAP first computes bounded manifolds by constructing fuzzy simplicial sets from each provided dataset, i.e., system state, by applying the UMAP algorithm (FuzzySimplicialSet, Algorithm 2). Then, pairwise directed nearest neighbor (NN) queries between bounded manifolds are computed in the ambient space of the cobordism (DirectedGeodesics, Algorithm 2). The directed NN queries between bounded manifolds are weighted by the native implementation of UMAP, for which the reader is referred to Equations 5 and 6. The resulting directed NN graphs between UMAP submanifolds are weighted, which reflect incompatible Riemannian metrics; i.e., they cannot simply be added or multiplied to integrate their weights. Therefore, we make the directed NN queries adaptable by stitching the cobordism metric and applying a fuzzy simplicial intersection, resulting in a weighted symmetric graph (FuzzyIntersection, Algorithm 2). The final cobordism generated by PatchMAP is obtained by taking a fuzzy union over the genus of all fuzzy simplicial sets (FuzzyUnion, Algorithm 2). To represent the connections between bounded manifolds in the PatchMAP cobordism projection, we implemented the Hammer edge bundling algorithm in the Datashader Python library. The pseudocode outlining the PatchMAP algorithm is given below.

[0171] Algorithm 2: PatchMAP Input: Dataset(D1,D 2… D n ), the bounded manifold ambient metric (g f ), cobordism ambient metric (g) Output: Cobordism (W) Function Stitch TIFF2024537615000121.tif64144

[0172] Domain / information transfer (i-PatchMAP). rq are the reference and query datasets, respectively. r and M. q Specifically, we define M rq is a matrix whose rows represent points in the reference bounded manifold and whose columns represent the nearest neighbors of the reference manifold points to the query factor manifold under a user-defined metric, and the i,j-th entry is TIFF2024537615000122.tif4128 and Points such as TIFF2024537615000123.tif4128 We multiply the transferred feature matrix F by the prediction P q where the weight matrix W rq The transpose of M rq of TIFF2024537615000125.tif4128 Normalization gives: TIFF2024537615000126.tif5128. In this situation, the matrix W rq is the state p derived from the geodesic distance on the cobordism. i andp j can be interpreted as a single-step transition matrix of the Markov chain between

[0173] Biological Methods. All patient tissue samples were obtained with approval from the Institutional Review Boards (IRB) of Massachusetts General Hospital (protocol #2005P000774) and Beth Israel Deaconess Medical Center (protocol #2018P000581).

[0174] Generation of imaging mass cytometry data. Frozen tissues were serially sectioned at 10 μm thickness using a Microm HM550 cryostat (Thermo Scientific) and thaw mounted onto SuperFrost™ Plus Gold charged microscopy slides (Fisher Scientific). After temperature equilibration to room temperature, tissue sections were fixed in 4% paraformaldehyde (Ted Pella) for 10 min and then rinsed three times with cytometry grade phosphate buffered saline (PBS) (Fluidigm). Nonspecific binding sites were blocked with 5% bovine serum albumin (BSA) (Sigma Aldrich) in PBS containing 0.3% Triton X-100 (Thermo Scientific) for 1 h at room temperature. Appropriate titrated concentrations of metal-conjugated primary antibodies (Fluidigm) were mixed with 0.5% BSA in DPBS and applied overnight at 4° C. in a humidified chamber. Sections were then washed twice with PBS containing 0.1% Triton X-100 and counterstained with iridium (Ir) intercalator (Fluidigm) at 1:400 in PBS for 30 min at room temperature. Slides were rinsed in cytometry grade water (Fluidigm) for 5 min and air-dried. Data acquisition was performed using a Hyperion imaging system (Fluidigm) and CyTOF software (Fluidigm) with 33 channels, a frequency of 200 pixels / sec, and a spatial resolution of 1 μm. Images were visualized with MCD Viewer software (Fluidigm) and then data were exported as text files for further analysis. After imaging, slides were quickly stained with 0.1% toluidine blue solution (Electron Microscopy Sciences) to reveal gross morphology. Slides were digitized using a digital camera at a resolution of approximately 2.75 μm / pixel.

[0175] Generation of mass spectrometry imaging data. Pairs of 10-μm-thick sections from the same tissue blocks used for imaging mass cytometry were melt-mounted onto indium tin oxide (ITO)-coated glass slides (Bruker Daltonics). Tissue sections were coated with 2.5-dihydroxybenzoic acid (40 mg / mL in 50:50 acetonitrile:water containing 0.1% TFA) with an automated matrix applicator (TM-Sprayer, HTX Imaging). Mass spectrometry imaging of the sections was performed using a rapifleX MALDI Tissuetyper (Bruker Daltonics, Billerica, MA). Data acquisition was performed using FlexControl software (Bruker Daltonics, version 4.0) with the following parameters: positive ion polarity, molecular weight scan range (m / z) 300–1000, 1.25 GHz digitizer, 50 μm spatial resolution, 100 shots per pixel, and 10 kHz laser frequency. Regions of interest for data acquisition were defined using FlexImaging software (Bruker Daltonics, version 5.0), and individual images were visualized using both FlexImaging and SCiLS Lab (Bruker Daltonics). After data acquisition, sections were washed with PBS and subjected to standard hematoxylin-eosin histological staining, followed by dehydration in graded concentrations of alcohol and xylene. Stained tissues were digitized at a resolution of 0.5 μm / pixel using an Aperio ScanScope XT brightfield scanner (Leica Biosystems).

[0176] Mass spectrometry imaging data preprocessing. Data were processed in SCiLS LAB 2018 using total ion number normalization to the average spectrum and peak centroiding with interval width of ±25 mDa. For all analyses, a peak range of m / z 400–1,000 was used after peak centroiding, which resulted in 9,753 m / z peaks. No peak picking was performed on the presented data unless explicitly mentioned. Data were exported from SCiLS Lab as imzML files for further analysis and processing.

[0177] Single-cell segmentation. To quantify single-cell parameters in the IMC and registered MSI data in the DFU dataset, we performed cell segmentation on the IMC ROIs using the pixel classification module in Ilastik (version 1.3.2)

[38] , which utilizes a random forest classifier for semantic segmentation. For each ROI, two 250 μm × 250 μm regions were cropped from the IMC data and exported in HDF5 format for use in supervised learning. To ensure that each cropped region was a representative learning sample, a respective global threshold was created using Otsu thresholding for iridium (nuclear) staining in the Scikit-image Python library. The cropped regions were required to contain more than 30% of the pixels for their respective thresholds.

[0178] Training regions were annotated for "background", "membrane", "nucleus" and "noise". Random forest classification incorporated Gaussian smoothing features, edge features (including Gaussian Laplacian, Gaussian gradient magnitude and difference of Gaussian features) and texture features (including structure tensor eigenvalues ​​and Gaussian Hessian eigenvalues). The trained classifier was used to predict the assignment probability of each pixel to the four classes in the full image and the predictions were exported as 16-bit TIFF stacks. To remove artifacts in the cell staining, the noise prediction channel was Gaussian blurred with sigma 2 and Otsu thresholding with a correction factor of 1.3 was applied, from which a binary mask was created that separates the foreground (high pixel probability is noise) from the background (low pixel probability is noise). Using the noise mask, a zero value in the other three probability channels from Ilastik (nucleus, membrane, background) was assigned to all pixels considered as foreground in the noise channel. Denoised three-channel probability images of nuclei, membranes and background were used for single-cell segmentation within CellProfiler (version 3.1.8)

[59] .

[0179] Single-cell parameter quantification. Single-cell parameter quantification was performed on IMC and MSI data using an in-house modified quantification (MCQuant) module within Multiple Choice Microscopy Software (MCMICRO)

[60] , which receives the NIfFTI-1 files after cell segmentation. IMC single-cell measurements were transformed using 99th percentile quantile normalization before downstream analysis.

[0180] Imaging Mass Cytometry Cluster Analysis. Cluster analysis was performed in Python using the Leiden community detection algorithm with the leidenalg Python package. A simplex set of UMAP (weighted undirected graph) created with 15 nearest neighbors and Euclidean metrics was used as input for community detection.

[0181] Microenvironment correlation network analysis. To calculate the association between MSI and IMC modalities, we used Spearman's correlation coefficient in Python Scipy library. M / z peaks from MSI data that were not correlated with IMC data with Bonferroni corrected P value above 0.001 were removed from the analysis. Correlation modules were formed with hierarchical Louvain community detection using Scikit network package. Resolution parameters used for community detection were selected based on the elbow point of the graph plotting resolution vs. modularity of the community detection results. Simplex sets of UMAPs created with 5 nearest neighbor and Euclidean metrics were used as inputs for community detection after inverse cosine transformation of Spearman's correlation coefficient to form the metric distance. Visualization of MSI correlation module trends to IMC parameters was calculated using exponentially weighted moving averages in Python's Pandas library after standard scaling of IMC and MSI single-cell data. For plotting, MSI moving averages were additionally minimax scaled to the range 0-1. The differential correlation of the variables u from the MSI data and v from the IMC data between conditions a and b was quantified and expressed as: TIFF2024537615000127.tif5128 is used to rank the correlation coefficients, where the change in correlation coefficient for each pair u,v between conditions is weighted by the maximum absolute correlation coefficient between both conditions. TIFF2024537615000128.tif4128Significance was calculated using one-sided Bonferroni corrected z-statistics after Fisher transformation.

[0182] Benchmarking dimensionality reduction algorithms. We outline the method we use to benchmark dimensionality reduction algorithms in HDIPrep dimensionality reduction validation in Supplementary Note 3.

[0183] Benchmarking spatial subsampling. The default subsampling parameters in MIAAIM are based on experiments across IMC data from DFU, tonsil and prostate cancer tissues, recording the Procrustes transformed sum-squared error between the subsampled UMAP embedding, the subsequent projection of out-of-sample pixels, and the full UMAP embedding using all pixels. We performed benchmarking of spatial subsampling over a wide range of subsampling ratios.

[0184] Benchmarking of spectral landmarks. Subsampling rates and dimensions for validating landmark-based steady-state UMAP dimensions were case-by-case determined from empirical studies and aligned with those used for the presented registration data. Parameters were selected due to the computational burden of calculating cross entropy on large data. To compare landmark steady-state dimension selection with subsampled data, we compared the shape of the exponential regression fits from both data sets using the sum of squares error. The sum of squares error was calculated over a wide range of resulting landmarks.

[0185] Submanifold stitching simulations. Simulations were performed with the MNIST digit dataset in the Python Scikit-learn library with default parameters for BKNN, Seurat v3, Scanorama, and PatchMAP across a wide range of nearest neighbor values. Data points were partitioned by digit label and stitched together using each method. The combined data from each test method except PatchMAP was then visualized in UMAP. The quality of the submanifold stitching for each algorithm was quantified using the Silhouette Coefficient, implemented in Python with the Scikit-learn library, in the UMAP embedding space. The Silhouette Coefficient is a dispersion measure for the partition of the dataset. High values ​​indicate that data from the same label / type are grouped together tightly, while low values ​​indicate that data from different types are grouped together. The Silhouette Coefficient (SC) is the average silhouette score s calculated across each data point in the dataset, as follows: TIFF2024537615000129.tif22128, where a(i) is the average distance from data point i to all points with the label, and b(i) is the average distance from point i to all other data that do not have the same label.

[0186] Transfer of CBMC CITE-seq data. CBMC CITE-seq data were preprocessed by the vignette provided by Satija lab at https: / / satijalab.org / seurat / articles / multimodal_vignette.html. RNA profiles were log-transformed and ADT abundances were normalized using centered log-ratio transformation. Various RNA features were then identified and principal component analysis was used to reduce the dimensionality of the RNA profiles of cells. The first 30 principal components of single-cell RNA profiles were used to predict ADT abundances in single cells. The CBMC dataset was randomly split into 15 evaluation cases with 75% training data and 25% test data. The training data was used to predict the test data measures. The prediction quality was quantified using Pearson's correlation coefficient between true and predicted ADT abundances. Correlations were calculated using Python library (Scipy). After checking the transfer anchors using default parameters, Seurat was implemented with the TransferData function (FindTransferAnchors function). PatchMAP and UMAP+ were applied with 80 nearest neighbors and Euclidean metrics in PCA space.

[0187] Transfer of spatially resolved image data. To benchmark information transfer from MSI to IMC, we performed leave-one-out cross-validation using segmented single cells from 23 image tiles (each with ~100 to ~500 cells) from the DFU dataset. The IMC ROI was divided into four equally sized quadrants to create 24 tiles. One tile was removed due to lack of cell content. Prior to information transfer, data were transformed using principal component analysis with 15 components using the Scikit-learn library. After checking the transfer anchors with default parameters and 15 principal components, Seurat was implemented using the TransferData function (FindTransferAnchors function). PatchMAP and UMAP+ were implemented using 80 nearest neighbors and Euclidean metrics in PCA space. Information transfer quality was calculated for each predicted IMC parameter by calculating Pearson's correlation between the Moran's autocorrelation index of the ground truth data and the predicted data in Python using the Scipy library. Moran's autocorrelation index (I) is: Given by TIFF2024537615000130.tif11128 13 where N is the number of spatial dimensions in the data (for our purposes, 2), χ is the abundance of the protein of interest, TIFF2024537615000131.tif3128 is the average abundance of protein χ, w ij is the spatial weight matrix, and W is the weight matrix for all w ij It is the sum of.

[0188] Additional Notes Supplementary Note 1. Combination of MIAAIM and existing bioimaging analysis software MIAAIM's core functionality allows for cross-technique and cross-tissue comparisons. As shown in our proof-of-principle example, its functionality has broad applications that can be configured and performed with other software applications. We anticipate that a challenge for many users will be performing successive registrations and analyses between various software. MIAAIM's multiple output data formats directly interface with numerous tools for visualization, cell segmentation and single-cell analysis (Table 2), paving the way for continued investigation of multimodal tissue portraits in a variety of contexts.

[0189] [Table 2]

[0190] Supplementary Note 2. A note on expected performance of the HDIreg workflow A fundamental assumption in intensity-based image registration is that a quantifiable relationship exists between modalities - this is often met in practice, as shown in our proof-of-principle application. However, this assumption can be undermined by artifacts such as folding, tearing, and, in the case of serial sectioning, nonlinear deformation. In our experience, glandular tissues, e.g., those derived from the prostate, are likely to exhibit high architectural variability over short distances, making alignment of images from separate sections difficult. Manual landmark guidance can be used in challenging use cases such as those created by serial tissue sectioning. By using the Elastix library, HDIreg also offers numerous similarity measures for single-channel registration in addition to the manifold alignment scheme used for multi-channel registration. We note that histogram-based mutual information outperforms KNN α-MI in these single-channel registration situations. 16 , which was used in our benchmark studies.

[0191] Supplementary Note 3. HDIprep dimensionality reduction validation Benchmarking dimensionality reduction algorithms (Figures 18A–18J, 19A–19H, and 20A–20H) Our investigation encompasses a wide range of dimensionality reduction methods, ranging from locally nonlinear to globally linear methods, including t-SNE, Uniform Manifold Approximation and Projection (UMAP), potential of heat diffusion for affinity-based transition embedding (PHATE), Isometries (Isomap), Nonnegative Matrix Factorization (NMF), and Principal Component Analysis (PCA).

[0192] To evaluate each method's ability to provide an appropriate data representation while enabling multimodal support, we measured (i.) the ability to generalize to any number of features or necessary degrees of freedom to accurately represent data modalities, (ii.) the ability to succinctly capture data complexity, (iii.) the ability to maximize the information content shared between imaging modalities, (iv.) robustness to noise, and (iv.) computational efficiency.

[0193] (i.-ii.) Estimation of intrinsic data dimensionality. To identify suitable methods for reducing the complexity of mass spectrometry-based image datasets, we hypothesized that introducing more degrees of freedom into the coordinates of the embedded data (i.e., increasing the dimensionality of the embedding) would result in increased similarity between each method's embedding and its higher dimensional counterpart with respect to each algorithm's objective function. Therefore, we examined each algorithm separately with its distinct objective function and identified for each method an appropriate target dimensionality for embedding the data by analyzing the objective function error generated by each method after embedding the data in increased dimensionality. To do this, we used a Euclidean n-space A fit score to estimate the error associated with embedding MSI data into TIFF2024537615000133.tif3128 was generated for each dimensionality reduction method across tissue types and increasing embedding dimensions. For this analysis, we focused on MSI data rather than IMC data, from which we found that application of most dimensionality reduction methods was not feasible due to the data size (number of pixels / high resolution).

[0194] To determine the intrinsic dimensionality of the dataset estimated by each method, we identified points in the error graph of each method where increasing the dimensionality does not reduce the embedding error. To do this, we compare the increase in dimensionality with an exponential increase in the potential locations of a point (i.e., a copy of the real number line). We observed the natural increase in real-valued data dimensionality by modeling the error curves as increasing linearity (increasing linearity). Therefore, we fitted a least-squares exponential regression to the data embedding error curves and constructed 95% confidence intervals (CIs) by modeling a Gaussian residual process. The optimal embedding dimension for each method was selected by simulating samples along the expected value of the fitted curve and identifying the first example of an integer value that falls within the 95% CI of the exponential asymptote. In this way, the minimum degrees of freedom required to capture data complexity was identified. The average error curves for each method across the five random initializations of each algorithm across each MSI dataset are shown in Figures 18A-18J, 19A-19H, and 20A-20H. The methods and rationale used to calculate the embedding error for each method are outlined below.

[0195] UMAP. The UMAP algorithm belongs to the category of manifold learning techniques and aims to optimize the embedding of fuzzy simplicial set representations of high-dimensional data into lower-dimensional Euclidean spaces. In practice, low-dimensional fuzzy simplicial sets are optimized such that the fuzzy set cross-entropy between their high-dimensional counterparts is minimized. The fuzzy set cross-entropy is explicitly defined in terms of definition 1 given by McInnes and Healy

[15] .

[0196] While the theoretical underpinnings of UMAP are based on category theory, practical implementations of UMAP result in weighted graphs. To provide an estimate of the intrinsic dimensionality of the data determined by UMAP, we use an open-source implementation in Python with 15 nearest neighbors, a value of 0.1 for the minimum distance in the resulting embedding, and we allow the algorithm to optimize the embedding with default values ​​of 200 iterations per dimension. We used the Python conversion module of the MATLAB UMAP implementation to calculate the cross-entropy per dimension between the high-dimensional fuzzy simplicial sets and their low-dimensional counterparts.

[0197] t-SNE. t-SNE is a manifold-based dimensionality reduction method that aims to preserve the local structure of a dataset for visualization. To achieve this, t-SNE minimizes the difference between distributions that represent the local similarity between points in the original high-dimensional ambient space and the respective low-dimensional embeddings. The difference between these two distributions is determined by the Kullback-Leibler (KL) divergence between them. As a result, we report the final value of the KL divergence at embedding time as a means to estimate the error associated with t-SNE embedding in each dimension. For all t-SNE calculations, we use an open-source multicore implementation with default parameters (perplexity 30).

[0198] Isomap. Isomap is a manifold-based dimensionality reduction method that uses classical multidimensional scaling (MDS) to preserve inter-point geodesic distances. To do this, the geodesic distances between points are determined by the shortest-path graph distance using a Euclidean metric. The pairwise distance matrix represented by this graph is then embedded into n-dimensional Euclidean space via classical MDS, a metric-preserving technique that finds an optimal transformation for preserving the inter-point Euclidean metric. As a consequence of the implicit linearity in classical MDS, we show that R can be 1-R, the standard linear correlation coefficient between the geodesic distance matrix and the pairwise Euclidean distance matrix in TIFF2024537615000135.tif3128 2 We estimate the intrinsic dimensionality of the data by computing the reconstruction error in each dimension using: For all computations, we select the 15 nearest neighbors for determining the shortest path graph distance, and use the Minkowski metric of the squared norm of the difference. TIFF2024537615000136.tif5128 was selected. All Isomap calculations were performed using Scikit-learn.

[0199] PHATE. PHATE is a manifold-based dimensionality reduction technique developed for data visualization that captures both global and local features of a dataset. PHATE achieves dimensionality reduction by modeling the relationships between data points as t-step random walk diffusion probabilities, and then computing latent distances between data points through a comparison of each pair of points' respective diffusion distributions with all others in the dataset. These latent distances are then embedded into an n-dimensional space using classical MDS followed by metric MDS. Metric MDS is suitable for embedding points with no similarity given by any metric, and reduces the Euclidean constraints imposed by classical MDS to the following stress function S: TIFF2024537615000137.tif17128, where D is the point χ 1… χ n is a metric defined over TIFF2024537615000138.tif4128 are the corresponding embedded data points in dimension n. This stress function spans a least-squares optimization problem. In the scalable form of PHATE used for large data sets, landmarks instead of points are embedded into n-dimensional Euclidean space based on their pairwise latent distances using the above stress function. Out-of-sample embedding of all data points is performed by computing a linear combination of t-step transition matrices from points to landmarks using the embedded landmark coordinates as weights. If the stress function for the metric MDS is zero, then the dimensionality reduction process can fully embed and capture the inter-point distances of the data. This is believed to provide an error estimate used to analyze the intrinsic data dimensionality for the full data set and the full PHATE algorithm; however, in landmark-based computations, not all points are embedded using the metric MDS. Given the linear interpolation scheme and initialization of the scalable PHATE using classical MDS for landmark latent distances, we show that R is 1-R, the linear correlation coefficient between the point-landmark transition matrix and the pairwise Euclidean distance matrix in TIFF2024537615000139.tif3128 2 We assumed that the reconstruction error, given by,provides an estimate of the error associated with embedding the full dataset.,All PHATE calculations were performed in Python using a default,number of 15 nearest neighbors and 2,000 landmark points.

[0200] NMF. Nonnegative matrix factorization (NMF) is a linear dimensionality reduction technique that aims to minimize the divergence between an input matrix X and its reconstruction WH obtained through a matrix factorization. Through this factorization, a linear combination of the columns of W is generated with weights from H. Using the Frobenius norm between X and WH in our calculation, the divergence between the two is TIFF2024537615000140.tif6128. We therefore plotted this divergence or reconstruction error to estimate the error associated with each embedding dimension. For all calculations, we minimax rescaled each channel in the dataset to the range 0-1 so that only positive elements were included in X. All calculations were performed using Scikit-learn.

[0201] PCA. Principal component analysis (PCA) is a linear dimensionality reduction method that aims to capture the principal axes of variation in the data at a global level. To determine the intrinsic dimensionality of a dataset estimated by PCA, the cumulative percentage of residual variance remaining after dimensionality reduction is plotted for each component. Given components 1 ≤ d ≤ n-1, where n is the dimensionality of the original dataset, the percentage of variance explained by the embedding in dimension d is determined by summing the d largest eigenvalues ​​of the covariance matrix of the full dataset. For all calculations, each channel in the dataset was standardized by removing the mean and scaling to unit variance. Standardization is used to ensure that no features dominate the PCA objective function. All calculations were performed using Scikit-Learn.

[0202] (iii.) Evaluation of information content for H&E tissue morphology. In order to have an unbiased evaluation of inter-image information content between the embedded data generated from each dimensionality reduction method and the corresponding H&E stained tissue biopsy sections, three channels were carefully chosen from the MSI data as representative peaks highlighting the morphological characteristics of the tissues (m / z peaks 782.399, 725.373, 566.770 for diabetic foot ulcer, prostate and tonsil), and hyperspectral images were generated, converted to grayscale and registered to the corresponding grayscale converted H&E images (Figures 18A, 19A, and 20A).

[0203] To ensure proper alignment between the manually selected grayscale MSI and grayscale H&E images of diabetic foot ulcers, the mutual information of the registration between the two images and the Dice score of seven pairs of ROIs were evaluated for the initial affine registration and the subsequent nonlinear registration across a hyperparameter grid (Figure 18C). For prostate and tonsil tissues, we optimized only the mutual information (Figures 19C and 20C). The results were then analyzed across a hyperparameter grid to select the optimal parameters for each step in the registration scheme.

[0204] For affine registration, the selected number of resolutions was obtained in a multiresolution pyramidal hierarchy by hyperparameter search. For nonlinear registration, both the number of resolutions and the final uniform grid spacing for the B-spline control points were determined by hyperparameter grid search. In both registrations, the number of resolutions either improved the registration result or the registration remained unchanged. However, during nonlinear registration, finer control point grid spacing schedules led to improved registration as indicated by the mutual information, but they led to regions with unrealistic distortions even with the addition of regularization with deformation bending energy penalty. A value of 300 was chosen for the final grid spacing as a balance between improved registration and increased distortion as indicated by the cost function.

[0205] The resulting deformation fields were then applied to the grayscale hyperspectral images produced from each dimensionality reduction algorithm, uniformly spatially aligning them with the H&E images of each tissue. Before calculating the mutual information between the H&E and embedded MSI images, a non-zero intersection was applied to the pair of images. The non-zero intersection is used to account for any edge effects introduced into the registration by using three manually picked MSI peaks, which would have adversely affected the registration and mutual information calculation if they were not well represented at all locations of the images in our analysis. The mutual information between each registered dimensionality reduced image (n=5 per method) was then calculated using the Parzen window-based method in SimpleITK (Figures 18B, 19B, and 20B).

[0206] (iv.) Evaluation of the algorithms' robustness to noise. Through evaluation of the intrinsic dimensionality of the data, we learned that both high-dimensional imaging modalities (MSI and IMC) follow a manifold structure, where the dimensionality of the data can be approximated with fewer degrees of freedom than the number of parameters initially given in the ambient space. Using this information, in addition to the visual quality of each method's subsequent space mapping back onto tissue, as evidence for the correctness of the assumption of such a manifold structure, we then decided to collate each algorithm's ability to preserve geodesic distances in low-dimensional embeddings with and without "noisy" peaks and / or technical variations.

[0207] To do this, we utilized the denoised manifold preserving (DEMaP) metric. By calculating the DEMaP metric (Spearman's rank correlation coefficient) between the geodesic distances in the ambient space of the peak-picked MSI dataset and the pairwise embedded Euclidean distances between data points from the corresponding non-peak-picked dataset, we evaluated the ability of each algorithm to preserve the manifold structure of the dataset in the presence of noise. All algorithms used were calculated using the Euclidean metric with 15 nearest neighbors, or, since they inherently assume a Euclidean structure, we calculated the geodesic distances in the peak-picked MSI dataset using the 15 nearest neighbors with the Euclidean metric. Peak picking was performed at SCiLS Lab 2018b using orthogonal matching pursuit with a maximum number of peaks of 1,000. The DEMaP scores for each method across five random initializations of each algorithm per MSI dataset are shown in Figures 18I, 19G, and 20G.

[0208] (v.) Evaluation of computational run times. Computational run times for all methods were obtained across diabetic foot ulcer, prostate cancer and tonsil tissue biopsy MSI data, with 1- to 10-dimensional embeddings, over five randomly initialized runs per algorithm (Figures 18J, 19H, and 20H).

[0209] Example 10. Differential diagnosis of diabetic foot ulcer tissue using multimodal imaging and MIAAIM analysis to aid clinical decision making Diabetic foot ulcer (DFU) biopsies collected and stored during the course of clinical care can be analyzed through our method to identify patients at high risk of developing chronic non-healing ulcers and those who heal naturally. Obtained DFU tissue biopsies are imaged using three or more modalities (e.g., H&E, MSI, IMC, IHC, RNAscope, or equivalent imaging methods) to quantify the abundance and spatial distribution of cells, tissue structures, and molecular analytes (e.g., proteins, nucleic acids, lipids, metabolites, carbohydrates, or therapeutic compounds). The resulting images and datasets from all imaging modalities are processed using the MIAAIM analysis pipeline by first processing and extracting pixel-level image data to identify regions, structures, and / or cell populations of interest (e.g., through image segmentation calculations via partitioning based on watershed, Ilastik, UNet, or similar classifications). The resulting processed images and underlying data from each imaging modality are spatially aligned and combined in the MIAAIM method as described above. This involves dimensionality reduction (using UMAP, tSNE, PCA, or similar methods) and clustering of the high-dimensional graph prior to reduction (embedding) of the actual data dimensionality. The resulting combined spatially aligned datasets from three or more image modalities are analyzed to generate a multidimensional signature of the biopsy microenvironment. The signature may include the abundance and distribution of individual cells, tissue structures, or analytes (as defined above), as well as the spatial relationships between two or more such elements (e.g., median distance of an immune cell population from a gradient of metabolites most enriched at tissue boundaries). The resulting multidimensional signatures, which correlate with each of these clinical information, if possible, are then compared and contrasted with existing and newly generated datasets using statistical means to assess the likelihood of wound status and clinical outcomes (e.g., chronic vs. healing) that may aid in clinical decision making.For example, chronic non-healing wounds may be significantly correlated with a signature of a median distance of NK cells from suppressive macrophages less than 20 uM, an elevated abundance of mature B cells compared to adjacent healthy tissue, and the presence of elevated levels of mass spectrometry analytes corresponding to bacterial-associated complement proteins, lipoproteins, and metabolites compared to naturally healing wounds. Based on these outputs identified using the MIAAIM signatures, global or specific associations with clinical outcomes may be presented, after which clinicians may be able to adopt or modify treatment strategies to improve patient care (e.g., by immediately using a more aggressive wound care regimen).

[0210] Example 11. Prognostic evaluation of prostate biopsies for clinical decision support using multimodal imaging and MIAAIM analysis Prostate tissue obtained during diagnostic biopsy or prostatectomy can be analyzed through our method to identify patients at increased risk of aggressive disease or recurrence, and also to guide the evaluation of additional follow-up monitoring and treatment options. Prostate tissue biopsies are imaged using three or more modalities (e.g., H&E, MSI, IMC, IHC, RNAscope, or equivalent imaging methods) to quantify the abundance and spatial distribution of cells, tissue structures, and molecular analytes (e.g., proteins, nucleic acids, lipids, metabolites, carbohydrates, or therapeutic compounds). The resulting images and datasets from all imaging modalities are processed using the MIAAIM analysis pipeline by first processing and extracting pixel-level image data to identify regions, structures, and / or cell populations of interest (e.g., through image segmentation calculations, e.g., via partitioning based on watershed, Ilastik, UNet, or similar classifications). The processed images and the underlying data from each imaging modality are then spatially aligned and combined in the MIAAIM method as described above. This involves dimensionality reduction (UMAP, tSNE, PCA, or similar methods) to cluster high-dimensional graphs before reducing (embedding) the actual data dimensionality. The combined spatially aligned datasets from two or more imaging modalities are analyzed to generate a multidimensional signature of the biopsy microenvironment. The signature may include the abundance and distribution of individual cells, tissue structures, or analytes (as defined above), as well as the spatial relationships between two or more such elements (e.g., the median distance of an immune cell population from a gradient of metabolites most enriched at tissue boundaries). The resulting multidimensional signatures correlating to each of these clinical information are then compared and contrasted with existing and newly generated datasets using statistical means to evaluate known clinical correlates in a novel, integrative manner to identify novel features and / or signatures with prognostic or theranostic utility.In one hypothetical example, rather than the current practice of imaging clinically validated antibody targets for histopathological evaluation either individually or in pairs, the method can collate many targets at once in a highly multiplexed manner (>20 antibodies simultaneously). The resulting data provides a detailed and comprehensive profile of all standard clinical antibodies, including quantification of overall abundance and distribution within the tissue, subcellular distribution, and relative spatial relationships between one or more targets labeled with each individual antibody (e.g., median distance between cell subsets defined using antibody labeling or intensity ratios of spatially matched antibodies). Multimodal imaging data can be collated together with data from matched H&E images and clinical information to generate multimodal signatures that identify risk of progression or recurrence both between and within tumor grade / stage groups. In a second hypothetical example, multimodal imaging of prostate biopsy tissue can be collated to identify signatures associated with responsiveness to treatment. The abundance and distribution of proteins and analytes associated with immune activity and genomic instability can be used to identify spatial relationships that correlate with positive or negative outcomes following treatment with immune-modulating or anti-cancer therapies, and to identify patients most likely to benefit from specific interventions. Based on these outputs identified using the MIAAIM signature, global or specific associations with clinical outcomes can be presented, and clinicians may then improve patient care by assessing the possible utility of additional clinical trials, choosing more frequent follow-up monitoring plans, assessing the risks / benefits of radical prostatectomy, and selecting treatment strategies that reduce the risk of recurrence or metastasis.

[0211] References TIFF2024537615000141.tif200155TIFF2024537615000142.tif240154TIFF2024537615000143.tif22315 7TIFF2024537615000144.tif241157TIFF2024537615000145.tif223157TIFF2024537615000146.tif73154

[0212] Other Aspects Various modifications and variations of the described invention will be apparent to those skilled in the art without departing from the scope and spirit of the invention. Although the invention has been described in connection with specific embodiments, it should be understood that the invention as claimed should not be unduly limited to such specific embodiments. Indeed, various modifications of the described modes for carrying out the invention that are obvious to those skilled in the art are intended to be within the scope of the present invention.

[0213] This application is related to U.S. Patent Application No. 63 / 073,816, the contents of which are incorporated herein in their entirety.

[0214] Other aspects are within the scope of the claims.

Claims

1. 1. A method for identifying cross-modal features from two or more spatially resolved data sets, comprising: (a) registering the two or more spatially resolved data sets to generate an aligned feature image comprising the two or more spatially resolved data sets that are spatially aligned; and (b) extracting the cross-modal features from the aligned feature images. The method comprising:

2. 2. The method of claim 1, wherein step (a) comprises dimensionality reduction for each of the two or more datasets.

3. 3. The method of claim 2, wherein the dimensionality reduction is performed by Uniform Manifold Approximation and Projection (UMAP), Isometry Mapping (Isomap), t-SNE (t-SNE), potential of heat diffusion for affinity-based transition embedding (PHATE), Principal Component Analysis (PCA), Diffusion Map, or Nonnegative Matrix Factorization (NMF).

4. The method of claim 3 , wherein the dimensionality reduction is performed by Uniform Manifold Approximation and Projection (UMAP).

5. The method of any one of claims 1 to 4, wherein step (a) comprises optimizing a global spatial alignment in the aligned feature images.

6. The method of any one of claims 1 to 5, wherein step (a) comprises optimizing a local alignment in the aligned feature images.

7. 7. The method of claim 1, further comprising the step of clustering said two or more spatially resolved datasets to complement said datasets with a similarity matrix representing inter-data point similarities.

8. The method of claim 7 , wherein the clustering step comprises extracting a high-dimensional graph from the aligned feature images.

9. 9. The method of claim 8, wherein the clustering is performed by the Leiden algorithm, the Louvain algorithm, random walk graph partitioning, spectral clustering, or affinity propagation.

10. 10. The method of any one of claims 7 to 9, comprising predicting cluster assignments to unseen data.

11. 11. The method of any one of claims 7 to 10, comprising modelling cluster-cluster spatial interactions.

12. 11. The method of any one of claims 7 to 10, comprising an intensity-based analysis.

13. The method of any one of claims 7 to 10, comprising analysis of cell type abundance or heterogeneity of a given region in said data.

14. 11. The method of any one of claims 7 to 10, comprising analysis of spatial interactions between objects.

15. The method of any one of claims 7 to 10, comprising the analysis of type-specific neighbourhood interactions.

16. 11. The method of any one of claims 7 to 10, comprising the analysis of higher order spatial interactions.

17. The method of any one of claims 7 to 10, comprising the analysis of spatial niche predictions.

18. The method of any one of claims 1 to 17, further comprising the step of classifying said data.

19. 20. The method of claim 18, wherein the classifying step is performed by a hard classifier, a soft classifier, or a fuzzy classifier.

20. The method of any one of claims 1 to 19, further comprising defining one or more spatially resolved objects in the aligned feature images.

21. 21. The method of claim 20, further comprising analyzing the spatially resolved object.

22. 22. The method of claim 21, wherein said step of analyzing the spatially resolved object comprises segmentation.

23. The method of any one of claims 1 to 22, further comprising inputting one or more landmarks into the aligned feature image.

24. The method of any one of claims 1 to 23, wherein step (b) comprises a permutation test for enrichment or depletion of cross-modal features.

25. 25. The method of claim 24, wherein the permutation test generates a list of p-values ​​and / or identities of enriched or depleted factors.

26. 26. The method of claim 24 or 25, wherein the permutation test is performed by a mean permutation test.

27. 27. The method of any one of claims 1 to 26, wherein step (b) comprises multi-domain transformation.

28. 28. The method of claim 27, wherein the multi-domain transformation generates a trained model or a predicted output based on the cross-modal features.

29. 29. The method of claim 27 or 28, wherein the multi-domain transformation is performed by a generative adversarial network or an adversarial autoencoder.

30. 30. The method of any one of claims 1 to 29, wherein at least one of said two or more spatially resolved data sets is an image from immunohistochemistry, imaging mass cytometry, multiplexed ion beam imaging, mass spectrometry imaging, cell staining, RNA-ISH, spatial transcriptome analysis, or co-detection by index imaging.

31. 31. The method of claim 30, wherein at least one of the spatially resolved measurement modalities is immunofluorescence imaging.

32. The method of claim 30 or 31, wherein at least one of the spatially resolved measurement modalities is imaging mass cytometry.

33. The method of any one of claims 30 to 32, wherein at least one of the spatially resolved measurement modalities is multiplexed ion beam imaging.

34. At least one of the spatially resolved measurement modalities is Mass spectrometry imaging, which may be MALDI imaging, DESI imaging, or SIMS imaging The method of any one of claims 30 to 33, wherein

35. At least one of the spatially resolved measurement modalities is Cell staining, which may be H&E, toluidine blue, or fluorescent stain The method of any one of claims 30 to 34, wherein

36. The method of any one of claims 30 to 35, wherein at least one of the spatially resolved measurement modalities is RNA-ISH, which is RNAScope.

37. The method of any one of claims 30 to 36, wherein at least one of the spatially resolved measurement modalities is spatial transcriptome analysis.

38. The method of any one of claims 30 to 37, wherein at least one of the spatially resolved measurement modalities is co-detection with index imaging.

39. Identifying a parameter of interest in the plurality of aligned feature images identified by the method of any one of claims 1 to 38; and comparing the parameters of interest among the plurality of aligned feature images to identify trends; a plurality of aligned feature images for identifying a trend of a parameter of interest, the plurality of aligned feature images comprising:

40. 1. A computer-readable storage medium, comprising: A computer readable storage medium having stored thereon a computer program comprising a routine set of instructions for causing a computer to carry out the steps of the method of any one of claims 1 to 38.

41. 1. A computer-readable storage medium, comprising: A computer program for identifying trends in a parameter of interest within a plurality of aligned feature images identified by the method of any one of claims 1 to 38 is stored on the computer readable storage medium, The computer program comprises a routine set of instructions for causing a computer to carry out the steps of the method of claim 39. The computer-readable storage medium.

Images