A method for identifying cross-modal features from spatially decomposed 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. TIFF2023539830000146.tif47165

Description

[Technical Field]

[0001] Field of Invention This application relates to a method and system for identifying diagnostic methods, prognoses, or theranostics from one or more correlates identified from aligned spatially decomposed datasets. [Background technology]

[0002] background The development of spatially resolved detection modalities has revolutionized diagnostics, prognosis, and theranostics. However, because each modality is generally analyzed independently of the others, the potential for their multimodal applications remains largely unrealized.

[0003] There is a need for new methods that utilize multiple spatially resolved detection modalities to identify multimodal diagnostics, prognoses, and theranostics. [Overview of the Initiative]

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

[0005] In some embodiments, step (a) includes dimensionality reduction for each of two or more datasets. In some embodiments, dimensionality reduction is performed by homogeneous manifold approximation and projection (UMAP), isometric mapping (Isomap), t-distributed stochastic neighbor embedding (t-SNE), PHATE (potential of heat diffusion for affinity-based transition embedding), principal component analysis (PCA), diffusion mapping, or non-negative matrix factorization (NMF). In some embodiments, dimensionality reduction is performed by homogeneous manifold approximation and projection (UMAP). In some embodiments, step (a) includes optimizing global spatial alignment in the aligned feature images. In some embodiments, step (a) includes optimizing local alignment in the aligned feature images.

[0006] In some embodiments, the method further includes a step of clustering a dataset to complement two or more spatially decomposed datasets with a similarity matrix representing the similarity between data points. In some embodiments, the clustering step includes extracting a high-dimensional graph from aligned feature images. In some embodiments, clustering is performed by the Leiden algorithm, the Leuven algorithm, random walk graph decomposition, spectral clustering, or affinity propagation. In some embodiments, the method includes predicting cluster assignments to unseen data. In some embodiments, the method includes a step of modeling cluster-cluster spatial interactions. In some embodiments, the method includes an intensity-based analysis. In some embodiments, the method includes an analysis of the abundance of cell types in the data or heterogeneity in a given region. In some embodiments, the method includes an analysis of spatial interactions between objects. In some embodiments, the method includes an analysis of type-specific neighbor interactions. In some embodiments, the method includes an analysis of higher-order spatial interactions. In some embodiments, the method includes an analysis of predictions of spatial niches.

[0007] In some embodiments, the method further includes a step of classifying data. In some embodiments, the classification step is performed by a hard classifier, a soft classifier, or a fuzzy classifier.

[0008] In some embodiments, the method further includes the step of defining one or more spatially resolved objects in an aligned feature image. In some embodiments, the method further includes the step of analyzing the spatially resolved objects. In some embodiments, the step of analyzing the spatially resolved objects includes segmentation. In some embodiments, the method further includes the step of inputting one or more landmarks into an aligned feature image.

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

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

[0011] In some embodiments, at least one of two or more spatially resolved datasets is an image from co-detection by immunohistochemical analysis, imaging mass cytometry, multiplex ion beam imaging, mass spectrometry imaging, cell staining, RNA-ISH, spatial transcriptome analysis, or 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 multiplex ion beam imaging. In some embodiments, at least one of the spatially resolved measurement modalities is mass spectrometry imaging, such as MALDI imaging, DESI imaging, or SIMS imaging. In some embodiments, at least one of the spatially resolved measurement modalities is cell staining, such as H&E, toluidine blue, or fluorescence staining. In some embodiments, at least one of the spatially resolved measurement modalities is RNA-ISH, such as 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 by index imaging.

[0012] In another aspect, the present invention provides a method for identifying diagnostic methods, prognoses, or theranostics relating to a disease state from two or more imaging modalities, the method comprising the step of comparing a plurality of cross-modal features to identify diagnostic methods, prognoses, or theranostics by identifying a correlation between at least one cross-modal feature parameter and a disease state, the plurality of cross-modal features being identified by the method described herein, each cross-modal feature comprising a cross-modal feature parameter, and the two or more spatially decomposed datasets are outputs from the corresponding imaging modal selected from the group consisting of two or more imaging modalities.

[0013] In some embodiments, cross-modal feature parameters are molecular signatures, single-molecule markers, or marker abundances. In some embodiments, diagnostics, prognoses, or theranostics are individualized to individuals that are the source of two or more spatially decomposed datasets. In some embodiments, diagnostics, prognoses, or theranostics are population-level diagnostics, prognoses, or theranostics.

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

[0015] In yet another aspect, the present invention provides a computer-readable storage medium in which a computer program for identifying cross-modal features from two or more spatially decomposed datasets is stored, the computer program including a routine set of instructions for causing a computer to perform steps of the method described herein.

[0016] In a further aspect, the present invention provides a computer-readable storage medium in which a computer program for identifying a diagnostic method, prognosis, or theranostics relating to a pathological condition from two or more imaging modalities is stored, and the computer program includes a routine set of instructions for causing a computer to perform steps of the method described herein.

[0017] In a further aspect, the present invention provides a computer-readable storage medium on which a computer program for identifying trends of parameters of interest within a plurality of aligned feature images identified by the method described herein is stored, and the computer program includes a routine set of instructions for causing a computer to perform the steps of the method described herein.

[0018] In a further aspect, the present invention provides a method for identifying a vaccine, the method comprising: (a) providing a first dataset of cytometry markers relating to a disease-naive population; (b) providing a second dataset of cytometry markers relating to a population suffering from the disease; (c) identifying one or more markers from the first and second datasets that correlate with clinical or phenotypic scales of the disease; and (d) identifying a composition as a vaccine that (1) induces one or more markers directly correlated with a positive clinical or phenotypic scale of the disease; or (2) identifying a composition as a vaccine that suppresses one or more markers directly correlated with a negative clinical or phenotypic scale of the disease. [Invention 1001] A method for identifying cross-modal features from two or more spatially decomposed datasets, (a) registering the two or more spatially decomposed datasets in order to generate an aligned feature image that includes the two or more spatially aligned spatially decomposed datasets; and (b) Steps to extract the cross-modal features from the aligned feature images. The method, including the method described above. [Invention 1002] The method of the present invention 1001, wherein step (a) includes dimensionality reduction for each of the two or more datasets. [Invention 1003] The method of the present invention 1002, wherein the dimensionality reduction is performed by homogeneous manifold approximation and projection (UMAP), isometric mapping (Isomap), t-distribution type stochastic neighborhood embedding (t-SNE), PHATE (potential of heat diffusion for affinity-based transition embedding), principal component analysis (PCA), diffusion mapping, or non-negative matrix factorization (NMF). [Invention 1004] The method of the present invention 1003, wherein the dimensionality reduction is performed by homogeneous manifold approximation and projection (UMAP). [Invention 1005] A method according to any one of the present invention 1001 to 1004, wherein step (a) includes optimizing global spatial alignment in the aligned feature images. [Invention 1006] A method according to any one of the present invention 1001 to 1005, wherein step (a) includes optimizing local alignment in the aligned feature image. [Invention 1007] Any method 1001 to 1006 of the present invention further comprises the step of clustering the two or more spatially decomposed datasets in order to complement the two or more spatially decomposed datasets with a similarity matrix representing the similarity between data points. [Invention 1008] The method of the present invention 1007, wherein the clustering step includes extracting a high-dimensional graph from the aligned feature images. [Invention 1009] The method of the present invention 1008, wherein clustering is performed by the Leiden algorithm, the Leuven algorithm, random walk graph partitioning, spectral clustering, or affinity propagation. [Invention 1010] Any method of the present invention 1007 to 1009, including prediction of cluster assignment for unseen data. [Invention 1011] A method according to any one of the present invention 1007 to 1010, comprising the step of modeling cluster-cluster spatial interactions. [Invention 1012] Any method of the present invention 1007 to 1010, including analysis based on intensity. [Invention 1013] A method of the present invention, any one of items 1007 to 1010, comprising analysis of the abundance of cell types or heterogeneity in a predetermined region in the aforementioned data. [Invention 1014] A method according to any one of the present invention 1007 to 1010, including the analysis of spatial interactions between objects. [Invention 1015] Any method of the present invention 1007 to 1010, including the analysis of type-specific neighbor interactions. [Invention 1016] Any method of the present invention 1007 to 1010, including the analysis of higher-order spatial interactions. [Invention 1017] Any method of the present invention 1007 to 1010, including analysis of prediction of spatial niches. [Invention 1018] Any method of the present invention 1001 to 1017, further comprising the step of classifying the aforementioned data. [Invention 1020] The method of the present invention 1018, wherein the classification step is performed by a hard classifier, a soft classifier, or a fuzzy classifier. [Invention 1021] Any method of the present invention 1001 to 1020, further comprising the step of defining one or more spatially resolved objects in the aligned feature images. [Invention 1022] The method of the present invention 1032, further comprising the step of analyzing spatially resolved objects. [Invention 1023] The method of the present invention 1033, wherein the step of analyzing spatially resolved objects includes segmentation. [Invention 1024] Any method 1001 to 1023 of the present invention further comprises the step of inputting one or more landmarks into the aligned feature image. [Invention 1025] A method according to any of the 1001-1024 of the present invention, wherein step (b) includes a sorting test relating to the enrichment or depletion of cross-modal features. [Invention 1026] The method of the present invention 1025, wherein the sorting test generates a list of p-values ​​and / or identities for enriched or depleted factors. [Invention 1027] The method of the present invention 1025 or 1026, wherein the sorting test is performed by mean sorting test. [Invention 1028] A method according to any of the 1001 to 1027 of the present invention, wherein step (b) includes multi-domain transformation. [Invention 1029] The method of the present invention 1028, wherein the multi-domain transformation generates a trained model or predictive output based on the cross-modal features. [Invention 1030] The method of the present invention 1028 or 1029, wherein the multi-domain transformation is performed by a generative adversarial network or an autoencoder adversarial. [Invention 1031] Any method of the present invention 1001 to 1030, wherein at least one of the two or more spatially resolved datasets is an image from co-detection by immunohistochemical examination, imaging mass cytometry, multiplex ion beam imaging, mass spectrometry imaging, cell staining, RNA-ISH, spatial transcriptome analysis, or index imaging. [Invention 1032] The method of the present invention 1031, wherein at least one of the spatially resolved measurement modalities is immunofluorescence imaging. [Invention 1033] The method of the present invention 1031 or 1032, wherein at least one of the spatially resolved measurement modalities is imaging mass cytometry. [Invention 1034] A method according to any one of the invention 1031 to 1033, wherein at least one of the spatially resolved measurement modalities is multiplex ion beam imaging. [Invention 1035] At least one of the spatially resolved measurement modalities, Mass spectrometry imaging, which includes MALDI imaging, DESI imaging, or SIMS imaging. The present invention relates to any of the methods described in 1031 to 1034. [Invention 1036] At least one of the spatially resolved measurement modalities, Cell staining, such as H&E, toluidine blue, or fluorescent staining. The present invention relates to any of the methods described in 1031 to 1035. [Invention 1037] A method according to any one of the invention 1031 to 1036, wherein at least one of the spatially resolved measurement modalities is RNA-ISH which is RNAScope. [Invention 1038] A method according to any one of the invention 1031 to 1037, wherein at least one of the spatially resolved measurement modalities is spatial transcriptome analysis. [Invention 1039] A method according to any one of the invention 1031 to 1038, wherein at least one of the spatially resolved measurement modalities is co-detection by index imaging. [Invention 1040] A method for identifying diagnostic methods, prognoses, or theranostics related to a disease state from two or more imaging modalities, The method comprises a step of comparing a plurality of cross-modal features in order to identify a diagnostic method, prognosis, or theranostics by identifying a correlation between at least one cross-modal feature parameter and the pathological condition, wherein the plurality of cross-modal features are identified by any one of methods 1 to 39, each cross-modal feature includes a cross-modal feature parameter, and two or more spatially decomposed datasets are outputs from corresponding imaging modalities selected from the group consisting of the two or more imaging modalities. [Invention 1041] The method of the present invention 1040, wherein the cross-modal feature parameter is a molecular signature, a single-molecule marker, or the abundance of a marker. [Invention 1042] The method of the present invention 1040 or 1041, wherein the diagnostic method, prognosis, or theranostics is individualized to the individual that is the source of the two or more spatially decomposed datasets. [Invention 1043] The method of the present invention 1040 or 1041, wherein the diagnostic method, prognosis, or theranostics is a population-level diagnostic method, prognosis, or theranostics. [Invention 1044] A step of identifying a parameter of interest in a plurality of aligned feature images identified by any of the methods 1001 to 1039 of the present invention, and a step of comparing the parameter of interest among the plurality of aligned feature images to identify a trend. A method for identifying trends in parameters of interest within the plurality of aligned feature images, including the above. [Invention 1045] A computer-readable storage medium, A computer program for identifying cross-modal features from two or more spatially decomposed datasets is stored in the computer-readable storage medium. The computer program includes a routine set of instructions for causing the computer to perform a step of any of the methods 1001 to 1039 of the present invention. The aforementioned computer-readable storage medium. [Invention 1046] A computer-readable storage medium, A computer program for identifying diagnostic methods, prognoses, or theranostics related to a disease state from two or more imaging modalities is stored in the computer-readable storage medium. The computer program includes a routine set of instructions for causing the computer to perform a step of any of the methods 1040 to 1043 of the present invention. The aforementioned computer-readable storage medium. [Invention 1047] A computer-readable storage medium, A computer program for identifying trends in parameters of interest within a plurality of aligned feature images identified by any of the methods 1001 to 1039 of the present invention is stored in the computer-readable storage medium. The computer program includes a routine set of instructions for causing the computer to perform the steps of the method of the present invention 1044, The aforementioned computer-readable storage medium. [Invention 1048] (a) A step of providing a first dataset of cytometry markers for a disease-naive population; (b) Providing a second dataset of cytometry markers for a population suffering from the disease; (c) the step of identifying one or more markers from the first and second datasets that correlate with a clinical scale or phenotypic scale of the disease; and (d)(1) Identifying a vaccine composition capable of inducing one or more markers that directly correlate with a positive clinical or phenotypic scale of the disease; or (2) A step of identifying a composition as a vaccine that can suppress one or more markers that are directly correlated with a negative clinical scale or phenotypic scale of the disease. Methods for identifying vaccines, including those mentioned above. [Brief explanation of the drawing]

[0019] [Figure 1] This diagram illustrates the process of imaging diabetic foot ulcer (DFU) biopsy tissue using multiple modalities, such as H&E staining, mass spectrometry imaging (MSI), and imaging mass cytometry (IMC), followed by the processing and analysis of a multimodal image dataset using an integrated analysis pipeline. [Figure 2]Figure 2A is a high-resolution scan image showing a DFU biopsy tissue section on a microscopic glass slide. Figure 2B is a schematic diagram showing a DFU biopsy tissue section on a glass slide before treatment with a spray matrix solution (optimized for each type of analyte) containing 40% 2,5-dihydroxybenzoic acid (DHB) in a 50:50 v / v acetonitrile:0.1% TFA aqueous solution. Figure 2C is a schematic diagram showing a DFU biopsy tissue section on a glass slide after treatment with a spray matrix solution (optimized for each type of analyte) containing 40% 2,5-dihydroxybenzoic acid (DHB) in a 50:50 v / v acetonitrile:0.1% TFA aqueous solution. Figure 2D is a graph showing the mass-versus-charge-average spectrum of a region of the DFU tissue obtained after characterization using laser desorption, ionization, and mass spectrometry. [Figure 3] This diagram illustrates the process underlying imaging of DFU biopsy tissue or cell lines using IMC. Following sample preparation, staining with metal-labeled antibodies is performed. Laser ablation of the sample produces aerosolized droplets, which are directly transported to the instrument's inductively coupled plasma torch, producing atomized and ionized sample components. Unwanted components are filtered out in a quadrupole ion deflector, where low-mass ions and photons are removed by filtration. High-mass ions, primarily representing metal ions associated with the labeled antibodies, are further pushed to a time-of-flight (TOF) detector, where the time of flight 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. Subsequently, each isotope-labeled sample component is represented by an isotope 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] This flowchart summarizes the multiple steps involved in acquiring a multimodal image dataset and extracting molecular signatures from that dataset. [Figure 5]Figures 5A–5F are a series of graphs showing the estimation of the essential dimension of an MSI dataset using dimensionality reduction methods: t-distribution stochastic neighborhood embedding (t-SNE), homogeneous manifold approximation and projection (UMAP), PHATE (potential of heat diffusion for affinity-based transition embedding), isometric mapping (Isomap), non-negative matrix factorization (NMF), and principal component analysis (PCA). Convergence to the embedding error value indicates that as the dimensionality of the resulting embedding increases, the algorithm's ability to accommodate data complexity no longer improves. Nonlinear dimensionality reduction methods, e.g., t-SNE, UMAP, PHATE, and Isomap, converged to a much lower essential dimension than linear methods, e.g., NMF and PCA, indicating that far fewer dimensions are needed to accurately represent the dataset. [Figure 6] Figures 6A and 6B are graphs showing the computational execution time for each algorithm across 1 to 10 embedding dimensions. The mean and standard deviation (n=5) are plotted for each method across each dimension. The results show that the nonlinear methods t-SNE and Isomap require longer execution times than the nonlinear methods PHATE and UMAP. [Figure 7A] This graph compares the amount of mutual information captured by each dimensionality reduction method tested between the 3D embedding of the grayscale version of the MSI data and the corresponding H&E-stained tissue sections. Mutual information is defined as greater than or equal to zero, and a negative value coincides with the minimization of the cost function in the registration process. The results show that Isomap and UMAP consistently share more information with H&E images than the other tested methods. [Figure 7B]This scheme illustrates the key technical steps of the analysis described herein. Using both complete (noisy) and denoised (peak-picked) datasets, we evaluated the ability of each dimensionality reduction method tested to reconstruct data connectivity (manifold structure). We then calculated a denoised manifold-preserving (DeMaP) metric

