Transcriptome data from GSE49710 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE49710) and E-MTAB-8248 cohort (https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-8248/) in GEO and ArrayExpress databases was included in this study. The two cohorts included 498 and 223 patients respectively, and the main clinical information included age, sex, COG risk, MYCN status, INSS status, survival time, and survival status. The transcriptome data of 39 NB cell lines in GSE89413 cohort (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE89413) was also obtained. After data correction and probe annotation, we finally got the transcriptome matrix data with row names of gene and column names of patient number/cell number. In order to obtain more accurate anoikis-related subtypes and build an anoikis-related prognostic model, we merged the transcriptome matrix data of GSE49710 and E-MTAB-8248, standardized them through the “normalizeBetweenArray ()” function, and finally got the merged matrix data.

Args were obtained by searching key word “anoikis” in Genecards (https://www.genecards.org/, score>0.4) and Harmonizome (https://maayanlab.cloud/Harmonizome/) databases. Through univariate cox regression analysis of the combined cohort, we screened Args closely related to survival status of children with NB by setting the standard of “HR≠1 and p<0.05”.

We conducted consensus clustering analysis in 721 children by analyzing the expression levels of Args closely related to the survival status of NB. We finally obtained the best k-value for clustering by using the “CONSENSUSClusterPlus ()” function and setting the parameters reps=50, pItem=0.8, pFeature=1, clusterAlg=“km”, distance=“euclidean”, seed=123, and used the “prcomp ()” function to perform PCA analysis on all samples. Kaplan-Meier (KM) survival analysis of the anoikis-related subtypes was conducted by the “survdiff ()” function, and Args (|LogFC|>0.5 and p<0.05) differentially expressed in anoikis-related subtypes were screened through the “Lima” package. “Ggboxplot () “ and “ Pheeatmap () “ function were used to draw the boxplot and heatmap of differently expressed Args between different subtypes, and “ Pheeatmap () “ function was also used to observe the distribution of clinical indicators such as age, sex, COG risk, MYCN status and INSS stage between different subtypes.

We obtained the hallmark gene sets (c5.go.v7.4.symbols.gmt and c2.cp.Kegg.V7.4.Symbols.gmt) from MSigDB database, conducted GSVA and GSEA enrichment analyses between different subtypes through “GSVA” R package and “GSEA ()” function, and visualize the data through “pheatmap ()” and “gseaplot ()” function respectively.

Single-sample gene set enrichment analysis (ssGSEA) was used to analyze the differences of immune cell infiltration between different anoikis-related subtypes. We firstly scored the infiltration levels of 28 kinds of immune cells in each NB sample by using the “GSVA ()” function, then analyzed the differences of immune cell infiltration in different subtypes, and finally visualized the data by using the “ggboxplot ()” function. We used “CIBERSORT ()” function to score immune cell infiltration level for each sample in different risk groups, analyzed the differences in immune cell infiltration level between different risk groups divided according to the prognostic model, and visualized the data by “vioplot ()” function. Finally, we used the “corrplot ()” function to visualize the correlation between immune cells. The “Cor.test()” function was used to analyze the correlation between immune cells and riskscores.

By setting the standard of Height=2000, we used the “hclust ()” function to cluster all the samples and draw a clustering tree after filtering out the outliers. Based on the criterion of approximate scale-free network (R²>0.85), we used the “pick soft threshold ()” function to select an appropriate soft threshold to make the network more in line with the power law distribution, thus reducing the error and making the results more characteristic of biological data. We set 50 as the minimum number of genes in each module and 0.25 as the threshold of cutting height, used the dynamic cutting tree method to identify modules, combined the modules with high similarity and finally obtained multiple important modules. In order to get key modules, we evaluated the correlation between genes and samples by calculating gene significance (GS), defined module significance as the average GS of all genes in the module, calculated the correlation between the modules and anoikis, and selected the module with the highest correlation as the key module. By calculating the GS and module membership (MM) of each gene and setting the interception criteria of GS>0.5 and |MM|>0.80, we finally got the core genes of the key module. By overlapping the prognosis-related Args between two subtypes and the core genes screened by WGCNA, we finally screened the hub genes closely related to anoikis in NB.

All samples were randomly divided into training group and test group with equal sample size. Lasso regression analysis was performed on all Args related to the survival state in the training group through the “Glmnet ()” function. The minimum regularization parameter lambda (λ) and genes highly associated with anoikis-related subtypes were obtained by cross-validation (alpha=1). Multivariate cox regression analysis was then used to screen the core genes and the weight of each core gene through the “coxph ()” function. Riskscore of each sample was obtained according to the expression level and corresponding weight of core genes. The total sample group, training group and test group were divided into high-risk group and low-risk group according to the median riskscore value respectively. KM survival analysis was performed for the three groups by “survdiff ()” function. “pROC ()” function was used to calculate area under the curve (AUC) in three groups and evaluated the predictive capacity of this model. We used the “coxp ()” function to conduct multivariate cox regression analysis on the clinical information that might affect the prognoses of children with NB, and drew the nomogram through the “replot ()” function. The “calibrate ()” function and “plot ()” function were used to analyze the parameters of calibration curve and plot calibration curve respectively. The “survfit ()” function and “ggsurvplot ()” function were used to analyze the parameters of the cumulative hazard curve and plot the cumulative hazard curve respectively. The parameters of the decision curve were analyzed by “dca ()” function, and the decision curve was plotted by “ggplot ()” function.

We obtained the training set gene expression matrix file (GDSC2_Expr.rds) and drug processing information file (GDSC2_Res.rds) from Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/). The “calcPhenotype ()” function in the “oncoPredict” R package was adopted, and the criteria of removeLowVaryingGenes = 0.2 and minNumSamples = 10 were set for the prediction of drug susceptibility.

