Single-cell rna sequencing pathway activity scoring method, system and device based on pagerank and probability model

Abstract

The application discloses a single-cell RNA sequencing pathway activity scoring method, system and device based on PageRank and a probability model, and the method comprises the following steps: acquiring a gene set of a specified pathway and an interaction relationship between genes, and constructing a gene interaction directed network graph; the importance score of each gene is calculated by using a PageRank algorithm, and normalization processing is performed to obtain the weight value of each gene; a statistical model is fitted to the expression count of each gene, model parameters are estimated, and the conditional mean function probability value of each gene in each cell is calculated based on the model parameters; the weight value and the conditional mean function probability value are weighted and summed to calculate the activity score of each pathway in each cell. The method of the application can more accurately quantify the activity level of biological pathways in single cells, reduce the limitations in weight differentiation, noise processing and network structure utilization, and provide a more reliable calculation tool for single-cell function analysis.

Description

Technical Field

[0001] This invention relates to the fields of bioinformatics and computational biology, and in particular to a method, system, and apparatus for scoring single-cell RNA sequencing pathway activity based on PageRank and probability models. Background Technology

[0002] The rapid development of single-cell RNA sequencing technology has enabled researchers to precisely analyze cellular heterogeneity and gene expression dynamics at single-cell resolution, providing unprecedented insights into development, disease, and life processes. In such studies, to deeply reveal the biological functions and states of cells, researchers often need to assess the activity levels of specific biological pathways, such as signal transduction, metabolic, or regulatory pathways, in single cells, thereby transforming gene expression data into biologically meaningful pathway characteristics.

[0003] Currently, existing pathway activity scoring methods suffer from the following main shortcomings, which limit their accuracy and reliability in single-cell data analysis: 1. Simple Averaging Method: This method typically uses a simple arithmetic average or summation of the expression levels of all genes within a pathway to represent pathway activity. This method ignores the differences in the contributions of different genes to pathway function and fails to reflect the regulatory hierarchy and synergistic relationships between genes, resulting in overly coarse scoring results that do not accurately reflect the biological activity of the pathway under specific cellular conditions.

[0004] 2. Equal-weighting method: Many existing methods treat all genes in a pathway as having equal importance weight. However, in actual biological processes, pathways often contain core driver genes, key regulatory genes, and execution or helper genes, and their influence on the overall pathway activity is not equal. Equal-weighting may underestimate the role of core genes and may amplify the noise from fluctuations in the expression of non-key genes, thereby reducing the discriminative power and biological interpretability of the scoring results.

[0005] 3. Lack of statistical modeling for single-cell data: Single-cell RNA-seq data exhibit significant technical variability, such as frequent dropout events and batch effects. Most existing pathway scoring methods fail to incorporate these data characteristics into their statistical frameworks, directly applying models based on population RNA-seq assumptions. This results in scoring methods that are sensitive to technical noise, lack stability and reproducibility, and are particularly prone to bias in low-expression or sparse genes.

[0006] 4. Ignore the network topology within the pathway: Biological pathways are essentially functional networks formed by the interactions between genes and proteins. Most existing methods do not integrate known gene interaction information in pathways, such as protein-protein interactions and regulatory relationships. Therefore, they cannot distinguish between genes at key network nodes, such as highly connected regulatory factors, and genes at the periphery. They also fail to utilize network topology information to weight gene contributions, thus limiting the ability to accurately infer pathway activity at the system level.

[0007] To address the shortcomings of existing technologies, this invention proposes a pathway activity scoring method that integrates the PageRank algorithm and probabilistic modeling. This method aims to more accurately and robustly quantify the activity level of biological pathways in single cells, overcoming the limitations of the aforementioned methods in terms of weight differentiation, noise processing, and network structure utilization, and providing a more reliable computational tool for single-cell functional analysis. Summary of the Invention

[0008] This invention aims to address at least one of the technical problems existing in the prior art. To this end, this invention proposes a single-cell RNA sequencing pathway activity scoring method based on PageRank and a probabilistic model. This method can more accurately quantify the activity level of biological pathways in single cells, reducing limitations in weight differentiation, noise processing, and network structure utilization, and providing a more reliable computational tool for single-cell functional analysis.

[0009] This invention also proposes a single-cell RNA sequencing pathway activity scoring system based on PageRank and a probability model.

[0010] This invention also proposes a single-cell RNA sequencing pathway activity scoring device based on PageRank and a probability model.

