A comprehensive performance evaluation method for single-cell RNA sequencing simulator

Abstract

The present application relates to a kind of comprehensive performance evaluation method for single cell RNA sequencing simulator, belong to bioinformatics technical field.The present application is to solve the problem of single cell RNA sequencing simulator evaluation dimension single, missing systematic bias verification.The steps of the present application include: S1: obtaining multiple public single cell RNA sequencing real data set and the simulated data set generated by multiple simulation methods.S2: the data set is preprocessed, and the data format is unified.S3: the comprehensive evaluation index system including data attribute dimension and functional dimension is constructed.S4: the simulated data generated by each method is quantitatively evaluated on data attribute and functional dimension, and the comprehensive score of each method on each dimension is calculated.S5: based on the evaluation result, interactive online platform is constructed, supports user to upload custom data evaluation result and carries out normalization comparative analysis.

Description

Technical Field

[0001] This invention relates to a comprehensive performance evaluation method for a single-cell RNA sequencing simulator, belonging to the field of bioinformatics technology. Background Technology

[0002] Single-cell RNA sequencing (scRNA-seq), a revolutionary high-throughput sequencing technology, can resolve transcriptional heterogeneity in cell populations at single-cell resolution, overcoming the limitations of traditional bulk RNA sequencing and providing a new perspective for life science research. Currently, the number of publicly available scRNA-seq datasets is growing exponentially, and they are becoming increasingly diverse in terms of data type, tissue origin, species coverage, and sequencing depth. To effectively analyze and utilize scRNA-seq data, numerous tools and publications for its analysis have been developed. However, real-world data often faces limitations such as sample scarcity, high acquisition costs, difficulty in completely eliminating batch effects, and lack of cell type labels. Therefore, many methods for simulating single-cell transcriptome data have been developed to compensate for these shortcomings. These simulation methods replicate the key attributes of single-cell transcriptome data through deep generative and statistical models, providing controllable and defined basic facts for the development and validation of scRNA-seq analysis tools.

[0003] These tools typically include evaluations of simulation quality upon publication. However, these results are mostly non-neutral, focusing on certain attributes of the simulated data. These attributes include data properties such as the mean and variance of the log2 (CPM) estimates, library size, and distribution of mean gene values; and functional properties such as the ability to simulate differentially expressed genes and the ability to preserve the true cell population structure. In reality, due to differences in design logic and technical characteristics, the generated data from different simulation tools have varying emphases across different dimensions, and different downstream analysis tasks have different requirements for the simulated data. Therefore, a comprehensive evaluation of the quality of transcriptome data generated by different simulation methods is crucial. Several evaluation studies and platforms have attempted to address this issue, but most suffer from limitations such as a single evaluation dimension and a lack of functional validation. For example, CountsimQC focuses only on the evaluation of data attributes, SimBench evaluates data attributes and the ability of genes to maintain different expression patterns, Crowell et al.'s research considered data attributes and the function of methods for generating different cell batches and clusters, and Duo et al.'s research was more comprehensive in terms of evaluation dimensions, but it did not directly provide a complete data attribute evaluation process, nor did it deeply analyze a systematic bias that is common in data generated by simulation methods, namely, the systematic deviation of signal intensity in some methods from the original data in terms of cell population structure maintenance and differentially expressed gene restoration. Summary of the Invention

[0004] To address the issues of single-cell RNA sequencing simulator evaluation having limited dimensions and lacking systematic bias verification, this invention proposes a comprehensive performance evaluation method for single-cell RNA sequencing simulators.

[0005] The technical solution adopted by the present invention to solve the above problems is as follows: The present invention includes the following steps: Step 1: Obtain multiple real single-cell RNA sequencing datasets and simulated datasets generated by various simulation methods; Step 2: Based on real and simulated datasets, construct a comprehensive evaluation index system that includes data attribute dimensions and functional dimensions; Step 3: Based on the comprehensive evaluation index system, quantitatively evaluate the simulated dataset in terms of data attribute dimension and functional dimension to obtain the evaluation scores of each dimension index; Step 4: Calculate the comprehensive evaluation score of each simulation method, and build an interactive online platform based on the comprehensive evaluation score to support the comparative analysis of the evaluation results.

[0006] The beneficial effects of this invention are: This invention proposes countSimEval as a systematic framework for assessing the simulation quality of single-cell RNA sequencing data. Compared to existing assessment tools, this framework integrates 13 data attribute dimensions and 3 functional dimensions (preserving the true cell population structure, maintaining differentially expressed genes, and simulating batch effects), forming a more comprehensive and detailed assessment system. It can verify the simulation fidelity under different dimensionality reduction strategies and assess the performance of maintaining genes with different expression trends. Furthermore, this invention reveals for the first time the systematic "enhancement" or "weakening" biases of some simulation methods when simulating cell population structure and differentially expressed genes, providing crucial empirical evidence for method developers to identify technical shortcomings and for researchers to accurately select simulation tools. In addition, the established interactive online platform allows users to upload evaluation results of custom data and perform normalized comparative analysis, significantly improving the practical value and scalability of the assessment framework. Attached Figure Description

[0007] Figure 1 This is a flowchart of a comprehensive performance evaluation method for a single-cell RNA sequencing simulator; Figure 2 For the countSimEval model diagram; Figure 3 Scores for 24 single-cell transcriptome simulation methods evaluated by countSimEval in terms of data attributes and various functions, (a) using a statistical model method, (b) using a deep generative model method; Figure 4 Detailed results of the scores of 24 single-cell transcriptome simulation methods evaluated for countSimEval on data attributes; (a) heatmap of the simulation method scores on each sub-attribute, (b) box plot of the method scores on each sub-attribute, (c, d) mean scores of the simulation methods at the cellular & gene level. Figure 5 The data attribute performance of different models and the correlation between methods are shown below; (a) scatter plots of data attribute scores and stability performance of different models on multiple datasets; (b) grouped bar charts of average model performance, with the gray dashed line representing the overall average of all methods; (c) Spearman rank(r) correlation heatmap of different methods. Figure 6Functional performance of 24 single-cell transcriptome simulation methods evaluated by countSimEval; (a, b, c) Bullet plots of mean scores for methods on cell population structure, differentially expressed genes, and batch effect; (d, e, f) Box plots of scores for methods on cell population structure, differentially expressed genes, and batch effect; (g) Correlation between data attribute scores and cell population structure scores (Pearson = 0.82); (h) Correlation between data attribute scores and differentially expressed gene scores (Pearson = 0.42); (i) Correlation between cell population structure scores and differentially expressed gene scores (Pearson = 0.5). Figure 7 To illustrate the strength bias of simulated cell population structure and differentially expressed gene function across different methods, and the Spearman correlation between methods, (a) a box plot of the difference between simulated data clustering score and real data clustering score (positive = over-enhanced cell population structure, negative = weakened cell population structure), categorized according to different dimensionality reduction algorithms. (b) a box plot of the difference between simulated data DEG number and real data DEG number (positive = high DEG number, negative = low DEG number), categorized according to up-regulated and down-regulated genes. (c) a heatmap of Spearman rank(r) correlation between simulated cell population scores from different methods. (d) a heatmap of Spearman rank(r) correlation between scores for maintaining differentially expressed genes from different methods. Figure 8 This is a schematic diagram of the online platform's functions. Detailed Implementation

