Method and system for differential detection and typing of immune microenvironment of medulloblastoma

Abstract

The application belongs to the technical field of biomedicine, and specifically discloses a method and system for detecting and typing the immune microenvironment of medulloblastoma. The method collects gene expression data of a sample suffering from medulloblastoma, introduces a protein interaction network, combines with pseudo-single-cell expression data, performs unsupervised clustering through similarity network fusion analysis, determines different immune microenvironment subtypes in medulloblastoma, and then performs machine learning through lasso to screen biological characteristics of different immune microenvironment subtypes. Finally, single-cell omics analysis is performed. Through the verification of immune microenvironment typing and marker screening, the technical solution can obtain the immune characteristics and potential therapeutic targets of medulloblastoma, and provide new ideas and methods for the treatment and prognosis evaluation of diseases.

Description

Technical Field

[0001] This invention belongs to the field of biomedical technology and relates to a method and system for detecting and classifying immune microenvironment differences in medulloblastoma. Background Technology

[0002] Medulloblastoma (MB) is the most common childhood brain tumor in the posterior fossa, accounting for approximately 25% of all childhood brain tumors, and is highly aggressive. Transcriptomic, genomic, epigenomic, and proteomic analyses have shown that human MB exhibits high heterogeneity.

[0003] These tumors were further subdivided into four main subgroups with different outcomes and molecular characteristics: WNT (Wingless-related integration site, the WNT protein family is named after the "Wingless" gene in fruit flies and the "Int-1" gene in mice, which play important roles in morphogenesis, cell proliferation and differentiation), Sonic Hedgehog (SHH), Group 3 and Group 4. WNT and SHH are usually mutated in the WNT and SHH pathway genes, respectively. GP3 subgroup tumors usually show MYC (Myelocytomatosis viral oncogene homolog, MYC is an important proto-oncogene that plays a key role in cell proliferation, differentiation and apoptosis) amplification and have the worst clinical outcomes. The genetic drivers of GP4 are still unclear.

[0004] While traditional treatments have made some progress, the outcomes remain unsatisfactory for some high-risk patients, with survival rates and prognoses still being unfavorable. Therefore, exploring new treatment strategies and predictive biomarkers is particularly important.

[0005] Tumor heterogeneity refers to the composition of tumor tissue by cellular subpopulations with different expression profiles or biological functions, leading to differences in growth rate, invasiveness, metastasis, and drug sensitivity. This is a prominent characteristic of malignant tumors. This heterogeneity is a major cause of tumor recurrence, metastasis, and drug resistance, directly impacting clinical treatment outcomes. Therefore, analyzing tumor heterogeneity has significant clinical diagnostic value, and in-depth research into the formation and regulatory mechanisms of tumor heterogeneity will provide theoretical support for precision targeted therapy of tumors.

[0006] Currently, there are two main theories regarding the causes of tumor heterogeneity: the clonal selection theory and the tumor stem cell theory. The clonal selection theory posits that tumor heterogeneity originates from the continued mutation of individual tumor cell populations during their development. In the process of adapting to the environment, only those cell populations capable of survival and reproduction are preserved, ultimately leading to tumor heterogeneity. The tumor stem cell theory, on the other hand, suggests that tumors are composed of a small group of self-renewing tumor stem cells and cells with varying degrees of differentiation, thus exhibiting heterogeneity in different cell types. Although the two theories offer different explanations, both emphasize the crucial role of the tumor microenvironment (TME) in the formation of tumor heterogeneity.

[0007] The tumor microenvironment refers to the internal environment in which tumor cells develop and progress, including the tumor cells themselves, immune cells, stromal cells, and the surrounding intercellular matrix, microvessels, and infiltrating biomolecules. In recent years, increasing research has shown that the tumor immune microenvironment plays a crucial role in tumor occurrence, development, and treatment response. The tumor microenvironment not only influences tumor cell growth and metastasis but is also closely related to treatment efficacy and patient prognosis. Personalized treatment based on the tumor microenvironment is receiving growing attention, especially immunotherapy, which has been widely used clinically to treat various cancers with significant results. However, our understanding of the tumor microenvironment in medulloblastoma remains insufficient.

[0008] Therefore, a deeper understanding of the cellular composition of the medulloblastoma microenvironment, the interaction between medulloblastoma cancer cells and other cells in the microenvironment, and the impact of medulloblastoma cancer cells on the microenvironment will help discover biomarkers and potential therapeutic targets for the development of medulloblastoma. This is of great significance for understanding the pathogenesis of this tumor, predicting the progression of the disease in patients, and developing individualized treatment plans.

[0009] However, research on the immune microenvironment of medulloblastoma is still in its early stages, especially regarding the immune characteristics and clinical significance of different immune subtypes of medulloblastoma. Summary of the Invention

[0010] The purpose of this invention is to provide a method and system for detecting and classifying the immune microenvironment differences in medulloblastoma, to obtain the immune characteristics and potential therapeutic targets of medulloblastoma, and to provide new ideas and methods for the treatment and prognosis assessment of the disease.

[0011] To achieve the above objectives, the basic solution of the present invention is: a method for detecting and classifying immune microenvironment differences in medulloblastoma, comprising the following steps:

[0012] S1: Collect gene expression data from samples of patients with medulloblastoma and perform pseudo-single-cell calculations.

[0013] S2, using Netclass to generate a protein interaction network as an adjacency matrix, introduces gene expression data into the protein interaction network, and then combines pseudo-single-cell gene expression data, and uses similarity network fusion analysis to perform unsupervised clustering to determine different immune microenvironment subtypes in medulloblastoma;

[0014] S3 uses lasso for machine learning to screen biological features of different immune microenvironment subtypes;

[0015] S4 validates the immune microenvironment typing and biological feature screening through single-cell omics analysis.

[0016] The working principle and beneficial effects of this basic scheme are as follows: This technical scheme obtains the heterogeneity of the tumor immune microenvironment of different subtypes of medulloblastoma by classifying the microenvironment of different subtypes of MB patients, and conducts in-depth research on the differences in their microenvironment landscape by combining single-cell sequencing data visualization.

[0017] Medulloblastoma (MB) was classified into different subtypes with distinct biological and clinical characteristics, and genes showing prognostic differences among these subtypes were screened. Different subtypes of medulloblastoma possess distinct immune microenvironment subtypes, each with unique markers and functions, playing a unique role in tumor development and progression, and leading to drastically different outcomes for patients. Through immune microenvironment typing and marker screening, the immune characteristics and potential therapeutic targets of medulloblastoma were revealed, providing new insights and methods for the treatment and prognostic assessment of this disease.

[0018] Furthermore, step S1 specifically includes:

[0019] Based on the ssGSEA pseudo-single-cell analysis method, gene expression data are compared with a given gene set to calculate the enrichment score of the gene set in each sample, simulating the composition of immune cells in the tumor immune microenvironment.

[0020] Based on the netClass analysis method, gene expression data and protein network topology are integrated using information from protein-protein interaction (PPI) networks.

[0021] The pseudo-single-cell results obtained by the ssGSEA method were combined with gene expression data to perform similar network fusion (SNF) analysis on the samples.

[0022] The survival, ggplot2, ESTIMATE, MCPcounter, and msigdbr packages in R were used to visualize and analyze the survival data of different immune subtypes of MB patients, evaluate the abundance of the tumor microenvironment and different cell types in the tumor microenvironment, and evaluate the gene set data related to biological processes.

[0023] Use the glmnet function in the glmnet package to perform multi-class lasso regression analysis on gene expression data;

[0024] Based on gene expression data, DEG analysis of differentially expressed genes between groups was performed using the R package limma (3.56.2) to obtain the differences in gene expression among sample groups; bubble charts were drawn using the R package GOplot (R 4.0.2) to observe pathway enrichment.

[0025] To identify different immune microenvironment subtypes and their corresponding biological characteristics in medulloblastoma for subsequent analysis and processing.

[0026] Furthermore, based on the netClass analysis method, gene expression data and protein network topology are integrated using information from protein-protein interaction networks. The specific steps are as follows:

[0027] Data on protein-protein interaction networks are obtained from public databases or literature, and a network containing the interaction relationships between proteins is constructed as an adjacency matrix.