[0011] The single-cell RNA sequencing pathway activity scoring method based on PageRank and probability models according to the first aspect of the present invention includes the following steps: S1. Pathway network construction and gene weight calculation: Obtain the gene set and gene interaction relationships of the specified pathway, construct a directed gene interaction network graph; use the PageRank algorithm to calculate the importance score of each gene in the pathway, and perform normalization to obtain the weight value of each gene. S2. Probability Model Fitting and Conditional Mean Function Probability Calculation: Fit a statistical model to the expression count of each gene, estimate the model parameters, and calculate the conditional mean function probability value of each gene in each cell based on the model parameters. The conditional mean function probability is used to characterize the activation probability of the gene under the current expression count. S3. Pathway activity score calculation: The weighted sum of the weighted values ​​and the probability values ​​of the conditional mean function is calculated to obtain the activity score of each pathway in each cell.

[0012] Optionally, the PageRank algorithm iteratively calculates the importance score of a gene according to the following formula: in, Indicates gene PageRank score, Indicates gene V j Pagerank score, This represents the total number of genes in the pathway. The damping coefficient is... Indicates pointing to a gene The entire set of genes, Indicates gene The degree of departure.

[0013] Optionally, the statistical model is a negative binomial distribution model or a zero-inflated negative binomial distribution model. When the proportion of zero values ​​in the data is lower than a preset threshold, the negative binomial distribution model is selected; when the proportion of zero values ​​in the data is higher than the preset threshold, the zero-inflated negative binomial distribution model is selected.

[0014] Optionally, when the NB model is selected, the probability mass function of the NB model is: When the ZINB model is selected, the probability mass function of the ZINB model is: Where Y represents gene expression count, μ represents conditional mean, θ represents discrete parameter, Γ represents Gamma function, and π represents zero inflation probability.

[0015] Optionally, the regression formula of the statistical model is: in, The intercept term of gene g is represented. This represents the batch effect coefficient of gene g. This indicates the batch to which cell c belongs. This represents the offset of cell c.

[0016] Optionally, the probability of the conditional mean function is calculated as follows: in, This represents the CMF probability of gene g in cell c. The cumulative distribution function (CDF) represents the negative binomial distribution. This indicates the adjusted count. This represents the expected expression level of gene g in cell c. The parameter representing the dispersion of gene g.

[0017] Optionally, the formula for calculating the pathway activity score is: in, This indicates the activity fraction of pathway p in cell c. This represents the set of all genes contained in pathway p. This represents the normalized PageRank weight of gene g in pathway p. This represents the CMF probability of gene g in cell c.

[0018] Optionally, a parallel computing method is used for model fitting for each gene. When the model fails to fit, the gene is automatically skipped, and only the successfully fitted genes and their renormalized weights are used to calculate the pathway score.

[0019] According to a second aspect of the present invention, a single-cell RNA sequencing pathway activity scoring system based on PageRank and a probability model is used to perform the above-described method, the system comprising: The data input module is used to receive and validate single-cell RNA sequencing data. The pathway weight calculation module is used to obtain pathway information from the KEGG database and calculate gene weights; The model fitting module is used to fit statistical models in parallel and calculate the probability of the conditional mean function. The pathway scoring module is used to calculate the pathway activity scoring matrix; The results output module is used to output the scoring results and perform quality control.

[0020] According to a third aspect of the present invention, a single-cell RNA sequencing pathway activity scoring device based on PageRank and a probability model includes a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the steps of the above-described method.

[0021] The present invention has at least the following beneficial effects: Improved accuracy: By using PageRank weights to differentiate gene importance, the shortcomings of simple averaging are avoided, and the scoring results more accurately reflect the true activity of the pathway.

[0022] Enhanced robustness: Statistical modeling methods can effectively handle technical noise, dropout events, and batch effects in single-cell data, resulting in more robust and reliable scoring results.

[0023] The biological significance is clear: the weight allocation based on network topology is more in line with biological principles, the core genes contribute more to the score, and the results are more biologically interpretable.

[0024] High computational efficiency: Supports multi-core parallel computing, enabling efficient processing of large-scale single-cell datasets.

[0025] Highly flexible: It supports multiple pathway categories and can be flexibly configured according to research needs, making it suitable for different biological research scenarios.

[0026] Good scalability: The algorithm framework can be extended to other pathway databases (such as Reactome, GO, etc.) and other types of omics data.

[0027] Additional aspects and advantages of the invention will be set forth in part in the description which follows, and in part will be obvious from the description, or may be learned by practice of the invention. Attached Figure Description

[0028] To more clearly illustrate the technical solutions of the embodiments of this application, the accompanying drawings used in the embodiments will be briefly introduced below. It should be understood that the following drawings only show some embodiments of this application and should not be regarded as a limitation of the scope. For those skilled in the art, other related drawings can be obtained based on these drawings without creative effort.

