An integrated method for tumor single-cell transcriptome metaprogram identification and functional annotation

Abstract

This invention discloses an integrated method for identifying and functionally annotating metaprograms in tumor single-cell transcriptomes, belonging to the fields of bioinformatics, tumor biology, and single-cell transcriptome sequencing. This invention aims to address the problems of unstable metaprogram identification, significant technical artifact interference, and fragmented functional annotation in existing technologies. Based on multi-rank nonnegative matrix factorization, it systematically extracts co-expression modules at different scales from the tumor single-cell expression matrix. Robustness screening ensures the consistency and reliability of the obtained metaprograms under multiple factorizations and cross-sample conditions. Furthermore, it performs systematic functional analysis of the metaprograms to reveal the multi-dimensional functional states of tumor cells. This method improves the stability and biological interpretability of tumor metaprogram identification, providing a new technical means for tumor heterogeneity research and the discovery of precision therapeutic targets.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of single-cell transcriptome sequencing, tumor biology and bioinformatics, and particularly relates to an integrated method for metaprogram identification and functional annotation of tumor single-cell transcriptome data based on non-negative matrix factorization. BACKGROUND

[0002] Tumor tissues are driven by multiple factors such as genetic mutations, epigenetic abnormalities and microenvironmental pressure during tumorigenesis and progression, and the internal tumor cells exhibit significant state diversity and expression heterogeneity. Different tumor cell subpopulations have obvious differences in proliferation activity, differentiation degree, metabolic demand, stress response and invasion ability, and these differences not only affect the biological behavior of tumors, but also directly determine the treatment sensitivity and drug resistance mechanism of patients. Therefore, systematic analysis of the functional state structure of tumor cells has important scientific significance and clinical value for understanding the rules of tumor occurrence and development, identifying key regulatory modules and discovering potential therapeutic targets.

[0003] With the development of single-cell RNA sequencing (scRNA-seq) technology, tumor cell populations can be isolated and accurately identified from complex tumor tissues, and high-resolution gene expression matrices can be obtained. By analyzing tumor cell subsets independently, the expression interference of tissue microenvironment cells can be avoided, and the regulatory characteristics of tumor cells themselves can be more accurately characterized. With the help of this technology, researchers can analyze the functional state differences within the tumor at the cellular level and further explore the potential co-expression modules that support these states. However, in multi-source tumor cell data across patients and samples, how to stably identify functional modules with biological consistency and interpretability remains a key challenge for existing research.

[0004] Existing analytical strategies for identifying functional co-expression modules in tumor cell expression matrices primarily rely on decompositional modeling of statistical relationships between genes and potential expression structures. Traditional methods based on statistical features or correlations, such as highly variable genes (HVGs), differentially expressed genes (DEGs), and weighted gene co-expression network analysis (WGCNA), acquire potential modules through gene variance or correlation structures. However, the prevalent high noise, high sparsity, and sequencing depth variations in tumor cell expression matrices lead to highly unstable correlation structures between genes. Differences in tumor cell expression backgrounds among different patients further weaken the reproducibility of co-expression structures, making it difficult for these methods to reliably extract conserved functional modules across samples. Furthermore, statistical correlations alone are insufficient to reflect the multi-dimensional and multi-level state structure within tumor cells, making it difficult to systematically characterize the core gene combinations driving tumor cell behavior.

[0005] In contrast, matrix factorization methods attempt to extract latent structural factors from the overall expression matrix, representing more complex co-expression patterns. Principal Component Analysis (PCA) and Independent Component Analysis (ICA) typically produce linear or independent decomposition results, lacking clear biological interpretability and significantly affected by noise and sparsity, making it difficult to obtain stable functional modules. Non-negative matrix factorization (NMF), due to its alignment with the non-negative nature of gene expression data, is considered more suitable for extracting interpretable co-expression patterns. However, existing NMF-based studies often employ fixed decomposition conditions with a single rank value, resulting in the capture of expression structures at only specific scales; the decomposition results are sensitive to initialization and parameter changes, lacking a systematic robustness assessment mechanism. Furthermore, technical factors such as differences in tumor cell sequencing depth and cellular complexity bias can easily create systematic artifacts during decomposition, interfering with the identification of true biological modules. The lack of coherent integration between analytical steps such as module extraction, artifact removal, and annotation further limits the interpretability and cross-sample reproducibility of modules.