All data analyses in this study were completed by R 4.2.1 version, and p < 0.05 was considered statistically significant. The t test and one-way ANOVA were used to analyze the parametric data. The Wilcoxon test and Kruskal-Wallis test were used to analyze the nonparametric data of two and multiple independent samples respectively. Related R packages described above were downloaded from Bioconductor packages or R packages.

In order to obtain relatively complete clinical data and transcriptome data with a large sample size, we selected GSE49710 and E-MTAB-8248 cohorts from GEO and ArrayExpress databases. The two cohorts included 498 and 223 clinical samples respectively. Both groups contained complete transcriptome array data and Log2 conversion was performed. The GSE49710 cohort contained relatively complete information of age (> 18 months n=193, ≤18 months n=305), sex (boy n=287, girl n=211), COG risk (High risk n=176, Low risk n=322), MYCN status (amplified n=92, non amplified n=401), INSS stage (Stage I n =121, Stage II n=78, Stage III n=63, Stage IV n=183, Stage IVs n=53), survival state (alive n=393, dead n=105) and survival time. E-MTAB-8248 had complete information of age (> 18 months n=119, ≤18 months n=104), MYCN status (amplified n=46, non amplified n=177), INSS stage (Stage I n =29, Stage II n=39, Stage III n=36, Stage IV n=89, Stage IVs n=30), alive state (alive n=181, dead n=42) and survival time. By combining and standardizing transcriptome data in two cohorts, we obtained transcriptome microarray matrix data with 721 samples. The transcriptome data of 39 NB cell lines in the GSE89413 cohort were obtained and Log2 conversion was also performed for further study.