[18] between the Euclidean distance of the embedding corresponding to the resulting unpeak-picked data and the geodesic distance of the ambient space (not dimensionally reduced after peak-picking) of the corresponding peak-picked data. [Figure 7C] This graph shows the mean and standard deviation DeMaP metric (Spearman's low correlation coefficient) for all dimensionality reduction methods tested (n=5). This 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 and have a consistent correlation of over 0.85 across 2 to 10 dimensions. [Figure 8] This is a schematic flowchart illustrating the process from mass spectrometry data and image reconstruction to dimensionality reduction using UMAP, and the process of data visualization through a pixelated embedding representation of the mass spectrometry data. [Figure 9] This diagram illustrates the mapping of the 3D embedding of the MSI data after dimensionality reduction by UMAP to the original DFU tissue section, where each of the three UMAP dimensions is colored red (U1), green (U2), or blue (U3). The combined image (RGB image) contains the overlap of all three pseudocolor images. The conversion from the RGB image to grayscale is achieved by adding pixel intensity for each of the three pseudocolor channels, as shown in the equation. For visualization, weighting coefficients can be added to each channel (x1, x2, x3) to adjust the signal contribution for each channel. A representative grayscale image of the dataset is shown in the pseudocolor image. [Figure 10] This is a series of grayscale images of DFU biopsy tissue samples, illustrating a comparison of various linear and nonlinear dimensionality reduction methods. [Figure 11] These are images of DFU biopsy tissue acquired by bright-field microscopy (H&E), MSI, and IMC. The spatial resolution of the three imaging modalities is displayed to illustrate the differences in imaging resolution between bright-field, MSI, and IMC images. [Figure 12] This is a flowchart of a typical grayscale DFU biopsy tissue image illustrating the image registration process between imaging modalities. [Figure 13] This flowchart illustrates the process of aligning multimodal images using the Local Region of Interest (ROI) approach. [Figure 14] This is a flowchart of typical grayscale DFU biopsy tissue images illustrating the process of fine-tuning registration at a local scale. Regions of interest within the toluidine blue images corresponding to each MSI image were selected for local scale registration. [Figure 15] These are 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 DFU biopsy tissue sections. Single-cell coordinates on each ROI were identified by segmentation using IMC parameters, and cell types (cell types 1-12) were defined using clustering analysis of the IMC profiles of the extracted single-cell measurements. Corresponding MSI data were extracted using these single-cell coordinates. Panels A, B, and C show the spatial distribution of MSI parameters identified through sorting tests. Panels A', B', and C' show the spatial distribution of IMC markers of interest before single-cell segmentation. Panels A'', B'', and C'' show panel overlaps A+A', B+B', and C+C'. Panels A''', B''', and C''' show 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 IMC parameters. [Figure 16] This image illustrates an exemplary workflow for integrating image modalities (box labeled (C)) and modeling composite tissue states using MIAAIM. Inputs and outputs (box labeled (A)) are connected to key modules (shaded boxes) or exploration and analysis modules (dashed arrows) via MIAAIM's Nextflow implementation (solid arrows). MIAAIM-specific algorithms (box labeled (D)) are detailed in the corresponding diagrams (bold black text). This includes methods incorporated for applications to single-channel image data types and external software tools that interface with MIAAIM (white boxes). [Figure 17-1] Figures 17A and 17B illustrate HDIprep compression and HDIreg manifold alignment, respectively. The HDIprep compression process may include: (i) a high-dimensional modality, (ii) subsampling, and (iii) a data manifold. The 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 a wide range of dimensions, and exponential regression identifies the steady-state dimensions. The compressed image is reconstructed using the pixel positions. HDIreg manifold alignment may include: (i) Spatial transformations are optimized to align the video to a static image. α-MI is calculated using the KNN graph length between the resampled points (yellow). The edge length distribution panel shows the Shannon MI between the distribution of in-graph edge lengths 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 aligned images. The KNN graph join shows the agreement between modalities. (ii) The optimized transformation aligns the images. The result of the transformation from H&E images (green) to IMC (red) is shown. [Figure 17-2]Figure 17C demonstrates exemplary alignment: (i) Registration from complete tissue MSI to H&E generates T0. (ii) H&E is transformed to IMC complete tissue criterion, generating T1. (iii) ROI coordinates extract the underlying MSI and IMC data in the IMC criterion space. (iv) H&E ROI is corrected in the IMC domain by the transformation, generating T2. Final alignment applies modality-specific transformations. Results for IMC ROI are shown. [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 then 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, represented by the determinant of the spatial Jacobian matrix (center), was saved for downstream use as a control registration. Subsequently, a 3D Euclidean embedding of the MSI data was created using random initialization of each dimensionality reduction algorithm (bottom). Then, using these embeddings, an RGB image was created according to the procedure described above. Finally, the spatial transformation created by registering the manually identified peaks together with the H&E image was applied to the dimensionality-reduced grayscale image, and each was aligned to the grayscale H&E image. Figure 18B: The mutual information between each aligned grayscale embedding image (n=5 per method) and the grayscale H&E image was calculated using Partzen window histogram density estimation with a histogram bin width of 64. The plot was oriented so that the results matched the concept of a "cost function" in the optimization situation, where the objective is to minimize cost. Therefore, larger negative values ​​represent higher mutual information. UMAP consistently captures multimodal informational content regarding the H&E data. Figure 18C: Optimization of image registration using mutual information as a cost function between manually identified mass spectrometry peaks in the grayscale version and the grayscale H&E image (Figure 18A, top), and external validation using dice scores for seven manually annotated regions. The registration parameters used in the final registration used in Figure 18A are shown by dashed lines. Registration was first performed by aligning the images with multi-resolution affine registration (left).Subsequently, manually identified mass spectrometry peaks from the converted grayscale version were registered to the grayscale H&E image using nonlinear multiresolution registration. Figure 18D: Average neighborhood entropy of each pixel (n=5) calculated within a 10-pixel disk across dimensionality reduction algorithms. The results show that UMAP has the ability to enhance structure in tissue sections. Figure 18E: Manual annotation of the grayscale H&E image used to verify registration quality by the controlled deformation field used for mutual information calculation in Figure 18B. Figure 18F: Cropped region using the same spatial coordinates as Figure 18E for the manually annotated region used to calculate the dice score in Figure 18C. The results show good spatial overlap across different annotations. Figure 18G: Radar plot showing a comparison of the performance of dimensionality reduction algorithms across a wide range of data representations—linear, nonlinear, local, and global data structure preservation (t-SNE, UMAP, PHATE, Isomap, NMF, PCA). The graph shows the average algorithm execution time (n=5) (top, logarithmic transformation), estimated stationary-state manifold embedding dimension (right), noise tolerance (bottom), and multimodal mutual information with DFU MSI data (left). All plots are oriented so that larger values ​​represent better algorithm performance. The results demonstrate 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 histological images. Figure 18H: Intrinsic dimensions of MSI data estimated by each dimensionality reduction method. Embedding error (y-axis) is not comparable between plots. Mean and standard deviation (n=5) embedding error are plotted over embedding dimensions from 1 to 10. Convergence to the y-axis indicates that the algorithm's ability to capture data complexity no longer improves as the dimensionality of the resulting embedding increases.The results show that the intrinsic dimensions estimated by nonlinear methods (t-SNE, UMAP, PHATE, Isomap) are far fewer than those of linear methods (NMF, PCA), meaning that fewer dimensions are needed to accurately represent the dataset. Figure 18I: Denoised Manifold Preservation (DeMaP) metric between the resulting Euclidean distance of the embedding corresponding to the unpeak-picked data and the geodesic distance of the corresponding ambient space (undimensionally unreduced after peak-picking) of the peak-picked data. Results showing the mean and standard deviation DeMaP metrics (Spearman's low 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 and have consistent correlations of over 0.85 across 2 to 10 dimensions. Figure 18J: Computation execution time per algorithm across 1 to 10 embedding dimensions. Mean and standard deviation (n=5) plotted for each method across each number of dimensions. Nonlinear methods such as t-SNE and Isomap require longer execution times than nonlinear methods such as PHATE and UMAP. Linear methods require minimal execution time; however, these methods cannot concisely incorporate data complexity. [Figure 18-2] See the explanation in Figure 18-1. [Figure 18-3] See the explanation in Figure 18-1. [Figure 18-4] See the explanation in Figure 18-1. [Figure 19-1]Figures 19A–19H provide an overview of the performance of a dimensionality reduction algorithm for summarizing mass spectrometry imaging data of prostate cancer. Figure 19A: Same as Figure 18A except for the prostate cancer tissue biopsy. Figure 18B: Same as Figure 18B except for the prostate cancer tissue biopsy. Figure 19C: Optimization of image registration using mutual information as a cost function between manually identified mass spectrometry peaks in the grayscale version and the grayscale H&E image (Figure 19A, top). The registration parameters used in the final registration used in Figure 19A are shown by the dashed line. Registration was first performed by aligning the images with multi-resolution affine registration (left). Then, the manually identified mass spectrometry peaks in the transformed grayscale version were registered to the grayscale H&E image using nonlinear multi-resolution registration. Figure 19D: Same as Figure 18D except for the prostate cancer tissue biopsy. Figure 19E: Same as Figure 18G except for the prostate cancer tissue biopsy. Figure 19F: Same as Figure 18H except for the prostate cancer tissue biopsy. Figure 19G: Same as Figure 18I except for the prostate cancer tissue biopsy. The nonlinear methods Isomap, PHATE, and UMAP all consistently preserve the manifold structure without pre-filtering of the data and have consistent correlations of over 0.75 across 2 to 10 dimensions. Figure 19H: Results showing computational execution time for each algorithm across 1 to 10 embedding dimensions. Mean and standard deviation (n=5) are plotted for each method across each dimension. The results show that the nonlinear methods t-SNE, PHATE, and Isomap require longer execution times than UMAP. Linear methods require minimal execution time; however, these methods cannot succinctly capture data complexity and are not noise-tolerant. [Figure 19-2] See the explanation in Figure 19-1. [Figure 19-3] See the explanation in Figure 19-1. [Figure 19-4] See the explanation in Figure 19-1. [Figure 20-1]Figures 20A–20H provide an overview of the performance of dimensionality reduction algorithms for summarizing tonsil mass spectrometry imaging data. Figure 20A: Same as Figure 18A except for tonsil biopsy. Figure 20B: Same as Figure 18B except for tonsil biopsy. Isomap and NMF consistently capture multimodal informational content regarding H&E data. Figure 20C: Same as Figure 19C except for tonsil biopsy. Figure 20D: Same as Figure 18D except for tonsil biopsy. Figure 20E: Same as Figure 18G except for tonsil biopsy. Figure 30F: Same as Figure 18H except for tonsil biopsy. Figure 20G: Same as Figure 18I except for tonsil biopsy. Figure 20H: Same as Figure 18J except for tonsil biopsy. [Figure 20-2] See the explanation in Figure 20-1. [Figure 20-3] See the explanation in Figure 20-1. [Figure 20-4] See the explanation in Figure 20-1. [Figure 21A] We demonstrate that spectral central landmarks reproduce the steady-state manifold embedding dimension across tissue types and imaging techniques. Figure 21A: The sum of squared errors of exponential regression fits the steady-state embedding dimension selection from spectral landmarks compared to the complete mass spectrometry imaging dataset across tissue types. The difference between the exponential regressions fitted to the cross-entropy of landmark central embeddings and complete dataset embeddings approaches zero as the number of landmarks increases. The dashed line shows MIAAIM's default selection of 3,000 landmarks for calculating the steady-state manifold embedding dimension. [Figure 21B] We demonstrate that spectral central landmarks reproduce the stationary manifold embedding dimension across tissue types and imaging techniques. Figure 21B: Same as Figure 21A except that pixels are subsampled in the imaging mass cytometry region of interest. [Figure 22A]Figures 22A and 22B demonstrate that UMAP embedding of spatially subsampled imaging mass cytometry data by extrasample projection reproduces full data embedding (Figure 22B) while reducing execution time (Figure 22A) in a diabetic foot ulcer sample. [Figure 22B] See the explanation in Figure 22A. [Figure 23A] Figures 23A and 23B demonstrate that UMAP embedding of spatially subsampled imaging mass cytometry data by extrasample projection reproduces full data embedding (Figure 23B) in a prostate cancer sample while reducing execution time (Figure 23A). [Figure 23B] See the explanation in Figure 23A. [Figure 24A] Figures 24A and 24B demonstrate that UMAP embedding of spatially subsampled imaging mass cytometry data by extrasample projection reproduces full data embedding (Figure 24B) in tonsil samples while reducing execution time (Figure 24A). [Figure 24B] See the explanation in Figure 24A. [Figure 25] Figures 25A and 25B demonstrate the scaling of MIAAIM image compression to wide-field and high-resolution multiplexed image datasets by incorporating parametric UMAP. Figure 25A: Multiplexed CyCIF images of lung adenocarcinoma metastases to lymph nodes (n = approximately 100 million pixels, 0.65 μm / pixel resolution, 44 channels, 27 antibodies) with corresponding steady-state UMAP embedding and spatial reconstruction (showing 3 UMAP channels with 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 = approximately 256 million pixels, 0.65 μm / pixel resolution). [Figure 26-1]Figures 26A–26I show that microenvironment correlation network analysis (MCNA) correlates protein expression with 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 ionic intensities for the top 5 proteins showing positive and negative correlations. The color scheme indicates module assignment. The heatmap (right) shows Spearman's Rho. Figure 26C: Exponentially weighted moving average of normalized mean ionic intensities per module, aligned away from the center of the DFU wound. Figure 26D: Raw IMC nuclei (Ir) and CD3 staining in the ROI (left) (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 and top molecular correlates (right). Figure 26E: Same as Figure 26D except for a different ROI. Figure 26F: Unsupervised phenotyping. Shaded boxes represent the CD3+ population. Heatmap shows normalized protein expression. Figure 26G: MCNA UMAP colored to reflect the correlation between ions and Ki-67 in the CD3+ and CD3- populations. The color scheme represents Spearman's Rhodes, and the size of the points represents the negative logarithmically transformed Benjamin-Hockberg corrected P-value of the correlation. Figure 26H: Tornado plot comparing the top 5 ions showing differential negative and positive correlations between Ki-67 and CD3+ with the CD3- cell population. The X-axis shows CD3+-specific Ki-67 values. The color scheme of each bar shows the change in correlation from the CD3- population to the CD3+ population. Figure 26I: Box plot showing the ionic intensities of top ions (top, positive; bottom; negative) differentially correlated with CD3+-specific Ki-67 expression across ROIs on DFU. The tissue map shows differentially related top CD3+ Ki-67 correlates (upper, positive; lower; negative), with ROIs on tissues containing CD3+ cells indicated by boxes (white). [Figure 26-2] See the explanation in Figure 26-1. [Figure 26-3] See the explanation in Figure 26-1. [Figure 26-4] See the explanation in Figure 26-1. [Figure 27-1] Figures 27A–27H illustrate cobordist projection and domain transfer using (i-)PatchMAP. Figure 27A: Diagram showing PatchMAP stitching between bounded manifolds (reference and query data) to form a cobordist (gray), information transfer across cobordist geodesics (top), and cobordist projection visualization (bottom). Figure 27B: Bounded manifold stitching simulation. PatchMAP projection (manually drawn dashed lines indicate stitching) and UMAP projection of the integrated data are shown with NN values ​​that maximize SC for each method. Figure 27C: Data transfer from MSI to IMC using i-PatchMAP. The line plot shows Spearman's Rhoad 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), tonsils (green), and prostate (orange) pixels based on MSI profiles. Individual plots show IMC representation for DFU single cells (right). Figure 27G: Data transfer from MSI to IMC from DFU single cells to complete tissue. Figure 27H: Data transfer from MSI to IMC from DFU single cells to tonsil tissue. [Figure 27-2] See the explanation in Figure 27-1. [Figure 27-3] See the explanation in Figure 27-1. [Figure 27-4] See the explanation in Figure 27-1. [Figure 28]Figures 28A and 28B demonstrate that PatchMAP preserves the structure of bounded manifolds while accurately embedding the relationships between bounded manifolds into cobordism. Figure 28A: PatchMAP embedding of the MNIST digit dataset (n=70,000) randomly partitioned into two bounded manifolds of the same size. When the nearest neighbor values ​​are small, PatchMAP is similar to the UMAP embedding because open neighbors are preserved after intersecting pairwise nearest neighbor queries. Under these conditions, the intersection operation is similar to the fuzzy union performed by UMAP. When the nearest neighbor values ​​are large, PatchMAP incorporates the relationships between manifolds into cobordism while preserving the structure of the bounded manifold. Here, PatchMAP aligns the bounded manifolds along the principal axes to produce something close to a mirror image. This results in the data being equally partitioned in half, which is then incorporated using cobordist geodesic distances. Figure 28B: Verification of Figure 27B against the complete MNIST digit dataset. Each digit in the dataset is considered a bounded manifold. When the nearest neighbor values ​​are smaller, PatchMAP is similar to UMAP embeddings, and when the nearest neighbor values ​​are larger, PatchMAP can accurately model cobordist geodesic distances. [Modes for carrying out the invention]

[0020] Detailed explanation Overall, the present invention provides a method and a computer-readable storage medium for processing two or more spatially decomposed datasets for identifying cross-modal features, identifying diagnostic methods, prognoses, or theranostics related to disease conditions, or identifying trends in parameters of interest.

[0021] As used herein, the term "theranostics" refers to diagnostic treatments. For example, theranostic approaches may be used in personalized medicine.

[0022] This method is designed as a general framework for identifying cross-modal features that can be used as high-value or actionable indicators (e.g., biomarkers or prognoses) consisting of one or more parameters uniquely revealed through the creation and analysis of multidimensional maps, by correlating spatially resolved datasets from a wide range of diverse origins (e.g., laboratory samples, various imaging modalities, geographic information system data) with other aligned data.

[0023] The present invention may be a method for identifying cross-modal features from two or more spatially decomposed datasets, comprising the steps of (a) registering two or more spatially decomposed datasets to generate an aligned feature image containing two or more spatially aligned feature images; and (b) extracting cross-modal features from the aligned feature image.

[0024] The present invention may be a method for identifying diagnostic methods, prognoses, or theranostics related to a disease state from two or more imaging modalities. The method includes the step of comparing multiple cross-modal features to identify diagnostic methods, prognoses, or theranostics by identifying a correlation between at least one cross-modal feature parameter and a disease state. The multiple cross-modal features may be identified as described herein. In the method described herein, each cross-modal feature includes a cross-modal feature parameter. Two or more spatially resolved datasets are outputs from corresponding imaging modalities selected from the group consisting of two or more imaging modalities described herein.

[0025] The present invention may be a method for identifying trends in a parameter of interest within a plurality of aligned feature images identified by the method described herein. The method includes the steps of identifying a parameter of interest in a plurality of aligned feature images, and comparing the parameter of interest among the plurality of aligned feature images to identify trends.

[0026] Figure 4 summarizes the necessary 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 the expression / abundance of parameters of interest for modeling and predicting biological processes at multiple scales: cellular niches (fine local contexts), local tissue heterogeneity (local population contexts), tissue-wide heterogeneity and trending features (global contexts), and disease / tissue states (combinations of local and global tissue contexts).

[0027] For data derived from biological samples relevant to biomedical and research applications, this 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], co-detection by index imaging (CODEX)[7], and imaging mass cytometry (IMC)[8].

[0028] The present invention also provides a computer-readable storage medium. A computer program for identifying cross-modal features from two or more spatially resolved datasets may be stored on the computer-readable storage medium as described herein, and the computer program includes a routine set of instructions for causing a computer to perform steps of a method for identifying cross-modal features from two or more spatially resolved datasets. A computer program for identifying diagnostic methods, prognoses, or theranostics relating to a pathological condition from two or more imaging modalities may be stored on the computer-readable storage medium, and the computer program includes a routine set of instructions for causing a computer to perform steps of the corresponding method described herein. A computer program for identifying trends of parameters of interest within a plurality of aligned feature images identified by the corresponding method described herein may be stored on the computer-readable storage medium, and the computer program includes a routine set of instructions for causing a computer to perform steps of the corresponding method described herein.

[0029] All computer-readable storage media described herein exclude any temporary media (e.g., volatile memory, carriers, such as data signals integrated into carriers in networks, such as the Internet). Examples of computer-readable storage media include non-volatile memory media, such as magnetic storage devices (e.g., conventional "hard drives," RAID arrays, floppy disks), optical storage devices (e.g., compact discs (CDs) or digital video discs (DVDs)), or integrated circuit devices, such as solid-state drives (SSDs) or USB flash drives.

[0030] Registration of spatially decomposed datasets Integrating spatially resolved datasets (e.g., spatially resolved datasets of high parameters from various imaging modalities) presents challenges due to the existence of uncertain statistical relationships between different modalities, considering the differing spatial resolutions between modalities, the possibility of spatial deformation and misalignment, technical variations within modalities, and the goal of discovering new relationships. Therefore, 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.

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

[0032] Image registration is often perceived as a fitting problem in which a quality function is iteratively optimized by applying transformations to images to spatially align one or more images. In practice, image registration frameworks typically consist of sequential pairwise registration or groupwise registration to selected reference images; the latter has been proposed as a way to register multiple images in a single optimization step and eliminate the biases incurred by selecting reference images and, consequently, reference modalities [9,10]. More recently, both of these frameworks have been extended to learning-based registration that can handle large datasets using spatial transformer networks [11,12,13,14]. In our investigation of appropriate registration pipelines, we have recognized the usability of groupwise registration schemes and learning-based models, respectively, in situations where tissue morphology changes significantly between adjacent sections (as in glandular prostatic tissue) or when dealing with large amounts of data.

[0033] The method disclosed herein is central to a sequential pairwise registration scheme that can be guided and optimized at each step. Thus, the method disclosed herein provides a platform for registration of multiple samples in datasets across acquisition techniques and tissue types, as well as for single-image registration.

[0034] Image registration Dimensionality Reduction High-parameter datasets, often congested by technological fluctuations and noise, present challenges in their analysis and integration. Spatial integration of each modality currently requires the presentation of representative images that enable statistical correspondence with other modalities within the image registration scheme. In the datasets under consideration, manual identification of such images becomes immediately difficult due to the number of parameters acquired and the complex relationships between these parameters.

[0035] The method of the present invention includes the step of registering two or more spatially decomposed datasets in order to generate a feature image containing two or more spatially aligned spatially decomposed datasets. Automatic definition of image features can be achieved using techniques that embed data in a space having a metric adapted to the construction of an entropy-wide graph. Such techniques include dimensionality reduction and compression techniques that embed high-dimensional data points (e.g., pixels) in Euclidean space. Non-limiting examples of dimensionality reduction techniques include homogeneous manifold approximation and projection (UMAP)

[15] , isometric mapping (Isomap)

[16] , t-distribution stochastic neighborhood embedding (t-SNE)

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

[18] , principal component analysis (PCA)

[19] , diffusion mapping

[20] , and non-negative matrix factorization (NMF)

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

[0036] Homogeneous Manifold Approximation and Projection (UMAP) is a machine learning technique for dimensionality reduction. UMAP is constructed from a theoretical framework based on Riemannian geometry and algebraic topology. It is, as a result, a practical and extensible algorithm applicable 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 runtime performance. Furthermore, UMAP has no computational constraints on embedding dimensions, making it a promising dimensionality reduction technique for a wide range of machine learning applications.

[0037] An isometric mapping (Isomap) is a nonlinear dimensionality reduction method. It is used to compute a quasi-isometric low-dimensional embedding of a set of high-dimensional data points. The method allows for the estimation of the internal geometry of a data manifold based on an approximation of the neighbors of each data point on the manifold.

[0038] t-distribution stochastic neighbor embedding (t-SNE) is a machine learning algorithm for nonlinear dimensionality reduction that enables the representation of high-dimensional data in a 2D or 3D lower-dimensional space for better visualization. Specifically, it models each high-dimensional object with 2D or 3D points such that, with high probability, similar objects are modeled by neighboring points and dissimilar objects are modeled by distant points.

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

[0040] Principal component analysis (PCA) is a technique for reducing the dimensionality of large datasets by creating new uncorrelated variables that sequentially maximize the variance.

[0041] Diffusion mapping is a dimensionality reduction or feature extraction method that computes a family of embeddings of a dataset into Euclidean space (often low-dimensional), whose coordinates can be calculated from the eigenvectors and eigenvalues ​​of the diffusion operator on the data. The Euclidean distance between points in the embedded space is equivalent to the diffusion distance between 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.

[0042] Non-negative matrix factorization (NMF) is a dimensionality reduction method that decomposes a non-negative matrix into the product of two non-negative matrices.

[0043] This dimensionality reduction process is often data-dependent, and observing the performance of the chosen algorithm is necessary for a proper representation of the dataset. For an exemplary dataset, the dimensionality reduction method chosen by the inventors is the homogeneous manifold approximation and projection (UMAP) algorithm

[17] . The inventors' results (Figures 5, 6, 7A, 7B, and 7C) show that this manifold-based nonlinear technique provides the best MSI data representation among the methods considered in a multimodal comparison with H&E, based on standard image registration and computational complexity testing, robustness to noise, and the ability to capture information with low-dimensional embeddings. The dimensionality reduction processes listed above are applicable to all datasets under consideration, but manual curation of representative modality features is possible and they are considered “induced” dimensionality reduction.

[0044] 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 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 there is no guarantee that pixels / data points are distinguishable from the image background (typically given a zero value). To avoid this, each channel is linearly rescaled to the range of zero to 1 according to the process in

[23] , making it possible to distinguish between the foreground (spatial location containing acquired data) and the background (spatial location with no information).

[0045] Enter landmarks The image registration process may include, for example, user-initiated input of landmarks. User-initiated input of landmarks is not a required step to complete image registration. Instead, this step may be included to improve the quality of the results, for example, when 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. User-defined landmarks may be input before the optimization of registration parameters.

[0046] In a particular preferred embodiment, user input is incorporated after dimensionality reduction. Alternatively, user input may be incorporated before dimensionality reduction by using the spatial coordinates of the raw data. In practice, user-defined landmarks may be stored within image visualization software (e.g., Image J, available from imagej.nih.gov).

[0047] Optimization of registration parameters By selecting features for modality registration by dimensionality reduction, the parameters of the alignment process can be optimized semi-automatically by hyperparameter grid search and, for example, by manual verification. The calculations for the registration procedure in the current implementation (separated from the dimensionality reduction process) may be performed, for example, by open-source Elastix software

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

[0048] Optimization of global spatial alignment The image registration process may include optimizing the global spatial alignment of registration parameters. This optimization of global spatial alignment may be performed for two or more datasets after reducing their dimensions.

[0049] Using hyperparameter grid search, registration parameters can be optimized to ensure proper alignment of each modality at the whole tissue scale for coarse-grained analysis (e.g., calculating the gradient of markers of interest across the tissue, marker / cell heterogeneity across the tissue, and identifying regions of interest (ROIs) for further examination). In some embodiments, spatial alignment of datasets can be performed in a propagation manner by registering complete tissue sections (e.g., MSI, H&E, and toluidine blue stained images) for each dataset. Subsequently, spatial coordinates for ROIs (e.g., IMC ROIs taken from toluidine blue stained images) can be used to correct for any local deformations that require further adjustment for finer-grained analysis (Figures 14 and 15).

[0050] In the exemplary datasets described herein, the spatial resolution for each modality was as follows: MSI approximately 50 μm, H&E approximately 0.2 μm, and IMC approximately 1 μm.

[0051] The methods described herein can preserve the spatial coordinates of high-dimensional, high-resolution structures and tissue morphologies. Therefore, in some of the methods described herein, higher-resolution ROIs may remain unchanged at each step of the registration scheme (e.g., the exemplary registration scheme described herein). Such higher-resolution ROIs may serve, for example, as a final reference image against which all other images are aligned. MSI data have been shown to reflect tissue morphologies present in traditional tissue staining

[24] . Considering this correspondence, combined with the ability of tissue (H&E) staining to capture the spatial configuration of cells, we choose to examine H&E images as a medium between MSI and IMC datasets, and as a key element for spatially aligning all modalities. Due to computational resource limitations, a resolution of approximately 1.2 μm per pixel is used for H&E images in the registration process.

[0052] However, the use of a hierarchical multi-resolution registration scheme similar to the one implemented in our dataset also has the potential to register datasets of arbitrary resolutions.

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

[0054] In the exemplary datasets described herein, the lack of morphological information currently available at the whole tissue scale for acquired IMC images, a consequence of the destructive nature of IMC techniques, necessitates this additional step to compensate for local deformations occurring within each ROI. For this purpose, single-cell multiplex imaging techniques that enable the acquisition of complete tissue data, such as tissue-based cyclic immunofluorescence (t-CyCIF)[4] and simultaneous detection by indexing (CODEX)[7], provide both a large-scale coarse analysis of specimen heterogeneity and a local analysis of ROIs; however, the dilution of single-cell relationships resulting from the overall tissue heterogeneity, combined with the potential encounter of artifacts on the edges of complete tissue specimens, often necessitates a finer analysis of regions of interest (ROIs) within the complete tissue. Consequently, with complete tissue specimens, slides are often scanned for coarse morphological characteristics using low-magnification fields before obtaining a finer, cell-level analysis at higher magnifications.

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

[25] . The control point grid spacing in a nonlinear deformation model represents the spacing between noses that fix the deformation plane of the transformed image. When used with a multiplex resolution registration approach, uniform control point spacing for nonlinear deformations is often scaled with image resolution. Thus, coarse nonlinear deformations are corrected before finer, higher-resolution registration at the local scale. Our pyramidal approach to complete tissue registration attempts to mitigate misalignment due to very fine or coarse grid spacing. However, we ultimately choose to ensure fine-grained structural registration for each ROI by reducing the sampling space and shifting the cost function from a global total tissue cost to one centered on each ROI.

[0056] Final registration proceeds by constructing the resulting transformations in a propagation scheme after dimensionality reduction, global spatial alignment optimization, and local alignment optimization. Subsequently, the original data corresponding to each modality is spatially aligned with all others by applying their respective transformation sequences to each channel.

[0057] 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 predefined spatially resolved object analysis). 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 (a similar problem arises when selecting feature images for registration). Clustering is a method of grouping similar data points (e.g., pixels, cells) together for the purpose of reducing data complexity and preserving the overall data structure. Through this approach, individual pixels in an image can be grouped together to aggregate homogeneous areas of tissue, providing a more interpretable discretized version of the complete image, reducing the complexity of analysis from millions of individual pixels to a predetermined number of clusters (e.g., tens to hundreds). When used in conjunction with heatmaps or other forms of data visualization, an overview of each cluster or tissue region can be visualized in a single image, helping to quickly interpret the profile of each region.