[0006] In summary, existing technologies are insufficient to simultaneously meet the needs of multi-scale mining of tumor cell metaprograms, cross-sample robustness verification, artifact removal, and systematic functional annotation. An effective method for integrated metaprogram identification and functional analysis targeting single-cell transcriptome expression matrices in tumor cells is still lacking. Therefore, there is an urgent need to develop new analytical strategies to comprehensively and stably reveal the functional state structure within tumor cells and to uncover core regulatory modules that are conserved across patients. Summary of the Invention

[0007] This invention addresses the problems of unstable metaprogram identification, severe technical artifact interference, insufficient cross-sample reproducibility, and fragmented functional annotation processes in existing tumor single-cell transcriptome data. It provides an integrated method for metaprogram identification and functional annotation in tumor single-cell transcriptome data. Through multi-rank matrix factorization, robustness screening, artifact factor removal, and systematic annotation, this invention achieves multi-level structural analysis of the tumor cell expression matrix, stably extracting biologically significant core metaprograms. This provides effective technical support for tumor cell functional status research and potential target discovery.

[0008] This invention provides an integrated method for tumor single-cell transcriptome meta-programmatic identification and functional annotation, comprising the following steps:

[0009] 1) Construction and data preprocessing of tumor cell expression matrix

[0010] Tumor cells were screened from single-cell transcriptome data, and a tumor cell gene expression matrix was constructed. The expression matrix underwent standardization, quality control, and low-expression gene filtering to eliminate the influence of sequencing bias and cell quality differences, providing high-quality input for subsequent meta-process analysis.

[0011] 2) Nonnegative matrix decomposition based on multi-rank values

[0012] Multiple rank values ​​were set on the preprocessed tumor cell expression matrix for nonnegative matrix factorization to obtain candidate gene programs at different scales. Multi-rank factorization can comprehensively cover the coarse-to-fine granular co-expression structure of tumor cells, which is an important prerequisite for obtaining multi-scale metaprograms.

[0013] 3) Robustness screening and consistency assessment of candidate meta-procedures

[0014] For candidate gene programs obtained under various rank conditions, their reproducibility under different initialization conditions, different samples, and different decomposition parameters is evaluated from multiple perspectives. Programs that appear stably under multiple conditions are screened to form a robust metaprogram candidate set.

[0015] 4) Construction and integration of cross-sample common meta-routines

[0016] Cross-sample consistency analysis was performed on real metaroutines screened from different patients and samples to identify common metaroutines that could be repeated and stably appear in a multi-patient context, and these metaroutines were integrated to construct a tumor cell metaroutines library.

[0017] 5) Identification and removal of technical artifacts

[0018] An artifact identification system was constructed to address technical factors such as sequencing depth differences, cell complexity shifts, and expression sparsity, and artifact factors were detected in candidate metaprograms. Programs significantly related to technical factors were removed, retaining only biological programs that truly reflect the functional state of tumor cells.

[0019] 6) Systematic functional annotation and biological analysis of metaprograms

[0020] The similarity of the procedures obtained after step 5) was recalculated and clustered, retaining meta-procedures that represent at least two patients. Enrichment analysis, pathway annotation, and biological feature analysis were performed on the final meta-procedures to identify their corresponding functional systems. Combining tumor biology, the potential roles of the meta-procedures in tumorigenesis, progression, and treatment response were further inferred.

[0021] Step 1) above, the construction and data preprocessing of the tumor cell expression matrix, specifically includes:

[0022] ① Input data requirements: The input data is single-cell transcriptome data, in the format of a Seurat object, containing a unique molecular identifier (UMI) counting matrix, i.e., the counts matrix.

[0023] ② Sample screening requirements: Ensure that each sample contains no fewer than 50 cells to guarantee the reliability of the analysis results.

[0024] ③ Data standardization: Divide the gene UMI count by the total number of UMIs in the cell, multiply by 100,000, add 1, and finally perform log2 transformation.

[0025] ④ Gene filtering: Only genes with a log2 expression level greater than 3.5 in at least 2% of cells are retained, and the top 7000 genes with the highest expression levels are selected.

[0026] ⑤ Data digitization: For each gene, subtract the average expression value of that gene in all cells from the gene expression value in each cell, and set all negative values ​​to 0.

[0027] Step 2 above is based on multi-rank nonnegative matrix factorization (NMF), where:

[0028] The snmf / r factorization algorithm from the NMFR package was used, with a rank value ranging from 4 to 9. The algorithm was run at least 100 times for each sample, and the result with the smallest error was retained to obtain the gene feature matrix (w_basis) and the cell feature matrix (h_coef).

[0029] Step 3) above, the robustness screening and consistency evaluation of candidate meta-procedures, specifically includes:

[0030] ① Robustness within the tumor: If a program repeatedly appears in multiple rank values ​​in the same tumor sample, and at least 70% (≥35) of the first 50 characteristic genes of any two decompositions are consistent, then the program is considered to be stable within the tumor.

[0031] ② Cross-sample consistency: If a program has at least 20% (≥10 genes) overlap in the first 50 characteristic genes with programs of other tumor samples, then the program is considered to have cross-sample reproducibility.

[0032] ③ Non-redundancy within the tumor: In the same tumor sample, programs that have passed the robustness screening are sorted from high to low according to their gene overlap with other sample programs, and the most representative programs are retained first; other programs with gene overlap of more than 20% (≥10 genes) with the retained programs are removed to ensure that only non-redundant metaprograms are ultimately retained for the tumor.

[0033] Step 4) above, the construction and integration of cross-sample common meta-routines, includes:

[0034] ①Program similarity analysis: Calculate the Jaccard similarity of the first 50 characteristic genes of each NMF program to generate a similarity matrix between NMF programs.

[0035] ② Metaprogram Clustering and Construction: Based on the similarity matrix between NMF programs, select an appropriate hierarchical clustering method and divide the clustering tree into an appropriate number of metaprograms. The clustering method can be flexibly selected from Ward.D2 or Average depending on the characteristics of the data.

[0036] ③ Metaprogram heatmap display: Draw a similarity heatmap with metaprogram annotations to intuitively display the program similarity in different samples and the distribution of common metaprograms across samples.

[0037] ④ Screening of the top 50 genes in each metaprogram: Statistically analyze the frequency of genes in each metaprogram, sort the genes according to their average ranking in different programs, and screen the top 50 core genes in each metaprogram based on frequency and ranking.

[0038] Step 5) above, the identification and removal of technical artifacts, includes:

[0039] ① Cell complexity calculation: For the expression matrix of each tumor sample, the number of non-zero expressed genes in a single cell is used as the cell complexity index, and the complexity vector of each cell is calculated column by column.

[0040] ②Program complexity correlation assessment: For each robust NMF program, firstly, extract the corresponding tumor sample and the cell score vector of the program in the sample according to the program name, and then calculate the Pearson correlation coefficient with the complexity vector of each cell in the same sample to obtain the program cell complexity correlation matrix.

[0041] ③ Visual analysis of correlation trends: Arrange the correlation coefficients of all programs in the order of the programs in the heatmap, and plot a scatter plot with the program number in the heatmap in step 4) as the horizontal axis and the correlation coefficient as the vertical axis. Use local weighted regression (LOESS) to fit a smooth curve to intuitively show the overall correlation trend and abnormal patterns between different programs and cell complexity.

[0042] ④ Identification and Removal of Technical Artifacts: Based on the strength of the correlation between the program and cell complexity, the shape of the correlation trend curve, and the characteristics of the genes enriched by the program, abnormal programs suspected of being driven by sequencing depth bias, cell complexity deviation, or other technical factors are manually identified. Simultaneously, the first 50 core genes of each metaprogram are checked for excessive ribosome and mitochondrial genes; excessive enrichment indicates technical noise. Combining the above analysis, abnormal programs suspected of being driven by technical factors are manually removed, retaining metaprograms that are not significantly correlated with cell complexity and can accurately reflect the biological state of tumor cells.

[0043] Step 6) above, systematic functional annotation and biological analysis of the meta-program, wherein:

[0044] ① Manual verification of single-gene functions: Combining existing literature and the GeneCards database, the first 50 genes of each metaprogram obtained in step 5 were manually verified. By reviewing the functions of the genes, their known roles in tumor biology, cell function, etc., were confirmed to ensure that the gene characteristics of the metaprogram were consistent with their biological significance.

[0045] ② Pathway enrichment analysis: For the first 50 genes of each metaprogram, pathway enrichment analysis was performed using databases such as KEGG, GO, Hallmark50, and Reactome. This helps identify the biological pathways, functional categories, and regulatory mechanisms involved in the genes, providing detailed biological annotations for the metaprogram and revealing their associations with specific biological processes, cell types, or disease states.