We searched the key word “anoikis” in Genecards and Harmonizome databases, and a total of 640 Args were obtained. Univariate cox regression analysis was conducted on the expression level of 640 genes and the survival state of the children in this study. A total of 283 Args closely related to the survival state were obtained according to the standard of HR≠1 and p<0.05. Ten Args with the greatest impact on the prognosis were: PARP1 (HR=8.51), ELAVL1 (HR=6.49), XRCC5 (HR=5.62), CSNK2A1 (HR=5.60), DAP3 (HR=4.96), MYO5A (HR=0.20), CDC42 (HR=0.24), ARHGEF7 (HR=0.25), PDPK1 (HR=0.27), ATF2 (HR=0.27). Ten Args with the most statistical difference were: NTRK1 (p=1.58e-35), ETV4 (p=2.80e-31), PARP1 (p=4.19e-30), CDC42 (p=7.44e-30), SKP2 (p=3.00e-29), CHEK2 (p=6.68e-29), CDKN3 (p=7.78e-29), HK2 (p=8.59e-27), ARHGEF7 (p=1.62e-26), KIF18A (p=3.48e-26). We showed 28 genes with the lowest p-value through the forest map ( Figure 2A ). According to the expression levels of these 283 genes, 721 samples were clustered into different anoikis-related subtypes by consistent cluster analysis. Through the clustering heatmap and consistent cumulative distribution function (CDF) graph, we found that the CDF slope was the smallest and the clustering effect was the best when the number of clusters k was equal to 2 ( Figures 2B, C ). Through principal component analysis (PCA), we found that the clustering method could well distinguish the samples into two subtypes ( Figure 2D ) and KM survival analysis suggested a worse prognosis for ARCluster B ( Figure 2E ). These results indicated that the two anoikis-related subtypes had significant differences in the prognosis, and aroused strong interests of us in the differential expression of Args and the differential distribution of clinical information between ARCluster B and ARCluster A. We further analyzed the differences of Args between the two subtypes and found 119 differentially expressed Args. The 10 Args with the largest value of |LogFC| were NTRK1 (LogFC=3.66), LMO3 (LogFC=2.47), TWIST1 (LogFC=2.46), BIRC3 (LogFC=2.02), HK2 (LogFC=1.94), CCR7(LogFC=-1.88), ITGA8(LogFC=-1.86), KIF18A(LogFC=1.82), E2F1(LogFC=1.78) and UBE2C(LogFC=1.77). 32 genes with the lowest p-value were shown by boxplot ( Figure 2F ). It could also be found that the expression levels of Args were different between ARCluster A and ARCluster B according to the heatmap ( Figure 2G ). The results also showed that ARCluster B with a poor prognosis contained more samples with age > 18 months(p<0.001), INSS stage III-IV(p<0.001), MYCN amplified(p<0.001) and high COG risk(p<0.001) (33, 34). However, ARCluster A with a better prognosis mainly included samples with age ≤18 months(p<0.001), INSS stage I/II/IVs(p<0.001), MYCN non amplified(p<0.001) and low COG risk(p<0.001)( Figure 2G ). The above results fully indicated that children with NB could be clustered into ARCluster A with a good prognosis or ARCluster B with a poor prognosis through univariate cox regression analysis and consistent cluster analysis, which was consistent with clinical indicators such as age, COG risk, MYCN status and INSS stage (33, 34). These results made us interested in the study of the important pathways and key genes between the two anoikis-related subtypes.

