A single-cell transcriptional program screening method based on donor stability

Abstract

The application discloses a single-cell transcriptional program screening method based on donor stability, and the method comprises the following steps: obtaining single-cell transcriptional expression data set, donor identification information and pathological staging information, establishing the corresponding relationship between cells and donors and screening a donor subset; extracting a target cell set in the donor subset, calculating the continuous state coordinates of each cell based on a quasi-time sequence algorithm, and constructing a cell continuous state representation; modeling the gene expression dynamics based on the continuous state representation to obtain a trend vector, clustering the genes according to the trend similarity, and forming a candidate transcriptional program set; performing donor resampling in different stages and layers for the candidate transcriptional program set, aggregating the program values in the donor layer, establishing a fitting relationship with the pathological staging information to obtain an evaluation result, calculating the effect stability index, the direction consistency index and the donor coverage index of cross-resampling, constructing a donor layer scoring index, and evaluating the correlation and the distinguishing ability in different source verification data sets according to the donor layer.

Description

Technical Field

[0001] This invention relates to the field of bioinformatics and medical data processing technology, specifically to a method for screening transcription programs and constructing donor layer scores based on donor layer resampling stability constraints for multi-donor single-cell / monocyte transcriptome data, under the background of pseudo-temporal continuous state modeling. Background Technology

[0002] The development of chronic diseases such as liver fibrosis exhibits a continuous evolutionary pattern. However, the pathological staging system commonly used in clinical practice usually expresses this continuous process at discrete levels, making it difficult to fully analyze subtle transcriptional changes near stage boundaries, early initiation events, and key states in the intermediate transition phase. With the development of single-cell or single-nuclear transcriptome sequencing technology, researchers can characterize cellular state profiles at cellular resolution and combine pseudo-temporal or trajectory inference methods to model continuous dynamic processes. This allows them to attempt to screen dynamic gene sets or transcriptional programs related to disease progression to explain the continuous activation patterns of specific cell types in the fibrotic process.

[0003] However, in real-world clinical cohorts, single-cell transcriptional expression datasets often originate from multiple donors. These donors differ in etiological composition, tissue sampling and processing procedures, sequencing depth, cellular composition, and microenvironment. These donor-level differences mean that the signal from candidate transcription programs may be contributed by a few donors or exhibit inconsistent effect directions among different donors. Consequently, while these programs may be detectable in discovery datasets, consistent conclusions may not be reached in cohorts from different sources or in repeated experiments, thus limiting the usability of candidate transcription programs for mechanistic explanation and biomarker construction.

[0004] Meanwhile, the basic observation unit of single-cell data is the cell, while disease staging or clinical phenotype is usually labeled at the donor level. If statistical inferences about staging correlations or differences are made directly at the cell level, it is easy to treat a large number of cells in the same donor as independent samples, introducing spurious duplication bias. This results in an overestimation of the effective sample size in statistical inference and an effect size assessment that deviates from the true level at the donor level. This type of bias is more likely to occur in scenarios with a large number of cells and a relatively limited number of donors, further exacerbating the problem of insufficient cross-donor consistency of candidate transcription programs.

[0005] To alleviate the above problems, existing technologies have developed expression aggregation analysis approaches based on donor units. For example, expressions from the same donor and the same cell population are aggregated to form pseudo-batch data, and the donor is used as an independent statistical unit to conduct differential analysis, thereby reducing the impact of pseudo-repetition. However, such methods are mostly designed for discrete grouping and are difficult to couple naturally with continuous variables obtained from pseudo-time series or trajectory inference. Furthermore, they rely heavily on single-fit or single-test results and lack an objective quantitative mechanism for constraining and screening candidate transcription programs in multi-donor differential scenarios. It is difficult to simultaneously characterize the effect stability, effect direction consistency, and calculable effect size coverage of candidate transcription programs under different donor extraction conditions, and it is also difficult to form a unified basis for screening and ranking.

[0006] Therefore, it is still necessary to provide a single-cell transcription program screening method based on donor stability in the context of continuous trajectory modeling. This method involves resampling the donor set and aggregating program values ​​at the donor layer. Then, based on the effect sequences obtained across resampling, reproducibility indicators such as effect stability, directional consistency, and donor coverage are calculated. Based on these indicators, a unified scoring criterion for ranking and screening candidate programs is formed, thereby outputting a core transcription program set that can be used for subsequent scoring construction and cross-source evaluation, thus improving the reproducibility and cross-cohort generalization ability of transcription programs. Summary of the Invention

[0007] The purpose of this invention is to provide a single-cell transcription program screening method based on donor stability to address the shortcomings of the prior art.

[0008] To achieve the above objectives, the present invention provides the following technical solution: a single-cell transcription program screening method based on donor stability, comprising:

[0009] Step S101: Obtain single-cell transcriptional expression dataset, donor identification information and pathological stage information, establish the correspondence between cells and donors based on donor identification information, and screen donors with complete stage annotations to form a donor subset;

[0010] Step S102: Extract the target cell set based on the donor subset and construct a continuous state representation of each cell;

[0011] Step S103: Model gene expression dynamics based on continuous state representation, and cluster genes with similar dynamic characteristics to obtain a set of candidate transcription programs;

[0012] Step S104: Perform staged and stratified donor resampling on the candidate transcription program set. In the donor layer, express and aggregate the candidate transcription programs to obtain the quantification value of the donor layer program, and establish a smooth fitting relationship with the pathological stage information. Calculate the effect stability index, direction consistency index, and donor coverage index on the results of multiple resampling, and combine the three indices to calculate the SASI score of the candidate transcription programs for ranking.

[0013] Step S105: Select the core transcription program set according to the sorting results, and construct the donor layer scoring index from the core transcription program set; when external validation is required, obtain a validation single-cell transcription expression dataset from a different source than the single-cell transcription expression dataset, and calculate the donor layer scoring index according to the donor layer caliber while keeping the core transcription program set and its weights unchanged, so as to evaluate the correlation and discrimination ability between the donor layer scoring index and pathological staging information.

[0014] Furthermore, methods for filtering donors with complete phase annotations to form donor subsets include:

[0015] Obtain the stage-related expression dataset and extract the pathological stage label corresponding to each donor identifier based on the stage-related expression dataset; perform missing pathological stage label filtering on the pathological stage label corresponding to each donor identifier, and retain donor identifiers with no missing pathological stage labels; summarize the retained donor identifiers to form a donor subset, and extract the corresponding cell records from the stage-related expression dataset based on the donor subset.

[0016] Furthermore, methods for extracting target cell sets based on donor subsets and constructing continuous state representations of each cell include:

[0017] Obtain the sub-expression dataset corresponding to the donor subset, extract the cell type annotation field in the sub-expression dataset and compare it with the preset target cell type label set, and retain the cell records whose cell type annotation belongs to the preset target cell type label set to form the target cell set;

[0018] Gene expression values ​​of the target cell set are converted into normalized expression values ​​according to a preset normalization rule, and feature gene sets are selected according to a preset feature screening rule to construct a feature expression matrix;

[0019] The distance between cells is calculated based on the feature expression matrix and adjacency relationships are formed according to a preset number of nearest neighbors. The proportion of non-target cells is calculated based on a preset set of non-target marker genes and compared with a preset proportion threshold to remove cell records and obtain the processed target cell set. The corresponding adjacency relationships are regenerated based on the processed target cell set.

[0020] A random walk transition matrix is ​​constructed based on adjacency relationships. The initial feature score is calculated and compared with a preset initial threshold to determine the initial candidate. The initial cell is determined according to a preset selection rule. The diffusion distance of each cell relative to the initial cell is calculated and the diffusion distance or its order is used as a continuous state representation. The continuous state representation is written back to the cell record corresponding to the processed target cell set to form a continuous state association expression dataset.

[0021] Furthermore, methods for modeling gene expression dynamics based on continuous state representation include:

[0022] Obtain a continuous state association expression dataset, obtain cell identifiers, continuous state representations and gene expression values, and establish gene-by-gene modeling inputs with continuous state representations as independent variables and gene expression values ​​as dependent variables;

[0023] In the continuous state associated expression dataset, the number of cells with non-zero gene expression values ​​is counted. The cell count is compared with the preset detection cell count threshold to retain genes that meet the threshold condition. The first preset number of genes are selected by sorting by cell count to form a candidate modeling gene set.

[0024] For each gene in the candidate modeling gene set, a smooth fitting model with continuous state representation as input is constructed to output the fitted expression curve. Based on the fitted expression curve, dynamic statistics are calculated, and a preset number of genes are selected according to the dynamic statistics to form a dynamic gene set.

[0025] A uniform continuous state grid is generated within the range of values ​​represented by the continuous state. The dynamic gene set is input into the smooth fitting model, and the predicted value sequence is output as a trend vector on the uniform continuous state grid. The predicted value sequence is then standardized to meet the preset mean and variance constraints.

[0026] Furthermore, methods for obtaining a set of candidate transcription programs include:

[0027] Obtain the trend vector of each gene on a uniform continuous state grid, and use the trend vector as the dynamic feature vector of the gene.

[0028] The dynamic distance is calculated based on the point-by-point difference of the acquisition trend vectors of any two genes on the uniform continuous state grid according to the preset aggregation rules, and a preset number of nearest neighbors is selected for each gene in ascending order of the acquisition dynamic distance to obtain the set of nearest neighbor genes.

[0029] A gene adjacency graph is constructed based on the set of nearest neighbor genes obtained for each gene, where nodes are gene identifiers, edge connections are determined by the set of nearest neighbor genes obtained, and edge weights are obtained by the dynamic distance obtained according to a preset weight conversion rule;

[0030] The obtained gene adjacency graph is subjected to graph partitioning and clustering. Multiple gene clusters are output under the constraints of preset resolution parameters, and each gene cluster is defined as a candidate transcription program to obtain a set of candidate transcription programs.

[0031] Furthermore, methods based on the donor set obtained through resampling include:

[0032] Obtain a subset of donors and their pathological staging information, and divide the obtained subset of donors into multiple staging layers based on the obtained pathological staging information;

[0033] Set the number of resampling iterations and the extraction rules for each iteration. In each iteration, extract donors according to the extraction rules in each acquisition stage and merge them to obtain the donor set obtained by resampling in that iteration.

[0034] Furthermore, methods for summarizing the cellular expression of candidate transcription programs within the same donor according to preset aggregation rules to obtain donor layer program quantification values ​​include:

[0035] For any candidate transcription program, obtain the gene set corresponding to the candidate transcription program from the candidate transcription program set, and obtain the gene expression values, cell identifiers and donor identifiers corresponding to the obtained gene set from the continuous state associated expression dataset;

[0036] For any donor identifier in the donor set obtained by resampling, all cell records corresponding to the donor identifier are extracted, and the gene expression values ​​of the obtained gene set in the obtained donor are summarized according to the preset aggregation rules to obtain the donor layer program quantification value of the obtained donor identifier on the obtained candidate transcription program.

[0037] Furthermore, methods for calculating the SASI score of candidate transcription programs include:

[0038] Within the donor set obtained by resampling in each iteration, a fitting relationship is established using smooth splines with pathological staging information as the independent variable and donor layer procedural quantification value as the dependent variable, and the evaluation result is output. The evaluation result includes at least the predicted value of the fitting curve at the staging extreme value.

[0039] Based on the evaluation results, the effect size is calculated and the effect direction is recorded. The effect size is the predicted difference between the highest and lowest stages of the fitted curve. The effect size sequence and effect direction are summarized across all iterations, and the effect stability index, direction consistency index, and donor coverage index are calculated. The indexes are combined to calculate the SASI score of the candidate transcription program, which is used to rank the candidate transcription programs.

[0040] Optionally, the core transcription program set is determined based on the ranking results of the SASI scores.

[0041] Furthermore, methods for ranking candidate transcription programs and selecting a core set of transcription programs based on SASI scores include:

[0042] Obtain the candidate transcription program set and the corresponding SASI score for each candidate transcription program. Sort the programs by SASI score from high to low and select the core transcription program set according to the preset output quantity.

[0043] Obtain the donor layer program quantization value and summarize the donor layer program quantization value of each core transcription program in the core transcription program set according to the donor identifier, forming the core program quantization vector corresponding to the donor identifier.

[0044] Based on the global effect obtained from the fitting relationship in the discovery queue, the weight value is determined for each core transcription program in the core transcription program set and associated with the core transcription program identifier to form a weight set, which is then fixed and stored.

[0045] For any donor identifier, the donor layer scoring index is obtained by combining the core program quantization vector and the obtained weight set according to a preset weighted aggregation rule.

[0046] Furthermore, methods for assessing the correlation and discriminative power between donor layer scoring indicators and pathological staging information include:

[0047] Obtain a validation single-cell transcription expression dataset from a different source than the single-cell transcription expression dataset, and simultaneously obtain the donor identification information and pathological stage information corresponding to the validation single-cell transcription expression dataset, so that the cell identifier and the donor identifier are associated and mapped to the pathological stage information;

[0048] The core transcription program set and its weight set that have been identified and solidified in the discovery dataset are invoked, and the core transcription program is expressed and aggregated according to the donor layer in the validation single-cell transcription expression dataset to form a core program quantization vector.

[0049] While keeping the weight set unchanged, calculate the donor layer score index for each donor;

[0050] The correlation statistics between donor layer scoring indicators and pathological staging information were calculated under the donor layer caliber, and the area under the curve statistics used to distinguish different staging states were calculated to output the correlation and discrimination ability assessment results.

[0051] The technical effects and advantages of the single-cell transcription program screening method based on donor stability provided by this invention are as follows:

[0052] This method employs an analysis workflow that uses the donor as the smallest statistical unit to organize a single-cell transcriptional expression dataset. It extracts a target cell set within a subset of the donor cells, constructs a continuous state representation for each cell, and uses this continuous state representation as an independent variable to model gene expression dynamics, obtaining a trend vector on a uniform continuous state grid. Furthermore, it defines dynamic distances using the trend vectors and clusters them to obtain a set of candidate transcription programs. This method explicitly maps the continuous changes in disease progression to comparable continuous coordinates, allowing dynamic features to be characterized by a sequence of predicted values ​​from a continuous grid. This avoids the problem of difficulty in analyzing subtle changes near stage boundaries and transitional states caused by relying solely on discrete pathological staging groups, thus improving the ability to characterize continuous dynamic processes.

[0053] In scenarios with multiple donor differences, this invention aggregates the cell expression of candidate transcription programs within the same donor according to donor-layer aggregation rules into donor-layer program quantification values. A smooth fitting relationship is established between pathological staging information and these donor-layer program quantification values. Predicted values ​​at staging extremes are output, and effect size and effect direction are calculated. Simultaneously, donor subsets are divided into staging layers according to pathological staging information, and staging-layered donor resampling is performed. A controlled-perturbation donor set sequence is formed through multiple iterations, extending the single-step fitting into a cross-iteration effect size sequence. Based on the effect size sequence and effect direction, a SASI score is formed by combining effect stability, directional consistency, and donor coverage indices. The donor is used as the smallest statistical unit to rank candidate transcription programs, thereby constraining the dispersion, directional consistency, and computable coverage of candidate transcription programs under different donor extraction conditions with computable indices. This suppresses accidental conclusions caused by contributions from a few donors or directional inconsistencies between donors. Furthermore, donor-layer aggregation reduces spurious duplication bias introduced by treating multiple cell records within the same donor as independent samples, making effect size assessments more closely aligned with the donor layer where the pathological staging is located.

[0054] In terms of output and cross-source assessment, a core set of transcription programs is selected based on the SASI score. Weights are assigned to the core transcription programs based on the global effect obtained from the fitting relationship and stored in a fixed manner. Then, the quantified vectors of the core programs corresponding to the donor identifiers are weighted and summarized according to the weight set to construct a donor-level scoring index. This compresses the results of multiple programs into a single continuous value at the donor granularity, facilitating direct assessment of the correlation with pathological staging information and the ability to distinguish late stages at the donor level. When external validation is required, a validation single-cell transcription expression dataset from a different source than the single-cell transcription expression dataset is obtained, along with the corresponding donor identifier information and pathological staging information. The donor-level scoring index is calculated according to the same aggregation and weighting rules, while keeping the weight set unchanged. Conclusions are output by comparing the correlation statistics based on donor ranking consistency and the curve area statistics based on threshold scanning with preset thresholds, respectively. This allows the same core set of transcription programs and the same weight set to be reused under different donor composition conditions, maintaining consistent calculation methods and forming a quantifiable cross-source assessment basis, thus improving the verifiability and cross-cohort generalization ability of the obtained transcription program set and scoring index. Attached Figure Description

[0055] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

[0056] Figure 1 This is a schematic diagram of a single-cell transcriptional screening method based on donor stability according to the present invention.

[0057] Figure 2 A schematic diagram of the quantization framework for the continuous dynamics of HSC and the procedure transition at the donor level;

[0058] Figure 3 A schematic diagram of a transcription program screening framework based on donor stability constraints;

[0059] Figure 4 A schematic diagram of the human liver mononuclear transcriptome across fibrosis stages;

[0060] Figure 5 A schematic diagram illustrating the reconstruction of the continuous activation trajectory of HSCs and the construction of a dynamic transcription program;

[0061] Figure 6 A schematic diagram illustrating the screening effectiveness and core program dynamics and functional characteristics of SASI. Detailed Implementation

[0062] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, 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.

[0063] Example 1

[0064] like Figures 1 to 6 As shown, the single-cell transcription program screening method based on donor stability in this embodiment uses... Figure 1 The overall process shown is the main thread.

[0065] Figure 4 The diagram shows a human liver mononuclear transcriptome atlas across fibrosis stages and its donor and stage annotations (A: Harmony batch-corrected integrated snRNA-seq dataset UMAP; B: Classical marker gene dot plot for cell type annotation; C: UMAP colored by fibrosis stage; D: Cell type composition of different fibrosis stages), corresponding to step S101 for obtaining single-cell transcriptional expression datasets, donor identification information, and pathological stage information, and establishing the correspondence between cells and donors.

[0066] Figure 2 The framework for HSC continuous dynamics and donor-level program transition quantification is shown (specifically, the continuous state trajectory of hepatic stellate cells is constructed based on DiffusionPseudotime, program activity changes are calculated along pseudotime bins and major transition windows are identified, and the proportion of HSCs within the window is quantified at the donor level to characterize the stage-related continuous dynamic progression). This corresponds to the target cell set extraction and continuous state characterization construction in step S102, and the aggregation of candidate transcription program values ​​under the donor layer caliber to form donor layer program quantification values ​​in step S104.

[0067] Figure 5The process of reconstructing the continuous activation trajectory of HSCs and constructing dynamic transcription programs is illustrated (A: UMAP embedding of HSC subsets, colored according to diffusion pseudo-time DPT; B: Violin plot showing the distribution of different HSC subsets along the pseudo-time axis; C: Violin plot showing the donor level score distribution of ECM remodeling and TGF-β signaling related modules at different fibrosis stages; D: Heatmap of dynamic gene smoothing expression trends predicted by tradeSeq, Z-score normalized and sorted by trend similarity; E: Average normalized trend curve of candidate transcription programs along the pseudo-time axis, with dashed lines representing major state transition windows), corresponding to the generation of continuous state representation in step S102 and the dynamic modeling of gene expression based on continuous state representation and the formation of candidate transcription program sets in step S103.

[0068] Figure 3 The framework for screening transcription programs based on donor stability constraints is shown (specifically, based on the program activity input at the donor level, the stage-related effects of candidate programs are estimated multiple times through stage-stratified resampling and generalized additive model fitting, and their reproducibility is comprehensively quantified from three dimensions: effect stability, directional consistency, and coverage, and finally a SASI score for program ranking is generated), corresponding to the stage-stratified donor resampling, donor layer fitting evaluation, and the calculation and screening of effect stability index, directional consistency index, and donor coverage index in step S104.

[0069] Figure 6 The results show the effectiveness of SASI screening and the dynamic and functional characteristics of the core procedures (A: Ranking distribution of candidate procedure SASI scores; B: Distribution of candidate procedures in three dimensions: effect stability, directional consistency, and donor coverage, with the dot color representing the SASI score; C: Comparison of pseudo-temporal dynamic trajectories of donors with representative high SASI and low SASI procedures; thin lines represent individual donors, and thick lines represent the population mean; D: GO and KEGG functional enrichment analysis of the core procedures identified by SASI; E: Heatmap showing the dynamic expression of representative genes of the core procedures along the pseudo-time axis), corresponding to the quantitative representation of the stability screening effect in step S104, and the results of step S105, which outputs the core transcription procedure set based on SASI scores, constructs donor layer scoring indicators, and evaluates the donor layer in an external validation dataset.

[0070] Please see Figure 1 As shown, this embodiment provides a single-cell transcription program screening method based on donor stability, including:

[0071] Step S101: Obtain single-cell transcriptional expression datasets, donor identification information, and pathological staging information. Establish a correspondence between cells and donors based on the donor identification information, and select donors with complete staging annotations to form a donor subset (e.g., ...). Figure 4As shown, the single-cell transcriptional expression dataset is a mononuclear transcriptome map spanning fibrosis stages, carrying donor identifiers and pathological stage annotations to support subsequent summarization and comparison with the donor as the smallest statistical unit.

[0072] Methods for establishing a correspondence between cells and donors based on donor identification information include:

[0073] Obtain a single-cell transcriptional expression dataset, which includes at least cell identifiers, gene identifiers, and corresponding expression values. Simultaneously obtain donor identifier information and pathological stage information corresponding to the single-cell transcriptional expression dataset. The donor identifier information includes at least the correspondence between the donor identifier and the sample source field, and the pathological stage information includes at least the correspondence between the donor identifier and the pathological stage label.

[0074] Donor writing is performed on single-cell transcriptional expression datasets based on donor identifier information: when the cell records of the single-cell transcriptional expression dataset carry a sample source field, the donor identifier is written to each cell record according to the sample source field-donor identifier mapping; when the single-cell transcriptional expression dataset directly carries a donor field, the donor identifier is written according to the donor field; after the writing is completed, the associated expression dataset is obtained, which is a set of expression data including cell identifier and donor identifier.

[0075] A donor identifier-pathological stage label mapping table is constructed based on pathological stage information. The pathological stage label is written into all cell records corresponding to the same donor identifier in the associated expression dataset according to the donor identifier, so as to obtain the stage-associated expression dataset. The purpose of this step is to establish the correspondence between cells and donors in the same dataset, so that subsequent processing can be directly summarized and calculated on a donor-by-donor basis.

[0076] For ease of understanding, the following example is provided:

[0077] Acquire data and three types of information. Obtain single-cell transcriptional expression dataset A. This dataset A includes 6 cell records, labeled cell 1, cell 2, cell 3, cell 4, cell 5, and cell 6. Each cell record includes a gene identifier and its corresponding expression value. For example, gene A has an expression value of 12 in cell 1, gene B has an expression value of 3 in cell 1, gene A has an expression value of 8 in cell 2, and gene B has an expression value of 4 in cell 2. Simultaneously acquire donor identifier information B corresponding to dataset A. Donor identifier information B includes the correspondence between the donor identifier and the sample source field. For example, sample A1 corresponds to donor A, and sample B1 corresponds to donor B. Simultaneously acquire pathological staging information C corresponding to dataset A. Pathological staging information C includes the correspondence between the donor identifier and the pathological stage label. For example, donor A corresponds to stage two, and donor B corresponds to stage three.

[0078] Based on donor identification information, perform donor writing and form an associated representation dataset:

[0079] Scenario 1: Cell records in dataset A carry a sample source field. Assume the sample source field for cells 1, 2, and 3 is sample A1, and for cells 4, 5, and 6 it is sample B1. Based on the correspondence between the sample source field and the donor identifier in donor identifier information B, a donor identifier is written to each cell record: donor A for cell 1, donor A for cell 2, and donor A for cell 3; donor B for cell 4, donor B for cell 5, and donor B for cell 6. After writing, the associated expression dataset D is obtained. Associative expression dataset D still includes cell identifiers, gene identifiers, and expression values, while adding a new donor identifier field. For example, the donor identifier for cell 1 is donor A, and its gene A expression value of 12 and gene B expression value of 3 are retained.

[0080] Scenario 2: Cell records in dataset A directly carry the donor field. Assuming the donor field for cells 1 to 3 is donor A, and the donor field for cells 4 to 6 is donor B, then the donor identifier is directly written to each cell record according to the donor field, and the associated expression dataset D is obtained.

[0081] Write pathological stage labels and form a stage-related expression dataset:

[0082] A mapping table between donor identifiers and pathological stage labels is constructed based on pathological staging information C. The mapping table includes stage two corresponding to donor A and stage three corresponding to donor B. The pathological stage labels are written into all cell records corresponding to the same donor identifier within the associated expression dataset D: stage two is written into cells 1, 2, and 3 corresponding to donor A; stage three is written into cells 4, 5, and 6 corresponding to donor B. After writing, the stage-associated expression dataset E is obtained. At this point, any cell record in the stage-associated expression dataset E simultaneously possesses a cell identifier, gene identifier, expression value, donor identifier, and pathological stage label. For example, cell 2 includes the donor identifier "donor A" and the pathological stage label "stage two," and retains its corresponding gene expression value. Through this stage-associated expression dataset E, subsequent processing can summarize cell records according to donor identifiers and perform comparative calculations according to pathological stage labels, thereby achieving the purpose of direct summarization and calculation by donor.

[0083] Methods for selecting donor subsets with complete phase annotations include:

[0084] Using the stage-related expression dataset as input, a donor summary table is obtained by summarizing the data according to the donor identifier. The donor summary table includes at least the donor identifier, the corresponding pathological stage label, and the number of cell records.

[0085] The completeness of staging annotations is determined in the donor summary table. When the pathological staging label is empty or has a missing code, it is marked as incomplete staging. When the pathological staging label does not belong to the preset staging label set, it is marked as incomplete staging. When the same donor identifier corresponds to multiple different pathological staging labels, it is determined as a staging conflict and marked as incomplete staging according to the mutual exclusion rule. The mutual exclusion rule includes at least one pathological staging label allowed for the same donor.

[0086] Remove donor identifiers marked as incomplete phases from the donor summary table, retain the remaining donor identifiers to form a donor set with complete phase annotations, and define the donor set with complete phase annotations as a donor subset.

[0087] The staging-related expression dataset is screened based on the donor subset, retaining only cell records whose donor identifiers belong to the donor subset, thus obtaining the sub-expression dataset corresponding to the donor subset. The purpose is to limit all subsequent calculations involving pathological staging to the donor with complete staging annotation, avoid interference from missing or conflicting staging on donor layer statistics, and provide a consistent input set for subsequent fitting and comparison with the donor as the smallest statistical unit.

