Method for generating single cell sequencing data

Abstract

The invention provides a single cell sequencing data generation method. The method comprises the following steps: S1, constructing a deep neural network model; s2, an scRNA-seq comprehensive data set is obtained; s3, the scRNA-seq comprehensive data set is subjected to preprocessing, and the scRNA-seq comprehensive data set is obtained; s4, training a deep neural network model on the basis of the preprocessed scRNA-seq comprehensive data set; and S5, based on the trained deep neural network model, carrying out noise addition, denoising and reconstruction on the features of the to-be-detected data to generate single cell sequencing data. According to the method, high-quality scRNA-seq data can be generated by the deep neural network model on data sets of various sequencing platforms. Compared with other latest data generation models, the deep neural network model shows more excellent performance on a plurality of evaluation indexes. In addition, the deep neural network model can simulate single cell data under a pseudo time scale, and provides high-quality data support for tracking cell differentiation and development trajectories, analyzing communication between cells, revealing cell heterogeneity and other analysis.

Description

Technical Field

[0001] The present invention relates to a method for generating single-cell sequencing data, and belongs to the technical field of bioinformatics. Background Art

[0002] Cells are the basic structural and functional units of organisms. In order to study complex biological systems such as developing embryos, tumors, and microorganisms, it is necessary to understand the behavior of individual cells that constitute the system, cell heterogeneity, and cell-to-cell interactions. With the rapid change of technology, a large number of single-cell sequencing technologies dominated by single-cell transcriptomes have emerged. Traditional sequencing (BulkRNA-seq) obtains the average transcriptome data of a group of cells and cannot study single cells. The currently developed single-cell sequencing (scRNA-seq) is a technology that performs sequencing and analysis of the genome, transcriptome, and epigenome at the level of a single cell. Single-cell multi-omics sequencing is used to simultaneously analyze the composition and characteristics of various biological macromolecules such as the genome, transcriptome, epigenome, and proteome of a single cell. The development of this technology enables researchers to gain a deep understanding of the function and diversity of individual cells. The data obtained through single-cell sequencing technology and the analysis of single-cell data can reveal the mechanism and regulatory network of biological processes. Single-cell data has high dimensionality, information redundancy, and contains a lot of noise, which makes data analysis difficult. Appropriate dimensionality reduction and feature selection of single-cell data can effectively remove noise and support biological downstream analysis such as visualization, clustering, and trajectory inference. Single-cell data analysis mainly involves data preprocessing (dimensionality reduction, data interpolation, batch correction, etc.) and analysis of biological downstream tasks (clustering, communication between cells, trajectory analysis, etc.).

[0003] Single-cell technology can study the genomics, transcriptomics, and multi-omics of single cells in tissues. Single-cell RNA sequencing (scRNA-seq) detects information such as gene expression at the single-cell level. The single-cell data obtained using single-cell sequencing technology cannot be directly applied to the analysis of downstream tasks, but these data need to be reduced in dimension, denoised, and filled with missing values. There are a large number of zero values ​​in high-dimensional single-cell data, which include biomass zero (i.e., the gene is not expressed in the cell) and technical errors (failed to detect due to low expression or other factors). A large number of missing values ​​in the data will affect biological signals. Single-cell data usually come from multiple experiments, with different capture times, processing personnel, reagent batches, equipment, and even technical platforms. These differences can lead to data anomalies or batch effects, which are likely to confuse cell specificity in the process of merging data. With the continued growth of scRNA-seq data, it is crucial to achieve effective data batch integration.

[0004] scRNA-seq is a high-throughput sequencing technology that analyzes thousands to millions of cells in a single experiment. It can identify heterogeneity in cell types, states, and populations, and provides gene expression information but does not provide spatial location information of cells. Transcriptome (ST) data is a low-throughput sequencing technology that is good at capturing spatial information of gene expression. It can understand the distribution of gene expression in tissue sections and obtain the average expression value of multiple cells, but it cannot obtain the expression value of a single cell. ST provides complementary information about gene expression and spatial location for scRNA-seq of dissociated cells. Since a large amount of scRNA-seq data has been generated, it is hoped that the cell location information of scRNA-seq can be restored by utilizing the relationship between gene expression and spatial location learned from ST. The restoration of location information of scRNA-seq data helps downstream task analysis, such as cell-to-cell communication.

[0005] The study of scRNA-seq data plays a vital role in many aspects. scRNA-Seq can help identify and describe different cell types in heterogeneous cell populations, recognize the heterogeneity of cancer cells and develop personalized treatments, understand the complexity of immune cell populations and their responses to stimuli, analyze the differentiation and regulation of stem cells and discover stem cell markers, and deeply understand the regulatory networks that control gene expression.

[0006] When performing downstream task analysis on scRNA-seq data, sufficient high-quality scRNA-seq data can help researchers draw relatively correct conclusions. However, factors such as the high cost of sequencing, low measurement accuracy of rare cell types, and data privacy will affect the quantity and quality of sequencing data, which is likely to affect the accuracy of subsequent analysis. To solve this problem, researchers can use different sequencing technologies to conduct multiple experiments to collect more samples. Doing so will not only increase the cost of the experiment, but the integration of data from different sequencing platforms will also introduce batch effects and lead to a decrease in data quality. On the other hand, due to privacy protection and other reasons, it is difficult for researchers to expand existing data. Summary of the invention

[0007] The present invention aims to solve the technical problems of high cost and difficulty in expanding existing data in scRNA-seq data sequencing experiments in the prior art, and further proposes a method for generating single-cell sequencing data;

[0008] The technical solution adopted by the present invention to solve the above problems is: the present invention proposes a method for generating single-cell sequencing data, comprising:

[0009] S1: Build a deep neural network model;

[0010] S2: Obtaining comprehensive scRNA-seq datasets;

[0011] S3, preprocessing the scRNA-seq comprehensive dataset;

[0012] S4: Training of deep neural network models based on preprocessed scRNA-seq comprehensive datasets;

[0013] S5: Based on the trained deep neural network model, the features of the test data are denoised, denoised and reconstructed to generate single-cell sequencing data.