In order to explore whether there were differences in the pathways between ARCluster B and ARCluster A, we performed GSVA and GSEA enrichment analyses for the two subtypes. It was found that there were obvious differences in the pathways between ARCluster B and ARCluster A through GSVA enrichment analysis. According to “c5.go.v7.4.symbols.gmt”, the five pathways with the most obvious differences were as follows: GOCC_METHYLOSOME, GOBP_ACTIN_FILAMENT_ORGANIZATION, GOBP_DNA_STRAND_ELONGATION, GOBP_REGULATION_OF_ACTIN_FILAMENT_ORGANIZATION, GOCC_SPLICEOSOMAL_TRI_SNRNP_COMPLEX ( Figure 3A ). We found that the expression level of GOBP_REGULATION_OF_CELL_MATRIX_ADHESION pathway decreased significantly in ARCluster B. According to “c2.cp.kegg.v7.4.symbols.gmt”, the most obvious five pathways differentially expressed between two subtypes were: KEGG_BASE_EXCISION_REPAIR, KEGG_SPLICEOSOME, KEGG_NUCLEOTIDE_EXCISION_REPAIR, KEGG_DNA_REPLICATION, KEGG_RNA_DEGRADATION ( Figure 3B ). We also found that the expression level of “KEGG_CELL_ADHESION_MOLECULES_CAMS” pathway related to cell adhesion was significantly reduced in ARCluster B. GSEA enrichment analysis was further performed. According to “c5.go.v7.4.symbols.gmt”, the five pathways with the most obvious expression differences were GOBP_ADAPTIVE_IMMUNE_RESPONSE, GOBP_IMMUNE_RESPONSE_REGULATING_CELL_SURFACE_RECEPTOR_SIGNALING_PATHWAY, GOBP_POSITIVE_REGULATION_OF_LEUKOCYTE_CELL_CELL_ADHESION, GOBP_LEUKOCYTE_MEDIATED_IMMUNITY, GOBP_LYMPHOCYTE_MEDIATED_IMMUNITY ( Figure 3C ). According to “c2.cp.kegg.v7.4.symbols.gmt”, the five pathways contained KEGG_HEMATOPOIETIC_CELL_LINEAGE, KEGG_GRAFT_VERSUS_HOST_DISEASE, KEGG_CELL_ADHESION_MOLECULES_CAMS, KEGG_CELL_CYCLE, KEGG_CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION ( Figure 3D ). Through GSVA and GSEA enrichment analyses, it was found that there were A large number of differentially expressed pathways between ARCluster A and ARCluster B. Moreover, the results showed that the expression levels of the pathways related to cell adhesion were all reduced in ARCluster B, which was consistent with the phenomenon that malignant tumors with AR were more likely to break down the barriers of ECM and form metastasis. We analyzed the differences of the infiltration levels of immune cells between ARCluster A and ARCluster B through ssGSEA. Except for the increased level of Activated. CD4 T cell and Type2 T helper cell, the infiltration level of other immune cells in ARCluster B decreased significantly ( Figure 3E ). These results fully showed that there were obvious differences in the infiltration of immune cells between different anoikis-related subtypes of NB, and the changes of TME might play an important role in breaking down the barriers of ECM and forming AR.

We further screened the hub genes that might play important roles in anoikis-related subtypes through WGCNA. The remaining samples were clustered after removing outliers ( Figure 4A ). When the threshold value was 7, R^{2 =} 0.92, the gene expression was the most consistent with the scale-free network (35) ( Figure 4B ). Through the construction of scale-free network, we obtained a total of 18 modules closely related to anoikis and further combined the 18 modules into 14 modules according to the standard of Height=0.25 ( Figures 4C, D ). Among the 14 modules, MElightcyan, MEcyan, MEpink, MEyellow, and MEgrey were positively related to anoikis, while MEturquoise, MEblue, and MEmagenta were negatively related to anoikis. Among all modules, MEturquoise module had the largest number of genes and was the most closely related to anoikis (correlation index=-0.77) ( Figure 4E ). According to the screening criteria of GS>0,5 and |MM|>0.8, 93 core genes were screened from the MEturquoise module ( Figure 4F ). We further crossed the prognosis-related Args with the 93 core genes and finally obtained 10 anoikis-related hub genes ( Figure 4G ): BIRC5 (GS=0.62, MM=-0.81), BUB1 (GS=0.60, MM=-0.80), CDKN3 (GS=0.64, MM=-0.82), CENPF (GS=0.61, MM=-0.82), CHEK2 (GS=0.66, MM=-0.85), E2F1 (GS=0.64, MM=-0.83), MAD2L1 (GS=0.64, MM=-0.81), NTRK1 (GS=-0.67, MM=0.82), PLK1 (GS=0.65, MM=-0.85) and UBE2C (GS=0.60, MM=-0.81).The above results indicated that there were multiple modules affecting anoikis between two anoikis-related subtypes of NB, and the MEturquoise module had the most obvious inhibitory effect on the anoikis of NB. The screening of the hub genes in this module would play an important role in our subsequent research on the mechanism of AR in NB.

