Immune Microenvironment Analysis and Establishment of Response Prediction Model for PD-1 Blockade Immunotherapy in Glioblastoma Based on Transcriptome Deconvolution

doi:10.21203/rs.3.rs-2951425/v1

Download PDF

Research Article

Immune Microenvironment Analysis and Establishment of Response Prediction Model for PD-1 Blockade Immunotherapy in Glioblastoma Based on Transcriptome Deconvolution

https://doi.org/10.21203/rs.3.rs-2951425/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 05 Jul, 2023

Read the published version in Journal of Cancer Research and Clinical Oncology →

You are reading this latest preprint version

Purpose

Only a small proportion of GBM patients obtain survival benefit from PD-1 blockade immunotherapy due to the highly heterogeneous and suppressed immune micro-environment of GBM. We aimed at revealing the characteristics of tumor micro-environment (TME) of GBM related to response to PD-1 inhibitors and constructing a response prediction model for screening patients possibly benefit from PD-1 inhibitors.

Methods

Based on the composition and expression profiles of cell subpopulations calculated by deconvoluting the GBM bulk RNA-seq, differentially expressed genes analysis and gene set enrichment analysis (GSEA) were performed to explore genes and pathways related to response to PD-1 inhibitors. Further by combining least absolute shrinkage and selection operator (LASSO) regression and expression correlation with PD-L1, the response prediction genes of PD-1 inhibitors were identified and the response model was constructed through binary logistic regression.

Results

The comparison of abundances of infiltrating immune cells showed that the abundance of M0 macrophages of responders was lower while the abundance of activated DCs was higher before PD-1 inhibitors treatment; the abundances of plasma cells and M0 macrophages of responders were lower after PD-1 inhibitors treatment. In addition, GSEA showed that the main up-regulation pathways in the tumor microenvironment of responders before PD-1 inhibitors treatment included the regulation of T-helper 1 type immune response, the positive regulation of natural killer cell mediated cytotoxicity, p53 signaling pathway, homotypic cell-cell adhesion and etc., the main down-regulation pathways include the activation pathways of microglia and myeloid leukocytes, Ras signaling pathway and etc. Afterwards, ITGAX, LRRFIP1 and FMN1 were identified as the key response prediction genes of PD-1 inhibitors and the response prediction model based on them showed good predictive performance with potential value of clinical application in its validation and verification.

Conclusions

ITGAX, LRRFIP1 and FMN1 were identified as the response prediction genes of PD-1 inhibitors and the response prediction model based on them was proved to have potential clinical value.

Glioblastoma

PD-1 inhibitors

Transcriptome analysis

Tumor microenvironment

Drug response prediction model

Glioblastoma, the most common tumor in central nervous system with an annual incidence rate of 3.19/100,000 (Ostrom et al., 2013), also the most malignant glioma classified as WHO grade IV, is almost an incurable cancer with a median overall survival down to 14.6 months (Stupp et al., 2005). In recent years, PD-1/PD-L1 inhibitors showed incomparable curative effect in melanoma, lung cancer, primary hepatocellular carcinoma and etc., while the response rate of PD-1 inhibitors is only approximately 8% in GBM (Ott et al., 2019). Nevertheless, recurrent patients responded to PD-1 inhibitors showed a significant improvement in median overall survival, reaching 14.3 months according a retrospective clinical study (Zhao et al., 2019). Besides, 2 newly diagnosed patients involved in the research of Kurt A. Schalper (Schalper et al., 2019) survived over 33 months and 28 months respectively following neoadjuvant treatment with nivolumab.

