In Silico Identification of Therapeutic Targets and Novel Drug Candidates for Malignant Peripheral Nerve Sheath Tumors

Background : Malignant peripheral nerve sheath tumors (MPNSTs) are an aggressive form of sarcomas with a poor prognosis and limited treatment options. Therefore, new therapeutic targets are urgently needed to identify novel drugs. Methods : Based on the Gene Expression Omnibus database, an integrated analysis was performed to identify differentially expressed genes (DEGs) in MPNSTs compared to neurofibromas (NFs). Then functional enrichment analyses, protein-protein interaction (PPI) network construction


Introduction
Malignant peripheral nerve sheath tumors (MPNSTs) are highly invasive cancers that account for approximately 10% of all soft tissue sarcoma; 50% of such tumors are associated with neurofibromatosis type 1 (NF1).These sarcomas generally have poor clinical outcomes, and are the leading cause of mortality and morbidity in adults with NF1.Complete surgical resection is still the main treatment method for MPNSTs due to the limited effectiveness of both chemotherapy and radiotherapy.However, a full resection may not be possible due to tumor size and location in many patients [1,2].
Several amplified genes that are important for the pathogenesis of MPNST have been identified.For example, the expression of epidermal growth factor receptor (EGFR) and erbB2 is stronger in MPNSTs than in neurofibromas (NFs).Holtkamp et al. [3] valuated the effects of erlotinib and trastuzumab drugs targeting EGFR and erbB2 in MPNST cell lines, and EGFR and erbB2 have emerged as potential targets for the treatment of patients with MP-NSTs.The progressive amplification of MET was found in the case of MPNST.NF1 ablated mice exhibit a strong MPNST phenotype without additional mutations, which are consistently sensitive to the highly selective MET inhibitor capmatinib [4].Frequent platelet-derived growth factor receptor alpha (PDGFRα) mutations in MPNST are typically associated with the co-amplification of KIT Proto-Oncogene (c-KIT).The tyrosine kinase inhibitor imatinib can inhibit the proliferation of MPNST cell lines in vitro, and it is known that imatinib can target PDGFRα and c-KIT.PDGFRα is a candidate for targeted therapy of MPNST [5].Furthermore, the insulin-like growth factor 1 receptor (IGF1R) pathway is also a potential therapeutic target for MPNST patients.Inhibition of IGF1R in MPNST cell lines using small interference RNAs or IGF1R inhibitors leads to a significant decrease in cell proliferation, invasion, and migration [6].Although these targeted drugs have proven to be effective in treating MPNSTs, there are currently no U.S. Food and Drug Administration (FDA)-approved drugs for the treatment of MPNSTs.Therefore, a better understanding of the initiation and progression is needed to identify new therapeutic strategies.
In this study, we screened druggable targets via integrated bioinformatics analysis and predicted therapeutic drugs for MPNSTs using the Library of Integrated Network-Based Cellular Signatures (LINCS) database, which identi-fies the connection between gene expression profiles and FDA-approved drugs.Furthermore, molecular docking was used in the search for matching predicted potential drugs and screened targets.Our finding expands the current understanding of MPNST progression and identifies potential therapeutic drugs for NF1-related MPNSTs.