[0058] In the exemplary datasets described herein (Figures 7B and 7C), the UMAP algorithm, which has proven robust to noisy (variable) features, and its computational efficiency, enabled iterative data partitioning over reasonable periods. As a result of UMAP's noise robustness and ability to incorporate complexity, we found that the algorithm is best suited to construct mathematical representations of very high-dimensional data, such as data obtained from MSI or similar methods where hundreds to thousands of channels are available per image.

[0059] The dimensionality reduction unit of the UMAP algorithm operates by maximizing the information content contained in the lower-dimensional graph representation of the dataset compared to its higher-dimensional counterpart

[15] . In a particular preferred embodiment, the dimensionality reduction optimization scheme can encompass the higher-dimensional graph itself. As a result, we extract the higher-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

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

[15] , Isomap

[16] , PHATE

[18] , etc.). The methods described herein perform clustering of the higher-dimensional graph before actual data dimensionality reduction (embedding) so that clusters are formed based on configurations representing a global manifold structure. The exemplary clustering approaches used herein preserve the global features of the data, in contrast to embeddings generated by local dimensionality reduction using methods such as t-SNE or UMAP (preferably t-SNE)

[18]

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

[31] , the approach employed in our exemplary dataset reduces the burden of identifying principal components from the raw data before clustering, although we have 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).

[0060] Regardless of the choice of clustering algorithm, the simplified representation of data through the process described above subsequently enables the execution of numerous analyses, ranging from predicting cluster assignments to unseen data, directly modeling cluster-cluster spatial interactions, to performing traditional intensity-based analyses independent of spatial context. The choice of analysis depends on whether the focus is on features outside the scope of the current study and / or task-spatial context (e.g., cell type abundance, heterogeneity of a given region in the data), or on spatial interactions between objects (e.g., type-specific neighbor interactions

[26] , higher-order spatial interactions-first-order interactions[7], prediction of spatial niches

[27] ). The resulting analyses and predictions can then be used as features to prove disease diagnosis and prediction, or as indicators of biological processes of interest for purely scientific reasons.

[0061] Clustering allows for the matching of data in an unsupervised manner. However, it is equally easy to manually annotate pixels on an image to identify a set of features corresponding to annotations of interest. In our UMAP embedding representation of an exemplary dataset from diabetic foot ulcer biopsy tissue, for example, two diametrically opposed extremes of tissue health can be easily identified. These tissue states can be labeled and then combined to provide the same analysis listed above. In either case, annotations and cluster identities function as a set of discretized labels that can be further analyzed.

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

[0063] Using soft classifiers (e.g., conditional probabilities generated by random forests, neural networks with sigmoid final activation functions, etc.) can lead to the creation of further probability maps, for example, regarding disease / healthy tissue regions and diagnoses. The concept of these probability maps is best illustrated by the pixel classification workflow in the image analysis software Ilastik

[38] . After classification using a random forest classifier, relevant features can then be extracted and used to make predictions easier to understand. For example, we identified prominent features between tissue regions using the MSI parameter that had the greatest impact on the cluster conditional probabilities in our random forest classification.

[0064] On the other hand, hard classifiers allow for clear assignment of classes to data and are therefore useful in imposing constraints when clear categorical assignments (decisions) are required. In our exemplary dataset, we clustered the MSI dataset at the pixel level using the UMAP-based method described above and extended the cluster assignments to new pixels by assigning pixels to the highest probability cluster using a random forest classifier (hard classification). This instruction was made due to computational constraints and computational efficiency, in addition to the ability to identify nonlinear decision boundaries generated by our manifold clustering scheme, which is robust to parameter selection

[37] .