Studies further showed that it is the highly heterogeneous and suppressed immune micro-environment of GBM that mainly influences the efficacy of PD-1/PD-L1-based immunotherapy (Bikfalvi et al., 2023) which indicates the vital importance of revealing the TME of GBM. RNA-seq has become an effective method to analyze TME at transcription level for its advantages of high throughput and low requirements of samples. However, bulk RNA-seq shows only the average expression of mixed cells in bulk tissue and masks heterogeneous expression among cells and subpopulations, scRNA-seq can precisely reveal cellular heterogeneity of TME at single cell level, its application is restricted by high cost though. Deconvolution algorithms represented by Cibersortx (Newman et al., 2019), developed to calculate abundances and expression of cell subpopulations within bulk RNA-seq, compensated for the limitations of bulk RNA-seq to a certain degree and provided a new approach to analyzing characteristics of GBM immune micro-environment related to response to PD-1 inhibitors. Cibersortx has been wildly used in the analyses of micro-environment of various cancers and has been proved to be one of the most advanced transcriptome deconvolution algorithms (Barroso-Sousa et al., 2022; Newman et al., 2019; Zhang et al., 2022).

Screening patients possibly benefit from ICIs by pre-treatment evaluation based on biomarkers is of great importance for improving clinical efficiency of immune checkpoint blockade therapy, while there is no reliable ICI response biomarker for GBM currently. In the present study, we deconvoluted the bulk RNA-seq of GBM tissue from patients treated with PD-1 inhibitors. On this basis, we sought the characteristics of GBM immune micro-environment that might influence the efficacy of PD-1 inhibitor treatment, screened response prediction genes and constructed a prediction model.

Data acquisition

The bulk RNA-seq data of GBM from 17 patients receiving PD-1 inhibitors and related clinical information used in this study were provided by Junfei Zhao (Zhao et al., 2019), which included 16 pre-treatment samples and 9 post-treatment samples.

The scRNA-seq data of GBM and related metadata used in this study was obtained from Gene Expression Omnibus database (GEO) by “GSE84465”.

The bulk RNA-seq data of GBM and related clinical information of the external cohort in the survival analysis of this study were provided by Robert M Prins (Cloughesy et al., 2019).

scRNA-seq analysis

This part of the study was mainly conducted by the “Seurat” R package V4.1.0. Genes and cells were firstly filtrated based on the standards–genes detected in less than 3 cells and cells with less than 50 detected genes. Then, further filtration was performed to remove genes and cells with outlier of nCount RNA, nFeature RNA, percent.MT or percent.ERCC by the standards–100 < nFeatureRNA < 7500, 10000 < nCount RNA < 1300000 and percent.ERCC < 25%, which were determined by plotting the distribution map of these parameters. Then, the gene expression of the remaining cells was normalized by “LogNormalize” method to eliminate the influence of the sequencing depth among different cells. Next, the top 1500 highly variable genes were identified by “vst” method to perform principal component analysis (PCA) following Z-score conversion and 20 initial PCs were determined to be involved in the dimensionality reduction of t-distributed stochastic neighbor embedding (tSNE), as a result, all the cells were classified into 13 clusters. To identify the cell type of each cluster, the differentially expressed genes analysis among all the genes within clusters was performed by Wilcoxon test with the standard–adjusted P value < 0.05 and |log2 (Fold Change)| > 0.5. On this basis, cell type of each cluster was identified by comparing top differentially expressed genes of each cluster with cell markers from the CellMarker database (Hu et al., 2023).

Deconvolution analysis of bulk RNA-seq

Deconvolution analysis in this study was conducted by Cibersortx (Newman et al., 2019), which aimed to analyze bulk RNA-seq (Mixture) into composition and expression of cell subpopulations taking a transcriptome map (Signature Matrix) of corresponding type of cells as reference.

The expression profile of the GBM scRNA-seq was annotated with cell types based on the previous analysis and uploaded to Cibersortx to create Signature Matrix. Meanwhile, the expression profile of the bulk RNA-seq based on the GBM tissue from the patients with PD-1 inhibitor therapy was firstly normalized from Count to CPM (Counts per million) and then uploaded to Cibersortx as Mixture. Afterwards, the composition and expression of cell subpopulations of the Mixture were calculated using Cibersortx taking the created Signature Matrix of GBM and the transcriptome map of lymphocytes and myeloid cells merged from transcriptome profile of 22 types of human immune cells as reference respectively.