[0014] Optionally, the deep neural network model in S1 consists of an encoder, a diffusion module, and a decoder;

[0015] The encoder and decoder form an autoencoder;

[0016] The encoder is used to extract gene expression data features of input data;

[0017] The diffusion module is used to add noise and denoise the extracted gene expression data features;

[0018] The decoder is used to reconstruct the denoised gene expression data features to generate single-cell sequencing data.

[0019] Optionally, the comprehensive scRNA-seq dataset in S2 includes scRNA-seq datasets under multiple sequencing technologies, scRNA-seq datasets with multiple attributes, and scRNA-seq datasets of cell growth and development at pseudo time nodes.

[0020] Optionally, the steps for obtaining scRNA-seq datasets using multiple sequencing technologies include:

[0021] S20101: The number of jump steps of the deep neural network model in the inference phase is set to 5, and compared with the existing four generative models: scDiffusion, cscGAN, LSH-GAN, and sciGAN;

[0022] S20102: Use SCC, PCC, Wasserstein, MMD, and ILISI metrics to evaluate the quality of single-cell gene expression data generated by these models and use UMAP to visualize single-cell gene expression data;

[0023] S20103: Use different sequencing technologies and combine single-cell gene expression data to obtain scRNA-seq datasets using multiple sequencing technologies.

[0024] Optionally, the steps for obtaining a scRNA-seq dataset with multiple attributes include:

[0025] S20201: Divide the single-cell data in the muris_mam_spl_T_B dataset into organ types and cell types. Organ types include breast and spleen, and cell types include B cells and T cells.

[0026] S20202: Construct training labels based on the divided organ type and cell type single-cell data. The training labels include breast T cells, breast B cells, spleen T cells, and spleen B cells.

[0027] S20203: Train the diffusion module based on the training labels and guide the diffusion module to generate scRNA-seq datasets with multiple attributes.

[0028] Optionally, the steps for obtaining a scRNA-seq dataset of cell growth and development at a pseudo-time node include:

[0029] S20301: Obtain Wadding-OT data, select the cell states of the first 15 time nodes in the Wadding-OT data, and use the time nodes as data labels;

[0030] S20302: Set the number of jump steps of the deep neural network model at the inference node to 5, use the cell state data and data labels of the first 15 time nodes in the Wadding-OT data to train the deep neural network model, and obtain the scRNA-seq data set of cell growth and development at the pseudo time nodes.

[0031] In the present invention, cfDiffusion can generate high-quality scRNA-seq data on data sets of various sequencing platforms, thereby ensuring the generation of higher-quality scRNA-seq data with multiple attributes, while greatly reducing the training cost and inference time of the model.

[0032] Optionally, the steps for preprocessing the scRNA-seq comprehensive dataset in S3 include:

[0033] S301: Convert scRNA-seq comprehensive dataset to gene expression matrix X n×m ;

[0034] S302: Gene expression matrix X n×m Perform normalization operation and transform the gene expression matrix X n×m The total technical scale of cells in is 10000. An offset of 1 is added and the logarithm is taken to obtain the normalized gene expression matrix S ori ;

[0035] The expression of the normalization operation is:

[0036]

[0037] In formula (1), the rows of X represent different cells, there are n cells in total, and the columns of X represent genes, there are m genes in total for each cell.

[0038] Optionally, the steps of training the deep neural network model in S4 include:

[0039] S401: The normalized gene expression matrix S ori Input encoder, the encoder uses 2 MLPs, S ori After the encoder, we get a 128-dimensional embedding layer x 0 ;

[0040] S402: x 0 The fully connected network in the input diffusion module performs continuous noise addition, and continuous denoising is performed through reverse diffusion, and the denoised gene expression matrix S is output. rec ;

[0041] S403: The denoised gene expression matrix S rec Input decoder, the decoder uses 3 MLPs, x 0 After the decoder, we get the same matrix S as the gene expression matrix ori The gene expression matrix S of the same dimension rec , and the gene expression matrix S rec Reconstruct and generate single-cell sequencing data;

[0042] Embedding layer x 0 The calculation formula is:

[0043] x 0 =Encoder(s ori ) (2);

[0044] In formula (2), Encoder is the encoder;

[0045] Gene expression matrix S rec The calculation formula is:

[0046] S rec =Decoder(x 0 ) (3);

[0047] In formula (3), Decoder is a decoder;

[0048] The denoised gene expression matrix S rec The reconstructed expression is:

[0049] Loss = MSE (S rec ,S ori ) (4).

[0050] Optionally, the diffusion module in S402 includes D 1 Layer, D 2 Layer, U 1 Layer, U 2 Layer, FG 1 Layer and FG 2 layer;

[0051] Output the denoised gene expression matrix S rec The steps include:

[0052] S40201: Get the embedding x at t time steps t , the time information t and label information y are integrated into the Chunk in the fully connected network, and each Chunk is used to identify x 0 Update and get the noised x 0 ;

[0053] The deep neural network model constructed by the present invention can simulate single-cell data under a pseudo-time scale, providing high-quality data support for tracking cell differentiation and development trajectories, analyzing cell-to-cell communication, revealing cell heterogeneity, and other analyses.

[0054] S40202: Train the diffusion module to obtain a reverse diffusion diffusion module, set the number of jumps to 3, and in the Tth time step, change x 0 x in t Input the trained fully connected network to get the denoised embedding x t-1 ;

[0055] S40203: In the T-1th time step, x t-1 After passing through the D1 layer, it directly passes through the U2 layer to calculate the denoised x t-2 ;

[0056] S40204: Repeat S40202-S40203 to update the noise sequence and jump process until the denoised x is obtained. 0 ;

[0057] The present invention adopts a jump calculation method to accelerate the reasoning speed of the deep neural network model, thereby increasing the degree of fusion of the data generated by the deep neural network model and the original data.

[0058] After denoising, x t-2 The calculation formula is:

[0059] ∈ t-2 =FC 2 {FC 1 [D 1 (x t-1 ,t-1,y)+U 2 (xsave ,t-1,y)]} (5);

[0060] x t-2 =Denoise(x t-1 ,∈ t-2 ) (6);