[0065] Analysis of predefined spatially decomposed objects (cells, tissue structures, etc.) Tissue samples contain many objects of interest, which are cells or other morphological features (e.g., blood vessels, nerves, extracellular matrix; or whole structures such as hair follicles or tumors). Identifying the spatial coordinates of these objects is crucial for understanding the imaging dataset at a level higher than the pixel level. In the exemplary dataset described herein, the IMC modality contains single-cell resolution data, and the goal of the analysis is to connect this single-cell information to the parameters of other modalities. Single-cell multiplex imaging analysis can apply computer vision and / or machine learning techniques to pinpoint the coordinates of cells in the image, use those coordinates to extract aggregated pixel-level data, and then 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, including 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 for data aggregation for the application of any of the above analyses (e.g., clustering, spatial analysis, etc.). Importantly, for our application, this segmentation enables the aggregation of pixel-level data for each single cell and enables cell clustering regardless of spatial location. This process allows for the formation of cell identity based on traditional surface or activation marker staining using the IMC modality alone. A similar approach can be applied to any object under conditions where the analysis and aggregation of pixel-level data are guaranteed.

[0066] Feature extraction and analysis of multimodal data The methods described herein may include, for example, a step of comparing data from different modalities concerning a spatially resolved object by using their spatial coordinates. The image registration process spatially aligns all imaging modalities so that, as a result, an object can be defined in any one of the modalities used, while relevant features are still accurately maintained across all modalities. In the example described herein (Figure 15), the IMC dataset was used to identify single-cell coordinates, which were then used to extract single-cell features from both the aligned MSI pixel-level data and the IMC pixel-level data itself. Subsequently, the data were clustered based on single-cell measurements for the IMC modality alone and for the MSI modality alone. Clustering of IMC single-cell measurements can be used to determine cell types. This ability to integrate multiple imaging modalities made it possible to perform sorting tests as a function of the corresponding cell types defined in the IMC dataset regarding the enrichment or depletion of certain features in the MSI modality. Alternatively, the methods described herein can also identify which IMC features are depleted or enriched based on the cell type defined by the MSI modality. This type of cross-modal analysis can be extended to any number of parameters and any number of modalities. The sorting test evaluates the randomized mean of each parameter for its observations, regardless of modality, thereby enabling a one-to-all comparison, in which case the measured values ​​being evaluated are aggregated against a single modality by label. Furthermore, instead of testing enrichment or depletion using a statistical significance cutoff, one can also ask how parameters from other modalities affect or correlate with the values ​​obtained in the current modality of interest. To address this question, correlation analyses can be performed across modalities, and models that take multiple modalities into account can be constructed. For this purpose, the aforementioned tools, such as random forest classifiers, can be used for the task of predictive modeling of objects based on multimodal portraits.As mentioned above, the subsequent classifier weight fractions can also be extracted later to understand the relative influence of each parameter in each modality on the immediate predictive task.

[0067] Conveniently, the integration of these spatially resolved imaging datasets provides flexibility in analysis. The analysis pipeline can be drawn from or used for many of the independently listed imaging modalities. Considering cross-modal analysis from this perspective reveals an exciting opportunity to validate new multimodal analysis techniques, in addition to demonstrating their usefulness along with new insights.

[0068] Additional anticipated applications Pixel level calculation and analysis Most of the analyses described above focus on either identifying spatially resolved objects or clustering pixel-level data for discretization of the dataset for summarization. Alternatively, if we want to analyze registered images when they are at the pixel level, we can also collect trends of parameters of interest across organizations or regions of interest within organizations. For example, the gradient of a parameter of interest across an entire image can be visualized by calculating a parameter density estimate. The resulting smoothed representation of the pixel-level data resembles a continuous gradient and can be visualized as a contour map or heatmap. In our exemplary dataset, we achieved this visualization by calculating smoothed versions of markers of interest in the IMC data relative to each other, demonstrating 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 registration processes and spatial alignment across modalities, we can also calculate gradients across modalities. These continuous representations, when formally implemented in spatial gradient models such as in

[49] , can also be used to provide numerical solutions for the attractive and repulsive forces that intramodality or intermodality parameters have on each other. When used in conjunction with time-dependent analysis, these numerical solutions and equations offer the possibility of developing cross-modal simulation models at the tissue level. For example, if the data acquisition sensitivity required to identify single molecules with high confidence is provided in MSI, our dataset can also be combined with known attractive and repulsive forces between single molecules to simulate biological processes in tissues.

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

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

[0070] Multi-domain conversion Recently, advances in computer vision and artificial intelligence algorithms have continued in both classification tasks and generative modeling. Of particular note are models capable of learning and generating underlying distributions to construct / represent immediate datasets, such as adversarial generative networks

[44] and adversarial autoencoders

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

[46] . From this perspective, any modality can be transformed from one to the other, provided that a training relationship exists between them. This process, which we consider to be unantibody labeling, is an extension of the application of generative modeling in biological image prediction, such as when generating IMC images from a trained generative model on another modality[47,48].

[0071] The following examples are intended to illustrate the present invention. They are not intended to limit the present invention in any way. [Examples]

[0072] Example 1. Multimodal imaging and analysis of diabetic foot ulcer tissue. Diabetic foot ulcer (DFU) biopsies, including completely excised ulcers and surrounding healthy margin tissue, were performed, followed by tissue processing in preparation for multimodal imaging. Serial sections of the DFU biopsies were imaged using matrix-assisted laser desorption / ionization (MALDI) mass spectrometry (MSI), imaging mass cytometry (IMC), and optical 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 thin section of the DFU biopsy was stained with hematoxylin and eosin (H&E) and imaged using bright-field microscopy scanning. To prepare DFU biopsy thin sections for MSI (Figure 2A), the thin sections were sprayed with matrix solutions (optimized for each type of analyte). In this example, a matrix containing 40% 2,5-dihydroxybenzoic acid (DHB) in a 50:50 v / v acetonitrile:0.1% TFA aqueous solution was used (Figures 2B and 2C) to preferentially image low molecular weight molecules and lipids. 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 thin section (Figure 2D). Imaging mass cytometry was performed on regions of interest within the DFU biopsy thin section imaged by H&E staining and MSI. Following tissue or cell culture pretreatment, samples were stained with metal-labeled antibodies (Figure 3). Subsequently, labeled molecular markers in the samples were ablated using an ultraviolet laser connected to the mass cytometer system (Figure 3). In a mass cytometer, sample cells were vaporized, atomized, ionized, and filtered through a quadrupole ion filter. 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).

[0073] Example 2. Processing and analysis of multimodal and high-dimensional data Multimodal imaging data acquired using any combination of modalities, including MSI, IMC, immunohistochemical testing (IHC), and H&E staining, were processed using an integrated analysis pipeline (Figure 4). The analysis pipeline was designed as a generalizable framework for identifying high-value or actionable indicators (e.g., biomarkers or prognoses) consisting of one or more parameters uniquely revealed through the creation and analysis of multidimensional maps by correlating spatially resolved datasets from a wide range of diverse origins (e.g., laboratory samples, various imaging modalities, geographic information system data) with other aligned data. A series of steps were employed to process the multimodal imaging data in order to create such multidimensional maps. First, spatial alignment of all modalities was performed in a process called image registration (Figure 4). Steps 2-4, (2) image segmentation, (3) clustering and annotation based on manifolds at the pixel level, and (4) extraction and analysis of multimodal data features, were performed in parallel and were complementary approaches used to identify trends in the expression or abundance of parameters of interest for modeling and predicting biological processes at multiple scales: cellular niches (fine local contexts), local tissue heterogeneity (local population contexts), features of overall tissue heterogeneity and trending (broad contexts), and disease / tissue states (combinations of local and broad tissue contexts).

[0074] Example 3. Comparison of execution times and estimation of data dimensions using multiple dimensionality reduction methods. To develop a rapid and accurate method for (1) identifying non-healing diabetic foot ulcers (DFUs) at the time of presentation, and (2) evaluating the effectiveness of debridement procedures in DFU wound healing, we characterized the execution time of multiple dimensionality reduction methods on multimodal and high-dimensional imaging MSI datasets. To condense the dimensionality of the MSI datasets, we used the following dimensionality reduction techniques: homogeneous manifold approximation and projection (UMAP), isometric mapping (Isomap), t-distribution stochastic neighborhood embedding (t-SNE), PHATE (potential of heat diffusion for affinity-based transition embedding), principal component analysis (PCA), and non-negative matrix factorization (NMF) (Figure 5). The intrinsic dimensionality of the MSI data was estimated for each dimensionality reduction method (Figure 5). Embedding errors as mean and standard deviation (n=5) were plotted as 1- to 10-dimensional functions for all methods. Convergence to the embedding error values ​​indicated that the algorithm's ability to incorporate data complexity no longer improved as the dimensionality of the resulting embeddings increased. The inventors observed that nonlinear dimensionality reduction methods, such as t-SNE, UMAP, PHATE, and Isomap, converged to a much lower intrinsic dimension than linear methods, such as NMF and PCA, indicating that far fewer dimensions are needed to accurately represent the dataset. The computational execution time for each algorithm was measured and plotted as mean execution time and standard deviation over each dimension (Figure 6). The nonlinear methods, t-SNE and Isomap, required longer execution times than the nonlinear methods PHATE and UMAP. While the linear methods required minimal execution time, they also failed to succinctly capture data complexity. The results demonstrate that the UMAP algorithm, a manifold-based nonlinear technique, provides the best MSI data representation compared to other methods, based on standard image registration and computational complexity testing, robustness to noise, and its ability to capture information with low-dimensional embeddings.

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

[0076] Example 5. Dimensionality Reduction Process Pipeline Dimensionality reduction using UMAP was performed on a DFU biopsy MSI dataset (Figures 8-9). Each UMAP dimension in the 3D embedding was pseudo-colorized, for example, dimension U1 in red, dimension U2 in green, and dimension U3 in blue (Figure 9). The overlap of the three channels resulted in a synthetic grayscale image 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 to 3D in this example (center panel). Assigning embedding dimensions to arbitrary colors allows for better visualization of the data projection along the 3D plane. After UMAP 3D embedding, each pixel of the dataset, now color-coded by its UMAP dimension, can be mapped back to its original position on the DFU image (right panel). This enables the visualization of arbitrary structures related to collected tissue sections in high-dimensional datasets.

[0077] Example 6. Relative evaluation of the noise robustness of selected dimensionality reduction methods. Linear dimensionality reduction methods, such as NMF and PCA, suffer from the problem of overestimating the intrinsic dimension of the data and are susceptible to noisy channels. Dimensionality reduction was performed using linear and nonlinear methods, and the first two dimensions of the 4-dimensional embedding for each method were visualized (Figure 10). Linear methods require a larger number of features to capture the complexity of the dataset, and often the captured features are confused by noise, with some features dedicated solely to representing the noise. To further evaluate the noise-tolerant properties of dimensionality reduction methods—nonlinear, e.g., t-SNE, UMAP, PHATE, and Isomap, as well as linear, e.g., NMF and PCA—the manifold structure of full mass spectrometry imaging (MSI) data (noisy) and denoised MSI data (peak-picked) was characterized using the Denoised Manifold Preservation (DeMaP) metric. The Euclidean distance of the obtained embedding corresponding to the noisy MSI data and the geodesic distance of the corresponding peak-picked data, along with the DeMaP metric between them, were calculated. The mean and standard deviation DeMaP metric for all dimensionality reduction methods tested were plotted over 1 to 10 dimensions (Figure 7C).

[0078] Example 7. Multiscale Image Registration Pipeline We first performed a multiscale iterative registration approach (referred to as global registration) to spatially align the multimodal image dataset at the whole tissue level, followed by higher-resolution registration in a subset region (ROI) of interest (referred to as local registration). The spatial resolution of the imaging modalities varied widely among them. For example, MSI resolution was approximately 50 μm, H&E and toluidine blue resolution approximately 0.2 μm, and IMC resolution approximately 1.0 μm (Figure 11). To preserve the spatial coordinates of high-dimensional, high-resolution structures and tissue morphologies during multimodal image registration, we maintained a high-resolution image that served as a constant reference image in each step of the registration scheme, and aligned all other images against it.

[0079] Global grayscale image registration of DFU biopsy tissues imaged with MSI, H&E staining, and toluidine blue staining was performed using a multi-step process with the Elastix registration toolkit (Figure 12). First, the MSI images were processed for dimensionality reduction using UMAP. The resulting MSI images after UMAP dimensionality reduction (referred to as MSI0) were registered to their corresponding H&E0 images to generate transformed MSI1 images (Figures 12 and 13). This transformation (T1) warps the MSI images while maintaining the fixed H&E images. As a result, transformed MSI images (MSI1) are obtained, which are aligned to the H&E images. In parallel, the H&E0 images were registered to their corresponding toluidine blue0 images, which were separate adjacent tissue sections from the same DFU biopsy, and these were used for IMC imaging. Toluidine blue0 contained the spatial coordinates of the IMC region of interest, which served as reference coordinates for subsequent local transformations of the images. This transformation (T2) warps the H&E image while maintaining the fixed toluidine blue image. Finally, transformation T2 is applied to the already transformed MSI1 to obtain an MSI image (MSI2) registered with toluidine blue 0. This process can be summarized in two equations: T MSI-fT1 is 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 However, T2 is the final transformed H&E image used for downstream analysis, and as mentioned above, T2 is the registration transformation from the H&E image to the toluidine blue (IMC) image. H&E-f =T2.

[0080] After spatially aligning images from all modalities at a global level, the inventors incorporated a secondary fine-tuning step for image alignment for smaller ROIs (Figure 13). Due to the destructive nature of IMC imaging, it is necessary to add spatial information about the imaged sample using a reference image of the same sample acquired before IMC. The reference image is obtained from a toluidine blue-stained image and provides the ability to compensate for local deformation occurring in the tissue sample within each ROI. Fine-tuning of registration at a local scale was performed by selecting regions of interest in the toluidine blue image corresponding to each MSI image. Overall registration for a single ROI is performed by appropriate (modality-dependent) sequential transformations, first at a global level and then by local transformations (Figure 14).

[0081] Example 8. Feature extraction and analysis of multimodal data Spatially aligned images from multimodal datasets were analyzed to identify objects in a process called segmentation. Once the spatially resolved objects were identified, the inventors began comparing data from different modalities concerning those objects using their spatial coordinates. The inventors compared features from registration images containing data from IMC datasets (used to identify single-cell coordinates due to their relatively high spatial resolution) and MSI datasets (images A-C and A''-C'' in Figure 15). Subsequently, data were clustered based on single-cell measurements for both the IMC modality alone and the MSI modality alone. Cell types were determined using the clustering of IMC single-cell measurements (images A'-C' and A'''-C''' in Figure 15). The ability to integrate multiple imaging modalities allowed for sorting tests to be performed as a function of the corresponding cell types defined in the IMC dataset regarding the enrichment or depletion of certain features in the MSI modality.

[0082] Example 9. Multi-omics image alignment and analysis using information manifolds (MIAAIM) MIAAIM is a sequential workflow aimed at providing a comprehensive portrait of tissue states. It includes four processing stages: (i) image preprocessing using a high-dimensional image creation (HDIprep) workflow, (ii) image registration using a high-dimensional image registration (HDIreg) workflow, (iii) tissue state transition modeling using cobordism approximation and projection (PatchMAP), and (iv) cross-modality information transfer using 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 standardization format of the assembled images vary depending on the technology. For example, methods based on cycle fluorescence (e.g., CODEX, CyCIF) assemble BioFormats / OME-compliant 20-60 plex whole tissue mosaic images after illumination uniformity correction (e.g., BaSiC) and tile stitching (e.g., ASHLAR); other methods acquire 20-100 plex data directly in a region of interest (ROI) (e.g., MIBI, IMC). Further methods quantify thousands of parameters at rasterized locations on the whole tissue or ROI and do not store them in a BioFormats / OME-compliant format. For example, the imzML format, which relies on the mzML format used by Human Proteome Organization, often stores MSI data.

[0083] Regardless of the technique, assembled images contain numerous heterogeneously distributed parameters, which makes comprehensive, manually guided image alignment impossible. In addition, high-dimensional imaging generates a large feature space that makes methods typically used in unsupervised settings difficult. The HDIprep workflow in MIAAIM produces compressed images that preserve multiple prominent features, enabling inter-technical statistical comparisons while minimizing computational complexity (HDIprep, Figure 16). For images acquired from tissue staining, HDIprep provides parallelized smoothing and morphological operations that can be applied sequentially to preprocessing. Image registration with HDIreg generates transformations to combine modalities within the same spatial domain (HDIreg, Figure 16). HDIreg uses Elastix, a parallelized image registration library for calculating transformations, and is optimized to transform large multi-channel images with minimal memory usage while assisting tissue staining. HDIreg automatically resizes, pads, and trims boundaries before applying image transformations.

[0084] Aligned data are well-suited for established single-cell and spatial neighborhood analyses—they can also be segmented to incorporate multimodal single-cell measurements (levels 3 and 4 data), such as mean protein expression or spatial features of cells, and can be analyzed at the pixel level. However, a common goal in pathology is to map the transition from health to disease using synthetic tissue portraits. Similarity between systems—at the tissue state level—can be visualized in the PatchMAP workflow (PatchMAP, Figure 16). PatchMAP models tissue states as smooth manifolds that are stitched together to form a higher-order manifold called a cobordism. It is, as a result, a nested model that incorporates nonlinear intra-system states and inter-system continuities. Applying this paradigm as a tissue-based atlas mapping tool, information can be transferred between modalities by i-PatchMAP (i-PatchMAP, Figure 16).

[0085] The MIAAIM workflow is nonparametric, using manifold-supported probability distributions rather than training data models. Therefore, MIAAIM is technology-independent and generalizable to multiple imaging systems (Table 1). However, nonparametric image registration is often an iterative parameter tuning process rather than a “black box” solution. This presents a significant challenge to reproducible data integration across institution and computing architectures. Therefore, we containerized the MIAAIM data integration workflow with Docker, developed a Nextflow implementation, documented the human-in-the-loop process, and removed language-specific dependencies in accordance with FAIR (Searchable, Accessible, Interoperable, and Reusable) data stewardship principles.