[0046] ③ Gene overlap analysis: The first 50 genes of each metaprogram were compared with 41 metaprograms published by Gavish et al. (Gavish, A., et al. Hallmarks of transcriptional intratumour heterogeneity across a thousand tumours. Nature 618, 598–606 (2023). https: / / doi.org / 10.1038 / s41586-023-06130-4) using Jaccard similarity calculation to assess gene overlap. Hypergeometric and permutation tests were used to evaluate the statistical significance of gene overlap, verify the biological relevance of the metaprograms to those in known literature, and ensure their consistency.

[0047] Compared with existing methods, the present invention has the following advantages:

[0048] 1) The multi-rank decomposition and robustness screening system established in this invention significantly improves the stability and cross-sample reproducibility of tumor cell meta-program identification.

[0049] 2) A systematic technical artifact identification and removal mechanism effectively reduces the interference of sequencing bias on the meta-program, making the extraction results more biologically reliable.

[0050] 3) This invention integrates metaprogram extraction, artifact removal, and functional annotation into a unified process, which improves parsing efficiency and enhances the coherence and interpretability of module parsing.

[0051] 4) The obtained meta-program can stably characterize the functional state structure inside tumor cells, providing a reliable technical foundation for studying tumor heterogeneity, identifying core regulatory modules, and exploring targets related to precision diagnosis and treatment of tumors. Attached Figure Description

[0052] Figure 1 A flowchart of the integrated method for tumor single-cell transcriptome metaprogram identification and functional annotation of the present invention.

[0053] Figure 2 The gene feature matrix (top figure) and cell feature matrix (bottom figure) obtained by nonnegative matrix decomposition in this embodiment of the invention.

[0054] Figure 3 This invention provides a hierarchical clustering dendrogram based on inter-program Jaccard similarity.

[0055] Figure 4 The heatmap in this embodiment of the invention is based on the Jaccard similarity between programs.

[0056] Figure 5 The scatter plot in this embodiment of the invention is based on the correlation between program and cell complexity.

[0057] Figure 6 Gene overlap analysis diagram between tumor metaprograms and known metaprograms in this embodiment of the invention.

[0058] Figure 7 A tumor meta-procedure heatmap with robust biological significance in this embodiment of the invention. Detailed Implementation

[0059] The following is a more detailed description of the present invention. The parameters and specific implementation details are used to explain the feasibility and implementation effects of the present invention, and do not constitute a limitation of the present invention.

[0060] This example uses single-cell transcriptome data from 21 patients with lung adenocarcinoma for demonstration. Figure 1 The flowchart of the method of the present invention is shown below, and the experimental steps and results are described in detail below:

[0061] 1) Construction and data preprocessing of tumor cell expression matrix

[0062] First, single-cell transcriptome data from lung adenocarcinoma patients were retrieved. The data contained only tumor cells identified using the inferCNV package and included a tumor cell count matrix. Each sample was screened to ensure it contained at least 50 cells to improve the reliability of the analysis results. The total UMI count for each cell was calculated, and the UMI count for each gene was divided by the total UMI count for that cell to obtain its expression percentage. This percentage was then multiplied by 100,000 for standardization. Finally, the standardized data underwent a log2 transformation (adding 1 to avoid zero values) to obtain the standardized and transformed gene expression matrix. Genes with expression values ​​greater than 3.5 in at least 2% of cells were further screened, and the top 7000 genes by expression level were retained. The expression value of each gene was then centered by subtracting its average expression value across all cells. All negative values ​​were set to 0 to ensure the data met NMF input requirements.

[0063] 2) Nonnegative matrix decomposition based on multi-rank values

[0064] Using the snmf / r factorization algorithm from the NMFR package, each sample is factored with a rank ranging from 4 to 9. Each sample will ultimately yield two matrices: one is a gene feature matrix (w_basis) containing 7000 genes and 39 metaprograms. Figure 2 The upper middle figure), and the other is the cell feature matrix (h_coef) containing the number of samples and 39 meta-programs ( Figure 2(See the lower half of the image). These two matrices together constitute the expression characteristics of tumor cells. Figure 2 As shown, for 21 samples, each sample is decomposed into 39 programs, resulting in a total of 819 programs.

[0065] 3) Robustness screening and consistency assessment of candidate meta-procedures

