Single cell sequencing data analysis method based on Bayesian modeling

Abstract

The invention discloses a single cell sequencing data analysis method based on Bayesian modeling, which is suitable for low expression gene identification and biological heterogeneity detection. The method comprises the following steps: S1, data acquisition and preprocessing; s2, establishing a Bayesian hierarchical model, and modeling technical noise and biological heterogeneity at the same time; s3, introducing a Poisson model and Poisson-negative binomial model combined strategy, and fitting and expressing a mean value and variance relationship by using a Gaussian radial basis function so as to accurately describe gene variation; s4, carrying out posterior inference by utilizing Markov chain Monte Carlo sampling to identify HVG and LVG genes; and S5, under the condition of lack of the spike-in gene, introducing a batch effect and a cross-batch sharing prior, and constructing a non-spike-in expansion model. Compared with a traditional method, the method has the advantages that the accuracy and robustness of expression variation recognition are improved, and the method is suitable for single cell data analysis of complex biological systems such as tumors and immune microenvironments.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of biological information, and particularly relates to a single-cell sequencing data analysis method based on Bayesian modeling, and more particularly relates to a gene expression heterogeneity detection technology based on a Bayesian hierarchical model. BACKGROUND

[0002] Single-cell RNA sequencing (scRNA-seq) technology is an emerging gene sequencing technology that can reveal intercellular heterogeneity of gene expression levels, providing new tools for cancer research, immunology research, and stem cell research. However, single-cell sequencing data is usually affected by high-level technical noise, making it difficult to identify truly differentially expressed genes in the population. Traditional methods (such as Seurat and SCran) are mostly based on fixed variance thresholds or simple linear models, which cannot effectively separate mean-dependent biological variation from technical noise, and are difficult to correct batch effects without spike-in genes, resulting in low accuracy of HVG / LVG identification.

[0003] Currently, in the analysis of single-cell sequencing data, Bayesian methods are widely used because they can combine prior information and experimental data to provide robust statistical inference. However, existing Bayesian models do not fully consider the characteristics of single-cell sequencing data, such as intercellular differences in mRNA content, measurement errors of spike-in genes, etc. Therefore, there is an urgent need for an improved Bayesian analysis method to accurately quantify intercellular gene expression heterogeneity while reducing the impact of technical noise. SUMMARY

[0004] The application provides a Bayesian analysis of single-cell sequencing data (BASiCS) based on Bayesian modeling, aiming to solve the problems of inaccurate gene expression difference analysis, low sensitivity, and high misjudgment rate caused by high technical noise interference in existing single-cell RNA sequencing technology.

[0005] To achieve the above-mentioned purpose, the technical solution adopted by the application is as follows:

[0006] S1, data collection and preprocessing: Obtain single-cell RNA sequencing data sets containing biological genes and spike-in genes (exogenous reference genes) from the Gene Expression Omnibus (GEO) of the public database NCBI. Perform quality control on the data to ensure the reliability of the data; standardize the data to construct a unified gene expression matrix, where the rows represent genes i (i = 1, 2, …, q, where q0 is the number of biological genes, and q-q0 is the number of spike-in genes), and the columns represent cells j (j = 1, 2, …, n);

[0007] S2, Bayesian hierarchical model construction: Under the framework of Bayesian Analysis of Single-Cell Sequencing (BASiCS), a Bayesian hierarchical model is established to model technical noise and biological heterogeneity, where technical noise is modeled by spike-in genes, and biological heterogeneity is described by hyperparameters;

[0008] S3, differential expression correction and residual modeling: Introduce a combination strategy of Poisson model and Poisson-negative binomial model to correct the confusion effect between mean expression and variance. Use Gaussian radial basis function to fit the relationship between expression mean and variance, and extract the residual overdispersion parameter to accurately describe gene variation;

[0009] S4, posterior inference and gene identification: Use Markov Chain Monte Carlo (MCMC) method for parameter estimation and posterior distribution inference, and use Gibbs Sampling and AM algorithm (Adaptive Metropolis) to realize joint sampling of normalization constant, technical noise parameter, gene expression mean, and overdispersion parameter. Identify highly variable genes (HVG) and lowly variable genes (LVG) through tail posterior probability, and optimize threshold selection by combining expected false discovery rate (EFDR) and expected false negative rate (EFNR). The output results include corrected gene expression mean and residual overdispersion parameter, HVG / LVG gene list and its statistical significance;

[0010] S5, spike-in-free modeling extension: In the absence of spike-in genes, introduce batch effects and cross-batch shared priors to realize stable modeling and biological effect identification under multiple batches of data.

[0011] Preferably, in S1, the collected single-cell RNA sequencing data includes biological genes and spike-in genes. The specific preprocessing method involves removing cells with fewer than 300,000 biological genes or fewer than 600,000 spike-in genes, deleting biological genes with an average count less than 1, and deleting undetected spike-in genes.

