Single-cell multi-omics data integration method and system based on domain self-adaption

Abstract

The invention discloses a single-cell multi-omics data integration method and system based on domain self-adaption, and mainly solves the problem of excessive integration in the prior art. According to the implementation scheme, the method comprises the steps that original chromatin accessibility data scATAC-seq is subjected to preprocessing of conversion, labeling and splicing in sequence; constructing an unsupervised data integration model based on a variational auto-encoder and a modal discriminator, taking a relative entropy KL and a binary cross entropy BCE as a loss function L, and introducing a domain specific batch normalization layer DSBN in the variational auto-encoder; inputting the preprocessed data into an unsupervised data integration model for training; and inputting the preprocessed data into a trained unsupervised data integration model to obtain a unified low-dimensional embedded vector integrated by two different omics data. According to the method, the homogeneity between different omics can be learned in the integration process, the heterogeneity between different omics can be effectively reserved, the efficiency of large-scale data set integration is improved, and the method can be used for cell type label migration.

Description

Technical Field

[0001] The present invention belongs to the technical field of information processing, and further relates to a method and system for integrating single-cell multi-omics data, which can be used for cell type label transfer. Background Art

[0002] The rapid development of single-cell multi-omics sequencing technology has generated various single-cell multi-omics data: single-cell transcriptome data (scRNA-seq) can determine the gene expression levels of individual cells, single-cell chromatin accessibility data (scATAC-seq) can evaluate the open state and regulatory regions of single-cell chromatin, and single-cell spatial transcriptome data (MERFISH, seqFISH, etc.) can reveal the spatial distribution of single cells in tissues. These data can reveal the characteristics of single cells from different perspectives, fundamentally changing the study of complex biological processes such as the immune system, embryonic development, and cancer. Integrating single-cell multi-omics data can provide a deeper understanding of biological systems than single-omics data. For example, a study found that through integration, researchers successfully refined the cell clusters of kidney cell types labeled by scRNA-seq using scATAC-seq, revealing more cell clusters with potential clinical significance. In another study, by integrating scRNA-seq and scATAC-seq, researchers could infer upstream regulatory factors from cis-element motif analysis, thereby identifying cell type-specific regulatory networks. In this study, the cis-regulatory elements identified using scATAC-seq overcame the difficulty of detecting regulatory genes, such as transcription factors, due to low abundance in transcriptome data. These studies have all shown that integrating multiple modalities of omics information, namely single-cell multi-omics data, can achieve combinatorial inferences that cannot be obtained from any single omics. Therefore, the integration of single-cell multi-omics data has the potential to obtain more information compared to single-omics data.

[0003] Cao et al. proposed a method to use optimal transport theory to close the distance between different omics data in low-dimensional space in their paper “A unified computational framework for single-cell data integration with optimal transport” (Nature Communications, 2022, 13(1): 7419.). The implementation steps of this method are as follows: first, the data is screened, the top 2000 dimensional highly variable features in different omics are retained as the private features of the omics, and the top 2000 dimensional highly variable features shared between different omics are input into the model as the shared features. At the decoder, the shared features and private features of the omics are simultaneously reconstructed as part of the loss function, and the similarity between the two omics distributions is calculated based on the optimal transport theory in the low-dimensional embedding space as the other part of the loss function. After the model converges, the low-dimensional embedding space of multi-omics data is obtained. Although this method performs well in modal mixing, it overemphasizes the unique cell subtypes and different biological conditions between different omics because it over-emphasizes fitting the distributions between different omics, which not only leads to over-integration, but also makes it difficult to extend to multi-omics data with strong heterogeneity, and has poor generalization.

[0004] Cao et al. proposed a data integration method based on variational autoencoders and graph neural networks in their paper “Multi-omics single-cell data integration and regulatory inference with graph-linked embedding” (Nature Biotechnology, 2022, 40(10): 1458-1466). The implementation steps of this method are to connect the characteristic peaks of scATAC-seq and the characteristic genes of scRNA-seq according to the proximity of the genome to construct a gene regulatory prior graph. Finally, the variational autoencoder VAE and the graph neural network GNN are used to simultaneously reconstruct the expression matrix and the regulatory graph, learn the representation of cells from different modal data, and align them through an iterative process. This method has two shortcomings. First, it relies too much on the construction of the prior regulatory graph. If the quality of the prior regulatory graph is not high, the integration effect will be seriously affected. It is only applicable to multimodal data and performs poorly on multi-omics data with strong heterogeneity. Second, there is a problem of over-integration, and training is slow when the data volume is large, which limits its ability to be expanded to large-scale data sets. Summary of the Invention