[0086] [Table 1]

[0087] High-dimensional image compression with HDIprep. To compress high-parameter images, HDIprep performs dimensionality reduction on pixels using homogeneous manifold approximation and projection (UMAP) (Figure 17A). The inventors conducted a rigorous comparison using novel imaging datasets of diverse tissue biopsies with high cellular complexity, 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 preservation algorithms in terms of its robustness to noise and its ability to efficiently preserve data complexity while incorporating morphological structure (Figures 18A-18J, 19A-19H, and 20A-20H).

[0088] HDIprep preserves global data complexity with the fewest degrees of freedom required by detecting stationary state manifold embeddings. To determine the stationary state dimension, information captured by UMAP pixel embeddings is computed over a wide range of embedding dimensions (cross-entropy, definition 1, method), and the first dimension in which the observed cross-entropy approaches the asymptote of the exponential regression fit is selected. The calculation of the stationary state embeddings corresponds to the number of pixels and their squares, and thus HDIprep embeds spectral landmarks into a pixel manifold representing their global structure (Figures 21A and 21B).

[0089] Pixel-level dimensionality reduction is computationally expensive for large images, i.e., high resolution (e.g., 1 μm / pixel). To reduce compression time while preserving quality, we developed a subsampling scheme in which a pixel subset of the spatial representation is embedded before spectral landmark selection, and out-of-sample pixels are projected onto the embedding (Figures 22A, 22B, 23A, 23B, 24A, and 24B). HDIprep also scales all optimizations to full tissue images, combined with a recent neural network UMAP implementation. We have demonstrated its effectiveness on publicly available 44-channel CyCIF images containing approximately 100 million and 256 million pixels (Figure 25). Thus, HDIprep presents an objective pixel-level compression method applicable to multiple modalities (Algorithm 1, Method).

[0090] High-dimensional image registration (HDIreg). MIAAIM connects the HDIprep and HDIreg workflows with a manifold alignment scheme parameterized by spatial transformations. The inventors developed a theory for calculating manifold α-entropy using an entropy graph for UMAP embeddings and applied it to image registration using Reny α-mutual information (α-MI) based on the entropy graph (HDIreg, methods). HDIreg generates a transformation that maximizes the α-MI from image to image (manifold to manifold) (Figure 17B). This image similarity measure is generalized to Euclidean embeddings of arbitrary dimensions by considering the distribution of k-nearest neighbor (KNN) graph lengths of compressed pixels, rather than directly comparing the pixels themselves. Combining HDIprep compression with KNN α-MI extends intensity-based registration to complex images without the need for corresponding contrast staining between techniques.

[0091] Proof of Principle 1: MIAAIM generates information on cellular phenotype, molecular ion distribution, and tissue state at the full scale. To highlight the usefulness of high-dimensional image integration, the inventors applied the HDIprep and HDIreg workflows to MALDI-TOF MSI, H&E, and IMC data from DFU tissue biopsies containing a wide range of tissue states, from necrotic centers to healthy margins of ulcers. Image acquisition was performed at 1.2 cm for H&E and MSI data. 2 This covered the following areas. MSI molecular imaging enabled non-targeted mapping of lipids and small metabolites across the entire sample in the 400–1000 m / z range with a resolution of 50 μm / pixel. Tissue morphology was acquired by H&E at 0.2 μm / pixel, and 27 plex IMC data were obtained from seven ROIs on adjacent sections at a resolution of 1 μm / pixel.

[0092] Cross-modality alignment was performed in a global to local manner (Figure 17C). The inventors used HDIprep compression for high-parameter data and HDIreg manifold alignment for registration of the compressed images. Destructive properties of IMC acquisition in small ROIs.2 Consequently, the inventors first aligned complete tissue data from MSI, H&E (downsampled to approximately 3.5 μm / pixel), and IMC reference images. Local deformations not captured at the full tissue scale within each ROI were corrected using manual landmark guidance. Serial sectioning deformations were considered by nonlinear transformations. Registration was initialized by performing an affine transformation on the coarse alignment before nonlinear correction. Resolution differences were considered by a multi-resolution smoothing scheme. The final alignment was carried out by configuring both modality and ROI-specific transformations.

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

[0094] Proof of Principle 2: Identification of molecular microenvironment niches correlated with cells and pathological conditions through multiple omics networking. The inventors confirmed the existence of cross-modal relationships from Proof of Principle 1 by performing microenvironment correlation network analysis (MCNA) on registered IMC and MSI data (Figures 26A-26I). The inventors performed community detection (i.e., clustering) on ​​MSI analytes (m / z peaks) based on their correlation with single-cell protein measurements and defined microenvironment correlation network modules (MCNM; different colors in Figure 26A). Examination of MCNMs with high correlations to protein levels identified in IMC revealed that molecular aggregates, rather than individual peaks, are related to cellular protein expression (Figure 26B). MCNMs were grouped on separate axes based on whether they showed a moderate positive correlation with cellular markers indicating inflammation and cell death (CD68, activated caspase-3), or with immunomodulatory markers (CD163, CD4, FoxP3) and vascular markers (CD31). Several proteins, such as CD14 (myeloid cell marker) and the cell proliferation marker Ki-67, did not show a strong correlation with any m / z peak in any cell type.

[0095] To gain insights into the relationship between molecular distribution and tissue health, the inventors plotted the ion intensity distribution of MCNM against proximity to the ulcer center (Figure 26C). This analysis revealed a shift in the molecular profile approximately 6 mm from the ulcer center as the tissue state progressed from healthy to damaged. The inventors validated the performance of HDIreg in aligning their observations and microstructures by visualizing the distribution of top correlated ions within the cellular microenvironment (Figures 26D and 26E).

[0096] The advantage of our analysis lies in its potential to identify molecular variations in one modality (here, MSI) that correlate with cellular states identified using different modalities (here, IMC). We investigated whether the m / z peak was distinctly associated with cell proliferation (the Ki-67 marker in IMC). We performed unsupervised clustering on segmented cellular-level expression patterns in IMC to identify cellular phenotypes (Figure 26F) and performed differential correlation network analysis between phenotypes within well-isolated CD3+ clusters (likely to identify infiltrating T cells at the wound site) and CD3- cell populations (Figure 26G). Interestingly, we found that the correlation with Ki-67 expression shifted by near significance (2σ) for several m / z peaks between the CD3- and CD3+ populations (Fischer-transformed one-sided z-statistic; Bonferroni-corrected P-value) (Figure 26H).

[0097] Subsequently, the inventors, utilizing the spatial context preserved by MIAAIM, observed that the ionic 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 primarily occurs near the healthy margin of the DFU, confirming that molecular correlates of T-cell proliferation can be identified through this unbiased analysis. In summary, these results provide insights into the molecular microenvironments associated with the different functional and metabolic states of specific cell subtypes, as well as insights into how these microenvironments are distributed in spatial context along the gradient from damaged tissue to healthy tissue.

[0098] Mapping of tissue state transitions via cobordism approximation and projection (PatchMAP). To model transitions between tissue states such as healthy or damaged, the inventors generalized manifold learning and dimensionality reduction to higher-order situations by developing a novel algorithm called PatchMAP, which integrates mutual neighbor calculation with UMAP (Figure 27A and Algorithm 2, Method). The inventors hypothesized that the topological space of system-level transitions is nonlinear and can be consistently parameterized by manifold learning. Therefore, PatchMAP represents the non-combined manifolds (i.e., system states) as boundaries of higher-dimensional manifolds called cobordisms (i.e., state transitions). Overlapping patches are connected by pairwise dependent neighbor queries representing cobordism geodesics between bounded manifolds, and their metrics become adaptable by stitching them together using the t-norm. Interpreting and executing PatchMAP embeddings is similar to existing dimensionality reduction algorithms—similar data within or between bounded manifolds are located close to each other, while heterogeneous data are located further apart. PatchMAP incorporates both the topological structure of bounded manifolds and the continuity across the entire bounded manifold to generate cobordism.

[0099] Currently, there is no method for forming cobordisms—the closest approach to achieving this is a 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 a stitching simulation using the “digit” machine learning method development dataset (Figure 27B). We partitioned the data by labeling bounded manifolds and then re-stitched the complete dataset using each method. In this task, complete stitching resulted in complete separation of projected bounded manifolds, which we quantified using silhouette coefficients (SC). For control of visualization, UMAP was used after data integration to provide PatchMAP with embeddings similar to all benchmark methods.

[0100] PatchMAP is robust to overlapping bounded manifolds and outperforms data integration methods at higher neighborhood (NN) counts. All other methods inaccurately blend bounded manifolds when there was no overlap, which is expected if we consider that their assumptions are broken by the lack of connections between manifolds. PatchMAP's stitching, on the other hand, uses fuzzy intersection, which strongly weights accurate connections while discarding inaccurately connected data across the manifold. We also verified whether PatchMAP preserves bounded manifold structure while embedding higher-order structures 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 annotation is strongly weighted, resulting in less blending and better manifold separation.

[0101] Imaging technology and inter-tissue information transfer (i-PatchMAP). The inventors assumed that information transfer between biological states should similarly describe continuous transitions and be robust to the strength (including the absence) of manifold connections. Therefore, the i-PatchMAP workflow uses PatchMAP as a paired method of domain transfer and quality control visualization to propagate information between different samples (information transfer, Figure 27A). To do this, i-PatchMAP first normalizes the connections between bounded manifolds of "reference" and "query" data to define local one-step Markov chain transition probabilities (transition probabilities, Figure 27A), and then linearly interpolates the measurements from the reference data to the query data (information transfer, Figure 27A). Quality control of i-PatchMAP can be performed by visualizing the connections between bounded manifolds in the PatchMAP embedding (visualization of manifold connections, Figure 27A).

[0102] To benchmark i-PatchMAP, the inventors compared i-PatchMAP with other nonparametric domain transfer tools, Seurat v3, and a modified version of UMAP (UMAP+) that incorporates transition probability-based interpolation similar to i-PatchMAP, as well as data from Proof of Principle 1 and the publicly available umbilical cord blood mononuclear cell (CBMC) CITE-seq dataset. 11This was compared to the methods used. UMAP+ utilizes a directed NN graph from query data to reference data for data interpolation, rather than the metric-fit stitching of PatchMAP. Therefore, it acts as a control mechanism for PatchMAP. The inventors constructed 23 evaluation examples by tiling ROIs from Proof of Principle 1 and predicted IMC protein expression by performing one-out cross-validation using single-cell MSI profiles. The inventors assessed the accuracy for each parameter using Spearman correlation between predicted spatial autocorrelation (Moran's I) and true autocorrelation, taking into account differences in resolution between imaging modalities. i-PatchMAP was superior to the methods tested in its ability to transfer IMC measurements to query data based on MSI profiles (Figure 27B), however, all methods performed were consistently insufficient for parameters (TGF-β, FoxP3, CD163) where the original spatial autocorrelation was not present in the tile. For the CITE-seq dataset, the inventors prepared 15 evaluation cases and predicted the abundance of antibody-derived tags (ADTs) using single-cell RNA profiles. The inventors quantified performance using Pearson correlation between true and predicted ADT values ​​(Figure 27C) and found that the i-PatchMAP method performed better or slightly better than other testing methods for all parameters.

[0103] Proof of Principle 3: i-PatchMAP transfers multi-protein distributions between tissues based on molecular microenvironment profiles. To evaluate whether molecular signature information can be transferred between imaging modalities and even between different tissue samples using i-PatchMAP, we extrapolated IMC information, based on the MSI profile, to complete DFU samples and similarly to separate prostate tumor and tonsil specimens, using single-cell IMC / MSI protein measurements (see Proof of Principle 1). PatchMAP embedding of single cells and individual pixels into DFU ROIs across the tissue based on MSI parameters revealed that the single-cell molecular microenvironment in the DFU ROI provides a good representation of the overall DFU molecular profile (Figure 27F). Therefore, we used i-PatchMAP to transfer DFU single-cell protein measurements to complete DFU tissues based on molecular similarity. i-PatchMAP predicted that the wound areas of the DFU tissue would show high expression levels for CD68 (a marker of pro-inflammatory macrophages) and activated caspase-3 (a marker of apoptotic cell death). On the other hand, the healthy margin of DFU biopsies was predicted to contain higher levels of CD4 (indicating invasive T cells) and the cell proliferation marker Ki-67. Interestingly, PatchMAP visualization revealed that the molecular microenvironment corresponding to specific single-cell scales (e.g., CD4) in DFUs had a strong correlation with MSI pixels in tonsil tissue (Figure 27F). In tonsil tissue, which is lymphocyte-rich, i-PatchMAP predictions for CD4 agreed well with lymphocyte structure, and regions lacking cellular contents were accurately predicted to not contain CD4. On the other hand, there was no strong correlation between the molecular profiles of prostate cancer samples and DFU biopsies. Therefore, in the current dataset, strong cellular and molecular correlations between samples are thought to be supported by the common presence of specific immune cell populations. Indeed, IMC examination of the prostate biopsies used here showed insufficient immune cell infiltration.

[0104] method MIAAIM implementation. By implementing the MIAAIM workflow in Python and connecting it via the Nextflow pipeline language, automated result caching, dynamic resumption of processing after changes in workflow parameters, and efficient parallel processing of multiple images are enabled. MIAAIM is also available as a Python package. Each data integration workflow is containerized to enable a reproducible environment and eliminate arbitrary language-specific dependencies. MIAAIM output interfaces with numerous existing image analysis software tools (see Appendix 1, Combinations of MIAAIM with Existing Bioimaging Software). Therefore, MIAAIM complements existing tools rather than replaces them.

[0105] High-dimensional image compression and preprocessing (HDIprep). HDIprep is implemented by specifying a series of 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 the Python NiBabel library. NIfTI-1 was chosen as the default file format for many MIAAIM operations due to its compatibility with Elastix, ImageJ for visualization, and its memory mapping capabilities in Python.

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

[0107] Definition 1. Given a reference set A, membership functions u: A → [0, 1], v: A → [0, 1], the fuzzy set intersection entropy C of (A, u) and (A, v) is defined as TIFF0007881548000002.tif12128.

[0108] The fuzzy set intersection entropy is a global measure of the agreement between simplicial sets, aggregated among members of the reference set A (here, graph edges). Calculating its exact value scales with the number of data points and is of the order of the square, limiting its use for large datasets. Therefore, the current implementation of UMAP cannot compute the exact intersection entropy during that optimization of the low-dimensional embedding. Instead, this relies on probabilistic edge sampling and negative sampling to reduce the runtime for large datasets. 35 . Together, to identify the steady-state embedding dimension, the inventors compute patches on the data manifold representing its global structure and then use these, after projecting them across a wide range of dimensions using UMAP, for the calculation of the exact intersection entropy. As a result, a global estimate of the dimension required to accurately capture the manifold complexity is obtained.

[0109] To identify patches of global representation on a data manifold, the inventors subjected a fuzzy simplicial set to a variation of spectral clustering. Using UMAP, the inventors repeatedly projected spectral centers onto a dimensionally increasing Euclidean space, calculated the fuzzy set cross-entropy for each case, and then minimax-normalized the obtained values. To identify the stationary embedding dimension, the inventors fitted least-squares exponential regression to the normalized cross-entropy depending on the dimension, and then simulated the sample along the regression line to confirm a first dimension that falls within the 95% confidence interval of the exponential asymptote. Subsampled data is embedded in the stationary dimension, and out-of-sample pixels are projected onto this embedding using a method based on UMAP's native neighborhood (transform function). Finally, all pixels are mapped back to their original spatial coordinates, and a compressed image is constructed with a number of channels equal to the stationary embedding dimension. These steps are summarized in the following pseudocode.

[0110] Algorithm 1: Image Compression Input: Multichannel image (X), SVD dimension (b), k-mean clusters (k), embedding dimension (n) Output: Compressed image (I) Execution of the compression function TIFF0007881548000003.tif169149

[0111] 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 mask standards for sampling regions, which can be useful for very large datasets.

[0112] By default, images with fewer than 50,000 pixels are not subsampled; images with 50,000 to 100,000 pixels are subsampled using 55% pseudo-random sampling starting from a 2x2 grid with uniform spacing; images with 100,000 to 150,000 pixels are subsampled using 15% pseudo-random sampling starting from a 3x3 grid; and images with more than 150,000 pixels are subsampled using a 3x3 grid. These default values ​​are based on empirical studies (Figures 22A, 22B, 23A, 23B, 24A, and 24B).

[0113] Subsampling was not used for the presented MSI data. The subsampling rate used for the presented IMC data was determined on a case-by-case basis from empirical studies and matched that used in spectral landmark sampling experiments. Subsampling using a 10x10 pixel grid with uniform spacing was used for CyCIF data compression.

[0114] Generation of Fuzzy Simplices. To construct a pixel-level data manifold, the inventors 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). The inventors then implement the UMAP algorithm to extract a set of fuzzy simplices representing the manifold structure of these resulting d-dimensional points. For all presented results, the inventors used default UMAP parameters to generate this manifold: 15 nearest neighbors and Euclidean metric.

[0115] Manifold landmark selection by spectral clustering. Spectral landmarks are identified using a variation of spectral clustering. The inventors scale spectral clustering to a large dataset using randomized singular value decomposition (SVD) and then minibatch k-means, following the procedure introduced in the PHATE (potential of heat diffusion for affinity-based transition embedding) algorithm. d-dimensional space Given a symmetric adjacency matrix A representing the pairwise similarity between nodes (in this case, pixels) originating from TIFF0007881548000004.tif3128, the inventors first calculate the eigenvectors corresponding to the largest eigenvalue k of A. Then, using these eigenvectors k as functions, the inventors perform a mini-batch k-average on the nodes of A. Subsequently, spectral landmarks are defined as the d-dimensional centers of the obtained clusters.

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

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