[0066] All 819 programs were screened to retain those that were stable across multiple samples and different rank values. First, the top 50 characteristic genes of each program were extracted. The overlap of the top 50 characteristic genes of the same program in any two decompositions was calculated. If the number of overlapping genes reached 70% (i.e., at least 35 genes were identical), the program was considered stable in that tumor sample. Further checks were made on the consistency of these robust programs across samples by calculating their gene overlap with corresponding programs in other samples. If a program had at least 20% gene overlap (i.e., at least 10 genes) across different samples, it was considered to have good reproducibility across different samples. Finally, redundancy was eliminated from programs that had passed the robustness screening within the same tumor sample. If a program had more than 20% gene overlap with a retained program (i.e., at least 10 genes overlap), it was removed to ensure that the final retained programs were non-redundant. Ultimately, we obtained 131 robust programs for subsequent analysis.

[0067] 4) Construction and integration of cross-sample common meta-routines

[0068] For all 131 robust programs, the top 50 feature genes were extracted. A similarity matrix between programs was generated by calculating the Jaccard similarity of the top 50 feature genes for each program. The Average method was used to calculate the average distance between each pair of elements for clustering, and then... Figure 3 The tree diagram shown indicates that 15 metaroutes are a suitable number. To visually demonstrate the similarity of programs across different samples and the distribution of common metaroutes across samples, we used the ComplexHeatmap package to create a similarity heatmap with metaroute annotations, as shown below. Figure 4 As shown.

[0069] 5) Identification and removal of technical artifacts

[0070] The number of genes with expression values ​​greater than zero in each cell is counted to obtain the complexity of each cell. Next, for each robust program, the cell feature matrix (h_coef) of the program in each sample is extracted, and Pearson correlation analysis is performed between the program's score vector and the complexity vector of each cell in the sample to obtain the program cell complexity correlation matrix.

[0071] Based on the order of the programs in the similarity heatmap, the correlation coefficients of all programs are arranged, with the program number on the horizontal axis and the correlation coefficient on the vertical axis. A scatter plot is then drawn between all programs and cell complexity, and Loess is used to smooth the scatter plot to generate a trend curve. Programs with abnormally fluctuating curves are marked, such as... Figure 5 As shown, it is an abnormal program driven by technical deviation.

[0072] The first 50 characteristic genes of each obtained metaprogram are counted. If there are more than 3 ribosomal genes or more than 3 mitochondrial genes, it is considered that the metaprogram is closely related to technical factors such as cell quality and sequencing depth, reflecting technical noise.

[0073] The two types of anomalous programs and metaprograms were manually removed, resulting in 72 robust and biologically significant programs.

[0074] 6) Systematic functional annotations and biological analysis of meta-programs, including:

[0075] Step 4) was repeated for 72 robust and biologically significant procedures to determine that 10 metaroutines were appropriate. Metaroutines representing at least two patients were annotated, and the first 50 genes of each metaroutine were functionally queried using GeneCards. Jaccard similarity was calculated against 41 known metaroutines published by Gavish et al. (Gavish, A., et al. Hallmarks of transcriptional intratumour heterogeneity across a thousandtumours. Nature 618, 598–606 (2023). https: / / doi.org / 10.1038 / s41586-023-06130-4). The statistical significance of gene overlap was assessed using hypergeometric and permutation tests, such as... Figure 6 As shown. The final result is as follows: Figure 7 The six tumor metaroutines shown have the following biological meanings: oxidative stress, secretory cell type I, secretory cell type II, mucin epithelium, senescent epithelium, and mesenchymal transition epithelium.