[0029] Figure 1 This is a flowchart of a single-cell RNA sequencing pathway activity scoring method based on PageRank and probability models according to an embodiment of the present invention. Figure 2 This is a schematic diagram of the experimental steps using a pathway activity scoring method according to an embodiment of the present invention; Figure 3 The results are a comparison of the distribution of activity scores on the KEGG pathway “Nicotinate and nicotinamide metabolism” after normalization from 0 to 1 using a pathway activity scoring method according to an embodiment of the present invention and other methods. Figure 4 This is a schematic diagram illustrating the scoring and intergroup comparison results of the pathway activity scoring method according to an embodiment of the present invention and other methods for the KEGG pathway “alpha-linolenic acid metabolism”. Figure 5 This is a schematic diagram illustrating the results of scoring and comparing the KEGG pathway “Nicotinate and nicotinamide metabolism” using a pathway activity scoring method according to an embodiment of the present invention and other methods. Detailed Implementation

[0030] Embodiments of the present invention are described in detail below, examples of which are illustrated in the accompanying drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to the accompanying drawings are exemplary and intended to explain the present invention, and should not be construed as limiting the present invention.

[0031] The following disclosure provides numerous different embodiments or examples for implementing various structures of the invention. To simplify the disclosure, specific examples of components and arrangements are described below. These are merely examples and are not intended to limit the invention. Furthermore, reference numerals and / or letters may be repeated in different examples. Such repetition is for simplification and clarity and does not in itself indicate a relationship between the various embodiments and / or arrangements discussed. Additionally, examples of various specific processes and materials are provided in this invention; however, those skilled in the art will recognize the applicability of other processes and / or the use of other materials.

[0032] The single-cell RNA sequencing pathway activity scoring method based on PageRank and probability models according to the first aspect of the present invention includes the following steps: S1. Pathway network construction and gene weight calculation: Obtain the gene set and gene interaction relationships of the specified pathway, construct a directed gene interaction network graph; use the PageRank algorithm to calculate the importance score of each gene in the pathway, and perform normalization to obtain the weight value of each gene. S2. Probability Model Fitting and Conditional Mean Function Probability Calculation: Fit a statistical model to the expression count of each gene, estimate the model parameters, and calculate the conditional mean function probability value of each gene in each cell based on the model parameters. The conditional mean function probability is used to characterize the activation probability of the gene under the current expression count. S3. Pathway activity score calculation: The weighted sum of the weighted values ​​and the probability values ​​of the conditional mean function is calculated to obtain the activity score of each pathway in each cell.

[0033] In some embodiments, step S1 retrieves all human pathways from the KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway database, including pathway IDs, gene lists, and interaction relationships, and obtains their hierarchical information. The hierarchical information is a three-level classification structure (A / B / C), divided into major categories, subcategories, and specific pathways. The gene expression count matrix obtained by single-cell sequencing technology is input as the first parameter of the function, and the KEGG pathway category to be analyzed, the statistical model, and the number of parallel computation cores are selected.

[0034] In some embodiments, step S1 further includes preprocessing and filtering the gene expression data. The input data format is: a gene × cell expression count matrix, with row names being gene names and column names being cell identifiers.

[0035] In some embodiments, step S1 constructs a directed gene interaction network graph. Gene information for a specified set of pathways is obtained from the KEGG (Kyoto Encyclopedia of Genes and Genomes) database. Specifically, the steps include: pathway classification selection: one or more of the following categories can be selected: "Metabolism", "Genetic Information Processing", "Environmental Information Processing", "Cellular Processes", "Organismal Systems", "Human Diseases", and "Drug Development".

[0036] Furthermore, for each selected pathway, a directed gene interaction network graph G=(V,E) is constructed, where: V represents the set of nodes, each node corresponding to a gene in the pathway; E represents the set of edges, each directed edge representing a regulatory or metabolic interaction between genes. The direction of the edge indicates the direction of regulation (e.g., gene A activates gene B). The number of pathway categories is approximately 1-7, and the minimum number of genes in each pathway is ≥3. The network graph is constructed using the pathwayGraph function based on the KEGG database.

[0037] In some embodiments, step S1 uses the PageRank algorithm to calculate the importance score of each gene in the pathway and performs normalization to obtain the weight value of each gene.