[0118] Tissue image preprocessing. HDIprep processing options for hematoxylin-eosin (H&E) stained tissues and other low-channel tissue stains include image filtering (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, and then a binary mask representing the foreground was created by Otsu thresholding. Subsequently, sequential morphological operations were applied to the mask, including morphological opening to remove small bound foreground components, morphological closing to fill small holes in the foreground, and filling to close large holes in the foreground.

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

[0120] High-Dimensional Image Registration (HDIreg). HDIreg is a containerized workflow that implements open-source Elastix software along with a custom Python module for automating image resizing, padding, and cropping, which are often applied before registration. HDIreg incorporates several different registration parameters, cost functions, and deformation models, and additionally allows for manual definition of point correspondences for difficult problems, as well as configuration of transformations for fine tuning (see Supplementary Document 2, Notes on Predictive Performance of the HDIreg Workflow).

[0121] High-parameter images are registered using a manifold alignment scheme parameterized by a spatial transformation aimed at maximizing image similarity. Formally, the inventors have optimized the registration process as follows: 40 Let's consider it as such.

[0122] Domain Ω F A fixed d-dimensional image I F : TIFF0007881548000005.tif4128 and Domain Ω m A dynamic q-dimensional image I M : Given TIFF0007881548000006.tif3128, the inventors, The goal is to optimize TIFF0007881548000007.tif8128. In the formula, T μ : TIFF0007881548000008.tif4128 is a parameter A smoothing transformation defined by the vector TIFF0007881548000009.tif4128, where S is TIFF0007881548000010.tif5128 and I F This is the similarity measure that yields the maximum result when the elements are aligned.

[0123] Differential Geometry and Manifold Learning: The manifold alignment scheme of MIAAIM uses the Renyi α - mutual information (α - MI) based on the entropy graph as the similarity measure S in Equation 1, which is extended to the manifold representation of images (i.e., compressed images) embedded in Euclidean spaces with potentially different dimensions. This measure is justified in the HDIreg manifold alignment scheme through the concept of intrinsic manifold information (i.e., entropy). Below, the inventors introduce basic differential geometric concepts that can extend the existing basis of intrinsic manifold entropy estimation to the UMAP algorithm.

[0124] Definition 2: Let X and Y be topological spaces. A function f: X → Y is at each point for each open neighborhood N of TIFF0007881548000011.tif3128 and TIFF0007881548000012.tif4128, if the set f -1 (N) is an open neighborhood of TIFF0007881548000013.tif3128, then it is continuous. A function ∫: X → Y is a homeomorphism if it is one - to - one, onto, continuous, and has a continuous inverse. If there exists a homeomorphism between spaces X and Y, they are called homeomorphic spaces.

[0125] Definition 3. A manifold M of dimension n (i.e., an n - manifold) is a second - countable Hausdorff space, and each of its points has an open neighborhood that is homeomorphic to an n - dimensional Euclidean space TIFF0007881548000014.tif3128. For any open set TIFF0007881548000015.tif3128, the inventors can define a chart (φ, U), where TIFF0007881548000016.tif4128 is a homeomorphism. (φ, U) can be said to act as a local coordinate system for M, and the inventors When TIFF0007881548000017.tif3128 is non-empty, the transition between two charts (φ, U) and (ω, V) can be defined as TIFF0007881548000018.tif4128.

[0126] Definition 4. A smooth manifold is a manifold in which there exists a smooth transition mapping between each chart of M. The Riemannian metric g is a mapping that determines the inner product g TIFF0007881548000019.tif4128 at each point y between the tangent vectors to M. y ( . , . ) for the tangent vectors of y, denoted as T y M. A Riemannian manifold written as (M, g) is a smooth manifold M with the Riemannian metric g. Given a Riemannian manifold, the Riemannian volume element provides a means of incorporating a function related to volume in local coordinates. Given (M, g), the inventors can express the volume element ω in terms of the metric g at the local coordinates χ = χ1,…,χ n as TIFF0007881548000020.tif4128, where g(χ)>0 and ∧ denotes the wedge product. The volume of M in this volume form is Vol(M)=∫ M ω.

[0127] Definition 5. An embedding of a smooth n-manifold M into N is a differentiable mapping Ψ:M→N such that dΨ p :T p M→T Ψ(p) N is injective at all points TIFF0007881548000021.tif4128. Therefore, Ψ is an embedding if its derivative is injective everywhere.

[0128] Definition 6. An embedding between smooth manifolds M and N is a smooth function f: M → N, where f is the embedding, and its continuous function is an embedding in a topological space (i.e., an injective homeomorphism). A closed embedding between M and N is: TIFF0007881548000022.tif4128 is a closed embed.

[0129] (M,g) is an ambient where n≪d Let TIFF0007881548000023.tif4128 be a compact n-dimensional Riemannian manifold, Let TIFF0007881548000024.tif5128 be the set of independent, identically distributed random vectors whose values ​​are derived from the distribution supported by M. TIFF0007881548000025.tif5128 We define it as an open neighborhood of elements in TIFF0007881548000026.tif5128. The purpose of manifold learning is to determine that the measure of skewness D is TIFF0007881548000027.tif4128 and The goal is to approximate the embedding f such that it is minimized between TIFF0007881548000028.tif5128 and the original TIFF0007881548000028.tif5128. Therefore, the manifold learning problem is: It can be written out as TIFF0007881548000029.tif8128, and in the formula, TIFF0007881548000030.tif4128 is, Represents a family of measurable functions that can take the form TIFF0007881548000031.tif4128. In a machine learning setting, this is a vector. Open neighbors for TIFF0007881548000032.tif4128 TIFF0007881548000033.tif4128 is often defined as a geodesic distance (or its stochastic coding) approximated using a positive-definite kernel, enabling the calculation of the inner product in the Riemann framework (compared to the pseudo-Riemann framework, which requires non-positive-definite values). The measure of distortion varies depending on the algorithm (see Appendix 3, e.g., HDIprep dimensionality reduction verification). Our interest in this explanation lies in the measures induced by the embedded geodesics through the volume element of the glued-together coordinates. These provide the components necessary to quantify the intrinsic Reny-α-entropy of the embedded data manifold.

[0130] Entropy graph estimator. Lebesgue density f, and Identical distributed random vectors X1, ..., X have values ​​in the compact subset of TIFF0007881548000034.tif3128. n Given that, the external Reynolds α-entropy of f is, Given by TIFF0007881548000035.tif10128, Here, The filename is TIFF0007881548000036.tif4128.

[0131] Definition 7 (Source: Costa and Heroes) 38 ). TIFF0007881548000037.tif4128 If we consider a random vector with the same distribution whose values ​​are in a compact subset of TIFF0007881548000038.tif3128, then under the Euclidean metric... The nearest neighbor to TIFF0007881548000039.tif4128 is: This is given by TIFF0007881548000040.tif8128.

[0132] The k-nearest neighbor (KNN) graph is, This represents the edges between TIFF0007881548000041.tif4128 and its k-nearest neighbors. Let TIFF0007881548000042.tif4128 be the set of k-nearest neighbors of TIFF0007881548000043.tif4128. Then, the total edge length of the KNN graph for TIFF0007881548000044.tif4128 is given by TIFF0007881548000045.tif15128, where γ > 0 is the power weighting constant.

[0133] 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 a set of random vectors having values in a compact subset of TIFF0007881548000046.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 TIFF0007881548000046.tif3128.

[0134] Beardwood-Halton-Hammersley (BHH) theorem. Let TIFF0007881548000047.tif4128 be a compact Riemannian m-manifold embedded in TIFF0007881548000048.tif4128. Let TIFF00078815480OOO048.tif41Z8 be a set of identically distributed random vectors having values in a compact subset of TIFF0007881548000049.tif4128, and consider the Lebesgue density to be f. Assume d ≥ 2, 1 ≤ γ < d, and define TIFF0007881548000050.tif6128. Then, with probability 1, the following holds. The value that determines the right side of the limit line in Equation 6 of TIFF0007881548000051 is the extrinsic Renyi α - entropy given by Equation 4. When the i.i.d. random vectors are restricted to a compact and smooth m - manifold M in ambient In TIFF0007881548000052, when restricted to a compact and smooth m - manifold M in 4128, the BHH theorem is generalized to The intrinsic Renyi α - entropy of the multivariate density f with respect to M defined by 10128 in TIFF0007881548000053 The estimation of 4128 in TIFF0007881548000054 becomes possible by incorporating a measure μ that is naturally induced by the Riemannian metric through the Riemannian volume element. This is formalized by the following given by Costa and Hero. g

[0135] Theorem 1 (Costa and Hero): Let (M,g) be a compact Riemannian m - manifold embedded in ambient In TIFF0007881548000055, consider it as a compact Riemannian m - manifold embedded in 4128. In TIFF0007881548000056, consider 4128 as an i.i.d. random vector on M with a bounded density f relative to the differential volume element μ induced by the metric g. Assume m≧2, 1≦γ<m, and g Define 6128 in TIFF0007881548000057. Then, with probability 1, the following holds. In Equation 23128 of TIFF0007881548000058, where β m,γ,k is a constant independent of f and (M,g). Similarly, the expected value In TIFF0007881548000059 converges to the same limit line.

[0136] ​​The quantity that determines the limit line when d'=m is the intrinsic Leini alpha entropy of f, given by equation 7. Theorem 1 is used in conjunction with the manifold learning algorithms Isomap and modified C-Isomap to estimate the essential dimension of the embedded manifold. 39 In contrast to these results, which use all pairwise geodesic approximations for each point in the dataset to estimate α-entropy, we aim to provide an analogous formula that utilizes local information contained in the data manifold, following the results of our dimensionality reduction benchmark, which show 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). In fact, the information density of the volume of a contiguous region of a family of models (i.e., the output embedding space or the family of input points) is recognized in the definition of information geometry for statistical manifold learning.

[0137] Entropy graph estimator of local information of embedded manifolds: Hereinafter, we utilize two concepts to show that the intrinsic information of multivariate probability distributions supported by embedded manifolds in Euclidean space using the UMAP algorithm can be approximated using the BHH theorem: (i) the compactness of the constructed manifold and (ii) the preservation of Riemann volume elements. We address (i) with a simple proof, and we address (ii) by providing an example of the motivation for preserving volume elements using UMAP.

[0138] Definition 8. A topological space X is compact if every open cover A of X contains a finite family of subsets that also cover X. An open cover means that the elements of A are open, and the union of the elements of A is equal to X: TIFF0007881548000060.tif4128.

[0139] Proposition 1. Let n > d, and M be ambient We consider TIFF0007881548000061.tif3128 to be a compact manifold of dimension γ (γ≦d). As a result, the projection f: The image f(M) of M below in TIFF0007881548000062.tif4128 is compact. Proof. (M,g) is ambient Let TIFF0007881548000063.tif3128 be a compact Riemannian manifold (e.g., a manifold constructed with UMAP) with metric g, and let f be from M. Let f be a projection up to TIFF0007881548000064.tif3128. Since f is a projection, it is continuous and maps compact sets to compact sets.

[0140] Proposition 1 states that the d-dimensional Euclidean projection of a compact Riemannian manifold is, We have shown that we can map values ​​into a compact subset of TIFF0007881548000065.tif3128, which is a sufficient condition for the BHH theorem. The UMAP algorithm considers a set of fuzzy simplices, i.e., manifolds, constructed from a finite extended pseudometric space (a finite fuzzy realization functor, see Definition 7). Finite means that these extended pseudometric spaces are constructed from a finite family of points. Given this finiteness condition, the compactness of a UMAP manifold can be naturally derived from Definition 8 - given an open cover on the manifold, a finite subcover can be found.

[0141] Therefore, the UMAP projection is compact according to Proposition 1. To extend the BHH theorem to the calculation of the intrinsic α-entropy of the UMAP embedding, as in Equation 7, we must show that the volume element induced through the embedding is sufficiently approximated. Note that these results are applicable to any dimensionality reduction algorithm that can clearly preserve distance in open neighborhoods when embedding a compact manifold into Euclidean space. Below, we do not provide proof that UMAP preserves distance in open neighborhoods around a point, although this would be an ideal scenario. Rather, we assume that this ideal scenario exists and describe how to find the optimal dimension for projecting data to satisfy this premise.

[0142] In contrast to global data preservation algorithms such as Isomap, which calculate all pairwise geodesic distances or approximations thereof using a landmark-based approach, UMAP approximates geodesic distances with open neighborhoods close to each point (see Lemma 2 below). Y is a random vector with values ​​constrained to be located on a compact, uniformly distributed Riemannian manifold M, where μ is the Lebesgue density. 1,…, Y n Given, sample Y derived from μ i and Y j The geodesics between [point 1] and [point 2] are probabilistically encoded using UMAP and scaled to an exponential distribution: This represents TIFF0007881548000066.tif23128, and in the formula, ρ i is the vector Y i It is the distance from to its nearest neighbor, σ i is an adaptively selected normalization coefficient. Using the terminology of Equation 2, the objective of the embedding in UMAP is given by minimizing the fuzzy simplicial set cross-entropy representing the strain D (Definition 1). Formally, sample Y i and Y j Probability distribution P that encodes the geodesic between [point 1] and [point 2]. ij and sample f (Y i ) and f (Y jA probability distribution Q that encodes the distance between ) ij Given, the inventors of the present invention use the cross-entropy loss used by UMAP, It can be represented as TIFF0007881548000067.tif12128, and in the formula, Q ij teeth, The embedded vector f (Y) is obtained by TIFF0007881548000068.tif10128. i ) and f (Y j This is a probability distribution formed from the low-dimensional positions of ), where a and b are user-defined parameters for controlling the embedding spread.

[0143] Minimizing equation 11 is not, in general, a convex optimization problem. Optimization of the family F from equation 2 is restricted to a subset rather than the complete family, and therefore, in the best case, represents a local optimum. As outlined in the HDIprep workflow for a “pseudo-global” optimization procedure, we include a larger family of measurable functions to more accurately approach the optimal embedding of geodesic distances within the open neighborhood of a vector through the identification of the stationary-state embedding dimension.

[0144] Using the notation of Equation 2, The volume element induced by TIFF0007881548000069.tif5128 is similar, σ=g(χ)dχ which is TIFF0007881548000070.tif4128 Local coordinates χ 1,…, χ n The method that minimizes the strain between the volume elements of M given, and Local coordinates by TIFF0007881548000071.tif4128 Given TIFF0007881548000072.tif4128 Similar to TIFF0007881548000073.tif4128, from which we assume that the dimension n of the local coordinates is conserved. The global optimal solution for equation 7 (i.e., P ij =Qij Under the condition ( ), the volume element is preserved: TIFF0007881548000074.tif5128. Moser has proven that for a compact manifold under consideration, there exists a volume that preserves the diffeomorphism. 53 .

[0145] To identify a manifold embedding that minimizes the distortion of the volume element, the inventors have found that increasing the dimension is equivalent to an exponential increase in the potential position of a point (i.e., a copy of the real number line). By modeling TIFF0007881548000075.tif3128 as increasing, we observed a natural increase in the dimensionality of real-valued data. The relationship between the radius of the open neighborhood of a manifold in the learning environment and its volume and density was examined by Narayan et al. 52 This is formalized in applications where density is preserved at a given embedding dimension by changing the open neighborhood radius; however, the inventors can easily extend this to volume-preserving scenarios adapted under a fixed radius to infer dimensions that satisfy the geodesic distance preservation premise. Consider the following example.

[0146] TIFF0007881548000076.tif4128, ambient Let Y be a vector of the manifold M embedded in TIFF0007881548000077.tif3128. i The k-nearest neighbor of is with radius r d sphere In TIFF0007881548000078.tif4128, proportional volume Assume that the data is uniformly distributed in TIFF0007881548000079.tif4128. The mapping f is Y i Open neighborhood TIFF0007881548000080.tif4128, radius r m m-dimensional sphere of manifold N When copied to TIFF0007881548000081.tif4128, the uniformly distributed and induced Riemann volume element V mAssume that the structure containing it is preserved. Narayan et al. 52 Accordingly, the inventors have determined that proportional Using TIFF0007881548000082.tif4128, the local radius r of the embedding space m and The original radius and the relationship between them in TIFF0007881548000083.tif4128 are related by a power law. It can be inferred that TIFF0007881548000084.tif4128 exists.

[0147] To extend this example, m and r d It is assumed that the radius is fixed in the original space, and that the radius in the embedded space is Q in equation 11. ij (It can be assumed that this is controlled by the a and b parameters that affect it). The ambient metrics of the embedding space and the original space are the same, and these can be obtained using the native UMAP method. TIFF0007881548000085.tif4128 and Assume that a geodetic distance is generated within TIFF0007881548000086.tif4128 (Lemma 2 below). Ambient metric and radius are preserved, Since it is TIFF0007881548000087.tif5128, this is, TIFF0007881548000088.tif4128 and Geodetic distance δ between points in TIFF0007881548000089.tif4128 m and δ d The relationship of the power law This means that TIFF0007881548000090.tif4128 is shown. The inventors, If we further assume TIFF0007881548000091.tif4128 (i.e., an ideal scenario in which geodesic distances between open neighborhood points are preserved), then we can solve the relationship between geodesics and dimension m. Specifically, Let's consider TIFF0007881548000092.tif4128. TIFF0007881548000093.tif4128 is δ m When substituted with, the inventors of this invention, We can verify that it is TIFF0007881548000094.tif38128, which means that the geodesics in the open neighborhood of a point in the original space are exponentially related to dimension m.

[0148] Using the power law relationship between the volumes of open neighborhoods in UMAP manifolds and their embedded counterparts, the inventors can attempt to identify a dimension m such that geodesics within open neighborhoods are conserved by exponential regression. It should be noted that geodesic distance conservation is stronger than volume conservation; volume conservation is implicitly included. As a result, the stationary-state manifold embedding provides a Euclidean dimension for approximating the manifold geodesics of vectors sampled from M, and thus the volume element of M. The KNN graph function computed in the stationary-state embedding space provides the necessary machinery for calculating the intrinsic α-entropy of the data manifold embedded in MIAAIM by applying the BHH theorem with the derived scale used in the gluing of all coordinates as in Theorem 1.

[0149] However, given the assumptions introduced in our example, it should be noted that there is no guarantee that the distance outside the open neighborhood of a point can be accurately modeled in the embedding space using UMAP. Therefore, even if Theorem 1 is applied together with UMAP by replacing the KNN graph length with one obtained using the length function of a geodesic minimum spanning tree (GMST), which is another type of entropy graph, Costa and Hero 39One should not expect to reproduce the original endogenous entropy reported by []. The main contribution of the inventors here is to combine the KNN endogenous entropy estimator with a dimensionality reduction algorithm that preserves local information. The inventors hope to extend these results to settings where two such manifolds are compared, which forms the basis of the inventors' image registration application. The estimator based on the entropy graph of α-MI in the image registration setting is described below 40 .

[0150] Let \(z(\chi\) i ) = [z i (\chi\) i ), … , z d (\chi\) i )] be a d-dimensional vector encoding the features of the point \(\chi\) i . Let \(Z f (\chi)=\{z f (\chi_1), … , z f (\chi\) N )\} be the set of features of the fixed image, and \(Z m (T μ (\chi))=\{z m (T μ (\chi_1)), … , z m (T μ (\chi\) N ))\} be the set of features of the transformed video at the point \(T μ (\chi)\). Let TIFF0007881548000095.tif6128 be the concatenation of the feature vectors of the fixed and transformed videos at \(\chi\) i . As a result, TIFF0007881548000096.tif21132 is an estimator based on the graph for α-MI, where \(\gamma = d(1 - \alpha)\), \(0 < \alpha < 1\), and three graphs TIFF0007881548000097.tif55128 is the Euclidean graph function (length) from the feature vector \(z\) to its \(p\)-th nearest neighbor over the \(k\) considered nearest neighbors.