[0061] The expressions for updating the noise sequence and the jump process are:

[0062]

[0063] In formula (7), is the high-level feature numbered b-1 at time step t.

[0064] Optionally, the diffusion module for obtaining reverse diffusion in S40202 specifically includes:

[0065] The diffusion module is parameterized as ∈ θ , and add label information y during training to obtain the learned noise ∈ θ , and after adding noise x 0 Subtract the learned noise ∈ θ Afterwards, the diffusion module of reverse diffusion is obtained by removing the logarithm from the Bayesian formula and introducing the hyperparameter k;

[0066] After denoising, x t for:

[0067]

[0068] The Bayesian formula and the expressions for taking logarithms are:

[0069]

[0070] Introducing a hyperparameter k:

[0071]

[0072] In formula (10), k is used to adjust the fidelity and diversity of generated data, and its value range is (1,2).

[0073] The beneficial effects of the present invention are:

[0074] 1. The present invention proposes a cfDiffusion model. Based on the Diffusion model, by introducing a Classifier-Free method and a high-level feature caching mechanism, the training cost and inference time of the model are greatly reduced while ensuring the generation of higher quality scRNA-seq data with multiple attributes.

[0075] 2. The cfDiffusion model of the present invention can generate high-quality scRNA-seq data on data sets of multiple sequencing platforms. By comparing with other latest data generation models, the cfDiffusion model shows better performance in multiple evaluation indicators.

[0076] 3. The cfDiffusion model constructed in the present invention can also simulate single-cell data under a pseudo-time scale, providing high-quality data support for tracking cell differentiation and development trajectories, analyzing cell-to-cell communication, revealing cell heterogeneity, and other analyses. BRIEF DESCRIPTION OF THE DRAWINGS

[0077] Figure 1 A flow chart of a method for generating single-cell sequencing data provided by the present invention;

[0078] Figure 2 A structural diagram of the deep neural network model provided by the present invention;

[0079] Figure 3 Comparison diagram of cfDiffusion and scDiffusion simulation data quality assessment under multiple conditions provided by the present invention

[0080] Figure 4 A quality assessment diagram of cfDiffusion and scDiffusion simulation data under combined conditions provided by the present invention;

[0081] Figure 5 Distribution diagrams of Alles, Baron_Human, Marshall, and Mizrak data simulated by five UMAP visualization methods provided by the present invention;

[0082] Figure 6 The UMAP visualization method provided by the present invention simulates the distribution of WOT data and the random forest ROC curve diagram;

[0083] Figure 7 The UMAP visualization provided by the present invention is a distribution diagram of the real data and generated data of muris, pbmc68k, and human_lung in two-dimensional space. DETAILED DESCRIPTION

[0084] Specific implementation method 1: In this specific implementation method, the cfDiffusion model is a deep neural network model. Figure 1-Figure 7 To illustrate this embodiment, Figure 1 As shown, the method for generating single-cell sequencing data described in this embodiment includes the following steps:

[0085] S1. Obtain scRNA-seq dataset;

[0086] This implementation compares cfDiffusion (setting the jump step number of cfDiffusion to 5 in the inference stage) with the existing four generative models of scDiffusion, cscGAN, LSH-GAN, and sciGAN, and uses SCC, PCC, Wasserstein, MMD, and ILISI indicators to evaluate the quality of single-cell gene expression data generated by these models, and uses UMAP to visualize single-cell gene expression data. In this experiment, scRNA-seq datasets were selected, which include Alles, Baron_Human, Marshall, and Mizrak. Alles and Mizrak use Drop-seq sequencing technology, Baron uses InDrop sequencing technology, and Marshall uses 10X sequencing technology. The cell types of these four datasets are different. The reason for selecting datasets with different characteristics is to evaluate the robustness of each model in generating single-cell gene expression data in different scenarios.

[0087] The scRNA-seq data set obtained in this embodiment is as follows:

[0088] The single-cell data in the muris_mam_spl_T_B dataset has two attributes: organ type and cell type, and is composed of breast T cells, breast B cells, spleen T cells, and spleen B cells. When training cfDiffusion, this implementation uses organ type (breast, spleen) and cell type (B cells, T cells) as two labels to pair single-cell gene expression data, and trains Diffusion together. In the inference stage (the number of jump steps of cfDiffusion is set to 5), the organ type and cell type are combined into four labels: breast T cells, breast B cells, spleen T cells, and spleen B cells, and the label information is used to guide Diffusion to generate specific single-cell gene expression data. According to Figure 3 As shown in A, most of the generated data is overlaid on the real data, and some of the generated data is scattered around the real data. The model can learn the attribute characteristics of organs and cell types. This embodiment selects the marker gene Cd74 of breast B cells, and plots the expression distribution of gene Cd74 in real cells other than breast B cells, real breast B cells, and generated breast B cells into a box plot as shown in Figure 3As shown in .CD, the Wilcoxon rank sum test calculated the p-value of Cd74 gene expression in the breast B cell data simulated by cfDiffusion and the real breast B cell data to be 0.07618, indicating that there is no significant difference between the two distributions. However, the p-value of Cd74 gene expression in the breast B cell data simulated by scDiffusion and the real breast B cell data is 0.00017, indicating that there is a significant difference between the two distributions.

[0089] Wadding-OT includes the differentiation state of pluripotent stem cells (iPSCs) for 18 days, with a recording interval of 0.5 days. A total of 254,197 cells were sampled for differentiation state, and 19,423 genes were retained. Since the data set is too large, this embodiment selects the cell states of 15 time nodes, including Day 0, Day 0.5, Day 1, Day 1.5, Day 2, Day 2.5, Day 3, Day 4.5, Day 5, Day 5.5, Day 6, Day 6.5, Day 7, Day 7.5, Day 8, a total of 82,920 cells. This embodiment uses the time nodes as labels for these gene expression data, and uses the paired gene expression data and labels to train cfDiffusion. The number of jump steps of cfDiffusion is set to 5 in the inference stage. This embodiment generates 5,600 cell states for each time node, as shown in Tables 13 and 14, and uses various indicators to evaluate the data generated by cfDiffusion. In UMAP's dimensionality reduction and visualization of data, Figure 6 As shown, in this embodiment, it can be seen that the generated data and the real data distribution are roughly the same, and the SCC and PCC values ​​are both greater than 0.99, which indicates that the gene expression pattern of the generated single-cell data is extremely similar to the gene expression pattern of the real single-cell data. Figure 6 As shown, the value of ILISI is 0.8198, which represents the degree of mixing between the generated single-cell data and the real single-cell data, and the values ​​of the random forest in training, testing and AUC indicators show that the trained random forest model is not overfitting. Since the generated data and the real data are similar, it is largely impossible to distinguish whether the single-cell data is generated by the model. This embodiment uses real data and model-generated data at different time nodes to train the K nearest neighbor classifier, and the average accuracy of the prediction at each time node is about 0.5, which can also reflect that the quality of the single-cell data generated by cfDiffusion at each time node is very high.