Differentially expressed genes analysis

Differentially expressed genes analysis based on deconvolution was performed with the “Limma” R package V3.50.3. The selection criteria for differentially expressed genes is |log2 FoldChange|> 1 and p < 0.05.

Gene set enrichment analysis (GSEA)

GSEA was performed by the “clusterProfiler” R package V4.2.2. Genes sorted in descending order by log2 FoldChange based on differentially expressed genes analysis

were enriched taking the gene sets of GO (BP) and KEGG as reference with the following criterion: p.adj (Benjamin Hochberg) < 0.25.

Selection of response prediction genes

The response prediction genes of PD-1 inhibitors were firstly screened by least absolute shrinkage and selection operator (LASSO) regression with threefold cross-validation to determine the optimal regularization parameter (λ.min). Then, the response prediction genes were further selected by the expression correlation analysis with PD-L1 on the online analysis platform of GEPIA2 (Tang et al., 2019) with the following criteria: |Spearman correlation coefficient|<0.3 and Spearman correlation p value < 0.05. Next, the expression correlation between response prediction genes of the immune cells in GBM was conducted to remove gene with multi-collinearity under the same criteria.

Construction of the response prediction model

The response prediction model was constructed through binary logistic regression using the “glm” function of R software with the “family” parameter of “binomial”. Variance inflation factor (VIF) of each gene in the regression model was calculated using the “car” R package V3.0-13 to detect genes with multi-collinearity. The nomogram of the model was plotted using the “regplot” R package V1.1.

Validation and verification of the response prediction model

The ROC, Hosmer-Lemeshow test and calibration curve were calculated respectively using the R package “pROC” V1.18.0, “ResourceSelection” V0.3-5 and “riskRegression” V2022.09.23, tested by leave-one-out cross-validation (LOOCV). The decision curve and clinical impact curve were calculated using the “rmda” R package V1.6, tested by LOOCV. The accuracy value and Kappa value were calculated using the “caret” R package V6.0-9.2, tested by LOOCV.

Kaplan-Meier curves of the two groups divided by the median response probability according to the model were plotted using the the R package “survival” V3.3-1 and “survminer” V0.4.9. Logrank test was performed to compare the differences of survival distributions between the two groups.

Difference in abundances of infiltrating immune cells in GBM between PD-1 inhibitor responders and non-responders

To explore the relationship between abundances of infiltrating immune cells in GBM and response to PD-1 inhibitors, the component fractions of 22 types of immune cells in the TME of each patient before and after anti-PD-1 treatment were calculated by deconvoluting the bulk RNA-seq taking the transcriptome map of 22 types of human immune cells (Newman et al., 2019) as reference.

The comparison of abundances of infiltrating immune cells between PD-1 inhibitor responders (n = 9) and non-responders (n = 7) before PD-1 inhibitor treatment showed that the abundance of M0 macrophages of responders was lower while the abundance of activated DCs was higher (Fig. 1A; Wilcoxon test; M0 macrophages, p = 0.047; activated DCs, p = 0.029), which suggest that low abundance of M0 macrophages and high abundance of activated DCs in the TME of GBM might be predictive factors for response to PD-1 inhibitors.

The comparison of abundances of infiltrating immune cells between PD-1 inhibitor responders (n = 7) and non-responders (n = 2) after PD-1 inhibitor treatment showed that the abundances of plasma cells and M0 macrophages of responders was lower (Fig. 1B; Wilcoxon test; plasma cells, p = 0.024; M0 macrophages, p = 0.037), which suggest that macrophages and antibodies might play a certain role in the treatment of GBM with PD-1 inhibitors.

Genes expression and pathway regulations related to response to PD-1 inhibitors

Differentially expressed genes analysis was considered to further excavate genes that are predictive for response to PD-1 inhibitors. Taking into consideration of the high heterogeneity of the micro-environment of GBM, we first constructed the TME transcriptome map of 9 types of cells by analyzing a scRNA-seq data from 4 GBM patients (Darmanis et al., 2017) (Fig. 2).

