A single-cell multi-omics clustering method based on multi-subspace contrastive learning

Abstract

The application discloses a single-cell multi-omics clustering method based on multi-subspace contrast learning, comprising the following steps: extracting important individual features of transcriptome and chromatin accessibility of different omics based on an automatic encoder; extracting common features of the transcriptome and the chromatin accessibility by using multi-subspace contrast learning, so as to effectively reveal deep correlations between omics data; reconstructing data by using the automatic encoder for specific omics chromatin accessibility and transcriptome, so as to better denoise and reconstruct features conforming to omics information; and fusing individual information of specific omics and common information of multiple omics to obtain a common embedding more conducive to clustering. The application can make up for the deficiency of single-space contrast, and construct a more comprehensive and larger information representation for fusion and clustering of single-cell multi-omics data.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of single-cell multi-omics data analysis, in particular to a single-cell multi-omics clustering method based on multi-subspace contrastive learning. BACKGROUND

[0002] Single-cell multi-omics joint sequencing technologies are developing rapidly, such as sci-CAR, SNARE-seq, Paired-seq and 10XGenomics Multiome. These sequencing methods simultaneously sequence gene expression and chromatin accessibility in the same cell, providing valuable data resources for studying single-cell multi-omics clustering methods. However, the integration of single-cell multi-omics data remains a great challenge. This is mainly due to the inherent high sparsity of the data, the huge heterogeneity caused by measurement noise, the huge dimension difference between single-cell chromatin accessibility sequencing scATAC-seq and single-cell RNA sequencing scRNA-seq data, about 10-20 times, and the imbalance of information between omics, which together bring great computational and analytical challenges.

[0003] The integration of single-cell multi-omics data can promote the study of complex biological information. Researchers have been developing single-cell multi-omics integration methods using machine learning and bioinformatics techniques. Some methods are based on non-negative matrix factorization or principal component analysis to integrate single-cell multi-omics data and address cell heterogeneity. At the same time, with the development of deep learning and its ability to extract expression features, deep learning methods have gradually become the mainstream technology for single-cell data analysis. Furthermore, graph-based methods use graphs to represent the relationship between different cells, where cells are treated as nodes, and then the integration is completed by manipulating this graph. This method emphasizes neighborhood structure and sometimes aggregates information between neighborhoods, resulting in additional robustness to noise.

[0004] Although the above methods are effective, they still have some problems. First, most of them focus on features common to different omics, but ignore features unique to each omic, making it difficult to better integrate biological features at different levels to learn more distinct cell clusters. Second, most of them ignore the high-dimensional sparsity and noise inherent in multi-omics data, which limits feature extraction. SUMMARY

[0005] The single-cell multi-omics clustering method based on multi-subspace contrastive learning provided by the present application constructs a more comprehensive and information-rich representation for the fusion and clustering of single-cell multi-omics data, and can at least solve one of the above technical problems.

[0006] To solve the above technical problems, the present application adopts the following technical solutions:

[0007] A single-cell multi-omics clustering method based on multi-subspace contrastive learning, comprising the following steps:

[0008] S1, based on an autoencoder, important individual features are extracted for transcriptome and chromatin accessibility of different omics;

[0009] S2, multi-subspace contrastive learning is used to extract common features of transcriptome and chromatin accessibility to fully mine the common features of transcriptome and chromatin accessibility and effectively reveal the deep relationship between omics data;

[0010] S3: Use the autoencoder for specific omics chromatin accessibility and transcriptome to reconstruct the data to better denoise and reconstruct the features conforming to the omics information;

[0011] S4: The individual information of the specific omics of S3 and the common information of multiple omics of S2 are fused to obtain a common embedding more conducive to clustering.

[0012] Further, in S1, the single-cell multi-omics data is respectively mapped to the respective nonlinear embedding space by using the autoencoder, the individual features are retained, and the robustness to noise and outliers is enhanced;

[0013] Specifically, let , be the standardized scRNA-seq data and scATAC-seq data, respectively, where N is the number of samples, Dr and Da are the number of features, and two independent encoders and are used to learn the corresponding d-dimensional feature representation , } , where d is the dimension of the embedding space, and the expressions are respectively:

[0014] ;

[0015] In these two expressions, is the latent low-dimensional feature representation of cells and genes in scRNA-seq data, is the latent low-dimensional feature representation between cells and peaks in scATAC-seq data.

[0016] Further, S2 further comprises:

[0017] S2.1, define two views u and r, from the feature perspective, input the view u and the view r into the multiple subspace encoders MSEncoder to generate the node embedding of multiple subspaces, and the process is represented as:

[0018]

[0019]

[0020] where l denotes the l-th subspace, p denotes the number of subspace, , and are the node embeddings of the i-th subspace in view u and view r respectively, d denotes the dimension of the embedding space, and U and R are the feature representations obtained by multiple subspace encoders, and are the common features of the same omics in view u and view r respectively, are the latent low-dimensional feature representations of cells and genes in scRNA-seq data, are the latent low-dimensional feature representations between cells and peaks in scATAC-seq data;

[0021] The subspace encoder MSEncoder maps the feature matrix X into the key embedding K and the value embedding V using two independent linear projections for the input feature matrix X, and the corresponding expressions are respectively:

[0022]

[0023]

[0024] where, , and denote the different projection matrices of keys and values respectively;

[0025] S2.2, in order to make the features captured by multi-subspace contrastive learning more diverse and focus on more important node information, make full use of the degree information obtained from the graph structure, integrate the graph structure information into the attention mechanism to enhance the discriminability of the feature representation;

[0026] The initialization expression of the degree embedding matrix S is:

[0027] ;

[0028] where, denotes the projection matrix corresponding to the degree, and D is the degree matrix;

[0029] The expression of the attention weight matrix is:

[0030] ;

[0031] where, Att is the output of the entire expression, denoting the "attention weight matrix", the Softmax(·) function is a function that "compresses" an arbitrary real-valued vector into a probability distribution vector, Q is the query matrix, and a is a learnable or settable hyperparameter;