[0005] The object of the present invention is to propose a method and system for integrating single-cell multi-omics data based on domain adaptation in view of the deficiencies of the above-mentioned existing technologies based on the domain adaptation theory, so as to avoid over-integration and introduction of incorrect information in data integration, improve the training speed, and enhance the efficiency of integrating large-scale data sets.

[0006] The idea to achieve the object of the present invention is as follows: By introducing a domain-specific batch normalization (DSBN) layer at the decoder, and performing independent batch normalization (BN) processing on data with different omics labels when using the decoder to reconstruct the expression matrix, the model can retain the private features between different omics during the integration process, avoid the over-integration problem, and effectively identify and retain the heterogeneity between different omics such as cell subtypes and biological condition differences; at the same time, by dynamically adjusting the intensity α value of the modality discriminator during the training process, the training speed is increased, and the efficiency of integrating large-scale data sets is enhanced.

[0007] Based on the above idea, the technical solution of the present invention is as follows:

[0008] 1. A method for integrating single-cell multi-omics data based on domain adaptation, characterized by comprising:

[0009] 1) Preprocess the original chromatin accessibility data scATAC-seq by sequentially performing transformation, labeling, and splicing.

[0010] 2) Construct an unsupervised data integration model including a variational autoencoder and a modality discriminator, use relative entropy KL and binary cross-entropy BCE as its loss function L, and introduce a domain-specific batch normalization layer DSBN in the variational autoencoder to retain the specificity between omics.

[0011] The variational autoencoder is used to learn a low-dimensional embedding vector that can represent the input data.

[0012] The modality discriminator is used to discriminate the data type from which the low-dimensional embedding vector comes.

[0013] 3) Input the preprocessed data into the unsupervised data integration model for training, iteratively update the network parameters until the loss function of the network converges.

[0014] 4) Input the preprocessed data into the trained unsupervised data integration model to obtain a low-dimensional embedding vector.

[0015] Furthermore, the variational autoencoder includes an encoder and a decoder.

[0016] The encoder consists of two cascaded composite encoding layers. Each composite encoding layer includes a linear layer, a batch normalization (BN) layer, and an activation layer. The input and output dimensions of the two composite encoding layers are set to 2000-1024 and 1024-64 respectively, for obtaining a 64-dimensional low-dimensional embedding vector.

[0017] The decoder consists of a linear layer, a domain-specific batch normalization (DSBN) layer, and an activation layer in series. The input and output dimensions of the decoder are 64-2000, for reconstructing the expression profile information.

[0018] Furthermore, the modality discriminator consists of five cascaded composite discrimination layers. Each composite discrimination layer includes a linear layer, an activation layer, and a dropout layer. The input and output dimensions of the five composite discrimination layers are set to 64-64, 64-64, 64-32, 32-16, and 16-1 respectively, for obtaining a data type prediction vector.

[0019] 2. A single-cell multi-omics data integration system based on domain adaptation, characterized by comprising:

[0020] A preprocessing module, used for converting the original chromatin accessibility data scATAC-seq into a gene activity matrix; labeling different tags for scRNA-seq and scATAC-seq data; splicing according to the overlapping genes between the gene expression data scRNA-seq and scATAC-seq data to obtain preprocessed data.

[0021] A variational autoencoder, used for first encoding and mapping an input vector into a low-dimensional embedding vector, and then decoding the low-dimensional embedding vector back to the original input vector to obtain a low-dimensional embedding vector that can represent the input data.

[0022] A modality discriminator, used for mapping the low-dimensional embedding vector into a data type prediction vector and judging which data type the low-dimensional embedding vector comes from.

[0023] An integration module, used for constructing an unsupervised data integration model including a variational autoencoder and a modality discriminator, and connecting the modality discriminator to the middle layer of the variational autoencoder to construct an unsupervised data integration model.

[0024] A training module, used for iteratively training the unsupervised data integration model, so that the variational autoencoder can obtain a low-dimensional embedding vector that represents the input data, while the modality discriminator cannot identify which data type the input low-dimensional embedding vector is, and obtaining a variational autoencoder network that can extract potential features.

[0025] A visualization module, used for evaluating the quality of the integration result, reducing the low-dimensional embedding vector to a two-dimensional vector and visualizing it.

[0026] The present invention has the following advantages compared with the prior art:

[0027] First, the present invention introduces a domain-specific batch normalization (DSBN) layer in the decoder part of the variational autoencoder, which enables the decoder to perform independent batch normalization (BN) processing on data with different omics tags when reconstructing the expression matrix. At the same time, during the training process, the intensity α value of the modality discriminator can be dynamically adjusted, enabling the model to retain the private features between different omics during the integration process, overcoming the over-integration problem of the prior art, so as to effectively identify and retain the heterogeneity between different omics, such as cell subtypes and biological condition differences, during the data integration process.

