Advertisement

Identification of key pathways and genes in nasopharyngeal carcinoma based on WGCNA

  • Author Footnotes
    1 Yongmei Dai and Wenhan Chen have contributed equally to this work.
    Yongmei Dai
    Correspondence
    Corresponding author at: Departments of Oncology, Shengli Clinical Medical College of Fujian Medical University & Fujian Provincial Hospital, No.134 East Street, Gulou District, Fuzhou 350001, China.
    Footnotes
    1 Yongmei Dai and Wenhan Chen have contributed equally to this work.
    Affiliations
    Departments of Oncology, Shengli Clinical Medical College of Fujian Medical University & Fujian Provincial Hospital, Fuzhou 350001, China
    Search for articles by this author
  • Author Footnotes
    1 Yongmei Dai and Wenhan Chen have contributed equally to this work.
    Wenhan Chen
    Footnotes
    1 Yongmei Dai and Wenhan Chen have contributed equally to this work.
    Affiliations
    The Second Clinical Medical College of Fujian Medical University, Fujian 362000, China

    Department of Clinical Medicine, Fujian Medical University, Fujian 350122, China
    Search for articles by this author
  • Junpeng Huang
    Affiliations
    Departments of Oncology, Shengli Clinical Medical College of Fujian Medical University & Fujian Provincial Hospital, Fuzhou 350001, China
    Search for articles by this author
  • Li Xie
    Affiliations
    Departments of Oncology, Shengli Clinical Medical College of Fujian Medical University & Fujian Provincial Hospital, Fuzhou 350001, China
    Search for articles by this author
  • Jianfang Lin
    Affiliations
    Departments of Oncology, Shengli Clinical Medical College of Fujian Medical University & Fujian Provincial Hospital, Fuzhou 350001, China
    Search for articles by this author
  • Qianshun Chen
    Affiliations
    Department of Thoracic Surgery, Fujian Provincial Hospital, Shengli Clinical College of Fujian Medical University, Fuzhou 350001, China
    Search for articles by this author
  • Guicheng Jiang
    Affiliations
    Departments of Oncology, Shengli Clinical Medical College of Fujian Medical University & Fujian Provincial Hospital, Fuzhou 350001, China
    Search for articles by this author
  • Chen Huang
    Affiliations
    Department of Thoracic Surgery, Fujian Provincial Hospital, Shengli Clinical College of Fujian Medical University, Fuzhou 350001, China
    Search for articles by this author
  • Author Footnotes
    1 Yongmei Dai and Wenhan Chen have contributed equally to this work.