[0008] like Figure 1 As shown, the steps of the comprehensive performance evaluation method for a single-cell RNA sequencing simulator described in this embodiment include: S1: Obtain multiple publicly available real datasets of single-cell RNA sequencing and simulated datasets generated by various simulation methods; We obtained multiple publicly available real single-cell RNA sequencing datasets and simulated datasets generated by various simulation methods. The real datasets came from multiple public databases, including the 10x Genomics website, NCBI Gene Expression Omnibus (GEO), Figshare platform, and Broad Single Cell Portal database, covering 9 different sequencing technologies and a total of 22 datasets. The simulation methods included 24 methods, divided into statistical model-based methods and deep generative model-based methods, which were used to simulate and generate real data.

[0009] S2: Preprocess the dataset to standardize the data format; Standardized preprocessing includes: basic filtering of raw data to retain genes expressed in at least 3 cells and cells with at least 200 genes detected; standardizing the data format for statistical models to SingleCellExperiment objects and standardizing the data format for deep generative models to AnnData structures.

[0010] S3: Construct a comprehensive evaluation index system that includes data attribute dimensions and functional dimensions; The data attributes cover 13 dimensions, including: AveLogCPM_Tagwise: The mean log(CPM) and tagwise dispersion of gene expression, which is a two-dimensional attribute.

[0011] Mean_Dispersion: The baseline mean of gene expression and the gene-specific dispersion, which is a two-dimensional attribute.

[0012] Mean_var: The mean and variance of the gene log(CPM) estimate, which is a two-dimensional attribute.

[0013] Lib: The sum of the expression levels of all genes in each cell (library size), which is a one-dimensional attribute.

[0014] TMM: Mean M value after pruning (standardization factor), which is a one-dimensional attribute.

[0015] EffLibsize: The total cell count multiplied by its corresponding TMM normalization factor (effective library size), is a one-dimensional attribute.

[0016] AveLogCPM: The average abundance value of genes. Abundance is represented by the log(CPM) value calculated by edgeR, which is a one-dimensional attribute.

[0017] sampleFraczero: The zero-proportion observation for each cell (column) in the counting matrix, which is a one-dimensional attribute.

[0018] featureFraczero: The zero-proportion observation for each gene (row) in the counting matrix, which is a one-dimensional attribute.

[0019] sampleCorrelation: The log(CPM) value calculated based on edgeR, representing the distribution of Spearman correlation coefficients between cell pairs, and is a one-dimensional attribute.

[0020] featureCorrelation: The log(CPM) value calculated based on edgeR, representing the distribution of Spearman correlation coefficients between gene pairs, is a one-dimensional attribute.

[0021] Lib_Fraczero: The correlation between the size of the cell library and the number of zeros observed per cell is a two-dimensional attribute.

[0022] AveLogCPM_Fraczero: The correlation between the average abundance and the zero fraction observed for each gene. Abundance is represented by the log(CPM) value calculated by edgeR and is a two-dimensional attribute.

[0023] One-dimensional attributes are those that describe and evaluate the distribution and characteristics of a single variable and do not involve the relationship between variables; two-dimensional attributes involve the relationship between two variables and focus on the collaborative change or dependency pattern between variables, rather than the independent characteristics of a single variable.

[0024] The three dimensions encompassed by functionality include: Cell population structure preservation capability: assessing whether simulated data can accurately reproduce the cell population structure in real data.

[0025] Capacity to maintain differentially expressed genes: assessing whether simulated data can preserve differentially expressed signals between cell populations in real data.

[0026] Batch effect simulation capability: Evaluating whether simulated data can reproduce the batch effect structure in real data.

[0027] S4: Quantitatively evaluate the simulation data generated by each method in terms of data attributes and functionality, and calculate the comprehensive score of each method in each dimension; The comprehensive score calculation for data attributes includes: S30101: Use the following 7 indicators to evaluate the quality of data attributes, where 1, 2, 3, 4, and 5 are used to evaluate one-dimensional attributes, and 4, 5, 6, and 7 are used to evaluate two-dimensional attributes.

[0028] (1) Kolmogorov-Smirnov tests: The KS test is used to measure the maximum difference between two sample distributions. The difference is calculated using the following formula: For the sample X = {x1, x2... }, its empirical distribution function Defined as: (1); For another sample Y={y1,y2... }, its empirical distribution function for: (2); The KS statistic is the maximum vertical distance between two empirical distribution functions: (3); Where sup represents taking the supremum (i.e. the maximum difference) for all possible x, D∈[0,1], D=0 means that the two distributions are completely identical, and D=1 means that they do not overlap at all.

[0029] (2) Scaled area between eCDFs: The empirical distribution function of the two distributions is and The scaled area between eCDFs (SAE) is defined as follows: (4); The smaller the SAE, the more similar the two distributions are; the larger the SAE, the greater the difference between the two distributions.

[0030] (3) Wald-Wolfowitz runs test: This test is a nonparametric test used to examine whether two independent samples come from the same distribution. The basic idea is to merge and sort the values ​​of the two samples, and then convert them into a binary sequence based on which sample they came from. Next, the number of runs in this binary sequence is calculated, where a run is a sequence of consecutive identical elements (0 or 1). Under the null hypothesis, i.e., the two samples come from the same distribution, the number of runs should be close to a certain expected value. If the number of runs deviates significantly from this expected value, then this invention has reason to reject the null hypothesis. The specific steps are as follows: The values ​​of the two samples are merged and sorted in ascending order. Based on which sample the value comes from, the sorted values ​​are converted into a binary sequence and the number of runs in the binary sequence is calculated.

[0031] Calculate the expected value and standard deviation of the number of runs. For a length of... The first sample and its length are The second sample, the expected value E(R) and standard deviation of the number of runs. They are respectively: (5); (6); Calculate the Wald-Wolfowitz runs test statistic Z, which is the difference between the number of runs and the expected value divided by the standard deviation: (7); Where R is the actual number of runs. A smaller Z value indicates that the same set of data tends to cluster together after sorting, indicating a poor degree of mixing; a larger Z value indicates that the two sets of data are mixed more evenly after sorting, indicating a good degree of mixing.

[0032] (4) k-Nearest Neighbors Rejection Fraction: This statistic is used to assess the degree of mixing between two datasets. The calculation process is as follows: Dataset: X (real data, n samples), Y (simulated data, m samples) Total sample size: N = n + m Reference subset size: R = min(500, N) (Randomly select R samples for testing) Nearest neighbor count: k = max(5, 0.01N) (Number of nearest neighbors for each sample) Expected proportion (under the null hypothesis, the composition of nearest neighbors should be consistent with the global composition): ,

[0033] For each sample Perform a chi-square test: calculate The distance to all samples is used to select the k nearest samples. Let the number of samples from X in the nearest neighbors be .

[0034] The number of samples from Y is , + =k, the formula for the chi-square test statistic is: (8); Among them, the expected value =k· , =k·

[0035] like (With 1 degree of freedom and a significance level of 0.05), then reject.

[0036] NN Rejection Fraction = The higher the NN Rejection Fraction, the worse the mixing effect of the two datasets; the lower the NN Rejection Fraction, the better the mixing effect of the two datasets.

[0037] (5) Average silhouette width: Randomly select R samples, and for each sample, calculate the Euclidean distance with all other samples, and calculate the silhouette width based on these distances.