[0028] Second, since the present invention uses a common encoder in the variational autoencoder to extract the common features between different data sources and introduces a modality discriminator to accelerate the model convergence, it overcomes the problem of slow training in the prior art and significantly improves the efficiency of integrating large-scale datasets. BRIEF DESCRIPTION OF THE DRAWINGS

[0029] Figure 1 It is a flowchart of the implementation of the single-cell multi-omics data integration method provided in the first embodiment of the present invention;

[0030] Figure 2 It is a block diagram of the structure of the single-cell multi-omics data integration system provided in the second embodiment of the present invention;

[0031] Figure 3 It is a comparison chart of the integration results of the present invention and several existing data integration methods on multi-omics data of heterogeneous mouse spleens;

[0032] Figure 4 For Figure 3 It is an index evaluation chart of the integration results;

[0033] Figure 5 It is a comparison chart of the integration results of the method of the present invention and several existing data integration methods on a dataset with biological condition differences;

[0034] Figure 6 For Figure 5 It is an index evaluation chart of the integration results. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0035] In order to enable those skilled in the art to better understand the solution of the present invention, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0036] It should be noted that the step numbers in the specification and claims of the present invention are only for clearly describing the implementation embodiments of the present invention for easy understanding, and their sequence numbers are not limited.

[0037] Example 1, A Single-Cell Multi-Omics Data Integration Method Based on Domain Adaptation

[0038] Refer to Figure 1 , the implementation steps of this example are as follows:

[0039] Step 1, Obtain single-cell multi-omics data and perform preprocessing.

[0040] 1.1) Download single-cell multi-omics data from an existing database, which includes gene expression data scRNA-seq and chromatin accessibility data scATAC-seq. Among them, the scRNA-seq is in the form of a two-dimensional matrix, with the rows being cell names cell and the columns being gene names gene. The scATAC-seq is in the form of a two-dimensional matrix, with the rows being cell names cell and the columns being genomic region names peak;

[0041] 1.2) Perform preprocessing on the single-cell multi-omics data:

[0042] 1.2.1) Retain the genes that are expressed in at least 3 cells and the cells with at least 200 genes expressed in the dataset, and convert the chromatin accessibility scATAC-seq data into a gene activity matrix;

[0043] 1.2.2) Label different data types, that is, set the label of the gene expression data scRNA-seq to 0 and the label of scATAC-seq to 1;

[0044] 1.2.3) Select the common genes in the scRNA-seq and scATAC-seq data, and splice the scRNA-seq and scATAC-seq levels together according to the selected common genes to obtain the preprocessed data.

[0045] Step 2, Construct an unsupervised data integration model.

[0046] 2.1) Establish a variational autoencoder composed of a series connection of an encoder and a decoder:

[0047] The encoder is composed of two cascaded composite encoding layers, with its input and output dimensions set to 2000-1024 and 1024-64 respectively, and is used to obtain a 64-dimensional low-dimensional embedding vector; each composite encoding layer includes a linear layer, a batch normalization (BN) layer, and an activation layer; the linear layer is implemented by the torch.nn.Linear function and is used to extract the low-dimensional embedding vector; the batch normalization (BN) layer is implemented by torch.nn.BatchNorm1d and is used to stabilize the values of the hidden layer during training and accelerate the training speed; the activation layer is implemented using the torch.nn.relu function;

[0048] The decoder is composed of a linear layer, a domain-specific batch normalization (DSBN) layer, and an activation layer in series, with its input and output dimensions being 64-2000 and is used to reconstruct the expression profile information; the linear layer is implemented by the torch.nn.Linear function and is used to reconstruct the expression profile information; the domain-specific batch normalization (DSBN) layer is based on torch.nn.BatchNorm1d and is used to retain the specificity between omics; the activation layer is implemented using torch.nn.Sigmoid;

[0049] 2.2) Establish a modality discriminator composed of five cascaded composite discriminant layers:

[0050] Each composite discriminant layer includes a linear layer, an activation layer, and a dropout layer. The input and output dimensions of the five composite discriminant layers are set to 64-64, 64-64, 64-32, 32-16, and 16-1 respectively, and are used to obtain a data type prediction vector, so that the modality discriminator cannot identify which data type the input low-dimensional embedding vector is, thereby accelerating the capture of shared features between different data; the linear layer is implemented by the torch.nn.Linear function and is used to extract features; the activation layer is implemented using torch.nn.LeakyReLU; the dropout layer is implemented by the torch.nn.dropout function, where the parameter dropout_rate is set to 0.5;

[0051] 2.3) Connect the modality discriminator to the middle layer of the variational autoencoder to construct an unsupervised data integration model;