[0038] The PageRank algorithm was used to calculate the importance score of each gene in the pathway: The PageRank algorithm was applied to the directed network graph of gene interactions in each pathway to calculate the importance weight of each gene in the pathway. For the pathway network graph G, the PageRank algorithm iteratively calculated the importance score of each gene: in: Indicates gene PageRank score, Indicates gene V j Pagerank score, This represents the total number of genes in the pathway. The damping coefficient has a value ranging from 0.85 to 0.95, with 0.85 being preferred. Indicates pointing to a gene The set of all genes, i.e. Incoming neighbor nodes, Indicates gene out-degree (i.e. from) (Number of edges pointing to other genes) The gene weights within each pathway are normalized so that the sum of the weights is 1: in, Indicates gene Normalized weights. Parameter control range: damping coefficient d: 0.80-0.99, preferably 0.85; iterative convergence threshold: Maximum number of iterations: 50-200 times, with 100 times being the optimal number.

[0039] In some embodiments, step S1 is followed by preprocessing of the input single-cell RNA-seq expression matrix. Only genes contained in the KEGG pathway are retained, and low-expression genes expressed in fewer than five cells are filtered out to reduce computation and improve model stability. ENTREZIDs in the KEGG database are converted to gene symbols (SYMBOL) to handle one-to-many mappings. Preprocessing transforms the original data matrix into the final input matrix.

[0040] In some embodiments, step S2 fits a statistical model to the expression count of each gene, estimates the model parameters, and calculates the conditional mean function probability value of each gene in each cell based on the model parameters. The conditional mean function probability is used to characterize the activation probability of the gene under the current expression count. The statistical model is either a negative binomial distribution model or a zero-inflated negative binomial distribution model. When the proportion of zero values ​​in the data is lower than a preset threshold, the negative binomial distribution model is selected; when the proportion of zero values ​​in the data is higher than the preset threshold, the zero-inflated negative binomial distribution model is selected.

[0041] When the NB model is selected, the probability mass function of the NB model is: When the ZINB model is selected, the probability mass function of the ZINB model is: Where Y represents gene expression count, μ represents conditional mean, θ represents discrete parameter, Γ represents Gamma function, and π represents zero inflation probability.

[0042] For example, a statistical model is fitted to the expression count of each gene, using a negative binomial distribution or a zero-inflated negative binomial distribution model to fit the gene expression data. For genes Expression count in cell c The NB model assumes: in: This represents the expected expression level of gene g in cell c (mean parameter). This represents the dispersion parameter of gene g, which controls the variance of the distribution.

[0043] The mean parameter is modeled using a logarithmic link function: in: The intercept term of gene g is represented. This represents the batch effect coefficient of gene g. This indicates the batch to which cell c belongs. This represents the offset of cell c, which is usually the total sequencing depth of the cell or the normalized library size.

[0044] For sparse data with a large number of zero values, use the ZINB model: in: This represents the zero inflation probability of gene g. This represents the probability mass of the negative binomial distribution at point k.

[0045] Furthermore, the zero-inflation probability is modeled using the logit link function: In model fitting, the glmmTMB package is used for model fitting, and the BFGS optimization algorithm is adopted. The maximum number of iterations is in the range of 20-100, preferably 30-50. Multi-core parallelism is supported, with the number of cores ncores=2-32, preferably 4-8 cores.

[0046] In some embodiments, a parallel computation method is used for model fitting for each gene. When the model fails to fit, the gene is automatically skipped, and only the successfully fitted genes and their renormalized weights are used to calculate the pathway score.

[0047] Further, set the parameter control range. Choose either the NB or ZINB model based on data sparsity, and set maxit (maximum number of iterations) and ncores (number of parallel cores). Perform a convergence check: gradient norm < Or relative change .

[0048] In some embodiments, step S2 calculates the conditional mean function probability value of each gene in each cell based on the model parameters: based on the fitted model parameters, the activation probability of each gene in each cell is calculated, i.e., the conditional mean function (CMF) probability.

[0049] Counting the expression of gene g in cell c The CMF probability is defined as: in: This represents the CMF probability of gene g in cell c. The cumulative distribution function (CDF) represents the negative binomial distribution. This indicates the adjusted count (for) In this situation, (set to 0) This represents the expected expression level of gene g in cell c. The parameter representing the dispersion of gene g.

[0050] Special case handling: when Force setting This indicates that the gene is not activated in this cell. When a NaN is encountered in the calculation, it is replaced with 0.

[0051] In one embodiment, step S3 involves weighting and summing the weight values ​​and the probability values ​​of the conditional mean function to calculate the activity score of each pathway in each cell.

[0052] Formula for calculating pathway activity fraction: in: This indicates the activity fraction of pathway p in cell c. This represents the set of all genes contained in pathway p. This represents the normalized PageRank weight of gene g in pathway p. This represents the CMF probability of gene g in cell c.

