Skip to main content

ORIGINAL RESEARCH article

Front. Immunol., 12 October 2022
Sec. Cancer Immunity and Immunotherapy

Identification of hub genes and their correlation with infiltration of immune cells in MYCN positive neuroblastoma based on WGCNA and LASSO algorithm

Ji Chen&#x;Ji Chen1†Mengjiao Sun&#x;Mengjiao Sun2†Chuqin Chen&#x;Chuqin Chen2†Bin Jiang*Bin Jiang1*Yongjun Fang*Yongjun Fang2*
  • 1Department of General Surgery, Children’s Hospital of Nanjing Medical University, Nanjing, China
  • 2Department of Hematology and Oncology, Children’s Hospital of Nanjing Medical University, Nanjing, China

Background: The prognosis of MYCN positive NB is poor, and there is no targeted drug for N-myc at present. This study aims to screen out hub genes closely related to MYCN, analyze the relationship between hub genes and NB microenvironment, and provide basis for molecular targeted therapy of MYCN positive NB.

Methods: We combined the microarray data of GSE45547 (n=649) and GSE49710 (n=498), screened the DEGs between MYCN positive (n=185) and MYCN negative NB (n=951), performed WGCNA, Lasso regression and Roc analyses on the merged matrix, and obtained the hub genes related to MYCN in the training group. We performed ssGSEA on the experimental group to calculate the infiltration level of 28 kinds of immune cells in each sample, compared the differences of immune cell infiltration between MYCN positive and MYCN negative group. The influences of hub genes on the distribution of each immune cell were also analyzed by ssGSEA. The expression differences of the three hub genes were verified in the E-MTAB-8248 cohort (n=223), and the correlation between hub genes and prognosis of NB was calculated by Kaplan-Meier method in GSE62564 (n=498) and the validation group. We also verified the expression differences of hub genes by qRT-PCR in SK-N-BE(2), SKNDZ, Kelly and SH-SY5Y cell lines.

Results: Here were 880 DEGs including 420 upregulated and 460 downregulated genes in MYCN positive NB in the training group. Overlap of the DEGs and WGCNA networks identified four shared genes, namely, ZNF695, CHEK1, C15ORF42 and EXO1, as candidate hub genes in MYCN positive NB. Three core genes, ZNF695, CHEK1 and C15ORF42, were finally identified by Lasso regression and Roc analyses. ZNF695, CHEK1 and C15ORF42 were highly expressed in MYCN positive NB tissues and cell lines. These three genes were closely related to the prognosis of children with NB. Except that Activated CD4 T cell and Type2 T helper cell increased, the infiltration levels of the other 26 cells decreased significantly in MYCN positive NB tissues. The infiltration levels of Type2 T helper cell and Activated CD4 T cell were also significantly positively correlated with the expression levels of the three hub genes.

Conclusion: ZNF695, CHEK1 and C15ORF42 are highly expressed in MYCN positive NB, and their expression levels are negatively correlated with the prognosis of children with NB. The infiltration levels of Activated CD4 T cell and Type2 T helper cell increased in the microenvironment of MYCN positive NB and were significantly positively correlated with the expression levels of the three hub genes. The results of this study provide that ZNF695, CHEK1 and C15ORF42 may be potential prognostic markers and immunotherapy targets for MYCN positive NB.

Introduction

Neuroblastoma (NB) is one of the most common solid extracranial malignancies in infants (1). Surgery combined with chemotherapy can achieve a good prognosis for children with low-risk (2) or intermediate-risk NB (3, 4). For high-risk NB, although combined surgery (5, 6), chemotherapy (5, 7), radiotherapy (5), immunotherapy (811), bone marrow transplantation (5) and other treatments, the prognoses of these patients remain poor (8). According to the revised NB risk classification system promulgated by Children’s Oncology Group in 2021, MYCN amplification is the most important biological factor affecting the risk staging of NB (12).

As an important member of MYC oncogene family, MYCN has highly conserved amino acid sequences with other members of MYC family, and these highly conserved amino acid sequences determine that members of MYC family have many identical biological properties (13). As an important transcription factor, N-myc protein formed by MYCN translation promotes the transcription of many genes related to NB metabolism (14), immune escape (15), apoptosis (16), ferroptosis (17) and other fields. Especially in the aspect of immune surveillance, recent studies have found that MYC family regulates the expression and production of a variety of immune ligands, receptors and immune effector molecules (1820), thus affecting the biological activities of CD4+T (21), Macrophages (22), NK cells (23), B cells (24) and other cells, exerting an important influence on the microenvironment in NB (15). Although there have been attempts to use small molecules or low molecular weight compounds such as NY2267 (25), IIA6B17 (26), Tz-1 (27) and 10058-F4 (28) to target MYC family, no clinically available N-myc inhibitors can accomplish this task, and N-myc has long been thought to be an “untreatable” proto-oncoprotein (29). It has forced researchers to continue to innovate in targeted therapies for MYCN-driven malignancies.

Weighted gene Co-expression Network analysis (WGCNA) is a bioinformatics analysis method that is often used to effectively explore the relationship between genes and phenotypes (30). The significant advantage of WGCNA lies in its ability to aggregate genes into co-expression modules, bridging the gap between sample characteristics and gene expression changes. WGCNA analyzes thousands of genes to identify gene modules associated with clinical features, identifies key genes in disease pathways for further validation, and ultimately provides a system-level insight into signal networks that may be associated with a phenotype of interest.

In this study, differentially expressed genes (DEGs) of MYCN positive NB were screened out by GEO and ArrayExpress databases, and gene co-expression network was constructed by WGCNA, Lasso regression and Roc analyses to obtain modules and hub genes that play important biological roles in MYCN positive NB. The immunocyte correlation analyses of MYCN positive NB, MYCN negative NB and hub genes were also carried out. Figure 1 illustrates the workflow chart of data preparation, processing, analysis, and validation. The views provided in this paper will provide new insights into the immunotherapy of MYCN positive NB.

FIGURE 1
www.frontiersin.org

Figure 1 Flow chart of this study.

Methods

Cell culture

The NB cell lines used for the experiments including SK-N-BE (2), SKNDZ and SH-SY5Y were obtained from American Type Culture Collection (ATCC). Kelly cell line was kindly provided by Dr Guoliang Qing from Department of Pathophysiology, School of Basic Medical Sciences, Wuhan University. MYCN is highly expressed in SK-N-BE (2), SKNDZ and Kelly cell lines and low expressed in SH-SY5Y. All four cell lines are representative, widely used NB cell lines with relatively high feasibility. SK-N-BE (2), SKNDZ and SH-SY5Y cell lines were routinely cultured with DMEM (C11995500BT, Gibco, USA) supplemented with 10% fetal bovine serum (10099-141C, Gibco, USA) and1% Penicillin/Streptomycin (C100C5, New cell & Molecular Biotech, China), Kelly cell line was grown in complete RPMI-1640 (C11875500BT, Gibco, USA) supplemented with 10% FBS and 1% penicillin/streptomycin. All short tandem repeat (STR) -certified cell lines were cultured at 37C and 5% CO2 for up to 6 months after resuscitation, with mycoplasma contamination detected regularly using MycoAlert (LT07-710, Lonza, Switzerland).