[0052] 2.4) Use relative entropy KL and binary cross-entropy BCE as the loss function L of the unsupervised data integration model, which is expressed as follows:

[0053]

[0054] where D KL(q(z|x)||p(z)) represents the difference between the posterior distribution q(z|x) and the prior distribution p(z). Here, x represents the original expression matrix, z represents the low-dimensional embedding vector, q(z|x) represents the posterior distribution of the low-dimensional embedding vector z obtained from the input x, and p(z) represents the Gaussian distribution;

[0055] represents the difference between the predicted value and the actual value, represents the predicted value, and y represents the actual value.

[0056] Step 3: Iteratively train the unsupervised data integration model.

[0057] 3.1) Obtain the low-dimensional embedding vector:

[0058] Randomly sample the preprocessed data and input it in batches into the encoder module of the variational autoencoder. Map the input data to the mean μ and variance σ. Obtain the low-dimensional embedding vector by sampling the mean μ and variance σ, and calculate the current relative entropy KL’ between N(μ, σ) and the Gaussian distribution N(0, 1). The formula is as follows:

[0059]

[0060] where q(z|x) represents the posterior distribution N(μ, σ 2 I) of the low-dimensional embedding vector z obtained from the original expression matrix x, and p(z) represents the Gaussian distribution N(0, I), aiming to calculate the difference between the posterior distribution q(z|x) and the prior distribution p(z).

[0061] 3.2) Data type prediction:

[0062] Input the low-dimensional embedding vector into the modality discriminator, and map it to a data type prediction vector through the modality discriminator. Calculate the current predicted binary cross-entropy BCE1’ between the data type prediction vector and the true label:

[0063]

[0064] where, X r represents the gene expression value of scRNA-seq, X a represents the gene activity matrix expression value of scATAC-seq, E represents the encoder, and C represents the modality discriminator, aiming to calculate which data type the input low-dimensional embedding vector comes from;

[0065] 3.3) Reconstruct the expression matrix:

[0066] Input the low-dimensional embedding vectors into a decoder with a DSBN layer to perform batch normalization on the low-dimensional embedding vectors of different labels using different BN layers, and reconstruct the normalized low-dimensional embedding vectors into a reconstructed expression matrix with the same dimension as the original input data, and calculate the current reconstruction binary cross-entropy BCE2' between the reconstructed expression matrix and the original expression matrix:

[0067]

[0068] where x represents the original expression matrix, represents the reconstructed expression matrix, and N represents the number of samples, aiming to calculate the difference between the reconstructed expression matrix and the original expression matrix;

[0069] 3.4) Calculate the loss and update the network parameters:

[0070] 3.4.1) Calculate the total loss function L' of the loss through forward propagation:

[0071] L' = λ1 * BCE1' + λ2 * KL' + λ3 * BCE2'

[0072] where λ1, λ2, and λ3 respectively represent the weight parameters of each part, and the default settings are 55, 0.5, and 2000;

[0073] 3.4.2) Use the Adaptive Moment Estimation (ADAM) optimizer to calculate its gradient and perform backpropagation of the gradient to update the network parameters;

[0074] 3.5) Repeat steps 3.1)-3.4) until the loss function of the network converges to obtain a trained variational autoencoder network.

[0075] Step 4, input the preprocessed single-cell multi-omics data into the trained variational autoencoder, extract the shared features between the two omics through encoding, and obtain a unified low-dimensional embedding vector integrating the two different omics datasets.

[0076] Example 2, A Single-Cell Multi-Omics Data Integration System Based on Domain Adaptation

[0077] Refer to Figure 2 , this example includes: a preprocessing module 1, a variational autoencoder 2, a modality discriminator 3, an integration module 4, a training module 5, and a visualization module 6. Among them:

[0078] A preprocessing module 1, which is used to convert the original chromatin accessibility data scATAC-seq into a gene activity matrix; label different tags for scRNA-seq and scATAC-seq data; splice according to the overlapping genes between the gene expression data scRNA-seq and scATAC-seq data to obtain preprocessed data, and input it to the variational autoencoder 2;

[0079] The variational autoencoder 2 encodes and maps the input vector into a low-dimensional embedding vector, and then decodes the low-dimensional embedding vector back to the original input vector to obtain a low-dimensional embedding vector that can represent the input data, and inputs it to the modality discriminator 3;

[0080] The modality discriminator 3 is used to map the low-dimensional embedding vector into a data type prediction vector and determine which data type the low-dimensional embedding vector comes from;

[0081] The integration module 4 is used to construct an unsupervised data integration model including the variational autoencoder 2 and the modality discriminator 3, and connect the modality discriminator 3 to the middle layer of the variational autoencoder 2 to form an unsupervised data integration model;