We had found that there were two anoikis-related subtypes with different prognoses in NB through the above cluster analysis. Therefore, we hoped to construct a prognostic model that could classify NB into anoikis-related subtypes with different prognoses through their transcriptome data. Firstly, we randomly divided 721 cases of NB into the training group and the test group with the same sample size, and constructed lasso regression analysis based on the expression levels of 283 Args and the prognoses of these children. We got a total of 26 key genes through further cross-validation (alpha=1). The multivariate cox regression analysis was carried out between the expression of these 26 key genes and the prognoses of these patients. Finally, we obtained 15 key Args for building the model and listed the location of each gene ( Figures 5A, B ): STK11, CCND1, PIK3CG, DAPK1, FADD, CPT1A, TWIST1, BAG1, RAD9A, MALAT1, NTRK3, MMP3, HTRA2, KIF18A, DAP3. The weight of each gene in this model was stated in Table 1 . By calculating the riskscore of each sample, we divided the samples of training group and test group into high-risk group and low-risk group according to the median value of each group. The results suggested that the expression levels of 15 Args participating in the construction of this model were different between high and low-risk groups ( Figure 5C ), and the riskscores in ARCluster B were significantly higher than those in ARCluster A ( Figure 5D ). Almost all of the samples in ARCluster B were classified into the high-risk group, and most of the samples in ARCluster A were classified into the low-risk group. Almost all the dead patients came from the high-risk group, while the majority of the surviving children belonged to the low-risk group ( Figure 5E ). We found that Args differentially expressed in ARCluster A and ARCluster B groups also had consistent differences between high-risk and low-risk groups, and the clinical indicators related to poor prognosis including ≥ 18 months, MYCN amplified, high COG risk, INSS stage III-IV were mainly enriched in the anoikis-related high-risk group ( Figure 5F ), which was consistent with our previous conclusions. According to the commonly used clinical indicators: Age, MYCN status, COG risk, INSS stage and the risk level calculated by anoikis-related model, we drew the nomogram to comprehensively evaluate the prognoses of children with NB. Different clinical indicators were given different points in the nomogram. The probabilities of the 1-, 3-, and 5-year survival were predicted by the accumulated scores. It was found that the scores of the two opposite states of other indicators were very close except for the risk level. The nomogram showed that the risk level was an independent diagnostic factor and synergistically predicted the survival probability of NB patients. ( Figure 5G ). KM survival analysis indicated that the prognoses of patients with high riskscores were poor (p<0.001 in three groups) whether in the total sample group, training group or test group ( Figure 5H ). AUC values at 1-year, 3-years, 5-years in the total sample group (0.947, 0.900, 0.917), training group (0.947, 0.933, 0.958), and test group (0.942, 0.873, 0.888) indicated that the model we built could well predict the prognoses of children with NB ( Figure 5I ). By drawing the decision curve ( Figure 5J ) and calibration curve ( Figure 5K ), we could also find that the anoikis-related model we built could well predict the 1,3,5-year survival rate of children with NB. The cumulative hazard curve showed that the overall survival (OS) risk of children in high-risk group increased gradually with the prolongation of survival time ( Figure 5L ). The above results suggested that this anoikis-related model could accurately reflect the prognoses of children with NB, which laid the foundation for the popularization of this model and the researches on the potential mechanisms and immune microenvironment of NB in high-risk group.

We found that the model constructed in this study could effectively divide the samples of GSE49710 and E-MTAB-8248 cohorts into high-risk group and low-risk group ( Figures 6A, B ), and the prognoses of children in high-risk group were worse than that in low-risk group ( Figures 6C, D ). Args differentially expressed between ARCluster A and ARCluster B were also differentially expressed between high and low-risk groups of the GSE49710 and E-MTAB-8248 cohorts according to the heatmap ( Figures 6A, B ). Clinical indicators closely related to poor prognoses of children: ≥18 months, MYCN amplified, high COG risk, INSS stage III-IV were also enriched in high-risk group of both cohorts ( Figures 6A, B ). The results of independent prognostic analysis suggested that sex, MYCN status, COG risk, INSS stage and anoikis-related riskscore could be used as independent prognostic indicators in the GSE49710 cohort ( Figure 6E ). Age, MYCN status and anoikis-related riskscore were independent prognostic indicators in the E-MTAB-8248 cohort ( Figure 6F ). 39 NB cell lines in GSE89413 cohort were also scored according to the constructed model. RD, NB-SD, COG-N-519, CHP-212, COG-N-561, IMR-05, SK-N-FI, COG-N-440, LA-N-6, RD, NB-SD, COG-N-519, CHP-212, COG-N-561, IMR-05, SK-N-Fi, COG-N-440, LA-N-6 and Nb-1 were the top ten NB cell lines with the highest riskscore, while COG-N-557nb, SK-N-SH, NB-1643, IMR-32, COG-N-471nb, NB-69, NMB, SMS-KAN, SK-N-BE (2) and SH-SY5Y were the top ten cell lines with the lowest riskscore ( Figure 6G ). The above results fully indicated that the anoikis-related model constructed by us could generally predict the prognoses of children with NB, and the scores obtained by this model could be used as an independent indicator to better predict the prognosis. The above results also suggested that there were generally two anoikis-related subtypes in the population and cell lines of NB, among which the subtype with high riskscore and poor prognosis was likely to be closely related to the AR. The classification of NB cell lines would help us study the mechanism of AR in NB.

