Single cell time sequence data simulation and completion method and system

Abstract

The invention belongs to the technical field of single cell data processing, and particularly relates to a single cell time sequence data simulation and completion method and system. The method comprises the following steps that interpolation is carried out on hidden space vector representation in a time starting point and end point single cell matrix variational auto-encoder (VAE) model, then interpolation is carried out on the two vectors in a hidden space, and a middle time point hidden space vector is generated; finally, the vectors are reduced into a gene expression matrix through a decoder, so that missing or unobserved single cell data at the middle moment in the time sequence are simulated, the single cell data in the time sequence can be simulated through the method, and therefore the continuous change of the cell state can be deduced and complemented; the method is suitable for single cell level pedigree reconstruction and trajectory analysis in dynamic biological processes such as development, aging and disease progression, and has high biological interpretation value and practical application significance.

Description

Technical Field

[0001] This invention belongs to the field of single-cell data processing technology, and particularly relates to a method and system for simulating and completing single-cell time series data. Background Technology

[0002] Single-cell RNA sequencing (SCRNA) technology enables quantitative analysis of RNA within cells at single-cell resolution. It reveals cell population heterogeneity and elucidates gene regulatory networks at the RNA level. Because it can essentially capture the entire gene expression profile of a single cell, SCRNA data is highly interpretable, offering unique advantages for exploring key gene events and identifying rare cell types. It greatly advances the understanding of complex events in life development and is an important tool in life science and medical research. Commonly used single-cell sequencing technologies include 10× Genomics, SMART-seq2, and MERFISH. These technologies are based on different principles (such as microfluidic droplets, full-length cDNA amplification, and multiplex fluorescence in situ hybridization) to perform single-cell-level gene expression analysis on tissue samples. Among them, 10× Genomics technology, due to its high-throughput characteristics, can efficiently map large-scale cell atlases and has become the most widely used platform.

[0003] However, current single-cell sequencing technology still has certain limitations. A prominent problem is that mainstream methods are usually endpoint detection methods: to obtain a single-cell suspension, tissue dissociation is required, and the sample cannot be recovered after sequencing; or the sample needs to be fixed and stained, leading to loss of cell viability. This means that to analyze continuous dynamic changes in cell state (such as development or disease progression), multiple independent sequencing operations must be performed at different time points. This not only significantly increases experimental costs but also introduces significant batch effects and noise due to technical variations between samples, posing challenges to data integration and causal inference. Summary of the Invention

[0004] This invention addresses the problem of multiple sequencing attempts or sample loss during continuous dynamic change studies using single-cell sequencing technology. It provides a method and system for simulating and completing single-cell time-series data. In studies spanning time and physiological states, this method uses single-cell transcriptomes collected at two different time points or physiological states as references to simulate gene expression profiles at intermediate time points or transitional physiological states, thereby obtaining continuous single-cell expression profiles. This enables the simulation and completion of single-cell data during time-series studies and provides good biological interpretability.

[0005] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0006] This invention provides a method for simulating and completing single-cell time series data, comprising the following steps:

[0007] (1) Obtain single-cell sequencing data at the start and end time points or physiological states, use cell annotation tools to annotate the data by cell type, obtain two corresponding annotation data objects, integrate the two annotation data objects and perform sequencing depth standardization to obtain a standardized expression matrix, and use it as an input matrix to train and optimize the variational autoencoder model.

[0008] (2) Input the standardized expression matrix into the optimized variational autoencoder model, reduce the dimension and denoise the biological information in the expression matrix to obtain a continuous single-cell level latent space vector representation. Based on the distance minimization criterion, establish the correspondence between cells at the start and end time points in the latent space. For each pair of corresponding cells, use the interpolation method to calculate the latent space vector of the intermediate time point cell. Decode the latent space vector of the intermediate time point through the decoder of the variational autoencoder model to obtain the simulated intermediate time point expression matrix.

[0009] (3) Based on the single-cell sequencing data of the start and end time points or physiological states, the expression level of the obtained intermediate time point expression matrix is ​​corrected, and finally the simulated intermediate state single-cell expression matrix is ​​output.