[0038] For the i-th sample point, the silhouette coefficient is defined as: (9); in a(i): The average distance between sample i and all other samples within the same source data. b(i): The average distance between sample i and all other samples from different sources. The average profile width is the average of the profile coefficients for all data points, i.e.: (10); Where n is the total number of data points, and the contour coefficient ranges from [-1, 1]. If the contour coefficient is closer to ±1, it means that the simulated data and the original data are not mixed, and the simulation effect is not good; if the contour coefficient is closer to 0, it means that the simulated data and the original data overlap, and the simulation effect is better.

[0039] (6) Multivariate KS statistic: The multivariate Kolmogorov-Smirnov test used in this invention is a generalization of the classic KS test proposed by Fasano and Franceschini in two-dimensional and three-dimensional space, used to test whether two independent samples come from the same distribution. The basic idea is to compare two samples X (of size m) and Y (of size n) in multidimensional space d. The maximum difference in the cumulative distribution function (eCDF) in section 2 is used to determine whether the distributions are the same. The specific calculation process is as follows: Let sample X = { , ,..., } and Y={ , ,..., } is a two-dimensional array, and after merging, a coordinate sequence is generated. ), calculate its cumulative distribution difference in the four quadrants: (x < ,y< ), (x< ,y> ), (x> ,y< ), (x> ,y> ).

[0040] Calculate the eCDF of sample X in quadrant Q: (11); Calculate the eCDF of sample Y in quadrant Q: (12); Where 1{ } is an indicator function.

[0041] Take the maximum absolute difference in eCDF between two samples across all quadrants: (13); Standardized statistics: (14); When Z is smaller, the difference in the cumulative distribution function (eCDF) between the two samples is smaller, and the data are similar; when Z is larger, the difference in the cumulative distribution function (eCDF) between the two samples is more significant, and the data are not similar.

[0042] (7) Kernel density based global two-sample comparison test: This test is a non-parametric statistical method used to estimate the shape of the probability density function. It estimates the density function of the entire dataset by placing a "kernel" (usually a Gaussian function) around each data point and then summing these kernel functions. It tests whether two samples come from the same distribution by comparing the global differences in their density functions. The basic calculation steps are as follows: For two datasets X={ , ,..., } and Y={ , ,..., }, calculate their kernel density respectively Degree estimation The formula for estimating kernel density is: (15); Where K is the kernel function and H is the bandwidth matrix.

[0043] The test statistic T is defined as the integral of the squares of the differences between the two density functions: (16); This statistic measures the difference between the two density functions.

[0044] Substituting the kernel density estimate into the formula for the test statistic, we obtain the estimated value of the test statistic: (17); in: (18); (19); (20); (twenty one); Under the null hypothesis : Lower test statistic The distribution is close to normal. Specifically, when , At that time, there were: (twenty two); Mean correction term: (twenty three); Variance estimation: (twenty four); The standardized statistic Z is: (25); When Z is smaller, it indicates that the difference in the density function between the two samples is close to the mean correction term and the difference is not significant; when Z is larger, it indicates that the difference in the density function between the two samples is far from the mean correction term and the difference is significant.

[0045] S30102: Standardized Indicator Score Calculation For each evaluation metric, the difference between the real data and the simulated data on that attribute is calculated. If the difference is in the range [0,1], no processing is performed. Otherwise, the scores of the same attribute from multiple methods from a dataset are first Z-score standardized, and then min-max normalized to the range [0,1] to ensure that the higher the score of each metric, the more similar it is to the original data.

[0046] S30103: Multi-indicator aggregation The arithmetic mean of the scores for all corresponding evaluation metrics under this attribute is taken: (26); Where m is the number of evaluation indicators applicable to this attribute.

[0047] S30104: Calculation of Overall Score Comprehensive score of data attributes: The average score of the 13 attributes is taken. (27); The quantitative assessment of the functional aspects of maintaining the original cell population structure includes the following steps: S30201: First, map the high-dimensional gene expression matrix to the low-dimensional feature space using PCA or t-SNE, and then calculate the Euclidean distance between cells.

[0048] S30202: Use the following five indicators to assess cell population structure quality: (1) Average silhouette width: Used to evaluate the quality of clustering results, measuring the similarity of each data point to its own cluster and other clusters. The calculation process is as follows: For the i-th sample point, the silhouette coefficient is defined as: (28); In formula (28), a(i) is the average distance between sample i and all other samples in the same cluster, and b(i) is the average distance between sample i and all other samples in different clusters. The average profile width is the average of the profile coefficients for all data points, i.e.: (29); Where n is the total number of data points, and the silhouette coefficient ranges from [-1, 1]. The closer the silhouette coefficient is to 1, the more correctly the data points are assigned to their respective clusters, and the better the separation between clusters. If the silhouette coefficient is close to a negative value or 0, it indicates that the data points may have been assigned to the wrong clusters, or the separation between clusters is insufficient.

[0049] (2) Dunn Index: This is an internal validity metric used to evaluate cluster quality. Its core idea is to measure the balance between inter-cluster separation and intra-cluster compactness. The specific calculation steps are as follows: Find the minimum distance among all cluster pairs: (30); in and It can be any two different clusters. It is the minimum distance between all sample point pairs between them.

[0050] Take the maximum value of the maximum diameter of all clusters: (31); in Represents the k-th cluster, This represents the maximum distance between all pairs of sample points within it, i.e., the diameter of the cluster.

[0051] The Dunn index is defined as the ratio of the minimum inter-cluster distance to the maximum intra-cluster diameter: (32); A higher Dunn value indicates better inter-cluster separation, higher intra-cluster density, and better clustering results; a lower Dunn value indicates overlap between clusters or loose clustering within clusters, resulting in poor clustering performance.

[0052] (3) Connectivity: Used to evaluate clustering effectiveness, it measures the density of similar samples within their neighborhood, reflecting the local consistency of the clustering results. The specific calculation process is as follows: For dataset X={ } and clustering results C={ },set up For point The set of m nearest neighbors, for each point Count the number of its m nearest neighbors that do not belong to the same cluster: (33); in, for The cluster to which it belongs, 1( ) is an indicator function.

[0053] Take the average connectivity of all points: (34); The range of connectivity values ​​is [0, +]. When connectivity is 0, it represents perfect clustering, meaning all nearest neighbors are in the same cluster. The larger the value, the worse the clustering density and the more dispersed the clustering effect.

[0054] (4) Calinski-Harabasz Index: Used to evaluate clustering performance. Its design logic is based on the principle that "the essence of clustering is to maximize inter-cluster variance and minimize intra-cluster variance." It quantifies the global clustering quality by calculating the ratio of inter-cluster variance to intra-cluster variance. The formula is as follows: (35); Where k is the number of cell groups (clusters), and n is the total number of cells. The sum of squared differences between clusters is calculated by multiplying the squared distance between the center point of each cluster and the global center point by the number of cells in that cluster, and then summing the results. The sum of squared deviations from the mean within a cluster is calculated by squared the distance between each cell and the center of its cluster, and then summing the results over all cells.

[0055] The CHI value ranges from (0, +∞). A larger CHI value indicates a greater contribution of inter-cluster differences to the total variance, smaller intra-cluster noise, and better clustering quality.