Combing that and the transcriptome map of 22 types of human immune cells as reference, the expression profiles of neoplastic cells, lymphocytes, and myeloid cells within the bulk RNA-seq of GBM treated with PD-1 inhibitors was calculated by deconvolution to reduce the interference caused by the heterogeneity between cell subpopulations. With that, the differentially expressed genes analysis was performed and 384 differentially expressed genes of neoplastic cells (Fig. 3A), 81 differentially expressed genes of lymphocytes (Fig. 2B) and 205 differentially expressed genes of myeloid cells (Fig. 3C) was identified.

Gene set enrichment analysis (GSEA) was performed to explore changes of pathways in the micro-environment of responders compared to non-responders before PD-1 blockade treatment based on the differential gene expression analysis. GSEA of neoplastic cells showed that the up-regulation pathways of responders mainly included the regulation of T-helper 1 type immune response, the positive regulation of natural killer cell mediated cytotoxicity and p53 signaling pathway; the down-regulation pathways mainly included the activation of microglia and myeloid leukocytes and Ras signaling pathway (Fig. 4A). GSEA of lymphocytes showed that the up-regulation pathways in responders mainly included the positive regulation of cell population proliferation and the regulation of cellular amide metabolic process (Fig. 4B). GSEA of myeloid cells showed that the up-regulation pathways in responders mainly included cell-substrate junction, integrin mediated signaling pathway, homotypic cell-cell adhesion and leukocyte mediated cytotoxicity (Fig. 4C).

Identification of predictive genes of response to PD-1 inhibitors in GBM

With the expression profiles of immune cells of 27 pre-treatment samples from 17 patients (all pre-samples of the same patient were included to expand the sample size) based on deconvoluting the GBM bulk RNA-seq and the differentially expressed genes of the immune cells, 15 predictive genes of response to PD-1 inhibitors with nonzero coefficients, ASS1, NR1H3, ITGAX, CPEB1, LRRFIP1, FKBP5, DNAJC5B, PI4K2A, TESK2, USB1, HIAT1, NCOA7, FMN1, GALC and SPNS3, were screened by least absolute shrinkage and selection operator (LASSO) regression (Fig. 5).

Preclinical and clinical evidences suggested a strong expression correlation between PD-L1 and the response to PD-1 inhibitors in multiple malignant tumors (Davis and Patel, 2019; Lantuejoul et al., 2020), according to which four predictive genes highly correlated to PD-L1 at transcription level were further screened by expression correlation analysis (Table 1). The four predictive genes were ITGAX, LRRFIP1, PI4K2A and FMN1, PI4K2A of which was next removed in order to select genes as independent factors of response to PD-1 inhibitors due to the multi-collinearity between PI4K2A and LRRFIP1 that was found by calculating the variance inflation factor (ITGAX, VIF = 1.482; LRRFIP1, VIF = 5.879 > 5; FMN1, VIF = 2.147; PI4K2A, VIF = 8.122 > 5) of each gene in the linear regression and analyzing the expression correlation between them (Fig. 6).

Table 1

Expression correlation between 15 candidate response prediction genes and PD-L1 in GBM
Gene	Spearman correlation coefficient (ρ)	Spearman correlation P value
LRRFIP1	0.47	$\text{2×}{\text{10}}^{\text{-10}}$
ITGAX	0.36	$\text{3.3×}{\text{10}}^{\text{-6}}$
FMN1	0.35	$\text{4.2×}{\text{10}}^{\text{-6}}$
PI4K2A	0.32	$\text{3.7×}{\text{10}}^{\text{-5}}$
GALC	0.3	$\text{8.2×}{\text{10}}^{\text{-5}}$
HIAT1	0.3	$\text{1×}{\text{10}}^{\text{-4}}$
FKBP5	0.29	0.00018
USB1	0.28	0.00029
NR1H3	0.24	0.0018
NCOA7	0.23	0.0037
CPEB1	0.22	0.0047
SPNS3	0.21	0.0082
DNAJC5B	0.051	0.52
ASS1	-0.0073	0.93
TESK2	-0.0066	0.93