[0010] Further, in step (1), the cell annotation tool is the Scanpy tool; the annotation data object is an AnnData object, which includes cell annotation information stored in the AnnData observation structure and gene information stored in the AnnData feature structure.

[0011] Furthermore, the cell annotation information includes the sample source, collection time point, and physiological state, and the gene information includes the gene name and the total number of gene reads. Further, in step (1), the sequencing depth normalization specifically involves normalizing the total number of counts for each cell to 10,000.

[0012] Further, in step (1), the variational autoencoder model is trained and optimized using a composite loss function, which is composed of a weighted sum of a mean squared error loss function and a KL divergence loss function, wherein the weight coefficient of the KL divergence loss function is 0.5.

[0013] Further, in step (2), the minimization criterion is: using the minimum Euclidean distance between latent space vectors as the criterion, a mapping relationship between cells from the endpoint to the starting time point is established.

[0014] Further, in step (2), the interpolation method specifically involves: finding the midpoint for each pair of corresponding cells to obtain the single-cell level latent space vector representation of the intermediate state, and decoding it to obtain the standardized expression matrix of the sequencing depth at the intermediate time point. Further, step (3) specifically involves: evaluating the gene dropout rate of single-cell sequencing data at the start and end time points or physiological states; based on the evaluation results, scaling the gene expression matrix of the simulated intermediate time point matrix and setting low-expressed genes to zero to correct the expression level, thereby obtaining the simulated single-cell expression matrix of the intermediate state.

[0015] This invention also provides a single-cell time series data simulation and completion system, the system comprising:

[0016] The data preprocessing and training module acquires and integrates single-cell sequencing data from multiple time points, performs cell annotation and expression matrix standardization, and uses the standardized data to train a variational autoencoder model.

[0017] The latent space simulation module, connected to the data preprocessing and training module, encodes the standardized expression matrix into the latent space, establishes cell mapping relationships, generates latent space vectors at intermediate time points through interpolation, and decodes them into standardized expression matrices.

[0018] The expression correction and output module is connected to the latent space simulation module. It performs correction processing on the decoded simulated expression matrix based on the statistical features of real data and outputs the final simulated single-cell gene expression matrix.

[0019] Compared with the prior art, the present invention has the following beneficial effects:

[0020] 1. For time-series single-cell data, it can combine the start and end point sequencing matrices to perform transition state simulation in high-dimensional space and restore it to the expression matrix, which has biological interpretability;

[0021] 2. It generates data quickly and robustly, and also performs efficiently in large-scale simulation generation.

[0022] 3. Simulating biological processes can, to some extent, solve the problem that single-cell sequencing technology cannot achieve continuous observation and can be applied to time-series research. Attached Figure Description

[0023] Figure 1 The flowchart illustrates the single-cell time series data simulation and completion method provided by this invention.

[0024] Figure 2 The data is from the Zygote-CS16 developmental lineage of human embryos, with the CS8 time point representing single-cell data generated by this method through simulation and completion.

[0025] Figure 3 The data presented here are real data from the CS8 stage of human gastrulation and data simulated and completed using the methods described in this invention, along with scDesign3, scVAEDer, Splatter, and Splatter-simple methods. Cell type annotation results are shown in the bubble chart, which displays the marker genes for different cell types during annotation. The UMAP diagram shows the cell clustering results. Figure 3 In this context, A represents the annotation result of the actual cell type data. Figure 3 In this paper, B represents the cell type annotation result of the method described in this invention. Figure 3 In this context, C represents the cell type annotation result from the scDesign3 method. Figure 3 In this context, D represents the cell type annotation result from the scVAEDer method. Figure 3 In this context, E represents the cell type annotation result from the Splatter method. Figure 3 F in the figure represents the annotation result of the Splatter-simple method.

[0026] Figure 4 The data includes real data from the E14 stage of pre-gastrulation in humans, and data simulated and completed using the methods described in this invention, along with scDesign3, scVAEDer, Splatter, and Splatter-simple methods. Cell type annotation results are shown in the bubble chart (marker genes for different cell types during annotation) and the UMAP diagram (cell clustering results). Figure 4 In this context, A represents the annotation result of the actual cell type data. Figure 4 In this paper, B represents the cell type annotation result of the method described in this invention. Figure 4 In this context, C represents the cell type annotation result from the scDesign3 method. Figure 4 In this context, D represents the cell type annotation result from the scVAEDer method. Figure 4 In this context, E represents the cell type annotation result from the Splatter method. Figure 4 F in the figure represents the annotation result of the Splatter-simple method.