[0056] (5) Davies-Bouldin Index: The core idea is to measure the average similarity between each cluster and the other most similar clusters. Similarity is the ratio of intra-cluster dispersion to inter-cluster distance, as shown in the following formula: (36); k: Number of cell groups (clusters) Cluster The internal dispersion, i.e., the average Euclidean distance from all cells within a cluster to its center point. Cluster and Inter-cluster Euclidean distance The DBI value ranges from (0, +∞). The smaller the DBI value, the better the clustering quality and the clearer the cell population structure.

[0057] S30203: Standardized Indicator Score Calculation For each evaluation metric, the difference between the real data and the simulated data is calculated. If the difference is within the range of [0,1], no processing is performed. Otherwise, the scores of the same attribute from multiple methods in a dataset are first Z-score standardized, and then min-max normalized to the range of [0,1] to ensure that the higher the score of each metric, the more similar it is to the original data.

[0058] S30204: Calculation of Overall Score Cell population structure comprehensive score: The average score of the five assessment indicators is used. (37); Quantitative assessment of the functional maintenance of differentially expressed genes includes the following steps: S30301: Prepare real single-cell expression data and simulated single-cell expression data, ensuring that the cell population annotations of the two types of data are completely consistent.

[0059] S30302: Use the edgeR package to perform differential expression analysis on real and simulated data, and perform pairwise comparisons on different cell populations in the data to identify upregulated and downregulated genes.

[0060] S30303: Calculate the proportion of differentially expressed genes in the corresponding cell populations of real and simulated data for upregulated and downregulated genes respectively, using the following formulas: (38); In formula (38), These are the corresponding upregulated / downregulated genes in the simulated data. These are the upregulated / downregulated genes in the original data.

[0061] S30303: If there are more than two cell populations, the proportions of the different population pairs are summed, and the final score is used as the score of the method for maintaining differentially expressed genes on this dataset.

[0062] The quantitative assessment of batch effects in functionality includes the following steps: S30401: First, map the high-dimensional gene expression matrix to the low-dimensional feature space using PCA or t-SNE, and then calculate the Euclidean distance between cells.

[0063] S30402: Use the following five indicators to assess batch-effect quality of cell distance: (1) k-nearest neighbor batch effect test: used to evaluate batch effect. The basic idea is that if there is no batch effect, then for any cell, its k nearest neighbors should be like a small "random sampling", where the proportion of each batch is roughly the same as the proportion of batches in the whole dataset.

[0064] If batch effects exist, then a cell's neighbors will primarily come from its own batch, because technological biases cause cells from the same batch to cluster more closely together in space.

[0065] The k-nearest neighbor batch effect test in this invention is implemented using the 'kBET'R package.

[0066] (2) Average silhouette width: Similar to the Average silhouette width used to evaluate clustering above, batch labels are treated as cluster labels in this section. The closer the ASW value is to 1, the more significant the batch effect of the dataset. If the silhouette coefficient value is close to negative or 0, the batch effect is more ambiguous.

[0067] (3) Local Inverse Simpson's Index: An index for judging the degree of mixing of different batches of cells in a local region. For each cell i in the dataset, select the k nearest cells in the low-dimensional feature space to form the local nearest neighbor set of that cell. The local Simpson index is defined as: (39); m: Batch quantity Local nearest neighbor set The proportion of cells in category c Take the reciprocal of the local Simpson index: (40); Calculate the mean LISI of all cells as an indicator of the overall dataset. (41); The mean LISI ranges from (1, m). A higher value indicates a more ambiguous batch effect, while a lower value indicates a more significant batch effect.

[0068] (4) Adjusted Rand Index: To quantify the batch effect, this invention first uses the K-means clustering algorithm to group cells, obtain cluster labels, and compares the obtained cluster labels with the true batch labels to calculate the Adjusted Rand Index (ARI). The Adjusted Rand Index is calculated by comparing the distribution of all sample pairs in two partitions (clustering result U and batch label V). Its calculation formula is as follows: (42); in, The number of samples that simultaneously belong to both cluster Ui and batch Vj. = For clustering The total number of samples in the sample, = For batch The total number of samples in the sample, where n is the total number of samples. () represents the number of combinations of selecting 2 from n samples. The range of ARI is [-1, 1]. The closer the ARI is to 1, the more significant the batch effect; the closer the ARI is to 0, the less significant the batch effect.

[0069] S30403: For each simulated data, calculate the evaluation index. If the index is in the range of [0,1], no processing is performed. Otherwise, the scores of the same attribute from multiple methods in a dataset are first Z-score standardized, and then min-max normalized to the range of [0,1] to ensure that the higher the score of each index, the more similar it is to the original data.

[0070] S30404: Calculation of Overall Score Batch effect composite score: The average score of the four evaluation indicators: (43).

[0071] S5: Based on the evaluation results, build an interactive online platform that allows users to upload custom data evaluation results and perform normalized comparative analysis.

[0072] This invention utilizes JavaWeb in conjunction with Rserve to establish an online platform. This platform clearly displays the detailed performance of 24 simulation methods on 22 corresponding real datasets, focusing on data attributes (including 6 at the cellular level and 7 at the gene level, totaling 13 sub-attributes) and functionality (including cell population structure, biological signals, and batch effects). The platform includes scores for each dimension, method mean, and ranking. The scores of the 24 single-cell transcriptome simulation methods evaluated by countSimEval for data attributes and functionality are shown below. Figure 3As shown, it supports users filtering datasets by cell count for comparison. The countSimEval model is as follows. Figure 2 As shown, the detailed results of the data attribute scores for 24 single-cell transcriptome simulation methods evaluated by countSimEval are as follows: Figure 4 As shown, the platform also supports users in expanding analysis dimensions. Researchers can upload scores for data attribute dimensions that they generate themselves and that have been evaluated by the countSimEval tool to the online platform. After uploading, the backend will call the built-in normalization algorithm to normalize these external data against the benchmark scores of all 24 methods in this invention.

[0073] Example This embodiment mainly conducted experiments on 24 simulation methods in 22 real scRNA-seq datasets derived from 9 different sequencing technologies and evaluated them from the following dimensions: (1) Data attributes: basic statistical attributes and distribution of each data, including 6 at the cellular level and 7 at the gene level, totaling 13 sub-attributes. (2) Functionality: the ability of the simulation data to support downstream analysis tasks, including cell population structure (to what extent the simulation data can maintain the cell population structure of the original data), differentially expressed genes (the ability of the simulation data to maintain the biological signals of the original data), and batch effects (to what extent the simulation data generated by the simulation method can simulate batch effects).

[0074] The results show significant performance heterogeneity among different simulation methods in terms of data attributes. ZINB-WaVE achieved the highest data attribute score. Although its initial design goal was to extract low-dimensional signals from scRNA-seq data rather than for simulation generation, it demonstrated the highest accuracy in this evaluation and performed well across all sub-attributes. Similarly, this method has been proven in other studies to generate excellent simulated data. SPARSim, scDesign3, and scDesign2 also showed relatively good performance across all metrics. In contrast, hierarchicell, SymSim, and scDD ranked lower across all metrics. Splat provides a very popular simulation method, with over 600 related publications. The Splat framework, Splatter, integrates multiple simulation methods and defines a statistical model generation process that includes parameter estimation and data generation. However, Splat's performance in most data attributes was only at a mid-level. Overall, deep generative model-based methods perform only moderately well in simulating data attributes compared to statistical model-based methods. While the proliferation of deep generative models in recent years has led to a surge in their application, it must be acknowledged that their performance in simulating data attributes is only moderate due to the large amount of data required for training deep neural networks and their sensitivity to variance. Among data attributes, those related to library size are mostly simulated well (Lib, TMM, EffLibsize, Lib_Fraczero), resulting in higher cell-level accuracy scores than gene-level accuracy scores for most methods. Regarding attributes like sampleCorrelation, Mean_Var, and AveLogCPM_Fraczero, most methods do not specifically consider these attributes, meaning only a few methods with high overall scores can capture these attributes well.