[0088] Step S102: Extract the target cell set based on the donor subset and construct a continuous state representation of each cell (e.g., Figure 2 and Figure 5As shown, continuous state representation is used to characterize the relative position of the target cell in a continuous dynamic process, and serves as a continuous independent variable for subsequent gene-by-gene dynamic modeling and trajectory-related feature extraction.

[0089] Obtain the sub-expression dataset corresponding to the donor subset. The sub-expression dataset includes at least cell identifiers, donor identifiers, pathological stage labels, and gene expression values. Obtain cell type annotation fields from the sub-expression dataset and compare them with a preset target cell type label set. Cell records whose cell type annotations belong to the preset target cell type label set are retained to form a target cell set. The target cell set is a set of target cell types limited to the donor subset to ensure that the computational objects for subsequent continuous state representation are consistent with the donor subset.

[0090] Expression value standardization is performed on the target cell set to form an input matrix that can be used for continuous state representation computation: For each cell record, the gene expression value of the cell record is converted into a normalized expression value according to a preset normalization rule. The preset normalization rule includes at least scaling the expression value proportionally to the total expression value of the cell record and performing a logarithmic transformation. Based on the normalized expression value, a set of feature genes is selected according to a preset feature selection rule. The preset feature selection rule includes at least calculating the variance or dispersion of each gene in the target cell set and comparing it with a preset feature threshold. Genes that meet the threshold condition are used as the set of feature genes. Based on this, a feature expression matrix of cell identifier-feature gene expression vector is constructed to provide a unified input for the next step of adjacency relationship construction. Normalization can be performed using the SCTransform method.

[0091] Cell adjacency relationships are constructed based on the feature expression matrix, and cell records that do not meet the consistency requirements of the target cells are removed by comparison with a threshold. For any cell, the distance between it and other cells on the feature expression matrix is ​​calculated, and a set of nearest neighbors is selected according to a preset number of nearest neighbors to form an adjacency relationship. At the same time, the non-target percentage of each cell is calculated based on a preset set of non-target marker genes. The non-target percentage is the ratio of the sum of the expression values ​​of the cell on the preset set of non-target marker genes to the sum of the expression values ​​of the cell's whole genes. The non-target percentage is compared with a preset percentage threshold. Cell records with a non-target percentage not less than the preset percentage threshold are removed from the target cell set to obtain the processed target cell set. The corresponding adjacency relationship is regenerated based on the processed target cell set to ensure that the continuous state representation calculation is based on a consistent cell set and adjacency relationship.

[0092] Construct a continuous state representation for each cell based on the processed target cell set and its adjacency relationships:

[0093] A random walk transition matrix is ​​constructed based on adjacency relations. The elements of the random walk transition matrix are calculated by adjacency relations and distance weights according to a preset weighting rule.

[0094] A continuous state starting cell set is determined through threshold comparison. This includes calculating a starting feature score for each cell, which is the sum of expression values ​​of that cell on a preset starting feature gene set. The starting feature score is compared with a preset starting threshold, and cells with a starting feature score not less than the preset starting threshold are selected as starting candidates. Next, starting cells are selected from the starting candidates according to preset selection rules, which at least include selecting the cell with the smallest average transfer distance from other cells as the starting cell. Then, starting from the starting cell, the diffusion distance of each cell relative to the starting cell is calculated based on the random walk transition matrix, and the cells are sorted in ascending order of diffusion distance. The sorting position or diffusion distance value is used as the continuous state representation of that cell, thereby generating a comparable one-dimensional continuous state representation for each cell.

[0095] The continuous state representation is written back into the cell record corresponding to the processed target cell set, and the same recording granularity is maintained with the cell identifier, donor identifier and pathological stage label to form a continuous state associated expression dataset. The continuous state associated expression dataset serves as the input for subsequent modeling of gene expression dynamics and forming a candidate transcription program set, so that subsequent steps can directly characterize tissue expression changes according to the continuous state and perform summary calculations according to the donor identifier.

[0096] Through the above steps, the purpose of step S102 is to obtain the target cell set within the donor subset and generate a continuous state representation for each cell in the target cell set, so that cell state changes are expressed in the form of sortable continuous variables. This provides a unified input for subsequent gene expression dynamic modeling based on continuous variables and donor layer statistics. For ease of understanding, a specific example is provided below:

[0097] The target cell set is extracted based on the donor subset; a stage-associated expression dataset E is used, where donor A corresponds to stage two and donor B corresponds to stage three, and the donor subset remains donor A and donor B. The sub-expression dataset F corresponding to the donor subset includes cells 1 to 6, and includes cell identifiers, donor identifiers, pathological stage labels, and gene expression values; the sub-expression dataset F also includes cell type annotation fields, for example, cell 1 is annotated as target cell type one, cell 2 is annotated as target cell type one, cell 3 is annotated as non-target cell type, cell 4 is annotated as target cell type one, cell 5 is annotated as target cell type one, and cell 6 is annotated as non-target cell type; the preset target cell type label set is set as target cell type one; after comparison, cells 1, 2, 4, and 5 are retained to form the target cell set G; where cells 1 and 2 belong to donor A, and cells 4 and 5 belong to donor B, thus the target cell set G is still limited to the donor subset.

[0098] Expression value uniformity and feature expression matrix construction; normalization is performed on each cell record in the target cell set G. Assume the target cell set G contains only three gene identifiers: gene A, gene B, and gene C; the original gene expression values ​​are as follows:

[0099] In cell 1, gene A has 12, gene B has 3, and gene C has 5, with a total expression level of 20; in cell 2, gene A has 8, gene B has 4, and gene C has 8, with a total expression level of 20; in cell 4, gene A has 2, gene B has 10, and gene C has 8, with a total expression level of 20; in cell 5, gene A has 1, gene B has 9, and gene C has 10, with a total expression level of 20. When scaling by the total expression level, since the total expression levels of the four cells are the same, the scaling does not change the relative proportions. Then, a logarithmic transformation is performed on the expression values ​​of each gene to obtain normalized expression values; subsequently, calculations are performed... Calculate the variance or dispersion of each gene in the target cell set G. For example, the values ​​of gene A and gene B differ significantly across the four cells, while the values ​​of gene C differ less. After comparing the variance or dispersion with a preset feature threshold, assume that gene A and gene B are retained as the feature gene set. Construct a feature expression matrix H, where each row corresponds to a cell identifier and each column corresponds to a feature gene. The feature gene expression vector of cell 1 is formed by the normalized values ​​of gene A and gene B. The same applies to cells 2, 4, and 5.

[0100] Adjacency relationship construction and threshold-based elimination; calculation of inter-cell distance based on feature expression matrix H; taking cell 1 as an example, calculating the distance between cell 1 and cell 2, cell 4, and cell 5 on the feature gene expression vector; with a preset nearest neighbor number of 2, the nearest neighbor set of cell 1 selects the two cells with the smallest distance, such as cell 2 and cell 4, forming the adjacency relationship between cell 1 and cell 2, and cell 1 and cell 4; the remaining cells form the adjacency relationship set in the same way.

[0101] Simultaneously, the proportion of non-target genes is calculated based on a pre-defined set of non-target marker genes. Using the same dataset, the pre-defined set of non-target marker genes is assumed to include only gene D. Assuming that gene D expression in cell 1 is 1, and the sum of all gene expression values ​​in cell 1 is 21, then the proportion of non-target genes is 1 divided by 21. In cell 2, gene D is 0, so the proportion of non-target genes is 0. In cell 4, gene D is 6, and the sum of all gene expression values ​​in cell 4 is 26, so the proportion of non-target genes is 6 divided by 26. In cell 5, gene D is 1, and the sum of all gene expression values ​​in cell 5 is 2... 1. The non-target percentage is 1 divided by 21. The non-target percentage is compared with the preset percentage threshold. Assuming the preset percentage threshold is set to 0.15, the non-target percentage of cell 4 is about 0.23, which is not less than the preset percentage threshold. Cell 4 is removed from the target cell set G, and the processed target cell set I is cell 1, cell 2, and cell 5. The distance and nearest neighbor set are recalculated based on the processed target cell set I, and the adjacency relationship set is regenerated to ensure that the subsequent continuous state representation is based on the consistent cell set and adjacency relationship.

[0102] Continuous state representation construction; random walk transition matrix is ​​constructed based on the processed target cell set I and its adjacency relationship; taking cell 1 as an example, its nearest neighbor set is assumed to be cell 2 and cell 5, then the transition weights are assigned to cell 1 to cell 2 and cell 1 to cell 5 according to the distance weight rule, and the sum of the transition weights of the row where cell 1 is located is 1, forming the corresponding row of the random walk transition matrix.

[0103] Determine the set of starting cells in a continuous state. Assume a preset set of starting characteristic genes includes gene A, and the starting characteristic score is the normalized expression value of that cell on gene A. Calculate the starting characteristic scores for cells 1, 2, and 5 respectively, and compare them with a preset starting threshold. If cell 5's starting characteristic score is lower than the preset starting threshold, while cells 1 and 2 are not lower, then cells 1 and 2 are the starting candidates. Next, select the starting cell according to a preset selection rule: the cell with the smallest average transfer distance to other cells is selected as the starting cell. For example, if cell 1's average transfer distance is less than cell 2's, then cell 1 is the starting cell.

[0104] Starting with cell 1, calculate the diffusion distances of cells 2 and 5 relative to cell 1 based on the random walk transition matrix, and sort them in ascending order of diffusion distance. Assuming the diffusion distance results are 0 for cell 1, 0.3 for cell 2, and 0.8 for cell 5, the sorting order is cell 1, cell 2, and cell 5. Use the sorting position or diffusion distance value as a continuous state representation. For example, if the continuous state representation is the diffusion distance value, then cell 1 is 0, cell 2 is 0.3, and cell 5 is 0.8.

[0105] A continuous state-associated expression dataset is formed. The continuous state representations are written back to the cell records corresponding to the processed target cell set I, maintaining the same record granularity as the cell identifier, donor identifier, and pathological stage label, forming a continuous state-associated expression dataset J. For example, the cell 1 record includes the donor identifier "Donor A", the pathological stage label "Stage II", and the continuous state representation "0". The cell 2 record includes the donor identifier "Donor A", the pathological stage label "Stage II", and the continuous state representation "0.3". The cell 5 record includes the donor identifier "Donor B", the pathological stage label "Stage III", and the continuous state representation "0.8". Subsequent steps can directly organize the expression changes based on the continuous state representations and summarize and calculate the cell records according to the donor identifier.

[0106] Step S103: Model gene expression dynamics based on continuous state representation, and cluster genes with similar dynamic characteristics to obtain a set of candidate transcription programs (e.g., ...). Figure 5 As shown, the process of constructing candidate transcription programs is used to organize the dynamic changes of genes on a continuous trajectory into a set of programs that can be used for subsequent screening.

[0107] Methods for modeling gene expression dynamics based on continuous state representation:

[0108] Optionally, the construction of candidate transcription programs can be completed collaboratively based on pseudo-time inference, trend modeling, dynamic gene screening, trend clustering, and other steps. The key steps and parameter settings are shown in Table 1.

[0109] Table 1. Examples of steps and parameters for constructing candidate transcription programs based on trajectory information.

[0110]

[0111] Obtain a continuous state association expression dataset, acquire cell identifiers, continuous state representations, and gene expression values, and establish gene-by-gene modeling inputs with continuous state representations as independent variables and gene expression values ​​as dependent variables.

[0112] Availability screening is performed on the gene set. The number of cells with non-zero expression values ​​for each gene in the continuous state associated expression dataset is counted. The number of cells is compared with the preset detection cell number threshold. Genes that meet the detection cell number threshold are retained. After sorting by cell number, the top preset number are selected to form a candidate modeling gene set, so as to reduce the impact of low expression noise and computational burden.

[0113] Dynamic discrimination is performed on the candidate modeling gene set. For each gene, a smooth fitting model is constructed with continuous state representation as input. The smooth fitting model outputs the fitted expression curve of the gene on the continuous state representation. Dynamic statistics are calculated based on the fitting model. Dynamic statistics are used to measure the degree of deviation of the fitting curve from the constant level. The first preset number of genes are selected according to the dynamic statistics to form a dynamic gene set.

[0114] Dynamic features for clustering are generated. A uniform continuous state grid is generated within the value range of the continuous state representation. The dynamic gene set is input one by one into the smooth fitting model, and the predicted value sequence of the gene is output as a trend vector on the uniform continuous state grid. Each trend vector is standardized by gene to ensure that the gene trend vector meets the preset mean and variance constraints in the grid dimension, thereby highlighting the differences in trend shape and weakening the differences in absolute magnitude. A specific example is as follows:

[0115] Obtain a continuous-state association expression dataset and establish gene-by-gene modeling input. The continuous-state association expression dataset J contains cells 1, 2, and 5. Cells 1 and 2 belong to donor A and are labeled as stage 2, while cell 5 belongs to donor B and is labeled as stage 3. The continuous-state representations are 0 for cell 1, 0.3 for cell 2, and 0.8 for cell 5. Gene expression values ​​are taken from genes A, B, and C, with gene D added as a subsequent screening example. Establish gene-by-gene modeling... When inputting the model, the continuous state representation is used as the independent variable sequence 0, 0.3, 0.8, and the expression values ​​of each gene in cell 1, cell 2, and cell 5 are used as the dependent variable sequence. For example, the dependent variable sequence of gene A is 12 in cell 1, 8 in cell 2, and 1 in cell 5; the dependent variable sequence of gene B is 3 in cell 1, 4 in cell 2, and 9 in cell 5; the dependent variable sequence of gene C is 5 in cell 1, 8 in cell 2, and 10 in cell 5; and the dependent variable sequence of gene D is 1 in cell 1, 0 in cell 2, and 1 in cell 5.

[0116] Availability screening is performed on the gene set to form a candidate modeling gene set. The number of cells with non-zero expression values ​​for each gene in the continuous state associated expression dataset J is counted. Gene A has non-zero values ​​in cells 1, 2, and 5, with a cell count of 3; Gene B has non-zero values ​​in cells 1, 2, and 5, with a cell count of 3; Gene C has non-zero values ​​in cells 1, 2, and 5, with a cell count of 3; Gene D has non-zero values ​​in cells 1 and 5, but zero values ​​in cell 2, with a cell count of 2. The cell counts are compared with a preset detection cell count threshold. Assuming the preset detection cell count threshold is set to 3, Gene A, Gene B, and Gene C meet the threshold condition and are retained, while Gene D does not meet the threshold condition and is not retained. When the retained genes are sorted by cell count, the three are the same. The top 3 are then selected according to the preset number to form the candidate modeling gene set of Gene A, Gene B, and Gene C.

[0117] Dynamic discrimination is performed on the candidate modeling gene set to form a dynamic gene set. For each gene in the candidate modeling gene set, a smooth fitting model (such as using a generalized additive model) is constructed with continuous state representation as input. Taking gene A as an example, the inputs are 0, 0.3, and 0.8, corresponding to expression values ​​of 12, 8, and 1. The smooth fitting model outputs the fitted expression curve of gene A on the continuous state representation. The fitted expression curves of genes B and C are obtained in the same way. Dynamic statistics are calculated based on the fitting model. Dynamic statistics are used to measure the degree of deviation of the fitted expression curve from the constant level. For example, the constant level is taken as the average expression value of the gene in three cells. The average value of gene A is 7, with a large deviation at the three points; the average value of gene B is 5.33, with a moderate deviation; and the average value of gene C is 7.67, with a moderate deviation. After sorting the dynamic statistics of the three genes, assuming the order is gene A first, gene C second, and gene B third; and then selecting the first two according to a preset number, the dynamic gene set is gene A and gene C.

[0118] A trend vector for clustering is generated and gene-based normalization is performed. A uniform continuous state grid is generated within the value range of 0 to 0.8 for continuous state representation. Assume the uniform continuous state grid has 5 points: 0, 0.2, 0.4, 0.6, and 0.8. Gene A from the dynamic gene set is input into its smoothing model, and the predicted value sequence is output at the above 5 grid points to obtain the trend vector of gene A. For example, it is predicted as 12 at 0, 10 at 0.2, 7 at 0.4, 4 at 0.6, and 1 at 0.8. Gene C is input into its smoothing model, and the predicted value sequence is output at the same grid points to obtain the trend vector of gene C. For example, it is predicted as 5 at 0, 6 at 0.2, 8 at 0.4, 9 at 0.6, and 10 at 0.8.

[0119] For each trend vector, standardization is performed on the gene to ensure that the trend vector satisfies preset mean and variance constraints in the grid dimension. For example, the preset mean constraint is that the standardized mean is 0, and the preset variance constraint is that the standardized variance is 1. Then, for gene A, the five predicted values ​​are subtracted from its mean and divided by its standard deviation to obtain the standardized trend vector of gene A. The same applies to gene C to obtain the standardized trend vector of gene C. After completion, a dynamic feature set for clustering is obtained, which can directly compare the differences in trend shape on the same uniform continuous state grid and reduce the impact of absolute amplitude differences on subsequent clustering.

[0120] Methods for clustering genes with similar dynamic characteristics to obtain a set of candidate transcription programs include:

[0121] Using trend vectors as the dynamic feature vectors of genes, the dynamic distance between any two genes is defined. The dynamic distance is obtained by the point-by-point difference between the two trend vectors on a uniform continuous state grid according to a preset aggregation rule. For each gene, a preset number of nearest neighbors is selected in ascending order of dynamic distance to obtain the set of nearest neighbor genes for that gene.

[0122] Gene adjacency graphs are constructed based on the set of nearest neighbor genes for each gene. Nodes in the graph are gene identifiers, and edge connections are determined by the nearest neighbor relationship. Edge weights are obtained by dynamic distance according to a preset weight conversion rule, so that the clustering input is completely derived from the trend vector.

[0123] The gene adjacency graph is divided and clustered, and multiple gene clusters are output under the constraints of preset resolution parameters. Each gene cluster is defined as a candidate transcription program, and all candidate transcription programs are summarized to obtain a candidate transcription program set.

[0124] To facilitate subsequent program quantification and visualization, the trend amplitude is calculated for genes within each candidate transcription program. The trend amplitude is the difference between the maximum and minimum predicted values ​​of the gene's trend vector. A representative gene list for each candidate transcription program is formed by sorting genes by trend amplitude and selecting the top preset number of genes, and this list is stored in association with the candidate transcription program identifier. A specific example is shown below:

[0125] Calculate the dynamic distance and obtain the nearest neighbor gene set, constructing a gene adjacency graph with a uniform continuous grid of 0, 0.2, 0.4, 0.6, and 0.8. Using genes A and C from the dynamic gene set above, add gene E, whose trend vector is obtained through smoothing and standardization, to demonstrate nearest neighbor selection and graph construction. Assume that gene E's trend vector on the above grid is predicted as 11 at 0, 9 at 0.2, 6 at 0.4, 3 at 0.6, and 1 at 0.8, and obtain the standardized trend vector after gene standardization.

[0126] When defining the dynamic distance between any two genes, the point-by-point difference between the two trend vectors at each grid point is used as the basis, and the results are summarized according to a preset summation rule. For example, the preset summation rule takes the square root of the sum of the squares of the point-by-point differences. Therefore, when calculating the dynamic distance between gene A and gene C, the difference between their predicted values ​​is calculated for each of the five grid points, the squares are summed, and the square root is taken to obtain a distance value. The calculation of the dynamic distance between gene A and gene E is similar. The calculation of the dynamic distance between gene C and gene E is similar. If the calculation result shows that the dynamic distance between gene A and gene E is less than the dynamic distance between gene A and gene C, and the dynamic distance between gene C and gene E is greater than the dynamic distance between gene A and gene C, then when sorting gene A by dynamic distance from smallest to largest, gene E is the nearest neighbor with the smallest distance, and gene C is the second nearest neighbor with the second smallest distance.

[0127] With a preset nearest neighbor count of 1, the nearest neighbor set of gene A is assumed to be gene E, the nearest neighbor set of gene C is assumed to be gene A, and the nearest neighbor set of gene E is assumed to be gene A. When constructing the gene adjacency graph based on the nearest neighbor sets of each gene, the nodes in the graph are gene A, gene C, and gene E. Edge connections are determined by the nearest neighbor relationships; therefore, there are edges connecting gene A and gene E, genes C and gene A, and gene E and gene A. Edge weights are obtained by calculating the dynamic distance according to a preset weight conversion rule. For example, the preset weight conversion rule takes the edge weight as 1 divided by 1 plus the dynamic distance. Thus, the smaller the dynamic distance, the larger the edge weight, thereby ensuring that the clustering input is completely derived from the trend vector.

[0128] Graph partitioning and clustering are performed on the gene adjacency graph to obtain a set of candidate transcription programs. The clustering process outputs multiple gene clusters under the constraint of a preset resolution parameter. For example, when the preset resolution parameter is set to a lower value, gene A and gene E are assigned to the same gene cluster due to their larger edge weights, while gene C is assigned to a different gene cluster due to its smaller edge weight to gene A. Therefore, two gene clusters are output: the first gene cluster consists of gene A and gene E, and the second gene cluster consists of gene C. Each gene cluster is defined as a candidate transcription program. The set of candidate transcription programs includes candidate transcription program one (corresponding to gene A and gene E) and candidate transcription program two (corresponding to gene C).

[0129] The trend amplitude is calculated, and a representative gene list is generated and stored in association. To facilitate subsequent program quantification and display, the trend amplitude is calculated for each gene in the candidate transcription program. The trend amplitude is the difference between the maximum and minimum predicted values ​​of the trend vector for that gene. Using the unstandardized predicted value sequence for illustration, the maximum predicted value of the trend vector for gene A is 12, the minimum predicted value is 1, and the trend amplitude is 11; the maximum predicted value of the trend vector for gene E is 11, the minimum predicted value is 1, and the trend amplitude is 10; the maximum predicted value of the trend vector for gene C is 10, the minimum predicted value is 5, and the trend amplitude is 5.

[0130] Genes A and E within candidate transcription program one are sorted by trend amplitude, with gene A being greater than gene E. The preset quantity is 1, so the representative gene list of candidate transcription program one is gene A. Candidate transcription program two only includes gene C, so its representative gene list is gene C. The representative gene list is associated and stored with the identifiers of candidate transcription program one and candidate transcription program two, respectively, to obtain the displayable result of the candidate transcription program set, which is consistent with the clustering result derived from the trend vector above.

[0131] Step S104: Perform donor resampling on the candidate transcription program set (e.g., Figure 3 As shown, staged and stratified donor resampling performs coefficient estimation and stability quantification of candidate transcription programs at the donor layer level, and generates a SASI score based on stability screening. Based on the donor set obtained from resampling: the cell expression of candidate transcription programs in the same donor is summarized according to a preset aggregation rule to obtain the donor layer program quantification value; a fitting relationship is established between pathological stage information and the donor layer program quantification value to obtain the evaluation result; the effect stability index, directional consistency index, and donor coverage index across resampling are calculated, and the SASI score is obtained by combining the above indices (further as follows). Figure 6 As shown, the joint screening based on effect stability index, orientation consistency index and donor coverage index is used to suppress sporadic candidates driven by a few donors or caused by orientation instability, thereby improving the stability of the screening results under different donor extraction conditions.

[0132] Optionally, to facilitate the description of the donor layer fitting and stability quantification process, the key quantities can be defined using the following mathematical expressions:

[0133] (1) Donor layer fitting relationship: ,in For candidate transcription programs The donor layer procedure quantization value, Let f(·) be the pathological stage variable, and f(·) be the smoothing spline function. This is the residual term.

[0134] (2) Effect size of the b-th resampling: ,in and These are the highest and lowest installment periods, respectively. (·) represents the predicted value of the fitted curve at the corresponding stage.

[0135] (3) Indicators of effect stability: ,in (·) represents the standard deviation. It is a very small positive number.

[0136] (4) Directional consistency index:

[0137] ,in, In the first On the resampled data, for the object The estimated effect parameters or regression coefficients are represented by their positive or negative signs, indicating the direction of influence. The 1(·) is an indicator function, taking the value 1 if the condition within the parentheses is true, and 0 otherwise. This represents the number of resampling iterations.

[0138] (5) Donor coverage index: ,in, Used for statistics In the second resampling iteration, the candidate transcription program The number of iterations corresponding to a positive coefficient estimation result.

[0139] (6) SASI score: .

[0140] (7) Donor layer scoring indicators: ,in As a core set of transcription programs, For candidate transcription programs The donor layer procedure quantification value or its derived score, The weights are determined and solidified by the global effect.

[0141] Methods for resampling donors from a set of candidate transcription programs, and methods based on the resampled donor set, include:

[0142] Obtain donor subsets and their pathological staging information, and divide the donor subsets into multiple staging layers according to the pathological staging information. The staging layers are used to limit the source of donor extraction.

[0143] Set the number of resampling iterations and the extraction rules for each iteration. In each iteration, donors are extracted and merged according to the extraction rules within each phase layer to obtain the donor set obtained by resampling in that iteration. This donor set is then stored in association with the iteration identifier. The purpose is to form a phase-composition controlled donor set sequence over multiple iterations, so that subsequent evaluations can be repeatedly calculated under different donor composition perturbations. A specific example is as follows:

[0144] Continuing with the donor subsets from the previous text, donor A and donor B are identified, and the pathological staging information is staging two for donor A and staging three for donor B. When the donor subset is divided into multiple staging layers based on the pathological staging information, staging layer two includes donor A, and staging layer three includes donor B. The staging layers two and three are used to limit the source of subsequent donor extraction, that is, the extraction is carried out within staging layers two and three, respectively.

[0145] The resampling iteration is set to 3 times, and the sampling rule for each iteration is to extract one donor from each phase layer and merge them. In the first iteration, donor A is extracted from phase two layer and donor B is extracted from phase three layer. The two are merged to obtain the donor set obtained by resampling in this iteration, which is donor A and donor B. This donor set is associated with the iteration identifier Iteration 1 and stored. In the second iteration, the same sampling rule is used to extract donor A from phase two layer and donor B from phase three layer. The resampling donor set is obtained as donor A and donor B and associated with the iteration identifier Iteration 2 and stored. The third iteration is similar, obtaining the resampling donor set as donor A and donor B and associated with the iteration identifier Iteration 3 and stored. Through the above settings, a phased controlled donor set sequence is formed in multiple iterations. Subsequent evaluations can be performed under different iteration identifiers corresponding to iteration 1, iteration 2, and iteration 3, thereby repeating the subsequent evaluation process on the donor sets of different iterations.

[0146] Methods for summarizing the cellular expression of candidate transcription programs in the same donor according to a preset aggregation rule to obtain the donor layer program quantification value include:

[0147] For any candidate transcription program, obtain the gene set corresponding to that candidate transcription program in the candidate transcription program set, and obtain the gene expression values, cell identifiers, and donor identifiers corresponding to the gene set in the continuous state associated expression dataset.

[0148] For any donor identifier in the donor set obtained by resampling, all cell records corresponding to that donor identifier are extracted, and the gene set of the candidate transcription program is summarized by a preset aggregation rule within that donor to obtain the donor layer program quantization value on the candidate transcription program; wherein, the donor layer program quantization value is the program value after summarizing the cell expression within the donor with the donor as the statistical unit, in order to replace the scope of direct modeling at the cell granularity.

[0149] The donor layer program quantization value is associated with the donor identifier, candidate transcription program identifier, and iteration identifier to form an in-iteration donor layer program quantization table. The purpose is to ensure that subsequent fitting is performed with the donor as the smallest statistical unit, avoiding inference bias caused by treating multiple cell records from the same donor as independent samples. A specific example is as follows:

[0150] Using the continuous state association expression dataset J from the previous section, cells 1 and 2 belong to donor A and are labeled as stage 2, while cell 5 belongs to donor B and is labeled as stage 3. The donor set obtained from the resampling in the first iteration is also used, designated as donor A and donor B, with the iteration identifier being iteration 1. The candidate transcription program set follows the clustering example from the previous section, where candidate transcription program 1 corresponds to gene set A and gene E, and candidate transcription program 2 corresponds to gene set C. Candidate transcription program 1 will now be used as any candidate transcription program for illustration.

[0151] Obtain the gene set and expression values, cell identifier and donor identifier; obtain the gene set corresponding to candidate transcription program one in the candidate transcription program set, and obtain gene A and gene E; then obtain the gene expression values ​​corresponding to gene A and gene E from the continuous state associated expression dataset J, and at the same time obtain the cell identifier and donor identifier in the corresponding record; for example, the donor identifier of cell 1 is donor A, the gene A expression value of cell 1 is 12, and the gene E expression value of cell 1 is 11; the donor identifier of cell 2 is donor A, the gene A expression value of cell 2 is 8, and the gene E expression value of cell 2 is 9; the donor identifier of cell 5 is donor B, the gene A expression value of cell 5 is 1, and the gene E expression value of cell 5 is 1.

[0152] Within the same donor, the donor layer program quantization value is obtained by summarizing according to the preset aggregation rules. For donor A in the resampled donor set, all cells corresponding to donor A are extracted and recorded as cell 1 and cell 2. Within donor A, the gene sets gene A and gene E are summed according to the preset aggregation rules. For example, the preset aggregation rule is to first calculate the average expression value of the same gene in all cells within the donor, and then calculate the average gene value within the program. Thus, the cell mean of gene A in donor A is the average of 12 in cell 1 and 8 in cell 2, which is 10. The cell mean of gene E in donor A is the average of 11 in cell 1 and 9 in cell 2, which is 10. The average of the two gene values ​​is then calculated within the program to obtain the donor layer program quantization value of donor A in candidate transcription program one, which is 10.

[0153] For donor B in the resampled donor set, all cells corresponding to donor B are extracted and recorded as cell 5. Within donor B, the cells are aggregated according to the same preset aggregation rules. Since donor B only includes cell 5, the cell mean value of gene A in donor B is 1, and the cell mean value of gene E in donor B is 1. The program calculates the average of the two gene mean values ​​to obtain the donor layer program quantization value of donor B on candidate transcription program one as 1. Thus, the donor layer program quantization value is obtained with the donor as the statistical unit, providing the program value caliber of donor granularity for subsequent fitting.

[0154] The donor layer program quantization value 10 of donor A is associated with the donor identifier donor A, the candidate transcription program identifier candidate transcription program one, and the iteration identifier iteration one to form one record; the donor layer program quantization value 1 of donor B is associated with the donor identifier donor B, the candidate transcription program identifier candidate transcription program one, and the iteration identifier iteration one to form another record; the two records together constitute the iterative donor layer program quantization table of iteration one; this iterative donor layer program quantization table is directly used in subsequent fitting with the donor as the smallest statistical unit, thereby avoiding inference bias caused by treating cell 1 and cell 2 in donor A as independent samples participating in fitting.

[0155] The evaluation results are obtained by establishing a fitting relationship between pathological staging information and donor layer procedural quantification values; the method for calculating the SASI score by combining the effect stability index, directional consistency index, and donor coverage index across resampling includes:

[0156] Within the donor set obtained by resampling in each iteration, the pathological stage information and donor layer program quantization value of that iteration are obtained. A fitting relationship is established with the pathological stage information as the independent variable and the donor layer program quantization value as the dependent variable. The fitting relationship adopts the form of smooth splines to characterize nonlinear changes, and the evaluation result of the candidate transcription program is output. The evaluation result includes at least the predicted value of the fitting curve at the stage extreme value.

[0157] The effect size of the iteration is calculated based on the evaluation results. The effect size is defined as the difference between the predictions of the fitted curve at the highest and lowest stages, and the sign of the difference is recorded as the direction of the effect.

[0158] The effect size series is aggregated across all iterations, and the effect stability index is calculated according to the preset aggregation rules. The effect stability index consists of the standard deviation of the effect size series and the smallest positive number, which is used to reflect the dispersion of the effect size under multiple resampling.

[0159] The directional consistency index is calculated by counting the number of times the statistical effect occurs across all iterations. The directional consistency index is the proportion of the larger of the positive and negative counts to the total number of iterations.

[0160] The donor coverage index is obtained by calculating the percentage of iterations across all iterations that yield statistical effect sizes. This index reflects whether the candidate transcription program can consistently output effect sizes under different resampled donor sets.

[0161] Based on the effect size sequence obtained through resampling, the effect stability index, orientation consistency index, and donor coverage index of candidate transcription programs are calculated. These three indices are then combined to calculate the SASI score of the candidate transcription program, which is then associated with and stored in relation to the corresponding candidate transcription program identifier. The aim is to extend single-step fitting to effect sequences obtained through resampling iterations and to impose constraints on candidate transcription programs using three types of computable indices. This ensures that subsequent outputs preferentially come from candidate transcription programs that can consistently produce consistent orientation and computable effects under multiple donor subset perturbations. A specific example is as follows:

[0162] A fitting relationship is established and the evaluation results are output. The donor set obtained from the resampling in Iteration 1 is donor A and donor B. The pathological staging information is stage 2 for donor A and stage 3 for donor B. The donor layer program quantization value follows the result of candidate transcription program 1 above, with donor A set to 10 and donor B set to 1. After obtaining the pathological staging information and donor layer program quantization value for this iteration, a fitting relationship is established using the pathological staging information as the independent variable and the donor layer program quantization value as the dependent variable. The fitting relationship uses a smooth spline form, with input points set to 10 for stage 2 and 1 for stage 3. The evaluation results of candidate transcription program 1 in Iteration 1 are output. The evaluation results include at least the predicted values ​​of the fitted curve at the stage extremes. Here, the lowest stage is stage 2, and the highest stage is stage 3. The predicted value of the fitted curve at stage 2 is 10, and the predicted value at stage 3 is 1.

[0163] Calculate the effect size and record the effect direction. The effect size is calculated based on the evaluation results of iteration one. The effect size is defined as the difference between the predictions of the fitted curve at the highest and lowest stages, that is, the predicted value of stage three is 1 minus the predicted value of stage two is 10, and the effect size is negative. Record the sign of the difference as the effect direction, which is negative in this case.

[0164] The effect size sequences are summarized and the effect stability index is calculated. To meet the requirements of cross-iteration summarization, the resampling settings above are used for iteration 1, iteration 2, and iteration 3. It is assumed that the donor layer program quantization value changes due to differences in the cell composition in the donor body in iteration 2 and iteration 3. For candidate transcription program 1, the following evaluation results are obtained: In iteration 2, the predicted value at stage 2 is 9, and the predicted value at stage 3 is 2, so the effect size is 2 minus 9, which is negative; In iteration 3, the predicted value at stage 2 is 11, and the predicted value at stage 3 is 1, so the effect size is 1 minus 11, which is negative. Thus, the effect size sequences across all iterations are summarized as the negative difference of iteration 1, the negative difference of iteration 2, and the negative difference of iteration 3. The effect stability index is calculated according to the preset summarization rules. The preset summarization rules are to calculate the standard deviation of the effect size sequence and add it to the smallest positive number to obtain the effect stability index, which is used to reflect the dispersion of the effect size under three resampling.

[0165] Calculate the directional consistency index by counting the number of times the effect direction appears in the three iterations. If the effect direction is negative in iteration 1, iteration 2, and iteration 3, then the number of negative effects is 3, and the number of positive effects is 0. The directional consistency index is the ratio of the larger of the number of positive effects and the number of negative effects to the total number of iterations. Therefore, the directional consistency index is 3 divided by 3, which equals 1.

[0166] To calculate the donor coverage index, the percentage of iterations in which the statistical effect size can be calculated is assumed to be 3. Assuming that iterations 1, 2, and 3 can all establish a fitting relationship and obtain the predicted values ​​at the extreme values ​​of each period, and thus the effect size can be calculated, the number of iterations that can be calculated is 3, which accounts for 1 of the 3 iterations, and the donor coverage index is 1.

[0167] Based on the effect size sequence obtained through cross-resampling, the effect stability index, orientation consistency index, and donor coverage index of candidate transcription programs are calculated, and the SASI score of the candidate transcription program is calculated by combining the three indices. For example, in a specific implementation, the effect size sequence of a candidate transcription program is obtained after multiple staged and stratified resampling. The effect stability index, orientation consistency index, and donor coverage index are calculated based on the effect size sequence, and the SASI score of the candidate transcription program is calculated by combining the three indices. After calculating the corresponding SASI scores for all candidate transcription programs, they are sorted from high to low according to the SASI scores, and the candidate transcription programs with the highest ranking are selected as the core transcription program set. In this way, the single-fit result is extended to the cross-iteration robustness evaluation result, so that the subsequent output comes first from the candidate transcription programs that still show stable effects and orientation consistency under multiple iterations.

[0168] Step S105: Based on the ranking results of the SASI scores, a core transcription program set is selected, and a donor layer scoring index is constructed from the core transcription program set. When external validation is required, a validation single-cell transcription expression dataset from a different source than the single-cell transcription expression dataset is obtained. While keeping the core transcription program set and its weights unchanged, the donor layer scoring index is calculated according to the donor layer caliber to evaluate the correlation and discriminative ability of the donor layer scoring index with pathological staging information (e.g., Figure 6 (As shown).

[0169] Methods for selecting a core set of transcription programs based on SASI score ranking results and constructing donor layer scoring indicators from the core set of transcription programs include:

[0170] Obtain the candidate transcription program set and the corresponding SASI score for each candidate transcription program. Sort the programs by SASI score from high to low and select the core transcription program set according to the preset output quantity.

[0171] Obtain the donor layer program quantization value, and summarize the donor layer program quantization value of each core transcription program in the core transcription program set according to the donor identifier to form the core program quantization vector corresponding to the donor identifier.

[0172] The global effect obtained based on the fitting relationship determines the weight value for each core transcription program in the core transcription program set, and the weight values ​​are associated with the core transcription program identifier to form a weight set and then stored in a fixed manner; the weight set is a fixed input for subsequent cross-source data reuse.

[0173] For any donor identifier, the donor's core program quantization vector and weight set are synthesized according to a preset weighted aggregation rule to obtain the donor-level score index. This score index is then stored in association with the donor identifier and pathological stage information. The donor-level score index is a continuous value at the donor granularity, used to characterize the overall activity level of the core transcription program set within the donor. The purpose is to convert multi-program results into a single continuous score at the donor level, enabling subsequent analysis to be conducted directly at the donor level and maintaining the verifiability of the score construction process. A specific example is as follows:

[0174] Candidate core transcription programs are ranked based on SASI scores. Following the previous approach, the candidate transcription program set includes candidate transcription program one and candidate transcription program two. Candidate transcription program one corresponds to gene A and gene E, while candidate transcription program two corresponds to gene C. The SASI score for candidate transcription program one is set to 0.85, and the SASI score for candidate transcription program two is set to 0.40. The programs are sorted from highest to lowest SASI score, and the core transcription program set is selected based on a preset output quantity. When the preset output quantity is 1, since the initial core transcription program set does not exceed the preset output quantity, the core transcription program set is directly obtained as candidate transcription program one.

[0175] The donor layer program quantization values ​​are summarized to obtain the core program quantization vector. The donor layer program quantization table in iteration one above is used, where the donor layer program quantization value of candidate transcription program one is 10 for donor A and 1 for donor B. After obtaining the donor layer program quantization values, the donor layer program quantization values ​​of each core transcription program in the core transcription program set are summarized according to the donor identifier. Since the core transcription program set only includes candidate transcription program one, the core program quantization vector corresponding to donor A is a vector containing only the value 10, and the core program quantization vector corresponding to donor B is a vector containing only the value 1.

[0176] The weights are determined by the global effect and stored in a fixed manner. Following the effect size example of the fitting relationship above, the effect size of candidate transcription program 1 is negative across iterations, and its absolute effect size is relatively large in three iterations. The global effect obtained based on the fitting relationship determines the weight value of candidate transcription program 1. For example, the weight value is the result of normalizing the absolute value of the global effect. Since the core transcription program set only includes candidate transcription program 1, its weight value can be set to 1. The weight value of 1 is associated with the core transcription program identifier candidate transcription program 1 to form a weight set and stored in a fixed manner as a fixed input for subsequent cross-source data reuse.

[0177] The donor layer scoring index is obtained through weighted summation. For any donor identifier, the donor's core program quantization vector and weight set are synthesized according to a preset weighted summation rule. For example, the preset weighted summation rule is to sum the core program quantization vectors according to their weight values. Thus, the donor layer scoring index for donor A is 10 multiplied by 1, and the donor layer scoring index for donor B is 1 multiplied by 1, resulting in 1. The donor layer scoring index 10 for donor A is associated and stored with the donor identifier donor A and the pathological stage information stage 2. The donor layer scoring index 1 for donor B is associated and stored with the donor identifier donor B and the pathological stage information stage 3. In this way, the multidimensional expression of the core transcription program set in the donor is converted into a single continuous score at the donor layer, and subsequent analysis can directly use this scoring index for calculation at the donor layer level.

[0178] When external validation is required, methods include obtaining validation single-cell transcriptional expression datasets from sources different from the single-cell transcriptional expression datasets, and calculating donor layer scoring indices according to donor layer caliber to assess the correlation and discriminative ability of donor layer scoring indices with pathological staging information, such as:

[0179] When an external validation requirement is triggered, a validation single-cell transcriptional expression dataset from a different source than the single-cell transcriptional expression dataset is obtained. Simultaneously, the donor identification information and pathological stage information corresponding to the validation single-cell transcriptional expression dataset are obtained, so that the validation single-cell transcriptional expression dataset meets the conditions for establishing the correspondence between cell identifiers and donor identifiers, and the donor identifiers can be mapped to pathological stage information.

[0180] The set of core transcription programs and their weights that have been identified and solidified in the discovery dataset are invoked. The core transcription programs are then expressed and aggregated according to the donor layer in the validation single-cell transcription expression dataset to form the core program quantization vector of the validation data.

[0181] While keeping the weight set unchanged, each donor calculates a donor layer score index.

[0182] The correlation between donor layer scoring indicators and pathological staging information is calculated under the donor layer caliber. The correlation includes at least the correlation statistics based on donor ranking consistency, and the correlation statistics are compared with the preset correlation threshold to output the correlation conclusion.

[0183] The document calculates the discriminative power of donor layer scoring indicators for late and non-late stages under donor layer caliber. Discriminative power includes at least the area under the curve (AUC) obtained from threshold scanning, and the AUC is compared with a preset discriminative threshold to output the discriminative power conclusion. Specifically, the document clarifies the caliber for comparing stage association and late stage discriminative power in independent external cohorts. The aim is to reuse the same core transcription program set and the same weight set under different donor source compositions, outputting association and discriminative power results according to the donor layer caliber, thereby forming a quantifiable cross-source assessment basis. A specific example is shown below.

[0184] We acquire validation single-cell transcriptional expression datasets from different sources and establish donor mappings. Following the previous approach, the donor subsets in the discovery dataset are donor A and donor B, the core transcription program set is candidate transcription program one, and its gene set is gene A and gene E. In the fixed-store weight set, the weight of candidate transcription program one is set to 1. Now, when external validation is triggered, we acquire validation single-cell transcriptional expression dataset K from a different source than the discovery dataset. Validation single-cell transcriptional expression dataset K includes cell records from donor C and donor D, and simultaneously acquires donor identification information and pathological stage information. The donor identification information indicates that cells 7 and 8 come from donor C, and cells 9 and 10 come from donor D. The pathological stage information indicates that donor C corresponds to stage two, and donor D corresponds to stage three. Therefore, validation single-cell transcriptional expression dataset K satisfies the conditions for establishing the correspondence between cell identifiers and donor identifiers, and donor C and donor D can be mapped to stages two and three, respectively.

[0185] According to the preset aggregation rules, the donor layer program quantization value and core program quantization vector of the validation data are obtained. The core transcription program set and its gene set are obtained from the validation single-cell transcription expression dataset K, i.e., candidate transcription program one and its gene set genes A and E are obtained. The core transcription program set is then summarized within the same donor cell according to the preset aggregation rules. Following the preset aggregation rules mentioned above, the average expression value of the same gene in all cells within the donor cell is calculated first, and then the average gene expression value is calculated within the program. Assuming that the gene A expression values ​​in cells 7 and 8 of donor C in the validation data are 9 and 7 respectively, and the gene E expression values ​​are 10 and 8 respectively, then the gene A cells within donor C... The mean value is 8, the mean value of gene E cells in donor C is 9, and the mean value within the program is 8.5, resulting in a donor layer program quantization value of 8.5 for donor C in candidate transcription program one. Assuming that the expression values ​​of gene A in cells 9 and 10 of donor D are 2 and 1 respectively, and the expression values ​​of gene E are 2 and 1 respectively, then the mean value of gene A cells in donor D is 1.5, the mean value of gene E cells in donor D is 1.5, and the mean value within the program is 1.5, resulting in a donor layer program quantization value of 1.5 for donor D in candidate transcription program one. Thus, the core program quantization vector of the validation data is formed according to the donor identifier. The vector corresponding to donor C is a vector that only includes 8.5, and the vector corresponding to donor D is a vector that only includes 1.5.

[0186] The donor layer score index is calculated by calling the fixed weight set without updating the weights. The fixed weight set is called, and the weight of candidate transcription program one is set to 1. The donor layer score index is calculated for each donor identifier in the validation data according to the same preset weighted summation rule. The preset weighted summation rule above is to sum by weight, so the donor layer score index of donor C is 8.5 multiplied by 1 to get 8.5, and the donor layer score index of donor D is 1.5 multiplied by 1 to get 1.5. During this calculation process, the weight set is kept unchanged, so as to ensure that the validation data and the discovery data use the same weight input and the same weighted summation rule.

[0187] The correlation between donor layer scoring indicators and pathological staging information is calculated under the donor layer caliber. The correlation includes at least a correlation statistic based on the consistency of donor ranking. For example, donor layer scoring indicators are sorted from smallest to largest as donor D first, followed by donor C; pathological staging is sorted from lowest to highest as stage II first, followed by stage III, corresponding to donor C first, followed by donor D. The correlation statistic is calculated based on the consistency of the two rankings to quantify the correlation between donor layer scoring indicators and pathological staging information. This correlation statistic is then compared with a preset correlation threshold to output a correlation conclusion. If the comparison result meets the pass condition corresponding to the preset correlation threshold, the correlation conclusion is output as pass; otherwise, the correlation conclusion is output as fail.

[0188] The ability of donor-level scoring indicators to distinguish between late and non-late stages is calculated under the donor-level caliber. For example, if stage 3 is defined as late stage and stage 2 as non-late stage, then donor D is late stage with a score of 1.5, and donor C is non-late stage with a score of 8.5. A threshold scan is performed on the donor-level scoring indicators, for example, using thresholds of 2, 4, 6, and 8 respectively. For each threshold, the distinction between late and non-late stages is calculated, resulting in the area under the curve (AUC). The AUC is compared with a preset discrimination threshold to output the discrimination ability conclusion. If the AUC meets the passing condition corresponding to the preset discrimination threshold, the discrimination ability conclusion is "passed"; otherwise, it is "failed." The above process corresponds to the comparative evaluation caliber of stage association and late stage discrimination ability in an independent external queue.

[0189] The above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A single-cell transcriptional screening method based on donor stability, characterized in that, include: Step S101: Obtain single-cell transcriptional expression dataset, donor identification information and pathological stage information, establish the correspondence between cells and donors based on donor identification information, and screen donors with complete stage annotations to form a donor subset; Step S102: Extract the target cell set based on the donor subset and construct a continuous state representation of each cell; Step S103: Model gene expression dynamics based on continuous state representation, and cluster genes with similar dynamic characteristics to obtain a set of candidate transcription programs; Step S104: Perform staged and stratified sampling of candidate transcription programs; Donor layer aggregation yields donor layer program quantification values; based on the donor set obtained by resampling: the cell expression of candidate transcription programs in the same donor is summarized according to a preset aggregation rule to obtain donor layer program quantification values; a fitting relationship is established between pathological staging information and donor layer program quantification values ​​to obtain evaluation results; effect stability index, directional consistency index, and donor coverage index are calculated, and the SASI score of candidate transcription programs is calculated based on the three indices and the candidate transcription programs are ranked; Step S105: Output the core transcription program set based on the SASI score, and determine the weights by the core transcription program set and its global effect to construct the donor layer scoring index; when external validation is required, obtain a validation single-cell transcription expression dataset from a different source than the single-cell transcription expression dataset, and calculate the donor layer scoring index according to the donor layer caliber to evaluate the correlation and distinguishing ability of the donor layer scoring index with pathological staging information.

2. The single-cell transcription program screening method based on donor stability according to claim 1, characterized in that, Methods for filtering donors with complete phase annotations to form donor subset correspondences include: Single-cell transcriptional expression data and its staging annotation information are obtained, and the pathological stage labels corresponding to each donor identifier are extracted. Missing pathological stage labels are filtered for each donor identifier, and donor identifiers with complete pathological stage labels are retained. The retained donor identifiers are summarized to form a donor subset, and the corresponding cell records are extracted from the staging-associated expression dataset based on the donor subset.

3. The single-cell transcription program screening method based on donor stability according to claim 2, characterized in that, Methods for extracting target cell sets based on donor subsets and constructing continuous state representations of each cell include: Obtain the sub-expression dataset corresponding to the donor subset, extract the cell type annotation field in the sub-expression dataset and compare it with the preset target cell type label set, and retain the cell records whose cell type annotation belongs to the preset target cell type label set to form the target cell set; Gene expression values ​​of the target cell set are converted into normalized expression values ​​according to a preset normalization rule, and feature gene sets are selected according to a preset feature screening rule to construct a feature expression matrix; The distance between cells is calculated based on the feature expression matrix and adjacency relationships are formed according to a preset number of nearest neighbors. The proportion of non-target cells is calculated based on a preset set of non-target marker genes and compared with a preset proportion threshold to remove cell records and obtain the processed target cell set. The corresponding adjacency relationships are regenerated based on the processed target cell set. A random walk transition matrix is ​​constructed based on adjacency relationships. The initial feature score is calculated and compared with a preset initial threshold to determine the initial candidate. The initial cell is determined according to a preset selection rule. The diffusion distance of each cell relative to the initial cell is calculated and the diffusion distance or its order is used as a continuous state representation. The continuous state representation is written back to the cell record corresponding to the processed target cell set to form a continuous state association expression dataset.

4. The single-cell transcription program screening method based on donor stability according to claim 3, characterized in that, Methods for modeling gene expression dynamics based on continuous state representation include: Obtain a continuous state association expression dataset, obtain cell identifiers, continuous state representations and gene expression values, and establish gene-by-gene modeling inputs with continuous state representations as independent variables and gene expression values ​​as dependent variables; In the continuous state associated expression dataset, the number of cells with non-zero gene expression values ​​is counted. The cell count is compared with the preset detection cell count threshold to retain genes that meet the threshold condition. The first preset number of cells are selected according to the cell count to form a candidate modeling gene set. For each gene in the candidate modeling gene set, a smooth fitting model with continuous state representation as input is constructed to output the fitted expression curve. Based on the fitted expression curve, dynamic statistics are calculated, and a preset number of genes are selected according to the dynamic statistics to form a dynamic gene set. A uniform continuous state grid is generated within the range of values ​​represented by the continuous state. The dynamic gene set is input into the smooth fitting model, and the predicted value sequence is output as a trend vector on the uniform continuous state grid. The predicted value sequence is then standardized to meet the preset mean and variance constraints.

5. The single-cell transcription program screening method based on donor stability according to claim 4, characterized in that, Methods for obtaining a set of candidate transcription programs include: Obtain the trend vector of each gene on a uniform continuous state grid, and use the trend vector as the dynamic feature vector of the gene. The dynamic distance is calculated based on the point-by-point difference of the acquisition trend vectors of any two genes on the uniform continuous state grid according to the preset aggregation rules, and a preset number of nearest neighbors is selected for each gene in ascending order of dynamic distance to obtain the set of nearest neighbor genes. Gene adjacency graphs are constructed based on the nearest neighbor gene sets of each gene, where nodes are gene identifiers, edge connections are determined by the nearest neighbor gene sets, and edge weights are obtained by dynamic distance according to preset weight conversion rules. The gene adjacency graph is subjected to graph partitioning and clustering. Multiple gene clusters are output under the constraints of preset resolution parameters, and each gene cluster is defined as a candidate transcription program to obtain a set of candidate transcription programs.

6. The single-cell transcription program screening method based on donor stability according to claim 5, characterized in that, Methods based on donor sets obtained through resampling include: Obtain donor subsets and their pathological staging information, and divide the donor subsets into multiple staging layers based on the pathological staging information; Set the number of resampling iterations and the extraction rules for each iteration. In each iteration, extract donors according to the extraction rules in each acquisition stage and merge them to obtain the donor set obtained by resampling in that iteration.

7. The single-cell transcription program screening method based on donor stability according to claim 6, characterized in that, Methods for summarizing the cellular expression of candidate transcription programs in the same donor according to a preset aggregation rule to obtain the donor layer program quantification value include: For any candidate transcription program, obtain the gene set corresponding to the candidate transcription program in the candidate transcription program set, and obtain the gene expression values, cell identifiers and donor identifiers corresponding to the gene set in the continuous state associated expression dataset; For any donor identifier in the donor set obtained by resampling, all cell records corresponding to the donor identifier are extracted, and the gene expression values ​​of the gene set in the donor are summarized according to the preset aggregation rules to obtain the donor layer program quantification value of the donor identifier on the candidate transcription program.

8. The single-cell transcription program screening method based on donor stability according to claim 7, characterized in that, The method for calculating the SASI score based on multiple resampling includes: Within the donor set obtained in each resampling, with pathological staging information as the independent variable and donor layer procedural quantization value as the dependent variable, a fitting relationship is established using a smooth spline function to obtain the fitting curve under that resampling, and the predicted values ​​of the fitting curve at the highest and lowest stages are obtained. Based on the predicted values, the effect size of the candidate transcription program in this resampling is calculated and the effect direction is recorded, where the effect size is the difference between the predicted values ​​at the highest stage and the lowest stage; the effect size sequence and effect direction are summarized across all resampling iterations, and the effect stability index composed of the standard deviation of the effect size, the direction consistency index composed of the proportion of the dominant effect direction, and the donor coverage index composed of the iterable proportion of the effect size are calculated. The three scores are combined to calculate the SASI score of the candidate transcription program, which is used to rank the candidate transcription programs.

9. The single-cell transcription program screening method based on donor stability according to claim 8, characterized in that, Methods that output a core set of transcription programs based on the SASI score and construct donor layer scoring indicators from the core set of transcription programs include: Obtain the candidate transcription program set and the corresponding SASI score for each candidate transcription program. Sort the programs by SASI score from high to low and select the core transcription program set according to the preset output quantity. Obtain the donor layer program quantization value and summarize the donor layer program quantization value of each core transcription program in the core transcription program set according to the donor identifier to form the core program quantization vector corresponding to the donor identifier. The global effect obtained based on the fitting relationship is used to determine the weight value of each core transcription program in the core transcription program set and associate it with the core transcription program identifier to form a weight set, which is then stored in a fixed manner. For any donor identifier, the core program quantization vector and weight set are synthesized according to a preset weighted aggregation rule to obtain the donor layer scoring index.

10. The single-cell transcription program screening method based on donor stability according to claim 9, characterized in that, Methods for assessing the correlation and discriminative power between donor layer scoring indicators and pathological staging information include: Obtain a validation single-cell transcription expression dataset from a different source than the single-cell transcription expression dataset, and simultaneously obtain the donor identification information and pathological stage information corresponding to the validation single-cell transcription expression dataset, so that the cell identifier and the donor identifier are associated and mapped to the pathological stage information; The set of core transcription programs and their corresponding weights that have been identified and solidified in the discovery dataset are invoked. The core transcription programs are aggregated by donor layer in the validation single-cell transcription expression dataset to obtain the core program quantization vector of each donor. Under the condition that the weights are not updated, the donor layer score index is calculated for each donor. Under the donor layer caliber, the correlation statistics between the donor layer scoring index and pathological staging information are calculated and compared with the preset correlation threshold to output the correlation conclusion. The area under the curve statistics of the donor layer scoring index are also calculated and compared with the preset discrimination threshold to output the discrimination ability conclusion.

Images