[0028] Calculate the diffusion kernel: Accepts the adjacency matrix as input and calculates the diffusion kernel matrix according to the specified parameters (p and a). The diffusion kernel matrix is ​​generated through iterative calculation. When calculating the diffusion kernel matrix, the parameter is.adjacency is used to determine whether the input is an adjacency matrix or a Laplacian matrix, and different calculation paths are selected accordingly.

[0029] Use the igraph library to perform graph operations, including creating undirected graphs and calculating the Laplacian matrix;

[0030] When calculating the diffusion kernel matrix and subsequent eigenvector multiplication, matrix operations are used to filter out the intersection with another dataset, and then the diffusion kernel is calculated using a portion of the data from the intersection.

[0031] netClass is used to analyze gene expression data and, in conjunction with the topology of protein-protein interaction networks, to identify biologically significant gene patterns.

[0032] Furthermore, the method for performing SNF analysis on the samples using the pseudo-single-cell results obtained through the ssGSEA method, combined with gene expression data, is as follows:

[0033] Use the complete expression matrix after removing pseudogenes and the ssGSEA matrix as input;

[0034] Using the SNFtool R package (v2.2.0), with the number of neighbors K=30, Gaussian kernel parameter alpha=0.5, and number of iterations T=10, the spectral clustering implemented by the SNF tool package is run on the SNF fused similarity matrix to obtain the most suitable grouping;

[0035] SNF analysis was performed in different subtypes: In the SHH (sample size n=233) subtype, three immune microenvironment subgroups were obtained when k=3;

[0036] In the GP3 (n=144) subtype, two immune microenvironment subpopulations were obtained when k=2;

[0037] In the GP4 (n=326) subtype, three immune microenvironment subgroups were obtained when k=3.

[0038] Netclass analysis combined with SNF analysis can classify the tumor microenvironment of different subtypes of MB tumors, and the operation is simple.

[0039] Furthermore, using the `survival`, `ggplot2`, `ESTIMATE`, `msigdbr`, and `MCPcounter` packages in R, we performed visualization analysis on the survival data of MB patients, assessed the tumor microenvironment based on gene expression data, obtained gene set data related to biological processes from the species' genomics database, and evaluated the abundance of different cell types in the tumor microenvironment. The steps were as follows:

[0040] S11, Prepare the data required for survival analysis, including the observation time (futime) and event status (fustat) for each sample, as well as the group variable for grouping;

[0041] S12, use the survfit() function to fit the survival curve, calculate the survival curve based on the specified observation time, event state and grouping variable group;

[0042] S13, use the ggsurvplot() function to create a survival curve chart. In the chart, set the following parameters: do not display confidence intervals (conf.int=F), display risk table (risk.table=T), do not add total patient survival curves (add.all=F), customize the color palette (palette="Dark2"), and adjust the title, axis titles, font size and style, as well as the range of the axes.

[0043] S21. Prepare the gene expression data file. Use the outputGCT() function to convert the data into the input format required by the ESTIMATE package and save it as a new file (in.gct.file).

[0044] S22, use the filterCommonGenes() function to filter out common genes from the original data and save them as a new input file, use the estimateScore() function to calculate the ESTIMATE score of the tumor sample, and save the score results as a GCT format file (out.score.file);

[0045] S23, then use the plotPurity() function to visualize the purity score of the sample, use the read.table() function to read the score result file and convert it into a data frame format, and finally save the result as a text file (ESTIMATE_score.txt);

[0046] S31, use the msigdbr() function to retrieve the gene set data of C5 category from the genomics database of the species and save it in the GO_df_all data frame;

[0047] The dplyr::select() function was used to select the desired columns (gs_name, gene_symbol, gs_exact_source, and gs_subcat), and unwanted data was filtered out according to the subcategory (gs_subcat) to obtain the final GO gene set data (GO_df);

[0048] S32, use the split() function to group the gene set data according to the GO entry name gs_name to obtain a gene list (go_list);

[0049] The gene expression data (rt1) to be analyzed is subjected to gene set variability analysis (GSVA) ​​using the gsva() function. The gene set variability score (gsva_mat) of each GO entry is calculated and the results are saved in a text file (gsva_hallmark.txt).

[0050] S33, set parameters: use Gaussian kernel density estimation function (kcdf = "Gaussian"), and call all available kernels (parallel.sz = parallel::detectCores());

[0051] S41, Read the data required for MCPcounter evaluation, including cell type gene expression data (MCP_counter_tz.txt) and gene information data (MCP_counter_gene.txt);

[0052] S42, use the MCPcounter.estimate() function to estimate the abundance of cell types in the tumor sample, specifying the gene expression data of cell types (probesets parameter), gene information of cell types (genes parameter), and feature type (featuresType parameter);

[0053] S43 saves the evaluation results in the results object for visualization, allowing you to understand the abundance of different cells in the immune microenvironment.

[0054] We use multiple packages in the R language for data analysis, which makes the process easier.

[0055] Furthermore, in step S3, the method for performing multi-class lasso regression analysis on gene expression data using the glmnet function from the glmnet package is as follows:

[0056] Use the glmnet function to fit a multinomial logistic regression model, where the independent variable x is the gene expression of the sample, the dependent variable y is the classification label of the sample, the family parameter specifies the distribution type of the multinomial as multinomial, and the type.multinomial parameter sets the combination method.

[0057] Use the plot function to visualize the fitted model, setting the x-axis to the coefficient lambda of the regularization penalty term and labeling the value of each coefficient;

[0058] The cv.glmnet function is used to perform cross-validation of the lasso regression model. The optimal lambda value is selected. Based on the optimal lambda value, the glmnet function is used to refit the lasso regression model. The alpha parameter is set to 1 to indicate that lasso regression is used, and the lambda parameter is set to obtain the optimal value through cross-validation. The family parameter and the type.multinomial parameter are the same as before.

[0059] Extract genes with non-zero coefficients from the lasso regression model and save them as gene_min.

[0060] Use the glmnet function from the glmnet package to perform multi-class lasso regression analysis on gene expression data.

[0061] Furthermore, based on gene expression data, DEG analysis of differentially expressed genes between groups was performed using the R package limma (3.56.2) to obtain the differences in gene expression among sample groups; bubble charts were then drawn using the R package GOplot (R 4.0.2) to observe pathway enrichment. The specific steps are as follows:

[0062] Differential gene expression analysis using the limma package:

[0063] Design matrix creation: Use the model.matrix function to create a design matrix, which describes the treatment groups in the experimental design;

[0064] Control matrix creation: Use the makeContrasts function to create a control matrix (contr.matrix), defining the comparison between two control groups (C1VSC2), where C1 and C2 are the two control groups in the design matrix;

[0065] Linear model fitting: Use the lmFit function to fit a linear model (vfit) that fits the expression data (y) and the design matrix (design);

[0066] Contrasts.fit: Use the contrasts.fit function to perform a contrast fitting on the linear model, passing in the contrast matrix (contr.matrix);

[0067] Bayesian estimation: The Bayesian function is used to perform Bayesian estimation on the fitted model to obtain the significance of genes;

[0068] Diagnostic plots were plotted using the plotSA function to assess the model's suitability and stability.

[0069] Significance test: Use the decideTests function to perform a significance test on the model and output a summary of the significance test results for the genes;

[0070] Differential expression gene analysis: The topTable function is used to obtain the differential expression of all genes based on the Bayesian estimation results, and then genes with significant differential expression are selected according to the given thresholds (padj and foldChange);

[0071] Results output: Write the screened differentially expressed genes into a CSV file;

[0072] DEGs functional annotation and enrichment analysis: GO functional annotation and KEGG enrichment analysis were performed on the screened differentially expressed genes (DEGs) using the R package clusterProfiler. Enrichment results were considered significant only when the p-value was less than 0.05. Specifically:

[0073] Gene ID conversion: Use the bitr function to convert the gene symbol from SYMBOL to ENTREZID;

[0074] GO analysis: GO enrichment analysis is performed using the enrichGO function, specifying the gene ID type, p-value, and q-value threshold parameters. The results are saved as a CCGO.rdata file.