[0012] Preferably, in S2, technical variability is quantified using spike-in genes, which are unaffected by biological factors and can directly reflect technical variations in the experimental process; intercellular biological heterogeneity is quantified using the hyperparameter δ. i The description is used to represent the natural differences in gene expression among cells. In the hierarchical modeling process, modeling is first performed on the expression count of each gene in a single cell. The expression count of biological genes is regulated by cell-specific factors, while the expression count of spike-in genes is determined solely by technological factors. Subsequently, cell-specific random effects v are introduced. j The former, along with the technology variation parameter θ, is used to capture technological differences between cells that cannot be explained by experimental design, while the latter revolves around the capture efficiency normalization parameter s. j Fluctuations; the latter is related to gene expression heterogeneity among cells, as measured by the hyperparameter δ. i Quantify the level of biological heterogeneity among cells. Ultimately, analyze cell-specific factors φ. j Acquisition efficiency normalization parameter s j Hyperparameters (δ) i Prior distributions are set for key parameters such as θ, etc., to reflect the initial assumptions about each parameter before the experimental data observation, thus providing a foundation for subsequent Bayesian inference. The hierarchical model in the BASiCS framework is as follows:

[0013]

[0014] v j |s j ,θ~Gamma(1 / θ,1 / (s j θ) and ρ ij |δ i ~Gamma(1 / δ i ,1 / δ i )

[0015] Among them, X ij q0 is a random variable representing the number of genes i (i∈{1,…,n}) in cell j (j∈{1,…,n}), where q0 is the number of biological genes and q-q0 is the number of spike-in genes. μ i φ represents the average expression of a gene in a cell. j It is a cell-specific factor that corrects for differences in mRNA levels between cells. jThis represents cell-specific random effects, used to capture unexplained technical noise. θ is a technical variation parameter. j To capture efficiency normalization parameters. Random effects ρ ij (E(ρ ij |δ i )=1,Var(ρ ij |δ i )=δ i It is related to the heterogeneity of gene expression among cells.

[0016] Preferably, in S3, a Poisson-negative binomial model is introduced to correct for the confounding effect of mean and variance. The model is constructed to separate biological heterogeneity from mean dependence. Specifically, this is achieved by introducing a gene-specific residual overdispersion parameter ∈ i This eliminates the confounding effect of mean expression on heterogeneity, enabling accurate assessment of gene expression variability. The Poisson-negative binomial model in the BASiCS framework is as follows:

[0017]

[0018] in, For shape parameters, the residual over-dispersion parameter ∈ i A nonlinear regression relationship is established with the mean, specifically as follows:

[0019]

[0020] Where g1(·),…,g l (·) denotes a set of Gaussian radial basis function kernels, α0, α1, β1, ..., β l For regression coefficients, ∈ i ~t η (0,σ 2 ) represents the residual term independent of the mean, used to capture gene-specific variations that exceed the mean-variance global trend. Here, η is the degree of freedom, 0 is the location parameter, and σ... 2 The scaling parameter is used. A nonlinear model is constructed using Gaussian radial basis function kernels and regression coefficients to fit the relationship between the mean and variance, capturing complex mean-variance dependencies. The gene-specific residual overdispersion parameter ∈ i It is used to quantify the portion of gene expression variation that exceeds the mean-variance relationship, thereby eliminating the confounding effect of mean expression level on heterogeneity and making the assessment of gene expression variability more accurate.

[0021] Preferably, in S4, the identification criteria for highly variable genes are:

[0022]

[0023] Where, (φs) *The technical noise is represented by the median of the normalized constant, expressed as: γ H This is the variance contribution threshold for highly variable genes. Specifically, it's the threshold where the total variation rate of gene i expression counts in the reference cell exceeds a variance contribution threshold γ. H This indicates the existence of biological heterogeneity among cells. The BASiCS model quantifies the evidence of highly variable genes based on the upper tail posterior probability and labels these genes as HVGs. For a given evidence threshold α... H have:

[0024]

[0025] in, This represents the upper tail posterior probability of a highly variable gene. Similarly, for a given low-variable gene, the variance contribution threshold γ is used. L and evidence threshold α L Low-variability genes can be represented as:

[0026]

[0027] This represents the upper tail posterior probability of a low-variable gene.

[0028] Preferably, in S4, the posterior probability calculation is achieved through sampling using Gibbs and AM algorithms, and the Markov Chain Monte Carlo (MCMC) method is used for parameter estimation and posterior distribution inference, including the mean expression μ. i Over-dispersion parameter δ i Technical noise θ and normalization constant φ j and s j .

[0029] Preferably, in S4, the selection of thresholds is further optimized by considering the expected false detection rate (EFDR) and expected false negative rate (EFNR) for highly variable genes (HVG), specifically using the following formula:

[0030]

[0031] Wherein, if A is true, then I(A) = 1, otherwise it is 0;

[0032] The specific formula for optimizing the threshold selection of expected false detection rate and expected false negative rate for low variable gene (LVG) is as follows:

[0033]

[0034] Generally speaking Methods for detecting highly (low) variable genes require the selection of a variance contribution threshold γ. H and γL and the evidence threshold α H and α L If at a specific γ H and γ L If there is a biological motivation behind the value, it can be fixed before analysis. However, α H and α L Uncertainties related to the detection of HVG and LVG can be quantified. For γ H and γ L The fixed value can be α. H and α L Choose the optimal value. In this invention, r H =0.79, r L =0.41, evidence threshold α H =0.7925, α H =0.7650.

[0035] Preferably, S5 also includes extended modeling of spike-in-free genes, utilizing information from multiple batches of sequenced cells to estimate biological effects unaffected by technological differences. Specifically, biological effect estimation and subsequent processing are achieved through the following method: For the spike-in-free gene dataset, batch effect modeling is introduced, separating technological variation from biological heterogeneity by integrating cross-batch data. Assume the experimental design contains k independent batches, and the expression count of gene i in cell j of batch k is X. ijk Then construct the following model:

[0036] X ijk ~Poisson(μ i v jk ρ ijk ), v jk ~Gamma(1 / θ) k ,1 / (s jk θ k )),

[0037] ρ ijk ~Gamma(1 / δ i ,1 / δ i )

[0038] Where, μ i This represents the global average expression of gene i, shared across all batches, reflecting the true biological signal. jk ρ represents the random effect of technical variation in cells j within batch k. ijk θ represents the random effect of gene i on the biological heterogeneity of cell j in batch k. k Let s be the technical variation intensity parameter for batch k. jkδ is a normalization constant for cell j in batch k, used to correct for differences in cell capture efficiency within the batch. i This is a biological heterogeneity hyperparameter for gene i, shared across all batches, quantifying the biological origins of intercellular expression differences.

[0039] Compared with the prior art, the beneficial effects of this invention are as follows:

[0040] The above-mentioned optimized scheme further ensures the high adaptability and robustness of this method to different sequencing depths, experimental batches, and platform standards. It is suitable for single-cell transcriptomics research in actual clinical samples, public database data, and scenarios without spike-in genes, providing a solid data analysis foundation for subsequent disease mechanism exploration, target screening, and personalized treatment strategies. Compared with advanced methods, the advantages of this invention are: (1) This invention can simultaneously model expression mean, technical noise and biological heterogeneity, breaking through the limitation of traditional differential analysis methods that cannot distinguish between "technical fluctuations" and "real variations", thus improving the interpretability of single-cell data; (2) This invention supports modeling based on spike-in genes and extended modeling under conditions without spike-in genes, which has higher adaptability and generalizability than existing statistical frameworks that rely on exogenous controls; (3) This invention solves the problem of systematic interference of expression mean on variability estimation by introducing Gaussian radial basis functions, which is superior to traditional variance or coefficient of variation indicators and effectively reduces the misjudgment rate of low-expression genes; (4) This invention adopts the Bayesian posterior significance test method, which can flexibly set significance control strategies (EFDR, EFNR), which is more robust and biologically interpretable than the classic hypothesis testing method; (5) This invention shows higher detection sensitivity and model stability in various datasets such as real mouse embryonic stem cell heterogeneity modeling, especially in scenarios with limited sample size or significant batch bias. This invention has good scalability and openness, and can be seamlessly integrated with new data types such as spatial transcriptomics and time-series single-cell RNA sequencing data, which facilitates further integration and application of downstream tasks such as cell type identification, trajectory inference and therapeutic target screening. Attached Figure Description

[0041] Figures 1 to 6 These are the parameters μ in the BASiCS model. i δ i v j s j φ j The θ iteration trajectory and autocorrelation plot show that the posterior distribution of the parameters tends to stabilize with the number of iterations, while the autocorrelation plot on the right indicates that, except for cell-specific random effects v... j In addition, the autocorrelation coefficients of the other parameters rapidly decreased to below 0.2, which was used to verify the convergence of Gibbs sampling and the AM algorithm, indicating the reliability of the model parameter estimation.

[0042] Figure 7 It is a cell-specific factor φ j and normalization constant s j The posterior distribution feature map, where (a) is φ j (a) High posterior density interval plot, reflecting the distribution of correction parameters for differences in mRNA content among cells; the narrower the interval, the more accurate the parameter estimation; (b) shows the s in each cell. j The posterior median distribution plot reveals the heterogeneity of cell capture efficiency, and the concentrated median region suggests that most cells have similar technical stability; (c) shows the s of the first cell (s1). j Posterior distribution morphology diagram, (d) is the s of the second cell (s2). j The posterior distribution pattern verifies the model's ability to accurately characterize single-cell technology noise, with the peak position corresponding to the optimal capture efficiency estimate.

[0043] Figure 8 This quantifies the spatial distribution and global intensity of technical noise, where (a) represents the cell-specific random effect v. j (a) The nuclear density estimation plot shows the heterogeneity of technical noise among cells, with the peak concentrated in the range of 0.4-0.6, indicating that most cells are subject to similar levels of technical perturbation; (b) is the posterior distribution histogram of the global technical variation parameter θ, with a peak of about 0.49, proving that the model effectively captures the overall technical error level of the experiment and provides a benchmark for distinguishing between biological variation and technical noise.

[0044] Figure 9 It is the biological heterogeneity parameter δ i The posterior median and the mean gene expression μ i The joint distribution map, where (a) is the biological heterogeneity parameter δ i The posterior median distribution plot, δ i The posterior median varies with μ i The increase showed a trend of first rising and then falling, verifying the necessity of Gaussian radial basis function (GRBF) correction for mean dependence, indicating that the heterogeneity of highly expressed genes may be affected by differences in transcriptional regulatory mechanisms; (b) is the average gene expression μ i The distribution map shows that the peak is concentrated in the low expression region (μ). i <1), reflecting the typical characteristic of a high proportion of low-expression genes in single-cell data, providing a data foundation for residual modeling.