[0032] To compute the embedding representation of nodes within each subspace, the value embedding V is broadcast to each subspace as , obtaining the node embeddings of p subspace, the expression is:

[0033] ;

[0034] ;

[0035] S2.3, by calculating and The mutual information between the node embeddings across subspace is maximized and the mutual information within the subspace is minimized, respectively as the respective subspace contrast loss function;

[0036] S2.4, the final objective function of subspace contrast learning is the weighted sum of mutual information maximization target and mutual information minimization target, the formula is:

[0037]

[0038] In this formula, L contra represents the subspace contrast loss, θ' and φ, are the parameters of the encoder in view u and view r branch respectively, γ is the weight coefficient balancing the influence of intra-space and inter-space target, MI(˙) is the mutual information, k represents the kth subspace, and respectively represent the node embedding of the kth and lth subspace in view u, represents the node embedding of the kth subspace in view r;

[0039] S2.5, the view U and view R are trained and spliced, and finally the latent feature representation .

[0040] Further, in the S2.3, the feature representation of view u and view r is output to p subspace, that is In order to achieve the goal of optimizing the subspace encoder MSEncoder, making the common features of the same omics in view r and view r, that is, the view U and the view R have diversified feature representation, a pair of loss functions is used, which focuses on maximizing the mutual information between subspace representations and minimizing the mutual information between different subspace representations. For the mutual information between subspace representations, the target is formulated as:

[0041] ;

[0042] where θ' and φ are the parameters of the encoders in view u and view r branches respectively, k denotes the k-th subspace, p denotes the number of subspace, and MI(·) is the mutual information, denotes the node embedding of the k-th subspace in view u, denotes the node embedding of the k-th subspace in view r.

[0043] Further, since and are unknown, the mutual information MI cannot be calculated, in order to estimate the mutual information MI, the sampling method is used to replace the maximum mutual information, to maximize the lower bound of the mutual information MI, and the contrastive log ratio upper bound CLUB is used to minimize the upper bound of the mutual information MI, further:

[0044] Regarding the maximization of the lower bound of the mutual information MI, the Jensen-Shannon estimation is used instead of the Kullback-Leibler estimation, and the calculation formula is:

[0045]

[0046] In the formula, is a discriminator, which takes the representations of the two views as input and scores the consistency between them, here the discriminator is simply instantiated as the dot product between the two representations, that is , SP(·) is an activation function, is the expectation of the joint distribution of the two views in the same subspace k, is the product of the expectation of the independent distribution of the two views in the same subspace k;

[0047] In addition to the objective in each space, constraints are also imposed on the relationship between different subspaces in the same view to increase the diversity of the space, for this purpose an optimization objective is proposed to maximize the mutual information, which can increase the key features in view u, and since , the expression of the optimization objective is:

[0048] ;

[0049] In the formula, l denotes the l-th subspace, denotes the node embedding of the l-th subspace in view u;

[0050] Regarding the minimization of the upper bound of the mutual information MI, the calculation formula of the contrastive log ratio upper bound CLUB is:

[0051] ;

[0052] In the formula, ICLUB Let E represent the expectation of CLUB, P(x,y) represent the joint distribution of x and y, P(x) and P(y) represent the marginal distributions of x and y respectively, log is the natural logarithm, and P(y|x) is the conditional probability of y given x. The key to calculating CLUB is the distribution... In this case, it is assumed that each dimension of x and y is correlated;

[0053] give Assuming that each dimension of y depends only on the corresponding dimension of x, then using and Substituting x and y respectively, the relevant expression is:

[0054] ;

[0055] In this formula, P represents the probability density. Represents the covariance matrix. ˙) represents a multivariate Gaussian distribution, and I represents the identity matrix. This represents the variance shared across all dimensions. and They are and The node embedding in the i-th subspace, i The multiplication symbol is for all elements or subspaces i, representing the product of a high-dimensional Gaussian distribution into a product of independent one-dimensional Gaussian distributions.

[0056] Furthermore, in S3, a decoder network based on the zero-inflated negative binomial regression model ZINB explores the global probability structure of scRNA-seq data. The ZINB distribution is derived from the mean of the negative binomial distribution. Deviation parameters and coefficients used to describe the probability of interpolated events. The definition and expression are:

[0057]

[0058]

[0059] In these two formulas, Indicates the parameter is u r The negative binomial probability mass function of x and θ. r It is a vector from the raw scRNA-seq data, representing a random variable taking non-negative integer values, u r Let θ be the mean parameter of the distribution, and θ be the deviation parameter. It is the gamma function, ZINB(˙) represents the zero-inflation negative binomial distribution, and π represents the zero-inflation probability. is x r is an indicator function when

[0060] In particular, the ZINB-based decoder estimates the parameters :

[0061] ;

[0062] ;

[0063] ;

[0064] Among them, is a matrix form of is a zero inflation probability parameter that determines the proportion of "extra zeros", is a negative binomial mean parameter that represents the expected count of non-zero parts, usually output by the decoder network, is a negative binomial dispersion parameter that controls the size of the variance, usually output by the decoder network, since the activation function of and is sigmoid(·), so the exit probability is between 0 and 1, and in addition, since the mean parameter and the dispersion parameter are non-negative, the exponential function exp is selected as the activation function of is a decoder with a fully connected layer, , and θ are three learnable parameter matrices;

[0065] Unlike traditional autoencoders based on mean square error loss, the loss function of the ZINB-based decoder network is the negative logarithm of the ZINB likelihood, expressed as:

[0066] ;

[0067] Among them, is a zero-inflated negative binomial distribution used to model sparse and over-dispersed single-cell count data, and -log(˙) represents taking the probability density and then taking the negative logarithm to get the negative log-likelihood loss, the smaller the better;

[0068] Considering the extremely sparse and almost binary nature of scATAC-seq data, a decoder network based on Bernoulli distribution Ber is selected to model scATAC-seq data, expressed as:

[0069] ;

[0070] In this formula, is a vector from the original scATAC-seq data, is the average parameter of Ber;

[0071] The decoder based on the Bernoulli distribution Ber estimates based on the input latent variable Z' through a fully connected layer with sigmoid(·) as the activation function, and the expression is:

[0072] ;

[0073] In this formula, is a in matrix form, sigmoid(·) is an activation function that compresses real values to between 0 and 1, is a weight parameter matrix, denotes the a-th decoder mapping;

[0074] Finally, the decoder based on the Bernoulli distribution Ber is optimized by the cross-entropy loss, and the expression is:

[0075] ;

[0076] In this formula, L Ber denotes the Bernoulli loss, CE(·) is the cross-entropy function, is the true observed Bernoulli variable.

[0077] Further, in the S4, by learning two latent representations and , the characteristics of individuals in different omics data are encoded, and the commonalities among the above characteristics are captured by the latent feature representation , and in this process, scale parameters and are introduced to element-wise weighted sum the two groups of latent representations and to achieve effective aggregation between them, and the formal expression is:

[0078] ;

[0079] In this formula, Z represents the co-embedding representation

[0080] Based on the generated more discriminative co-embedding representation Z, the accuracy and reliability of clustering analysis can be improved.

[0081] Further, in the S4, in order to pursue a more discriminative and informative co-embedding representation Z, the total loss The individual information and common information of the multi-omics data are combined, and the input targets of the scRNA-seq data and the scATAC-seq data are unified, and the total loss expression is:

[0082]

[0083] In the formula, represents taking the minimum value, represents the network parameters, and alpha1, alpha2 and gamma all represent weight coefficients, is a zero-inflated negative binomial distribution used to model sparse and over-dispersed single-cell count data, is a zero-inflated probability parameter that determines the proportion of "extra zeros", is a negative binomial mean parameter, representing the expected count of the non-zero part, which is usually output by the decoder network, is a negative binomial dispersion parameter that controls the size of the variance, which is usually output by the decoder network, -log(˙) represents taking the probability density and then taking the negative logarithm to get the negative log-likelihood loss, the smaller the value is, the better, and CE(·) is the cross-entropy function, is a real observed Bernoulli variable, is a matrix form of u a , and u a is the mean parameter of the Bernoulli distribution Ber, and theta' and are parameters of the encoder in the view u and the view r branch respectively, MI(·) is the mutual information, which is used to measure the degree of non-linear dependence between two random variables, k represents the kth subspace, and p represents the number of subspaces, and respectively represent the node embedding of the kth and lth subspaces in the view u, represents the node embedding of the kth subspace in the view r.

[0084] The total loss function is continuously minimized, and 500 epochs are trained for convergence optimization, to learn individual feature representation and shared feature representation from multi-omics data, and to combine the individual feature representation and the shared feature representation into the clustering information co-embedding representation Z.

[0085] The beneficial effects of the present application are embodied in:

[0086] The present application can capture the complementary information between different feature space omics data and the internal relationship between the omics data in the same feature space through multi-subspace contrast learning, thereby making up for the deficiency of single-space contrast, and constructing a more comprehensive and larger information quantity representation for the fusion and clustering of single-cell multi-omics data. BRIEF DESCRIPTION OF DRAWINGS

[0087] The accompanying drawings, which are included to provide a further understanding of the application and are incorporated in and constitute a part of this application, illustrate embodiments of the application and together with the description serve to explain the application. The drawings illustrate embodiments of the application and, together with the description, serve to explain the principles of the application. They must not be considered in a limiting manner in respect of the application.

[0088] Figure 1 is a schematic diagram of the overall process of the single-cell multi-omics clustering method of the embodiments of the application.

[0089] Figure 2 is a detailed process block diagram of the single-cell multi-omics clustering method of the embodiments of the application.

[0090] Figure 3 is a schematic diagram of the overall process of the single-cell multi-omics clustering method of the embodiments of the application. Figure 2 is an enlarged diagram of the working principle of the sub-space encoder MSEncoder.

[0091] Figure 4 is an enlarged diagram of the working principle of the sub-space encoder MSEncoder. Figure 2 is an enlarged diagram of the working principle of the sub-space encoder MSEncoder.

[0092] Figure 5 is a clustering UMAP plot of the scMVAE-PoE method on the dataset CellMix.

[0093] Figure 6 is a clustering UMAP plot of the scMVAE-NN method on the dataset CellMix.

[0094] Figure 7 is a clustering UMAP plot of the scMVAE-Direct method on the dataset CellMix.

[0095] Figure 8 is a clustering UMAP plot of the DCCA method on the dataset CellMix.

[0096] Figure 9 is a clustering UMAP plot of the scMCs method on the dataset CellMix.

[0097] Figure 10 is a clustering UMAP plot of the scEMC method on the dataset CellMix.

[0098] Figure 11 is a clustering UMAP plot of the scMSCL method of the application on the dataset CellMix.

[0099] Figure 12 is a structural block diagram of the computer device of the embodiments of the application. DETAILED DESCRIPTION

[0100] With reference to the drawings of the embodiments of the present application, the technical solutions in the embodiments of the present application will be clearly and completely described. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments of the present application. In the case of no conflict, the embodiments in the present application and the features in the embodiments can be combined with each other. Based on the embodiments in the present application, all other embodiments obtained by a person of ordinary skill in the art without creative work are within the protection scope of the present application.

[0101] It should be noted that the meaning of "and / or" appearing throughout the text includes three parallel solutions. Taking "A and / or B" as an example, it includes the A solution, or the B solution, or the solution of A and B being satisfied at the same time. In addition, "multiple" means more than two. In addition, the technical solutions of each embodiment can be combined with each other, but it must be based on the fact that a person of ordinary skill in the art can realize it. When the combination of technical solutions appears contradictory or unachievable, it should be considered that the combination of technical solutions does not exist and is not within the protection scope required by the present application.