[0075] In the field of single-cell transcriptome data generation, a wealth of statistical probability distribution models and deep generation models exist for constructing generation methods. These different underlying models are key factors affecting their ability to simulate data attributes. Furthermore, even when using the same underlying model, differences in technical approaches and engineering details can also impact generation quality. Therefore, this invention categorizes methods according to different models and evaluates their scores and stability on multiple datasets. Figure 5In this study, it was found that ZINB-WaVE (0.87, 25.225), using a zero-inflated model, combines good performance and stability. However, POWSC (0.619, 10.843), also a zero-inflated model, only ranks in the middle. Similarly, SPARSMim (0.855, 21.864), using GMH, performs excellently, while MOSMim (0.778, 11.427), also using GMH, performs well in accuracy but its stability is only average among all methods. This indicates that even when sharing the same probability distribution assumptions or models, different methods may still exhibit significant differences in accuracy and stability. This difference does not necessarily stem from the model category itself, but rather from the overall technical approach and implementation details. scDesign3 (0.849, 24.978) using GAMLSS and scDesign2 (0.836, 14.82) using Optimal fitting both demonstrate excellent performance in terms of accuracy and stability. Both employ a technique of fitting gene marginal distributions (using different models to fit the distribution for each gene) combined with gene correlation modeling. Simulation methods using GP models, such as ESCO, splat, SCRIP-GP-commonBCV, and SCRIP-GP-trendedBCV, are all around average in terms of accuracy and stability. This can be attributed to the fact that they not only use GP models as the underlying distribution assumptions but also share similar data generation techniques.

[0076] To further explore the correlation between methods, this invention calculated the Spearman correlation coefficients between different methods and plotted a correlation heatmap. Figure 5 c) The results show that SCRIP-BP, splat, SCRIP-GP-commonBCV, SCRIP-BGP-commonBCV, and ESCO exhibit a certain degree of correlation (mean Spearman = 0.517). This is mainly due to their similar simulation steps and shared use of GP and its variant models. Similarly, scDesign3 and scDesign2, based on similar strategies, also showed a moderate degree of correlation (Spearman = 0.612), and LSH-GAN and scGAN, based on the same deep generative model GAN, also showed a certain degree of correlation (Spearman = 0.543).

[0077] The significant value of simulated data lies in supporting downstream biological analysis tasks. This invention systematically evaluates and analyzes three functional requirements of the method: cell population structure preservation, maintenance of differentially expressed genes, and ability to simulate batch effects. Figure 6 In the scoring of simulated cell population structure, scDesign2 (0.937) received the highest score and had a small range (0.099), indicating that scDesign2 is not only highly accurate in simulating cell population structure but also performs stably and reliably. Following closely behind were scDesign3 (0.828), POWSC (0.825), ZINB-WaVE (0.812), and GLMsim (0.790), all of which performed well. The scores for maintaining differentially expressed genes showed a wide range, with scVI (0.957), based on a deep generative model, receiving the highest score. ZINB-WaVE (0.928), GLMsim (0.913), scDesign2 (0.913), and scDesign3 (0.9) also performed well, while ESCO only scored (0.209). In terms of batch effect simulation, SCRIP-BP achieved the highest score (0.896) and had a small range (0.195), indicating that SCRIP-BP has a strong batch effect simulation capability and stable performance.

[0078] This invention reveals the intrinsic relationships among the various assessment dimensions through further correlation analysis. It found a strong positive correlation between data attribute scores and cell population structure scores (Pearson r = 0.82), while the positive correlations between data attributes and differentially expressed gene scores (Pearson r = 0.42) and between cell population structure scores and differentially expressed gene scores (Pearson r = 0.5) were relatively low. This indicates that the simulation accuracy of data attributes is fundamental to high-quality simulation of cell population structure, but its impact on differentially expressed genes is relatively limited. Furthermore, no correlation was found between batch effect scores and data attribute scores (Pearson r = -0.004, P = 0.992). Although batch effect scores also showed some correlation with cell population structure (Pearson r = 0.764, P = 0.236) and differentially expressed gene scores (Pearson r = 0.686, P = 0.314), the relatively high P-values ​​indicate that a stable association could not be confirmed.

[0079] In single-cell data simulations, accurately reproducing the cell population structure and differentially expressed genes from the original data is crucial for the reliability of downstream analyses. However, most existing studies focus on the overall similarity between simulated and real data, neglecting potential systematic biases in functional strength. These biases manifest as follows: some methods may overemphasize cell population structure (producing overly segregated, denser clusters) or differentially expressed signals (generating excessive DEGs), while others may weaken these key features. Although such strength biases are not easily detected in routine assessments, they can affect downstream analytical tasks and further mislead biological discoveries based on simulated data.

[0080] To this end, this invention conducts an in-depth analysis from two key dimensions: First, by quantifying the differences between simulated and original data in the strength of cell population clustering structure, it identifies the "enhancement" or "weakening" trends and their degrees in the simulation process of each method; second, by comparing the relative proportions of the number of differentially expressed genes (DEGs) in simulated and original data, it assesses whether there are systematic biases of "excessive" or "insufficient" DEG signal restoration in each method. Based on this, this invention further quantifies the consistency of different methods in the aforementioned functional performance to reveal the performance correlation patterns between methods.

[0081] To quantify the deviation of the simulated method from the original data's cell population clustering strength, the five simulated cell population evaluation metrics of this invention were first used to calculate the scores of the original and simulated data on PCA and TSNE. For each metric, the method's score was first z-score standardized and then normalized to the [-1, 1] interval. Then, the clustering score of the dataset was calculated using the five evaluation metrics: score = (Average silhouette width + Dunn Index - Connectivity + Calinski-Harabasz Index - Davies-Bouldin Index) / 5. If the difference between the simulated data clustering score and the original data score is positive, it indicates that the simulated data excessively strengthens the cell population structure; conversely, if the difference is negative, it indicates that the simulated cell population structure is weaker than the original data. Figure 7As can be seen in section a, the overall difference values ​​for methods such as scDesign3, SPsimSeq, and ESCO are negative, indicating that the simulated cell population structure strength is generally weaker than the original data. Even though some of these methods score highly in simulating cell populations, it still shows that these methods have a good ability to simulate cell populations, but they generally exhibit a weaker simulated cell population structure strength. Similarly, the overall difference values ​​for methods such as POWSC, ZINB-WaVE, GLMsim, and SCRIP-GP-trendedBCV are positive, indicating that the simulated cell population structure strength is generally stronger than the original data.