[0045] Figure 10 These are the detection results of HVG and LVG, where (a) is the HVG detection result, with a threshold α. H Genes with a variability score above 0.7925 were marked as highly variable, and their variability contribution index σ was used to indicate this. i The value is significantly higher than the variance contribution threshold γH =0.79, proving that the model can effectively screen genes with significant intercellular heterogeneity; (b) is the LVG detection result diagram, with threshold α L Genes with a σ value above 0.7650 are classified as low variability, and their σ i Values ​​below γ L =0.41, indicating that the expression differences of this type of gene are mainly driven by technical noise.

[0046] Figure 11 These are the estimation results of spike-in-free genes under the BASiCS model, where (A) represents the average expression parameter μ among all spike-in-free genes. i The posterior estimate of the mean expression parameter μ of the spike-in gene is... i The prior estimate highly overlaps with the results of spike-in data, demonstrating that cross-batch prior sharing effectively eliminates batch effect interference; (B) represents the overdispersion parameter δ among all non-spike-in genes. i The posterior estimate of the overdispersion parameter δ i The posterior distributions are consistent across both classes of data, indicating that the model can still accurately quantify biological heterogeneity even without external references, thus expanding the applicability of the method in real clinical samples. Detailed Implementation

[0047] It should be noted that the following detailed descriptions are illustrative and intended to provide further explanation of this application. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this application pertains.