[0082] The training module 5 is used to iteratively train the unsupervised data integration model, so that the variational autoencoder 2 can obtain a low-dimensional embedding vector that represents the input data, while the modality discriminator 3 cannot identify which data type the input low-dimensional embedding vector is, and obtain a trained variational autoencoder network that can extract latent features;

[0083] The visualization module 6 is used to input the preprocessed data into the trained variational autoencoder network for mapping, obtain two different omics datasets with the low-dimensional embedding vector reduced to a two-dimensional space integrated into a unified low-dimensional embedding vector, and perform visualization.

[0084] The effect of the present invention can be further demonstrated by the following simulation:

[0085] 1. Simulation experiment conditions.

[0086] The software platform for the simulation experiment of the present invention is: WINDOWS10 operating system, Python 3.8 and Rstudio integrated development software.

[0087] The data used in Simulation Experiment 1 of the present invention is: the mouse spleen single-cell multi-omics dataset provided by the paper "MultiMAP: dimensionality reduction and integration of multimodal data". Among them, scRNA-seq contains 4,382 cells and 13,575 genes; after scATAC-seq is converted into a gene activity matrix, 3,166 cells and 19,410 genes are retained; scRNA-seq has 12 cell types, scATAC-seq has 12 cell types, and there are 8 cell types in common between the two omics.

[0088] The data used in Simulation Experiment 2 of the present invention is: the human PBMC dataset GSE156478 under different biological condition stimulations. scRNA-seq contains 4,464 cells and 17,441 genes; after scATAC-seq is converted into a gene activity matrix, 3,574 cells and 17,441 genes are retained; among them, scRNA-seq and scATAC-seq have 6 identical cell types, and the gene expression levels of two cell types between the two omics change due to biological condition stimulations. During the integration process, the Naive CD4+ T cells and Naive CD8+ T cells corresponding to the two omics should not be mixed.

[0089] 2. Simulation content and result analysis.

[0090] In Simulation Experiment 1, the data integration was performed on the dataset described in Simulation Dataset 1 above using 9 integration methods of the present invention and the prior art, namely MaxFuse, GLUE, Seurat, Harmony, SCALEX, uniPort, MultiMAP, Scanorama, and scCorrector, to obtain the visual comparison results of the cell type and modality mixing degree as Figure 3 shown, and the index comparison results as Figure 4 shown. Among them:

[0091] Figure 3 a is the visualization result graph of the cell type and modality mixing degree of the original data;

[0092] Figure 3 b is the visualization result graph of the cell type and modality mixing degree of the present invention;

[0093] Figure 3 c is the visualization result graph of the cell type and modality mixing degree of the existing SCALEX method;

[0094] Figure 3d is a visualization result graph of the cell type and modality mixing degree of the existing scCorrector method;

[0095] Figure 3 e is a visualization result graph of the cell type and modality mixing degree of the existing uniPort method;

[0096] Figure 3 f is a visualization result graph of the cell type and modality mixing degree of the existing MultiMAP method;

[0097] Figure 3 g is a visualization result graph of the cell type and modality mixing degree of the existing MaxFuse method;

[0098] Figure 3 h is a visualization result graph of the cell type and modality mixing degree of the existing GLUE method;

[0099] Figure 3 i is a visualization result graph of the cell type and modality mixing degree of the existing Seurat method;

[0100] Figure 3 j is a visualization result graph of the cell type and modality mixing degree of the existing Harmony method;

[0101] Figure 3 k is a visualization result graph of the cell type and modality mixing degree of the existing Scanorama method;

[0102] Figure 4 a is a comparison graph of the adjusted Rand index ARI, mutual information NMI, and F1-score metrics of the present invention and the existing technologies used;

[0103] Figure 4 b is a comparison graph of the silhouette coefficient Silhouette and batch entropy Batch Entropy metrics of the present invention and the existing technologies used;

[0104] Figure 3 For each method in, the two subgraphs, the left subgraph represents the visualization result colored by cell type, each point represents a cell, and each color corresponds to a cell type; the right subgraph represents the visualization result colored by omics label, each point represents a cell, and each color corresponds to a data type.

[0105] From Figure 3It can be seen that: the present invention can not only completely integrate the common cell types together, but also retain the different cell types between different omics without mixing; while other methods all have over-integration during the integration process, that is, the different cell types between different omics are mixed together, or under-integration, so that the common cell types between different omics are not completely integrated.

[0106] From Figure 4 (a) It can be seen that: among all the comparison methods, the adjusted Rand index ARI, mutual information NMI, and F1 score of the present invention are all higher than those of the existing comparison methods, proving that the present invention can accurately integrate single-cell multi-omics data;