[0102] Referring to Figure 1 The embodiment of the present application provides a single-cell multi-omics clustering method based on multi-subspace contrast learning, comprising the following steps:

[0103] S1, important individual features of transcriptome and chromatin accessibility of different omics are extracted based on an automatic encoder;

[0104] S2, multi-subspace contrast learning is used to extract common features of transcriptome and chromatin accessibility, so as to effectively reveal deep correlations between omics data;

[0105] S3: An automatic encoder for specific omics chromatin accessibility and transcriptome is used to reconstruct data, so as to better denoise and reconstruct features conforming to omics information;

[0106] S4: The individual information of the specific omics in S3 and the common information of the multiple omics in S2 are fused to obtain a common embedding more conducive to clustering.

[0107] Referring to Figures 2-4 In the embodiment, in S1, the automatic encoder, as a classical neural network, can map high-dimensional data to a low-dimensional representation space, effectively ignoring noise and outliers. Based on this characteristic, the single-cell multi-omics data is respectively mapped to the respective nonlinear embedding space by using the automatic encoder, so as to retain individual features and enhance the robustness to noise and outliers;

[0108] Specifically, let , are the normalized scRNA-seq data and scATAC-seq data, respectively, where N is the number of samples, Dr and Da are the number of features, and two independent encoders are adopted and to learn the corresponding d-dimensional feature representations , } where d is the dimension of the embedding space, and the expressions are

[0109] ;

[0110] In the two expressions, is the latent low-dimensional feature representation of cells and genes in the scRNA-seq data, is the latent low-dimensional feature representation between cells and peaks in the scATAC-seq data.

[0111] Referring to Figures 2-4 , attention mechanisms may cause the model to overemphasize individual features or noise of each omics, which is not conducive to data fusion and cell clustering, in addition, it is difficult to comprehensively represent the complementarity between omics from individual features alone, while information shared between omics can reflect their commonality, which is crucial for achieving high-quality consistent clustering, therefore, in the embodiment, the S2 further includes:

[0112] S2.1, define two views u and r, unlike data transformation in traditional contrastive learning, from the perspective of features, input the view u and the view r into multiple subspace encoders MSEncoder to generate node embeddings of multiple subspaces, the process is represented as:

[0113]

[0114]

[0115] In the two expressions, l represents the lth subspace, p represents the number of subspaces, , and are the node embeddings of the ith subspace in the view u and the view r, respectively, d represents the dimension of the embedding space, and the views U and R are the feature representations obtained by the multiple subspace encoders, and are the common features of the same omics in the view u and the view r, respectively, is the latent low-dimensional feature representation of cells and genes in the scRNA-seq data, is the latent low-dimensional feature representation between cells and peaks in the scATAC-seq data;

[0116] The subspace encoder MSEncoder is a multi-feature generator based on self-attention. For the input feature matrix X, the feature matrix X is mapped to the key embedding K and the value embedding V using two independent linear projections, and the corresponding expressions are respectively:

[0117]

[0118]

[0119] In the two expressions, , respectively represent the different projection matrices of the key and the value;

[0120] S2.2, in order to make the features captured by the multi-subspace contrast learning more diverse and focus on more important node information, make full use of the degree information obtained from the graph structure, integrate the graph structure information into the attention mechanism, and enhance the discriminability of the feature representation;

[0121] The degree embedding matrix S is initialized as:

[0122]

[0123] In the expression, represents the projection matrix corresponding to the degree, and D is the degree matrix;

[0124] The expression of the attention weight matrix is obtained as:

[0125]

[0126] In the expression, Att is the output of the entire expression, representing the "attention weight matrix", the Softmax(·) function is a function that "compresses" an arbitrary real-valued vector into a probability distribution vector, Q is the query matrix, and α is a learnable or settable hyperparameter;

[0127] To calculate the embedding representation of each node in each subspace, the value embedding V is broadcast to each subspace as , and the node embedding of the p subspaces is obtained, and the expression is:

[0128]

[0129]

[0130] S2.3, by calculating and the mutual information between the node embeddings across subspaces is maximized, and the mutual information within the subspace is minimized, respectively as the respective subspace contrast loss function;

[0131] ​​​​S2.4, the final objective function of subspace contrast learning is the weighted sum of the mutual information maximization objective and the mutual information minimization objective, which is formulated as:

[0132]

[0133] In the formula, L contra represents the subspace contrast loss, θ' and φ are the parameters of the encoders in the view u and view r branches respectively, γ is the weight coefficient balancing the influence of the intra-space and inter-space objectives, MI(·) is the mutual information, k represents the kth subspace, and respectively represent the node embedding of the kth and lth subspaces in the view u, represent the node embedding of the kth subspace in the view r.

[0134] S2.5, the view U and the view R are trained and spliced, and finally the latent feature representation is obtained.

[0135] Referring to Figures 2-4 , in the S2.3, the feature representations of the view u and the view r are output to p subspaces, that is, In order to achieve the goal of optimizing the subspace encoder MSEncoder to make the common features of the same omics in the view r and the view r, that is, the view U and the view R have diversified feature representations, a pair of loss functions is used, which focuses on maximizing the mutual information between the subspace representations and minimizing the mutual information between different subspace representations. For the mutual information between the subspace representations, the objective is formulated as:

[0136] ;

[0137] In the formula, θ' and φ are the parameters of the encoders in the view u and view r branches respectively, k represents the kth subspace, p represents the number of subspaces, MI(·) is the mutual information, represents the node embedding of the kth subspace in the view u, represents the node embedding of the kth subspace in the view r.

[0138] Referring to Figures 2-4 , in the embodiment, since and are unknown, the mutual information MI cannot be calculated. In order to estimate the mutual information MI, the maximum mutual information is replaced based on the sampling method to maximize the lower bound estimate of the mutual information MI, and the contrast log ratio upper bound CLUB is used to minimize the upper bound estimate of the mutual information MI, and further:

[0139] For the lower bound of MI maximization, the Jensen-Shannon estimation is used instead of the Kullback-Leibler estimation, and the formula is:

[0140]

[0141] In the formula, is a discriminator that takes as input the representations of the two views and scores the consistency between them, here the discriminator is simply instantiated as the dot product between the two representations, i.e. , SP(·) is an activation function, is the expectation of the joint distribution of the two views within the same subspace k, is the product of the expectation of the independent distribution of the two views within the same subspace k;

[0142] In addition to the objectives within each space, constraints are also imposed on the relationship between different subspace in the same view to increase the diversity of the space, for this purpose an optimization objective is proposed to maximize the mutual information, which can increase the key features within the view u, and since , the expression of the optimization objective is:

[0143] ;

[0144] In the formula, l represents the lth subspace, represents the node embedding of the lth subspace in the view u;

[0145] For the upper bound of MI minimization, the formula of the contrastive log-ratio upper bound (CLUB) is:

[0146] ;

[0147] In the formula, I CLUB represents the CLUB estimation, E represents the expectation, P(x,y) represents the joint distribution of x and y, P(x) and P(y) represent the marginal distribution of x and y respectively, log is the natural logarithm, P(y|x) is the conditional probability of y given x, and the key to calculating the CLUB is to assume that each dimension of x and y is related under the distribution ;

[0148] Given the assumption that , each dimension of y only depends on the corresponding dimension of x, and by using and to replace x and y respectively, the relevant expression is

[0149] ;

[0150] In this formula, P represents the probability density. Represents the covariance matrix. ˙) represents a multivariate Gaussian distribution, and I represents the identity matrix. This represents the variance shared across all dimensions. and They are and The node embedding in the i-th subspace, i The multiplication symbol is for all elements or subspaces i, representing the product of a high-dimensional Gaussian distribution into a product of independent one-dimensional Gaussian distributions.

[0151] See Figures 2-4 Previous studies have shown that scRNA-seq data typically exhibits characteristics such as dispersion, variance greater than the mean, and high sparsity. Nevertheless, some research reports indicate that the zero-inflated negative binomial (ZINB) distribution can effectively describe these data characteristics. Therefore, in this embodiment, in step S3, a decoder network based on the zero-inflated negative binomial regression model ZINB explores the global probability structure of scRNA-seq data. The ZINB distribution is derived from the mean of the negative binomial distribution. Deviation parameters and coefficients used to describe the probability of interpolated events. The definition and expression are:

[0152]

[0153]

[0154] In these two formulas, Indicates the parameter is u r The negative binomial probability mass function of x and θ. r It is a vector from the raw scRNA-seq data, representing a random variable taking non-negative integer values, u r Let θ be the mean parameter of the distribution, and θ be the deviation parameter. It is the gamma function, ZINB(˙) represents the zero-inflation negative binomial distribution, and π represents the zero-inflation probability. It is x r Indicator function when =0;

[0155] Specifically, the ZINB-based decoder estimates parameters based on the co-embedding representation Z through three different fully connected layers. :

[0156] ;

[0157] ;

[0158] ;

[0159] Among the three, is in the form of a matrix, is a zero inflation probability parameter determining the proportion of "extra zeros", is a negative binomial mean parameter representing the expected count of non-zero parts, usually output by the decoder network, is a negative binomial dispersion parameter controlling the size of the variance, usually output by the decoder network, since The activation function of is sigmoid(·), so the exit probability is between 0 and 1, and in addition, since the mean parameter and the dispersion parameter are non-negative, the exponential function exp is selected as the activation function of and , is a decoder with a fully connected layer, , and θ are three learnable parameter matrices;

[0160] Unlike traditional autoencoders based on mean square error loss, the loss function of the ZINB-based decoder network is the negative logarithm of the ZINB likelihood, expressed as:

[0161] ;

[0162] Among the three, is a zero-inflated negative binomial distribution used to model sparse and over-dispersed single-cell count data, and -log(˙) represents taking the probability density and then taking the negative logarithm to get the negative log-likelihood loss, which is better the smaller the value is;

[0163] Considering the extremely sparse and almost binary nature of scATAC-seq data, a decoder network based on Bernoulli distribution Ber is selected to model scATAC-seq data, expressed as:

[0164] ;

[0165] Among the three, is a vector from the original scATAC-seq data, is the average parameter of Ber;

[0166] The decoder based on Bernoulli distribution Ber estimates based on the input latent variable Z' through a fully connected layer with sigmoid(·) as the activation function, expressed as:

[0167] ​ ;

[0168] In this formula, is a a matrix form, sigmoid(·) is an activation function that compresses real values to between 0 and 1, is a weight parameter matrix, denotes the a-th decoder mapping;

[0169] Finally, the decoder based on the Bernoulli distribution Ber is optimized by the cross-entropy loss, and the expression is:

[0170] ;

[0171] In this formula, L Ber denotes the Bernoulli loss, CE(·) is the cross-entropy function, is the true observed Bernoulli variable.

[0172] Referring to Figures 2-4 , in this embodiment, in the S4, by learning two latent representations and , the characteristics of individuals in different omics data are encoded respectively, and the commonalities among the above characteristics are also captured by the latent feature representation , and in this process, scale parameters and are introduced to element-wise weighted sum the two groups of latent representations and to achieve effective aggregation between them, and the formulaic expression is:

[0173] ;

[0174] In this formula, Z denotes the co-embedding representation

[0175] Based on the generated more discriminative co-embedding representation Z, the accuracy and reliability of clustering analysis can be improved.

[0176] Referring to Figures 2-4 , in this embodiment, in the S4, in order to pursue a more discriminative and informative co-embedding representation Z, the total loss combines the individual information and common information of multi-omics data, and unifies the input targets of scRNA-seq data and scATAC-seq data, and the total loss expression is:

[0177]