[0048] It should be noted that the terminology used herein is for the purpose of describing particular embodiments only and is not intended to limit the exemplary embodiments according to this application. As used herein, the singular form includes the plural form unless the context clearly indicates otherwise. Furthermore, it should be understood that when the terms "comprising" and / or "including" are used in this specification, they indicate the presence of features, steps, operations, devices, components, and / or combinations thereof.

[0049] This invention provides a method for analyzing single-cell sequencing data based on Bayesian modeling, comprising the following steps:

[0050] S1. Collect raw RNA count data from single cells. Obtain raw sequencing data containing biological genes and spike-in genes from the public database NCBI Gene Expression Omnibus (GEO). After removing low-quality cells and invalid genes, normalize the data and construct an expression matrix.

[0051] S2. Introduce the BASiCS framework and construct a Bayesian hierarchical structure including cell-specific factors, technical noise parameters, and biological heterogeneity parameters. Under the BASiCS framework, jointly model biological heterogeneity and technological variability, where technological variability is estimated using spike-in genes, and biological variability is modeled using differences in gene expression among cells.

[0052] S3. To correct the systematic influence of expression mean on variability, a Gaussian radial basis function (GRBF) is introduced to establish a nonlinear relationship between expression mean and variance, thereby extracting the residual over-dispersion parameter to truly reflect the heterogeneity among genes.

[0053] S4. Parameter estimation and posterior distribution inference are performed using the Markov Chain Monte Carlo (MCMC) method. Gibbs sampling and the AM algorithm are used to jointly sample the normalization constant, technical noise parameter, gene expression mean, and overdispersion parameter. The posterior distribution is used to calculate the proportion of biological variation for each gene, and then high-variable genes (HVG) and low-variable genes (LVG) are identified based on the Bayesian significance threshold. This step allows for flexible setting of EFDR and EFNR to control the risk of false positives and false negatives.

[0054] S5. To adapt to real-world sample scenarios without spike-in genes, the model is further extended to support cross-batch modeling. By introducing batch effect factors and cross-batch shared prior mechanisms, the model achieves accurate decomposition of expression differences and stable identification of biological heterogeneity without relying on spike-in genes.

[0055] Preferably, in step S1, the raw data primarily comes from 25,914 genes from 96 publicly available mouse embryonic stem cells, and 92 spike-in genes are introduced. To ensure data quality, preprocessing is used to exclude genes with low expression rates and very weak biological relevance, remove cells with fewer than 300,000 biological genes or fewer than 600,000 spike-in genes, delete biological genes with an average count less than 1, and delete undetected spike-in genes. Finally, 7,941 genes (7,895 biological genes and 46 spike-in genes) from 41 mouse embryonic stem cells are retained, ensuring data stability and representativeness.

[0056] Preferably, in S2, the Bayesian hierarchical model is based on the Poisson noise model, and the biological expression level X ij Satisfy the following formula:

[0057] X ij |μ i ,φ j ,v j ,ρ ij ~Poisson(φj v j μ i ρ ij ), i=1,…,q0,j=1,…,n.

[0058] Where i represents a gene, j represents a cell, and q0 is the number of biological genes. μ i φ represents the average expression of a gene in a cell. j It is a cell characteristic factor that corrects for differences in mRNA levels between cells. j For random effects of technological variation, used to capture unexplained technological noise, v j |s j ,θ~Gamma(1 / θ,1 / (s j θ), where θ is the global technical variation parameter, s j This is a capture efficiency normalization constant. Random effects ρ ij (E(ρ ij |δ i )=1,Var(ρ ij |δ i )=δ i It is related to the heterogeneity of gene expression among cells, and the hyperparameter δ i Used to describe biological heterogeneity between cells. For spike-in genes, the deviation from the Poisson formula is solely due to unexplained technical noise, which is assumed here to depend on cell-specific characteristics and, for a given cell, affects the expression count of all genes in the same way, with spike-in gene expression level X. ij Satisfy the following formula:

[0059] X ij |μ i ,v j ~Poisson(v j μ i ), i = q0 + 1, ..., q, j = 1, ..., n.

[0060] Where q-q0 represents the number of spike-in genes.

[0061] Preferably, in step S3, to eliminate the interference of the expression mean on variability, a Gaussian radial basis function (GRBF) is used to apply the log(δ) method. i ) and log(μ i The relationship between ) is modeled using nonlinear regression to extract the residuals ∈ i Characterizing the true overdispersion effect after removing mean dependence:

[0062] log(δi )=f(μ i )+∈ i

[0063] Where, f(μ) i ) indicates that the average expression μ has been given. i The global trend of over-dispersion predicted (by all genes) ∈ i The excessive dispersion parameter for potential gene-specific residues is used to capture the overall trend of deviation, g1(·),…,g l (·) denotes a set of Gaussian radial basis function kernels, α0, α1, β1, ..., β l For regression coefficients, the GRBF kernel can be defined as:

[0064]

[0065] Where, m l and h l This represents the location and scale hyperparameters of the GRBF kernel.

[0066] Preferably, in step S4, based on the above model structure, the Markov Chain Monte Carlo (MCMC) method is used for parameter estimation and posterior distribution inference, and Gibbs sampling and the AM algorithm are used to estimate the parameter μ. i , θ, φ j and s j Posterior estimation is performed to calculate the proportion of biological variation for each gene, thereby identifying highly variable genes (HVGs) and low-variable genes (LVGs). Highly variable genes (HVGs) are those whose total expression variability can be largely explained by intercellular heterogeneity. A significance threshold and a variance contribution threshold γ are set. H =0.79, r L =0.41 and the evidence threshold α H =0.7925, α L =0.7650, if

[0067]

[0068] Where, σ i The contribution index to the variability of gene i is used to measure the extent to which the intercellular heterogeneity of this gene explains the total expression variation; a higher value indicates higher gene variability. (φs) * For technical noise, Represents φ in n cells j s j Take the median. The total variability in gene i expression counts in cells exceeds the variance contribution threshold γ. HWhen this indicates the presence of biological heterogeneity among cells, the gene is classified as a highly variable gene (HVG). The BASiCS model quantifies the evidence of high gene variability based on the upper tail posterior probability and labels these genes as HVGs. For a given evidence threshold α... H The following expression exists:

[0069]

[0070] in, This represents the upper tail posterior probability of highly variable genes. Similarly, for a given variance contribution threshold γ of a low-variable gene... L and evidence threshold α L Low-variable genes (LVG) can be represented as:

[0071]

[0072] This represents the upper tail posterior probability of low-variable genes. The method for detecting high (and low) variable genes in this invention requires selecting a variance contribution threshold γ. H and γ L and the evidence threshold α H and α L If γ H or γ L The specific values ​​of α have a biological basis and can be preset before analysis. However, α H and α L There are uncertainties when using it to quantify highly variable genes (HVG) and low variable genes (LVG). For fixed γ H and γ L Value, the present invention will α H and α L The optimal value is chosen when the expected false discovery rate (EFDR) equals the expected false negative rate (EFNR). Generally speaking, Therefore, the expected false discovery rate of highly variable genes and expected false negative rate They are respectively:

[0073]

[0074] Here, I(A) is an indicator function; if A is true, I(A) = 1, otherwise it is 0. Similarly, for low-variable genes, the expected false discovery rate... and expected false negative rate for

[0075]

[0076] Preferably, in S5, in the absence of spike-in gene data, it is assumed that the biological effects are shared across batches, and that technical differences will be reflected through spurious differences between cells and batches. The extended model is as follows:

[0077] X ijk |μ i ,v jk ,ρ ijk ~Poisson(μ i v jk ρ ijk )