[0075] KEGG analysis: KEGG enrichment analysis is performed using the enrichKEGG function. Parameters such as species, p-value, and q-value thresholds are specified, and the results are saved as a CCKEGG.rdata file.

[0076] GO / KEGG enrichment visualization: Visualize the GO and KEGG enrichment results using the barplot and dotplot functions, and save the images as PDF files;

[0077] Gene set comparison analysis: The enrichKEGG function is used to perform KEGG enrichment analysis on genes in different groups, compare the KEGG enrichment results of different groups, and visualize the results.

[0078] Based on gene expression data, DEG analysis of differentially expressed genes between groups was performed using the R package limma (3.56.2) to obtain the differences in gene expression among sample groups. Bubble plots were then generated using the R package GOplot (R 4.0.2) to observe pathway enrichment.

[0079] Furthermore, the method for validating step S4, which involves immune microenvironment typing and biological characteristic screening through single-cell omics analysis, is as follows:

[0080] Single-cell data were collected from samples of patients with medulloblastoma, and Seurat was used for quality control and dimensionality reduction clustering of the single-cell data.

[0081] Based on the collected annotation information, the cell subpopulations of dimensionality reduction clustering are annotated;

[0082] The R package scRNAtoolVis was used to obtain the distribution of specific markers in single-cell subsets and the differences in single-cell levels among medulloblastoma patients with different immune microenvironment subtypes.

[0083] The method for collecting single-cell data from samples of medulloblastoma patients and using Seurat for quality control and dimensionality reduction clustering of single-cell data is as follows:

[0084] For each sample, the number of UMIs, the total number of genes, the number of mitochondrial genes, and the number of ribosomal genes were counted. Cells with a total number of genes greater than 6000 or a number of genes less than 200, as well as cells with a mitochondrial gene ratio greater than 30%, were filtered out.

[0085] For each sample, the top 2000 genes with the greatest variation were identified based on the mean and dispersion of all genes for integrated analysis to eliminate batch effects between samples. Principal component analysis (PCA) was then performed on the integrated data to reduce dimensionality and retain the top 20 principal components to capture the main changes in the data.

[0086] HGV was detected using Seurat's pipeline. The average expression and dispersion of each gene were calculated. Genes were placed into bins, and the z-score for dispersion within each bin was calculated. A z-score of 0.5 was used as the cutoff value for dispersion, and a lower cutoff of 0.0125 and a higher cutoff of 3.0 were used as the average expression. Principal component analysis (PCA) was used for linear dimensionality reduction. Elbow Plot and Jackstraw methods were used to select statistically significant principal components. Clustering was visualized using Seurat based on UMAP. The specific steps are as follows:

[0087] JackStraw: Use the JackStraw method to evaluate the results of PCA and identify significant principal components;

[0088] ScoreJackStraw: Scores the results of the JackStraw method to identify the principal components that are statistically significant;

[0089] JackStrawPlot and ElbowPlot: Visualize the results of the JackStraw method and select the number of principal components to retain. JackStrawPlot is used to observe the p-value of each principal component, while ElbowPlot is used to observe the changes in the variance explained. Typically, the number of principal components at the "elbow" is selected.

[0090] Cell clustering:

[0091] FindNeighbors: Calculates the proximity relationships between cells based on PCA results;

[0092] FindClusters: Based on the K-nearest neighbor graph, a clustering algorithm is used to divide cells into different cell clusters;

[0093] Idents and table: View the clustering results. Use Idents to see the cluster to which each cell belongs, and use table to count the number of cells in each cluster.

[0094] UMAP dimensionality reduction and visualization:

[0095] RunUMAP: Uses the UMAP algorithm to reduce the dimensionality of data. UMAP uses the stochastic gradient descent optimization algorithm, based on the minimum spanning tree and Gaussian mixture model approximation method, and adopts a distance-based weighting function to adjust the similarity weights between different data points, mapping high-dimensional data to two-dimensional or three-dimensional space.

[0096] Minimize the loss function between the distances between data points in high-dimensional space and the distances between corresponding points in low-dimensional space to optimize the data representation in low-dimensional space;

[0097] DimPlot: Based on UMAP dimensionality reduction, it plots the distribution of cells, visually demonstrating the relationships and distribution of different cell clusters;

[0098] UMAP clustering visualization results showed that individual cells from 26 samples from SHH, GP3, and GP4 were clustered into 21 (SHH), 18 (GP3), and 15 (GP4) subgroups, respectively.

[0099] Single-cell data were collected from samples of patients with medulloblastoma, and Seurat was used for quality control and dimensionality reduction clustering of the single-cell data for ease of use.

[0100] Furthermore, based on the collected annotation information, the method for annotating the cell subpopulations of dimensionality reduction clustering is as follows:

[0101] scRNA sequencing was performed on MB patient samples and single matched samples from relapsed patients. Quality-controlled cells were projected into a two-dimensional UMAP map, and cell subpopulations in this experiment were manually annotated based on marker genes from the Cell Marker database and other brain tumor-related literature.

[0102] Based on the collected annotation information, the cell subpopulations of dimensionality reduction clustering are annotated.

[0103] The present invention also provides a system for detecting and classifying differences in the immune microenvironment of medulloblastoma, comprising a processing unit, wherein the processing unit performs the method described in the present invention to detect and classify differences in the immune microenvironment of medulloblastoma.

[0104] This system, through its processing unit, acquires the distribution of specific markers in single-cell subsets and the differences in single-cell levels among medulloblastoma patients with different immune microenvironment subtypes, providing new ideas and methods for disease treatment and prognostic assessment. Attached Figure Description

[0105] Figure 1 This is a schematic flowchart of the method for detecting and classifying immune microenvironment differences in medulloblastoma according to the present invention. Detailed Implementation

[0106] Embodiments of the present invention are described in detail below. Examples of these embodiments are shown in the accompanying drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to the accompanying drawings are exemplary and are only used to explain the present invention, and should not be construed as limiting the present invention.

[0107] In the description of this invention, it should be understood that the terms "longitudinal", "lateral", "up", "down", "front", "rear", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "outer", etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings. They are only for the convenience of describing this invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limitations on this invention.

[0108] In the description of this invention, unless otherwise specified and limited, it should be noted that the terms "installation", "connection" and "linking" should be interpreted broadly. For example, they can refer to mechanical or electrical connections, or internal connections between two components. They can be direct connections or indirect connections through an intermediate medium. Those skilled in the art can understand the specific meaning of the above terms according to the specific circumstances.

[0109] This invention discloses a method for detecting and classifying the immune microenvironment differences in medulloblastoma, exploring the immune microenvironment of medulloblastoma and classifying it, aiming to provide strong support for the formulation of personalized treatment strategies. This method can be used for early diagnosis, prognosis prediction, development of specific therapies, and monitoring of treatment effects. Figure 1 As shown, the method for detecting and classifying differences in the immune microenvironment of medulloblastoma includes the following steps:

[0110] S1: Collect gene expression data from samples of patients with medulloblastoma and perform pseudo-single-cell calculations.

[0111] S2 uses Netclass (NetClass is a network-based classification method that transforms network structure information into feature vectors and then uses a traditional machine learning classifier for prediction. Its main goal is to improve classification performance by utilizing the topological properties of the network. The core idea of ​​the NetClass algorithm is to extract structural information from network data into feature vectors and then use these feature vectors for classification) to generate a protein-protein interaction (PPI) network as an adjacency matrix. Gene expression data is introduced into the protein-protein interaction network, and then combined with pseudo-single-cell gene expression data. Unsupervised clustering is performed through similarity network fusion analysis to determine different immune microenvironment subtypes in medulloblastoma.

[0112] S3 uses lasso for machine learning to screen biological features of different immune microenvironment subtypes;

[0113] S4 validates the immune microenvironment typing and biological feature screening through single-cell omics analysis.

[0114] In a preferred embodiment of the present invention, step S1 specifically comprises:

[0115] Based on the ssGSEA (Single-sample Gene Set Enrichment Analysis) pseudo-single-cell analysis method, gene expression data are compared with a given gene set to calculate the enrichment score of the gene set in each sample, simulating the composition of immune cells in the tumor immune microenvironment. Specifically:

[0116] The gene expression data for each sample are sorted from high to low (sorting the gene expression data for each sample is usually done by sorting the gene expression level in the sample from high to low. This sorting reflects the expression intensity of the gene in the sample), and the cumulative distribution function (CDF) of each gene set is calculated. CDF is an important concept in probability theory and statistics. For a random variable X, its cumulative distribution function F(x) is defined as the probability that the random variable X is less than or equal to x, i.e., F(x) = P(X≤x)). This includes the gene set CDF (calculated based on the position of genes in the sorted list) and the background gene CDF (the cumulative distribution calculated based on the expression levels of all genes).

[0117] Enrichment scores are calculated by comparing the CDF of the gene set with that of the background genes, where N is the total number of genes and CDF represents the cumulative distribution function.

[0118] Obtain expression data for the gene set and background genes: This can be gene chip data, RNA-seq data, or the results of other gene expression measurements, ensuring that the expression values ​​of the gene set and background genes are ready.

[0119] Calculate the CDF of the gene set and background genes: Sort the gene set and background genes in ascending order of expression value. Then, calculate the cumulative distribution value of each gene, which is the number of genes with an expression value less than or equal to that gene divided by the total number of genes. This will give you the CDF curves for the gene set and background genes.

[0120] Comparing the CDF of the gene set with that of the background gene set: Visualization methods, such as plotting CDF curves, can be used to visually compare the distribution of the gene set and the background genes. Observe whether the gene set CDF is above or below the background gene CDF, and the degree of difference between them.

[0121] Calculating the enrichment fraction: Calculate the area difference between the CDF of the gene set and the CDF of the background gene set. The specific calculation method can be selected according to the specific situation. A simple method is to calculate the integral area between the two CDF curves as the enrichment fraction. Numerical integration or other appropriate methods can be used to calculate the area.

[0122] The enrichment scores are standardized and compared across different gene sets, assigning them to the most likely cells in order to identify the cell type most likely to be represented by each cell at the single-cell level; here, "most likely" means that the cell has the highest enrichment score on a certain gene set, and is therefore most likely to represent the cell type associated with that gene set.

[0123] Based on the netClass analysis method, gene expression data and protein network topology are integrated using information from protein-protein interaction (PPI) networks.

[0124] The pseudo-single-cell results obtained by the ssGSEA method were combined with gene expression data to perform similarity network fusion (SNF) analysis on the samples.