[0107] From Figure 4 (b) It can be seen that: the silhouette score of the present invention is the highest, indicating that the present invention can maintain a good clustering effect after integration, that is, tight within clusters and separated between clusters. The reason why the batch entropy score of the present invention is not high is that the heterogeneity between different omics is successfully retained, proving that the present invention can not only capture the homogeneity between different omics during integration, but also retain the heterogeneity between different omics.

[0108] In simulation experiment 2, the present invention and 7 existing integration methods, namely Seurat, Harmony, SCALEX, uniPort, MultiMAP, Scanorama, and scCorrector, are used to perform data integration on the dataset described in the above simulation dataset 2, and the visual comparison results of the cell type and modality mixing degree are as Figure 5 shown, and the index comparison results are as Figure 6 shown. Among them:

[0109] Figure 5 a is the visualization result diagram of the cell type and modality mixing degree of the original data;

[0110] Figure 5 b is the visualization result diagram of the cell type and modality mixing degree of the present invention;

[0111] Figure 5 c is the visualization result diagram of the cell type and modality mixing degree of SCALEX;

[0112] Figure 5 d is the visualization result diagram of the cell type and modality mixing degree of scCorrector;

[0113] Figure 5 e is the visualization result diagram of the cell type and modality mixing degree of uniPort;

[0114] Figure 5f is a visualization result graph of the cell type and modality mixing degree of MultiMAP;

[0115] Figure 5 g is a visualization result graph of the cell type and modality mixing degree of Seurat;

[0116] Figure 5 h is a visualization result graph of the cell type and modality mixing degree of Harmony;

[0117] Figure 5 i is a visualization result graph of the cell type and modality mixing degree of Scanorama;

[0118] Figure 6 a is a comparison graph of the adjusted Rand index ARI, mutual information NMI, and F1-score metrics between the present invention and the existing technologies used;

[0119] Figure 6 b is a comparison graph of the Silhouette coefficient and Batch Entropy metrics between the present invention and the existing technologies used;

[0120] From Figure 5 it can be seen that: the present invention has successfully integrated four identical cell types in different omics. The Naive CD4+T and Naive CD8+T cells circled by red and blue frames in the figure have obvious differences after integration, while other technologies have problems of over-integration and under-integration, proving that the present invention can not only capture the homogeneity between different omics during integration, but also retain the heterogeneity between different omics.

[0121] From Figure 6 (a) it can be seen that: among all the comparison technologies, the adjusted Rand index ARI, mutual information NMI, and F1-score of the present invention are higher than those of the existing comparison methods, proving that the present invention can accurately integrate single-cell multi-omics data of different biological conditions;

[0122] From Figure 6 (b) it can be seen that: the Silhouette coefficient and Batch Entropy of the present invention both have great advantages, proving that the present invention is also applicable to the data integration of different biological conditions, making the visualization results after integration more biologically meaningful.

[0123] The above simulation experiments show that: the present invention can accurately integrate single-cell multi-omics data, and can retain the heterogeneity between different omics while capturing the homogeneity between different omics.

[0124] The sources of the existing technologies in the above simulation experiments are as follows:

[0125] MaxFuse is a method proposed by Chen et al. in "Integration of spatial and single-cell data across modalities with weakly linked features, in: Nature biotechnology, 2024, 42(7): 1096-1106."

[0126] scCorrector is a method proposed by Guo et al. in "scCorrector: a robust method for integrating multi-study single-cell data, in: Briefings in Bioinformatics, 2024, 25(2): bbad525."

[0127] uniPort is a method proposed by Cao et al. in "A unified computational framework for single-cell data integration with optimal transport, in: Nature Communications, 2022, 13(1): 7419."

[0128] Seurat is a method proposed by Stuart T et al. in "Comprehensive integration of single-cell data, in: cell, 2019, 177(7): 1888-1902.e21."

[0129] GLUE is a method proposed by Cao et al. in "Multi-omics single-cell data integration and regulatory inference with graph-linked embedding, in: Nature Biotechnology, 2022, 40(10): 1458-1466."

[0130] SCALEX is a method proposed by Xiong et al. in "Online single-cell data integration through projecting heterogeneous datasets into a common cell-embedding space, in: Nature Communications, 2022, 13(1): 6118."

[0131] Harmony is a method proposed by Korsunsky et al. in "Fast, sensitive and accurate integration of single-cell data with Harmony, in: Nature methods, 2019, 16(12): 1289-1296."

[0132] MultiMAP is a method proposed by Jain et al. in "MultiMAP: dimensionality reduction and integration of multimodal data, in: Genome biology, 2021, 22: 1-26."

