Crosstalk Between Immunity and Oncogenes Within the Tumor Microenviron

Introduction

Cervical cancer remains a major global health concern and is the second most common gynecologic malignancy, particularly among women with persistent high-risk human papillomavirus (HPV) infections.¹ Approximately 93.8% of CSCCs are linked to high-risk human papillomavirus (hrHPV) infections.^2–4 Globally, an estimated 604,127 new cases of cervical cancer are diagnosed each year, resulting in over 350,000 deaths.^3,5 Globally, nearly 291 million women are carriers of hrHPV DNA (Cervical cancer (who.int). Persistent hrHPV infection can disrupt the local immune microenvironment, creating a proinflammatory and immunosuppressive tumor microenvironment (TME) that promotes CSCC development and alters therapeutic outcomes.⁶ The interplay between oncogenes, tumor suppressor genes, and tumor-promoting inflammation is crucial in the etiology and progression of cancer. Inflammation serves as a double-edged sword, awakening a favorable tumor microenvironment (TME) through angiogenesis and immune surveillance inhibition while perpetuating cancer initiation via processes like DNA damage and inducing genetic mutations. Cancers can also evoke chronic inflammation and an immunosuppressive microenvironment, supporting tumor progression and metastasis.⁷

This inflammatory milieu releases cytokines, chemokines, and growth factors from stromal and immune cells within the TME, fueling tumor-promoting pathways.^7,8 The TME is a highly complex and heterogeneous network composed of cancer-associated fibroblasts (CAFs), diverse immune cell populations (e.g., T cells, macrophages, dendritic cells), and non-cellular components such as the extracellular matrix.⁹ Previous data revealed the importance of immune and tumor-related genes in CSCC, particularly in their prognostic impact and therapeutic potential. For instance, the expression profiles of genes associated with immune function have been linked to improved clinical outcomes, indicating their potential as important biomarkers for predicting prognosis and survival.² The dynamic interactions and continuous crosstalk between tumor cells and their surrounding microenvironment can significantly influence tumor behavior, either fostering or inhibiting growth, progression, and metastasis, as well as the production of inflammatory signals, the attraction of immune cells, and the regulation of gene activity.^10–13

Targeted transcriptomic profiling using RT² Profiler PCR arrays offers a robust, scalable, and clinically relevant alternative for interrogating immune-oncogene interactions in cancer patient cohorts. This method allows for the focused quantification of immune response genes, oncogenes, and tumor suppressors with high sensitivity and reproducibility. Spatial transcriptomics and single-cell RNA sequencing (scRNA-seq) have significantly advanced our understanding of the tumor microenvironment by offering high-resolution spatial and cellular context.^14,15 These technologies are currently limited by cost, complexity, and clinical scalability. In contrast, targeted PCR-array platforms remain a practical and robust tool for translational research, enabling high-throughput quantification of clinically relevant genes across well-characterized patient cohorts. Our study leverages this approach to generate actionable insights into immune-oncogene dynamics in HPV-positive CSCC. There remains an incomplete understanding of the HPV-associated CSCC immune microenvironment, particularly concerning the interactions between inflammation, immunity, oncogenes, and tumor suppressor genes, as well as disease progression.

The objective of this study was to investigate the transcriptomic landscape of the HPV-associated CSCC tumor microenvironment by examining the expression of 168 functionally annotated genes involved in inflammation, immune regulation, signal transduction, transcriptional control, oncogenesis, tumor suppression, angiogenesis, and apoptosis. Using RT² profiler PCR profiling, immune infiltration assessment, survival correlation, pathway enrichment, and interaction network analysis, we aim to identify prognostic biomarkers and uncover potential therapeutic targets relevant to the immune-oncogenic axis in HPV-associated CSCC.

Materials and Methods

Biological Specimen Collection and Processing and Storage

Formalin-fixed paraffin-embedded (FFPE) tissues were obtained from CSCC and nonmalignant cervical tissue (37 CSCC and 13 nonmalignant tissues). The FFPE blocks were stored in a dry place without humidity at 15°C and protected from light. This study was conducted at the Research Center, King Fahad Medical City. The CSCC cases were confirmed by two independent pathologists. From each paraffin block, a series of 10-μm-thick sections were prepared and stored at 4°C until used. We included participants who have: (1) Confirmed cases of CSCC; (2) aged 18 years or older; (3) complete patient clinical and demographic information available. Exclusion criteria included: (1) Patients who had received chemotherapy, immunotherapy, or radiation treatment; (2) Patients with a history of other malignancies. All tissue samples were screened for HPV using quantitative real-time PCR (qRT-PCR).

Nucleic Acid Extraction and First-Strand cDNA Synthesis

RNA was extracted from FFPE samples using the Qiagen AllPrep FFPE Kit (Hilden, Germany) in accordance with the provided protocol. Subsequent reverse transcription and genomic DNA removal were performed using the RT² First Strand Kit (Qiagen, Germantown, MD, USA) to generate high-quality complementary DNA (cDNA) for downstream applications. To assess quality and integrity of the extracted nucleic acids RNA, all samples were evaluated using a NanoDrop 2000 (Thermo Scientific) for concentration and purity, and two internal controls hydroxymethylbilane synthase (HMBS) and glyceraldehyde 3-phosphate dehydrogenase (GAPDH) were used to evaluate the quality and integrity of the extracted nucleic acids. To avoid contamination, rigorous protocols were followed during RNA handling, providing RNase-free conditions and working in a biological safety cabinet.

Gene Expression Profiling of Immunity and Tumor-Related Genes

mRNA expression was analyzed through a series of standardized procedures. Complementary DNA (cDNA) was synthesized, and genomic DNA was removed using the RT² First Strand Synthesis Kit (Qiagen, Germantown, MD, USA), in accordance with the manufacturer’s guidelines. The cDNA was analyzed by RT2-PCR using the two panels Cancer Inflammation and Immunity Crosstalk Array (PAHS-181Z) and the Oncogenes and Tumor Suppressor Genes Array (PAHS-502ZC-12) from SABiosciences (Frederick, MD, USA). Expression levels were evaluated for a total of 168 genes involved in cancer cell-cell interaction, inflammatory mediators, oncogenesis, tumor suppression, transcription factors, angiogenesis, and apoptosis. A complete listing of RT2 Profiler PCR Array Kits genes is available at https://geneglobe.qiagen.com/us/product-groups/rt2-profiler-pcr-arrays). The qPCR was executed using the RT2 SYBR Green/ROX qPCR Master Mix (Qiagen), as described previously.^4,16 The final volume for each 96-well plate was 2700 µL, which contained 1350 µL of 2× RT2 SYBR Green Master Mix, 102 µL of cDNA and 1248 µL of RNase-free water. 25 µL of the qPCR component mix was dispensed into each well of the PCR array. The 96-well array was subsequently introduced into the ABI 7500 Fast instrument (Applied Biosystems, Waltham, MA, USA) for amplification. The PCR cycling protocol consisted of one cycle of 10 minutes of denaturation at 95 °C, followed by 40 cycles of 15 seconds at 95 °C for denaturation and 1 minute at 60 °C for annealing, during which extension and fluorescence data were also acquired. Relative gene expression data were analyzed by Qiagen’s online RT2-PCR Profiler data analysis software (www. qiagen.com). The relative expression was quantified by the ΔΔCt method and normalized using five housekeeping genes.

Analysis of Real-Time RT² Profiler PCR Array

PCR Array data was examined using the Qiagen Online Analysis Software (https://geneglobe.qiagen.com/us/analyze). Each plate included 3 positive PCR controls and 3 RT controls to ensure assay reliability. Contamination levels of human genomic DNA were assessed following the manufacturer’s guidelines. The cycle threshold (Ct) values acquired during the reaction were utilized to calculate fold changes in mRNA expression, employing the ²⁻ΔΔCt method. To enhance accuracy, ΔCt values were normalized against the average of 5 housekeeping genes.

Bioinformatics Analysis of Experimental-Derived Data

Expression and Patient Prognosis Analysis

The correlation and clinical prognostic value of the expression levels of the 10 most co-upregulated genes in CSCC, CASP8, ZHX2, BCL2L1, RB1, BAX, CD274, CCL20, CSF2, FOXP3, and CCL18 were analyzed with OS and DFS. The Kaplan–Meier estimator was employed to evaluate OS and DFS in CSCC patients, using the median (50%) expression level as the cut-off to define high and low expression groups. Statistical significance between these groups was determined using Log rank tests, with a p-value of < 0.05 considered significant.

This analysis utilized the GEPIA2 (Gene Expression Profiling Interactive Analysis) database (http://gepia2.cancer-pku.cn/#index u.cn), which provides comprehensive data sourced from The Cancer Genome Atlas (TCGA) and Genotype-Tissue Expression (GTEx).

Immune Cell Infiltration in the Tumor Microenvironment

The Tumor IMmune Estimation Resource (TIMER) assessed immune cell infiltration within tumor tissues. The “Immune” module under the “Gene” section of the TIMER platform (https://cistrome.shinyapps.io/timer/) was utilized to examine the association between immune infiltration and the 10 expressed genes in CSCC. A p-value of less than 0.05 was considered statistically significant. Cancer-associated fibroblasts (CAFs) and regulatory T cells (Tregs) were also selected for infiltration assessment. Immune infiltration levels were estimated using the EPIC, MCPCOUNTER, XCELL, and TIDE algorithms, which facilitate the quantification of various immune cell types within the TME. The gene expression data and corresponding infiltration levels were visualized using scatter plots to illustrate the relationships. Purity and immune infiltration p-values were calculated using the Spearman’s rank correlation test, adjusted for tumor purity. The most upregulated genes in various cancer types were performed by TIMER Analysis. We further investigated the immune microenvironment of CESC and the associated immune cell infiltrates via GSCA: Gene Set Cancer Analysis (wchscu.cn) database. We analyzed the correlation between immune cell infiltrates and the Gene Set Variation Analysis (GSVA) enrichment scores.

Correlation and Multiple Gene Analysis

Differential gene expression analyses across tumor, normal, and metastatic tissues as well as multiple gene analysis, correlation matrix generation, and gene signature evaluation were performed using the TNMplot tool (https://tnmplot.com/). TNMplot offers comprehensive in silico analysis by integrating data from TCGA, GEO, and GTEx to compare gene expression across tissue types.

Pathways, Gene Ontology and Enrichment Analysis

Metascape (http://metascape.org) tool was used for gene function enrichment and classification, facilitating the analysis of gene sets. The Gene Ontology (GO) framework encompasses three primary components: biological processes, cellular components, and molecular functions, which collectively provide a comprehensive annotation of gene roles and their interrelationships. Genes were analyzed using the GO function and KEGG pathway enrichment network available on the Metascape platform. A minimum overlap threshold of 1.5 and a minimum enrichment factor of 1 were established, with results filtered by a significance cutoff of P < 0.05. The Molecular Complex Detection (MCODE) algorithm was employed to identify interconnected hub genes within the networks. Each identified MCODE network was color-coded to visually represent the strength of interactions between the molecules. For the interaction enrichment analysis, default parameters for minimum and maximum network sizes were set at 2 and 500, respectively. To refine the results, only those gene sets with a false discovery rate (FDR) of less than 0.05 were included in the final output. The Kyoto Encyclopedia of Genes and Genomes (KEGG) served as a robust database for pathway annotation, offering insights into biological pathways.

Pathway Enrichment Analysis and Data Visualization

Pathway enrichment analysis was employed separately to obtain a complete picture of 45 genes that exhibited relatively similar expression levels and were selected for further investigation using the g: Profiler platform. This analysis focused on the ontology biological process (GO BP), gene ontology molecular function (GO MF), gene ontology cellular compartment (GO CC), and KEGG pathways.

Functional Analysis of Gene Expression

NetworkAnalyst (http://www.networkanalyst.ca) is an extensive online platform that facilitates both basic and advanced meta-analyses of gene expression data through a user-friendly interface. The NetworkAnalyst platform simplifies the complex molecular landscapes of the gene expression data. The ExpressAnalyst module facilitates the visualization, statistical analysis, and functional interpretation of gene expression data uploaded in list or table formats. It features three distinct modules designed for statistical analysis-based data input. The network-based representation enables the identification of key signaling hubs, crosstalk between pathways, and potential functional relationships that may underlie various biological and pathological states.

Drug Sensitivity Analysis Using Genomics of Drug Sensitivity in Cancer (GDSC) and Cancer Therapeutics Response Portal (CTRP)

Gene Set Cancer Analysis (GSCA) (https://guolab.wchscu.cn/GSCA/#/) is a web server that integrates genomics and immunomics data derived from the TCGA database. The relationship between the expression of the 10 upregulated genes and drug sensitivity was examined using data from the Genomics of Drug Sensitivity in Cancer (GDSC) and the Cancer Therapeutics Response Portal (CTRP). This analysis was conducted across a range of cancer types, with a particular focus on the 30 most relevant types.

Physical Interaction Network of Host-Virus Interactions

The Enrichr (http://amp.pharm.mssm.edu/Enrichr) platform was used for the biological properties analysis. Enrichr is an open-source web server and contains about 100 gene set libraries with more than 180,000 gene sets. The study focused on “Diseases/Drugs” and used VirusMINT to explore overexpressed genes and their physical interaction with viral proteins. Viruses can significantly alter cellular regulatory mechanisms, affecting critical processes such as gene expression, cell growth, and differentiation at both the transcriptomic and proteomic levels. To elucidate the molecular mechanisms that underlie this disruption of cellular physiology, it is essential to investigate how viral proteins perturb these networks.

Statistical Analysis

Data were analyzed using GraphPad Prism version 10.0 (GraphPad Software, San Diego, CA, USA). For comparisons between two groups, Student’s t-test was used. For analyses involving more than two groups, one-way or two-way ANOVA was performed, followed by Sidak’s or Tukey’s post hoc tests for multiple comparisons, as appropriate. Correlation analyses were conducted using Pearson or Spearman correlation coefficients, and a correlation matrix was generated to visualize gene-to-gene relationships. For RT² Profiler PCR array data, the ΔΔCt method was applied using the Qiagen GeneGlobe Data Analysis Center. All primer assays had a PCR efficiency greater than 90%. Each array included quality control wells to assess genomic DNA contamination, reverse transcription efficiency, and overall PCR performance. Data normalization was performed using the geometric mean of at least six validated housekeeping genes. To determine differential gene expression, genes with greater than a 1.5-fold change and a raw p-value < 0.05 were initially considered significant. To account for multiple hypothesis testing across the high-throughput dataset, the false discovery rate (FDR) correction was used to adjust p-values in the context of differential gene expression, correlation matrices, and pathway enrichment analyses. Kaplan–Meier survival curves were generated to assess associations between gene expression and overall or disease-free survival. The log-rank (Mantel–Cox) test was used to compare survival distributions. Only adjusted p-values (q-values) < 0.05 were considered statistically significant.

Results

Analysis of Gene Expression Patterns of Tumor Inflammation, Immunity Crosstalk, Oncogenes and Tumor Suppressor

The interaction between cancer, inflammation, and the regulation of oncogenes and tumor suppressor genes plays a critical role in the initiation and progression of cancer. Inflammation has been shown to initiate cancer and promote angiogenesis and hinder immune surveillance. To explore the relationships between inflammation and gene regulation in cancer, the mRNA expression profiles of 168 immune, oncogene, and tumor suppressor-related genes within the TME of CSCC associated with HPV were analyzed. Our analysis revealed that 94 genes were differentially expressed in HPV-positive CSCC, with a fold change greater than 1.5. Notably, 55 of these genes were upregulated (Table 1). Of the upregulated genes, we identified ten that exhibited particularly high expression levels: CASP8, ZHX2, BCL2L1, CTNNB1, RB1, BAX, CD274, CCL20, FOXP3, and CCL18. CASP8 demonstrated the most significant increase, with a fold change of 22.47, followed by ZHX2 and BCL2L1, which showed substantial upregulations of 17.57 and 16.70, respectively. Other noteworthy genes included RB1 (fold change 12.12), BAX (fold change 9.33), and CCL20 (fold change 7.08). Additionally, genes such as CTNNB1 (fold change 8.43), CXCL18 (fold change 5.10), and FOXP3 (fold change 5.22) were also significantly expressed (Table 1). These findings showed a relationship between inflammation and oncogenic pathways in HPV-positive CSCC, highlighting potential targets for therapeutic intervention.

Table 1 Up-Regulated Gene Expression

Downregulation Patterns of Tumor Inflammation, Immunity Crosstalk, Oncogenes, and Tumor Suppressor

The transcriptome analysis showed that there was a significant downregulation of 39 genes associated with tumor biology, where 13 of the 39 genes were downexpressed in the study group when compared with the control group (Table 2). HRAS also had the most significant downregulation with a fold change of −30.18, suggesting that cell signaling whose modulation is responsible for controlling growth may be interrupted. In the same way, CCND1 (fold change −28.92) and ATM (fold change −17.90) indicated important effects on cell cycle regulation and DNA damage response, respectively. Notably, other pivotal genes like RUNX3 (fold change −17.11), E2F1 (fold change −15.32), CXCR1 (fold change −14.16) and MIF (fold change −12.45) were marked downregulated (Table 2).

Table 2 Downregulated Gene Expression

Gene-Gene Interactions and Correlation Analysis of Gene Expression

A correlation analysis was conducted to examine the relationships between various genes (Figure 1). The correlation analysis revealed several strong positive and negative correlations. There was a direct correlation between RB1 and CASP8 (r = 0.56). CCL20 showed a positive correlation with CCL18 (r=0.55) and BAX (r=0.39). CCL18 and CCL20 (r = 0.41) exhibited a positive correlation. BCL2L1 and FOXP3 (r = 0.49) showed a positive correlation. Moderate to weak positive association was identified between CSF2 and BAX (r = 0.22), CSF2 and BCL2L1 (r = 0.25), CSF2 and FOXP3 (r = 0.26), CASP8 and ZHX2 (r = 0.24) and CASP8 and BAX (r = 0.24). In contrast, RB1 and BCL2L1 showed negative correlations (r = −0.48) and RB1 and FOXP3 (r = −0.35) revealed a negative correlation. Other negative correlations included CASP8 with BCL2L1 (r = −0.28), CTNNB1 with FOXP3 (r= −0.42), and FOXP3 (r = −0.34) (Figure 1).

Figure 1 Heatmap Gene-Gene Interactions and Correlation Analysis of Gene Expression. Correlations are color-coded, with red indicating positive and blue indicating negative associations. The intensity of the color corresponds to the magnitude of the correlation coefficient. Each cell in the matrix displays the correlation coefficient between the corresponding gene pairs.

Gene Expression Levels in Tumor, Normal and Metastatic Cervical Tissues

To elucidate the potential biological roles of the most highly expressed genes in CSCC, gene expression data from the TNMplot platform was utilized. A differential expression analysis of these genes was conducted across tumor, normal, and metastatic cervical tissues. The results demonstrated that cervical and metastatic tumor samples exhibited significantly higher expression levels than normal samples (Figure 2). Notably, CASP8 and BAX expression levels were particularly elevated in metastatic cervical tumors (Figure 2). The upregulation of CD274 can promote tumors from immune attack.

Figure 2 Boxplots of the Most Highly Upregulated Genes, illustrating the expression levels in tumor, normal and metastatic cervical tissues. The Kruskal. Wallis’ p-value was utilized to assess statistical significance. The bars illustrate the proportion of tumor and metastatic tissue samples with elevated gene expression compared to normal tissues, evaluated at multiple quantile thresholds (minimum, 1st quartile, median, 3rd quartile, and maximum).

Highly Expressed Genes are Associated with Both the Prognosis and Long-Term Survival of Patients with CSCC

The study then aimed to evaluate the prognostic impact of upregulated genes. The prognostic significance of the ten most upregulated genes in CSCC was assessed using the GEPIA2 tool. Of these 10 genes, there were 6 genes (CASP8, BAX, RB1, CD274, FOXP3, and CCL18) whose expressions were positively correlated with the lengths of the OS and DFS (Figure 3). Patients with CASP8, BAX, RB1, CD274, FOXP3, and CCL18 overexpression had significantly longer survival and better prognosis than patients with low expression. The expressions of the remaining 4 genes exhibited similar correlations with the lengths of the OS and PFS (Figure 3). These findings suggest that the upregulation of these genes correlates with a better prognosis, highlighting their potential role as prognostic biomarkers.

Figure 3 The prognostic value of highly upregulated genes. (A-J) Kaplan-Meier curves of survival outcomes for each gene by GEPIA2. The p-values were calculated using the Mantel-Cox test to assess statistical significance. Hazard ratios (HR) were calculated from the Cox proportional hazards model. HR >1 indicates that high gene expression is associated with a reduced survival rate, whereas HR < 1 suggests that high gene expression correlates with an increased survival rate.

Cancer-Related Pathways Activity

Gene Set Cancer Analysis (GSCA) was employed to evaluate the potential cancer pathway activities. The GSVA results indicated that cancer-related pathways associated with different activities were significantly linked to the most highly upregulated genes. Gene expression showed a positive association with the activation of key biological processes, including epithelial-mesenchymal transition (EMT), DNA damage response, apoptosis, cell cycle regulation, and androgen receptor signaling (hormone AR), PI3KAKT, RASMAPK, RTK and TSCmTOR pathways (Figure 4). This suggests that RB1 expression is positively associated with apoptosis but negatively associated with pathways like mTOR and RTK, which are critical for cell growth and survival.

Figure 4 Gene set variation analysis (GSVA) of the most upregulated gene. The heatmap presents the relationship between GSVA scores and the activity of cancer-associated pathways in CESC. The GSVA score reflects the overall expression level of a gene set and is directly correlated with individual gene expression levels. *P value ≤ 0.05; #FDR ≤ 0.05.

Pathway and Process Enrichment Analysis

The Metascape database was used to analyze and visualize gene enrichment across Gene Ontology (GO) terms, including “biological processes” “cellular components” and ‘molecular functions. The pathway and process enrichment analysis revealed deeper understanding of the biological functions associated with the 10 upregulated genes, with Gene Ontology (GO) analysis identifying enrichment in 13 distinct biological processes, including multicellular organismal process, biological processes involved in interspecies interaction, locomotion, reproductive process, homeostatic process, immune system process, developmental process, response to stimulus, regulation of biological processes, both downregulation and Upregulation of biological processes, cellular processes, and biological regulation (Figure 5A). The network of enriched terms involved small-cell lung cancer, regulation of mononuclear cell proliferation, head and neck squamous cell carcinoma, negative regulation of cell population proliferation, camera-type eye morphogenesis and elevated expression of cell migration (Figure 5B).

Figure 5 Top 6 Enriched Ontology Clusters with their representative enriched terms analyzed by Metascape. (A) Enriched ontology clusters were color-coded based on their cluster IDs, with nodes belonging to the same cluster generally positioned in close proximity. (B) The thickness of the edge represents the similarity score.

Enrichment Analysis of Biological Pathways

Analysis of gene lists with similar expression levels was conducted to identify enriched Gene Ontology (GO) categories and biological pathways. The results were classified into distinct functional categories based on GO term enrichment for biological processes (BP), molecular functions (MF), cellular components (CC), and KEGG pathways (Figure 6). The GO BP analysis revealed significant enrichments in immune system processes, intracellular signaling cascades, positive regulation of developmental processes, inflammatory responses, cellular responses to cytokine stimuli, cell migration, and locomotion (Figure 6A). Regarding molecular functions, the most significantly enriched activities included C-C chemokine receptor activity, G protein-coupled chemoattractant receptor activity, chemokine receptor activity, cytokine binding, and chemokine binding (Figure 6B). According to the KEGG database, the most significantly enriched pathways included cytokine-cytokine receptor interactions, viral protein interactions with cytokines and cytokine receptors, pathways in cancer, and the MAPK signaling pathway (Figure 6C). The CC enriched categories were primarily associated with the external side of the plasma membrane, cell surface, membrane side, interleukin-23 complex, plasma membrane, and cell periphery. REACTOME pathway analysis revealed that the most significantly enriched pathways involved chemokine receptor interactions, peptide ligand-binding receptors, IL-4 and IL-13 signaling, as well as multiple other interleukin- and chemokine-related pathways. The GO was enriched in several pathways: inflammatory response, interleukin-4 and interleukin-13 signaling, signaling by interleukins, cancer pathways, PID telomerase pathway and cellular response to cytokine stimulus G.

Figure 6  Functional enrichment analysis of 21 with similar expression levels in CSCC using g: Profiler. Graphs (A–C) showed the results of the enriched terms with the statistical significance (p-value) according to the (A) gene ontology biological processes (GO BP), (B) gene ontology molecular functions (GO MF), and (C) KEGG.

Interaction Network Analysis

To investigate the functional relationships among upregulated genes identified in HPV-associated CSCC, we constructed a protein–protein interaction (PPI) network using NetworkAnalyst. The PPI network revealed both direct and indirect interactions between upregulated genes and their broader regulatory context as well as shared biological pathways. The network encompassed a wide range of biological entities, including proteins, protein families, complexes, chemicals, small molecules, microRNAs (miRNAs), phenotypes, and stimuli (Figure 7). A total of 859 nodes and 1,196 edges were identified, where nodes represent proteins or associated biological entities, and edges indicate known or predicted interactions between them (Figure 7). The analysis revealed several hub proteins with high connectivity, suggesting their potential central roles in HPV-associated CSCC biological processes. Notably, the presence of regulatory elements such as miRNAs and chemical stimuli in the network highlights the multifaceted control of upregulated genes signaling in HPV-associated CSCC.

Figure 7 Topology network analysis illustrating protein-protein interactions (PPI), networks encompassing various biological entities such as chemical, complex, miRNA, phenotype, protein, protein family, small molecule, and stimulus using NetworkAnalyst. The central panel displays the network with distinct modules highlighted in various colors. MicroRNAs are represented by blue squares, while genes are depicted as red circles. The size of each circle corresponds to the degree of the node within the network. This analysis comprises a total of 859 nodes and 1196 edges. Nodes in the network represent proteins and biological entities, while edges represent the interactions between them.

Analysis of Drug Sensitivity

Patient sensitivity to chemotherapy and targeted therapies can be significantly influenced by genomic abnormalities. The GDSC and CTRP database allows researchers to evaluate how cancer types respond to various therapeutic agents. This analysis is important to analyze the effectiveness of drugs in treating specific cancers. This study investigated the role of the ten most highly co-upregulated genes in CSCC in mediating patient responses to cancer treatments. The expression levels of CASP8, ZHX2, BCL2L1, RB1, BAX, CD274, CCL20, CSF2, FOXP3, and CCL18 were positively correlated with drug sensitivity to several compounds (Figure 8). This positive correlation suggests that higher expression of these genes may contribute to drug resistance. The expression of RB1 and BAX was found to be negatively associated with drug sensitivity to some compounds, indicating that increased expression of these genes may enhance the sensitivity of the drug (Figure 8).

Figure 8 Correlation analysis between (A) GDSC and (B) CTRP drug sensitivity datasets and the mRNA levels of the ten most highly co-upregulated genes across various cancers. Spearman correlation indicates the relationship between gene expression and drug response. Positive correlation suggests that elevated gene expression is associated with increased drug resistance, whereas a negative correlation implies that higher expression corresponds to greater drug sensitivity. The figure summarizes the correlation between the ten most co-upregulated gene expression levels and the sensitivity of the top 50 drugs in pan-cancer.

Tumor Microenvironment and Immune-Infiltration Analysis

A comprehensive analysis was conducted using the TIMER database to examine the immunological landscape associated with the ten most highly upregulated genes in cervical squamous cell carcinoma (CSCC). Immune infiltration analysis demonstrated associations between the expression levels of these genes and various immune cell types, including B cells, CD8⁺ T cells, CD4⁺ T cells, macrophages, neutrophils, and dendritic cells (Figure S1). Correlation analysis revealed both positive and negative relationships between specific genes and the infiltration levels of immune cells such as B cells, macrophages, and dendritic cells (Figure S2).

Further analysis using the TIMER-derived heatmap demonstrated Spearman correlation coefficients between gene expression and cancer-associated fibroblasts (CAFs), highlighting a positive correlation between certain genes and CAFs. The EPIC, MCPCOUNTER, XCELL, and TIDE algorithms were applied to assess immune infiltration across multiple TME components. A Spearman’s Rho value below 0.2 suggested no significant correlation between tumor purity and the expression levels of the top ten genes.

Gene Set Variation Analysis (GSVA) was performed to evaluate immune cell enrichment scores. The resulting heatmap revealed that the expression of the ten most upregulated genes was positively associated with a broad range of immune cell types, including B cells, CD8-naive cells, monocytes, macrophages, dendritic cells, natural regulatory T cells (nTregs), induced regulatory T cells (iTregs), natural killer (NK) cells, NKT cells, γδ T cells, exhausted T cells, and CD4⁺ T cells (Figure S3). In contrast, several genes showed negative associations with specific immune cell infiltrates (Figure S3).

Signaling Pathway Analysis

To understand the underlying biological mechanisms and pathways of all upregulated genes, EpressAnalyst was used to conduct a comprehensive gene expression analysis, visualizing the enrichment of functional Gene Ontology (GO) biological processes (BP) and KEGG pathway networks. The gene set network alongside the functional categories enriched in the GO analysis showed various functional categories that are significantly represented within the network (Figure S4A). Similarly, KEGG pathway enrichment analysis showed the relationships between different pathways (Figure S4B). The connections between the nodes showed the complex interactions and relationships between these biological entities. This comprehensive network map provides a general view of the intricate interplay of diverse signaling pathways and their associations with a wide range of diseases, conditions, and cellular phenomena.

Interaction Network Analysis

Using the GeneMANIA software, we constructed an extensive protein-protein interaction (PPI) network. This analysis allowed us to visualize the genetic interactions and pathways linking these genes (Figure S5). The PPI network revealed both direct and indirect interactions, including physical associations, co-expression relationships, and shared biological pathways. The analysis highlighted several co-expressive genes that exhibited strong connectivity with multiple hub genes (Figure S5).

Physical Interaction Network of Host-Virus Interactions

Viruses can significantly alter cellular regulatory mechanisms, affecting critical processes such as gene expression, cell growth, and differentiation at both the transcriptomic and proteomic levels. To elucidate the molecular mechanisms underlying this disruption of cellular physiology, it is essential to investigate how viral proteins perturb these networks. The results of this study indicated that upregulated genes interact with proteins encoded by several clinically significant viruses, including HPV-1a, HPV-11, HPV-16, HPV-18, HPV-31, and HPV-6b, as well as human adenoviruses (types 2, 5, and 12), hepatitis C virus, and vaccinia virus strains (Figure S6). The enrichment biological properties analysis enables us to identify the most interactions with clinically significant viruses where our input genes interacted with and over-represented, playing a crucial role.

Discussion

The Tumor Microenvironment (TME) is critical in initiating and progressing various cancer types, leading to it becoming a prominent research focus in recent years. Understanding immune-oncogenes and tumor suppressor genes in the TME is still highly challenging.¹⁶

To study gene expression in the specific TME of HPV-associated conditions, the expression of 168 genes involved in tumor cell interaction with cellular mediators of inflammation and immunity, oncogenesis, tumor suppression, transcription regulation, angiogenesis and apoptosis were profiled. Bioinformatics methods and clinical data were combined. The analysis identified 94 genes with differential expression in HPV-positive CSCC. Importantly, 55 of these genes were upregulated. These included the high expression signatures of CASP 8, ZHX 2, BCL 2 L 1, CTNNB 1, RB 1, BAX, CD 274, CCL 20, FOXP 3 and CCL 18 of the complete upregulated genes. Patients with higher expression levels of CASP8, BAX, RB1, CD274, FOXP3, and CCL18 survived longer than those with lower expression levels and had a better prognosis. The study found crosstalk of inflammation, oncogenic pathways, and immune response in HPV- positive CSCC. The upregulated expression of apoptosis regulator (CASP8) and pro-apoptotic protein (BAX) may play a role in increased apoptosis, inhibition of tumor growth, and repression of tumor cell growth.¹⁷ In contrast, the downregulation of HRAS, CCND1, and ATM suggested dysregulation of DNA repair pathways and cell cycle regulation. This dysregulation could lead to genomic instability and possibly even tumorigenesis and subsequent malignant transformation. The positive correlations between immune and tumor suppressor-related genes identified in the analysis indicated the interplay of gene-gene interactions. The prognostic significance of our findings is considerable. Common among the most upregulated genes were CASP8, CCL18, FOXP3, and BAX. These results provide exciting opportunities for therapy and biomarker development with improved patient survival.

The tumor and immune-related genes identified in this study align with those previously observed in studies that describe the critical role that the immune microenvironment plays in HPV-positive cancers.⁴ High expression of CASP8 has demonstrated a good prognostic value for gastric cancer patients, which is associated with prolonged survival.^18,19

CASP8 acts as an important mediator of apoptosis in the cell death pathway, responsible for eliminating damaged or cancer cells. Elevated CASP8 expression indicates that tumors are more susceptible to apoptogenic signals and, therefore, undergo cell death more frequently with the decrease in tumor mass.^18,19 Furthermore, patients with elevated CASP8 expression may benefit from improved treatment outcomes, as therapies that induce apoptosis can be more effective.²⁰ Thus, assessing CASP8 expression could serve as a valuable biomarker for predicting patient outcomes and tailoring therapeutic strategies in CSCC, highlighting its potential role in improving prognosis.²¹ High mRNA expression levels of CASP8 were found in this study. Both CASP8 and BAX are included in apoptotic pathways. Their co-overexpression may suggest an increased tumor apoptosis activity, which could be associated with better responses to therapy and longer survival. Accumulated data have shown that different upregulated genes can be utilized as prognostic cancer biomarkers or therapy prognosis response, supporting the implications of the findings in this study.^17,22,23

FOXP3, a member of the FOX protein family, is a transcription factor originally identified in regulatory T cells, where it plays a crucial role in their differentiation. Growing evidence indicates that FOXP3 is also expressed across a range of cancer types.²⁴ Several studies have revealed the upregulation of FOXP3 in many cancers, such as tonsillar squamous cell carcinoma, thyroid, breast, lung and gastric cancers.^16,25

FOXP3 plays a dual and broader role in tumor biology. However, many studies have shown that the role and expression of FOXP3 in cancer is inconsistent and even the opposite. For example, the overexpression of FOXP3 in hepatocellular carcinoma suppresses tumor progression, emphasizing FOXP3’s potential as both a prognostic biomarker and a novel therapeutic target.²⁶ In breast cancer, FOXP3 was shown to be a marker for longer survival and treatment response.^27,28 FOXP3 expression is upregulated in certain cancers and has been implicated in promoting tumor progression, often correlating with poor prognosis. In non-small cell lung cancer, the overexpression of FOXP3 has been linked with tumor growth and metastasis.²⁹ The exact reason for the contradictory results is undetermined. The biological significance of these results warrants further in vitro and in vivo studies.

In this study, high levels of CD274 expression were correlated with better survival rates. CD274 plays critical roles in tumor development and immune response, making it a useful biomarker for predicting patient prognosis and understanding the cancer immune landscape.^25,30

CD274 (PD-L1) plays a well-characterized role in immune checkpoint regulation, inhibiting T-cell activation when engaged with PD-1. While high PD-L1 expression is often associated with immune evasion and poor outcomes, it can also indicate a T cell-inflamed TME that is more responsive to checkpoint blockade therapy.³¹ The association of CD274 with improved survival in our cohort may reflect an “immune-hot” tumor state, where PD-L1 is upregulated as part of an active anti-tumor immune response. These findings emphasize the importance of cellular context, tumor-immune dynamics, and HPV-related immunomodulation in interpreting gene function. Active immune response was associated with high levels of PD-L1 expression. The overexpression of PD-L1 in cancer patients may benefit from immunotherapies, leading to better clinical outcomes and longer survival. Its expression in various tumors has been a biomarker of prognosis in many cancer types.^32,33 Many studies also report upregulation of BCL2L1 in HPV-positive cancers, suggesting common pathways involved in tumor progression and immune modulation and chemotherapy or radiation sensitivity.^34,35 Another study showed significantly higher BCL2 and lower p53 expression in HPV-positive tumors than HPV-negative tumors, suggesting a possible role of HR-HPV in cancer carcinogenesis, in which BCL2 and p53 are involved.^36,37 These genes show a significant association with improved OS in CSCC patients, with higher expression levels correlating with favorable clinical outcomes, highlighting their potential as prognostic biomarkers.³⁸ Given their involvement in apoptosis, immune regulation, and cell cycle control, these genes represent promising candidates for inclusion in multiplex diagnostic or prognostic panels, especially in the context of immune checkpoint blockade and other immunotherapy strategies. Their expression profiles, when validated in larger cohorts, may aid in risk stratification, inform personalized treatment decisions, and potentially serve as companion biomarkers for patient selection in clinical trials targeting the tumor microenvironment.

Targeting the identified pathways in HPV-positive CSCC presents several challenges. One major challenge is that the heterogeneity of tumors can result in the variable expression of these targets, diminishing the efficacy of treatments. Additionally, the complexity of signaling networks that disrupt one pathway may activate anti-pathways, potentially leading to treatment failure. Lastly, there is a lack of comprehensive biomarkers to predict patient response to these targeted therapies, which complicates patient selection and treatment personalization.

Combining therapies that target different pathways may enhance treatment efficacy and overcome resistance. For instance, combining a BCL2L1 inhibitor with a CASP8 activator could induce robust apoptosis in cancer cells.³⁹ Similarly, pairing a ZHX2 modulator with immune checkpoint inhibitors may not only inhibit tumor progression but also enhance anti-tumor immunity.^40,41 Additionally, combining RB1 restoration strategies with CTNNB1 inhibitors could synergistically prevent tumor growth by targeting multiple aspects of tumorigenesis.^42–44 The integration of agents that enhance BAX activity with CCL20-targeted therapies could further improve immune responses while promoting apoptosis.^45–47

This study hypothesized that studies exploring therapies targeting these genes could lead to promising treatment strategies. BCL2L1 is an antiapoptotic gene essential for the survival of many cell types that can be targeted for antitumor therapy. The promising effects of BCL2L1 inhibitors for many cancers provide a potential avenue for CSCC.^48,49 Harnessing BCL2L1 inhibitors with cancer therapies (eg, chemotherapy, immunotherapy, or radiotherapy) may potentiate treatment efficacy. While specific studies targeting RB1 and CTNNB1 in HPV-positive CSCC may be limited, the experience of targeted therapies for other cancers suggests potential opportunities for future research.

The pathway and process enrichment analysis reveals significant findings concerning the biological functions associated with the 10 upregulated genes, with Gene Ontology (GO) analysis identifying enrichment in 13 distinct biological processes, including multicellular organismal process, biological processes involved in interspecies interaction, locomotion, reproductive process, homeostatic process, immune system process, developmental process, response to stimulus, regulation of biological processes, and both negative and positive regulation of biological processes, cellular processes, and biological regulation. This broad range of enriched processes indicates a multifaceted role of the upregulated genes in various physiological and pathological contexts. Additionally, the enrichment analysis is associated with small-cell lung cancer, head and neck squamous cell carcinoma, regulation of mononuclear cell proliferation, head and neck squamous cell carcinoma, negative regulation of cell population proliferation, camera-type eye morphogenesis, and positive regulation of cell migration.^50,51 The enrichment of terms related to cellular processes and the regulation of cell migration aligns with literature indicating that TME notably impacts patients’ clinical outcomes.^52,53 The involvement of negative and positive regulation of cell population proliferation indicates dysregulation of cell cycle control.⁵² Studying the molecular landscapes linked with these enriched terms might lead to new prognostic biomarkers.

In this study, the analysis of the immune landscape associated with the top upregulated genes in HPV-associated CSCC highlighted their association with various immune cell infiltrates and TME components. Our findings align with previous studies emphasizing the importance of immune cell infiltration in cervical cancer prognosis and immunotherapy response. The observed positive associations between gene expression and infiltrating immune cell types, such as CD8⁺ T cells, CD4⁺ T cells, dendritic cells (DCs), macrophages, and B cells, underscore the immunogenic nature of the TME in CSCC. Similar trends have been reported in prior studies. For instance, one study demonstrated that increased infiltration of CD8⁺ T cells and DCs was associated with better survival outcomes in cervical cancer patients, supporting our observations of their abundance in high-gene-expression patients.⁵⁴ Additionally, immune infiltration patterns are increasingly recognized as predictors of response to checkpoint blockade therapies in cervical and other squamous malignancies.^55,56

Our analysis also identified regulatory immune cells, such as natural and induced T regulatory cells (nTregs and iTregs), in association with highly expressed redox genes. The presence of these immunosuppressive populations has been implicated in immune evasion in cervical cancer, contributing to tumor progression.⁵⁷ Moreover, the involvement of exhausted T cells and γδ T cells suggests a TME undergoing immune exhaustion, a feature commonly seen in advanced-stage cancers with persistent antigenic stimulation.⁵⁸

The stromal landscape also appeared to shape the immune contexture. Our findings of positive correlations between upregulated genes and cancer-associated fibroblasts (CAFs) are consistent with previous studies that showed CAFs can modulate immune exclusion by altering extracellular matrix dynamics and cytokine profiles.^59,60 This stromal–immune interaction may contribute to immunotherapy resistance by limiting effector cell infiltration into tumor nests.

Additionally, the lack of significant correlation between tumor purity and gene expression supports the robustness of our results. Similarly reported that accounting for tumor purity is critical in interpreting bulk tumor transcriptomic data, particularly when studying immune infiltration.⁶¹

Our use of GSVA further allowed us to quantify immune cell enrichment scores and assess the coordinated regulation of immune pathways. This integrative approach provides a broader view of the immune landscape and supports the potential functional relevance of these genes in shaping immune responses. Comparable analyses in breast and colorectal cancer have shown that immune metagenes strongly correlate with prognosis and therapy responsiveness.⁶¹ Overall, our results suggest that the most highly upregulated genes in HPV-associated CSCC are closely linked to immune modulation in CSCC, interacting with both effector and suppressive immune elements as well as stromal regulators like CAFs. These genes may serve as potential therapeutic targets aimed at remodeling the immune microenvironment. Further functional validation and integration with clinical datasets may pave the way for their application in immunotherapy approaches.

Limitations

This study focused on a subset of protein-coding genes, which may overlook other critical molecular contributors to tumor progression and the TME. A more comprehensive investigation incorporating non-coding RNAs, epigenetic modifications, and DNA methylation patterns is needed to achieve a holistic understanding of the regulatory networks involved in HPV-associated CSCC. In addition, the number of patients included in this study was relatively small, which may limit the generalizability of the findings and the statistical power of the analyses. Consequently, findings of this study should be interpreted as preliminary, with further validation needed in larger cohorts. Also, the study was cross-sectional in design, capturing static gene expression profiles at a single time point. While suitable for exploratory transcriptomic analyses, this approach does not permit evaluation of temporal dynamics in gene expression during disease progression or therapeutic response. Future longitudinal studies will be necessary to validate and expand upon these findings in the context of clinical outcomes.

This study also did not include in vitro or in vivo functional validation experiments using CRISPR-Cas9 knockout or siRNA-mediated silencing in cervical cancer cell lines and xenograft models. Such follow-up studies will be essential to investigate the mechanistic roles of the identified genes and to evaluate their potential as therapeutic targets or prognostic biomarkers in preclinical models. Lastly, although RNA quality was assessed using spectrophotometric and gel-based methods before analysis, RNA Integrity Number (RIN) scores were not documented. Incorporating RIN-based quality control measures in future studies will enhance methodological transparency and reproducibility.

Conclusion

This study revealed a complex interplay between gene expression, inflammation, and immune regulation within the tumor microenvironment (TME) of HPV-positive cervical squamous cell carcinoma (CSCC). By analyzing 168 immune-associated, oncogenic, and tumor suppressor genes, we identified a subset of significantly upregulated genes potentially involved in tumor progression, immune evasion, and modulation of the TME. Alongside, the downregulation of several genes suggests disruptions in critical signaling pathways that may contribute to unchecked tumor growth.

Notably, six of the upregulated genes demonstrated significant associations with improved overall survival (OS) and disease-free survival (DFS), underscoring their potential as prognostic biomarkers. Immune infiltration analysis further revealed correlations between gene expression and the presence of key immune cell types, highlighting the relevance of these genes in shaping immune dynamics within the tumor milieu.

Overall, the identified gene candidates and associated pathways represent promising targets for future therapeutic intervention and warrant further investigation in translational and clinical contexts.

Future Directions

To strengthen and expand upon our current findings, several avenues will be pursued. First, functional validation of the six survival-associated genes (CASP8, BAX, FOXP3, CD274, RB1, CCL18) will be performed using CRISPR-Cas9 gene knockout and siRNA-mediated silencing in HPV-positive cervical cancer cell lines. These experiments will evaluate effects on tumor cell proliferation, apoptosis, immune evasion, and cytokine signaling. Second, in vivo studies, including xenograft or syngeneic mouse models, will be used to assess how manipulation of these genes influences tumor growth and immune cell infiltration. Third, validation in larger, independent CSCC cohorts, ideally with paired immune profiling and clinical outcome data, will help confirm the prognostic value of the identified gene signature. Lastly, integrating single-cell RNA sequencing and spatial transcriptomics will allow mapping the cell-type specificity and spatial distribution of immune-oncogene interactions within the TME, offering deeper insights into their functional relevance and therapeutic potential.

Data Sharing Statement

The data presented in this study are available on request from the corresponding author.

Institutional Review Board Statement

The study received ethical approval from the Institutional Review Board (IRB) of King Fahad Medical City (IRB number 22-008). The study was conducted in accordance with the principles outlined in the Declaration of Helsinki. Patient consent for the use of retrospectively collected tissue samples and to review medical records was not required, as the Ethics Committee classified the study as minimal risk and granted a waiver of consent in line with national regulations and institutional policies for retrospective studies involving de-identified data. All participant data were treated with strict confidentiality; each patient was assigned a unique identifier, and the data were securely stored in a password-protected Excel file.

Acknowledgments

This research project was supported by a grant from the Ongoing Research Funding Program (ORF) number (ORF-2025-543), King Saud University, Riyadh, Saudi Arabia. The authors extend their appreciation to the Ongoing Research Funding Program number (ORF-2025-543), King Saud University, Riyadh, Saudi Arabia. This research project was supported by Princess Nourah bint Abdulrahman University Researchers Supporting Project number (PNURSP2025R299), Princess Nourah bint Abdulrahman University, Riyadh, Saudi Arabia.

Author Contributions

All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis and interpretation, or in all these areas; took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.

Funding

Princess Nourah bint Abdulrahman University Researchers Supporting Project number (PNURSP2025R299), Princess Nourah bint Abdulrahman University, Riyadh, Saudi Arabia. This research project wassupported by a grant from the Ongoing Research Funding Program (ORF) number (ORF-2025-543), King Saud University, Riyadh, Saudi Arabia. The authors extend their appreciation to the Ongoing Research Funding Program number (ORF-2025-543), King Saud University, Riyadh, Saudi Arabia.

Disclosure

The authors declare no conflicts of interest in this work.

References

1. Wang X, Huang X, Zhang Y. Involvement of human papillomaviruses in cervical cancer. Front Microbiol. 2018;9. doi:10.3389/fmicb.2018.02896.

2. Burmeister CA, Khan SF, Schäfer G, et al. Cervical cancer therapies: current challenges and future perspectives. Tumour Virus Res. 2022;13:200238. doi:10.1016/j.tvr.2022.200238

3. Lintao RCV, Cando LFT, Perias GAS, et al. Current status of human papillomavirus infection and cervical cancer in the Philippines. Front Med. 2022;9:929062. doi:10.3389/fmed.2022.929062

4. Alahmadi RM, Awadalla M, Marraiki N, et al. “profiling the tumor immune microenvironment of hpv-associated base of tongue squamous cell carcinoma,”. Onco Targets Ther. 2025;18:263–281. doi:10.2147/ott.S505376

5. Molina MA, Steenbergen RDM, Pumpe A, Kenyon AN, Melchers WJG. HPV integration and cervical cancer: a failed evolutionary viral trait. Trends Mol Med. 2024;30(9):890–902. doi:10.1016/j.molmed.2024.05.009

6. Yu X, et al. comprehensive analysis of the carcinogenic process, tumor microenvironment, and drug response in HPV-positive cancers. Front Oncol. 2022;12. doi:10.3389/fonc.2022.842060

7. Landskron G, De la Fuente M, Thuwajit P, Thuwajit C, Hermoso MA. Chronic inflammation and cytokines in the tumor microenvironment. J Immunol Res. 2014;2014:149185. doi:10.1155/2014/149185

8. Nengroo MA, Verma A, Datta D. Cytokine chemokine network in tumor microenvironment: impact on CSC properties and therapeutic applications. Cytokine. 2022,;156:155916. doi:10.1016/j.cyto.2022.155916

9. Trujillo-Cirilo L, Weiss-Steider B, Vargas-Angeles CA, Corona-Ortega MT, Rangel-Corona R. Immune microenvironment of cervical cancer and the role of IL-2 in tumor promotion. Cytokine. 2023,;170:156334. doi:10.1016/j.cyto.2023.156334

10. De Nola R, Loizzi V, Cicinelli E, Cormio G. Dynamic crosstalk within the tumor microenvironment of uterine cervical carcinoma: baseline network, iatrogenic alterations, and translational implications. Crit Rev Oncol/Hematol. 2021,;162:103343. doi:10.1016/j.critrevonc.2021.103343

11. Khaliulin M, Valiullina A, Petukhov A, Yuan Y, Spada S, Bulatov E. Breaking the shield of solid tumors: a combined approach for enhanced efficacy of CAR-T cells. Cancer Immunol Immunother. 2024,;74(1):3. doi:10.1007/s00262-024-03817-z

12. Abdul-Rahman T, Ghosh S, Badar SM, et al. The paradoxical role of cytokines and chemokines at the tumor microenvironment: a comprehensive review. Eur J Med Res. 2024,;29(1):124. doi:10.1186/s40001-024-01711-z

13. Kartikasari AER, Huertas CS, Mitchell A, Plebanski M. Tumor-induced inflammatory cytokines and the emerging diagnostic devices for cancer detection and prognosis. Front Oncol. 2021;11:692142. doi:10.3389/fonc.2021.692142

14. Zhu C, Preissl S, Ren B. Single-cell multimodal omics: the power of many. Nat Methods. 2020;17(1):11–14. doi:10.1038/s41592-019-0691-5

15. Moncada R, Barkley D, Wagner F, et al. Integrating microarray-based spatial transcriptomics and single-cell RNA-seq reveals tissue architecture in pancreatic ductal adenocarcinomas. Nat Biotechnol. 2020;38(3):333–342. doi:10.1038/s41587-019-0392-8

16. Alahmadi RM, Marraiki N, Alswayyed M, et al. comprehensive transcriptome analysis reveals the distinct gene expression patterns of tumor microenvironment in HPV-Associated and HPV-non associated tonsillar squamous cell carcinoma. Cancers. 2023;15(23). doi:10.3390/cancers15235548

17. Zhou Y, et al. Tumor biomarkers for diagnosis, prognosis and targeted therapy. Signal Transduction Targeted Ther. 2024,;9(1):132. doi:10.1038/s41392-024-01823-2

18. Carreras J, Kikuti YY, Roncador G, et al. High expression of caspase-8 associated with improved survival in diffuse large b-cell lymphoma: machine learning and artificial neural networks analyses. BioMedInformatics. 2021;1(1):18–46. doi:10.3390/biomedinformatics1010003

19. Gryko M, Łukaszewicz-Zając M, Guzińska-Ustymowicz K, Kucharewicz M, Mroczko B, Algirdas U. The caspase-8 and procaspase-3 expression in gastric cancer and non-cancer mucosa in relation to clinico-morphological factors and some apoptosis-associated proteins. Adv Med Sci. 2023,;68(1):94–100. doi:10.1016/j.advms.2023.02.001

20. Stupack DG. Caspase-8 as a therapeutic target in cancer. Cancer Lett. 2013;332(2):133–140. doi:10.1016/j.canlet.2010.07.022

21. Zhang W, Zhu C, Liao Y, Zhou M, Xu W, Zou Z. Caspase-8 in inflammatory diseases: a potential therapeutic target. Cell. Mol. Biol. Lett. 2024,;29(1):130. doi:10.1186/s11658-024-00646-x

22. Passaro A, Al Bakir M, Hamilton EG, et al. Cancer biomarkers: emerging trends and clinical implications for personalized treatment,” (in eng. Cell. 2024;187(7):1617–1635. doi:10.1016/j.cell.2024.02.041

23. Guan B, Xu M, Zheng R, Guan G, Xu B. Novel biomarkers to predict treatment response and prognosis in locally advanced rectal cancer undergoing neoadjuvant chemoradiotherapy. BMC Cancer. 2023,;23(1):1099. doi:10.1186/s12885-023-11354-8

24. Karanikas V, Speletas M, Zamanakou M, et al. Foxp3 expression in human cancer cells. J Transl Med. 2008;6:19. doi:10.1186/1479-5876-6-19

25. Zhou J, Liu Y, Zhang Y, et al. Comprehensive analysis of a novel subtype of immune microenvironment-derived HPV-infected colorectal cancer. Microb Infect. 2024,;26(4):105315. doi:10.1016/j.micinf.2024.105315

26. Shi JY, Ma L-J, Zhang J-W, et al. FOXP3 Is a HCC suppressor gene and Acts through regulating the TGF-β/Smad2/3 signaling pathway. BMC Cancer. 2017,;17(1):648. doi:10.1186/s12885-017-3633-6

27. Liu L, Xiao W, Zhang C, Fan P, Zeng J, Yi J. The potential of foxp3 in predicting survival and treatment response in breast cancer. Int J Gen Med. 2024;17:1233–1251. doi:10.2147/ijgm.S454421

28. Fiori Lopes L, Losi Guembarovski R, Guembarovski AL, et al. FOXP3 transcription factor: a candidate marker for susceptibility and prognosis in triple negative breast cancer. Biomed Res. Int. 2014;2014(1):341654. doi:10.1155/2014/341654

29. Yang S, Liu Y, Li M-Y, et al. FOXP3 promotes tumor growth and metastasis by activating Wnt/β-catenin signaling pathway and EMT in non-small cell lung cancer. Mol Cancer. 2017,;16(1):124. doi:10.1186/s12943-017-0700-1

30. Zhou X, Wang Y, Zheng J, et al. Integrative study reveals the prognostic and immunotherapeutic value of CD274 and PDCD1LG2 in pan-cancer. Front Genetics. 2022;13. doi:10.3389/fgene.2022.990301

31. Schalper KA, Brown J, Carvajal-Hausdorf D, et al. Objective measurement and clinical significance of TILs in non-small cell lung cancer. J Natl Cancer Inst. 2015;107(3). doi:10.1093/jnci/dju435

32. Gao G, Zhang X-D, Qu H, et al. A comprehensive pan-cancer analysis of CD274 gene amplification, tumor mutation burden, microsatellite instability, and PD-L1 expression in Chinese cancer patients. Ann Transl Med. 2021;9(8):677. doi:10.21037/atm-21-853

33. Huang CY, Chiang S-F, Ke T-W, et al. Clinical significance of programmed death 1 ligand-1 (CD274/PD-L1) and intra-tumoral CD8+ T-cell infiltration in stage II–III colorectal cancer. Sci Rep. 2018,;8(1):15658. doi:10.1038/s41598-018-33927-5

34. Lohavanichbutr P, Houck J, Fan W, et al. Genomewide gene expression profiles of HPV-positive and HPV-negative oropharyngeal cancer: potential implications for treatment choices. Arch Otolaryngol Head Neck Surg. 2009;135(2):180–188. doi:10.1001/archoto.2008.540

35. Sasa N, Kishikawa T, Mori M, et al. Intratumor heterogeneity of HPV integration in HPV-associated head and neck cancer. Nat Commun. 2025,;16(1):1052. doi:10.1038/s41467-025-56150-z

36. Wang YW, Zhang K, Zhao S, et al. HPV status and its correlation with BCL2, p21, p53, Rb, and survivin expression in breast cancer in a Chinese population. Biomed Res Int. 2017;2017:6315392. doi:10.1155/2017/6315392

37. Karstensen B, Poppelreuther S, Bonin M, Walter M, Iftner T, Stubenrauch F. Gene expression profiles reveal an upregulation of E2F and downregulation of interferon targets by HPV18 but no changes between keratinocytes with integrated or episomal viral genomes. Virology. 2006,;353(1):200–209. doi:10.1016/j.virol.2006.05.030

38. Chai J, Lei Y, Xiang X, Ye J, Zhao H, Yi L. High expression of caspase-8 as a predictive factor of poor prognosis in patients with esophageal cancer. Cancer Med. 2023;12(6):7651–7666. doi:10.1002/cam4.5496

39. Qian S, Wei Z, Yang W, Huang J, Yang Y, Wang J. The role of BCL-2 family proteins in regulating apoptosis and cancer therapy. Front Oncol. 2022;12:2022–October–12. doi:10.3389/fonc.2022.985363

40. Liang Y, Wang L, Ma P, Ju D, Zhao M, Shi Y. Enhancing anti-tumor immune responses through combination therapies: epigenetic drugs and immune checkpoint inhibitors. Front Immunol. 2023;14:1308264. doi:10.3389/fimmu.2023.1308264

41. Varayathu H, Sarathy V, Thomas BE, Mufti SS, Naik R. Combination strategies to augment immune check point inhibitors efficacy – implications for translational research. Front Oncol. 2021;11:559161. doi:10.3389/fonc.2021.559161

42. Knudsen ES, Nambiar R, Rosario SR, Smiraglia DJ, Goodrich DW, Witkiewicz AK. Pan-cancer molecular analysis of the RB tumor suppressor pathway. Commun. Biol. 2020,;3(1):158. doi:10.1038/s42003-020-0873-9

43. AmeliMojarad M, AmeliMojarad M, Cui X, Shariati P. Pan-cancer analysis of CTNNB1 with potential as a therapeutic target for human tumorigenesis. Inf Med Unlocked. 2023,;42:101331. doi:10.1016/j.imu.2023.101331

44. Katkoori VR, Suarez-Cuervo C, Shanmugam C, et al. Bax expression is a candidate prognostic and predictive marker of colorectal cancer. J Gastrointestinal Oncol. 2010;1(2):76–89. doi:10.3978/j.issn.2078-6891.2010.019

45. Kwantwi LB, Boafo JD, Egleh BE, Li M. CCL20 in the tumor microenvironment: implications for cancer progression and therapeutic approaches. Clin Transl Oncol. 2025,;27:3285–3292. doi:10.1007/s12094-025-03874-5

46. Muscella A, Vetrugno C, Marsigliante S. CCL20 promotes migration and invasiveness of human cancerous breast epithelial cells in primary culture. Molecular Carcinogenesis. 2017;56(11):2461–2473. doi:10.1002/mc.22693

47. Zhang L, Yu J, Park BH, Kinzler KW, Vogelstein B. Role of BAX in the Apoptotic Response to Anticancer Agents. Science. 2000;290(5493):989–992. doi:10.1126/science.290.5493.989

48. Loo LSW, Soetedjo AAP, Lau HH, et al. BCL-xL/BCL2L1 is a critical anti-apoptotic protein that promotes the survival of differentiating pancreatic cells from human pluripotent stem cells. Cell Death Dis. 2020,;11(5):378. doi:10.1038/s41419-020-2589-7

49. Chiou J-T, Wu -Y-Y, Lee Y-C, Chang L-S. BCL2L1 inhibitor A-1331852 inhibits MCL1 transcription and triggers apoptosis in acute myeloid leukemia cells. Biochem Pharmacol. 2023,;215:115738. doi:10.1016/j.bcp.2023.115738

50. Li Q, et al. Signaling pathways involved in colorectal cancer: pathogenesis and targeted therapy. Signal Transduction Targeted Ther. 2024,;9(1):266. doi:10.1038/s41392-024-01953-7

51. Koveitypour Z, Panahi F, Vakilian M, et al. Signaling pathways involved in colorectal cancer progression. Cell Biosci. 2019,;9(1):97. doi:10.1186/s13578-019-0361-4

52. Tian Y, Li Q, Yang Z, et al. Single-cell transcriptomic profiling reveals the tumor heterogeneity of small-cell lung cancer. Signal Transduction Targeted Ther. 2022,;7(1):346. doi:10.1038/s41392-022-01150-4

53. Liu D, Wen J, Chen J, et al. A comparative analysis of the gene expression profiles of small cell esophageal carcinoma, small cell lung cancer, and esophageal adeno/squamous carcinoma. Frontiers in Surg. 2021;8:2021–July–30. doi:10.3389/fsurg.2021.655159

54. Zhou R, Xie Y, Wang Z, et al. “Single-cell transcriptomic analysis reveals CD8 + T cell heterogeneity and identifies a prognostic signature in cervical cancer,”. BMC Cancer. 2025,;25(1):498. doi:10.1186/s12885-025-13901-x

55. Li X, Pan L, Li W, et al. Deciphering immune predictors of immunotherapy response: a multiomics approach at the pan-cancer level. Cell Rep Med. 2025;6(4):101992. doi:10.1016/j.xcrm.2025.101992

56. Chen F, Shen L, Wang Y, et al. Signatures of immune cell infiltration for predicting immune escape and immunotherapy in cervical cancer. Aging. 2023;15(5):1685–1698. doi:10.18632/aging.204583

57. Kim CW, Kim J, Park Y, et al. ERRATUM: prognostic implications of extranodal extension in relation to colorectal cancer location. Cancer Res Treat. 2021;53(3):893. doi:10.4143/crt.2018.392.E

58. Chen D, Guo Y, Jiang J, et al. γδ T cell exhaustion: opportunities for intervention. J Leukoc Biol. 2022;112(6):1669–1676. doi:10.1002/jlb.5mr0722-777r

59. Arima Y, Matsueda S, Saya H. Significance of cancer-associated fibroblasts in the interactions of cancer cells with the tumor microenvironment of heterogeneous tumor tissue. Cancers. 2023;15(9):2536. doi:10.3390/cancers15092536

60. Yang D, Liu J, Qian H, Zhuang Q. Cancer-associated fibroblasts: from basic science to anticancer therapy. Exp. Mol. Med. 2023,;55(7):1322–1332. doi:10.1038/s12276-023-01013-0

61. Haider S, Tyekucheva S, Prandi D, et al. Systematic assessment of tumor purity and its clinical implications. JCO Precis Oncol. 2020;4. doi:10.1200/po.20.00016