[0053] Furthermore, for each pathway, ensure that the gene weights used are normalized: Output data in matrix format, with row names as KEGG pathway IDs and column names as cell IDs. Scores close to 1 indicate a high score, representing highly active pathways in the cell; scores close to 0 indicate a low score, representing inactive or unactivated pathways; and scores between 0.3 and 0.7 indicate a moderate score, representing partial activation.

[0054] Furthermore, the output data can be further processed through standardization, dimensionality reduction visualization, cluster analysis, and other follow-up steps. The algorithm can provide feedback on model convergence, data quality, parameter validity, and batch effects.

[0055] According to a second aspect of the present invention, a single-cell RNA sequencing pathway activity scoring system based on PageRank and a probability model is used to perform the above-described method, the system comprising: The data input module is used to receive and validate single-cell RNA sequencing data. The pathway weight calculation module is used to obtain pathway information from the KEGG database and calculate gene weights; The model fitting module is used to fit statistical models in parallel and calculate the probability of the conditional mean function. The pathway scoring module is used to calculate the pathway activity scoring matrix; The results output module is used to output the scoring results and perform quality control.

[0056] According to a third aspect of the present invention, a single-cell RNA sequencing pathway activity scoring device based on PageRank and a probability model includes a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the steps of the above-described method.

[0057] The following specific examples demonstrate the application of the single-cell RNA sequencing pathway activity scoring method based on PageRank and a probability model of this invention. In these examples, the method is named the ProbScorer algorithm.

[0058] Example 1 like Figure 2 As shown, the single-cell RNA sequencing pathway activity scoring method based on PageRank and probability models of this invention, i.e., the ProbScorer algorithm, includes the following steps before scoring the pathways in a single-cell dataset: Preliminary Step 1: Data Quality Control First, import the single-cell dataset into RStudio and perform quality control. The data we used comes from healthy individuals, and the dataset information is shown in the table below:

[0059] After annotating the cells in the peripheral blood single-cell transcriptome data with immune cell types, the dataset contains multiple immune cells and related cell types, specifically including: CD4. + T cell subsets (NaiveCD4) + T, MemoryCD4 + T, Th17, Th2, Treg), CD8 + T cell subsets (NaiveCD8) + T、MemoryCD8 + T, EffectCD8 + T cells, γδT cells, NK cells; B cell subsets (Naive B, Memory B, Plasma); monocytes and myeloid cell subsets (CD14 Monocyte, CD16CD14 Monocyte, C1QAMonocyte, MDSC, ProMyeloid, Neutrophil); as well as dendritic cells (cDC, pDC), megakaryocytes, mast cells, and cycling cells in an active cell cycle state. These cell types collectively constitute the main immune cell composition in peripheral blood and are used for subsequent analysis.

[0060] The purpose of data quality control is to remove low-quality cells that may affect the analysis results. This process begins by calculating the percentage of mitochondrial gene expression in each cell; cells with a higher percentage of mitochondrial genes typically represent poorer cell quality or are already dead. We will use a violin plot to visually represent the distribution of cell quality by showing three quality control indicators: the number of genes detected in each cell, the total number of RNA molecules in each cell (UMI count), and the percentage of mitochondrial genes.

[0061] Next, two scatter plots were created: total RNA count vs. mitochondrial ratio, and total RNA count vs. number of detected genes. The former was used to check for the presence of low-quality cells with a high mitochondrial ratio, while the latter was used to detect the correlation between total RNA count and number of detected genes (which should normally be positively correlated). Based on these plots, cells of poor quality were filtered out. Specific filtering criteria included: removing cells with too few detected genes (which may be empty droplets), removing cells with too many detected genes (which may be double cells), and removing cells with an excessively high proportion of mitochondrial genes (which may be dead cells). These steps help preserve high-quality cell data for subsequent analysis.

[0062] Preliminary Step 2: Data Standardization After quality control was completed, the data was standardized using logarithmic normalization. The specific method was as follows: first, the raw cell count for each cell was divided by the total cell count for that cell, then multiplied by 10,000 to obtain the standardized data per 10,000 molecules. The normalization formula is: Adding 1 is to avoid taking the logarithm of 0 and ensure accurate calculation. This process eliminates the differences between different cells caused by varying sequencing depths, allowing for comparison of expression levels across different cells.

[0063] After completing the above preliminary steps, the ProbScorer algorithm is used to score the pathways.

[0064] First, the UMI count matrix for each cell was extracted, and then the data was modeled using the ZINB (zero-inflated negative binomial distribution) statistical model. The ZINB model takes into account zero values ​​in single-cell data, which come from two sources: biological zero values ​​and technical zero values ​​(dropout). This model can model these two types of zero values ​​separately and estimate pathway activity more accurately, rather than relying solely on the mean or median.