[0151] The Reny α-MI scale provides a quantitative measure of the correlation between the intrinsic structures of multiple manifold embeddings constructed with the UMAP algorithm. The Reny α-MI scale is extended to feature spaces of any dimension, from which MIAAIM, in conjunction with its image compression method, quantifies the similarity between stationary-state embeddings of image pixels in potentially different dimensions.

[0152] Proof of Concept study. Since data acquisition involved removing tissue details in the region of interest on the IMC complete tissue reference image, the inventors first aligned the complete tissue sections and then extracted data from all modalities for fine-tuning using the coordinates of the IMC region. Alignment was propagated between imaging modalities using a custom Python script. Manual landmark matching was used where unsupervised alignment proved suboptimal. Alignment errors around the IMC region were addressed following complete tissue registration by padding areas with extra pixels before cropping. All registrations, including MSI or IMC data, were performed using KNNα-MI with nearest neighbors α=0.99 and 15, as shown in equation 12. For rapid processing, all registrations were performed aligning low-channel slides (IMC reference toluidine blue images and H&E) using histogram-based MI after grayscale conversion.

[0153] For complete tissue images, a two-step registration process was implemented by first aligning the images using an affine model with respect to the parameter μ vector (Equation 1), and then aligning the images using a nonlinear model parameterized by B-splines. A hierarchical Gaussian smoothing pyramid was used to account for the resolution differences between image modalities, and a probabilistic gradient descent method with random coordinate sampling was employed for optimization. The inventors further optimized the final control point grid spacing and the number of hierarchical levels for the B-spline model, adding them individually to the pyramidal smoothing of each MSI dataset to H&E alignment (Figures 18A-18J, 19A-19H, and 20A-20H). A final control point spacing of 300 pixels was used for the nonlinear B-spline registration of MSI data to corresponding H&E data, balancing correct alignment with unrealistic distortion. The inventors then visually and by inspection identified spatial Jacobian matrices with values ​​substantially deviating from 1. H&E and IMC reference tissue registration utilized a final grid spacing of 5 pixels. Similar optimizations were applied to these data for the pyramid-level number of data points. All image-registered data was exported and saved as 32-bit NIfTI-1 images. IMC data was left unconverted and maintained in 16-bit OME-TIF(F) format.

[0154] Cobordism approximation and projection (PatchMAP). PatchMAP is an algorithm that constructs smooth manifolds by patching Riemannian manifolds to their boundaries, projecting higher-order manifolds onto lower-dimensional spaces for visualization. The higher-order manifolds produced by PatchMAP can be understood as cobordisms, which are described by the set of definitions below.

[0155] Definition 9. Family of sets Given TIFF0007881548000098.tif4128 and the index set I, The non-combined sum represented by TIFF0007881548000099.tif4128 is, for each S iFor an injective function φ i :S i →This is the set with S added. The non-compound union corresponds to the coproduct of the sets.

[0156] Definition 10. Two closed n-manifolds M and N are: The non-commutating union represented by TIFF0007881548000100.tif4128 is a coboldant if it is the boundary of some manifold W. The inventors call the manifold W a coboldism. The boundary of an n-manifold is the upper half This refers to the set of points on M that are topologically isomorphic to TIFF0007881548000101.tif4128. The inventors denote the boundary of W as ∂W.

[0157] PatchMAP addresses cobordism learning in a semi-supervised manner, assuming that the data follows a nonlinear cobordist structure. Our task is to generate a cobordism by gluing a lower-dimensional manifold to the boundary of a higher-dimensional manifold. Here, we want the coordinate transformation across the cobordism to have its own geometry, independent of the metric of the bounded manifold. In fact, this property allows us to explore data within the bounded manifold without depending on the specific geometry of the cobordism. Ultimately, the cobordist geodesic is a fundamental component of downstream applications such as the i-PatchMAP workflow. Furthermore, we want the cobordism to highlight bounded manifolds where points overlap with high confidence—such overlaps can give rise to interesting nonlinearities in higher-order spaces. A natural way to satisfy both of these conditions is to use the fuzzy set theoretical basis of the UMAP algorithm.

[0158] After that, the main goal of PatchMAP is to identify a smooth manifold whose boundary is the disjoint union of lower-dimensional smooth manifolds and which has a metric that is independent of the metrics of each of the manifolds with boundary that we have chosen to represent. We address this by a two-step algorithm that decouples the calculation of the manifolds with boundary from the cobordisms. First, we calculate the manifolds with boundary by applying the UMAP algorithm to each dataset with a user-provided metric. In fact, as a result of this step, we obtain a symmetric weighted graph that represents geodesics within each manifold with boundary. Our task is to construct a manifold (M, g) from a finite set F = {(M1, g1) … (M n , g n )} such that the metric g is F TIFF0007881548000102.tif4128 non-intersecting TIFF0007881548000102.tif4128. We approximate the geodesics of M F and TIFF0007881548000103.tif4128 for each element, we wish to approximate the tangent space T p M p F with the inner product g F

[0159] Lemma 2 (McInnes and Healy 35 ). Let (M, g) be a Riemannian manifold in ambient TIFF0007881548000104.tif3128 and TIFF0007881548000105.tif4128 be a point. If g is locally constant around p in an open neighborhood TIFF0007881548000106.tif3128 and as a result g is a constant diagonal matrix in ambient coordinates, then the volume TIFF0007881548000107.tif9128 of a sphere centered at p with respect to g TIFF0007881548000108.tif3128 is In TIFF0007881548000109.tif4128, the geodetic distance from p to any point is: The formula is TIFF0007881548000110.tif6128, where γ is the ambient space. This is the radius of the sphere in TIFF0007881548000111.tif3128. TIFF0007881548000112.tif4128 is a metric for ambient space.

[0160] Assuming that the inventors can compare data points across the entire bounded manifold with a suitable metric, they can use Lemma 2 in its non-commutating union to calculate geodesics between points on the projection of each bounded manifold under a user-provided ambient metric. i and M j Given, the inventors of this invention, TIFF0007881548000113.tif4128 and The fact that it is TIFF0007881548000114.tif4128 Calculate the pairwise geodetic distance between TIFF0007881548000115.tif4128, The components necessary to construct the metric for TIFF0007881548000116.tif4128 can be obtained. By extension, non-compounding To TIFF0007881548000117.tif4128, Bounded manifold M is TIFF0007881548000118.tif4128 i and M jThe concatenation of calculations of all pairwise distances between and provides components for constructing geodesics on a complete cobordism over all pairwise combinations of bounded manifolds. However, as a result of using Lemma 2 to approximate the projection of manifold geodesics, the inventors have a directed and incompatible view of geodesics to and from bounded manifolds across the cobordism. The inventors can interpret these directed geodesics as being defined on a directed cobordism. The inventors aim to encode directed geodesics and bounded manifold geodesics in a directed cobordism in a single data representation.

[0161] The inventors' subsequent goal is to construct an undirected cobordism, where the directed geodesics of a directed cobordism are decomposed into a single symmetric matrix representation. To this end, the inventors can transform each of the extended pseudometric spaces described above into a fuzzy simplicial set using a fuzzy single-set functor (see Definition 9), which incorporates both the topological representation of the directed cobordism and the underlying metric information. Inconsistencies in the directed cobordist geodesics can be resolved using a norm of the inventors' choice. A natural choice for the fuzzy set representation of the inventors' choice is the t-norm (also known as the fuzzy intersection). When the inventors probabilistically interpret the fuzzy simplicial set representation of a directed cobordism, its intersection corresponds to a joint distribution of the directed cobordist metric space, which highly emphasizes directed cobordist geodesics occurring in both directions.

[0162] The final step is to integrate the bounded manifold geodesics with the symmetric cobordist geodesics obtained from the fuzzy product. The inventors can do this by employing a fuzzy union (stochastic t-conorm) across the extended pseudometric space genus, as in the original UMAP implementation. As a result, in addition to the individual bounded manifolds containing their own geometry, we obtain a cobordism containing its own geometric structure incorporated into the cobordist geodesic.

[0163] The optimization of low-dimensional representations of cobordisms can be achieved in numerous ways; however, for consistency, we choose to optimize the embedding using fuzzy set cross-entropy (Definition 1), as in the original UMAP implementation. Note that since our algorithm generates symmetric matrices, PatchMAP can also be iteratively applied to construct hierarchically dimensional "nested" cobordisms.

[0164] PatchMAP implementation. To construct a cobordism, PatchMAP first computes bounded manifolds by constructing fuzzy simplicial sets from each provided dataset, i.e., system states, 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 details on how this is done, readers should refer to equations 5 and 6. The resulting directed NN graphs between UMAP submanifolds are weighted, and these reflect misfitted Riemannian metrics; that is, they cannot be simply added or multiplied to unify their weights. Therefore, the inventors stitch the cobordism metrics and make the directed NN queries fitable by applying the fuzzy simplicial intersection, resulting in a weighted symmetric graph (FuzzyIntersection, Algorithm 2). The final cobordism generated by PatchMAP is obtained by adopting a fuzzy union across all genus sets of fuzzy simplices (FuzzyUnion, Algorithm 2). To represent the connections between bounded manifolds in the PatchMAP cobordism projection, we implemented a hammer edge bundling algorithm in the Datashader Python library. Pseudocode outlining the PatchMAP algorithm is shown below.

[0165] Algorithm 2: PatchMAP Input: Dataset(D1,D 2… D n ), bounded manifold ambient metric (g f ), Cobordism ambient weighing (g) Output: Kobordism (W) Execution of the stitching function TIFF0007881548000119.tif67162

[0166] Domain / information transfer (i-PatchMAP). M rq These are bounded manifolds obtained from PatchMAP, which serve as the reference and query datasets, respectively, and M r and M q Let M be the geodesic in the cobordism between the points. Specifically, rq This is a matrix where the rows represent points in the bounded manifold of references and the columns represent the nearest neighbors of the reference manifold points of the query factor manifold under a user-defined metric, and the i, j-th entries are: TIFF0007881548000120.tif4128 and Points such as TIFF0007881548000121.tif4128 Represents geodesics between TIFF0007881548000122.tif4128. The inventors predict the query dataset P by multiplying by the transferred feature matrix F. q We compute a new feature matrix, where the weight matrix W rq The transpose of is M rq of TIFF0007881548000123.tif4128 obtained through normalization: TIFF0007881548000124.tif6128. In this situation, matrix W rq This is the state p derived from the geodetic distance on the cobordism. i and p j This can be interpreted as a single-step transition matrix of a Markov chain between [points].

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

[0168] Generation of imaging mass cytometry data. Frozen tissue was sequentially sectioned into 10 μm thick sections using a Microm HM550 cryostat (Thermo Scientific) and thawed and mounted on SuperFrost® Plus Gold charged microscopy slides (Fisher Scientific). After temperature equilibration to room temperature, the tissue sections were fixed in 4% paraformaldehyde (Ted Pella) for 10 minutes, and then rinsed three times with cytometry-grade phosphate buffer (PBS) (Fluidigm). Nonspecific binding sites were blocked at room temperature for 1 hour with 5% bovine serum albumin (BSA) (Sigma Aldrich) in PBS containing 0.3% Triton X-100 (Thermo Scientific). A metal conjugate primary antibody (Fluidigm) at an appropriately titrated concentration was mixed with 0.5% BSA in DPBS and applied overnight at 4°C in a humidified chamber. Subsequently, sections were washed twice with PBS containing 0.1% Triton X-100 and counterstained with iridium (Ir) intercalator (Fluidigm) at a ratio of 1:400 in PBS at room temperature for 30 minutes. Slides were rinsed in cytometry-grade water (Fluidigm) for 5 minutes and air-dried. Data acquisition was performed using the Hyperion imaging system (Fluidigm) and CyTOF software (Fluidigm) with 33 channels, a frequency of 200 pixels / second, and a spatial resolution of 1 μm. After visualizing the images with MCD Viewer software (Fluidigm), the data was exported as a text file for further analysis. After imaging, slides were quickly stained with 0.1% toluidine blue (Electron Microscopy Sciences) to reveal macroscopic morphology. Slides were digitized using a digital camera at a resolution of approximately 2.75 μm / pixel.

[0169] Generation of mass spectrometry imaging data. A pair of 10 μm thick sections from the same tissue block used for imaging mass cytometry were fused and mounted on indium tin oxide (ITO) coated glass slides (Bruker Daltonics). The tissue sections were coated with 2,5-dihydroxybenzoic acid (40 mg / mL in 50:50 acetonitrile:water containing 0.1% TFA) using 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: cation 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. The region of interest for data acquisition was 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 tissue 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 bright-field scanner (Leica Biosystems).

[0170] Preprocessing of mass spectrometry imaging data. Data was processed at SCiLS LAB 2018 using total ion number normalization to the mean spectrum and peak centroiding with an interval width of ±25 mDa. For all analyses, the peak range of m / z 400–1,000 was used after peak centroiding, yielding 9,753 m / z peaks. Peak picking was not performed on the presented data unless explicitly mentioned. Data was exported from SCiLS Lab as imzML files for further analysis and processing.

[0171] Single-cell segmentation. To quantify single-cell parameters in IMC and registered MSI data within the DFU dataset, we performed cell segmentation on 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 excised from the IMC data and exported in HDF5 format for supervised learning. To ensure that each excised region was a representative training sample, global thresholds were created for each region using Otsu thresholding for iridium (nucleus) staining with the Scikit-image Python library. The excised regions needed to contain more than 30% of the pixels for each of their respective thresholds.

[0172] The training region was annotated for "background," "membrane," "nucleus," and "noise." Gaussian smoothing features, edge features (including Gaussian Laplacian features, Gaussian gradient intensity features, and Gaussian feature differences), and texture features (including structural tensor eigenvalues ​​and Gaussian Hessian eigenvalues) were incorporated into the random forest classification. Using the trained classifier, the assignment probability of each pixel to one of the four classes in a complete image was predicted, and the predictions were exported as a 16-bit TIFF stack. To remove artifacts in cell staining, the noise prediction channel was Gaussian blurred by sigma 2, and Otsu thresholding with a correction factor of 1.3 was applied to create a binary mask that separated the foreground (high pixel probabilities are noise) from the background (low pixel probabilities are noise). Using the noise mask, zero values ​​from the other three probability channels (nucleus, membrane, background) from Ilastik were assigned to all pixels considered to be in the foreground in the noise channel. Probabilistic images of the three channels—nucleus, membrane, and background—with noise reduction were used for single-cell segmentation in CellProfiler (version 3.1.8)

[59] .

[0173] Single-cell parameter quantification. Single-cell parameter quantification was performed on IMC and MSI data using the in-house modified quantification (MCQuant) module in multi-limb selection microscopy software (MCMICRO)

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

[0174] Imaging mass cytometry cluster analysis. Cluster analysis was performed in Python using the Leiden community detection algorithm along with the leidenalg Python package. A set of simplices (weighted undirected graphs) of UMAPs constructed with 15 nearest neighbors and Euclidean metrics was used as input for community detection.

