>

Identification of Hub Genes and Key Pathways Associated with Human Cell Carcinoma Based on Gene Expression Papillomavirus Status in Cervical Squamous Profiling via Integrated Bioinformatics

 

Marah Massoud ⁽¹⁾ Majd Aljamali ^(1,2) Lama Youssef ^(1,3)

⁽¹⁾ Syrian Virtual University, ⁽²⁾ Dept. Biochemistry and Microbiology, Faculty of Pharmacy, Damascus University, ⁽³⁾ Dept. Pharmaceutics and Clinical Pharmacy, Faculty of Pharmacy, Damascus University

marahmass94@gmail.com

 

Abstract

Using integrated bioinformatics, it may be possible to identify the pathogenic mechanism underlying the development of the tumor by analyzing differentially expressed genes (DEGs) linked with the two HPV statuses (HPV positive and HPV negative) in cervical squamous cell carcinoma.

In addition, the discovery of significant differentially expressed genes, enrichment of their biological functions and key pathways, and visualization of the network of DEGs and hub genes will produce more accurate and trustworthy biomarkers for early diagnosis, individualized preventative measures, and improved therapeutic efficacy.

In this study, a series of analyses were conducted using R2 software of HPV status in squamous cervical carcinoma-related data from TCGA database to screen and identify prognostic biomarkers related to differentially expressed genes. Then, the up- and downregulated DEGs were classified into three groups according to Gene Ontology (GO) terms, whereas KEGG pathway enrichment analysis was conducted using DAVID website. To further study the potential relationships between the genes, protein-protein interaction (PPI) network was created using STRING database and reconstructed via Cytoscape software. The Cytoscape plug-in CytoHubba was used to identify the hub genes. We used Kaplan-Meier plotter to assess the survival rate of cancer patients and the prognostic value of the extracted key genes. The UALCAN online tool was then used to confirm the expression of these important genes in squamous cell cervical cancer.

On this basis, eleven hub genes (STAT1, CTNNB1, IRF9, EGFR, RSAD2, IRF7, MX1, IFIH1, IRF5, IRF1, DDX58) were defined.

Cancer Atlas investigation revealed that four proteins stained differently in tumor tissues compared to normal tissues (RSAD2, MX1, STAT1, DDX58). Finally, functional enrichment analysis was conducted to predict gene-drug interactions. Interestingly, we identified a number of medecations that could have significant potential in the future management of this cancer type.

 

Keywords: bioinformatics analysis, biomarkers, Cervical Squamous Cell Carcinoma, prognosis, differentially expressed genes

 

 

 

 

 

 

 

Introduction

With around 450,000 new cases diagnosed each year, cervical cancer continues to be the second-most prevalent cause of cancer-related deaths in women globally (1).

The prognosis and treatment plan can be determined by identifying the kind of cervical cancer. Yet squamous cell carcinomas are the most common type of cervical cancer. Squamous Cell Carcinoma is a type of cervical cancer that begins in the thin, flat cells (squamous cells) lining the outer part of the cervix, which projects into the vagina. In developing countries, cervical cancer is the most common cancer in women and may constitute up to 25% of all female cancers. In 1996, the World Health Association, along with the European Research Organization on Genital Infection and Neoplasia and the National Institutes of Health Consensus Conference on Cervical Cancer, recognized HPV as an important cause of cervical cancer. Additionally, HPV has been implicated in 99.7% of cervical squamous cell cancer cases worldwide (2) (3). Various strains of Human papillomavirus (HPV), a sexually transmitted infection, represents the most important risk factor for the development of cervical cancer (1). When exposed to HPV, the body's immune system typically prevents the virus from doing harm. In a small percentage of people, the virus survives for years, contributing to the process that causes some cervical cells to become cancer cells. Prevalence surveys, large case-control studies, and case series have unequivocally shown that HPV DNA can be detected in cervical cancer specimens in 90–100% of cases, compared with a prevalence of 5–20% in cervical specimens from women identified as suitable epidemiological controls (4). Although more than 100 distinct HPV genotypes have been described, and at least 20 are associated with cervical cancer, HPV types 16 and 18 are the most frequently detected in cervical cancer regardless of the geographical origin of the patients (4). Despite the fact that early stage cervical cancer can be cured by radical surgery or radiotherapy with equal effectiveness (5), pelvic radiation represents the standard therapy for the treatment of locally advanced disease. Despite technological advances, however, up to 35% of patients overall will develop advanced, metastatic disease, for which treatment results are poor. A deeper understanding of the molecular basis of cervical cancer has the potential to refine significantly the diagnosis and management of these tumors and may eventually lead to the development of novel, more specific and more effective treatments for prevention of disease progression following first-line therapy.

