Identification and validation of a prognostic model based on immune-related genes in ovarian carcinoma

Data extraction and preprocessing

UCSC Xena database (https://xena.ucsc.edu/) collected RNA-seq expression and clinical data from The Cancer Genome Atlas (TCGA) dataset and GTEx, including 379 tumor samples and 88 normal controls. In the present study, the TCGA ovarian cancer (OC) cohort and GTEx ovarian cohort were employed as the train cohort. The RNA-seq expression and clinical-pathological parameters of 85 OC patients were retrieved from International Cancer Genome Consortium (ICGC, https://dcc.icgc.org/). To operate the study data, the microarray probes with null gene test value were first eliminated and the remaining probes were mapped to human genes. The tumor samples absence of clinical information and survival data were excluded. Details of the steps were summarized in Fig. 1 through a flowchart.

Screening for immune gene signature related to prognosis

The immune-related genes were downloaded from GSEA database (http://www.gsea-msigdb.org/gsea/index.jsp). The search keyword was “immune” and organism was “Homo sapiens”. A total of 458 immune-related genes were included in the subsequent analysis (Systematic name: M41715). Batch effects of the train dataset were removed using the “sva” package in R software. We then intersected all differentially expressed genes (DEGs) with immune genes to arrive at all of the immune DEGs in OC. Differentially expressed immune-related genes between normal and tumor tissues of the train dataset were identified using the Wilcoxon test (FDR < 0.05, |Log2(fold-change)| > 1). We additionally screened the immune DEGs related to prognosis through Univariate Cox Regression analysis. For the course of analysis, we used the R package “survival”. The screening criteria was p-value < 0.05. Furthermore, the immune gene correlation was used for the train cohort to predict integrative clusters and molecular subtype classifications through the “ConsensusClusterPlus” R package.

Description of prognosis-related molecular subtypes

Consensus clustering was conducted with the R software ConsensusClusterPlus package to distinguish the subgroups of OC according to those IRGs possessing of the highest variability. Initially, partial features and items were subsampled in the light of the data matrix by this algorithm, and all subsamples were classified as k groups based on the k-means. Their consensus values were calculated after repeating for the number specified by the user, and cluster stability was evaluated by using several clustering runs of algorithm. Subsequently, the proportion of cluster runs that divided into two items were regarded as the values of pairwise consensus, which were measured and saved within the consensus matrix for every k. Then, for every k, the 1-consensus value distance was applied to the eventual hierarchical agglomerative consensus clustering, which was reduced into k groups. Specially, the stability of clustering result was determined with the methods of the particular clustering approach into random data subsets, which was named as the “consensus” clustering (Hu, Yang & Sang, 2020; Wu et al., 2017). To generate more specific OC categories, more categories were required. Commonly, the consensus values set as 0 (white) and 1 (dark blue) were characterized with the color gradient. Furthermore, there were the same cluster items placed near each other in the matrix. As a result, the excellent consensus matrix exhibited the distinct color-coded heatmap. This heatmap consisted of blue blocks along the diagonal with the white background.

The construction of prognostic model

With the help of the “sva” package in R software, batch effects of the TCGA and ICGC dataset were deleted. Univariate Cox regression analysis for screening the prognostic genes with the “survival” R package was employed to constitute a reliable prognosis signature. In order to make a feature selection and screen the importance of prognostic-related genes with p < 0.05, the classic random forest algorithm was utilized. Thereafter, based on the parameters of the least absolute shrinkage and selection operator (LASSO) COX regression analysis, we build the prognostic signature with 1,000-fold cross-validation with the “glmnet” package in R software. Then the risk score of every patients was determined according to the following formula = ∑Coef_gene × Exp_gene, where Coef_gene denotes the each prognostic gene coefficient, the Exp_gene denotes the expression of each gene (Dong & Xu, 2019). After obtaining the median risk score, all patients were divided into high-risk and low-risk groups.

Evaluation of prognostic model

The time-dependent receiver operating characteristic (Rocconi et al., 2021) analysis with the “timeROC” package in R software was chosen to evaluate the predictive accuracy of the signature. Furthermore, the risk curve was created by the “pheatmap” R package to evaluate the prognostic model. In order to test the feasibility of prognostic model as an independent predictor, the indicators containing clinical traits and risk values in train cohort and test cohort was measured using univariate and multivariate Cox regression analyses.

Analysis for GO and KEGG enrichment pathway

GO analysis including cellular component, molecular function, and biological process to illustrate the multiple cellular functions of certain genes.

KEGG is an analysis method for uncovering the specific biological pathways of certain enriched genes. Therefore, GO and KEGG pathway analysis were conducted with the “clusterProfiler”, “org.Hs.eg.db”, “enrichplot”, and “ggplot2” R package according to DEGs with high-risk and low-risk groups of OC patients. The significant pathways were enriched with q-value < 0.05 (the corrected p-value).

Drug sensitivity analysis

In order to examine the drug sensitivity and tolerance effects of genes involved in developing the prognostic model, the transcriptome and FDA-certified drug sensitivity related data were downloaded from the CellMiner database (https://discover.nci.nih.gov/cellminer/). The Pearson correlation analysis was utilized to demonstrate the relationships between gene expression and drug sensitivity.

Estimation of immune cell infiltration

To examine the connection between the prognostic genes and immune cell infiltration, we analyzed TCGA cohort using online tool of GEPIA (http://gepia.cancer-pku.cn/) and TIMER (https://cistrome.shinyapps.io/timer/). Median expression value of prognostic genes was set as critical point, and divided into high and low expression groups. To further evaluate immune cell infiltration, the single sample gene set enrichment analysis (ssGSEA) was selected to study the different infiltration degrees of immune cell subgroups between above two high and low expression groups using the R package “GSVA”.

Single cell sequencing data download and processing

The scRNA-seq dataset GSE118828 was downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/). Total number of 16 OV samples were included in the cohort in addition to benign and normal tissues. The R package Seurat (version 4.3.4) (Butler et al., 2018) was used to conduct the process of quality control (QC). The data filtering parameters were set as follows: 1. Genes expressed in at least three single cells; 2. RNA in cells was greater than 50; 3. Mitochondrial genes less than 5%. At the bottom of QC, the cell populations were tagged based on marker genes. These markers were recognized as major immune cell types and visualized with the dimensionality reduction algorithm uniform manifold approximation and projection (UMAP).

Establishment of specific ovarian cancer cell line

Tissue sample of primary OC was obtained from Tianjin Medical University Cancer Institute and Hospital after approval of the Institutional Ethics Committee (NO. bc2023167) and informed written consent from the patient’s family. This clinical specimen was donated from a 61-year-old Chinese female patient with advanced ovarian carcinoma. The willingness to provide verbal informed consent was abstained. Tumor specimen was washed with sterile phosphate-buffered saline (PBS) for three times, and then re-suspended with RPMI-1640 medium containing 10% FBS, 5–10% cell-free ascites, penicillin G (100 U/ml), streptomycin (100 µg/ml). Subsequently, it was transferred into a 10 cm petri dish at 37 °C and 5% CO2. After the confluence rate up to of 80–90%, these primary cultured cells were sub-cultured at a ratio of 1:1. This cell line was considered as a continuous cell line after 50 passages and named W038. The detail morphological properties of W038 were monitored with the help of an inverted microscope. Meanwhile, cells treated with hematoxylin and eosins (HE) were detected using a light microscope.

Transfection

The GV298 lentiviral particles consisting of GBP1P1 shRNA sequence (5′-GCCAAGTCTG GTCACTAAACT-3′) were designed and synthesized by Jikai Kiin Technology Co. LTD (Shanghai, China). Lentiviral particles containing 5 µg/mL polybrene were added into W038 cells. After the transfection for 48 h, targeting cells were selected with 2 µg/mL of puromycin for two months. Then the resistant colonies were pooled, expanded, and named sh-GBP1P1 W038. The negative control cells with disordered shRNA (5′-TTCTCCGAACGTGTCACGT-3′) were named sh-CTRL W038. The cells of GBP1P1 knockdown were identified by quantitative real-time polymerase chain reaction (qRT-PCR), and fluorescence imaging using fluorescence microscopy (Lecia, Wetzlar, Germany). The minimum information for publication of quantitative real-time PCR experiments (MIQE) was listed in the Supplemental Material.

Immunohistochemistry staining

Immunohistochemical studies were performed on paraffin section with precipitation from cell centrifugation. The sections (4 µm) were heated at 60 °C for 1 h, deparaffinized in xylene, and rehydrated with graded ethanol. Then samples were heated in citrate buffer (pH 6.0) for 2 min to retrieve antigens. Endogenous peroxidase activity was quenched using a methanol and hydrogen peroxide bath for 20 min. After pretreatment, slides were loaded with the primary antibody (1:200) overnight at 4 °C. Next, the biotinylated secondary antibody named goat anti-human IgG (Santa Cruz Biotechnology, USA) was labeled with streptavidin–horseradish peroxidase (HRP) using a DAB staining kit (Maixin Biotechnology, China). Negative controls were performed by omission of primary antibody. Positive cells showed brownish yellow in the cytoplasm or the cell membrane. These images were acquired with an Olympus BX51 microscope. Five high-power fields were chosen randomly for histological evaluation. The 2 pathologists independently reviewed the antibody stained slides along with the HE slides 1 to 3 months later.

The assays of cell proliferation, invasion, and migration

The CCK-8 kit was carried out with a density of 2,000 cell/well in 96-well plates. Cells were then treated with 10 µL of CCK-8 solution (Dojindo Molecular Technologies, Inc., Rockville, MD, USA) for 60 min. Microplate spectrophotometer was used to determine the optical density at 450 nm of cells at different time. The invasion assay was performed as follows: 4,000 cells suspended in 200 µL of serum-free RPMI 1640 media were seeded into the upper chambers, which were coated with matrigel (1:5 diluted solution; BD Biosciences); then dried at 37 °C for 2 h; 650 µL RPMI 1640 media with 10% FBS was next added into the lower chambers. After the incubation for 48 h, the filtered cells were fixed in 4% paraformaldehyde and loaded with 1% crystal violet at 37 °C and 5% CO₂; the stained cells from five selected randomly views were counted with a microscope at a 200× magnification. For the migration assay, the ibidi Culture-Insert system was used according to the manufacturer’s in structions (ibidi, Martinsried, Germany).

Flow cytometry apoptosis assay

After addition with 3µg/mL carboplatin for 24 h, trypsin was used to digest OC cells which then washed with PBS for twice. Cells were suspended with 1× binding buffer at the density of 1 × 10⁶ cells/mL in flow tubes, and volume of cell solution was 100 µL. Thereafter Annexin V-IF488 (5 µL) and 7-AAD (5 µL) solutions were added into each tube. After the reaction for 15 min, another 400 µL of 1× binding buffer was added to stop the reaction. Eventually, Flow cytometry was used for on-board detection within 1 h. The percentage of apoptotic cells in each group was evaluated by calculating the sum of Annexin V + 7-AAD + and Annexin V + 7-AAD − cells. The flow cytometry files could be found here: https://doi.org/10.6084/m9.figshare.25991899.v1.

Statistical analysis

R version 4.3.4 was used to conduct the statistical analysis of this research. A difference of p value < 0.05 indicated statistical significance.

Determining the consensus clustering result

According to the expression profile of IRGs, we employed a consensus clustering algorithm to categorize the TCGA OC patients into de novo groups. In comparison with traditional k-mean and hierarchical clustering algorithms, consensus clustering was illustrated to be more robust and insensitive to random starts, and had been broadly used to recognize biologically meaningful clusters (Monti et al., 2003). A total of 243 differentially expressed IRGs (TCGA tumor vs GETx normal tissues, FDR < 0.05, |Log2 (fold-change) | > 1) were used for consensus clustering to identify selection of adequate groups. As indicated in Fig. 2A and Fig. S1, the color-coded heat-map showed that TCGA OC samples were divided into 2 IRGs subgroups performed better than those clustered into more than two categories based on consensus clustering. The color gradients represented consensus values from 0 to 1 (white corresponded to 0 and dark blue to 1). To sum up, the dataset was finally partitioned into two subgroups: cluster 1 (N = 180) and cluster 2 (N = 155). As implied by results of Kaplan–Meier survival analysis, differences in the prognosis were statistically significant across those two de novo clusters (Fig. 2B, P < 0.001). Additionally, the 243 IRGs expression heat-map of all OC patients in the training cohort were shown in Fig. 2C. We could see from heat-map that there was no correlation between two subgroups and clinical characteristics such as age, stage, grade, tumor residue and race.

Construction and assessment of the prognostic model based on IRGs

In order to establish a reliable prognosis signature, batch effects of the TCGA, GETx and ICGC dataset were removed using the “sva” package in R software. Here, we assessed the prognostic value of the 243 differentially expressed IRGs in OC patients using univariate Cox regression analysis (Tables S1 and S2). Afterwards, 25 genes were taken to establish a prognostic risk model for categorizing OC patients into low-risk and high-risk group with discrete OS by LASSO Cox regression analysis. As indicated in Figs. 3A, 3B and S2, the 25 genes were appropriate for building the prognostic risk model for OC patients depending on the primary screen of the IRGs. The risk score was calculated as follows: risk score = 0.0361 × TNFAIP8L3 + 0.0062 × PI3 + 0.0053 × TMEM181 + (− 0.1056 × GBP1P1) + (− 0.0018 × STX18) + 0.0093 × KIF26B + (− 0.0066 × MRPS11) + 0.0135 × CACNA1C + (− 0.0124 × PACSIN3) + (− 0.0027 × GMPR) + (− 0.0040 × MANF) + 0.0003 × PYGB + (− 0.0008 × SNRPA1) + (− 0.0103 × ST7L) + (− 0.0004 × ZBP1) + (− 0.0114 × BMPR1B-DT) + 0.0030 × STAC2 + (− 0.0247 × LINC02585) + (− 0.0252 × LYPD6) + (− 0.0064 × NSG1) + (− 0.0059 × ACOT13) + 0.0037 × FAM120B + (− 0.0061 × LEFTY1) + (− 0.0513 × SULT1A2) + (− 0.0145 × FZD3). Furthermore, all the OC patients were divided into two groups based on the median risk scores calculated by the above formula. Patients in the high-risk group had significantly shorter survival time than the patients in the low-risk group via Kaplan–Meier survival analysis (p < 0.001, Fig. 3C). The ROC curve suggested the predictive ability of the risk score with regard to prognosis; the areas under the curves (AUC) of 1, 2 and 3 years were 0.806, 0.773 and 0.762, respectively (Fig. 3D). Subsequently, the risk curves were used to evaluate the prognostic models by setting the median risk score as the cut-off point (Figs. 3E and 3F). As the risk score rose, the death numbers increased, and the survival time decreased. Univariate and multivariate Cox regression analysis revealed the risk score as an independent prognosis predictor for OS in TCGA cohort (Table 1).

Validation of the prognostic model in ICGC cohort

A group of 85 ICGC OC patients with complete survival information served as an external validation cohort. The formula mentioned above was used to calculate the risk score, and patients were divided into two groups with the same cutoff. As indicated in Fig. 4A, the high-risk group in ICGC cohort also showed poor outcomes as demonstrated by Kaplan–Meier survival analysis (p < 0.05). The AUC obtained from the ROC curve for OS at 1, 2, and 3 years were 0.796, 0.743, and 0.716, respectively (Fig. 4B). Similar to the TCGA cohort, the risk curves were used to show the distribution of risk score and survival status. We could see from Figs. 4C and 4D that the death numbers increased and the survival time decreased as the risk score rose. Univariate and multivariate Cox regression analysis exposes the risk score as an independent prognosis predictor for OS in ICGC cohort (Table 1). These results further confirmed the hardiness of the signature in predicting overall survival of OC.

GO and KEGG pathway enrichment analysis of the different risk groups

To investigate the potential function and mechanisms of the low-risk and high-risk groups, we summarized differential genes in two risk groups to conduct GO and KEGG enrichment analysis in TCGA OC cohort. First, the heat-map of the 25 genes involved in the construction of the prognostic model in the different risk group was shown in Fig. 5A. We could see from the heat-map that the high-risk genes were TNFAIP8L3, PI3, TMEM181, KIF26B, CACNA1C, PYGB, STAC2 and FAM120B.

The top five GO enrichment biological processes were oxidative phosphorylation, mitochondrial gene expression, mitochondrial translation, translational elongation and mitochondrial respiratory chain complex assembly (Fig. 5B, q-value < 0.05). Additionally, KEGG analysis revealed that several immune-related pathways were enriched such as antigen processing and presentation, human T-cell leukemia virus 1 infection, and leukocyte transendothelial migration (Fig. 5C, q-value < 0.05).

Drug sensitivity analysis

CellMiner (http://discover.nci.nih.gov/cellminer) was web-based applications for mining publicly available genomic, molecular, and pharmacological datasets of human cancer cell lines (Reinhold et al., 2012; Reinhold et al., 2015). We obtained the top 16 drugs approved by the US Food and Drug Administration (FDA) with the most statistically significant differences, by performing a separate drug sensitivity analysis on genes involved in the construction of the prognostic model (Fig. 6). The results showed that the expression of NSG1 was positively correlated with the sensitivity of vemurafenib and dabrafenib. The expression of ZBP1 was positively correlated with the sensitivity of LDK-378 and alectinib. The expression of GBP1P1 was positively correlated with the sensitivity of idelalisib and IPI-145. The expression of KIF26B, STX18 and CACNA1C were positively correlated with the sensitivity of zoledronate, nelarabine and idelalisib, respectively. It was indicating that the higher the expression of NSG1, ZBP1, GBP1P1, KIF26B, STX18 and CACNA1C, the stronger the sensitivity to the abovementioned drugs. The expression of GMPR was positively correlated with the sensitivity of vemurafenib and dabrafenib, but it was negatively correlated with the sensitivity of brigatinib. The expression of PYGB was positively correlated with the sensitivity of vemurafenib, but it was negatively correlated with the sensitivity of oxaliplatin. In addition, the higher the expression of PACSIN3 in OC patients, the patient’ s drug resistance to carmustine and bendamustine was stronger. We also predicted the correlation between drug sensitivity of targeted drugs and commonly used chemotherapy agents and patterns of 25 IRGs expression involved in the construction of prognostic models for ovarian cancer (Fig. 6 and Fig. S3). We found that targeted drugs were more significantly associated with gene expression patterns than conventional chemotherapy drugs.

GBP1P1 (LOC400759) was involved in immune cell infiltration in OC

To find biological function for genes involved in building prognostic models, we analyzed the correlation between gene expression and abundance of immune infiltrates by online tool of GEPIA and TIMER. Among all the 25 prognostic genes, the result showed that significant correlation between GBP1P1 expression and immune cell infiltration was the most robust (Fig. 7A and Figs. S4–S7). These findings suggested that GBP1P1 might be involved in the immune cell infiltration in OC.

Cell Localization of GBP1P1

Based on scRNA-seq data of GSE118828, we obtained gene expression profiles from 16 OV samples in different cell types. We accomplished PCA using the top 1,500 variable genes to reduce the dimensionality, and 10 cell clusters were identified. Thereafter, the cell identity of each cluster was annotated using a reference dataset from the Human Primary Cell Atlas, and finally determined five cell types. T cells, endothelial cells, epithelial cells, smooth muscle cells and monocyte, respectively (Fig. 7B). Then, GBP1P1 was distributed on all five cell types, with the highest distribution on epithelial cells (red) as shown in Fig. 7C.

GBP1P1 knockdown inhibits malignant phenotype of ovarian cancer in vitro

We established the ovarian cancer cell line W038 for this study. To evaluate the specific role of GBP1P1 in ovarian cancer, we examined the effects of GBP1P1 knockdown on the proliferation, apoptosis, migration, and invasion of W038 cells in vitro. The W038 cell line was derived from ascites from patients with advanced ovarian cancer. After 50 passages, the cell morphology was homogeneous and cells with mesenchymal morphology were not observed (Fig. 8A). H&E staining showed that W038 cells were adenocarcinoma cells. Immunocytotochemical study revealed that W038 cell line was positive for cytokeratin (CK7), P53 and pair box gene8 protein (PAX-8), but negative for Wilms tumor protein (WT-1), estrogen receptors (ER), progesterone receptors, caudal type homeobox 2 (CDX-2) (Fig. 8B). To determine the efficiency of siRNA depletion, fluorescence images and qRT-PCR were performed. As shown in Fig. 8C, si-GBP1P1 and red fluorescent protein were co-expressed on the carrier, and the transfection efficiency was directly visualized by fluorescence, indicating that siRNA had been transferred into the W038 cell. qRT-PCR analysis for viral transcripts showed GBP1P1 was significantly reduced in sh-GBP1P1 W038 cells compared with sh-CTRL W038 cells (Fig. 8D). The results of CCK8 assays specified that GBP1P1 knockdown significantly inhibited the proliferation rate of W038 cells (Fig. 8E). The results of transwell assay demonstrated that the W038 cells exhibited markedly decreased invasion upon GBP1P1 knockdown (Figs. 8F and 8G). The results of wound healing assay revealed that the knockdown of GBP1P1 significantly decreased the migration of W038 cells (Figs. 8H and 8I). The apoptosis was observed by flow cytometry (IF488-Annexin V/7-AAD) and the results showed that the apoptosis rate of W038 cells increased significantly in the sh-GBP1P1 groups (Figs. 8J and 8K).