Establishment and validation of response prediction model for PD-1 inhibitors in GBM

Three genes, ITGAX, LRRFIP1 and FMN1, was finally determined as the key genes predictive of response to PD-1 inhibitors in GBM and the response prediction model was constructed then by binary multiple logistic regression (Table 2, Fig. 7).

The probability for GBM patients to respond to PD-1 inhibitors can be calculated by the following formula according to the model:

$$P=\frac{1}{1+{e}^{-(-114.867+18.560\times （ITGAX）+473.179\times （LRRFIP1）-36.500\times （FMN1）)}}$$

Transcriptome expression quality of ITGAX, LRRFIP1 and FMN1 in the formula can be obtained from the expression profiles of immune cells calculated by deconvoluting the bulk RNA-seq of GBM tumor tissue or directly from the bulk RNA-seq of immune cells sorted from GBM tumor tissue, following the max-min normalization as below according to the minimum value ${\text{m}\text{i}\text{n}\{x}_{i}\}$ and the maximum value ${\text{m}\text{a}\text{x}\{x}_{i}\}$ of ${x}_{i}$--the expression value of ITGAX, LRRFIP1 and FMN1:

$${x}_{i}^{*}=\frac{{x}_{i}-{\text{m}\text{i}\text{n}\{x}_{i}\}}{{\text{m}\text{a}\text{x}\{x}_{i}\}-{\text{m}\text{i}\text{n}\{x}_{i}\}}$$

Table 2

Results of the multivariate logistic regression
Variable	Estimated regression coefficient	Z value of Wald’s test for regression coefficient	P value of Wald’s test for regression coefficient	Odds ratio (OR)
ITGAX	18.560	2.182	0.029	1.047 (1.005–1.091)
LRRFIP1	473.179	2.382	0.017	3.215 (1.23–8.406)
FMN1	-36.500	-1.793	0.073	0.914 (0.828–1.008)

Then we performed leave-one-out cross validation (LOOCV) to test the model. The area under the curve (AUC) of receiver operating characteristic (ROC) was 0.922 (95% CI 0.791-1), suggesting a good discriminative ability (Fig. 8A). The P-value of the Hosmer-Lemeshow (H-L) test was 0.574 and the calibration curves demonstrated a high consistency between prediction and actual observation, indicating a high prediction accuracy (Fig. 8B).

Clinical utility of the response prediction model

The clinical utility of the model was evaluated by the decision curve analysis (DCA) and the clinical impact curve (CIC).

The DCA curve showed that carrying out PD-1 blockade immunotherapy according to the model when the threshold probability was lower than 82% provided greater net benefit than the “treat all” or “treat none” strategies, indicating the high clinical value of the model (Fig. 9A).

The CIC showed that the population responsive to PD-1 inhibitors identified by the model was highly consistent with the actual population responded to PD-1 inhibitors when the threshold probability was higher than 70%, confirming the high clinical efficiency of the model (Fig. 9B).

Survival differences between patients with PD-1 blockade treatment classified by the model

The ultimate significance of drug response prediction models lies in whether patients can benefit from their application in survival time and survival outcomes. Therefore, we performed survival analysis to evaluate the survival significance of the model. The probabilities of response to PD-1 inhibitors of each patient in the internal cohort and the external cohort were calculated according to the bulk RNA-seq of tumor tissue sampled before (internal cohort) or on the early stage (external cohort) of the treatment by the model. The cohorts were next respectively divided into low predictive response group and high predictive response group by the median response probability to plot survival curves using Kaplan-Meier method and compare differences of the survival distributions using Log-rank test. The survival analysis of the internal cohort (n = 16) showed that the survival rate of the high predictive response group was higher than that of the low predictive response group (Log rank test, p = 0.033) (Fig. 10A). The survival analysis of the external cohort (n = 26) showed that the survival rate of the high predictive response group was markedly higher than that of the low predictive response group (Log rank test, p = 0.005) with a median survival time two time over that of the low predictive response group (Fig. 10B).