RNA extraction and real-time quantitative PCR (qRT-PCR)

Total RNAs were extracted from NB cells by Trizol reagent (15596018, ambion, USA). Then cDNA was synthesized using HiScript@IIQ RT SuperMix (R222-01, Vazyme, China) and qRT-PCR was performed with AceQ qPCR SYBR Green Master Mix (Low ROX Premixed)(Q131, Vazyme, China) according to the manufacturer’s instructions. The reaction conditions were: pre-denaturation at 95 C for 5 min, 40 cycles, denaturation at 95 C for 10 s, annealing at 60.0 C for 30 s, extension at 95C for 15 s, 60C for 60 s and 95 C for 15 s. qRT-qPCR was performed on a Lightcyler@ 96 instrument (05815916001, Roche, USA). Relative RNA amount was calculated by 2-ΔΔCt (31)with the normalization to β-actin. β-actin is a cytoskeletal protein encoded by housekeeping gene, and its expression is relatively constant in various tissues and cell lines, so it is a suitable reference gene. The primer sequences (General Biol, China) were listed in Table S1.

Data collection and DEGs analysis

Through data retrieval of GEO and ArrayExpress databases, GSE45547 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE45547), GSE49710 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE49710), E-MTAB-8248 (https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-8248/) and GSE62564 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE62564) cohorts were selected as the objects of this research. The same patients were enrolled but different platforms were used in GSE62564 and GSE49710. All cohorts met the following criteria: 1) large sample size, 2) complete clinical information and microarray data, 3) fresh NB tissues for microarray analysis. Because of the same platform, GPL16876, GSE45547 and GSE49710 used, we integrated the data of the two cohorts into a large, fused expression matrix, after removing batch effects in order to obtain as many samples as possible to improve the statistic power. Meanwhile, E-MTAB-8248 cohort was selected as the validation group to fully verify the reliability of the results. GSE62564 was only used for survival analysis. After excluding the cases with unclear MYCN expression, 643 cases were included in GSE45547, including 93 MYCN positive cases and 550 MYCN negative cases. GSE49710 included 493 cases, including 92 MYCN positive cases and 401 MYCN negative cases. A total of 222 patients were included in the E-MTAB-8248 cohort, including 46 MYCN positive cases and 176 MYCN negative cases. Log2 transformation of transcriptome data was performed and probes were matched with genetic symbols according to the annotated documentation of the corresponding platform. Finally, the gene matrix with row name as sample name and column name as gene symbol was obtained for subsequent analysis. Data from GSE45547 and GSE49710 was averaged and combined into genomicMatrix as training group. By comparing the gene expression profiles of MYCN positive NB with those of MYCN negative NB, a list of DEGs with |log2fc|> 1.2 and p-value< 0.05 was obtained in the combined cohort.

WGCNA

Data were preprocessed by “WGCNA” R package and abnormal genes were detected and excluded by the “GoodSamplegenes()” function. Then we used the “hclust()” function to draw the cluster graph, select the appropriate threshold, eliminate the outlier samples, and obtain relatively consistent gene expression data. We performed automatic network topology analysis according to the “pickSoftThreshold()” function, determined the soft threshold parameter (β) according to the principle of scale-free network and finally constructed a scale-free weighted gene co-expression network. Topological overlap measurement (ToM) was used to identify highly co-expressed gene modules to reduce the sensitivity of the network to false connections. By cutting the cluster tree into branches, the genes with high absolute correlation were clustered into the same module. Only modules with more than 60 genes could be defined as valid modules and the modules with high correlation would be merged (MEDissThres = 0.25). Each module would be assigned a different color for visual analysis. Cluster analysis of modules was based on the eigenvector value of each module. The “hclust()” and “merge()” functions were used to draw hierarchical clustering trees and merge modules with high correlation respectively. Finally, we plotted the scatter plot between Gene Significance (GS) and Module Membership (MM) within each module through the “plot()” function to understand the significance of high-connectivity genes within the module.

Gene ontology (GO) enrichment analyses

GO function analysis was divided into three parts: biological process (BP), cell component (CC) and molecular function (MF). We performed GO enrichment analyses of the DEGs by “clusterProfiler” R package and drew bubble graph of GO enrichment analysis through “GGplot” R package.

Gene set enrichment analysis