[0090] Tabula Muris: RNA sequencing was performed on single cells from 20 tissues of 3-month-old mice on the Illumina NovaSeq 6000 platform, and single-cell data from 12 organs were selected, including Bladder, Heart_and_Aorta, Kidney, Limb_Muscle, Liver, Lung, Mammary_Gland, Marrow, Spleen, Thymus, Tongue, Trachea. There are a total of 57,004 cells, each with 18,996 genes.

[0091] PBMC68k: Peripheral blood mononuclear cells (PBMCs) are blood cells that are an important part of the immune system, used to fight infection and protect the body from harmful pathogens. In biomedical research, peripheral blood mononuclear cells are often used to study global immune responses to disease outbreaks and development, pathogen infection, vaccine development, and a variety of other clinical applications. This dataset contains single cell data for 11 cell types, including CD14+Monocyte, CD19+B, CD34+, CD4+THelper2, CD4+ / CD25 T Reg, CD4+ / CD45RA+ / CD25-NaiveT, CD4+ / CD45RO+Memory, CD56+NK, CD8+CytotoxicT, CD8+ / CD45RA+NaiveCytotoxic, and Dendritic, with a total of 68,579 cells, each with 17,789 genes.

[0092] Waddington-OT: A cell reprogramming dataset of mouse embryonic fibroblasts (MEFs) is used to study the effects of induced pluripotent stem cells (iPSCs) and growth differentiation factors (GDF9) on reprogramming. This dataset contains cells at different time stamps during the 18-day reprogramming process. Using data from three datasets, cells with expression counts less than 10 and genes expressed in less than 3 cells were filtered out, retaining 82,920 cells, each with 19,423 genes.

[0093] Homo Sapiens: measured on the Illumina NextSeq 500 platform, there are fibroblast, macrophage, memory bcell, There are six types of cells, nkcell, plasmcell, which come from Bladder, Blood, Spleen, Thymus, and Vasculature organs. The number of cells is 46177, and each cell has 28231 genes.

[0094] Human_PF_Lung: AT1, AT2, BCells, Basal, Ciliated, DifferentiatingCiliated, EndothelialCells, Fibroblasts, HAS1HighFibroblasts, KRT5- / KRT17+, Ly measured on human lung cells mphaticEndothelialCells,MUC5AC+High,MUC5B+,Macrophages,MastCells,MesothelialCells,Monocytes,Myofibroblasts,NKCells,PLIN2+Fibroblasts,PlasmaCells,ProliferatingEpithelialCel There are 31 types of cells, including ls, ProliferatingMacrophages, ProliferatingTCells, SCGB3A2+, SCGB3A2+SCGB1A1+, SmoothMuscleCells, TCells, TransitionalAT2, cDCs, and pDCs. The number of cells is 114,396, and each cell has 27,281 genes.

[0095] muris_mam_spl_T_B: Consists of mammalian mammary T cells, mammary B cells, splenic T cells, and splenic B cells. Cells with expression counts less than 10 and genes expressed in less than 3 cells are filtered out, thus retaining 11,330 cells, each with 14,652 genes.

[0096] Alles: The Alles dataset consists of 14 types of cells: LVM longitudinal visc.muscle, amnioserosa, developing midgut, epidermis, fat body, germ cells, head mesoderm / hemocyte differentiation, midgut, muscle, neurogenesis, neurons, undifferentiated cells, visceral muscle, yolk. Cells with expression counts less than 10 and genes expressed in less than 3 cells were filtered out, thus retaining 4614 cells, each with 14850 genes.

[0097] Baron Human: Using InDrop, a droplet-based single-cell RNA-seq method, we obtained the transcriptome of human pancreatic cells, including 14 types of cells: acinar, activated_stellate, alpha, beta, delta, ductal, endothelial, epsilon, gamma, macrophage, mast, quiescent_stellate, schwann, and t_cell. We filtered out cells with expression counts less than 10 and genes expressed in less than 3 cells, thus retaining 8569 cells, each with 16359 genes.

[0098] Marshall: The single-cell transcriptome data of Mus musculus liver were measured on the 10X platform. There were 16 types of cells, including B1, B2, B3, B4, B5, B6, B7, CD8 T, Cd4 T, Monocyte, Monocyte derived Macrophage, NK, Redpulp Macrophage, Treg, cDC, and pDC. Cells with expression counts less than 10 and genes expressed in less than 3 cells were filtered out, thus retaining 6022 cells, each with 16548 genes.

[0099] Mizrak: Drop-seq technology was used to analyze the V-SVZ in male and female mice aged 8-10 weeks. There were 17 types of cells, including Astrocyte, Astrocytes, COP, Doublet, Endothelial, Ependymal, Microglia, Mural Cell+Fibroblast, Mural Cell+Fibroblast, Mural+Fibroblast, Neuron, OPC, Oligodendrocyte, Oligodendrocytes, T Cell, aNSC+TAC+NB, and T cell, totaling 42,374 cells, each with 48,529 genes.

[0100] S2, build a deep neural network model, such as Figure 2 As shown, the cfDiffusion model framework constructed in this embodiment is used to generate single-cell gene expression data;