[0082] To quantify the deviation of DEG signal intensity between the simulated and original data, the number of upregulated and downregulated differentially expressed genes between each pair of cell populations in both the original and simulated data was first calculated. Then, the similarity between the upregulated and downregulated genes was calculated using the formula: Differential gene ratio = 1 - (Number of differentially expressed genes in simulated data - Number of differentially expressed genes in original data) / (Sum of both). If there were more than two cell populations, the ratios of the different population pairs were summed and the mean was calculated. A positive ratio indicates that the method has an enhancing effect on maintaining the differentially expressed gene, meaning that more differentially expressed genes of that type were generated than in the original data. A negative ratio indicates that the method weakened the differentially expressed signal, meaning that fewer differentially expressed genes of that type were generated than in the original data. Figure 7 As shown in b, while methods such as ZINB-WaVE, POWSC, and scDiffusion can maintain differentially expressed genes among cell populations relatively well, the overall proportion of differentially expressed genes is negative, indicating that the intensity of differentially expressed signals is weaker than the original data. In contrast, methods such as GLMsim, scDesign2, and MOSIM show an overall positive proportion of differentially expressed genes, indicating that they simulate excessively strong differentially expressed signals.

[0083] To further explore the functional correlations among the methods, this invention calculated the Spearman correlation coefficients between cell population scores and differentially expressed gene scores of different model methods and plotted a correlation heatmap. Figure 7 As can be seen from c, most methods show a positive correlation with each other in terms of simulated cell population scores, indicating that different platform methods have a certain similarity trend in their performance in simulating cell population structure. Muscat and SPARSMim (Spearman = 0.943), scDesign2 and GLMsim (Spearman = 0.886), and scDesign3 and POWSC (Spearman = 0.829) show strong Spearman correlations. Figure 7As can be seen in d, the Spearman correlation of differentially expressed gene scores is extremely strong among ESCO, GLMsim, muscat, and scDesign2 methods (mean Spearman = 0.929).

[0084] Furthermore, to facilitate researchers' intuitive understanding of the applicable scenarios and performance differences of various simulation methods, this invention establishes an online platform. This platform can clearly display the detailed performance of 24 simulation methods on 22 corresponding real datasets, focusing on data attributes (including 6 at the cellular level and 7 at the gene level, totaling 13 sub-attributes) and functionality (including cell population structure, biological signals, and batch effects). This includes scores for each dimension, method mean, and ranking. It also supports users filtering datasets by cell count for comparison. Figure 8 As shown, the platform also supports users in expanding analysis dimensions. Researchers can upload scores for data attribute dimensions generated by themselves and evaluated using the countSimEval tool to the online platform. After uploading, the backend calls the built-in normalization algorithm to normalize these external data against the benchmark scores of all 24 methods in this invention. This process ensures the direct comparability of evaluation results from different sources, ultimately generating a comprehensive cross-sectional comparison report for users that integrates built-in methods and custom data.

[0085] This invention utilizes countSimEval, a comprehensive single-cell transcriptome simulation data quality assessment platform, to systematically evaluate the performance of 24 mainstream single-cell RNA sequencing data simulation methods. The evaluation is based on 22 experimental datasets, covering 13 data attribute dimensions and three functional dimensions: preservation of realistic cell population structure, maintenance of differentially expressed genes, and batch effects. Furthermore, this invention establishes an interactive online platform where users can filter datasets of interest to view the performance of simulation methods. The platform also supports users uploading evaluation results for custom datasets, enabling convenient cross-sectional comparisons through normalization with existing data. Compared to existing evaluation tools, countSimEval offers a more comprehensive evaluation scope and a more detailed functional evaluation design. Combined with the support of an interactive platform, it not only enhances the flexibility of user-defined data analysis but also strengthens the practical value of the evaluation results.

[0086] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention in any way. Although the present invention has been disclosed above with reference to preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some modifications or alterations to the above-disclosed technical content to create equivalent embodiments without departing from the scope of the present invention. Any simple modifications, equivalent substitutions, and improvements made to the above embodiments without departing from the scope of the present invention, based on the technical essence of the present invention and within the spirit and principles of the present invention, shall still fall within the protection scope of the present invention.

Claims

1. A comprehensive performance evaluation method for a single-cell RNA sequencing simulator, characterized in that, include: Step 1: Obtain multiple real single-cell RNA sequencing datasets and simulated datasets generated by various simulation methods; Step 2: Based on the real dataset and the simulated dataset, construct a comprehensive evaluation index system that includes data attribute dimensions and functional dimensions; Step 3: Based on the comprehensive evaluation index system, quantitatively evaluate the simulated dataset in terms of data attribute dimension and functional dimension to obtain the evaluation scores of each dimension index; Step 4: Calculate the comprehensive evaluation score of each simulation method, and construct an interactive online platform based on the comprehensive evaluation score to support the comparative analysis of the evaluation results.

2. The comprehensive performance evaluation method for a single-cell RNA sequencing simulator according to claim 1, characterized in that, In step 1, real datasets of single-cell RNA sequencing are obtained by querying public databases. These real datasets cover 22 datasets of 9 different sequencing technologies. Simulated datasets are generated from the real datasets of single-cell RNA sequencing using statistical modeling and deep generative modeling methods. There are 24 simulation methods for these simulated datasets.

3. The comprehensive performance evaluation method for a single-cell RNA sequencing simulator according to claim 1, characterized in that, Step 2 specifically includes: Step 2.1: Preprocess the real RNA sequence dataset and the simulated dataset. The preprocessing includes basic filtering of the real RNA sequence dataset and the simulated dataset, retaining genes expressed in at least 3 cells and cells in which at least 200 genes are detected, and unifying the data format used for statistical models in the basic filtered real RNA sequence data and simulated data into SingleCellExperiment objects, and unifying the data format used for deep generative models into AnnData structures. Step 2.2: Establish a comprehensive evaluation index system that includes data attribute dimensions and functional dimensions based on the preprocessed real RNA sequence dataset and simulated dataset.

4. The comprehensive performance evaluation method for a single-cell RNA sequencing simulator according to claim 3, characterized in that, The data attribute dimensions include 13 dimensions, which are divided into one-dimensional attributes and two-dimensional attributes. One-dimensional attributes describe and evaluate only the distribution and characteristics of a single variable, without involving the correlation between variables; two-dimensional attributes involve the correlation between two variables, focusing on the attributes of collaborative changes or dependency patterns between variables. The two-dimensional attributes include: AveLogCPM_Tagwise, Mean_Dispersion, Mean_var, Lib_Fraczero, and AveLogCPM_Fraczero; AveLogCPM_Tagwise represents the mean log and tagwise dispersion of gene expression, Mean_Dispersion represents the baseline mean and gene-specific dispersion of gene expression, Mean_var represents the mean and variance of gene log estimates, Lib_Fraczero represents the correlation between cell library size and the observed zero fraction per cell, and AveLogCPM_Fraczero represents the correlation between mean abundance and the observed zero fraction per gene; One-dimensional attributes include: Lib, TMM, EffLibsize, AveLogCPM, sampleFraczero, featureFraczero, sampleCorrelation, and featureCorrelation; Lib represents the sum of all gene expression levels in each cell, TMM represents the pruned mean M value, EffLibsize represents the total count of each cell multiplied by its corresponding TMM normalization factor, AveLogCPM represents the mean abundance of genes, sampleFraczero represents the zero proportion of observations in each cell in the counting matrix, featureFraczero represents the zero proportion of observations in each gene in the counting matrix, sampleCorrelation represents the Spearman correlation coefficient distribution between log values ​​calculated based on edgeR and cell pairs, and featureCorrelation represents the Spearman correlation coefficient distribution between log values ​​calculated based on edgeR and gene pairs. The functional dimensions include three dimensions: cell population structure preservation capability, differentially expressed gene maintenance capability, and batch effect simulation capability. Cell population structure preservation capability is used to assess whether the simulated data can accurately reproduce the cell population structure in the real data; differentially expressed gene maintenance capability is used to assess whether the simulated data can retain the differential expression signals between cell populations in the real data; and batch effect simulation capability is used to assess whether the simulated data can reproduce the batch effect structure in the real data.