[0076] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made in accordance with the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. An integrated method for tumor single-cell transcriptome meta-programmatic identification and functional annotation, comprising the following steps: 1) Construction and data preprocessing of tumor cell expression matrix: Tumor cells were screened from single-cell transcriptome data, and a tumor cell gene expression matrix was constructed. Then, the expression matrix was standardized, quality controlled, and low-expression genes were filtered to eliminate the influence of sequencing bias and cell quality differences. 2) Non-negative matrix factorization based on multi-rank values: Multiple rank values ​​are set for the preprocessed tumor cell expression matrix, and non-negative matrix factorization is performed to obtain candidate gene programs at different scales. 3) Robustness screening and consistency assessment of candidate metaprograms: For candidate gene programs obtained under each rank condition, evaluate their reproducibility under different initialization conditions, different samples and different decomposition parameters, screen programs that appear stably under multiple conditions, and form a robust metaprogram candidate set. 4) Construction and integration of common metaroutines across samples: Perform cross-sample consistency analysis on real metaroutines screened from different patients and samples, identify common metaroutines that can be repeated and stably appear in a multi-patient context, and integrate them to construct a tumor cell metaroutines library. 5) Identification and elimination of technical artifacts: Construct an artifact identification system targeting technical factors including sequencing depth differences, cell complexity shifts and expression sparsity. Perform artifact factor detection on candidate metaprograms, eliminate programs that are significantly related to technical factors, and retain only programs that truly reflect the functional state of tumor cells. 6) Systematic functional annotation and biological analysis of meta-programs: The similarity of the programs obtained after step 5) is recalculated and clustered, retaining meta-programs that represent at least two patients; enrichment analysis, pathway annotation and biological feature analysis are performed on the finally obtained meta-programs to identify their corresponding functional systems.

2. The method as described in claim 1, characterized in that, Step 1) Screen tumor cells from single-cell transcriptome data and extract the UMI count matrix; retain no less than 50 cells in each sample, calculate the total UMI count for each cell, divide the gene UMI count by the total UMI count of the cells, multiply by 100,000, add 1, and then perform log2 transformation; screen for genes with a log2 expression value greater than 3.5 in at least 2% of cells, and retain the top 7000 genes with the highest expression levels; center each gene according to its average expression in all cells and set negative values ​​to 0 to construct a normalized expression matrix for non-negative matrix factorization.

3. The method as described in claim 2, characterized in that, The data processed in step 1) is single-cell transcriptome data in Seurat object format containing a UMI counting matrix.

4. The method as described in claim 2, characterized in that, Step 2) Use the snmf / r factorization algorithm of the NMFR package to perform at least 100 factorizations on each sample within the rank range of 4 to 9, select the result with the smallest error, and obtain a gene feature matrix containing a 7000-gene × 39 program and a cell feature matrix containing a cell number × 39 program.

5. The method as described in claim 4, characterized in that, In step 2), each sample obtained by nonnegative matrix factorization contains 39 programs, and n samples yield a total of n×39 programs.

6. The method as described in claim 1, characterized in that, Step 3) includes: (3a) Intratumor robustness screening: If a program appears repeatedly in the same sample under different rank values ​​and at least 70% of its first 50 genes are identical, it is determined to be a robust program; (3b) Cross-sample consistency screening: If a program has at least 20% overlap in the first 50 genes of other sample programs, it is considered to have cross-sample reproducibility. (3c) Non-redundancy screening: Robust programs in the same sample are sorted from high to low cross-sample overlap, and only the most representative programs are retained, while redundant programs with more than 20% overlap with their genes are removed.

7. The method as described in claim 1, characterized in that, In step 4), the top 50 genes of all programs selected in step 3) are extracted, and the Jaccard similarity between programs is calculated to form a similarity matrix; a clustering tree is constructed based on hierarchical clustering and the number of metaprograms is determined. The frequency and average ranking of genes within each metaprogram were statistically analyzed, and the top 50 core genes of each metaprogram were selected.

8. The method as described in claim 1, characterized in that, In step 5), the number of non-zero genes in each cell is calculated as the cell complexity vector; for each program, the Pearson correlation coefficient is calculated between its cell score vector and cell complexity vector; a scatter plot is plotted with the program number as the horizontal axis and the correlation coefficient as the vertical axis, and Loess smoothing is performed to identify abnormal programs driven by cell complexity; the number of ribosomal genes and mitochondrial genes in the top 50 genes of the program is counted, and if the number of each exceeds 3, it is determined to be a technical artifact; technical noise programs are eliminated by combining the above indicators.

9. The method as described in claim 1, characterized in that, In step 6), the Jaccard similarity of the programs retained in step 5) is recalculated and clustered to select metaprograms from at least two patients; GeneCards single-gene function checks are performed on the top 50 genes of each metaprogram, and pathway enrichment analysis is performed based on the KEGG, GO, Hallmark50 and Reactome databases; and Jaccard similarity calculation, hypergeometric test and permutation test are performed between the metaprogram genes and published metaprograms to assess biological consistency.

Images