[0101] S201: The cfDiffusion model framework consists of an encoder, a diffusion model, and a decoder. The encoder extracts important features of gene expression data, and Diffusion continuously adds noise and denoises these features, and finally reconstructs them to the decoder to obtain the generated gene expression data.

[0102] S202: Diffusion Backbone network structure. This is designed based on the UNet model, where D 1 ,D 2 ,U 1 ,U 2 ,FC 1 ,FC 2 ,Mid modules all use MLP. This implementation integrates time information t and label information y into D 1 ,D 2 ,U 1 ,U 2 ,In the Mid module, the right half of the figure shows how the time information, label information and feature information are fused,x t is the feature information of a certain time step, and z is the feature information passed in by the previous module, such as module D 2 Get D 1 The incoming feature information, module Mid obtains D 2 The incoming feature information.

[0103] S203: When the number of jump steps is 3 and the module for specifying cache information is U 2 The model accelerates the process of generating single-cell gene expression data in the inference phase. First, U 2 The feature information is cached and passed to module U in the next t-1, t-2, and t-3 time steps. 1 Middle,U 1 After completing the feature information processing, 1 The feature information of the module is fused. In the whole process of processing feature information, only module D 1 ,U 1 ,FC 1 ,FC 2 without going through the gray module (D 2 ,Mid,U 2 ). The operation process of the t-4 time step is consistent with the t time step, which realizes the update of the cache information.

[0104] S3. Deep neural network model training.

[0105] S301: Train the autoencoder. n×m S is obtained by normalization (scaling the total count of each cell to 10,000, adding a small offset of 1, and finally taking the logarithm) ori As the input of VAE, the rows of X represent different cells, there are n cells in total, and the columns of X represent genes, and each cell has m genes in total.

[0106]

[0107] This operation not only helps to reduce the difference in sequencing depth between cells and the impact of extreme values ​​on the analysis, but also stabilizes the variance of the data and converts the distribution of the original data into an approximate normal distribution, which is consistent with the Gaussian distribution used in the diffusion process, making it easier for Diffusion to learn the reverse process. The VAE encoder uses 2 MLPs, S ori After the encoder, we get a 128-dimensional embedding layer (Embedding) x 0 :

[0108] x 0 =Encoder(s ori ) (2);

[0109] The decoder of VAE uses 3 MLPs, x 0 After the decoder, we get the same matrix S as the original gene expression matrix ori Same dimensions:

[0110] S rec =Decoder(x 0 ) (3);

[0111] The goal of optimization is to maximize S rec Refactoring to S ori :

[0112] Loss = MSE (S rec ,S ori ) (4);

[0113] S302: Training Classifier-Free Diffusion. In the forward process of Diffusion, the embedding x is trained in a series of random time steps. 0 Adding noise gives the embedding x at the tth time step t :

[0114]

[0115] In formula (5), ∈ fuses t-1 Gaussian distributions. When t is large enough, x tApproximately follows Gaussian distribution. Due to the high sparseness and high dimensionality of the gene expression matrix, the direct use of U-Net convolution operation is unlikely to capture the correlation between two genes that are far apart in the vector, or the captured correlation is weak. MLP can solve this problem. Therefore, using a fully connected network instead of a convolutional layer can make Diffusion generate single-cell data with better performance. The improved U-Net is as follows: Figure 2 B. The time information t and label information y are imported into each Chunk. The first Chunk is used as an example. The specific operation is shown in formula (6):

[0116] x t =x t ⊙y+t (6);

[0117] In formula (6), ⊙ is the dot product operation of the matrix. After a Chunk operation on x t In the improved U-Net, let the number of Chunks be n and n be an odd number, and The information of a Chunk is directly connected to the low-level information via a jump connection. Chunks, as shown in formula (7):

[0118] U i+1 =U i (·)+D i (·) (7);

[0119] In formula (7), U i After Chunk, D i It is before Chunk, A Chunk is equivalent to the middle layer of the U-Net network. t The whole process of , t, and y passing through the entire U-Net is shown in the pseudo code Algorithm 1.

[0120]

[0121]

[0122] S5: Train Diffusion so that it learns from x in the reverse process t Noise reduction x t-1 , the process can be expressed as: p(x t-1 |x t ). The Diffusion model is parameterized as ∈ θ , and add label information y during the training process, assuming ∈θ (x t ,t,y) obeys the standard normal distribution, according to formula (8):

[0123]

[0124] In fact, the noise ∈ θ It is approximately subject to normal distribution. According to TweedieEstimator, we can further obtain formula (9):

[0125]

[0126] Then, using the Bayesian formula and taking the logarithm, we get:

[0127]

[0128] Introducing a hyperparameter k:

[0129]

[0130] In formula (11), k is used to adjust the fidelity and diversity of the generated data, and its value range is (1,2).

[0131] S501: Diffusion reasoning acceleration;

[0132] S50101: In the reverse diffusion process, x i Its neighboring feature x i-1 、x i-2 Similarly, in consecutive time steps, redundant calculations are eliminated by skipping some U-Net network layers. To better understand the entire accelerated reasoning process, this implementation sets the number of jumps to 3, such as Figure 2 As shown in C, each jump to the branch chunk 5 That is, upsampling U 1 Up,U i The features above are called high-level features, D i The features above are called low-level features, and a high-level feature will be saved for skip calculation. Figure 3 As shown in C, in this implementation, it can be observed that U 1 The features of the layer are saved as x save , x t After passing through the entire U-Net and denoising, we get the embedding x t-1 . In the T-1th time step, x t-1 Did not go through D 2 , M, U 1 layer, but through D 1 After calculation, it goes directly through U2 , the calculated ∈ t-2 As shown in formula (12):

[0133] ∈ t-2 =FC 2 {FC 1 [D 1 (x t-1 ,t-1,y)+U 2 (x save ,t-1,y)]} (12);

[0134] Finally, further denoising is obtained to obtain x t-2 :

[0135] x t-2 =Denoise(x t-1 ,∈ t-2 ) (13);

[0136] S50102: In the T-2 and T-3 time steps, this embodiment can obtain x according to the jump calculation method. t-3 、x t-4 , so far it has jumped 3 times. In the T-4th time step, x t-4 Need to go through the entire U-Net to update x save , and then use the same jump calculation method to get x t-5 、x t-6 、x t-7 And so on, until the implementation method obtains x 0 .