[0133] Scanorama is a method proposed by Hie et al. in "Efficient integration of heterogeneous single-cell transcriptomes using Scanorama, in: Nature biotechnology, 2019, 37(6): 685-691."

[0134] The calculation formulas of the evaluation metrics in the above simulation experiments are as follows:

[0135] a) Adjusted Rand Index ARI:

[0136]

[0137] where RI is the Rand Index, E[RI] is the expected value of the Rand Index, max(RI) is the maximum value of the Rand Index, a is the number of sample pairs that are in the same cluster in both clusterings, b is the number of sample pairs that are not in the same cluster in both clusterings, c is the number of sample pairs that are in the same cluster in the first clustering but not in the same cluster in the second clustering, and d is the number of sample pairs that are not in the same cluster in the first clustering but are in the same cluster in the second clustering.

[0138] b) Normalized Mutual Information NMI:

[0139]

[0140] H(U) = -∑ u∈U P(u) log P(u)

[0141] where U and V represent two different clustering results respectively, P(u, v) represents the probability that a sample belongs to clusters u and v simultaneously, and P(u) and P(v) are the probabilities that a sample belongs to clusters u and v respectively.

[0142] c) F1 score:

[0143]

[0144] where Precision represents the accuracy and recall represents the recall rate.

[0145] d) Silhouette coefficient

[0146]

[0147] where a is the average within-cluster distance and b is the average nearest-cluster distance (the distance to the nearest cluster that the sample does not belong to). The score is scaled between 0 and 1, with 1 indicating a perfect clustering and 0 indicating that the samples belong to different clusters.

[0148] e) Batch Entropy

[0149]

[0150] where pi is the proportion of cells from modality i in a given region, and Pi is the proportion of cells from modality i relative to the total number of cells.

[0151] In the above two simulation experiments, the visualization result of the integrated cell type and modality mixing degree refers to: inputting the preprocessed data into the trained variational autoencoder to obtain the low-dimensional embedding vectors of all cells, and using the standard process in scanpy to find neighbors and perform UMAP dimensionality reduction, that is, using the scanpy.pp.neighbors function to find neighbors, where the parameter n_neighbors is set to 30; using the scanpy.tl.umap function to perform dimensionality reduction, where the parameter min_dist is set to 0.1.

Claims

1. A method for integrating single-cell multi-omics data based on domain adaptation, characterized in that, including (1) Preprocess the original chromatin accessibility data scATAC-seq by sequentially performing transformation, tagging, and splicing (2) Construct an unsupervised data integration model including a variational autoencoder and a modality discriminator, using relative entropy KL and binary cross-entropy BCE as its loss function L, and introducing a domain-specific batch normalization layer DSBN in the variational autoencoder to preserve the specificity between omics; The variational autoencoder is used to learn a low-dimensional embedding vector that can represent the input data; The modality discriminator is used to discriminate the data type from which the low-dimensional embedding vector comes; (3) Input the preprocessed data into the unsupervised data integration model for training, iteratively update the network parameters until the loss function of the network converges; (4) Input the preprocessed data into the trained unsupervised data integration model to obtain a low-dimensional embedding vector.

2. The method according to claim 1, wherein In the above (1), the preprocessing of sequentially performing transformation, tagging, and splicing on the original chromatin accessibility data scATAC-seq includes: (1a) Select a 2kb region upstream and downstream of the transcription start site TSS, and assign peaks to the genes that may regulate or are closest based on the gene locus, that is, assign the characteristic peaks of scATAC-seq to the adjacent genes, to obtain a gene activity matrix converted from scATAC-seq data; (1b) Tag different data types, that is, set the tag of gene expression data scRNA-seq to 0 and the tag of scATAC-seq to 1; (1c) Select the common genes in scRNA-seq and scATAC-seq data, and splice scRNA-seq and scATAC-seq levels together according to the selected common genes.

3. The method according to claim 1, characterized in that, The variational autoencoder in the above (2) includes an encoder and a decoder; The encoder is composed of two cascaded composite encoding layers. Each composite encoding layer includes a linear layer, a batch normalization BN layer, and an activation layer. The input and output dimensions of the two composite encoding layers are set to 2000 - 1024 and 1024 - 64 respectively, for obtaining a 64-dimensional low-dimensional embedding vector; The decoder is composed of a linear layer, a domain-specific batch normalization DSBN layer, and an activation layer in series. The input and output dimensions of the decoder are 64 - 2000, for reconstructing the expression profile information.

4. The method according to claim 1, wherein The modality discriminator in the above (2) is composed of five cascaded composite discrimination layers. Each composite discrimination layer includes a linear layer, an activation layer, and a dropout layer. The input and output dimensions of the five composite discrimination layers are set to 64 - 64, 64 - 64, 64 - 32, 32 - 16, and 16 - 1 respectively, for obtaining a data type prediction vector.