[0027] Figure 5 The data presented here are real data from the CS10 stage of organogenesis and data simulated and completed using the methods described in this invention, along with scDesign3, scVAEDer, Splatter, and Splatter-simple methods. The bubble chart shows the marker genes for different cell types during annotation, and the UMAP diagram shows the cell clustering results. Figure 5 In this context, A represents the annotation result of the actual cell type data. Figure 5 In this paper, B represents the cell type annotation result of the method described in this invention. Figure 5 In this context, C represents the cell type annotation result from the scDesign3 method. Figure 5 In this context, D represents the cell type annotation result from the scVAEDer method. Figure 5 In this context, E represents the cell type annotation result from the Splatter method. Figure 5 F in the figure represents the annotation result of the Splatter-simple method.

[0028] Figure 6 The figure shows a comparison of the performance and biological information restoration accuracy of this method with scDesign3, scVAEDer, Splatter, and Splatter-simple simulation and completion methods. Figure 6 In the graph, A represents the runtime test results. Figure 6 B in the diagram represents the method. Figure 6 C in the figure represents the correlation results between cell populations assessed using the Silouette profile coefficient; Figure 6 In the figure, D represents the results of assessing the correlation between cell populations using the kBET parameter; Figure 6 E in the figure is a schematic diagram illustrating the principle of evaluating the spectral reduction degree of different matrix generation methods; Figure 6 F in the figure represents the result of using Pearson_start to evaluate the relationship between the "starting point" column of the moment and the ground truth. Figure 6 G in the figure represents the result of using the Pearson_simulated evaluation matrix to assess the relationship between the "transition state" column and the ground truth. Figure 6 H in the figure represents the result of using Pearson_end to evaluate the relationship between the "endpoint" column of the matrix and the ground truth. Figure 6 The "I" in the graph represents the result of using Pearson_relative to evaluate the relationship between the "relative distance" column of the matrix and the ground truth. Detailed Implementation

[0029] To make the technical problems, technical solutions and advantages of the present invention clearer, a detailed description will be given below in conjunction with the accompanying drawings and specific embodiments.

[0030] This invention provides a method for simulating and completing single-cell time series data, the flowchart of which is shown below. Figure 1 As shown, it includes the following steps:

[0031] Step 1: Data Preprocessing

[0032] (1.1) Obtain single-cell sequencing data at two time points or physiological states, defined as follows: and The Scanpy tool was used to annotate the cell types in the sequencing data, resulting in two constructed AnnData objects, adata1 and adata2. Cell annotation information, such as sample source, collection time point, and physiological state, is stored in the AnnData Observations (Obs) structure, while gene information, such as gene name, total number of reads, and whether it is hypervariable, is stored in the AnnData Variables (Var) structure.

[0033] Where adata1 is X start The annotation results contain the expression of X genes in M ​​cells, and adata2 is... The annotation results include the expression of Y genes in N cells.

[0034] (1.2) Integrate the adata1 and adata2 obtained from the annotations to obtain the total In cells Gene expression: The integrated expression matrix contains expression information of all cells obtained from single-cell sequencing data at two time points. The Obs structure is the union of two Obs structures, and the Var structure is the union of two Var structures.

[0035] (1.3) The sequencing depth of the integrated matrix is ​​normalized so that the total number of counts for each cell is 10,000, and this matrix is ​​used as the input matrix to train the VAE model.

[0036] Step 2: VAE Model Training and Output Simulation

[0037] A variational autoencoder (VAE) is used as the generative model framework to learn the probability distribution of single-cell expression data in the latent space of a deeply normalized matrix.

[0038] (2.1) Use a VAE model to learn the distribution of the input matrix (mean μ and variance σ) 2 );

[0039] The model's loss function (L) consists of the mean squared error loss function (MSE) and the KL divergence loss function. MSE is defined as follows:

[0040]