[0137] S502: Figure 7 As shown in Figure 1, in the inference phase of the model, it is slow to generate a specific type of single-cell gene expression data using no classifier. We use jump calculation to speed up the inference speed. Since the jump interval is too small, the comparison of experimental results is not very significant, so the jump interval is set to 5. Figure 7The figure shows the change in the quality of the generated single-cell gene expression data as the number of jump intervals increases. The most obvious trend is that the time overhead is gradually decreasing, as shown in Table 2. As shown in Tables 1 to 4, when the number of jump intervals increases from 5 to 20, the quality of the generated single-cell expression data decreases slightly, which shows that the features of adjacent time steps are similar and there is a certain redundancy. As shown in Table 5, when the number of jump intervals is 50, the quality loss of the generated single-cell expression data is more obvious, which shows that some important features are deleted. Although the data quality decreases within the jump interval of 5 to 20 steps, the quality of the data generated by cfDiffusion is not worse than that of scDiffusion from the overall reflection of the indicators. For example, in the experiment of the WOT dataset, it is more difficult to distinguish the data generated by cfDiffusion than the random forest of scDIffusion. The increase in the ILISI index also shows that the fusion degree of the data generated by cfDiffusion and the original data has increased.

[0138] Table 1

[0139]

[0140] Table 2

[0141]

[0142] Table 3

[0143]

[0144] Table 4

[0145]

[0146] Table 5

[0147]

[0148] like Figure 7 As shown, UMAP visualizes the distribution of real data and generated data of muris, pbmc68k, and human_lung in two-dimensional space, and the quality changes of the generated data when the model jumps to 5, 10, 20, and 50 steps. (A) The quality changes of muris generated by cfDiffusion when the jump steps are 5, 10, 20, and 50. (B) The quality changes of pbmc68k generated by cfDiffusion when the jump steps are 5, 10, 20, and 50. (C) The quality changes of human_lung generated by cfDiffusion when the jump steps are 5, 10, 20, and 50.

[0149] Now generalize all the conditions, first define a time step to update x save The sequence ξ={p 1 ,p 2 ,…,p m}, where p i ∈[0,999], m≤1000. If the jump step size is fixed, then |p i -p i+1 |=c, c is a constant; if the jump step size is dynamically changed, then P(|p i -p i+1 |=|p i+1 -p i+2 |)>0, where P(·) represents the probability of something happening. Specify a jump to any U b , b represents the number. Since the structure of U-Net is symmetrical, U b The corresponding chunk is D b . Update x save The jumping process is shown in formula (14):

[0150]

[0151] in represents the high-level feature numbered b-1 at time step t.

[0152] Example

[0153] As described in S1 of the specific implementation, this embodiment compares cfDiffusion (setting the jump step number of cfDiffusion to 5 in the inference stage) with the existing four generative models of scDiffusion, cscGAN, LSH-GAN, and sciGAN, and uses SCC, PCC, Wasserstein, MMD, and ILISI indicators to evaluate the quality of single-cell gene expression data generated by these models, and uses UMAP to visualize single-cell gene expression data. In this experiment, Alles, Baron_Human, Marshall, and Mizrak data sets were selected, among which Alles and Mizrak used Drop-seq sequencing technology, Baron used InDrop sequencing technology, and Marshall used 10X sequencing technology, and the cell types of these four data sets were different. The reason for selecting data sets with different characteristics is to evaluate the robustness of each model in generating single-cell gene expression data in different scenarios. Tables 6, 7, 8, and 9 respectively show the performance of the five models on the quality evaluation indicators of Alles data, Baron_Human data, Marshall data, and Mizrak data. The method cfDiffusion of this embodiment performs best on the Baron_Human, Marshall, and Mizrak data sets, and has the highest ILISI score on the Alles data set. The scores of the other indicators are very close to those of scDiffusion and can be almost ignored. In particular, the MMD and ILISI indicators show that the distribution difference between the single-cell gene expression data generated by cfDiffusion and the real single-cell gene expression data is small, and the degree of fusion is high. At the same time, this embodiment uses UMAP to visualize the highly variable genes in the single-cell gene expression data, such as Figure 5 As shown in Figure 3, the distribution of single-cell gene expression data generated by cfDiffusion and scDiffusion are close to the distribution of true single-cell gene expression data, while the distribution of single-cell data generated by other methods is far from the distribution of true single-cell gene expression data.

[0154] Table 6

[0155]

[0156] Table 7

[0157]

[0158] Table 8

[0159]

[0160] Table 9

[0161]

[0162] As described in S2 of the specific implementation manner, the single-cell data in the muris_mam_spl_T_B data set has two attributes, organ type and cell type, and is composed of breast T cells, breast B cells, spleen T cells, and spleen B cells. When training cfDiffusion, this embodiment uses organ type (breast, spleen) and cell type (B cells, T cells) as two labels for pairing single-cell gene expression data, and trains Diffusion together. In the inference stage (the number of jump steps of cfDiffusion is set to 5), the organ type and cell type are combined into four labels: breast T cells, breast B cells, spleen T cells, and spleen B cells, and the label information is used to guide Diffusion to generate specific single-cell gene expression data. According to Figure 3 As shown in A, most of the generated data is overlaid on the real data, and some generated data is scattered around the real data. The model can learn the attribute characteristics of organs and cell types. This example selects the marker gene Cd74 of breast B cells, and plots the expression distribution of gene Cd74 in real cells other than breast B cells, real breast B cells, and generated breast B cells into a box plot, as shown in Figure 3 As shown in C and D. The Wilcoxon rank sum test calculated that the p-value of the gene expression of Cd74 in the breast B cell data simulated by cfDiffusion and the real breast B cell data was 0.07618, indicating that there was no significant difference between the two distributions. However, the p-value of the gene expression of Cd74 in the breast B cell data simulated by scDiffusion and the real breast B cell data was 0.00017, indicating that there was a significant difference between the two distributions. When scDiffusion generates data on breast T cells, breast B cells, spleen T cells, and spleen B cells, it is necessary to train classifiers for two attributes, one for guiding the model to generate data for specific organs, and the other for guiding the model to generate data for specific cell types. cfDiffusion does not need to train any classifiers, which greatly saves training costs. Under multiple conditions, according to Figure 3 As shown in B and Table 10, Table 10 is a comparison of the quality of muris_T_B data simulated by cfDiffusion and scDiffusion. The distribution of single-cell data generated by scDiffusion is significantly different from that of the real data, and scDiffusion does not perform as well as cfDiffusion.