The difference of the distribution of survival time between the two groups was consistent with the prediction from the model both in the survival analyses of the internal cohort and the external cohort which indicated a certain predictive value for survival in GBM patients with PD-1 inhibitor treatment.

GBM is one of the most malignant refractory tumors with rapid growth, strong invasiveness and poor prognosis. PD-1 blockade immunotherapy was proved effective in extending the survival of GBM patients in multiple clinical researches while the response rate was unsatisfactory (Omuro et al., 2022; Ott et al., 2019; Schalper et al., 2019; Zhao et al., 2019). Searching for biomarkers predictive for response to PD-1 inhibitors is a major problem demanding prompt solution in the clinical practice of treating GBM. Further research showed that the efficacy of ICIs was associated with the special micro-environment of GBM (Bikfalvi et al., 2023), which was proved heterogeneous, immune-suppressed and complex, while relevant mechanisms were still remaining to be investigated. With the development of high-throughput sequencing technology, large amounts of transcriptome data were generated to study TME which greatly contributed to develop molecular markers and targets for tumor therapy, driving the growth and implementation of precision medicine.

Transcriptome deconvolution by Cibersortx (Newman et al., 2019), aimed to analyze bulk RNA-seq into composition and expression of cell subpopulations, showed excellent calculation performance in the analyses of microenvironment of multiple tumors, for instance, melanoma, liver cancer and breast cancer (Barroso-Sousa et al., 2022; Newman et al., 2019; Zhang et al., 2022). Moreover, the outstanding calculation accuracy and robustness of Cibersortx have been proved in several studies by comparing it with mainstream deconvolution algorithms and verifying it through flow cytometry and scRNA-seq in analyzing the cell compositions of peripheral blood mononuclear cell from healthy people, lymph nodes from patients with follicular lymphoma, micro-environment of colorectal cancer, squamous cell carcinoma of head and neck (HNSCC) and melanoma (Newman et al., 2015; Wang et al., 2023), which indirectly supported for the reliability of deconvolution calculation based on GBM scRNA-seq analysis in our research.

The response prediction model performed well in the validation and testing of discrimination, calibration, clinical utility and survival significance, showing potential clinical value. However, what should be noted is that the limitation of the model lies in lack of validation in large independent cohort which is key to improving its clinical predictive capability. With extensive development of GBM research, more and more RNA-seq of GBM is going to be performed and released, which will provide more powerful support for the improvement of the model.

In order to extend the application of the model, we evaluated it in melanoma and lung cancer cohorts with PD-1 blockade treatment, but unfortunately, the model showed no prediction effect (AUC = 0.5), which could be because the huge differences among these cancers, such as gene expression, protein function, tissue morphology and etc.

Cibersortx provided a new approach to analyze the immune micro-environment of GBM related to response to PD-1 inhibitors. Nevertheless, it still could not analyze bulk RNA-seq into single cell level comparing to scRNA-seq which was recognized as the best means to study heterogeneity of TME at transcription level. We expected that scRNA-seq of large cohort will depict the micro-environment of GBM more comprehensively and accurately in the future.

In the present study, the comparison of abundances of infiltrating immune cells between responders and non-responders suggested that high abundances of M1 macrophages and activated DCs in the micro-environment of GBM might be related to response to PD-1 inhibitors. In addition, the gene set enrichment analysis indicated that the potential mechanisms of PD-1 inhibitors in treating GBM might involve the down-regulation of Ras pathway, the up-regulation of p53 pathway and inhibiting tumor development by regulating tumor immune micro-environment and affecting adhesion between tumor cells and extracellular matrix mediated by integrin and cadherin, aside from activating T cells through binding PD-L1. Furthermore, ITGAX, LRRFIP1 and FMN1 were identified as key response prediction genes and the response prediction model based on them was validated and verified to have good predictive performance and potential clinical value.