[0041] In the formula, n is the dimension of the input data, x i This represents the original input data, x' i This represents the output reconstructed by the VAE decoder after sampling from the latent space.

[0042] The KL divergence loss is defined as:

[0043]

[0044] In the formula For the prior distribution, Let m be the posterior distribution of the model output, and m be the dimension of the latent space. , Let and represent the mean and variance of the j-th dimension of the latent space, respectively.

[0045] The formula for calculating the loss function (L) is as follows:

[0046] The ratio between the mean squared error loss function and the KL divergence loss function is adjusted by the parameter β (default β = 0.5).

[0047] (2.2) The integrated standardized expression matrix is ​​input into the VAE model, and the biological information in the expression matrix is ​​reduced in dimension and denoised to obtain continuous single-cell level latent space vector representations, which are respectively and The dimension is the number of cells × the dimension of the latent space. The latent space vector representation Z of cells at intermediate time points is calculated using interpolation. mid Then, decoding is performed to obtain the representation matrix of the intermediate time points in the simulation. Specifically, it includes the following steps:

[0048] Standardized expression and Input the VAE model trained in step (2.1) to obtain their respective latent space vector representations. and The expression is as follows:

[0049]

[0050] for For each cell j in the array, calculate its path to... The Euclidean distance D of each cell i is calculated using the following formula:

[0051]

[0052] In the formula, Represented as Each cell j in Represented as Each cell i in the array.

[0053] Find each cell based on this distance. Recent cells Specifically, each terminal cell Match the nearest starting cell i in the latent space:

[0054]

[0055] For matching cells and ,calculate and Linear interpolation of vectors, i.e., the latent space representation of cells at intermediate time points. The formula is as follows:

[0056]

[0057] In the formula, express Each terminal cell Match the nearest starting cell i in the latent space. for Each cell j in the process.

[0058] Using step (2.1), the VAE model decoder is trained to obtain the intermediate time point representation matrix of the simulation. The expression is as follows:

[0059] .

[0060] Step 3: Gene Expression Correction

[0061] For the output matrix Every cell in , its gene vector Sort the cells in descending order of expression level and normalize them to a value X, so that the total number of counts in each cell of the output matrix is ​​X:

[0062]

[0063] ;

[0064] Calculate the cutoff parameter k such that the sum of the expression levels of the first k genes in the sort is Y (the default value for parameter Y is 10000 and Y≤X):

[0065]

[0066] Gene expression after the k-position is set to 0 to simulate dropout in real sequencing. The correction results are evaluated using dropout visualization.

[0067]

[0068] The obtained corrected generator matrix It can be used for subsequent bioinformatics analysis.

[0069] Example 1

[0070] The method described in this invention is applied to simulate and complete human embryonic single-cell sequencing data.

[0071] Single-cell sequencing data for the CS7 stage of gastrulation (HumanGastrula Atlas, E-MTAB-9388) and the CS9 stage of gastrulation (Genome Sequence Archive, HRA013158) were obtained from the public database ArrayExpress as input matrices.

[0072] 1. Simulate and complete the single-cell expression profile data of gastrulation CS8 stage.

[0073] The generated CS8 dataset was integrated with the Zygote-CS7 and CS9-CS16 datasets, and the Seurat function `NormalizeData` was used to normalize each single-cell dataset. The `SelectIntegrationFeatures` function was used to select 4000 hypervariable genes, and the `ScaleData` function was used for scaling. After performing PCA dimensionality reduction on the integrated data using `RunPCA` (npcs=50), the 50 PCA features were corrected using the `fastMNN` algorithm with `IntegrateLayers`. The corrected `fastMNN` features were then used for UMAP embedding and visualization using `RunUMAP` (dims=1:50). The results are shown below. Figure 2 .

[0074] The “CS8_Generation” section in the figure represents the CS8 data generated by this method. The results show that the simulated CS8 matrix can be used for cell clustering and subsequent analysis, and has interpretability for key biomolecular events during development (morphogen expression and transcription factor expression). The portion of the data that is discontinuous between CS7 and CS9 in the dimensionality reduction plot can be used to simulate relevant events during embryonic development.