Gene set enrichment analysis (GSEA) was performed with “clusterProfiler” R package to calculate the gene sets with statistical differences between MYCN positive group and MYCN negative group. P-value < 0.05 was selected as the cutoff value. The annotated gene set entitled “immunesigdb” used in GSEA was downloaded from Molecular Signatures Database (MSigDB, https://www.gseamsigdb.org/).

Establishment and validation of lasso regression analyses

The intersecting genes were obtained by overlap of DEGs (|log2fc|> 2.0 and p-value < 0.05) and hub genes, and the Lasso regression analysis was constructed by “Glmnet()” function. The minimum regularization parameter lambda (λ) and genes highly associated with MYCN positive NB were obtained by cross-validation (alpha=1). We draw receiver operator characteristic curve (Roc) for each gene by running “pROC()” function, calculated area under the curve (AUC) and evaluated the degree of association between MYCN positive NB and the hub genes in the training group and validation group. Survival curves of core genes were plotted by running the “SURv()” function, and Kaplan-Meier analysis was conducted to estimate the correlation between hub genes and prognosis.

Analysis of immune cell infiltration

We performed single-sample gene set enrichment analysis (ssGSEA) on the training group using the “GSVA()” function to calculate the infiltration level of 28 kinds of immune cells in each sample, compared the difference of immune cell infiltration between MYCN positive and MYCN negative samples by Wilcoxon test and plotted the violin through the “Vioplot” R package. The correlation between 28 infiltrating immune cells and hub genes was further analyzed by “cor.test()” function, and the heatmap was drawn by “ggplot()” function.

Results

DEGs analysis

A total of 1136 cases and 19806 genes were included by filtering and integrating the microarray data of GSE45547 and GSE49710 cohort. We selected 880 DEGs, including 420 highly expressed genes and 460 low expressed genes in MYCN positive NB by analyzing the gene expression levels of the two groups (|log2fc| > 1.2 and p-value < 0.05) (Figures 2A, B). Among all DEGs, the 10 up-regulated genes with highest value of log2fc were SLCO4A1, LMO3, MGC16291, ABCA12, HOXD10, TWIST1, NPW, GABRA5, NMU and DUXAP10. The 10 down-regulated genes with minimal log2fc were NTRK1, KRT19, LOC388002, IL7, RGS9, CCL19, RP13-102H20.1, SLC18A2, ECEL1 and FAM70A.

FIGURE 2
www.frontiersin.org

Figure 2 Identification of DEGs (A) Heatmap of representative DEGs in MYCN positive and MYCN-negative NB tissues. (B) Volcano plot of all genes in NB tissues. The vertical dotted lines represent |log2fc|= 1.2, the horizontal dotted lines represent p-value= 0.05.

GO enrichment analyses

Three subontologies of DEGs including biological process (BP), cellular component (CC) and molecular function (MF) were analyzed by GO analysis. 639, 43 and 20 pathways were enriched in each subontology respectively (Figure 3A). For BP, DEGs were mainly enriched in T cell activation (GeneRatio: 54/661, adjust p-value: 4.08×10^-10), T cell differentiation (GeneRatio: 36/661, adjust p-value: 3.78×10^-9), positive regulation of leukocyte cell-cell adhesion (GeneRatio: 32/661, adjust p-value: 1.48×10^-7), positive regulation of T cell activation (GeneRatio: 30/661, adjust p-value:1.73×10^-7) and lymphocyte differentiation (GeneRatio: 40/661, adjust p-value: 4.13×10^-7). With respect to CC, DEGs were mainly concentrated in MHC class II protein complex (GeneRatio: 10/684, adjust p-value: 4.01×10^-11), chromosome-centromeric region (GeneRatio: 27/684, adjust p-value: 1.53×10^-9), clathrin-coated endocytic vesicle membrane (GeneRatio: 16/684, adjust p-value: 2.81×10^-9), MHC protein complex (GeneRatio: 10/684, adjust p-value: 5.23×10^-9) and condensed chromosome, centromeric region (GeneRatio: 22/684, adjust p-value: 5.99×10^-9). DEGs were mainly clustered in MHC class II protein complex binding (GeneRatio: 10/664, adjust p-value: 1.72×10^-8), MHC protein complex binding (GeneRatio: 11/664, adjust p-value: 3.34×10^-8), hormone activity (GeneRatio: 18/664, adjust p-value: 3.95×10^-7), neuropeptide hormone activity (GeneRatio: 9/664, adjust p-value: 7.23×10^-7) and MHC class II receptor activity (GeneRatio: 5/664, adjust p-value: 1.32×10^-5) in MF category. These results suggested that a large number of MYCN related DEGs were enriched into immune-related pathways, and these DEGs might play important roles in the microenvironment of MYCN positive NB.

FIGURE 3
www.frontiersin.org

Figure 3 GO and GSEA enrichment analyses of DEGs (A) Bubble diagram of GO enrichment analysis of DEGs in MYCN positive NB. (B) GSEA analysis of representative immune-related pathways for all genes in MYCN positive NB. (C) GSEA analysis of representative immune-related pathways for all genes in MYCN negative NB.

GSEA

To explore the potential immune regulatory mechanism in NB, we used the annotated gene set entitled “immunesigdb” from MsigDB as a reference for GSEA analysis in this study. According to the standard of adjust p-value<0.05, A total of 439 immune-related pathways were enriched in MYCN positive NB (Figure 3B). The top five enrichment score pathways among MYCN positive NB were GSE15750_DAY6_VS_DAY10_TRAF6KO_EFF_CD8_TCELL_UP (enrichment score: 0.7155, NSE: 0.3534, adjust p-value: 3.12×10^-9), GSE15750_DAY6_VS_DAY10_EFF_CD8_TCELL_UP (enrichment score: 0.7082, NSE: 2.9948, adjust p-value: 3.12×10^-9), GSE36476_CTRL_VS_TSST_ACT_72H_MEMORY_CD4_TCELL_YOUNG_DN (enrichment score: 0.6930, NSE: 2.9574, adjust p-value: 3.12×10^-9) and GSE39556_CD8A_DC_VS_NK_CELL_MOUSE_3H_POST_POLYIC_INJ_UP (enrichment score: 0.6681, NSE: 2.8580, adjust p-value: 3.12×10^-9), GSE18893_TCONV_VS_TREG_24H_TNF_STIM_UP (enrichment score: 0.6648, NSE: 2.8257, adjust p-value: 3.12×10^-9). 1904 immune-related pathways were enriched in MYCN negative NB (Figure 3C). The top five pathways enriched among MYCN negative NB were GSE7218_UNSTIM_VS_ANTIGEN_STIM_THROUGH_IGG_BCELL_DN (enrichment score: -0.7259, NSE: -2.9380, adjust p-value: 3.12×10^-9), GSE7218_IGM_VS_IGG_SIGNAL_THGOUGH_ANTIGEN_BCELL_DN (enrichment score: -0.6682, NSE: -2.7121, adjust p-value: 3.12×10^-9), GSE10325_LUPUS_CD4_TCELL_VS_LUPUS_BCELL_UP (enrichment score: -0.6428, NSE: -2.6871, adjust p-value: 3.12×10^-9), GSE7509_UNSTIM_VS_FCGRIIB_STIM_DC_DN (enrichment score: -0.6405, NSE: -2.6303, adjust p-value: 3.12×10^-9) and GSE7509_UNSTIM_VS_IFNA_STIM_IMMATURE_DC_DN (enrichment score: -0.6363, NSE: -2.6296, adjust p-value: 3.12×10^-9). The above results fully indicated that there were significant differences in immune-related pathways enriched in MYCN positive and MYCN negative NB, especially in the expression of CD4 T, CD8T, B cell and other immune cells. MYCN and MYCN related genes might play an important role in the immune regulation of NB.

Identification of clinically significant module and hub genes by WGCNA

The transcriptome expression data of 185 MYCN positive and 951 MYCN negative NB tissues were preprocessed, and duplicated genes and missing values were removed to obtain a combined matrix containing 19806 genes. To ensure the accuracy of the results, we performed cluster analysis on the samples after removing the outlier samples. The sample dendrogram and trait heatmap are shown in Figure 4A. We selected the minimum threshold 8 that could make the curve tend to be smooth to make the co-representation network to be closed to a scale-free network (Figures 4B, C). Nineteen modules were initially obtained and similar modules were further combined into 15 modules as shown in the hierarchical clustering tree (Figures 5A, B). We found that turquoise module contained the largest number of genes among all modules (Figure 5B), with a total of 1842 genes included. The genetic importance of each module was listed in Figure 5C, from which we could see that turquoise module had the greatest effect on the gene regulatory network in NB (Figure 5C). In order to further understand the correlation between each module and MYCN, we analyzed the clinical information of training group and calculated the p-value of correlation between each module and MYCN expression level (Figure 5D), so as to screen out the module most related to MYCN. We could find that except the black module, salmon module and cyan module, all other modules were correlated with the expression level of MYCN. Among them, magenta module, lightcyan module, grey module and turquoise module were positively correlated with the amplification of MYCN, while other modules are negatively correlated with the expression of MYCN. Finally, turquoise module with the strongest correlation with MYCN (correlation index: 0.69, p-value: 5×10^-163) was obtained as the object of our further study. According to the filtering conditions of gene significance (GS)>0.5 and module membership (MM)>0.8, a total of 27 core genes were screened in turquoise module (Figure 5E). We crossed 166 DEGs (|log2fc| > 2 and p-value < 0.05) with 27 core genes of turquoise module and obtained four hub genes: ZNF695 (GS: 0.6229 MM: 0.8559), CHEK1 (GS: 0.5655 MM: 0.8746), C15ORF42 (GS: 0.5479 MM: 0.9063), EXO1 (GS: 0.5051 MM: 0.8569) (Figure 5F).These results suggested that there were multiple gene clusters in the gene network of NB. Due to the differential expression of MYCN, these gene clusters played different roles in NB. Exploring hub genes in turquoise module would help us understand the important roles of MYCN and its related genes in the gene regulatory network of NB.

FIGURE 4
www.frontiersin.org

Figure 4 Establishment of scale-free network (A) Clustering dendrogram of 1136 samples. (B) Scale-free exponential analysis and average connectivity analysis of various soft-threshold powers (β). (C) Check the scale-free topology, the distribution approximately follows an approximate straight line, which is called approximate scale-free topology when β=8.

FIGURE 5
www.frontiersin.org

Figure 5 Identification of clinically significant module and hub genes by WGCNA (A) As shown in the hierarchical clustering tree, 19 modules were screened out and similar modules were further combined into 15 modules based on WGCNA (Red line: MEDissThres = 0.25). (B) Merged dynamic tree of genes with dissimilarity based on topological overlap and assigning module colors. (C) Diagram of correlation of NB and module’s color. The turquoise module was most closely related to MYCN. (D) Heatmap of the correlation between MYCN expression and module eigengenes. The turquoise module had the lowest p value and was selected as our next research object. (E) Scatterplot analysis of turquoise modules. core genes were filtered in the upper right region of GS > 0.5 and MM > 0.8. (F) Displays the Venn plot of 27 core genes in turquoise modules and 166 DEGs (|log2fc| > 2 and p-value < 0.05).

MYCN positive NB model constructed by Lasso regression and Roc analysis

We further constructed Lasso regression model for the four genes to obtain the characteristic genes of MYCN. Three characteristic genes (ZNF695, CHEK1, C15ORF42) highly associated with MYCN positive NB were obtained by cross-validation (alpha=1) (Figures 6A, B). Roc regression analysis was further performed to evaluate the association between the three genes and MYCN positive NB. In the training group, the AUC values of ZNF695, CHEK1 and C15ORF42 reached 0.956, 0.926 and 0.930 respectively (Figure 6C). In order to verify the results of the training group, we further calculated the AUC values of ZNF695 (0.954), CHEK1 (0.922) and C15ORF42 (0.921) in the validation group (Figure 6D). These results indicated that ZNF695, CHEK1 and C15ORF42 are closely related to MYCN positive NB, which is worthy of further study.

FIGURE 6
www.frontiersin.org

Figure 6 MYCN positive NB model constructed by Lasso and Roc regression analysis (A) Ten-fold cross-validation of tuning parameter selection in Lasso models. (B) Lasso coefficient curves of four MYCN closely related DEGs. ZNF695, CHEK1 and C15ORF42 were eventually screened out. (C) The Roc curves of the three hub genes in 1136 cases. (D) The Roc curves of the three hub genes in the validation set of E-MTAB-8248.

Validation of correlation between ZNF695, CHEK1, C15ORF42 and MYCN positive NB

According to the boxplot, we found that ZNF695, CHEK1 and C15ORF42 were significantly overexpressed in MYCN positive NB in the training group (Figure 7A), and the same results were obtained in the verification group (Figure 7B). We also confirmed differential expression of ZNF695, CHEK1, and C15ORF42 genes in MYCN positive NB cell lines SK-N-BE (2), Kelly, SKNDZ and MYCN negative cell line SH-SY5Y except for the expression of C15ORF42 in Kelly cell line by qRT-PCR (Figure 7C). Survival analysis showed that ZNF695, CHEK1 and C15ORF42 had significant effects on event-free survival rate (GSE62564 cohort; Figure 8A) and overall survival rate (E-MTAB-8248 cohort; Figure 8B) in NB. These results indicated that the three hub genes screened in this study had highly consistent expression trends in NB tissues and cell lines, and their expression levels were closely related to the prognosis of children with NB. This provided bases for us to carry out researches on the mechanism of these three genes in NB.

FIGURE 7
www.frontiersin.org

Figure 7 Validation of correlation between ZNF695, CHEK1, C15ORF42 and MYCN positive NB (A) The comparison of the expression levels of the hub genes between the MYCN positive and negative group with a total of 1136 NB cases. (B) The comparison of the expression levels of the hub genes between the MYCN positive and negative group in validation set. (C) The comparison of the expression levels of the hub genes between the MYCN positive (SKNDZ, Kelly, SK-N-BE(2)) and negative (SH-SY5Y) cell lines. (**p<0.01, ***p<0.001, ****p<0.0001).

FIGURE 8
www.frontiersin.org

Figure 8 Survival analysis of three identified hub genes (A) Event-free survival analysis of three hub genes in GSE62564 cohort. (B) Overall survival analysis of three hub genes in E-MTAB-8248 cohort.

Correlation between immune cell status and MYCN-related DEGs in NB

According to the microarray data, we further scored the immune cell infiltration status of 1136 NB cases in the experimental group, and a total of 28 immune-related cells were included in this study. As shown in the heatmap, the level of infiltration of immune cells was closely related to the expression level of MYCN (Figure 9A). It was found that the infiltration level of each immune cell showed an obvious opposite trend in MYCN positive and MYCN negative NB. Except that activated CD4 T cells and Type 2 T helper cells were significantly enriched in MYCN positive NB, the infiltration level of other immune cells showed an increasing trend in MYCN negative NB. In order to more intuitively reflect the distribution level of each immune cell in MYCN positive and MYCN negative NB, we drew the corresponding violin diagram, from which we could see that the distribution of 28 immune cells was different in MYCN positive and MYCN negative group (Figure 9B). Except that Activated CD4 T cell and Type2 T helper cells increased, the infiltration levels of the other 26 cells decreased significantly in MYCN positive NB tissues. This made us realize that immune cells might play an important role in MYCN positive NB, and the infiltration level of immune cells could also reflect the expression level of MYCN. Finally, we analyzed the correlation between the level of immune cell infiltration and three characteristic genes. The results showed that the infiltration levels of Type2 T helper cell and Activated CD4 T cell were significantly positively correlated with the expression levels of the three hub genes, while the infiltration levels of other 26 immune cells, especially Type17 T helper cell, Type1 T helper cell, T follicular helper cell, Regulatory T cell, Natural killer T cell, Natural killer cell, Monocyte, Immature dendritic cell, Effector memory CD8 T cell, CD56dim natural killer cell and CD56bright natural killer cell decreased significantly with the increased expression level of three hub genes (Figure 9C). These results suggested that MYCN and MYCN related genes have consistent effects on the immune cell infiltration in the NB microenvironment, and MYCN may regulate the microenvironment of MYCN positive NB partly through these hub genes.

FIGURE 9
www.frontiersin.org

Figure 9 Correlation between immune-related cells and MYCN-related hub genes in NB (A) Heatmap shows the compositions of infiltrated immune cells between MYCN positive and MYCN negative group with a total of 1136 NB cases. (B) The violin plot shows the abundance of 28 immune-related cells in 1136 cases, comparisons between immune cells in MYCN positive and MYCN negative NB. (C) Heatmap shows the correlations between 28 immune-related cells and three hub genes in 1136 cases.

Discussion

Tumor formation is a very complex process, which involves the abnormal expression of a large number of tumor related genes and the changes of tumor microenvironment (32). It is particularly important to deeply understand the tumor microenvironment and strengthen the research of drugs targeting the tumor microenvironment. With the development of high-throughput sequencing (33) and single cell sequencing technology (34), the understanding level of tumor microenvironment has been further improved. Therefore, based on GEO and Arrayexpress database, this paper screened out the core genes closely related to MYCN positive NB and revealed the relationship between these genes and NB microenvironment.

We screened MYCN related genes and found that even if we raised the screening criteria to |log2fc| > 1.2 and p-value < 0.05, 880 DEGs were screened out. As an important transcription factor, N-myc plays an important role in the tumor biological characteristics of MYCN positive NB. As the results showed, the target genes of N-myc such as Twist1 (35), GLDC (36), TP53 (37) were significantly increased in MYCN positive NB. Therefore, MYCN gene may regulate the expression of key genes in the process of tumor metabolism (14), cycle (38), apoptosis (16), immunity (15) and so on through the highly expressed N-myc protein.

By immune-related pathway enrichment of DEGs, we found that DEGs were significantly enriched in NB microenvironment, especially in T cells’ activation, differentiation, adhesion and formation of MHC proteins. GSEA enrichment analysis also showed that the regulation of CD8+T cells, CD4+T cells, NK cells and Treg cells was closely related to the expression level of MYCN. In the past few years, more attention has been paid to the research on MYC, there have also been studies focusing on the effect of age of onset on immune cell infiltration in NB (39). However, with a deeper understanding of MYCN, we have found that MYCN plays an important role in many fields of tumor microenvironment. For example, in the early years, we found that N-myc suppressed the expression of MHC I through enhancer inactivation (40). Subsequently, it was found that N-myc could mediate tumor immune escape by inhibiting NKT cell enrichment to the site of disease in NB (15). In recent years, it was found that N-myc could regulate the expression of PD-L1, thereby regulating tumor immune surveillance and immune escape (41). Through WGCNA and Lasso regression analysis, we finally screened three hub genes closely related to MYCN: ZNF695, CHEK1 and C15ORF42. Since MYCN is closely related to NB tumor immunity, the roles of these three genes in the regulation of the microenvironment of MYCN positive NB are worthy of further study.

At present, there are few studies on ZNF695, let alone related reports on ZNF695 in NB. According to the projections provided by Alliance of Genome Resources, ZNF695 may enable DNA-binding transcription factor activity, RNA polymerase II-specific and RNA polymerase II cis-regulatory region sequence-specific DNA binding activity, and is predicted to be involved in regulation of transcription (https://www.alliancegenome.org/gene/HGNC:30954#summary). Ke ZB et al. found that ZNF695 is an important marker of radiotherapy resistance in prostate cancer (42), and is closely related to the proliferation, invasion and migration of prostate cancer cells, which may be an important factor affecting the microenvironment of prostate cancer (42, 43).Takahashi T et al. found that methylation of ZNF695 is an important factor affecting the chemotherapeutic resistance of esophageal squamous cell carcinoma (44). Some studies have also found that ZNF695 is closely related to the prognosis of breast cancer (45) and ovarian cancer (46). These results fully indicated that ZNF695 is closely related to the prognosis of many tumors, and may have an important role in the microenvironment of malignant tumors.

CHK1 (CHEK1), a member of the CHEK family, is a serine/threonine-specific protein kinase that mediates cell cycle arrest in response to DNA damage (47). Previously, CHK1 was considered to be a tumor suppressor gene because of its role in regulating DNA damage (48). However, current studies have found that this gene is overexpressed in some solid tumors and affects the prognosis of these patients (49). Therefore, it is of great significance to study the regulatory mechanism of CHK1 in solid tumors. Previous studies have found that N-myc in MYCN positive tumor cells abnormally regulates G1/S checkpoint, inhibits G1 arrest after DNA damage, promotes tumor cells to enter the S phase with DNA damage repairing and experiencing chronic replication stress (50, 51). They are more dependent on ATR/CHK1 pathway to solve DNA damage during replication (51). CHK1 inhibitors can well inhibit the protective effect of ATR/CHK1 pathway on MYCN positive tumor cells (52). A number of studies have shown that CHK1 inhibitors or combined with Wee1 (53), R9-capep (54) and other drugs have significant killing effects on MYCN positive NB cell lines, and CHK1 inhibitors can induces regression of preclinical models of human NB (55, 56). However, there are currently no reports on the correlation between CHK1 and the microenvironment of MYCN positive NB.

C15ORF42 (TICRR), a human homolog of yeast sld3 protein (57), is an important DNA replication factor regulated by cyclin dependent kinases and DNA damage checkpoints. Yu Q et al. suggested that C15ORF42 participates in tumorigenesis by accelerating DNA replication, and a higher level of C15ORF42 is associated with poor overall survival and disease free survival in multiple cancer types (58). For example, the increased expression of C15ORF42 in papillary renal cell carcinoma may participate in tumorigenesis by regulating the cell cycle, and have an impact on prognosis (59). Wang Y et al. found that C15ORF42 is closely related to the prognosis of pediatric brain tumors, and may affect the prognosis of children by changing the tumor microenvironment (60). At present, there are no studies on the effects of C15ORF42 in the microenvironment of children with NB, and no researches on the correlation between MYCN and C15ORF42. We found that there are specific binding sites of N-myc and N-myc was highly enriched in the transcription initiation region of CHK1 and C15ORF42 through UCSC database (http://genome.ucsc.edu/) and chip-seq data from Zeid R’s research (61)(Figures 10A, B), which suggests that N-myc may be an important transcription factor of CHK1 and C15ORF42.

FIGURE 10
www.frontiersin.org

Figure 10 Prediction of the correlation between CHEK1, C15ORF42 and MYCN (A) The UCSC database predicts that N-myc is highly enriched near the CHEK1 transcription start site. (B) The UCSC database predicts that N-myc is highly enriched near the C15ORF42 transcription start site.

In conclusion, the precise mechanisms of ZNF695, CHEK1 and C15ORF42 in MYCN positive NB are still unclear. This study successfully predicts that these three genes, as key genes in the regulation of tumor biological characteristics in MYCN positive NB, are closely related to tumor microenvironment. The limitation of this study is that all the above results are based on microarray data and bioinformatics technology, and have not been fully verified by experiments. In the next step, we will further explore the mechanism of the correlation between these three hub genes and the microenvironment of MYCN positive NB.

Conclusion

The treatment of MYCN positive high-risk NB is still a difficulty in clinical work, and currently there is no targeted drug directly targeting N-myc. In this study, three hub genes closely related to MYCN were screened out by WGCNA and Lasso regression analysis through GEO and ArrayExpress databases. It is predicted that these three genes may be closely related to the microenvironment of MYCN positive NB. Subsequent mechanism studies may bring new ideas for the clinical diagnosis and treatment of MYCN positive NB and bring hope to children with NB.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author contributions

BJ and YF, conceptualization, methodology, writing-reviewing, and editing. JC and MS, investigation, data curation, and writing-original draft preparation. JC and CC, visualization, validation, supervision, and software. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by the General Project of Nanjing Medical University (NMUB20210060), National Natural Science Foundation of China (81903383), Natural Science Foundation of Jiangsu Province (BK20211009), Scientific Research Projects of Jiangsu Health Commission (ZDB2020018), China Postdoctoral Science Foundation funded project (2021M701764), Special Fund for Health Science and Technology Development in Nanjing (JQX19008), Nanjing Medical Science and Technology Development Project (YKK21149), Young Talent Support Project of Children’s Hospital of Nanjing Medical University (TJGC2020016, TJGC2020007, TJGC2020014).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.1016683/full#supplementary-material

References

1. Ward E, DeSantis C, Robbins A, Kohler B, Jemal A. Childhood and adolescent cancer statistics, 2014. CA Cancer J Clin (2014) 64(2):83–103. doi: 10.3322/caac.21219

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Strother DR, London WB, Schmidt ML, Brodeur GM, Shimada H, Thorner P, et al. Outcome after surgery alone or with restricted use of chemotherapy for patients with low-risk neuroblastoma: results of children's oncology group study P9641. J Clin Oncology: Off J Am Soc Clin Oncol (2012) 30(15):1842–8. doi: 10.1200/JCO.2011.37.9990

CrossRef Full Text | Google Scholar

3. Twist CJ, Schmidt ML, Naranjo A, London WB, Tenney SC, Marachelian A, et al. Maintaining outstanding outcomes using response- and biology-based therapy for intermediate-risk neuroblastoma: A report from the children's oncology group study ANBL0531. J Clin Oncology: Off J Am Soc Clin Oncol (2019) 37(34):3243–55. doi: 10.1200/JCO.19.00919

CrossRef Full Text | Google Scholar

4. Baker DL, Schmidt ML, Cohn SL, Maris JM, London WB, Buxton A, et al. Outcome after reduced chemotherapy for intermediate-risk neuroblastoma. New Engl J Med (2010) 363(14):1313–23. doi: 10.1056/NEJMoa1001527

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Matthay KK, Villablanca JG, Seeger RC, Stram DO, Harris RE, Ramsay NK, et al. Treatment of high-risk neuroblastoma with intensive chemotherapy, radiotherapy, autologous bone marrow transplantation, and 13-cis-retinoic acid. Children's Cancer Group New Engl J Med (1999) 341(16):1165–73. doi: 10.1056/NEJM199910143411601

CrossRef Full Text | Google Scholar

6. Holmes K, Pötschger U, Pearson ADJ, Sarnacki S, Cecchetto G, Gomez-Chacon J, et al. Influence of surgical excision on the survival of patients with stage 4 high-risk neuroblastoma: A report from the HR-NBL1/SIOPEN study. J Clin Oncology: Off J Am Soc Clin Oncol (2020) 38(25):2902–15. doi: 10.1200/JCO.19.03117

CrossRef Full Text | Google Scholar

7. Pearson ADJ, Pinkerton CR, Lewis IJ, Imeson J, Ellershaw C, Machin D. High-dose rapid and standard induction chemotherapy for patients aged over 1 year with stage 4 neuroblastoma: a randomised trial. Lancet Oncol (2008) 9(3):247–56. doi: 10.1016/S1470-2045(08)70069-X

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Yu AL, Gilman AL, Ozkaynak MF, London WB, Kreissman SG, Chen HX, et al. Anti-GD2 antibody with GM-CSF, interleukin-2, and isotretinoin for neuroblastoma. New Engl J Med (2010) 363(14):1324–34. doi: 10.1056/NEJMoa0911123

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Furman WL, McCarville B, Shulkin BL, Davidoff A, Krasin M, Hsu C-W, et al. Improved outcome in children with newly diagnosed high-risk neuroblastoma treated with chemoimmunotherapy: Updated results of a phase II study using hu14. 18K322A. J Clin Oncology: Off J Am Soc Clin Oncol (2022) 40(4):335–44. doi: 10.1200/JCO.21.01375

CrossRef Full Text | Google Scholar

10. Abbasi J. Mixed findings in pediatric neuroblastoma CAR-T therapy trial. JAMA (2021) 325(2):121. doi: 10.1001/jama.2020.26233

CrossRef Full Text | Google Scholar

11. Park JA, Cheung N-KV. Targets and antibody formats for immunotherapy of neuroblastoma. J Clin Oncology: Off J Am Soc Clin Oncol (2020) 38(16):1836–48. doi: 10.1200/JCO.19.01410

CrossRef Full Text | Google Scholar

12. Irwin MS, Naranjo A, Zhang FF, Cohn SL, London WB, Gastier-Foster JM, et al. Revised neuroblastoma risk classification system: A report from the children's oncology group. J Clin Oncology: Off J Am Soc Clin Oncol (2021) 39(29):3229–41. doi: 10.1200/JCO.21.00278

CrossRef Full Text | Google Scholar

13. Baluapuri A, Wolf E, Eilers M. Target gene-independent functions of MYC oncoproteins. Nat Rev Mol Cell Biol (2020) 21(5):255–67. doi: 10.1038/s41580-020-0215-2

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Guo Y-F, Duan J-J, Wang J, Li L, Wang D, Liu X-Z, et al. Inhibition of the ALDH18A1-MYCN positive feedback loop attenuates MYCN-amplified neuroblastoma growth Sci Trans Med (2020) 12(531):eaax8694. doi: 10.1126/scitranslmed.aax8694

CrossRef Full Text | Google Scholar

15. Song L, Ara T, Wu H-W, Woo C-W, Reynolds CP, Seeger RC, et al. Oncogene MYCN regulates localization of NKT cells to the site of disease in neuroblastoma. J Clin Invest (2007) 117(9):2702–12. doi: 10.1172/JCI30751

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Ham J, Costa C, Sano R, Lochmann TL, Sennott EM, Patel NU, et al. Exploitation of the apoptosis-primed state of MYCN-amplified neuroblastoma to develop a potent and specific targeted therapy combination. Cancer Cell (2016) 29(2):159–72. doi: 10.1016/j.ccell.2016.01.002

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Alborzinia H, Flórez AF, Kreth S, Brückner LM, Yildiz U, Gartlgruber M, et al. MYCN mediates cysteine addiction and sensitizes neuroblastoma to ferroptosis. Nat Cancer (2022) 3(4):471–85. doi: 10.1038/s43018-022-00355-4

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Xu Y, Poggio M, Jin HY, Shi Z, Forester CM, Wang Y, et al. Translation control of the immune checkpoint in cancer and its therapeutic targeting. Nat Med (2019) 25(2):301–11. doi: 10.1038/s41591-018-0321-2

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Maeda T, Hiraki M, Jin C, Rajabi H, Tagde A, Alam M, et al. MUC1-c induces PD-L1 and immune evasion in triple-negative breast cancer. Cancer Res (2018) 78(1):205–15. doi: 10.1158/0008-5472.CAN-17-1636

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Casey SC, Tong L, Li Y, Do R, Walz S, Fitzgerald KN, et al. MYC regulates the antitumor immune response through CD47 and PD-L1. Sci (New York NY) (2016) 352(6282):227–31. doi: 10.1126/science.aac9935

CrossRef Full Text | Google Scholar

21. Casey SC, Li Y, Fan AC, Felsher DW. Oncogene withdrawal engages the immune system to induce sustained cancer regression. J Immunother Cancer (2014) 2:24. doi: 10.1186/2051-1426-2-24

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Dhanasekaran R, Baylot V, Kim M, Kuruvilla S, Bellovin DI, Adeniji N, et al. MYC and Twist1 cooperate to drive metastasis by eliciting crosstalk between cancer and innate immunity. Elife (2020) 9:e50731. doi: 10.7554/eLife.50731

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Collins PL, Cella M, Porter SI, Li S, Gurewitz GL, Hong HS, et al. Gene regulatory programs conferring phenotypic identities to human NK cells. Cell (2019) 176(1-2):348-360.e12. doi: 10.1016/j.cell.2018.11.045

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Sodir NM, Kortlever RM, Barthet VJA, Campos T, Pellegrinet L, Kupczak S, et al. MYC instructs and maintains pancreatic adenocarcinoma phenotype. Cancer Discov (2020) 10(4):588–607. doi: 10.1158/2159-8290.CD-19-0435

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Xu Y, Shi J, Yamamoto N, Moss JA, Vogt PK, Janda KD. A credit-card library approach for disrupting protein-protein interactions. Bioorg Med Chem (2006) 14(8):2660–73. doi: 10.1016/j.bmc.2005.11.052

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Panda D, Saha P, Das T, Dash J. Target guided synthesis using DNA nano-templates for selectively assembling a G-quadruplex binding c-MYC inhibitor. Nat Commun (2017) 8:16103. doi: 10.1038/ncomms16103

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Lu X, Vogt PK, Boger DL, Lunec J. Disruption of the MYC transcriptional function by a small-molecule antagonist of MYC/MAX dimerization. Oncol Rep (2008) 19(3):825–30. doi: 10.3892/or.19.3.825

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Yin X, Giap C, Lazo JS, Prochownik EV. Low molecular weight inhibitors of myc-max interaction and function. Oncogene (2003) 22(40):6151–9. doi: 10.1038/sj.onc.1206641

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Huang M, Weiss WA. Neuroblastoma and MYCN. Cold Spring Harb Perspect Med (2013) 3(10):a014415. doi: 10.1101/cshperspect.a014415

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Langfelder P, Horvath S. WGCNA: an r package for weighted correlation network analysis. BMC Bioinf (2008) 9:559. doi: 10.1186/1471-2105-9-559

CrossRef Full Text | Google Scholar

31. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-delta delta C(T)) method. Methods (2001) 25(4):402–8. doi: 10.1006/meth.2001.1262

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Pitt JM, Marabelle A, Eggermont A, Soria JC, Kroemer G, Zitvogel L. Targeting the tumor microenvironment: removing obstruction to anticancer immune responses and immunotherapy. Ann Oncol (2016) 27(8):1482–92. doi: 10.1093/annonc/mdw168

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Rizvi H, Sanchez-Vega F, La K, Chatila W, Jonsson P, Halpenny D, et al. Molecular determinants of response to anti-programmed cell death (PD)-1 and anti-programmed death-ligand 1 (PD-L1) blockade in patients with non-Small-Cell lung cancer profiled with targeted next-generation sequencing. J Clin Oncol Off J Am Soc Clin Oncol (2018) 36(7):633–41. doi: 10.1200/JCO.2017.75.3384

CrossRef Full Text | Google Scholar

34. Finotello F, Rieder D, Hackl H, Trajanoski Z. Next-generation computational tools for interrogating cancer immunity. Nat Rev Genet (2019) 20(12):724–46. doi: 10.1038/s41576-019-0166-7

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Selmi A, de Saint-Jean M, Jallas A-C, Garin E, Hogarty MD, Bénard J, et al. TWIST1 is a direct transcriptional target of MYCN and MYC in neuroblastoma. Cancer Lett (2015) 357(1):412–8. doi: 10.1016/j.canlet.2014.11.056

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Alptekin A, Ye B, Yu Y, Poole CJ, van Riggelen J, Zha Y, et al. Glycine decarboxylase is a transcriptional target of MYCN required for neuroblastoma cell proliferation and tumorigenicity. Oncogene (2019) 38(50):7504–20. doi: 10.1038/s41388-019-0967-3

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Chen L, Iraci N, Gherardi S, Gamble LD, Wood KM, Perini G, et al. p53 is a direct transcriptional target of MYCN in neuroblastoma. Cancer Res (2010) 70(4):1377–88. doi: 10.1158/0008-5472.CAN-09-2598

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Hogarty MD. The requirement for evasion of programmed cell death in neuroblastomas with MYCN amplification. Cancer Lett (2003) 197(1-2):173–9. doi: 10.1016/S0304-3835(03)00103-4

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Asgharzadeh S, Salo JA, Ji L, Oberthuer A, Fischer M, Berthold F, et al. Clinical significance of tumor-associated inflammatory cells in metastatic neuroblastoma. J Clin Oncology: Off J Am Soc Clin Oncol (2012) 30(28):3525–32. doi: 10.1200/JCO.2011.40.9169

CrossRef Full Text | Google Scholar

40. Lenardo M, Rustgi AK, Schievella AR, Bernards R. Suppression of MHC class I gene expression by n-myc through enhancer inactivation. EMBO J (1989) 8(11):3351–5. doi: 10.1002/j.1460-2075.1989.tb08497.x

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Melaiu O, Mina M, Chierici M, Boldrini R, Jurman G, Romania P, et al. PD-L1 is a therapeutic target of the bromodomain inhibitor JQ1 and, combined with HLA class I, a promising prognostic biomarker in neuroblastoma. Clin Cancer Research: an Off J Am Assoc For Cancer Res (2017) 23(15):4462–72. doi: 10.1158/1078-0432.CCR-16-2601

CrossRef Full Text | Google Scholar

42. Ke Z-B, You Q, Chen J-Y, Sun J-B, Xue Y-T, Zhuang R-B, et al. A radiation resistance related index for biochemical recurrence and tumor immune environment in prostate cancer patients. Comput Biol Med (2022) 146:105711. doi: 10.1016/j.compbiomed.2022.105711

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Zhang L, Li Y, Wang X, Ping Y, Wang D, Cao Y, et al. Five-gene signature associating with Gleason score serve as novel biomarkers for identifying early recurring events and contributing to early diagnosis for prostate adenocarcinoma. J Cancer (2021) 12(12):3626–47. doi: 10.7150/jca.52170

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Takahashi T, Yamahsita S, Matsuda Y, Kishino T, Nakajima T, Kushima R, et al. ZNF695 methylation predicts a response of esophageal squamous cell carcinoma to definitive chemoradiotherapy. J Cancer Res Clin Oncol (2015) 141(3):453–63. doi: 10.1007/s00432-014-1841-x

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Shinden Y, Hirashima T, Nohata N, Toda H, Okada R, Asai S, et al. Molecular pathogenesis of breast cancer: impact of miR-99a-5p and miR-99a-3p regulation on oncogenic genes. J Hum Genet (2021) 66(5):519–34. doi: 10.1038/s10038-020-00865-y

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Juárez-Méndez S, Zentella-Dehesa A, Villegas-Ruíz V, Pérez-González OA, Salcedo M, López-Romero R, et al. Splice variants of zinc finger protein 695 mRNA associated to ovarian cancer. J Ovarian Res (2013) 6(1):61. doi: 10.1186/1757-2215-6-61

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Lunardi A, Varmeh S, Chen M, Taulli R, Guarnerio J, Ala U, et al. Suppression of CHK1 by ETS family members promotes DNA damage response bypass and tumorigenesis. Cancer Discov (2015) 5(5):550–63. doi: 10.1158/2159-8290.CD-13-1050

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Yarden RI, Pardo-Reoyo S, Sgagias M, Cowan KH, Brody LC. BRCA1 regulates the G2/M checkpoint by activating CHK1 kinase upon DNA damage. Nat Genet (2002) 30(3):285–9. doi: 10.1038/ng837

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Fadaka AO, Bakare OO, Sibuyi NRS, Klein A. Gene expression alterations and molecular analysis of CHEK1 in solid tumors. Cancers (Basel) (2020) 12(3):662. doi: 10.3390/cancers12030662

CrossRef Full Text | Google Scholar

50. Paffhausen T, Schwab M, Westermann F. Targeted MYCN expression affects cytotoxic potential of chemotherapeutic drugs in neuroblastoma cells. Cancer Lett (2007) 250(1):17–24. doi: 10.1016/j.canlet.2006.09.010

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Southgate HED, Chen L, Tweddle DA, Curtin NJ. ATR inhibition potentiates PARP inhibitor cytotoxicity in high risk neuroblastoma cell lines by multiple mechanisms. Cancers (Basel) (2020) 12(5):1095. doi: 10.3390/cancers12051095

CrossRef Full Text | Google Scholar

52. Di Giulio S, Colicchia V, Pastorino F, Pedretti F, Fabretti F, Nicolis di Robilant V, et al. A combination of PARP and CHK1 inhibitors efficiently antagonizes MYCN-driven tumors. Oncogene (2021) 40(43):6143–52. doi: 10.1038/s41388-021-02003-0

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Russell MR, Levin K, Rader J, Belcastro L, Li Y, Martinez D, et al. Combination therapy targeting the CHK1 and Wee1 kinases shows therapeutic efficacy in neuroblastoma. Cancer Res (2013) 73(2):776–84. doi: 10.1158/0008-5472.CAN-12-2669

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Gu L, Chu P, Lingeman R, McDaniel H, Kechichian S, Hickey RJ, et al. The mechanism by which MYCN amplification confers an enhanced sensitivity to a PCNA-derived cell permeable peptide in neuroblastoma cells. EBioMedicine (2015) 2(12):1923–31. doi: 10.1016/j.ebiom.2015.11.016

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Walton MI, Eve PD, Hayes A, Valenti MR, De Haven Brandon AK, Box G, et al. CCT244747 is a novel potent and selective CHK1 inhibitor with oral efficacy alone and in combination with genotoxic anticancer drugs. Clin Cancer Research: an Off J Am Assoc For Cancer Res (2012) 18(20):5650–61. doi: 10.1158/1078-0432.CCR-12-1322

CrossRef Full Text | Google Scholar

56. Lowery CD, VanWye AB, Dowless M, Blosser W, Falcon BL, Stewart J, et al. The checkpoint kinase 1 inhibitor prexasertib induces regression of preclinical models of human neuroblastoma. Clin Cancer Research: An Off J Am Assoc For Cancer Res (2017) 23(15):4354–63. doi: 10.1158/1078-0432.CCR-16-2876

CrossRef Full Text | Google Scholar

57. Sanchez-Pulido L, Diffley JFX, Ponting CP. Homology explains the functional similarities of Treslin/Ticrr and Sld3. Curr Biol (2010) 20(12):R509–R10. doi: 10.1016/j.cub.2010.05.021

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Yu Q, Pu S-Y, Wu H, Chen X-Q, Jiang J-J, Gu K-S, et al. TICRR contributes to tumorigenesis through accelerating DNA replication in cancers. Front Oncol (2019) 9:516. doi: 10.3389/fonc.2019.00516

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Xia S, Lin Y, Lin J, Li X, Tan X, Huang Z. Increased expression of predicts poor clinical outcomes: A potential therapeutic target for papillary renal cell carcinoma. Front Genet (2020) 11:605378. doi: 10.3389/fgene.2020.605378

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Wang Y, Zhou C, Luo H, Cao J, Ma C, Cheng L, et al. Prognostic implications of immune-related eight-gene signature in pediatric brain tumors. Braz J Med Biol Res (2021) 54(7):e10612. doi: 10.1590/1414-431x2020e10612

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Zeid R, Lawlor MA, Poon E, Reyes JM, Fulciniti M, Lopez MA, et al. Enhancer invasion shapes MYCN-dependent transcriptional amplification in neuroblastoma. Nat Genet (2018) 50(4):515–23. doi: 10.1038/s41588-018-0044-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: MYCN, neuroblastoma, WGCNA, LASSO, tumor microenvironment, immune cell infiltration

Citation: Chen J, Sun M, Chen C, Jiang B and Fang Y (2022) Identification of hub genes and their correlation with infiltration of immune cells in MYCN positive neuroblastoma based on WGCNA and LASSO algorithm. Front. Immunol. 13:1016683. doi: 10.3389/fimmu.2022.1016683

Received: 11 August 2022; Accepted: 28 September 2022;
Published: 12 October 2022.

Edited by:

Hui Zhao, The Chinese University of Hong Kong, China

Reviewed by:

Nevim Aygun, Ege University, Turkey
Simon Keane, University of Gothenburg, Sweden

Copyright © 2022 Chen, Sun, Chen, Jiang and Fang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Bin Jiang, jiangbin363@163.com; Yongjun Fang, fyj322@189.cn

These authors have contributed equally to this work and share first authorship

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.