Data Collection and DEG Acquisition
The gene expression profiles of the microarray datasets GSE41747, GSE66743, and GSE141439 were obtained for analysis from the Gene Expression Ominbus (GEO) database at the National Center for Biotechnology Information (http://www.ncbi.nlm.nih.gov/geo/).Bioinformatics analyses of the data from online databases were carried out by R 4.1.1software (R Foundation for Statistical Computing, Vienna, Austria).We performed differential expression analysis using the limma package (https://bi oconductor.org/packages/release/bioc/html/limma.html) to determine differentially expressed genes (DEGs) with the criteria of |log2FC| >1 and p < 0.05 and created volcano plots of the DEGs [7].A Venn diagram and heatmap of DEGs were generated using the online database SangerBox website (http://sangerbox.com/tool.html).

Protein-Protein Interaction Network Construction, Gene Ontology, and Kyoto Encyclopedia of Genes and Genomes Pathway Analysis of DEGs
Search Tool for the Retrieval of Interacting Genes (STRING, https://string-db.org/) is an online tool for evaluating the information on protein-protein interaction (PPI) [8].We use the STRING to find possible relations among the DEGs.The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of the DEGs were conducted with the online database SangerBox website.

PPI Network Module Analyses and Identification of Hub Genes
PPI network and module analyses were visualized using Cytoscape version 3.7.2(Institute for Systems Biology, Seattle, WA, USA), and subsequent hub gene identification of the modules was performed [9].First, the Cytoscape Molecular Complex Detection (MCODE) plug-in was used to identify the most important modules in the PPI networks.The criteria for selection were set as follows: MCODE score >3, number of nodes >4.Next, we use the Cyto-Hubba plug-in of Cytoscape to filter hub genes from the PPI network with the EcCentricity, BottleNeck, and Closeness ranking methods [8].The top 10 genes were considered hub genes.

Prognosis Analysis and Differential Expression of Twist Family bHLH Transcription Factor 1 in Pan-Cancer
We used the Kaplan-Meier method and Cox proportional hazards (CoxPH) regression model to explore the overall survival (OS) of patients with high and low Twist family bHLH transcription factor 1 (Twist1) expression in 39 cancer types.Hazard ratios (HRs) with 95% confidence intervals (CIs) and log-rank p values were reported.p < 0.05 was considered statistically significant.The online database SangerBox website (http://sangerbox.com/tool.html) was used to analyze the expression of Twist1 in different cancers and normal tissues [10].The expression profiles of Twist1 in different cancer and normal cell lines were obtained from the BioGPS database (http://biogps.org)[11].

Mutation Profiles
The c-BioPortal for Cancer Genomics (http://cbioport al.org) is a publicly available database for exploring multidimensional cancer genomics datasets [12].The genetic alteration of Twist1 in different cancers was analyzed by the c-BioPortal database [13].

MPNST-Associated Drug Prediction
We use The Connectivity Map, which is designed to identify biological pathways and expression modules by querying the gene lists of up-and down-regulated genes to rank small molecules, drugs and genes.A so-called connectivity score was used to evaluate the correlation between the drugs and input genes; a positive score indicates an inducement effect of the compound on the gene signatures, while a negative score reflects an inverse impact of a compound on the gene signatures.The L1000 database of the LINCS project includes 476,251 gene expression signatures gathered from 72 cell lines stimulated by 27,927 small molecule compounds.We identified the potential compounds based on the connectivity score (-2 to 2).

Molecular Docking Between MPNST Candidate Drugs and Hub Genes
The 3D structure files of the compounds were downloaded from PubChem (https://pubchem.ncbi.nlm.nih.gov/).The crystal structures of target proteins were retrieved from the RCSB Protein Data Bank (PDB) (https://www.rcsb.org/).AutoDockTools 1.5.7 software (Molecular Graphics Laborary, The Scripps Research Institute, La Jolla, CA, USA) was used to prepare the required files, converting ligands and receptors structures from PDB file format to Pdbqt format.All water molecules were detached for predocking preparation, and then Kolmman charges and polar hydrogen atoms were added to the compounds using Autodock-Tools.Molecular docking calculations were performed by AutoDockTools.The docking results were expressed as binding affinity values of the obtained receptor/ligand complex (kcal/mol) according to hydrogen bonds, electrostatic, and hydrophobic interactions.PyMOL 2.6.0 software (PyMOL Molecular Graphics System, Schrödinger, LLC, Shanghai, China) was used to visualize the 3D structures of the protein and major receptor/ligand interaction in terms of hydrogen bonds binding to the key amino acid residues.

Statistical Analysis
The primary statistical analyses were performed in R 4.1.1software.The normality test was performed with the Shapiro-Wilk test.The Kaplan-Meier curve was drawn to compare time to event distribution, and Univariate Cox regression analysis was used to evaluate the correlations between the DEGs and OS.All hypothetical tests were twosided, and p < 0.05 was considered statistically significant with the following thresholds: *p < 0.05, **p < 0.01, ***p < 0.001.

Functional Enrichment Analyses
The GO enrichment analysis showed that biological process (BP) terms were related to the cellular process, single-organism process, metabolic process, and biological regulation.GO_cellular component (GO_CC) analysis showed enrichment in macromolecular complexes, membranes, and membrane-enclosed lumen.GO_molecular function (GO_MF) analysis showed enrichment in catalytic activity, DNA binding transcription factor activity, and structural molecular activity (Fig. 2A).The KEGG pathways were significantly enriched in antifolate resistance, proteoglycans in cancer, viral carcinogenesis, and the hypoxia inducible factor 1 (HIF-1) signaling pathway, which involved in multipe mechanism of tumorgenesis and related to occurrence and development of MPNSTs (Fig. 2B,C).

PPI Network Construction, Module Selection, and Identification of Hub Genes
The PPI network was constructed on STRING, which contained 287 nodes and 848 edges.The most significant module was obtained in the PPI network of the upregulated and downregulated genes using Cytoscape MCODE plugin.It was enriched in viral carcinogenesis, proteoglycans in cancer, and necroptosis (Fig. 3A).

Prognostic Value and Differential Expression of Twist1 across Pan-Cancer
We analyzed Twist1 expression in various cancer cell lines and normal tissues extracted from the BioGPS database.The results revealed that Twist1 expression increased not only in almost all cancer cells but also significantly varied among different cancer cell lines (Fig. 3B).Subsequently, we examined the expression levels of Twist1 in pan-cancer using the online database SangerBox.The expression of Twist1 was significantly increased in adrenocortical carcinoma (ACC), lower grade glioma (LGG), cholan-giocarcinoma (CHOL), lung adenocarcinoma (LUAD), acute myeloid leukemia (LAML), stomach adenocarcinoma (STAD), head and neck squamous cell carcinoma (HNSC), testicular germ cell tumors (TGCT), kidney renal clear cell carcinoma (KIRC), lung squamous cell carcinoma (LUSC), rectum adenocarcinoma (READ), glioblastoma multiforme (GBM), pancreatic adenocarcinoma (PAAD), prostate adenocarcinoma (PRAD), and uterine carcinosarcoma (UCS) (Fig. 3C).These results suggested that Twist1 is overexpressed in tumor tissues.The univariate Cox regression analysis was used to investigate the correlation between Twist1 expression and prognosis, and a forest plot was drawn via R package "survival" (Fig. 3D ).Furthermore, high Twist1 expression was also associated with poor OS time in MPNSTs (Fig. 3E).Together, these data suggested that Twist1 expression was associated with prognosis in most cancers, including MPNSTs.

Twist1 Genomic Alterations in Different Cancer Groups
The results demonstrated that Twist1 genomic alterations occurred in 1.5% of patients (Fig. 4A).The frequency of Twist1 gene alterations was higher in different cancer types, especially in the type of amplification (Fig. 4B).

MPNST-Associated Drugs
Overlapping DEGs between MPNST and NF were used to search for drugs that could revert their expression signatures using the LINCS Query tool.Drug-gene combinations were ranked by search score.Drugs with low scores correspond to higher reversal potency and more significant potential for application.Finally, we selected five drugs (cytochalasin-d, cabozantinib, everolimus, refametinib, and BGT-226).

Molecular Docking Simulation
To know the interaction between predicted five drugs and Twist1, the molecular docking simulation study was performed by the AutoDockTools software.At the same time, we selected selumetinib, a classic mitogen-activated protein kinase (MAPK) inhibitor, as the reference drug for drug identification.The binding energy of cytochalasind, cabozantinib, everolimus, refametinib, BGT-226, and selumetinib with Twist1 was -6.61, -7.03, -7.73, -3.94, -7.07, and -6.22 kcal/mol, respectively (Table 2, Fig. 5).The molecular docking results revealed that everolimus interacted with Twist1 via two H-bonds in close proximity of 2.1 Å and 2.2 Å with ARG-423 and TYR-430 residues, respectively.

Discussion
Only 8-13% of NF1 patients develop MPNSTs.The genetic alteration from NFs to MPNST formation is unclear [14].Cyclin-dependent kinase inhibitor 2A (CDKN2A)/p16 and polycomb repressive complex 2 (PRC2) gene mutations drive the transformation from NFs to MPNST [15][16][17][18][19][20].p16 is a tumor suppressor protein encoded by the CDKN2A gene, which inhibits cyclin D-dependent protein kinases from entering the S phase of the cell cycle, thereby maintaining retinoblastoma protein (Rb) in a lowphosphorylated state and preventing its isolation from E2F transcription factor.Loss of p16 expression and homozygous deletion of the CDKN2A gene in NFs may be related to their transformation to MPNSTs [15].The core PRC2 complex comprises four components: EZH1/2, SUZ12, EED, and PbAp46/48 [21].Lee et al. [17] identified loss-offunction somatic alterations of the /textitPRC2 components (EED or SUZ12) in MPNSTs.MPNSTs with PRC2 inactivation have been shown consistent and complete loss of trimethylation at lysine 27 of histone H3 (H3K27me3).The introduction of the missing PRC2 component in a PRC2deficient MPNST cell line restored H3K27me3 levels and reduced cell proliferation.Some individuals with NF1 will develop MPNSTs that arise due to the malignant transformation of plexiform neurofibroma (pNF) [22].In the transition from NFs to MP-NSTs, a recently identified tumor is called Atypical neurofibromatosis neoplasms of uncertain biological potential (ANNUBP), which is hypothesized to be precursor lesions of MPNSTs [23].However, the molecular mechanism underlying the transformation from pNF to MPNSTs is still unknown.In this study, we first identified DEGs between MPNSTs and NFs using three independent GEO databases.In total, 34 upregulated DEGs and 55 downregulated DEGs were revealed.Bioinformatics analysis revealed that the upregulated genes are mainly involved in proteoglycans in cancer, and viral carcinogenesis, and the HIF-1 signaling pathway.To identify the hub genes involved in the progres-sion from NFs to MPNSTs, we constructed a PPI network, and Twist1 was identified as a hub gene in the progression of MPNSTs.Further analysis revealed that higher expression of Twist1 correlated with shorter survival time in MPNSTs as in pan-cancer.Thus it may be a potential biomarker for the prognosis of MPNSTs.
Cutaneous NFs are the most prevalent tumour in NF1 patients.Plexiform NFs, defined as NFs involving multiple peripheral nerve fascicles, are less frequent but lead to serious consequences.Malignant transformation develops in 5-10% of NF1 patients, typically in a plexiform NFs subtype [24].The neural crest transcription (SOX9) is found to be a biomarker of MPNSTs, which is possibly a therapeutic target in NF1.Interestingly, <30% of cells in most cutaneous NFs are SOX9 positive; plexiform NFs contain an average of 50% of SOX9 positive cells, while in MPN-STs, >70% of cells are SOX9 positive.This may indicate a potential evolutionary relationship among the three types of tumors [25,26].
It appears that plexiform NFs have a higher ability to evolve into MPNSTs than cutaneous NFs.But so far, it is still unknown whether there is a molecular difference between cutaneous NFs and plexiform NFs.In the GSE41747 dataset we used to analyse, transcripts that passed the ANOVA failed to distinguish the two NF subtypes [25,26].The explanation is that cutaneous and plexiform NFs cells may be essentially the same but exposed to different tumor microenvironments.In addition, the other two datasets we used (GSE66743 and GSE141439) did not specify whether NFs were cutaneous NFs or plexiform NFs [27,28], so our analysis did not distinguish the two types of NFs.Due to the different potential of developing MPNSTs of two types of NFs, it is necessary to conduct gene differential analysis between different types of NFs and MPNSTs in future studies.
Twist1, a transcription factor and member of the basic helix-loop-helix protein family, is known to stimulate epithelial-mesenchymal transition (EMT) and promote tumor invasion and metastasis.Dysregulation of Twist1 is associated with the loss of epithelial markers, such as Ecadherin, and alpha/gamma-catenins, and the activation of mesenchymal markers, such as N-cadherin, suggesting that Twist1 promotes the EMT [29].Twist1-induced EMT is clinically associated with distant cancer metastasis and poor prognosis in multiple tumor types.Twist1 transcription factors mediate EMT target genes involved in cell migration and invasion [30], multidrug resistance [31][32][33], cancer stems cell self-renewal [34,35], immune surveillance [36], and apoptosis [37].
The widespread involvement of Twist1 in human cancer suggests that inhibition of Twist1 is a promising therapeutic strategy for the prevention and treatment of cancer.We also showed that Twist1 expression was linked to clinical prognosis in multiple cancers.In addition, high expression levels of Twist1 had a significant association with poor OS in MPNSTs.Twist1 is highly expressed in MPNSTs and has low expression in normal adult tissues, making it an promising therapeutic target for drug development.Therefore, the application of Twist1 inhibitor has great clinical value for the treatment of MPNSTs, with few side effects.
We used LINCS datasets to identify novel MPNSTs drugs based on the overlapping DEGs.Finally, we selected five drugs (cytochalasin-d, cabozantinib, everolimus, refametinib, and BGT-226) with lower scores regarding their potential as drugs.To further predict the five drugs that might be direct Twist1 inhibitors, the molecular docking analysis was performed using the AutoDockTools and PyMOLsoftware.The results showed that everolimus interacting with Twist1 has lowest binding energy, making everolimus an optimal potential drug.
The mammalian target of rapamycin (mTOR), which involve in the growth and proliferation of tumor cells and play a role in the switch between cell catabolism and anabolism, is dysregulated in many types of tumors [67].This explains the potential use of mTOR inhibitors in oncology.Everolimus is a new mTOR inhibitor, which inhibits the growth and angiogenesis of tumor cells by continuously inhibiting mTOR.Everolimus blocking the bind-ing of the raptor to mTOR can restore control of the activated PI3K/AKT/mTOR signalling pathway.Everolimus has been approved for treating advanced breast cancer with positive hormone receptors and negative human epidermal growth factor receptor 2 [68], advanced renal cell carcinoma after targeted treatment of vascular endothelial growth factor [69], and advanced neuroendocrine tumors of pancreatic, gastrointestinal or pulmonary origin [70,71].In addition, it can be used to control the growth of subependymal giant cell astrocytoma, renal angiomyolipoma, and hamartoma in tuberous sclerosis complex (TSC) [72].It can also be used for immunosuppression after heart or kidney transplantation [73].
Interestingly, Twist1 promotes energy metabolism reprogramming (EMR) of glucose metabolism in breast cancer cells by activating the PI3K/AKT/mTOR pathway.EMR is crucial for the survival of cancer cells, as it can increase cancer proliferation, migration, and invasion [74].In addition, previous studies have confirmed that Twist1 promotes tumor invasion and metastasis by inducing EMT, and multiple signals that induce EMT often converge in the PI3K/AKT/mTOR signalling pathway.Xue et al. [75] found that complete EMT phenotype and morphological changes require PI3K/AKT mediated Twist1 phosphorylation in breast tumor cells.Twist1 may be a PI3K/AKT/mTOR signalling pathway member.Based on our research, everolimus could be a potential therapeutic drug for MPNSTs.The results of molecular docking show that everolimus has the lowest binding energy with Twist1, indicating that it may play an anti-tumor role by targeting Twist1 in addition to targeting mTOR.Another possible explanation is that it indirectly regulates Twist1 phosphorylation by targeting PI3K/AKT/mTOR pathway.However, our study results have not been experimentally validated and further experiments are needed.
This study is the first to systematically investigate the potential role of Twist1 in MPNSTs, which may be a key gene related to the occurrence and progression of MPNSTs.In addition, through molecular docking, everolimus, a potential target drug for MPNSTs, has a clear binding site and good binding ability with Twist1.It is expected to provide new ideas and insights for exploring the treatment strategies and therapeutic drugs of MPNSTs.

Conclusions
Overall, our results describe the importance of Twist1 in MPNST pathogenesis.everolimus was also found to be a potential therapeutic drug for MPNSTs.These findings require further experimental studies for confirmation.

Limitation
This study had limitations.First, we have not performed further experiment methods to verify the functional roles of Twist1.Therefore, functional experiments including cell proliferation, migration, gain-and loss-of-function assays, should be conducted in the future.Second, although we have identified the clinical significance of Twist1 in MP-NSTs in three independent datasets, verification in large clinical samples is still needed.

Fig. 1 .
Fig. 1.Identification of overlapping DEGs.(A-C) Differential expression of genes between MPNST tissues and neurofibroma tissues in the microarray datasets (GSE41747, GSE66743, and GSE141439).(D) Venn diagram showing DEGs across the three datasets, with an overlap of 34 upregulated gene.(E) Venn diagram showing DEGs across the three datasets, with an overlap of 55 downregulated genes.(F) The heat map of overlap upregulated and downregulated genes.The columns and rows in the heat map represent GEO datasets and genes.MPNST, malignant peripheral nerve sheath tumors; DEGs, differentially expressed genes; GEO, Gene Expression Ominbus.

Fig. 2 .
Fig. 2. Functional enrichment analyses of DEGs.(A) GO_MF analysis showed enrichment in catalytic activity, DNA binding transcription factor activity, and structural molecular activity.(B-C) The KEGG pathways were significantly enriched in antifolate resistance, proteoglycans in cancer, viral carcinogenesis, and the HIF-1 signaling pathway, which involved in multipe mechanism of tumorgenesis and related to occurrence and development of MPNSTs.GO_MF, GO_molecular function; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Fig. 3 .
Fig. 3. Prognostic value and differential expression of Twist1 across pan-cancer.(A) The protein-protein interaction network and hub genes.(B) Twist1 expression in different cancer cell lines extracted from the BioGPS database.(C) Twist1 expression in pan-cancers compared with normal tissues in the GEPIA database.(D) The correlation between the Twist1 expression and prognosis in pan-cancer.(E) High Twist1 expression was also associated with poor OS time in MPNST.

Fig. 4 .
Fig. 4. GZMA Genomic Alterations in different cancer groups analyzed by the cBioPortal database.(A) Twist1 genomic alterations occurred in 1.5% of patients.(B) The Twist1 gene alterations frequency was higher in different cancer types, especially in the type of Amplification.