[0065] During this process, sequencing depth correction is also performed. The offset term in the ZINB model acts similarly to a normalization factor, ensuring that differences in sequencing depth do not affect pathway activity scores. The mathematical principle is as follows: in, The pathway activity calculated in the model, This is a factor affecting sequencing depth. In this way, the ProbScorer algorithm can eliminate technical noise and ensure that pathway activity reflects the true biological differences in cells.

[0066] The ProbScorer algorithm uses the gene set from the KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway database as input. KEGG is a widely used biological database containing a large number of biological pathways, metabolic pathways, and signaling pathways, suitable for various biological research. The KEGG gene set is constructed based on experimental data and biological knowledge bases, covering multiple fields such as cell metabolism, immune responses, and cancer-related pathways. Next, the parameters of the ProbScorer algorithm need to be set. The main parameters here include: 1. Maximum number of iterations: To ensure algorithm convergence, we set a maximum number of iterations (50). If the model does not converge within this number of iterations, the algorithm will stop and report a non-convergence warning.

[0067] 2. Number of CPU cores for parallel computing: Due to the large size of the single-cell dataset and the significant computational load, we will choose parallel computing to accelerate the analysis process. Based on the size of the dataset, we select 4 CPU cores for parallel processing to improve computational efficiency.

[0068] By setting these parameters, ProbScorer can efficiently process large-scale single-cell data and calculate the activity score of each cell in each pathway. These scores can ultimately provide accurate input data for subsequent analyses (such as T-tests, differential analyses, etc.), thereby enabling in-depth analysis of the biological significance of the pathways.

[0069] Example 2 Based on the method in Example 1, we sought pathways for cell combinations with significant scores from the ProbScorer algorithm. In this process, we compared the scores from multiple algorithms to explore ProbScorer's performance in pathway scoring. First, we normalized the pathway scores so that each pathway's score fell between 0 and 1. This step ensured that scores from different algorithms could be compared under the same standard. We merged the scores from seven algorithms into a matrix, where each row represented the score of a pathway and each column represented the score of an algorithm. This resulted in a matrix containing the scores from the seven algorithms. After obtaining the normalized scoring matrix, we then filtered based on the differences between different cell combinations. We selected two typical cell types (denoted as cell A and cell B) for comparison and obtained cell indices for these two cell types. The selected cell combinations typically represent biologically significant cell types, and their differences can reveal changes in the activity of specific pathways.

[0070] To determine the differences in pathway activity among these cell combinations, we used a t-test to compare the scores of these pathways. Specifically, we performed statistical analysis on the scores of each pathway between cell A and cell B, and calculated the p-value for each pathway. A p-value less than 0.01 indicated a significant difference in pathway activity between the two cell combinations. Based on the results of the t-test, we categorized the differences in each pathway. Specifically, a p-value greater than or equal to 0.01 was marked as 0, indicating no significant difference; a p-value less than 0.01 was marked as +1 or -1. +1 indicates that the p-value of cell A is greater than that of cell B, and -1 indicates that the p-value of cell B is greater than that of cell A. Through this process, we obtained the significance of the differences in each pathway and categorized them. Finally, through the above t-test analysis, we screened out the pathways with significant differences in the ProbScorer algorithm. It is worth noting that the differences in these pathways were only significant in the ProbScorer scoring results, and did not show the same differences in other algorithms. This result indicates that the ProbScorer algorithm has unique advantages in the evaluation of certain pathways.

[0071] The table below shows the results obtained after performing the above operations on a single-cell dataset:

[0072] Example 3: Based on the method in Example 2, after performing a t-test on the pathway scoring results of the single-cell dataset, further visualization processing is performed. First, pathways showing only significant differences in the t-test results are extracted. These pathways were scored across different cell combinations using the ProbScorer algorithm and showed significant differences after the t-test. Through this screening step, we can determine which pathways exhibit strong biological differences between specific cell combinations. Selected pathways include alpha-linolenic acid metabolism, nicotinamide metabolism, etc.

[0073] Next, based on the T-test matrix from Example 2, cell combinations with significant differences were selected for further analysis. For example, in the alpha-linolenic acid metabolism pathway, the ProbScorer algorithm showed a significant difference in scoring results between cell A and cell B; therefore, these two cell combinations (cell A and cell B) were selected for subsequent analysis. Box plots were used to visualize the results for each algorithm. The pathway score distribution for each algorithm was presented using box plots. Box plots help us intuitively understand the score distribution of different algorithms on the same pathway, thereby assessing the performance differences between the algorithms. The box plots of all seven algorithms were combined into a single chart to directly compare the scoring results of different algorithms on the same pathway. This makes it clearer to see the superiority of the ProbScorer algorithm on certain pathways, especially its accuracy and consistency in pathways such as alpha-linolenic acid metabolism and nicotinamide metabolism. The visualization results are as follows: Figure 3 , 4 As shown in Figure 5.