[0075] Further annotations were provided on the CS8 simulation data generated by this method; the results are shown below. Figure 3 . Figure 3 This figure shows the annotation results of ground truth (Genome Sequence Archive, HRA005567) spatial transcriptome sequencing data from human gastrulation CS8, compared with the annotation results of this method and simulated data generated by different methods, including scDesign3, scVAEDer, Splatter, and Splatter-simple. The bubble chart shows the marker genes for different cell types during annotation, and the UMAP diagram shows the cell clustering.

[0076] Among them, scDesign3 is a statistical modeling framework based on Copula, scVAEDer is a combination of VAE and denoising diffusion model (DDM), and Splatter and Splatter-simple are statistical models based on Gamma-Poisson hierarchical distribution. These models were generated using the officially released usage methods and default parameters. In comparison with other methods, this method demonstrates better biological interpretability and can generate annotated single-cell expression matrices compared to published single-cell transcriptome data simulation methods.

[0077] 2. The single-cell expression profiles of human gastrulation during the E14 stage and organogenesis during the CS10 stage were simulated and completed.

[0078] The E14 data were simulated and completed using E12 data (Gene Expression Omnibus, GSE136447) and CS9 data (Genome Sequence Archive, HRA013158), and compared with the actual E14 data (Gene Expression Omnibus, GSE136447). Similarly, the CS10 data were simulated and completed using CS9 (Genome Sequence Archive, HRA013158) and CS11 data (Genome Sequence Archive, HRA013158), and compared with the actual CS10 data (Gene Expression Omnibus, GSE155121). This demonstrates that this method can complete single-cell data across multiple developmental time points, based on a single developmental start and end point.

[0079] The completion results in E14 are as follows: Figure 4 As shown, the completion result in CS10 is as follows: Figure 5 As shown in the diagram. The bubble chart represents the marker genes for different cell types during annotation, and the UMAP diagram shows the cell clustering. The results indicate that, compared with other methods, the method in this study can better simulate the biological characteristics of single-cell matrices, and the generated single-cell matrices have good annotability.

[0080] Accordingly, the present invention also provides a single-cell time series data simulation and completion system, the system comprising:

[0081] The data preprocessing and training module acquires and integrates single-cell sequencing data from multiple time points, performs cell annotation and expression matrix standardization, and uses the standardized data to train a variational autoencoder model.

[0082] The latent space simulation module, connected to the data preprocessing and training module, encodes the standardized expression matrix into the latent space, establishes cell mapping relationships, generates latent space vectors at intermediate time points through interpolation, and decodes them into standardized expression matrices.

[0083] The expression correction and output module is connected to the latent space simulation module. It performs correction processing on the decoded simulated expression matrix based on the statistical features of real data and outputs the final simulated single-cell gene expression matrix.

[0084] Example 2

[0085] The performance and biological information restoration accuracy of this method, as well as scDesign3, scVAEDer, Splatter, and Splatter-simple simulation and completion methods, were evaluated. The results are as follows: Figure 6 As shown.