[0078] Among them, X ijk μ represents the number of genes i (i∈{1,…,q}) in cell j (j∈{1,…,n}) of batch k (k∈{1,…,K}). i This represents the global average expression of gene i, shared across all batches, reflecting the true biological signal. jk For the random effect of technical variation in cells j in batch k, v jk |s jk ,θ~Gamma(1 / θ k ,1 / (s jk θ k )), θ k Let s be the technical variation intensity parameter for batch k. jk ρ is the normalization constant for cell j in batch k. ijk ρ represents the random effect of gene i on the biological heterogeneity of cell j in batch k. ijk |δ i ~Gamma(1 / δ i ,1 / δ i ), δ i This is a biological heterogeneity hyperparameter for gene i, shared across all batches, quantifying the biological origins of intercellular expression differences.

[0079] The technical solution of the present invention will be clearly and completely described below with reference to the embodiments. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0080] Example

[0081] This embodiment provides a single-cell sequencing data analysis method based on Bayesian modeling, including the following steps:

[0082] S1. Data Acquisition and Preprocessing.

[0083] In this embodiment, the dataset used in the experiment came from the GEO database, which contained the distribution of 25,914 genes in 96 mouse embryonic stem cells. To better identify genes exhibiting true heterogeneous expression in the cell population, 92 spike-in genes were introduced. During the quality control phase, cells with fewer than 300,000 biological genes or fewer than 600,000 spike-in genes were removed, along with biological genes with an average count less than 1 and undetected spike-in genes. Ultimately, 7,941 genes (7,895 biological genes and 46 spike-in genes) from 41 mouse embryonic stem cells were retained, completing the construction of the standardized expression matrix. This process ensured the representativeness of the analysis samples and the reliability of the results.

[0084] S2. Bayesian model construction and parameter estimation.

[0085] A BASiCS model was constructed, with all hyperparameters set to 1. The model utilizes the Markov Chain Monte Carlo (MCMC) method for parameter estimation and posterior distribution inference, employing a joint sampling approach using Gibbs sampling and the AM algorithm. The total number of iterations was 20,000 steps, the burn period was set to 10,000 steps, and the sampling interval was 10 steps. The sampling results show that the posterior trajectories of each parameter are stable. Figures 1 to 6 Parameter μ in the BASiCS model i δ i v j s j φ j The θ iteration trajectory and autocorrelation plot show that, apart from v j and φ j The autocorrelation plots of all external parameters quickly dropped to zero, indicating that the model had converged sufficiently. Figure 7 and Figure 8 The posterior inference results demonstrate the cell-specific technical noise characteristics and global technical variation level in single-cell RNA sequencing data processing. Figure 7 It is a cell-specific factor φ j and normalization constant s j Posterior inference, Figure 8 It is a cell-specific random effect v j The distribution plots of the technical noise parameter θ show that the posterior distributions of both types of parameters conform to the technical characteristics of single-cell sequencing data, indicating that the BASiCS model has reasonableness and reliability in modeling technical noise. The parameter estimation results of the BASiCS model are shown in Table 1 below:

[0086] Table 1. Parameter estimation results for the BASiCS model.

[0087]

[0088]

[0089] S3. Residual variability modeling and mean correction.