[0074] Figure 3 This paper presents the results of a comparison of the distribution of activity scores (normalized from 0 to 1) of CD8_naive cells and CD8_EM cells in the KEGG pathway "Nicotinate and nicotinamide metabolism" in the Monaco peripheral blood single-cell RNA sequencing dataset, using the ProbScorer algorithm of this invention and six comparison algorithms: VISION, AUCell, ssGSEA, GSVA, UCell, and singscore. The horizontal axis represents cell type, and the vertical axis represents pathway score; the p-value of the two-sided t-test between the two groups is labeled above each subplot. Figure 3As can be seen, the difference in scores between the two groups obtained by ProbScorer was statistically significant (p=0.0082), and the overall score of the CD8_naive group was higher than that of the CD8_EM group, with obvious box separation between the groups. Although VISION (p=0.021) and singscore (p=0.047) were marginally significant at the α=0.05 level, their significance was weaker than that of ProbScorer. The p values ​​of AUCell, ssGSEA, GSVA, and UCell were 0.95, 0.052, 0.13, and 0.52, respectively, and the distributions of the two groups overlapped significantly. The above results are consistent with the sign matrix of the T-test obtained in Example 2 for the comparison between CD8_naive and CD8_EM: under this pathway and the criterion of p<0.01, only ProbScorer was marked as significant. This demonstrates that, in cases where single-cell data exhibits sparse expression and zero inflation, the pathway scoring method of this invention, which combines PageRank gene weights with the probability of the conditional mean function of a negative binomial / zero-inflation negative binomial model, is more effective in identifying CD8 than conventional methods based on gene set enrichment or expression ranking. + Differences in metabolic pathway activity among T cell subsets.

[0075] Figure 4 This shows the scoring and inter-group comparison results of seven algorithms for the KEGG pathway "alpha-linolenic acid metabolism" under the same Monaco dataset and the same CD8_naive and CD8_EM cell comparison conditions. The plotting method is the same as... Figure 3 Same. With Figure 3 The nicotinic acid / nicotinamide metabolic pathways shown are different, and the activity score of this pathway in CD8_EM cells is higher than that in CD8_naive cells. Figure 4 As can be seen, ProbScorer detected the highest significance of differences between groups (p=0.0072), and its separation of the score distributions between the two groups was better than the other six algorithms. The p-values ​​of VISION, AUCell, ssGSEA, GSVA, UCell, and singscore were 0.71, 0.25, 0.84, 0.82, 0.18, and 0.31, respectively, none of which reached the significance level of α=0.05, indicating a high degree of overlap between the groups. This result is consistent with the "alpha-Linolenicacidmetabolism" entry in the representative summary table of Example 2: in the comparison between CD8_EM and CD8_naive, only ProbScorer was marked as +1 in the T-test sign matrix of the seven algorithms, while the other algorithms were marked as 0. This shows that ProbScorer is consistent with the differences in cell subpopulations and metabolic pathways. Figure 3Under different conditions, ProbScorer can still reliably detect biologically significant differences in pathway activity, while other comparison algorithms are susceptible to technical noise in small peripheral blood single-cell subsets and are unable to reflect the true inter-group differences caused by synergistic changes in core pathway genes.

[0076] Figure 5 This shows the results of scoring and inter-group comparisons of the KEGG pathway "Nicotinate and nicotinamide metabolism" using seven algorithms, extending cell contrast to B_naive cells and B_SM cells in the Monaco dataset; the mapping rules are as follows. Figure 3 , Figure 4 Same. By Figure 5 As can be seen, the ProbScorer score for the B_naive group was significantly higher than that for the B_SM group, with a p-value of 0.0046, the lowest among the seven subplots. VISION (p=0.04) and singscore (p=0.044) were only weakly significant at the α=0.05 level. The p-values ​​for ssGSEA, GSVA, AUCell, and UCell were 0.10, 0.083, 0.39, and 0.40, respectively, all below the conventional significance threshold, indicating significant overlap in group distributions. These visualization results are consistent with the T-test matrix and representative pathway summary table obtained in Example 2 regarding the comparison between B_naive and B_SM: under this pathway and the p<0.01 criterion, only ProbScorer was marked as a significant positive difference. Figures 3 to 5 It can be seen that the ProbScorer algorithm of this invention is effective in CD8. + The method provides more significant statistical test results and clearer intergroup score separation with a lower risk of false positives in T cell subsets and B cell subsets, different metabolic pathways, and different group differences. This verifies the applicability of the pathway activity scoring method of this invention in various immune cell types and KEGG metabolic pathway scenarios, as well as its superiority over existing technologies.