[0175] Micro-environment correlation network analysis. To calculate the relationships between MSI and IMC modalities, the inventors used Spearman's correlation coefficient with the Python Scipy library. M / z peaks from MSI data that did not correlate with IMC data and had Bonferroni-corrected P-values ​​greater than 0.001 were excluded from the analysis. Correlation modules were formed using hierarchical Louvain community detection with the Scikit network package. Resolution parameters used for community detection were selected based on the elbow point of a graph plotting resolution versus modularity of the community detection results. A set of simplices of nearest neighbors and Euclidean metrics, constructed with the nearest neighbor and Euclidean metrics, was used as input for community detection after the inverse cosine transform of Spearman's correlation coefficient to form metric distances. Visualization of the MSI correlation module trend to IMC parameters was calculated using exponentially weighted moving averages in the Python Pandas library after standard scaling of IMC and MSI single-cell data. For plotting, the MSI moving averages were further minimax-scaled to the range of 0-1. The differential correlation between the variable value u from MSI data and v from IMC data between conditions a and b was quantified using the formula: Ranking is performed using TIFF0007881548000125.tif6128, where the change in the correlation coefficient for each pair u,v between the conditions is weighted by the maximum absolute correlation coefficient between the two conditions. Differential correlation The significance of TIFF0007881548000126.tif4128 was calculated using a one-sided Bonferroni-corrected z-statistic after Fisher transformation.

[0176] Benchmarking of dimensionality reduction algorithms. The method used to benchmark dimensionality reduction algorithms is outlined in the HDIprep dimensionality reduction validation in Supplementary Note 3.

[0177] Benchmarking of spatial subsampling. The default subsampling parameters in MIAAIM are based on experiments across IMC data from DFU, tonsil, and prostate cancer tissue, recording the Procrustes transform sum of squared errors between subsampled UMAP embeddings, subsequent projection of out-of-sample pixels, and full UMAP embedding using all pixels. Benchmarking of spatial subsampling was performed across a wide range of subsampling rates.

[0178] Benchmarking of spectral landmarks. Subsampling rates and dimensions for validating the steady-state UMAP dimension based on landmarks were determined on a case-by-case basis from empirical studies and aligned with those used in the presented registration data. Parameters were selected based on the computational load of calculating cross-entropy on large datasets. To compare the landmark steady-state dimension selection with the subsampled data, the inventors compared the shape of exponential regression fits from both datasets using the sum of squared errors. The sum of squared errors were calculated over a wide range of landmarks obtained.

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

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

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

[0182] Supplementary note Supplementary Note 1. Combination of MIAAIM with existing bioimaging analysis software The core functionality of MIAAIM enables comparisons between technologies and tissues. As demonstrated in our proof-of-principle examples, its functionality has a wide range of applications that can be configured and executed using other software applications. We anticipate that a challenge for many users will be performing continuous registration and analysis across various software. MIAAIM's multiple output data formats interface directly with numerous tools for visualization, cell segmentation, and single-cell analysis (Table 2), paving the way for continuous investigation of multimodal tissue portraits under various circumstances.

[0183] [Table 2]

[0184] Supplementary Note 2. Notes on the expected performance of the HDIreg workflow. A fundamental premise in intensity-based image registration is the existence of quantifiable relationships between modalities—which are often met in practice, as demonstrated in the application of our proof-of-principle. However, this premise can be compromised by artifacts such as folding, tearing, and, in the case of serial sectioning, nonlinear deformation. In our experience, glandular tissue, such as that derived from the prostate, is likely to exhibit high structural 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 resulting from serial tissue sectioning. By using the Elastix library, HDIreg also provides numerous similarity measures for single-channel registration, in addition to the manifold alignment scheme used for multi-channel registration. We note that in these single-channel registration situations, histogram-based mutual information is superior to KNN α-MI. 16 This was used in the inventors' benchmark study.

[0185] Supplementary Note 3. Verification of HDIprep Dimensionality Reduction Benchmarking of dimensionality reduction algorithms (Figures 18A-18J, 19A-19H, and 20A-20H) The inventors' investigation encompasses a wide range of dimensionality reduction methods, from local nonlinear methods to global linear methods. Methods examined include t-distribution stochastic neighborhood embedding (t-SNE), homogeneous manifold approximation and projection (UMAP), PHATE (potential of heat diffusion for affinity-based transition embedding), isometric mapping (Isomap), non-negative matrix factorization (NMF), and principal component analysis (PCA).

[0186] To evaluate the ability of each method to provide appropriate data representation while enabling multimodal support, the inventors measured (i) the ability to accurately represent data modalities by generalizing any number of features or required degrees of freedom, (ii) the ability to concisely incorporate data complexity, (iii) the ability to maximize information content shared between imaging modalities, (iv) robustness to noise, and (iv) efficient computational capability.

[0187] (i.~ii.) Estimation of intrinsic data dimension. To identify a suitable method for reducing the complexity of image datasets based on mass spectrometry, the inventors hypothesized that introducing more degrees of freedom (i.e., increased embedding dimension) to the coordinates of the embedded data would result in increased similarity between the embedding of each method and its higher-dimensional counterpart with respect to the objective function of each algorithm. Therefore, the inventors identified a suitable target dimension for embedding the data for each method by examining each algorithm separately with its distinct objective function and analyzing the objective function error generated by each method after embedding the data in the increased dimension. To do this, the inventors used Euclidean n-space. Appropriate scores were generated for each dimensionality reduction method across organizational type and ascending embedding dimension to estimate errors associated with embedding MSI data into TIFF0007881548000131.tif3128. For this analysis, the inventors focused on MSI data rather than IMC data, from which they found that application to most dimensionality reduction methods was not feasible due to data size (number of pixels / resolution).

[0188] To determine the intrinsic dimension of the dataset estimated by each method, the inventors identified points in the error graph of each method where increasing the dimension does not reduce the embedding error. To do this, the inventors used an exponential increase in the potential position of the points (i.e., a copy of the real number line) to represent the increase in dimension. By modeling TIFF0007881548000132.tif3128 as increasing, we observed a natural increase in the dimensionality of real-valued data. Therefore, we constructed 95% confidence intervals (CIs) by fitting least-squares exponential regression to the error curve of the data embedding and modeling the Gaussian residual process. We selected the optimal embedding dimension for each method by simulating samples along the expected value of the fitted curve and identifying the first example of integer values ​​that fall within the 95% CI of the exponential asymptote. In this way, the minimum degrees of freedom required to capture data complexity were identified. The mean error curves for each method across 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.

[0189] UMAP. The UMAP algorithm falls into 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 space. In practice, the lower-dimensional fuzzy simplicial set is optimized so as to minimize the fuzzy set cross-entropy between its higher-dimensional counterparts. Fuzzy set cross-entropy is explicitly defined in the way given by McInnes and Healy

[15] .

[0190] While the theoretical basis of UMAP is based on category theory, the practical implementation of UMAP reduces to a weighted graph. To provide an estimate of the essential dimension of the data determined by UMAP, the inventors used an open-source implementation in Python with 15 nearest neighbors and a minimum distance of 0.1 in the resulting embeddings, and the inventors enabled the algorithm to optimize the embeddings with default values ​​of 200 iterations per dimension. Using the Python translation module of the MATLAB UMAP implementation, the dimensional cross-entropy between a set of high-dimensional fuzzy simplices and their lower-dimensional counterparts was calculated.

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

[0192] Isomap. Isomap is a manifold-based dimensionality reduction method that uses classical multidimensional scaling (MDS) to preserve point-to-point geodetic distances. To do this, point-to-point geodetic distances are determined by the shortest path graph distance using a Euclidean metric. The pairwise distance matrix represented by this graph is then embedded in n-dimensional Euclidean space via classical MDS, a metric-preservation technique that finds the optimal transformation for preserving point-to-point Euclidean metrics. As a result of the implicit linearity in classical MDS, the inventors have found that R is The standard linear correlation coefficient between the geodesic distance matrix and the pairwise Euclidean distance matrix in TIFF0007881548000133.tif3128 is 1-R. 2 The intrinsic dimension of the data is estimated by calculating the reconstruction error in each dimension using the following method. For all calculations, 15 nearest neighbors are selected to determine the shortest path graph distance, and the Minkowski metric of the square of the norm of the difference is used. TIFF0007881548000134.tif4128 was selected. All Isomap calculations were performed using Scikit-based learning.

[0193] 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 calculating the latent distances between data points by comparing each pair of diffusion distributions of points in the dataset with all others. These latent distances are then embedded in n-dimensional space using classical MDS, followed by metric MDS. Metric MDS is suitable for embedding points that do not have similarity given by any metric, and the Euclidean constraint imposed by classical MDS is overcome by the following stress function S: Relaxation is achieved by minimizing TIFF0007881548000135.tif18128, where D is the point χ of the original dataset. 1… χ n It is a metric defined over a certain period, TIFF0007881548000136.tif4128 is the corresponding embedded data point in dimension n. This stress function extends to a least-squares optimization problem. In the scalable form of PHATE used for large datasets, landmarks are embedded in n-dimensional Euclidean space instead of points, based on their pairwise latent distances, using the above stress function. Out-of-sample embedding of all data points is performed by calculating a linear combination of t-step transition matrices from point to landmark, using the embedded landmark coordinates as weights. If the stress function for metric MDS is zero, the dimensionality reduction process can fully embed and incorporate the distances between data points. This is thought to provide an error estimate used for analyzing the intrinsic data dimension for the complete dataset and the complete PHATE algorithm; however, in landmark-based calculations, not all points are embedded using metric MDS. Given a linear interpolation scheme and an initialization of scalable PHATE using classical MDS for landmark latent distances, the inventors can perform R The linear correlation coefficient between the point-landmark transition matrix and the pairwise Euclidean distance matrix in TIFF0007881548000137.tif3128 is 1-R. 2 We assumed that the reconstruction error given by provides an estimate of the error associated with embedding the complete dataset. All PHATE calculations were performed in Python using the default number of 15 nearest neighbors and 2,000 landmark points.

[0194] NMF. Non-negative matrix factorization (NMF) is a linear dimensionality reduction technique aimed at minimizing the divergence between an input matrix X and its reconstruction WH obtained through matrix factorization. Through this factorization, linear combinations of columns of W are generated using weights from H. Using the Frobenius norm between X and WH in our calculation, the divergence between the two is: The result is calculated as TIFF0007881548000138.tif6128. Therefore, this divergence or reconstruction error was plotted to estimate the error associated with each embedding dimension. For all calculations, each channel in the dataset was minimax rescaled to the range of 0 to 1 so that only positive elements were included in X. All calculations were performed using Scikit-Learning.

[0195] PCA. Principal Component Analysis (PCA) is a linear dimensionality reduction method that aims to capture the major axes of variability in data at a global level. To determine the intrinsic dimension of the dataset estimated by PCA, the cumulative rate of residual variance remaining after dimensionality reduction is plotted for each component. Given a component 1 ≤ d ≤ n-1 where n is the number of dimensions 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 complete dataset. For all calculations, each channel in the dataset was standardized by removing the mean and scaling to the unit variance. Standardization is used to ensure that there are no features that dominate the PCA objective function. All calculations were performed using Scikit-Learning.

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

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

[0198] For affine registration, the number of selected resolutions was obtained in a multi-resolution pyramidal hierarchy by hyperparameter search. For nonlinear registration, both the number of resolutions and the final uniform grid spacing for B-spline control points were determined by hyperparameter grid search. In both registrations, the number of resolutions either improved the registration results or did not change the registration at all. However, during nonlinear registration, finer control point grid spacing schedules resulted in improved registration as indicated by mutual information, but these also led to regions with unrealistic strains even with the addition of regularization using deformation bending energy penalties. A value of 300 was chosen for the final grid spacing to strike a balance between the improvement in registration and the increase in strain as indicated by the cost function.

[0199] Subsequently, the resulting deformation fields were applied to grayscale hyperspectral images prepared from each dimensionality reduction algorithm, and they were spatially aligned evenly with the H&E images of each tissue. Before calculating the mutual information between the H&E and the embedded MSI images, non-zero common areas were applied to pairs of images. Non-zero common areas are used to account for any edge effects introduced into the registration by using three manually selected MSI peaks, which are considered to have adverse effects on registration and mutual information calculation if they are not adequately represented at all locations in the image in our analysis. Subsequently, the mutual information between each registered dimensionality reduction image (n=5 per method) was calculated using the SimpleITK Partzen window-based method (Figures 18B, 19B, and 20B).

[0200] (iv.) Evaluation of the algorithm's robustness to noise. Through evaluation of the intrinsic dimension of the data, the inventors learned that both high-dimensional imaging modalities (MSI and IMC) follow a manifold structure, in which case the dimension of the data can be approximated with fewer degrees of freedom than the number of parameters initially given to the ambient space. Using this information, and as evidence of the validity of the assumption of such a manifold structure, in addition to the visual quality of the spatial map that each method returns to the subsequent organization, the inventors then decided to compare the ability of each algorithm to preserve geodesic distances in the low-dimensional embedding with and without "noisy" peaks and / or technical variations.

[0201] To do this, the inventors utilized a denoised manifold preservation (DEMaP) metric. By calculating the DEMaP metric (Spearman's rank correlation coefficient) between the ambient space geodesic distance of the peak-picked MSI dataset and the pairwise embedding Euclidean distance between data points from the corresponding non-peak-picked dataset, the inventors 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 along with 15 nearest neighbors, or, since they inherently assume a Euclidean structure, the inventors calculated the geodesic distance in the peak-picked MSI dataset using 15 nearest neighbors with the Euclidean metric. Peak picking was performed at SCiLS Lab 2018b using orthogonal matching tracking with a maximum of 1,000 peaks. The DEMaP scores for each algorithm across five random initializations for each MSI dataset are shown in Figures 18I, 19G, and 20G.

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

[0203] References TIFF0007881548000139.tif58164TIFF0007881548000140.tif241164TIFF0007881548000141.tif235164 TIFF0007881548000142.tif253165TIFF0007881548000143.tif247165TIFF0007881548000144.tif229165

[0204] Other embodiments 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 relation to specific embodiments, it should be understood that the claimed invention should not be unduly limited to such specific embodiments. In fact, various modifications of the described embodiments for carrying out the invention, which will be apparent to those skilled in the art, are intended to fall within the scope of the invention.

[0205] Other embodiments are described in the claims.

Claims

1. A method for identifying cross-modal features from two or more spatially decomposed datasets, (a) registering the two or more spatially decomposed datasets to generate an aligned feature image obtained from multi-omics image alignment, which includes the two or more spatially aligned spatially decomposed datasets; and (b) Step of extracting the cross-modal features from the aligned feature images. Includes, The method wherein at least two of the two or more spatially decomposed datasets are datasets from different modalities.

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

3. The method according to claim 2, wherein the dimensionality reduction is performed by homogeneous manifold approximation and projection (UMAP), isometric mapping (Isomap), t-distribution stochastic neighborhood embedding (t-SNE), PHATE (potential of heat diffusion for affinity-based transition embedding), principal component analysis (PCA), diffusion mapping, or non-negative matrix factorization (NMF).

4. The method according to claim 3, wherein the dimensionality reduction is performed by homogeneous manifold approximation and projection (UMAP).

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

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

7. The method according to any one of claims 1 to 6, further comprising the step of clustering the two or more spatially decomposed datasets in order to complement them with a similarity matrix representing the similarity between data points.

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

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

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

11. The method according to any one of claims 7 to 10, comprising the step of modeling cluster-cluster spatial interactions.

12. A method according to any one of claims 7 to 10, including analysis based on intensity.

13. The method according to any one of claims 7 to 10, comprising analysis of the abundance of cell types or heterogeneity of a predetermined region in the dataset.

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

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

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

17. A method according to any one of claims 7 to 10, comprising analysis of predictions of spatial niches.

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

19. The method according to claim 18, wherein the classification step is performed by a hard classifier, a soft classifier, or a fuzzy classifier.

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

21. The method according to claim 20, further comprising the step of analyzing spatially resolved objects.

22. The method according to claim 21, wherein the step of analyzing spatially resolved objects includes segmentation.

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

24. The method according to any one of claims 1 to 23, wherein step (b) includes a sorting test relating to the enrichment or depletion of cross-modal features.

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

26. The method according to claim 24 or 25, wherein the sorting test is performed by mean sorting test.

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

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

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

30. The method according to any one of claims 1 to 29, wherein at least one of the two or more spatially resolved datasets is an image from co-detection by immunohistochemical examination, imaging mass cytometry, multiplex ion beam imaging, mass spectrometry imaging, cell staining, RNA-ISH, spatial transcriptome analysis, or index imaging.

31. The method according to claim 30, wherein at least one of the spatially resolved datasets is an immunofluorescence image.

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

33. The method according to any one of claims 30 to 32, wherein at least one of the spatially resolved datasets is a multiplex ion beam image.

34. At least one of the spatially resolved datasets is a mass spectrometry image, The method according to any one of claims 30 to 33, wherein the mass spectrometry imaging is MALDI imaging, DESI imaging, or SIMS imaging.

35. At least one of the spatially resolved datasets is cell staining, The method according to any one of claims 30 to 34, wherein the cell staining is H&E staining, toluidine blue staining, or fluorescent staining.

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

37. The method according to any one of claims 30 to 36, wherein at least one of the spatially decomposed datasets is a spatial transcriptome analysis.

38. The method according to any one of claims 30 to 37, wherein at least one of the spatially decomposed datasets is co-detected by index imaging.

39. A computer-readable storage medium, A computer program for identifying cross-modal features from two or more spatially decomposed datasets is stored in the computer-readable storage medium. The computer program includes a routine set of instructions for causing the computer to perform the steps of the method described in any one of claims 1 to 38. The aforementioned computer-readable storage medium.

Images