Demographic and exposure differences between HPV-positive (HPV+) and negative (HPV−) Cervical Squamous Cell Carcinoma suggest that HPV+ tumors may constitute a subclass with different biology, whereas clinical differences have also been observed. Further investigation of differentially expressed genes may reveal the unique pathways in HPV+ tumors that may explain the different natural history and biological properties of these tumors. These properties may be exploited as a target of novel therapeutic agents in cervical squamous cell carcinoma treatment. The obtained hub genes and key pathways could be the therapeutic targets for the precise treatment of these two HPV status with different prognoses. using integrated bioinformatics to screen differentially expressed genes (DEGs) in cervical squamous cell cancer could benefit us for understanding the pathogenic mechanism underlying the tumor progression. Also to identify significantly upregulated genes in differentially expressed genes (DEGs) from both HPV positive and HPV negative status, which might be used as the biomarkers for early diagnosis and prevention of the disease.

 

 

Materials and Methods

Obtaining processed differentially expressed genes and data analysis

R2 platform (https://hgserver1.amc.nl/cgi-bin/r2/) was used to obtain the normalized RNA-seq data samples of patients with Cervical Squamous Cell Carcinoma derived from TCGA database (https://portal.gdc.cancer.gov/), data were collected from 305 samples including both HPV status (including HPV negative: 22 patients, and HPV positive: 281 patients). Analysis was performed by R2 to obtain differentially expressed genes (DEGs) and their expression levels in HPV positive and HPV negative samples (Figure 1), multiple testing correction applied is False. Discovery Rate(FDR), p-value cutoff =< 0.01 and |log2 fold change (FC)| >=1. The results were considered statistically significant.

HPV status data stands for:

1.  HPV positive for high risk HPV subtypes as they are associated with cervical cancer development: 16, 18, 35, 33, 31, 45, 52, and 58

2.  HPV negative: not associated with any HPV subtype infection.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(Figure 1) scatter plot of differentially expressed genes (DEGs) and their expression levels associated with HPV positive and HPV negative in Cervical Squamous Cell Carcinoma samples colored by p-value via R2 software. Each dot represents a different gene; the darkness of the color is proportional to the level of expression.

 

 

 

 

 

Functional Enrichment and Pathway Analysis of Differentially Expressed Genes (DEGs)

DAVID (6)(7)was utilized to conduct Gene ontology (GO) and Kyoto Gene and Genome Encyclopedia (KEGG) enrichment analysis (https://david.ncifcrf.gov/home.jsp).

The DAVID online analysis was used to conducted biological annotation of the identified common differentially expressed genes from integrated analysis of RNA-seq data in squamous cervical cancer. We obtained Gene ontology (GO) functional enrichments of up- and downregulated genes with a P-value<0.05

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

as a cutoff value. Three functional groups,

including molecular function (MFs), biological processes (BPs), and cell composition (CCs) and KEGG pathways of the enriched genes, were divided in GO analysis of the common DEGs (Table1). In the biological processes group, the identified DEGs were mainly enriched in:

defense response to virus, regulation of transcription from RNA polymerase II promoter, negative regulation of viral genome replication, and in the cellular composition group, genes were mainly categorized in cytoplasm, nucleus, nucleoplasm, plasma membrane. 

 

In the molecular function group, genes were mainly concentrated in protein binding, RNA polymerase II core promoter proximal region sequence-specific DNA binding, RNA polymerase II transcription factor activity, sequence-specific DNA binding, transcription factor activity, sequence-specific DNA binding, metal ion binding. Pathway enrichment analysis revealed the enrichment of the differentially expressed genes in many pathways: Herpes simplex virus 1 infection, Pathways in cancer and Epstein-Barr virus infection.

 

 

 

 

Construction of the Protein-Protein Interaction (PPI) Network for the DEGs

To further study the potential relationships between the genes, STRING online database (https://string-db.org/) was used to mine and provide the significant interactions among these genes (8) (9).

(Figure 2-A), and the degree of connectivity of each node of the network was calculated. The Cytoscape plug-in cytoHubba was used to identify the hub genes by finding the intersections of the first 30 genes from 12 topological analysis methods (Figures 2-B). Afterward, we used molecular complex detection (MCODE) to find clusters according to topology locating densely connected regions. MCODE was used to establish the important modules from the PPI network (Figure 2-C) with a degree cutoff = 2, node score cutoff = 0.2, k‐core = 2, and max depth = 100. Cluster 1 includes 27 nodes and 351 edges (Table 2).

Metascape tools (https://metascape.org/) were used to analyze the pathways and biological processes enrichment of hub genes (11). It was observed that key genes are enriched in defense response to virus, Interferon alpha/beta signaling, cellular response to cytokine stimulus, etc (Figure 3).

 

 

(Figure 2) The network of Differentially expressed genes associated with HPV positive and HPV negative status in Cervical Squamous Cell Carcinoma is shown in pink via Cytoscape software, hub genes network is highlighted in green generated via CytoHubba plug-in(A). Hub genes network of Differentially expressed genes associated with HPV positive and HPV negative status in Cervical Squamous Cell Carcinoma generated via CytoHubba plug-in(B). Best ranked cluster of the differentially expressed genes associated with HPV positive and HPV negative status in Cervical Squamous Cell Carcinoma according to MCODE app results(C).

 

 

 

 

 

 

 

 

(Figure 3) Network of enriched pathways and biological processes of differentially expressed hub genes (key genes that were extracted from CytoHubba in the previous step) associated with HPV positive and HPV negative status in Cervical Squamous Cell Carcinoma samples created by Metascape website: colored by cluster ID, where nodes that share the same cluster ID are typically close to each other.

 

Kaplan-Meier Survival Curve of the Differentially Expressed Genes and Screening of Prognostic Biomarkers

Gene expression data can be used through Kaplan-Meier Plotter(https://kmplot.com/) to estimate the survival rate of cancer patients. This tool's primary objective is to evaluate biomarkers using meta-analysis (12). In attempts to identify genes that can be employed as prognostic biomarkers for this disease, the effects of the hub genes are assessed on the prognosis of squamous cervical cancer patients using the Kaplan-Meier chart (Figure 4).

(Figure 4) Kaplan-Meier plots of differentially expressed hub genes associated with HPV positive and HPV negative status in Cervical Squamous Cell Carcinoma samples with

significant p-value. Kaplan-Meier Plotter use gene expression data to assess the survival rate.

Hub Gene Verification Through UALCAN

UALCAN is an online tool (http://ualcan.path.uab.edu/) with data from TCGA and GTEx(13) (14), was used to verify the expression of these key genes in squamous cell cervical cancer. In this study, according to the RNA sequence data from TCGA database, the mRNA expression levels of 30 hub genes were compared between the squamous cervical cancer samples and the adjacent normal tissues. Twelve genes were found to be highly expressed at the transcriptional level in cancer tissues compared with normal tissues with statistical significance (HLA-B, IFIH1, STAT1, STAT2, IRF7, RSAD2, IFIT3, IRF5, OAS2, EGFR, MX1, OASL). (Figure 5).

 

 

 

 

 

 

 

 

 

 

(Figure 5) Expression analysis of key mutated genes of HPV positive and HPV negative status in Cervical Squamous Cell Carcinoma samples with significant p-value via UALCAN. The figure demonstrates nine of twelve key genes that found to be highly expressed in Cervical Squamous Cell Carcinoma tissues compared with normal tissues.

Analysis of Cancer Genomics Data Through cBioPortal

The cBioPortal (6) for Cancer Genomics online platform (http://www.cbioportal.org/) provides resources for visualizing and analyzing multidimensional cancer genomics data. A graphical analysis of the genetic variation was done via OncoPrint tool to investigate the clinical significance of the hub genes. As shown in the figure, 30 key genes all showed a high mutation rate in squamous cell cervical cancer, with a rate of genome change ranging from 2.5% to 15% (Figure 6). OncoPrint from cBioPortal revealed that 68% of cases (187 out of 275) exhibited genetic alterations, including amplification, deep deletion, missense mutation, severe depletion, truncating mutation, and various mutations.

One limitation of the cBioportal analysis was that none of the cases were divided into HPV+ and HPV- groups. However, genetic alterations and up-or downregulation of these hub genes could be demonstrated.

 

 

 

 

 

(Figure 6) Graphic analysis of the genetic alteration of key genes in squamous cell cervical cancer via cBioPortal platform. OncoPrint from cBioPortal revealed that 68% of cases (187 out of 275) exhibited genetic alterations, including amplification, deep deletion, missense mutation, severe depletion, truncating mutation, and various mutations.

 

Screening and Survival Analysis of Pivotal Genes

Using the intersection of the first 30 genes in cytoHubba’s 12 algorithms, 11 key genes were identified: STAT1, CTNNB1 IRF9, EGFR, RSAD2, IRF7, MX1, IFIH1, IRF5, IRF1 and DDX58. DAVID software was used to analyze the pathway and biological process enrichment of these eleven hub genes. It was observed that key genes are enriched in the defense response to virus, sequence-specific DNA binding, Hepatitis C pathway. The cBioPortal online platform provided a graphic analysis of the genetic variation of the hub genes.

As shown in Figure 6, eleven key DEGs all showed a high mutation rate in SCC, with a rate of genome change ranging from 4% to 7%. To determine whether the selected hub genes have clinical correlations, Kaplan-Meier curves was used to analyze the univariate survival of these genes and found that the expression of those genes was correlated with prognosis (as shown in the plots above). Thus, these genes can be used as prognostic indicators of squamous cervical cancer. UALCAN, an online tool with data from TCGA and GTEx, was used to verify the expression of these key genes in squamous cervical cancer.

In this study, and according to the RNA sequence data from TCGA database, the mRNA expression levels of 11 genes were compared between the tumor samples and the adjacent normal tissues. These genes were found to be highly expressed at the transcriptional level, except one gene” CTNNB1” was found to be highly expressed in normal tissues than tumor tissues with statistical significance=7.31^E-01. P-Value=0.05 was chosen as a threshold.

Human Protein Atlas as a tool to facilitate clinical biomarker discovery

It has been widely acknowledged that protein expression patterns in a tumour may provide critical diagnostic and prognostic information, and that immunohistochemistry can offer an essential new level of information on top of morphology (16). Protein expression profiles were mapped using immunohistochemistry in normal tissues, cancerous tissues, and cell lines. For each antibody, altogether 708 spots of tissues and cells are analyzed and the resulting protein expression data, including underlying high-resolution images, are published on the free and publically available Human Protein Atlas portal (proteinatlas.org/).Taking the aforementioned information into consideration, the expression level of each protein of the eleven hub genes, that were obtained from the previous step, in different tissues was obtained using the Cancer Atlas, and the comparison was done by detecting the level of antibody staining in normal tissue verses cancer tissues. it has been realized that STAT1 and RSAD2 are weakly stained or negative in normal tissues, while the expression level is ranged from low, moderate to highly stained in malignant tissues. Those two proteins/genes might be a novel biomarker for cervical squamous cell carcinoma.

MX1 and DDX58 displayed a moderate staining in normal tissues, while malignant tissues showed strong to moderate immunoreactivity.

Computational-Based Drug Repurposing using WebGestalt tool

When analyzing the findings of high-throughput experiments, which usually produce a list of intriguing genes or proteins, functional enrichment analysis is crucial. One of the most popular gene set enrichment analysis tools, WebGestalt (WEB-based GEne SeT AnaLysis Toolkit), aids users in deriving biological insights from genes of interest (17). Herein, analysis was proceeded in a suggested functional database (ie, Drug Bank) through WebGestalt to predict gene-drug interactions and potential therapeutic options for candidate genes. To predict the main significant drugs for ten hub genes (STAT1, STAT2, IRF1, IRF9, OAS2, IRF7, IFIT5, OAS1, RSAD2, RNASEL), “drug” and “DrugBank” were selected as the functional database and enrichment category, respectively. We considered the analysis as the significance level of P value < 0.05, and over-representation analysis as the method of Interest. WebGestalt gene-drug analysis predicted ten existing drugs (Figure 7)

(Figure 7) Tabular form the enriched drugs and potential therapeutic options for ten hub genes associated with HPV positive and HPV negative status in Cervical Squamous Cell Carcinoma samples using DrugBank database as a functional database via WebGestalt website.

 

 

 

 

 

 

 

 

 

Discussion

 

 

In 2020, an estimated 604,000 women were diagnosed with cervical cancer worldwide and about 342,000 women died from the disease. The main cause of cervical cancer is persistent infection with high-risk types of human papillomavirus (HPV). High incidence rates and high mortality rates of cervical cancer occur mainly (~90% for both) in low-and middle-income countries (18). To improve the quality of life, prognosis of patients and to prolong their survival time, researchers must further identify key genes that affect the development of this disease. Therefore, bioinformatics analysis of dysregulated genes in clinical samples of cervical cancer may pave the way for development of better prognostic markers and therapeutic targets. In the present study, a series of analyses were conducted using R2 software to explore downregulated and upregulated genes associated with HPV status in squamous cervical carcinoma. The normalized RNA-seq data samples of patients with Cervical Squamous Cell Carcinoma were derived from TCGA database. Some genes were recognized as HPV positive or HPV negative-specific markers. For instance, the expressions of STAT1, IRF7, DDX58, MX1, BRCA1 and FOS were related to HPV positive according to R2 results. In a previous study, it was confirmed that the expression of antiviral genes (IFIT1 and MX1), genes involved in IFN signaling (STAT1), and proapoptotic genes (TRAIL and XAF1), are inhibited to similar extents by HPV16, -18, and -31(19). In other hand, R2 results indicated that the expression levels of BMP4, CTNNB1, NKD1, WNT11 and STK36 were related to HPV negative status. Furthermore, studies showed that Wnt/CTNNB1 pathway is constitutively activated in squamous cell carcinoma (SCC) cell lines compared to normal keratinocytes (20). The differentially expressed genes were mainly enriched in the nucleus and cytoplasm, mainly regulating the defense response to virus. Pathway analysis showed that differentially expressed genes in cervical carcinoma are mainly involved in the Herpes simplex virus 1 and Epstein-Barr virus infection and pathways in cancer. Researches clarified that HPV-induced cellular changes facilitate the establishment of a latent EBV infection, such latent infections would allow for long-term expression of EBV oncogenes and EBV-induced epigenetic reprogramming that contribute to the progression of HPV-positive oropharyngeal squamous cell carcinoma (21). Additionally, increasing evidence indicates that dysregulation of Wnt signaling pathway cascade is contributed to the development and progression of some solid tumors and hematological malignancies (22) (23) (24). Hence, the genes involved in these pathways might provide a new direction of research on the original basis. Using integrated bioinformatics including CytoHubba plug-in, UALCAN online analysis and Kaplan-Meier curves survival analysis, the expression of eleven genes in squamous cervical cancer tissues was found to be significantly higher than that in normal tissues and the effects of the hub genes were evaluated on the prognosis of squamous cervical cancer patients, among which only one was upregulated in the HPV negative status (CTNNB1), while the others were upregulated in the HPV positive status (STAT1, IRF9, EGFR, RSAD2, IRF7, MX1, IFIH1, IRF5, IRF1 and DDX58). Interestingly, previous studies revealed that the genetic variants in the CTNNB1 gene might contribute to the development of cervical cancer and it might be a therapeutic target in this disease (25). The eleven key DEGs all showed a high mutation rate in SCC, with a rate of genome change ranging from 4% to 7% according to cBioPortal results. The expression levels of each protein in different tissues were obtained using the Cancer Atlas part from Human Protein Atlas portal tool. It has been realized that four proteins were differently stained in normal and tumor tissues: STAT1, RSAD2, MX1 and DDX58. STAT1 and RSAD2 are weakly stained or negative in normal tissues, while the expression level is ranged from low, moderate to highly stained in malignant tissues. One study used bioinformatics analysis to identify latent biomarkers in connection with progression and prognosis in oral cancer(OC), showed that RSAD2 is one of the independent prognostic indicators of OC (26). Another study showed that high GLUT4 RNA expression in combination with low DDX58 RNA expression levels was significantly correlated with the worst head and neck squamous cell carcinoma patient survival (27). STAT1 has been identified to have prognostic value in patients with solid cancer (28), and a major transcriptional target of human papillomavirus type 31(29). Therefore, those four proteins/genes might be a novel biomarkers and therapeutic targets for the precise treatment of cervical squamous cell carcinoma associated with HPV infection. Functional enrichment analysis was conducted using WebGestalt server through a suggested functional database (DrugBank) to predict gene-drug interactions and potential therapeutic options for ten ranked key genes that were generated by CytoHubba with their scores (STAT1, STAT2, IRF1, IRF9, OAS2, IRF7, IFIT5, OAS1, RSAD2, RNASEL). WebGestalt gene-drug analysis predicted ten existing drugs that were known to treat several types of cancer like lung cancer and breast cancer. It is worth noting that most of them are tyrosine kinase inhibitor like Afatinib [Giotrif®], which is an orally administered irreversible inhibitor of the ErbB family of tyrosine kinases that provides an important first-line treatment option for advanced non-small cell lung cancer (NSCLC) with activating epidermal growth factor receptor (EGFR) mutations, and an additional treatment option for squamous NSCLC that has progressed following first-line platinum-based chemotherapy (30).

 

Conclusions

Taken together, eleven differentially expressed genes (STAT1, CTNNB1, IRF9, EGFR, RSAD2, IRF7, MX1, IFIH1, IRF5, IRF1, and DDX58) have been discovered to be considerably enriched in a variety of biological processes and key pathways via the application of bioinformatics and computational approaches. Besides, they may be key factors in the occurrence and prognosis of squamous cervical cancer and may have a significant role in many pathways related to tumor development. Therefore, on this basis, further studies should be performed to detect the polymorphic sites of these genes and explore their corresponding expression levels, which can be used to predict the prognosis of patients. Additionally, functional enrichment analysis used to predict gene-drug interactions and potential therapeutic options. To ascertain their precision and sensitivity in carcinogenesis and to inform the patient's unique treatment plan, these findings will need to be validated in experimental and clinical research. Therefore, the goal of this study is to influence future research on cervical cancer by offering new perspectives on clinical diagnosis and prognostic assessment through bioinformatics analysis in order to guide follow-up researches on cervical cancer.

 

References

1.    Santin AD, Zhan F, Bignotti E, et al. Gene expression profiles of primary HPV16- and HPV18-infected early stage cervical cancers and normal cervical epithelium: identification of novel candidate molecular markers for cervical cancer diagnosis and therapy. Virology. 2005 Jan 20;331(2):269-91.

2.    Okunade KS. Human papillomavirus and cervical cancer. J Obstet Gynaecol. 2020 Jul;40(5):602-608. doi: 10.1080/01443615. 2019.1634030. 2019 Sep 10. Erratum in: J Obstet Gynaecol. 2020 May;40 (4):590.

3.    Walboomers JM, Jacobs MV, Manos MM, et al. Human papillomavirus is a necessary cause of invasive cervical cancer worldwide. J Pathol. 1999 Sep;189(1):12-9.

4.    Bosch FX, Lorincz A, Muñoz N, et al. The causal relation between human papillomavirus and cervical cancer. J Clin Pathol. 2002 Apr;55(4):244-65.

5.    Landoni F, Maneo A, Colombo A, et al. Randomized study of radical surgery versus radiotherapy for stage Ib-IIa cervical cancer, The Lancet, Volume 350, Issue 9077,1997, Pages 535-540.

6.    Sherman BT, Hao M, Qiu J, et al. DAVID: a web server for functional enrichment analysis and functional annotation of gene lists. Nucleic Acids Res. 2022 Mar 23;50 (W1): W216–21.

7.    Huang daW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44-57.

8.    Szklarczyk D, Gable AL, Lyon D, et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2019 Jan 8;47(D1): D607-D613.

9.    Szklarczyk D, Gable AL, Nastou KC, et al. The STRING database in 2021: customizable protein-protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Res. 2021 Jan 8;49(D1): D605-D612.

10.              Shannon P, Markiel A, Ozier O, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003 Nov;13(11):2498-504.

11.              Zhou Y, Zhou B, Pache L, et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun. 2019 Apr 3;10(1):1523.

12.              Lánczky A, Győrffy B. Web-Based Survival Analysis Tool Tailored for Medical Research (KMplot): Development and Implementation. J Med Internet Res. 2021 Jul 26;23(7):e27633.

13.              Chandrashekar DS, Karthikeyan SK, Korla PK, et al. UALCAN: An update to the integrated cancer data analysis platform. Neoplasia. 2022 Mar; 25:18-27.

14.              Chandrashekar DS, Bashel B, Balasubramanya SAH, et al. UALCAN: A Portal for Facilitating Tumor Subgroup Gene Expression and Survival Analyses. Neoplasia. 2017 Aug;19(8):649-658.

15.              Gao J, Aksoy BA, Dogrusoz U, et al. Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci Signal. 2013 Apr 2;6(269):pl1.

16.              Uhlén M, Björling E, Agaton C, et al. A Human Protein Atlas for Normal and Cancer Tissues Based on Antibody Proteomics, Molecular & Cellular Proteomics, Volume 4, Issue 12, 2005, Pages 1920-1932.

17.              Liao Y, Wang J, Jaehnig EJ, et al. WebGestalt 2019: gene set analysis toolkit with revamped UIs and APIs, Nucleic Acids Research, 2019,Volume 47, Issue W1, 02 , Pages W199–W205.

18.              Sung H, Ferlay J, Siegel RL, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin. 2021 May;71(3):209-249.

19.              Reiser J, Hurst J, Voges M, et al. High-risk human papillomaviruses repress constitutive kappa interferon transcription via E6 to prevent pathogen recognition receptor and antiviral-gene expression. J Virol. 2011 Nov;85(21):11372-80.

20.              Zanfi ED, Fantini S, Lotti R, et al. Wnt/CTNNB1 Signal Transduction Pathway Inhibits the Expression of ZFP36 in Squamous Cell Carcinoma, by Inducing Transcriptional Repressors SNAI1, SLUG and TWIST. Int. J. Mol. Sci. 2020, 21, 5692.

21.              Guidry JT, Myers JE, Bienkowska-Haba M, et al. Inhibition of Epstein-Barr Virus Replication in Human Papillomavirus-Immortalized Keratinocytes. J Virol. 2019 Jan 4;93(2):e01216-18.

22.              Zhu H, Mazor M, Kawano Y, et al. Analysis of Wnt gene expression in prostate cancer: mutual inhibition by WNT11 and the androgen receptor. Cancer Res. 2004 Nov 1;64(21):7918-26.

23.              Zhang M, Weng W, Zhang Q, et al. The lncRNA NEAT1 activates Wntβ-catenin signaling and promotes colorectal cancer progression via interacting with DDX5. J Hematol Oncol. 2018;11(1):113.

24.              Ge X and Wang X. Role of Wnt canonical pathway in hematological malignancies. J Hematol Oncol. 2010 Sep 15;3:33.

25.              Wang B, Li X, Liu L, et al. β-Catenin: oncogenic role and therapeutic target in cervical cancer. Biol Res 53, 33 (2020).

26.              Reyimu A, Chen Y, Song X, et al. Identification of latent biomarkers in connection with progression and prognosis in oral cancer by comprehensive bioinformatics analysis. World J Surg Onc, 2021, 19, 240.

27.              Chang YC, Chi LH, Chang WM,  et al. Glucose transporter 4 promotes head and neck squamous cell carcinoma metastasis through the TRIM24-DDX58 axis. J Hematol Oncol, 2017, 10, 11.

28.              Zhang J, Wang F, Liu F, Xu G. Predicting STAT1 as a prognostic marker in patients with solid cancer. Ther Adv Med Oncol. 2020 May 8;12:1758835920917558.

29.              Chang YE, Laimins LA. Microarray analysis identifies interferon-inducible genes and Stat-1 as major transcriptional targets of human papillomavirus type 31. J Virol. 2000 May;74(9):4174-82.

30.              Deeks ED, Keating GM. Afatinib in advanced NSCLC: a profile of its use. Drugs Ther Perspect. 2018;34(3):89-98.