Previous studies had shown that TME was closely related to the metastasis of tumor cells, and AR was an important mechanism for metastasis. Therefore, understanding the TME characteristics of anoikis-related subtypes of NB was very important to immunotherapy for anoikis-related subtype of NB with a poor prognosis. CIBERSORT algorithm was used to analyze the immune cell infiltration level of NB samples from different anoikis-related risk groups. The results showed that B cells memory, Plasma cells, T cells follicular helper, NK cells resting, and Neutrophils were highly expressed, while T cells CD4 memory resting, T cells gamma delta, Macrophages M2 and Mast cells activated were weekly expressed in anoikis-related subtype with high riskscore ( Figure 7A ). We also found that there were correlations between immune cells, especially positive relationship between B cells naive and T cells follicular helper (correlation index=0.37). But the negative relationship between B cells naive and Macrophages M0 was more obvious (correlation index=-0.51) ( Figure 7B ). The estimate algorithm was used to further analyze the differences of TME between different risk groups. The results showed that the estimate score of high-risk group was lower, and the content of stromal cells and immune cells in tumor tissue was lower than that of low-risk group ( Figure 7C ). By analyzing the correlations between 15 Args participating in the model construction and immune cells, we found that except STK11 and CCND1, the other 13 Args were all related to the infiltration levels of immune cells. In particular, the expression level of PIK3CG was positively correlated with the infiltration levels of T cells regulatory (Tregs), T cells follicular helper, T cells CD8, T cells CD4 naive, T cells CD4 memory activated, Plasma cells, Monocells, Macrophages M1, B cells naive, and negatively correlated with the expression levels of NK cells activated, Mast cells activated, Macrophages M2, Macrophages M0 ( Figure 7D ). By further analyzing the correlation between the infiltration levels of 23 kinds of immune cells and the riskscores of patients, we found that a total of 6 kinds of immune cells were related to the riskscore. B cells memory (R=0.13, p=0.0079), Plasma cells (R=0.23, p=1.8E-6) and T cells follicular helper (R=0.12, p=0.0013) were positively correlated with the riskscore, while T cells CD4 memory resting (R=-0.27, p=1.2E-8), Macrophages M2 (R=-0.12, p=0.0012) and Mast cells activated (R=-0.12, p=0.0016) were negatively correlated with the riskscore ( Figure 7E ). The above results fully demonstrated that TME and immune cell infiltration level were different between anoikis-related subtypes. The key genes involved in the construction of anoikis-related model were also related to the infiltration levels of immune cells. The differences of the TME and immune cell infiltration level might be important factors leading to the AR of NB, and the key genes for constructing the model might play important roles in the above mechanisms.

According to the results calculated based on the GDSC database, we predicted that the anoikis-related high-risk group and low-risk group had different drug sensitivities to 158 kinds of chemotherapeutic drugs. The low-risk group had higher sensitivity to most of the predicted drugs, while the high-risk group had higher sensitivity to only17 kinds of drugs. The 8 kinds of chemotherapeutic drugs with the highest drug sensitivity in high-risk group were: JQ1, GSK269962A, Doramapimod, BMS-754807, AZD8055, KU-55933, AZD6482 and Axitinib ( Figure 8 ). The above results indicated that these selected chemotherapeutic drugs might have better therapeutic effects on NB in high-risk group, and might have potential targeted therapeutic effects on the AR of NB.