[0178] In this formula, denotes taking the minimum value, denote network parameters, and a1, a2 and g are weight coefficients, is a zero-inflated negative binomial distribution used to model sparse and over-dispersed single-cell count data, is a zero-inflation probability parameter that determines the proportion of "extra zeros", is a negative binomial mean parameter representing the expected count of non-zero parts, usually output by the decoder network, is a negative binomial dispersion parameter that controls the size of the variance, usually output by the decoder network, -log(˙) represents taking the probability density and then taking the negative logarithm to get the negative log-likelihood loss, the smaller the value is, the better, and CE(·) is the cross-entropy function, is a real observed Bernoulli variable, is a matrix form of u a a is the mean parameter of the Bernoulli distribution Ber, and and are the parameters of the encoder in the view u and view r branch respectively, MI(·) is the mutual information, which is used to measure the degree of non-linear dependence between two random variables, k represents the kth subspace, and p represents the number of subspaces, and represent the node embedding of the kth and lth subspaces in view u respectively,

[0179] The total loss function is continuously minimized, and 500 epochs are trained for convergence optimization, to learn individual feature representations and shared feature representations from multiple omics data, and to combine the individual feature representations and shared feature representations into the clustering information co-embedding representation Z.

[0180] In order to comprehensively evaluate the clustering performance of the scMSCL method proposed in the present application, the present application carries out experiments on four multi-modal single-cell data sets, including six competitive methods, specifically: scEMC proposes a SAN module combined with a Transformer structure to capture global structural relationships in different feature spaces; scMCs analyzes by mapping omics data from a high-dimensional space to a low-dimensional space through multiple autoencoders; scMVAE proposes three strategies-scMVAE-PoE, scMVAE-NN and scMVAE-Direct, for learning joint latent features for data fusion and clustering, wherein scMVAE-Direct directly connects the original features of each omics, scMVAE-NN fuses low-dimensional features extracted from different omics, and scMVAE-PoE estimates the joint posterior distribution through expert model multiplication; in addition, DCCA projects different omics into corresponding low-dimensional spaces and realizes effective fusion of multi-omics data through a "teacher-student" mechanism. ​

[0181] To verify the superiority of the scMSCL of the present application, NMI and ARI are used as evaluation indexes in this embodiment, and the higher the values are, the better the performance is. The specific evaluation structure is shown in Table 1 below:

[0182] Table 1. Evaluation index result table of six kinds of comparison algorithms

[0183]

[0184] As can be seen from the above table, the indexes NMI and ARI of the scMSCL of the present application both show the most superior performance. Meanwhile, referring to Figures 5-11 , the UMAP clustering visualization of the scMSCL of the present application is also the most compact, and the distance between each class is the largest, showing great advantages. UMAP (Uniform Manifold Approximation and Projection) is a nonlinear dimensionality reduction algorithm, mainly used to map high-dimensional data to low-dimensional space (such as 2D or 3D) while preserving the local and global structure of the data.

[0185] The embodiment of the present application also provides a computer readable storage medium, which stores a computer program. When the computer program is executed by a processor, the processor executes the steps of the single-cell multi-omics clustering method based on multi-subspace contrast learning as described above.

[0186] Referring to Figure 12 , the embodiment of the present application also provides a computer device, which includes a memory and a processor. The memory stores a computer program. When the computer program is executed by the processor, the processor executes the steps of the single-cell multi-omics clustering method based on multi-subspace contrast learning as described above.

[0187] The embodiment of the present application also provides a computer program product containing instructions, which, when running on a computer, causes the computer to execute the steps of the single-cell multi-omics clustering method based on multi-subspace contrast learning as described above.

[0188] It can be understood that the system, device and storage medium provided by the embodiment of the present application correspond to the method provided by the embodiment of the present application, and the explanation, examples and beneficial effects of the related content can be referred to the corresponding part in the single-cell multi-omics clustering method based on multi-subspace contrast learning as described above.

[0189] It should be understood that all or part of the steps in the embodiments of the present application can be implemented by software, hardware, firmware or any combination thereof. When implemented by hardware, all or part of the steps can be implemented in the form of a purchased standard component or a custom component. When implemented by software, all or part of the steps can be implemented in the form of a computer program product. The computer program product includes one or more computer instructions. When the computer program instructions are loaded and executed on a computer, all or part of the processes or functions described in the embodiments of the present application are generated. The computer can be a general purpose computer, a special purpose computer, a computer network, or other programmable device. The computer instructions can be stored in a computer readable storage medium or transmitted from one computer readable storage medium to another computer readable storage medium, for example, the computer instructions can be transmitted from one website, computer, server or data center to another website, computer, server or data center through wired (such as coaxial cable, optical fiber, digital subscriber line (DSL)) or wireless (such as infrared, wireless, microwave, etc.) mode. The computer readable storage medium can be any available medium accessible by a computer or a data storage device such as a server, data center, etc. containing one or more available media. The available media can be magnetic media (such as floppy disk, hard disk, magnetic tape), optical media (such as DVD), or semiconductor media (such as solid state disk (SSD)) and the like.

[0190] In summary, the present application provides a single-cell multi-omics clustering method based on multi-subspace contrast learning to make up for the deficiency of single-space contrast. The method can jointly model single-cell transcriptomic and epigenetic data. First, it can fully exploit the specificity of different omics and the potential common representation between different omics, and then fuse them into a co-embedding representation, and this co-embedding can analyze cell heterogeneity. Specifically, an information extraction and fusion module is designed to finely process individual features and common features of heterogeneous omics data. Through multi-subspace contrast learning, the scMSCL can capture the complementary information between different feature space omics data and the internal relationship between omics data in the same feature space, and build a more comprehensive and information-rich representation for the fusion and clustering of single-cell multi-omics data. A large number of experiments show that the scMSCL outperforms other benchmark methods in performance, demonstrating its superior clustering ability.

[0191] It should be understood that the examples and embodiments described herein are only for illustration and are not intended to limit the present application, and those skilled in the art can make various modifications or changes based on it, and any modifications, equivalent replacements, improvements, etc. made within the spirit and principles of the present application shall be included in the protection scope of the present application.

Claims