[0163] Table 10

[0164]

[0165] The single-cell data in the Sapiens dataset also has two attributes. Unlike the muris_mam_spl_T_B dataset, the organ types and cell types in Sapiens are not always in pairs. Therefore, multiple conditions are combined, such as combining the organ type Bladder and the cell type fibroblast into one condition `fibroblast_Bladder`. In this example, cfDiffusion and scDiffusion are trained separately in this way, and the simulated data are visualized by observing Figure 4 .A It can be seen that the distribution of the data simulated by the two methods is roughly similar to that of the real data, such as Figure 4 As shown in B, the trend of the five colored lines in the QQ-plot is roughly close to the reference line (dashed line), indicating that the expression distribution of the gene AKIRIN2 in the five organs (Bladder, Blood, Spleen, Thymus, Vasculate) in the data simulated by the two methods is similar to that in the real data. The data quality simulated by cfDiffusion and scDiffusion is shown in Tables 11 and 12. Table 11 is a comparison of the quality of Homo Sapiens data simulated by cfDiffusion and scDiffusion, and Table 12 is a measure of the authenticity of the two methods of simulating Homo Sapiens data using random forests. Under the evaluation of various indicators, cfDiffusion is almost completely superior to scDiffusion.

[0166] This experiment shows that scDiffusion is only suitable for simulating single-cell gene expression data under a single condition, and the data quality simulated under multiple conditions is far inferior to cfDiffusion, which is most likely caused by improper guidance of multiple classifiers. This experiment also shows that cfDiffusion has good performance under any conditions.

[0167] Table 11

[0168]

[0169] Table 12

[0170]

[0171] As described in S3 of the specific implementation method, Wadding-OT includes pluripotent stem cells (iPSCs) undergoing 18 days of differentiation state, with a recording interval of 0.5 days. A total of 254,197 cells were sampled for differentiation state, and 19,423 genes were retained. Since the data set is too large, this embodiment selects the cell states of 15 time nodes, including Day 0, Day 0.5, Day 1, Day 1.5, Day 2, Day 2.5, Day 3, Day 4.5, Day 5, Day 5.5, Day 6, Day 6.5, Day 7, Day 7.5, Day 8, a total of 82,920 cells. In this embodiment, the time nodes are used as labels of these gene expression data, and the paired gene expression data and labels are used to train cfDiffusion. The number of jump steps of cfDiffusion is set to 5 in the inference stage, as shown in Tables 13 and 14. Table 13 is a quality comparison of WOT data simulated under pseudo time scale by two methods, and Table 14 is a comparison of the authenticity of WOT data simulated under pseudo time scale by two methods using random forest. In this embodiment, 5600 cell states are generated for each time node, and various indicators are used to evaluate the data generated by cfDiffusion. Figure 6 As shown in the dimensionality reduction and visualization of data by UMAP, it can be seen that the generated data and the real data distribution are roughly the same, and the values ​​of SCC and PCC are both greater than 0.99, which indicates that the gene expression pattern of the generated single-cell data is extremely similar to the gene expression pattern of the real single-cell data. The value of ILISI is 0.8198, which represents the degree of mixing between the generated single-cell data and the real single-cell data. Figure 6 As shown, the values ​​of random forest in training, testing and AUC indicators show that the trained random forest model is not overfitting. Since the generated data and the real data are similar, it is largely impossible to distinguish whether the single-cell data is generated by the model. In this embodiment, the real data and model-generated data at different time nodes are used to train the K nearest neighbor classifier. The average accuracy of the prediction at each time node is about 0.5, which can also reflect that the quality of the single-cell data generated by cfDiffusion at each time node is very high.

[0172] In terms of various indicators, cfDiffusion and scDiffusion are obviously better than scDiffusion in MMD, ILIST, and random forest evaluation results, and the other indicators are numerically similar or equivalent. To a certain extent, it can be shown that cfDiffusion performs better than scDiffusion.

[0173] Table 13

[0174] Method cfDiffusion scDiffusion SCC↑ 0.98569 0.99261 PCC↑ 0.9933 0.9997 Wasserstein↓ 0.01637 0.00457 MMD↓ 2.6222 2.7052 ILISI↑ 0.81984 0.79245 KNN_AUC↓ 0.51 0.51 KNN_ACC↓ 0.49969 0.49969

[0175] Table 14

[0176] Method AUC↓ ACC_Train↓ ACC_Test↓ OOB_Score↓ cfDiffusion 0.679 0.635689018 0.628147612 0.629048078 scDiffusion 0.7444 0.679112397 0.685021708 0.67895964

[0177] The above is only a preferred embodiment of the present invention and does not limit the present invention in any form. Although the present invention has been disclosed as a preferred embodiment as above, it is not used to limit the present invention. Any technician familiar with this profession can make some changes or modify the technical contents disclosed above into equivalent embodiments without departing from the scope of the technical solution of the present invention. However, any simple modification, equivalent replacement and improvement made to the above embodiments without departing from the content of the technical solution of the present invention, based on the technical essence of the present invention, within the spirit and principles of the present invention, still fall within the protection scope of the technical solution of the present invention.

Claims

1. A method for generating single-cell sequencing data, characterized in that: It includes the following steps: S1: Build a deep neural network model; S2: Obtaining comprehensive scRNA-seq datasets; S3, preprocessing the scRNA-seq comprehensive dataset; S4: Training of deep neural network models based on preprocessed scRNA-seq comprehensive datasets; S5: Based on the trained deep neural network model, the features of the test data are denoised, denoised and reconstructed to generate single-cell sequencing data.