bulk RNA-seq: bulk RNA Sequencing

CIC: Clinical Impact Curve

DC: Dendritic Cell

DCA: Decision Curve Analysis

GAM: GBM Associated Microglia/Macrophage

GBM: Glioblastoma

GEO: Gene Expression Omnibus

GO: Gene Ontology

GSEA: Gene Set Enrichment Analysis

KEGG: Kyoto Encyclopedia of Genes and Genomes

LASSO: Least Absolute Shrinkage and Selection Operator

LOOCV: Leave-one-out Cross-validation

NK: Natural Killer Cell

NSCLC: Nonsmall-cell Lung Cancer

PC: Principal Component

PCA: Principal Component Analysis

PD-1: Programmed Cell Death Protein 1

PD-L1: Programmed Cell Death 1 Ligand 1

RNA-seq: RNA Sequencing

ROC: Receiver Operating Characteristic

scRNA-seq: Single Cell RNA Sequencing

Th: Helper T Cells

TME: Tumor Microenvironment

Treg: Regulatory T Cell

tSNE: t-distributed Stochastic Neighbor Embedding

VIF: Variance Inflation Factor

Data availability statement

The data of this study are available from the lead author upon reasonable request.

Barroso-Sousa R., Forman J., Collier K.et al (2022). Multidimensional molecular profiling of metastatic triple-negative breast cancer and immune checkpoint inhibitor benefit. JCO Precis Oncol 6: e2100413.
Bikfalvi A., Da Costa C.A., Avril T.et al (2023). Challenges in glioblastoma research: focus on the tumor microenvironment. Trends Cancer 9: 9-27.
Cloughesy T.F., Mochizuki A.Y., Orpilla J.R.et al (2019). Neoadjuvant anti-PD-1 immunotherapy promotes a survival benefit with intratumoral and systemic immune responses in recurrent glioblastoma. Nat Med 25: 477-486.
Darmanis S., Sloan S.A., Croote D.et al (2017). Single-cell RNA-seq analysis of infiltrating neoplastic cells at the migrating front of human glioblastoma. Cell Rep 21: 1399-1410.
Davis A.A., and Patel V.G. (2019). The role of PD-L1 expression as a predictive biomarker: an analysis of all US Food and Drug Administration (FDA) approvals of immune checkpoint inhibitors. J Immunother Cancer 7: 278.
Hu C., Li T., Xu Y.et al (2023). CellMarker 2.0: an updated database of manually curated cell markers in human/mouse and web tools based on scRNA-seq data. Nucleic Acids Res 51: D870-D876.
Lantuejoul S., Sound-Tsao M., Cooper W.A.et al (2020). PD-L1 testing for lung cancer in 2019: Perspective from the IASLC pathology committee. J Thorac Oncol 15: 499-519.
Newman A.M., Liu C.L., Green M.R.et al (2015). Robust enumeration of cell subsets from tissue expression profiles. Nat Methods 12: 453-457.
Newman A.M., Steen C.B., Liu C.L.et al (2019). Determining cell type abundance and expression from bulk tissues with digital cytometry. Nat Biotechnol 37: 773-782.
Omuro A., Reardon D.A., Sampson J.H.et al (2022). Nivolumab plus radiotherapy with or without temozolomide in newly diagnosed glioblastoma: Results from exploratory phase I cohorts of CheckMate 143. Neurooncol Adv 4: vdac025.
Ostrom Q.T., Gittleman H., Farah P.et al (2013). CBTRUS statistical report: Primary brain and central nervous system tumors diagnosed in the United States in 2006-2010. Neuro Oncol 15 Suppl 2: ii1-56.
Ott P.A., Bang Y.J., Piha-Paul S.A.et al (2019). T-cell-inflamed gene-expression profile, programmed death ligand 1 expression, and tumor mutational burden predict efficacy in patients treated with Pembrolizumab across 20 cancers: Keynote-028. J Clin Oncol 37: 318-327.
Schalper K.A., Rodriguez-Ruiz M.E., Diez-Valle R.et al (2019). Neoadjuvant nivolumab modifies the tumor immune microenvironment in resectable glioblastoma. Nat Med 25: 470-476.
Stupp R., Mason W.P., Van Den Bent M.J.et al (2005). Radiotherapy plus concomitant and adjuvant temozolomide for glioblastoma. N Engl J Med 352: 987-996.
Tang Z., Kang B., Li C.et al (2019). GEPIA2: an enhanced web server for large-scale expression profiling and interactive analysis. Nucleic Acids Res 47: W556-W560.
Wang X., Xie C., and Lin L. (2023). Development and validation of a cuproptosis-related lncRNA model correlated to the cancer-associated fibroblasts enable the prediction prognosis of patients with osteosarcoma. J Bone Oncol 38: 100463.
Zhang S., Bacon W., Peppelenbosch M.P.et al (2022). Deciphering tumour microenvironment of liver cancer through deconvolution of bulk RNA-seq data with single-cell atlas. Cancers (Basel) 15: 153.
Zhao J., Chen A.X., Gartrell R.D.et al (2019). Immune and genomic correlates of response to anti-PD-1 immunotherapy in glioblastoma. Nat Med 25: 462-469.