1. A single-cell multi-omics clustering method based on multi-subspace contrastive learning, characterized in that, Includes the following steps: S1. Based on the autoencoder, important individual features are extracted from the transcriptome and chromatin accessibility of different omics. S2. Multi-subspace contrastive learning is used to extract and fully explore the common features of transcriptome and chromatin accessibility, effectively revealing the deep-seated relationships between omics data. S3: Employ autoencoders targeting specific omics chromatin accessibility and transcriptomes to reconstruct data for better noise reduction and reconstruction of features consistent with omics information; S4: The individual information of the specific omics in S3 and the common information of multiple omics in S2 are fused to obtain a co-embedding that is more conducive to clustering. S2 further includes: S2.1 Define two views u and r. From the perspective of features, input views u and r into multiple subspace encoders MSEncoder to generate node embeddings for multiple subspaces. S2.2 To make the features captured by multi-subspace contrastive learning more diverse, we should focus on more important node information, make full use of the degree information obtained by the graph structure, and integrate the graph structure information into the attention mechanism to enhance the discriminativeness of the feature representation. To compute the embedding representation of a node within each subspace, the value embedding V is broadcast to each subspace as follows: We obtain the node embeddings of p subspaces, expressed as: ； ； Where l represents the l-th subspace, and p represents the number of subspaces. , and , respectively, are the node embeddings of the i-th subspace in view u and view r, view U and view R are the feature representations obtained by multiple subspace encoders, and are the common features of the same omics in view u and view r, respectively, Att represents the "attention weight matrix", N is the number of samples, and d represents the dimension of the embedding space; S2.3, through calculation and The mutual information maximization between node embeddings across subspaces and the mutual information minimization within subspaces are used as their respective subspace contrastive loss functions. S2.4 The final objective function of subspace contrastive learning is a weighted sum of the mutual information maximization objective and the mutual information minimization objective, which can be formally expressed as: In this formula, L contra Let represent the subspace contrast loss, θ' and φ be the encoder parameters in the view u and view r branches respectively, γ be the weighting coefficient balancing the influence of targets within and between spaces, MI(˙) be the mutual information, and k represent the k-th subspace. and Let represent the node embeddings of the k-th and l-th subspaces in view u, respectively. This represents the node embedding in the k-th subspace of view r; S2.

5. Train and concatenate views U and R to obtain the latent feature representation. .

2. The single-cell multi-omics clustering method based on multi-subspace contrastive learning as described in claim 1, characterized in that, In S1, an autoencoder is used to map single-cell multi-omics data to their respective nonlinear embedding spaces, preserving individual characteristics and enhancing robustness to noise and outliers. Specifically, let's set , These are standardized scRNA-seq and scATAC-seq data, respectively, where Dr and Da are the feature numbers, and two independent encoders are used. and To learn the corresponding d-dimensional feature representation { , } The expressions are as follows: ； In these two formulas, It represents the potential low-dimensional features of cells and genes in scRNA-seq data. It is a potential low-dimensional feature representation between cells and peaks in scATAC-seq data.

3. The single-cell multi-omics clustering method based on multi-subspace contrastive learning as described in claim 1, characterized in that, The process in S2.1 is represented as follows: In these two formulas, It represents the potential low-dimensional features of cells and genes in scRNA-seq data. It is a potential low-dimensional feature representation between cells and peaks in scATAC-seq data; The subspace encoder MSEncoder maps the input feature matrix X to a keyword embedding K and a value embedding V using two independent linear projections, with the corresponding expressions as follows: In these two formulas, , , representing different projection matrices of the key and value, respectively; In S2.2: The initialization expression for the degree embedding matrix S is: ； In this formula, Let D represent the projection matrix corresponding to the degree, where D is the degree matrix; The expression for the attention weight matrix is ​​as follows: ； In this formula, Att is the output of the entire expression, representing the "attention weight matrix", the Softmax(˙) function is a function that "compresses" a vector of arbitrary real values ​​into a probability distribution vector, Q is the query matrix, and α is a learnable or configurable hyperparameter.

4. The single-cell multi-omics clustering method based on multi-subspace contrastive learning as described in claim 3, characterized in that, In step S2.3, the feature representations of views u and r are output to p subspaces, i.e. To optimize the subspace encoder MSEncoder and ensure that views U and R both possess diverse feature representations of common features of the same omics, a pair of loss functions is used. These functions focus on maximizing the mutual information between subspace representations and minimizing the mutual information between different subspace representations. For the mutual information between subspace representations, the objective is formulated as follows: ； In this formula, θ' and φ are the encoder parameters in the branches of view u and view r, respectively, k represents the k-th subspace, p represents the number of subspaces, and MI(˙) is the mutual information. This represents the node embedding in the k-th subspace of view u. This represents the node embedding in the k-th subspace of view r.

5. The single-cell multi-omics clustering method based on multi-subspace contrastive learning as described in claim 4, characterized in that, because and The mutual information (MI) is unknown, making it impossible to calculate. To estimate MI, a sampling method is used instead of maximizing mutual information to maximize the lower bound estimate of MI, and a contrast log-ratio upper bound (CLUB) is employed to minimize the upper bound estimate of MI. Further: Regarding the estimation of the lower bound for maximizing mutual information (MI), the Jensen-Shannon estimation is used instead of the Kullback-Leibler KL divergence estimation. The calculation formula is as follows: In this formula, It is a discriminator that takes the representations of two views as input and scores the consistency between them. Here, the discriminator is simply instantiated as the dot product between two representations, i.e. SP(·) is the activation function. For the joint distribution expectation of two views within the same subspace k, The product of the independent distribution expectations of two views within the same subspace k; In addition to the objectives within each space, constraints are imposed on the relationships between different subspaces within the same view to increase spatial diversity. To this end, an optimization objective is proposed to minimize mutual information. This objective can increase the key features within view u, and because... Therefore, the expression for the optimization objective is: ； In this formula, l represents the l-th subspace. This represents the node embedding in the l-th subspace of view u; Regarding the estimation of the upper limit for minimizing mutual information (MI), the formula for calculating the upper limit of the logarithmic ratio (CLUB) is as follows: ； In this formula, I CLUB Let E represent the expectation of CLUB, P(x,y) represent the joint distribution of x and y, P(x) and P(y) represent the marginal distributions of x and y respectively, log is the natural logarithm, and P(y|x) is the conditional probability of y given x. The key to calculating CLUB is the distribution... In this case, it is assumed that each dimension of x and y is correlated; give Assuming that each dimension of y depends only on the corresponding dimension of x, then using and Substituting x and y respectively, the relevant expression is: ； In this formula, P represents the probability density. Let N(˙) denote the covariance matrix, N(˙) denote the multivariate Gaussian distribution, and I denote the identity matrix. This represents the variance shared across all dimensions. and They are and The node embedding in the i-th subspace, Π i The multiplication symbol is for all elements or subspaces i, representing the product of a high-dimensional Gaussian distribution into a product of independent one-dimensional Gaussian distributions.