5. The comprehensive performance evaluation method for a single-cell RNA sequencing simulator according to claim 1, characterized in that, Step 3 involves quantitative evaluation along the data attribute dimensions, including: Step 3.1: Establish indicators one through seven to evaluate the data indicators of the comprehensive evaluation indicator system. Indicators one through five are used to evaluate one-dimensional attributes, and indicators four through seven are used to evaluate two-dimensional attributes. Indicator 1 is used to measure the maximum difference between two sample distributions, including: Based on the sample X={x1,x2... } and sample Y={y1,y2... Calculate the KS statistic between samples using the empirical distribution function. To obtain the maximum difference between the distribution of simulated data and real data; KS statistic The calculation formula is: （1）； In formula (1), For all possible observation samples Take the maximum difference, D∈[0,1], where D=0 indicates that the two distributions are completely identical, and D=1 indicates that they do not overlap at all. Let X be the empirical distribution function of the sample. Let Y be the empirical distribution function of the sample. Indicator 2 is used to calculate the area of ​​overlap of the empirical distribution function of a one-dimensional attribute. The smaller the area of ​​overlap, the more similar the distributions of the simulated data and the real data are; the larger the area of ​​overlap, the greater the difference between the two distributions. The formula for calculating the overlapping area is: （2）； In formula (2), The empirical distribution function of the simulated data, The empirical distribution function of the real data; Indicator 3 is used to test the degree of mixing between samples, including: After merging the values ​​of the two samples, sort them in ascending order. Convert the sequence into a binary sequence according to the sample source, calculate the number of runs in the binary sequence, and obtain the standard deviation and expected value of the number of runs. Calculate the statistic Z based on the standard deviation and expected value. The smaller the statistic Z, the worse the mixing between the simulated data and the real data. The larger the statistic Z, the better the mixing between the simulated data and the real data. The formula for calculating the statistic Z is: （3）； In formula (3), R is the actual number of runs, and E(R) is the expected number of runs. The standard deviation of the number of runs; Metric 4 is used to evaluate the mixed effect between real and simulated datasets, including: Sample N is composed of real data and simulated data. For each sample within sample N, ... Perform a chi-square test to calculate the sample size. Calculate the chi-square test statistic based on the distances to all samples, select the k closest samples, and then calculate the chi-square test statistic. ,like If the sample is rejected, the NN Rejection Fraction value is calculated based on the number of rejected samples. The higher the NN Rejection Fraction, the worse the mixing effect between the real dataset and the simulated dataset. The lower the NN Rejection Fraction, the better the mixing effect between the real dataset and the simulated dataset. Chi-square test statistic The calculation formula is: （4）； In formula (4), and The expected values ​​of real data and simulated data are respectively, and the sample size is... The number of samples from X in the nearest neighbors is The number of samples from Y is ; Indicator 5 is used to evaluate the simulation performance between the simulated dataset and the real dataset, including: R samples are randomly selected from the simulated dataset. For each selected sample, the Euclidean distance between it and all real data is calculated, and the silhouette coefficient is calculated based on the calculated distance. According to the profile coefficient Calculate the average profile width Among them, the average profile width The average contour coefficient and average contour width for all simulated data. The closer the value is to ±1, the less the simulated data is mixed with the real data, resulting in a poor simulation effect; the closer the silhouette coefficient is to 0, the more the simulated data overlaps with the real data, resulting in a good simulation effect. Average profile width The calculation formula is: （5）； In formula (5), n The total number of data points, the silhouette coefficient. The range is [-1, 1]; Indicator 6 is used to calculate the difference in the cumulative distribution function between the simulated dataset and the real dataset, including: The simulated sample X and the real sample Y are merged to generate a coordinate sequence. Calculate the cumulative distribution difference of samples X in quadrant Q, calculate the eCDF of sample Y in quadrant Q, and calculate the maximum absolute difference of eCDF between the two samples in all quadrants. According to the maximum absolute difference Calculate standardized statistics Among them, when the standardized statistic The smaller the value, the smaller the difference between the cumulative distribution functions of the simulated sample and the real sample, and the more similar the data are. When the standardized statistic... The larger the value, the greater the difference between the cumulative distribution functions of the simulated sample and the real sample, indicating that the data are dissimilar. Standardized statistics in indicator six The calculation formula is: （6）； In formula (6), m is the size of the simulated sample X, and n is the size of the real sample Y; Indicator 7 is used to calculate the standardized statistics of the density functions of the simulated and real datasets, obtaining the differences between them, including: Based on the kernel function K and the bandwidth matrix H, the simulated sample X={ , ,..., } and the real sample Y={ , ,..., Kernel density estimation of} For kernel density estimation Integrate to obtain the test statistic T; substitute the kernel density estimate into the formula for calculating the test statistic T to obtain the estimated value of the test statistic. The mean correction term is calculated by combining the kernel function K and the bandwidth matrix H. and variance estimation Based on mean correction term and variance estimation Calculate the standardized statistic Z. The smaller the standardized statistic Z, the closer the difference in the density function of the two samples is to the mean correction term, and the difference is not significant. The larger Z is, the further the difference in the density function of the two samples is from the mean correction term, and the difference is significant. The formula for calculating the standardized statistic Z in Indicator 7 is as follows: （7）； Step 3.2: For indicators one through seven, calculate the difference measure between the real data and the simulated data on the corresponding indicators. If the final difference is within [0,1], no processing is performed; otherwise, the indicator scores from the same dataset on the same attribute are Z-score standardized and then min-max normalized to the range [0,1] to obtain the indicator scores. ; Step 3.3: Score all evaluation items under the corresponding attributes. Aggregate and obtain attribute scores ; Step 3.4: Take the mean of all attribute scores to obtain the comprehensive score of the data attributes. ; Attribute Score The calculation formula is: （8）； In formula (8), m The number of evaluation metrics applicable to this attribute. j For the first of this attribute j One evaluation indicator; Comprehensive score of data attributes The calculation formula is: （9）。 6. The comprehensive performance evaluation method for a single-cell RNA sequencing simulator according to claim 1, characterized in that, Step 3 involves a quantitative evaluation of the simulated dataset across functional dimensions, including: Establish a quality assessment index for cell population structure, and calculate the comprehensive score of cell population structure in simulated and real datasets. Differential expression analysis was performed on real and simulated data to quantitatively assess genes that maintain differential expression in functional groups. Establish evaluation indicators for differential expression analysis and calculate the overall score of batch effect.

7. The comprehensive performance evaluation method for a single-cell RNA sequencing simulator according to claim 6, characterized in that, Establish cell population structure quality assessment indicators, and calculate the comprehensive score of cell population structure in simulated and real datasets, including: The high-dimensional gene expression matrices of real and simulated data are mapped to a low-dimensional feature space to obtain several cell groups and calculate the Euclidean distance between cells. The quality of cell clustering is assessed based on the average silhouette width, a metric used to evaluate cell population structure quality. This includes: For the i-th sample point in the high-dimensional gene expression matrix, its silhouette coefficient is calculated based on the calculated Euclidean distance between cells. Based on the calculated profile coefficient Calculate the mean profile coefficient for all data points The mean of the profile coefficient A value closer to 1 indicates better separation between cell clusters; the mean silhouette coefficient... Values ​​close to negative or 0 indicate poor separation between cell clusters; Mean profile coefficient The calculation formula is: （10）； In formula (10), n The total number of data points. i For the first i One data point; The Dunn Index, a cell population structure quality assessment metric, measures the balance between the separation between cell clusters and the compactness within each cluster, including: Calculate the minimum value of the minimum distance between all pairs of cell clusters. And the maximum value of the maximum diameter among all cell clusters Based on the above parameters, the Dunn value is calculated. A higher Dunn value indicates better inter-cluster separation, higher intra-cluster compactness, and excellent clustering results; a lower Dunn value indicates overlap between clusters or loose clustering within clusters, and poor clustering effect. The formula for calculating the Dunn value is: （11）； Based on the cell population structure quality assessment metric Connectivity, the clustering effect of similar data samples is evaluated, including: For the real dataset / simulated dataset X={ } and clustering result C={ },set up For point The set of m nearest neighbors, for each point Count the number of its m nearest neighbors that do not belong to the same cluster. Calculate the mean connectivity of all data points. When the mean connectivity Time represents perfect clustering, connectivity mean The larger the value, the more dispersed the clustering effect; connectivity mean The calculation formula is: （12）； Global clustering quality is quantified based on the Calinski-Harabasz Index, a cell population structure quality assessment metric, including: Calculate the sum of squared differences within and between cell clusters, and calculate based on the sum of squared differences within and between cell clusters. The CHI value ranges from (0, +∞), and a higher CHI value indicates better clustering quality. The formula for calculating the CHI value is: （13）； In formula (13), k is the total number of cell clusters and n is the total number of cells; The sum of squared inter-cluster mean deviations is calculated by multiplying the squared distance between the center point of each cluster and the global center point by the number of cells in that cluster, and then summing the results. The sum of squares of the deviations from the mean within a cluster is calculated by squared the distances between each cell and the center of its cluster, and then summing the sums over all cells. The Davies-Bouldin Index, a cell population structure quality assessment metric, measures the average similarity between each cell cluster and other cell clusters, including: Calculate the average Euclidean distance from all cells within any cell cluster to the center point of other cell clusters, and the inter-cluster Euclidean distance between cell clusters, and then calculate... Value, when The smaller the value, the better the clustering quality; The formula for calculating the value is: （14）； In formula (14), k is the total number of cell clusters; For clusters The internal dispersion is calculated as the average Euclidean distance from all cells within the cluster to its center point. For clusters and The inter-cluster Euclidean distance; For each evaluation metric, the difference between the real data and the simulated data is calculated. If the difference metric falls within the range of [0,1], no processing is performed; otherwise, the scores of the same attribute for the cell population structure quality assessment metric from a dataset are first Z-score standardized, and then min-max normalized to the range of [0,1] to obtain the final score. score; Five cell population structure quality assessment indicators The average score is used to obtain the overall score of cell population structure. ; Overall score of cell population structure The calculation formula is: （15）。 8. The comprehensive performance evaluation method for a single-cell RNA sequencing simulator according to claim 6, characterized in that, Establish evaluation indicators for differential expression analysis and calculate the overall score of batch effects, including: The high-dimensional gene expression matrices of real and simulated data are mapped to a low-dimensional feature space to obtain several cell groups and calculate the Euclidean distance between cells. Based on differential expression analysis, the k-nearest neighbor batch effect test was used to evaluate the batch effect using the 'kBET'R package. The average silhouette width, an evaluation metric based on differential expression analysis, is used to assess the separation between cell clusters, including: For the i-th sample point in the high-dimensional gene expression matrix, its silhouette coefficient is calculated based on the calculated Euclidean distance between cells. Based on the calculated profile coefficient Calculate the mean profile coefficient for all data points The mean of the profile coefficient A value closer to 1 indicates better separation between cell clusters; the mean silhouette coefficient... Values ​​close to negative or 0 indicate poor separation between cell clusters; Local Inverse Simpson's Index, an indicator used in differential expression analysis to assess the degree of mixing of different batches of cells within a local region, includes: For each cell i in the real / simulated dataset, select the k nearest cells in the low-dimensional feature space to form the local nearest neighbor set of that cell. Calculate its local Simpson index The inverse of the Simpson index is taken and the mean LISI of all cells is calculated. The higher the mean LISI, the lower the degree of mixing between different batches of cells in the local area. The formula for calculating the LISI mean is: （16）； In formula (16), n is the total number of cells; The Adjusted Rand Index, an evaluation metric based on differential expression analysis, quantifies batch effects, including: The K-means clustering algorithm was used to group all cells and obtain cluster labels. The obtained cluster labels were compared with the actual batch labels to calculate the adjusted Rand index (ARI). The range of the adjusted Rand index (ARI) is [-1, 1]. The closer the adjusted Rand index (ARI) is to 1, the more obvious the batch effect is. The closer the ARI is to 0, the less significant the batch effect is. The formula for calculating the RAND Corporation Index (ARI) has been adjusted as follows: （17）； In formula (17), The number of samples that simultaneously belong to both cluster Ui and batch Vj; = For clustering The total number of samples in the sample; = For batch The total number of samples in the sample; n is the total number of samples; ) represents the number of combinations of selecting 2 from n samples. For each evaluation metric, the difference between the real data and the simulated data is calculated. If the difference metric falls within the range of [0,1], no processing is performed; otherwise, the scores of the same attribute for the cell population structure quality assessment metric from a dataset are first Z-score standardized, and then min-max normalized to the range of [0,1] to obtain the final score. score; Assessment of the structural quality of four cell populations The average score is used to obtain the batch effect composite score. ; Batch effect composite score The calculation formula is: （18）。 9. A comprehensive performance evaluation method for a single-cell RNA sequencing simulator according to claim 6, characterized in that, Differential expression analysis was performed on real and simulated data to quantitatively assess genes that maintain differential expression in functional groups, including: The edgeR package was used to perform differential expression analysis on real and simulated data, and pairwise comparisons were made between different cell populations in the data to identify upregulated and downregulated genes. The proportion of differentially expressed genes in corresponding cell populations in real and simulated data was calculated based on upregulated and downregulated genes. ; If there are more than two cell populations, the proportion of differentially expressed genes between the different population pairs should be considered. The scores are summed, and the final score is used as the score for the method to maintain differentially expressed genes on this dataset. Proportion of differentially expressed genes The calculation formula is: （19）； In formula (19), These are the corresponding upregulated / downregulated genes in the simulated data. These are the upregulated / downregulated genes in the original data.

10. The comprehensive performance evaluation method for a single-cell RNA sequencing simulator according to claim 1, characterized in that, Step 4 specifically includes: An online platform was built using JavaWeb and Rserve. The platform can display the scores of 24 simulation methods on 22 corresponding real datasets, including the scores of data attribute dimensions and functional dimensions, the mean and ranking of each simulation method, and supports filtering datasets by cell count for comparison. The platform can also accept externally input data dimension attribute scores, normalize the externally input data dimension attribute scores with the index scores generated by the online platform, and output a comprehensive performance evaluation score.

Images