[0077] Based on these boxplot results, we can conclude that the ProbScorer algorithm outperforms other algorithms on multiple pathways, especially in the specific case of cell A and cell B combinations. These intuitive visualizations demonstrate the ProbScorer algorithm's advantage in handling zero-inflation data, particularly in pathway activity scoring within a biological context.

[0078] In the description of this specification, the references to terms such as "one embodiment," "some embodiments," "example," "specific example," or "some examples," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of the present invention. In this specification, the illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples. Moreover, without contradiction, those skilled in the art can combine and integrate the different embodiments or examples described in this specification, as well as the features of different embodiments or examples.

[0079] Although embodiments of the invention have been shown and described, those skilled in the art will understand that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the claims and their equivalents.

Claims

1. A single-cell RNA sequencing pathway activity scoring method based on PageRank and a probability model, characterized in that, Includes the following steps: S1. Pathway network construction and gene weight calculation: Obtain the gene set and gene interaction relationships of the specified pathway, construct a directed gene interaction network graph; use the PageRank algorithm to calculate the importance score of each gene in the pathway, and perform normalization to obtain the weight value of each gene. S2. Probability Model Fitting and Conditional Mean Function Probability Calculation: Fit a statistical model to the expression count of each gene, estimate the model parameters, and calculate the conditional mean function probability value of each gene in each cell based on the model parameters. The conditional mean function probability is used to characterize the activation probability of the gene under the current expression count. S3. Pathway activity score calculation: The weighted sum of the weighted values ​​and the probability values ​​of the conditional mean function is calculated to obtain the activity score of each pathway in each cell.

2. The method according to claim 1, characterized in that, The PageRank algorithm iteratively calculates the importance score of a gene according to the following formula: in, Indicates gene PageRank score, Indicates gene V j Pagerank score, This represents the total number of genes in the pathway. The damping coefficient is... Indicates pointing to a gene The entire set of genes, Indicates gene The degree of departure.

3. The method according to claim 1, characterized in that, The statistical model is either a negative binomial distribution model or a zero-inflated negative binomial distribution model. When the proportion of zero values ​​in the data is lower than a preset threshold, the negative binomial distribution model is selected; when the proportion of zero values ​​in the data is higher than the preset threshold, the zero-inflated negative binomial distribution model is selected.

4. The method according to claim 3, characterized in that, When the NB model is selected, the probability mass function of the NB model is: When the ZINB model is selected, the probability mass function of the ZINB model is: Where Y represents gene expression count, μ represents conditional mean, θ represents discrete parameter, Γ represents Gamma function, and π represents zero inflation probability.

5. The method according to claim 3, characterized in that, The regression formula for the statistical model is: in, The intercept term of gene g is represented. This represents the batch effect coefficient of gene g. This indicates the batch to which cell c belongs. This represents the offset of cell c.

6. The method according to claim 1, characterized in that, The probability of the conditional mean function is calculated as follows: in, This represents the CMF probability of gene g in cell c. The cumulative distribution function (CDF) represents the negative binomial distribution. This indicates the adjusted count. This represents the expected expression level of gene g in cell c. The parameter representing the dispersion of gene g.

7. The method according to claim 1, characterized in that, The formula for calculating the pathway activity score is as follows: in, This indicates the activity fraction of pathway p in cell c. This represents the set of all genes contained in pathway p. This represents the normalized PageRank weight of gene g in pathway p. This represents the CMF probability of gene g in cell c.

8. The method according to claim 1, characterized in that, The model fitting for each gene is performed using a parallel computing method. When the model fails to fit, the gene is automatically skipped, and only the successfully fitted genes and their renormalized weights are used to calculate the pathway score.

9. A single-cell RNA sequencing pathway activity scoring system based on PageRank and a probability model, used to perform the method according to any one of claims 1 to 8, characterized in that, The system includes: The data input module is used to receive and validate single-cell RNA sequencing data. The pathway weight calculation module is used to obtain pathway information from the KEGG database and calculate gene weights; The model fitting module is used to fit statistical models in parallel and calculate the probability of the conditional mean function. The pathway scoring module is used to calculate the pathway activity scoring matrix; The results output module is used to output the scoring results and perform quality control.

10. A single-cell RNA sequencing pathway activity scoring device based on PageRank and a probability model, characterized in that, It includes a memory and a processor, the memory storing a computer program, and the processor executing the computer program to implement the steps of the method according to any one of claims 1 to 8.

Images