[0086] Generation time was evaluated using a matrix of 8000 genes with different cell numbers. Figure 6 (A) This method still maintains a relatively fast running speed on large datasets (running time less than 300 s for 20,000 cells × 8,000 genes). The difference between the generated matrix and the ground truth is compared. Figure 6 In section B), the matrices generated by different methods are integrated with the ground truth and their dimensionality reduced for UMAP visualization, using Silouette contour coefficients ( Figure 6 C), k-nearest neighbor batch effect test kBET ( Figure 6 In the D parameter of the UMAP graph, the distance between two groups is evaluated. A closer distance and a parameter closer to 1 indicate a better simulation result. This demonstrates that the method can largely represent the biological state of the matrix during the simulation process. The comparison of whether the generated matrix accurately reflects the actual biological development process between different time points is also important. Figure 6 The matrix generated by the method is integrated with the starting and ending matrices and its dimensionality reduced for UMAP visualization. The pairwise positional relationships between the three matrices are evaluated and a matrix is ​​constructed. The Pearson correlation coefficient is used to evaluate the relationship between each column of the matrix and the ground truth. Figure 6 The closer F–I in the model is to 1, the better the genealogical relationships in the ground truth are reproduced. It can be seen that this method has a good simulation effect and can reproduce the intermediate biological states in the development process of events to a large extent, and fill in the missing biological processes.

[0087] Obviously, those skilled in the art can make various modifications and variations to this disclosure without departing from its spirit and scope. Therefore, if such modifications and variations fall within the scope of the claims of this disclosure and their equivalents, this disclosure is also intended to include such modifications and variations.

Claims

1. A method for simulating and completing single-cell time series data, characterized in that, Includes the following steps: (1) Obtain single-cell sequencing data at the start and end time points or physiological states, use cell annotation tools to annotate the data by cell type, obtain two corresponding annotated data objects, integrate the two annotated data objects and perform sequencing depth standardization to obtain a standardized expression matrix, and use it as an input matrix to train and optimize the variational autoencoder model. (2) Input the standardized expression matrix into the optimized variational autoencoder model, reduce the dimension and denoise the biological information in the expression matrix to obtain a continuous single-cell level latent space vector representation. Based on the distance minimization criterion, establish the correspondence between cells at the start and end time points in the latent space. For each pair of corresponding cells, use the interpolation method to calculate the latent space vector of the intermediate time point cell. Decode the latent space vector of the intermediate time point through the decoder of the variational autoencoder model to obtain the simulated intermediate time point expression matrix. (3) Based on the single-cell sequencing data of the start and end time points or physiological states, the expression level of the obtained intermediate time point expression matrix is ​​corrected, and finally the simulated intermediate state single-cell expression matrix is ​​output.

2. The method for simulating and completing single-cell time series data according to claim 1, characterized in that, In step (1), the cell annotation tool is the Scanpy tool.

3. The method for simulating and completing single-cell time series data according to claim 1, characterized in that, In step (1), the annotation data object is an AnnData object, which includes cell annotation information stored in the AnnData observation structure and gene information stored in the AnnData feature structure.

4. The method for simulating and completing single-cell time series data according to claim 3, characterized in that, The cell annotation information includes the sample source, collection time point, and physiological state, and the gene information includes the gene name and the total number of gene reads.

5. The method for simulating and completing single-cell time series data according to claim 1, characterized in that, In step (1), the sequencing depth standardization specifically means: standardizing the total number of counts for each cell to 10,000.

6. The method for simulating and completing single-cell time series data according to claim 1, characterized in that, In step (1), the variational autoencoder model is trained and optimized using a composite loss function, which is composed of a weighted sum of a mean squared error loss function and a KL divergence loss function, wherein the weight coefficient of the KL divergence loss function is 0.

5.

7. The method for simulating and completing single-cell time series data according to claim 1, characterized in that, In step (2), the minimization criterion is: using the minimum Euclidean distance between latent space vectors as the criterion, a mapping relationship between cells from the endpoint to the starting time point is established.

8. The method for simulating and completing single-cell time series data according to claim 1, characterized in that, In step (2), the interpolation method is specifically as follows: the midpoint of each pair of corresponding cells is obtained to get the single-cell level latent space vector representation of the intermediate state, and the standardized expression matrix of the sequencing depth at the intermediate time point is decoded in this way.

9. The method for simulating and completing single-cell time series data according to claim 1, characterized in that, The specific steps (3) are as follows: the gene shedding rate of single-cell sequencing data at the start and end time points or physiological states is evaluated, and based on the evaluation results, the gene expression matrix of the simulated intermediate time point matrix is ​​scaled and the low-expression genes are set to zero to correct the expression level, so as to obtain the simulated intermediate state single-cell expression matrix.

10. A single-cell time series data simulation and completion system, characterized in that, The system includes: The data preprocessing and training module acquires and integrates single-cell sequencing data from multiple time points, performs cell annotation and expression matrix standardization, and uses the standardized data to train a variational autoencoder model. The latent space simulation module, connected to the data preprocessing and training module, encodes the standardized expression matrix into the latent space, establishes cell mapping relationships, generates latent space vectors at intermediate time points through interpolation, and decodes them into standardized expression matrices. The expression correction and output module is connected to the latent space simulation module. It performs correction processing on the decoded simulated expression matrix based on the statistical features of real data and outputs the final simulated single-cell gene expression matrix.

Images