No competing interests reported.

Download PDF

Journal Publication

published 05 Jul, 2023

Read the published version in Journal of Cancer Research and Clinical Oncology →

Editorial decision: Accepted
26 Jun, 2023
Reviewers agreed at journal
19 May, 2023
Reviewers invited by journal
19 May, 2023
Editor assigned by journal
19 May, 2023
Submission checks completed at journal
19 May, 2023
First submitted to journal
18 May, 2023

You are reading this latest preprint version

Immune Microenvironment Analysis and Establishment of Response Prediction Model for PD-1 Blockade Immunotherapy in Glioblastoma Based on Transcriptome Deconvolution

Status:

Journal Publication

Version 1

Abstract

Purpose

Methods

Results

Conclusions

Figures

Introduction

Materials and methods

Data acquisition

scRNA-seq analysis

Deconvolution analysis of bulk RNA-seq

Differentially expressed genes analysis

Gene set enrichment analysis (GSEA)

Selection of response prediction genes

Construction of the response prediction model

Validation and verification of the response prediction model

Results

Difference in abundances of infiltrating immune cells in GBM between PD-1 inhibitor responders and non-responders

Genes expression and pathway regulations related to response to PD-1 inhibitors

Identification of predictive genes of response to PD-1 inhibitors in GBM

Establishment and validation of response prediction model for PD-1 inhibitors in GBM

Clinical utility of the response prediction model

Survival differences between patients with PD-1 blockade treatment classified by the model

Discussion

Conclusions

Abbreviations

Declarations

References

Additional Declarations

Status:

Journal Publication

Version 1

Gene	Spearman correlation coefficient (ρ)	Spearman correlation P value
LRRFIP1	0.47	\(\text{2×}{\text{10}}^{\text{-10}}\)
ITGAX	0.36	\(\text{3.3×}{\text{10}}^{\text{-6}}\)
FMN1	0.35	\(\text{4.2×}{\text{10}}^{\text{-6}}\)
PI4K2A	0.32	\(\text{3.7×}{\text{10}}^{\text{-5}}\)
GALC	0.3	\(\text{8.2×}{\text{10}}^{\text{-5}}\)
HIAT1	0.3	\(\text{1×}{\text{10}}^{\text{-4}}\)
FKBP5	0.29	0.00018
USB1	0.28	0.00029
NR1H3	0.24	0.0018
NCOA7	0.23	0.0037
CPEB1	0.22	0.0047
SPNS3	0.21	0.0082
DNAJC5B	0.051	0.52
ASS1	-0.0073	0.93
TESK2	-0.0066	0.93