[0125] This study utilizes the following packages in R: `survival` (primarily used for survival analysis, a statistical method for processing time-of-occurrence data, commonly used in medical research (e.g., patient survival time), engineering reliability studies (e.g., equipment failure time), `ggplot2` (a powerful tool for creating high-quality data visualizations based on a "graphical syntax" that allows users to construct complex graphs hierarchically), and `ESTIMATE` (primarily used to assess stromal and immune cell infiltration levels in tumor tissue; based on gene expression data, it provides researchers with information on tumor purity, the level of stromal cells present, and the level of immune cells in the tumor tissue). The following packages were used to visualize and analyze survival data of different immune subtypes of medulloblastoma (MB) patients: the score of immune cell infiltration level, the score of MCPcounter (which can quantify the absolute abundance of 8 immune cells and 2 stromal cells in heterogeneous tissues using transcriptome data), and the score of msigdbr (which is mainly used to assess the level of stromal and immune cell infiltration in tumor tissues. It provides researchers with scores of tumor purity, the level of stromal cells present, and the level of immune cell infiltration in tumor tissues based on gene expression data).

[0126] Use the glmnet function in the glmnet package (the "glmnet" package is a language in R for fitting generalized linear models) to perform multi-class lasso regression analysis on gene expression data;

[0127] Based on gene expression data, the R package limma ("limma" is a powerful and widely used package in R language for gene expression data analysis) (3.56.2) was used to perform DEG analysis of differentially expressed genes between groups to obtain the differences in data expression between sample groups; the R package GOplot (R4.0.2) was used to draw bubble charts to observe the pathway enrichment.

[0128] In a preferred embodiment of the present invention, gene expression data and protein network topology are integrated using the netClass analysis method and information from protein-protein interaction networks. The specific steps are as follows:

[0129] Construction of protein interaction networks: Obtain protein interaction network data from public databases or literature, and construct a network containing the interaction relationships between proteins as an adjacency matrix;

[0130] Data processing, including data filtering and dimension adjustment, is performed to ensure the accuracy and effectiveness of subsequent calculations.

[0131] Calculate the diffusion kernel: This function accepts an adjacency matrix as input and calculates the diffusion kernel matrix based on specified parameters (p and a). The kernel matrix is ​​generated iteratively. During kernel calculation, the parameter `is.adjacency` (a user-defined parameter indicating whether a matrix is ​​an adjacency matrix) determines whether the input is an adjacency matrix or a Laplace matrix, selecting different computation paths accordingly. (If it's an adjacency matrix, calculations are performed based on the node connections represented by the adjacency matrix; if it's a Laplace matrix, its unique properties (such as the relationship between the degree matrix and the adjacency matrix) are utilized for Laplace matrix calculation.)

[0132] Use the igraph library to perform graph operations, including creating undirected graphs and calculating the Laplacian matrix;

[0133] When calculating the diffusion kernel matrix and subsequent eigenvector multiplication, matrix operations are used to filter out the intersection with another dataset, and then the diffusion kernel is calculated using a portion of the data from the intersection.

[0134] In a preferred embodiment of the present invention, the method for performing similar network fusion (SNF) analysis on samples by combining pseudo-single-cell results obtained through the ssGSEA method with gene expression data is as follows:

[0135] The SNF method does not require any prior feature selection, therefore it uses the complete expression matrix after removing pseudogenes and the ssGSEA matrix as input;

[0136] Using the SNFtool R package (v2.2.0), with the number of neighbors K=30, Gaussian kernel parameter alpha=0.5, and number of iterations T=10, the spectral clustering implemented by the SNF tool package is run on the SNF fused similarity matrix to obtain the most suitable grouping;

[0137] SNF analysis was performed in different subtypes: In the SHH (Sonic Hedgehog) subtype (sample size n=233, the optimal grouping recommended by SNF, this parameter is the cluster number parameter in the spectral clustering algorithm. During spectral clustering, the algorithm divides the data into K clusters), 3 subclusters of immune microenvironments were obtained when k=3;

[0138] In the GP3 (n=144) subtype, two immune microenvironment subpopulations were obtained when k=2;

[0139] In the GP4 (n=326) subtype, three immune microenvironment subgroups were obtained when k=3.

[0140] K is the number of neighbors for each node when constructing a similarity network. In SNF (Similarity Network Fusion), each sample node is connected to its K most similar sample nodes to build a K-nearest neighbors graph. Setting K=30 means that each sample is connected to its 30 most similar samples. Choosing an appropriate value for K is important because it affects the sparsity of the network and the final clustering results.

[0141] alpha is a parameter that controls the width of the Gaussian kernel. When calculating the similarity matrix, the Gaussian kernel function is used to measure the similarity between samples. Specifically, the elements of the similarity matrix are calculated using the Gaussian kernel function. Setting alpha = 0.5 means that 0.5 is used as the parameter for the Gaussian kernel function when calculating similarity. A larger alpha value will make the similarity values ​​more concentrated, while a smaller alpha value will make the similarity values ​​more dispersed.

[0142] T represents the number of iterations in the SNF algorithm. In SNF, information from different data sources is fused by iteratively updating the similarity matrix. In each iteration, the current similarity matrix is ​​fused and updated with similarity matrices from other data sources. Setting T=10 indicates 10 iterations to update and fuse the similarity matrix. Generally, more iterations result in a more stable and accurate fusion result, but also increase computation time.

[0143] In a preferred embodiment of the present invention, the survival, ggplot2, ESTIMATE, msigdbr, and MCPcounter packages in R language are used to visualize and analyze the survival data of MB patients, assess the tumor microenvironment by evaluating gene expression data, obtain gene set data related to biological processes from the species' genomics database, and assess the abundance of different cell types in the tumor microenvironment. The steps are as follows:

[0144] S11, Prepare the data required for survival analysis, including the observation time (futime) and event status (fustat) for each sample, as well as the group variable for grouping;

[0145] S12, use the survfit() function to fit the survival curve, calculate the survival curve based on the specified observation time, event state and grouping variable group;

[0146] S13. Use the ggsurvplot() function to create a survival curve chart. In the chart, set parameters such as not displaying confidence intervals (conf.int=F), displaying risk tables (risk.table=T), not adding overall patient survival curves (add.all=F), and customizing the color palette (palette="Dark2"). Adjust the title, axis titles, font size and style, and the range of the axes to help understand the survival situation between different groups and provide a basis for subsequent data interpretation and clinical research.

[0147] S21. Prepare the gene expression data file. Use the outputGCT() function to convert the data into the input format required by the ESTIMATE package and save it as a new file (in.gct.file).

[0148] S22, use the filterCommonGenes() function to filter out common genes from the original data and save them as a new input file, use the estimateScore() function to calculate the ESTIMATE score of the tumor sample, and save the score results as a GCT format file (out.score.file);

[0149] S23, then use the plotPurity() function to visualize the purity score of the sample, use the read.table() function to read the score result file and convert it into a data frame format, and finally save the result as a text file (ESTIMATE_score.txt) for subsequent analysis and visualization;

[0150] S31, use the msigdbr() function to retrieve the gene set data of C5 category from the genomics database of the species and save it in the GO_df_all data frame;

[0151] The dplyr::select() function was used to select the desired columns (gs_name, gene_symbol, gs_exact_source, and gs_subcat), and unwanted data was filtered out according to the subcategory (gs_subcat) to obtain the final GO gene set data (GO_df);

[0152] S32, use the split() function to group the gene set data according to gs_name (GO entry name) to obtain a gene list (go_list);

[0153] Gene set variability analysis (GSVA) ​​is performed on the gene expression data (rt1) to be analyzed using the gsva() function. The gene set variability score (gsva_mat) for each GO entry is calculated and the results are saved in a text file (gsva_hallmark.txt).

[0154] S33, set parameters: use the Gaussian kernel density estimation function (kcdf = "Gaussian") and call all available kernels (parallel.sz = parallel::detectCores()); these analysis steps can help understand the relative expression levels of gene sets in a sample and reveal the activity levels of different biological processes in the sample.

[0155] S41, Read the data required for MCPcounter (Microenvironment Cell Populations-counter) assessment, including gene expression data of cell types (MCP_counter_tz.txt) and gene information data (MCP_counter_gene.txt);

[0156] S42, use the MCPcounter.estimate() function to estimate the abundance of cell types in the tumor sample, specifying the gene expression data of cell types (probesets parameter), gene information of cell types (genes parameter), and feature type (featuresType parameter);

[0157] S43 saves the evaluation results in the results object for visualization, allowing you to understand the abundance of different cells in the immune microenvironment.

[0158] In a preferred embodiment of the present invention, step S3, using the glmnet function in the glmnet package to perform multi-class lasso regression analysis (Least Absolute Shrinkage and Selection Operator Regression, a regularization method for linear regression) on gene expression data, is as follows:

[0159] Use the `glmnet` function to fit a multinomial logistic regression model, where the independent variable `x` is the gene expression of the samples and the dependent variable `y` is the classification label of the samples. Use the `family` parameter (in the `glmnet` package of R, the `family` parameter specifies the type of regression model. For example, for Lasso regression, different distribution types can be selected to fit the corresponding model by setting the `family` parameter) to specify the multinomial distribution type as `multinomial` (when using the `glmnet` package for regression analysis, `multinomial` is usually used to specify the type of regression model as multinomial logistic regression). Use the `type.multinomial` parameter (the `type.multinomial` parameter specifies the type of treatment in multi-class problems. In the `glmnet` package, it has two optional values: "ungrouped" (default) and "grouped") to set the combination method.

[0160] Use the plot function to visualize the fitted model, setting the x-axis to the coefficient lambda of the regularization penalty term and labeling the value of each coefficient;

[0161] The cv.glmnet function is used to perform cross-validation of the lasso regression model. The optimal lambda value is selected. Based on the optimal lambda value, the glmnet function is used to refit the lasso regression model. The alpha parameter is set to 1 to indicate that lasso regression is used, and the lambda parameter is set to obtain the optimal value through cross-validation. The family parameter and the type.multinomial parameter are the same as before.

[0162] Genes with non-zero coefficients in the lasso regression model are extracted and saved as gene_min, thus identifying feature genes that have a significant impact on classification results in multi-classification problems.

[0163] In a preferred embodiment of the present invention, based on gene expression data, differentially expressed genes (DEG) analysis between groups is performed using the R software package limma (3.56.2) to obtain the differences in gene expression among sample groups; bubble charts are drawn using the R software package GOplot (R4.0.2) to observe pathway enrichment. The specific steps are as follows:

[0164] Differential gene expression analysis using the limma package:

[0165] Design matrix creation: Use the model.matrix function to create a design matrix, which describes the treatment groups in the experimental design;

[0166] Control matrix creation: Use the makeContrasts function to create a control matrix (contr.matrix), defining the comparison between two control groups (C1VSC2), where C1 and C2 are the two control groups in the design matrix;

[0167] Linear model fitting: Use the lmFit function to fit a linear model (vfit) that fits the expression data (y) and the design matrix (design);

[0168] Contrasts.fit: Use the contrasts.fit function to perform a contrast fitting on the linear model, passing in the contrast matrix (contr.matrix);

[0169] Bayesian estimation: The Bayesian function is used to perform Bayesian estimation on the fitted model to obtain the significance of genes;

[0170] Diagnostic plots were plotted using the plotSA function to assess the model's suitability and stability.

[0171] Significance test: Use the decideTests function to perform a significance test on the model and output a summary of the significance test results for the genes;

[0172] Differential expression gene analysis: The topTable function is used to obtain the differential expression of all genes based on the Bayesian estimation results, and then genes with significant differential expression are selected according to the given thresholds (padj and foldChange);

[0173] Results output: Write the screened differentially expressed genes into a CSV file;

[0174] DEGs (Different express genes) functional annotation and enrichment analysis: The R package clusterProfiler was used to perform GO functional annotation and KEGG enrichment analysis on the screened differentially expressed genes (DEGs). Enrichment results were considered significant only when the p-value was less than 0.05. Specifically:

[0175] Gene ID conversion: Use the bitr function to convert the gene symbol from SYMBOL to ENTREZID;

[0176] GO (Gene Ontology) analysis: GO enrichment analysis is performed using the enrichGO function, specifying the gene ID type, p-value, and q-value threshold parameters. The results are saved as a CCGO.rdata file.

[0177] KEGG (Kyoto Encyclopedia of Genes and Genomes) analysis: KEGG enrichment analysis was performed using the enrichKEGG function, specifying parameters such as species, p-value, and q-value thresholds, and the results were saved as a CCKEGG.rdata file;

[0178] GO / KEGG enrichment visualization: Visualize the GO and KEGG enrichment results using the barplot and dotplot functions, and save the images as PDF files;

[0179] Gene set comparison analysis: The enrichKEGG function is used to perform KEGG enrichment analysis on genes in different groups, compare the KEGG enrichment results of different groups, and visualize the results.

[0180] For example, using SNF, the tumor microenvironment of MB samples of the SHH subtype was further divided into three subtypes: C1, C2, and C3. The C2 subtype showed significant enrichment in activated CD8 T cells, γδ T cells, activated CD4 T cells, and type II helper T cells, with expression levels higher than the other two subtypes. Most immune cells, such as macrophages, natural killer cells, central memory CD4 T cells, type I helper T cells, activated dendritic cells, immunosuppressive myeloid cells, central memory CD8 T cells, and regulatory T cells, were enriched in both the C2 and C3 subtypes, but the expression in the C3 subtype was stronger than that in the C2 subtype. The C1 subtype showed almost no immune cells, with only slight expression in a small number of immune cells such as immature dendritic cells, CD56 bright-type natural killer cells, CD56 dark-type natural killer cells, and plasmacytoid dendritic cells, exhibiting an overall immune cold state.

[0181] Survival analysis revealed significant differences in survival rates among different subtypes, with subtype C2 exhibiting the worst survival and subtype C1 showing the highest overall survival. Subtype C1 patients were significantly older, predominantly adults; subtype C3 patients were predominantly infants and young children; and subtype C2 patients were predominantly children, but there was no significant difference in gender distribution.

[0182] The immune microenvironment of GP3 MB was divided into two subtypes, C1 and C2, using SNF. The C2 subtype showed significantly higher enrichment in γδ T cells, type II helper T cells, activated CD8 T cells, activated CD4 T cells, and memory B cells. Both subtypes expressed CD56 bright-type natural killer cells, effector memory CD4 T cells, plasmacytoid dendritic cells, immunosuppressive myeloid cells (immunosuppressive myeloid-derived cells), neutrophils, central memory CD4 T cells, natural killer cells, eosinophils, and macrophages, but the C2 subtype showed slightly stronger expression than the C1 subtype. Both subtypes also expressed activated B cells, monocytes, type 17 helper T cells, immature B cells, T follicular helper cells, mast cells, and type I helper T cells, but the C1 subtype showed stronger expression than the C2 subtype.

[0183] In the survival analysis, the C2 subtype had a significantly worse survival rate than the C1 subtype within 100 months.

[0184] The immune microenvironment of GP4 MB was divided into three subtypes: C1, C2, and C3 using SNF. Mast cells, type 17 T helper cells, activated B cells, monocytes, CD56 bright-type natural killer cells, CD56 dark-type natural killer cells, plasmacytoid dendritic cells, and effector memory CD4 T cells were enriched in subtypes C1 and C3 with similar enrichment levels, but were almost not expressed in subtype C2. Memory B cells, activated CD8 T cells, activated CD4 T cells, type II T helper cells, macrophages, T follicular helper cells, eosinophils, and immature dendritic cells were expressed in all three subtypes, but the expression was weakest in subtype C3, followed by subtype C2, and strongest in subtype C1. The remaining immune cell types, including natural killer T cells, type I T helper cells, immature B cells, neutrophil populations, activated dendritic cells, and effector memory CD8 T cells, were most strongly expressed in subtype C1, and almost not expressed in the other two subtypes. In the survival analysis, the C1 subtype had the worst survival rate over 100 months.

[0185] In a preferred embodiment of the present invention, the method for verifying the immune microenvironment typing and biological characteristic screening through single-cell omics analysis in step S4 is as follows:

[0186] Single-cell data were collected from samples of patients with medulloblastoma, and Seurat (a single-cell workflow analysis method) was used for quality control and dimensionality reduction clustering of the single-cell data.

[0187] Based on the collected annotation information, the cell subpopulations of dimensionality reduction clustering are annotated;

[0188] Using the R package scRNAtoolVis (“scRNAtoolVis” is an R package for single-cell RNA sequencing (scRNA-seq) data analysis and visualization. An R package is a collection of functions, datasets, and documents in the R language environment used to extend its functionality), we can obtain the distribution of specific markers in single-cell subsets and the differences in single-cell levels among medulloblastoma patients with different immune microenvironment subtypes.

[0189] The method for collecting single-cell data from samples of medulloblastoma patients and using Seurat for quality control and dimensionality reduction clustering of single-cell data is as follows:

[0190] For each sample, the number of UMIs, the total number of genes, the number of mitochondrial genes, and the number of ribosomal genes were counted. Cells with a total number of genes greater than 6000 or a number of genes less than 200, as well as cells with a mitochondrial gene ratio greater than 30%, were filtered out.

[0191] For each sample, the top 2000 genes with the greatest variation were identified based on the mean and dispersion of all genes for integrated analysis to eliminate batch effects between samples. Principal component analysis (PCA) was then performed on the integrated data to reduce dimensionality and retain the top 20 principal components to capture the main changes in the data.

[0192] HGV was detected using Seurat's pipeline. The average expression and dispersion of each gene were calculated. Genes were placed into bins, and the z-score of dispersion within each bin was calculated. A z-score of 0.5 was used as the cutoff value for dispersion, and a lower cutoff of 0.0125 and a higher cutoff of 3.0 were used as the average expression. Principal component analysis (PCA) was used for linear dimensionality reduction. Elbow Plot and Jackstraw methods were used to select statistically significant principal components. Clustering was visualized using Seurat based on UMAP (Uniform Manifold Approximation and Projection). The specific steps are as follows:

[0193] JackStraw: Use the JackStraw method to evaluate the results of PCA and identify significant principal components;

[0194] ScoreJackStraw: Scores the results of the JackStraw method to identify the principal components that are statistically significant;

[0195] JackStrawPlot and ElbowPlot: Visualize the results of the JackStraw method and select the number of principal components to retain. JackStrawPlot is used to observe the p-value of each principal component, while ElbowPlot is used to observe the changes in the variance explained. Typically, the number of principal components at the "elbow" is selected.

[0196] Cell clustering:

[0197] FindNeighbors: Calculates the proximity relationships between cells based on PCA results;

[0198] FindClusters: Based on the K-nearest neighbor graph, clustering algorithms (such as the Louvain algorithm) are used to divide cells into different cell clusters.

[0199] Idents and table: View the clustering results. Use Idents to see the cluster to which each cell belongs, and use table to count the number of cells in each cluster.

[0200] UMAP dimensionality reduction and visualization:

[0201] RunUMAP: Uses the UMAP algorithm to reduce the dimensionality of data. UMAP uses the stochastic gradient descent optimization algorithm, based on the minimum spanning tree and Gaussian mixture model approximation method, and adopts a distance-based weighting function to adjust the similarity weights between different data points, mapping high-dimensional data to two-dimensional or three-dimensional space.

[0202] Minimize the loss function between the distances between data points in high-dimensional space and the distances between corresponding points in low-dimensional space to optimize the data representation in low-dimensional space;

[0203] DimPlot: Based on UMAP dimensionality reduction, it plots the distribution of cells, visually demonstrating the relationships and distribution of different cell clusters;

[0204] UMAP clustering visualization results showed that individual cells from 26 samples from SHH, GP3, and GP4 were clustered into 21 (SHH), 18 (GP3), and 15 (GP4) subgroups, respectively.

[0205] In a preferred embodiment of the present invention, the method for annotating cell subpopulations of dimensionality reduction clustering by combining the collected annotation information is as follows:

[0206] scRNA sequencing was performed on MB patient samples and single matched samples from relapsed patients. Quality-controlled cells were projected into a two-dimensional UMAP map, and cell subpopulations in this experiment were manually annotated based on marker genes from the Cell Marker database and other brain tumor-related literature.

[0207] The first part of this invention involves pseudo-single-celling the gene expression data (SHH=223, GP3=144, GP4=326) of MB patients, generating a PPI network as an adjacency matrix using Netclass, and introducing the gene expression matrix (16116 genes) after removing pseudogenes to perform similarity network fusion (SNF) analysis.

[0208] (1) ssGSEA analysis was performed on MB data of different subtypes. Based on the pseudo-single cell results obtained, SNF analysis was performed in combination with gene expression data. Based on the differences in the tumor immune microenvironment, different immune microenvironment subtypes in medulloblastoma of the three subtypes SHH, GP3 and GP4 were identified.

[0209] (2) Combine other information, such as histopathological information and clinical information, to conduct bioinformatics analysis and study and discuss the biological characteristics of each subtype;

[0210] (3) Perform various immune infiltration analyses to verify the stability of SNF typing;

[0211] (4) Marker genes with significant survival differences in each subtype were screened, and functional differences between different subtypes were observed through gene enrichment analysis. MB were clearly divided into different subtypes with different biological and clinical characteristics, and genes with prognostic differences in different subtypes were screened. Each subtype has its own unique biological marker. The C2 subtype of SHH MB, which has the worst prognosis, specifically expresses PAWR; in GP3MB, the C1 subtype, which has the best prognosis, specifically expresses ZMYND8; and the GP4 MB of the C1 subtype, which has the worst prognosis, specifically expresses PLAT.

[0212] The second part of this invention uses single-cell data from medulloblastoma patients (SHH subtype = 9, GP3 subtype = 6, GP4 subtype = 12) for analysis:

[0213] (1) Use Seurat for quality control and dimensionality reduction clustering of single-cell data;

[0214] (2) Annotate the cell subpopulations of dimensionality reduction clustering based on the collected annotation information;

[0215] (3) Observe the distribution of specific maker in single-cell subsets and the differences in single-cell levels among MB patients with different immune microenvironment subtypes.

[0216] At the single-cell level, in SHH MB, compared to C1 and C3, C2 has significantly fewer myeloid cells and a large number of unipolar brush cells (UBCs); C1 has a large number of granule cell precursors (GCPs), which are almost absent in C2 and C3. In GP3 MB, the main difference between C1 and C2 remains UBCs. C2 contains a large number of UBCs, while C1 has almost none. In GP4 MB, the situation is completely reversed. The C3 subtype, with the best prognosis, expresses a large number of UBCs and has the fewest myeloid cells, while the C1 subtype, with the worst prognosis, is the opposite, highly expressing transitional cerebellar progenitors (TCPs). Low expression of myeloid cells is always accompanied by high expression of UBCs, suggesting that UBCs and myeloid cells may play a key role in promoting tumor development and progression, and that there may be some kind of relationship between them.

[0217] The three subtypes of medulloblastoma each possess distinct immune microenvironment subtypes. These subtypes have unique markers and functions, playing unique roles in tumor development and progression, and leading to completely different outcomes for patients. Through immune microenvironment typing and marker screening, the immune characteristics and potential therapeutic targets of medulloblastoma have been revealed, providing new ideas and methods for the treatment and prognostic assessment of this disease.

[0218] For example, by drawing violin diagrams of differentially expressed genes in different samples and cells, their distribution at the single-cell level can be observed. FBXO16 and HCAR1, which are highly expressed in the C1 subtype, are more abundant in oligodendrocytes and GCP cells, and less abundant in myeloid and endothelial cells; ZNF467 and MTHFD1L, which are highly expressed in the C3 subtype, are more abundant in myeloid cells; and PAWR, which is specifically expressed in the C2 subtype, is highly enriched in UBC cells.

[0219] The present invention also provides a system for detecting and classifying differences in the immune microenvironment of medulloblastoma, comprising a processing unit that executes the method described in the present invention to detect and classify differences in the immune microenvironment of medulloblastoma.

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

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

Claims

1. A method for detecting and classifying immune microenvironment differences in medulloblastoma, characterized in that, The steps include the following: S1: Collect gene expression data from samples of patients with medulloblastoma and perform pseudo-single-cell calculations. S2, using Netclass to generate a protein interaction network as an adjacency matrix, introduces gene expression data into the protein interaction network, and then combines pseudo-single-cell gene expression data, and uses similarity network fusion analysis to perform unsupervised clustering to determine different immune microenvironment subtypes in medulloblastoma; S3 uses lasso for machine learning to screen biological features of different immune microenvironment subtypes; S4, validated by single-cell omics analysis for immune microenvironment typing and biological feature screening; Step S1 is as follows: Based on the ssGSEA pseudo-single-cell analysis method, gene expression data are compared with a given gene set to calculate the enrichment score of the gene set in each sample, simulating the composition of immune cells in the tumor immune microenvironment. Based on the Netclass analysis method, gene expression data and protein network topology are integrated using information from protein-protein interaction (PPI) networks. The pseudo-single-cell results obtained by the ssGSEA method were combined with gene expression data to perform similar network fusion (SNF) analysis on the samples. The survival, ggplot2, ESTIMATE, MCPcounter, and msigdbr packages in R were used to visualize and analyze the survival data of different immune subtypes of MB patients, evaluate the abundance of the tumor microenvironment and different cell types in the tumor microenvironment, and evaluate the gene set data related to biological processes. Use the glmnet function in the glmnet package to perform multi-class lasso regression analysis on gene expression data; Based on gene expression data, DEG analysis of differentially expressed genes between groups was performed using the R package limma to obtain the differences in gene expression among sample groups; bubble charts were drawn using the R package GOplot to observe pathway enrichment.

2. The method for detecting and classifying immune microenvironment differences in medulloblastoma as described in claim 1, characterized in that, Based on the Netclass analysis method, gene expression data and protein network topology are integrated using information from protein-protein interaction networks. The specific steps are as follows: Data on protein-protein interaction networks are obtained from public databases or literature, and a network containing the interaction relationships between proteins is constructed as an adjacency matrix. Calculate the diffusion kernel: Accepts the adjacency matrix as input and calculates the diffusion kernel matrix according to the specified parameters p and a. The diffusion kernel matrix is ​​generated through iterative calculation. When calculating the diffusion kernel matrix, the parameter is.adjacency is used to determine whether the input is an adjacency matrix or a Laplacian matrix, and different calculation paths are selected accordingly. Use the igraph library to perform graph operations, including creating undirected graphs and calculating the Laplacian matrix; When calculating the diffusion kernel matrix and subsequent eigenvector multiplication, matrix operations are used to filter out the intersection with another dataset, and then the diffusion kernel is calculated using a portion of the data from the intersection.

3. The method for detecting and classifying immune microenvironment differences in medulloblastoma as described in claim 1, characterized in that, The method for performing SNF analysis on pseudo-single-cell results obtained through the ssGSEA method, combined with gene expression data, is as follows: Use the complete expression matrix after removing pseudogenes and the ssGSEA matrix as input; Using the SNFtool R package, with the number of neighbors K = 30, Gaussian kernel parameter alpha = 0.5, and number of iterations T = 10, the spectral clustering implemented by the SNF tool package is run on the SNF fused similarity matrix to obtain the most suitable grouping; SNF analysis was performed in different subtypes: In the SHH subtype, with a sample size of n=233, three immune microenvironment subpopulations were obtained when k=3; In the GP3 subtype, with n=144, two immune microenvironment subpopulations were obtained when k=2; In the GP4 subtype, n=326, and three subpopulations of the immune microenvironment were obtained when k=3.

4. The method for detecting and classifying immune microenvironment differences in medulloblastoma as described in claim 1, characterized in that, Using the `survival`, `ggplot2`, `ESTIMATE`, `msigdbr`, and `MCPcounter` packages in R, we performed visualization analysis on survival data of MB patients, assessed the tumor microenvironment using gene expression data, obtained gene set data related to biological processes from the species' genomics database, and evaluated the abundance of different cell types in the tumor microenvironment. The steps are as follows: S11, Prepare the data required for survival analysis, including the observation time (futime) and event status (fustat) for each sample, as well as the group variable for grouping; S12, use the survfit() function to fit the survival curve, calculate the survival curve based on the specified observation time, event state and grouping variable group; S13, use the ggsurvplot() function to create a survival curve chart, set the following settings in the chart: do not display confidence intervals conf.int = F, display risk table risk.table = T, do not add total patient survival curves add.all = F, customize the color palette palette = "Dark2", and adjust the title, axis titles, font size and style, as well as the range of the axes; S21. Prepare the gene expression data file. Use the outputGCT() function to convert the data into the input format required by the ESTIMATE package and save it as a new file in.gct.file. S22, use the filterCommonGenes() function to filter out common genes from the original data and save them as a new input file. Use the estimateScore() function to calculate the ESTIMATE score of the tumor samples and save the score results as a GCT format file out.score.file; S23, then use the plotPurity() function to visualize the purity score of the samples. Use the read.table() function to read the score result file and convert it to a data frame format. Finally, save the results as a text file ESTIMATE_score.txt; S31, use the msigdbr() function to retrieve the gene set data of C5 category from the genomics database of the species and save it in the GO_df_all data frame; The dplyr::select() function was used to select the required columns gs_name, gene_symbol, gs_exact_source, and gs_subcat, and the unwanted data was filtered according to the subcategory gs_subcat to obtain the final GO gene set data GO_df; S32, use the split() function to group the gene set data according to the GO entry name gs_name to obtain a gene list go_list; The gene expression data rt1 to be analyzed is subjected to gene set variability analysis (GSVA) ​​using the gsva() function. The gene set variability score gsva_mat for each GO entry is calculated and the results are saved in the text file gsva_hallmark.txt. S33, set parameters: use the Gaussian kernel density estimation function kcdf="Gaussian" and call all available kernels parallel.sz = parallel::detectCores(); S41, Read the data required for MCPcounter evaluation, including cell type gene expression data MCP_counter_tz.txt and gene information data MCP_counter_gene.txt; S42, use the MCPcounter.estimate() function to estimate the abundance of cell types in the tumor sample, specifying the probesets parameter for cell type gene expression data, the genes parameter for cell type gene information, and the featuresType parameter; S43 saves the evaluation results in the results object for visualization, allowing you to understand the abundance of different cells in the immune microenvironment.

5. The method for detecting and classifying immune microenvironment differences in medulloblastoma as described in claim 1, characterized in that, In step S3, the method for performing multi-class lasso regression analysis on gene expression data using the glmnet function from the glmnet package is as follows: Use the glmnet function to fit a multinomial logistic regression model, where the independent variable x is the gene expression of the sample, the dependent variable y is the classification label of the sample, the family parameter specifies the distribution type of the multinomial as multinomial, and the type.multinomial parameter sets the combination method. Use the plot function to visualize the fitted model, setting the x-axis to the coefficient lambda of the regularization penalty term and labeling the value of each coefficient; The cv.glmnet function is used to perform cross-validation of the lasso regression model. The optimal lambda value is selected. Based on the optimal lambda value, the glmnet function is used to refit the lasso regression model. The alpha parameter is set to 1 to indicate that lasso regression is used, and the lambda parameter is set to obtain the optimal value through cross-validation. The family parameter and the type.multinomial parameter are the same as before. Extract genes with non-zero coefficients from the lasso regression model and save them as gene_min.

6. The method for detecting and classifying immune microenvironment differences in medulloblastoma as described in claim 1, characterized in that, Based on gene expression data, DEG analysis of differentially expressed genes between groups was performed using the R package limma to obtain the differences in gene expression among sample groups; bubble charts were then drawn using the R package GOplot to observe pathway enrichment. The specific steps are as follows: Differential gene expression analysis using the limma package: Design matrix creation: Use the model.matrix function to create a design matrix, which describes the treatment groups in the experimental design; Control matrix creation: Use the makeContrasts function to create the control matrix contr.matrix, defining the comparison between two control groups C1VSC2, where C1 and C2 are the two control groups in the design matrix; Linear model fitting: Use the lmFit function to fit a linear model vfit, which is used to fit the expression data y and the design matrix design; Contrasts.fit: Use the contrasts.fit function to perform a contrast fitting on the linear model, passing in the contrast matrix contr.matrix; Bayesian estimation: The Bayesian function is used to perform Bayesian estimation on the fitted model to obtain the significance of genes; Diagnostic plots were plotted using the plotSA function to assess the model's suitability and stability. Significance test: Use the decideTests function to perform a significance test on the model and output a summary of the significance test results for the genes; Differential expression gene analysis: The topTable function is used to obtain the differential expression of all genes based on the Bayesian estimation results, and then genes with significant differential expression are selected based on the given thresholds padj and foldChange; Results output: Write the screened differentially expressed genes into a CSV file; DEGs functional annotation and enrichment analysis: GO functional annotation and KEGG enrichment analysis were performed on the screened differentially expressed genes (DEGs) using the R package clusterProfiler. Enrichment results were considered significant only when the p-value was less than 0.

05. Specifically: Gene ID conversion: Use the bitr function to convert the gene symbol from SYMBOL to ENTREZID; GO analysis: GO enrichment analysis is performed using the enrichGO function, specifying the gene ID type, p-value, and q-value threshold parameters. The results are saved as a CCGO.rdata file. KEGG analysis: KEGG enrichment analysis is performed using the enrichKEGG function, specifying the threshold parameters for species, p-value, and q-value, and the results are saved as a CCKEGG.rdata file; GO / KEGG enrichment visualization: Use the barplot and dotplot functions to visualize the GO and KEGG enrichment results and save the images as PDF files; Gene set comparison analysis: The enrichKEGG function is used to perform KEGG enrichment analysis on genes in different groups, compare the KEGG enrichment results of different groups, and visualize the results.

7. The method for detecting and classifying immune microenvironment differences in medulloblastoma as described in claim 1, characterized in that, The method for validating immune microenvironment typing and biological feature screening through single-cell omics analysis in step S4 is as follows: Single-cell data were collected from samples of patients with medulloblastoma, and Seurat was used for quality control and dimensionality reduction clustering of the single-cell data. Based on the collected annotation information, the cell subpopulations of dimensionality reduction clustering are annotated; The R package scRNAtoolVis was used to obtain the distribution of specific markers in single-cell subsets and the differences in single-cell levels among medulloblastoma patients with different immune microenvironment subtypes. The method for collecting single-cell data from samples of medulloblastoma patients and using Seurat for quality control and dimensionality reduction clustering of single-cell data is as follows: For each sample, the number of UMIs, the total number of genes, the number of mitochondrial genes, and the number of ribosomal genes were counted. Cells with a total number of genes greater than 6000 or a number of genes less than 200, as well as cells with a mitochondrial gene ratio greater than 30%, were filtered out. For each sample, the top 2000 genes with the greatest variation were identified based on the mean and dispersion of all genes for integrated analysis to eliminate batch effects between samples. Principal component analysis (PCA) was then performed on the integrated data to reduce dimensionality and retain the top 20 principal components to capture the main changes in the data. HGV was detected using Seurat's pipeline. The average expression and dispersion of each gene were calculated. Genes were placed into bins, and the z-score for dispersion within each bin was calculated. A z-score of 0.5 was used as the cutoff value for dispersion, and a lower cutoff of 0.0125 and a higher cutoff of 3.0 were used as the average expression. Principal component analysis (PCA) was used for linear dimensionality reduction. Elbow Plot and Jackstraw methods were used to select statistically significant principal components. Clustering was visualized using Seurat based on UMAP. The specific steps are as follows: JackStraw: Use the JackStraw method to evaluate the results of PCA and identify significant principal components; ScoreJackStraw: Scores the results of the JackStraw method to identify statistically significant principal components; JackStrawPlot and ElbowPlot: Visualize the results of the JackStraw method and select the number of principal components to retain. JackStrawPlot is used to observe the p-value of each principal component, while ElbowPlot is used to observe the changes in the variance explained. Select the number of principal components at the "elbow". Cell clustering: FindNeighbors: Calculates the proximity relationships between cells based on PCA results; FindClusters: Based on the K-nearest neighbor graph, a clustering algorithm is used to divide cells into different cell clusters; Idents and table: View the clustering results. Use Idents to see which cluster each cell belongs to, and use table to count the number of cells in each cluster. UMAP dimensionality reduction and visualization: RunUMAP: Uses the UMAP algorithm to reduce the dimensionality of data. UMAP uses the stochastic gradient descent optimization algorithm, based on the minimum spanning tree and Gaussian mixture model approximation method, and adopts a distance-based weighting function to adjust the similarity weights between different data points, mapping high-dimensional data to two-dimensional or three-dimensional space. Minimize the loss function between the distances between data points in high-dimensional space and the distances between corresponding points in low-dimensional space to optimize the data representation in low-dimensional space; DimPlot: Based on UMAP dimensionality reduction, it plots the distribution of cells, visually demonstrating the relationships and distribution of different cell clusters; UMAP clustering visualization results show that individual cells from 26 samples from SHH, GP3, and GP4 are clustered into 21SHH, 18GP3, and 15GP4 subpopulations, respectively.

8. The method for detecting and classifying immune microenvironment differences in medulloblastoma as described in claim 1, characterized in that, Based on the collected annotation information, the method for annotating the cell subpopulations of dimensionality reduction clustering is as follows: scRNA sequencing was performed on MB patient samples and single matched samples from relapsed patients. Quality-controlled cells were projected into a two-dimensional UMAP map, and cell subpopulations in this experiment were manually annotated based on marker genes from the Cell Marker database and other brain tumor-related literature.

9. A system for detecting and classifying differences in the immune microenvironment of medulloblastoma, characterized in that, The device includes a processing unit that performs the method described in any one of claims 1-8 to detect and classify the immune microenvironment differences in medulloblastoma.

Images