5. The method according to claim 1, characterized in that, The loss function L in the above (2) is expressed as follows: Among them, D KL (q(z|x)||p(z)) represents the difference between the posterior distribution q(z|x) and the prior distribution p(z). Here, x represents the original expression matrix, z represents the low-dimensional embedding vector, q(z|x) represents the posterior distribution of the low-dimensional embedding vector z obtained by inputting x, and p(z) represents a Gaussian distribution; Indicates the difference between the predicted value and the actual value, where represents the predicted value and y represents the actual value.

6. The method according to claim 1, characterized in that The training of the unsupervised data integration model in the above (3) is implemented as follows: 3a) Set the learning rate of the optimizer to 0.00022, the number of iterations to 330, and other parameters to default values; 3b) Input the preprocessed data into the encoder part of the variational autoencoder. Use the encoder to map the input data to the mean μ and variance σ, and calculate the current relative entropy KL’ between N(μ, σ) and the Gaussian distribution N(0, 1); 3c) Sample the mean and variance to obtain a low-dimensional embedding vector, and input the vector into the modality discriminator. Map it to a data type prediction vector through the modality discriminator, and calculate the current predicted binary cross-entropy BCE1’ between the data type prediction vector and the true label; 3d) Sample the mean and variance to obtain a low-dimensional embedding vector, and input it into the decoder with DSBN layers to perform batch normalization on the low-dimensional embedding vectors of different labels using different BN layers; 3e) Reconstruct the normalized low-dimensional embedding vector into a reconstructed expression matrix with the same dimension as the original input data, and calculate the current reconstructed binary cross-entropy BCE2’ between the reconstructed expression matrix and the original expression matrix; 3f) Calculate the total loss function L’ = BCE1’ + BCE2’ + KL’ through forward propagation, then use the Adaptive Moment Estimation (ADAM) optimizer to calculate its gradient, and perform backpropagation to update the network parameters; 3g) Repeat steps 3b)-3f) until the loss function of the network converges, and obtain the trained variational autoencoder network.

7. The method according to claim 6, characterized in that The current relative entropy KL’ between N(μ, σ) and the Gaussian distribution N(0, 1) calculated in step 3b) is as follows: Among them, q(z|x) represents the posterior distribution N(μ, σ 2 I) of the low-dimensional embedded vector z obtained from the original expression matrix x, and p(z) represents the Gaussian distribution N(0, I), aiming to calculate the difference between the posterior distribution q(z|x) and the prior distribution p(z).

8. The method according to claim 6, characterized in that, The current predicted binary cross-entropy BCE1’ between the data type prediction vector and the true label calculated in step 3c) is as follows: Among them, X r represents the gene expression value of scRNA-seq, and X a represents the gene activity matrix expression value of scATAC-seq. E represents the encoder, and C represents the modality discriminator, which aims to calculate which data type the input low-dimensional embedding vector comes from.

9. The method according to claim 6, characterized in that, The current reconstructed binary cross-entropy BCE2’ between the reconstructed expression matrix and the original expression matrix calculated in step 3e) is as follows: where x represents the original expression matrix, represents the reconstructed expression matrix, and N represents the number of samples, aiming to calculate the difference between the reconstructed expression matrix and the original expression matrix.

10. A single-cell multi-omics data integration system based on domain adaptation, characterized in that, It includes: A preprocessing module for converting the original chromatin accessibility data scATAC-seq into a gene activity matrix; Label different tags for scRNA-seq and scATAC-seq data; splice according to the overlapping genes between the gene expression data scRNA-seq and scATAC-seq data to obtain the preprocessed data; A variational autoencoder for first encoding and mapping the input vector into a low-dimensional embedding vector, and then decoding the low-dimensional embedding vector back to the original input vector to obtain a low-dimensional embedding vector that can represent the input data; A modality discriminator for mapping the low-dimensional embedding vector to a data type prediction vector and determining which data type the low-dimensional embedding vector comes from; An integration module for constructing an unsupervised data integration model including a variational autoencoder and a modality discriminator, and connecting the modality discriminator to the middle layer of the variational autoencoder to construct an unsupervised data integration model; A training module for iteratively training the unsupervised data integration model so that the variational autoencoder can obtain a low-dimensional embedding vector representing the input data, while the modality discriminator cannot identify which data type the input low-dimensional embedding vector is, and obtain a variational autoencoder network that can extract latent features; A visualization module for reducing the dimensionality of the low-dimensional embedding vectors output by the variational autoencoder network obtained by the training module to two-dimensional space and visualizing them.

Images