6. The single-cell multi-omics clustering method based on multi-subspace contrastive learning as described in claim 1, characterized in that, In S3, a decoder network based on the zero-inflated negative binomial regression model ZINB explores the global probability structure of scRNA-seq data. The ZINB distribution is derived from the mean of the negative binomial distribution. Deviation parameters and coefficients used to describe the probability of interpolated events. The definition and expression are: In these two formulas, Indicates the parameter is u r The negative binomial probability mass function of x and θ. r It is a vector from the raw scRNA-seq data, representing a random variable taking non-negative integer values, u r Let θ be the mean parameter of the distribution, and θ be the deviation parameter. It is the gamma function, ZINB(˙) represents the zero-inflation negative binomial distribution, and π represents the zero-inflation probability. It is x r Indicator function when =0; Specifically, the ZINB-based decoder estimates parameters based on the co-embedding representation Z through three different fully connected layers. : ； ； ； Of these three formulas, yes In matrix form, The proportion of "extra zeros" is determined by the zero-inflation probability parameter. It is the negative binomial mean parameter, representing the expected count of the non-zero portion, and is usually output by the decoder network. It is the negative binomial dispersion parameter, which controls the variance and is usually output by the decoder network. The activation function is sigmoid(·), so the exit probability is between 0 and 1. Furthermore, since the mean and discrete parameters are non-negative, the exponential function exp is chosen as the activation function. and Activation function, It is a decoder with a fully connected layer. , and W θ These are three learnable parameter matrices; Unlike traditional autoencoders based on mean squared error loss, the loss function of the ZINB-based decoder network is the negative logarithm of the ZINB likelihood, expressed as: ； In this formula, The zero-inflated negative binomial distribution is used to model sparse, overly discrete single-cell counting data. -log(˙) represents taking the probability density and then taking the negative log to obtain the negative log-likelihood loss. The smaller this value is, the better. Considering the extremely sparse and almost binary nature of scATAC-seq data, a decoder network based on the Bernoulli distribution (Ber) is chosen to model the scATAC-seq data, expressed as: ； In this formula, These are vectors from the original scATAC-seq data. It is Ber's average parameter; The decoder based on the Bernoulli distribution (Ber) estimates the input latent variable Z' through a fully connected layer using sigmoid (·) as the activation function. The expression is: ； In this formula, is u a In matrix form, sigmoid(·) is an activation function that compresses real values ​​to the range between 0 and 1. It is a weight parameter matrix. This represents the mapping of the a-th decoder; Finally, the decoder based on the Bernoulli distribution (Ber) is optimized using cross-entropy loss, expressed as follows: ； In this formula, L Ber This represents the Bernoulli loss, and CE(·) is the cross-entropy function. These are Bernoulli variables that are actually observed.

7. The single-cell multi-omics clustering method based on multi-subspace contrastive learning as described in claim 1, characterized in that, In S4, two latent representations are learned. and They encode the characteristics of individuals in different omics data, and also represent them through latent features. To capture the commonalities among these features, and in this process, a scale parameter is introduced. and For these two sets of potential representations and To achieve efficient aggregation between the two, we perform element-wise weighted summation, which can be formalized as follows: ； In this formula, Z represents co-embedding representation; The accuracy and reliability of cluster analysis can be improved by generating a more discriminative co-embedding representation Z.

8. The single-cell multi-omics clustering method based on multi-subspace contrastive learning as described in claim 7, characterized in that, In S4, in order to pursue a more discriminative and informative co-embedding representation Z, the total loss is... Combining individual and common information from multiple omics datasets, and unifying the input targets for scRNA-seq and scATAC-seq data, the total loss expression is as follows: ； In this formula, This indicates taking the minimum value. These represent network parameters, where α1, α2, and γ all represent weight coefficients. The zero-inflated negative binomial distribution is used to model sparse, overly discrete single-cell counting data. The proportion of "extra zeros" is determined by the zero-inflation probability parameter. It is the negative binomial mean parameter, representing the expected count of the non-zero portion, and is usually output by the decoder network. This is the negative binomial dispersion parameter, which controls the variance. It is usually output by the decoder network. -log(˙) indicates that the probability density is taken first, and then the negative logarithm is taken to obtain the negative log-likelihood loss. The smaller this value, the better. CE(·) is the cross-entropy function. For the actual observed Bernoulli variables, is u a The matrix form, u a The mean parameters of the Bernoulli distribution Ber are θ' and θ'. These are the encoder parameters in the branches of view u and view r, respectively. MI(˙) is the mutual information, used to measure the degree of nonlinear dependence between two random variables. k represents the k-th subspace, and p represents the number of subspaces. and Let represent the node embeddings of the k-th and l-th subspaces in view u, respectively. This represents the node embedding in the k-th subspace of view r; The total loss function is continuously minimized and trained for 500 epochs for convergence optimization. Individual feature representations and shared feature representations are learned from multi-omics data and then merged into the co-embedding representation Z for clustering information.

Images