[0090] To correct for the nonlinear correlation between expression mean and variability, a Gaussian radial basis function (GRBF) was used to adjust log(δ). i ) and log(μ i The relationship between ) is modeled using nonlinear regression, such as Figure 9 As shown, the fitting function f(μ) is obtained. i Then calculate the residual ∈ for each gene. i , that is ∈ i =log(δ) i )-f(μ i The residual reflects the true biological heterogeneity after controlling for the influence of the mean, which helps to more accurately distinguish between genes with high and low variability.

[0091] S4. Identification of highly variable genes (HVG) and low variable genes (LVG).

[0092] Highly variable genes (HVGs) and lowly variable genes (LVGs) were identified using a posterior probability thresholding method. A variance contribution threshold r was set. H =0.79, r L =0.41, evidence threshold α H =0.7925, α L =0.7650. The detection results showed that a total of 133 highly variable genes (HVG) and 589 low-variable genes (LVG) were identified, verifying that technical noise contributed 28% to the total variability. Figure 10 As shown, these highly variable genes exhibited significant expression differences among samples, and some are closely related to cell cycle, stemness maintenance, and early carcinogenesis pathways, possessing potential cancer biomarker functions. Furthermore, the expected false discovery rate (EFDR) and expected false negative rate (EFDR) were both controlled within 5%, validating the stability and high reliability of the method.

[0093] S5. Visual analysis and result verification.

[0094] pass Figures 1 to 10 It displays the posterior distribution, sampling trajectory, autocorrelation of key parameters, as well as the GRBF fitting curve and classification results. Figures 1 to 6 The sampling trajectory shows good convergence. Figures 7 to 8 They revealed the heterogeneity of cell-specific noise and technological perturbation, respectively. Figure 9 Showing δ i With μ i The nonlinear relationship Figure 10The significant differences in posterior probability distributions between HVG and LVG were clearly distinguished. The combined results of various analyses validated that this method possesses high stability, accuracy, and biological interpretability on this dataset.

[0095] S6, multi-batch data modeling under conditions without spike-in genes.

[0096] This experiment used mouse embryonic stem cells without the spike-in control gene. After initial quality control, batches were standardized, and batch labels were recorded as input variables for the model. Therefore, utilizing the principle of measurement error models, technical differences were quantified through replication in the absence of a gold standard feature. This horizontal data integration method is based on an experimental design where cells in the population are randomly assigned to multiple independent experimental replicates (referred to here as "batches"). In this experimental design, the spike-in-free implementation of BASiCS is assumed to have shared biological effects across batches, and technical differences will be reflected through spurious differences. Figure 11 As shown, through the comparison of parameter estimation in the absence of spike-in gene scenarios, it is intuitively demonstrated that the model of the present invention can still accurately quantify biological heterogeneity in the absence of exogenous references, solving the problem that traditional methods cannot model batch effects under spike-in conditions, and providing a more universal technical solution for single-cell transcriptome research (especially clinical sample analysis).

[0097] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely illustrative of the principles of the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the present invention as claimed. The scope of protection of this invention is defined by the appended claims and their equivalents.

Claims

1. A method for analyzing single-cell sequencing data based on Bayesian modeling, characterized in that, Includes the following steps: S1. Data Acquisition and Preprocessing: Obtain single-cell RNA sequencing datasets containing biological genes and spike-in genes from the publicly available NCBI Gene Expression Omnibus database; perform quality control on the data to ensure its reliability; standardize the data to construct a unified gene expression matrix, where rows represent genes i (i = 1, 2, ..., q, where q0 is the number of biological genes and q-q0 is the number of spike-in genes), and columns represent cells j (j = 1, 2, ..., n); S2. Bayesian hierarchical model construction: Under the framework of Bayesian Analysis of Single-Cell Sequencing, a Bayesian hierarchical model is established to model both technical noise and biological heterogeneity. Technical noise is modeled using spike-in genes, and biological heterogeneity is described using hyperparameters. S3. Differential Expression Correction and Residual Modeling: A strategy combining the Poisson model and the Poisson-negative binomial model is introduced to correct the confusion effect between mean expression and variance; the relationship between the expression mean and variance is fitted using the Gaussian radial basis function, and the residual over-dispersion parameter is extracted to accurately characterize gene variation; S4. Posterior Inference and Gene Identification: Parameter estimation and posterior distribution inference are performed using the Markov chain Monte Carlo method. Gibbs sampling and the AM algorithm are used to jointly sample the normalization constant, technical noise parameter, mean gene expression, and overdispersion parameter. High-variable genes (HVG) and low-variable genes (LVG) are identified through tail posterior probability, and threshold selection is optimized by combining the expected false detection rate and expected false negative rate. Output results include the corrected mean gene expression and residual overdispersion parameter, a list of HVG / LVG genes, and their statistical significance. S5. Spike-in-free modeling extension: In the absence of spike-in genes, batch effects and cross-batch shared priors are introduced to achieve stable modeling and biological effect identification under multi-batch data.

2. The single-cell sequencing data analysis method based on Bayesian modeling according to claim 1, characterized in that, In S1, the collected single-cell RNA sequencing data includes biological genes and spike-in genes; the specific preprocessing method is to remove cells with fewer than 300,000 biological genes or fewer than 600,000 spike-in genes, delete biological genes with an average count of less than 1 and undetected spike-in genes.

3. The single-cell sequencing data analysis method based on Bayesian modeling according to claim 1, characterized in that, In S2, technical variability is quantified using spike-in genes, which are unaffected by biological factors and directly reflect technical variations in the experimental procedure; intercellular biological heterogeneity is quantified using the hyperparameter δ. i The description is used to represent the natural differences in gene expression among cells; in the hierarchical modeling process, modeling is carried out on the expression count of each gene in a single cell. The expression count of biological genes is regulated by cell-specific factors, while the expression count of spike-in genes is determined only by technical factors; then, cell-specific random effects v are introduced. j The former, along with the technology variation parameter θ, is used to capture technological differences between cells that cannot be explained by experimental design, while the latter revolves around the capture efficiency normalization parameter s. j Fluctuations; the latter is related to gene expression heterogeneity among cells, as measured by the hyperparameter δ. i Quantifying the level of biological heterogeneity among cells; quantifying cell-specific factors φ j Acquisition efficiency normalization parameter s j Hyperparameters (δ) i The prior distribution of key parameters (θ) is set to reflect the initial assumptions about each parameter before the experimental data observation, providing a foundation for subsequent Bayesian inference; the hierarchical model in the BASiCS framework is as follows: v j |s j , θ ∼ Gamma(1 / θ, 1 / (sjθ)) and ρ ij |δ i ∼ Gamma(1 / δ i , 1 / δ i ) Among them, X ij A random variable represents the number of genes i (i∈{1,…,q}) in cell j (j∈{1,…,n}), where q0 is the number of biological genes and q-q0 is the number of spike-in genes; μ i φ represents the average expression of a gene in a cell; j It is a cell-specific factor that corrects for differences in mRNA levels between cells; v j Cell-specific random effects are used to capture unexplained technical noise; θ is a technical variation parameter; s j To normalize the capture efficiency parameter; random effect ρ ij (E(ρ ij |δ i )=1,Var(ρ ij |δ i )=δ i It is related to the heterogeneity of gene expression among cells.

4. The single-cell sequencing data analysis method based on Bayesian modeling according to claim 1, characterized in that, In S3, a Poisson-negative binomial model is introduced to correct for the confounding effects of mean and variance. The model is constructed to separate biological heterogeneity from mean dependence, specifically by introducing a gene-specific residual overdispersion parameter ∈ i To eliminate the confounding effect of mean expression on heterogeneity and achieve accurate assessment of gene expression variability; the Poisson-negative binomial model in the BASiCS framework is as follows: in, For shape parameters, the residual over-dispersion parameter ∈ i A nonlinear regression relationship is established with the mean, specifically as follows: Where g1(·),…,g l (·) denotes a set of Gaussian radial basis function kernels, α0, α1, β1, ..., β l For regression coefficients, ∈ i ~t η (0,σ 2 ) represents the residual term independent of the mean, used to capture gene-specific variations that exceed the mean-variance global trend; where η is the degree of freedom, 0 is the location parameter, and σ 2 The scaling parameter is used; a nonlinear model is constructed using Gaussian radial basis function kernels and regression coefficients to fit the relationship between the mean and variance, capturing complex mean-variance dependencies; the gene-specific residual overdispersion parameter ∈ i It is used to quantify the portion of gene expression variation that exceeds the mean-variance relationship, thereby eliminating the confounding effect of mean expression level on heterogeneity and making the assessment of gene expression variability more accurate.

5. The single-cell sequencing data analysis method based on Bayesian modeling according to claim 1, characterized in that, In S4, the criteria for identifying highly variable genes are: Where, (φs) * The technical noise is represented by the median of the normalized constant, expressed as: γ H The variance contribution threshold for highly variable genes; that is, when the total variation rate of gene i expression count in the reference cell exceeds a variance contribution threshold γ. H This indicates the existence of biological heterogeneity among cells; the BASiCS model quantifies the evidence of highly variable genes based on the upper tail posterior probability and labels these genes as HVG, for a given evidence threshold α. H have: in, Represents the upper tail posterior probability of a highly variable gene; γ is the variance contribution threshold for a given low-variable gene. L and evidence threshold α L Low-variable genes are represented as: This represents the upper tail posterior probability of a low-variable gene.

6. The single-cell sequencing data analysis method based on Bayesian modeling according to claim 1, characterized in that, In S4, the posterior probability is calculated through sampling using the Gibbs and AM algorithms. The Markov chain Monte Carlo method is employed for parameter estimation and posterior distribution inference, including the mean expression μ. i Over-dispersion parameter δ i Technical noise θ and normalization constant φ j and s j .

7. The single-cell sequencing data analysis method based on Bayesian modeling according to claim 1, characterized in that, S4 also includes optimizing the threshold selection based on the expected false detection rate and expected false negative rate of the highly variable gene HVG, with the specific formula as follows: Wherein, if A is true, then I(A) = 1, otherwise it is 0; For optimizing the threshold selection of the expected false detection rate and expected false negative rate for low variable gene LVG, the specific formula is as follows: Methods for detecting highly variable and low-variable genes require the selection of a variance contribution threshold γ. H and γ L and the evidence threshold α H and α L If in a specific γ H and γ L If a value has a biological motivation, it can be fixed before analysis; α H and α L Uncertainties related to HVG and LVG can be quantified and detected; for γ H and γ L The fixed value can be α. H and α L Choose the optimal value; r H =0.79, r L =0.41, evidence threshold α H =0.7925, α H =0.7650.

8. The single-cell sequencing data analysis method based on Bayesian modeling according to claim 1, characterized in that, S5 also includes extended modeling for spike-in-free genes, utilizing information from multiple batches of sequenced cells to estimate biological effects unaffected by technological differences. Specifically, biological effect estimation and subsequent processing are achieved through the following methods: For spike-in-free gene datasets, batch effect modeling is introduced, separating technological variation from biological heterogeneity by integrating cross-batch data; assuming the experimental design contains k independent batches, and the expression count of gene i in cells j of batch k is X. ijk Then construct the following model: X ijk ~Poisson(m i v jk r ijk ),v jk ~Gamma(1 / θ k ,1 / (s jk i k )), r ijk ~Gamma(1 / δ i ,1 / d i ) Where, μ i The global average expression of gene i, shared across all batches, reflects the true biological signal; v jk For the random effects of technical variation in cells j within batch k; ρ ijk θ represents the random effect of gene i on the biological heterogeneity of cell j in batch k; k Let s be the technical variation intensity parameter for batch k. jk δ is a normalization constant for cell j in batch k, used to correct for differences in cell capture efficiency within the batch; i This is a biological heterogeneity hyperparameter for gene i, shared across all batches, quantifying the biological origins of intercellular expression differences.

Images