2. The method for generating single-cell sequencing data according to claim 1, characterized in that: The deep neural network model in S1 consists of an encoder, a diffusion module, and a decoder; The encoder and decoder constitute an autoencoder; The encoder is used to extract gene expression data features of input data; The diffusion module is used to add noise and denoise the extracted gene expression data features; The decoder is used to reconstruct the features of the denoised gene expression data to generate single-cell sequencing data.

3. The method for generating single-cell sequencing data according to claim 1, characterized in that: The comprehensive scRNA-seq datasets in S2 include scRNA-seq datasets under various sequencing technologies, scRNA-seq datasets with multiple attributes, and scRNA-seq datasets of cell growth and development at pseudo-time nodes.

4. The method for generating single-cell sequencing data according to claim 3, characterized in that: The steps of obtaining the scRNA-seq data set under the multiple sequencing technologies include: S20101: setting the number of jump steps of the deep neural network model in the inference stage to 5, and comparing it with the existing four generative models of scDiffusion, cscGAN, LSH-GAN, and sciGAN; S20102: Use SCC, PCC, Wasserstein, MMD, and ILISI metrics to evaluate the quality of single-cell gene expression data generated by these models and use UMAP to visualize single-cell gene expression data; S20103: Use different sequencing technologies and combine single-cell gene expression data to obtain scRNA-seq datasets using multiple sequencing technologies.

5. The method for generating single-cell sequencing data according to claim 3, characterized in that: The steps of acquiring the scRNA-seq data set with multiple attributes include: S20201: Divide the single-cell data in the muris_mam_spl_T_B dataset into organ types and cell types, where the organ types include mammary gland and spleen, and the cell types include B cells and T cells; S20202: constructing training labels based on the divided organ type and cell type single-cell data, wherein the training labels include breast T cells, breast B cells, spleen T cells, and spleen B cells; S20203: Training a diffusion module based on the training labels, and guiding the diffusion module to generate scRNA-seq datasets with multiple attributes.

6. The method for generating single-cell sequencing data according to claim 3, characterized in that: The steps of acquiring the scRNA-seq data set of cell growth and development at the pseudo time node include: S20301: Obtain Wadding-OT data, select the cell states of the first 15 time nodes in the Wadding-OT data, and use the time nodes as data labels; S20302: The number of jump steps of the deep neural network model at the inference node is set to 5, and the deep neural network model is trained using the cell state data and data labels of the first 15 time nodes in the Wadding-OT data to obtain a scRNA-seq data set of cell growth and development at a pseudo time node.

7. The method for generating single-cell sequencing data according to claim 1, characterized in that: The steps for preprocessing the scRNA-seq comprehensive dataset in S3 include: S301: Convert the scRNA-seq comprehensive dataset into a gene expression matrix X n×m ; S302: Gene expression matrix X n×m Perform normalization operation and transform the gene expression matrix X n×m The total technical scale of cells in is 10000. An offset of 1 is added and the logarithm is taken to obtain the normalized gene expression matrix S ori ; The expression of the normalization operation is: In formula (1), the rows of X represent different cells, there are n cells in total, and the columns of X represent genes, there are m genes in total for each cell.

8. The method for generating single-cell sequencing data according to claim 1, characterized in that: The steps for training a deep neural network model in S4 include: S401: The normalized gene expression matrix S ori Input encoder, the encoder uses 2 MLPs, S ori After the encoder, a 128-dimensional embedding layer x0 is obtained; S402: Input x0 into the fully connected network in the diffusion module for continuous noise addition, perform continuous denoising through reverse diffusion, and output the denoised gene expression matrix S rec ; S403: The denoised gene expression matrix S rec Input decoder, the decoder uses 3 MLPs, x0 passes through the decoder to get the same as the gene expression matrix S ori The gene expression matrix S of the same dimension rec , and the gene expression matrix S rec Reconstruct and generate single-cell sequencing data; The calculation formula of the embedding layer x0 is: x0=Encoder(s ori ) (2); In formula (2), Encoder is the encoder; Gene expression matrix S rec The calculation formula is: S rec =Decoder(x0) (3); In formula (3), Decoder is a decoder; The denoised gene expression matrix S rec The reconstructed expression is: Loss=MSE(S rec ,S ori ) (4)。 9. The method for generating single-cell sequencing data according to claim 8, characterized in that: The diffusion module in S402 includes D1 layer, D2 layer, U1 layer, U2 layer, FG1 layer and FG2 layer; Output the denoised gene expression matrix S rec The steps include: S40201: Get the embedding x at t time steps t , the time information t and label information y are integrated into the Chunk in the fully connected network, and x0 is updated based on each Chunk to obtain the noisy x0; S40202: Train the diffusion module to obtain a reverse diffusion diffusion module, set the number of jumps to 3, and in the Tth time step, replace x in x0 t Input the trained fully connected network to get the denoised embedding x t-1 ; S40203: In the T-1th time step, x t-1 After passing through the D1 layer, it directly passes through the U2 layer to calculate the denoised x t-2 ; S40204: Repeat S40202-S40203 to iteratively update the noise sequence and the jump process until the denoised x0 is obtained; After denoising, x t-2 The calculation formula is: ∈ t-2 =FC2{FC1[D1(x t-1 ,t-1,y)+U2(x save ,t-1,y)]} (5); x t-2 =Denoise(x t-1 ,∈ t-2 ) (6); The expressions for updating the noise sequence and the jump process are: In formula (7), is the high-level feature numbered b-1 at time step t.

10. A method for generating single-cell sequencing data according to claim 9, characterized in that: The diffusion module for obtaining reverse diffusion in S40202 specifically includes: The diffusion module is parameterized as ∈ θ , and add label information y during training to obtain the learned noise ∈ θ , and subtract the learned noise ∈ from the noised x0 θ Afterwards, the diffusion module of reverse diffusion is obtained by taking the logarithm of the Bayesian formula and introducing the hyperparameter k; After denoising, x t for: The Bayesian formula and the expressions for taking logarithms are: Introducing a hyperparameter k: In formula (10), k is used to adjust the fidelity and diversity of the generated data, and its value range is (1,2).

Images