Abstract Background Basal cell carcinoma (BCC) is a frequent malignant tumor of skin cancers with high morbidity. The objective of this study was to identify critical genes and pathways related to the carcinogenesis of BCC and gain more insights into the underlying molecular mechanisms of BCC. Materials and methods The gene expression profiles of [43]GSE7553 and [44]GSE103439 were downloaded from the Gene Expression Omnibus database with 19 tumors and 6 normal skin tissues. Differentially expressed genes (DEGs) were screened between BCC samples and normal tissues, followed by gene ontology and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis. Subsequently, protein–protein interaction (PPI) network was constructed for these DEGs, and module analysis was performed. Results A total of 313 DEGs were obtained. Among them, 222 genes were upregulated and 91 genes were downregulated. Enrichment analysis indicated that the upregulated genes were significantly enriched in cell cycle and mitosis, while the downregulated genes were mainly associated with unsaturated fatty acid metabolic process and cell differentiation. In addition, TOP2A, CDK1, and CCNB1 were identified as the top three hub genes ranked by degrees in the PPI network. Meanwhile, three subnetworks were derived, which indicated that these DEGs were significantly enriched in pathways, including “cell cycle”, “extracellular matrix–receptor interaction”, “basal cell carcinoma”, and “hedgehog signaling pathway”. Conclusions The novel critical DEGs and pathways identified in this study may serve pivotal roles in the carcinogenesis of BCC and indicate more molecular targets for the treatment of BCC. Keywords: basal cell carcinoma, differentially expressed genes, enrichment analysis, bioinformatics analysis Introduction Cutaneous basal cell carcinoma (BCC) is recognized as a common subtype of nonmelanoma skin malignancies with high morbidity, which accounts for ~80% of newly diagnosed nonmelanoma skin carcinomas.[45]^1 In the last decade, there has been a substantial increase in the incidence of BCC.[46]^2 Due to the characteristics of slow-growing and locally aggressive, metastasis rarely occurred in patients with BCC, which resulted in a relatively good prognosis. As we all know, long-term exposure to sunlight, especially ultraviolet light, is considered as the main risk factor of skin cancers.[47]^3 However, the underlying molecular mechanisms for the development of BCC has not been completely illuminated. Meanwhile, the treatments of BCC are limited and drug resistance is ubiquitous in advanced or metastatic BCC patients. Therefore, an urgent need exists for further exploring the potential mechanisms of BCC and finding more effective molecular targets for the treatment of BCC. To date, several signaling pathways and molecules have been demonstrated to be involved in the tumorigenesis and progression of BCC at the molecular level, such as the hedgehog signaling pathway.[48]^4 Genes included in this pathway, such as the hedgehog receptors patched (PTCH1) or smoothened (SMO), have been extensively studied.[49]^5^,[50]^6 Mutations in these genes may cause constitutive hedgehog pathway activation, which promote the development of BCC. Recently, two new hedgehog pathway inhibitors, Vismodegib and Sonidegib, have been approved by the Food and Drug Administration for the targeted treatment of BCC.[51]^7^,[52]^8 However, the response rate of advanced or metastatic BCC is not promising and the secondary drug resistance may also occur. With the development of high-throughput technology, more and more new potential targets have been uncovered in BCC. In addition to canonical hedgehog pathway components, the transcription factor serum response factor was identified as a noncanonical hedgehog activator by multidimensional genomics analysis, which leads to the amplification of the hedgehog transcription factor glioma-associated oncogene family zinc finger-1 (GLI1).[53]^9 At the DNA level, Bonilla et al performed a genomic analysis of 293 BCC samples and revealed that mutations in other cancer-related genes also drove the initiation of BCC, including MYCN, PTPN14, and LATS1.[54]^10 Thus, much more molecular targets remain to be elucidated. Bioinformatics analysis of gene expression profiles or other high-throughput data are now playing a critical role in investigating the mechanisms of human disease, particularly in tumors. Accordingly, in the present study, we first time integratively reanalyzed the gene expression profiles of 19 BCC and 6 normal tissues deposited in two datasets by differentially expressed genes (DEGs) screening and functional and pathway enrichment analysis. By protein–protein interaction (PPI) network analysis, we identified top three hub genes (TOP2A, CDK1, and CCNB1). Finally, module analysis revealed that several critical pathways were mainly associated with the carcinogenesis of BCC, which might be used as molecular targets for the treatment of BCC. Materials and methods Microarray data Two datasets ([55]GSE7553 and [56]GSE103439) were respectively retrieved from Gene Expression Omnibus database ([57]http://www.ncbi.nlm.nih.gov/geo/), including 19 BCC and 6 normal tissues ([58]Table 1).[59]^11 These gene expression profiles were generated by [60]GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array) containing 54,675 probes. The latest annotation file of [61]GPL570 platform was downloaded from Affymetrix official website ([62]http://www.affymetrix.com/), in which 54,675 probes now mapped to 21,297 genes. Table 1. The basal information of two datasets in this study GEO datasets Platform Number of BCC Number of NS [63]GSE7553 [64]GPL570 15 4 [65]GSE103439 [66]GPL570 4 2 [67]Open in a new tab Abbreviations: BCC, basal cell carcinoma; GEO, Gene Expression Omnibus; NS, normal skin. Data preprocessing and DEGs screening The raw data files (.CEL files) of these 25 samples were processed by the R package “affy”.[68]^12 Background adjustment and normalization were performed using the Robust Multichip Average algorithm. Once multiple probes mapped to the same gene, the average value was finally selected to represent the gene expression value. DEGs were screened between BCC and normal tissues by the “limma” package in R.[69]^13 Then, hierarchical clustering analysis was applied to the DEGs by the “pheatmap” package in R based on the Euclidean distance. The criteria of DEGs was set as |log[2]fold change|>1 and false discovery rate (FDR) <0.05. Functional and pathway enrichment analysis Gene ontology (GO) analysis defines the functions of gene products covering three domains, including biological process, molecular function, and cellular component.[70]^14^,[71]^15 The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database is widely used to map large-scale datasets to pathway maps for higher-order functional information.[72]^16 The Database for Annotation, Visualization and Integrated Discovery (DAVID version 6.8, [73]http://david.abcc.ncifcrf.gov/) consists of an integrated biological knowledgebase and analytic tools, which can systematically extract biological meaning from large gene/protein lists.[74]^17 With the online DAVID tool, we performed functional and pathway enrichment analysis for these DEGs. P-value <0.05 was considered as significant. Construction of PPI network and module analysis Given the large number of DEGs, the “STRINGdb” package in R was used to investigate the potential interactions that existed in these DEGs.[75]^18 Briefly, 313 DEGs were mapped to their corresponding proteins in the Search Tool for the Retrieval of Interacting Genes/Proteins database. Only interactions with a combined score of >0.4 were imported into Cytoscape software to visualize the PPI network.[76]^19 Each node in the network represents one protein, and the degree of each protein was termed as the number of its interactions. Then, the Molecular Complex Detection (MCODE) plug-in was used to analyze the PPI network to identify significant modules.[77]^20 In addition, the functional and pathway enrichment analysis of genes in the subnetworks were performed. P-value <0.05 was set as the threshold. Results Identification of DEGs We screened DEGs in the two datasets ([78]GSE7553 and [79]GSE103439). Compared with normal skin tissues, 1,871 DEGs and 5,357 DEGs were obtained, respectively ([80]Figure 1A). Finally, a total of 313 aberrantly expressed genes (222 upregulated genes and 91 downregulated genes) were identified by integrated analysis ([81]Figure 1B and C). Strikingly, the number of upregulated genes were largely more than down-regulated genes ([82]Table S1). The heatmap of hierarchical clustering analysis showed that these DEGs could clearly distinguish BCC samples from the normal skin samples ([83]Figure 1D and E). Figure 1. [84]Figure 1 [85]Open in a new tab DEGs in the two datasets. Notes: (A) Common DEGs between [86]GSE7553 and [87]GSE103439. (B) Common upregulated DEGs between [88]GSE7553 and [89]GSE103439. (C) Common downregulated DEGs between [90]GSE7553 and [91]GSE103439. (D, E) Hierarchical clustering analysis of the DEGs in [92]GSE7553 and [93]GSE103439, respectively. Red and green indicate higher expression and lower expression, respectively. Abbreviation: DEGs, differentially expressed genes. GO and KEGG pathway enrichment analysis To further investigate the potential functions of these 313 DEGs, GO and KEGG pathways enrichment analysis was performed by the online DAVID tool. The results of GO analysis indicated that upregulated genes enriched in biological process were mainly involved in cell cycle and mitosis, such as the cell division (P=4.39×10^−11) and the mitotic nuclear division (P=5.90×10^−8) ([94]Table 2). Meanwhile, downregulated genes were significantly enriched in unsaturated fatty acid metabolic process (P=2.10×10^−3) and cell differentiation (P=6.76×10^−3) ([95]Table 3). With regard to pathway enrichment analysis, the most significant pathway of upregulated genes was cell cycle (P=4.75×10^−9) containing 13 genes. Interestingly, another five genes (LEF1, PTCH1, GLI2, FZD7, and GLI1) were enriched in the pathway named “basal cell carcinoma” (P=2.60×10^−3) ([96]Table 2), while downregulated genes were most significantly involved in the biosynthesis and metabolism of unsaturated fatty acids (P=5.26×10^−3) ([97]Table 3). Table 2. The top 10 GO terms and KEGG pathways of upregulated genes Term Count P-value Genes GO:0051301:cell division 24 4.39E-11 KIF14, CDK1, KIF11, NEK2, NUF2, KIF18B, NDC80, BIRC5, CDC20, CDC25C, MCM5, CCNE2, CCNB1, SPC25, MAD2L1, CCNB2, HMCN1, SGO2, SPAG5, NCAPG, NCAPG2, ZWINT, CENPW, BUB1B GO:0007067:mitotic nuclear division 17 5.90E-08 CDK1, KIF11, NEK2, KIF15, NUF2, BIRC5, NDC80, CDC20, PBK, CEP55, CDC25C, SPC25, CCNB2, NCAPG2, BUB1B, CENPW, ASPM GO:0000070:mitotic sister chromatid segregation 7 4.24E-07 MAD2L1, NEK2, SPAG5, ZWINT, NUSAP1, KIF18B, NDC80 GO:0007062:sister chromatid cohesion 11 4.75E-07 SPC25, MAD2L1, SGO2, ZWINT, KIF18A, NUF2, BUB1B, NDC80, BIRC5, CDC20, CENPK GO:0007052:mitotic spindle organization 7 1.35E-06 CCNB1, SPC25, KIF11, PCNT, TTK, NDC80, STMN1 GO:0007019:microtubule depolymerization 5 4.11E-06 KIF14, STMN3, KIF18A, KIF18B, STMN1 GO:0045893:positive regulation of transcription, DNA templated 21 4.56E-06 SOX11, PAX6, TGFB3, ATAD2, LEF1, TBX1, CREB5, SOX9, GLI2, MDK, FZD7, GLI1, MYCN, SMARCD3, LHX2, ZNF711, TFAP2B, CAND2, RFX3, PTCH1, SOX18 GO:0030574:collagen catabolic process 8 1.18E-05 MMP10, COL6A3, COL6A2, COL6A1, ADAMTS3, COL11A1, COL5A2, MMP12 GO:0007059:chromosome segregation 8 1.77E-05 SPC25, KIF11, NEK2, SPAG5, NUF2, CENPW, NDC80, TOP2A GO:0006260:DNA replication 11 1.92E-05 CDK1, GINS2, POLE2, DTL, RRM2, BRIP1, CDC25C, MCM5, FEN1, MCM6, NFIB hsa04110:cell cycle 13 4.75E-09 CCNE2, CCNB1, CDK1, MAD2L1, CCNB2, GADD45G, TGFB3, TTK, BUB1B, CDC20, CDC25C, MCM5, MCM6 hsa04115:p53 signaling pathway 6 6.65E-04 CCNB1, CCNE2, CDK1, CCNB2, RRM2, GADD45G hsa04974:protein digestion and absorption 6 0.002273 COL6A3, COL6A2, COL6A1, COL11A1, COL5A2, DPP4 hsa05217:basal cell carcinoma 5 0.002604 LEF1, PTCH1, GLI2, FZD7, GLI1 hsa03030:DNA replication 4 0.006291 POLE2, MCM5, FEN1, MCM6 hsa04512:ECM-receptor interaction 5 0.013198 COL6A3, COL6A2, COL6A1, COL11A1, COL5A2 hsa04914:progesterone-mediated oocyte maturation 5 0.013198 CCNB1, CDK1, MAD2L1, CCNB2, CDC25C hsa05200:pathways in cancer 10 0.021829 CCNE2, TGFB3, RUNX1T1, LEF1, BIRC5, PTCH1, GLI2, FZD7, GNG7, GLI1 hsa04114:oocyte meiosis 5 0.027764 CCNE2, CDK1, MAD2L1, CDC20, CDC25C hsa04340:hedgehog signaling pathway 3 0.032592 PTCH1, GLI2, GLI1 [98]Open in a new tab Abbreviations: ECM, extracellular matrix; GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes. Table 3. The top 10 GO terms and KEGG pathways of downregulated genes Term Count P-value Genes GO:0048704:embryonic skeletal system morphogenesis 4 7.52E-04 HOXB2, HOXB7, HOXA5, HOXA6 GO:0036109:alpha-linolenic acid metabolic process 3 0.001566736 ELOVL5, FADS1, FADS2 GO:0006636:unsaturated fatty acid biosynthetic process 3 0.002096572 ELOVL5, FADS1, FADS2 GO:0043651:linoleic acid metabolic process 3 0.002699483 ELOVL5, FADS1, FADS2 GO:0001558:regulation of cell growth 4 0.005913176 MELTF, BCAR1, NANOS1, CYR61 GO:0009952:anterior/posterior pattern specification 4 0.005913176 HOXB2, HOXB7, HOXA5, HOXA6 GO:0007267:cell–cell signaling 6 0.006210487 BMP2, ADRB2, FADS1, AREG, GDF15, CYR61 GO:0055007:cardiac muscle cell differentiation 3 0.006763636 BMP2, SIK1, NRG1 GO:0060325:face morphogenesis 3 0.008308263 DKK1, TIPARP, RRAS GO:2000726:negative regulation of cardiac muscle cell differentiation 2 0.013694378 BMP2, DKK1 hsa01040:biosynthesis of unsaturated fatty acids 3 0.005255803 ELOVL5, FADS1, FADS2 hsa01212:fatty acid metabolism 3 0.02175659 ELOVL5, FADS1, FADS2 hsa05230:central carbon metabolism in cancer 3 0.037090419 SLC1A5, HKDC1, MYC [99]Open in a new tab Abbreviations: GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes. PPI network analysis and module analysis After data of interactions imported into Cytoscape software, the PPI network with 202 nodes and 1,245 edges was constructed. Based on this network, TOP2A (degree =64), CDK1 (degree =59), and CCNB1 (degree =54) were screened as the top three hub genes due to the higher degrees ([100]Figure 2). Subsequently, we performed module analysis of the whole network by the MCODE plug-in. Three modules were identified and created as subnetworks. In addition, pathway enrichment analysis of genes included in each subnetwork was performed, which revealed that DEGs in modules 1–3 were mainly associated with “cell cycle”, “extracellular matrix (ECM)-receptor interaction”, “basal cell carcinoma”, and “hedgehog signaling pathway” ([101]Figure 3). Figure 2. [102]Figure 2 [103]Open in a new tab Histogram of degrees of the top 30 genes in the protein–protein interaction network. Note: The number displayed on each column is the degree of each gene. Figure 3. [104]Figure 3 [105]Open in a new tab Three subnetworks obtained from the whole protein–protein interaction network. Notes: (A, B) Module 1 and the pathway enrichment analysis of genes in module 1. (C, D) Module 2 and the pathway enrichment analysis of genes in module 2. (E, F) Module 3 and the pathway enrichment analysis of genes in module 3. Vertical axis represents GO or pathway terms. P-values are displayed by gradient colors. Abbreviations: ECM, extracellular matrix; KEGG, Kyoto Encyclopedia of Genes and Genomes. Discussion BCC, with low malignancy, is the most common skin cancer worldwide. Although rarely metastasize, BCC can cause substantial local tissue damage along with disfigurement and involve other adjacent areas of soft tissue, cartilage, and bone.[106]^7 Currently, the targeted treatments of BCC implicated in clinical practice mainly focus on the hedgehog signaling pathway.[107]^21 However, the issue of drug resistance and poor response rate cannot be ignored. In order to explore more potential therapeutic targets, the gene expression profiles of BCC need to be comprehensively studied. In our present study, a bioinformatics approach was conducted to reanalyze the gene expression profiles of 19 BCC and 6 normal skin tissues. A total of 313 DEGs were identified with 222 upregulated genes and 91 downregulated genes. Functional and pathway enrichment analysis indicated that these DEGs were significantly associated with mitosis, cell cycle, and unsaturated fatty acid metabolic process. By PPI network and module analysis, three critical genes and four pathways were finally identified, which may play a key role in the carcinogenesis of BCC. With regard to functional and pathway enrichment analysis, upregulated DEGs were mainly involved in the process of mitosis and cell cycle. Deregulation of cell cycle is a common feature in the initiation and progression of various cancers, which is often mediated by alterations in cyclin and cyclin-dependent kinase (CDK) activity.[108]^22 CDK1, as a mitotic CDK, is sufficient to drive the mammalian cell cycle without other interphase CDKs.[109]^23 Accumulating evidences indicated that dysregulation of CDK1 activity was participated in a variety of tumors, including lung cancer,[110]^24 prostate cancer,[111]^25 and colorectal cancer.[112]^26 Schmit et al also discovered that increased level of CDK1 and CCNB1 presented in nonmelanoma skin cancer cells (BCC and squamous cell carcinoma) compared with normal human epidermal keratinocytes growth.[113]^27 Moreover, patched1, the BCC-related protein, was found to be interacted with cyclin B1 to regulate cell-cycle progression in BCC.[114]^28^,[115]^29 Recently, targeting cyclin-dependent kinases has become a promising approach in cancer therapy. AZD5438, as a highly specific inhibitor of CDK1, 2, and 9, was discovered to enhance the radiosensitivity of non-small-cell lung cancer.[116]^30 In the present study, our results revealed that CDK1 was significantly upregulated in BCC samples and enriched in many cell cycle-related GO terms, which indicated the potential to be a therapeutic target in BCC. Topoisomerases have been considered as important therapeutic targets for human malignancies. TOP2A, the major isoform of topoisomerase II, is capable of resolving catenanes and supercoils during DNA metabolic processes and plays a critical role in condensation and segregation of chromosomes at mitosis. Accumulating studies highlighted that higher TOP2A expression level was correlated to advanced tumor stage and poor patients’ survival in human cancers. At the protein level, increased expression of topoisomerase IIα was demonstrated to be associated with elevated cell replication in BCC compared with squamous cell carcinoma.[117]^31 In our study, TOP2A was screened as the most significant gene with the highest degree and was up-regulated in BCC. Elevated expression of TOP2A was implicated in cell cycle, and targeting TOP2A was also considered as an important therapy for human cancers.[118]^32 Thus, TOP2A could be a critical target in BCC. COL6A1, COL6A2, COL6A3, COL5A2, and COL11A1 are members of the collagen family, and these five genes are enriched in the pathway of “ECM–receptor interaction”, which leads to a direct or indirect control of cellular activities such as adhesion, migration, differentiation, proliferation, and apoptosis. Accumulating evidence indicated that the “ECM–receptor interaction” pathway served as a critical role in the carcinogenesis and metastasis of human cancers, such as prostate cancer,[119]^33 breast cancer,[120]^34 and colorectal cancer.[121]^35 In this study, we also screened “ECM–receptor interaction” as an important pathway by module analysis, which indicated the potential role in the pathogenesis of BCC. Hedgehog signaling pathway, a highly conserved evolutionary pathway of signal transmission from the cell membrane to the nucleus, has been revealed to be associated with the development of cancers, especially in BCC.[122]^5 The main downstream genes of hedgehog signaling pathway include PTCH1, GLI1, and GLI2. In the module 3 analysis, these three genes were significantly enriched in “basal cell carcinoma”, “hedgehog signaling pathway”, and “pathways in cancer”. Currently, targeting the hedgehog signaling pathway has been an important strategy for cancer therapy, which has achieved a promising success in BCC.[123]^21 However, the targeted genes were restricted to two genes (PTCH1 and SMO). Therefore, the other critical genes in this pathway are expected to be studied. Of note, several limitations also existed in our work. First, the inclusive criteria for BCC patients and normal controls was not available due to lack of data from the public database. Second, the same as most previous studies, two relatively small patient cohorts were performed in this study. Third, there was a lack of validation in biological experiments or another dataset, which might increase the FDR in our results. In conclusion, we performed a comprehensive bioinformatics analysis of DEGs obtained from 19 BCC and 6 normal skin tissues. Three hub genes and four pathways were finally identified, which might play a critical role in BCC. Our results further revealed the potential molecular mechanisms during the initiation of BCC and laid the foundation of exploring effective molecular targets for the treatment of BCC. However, future biology experiments are required to confirm these findings. Supplementary materials Table S1. Differentially expressed genes between basal cell carcinoma and normal skin tissues Upregulated genes ADAMTS3 LHX2 CHGA LGR5 SOX11 S100A9 PTCH1 FBN3 FAT3 TMSB15A KCNE1 MYCN CRNN COL11A1 MMP10 BNC2 GAS2 TOX2 SPON2 TFAP2B GLI2 HEPH LMO3 ADAMTS17 VASH2 LINGO1 DIO2 CHST2 PCDHB2 PCDH8 NPNT SOX18 PITX2 UHRF1 TBX1 CREB5 ABI3BP LINC00865 EDIL3 GPC4 SHCBP1 SLC6A1 SERPINB4 GJB6 APCDD1L SOSTDC1 LRRN1 VCAN BGN FZD7 SFRP5 TNRC6C MARCH1 MUM1L1 ZNF711 SHOX2 LOC101929122 F2RL2 DTL TSPYL5 CASC15 PELI2 NRTN GLI1 SETBP1 FNDC1 MEGF6 RAD51AP1 PAPPA SOBP HUNK NINL UCP2 HIST1H4C ADGRL3 CHRDL2 NAP1L3 NTRK3 TOP2A SOX9 TSPAN18 H2BFXP DLGAP5 MAD2L1 S100A8 PLEKHG4B NUF2 GMPR NDC80 LRIG3 SLC7A2 CENPK KRT85 ALDH1A3 MMP12 BIRC5 PCNT KALRN KCNS3 SDC2 CYFIP2 KIF11 COL5A2 CNTN4 GBP6 BACH2 HS3ST3A1 LEF1 SGO2 GINS1 CDH11 TM4SF1 KIF14 MARCKSL1 STMN1 LOC440173 PCDHB10 NEK2 APOBEC3A SMARCD3 IFI27 SH3GL3 SLC6A6 NUDT10 MDK CMPK2 APELA SHANK2 GADD45G RUNX1T1 TTK CCNB1 TET1 LTBP1 PBK KRT13 CDC20 ABCC12 DENND2A BEND5 ASPM NUSAP1 RRM2 CENPW DMD BRIP1 SLAIN1 SYNPO PTPRN2 NRXN3 STMN3 MXRA5 FANCD2 DCHS1 MCM6 CDC25C CDH22 COL6A2 SORCS2 DPP4 TIGD1 SPAG5 MTFR2 KIF15 KIF18A RFX3 PRIMA1 CCNB2 STON1 CDC42EP4 MCM5 ZWINT SEMA6A ZNF367 CEP55 HMCN1 TMTC1 BUB1B PABPC4L FANCI ZNF566 CAND2 POLE2 CDK1 STC1 ALMS1 PLCE1 NFIB PCSK5 PLPPR1 COL6A1 NCAPG BHLHE41 CCNE2 PAX6 MIR3682 PRRT2 GINS2 NCAPG2 CFAP44 COL6A3 FEN1 TNFSF10 PLEKHO1 TNS3 KIF18B ATAD2 TMEM173 ZNF853 SLCO2A1 TBC1D1 GNG7 DEPDC1 SPC25 OSBPL7 TAGLN NTN1 TGFB3 BICC1 IFI44 MKI67 LUM Downregulated genes CRTAP ID3 CNTN1 NRG1 SLC22A15 TIPARP MMP28 IDH1-AS1 KLHDC8B XG SMIM21 HOXB2 NEFL CRELD1 HKDC1 C11orf70 MYADM LPAR3 IL17RC FAM110A BMP2 CBS NUDT16P1 PHLDB2 IMPACT HOXB7 RTN4RL1 SLC1A5 RRAS SNORA4 GDF15 C8orf88 CDKN2AIPNL ESPN ADGRF4 KLF6 PTPN20 BCAR1 PRSS21 MOCOS PLPP2 LURAP1L MINDY2 MAFF ERRFI1 DLK2 CTH GNAI1 HOXA6 CORO2A C1orf21 MST1R EVPLL PYROXD2 CYR61 SYBU ELOVL5 HOOK2 SERPINB2 FAS AREG FAM89A HIST1H2BD IL13RA2 QPRT AP1S1 BTBD16 ACP5 MFSD2A ANTXR2 RAB5C NANOS1 CCPG1 C11orf63 FADS1 HIST1H2AC EML2 TSC22D3 DKK1 FKBP5 MYC ATF3 HOXA5 EN1 SIK1 ECM2 CHMP4C RNF128 FADS2 MELTF ADRB2 [124]Open in a new tab Acknowledgments