Open AccessPublished:May 31, 2022DOI:https://doi.org/10.1016/j.anl.2022.05.013

      Abstract

      Objective

      We aim to identify the potential genes and signaling pathways associated with the nasopharyngeal carcinoma (NPC) prognosis using Weighted Gene Co-Expression Network Analysis (WGCNA).

      Methods

      Gene Expression Omnibus (GEO) query was utilized to download two NPC mRNA microarray data. WGCNA was conducted on differentially expressed genes (DEGs) to obtain tumor-associated gene modules. Genes in core modules were intersected with DEGs for gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analysis. GSE102349 dataset was devoted to identifying prognostic hub genes by survival analysis and the results were confirmed by quantitative polymerase chain reaction (qPCR).

      Results

      Co-expression networks were built, and we detected 12 gene modules. The Brown module and Magenta module were extremely associated with NPC samples. GO functional analysis and KEGG pathway analysis was carried out to the genes in the Brown and Magenta modules. Our data indicated that DEGs in Brown module and Magenta module were correlated with the biological regulation, metabolic process, reproduction, and cellular proliferation. Twenty-six hub genes were obtained and were considered to be closely related to NPC. GSE102349 dataset was devoted to identifying prognostic hub genes by survival analysis. The expression of IL33, MPP3 and SLC16A7 in GSE102349 dataset was significantly correlated with the progression-free survival (PFS). The results of qPCR indicated a strong correlation between SLC16A7 expression and the overall survival (OS).

      Conclusions

      WGCNA contributed to the detection of gene modules and identification of hub genes and crucial genes. These crucial genes might be potential targets for pharmaceutic therapies with potential clinical significance.

      Keywords

      1. Introduction

      Nasopharyngeal carcinoma (NPC) refers to an epithelial carcinoma arising from the nasopharyngeal mucosal lining [
      • Sung H.
      • Ferlay J.
      • Siegel R.L.
      • Laversanne M.
      • Soerjomataram I.
      • Jemal A.
      • et al.
      Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries.
      ]. Several factors have been considered to be associate with the pathogenesis of NPC including genetic susceptibility and environmental conditions [
      • Lin CT.
      Relationship between Epstein-Barr virus infection and nasopharyngeal carcinoma pathogenesis.
      ,
      • Stepan K.O.
      • Mazul A.L.
      • Skillington S.A.
      • Paniello R.C.
      • Rich J.T.
      • Zevallos J.P.
      • et al.
      The prognostic significance of race in nasopharyngeal carcinoma by histological subtype.
      ]. The treatment prognosis is closely related to local infiltration and lymphatic metastasis [
      • Zhang F.
      • Zhong L.Z.
      • Zhao X.
      • Dong D.
      • Yao J.J.
      • Wang S.Y.
      • et al.
      A deep-learning-based prognostic nomogram integrating microscopic digital pathology and macroscopic magnetic resonance images in nasopharyngeal carcinoma: a multi-cohort study.
      ]. However, there are large variances in the prognosis of patients with similar pathological types and staging even after the same combined-modality therapy [
      • Lyu Y.
      • Ni M.
      • Zhai R.
      • Kong F.
      • Du C.
      • Hu C.
      • et al.
      Clinical characteristics and prognosis of elderly nasopharyngeal carcinoma patients receiving intensity-modulated radiotherapy.
      ].
      Recently, extensive efforts have been made to investigate the efficiency of plasma Epstein-Barr virus (EBV) DNA copy number in predicting prognosis of NPC after treatment [
      • Sengar M.
      • Chorghe S.
      • Jadhav K.
      • Singh S.
      • Laskar S.G.
      • Pai P.
      • et al.
      Cell-free Epstein-Barr virus-DNA in patients with nasopharyngeal carcinoma: Plasma versus urine.
      ,
      • Wang H.
      • Liu W.
      • Luo B.
      The roles of miRNAs and lncRNAs in Epstein-Barr virus associated epithelial cell tumors.
      ]. For example, Lv et al. [
      • Lv J.
      • Chen Y.
      • Zhou G.
      • Qi Z.
      • Tan K.R.L.
      • Wang H.
      • et al.
      Liquid biopsy tracking during sequential chemo-radiotherapy identifies distinct prognostic phenotypes in nasopharyngeal carcinoma.
      ] revealed that real-time monitoring of liquid biopsy response added prognostic information. This suggested that pretreatment EBV-DNA may serve as an important index for establishing a treatment regimen. In contrast, post-treatment EBV-DNA rather than pre-treatment EBV-DNA was correlated with the overall survival (OS) among the NPC patients, which may serve as a predictable factor for the prognosis [
      • Lin J.C.
      • Wang W.Y.
      • Chen K.Y.
      • Wei Y.H.
      • Liang W.M.
      • Jan J.S.
      • et al.
      Quantification of plasma Epstein-Barr virus DNA in patients with advanced nasopharyngeal carcinoma.
      ]. Besides, several biological markers have been reported to be associated with the prognostic prediction of NPC, including epidermal growth factor receptor (EGFR) [
      • Huang F.
      • Liang X.
      • Min X.
      • Zhang Y.
      • Wang G.
      • Peng Z.
      • et al.
      Simultaneous Inhibition of EGFR and HER2 via Afatinib Augments the Radiosensitivity of Nasopharyngeal Carcinoma Cells.
      ], N-cadherin and β-catenin [
      • Sun H.
      • Liu M.
      • Wu X.
      • Yang C.
      • Zhang Y.
      • Xu Z.
      • et al.
      Overexpression of N-cadherin and β-catenin correlates with poor prognosis in patients with nasopharyngeal carcinoma.
      ], C-erbB2 [
      • Kim T.J.
      • Lee Y.S.
      • Kang J.H.
      • Kim Y.S.
      • Kang CS.
      Prognostic significance of expression of VEGF and Cox-2 in nasopharyngeal carcinoma and its association with expression of C-erbB2 and EGFR.
      ], p53 [
      • Cheng Y.
      • Wu S.
      • Tie X.
      • Huang X.
      • Cui L.
      Growth Inhibition of Nasopharyngeal Carcinoma Cells Mediated by p53 Gene-Containing Nanolipid Composites.
      ], interleukin-10 (IL-10) [
      • Savitri E.
      • Haryana MS.
      Expression of interleukin-8, interleukin-10 and Epstein-Barr viral-load as prognostic indicator in nasopharyngeal carcinoma.
      ], as well as vascular endothelial growth factor. Furthermore, Lam et al indicated that EBV DNA methylation profiles exhibited a disease-associated pattern [
      • Lam W.K.J.
      • Jiang P.
      • Chan K.C.A.
      • Peng W.
      • Shang H.
      • Heung M.M.S.
      • et al.
      Methylation analysis of plasma DNA informs etiologies of Epstein-Barr virus-associated diseases.
      ]. Nevertheless, there is still a lack of clinical trials focusing on the prediction efficiency of these biomarkers in NPC patients.
      High-throughput gene expression profiling techniques (e.g. microarray and high-throughput gene sequencing) provide new bioinformatic methods for identification of genes and pathological types that are associated with the onset of certain diseases. Meanwhile, they promote the identification of gene signatures related to the disease prognostic prediction and treatment. To our best knowledge, Weighted Gene Co-Expression Network Analysis (WGCNA) approach is proposed to provide a functional interpretation that is biologically significant and efficient in constructing global network structures [
      • Peña-Castillo L.
      • Mercer R.G.
      • Gurinovich A.
      • Callister S.J.
      • Wright A.T.
      • Westbye A.B.
      • et al.
      Gene co-expression network analysis in Rhodobacter capsulatus and application to comparative expression analysis of Rhodobacter sphaeroides.
      ,
      • Xu P.
      • Yang J.
      • Liu J.
      • Yang X.
      • Liao J.
      • Yuan F.
      • et al.
      Identification of glioblastoma gene prognosis modules based on weighted gene co-expression network analysis.
      ]. Therefore, it is considered as the most straightforward method for investigating the gene co-expression networks. In this study, WGCNA was utilized to identify the molecular markers associated with NPC prognosis.

      2. Materials and methods

      2.1 Tissues samples

      A total of 35 NPC patients admitted to our department between January 2014 and December 2015 were included in this study. All NPC patients had signed their written informed consent to participate this study. Clinical stages of NPC patients were classified by the TNM system of the American Joint Committee on Cancer. The patients with other malignancy and incomplete clinical or prognostic information were excluded from this study. The tissues samples were collected within 15 min after sample collection. The study protocols were approved by the Institutional Review Board of our institution.

      2.2 GEO download and pre-processing

      GEO query was utilized to download two NPC mRNA microarray datasets including GSE12452 and GSE53819 from GEO. The datasets were shown in Table 1. The study flow chart was added in Fig. 1. GSE53819 [
      • Bao Y.N.
      • Cao X.
      • Luo D.H.
      • Sun R.
      • Peng L.X.
      • Wang L.
      • et al.
      Urokinase-type plasminogen activator receptor signaling is critical in nasopharyngeal carcinoma cell growth and metastasis.
      ] dataset included 18 NPC tissue samples and 18 non-cancer tissue samples. One third of population was female. The samples were obtained from each subject before treatment. Agilent-014850 human genome microarray (4 × 44K G4112F, Agilent Tech, CA, USA) served as the sequencing platform. GSE12452 [
      • Dodd L.E.
      • Sengupta S.
      • Chen I.H.
      • den Boon J.A.
      • Cheng Y.J.
      • Westra W.
      • et al.
      Genes involved in DNA repair and nitrosamine metabolism and those located on chromosome 14q32 are dysregulated in nasopharyngeal carcinoma.
      ] was based on 10 normal nasopharyngeal tissues and 31 NPC tissues obtained before therapy. The sequencing platform was U133 Plus 2.0 array (Affymetrix, CA, USA). For the pre-processing of mRNA expression spectrum datasets, probes were projected to corresponding gene symbols. In cases of multiple probes on the same gene, expression level was calculated as the median, together with screening of mRNA dataset.
      Table 1Datasets in the NPC expressing spectrum of GEO.
      TypeGEO accessionPlatformTumorNormal
      mRNAGSE12452GPL5703110
      mRNAGSE53819GPL64801818

      2.3 Differentially expressed genes(DEGs)

      Limma package [
      • Ritchie M.E.
      • Phipson B.
      • Wu D.
      • Hu Y.
      • Law C.W.
      • Shi W.
      • et al.
      limma powers differential expression analyses for RNA-sequencing and microarray studies.
      ] was utilized for the differential analysis of GSE12452 and GSE53819. The screening for the mRNA expression spectrum was based on Benjamini-Hochberg adjusted P value of less than 0.05.

      2.4 WGCNA

      For the WGCNA, the gene expression was used to establish the similar matrix based on the absolute value of the Pearson correlation coefficient of both genes. The Pearson correlation coefficient was calculated according to the following formula:
      Sij=|1+cor(xi,yj)2|


      where i and j represented two genes, and xi and xj represented the expression values, respectively. Then the similar matrix was transmitted to the adjacency matrix (signed type) using the following formula:
      aij|Sij|β


      where i and j represented two genes and β represented the soft threshold.
      Afterwards, the adjacency matrix was transmitted to the topological matrix, and the topological overlap measure (TOM) was utilized to describe the correlation between two genes:
      TOM=uijaiuauj+aijmin(uaiu+uaju)+1aij


      where 1-TOM represented the variation between genei and genej.
      The 1-TOM served as the distance for the clustering analysis of the genes, and then the module recognition was performed using the dynamic splicing trees. The most representative gene in each module was the module eigengene (ME), which reflected the whole level of the gene expression in the modules. It was the first main constituent of each module, and the ME was calculated using the following formula:
      ME=princomp(xil(q))


      where i represented the gene in the module q, and 1 represented the chip sample in module q.
      The module membership (MM) was evaluated based on the expression spectrum of a certain gene in all the samples and the ME, using the formula as follows:
      MMiq=cor(xi,MEq)


      Where xi represented the expression spectrum of genei, MEq represented the ME of the module q. MMqi represented the membership of genei in the module q.
      In cases of MMqi=0, the genei was not in the module q; in the presence of a MMqi of approximately +1 or -1, there was a highly positive or negative correlation between genei and module q.
      The WGCNA package [
      • Pei G.
      • Chen L.
      • Zhang W.
      WGCNA Application to Proteomic and Metabolomic Data Analysis.
      ] was utilized to establish weighted gene co-expression network according to expression profile of the candidate genes. Then NPC-related modules were screened, and the interacted sub-networks were established using the co-expression data sets in the modules.

      2.5 Functional analysis of modules

      Candidate hub genes of specific target modules were subjected to KEGG pathway and GO functional enrichment analysis. KEGG pathways were analyzed using KEGG database (http://www.KEGG.jp) and GO functional enrichment analysis was carried out using GO database (http://geneontology.org/).

      2.6 Prognostic analysis based on hub genes

      For prognostic analysis, NPC dataset GSE102349 in GEO database was utilized. Disease progression and no progression in the last follow-up were considered as outcome of events. Median was used as boundary for evaluating high or low expression of genes. In the follow-up, Kaplan-Meier Plotter was utilized to depict progression-free survival (PFS) curvature. There were 113 samples in the GSE102349 dataset. The sequencing platform was GPL11154 Illumina HiSeq 2000 (Homo sapiens). After excluding the 25 samples involving data without event prognosis, 72 samples with no progression and 16 samples with disease progression were identified in the last follow-up.

      2.7 qRT-PCR verification

      Total RNA of tissues were extracted by TRIzol reagent (Invitrogen, Carlsbad, CA, USA). Approximately 1.2 μg RNA was reversely transcribed to cDNA by using Fast Quant RT Kit (TaKaRa, Dalian, China). Subsequently, qRT-PCR was carried out by SYBR Premix Ex Taq II kit (TaKaRa) and Roche Light Cycler® 480 System (Roche, Chicago, USA). Relative expression change of targets was analyzed by 2−ΔΔCt method [
      • Livak K.J.
      • Schmittgen TD.
      Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method.
      ], with GAPDH as the internal reference.

      2.8 Statistical analysis

      Statistical analysis was performed using GraphPad Prism v.7.0 (GraphPad Software Inc, San Diego, CA, USA). A P value of <0.05 was considered to be statistically significant.

      2.9 Patient and public involvement

      Patients and the public were not involved in the design, or conduct, or reporting, or dissemination plans of this research.

      3. Results

      3.1 DEGs screening

      Two mRNA microarray data of NPC including GSE12452 and GSE53819 were downloaded from GEO. The screening standard for mRNA expression spectrum was Benjamini-Hochberg adjusted P value of less than 0.05 and |logFC|>1). A total of 739 DEGs were identified in Supplementary File 1, including 457 that were up-regulated and 282 that were down-regulated in GSE12452 (Fig. 2a). For GSE53819, 1,756 DEGs (Supplementary File 2) were identified including 694 that were up-regulated and 1,062 that were down-regulated (Fig. 2b).
      Fig. 2
      Fig. 2Comparison of genes differentially expressed between GSE12452 (a) and GSE53819 (b) (volcano plot).

      3.2 Establishing of WGCNA based on DEGs

      WGCNA package [
      • Langfelder P.
      • Horvath S.
      WGCNA: an R package for weighted correlation network analysis.
      ] was utilized to establish WGCNA for candidate gene set (GSE53819). Co-expression network was in line with scale network, which demonstrated a negative correlation between the log(k) and log(P(k)). The correlation coefficient was larger than 0.9. To ensure scale-free network, we selected optimal β value (β=6) (Fig. 3a). Subsequently, the expression matrix was transmitted to the adjacency matrix, while adjacency matrix was transmitted to topological matrix. Based on TOM, average linkage clustering algorithm was used for clustering analysis of genes. Gene network module should contain at least 30 genes. Samples underwent analysis showed satisfactory consensus with no outlier. Therefore, the process for the outlier was not deleted.
      Fig. 3
      Fig. 3Correlation analyses of the modules. (a). Confirmation of soft threshold of β. There was a negative correlation between log(k) and log(P(k)), and the correlation coefficient was >0.9. (b). Co-expression module. There were 12 co-expression modules which were represented in different colors. These modules ranked in an order of the number of genes in each module. (c). Correlation between module and traits; the modules were represented by different abscissa and ordinates of various colors. Each rectangle displayed the correlation value of each module. (d). Co-expression TOM. Each module was indicated by different abscissa and ordinates of various colors. The brightness in the middle zone of a yellow color indicated the link between different modules.
      Adjacency modules were merged into new modules with a height set at 0.25, and then 12 modules were obtained (Fig. 3b). Gene list for each module was provided in Supplementary File 3. Correlation coefficient for the Brown, Magenta, and Turquoise modules was -0.87, 0.69, and -0.68, respectively. Both the positively and negatively highest absolute value for the correlation coefficients, namely Brown and Magenta, were analyzed. As shown in Fig. 3c, Brown module and Magenta module were extremely associated with NPC samples. Gene network with all genes were shown in Fig. 3d.

      3.3 Functional analysis of modules

      GO functional analysis was carried out to genes in the Brown module (Supplementary File 4) and Magenta module (Supplementary File 5). Brown module was associated with positive regulation of mitotic nuclear division, positive regulation of nuclear division, protein localization to microtubule cytoskeleton, endoplasmic reticulum to Golgi vesicle-mediated transport, cilium or flagellum dependent cell motility, cilium dependent cell motility, protein localization to cytoskeleton, protein localization to cell surface as well as cellular oxidant detoxification. Brown module and Magenta module were correlated with biological regulation, metabolic process, reproduction, as well as cellular proliferation (Fig. 4a and 4c).
      Fig. 4
      Fig. 4GO functional annotation and KEGG enrichment analysis of the Brown module and Magenta module. (a, c). GO functional annotation of the Brown module and Magenta module. (b, d). KEGG pathway enrichment analysis of Brown module and Magenta module.
      KEGG pathway analysis was given to genes in Brown module (Supplementary File 6) and Magenta module (Supplementary File 7). For Brown module, KEGG pathway was mainly enriched in systemic lupus erythematosus, alcoholism, glutathione metabolism, tuberculosis, drug metabolism-cytochrome P450, metabolism of xenobiotics by cytochrome P450, apoptosis-multiple species, Huntington disease, and viral carcinogenesis. For Magenta module, KEGG pathway was mainly enriched in breast cancer, basal cell carcinoma, gastric cancer, arrhythmogenic right ventricular cardiomyopathy (ARVC), Hippo signaling pathway, Cushing syndrome, Wnt signaling pathway, Hepatocellular carcinoma, as well as thyroid cancer (Fig. 4b and 4d).
      For Brown module and Magenta module, hub genes were further analyzed to identify their association with pathogenesis of cancer. Co-expression network was finally established based on selection of co-expression pairs with a weight of >0.4 in Brown module and a weight of >0.5 in Magenta module. Genes in two modules were subject to network visualization, followed by calculating degree of each gene node (Supplementary File 8). The genes with a screening degree ranked 10% among all the degrees were selected as the 26 hub genes for the modules, including FAM47E, C22orf15, WDR38, TSPAN1, AK7, IL33, RALA, DNAH12, PDSS1, TMEM190, EFCAB1, AKAP14, CLEC5A, GTF2H3, KIF6, TNFAIP6, OLR1, CELSR2, STAP2, RTKN, SPINT1, RARG, MPP3, SLC16A7, EPC1 and HOXC8.

      3.4 Prognostic analysis based on hub genes

      GSE102349 dataset was devoted to identifying prognostic hub genes by survival analysis. As shown in Fig. 5, PFS was significantly correlated with expression of IL33, MPP3 and SLC16A7. These indicated that these hub genes may be closely related to prognosis of NPC patients. Next, quantitative identification of hub genes by real-time fluorescence. Expression trend of SLC16A7 was strongly correlated with overall survival (OS, P=0.0366, Fig. 6).
      Fig. 5
      Fig. 5Kaplan–Meier survival plots of L33 (a), MPP3 (b) and SLC16A7 (c) indicated that the expression of these genes was significantly associated with the PFS for NPC patients.
      Fig. 6
      Fig. 6Association between OS and IL33 (a), MPP3 (b) and SLC16A7 (c). OS was significantly associated with SLC16A7 (P=0.0366).

      4. Discussion

      NPC is characterized by distinct geographical distribution and is particularly prevalent in east and southeast Asia [
      • Chen Y.P.
      • Chan A.T.C.
      • Le Q.T.
      • Blanchard P.
      • Sun Y.
      • Ma J.
      Nasopharyngeal carcinoma.
      ]. Nowadays, radiotherapy is the major treatment option for NPC. Meanwhile, concurrent cisplatin-based chemoradiotherapy is considered as the care standard for locally advanced NPC [
      • Liao W.
      • Huang J.
      • Wu Q.
      • Zhu G.
      • Wang X.
      • Wen F.
      • et al.
      Concurrent chemoradiotherapy with nedaplatin versus cisplatin in stage II-IVB nasopharyngeal carcinoma: A cost-effectiveness analysis.
      ,
      • Tang L.Q.
      • Chen D.P.
      • Guo L.
      • Mo H.Y.
      • Huang Y.
      • Guo S.S.
      • et al.
      Concurrent chemoradiotherapy with nedaplatin versus cisplatin in stage II-IVB nasopharyngeal carcinoma: an open-label, non-inferiority, randomised phase 3 trial.
      ]. To date, the prognosis of NPC is highly relied on disease stage. Generally, patients with localized lesion and/or early-stage NPC (i.e. stage I and II) were reported to present good prognosis, with a longer five-year survival than the other counterparts [
      • Zhou J.Y.
      • Chong V.F.
      • Khoo J.B.
      • Chan K.L.
      • Huang J.
      The relationship between nasopharyngeal carcinoma tumor volume and TNM T-classification: a quantitative analysis.
      ]. Thus, early screening and diagnosis of NPC could improve treatment prognosis.
      Recently, screening of diagnostic and prognostic genes among NPC patients play an important role in predicting unknown pathological features [
      • Dai Y.
      • Zhang Y.
      • Yang M.
      • Zhou L.
      • Pan H.
      • Xiao T.
      • et al.
      Radiosensitivity-Related Genes and Clinical Characteristics of Nasopharyngeal Carcinoma.
      ,
      • Zhang F.
      • Yin SC.
      Expression and clinical significance of Testin in nasopharyngeal carcinoma.
      ]. To date, there is indeed a harsh demand for the screening technique, however, few have been developed due to technical challenges. This requests investigation on development of specific biomarkers for diagnosis and screening of NPC. In this study, we aim to identify molecular markers associated with pathogenesis of NPC.
      WGCNA, an algorithm for distinguishing module information from chip, is a powerful tool for identifying possible regulatory mechanisms, which provides insights into possible pathways of certain diseases [
      • Li W.
      • Wang L.
      • Wu Y.
      • Yuan Z.
      • Zhou J.
      Weighted gene co‑expression network analysis to identify key modules and hub genes associated with atrial fibrillation.
      ,
      • Song Z.Y.
      • Chao F.
      • Zhuo Z.
      • Ma Z.
      • Li W.
      • Chen G.
      Identification of hub genes in prostate cancer using robust rank aggregation and weighted gene co-expression network analysis.
      ,
      • Yao Q.
      • Song Z.
      • Wang B.
      • Qin Q.
      • Zhang JA.
      Identifying Key Genes and Functionally Enriched Pathways in Sjögren's Syndrome by Weighted Gene Co-Expression Network Analysis.
      ]. For the genes always with similar expression patterns in different tissues or in a physiological process, they were defined as a module. On this basis, candidate hub genes and/or treatment targets were selected using module connectivity to relationship between phenotype and module. Unlike conventional gene screening techniques only focusing on differential expression, WGCNA involves near 10,000 of variable genes to screen modules of interest [
      • Pflieger L.T.
      • Dansithong W.
      • Paul S.
      • Scoles D.R.
      • Figueroa K.P.
      • Meera P.
      • et al.
      Gene co-expression network analysis for identifying modules and functionally enriched pathways in SCA2.
      ]. In this study, WGCNA was conducted to explore gene expression modules related to NPC. Our results indicated that 2,110 genes were constructed by co-expression network and finally 12 modules were generated, in which Brown and Magenta modules were extremely associated with NPC samples.
      KEGG pathway analysis contributes to screening of certain diseases based on gene structures and classifies pathways into 9 categories based on their targets. In this study, we screened candidate genes with biological functions in KEGG pathways and GO terms to investigate mechanisms of existing genes associated with NPC pathogenesis and prognosis, and mRNA interaction and future prediction of prognostic factors. Genes showed different enrichment levels in GO terms and KEGG pathways, implying diversity in biological function enrichment. In this study, KEGG pathway analysis was given to the genes in the Brown module and Magenta module, which indicated that Brown module and Magenta module were correlated with biological regulation, metabolic process, reproduction, and cellular proliferation.
      For Brown module and Magenta module, hub genes of the two modules were further analyzed to identify hub genes associated with pathogenesis of cancer. Twenty-six hub genes were obtained and were considered to be closely related to NPC. Then we verified the prognostic value of these biomarkers based on prognostic information data (GSE102349 dataset), which enhanced our understanding on NPC prognosis. PFS was significantly correlated with expression of IL33, MPP3 and SLC16A7. These indicated that these hub genes may be closely related to prognosis of NPC patients. Next, quantitative identification of hub genes by real-time fluorescence. The expression of SLC16A7 was strongly correlated with OS in NPC patients.
      SLC16A7 gene encoded protein is a member of monocarboxylate transporter family, which catalyzes proton-linked transport of monocarboxylates and has highest affinity for pyruvate. To our best knowledge, such protein was reported to be more highly expressed in prostate and colorectal cancer specimens. In a previous study, epigenetic and oncogenic regulation of SLC16A7 was associated with over-expression of SLC16A7 protein, which impacted in modulating signaling pathway and cellular phenotypes in prostate cancer patients [
      • Pertega-Gomes N.
      • Vizcaino J.R.
      • Felisbino S.
      • Warren A.Y.
      • Shaw G.
      • Kay J.
      • et al.
      Epigenetic and oncogenic regulation of SLC16A7 (MCT2) results in protein over-expression, impacting on signalling and cellular phenotypes in prostate cancer.
      ]. Meanwhile, in a study designed to assess effects of single-nucleotide polymorphisms of SLC16A7 gene on prognosis of colorectal cancer, SLC16A7 gene polymorphism may affect clinical outcome of patients, and could be used to predict responses among patients underwent adjuvant chemotherapy [
      • Fei F.
      • Guo X.
      • Chen Y.
      • Liu X.
      • Tu J.
      • Xing J.
      • et al.
      Polymorphisms of monocarboxylate transporter genes are associated with clinical outcomes in patients with colorectal cancer.
      ]. In this study, we demonstrated that SLC16A7 expression was associated with prognosis of NPC patients, which needed further approval by clinical trials.
      In this study, based on the analysis of WGCNA, we demonstrated that Brown module, Magenta module and 26 hub genes might be closely associated with the carcinogenesis of NPC, among which the expression trends of SLC16A7 were strongly correlated with OS in NPC. In the future, further studies are needed to verify the functions of the candidate core gene SLC16A7 in different level using methods such as western blot, and transgenic methods, in order to gain a deeper understanding of the relationship between the genes and modules identified in this study.

      Disclosure Statement

      The authors declare that there is no conflict of interest.

      Declaration of Competing Interest

      The authors declare that there is no conflict of interest.

      Acknowledgments

      Not applicable.

      Funding

      This work was supported by the Funding Project of Fujian Medical University College Student Innovation and Entrepreneurship Training Program [Grant No. C19071].

      Ethics approval and informed consent

      The research related to human use has been complied with all the relevant national regulations, institutional policies and in accordance the tenets of the Helsinki Declaration, and has been approved by the Ethics Committee of Fujian Provincial Hospital (K2021-04-039). All participants gave informed consent before taking part in the study.

      Consent for participate

      Patients and the public were not involved in the design, or conduct, or reporting, or dissemination plans of this research.

      Appendix. Supplementary materials

      References

        • Sung H.
        • Ferlay J.
        • Siegel R.L.
        • Laversanne M.
        • Soerjomataram I.
        • Jemal A.
        • et al.
        Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries.
        CA Cancer J Clin. 2021;
        • Lin CT.
        Relationship between Epstein-Barr virus infection and nasopharyngeal carcinoma pathogenesis.
        Ai Zheng. 2009; 28: 791-804
        • Stepan K.O.
        • Mazul A.L.
        • Skillington S.A.
        • Paniello R.C.
        • Rich J.T.
        • Zevallos J.P.
        • et al.
        The prognostic significance of race in nasopharyngeal carcinoma by histological subtype.
        Head Neck. 2021;
        • Zhang F.
        • Zhong L.Z.
        • Zhao X.
        • Dong D.
        • Yao J.J.
        • Wang S.Y.
        • et al.
        A deep-learning-based prognostic nomogram integrating microscopic digital pathology and macroscopic magnetic resonance images in nasopharyngeal carcinoma: a multi-cohort study.
        Ther Adv Med Oncol. 2020; 121758835920971416
        • Lyu Y.
        • Ni M.
        • Zhai R.
        • Kong F.
        • Du C.
        • Hu C.
        • et al.
        Clinical characteristics and prognosis of elderly nasopharyngeal carcinoma patients receiving intensity-modulated radiotherapy.
        Eur Arch Otorhinolaryngol. 2020;
        • Sengar M.
        • Chorghe S.
        • Jadhav K.
        • Singh S.
        • Laskar S.G.
        • Pai P.
        • et al.
        Cell-free Epstein-Barr virus-DNA in patients with nasopharyngeal carcinoma: Plasma versus urine.
        Head Neck. 2016; 38: E1666-E1673
        • Wang H.
        • Liu W.
        • Luo B.
        The roles of miRNAs and lncRNAs in Epstein-Barr virus associated epithelial cell tumors.
        Virus Res. 2021; 291198217
        • Lv J.
        • Chen Y.
        • Zhou G.
        • Qi Z.
        • Tan K.R.L.
        • Wang H.
        • et al.
        Liquid biopsy tracking during sequential chemo-radiotherapy identifies distinct prognostic phenotypes in nasopharyngeal carcinoma.
        Nat Commun. 2019; 10: 3941
        • Lin J.C.
        • Wang W.Y.
        • Chen K.Y.
        • Wei Y.H.
        • Liang W.M.
        • Jan J.S.
        • et al.
        Quantification of plasma Epstein-Barr virus DNA in patients with advanced nasopharyngeal carcinoma.
        N Engl J Med. 2004; 350: 2461-2470
        • Huang F.
        • Liang X.
        • Min X.
        • Zhang Y.
        • Wang G.
        • Peng Z.
        • et al.
        Simultaneous Inhibition of EGFR and HER2 via Afatinib Augments the Radiosensitivity of Nasopharyngeal Carcinoma Cells.
        J Cancer. 2019; 10: 2063-2073
        • Sun H.
        • Liu M.
        • Wu X.
        • Yang C.
        • Zhang Y.
        • Xu Z.
        • et al.
        Overexpression of N-cadherin and β-catenin correlates with poor prognosis in patients with nasopharyngeal carcinoma.
        Oncol Lett. 2017; 13: 1725-1730
        • Kim T.J.
        • Lee Y.S.
        • Kang J.H.
        • Kim Y.S.
        • Kang CS.
        Prognostic significance of expression of VEGF and Cox-2 in nasopharyngeal carcinoma and its association with expression of C-erbB2 and EGFR.
        J Surg Oncol. 2011; 103: 46-52
        • Cheng Y.
        • Wu S.
        • Tie X.
        • Huang X.
        • Cui L.
        Growth Inhibition of Nasopharyngeal Carcinoma Cells Mediated by p53 Gene-Containing Nanolipid Composites.
        J Nanosci Nanotechnol. 2020; 20: 6026-6032
        • Savitri E.
        • Haryana MS.
        Expression of interleukin-8, interleukin-10 and Epstein-Barr viral-load as prognostic indicator in nasopharyngeal carcinoma.
        Glob J Health Sci. 2015; 7: 364-372
        • Lam W.K.J.
        • Jiang P.
        • Chan K.C.A.
        • Peng W.
        • Shang H.
        • Heung M.M.S.
        • et al.
        Methylation analysis of plasma DNA informs etiologies of Epstein-Barr virus-associated diseases.
        Nat Commun. 2019; 10: 3256
        • Peña-Castillo L.
        • Mercer R.G.
        • Gurinovich A.
        • Callister S.J.
        • Wright A.T.
        • Westbye A.B.
        • et al.
        Gene co-expression network analysis in Rhodobacter capsulatus and application to comparative expression analysis of Rhodobacter sphaeroides.
        BMC Genomics. 2014; 15: 730
        • Xu P.
        • Yang J.
        • Liu J.
        • Yang X.
        • Liao J.
        • Yuan F.
        • et al.
        Identification of glioblastoma gene prognosis modules based on weighted gene co-expression network analysis.
        BMC Med Genomics. 2018; 11: 96
        • Bao Y.N.
        • Cao X.
        • Luo D.H.
        • Sun R.
        • Peng L.X.
        • Wang L.
        • et al.
        Urokinase-type plasminogen activator receptor signaling is critical in nasopharyngeal carcinoma cell growth and metastasis.
        Cell Cycle. 2014; 13: 1958-1969
        • Dodd L.E.
        • Sengupta S.
        • Chen I.H.
        • den Boon J.A.
        • Cheng Y.J.
        • Westra W.
        • et al.
        Genes involved in DNA repair and nitrosamine metabolism and those located on chromosome 14q32 are dysregulated in nasopharyngeal carcinoma.
        Cancer Epidemiol Biomarkers Prev. 2006; 15: 2216-2225
        • Ritchie M.E.
        • Phipson B.
        • Wu D.
        • Hu Y.
        • Law C.W.
        • Shi W.
        • et al.
        limma powers differential expression analyses for RNA-sequencing and microarray studies.
        Nucleic Acids Res. 2015; 43: e47
        • Pei G.
        • Chen L.
        • Zhang W.
        WGCNA Application to Proteomic and Metabolomic Data Analysis.
        Methods Enzymol. 2017; 585: 135-158
        • Livak K.J.
        • Schmittgen TD.
        Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method.
        Methods. 2001; 25: 402-408
        • Langfelder P.
        • Horvath S.
        WGCNA: an R package for weighted correlation network analysis.
        BMC Bioinformatics. 2008; 9: 559
        • Chen Y.P.
        • Chan A.T.C.
        • Le Q.T.
        • Blanchard P.
        • Sun Y.
        • Ma J.
        Nasopharyngeal carcinoma.
        Lancet. 2019; 394: 64-80
        • Liao W.
        • Huang J.
        • Wu Q.
        • Zhu G.
        • Wang X.
        • Wen F.
        • et al.
        Concurrent chemoradiotherapy with nedaplatin versus cisplatin in stage II-IVB nasopharyngeal carcinoma: A cost-effectiveness analysis.
        Oral Oncol. 2019; 93: 15-20
        • Tang L.Q.
        • Chen D.P.
        • Guo L.
        • Mo H.Y.
        • Huang Y.
        • Guo S.S.
        • et al.
        Concurrent chemoradiotherapy with nedaplatin versus cisplatin in stage II-IVB nasopharyngeal carcinoma: an open-label, non-inferiority, randomised phase 3 trial.
        Lancet Oncol. 2018; 19: 461-473
        • Zhou J.Y.
        • Chong V.F.
        • Khoo J.B.
        • Chan K.L.
        • Huang J.
        The relationship between nasopharyngeal carcinoma tumor volume and TNM T-classification: a quantitative analysis.
        Eur Arch Otorhinolaryngol. 2007; 264: 169-174
        • Dai Y.
        • Zhang Y.
        • Yang M.
        • Zhou L.
        • Pan H.
        • Xiao T.
        • et al.
        Radiosensitivity-Related Genes and Clinical Characteristics of Nasopharyngeal Carcinoma.
        Biomed Res Int. 2020; 20201705867
        • Zhang F.
        • Yin SC.
        Expression and clinical significance of Testin in nasopharyngeal carcinoma.
        Lin Chung Er Bi Yan Hou Tou Jing Wai Ke Za Zhi. 2016; 30: 982-985
        • Li W.
        • Wang L.
        • Wu Y.
        • Yuan Z.
        • Zhou J.
        Weighted gene co‑expression network analysis to identify key modules and hub genes associated with atrial fibrillation.
        Int J Mol Med. 2020; 45: 401-416
        • Song Z.Y.
        • Chao F.
        • Zhuo Z.
        • Ma Z.
        • Li W.
        • Chen G.
        Identification of hub genes in prostate cancer using robust rank aggregation and weighted gene co-expression network analysis.
        Aging (Albany NY). 2019; 11: 4736-4756
        • Yao Q.
        • Song Z.
        • Wang B.
        • Qin Q.
        • Zhang JA.
        Identifying Key Genes and Functionally Enriched Pathways in Sjögren's Syndrome by Weighted Gene Co-Expression Network Analysis.
        Front Genet. 2019; 10: 1142
        • Pflieger L.T.
        • Dansithong W.
        • Paul S.
        • Scoles D.R.
        • Figueroa K.P.
        • Meera P.
        • et al.
        Gene co-expression network analysis for identifying modules and functionally enriched pathways in SCA2.
        Hum Mol Genet. 2017; 26: 3069-3080
        • Pertega-Gomes N.
        • Vizcaino J.R.
        • Felisbino S.
        • Warren A.Y.
        • Shaw G.
        • Kay J.
        • et al.
        Epigenetic and oncogenic regulation of SLC16A7 (MCT2) results in protein over-expression, impacting on signalling and cellular phenotypes in prostate cancer.
        Oncotarget. 2015; 6: 21675-21684
        • Fei F.
        • Guo X.
        • Chen Y.
        • Liu X.
        • Tu J.
        • Xing J.
        • et al.
        Polymorphisms of monocarboxylate transporter genes are associated with clinical outcomes in patients with colorectal cancer.
        J Cancer Res Clin Oncol. 2015; 141: 1095-1102