{'fraction_non_alphanumeric': 0.04863544163747004, 'fraction_numerical': 0.03019546376544348, 'mean_word_length': 5.230040206777714, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 2, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A long non-coding RNAs expression signature to improve prognostic prediction of Wilms tumor in children

Background Wilms tumor (WT) is the most frequent malignancy of the kidney in children, and a subset of patients remains with a poor prognosis. This study aimed to identify key long non-coding RNAs (lncRNAs) related to prognosis and establish a genomic-clinicopathologic nomogram to predict survival in children with WT. Methods Clinical data of 124 WT patients and the relevant RNA sequencing data including lncRNAs expression signature of primary WT samples were obtained from the Therapeutically Applicable Research to Generate Effective Treatment (TARGET) Data Matrix. Then, lncRNAs associated with overall survival (OS) were identified through univariate Cox, least absolute shrinkage and selection operator (LASSO), and multivariate Cox regression analyses. The risk scores of 124 participants were calculated, and survival analyses were performed between low- and high-risk groups. A genomic-clinicopathologic nomogram was then developed and evaluated by time-dependent receiver operating characteristic (ROC) curves, including the area under the curve (AUC), calibration curve, and decision curve analysis. Subsequently, bioinformatics analyses were performed to explore the potential molecular mechanisms that affect the prognosis of WT. The package “DESeq2” was used to identify differentially expressed protein-coding genes (DEPCGs) between groups. Gene Set Enrichment Analysis (GSEA) was applied to explore the differences in pathways enrichment. The analytical tools CIBERSORTx and ESTIMATE were used to investigate the discrepancies of the immune microenvironment. Results A total of 10 lncRNAs were selected as independent predictors associated with OS (P<0.05). Participants in the high-risk group had a significantly worse OS and event-free survival (EFS) than those in the low-risk group (P<2E-16 and P=2.03E-04, respectively). The risk score and 3 clinicopathological features (gender, cooperative group protocol, and stage) were identified to construct the nomogram (combined model) (P=5.11E-17). The combined model (1-year AUC: 0.9272, 3-year AUC: 0.9428, 5-year AUC: 0.9259) and risk score model (1-year AUC: 0.9285, 3-year AUC: 0.9399, 5-year AUC: 0.9266) displayed higher predictive accuracy than that of the other models. Subsequently, 105 DEPCGs were identified. The GSEA revealed 4 significant pathways. Analysis with CIBERSORTx demonstrated that monocytes, macrophages M1, activated dendritic cells, and resting mast cells had significant infiltration differences between groups. Conclusions This study constructed a genomic-clinicopathologic nomogram, which might present a novel and efficient method for treating patients with WT.


Introduction
Wilms tumor (WT), also known as nephroblastoma, is the most frequent malignancy of the kidney in children, constituting nearly 90% of childhood renal tumors (1). Despite the development and advancement of various treatment approaches, 15% of WT patients undergo metastasis or relapse, suggesting poor prognosis after initial treatment (2). Therefore, it is of great clinical benefit to explore prognostic factors in children with different WT stages, which could provide valuable guidance on formulating an accurate therapeutic regime for each patient. Previous studies that were committed to identifying potential prognostic markers for WT patients relied on very few factors, mainly histology analysis and characterization of tumor stage, which has been deemed insufficient as WT in children was shown to be affected by a variety of factors, including genetic and epigenetic changes that underlie WT pathogenesis (1). The prognosis for WT children in clinical practice incorporates a diverse set of measurements, including age, tumor size, loss of heterozygosity (LOH) for chromosomes 16q and 1p, and sensitivity to medications (3)(4)(5)(6). Moreover, radiation (7), surgery, microscopic residual disease, diffuse anaplasia (8), lymph node involvement (9), and a combination of lymph node and LOH status (10) have been associated with the prognosis of WT in children. Therefore, it is important to build a prognostic model using multiple variables that underlie a broad spectrum of genetic and clinicopathological factors.
Long non-coding RNAs (lncRNAs) are characterized as non-coding transcripts longer than 200 nucleotides (11). Nextgeneration and high-throughput sequencing techniques have enabled a significant breakthrough in lncRNA identification and characterization in recent years. Gene expression of various lncRNAs is dysregulated in numerous cancers, and it has been reported to correlate with cancer recurrence, metastasis, and poor prognosis (12)(13)(14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24). Previous studies have suggested that several lncRNAs, such as LINC00473 (25), MIAT (26), SOX21-AS1 (27), and LINP1 (28), play important roles in the pathogenesis of WT. However, it remains largely unknown whether lncRNAs could be used as prognostic markers of WT in children. The present study sought to characterize lncRNA signature as a robust prediction for WT children's prognosis and construct a prognostic risk index based on comprehensive RNA-sequencing analysis. We established a prognosis model using RNA-sequencing data and subsequently developed and validated a genomic-clinicopathological nomogram that integrated risk scores and traditional clinicopathological factors. Moreover, Gene Set Enrichment Analysis (GSEA) revealed molecular pathways significantly enriched in different risk groups (29,30). Additionally, CIBERSORTx and ESTIMATE algorithms were applied to explore the immune microenvironment discrepancies between the high-risk and low-risk groups (31)(32)(33). Altogether, the current study presents a novel and efficient method for the prognosis of patients with WT and provides insight into the molecular mechanisms that affect WT's prognosis. We present the following article in accordance with the TRIPOD reporting checklist (available at http://dx.doi.org/10.21037/tp-20-318).

Study population and RNA-sequencing data processing
The level 3 RNAseq data of primary WT samples, including raw read counts and reads per kilobase million (RPKM), were obtained from the Illumina HiSeq 2000 platform (Illumina, San Diego, CA, USA), and relevant clinical information of the WT patients was downloaded from the Therapeutically Applicable Research to Generate Effective Treatment (TARGET) Data Matrix in April 2020. This study was conducted following the TARGET publication guidelines. The inclusion criteria of this study were as follows: (I) pathological diagnosis of primary WT; (II) clinical information including age at diagnosis, gender, cooperative group protocol, stage, histologic classification, the first event, event-free survival (EFS), vital status, and overall survival (OS) was documented; and (III) complete RNAseq data. The exclusion criteria were as follows: (I) patients with a history of other malignancies; (II) tumor tissues were not used for RNA sequencing analysis; (III) clinical information not available. Thus, 124 WT patients were selected and enrolled in this study (Tables S1,S2).
Participant clinicopathological characteristics are listed in Table 1. The geneset annotation GENCODE Human (GRCh38.p13) was used to annotate protein-coding genes and lncRNAs were selected for further study (34). The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Identification and validation of the prognostic lncRNA signature
The "survival" package in the R software (version 3.6.2; https://www.R-project.org/) was used to analyze RPKM data [transformed to log2(RPKM+1)] through univariate Cox regression analysis, and P<0.05 was used as the cut-off criterion. Next, through the "glmnet" package (version 3.0-2) in R, the least absolute shrinkage and selection operator (LASSO) Cox regression was performed to select key lncRNAs among the significant lncRNAs from univariate Cox analysis (35,36). We used 10 times cross-validations to determine the best penalty parameter lambda. The key lncRNAs selected by the LASSO method were considered as candidate variables for inclusion in the multivariate Cox analysis. Then, the independently prognostic lncRNAs were identified in the multivariate Cox analysis to obtain their coefficients (β values), with which the prognostic formula was constructed as follows: (β 1 × log 2 (RPKM+1) value of gene1) + (β2 × log 2 (RPKM+1) value of gene 2 ) + ... + (β n × log 2 (RPKM+1) value of genen). Subsequently, the risk score of each WT patient was calculated based on the prognostic formula. Then, the "surv_cutpoint" function in the "survminer" R package (version 0.4.6) was utilized to achieve the optimal cutoff values for risk scores. According to the cutoff values, we divided the 124 WT participants in the cohort into a low-risk group and a high-risk group. Between the low-and high-risk groups, differences in OS and EFS were examined by log-rank test, and Kaplan-Meier survival curves were constructed by the "survival" package in R software. A P value <0.05 was considered statistically significant.

Development and assessment of a predictive nomogram
Univariate and multivariate Cox regression analyses were conducted to ascertain whether the risk score of lncRNAs and the clinicopathological variables could be prognostic markers related to OS for WT patients, and a P<0.05 was deemed significant. Then, according to the results of Cox regression analyses, we constructed a genomic-clinicopathologic nomogram to predict 1-, 3-, and 5-year OS through the "rms" package (version 5.1-4) in R software. Using the same R package, we generated a calibration curve to verify the consistency between the nomogram-predicted probability of OS and the actual OS. Furthermore, to assess the accuracy of the predictive nomogram, we performed time-dependent receiver operating characteristic (ROC) curve analyses and calculated the area under the curve (AUC) through the "timeROC" package (version 0.4) in R (37). Additionally, we conducted a decision curve analysis (DCA) to examine the clinical value with the "stdca.R" statistical code (38) in R.

Identification of differentially expressed genes and construction of protein-protein interaction network
Using the "DESeq2" package (version 1.26.0) in R, the differentially expressed genes, including protein-coding genes (DEPCGs) and lncRNAs (DELs) were identified by comparing the low-risk group and high-risk group of primary WT, with the screening conditions of |log2 fold change| >1, and a false discovery rate (FDR) or adjusted P-value <0.05. These DEPCGs were then used to construct the protein-protein interaction (PPI) network through the Search Tool for the Retrieval of Interacting Genes (STRING) online database (39). Required interaction score >0.4 was used as the cut-off criterion. Then, the data downloaded from STRING were used to model the PPI network through Cytoscape software (version 3.7.1, https:// cytoscape.org) (40), and the densely connected regions were identified through the Molecular Complex Detection (MCODE) plug-in Cytoscape with the default criteria (degree cutoff =2, node score cutoff =0.2, Haircut = true, Fluff = false, K-score =2, and Max depth =100) (41).

GSEA
Before running GSEA, the normalization of the raw read counts of protein-coding genes was performed by "DESeq2". Then, GSEA (version 4.0.3) was conducted to demonstrate significant differences in protein-coding genes between the low-risk and high-risk groups (29,30 (42)(43)(44). The number of permutations was set at 1,000, and the nominal P value <0.01 was considered statistically significant.

Analysis of immune infiltration level
The analytical tool CIBERSORTx provides an estimation of member cell types' abundances in a mixed cell population, using gene expression data (31,32). We utilized CIBERSORTx to assess the proportions of 22 immune cell subtypes in the primary WT samples, using the RPKM data of protein-coding genes. The permutations for significance analysis were set at 1000. Then, we performed the 2-sided Wilcoxon test through the "ggpubr" package (version 0.2.5) in R, to compare differences in 22 immune cell subtypes between the low-and high-risk groups. Subsequently, through the "estimate" package (version 1.0.13), the immune and stromal scores were calculated based on the ESTIMATE algorithm (33), and Welch's t-test performed the comparisons of the stromal and immune scores between the 2 groups with "ggstatsplot" package (version 0.4.0) in R. A P value < 0.05 denoted statistical significance.

Statistical analysis
The R software version 3.6.2 and several R packages were used for statistical analyses, and a 2-tailed P value <0.05 denoted statistical significance. The univariate and multivariate Cox regression analyses were performed by the "survival" package. The LASSO Cox regression analyses were calculated by the "glmnet" package (35,36), and 10 times cross-validations were used to determine the best penalty parameter lambda. Kaplan-Meier analyses and survival curves were constructed through the "survival" package. The "rms" package generated the nomogram and the calibration curve. The time-dependent ROC curve analyses were performed by the "timeROC" package (37).

Identification and validation of the prognostic lncRNA signature
After filtering of low-expression lncRNAs, 2,927 lncRNAs were generated for further study. Subsequently, univariate Cox analysis was performed to obtain 244 lncRNAs associated with OS of WT patients (P<0.05). Among the 244 prognosis-associated lncRNAs, 24 were identified using the LASSO Cox method ( Figure 1). Then, through multivariate Cox analysis, ten lncRNAs (TENM3-AS1, AC022098.3, EMX2OS, AC099811.1, AL359710.1, ADAMTSL4-AS1, AC005944.1, MBNL1-AS1, AC002451.1, and AC120498.3) were selected from 24 candidates of lncRNAs as independent predictors associated with OS (P<0.05), and the coefficients of the 10 lncRNAs were obtained ( the risk score for each participant was calculated, and all participants were divided into high-risk (n=54) and low-risk (n=70) groups based on the optimal cutoff value (1.5069) for risk scores (Figure 2). Participants in the high-risk group had a significantly worse OS and EFS than those in the low-risk group (P<2E-16 and P=2.03E-04, respectively) ( Figure 3A,B, Table 3). Importantly, subgroup analyses based on the clinicopathological variables (i.e., gender, cooperative group protocol, stage, and histologic classification) indicated that the prognostic value of the risk scores for OS and EFS was independent of all these clinicopathological variables, except for EFS in the subgroup of stage I-II and favorable histology Wilms tumor (FHWT) (P=0.348 and P=0.326; Figure 3C,D,E,F,G,H,I,J,K,L,M,N,O,P, Table 3).

Development and assessment of a predictive nomogram
Based on the univariate Cox regression analysis results, the risk score of lncRNAs and 3 clinicopathological factors (i.e., gender, cooperative group protocol, and stage) were all significantly related to OS for WT patients ( Table 4). The subsequent multivariate Cox analysis (P value for Wald test: 5.11E-17) further demonstrated that the risk score remained a powerful and independent prognostic factor (P=9.22E-16) ( Table 4). Through integrating the 4 prognostic factors identified by univariate Cox analysis, a genomic-clinicopathologic nomogram was developed. This combined model is depicted in Figure 4A. The concordance index of the combined OS prediction model was 0.8780 (95% CI: 0.8407-0.9154). The nomograms' predictive accuracy was further evaluated by time-dependent ROC curves of 1-, 3-, and 5-year OS among the models. The results suggested that the combined model and the risk score model had a significantly greater AUC than that of the gender, cooperative group protocol, and stage models ( Figure 4B,C,D; Tables 5,6). Additionally, there was no significant difference in AUC between the combined model and the risk score model (1-year P=0.8610, 3-year P=0.7024, 5-year P=0.9354) ( Table 6). Furthermore, the calibration curve also demonstrated high consistency between the actual proportion of 500-day OS and the nomogram-predicted probability of 500-day OS ( Figure 4E). Meanwhile, DCA illustrated that the net benefit of the 500-day decision curve for the combined model was greater than that for the gender, cooperative group protocol, and stage model, further supporting the clinical value of the genomic-clinicopathologic nomogram ( Figure 4F).

Identification of differentially expressed genes and construction of PPI network
A total of 105 DEPCGs (24 upregulated and 81 downregulated) and 14 DELs (7 upregulated and 7 downregulated) were identified in the high-risk group ( Figure 5A,B). Furthermore, the PPI network of the 105 DEPCGs was generated ( Figure 5C). Following MCODE analysis, 2 densely connected regions in the PPI network were identified. The first region consisted of 15 target protein-coding genes, including MYL4, MYBPC2, ACTC1, MYLPF, DES, TRIM63, NEB, TTN, TNNI2,       Figure 6D), were significantly enriched in the low-risk group. Through GO gene sets, GSEA identified 11 GO terms consisted of 5 BP, 1 CC, and 5 MF, which were significantly enriched in the high-risk group (nominal P<0.01) ( Figure 6E), and 19 GO terms consisted of 16 BP, 2 CC, and 1 MF, which were significantly enriched in the lowrisk group (nominal P<0.01; Figure 6E).

Analysis of immune infiltration level
The CIBERSORTx analysis demonstrated that the proportion of monocytes activated dendritic cells and resting mast cells were significantly reduced in the high-risk group compared to that in the low-risk group (P=0.0464, P=0.0241, and P=0.00266), and the proportion of M1 macrophages (pro-inflammatory) in the low-risk group was substantially lower than that in the high-risk group (P=0.00944; Figure 7A). However, the stromal and immune scores showed no significance between groups (P=0.1339 and P=0. 8536; Figure 7B,C).

Discussion
In clinical practice, it is apparent that tumor progression and prognosis of patients are affected by multiple factors, yet many studies have only focused on single or few clinicopathologic factors for WT patients' prognosis (7)(8)(9). For instance, Tang et al. constructed nomograms to predict OS of children with WT based on 5 clinicopathologic factors, but AUCs of 3-and 5-year OS were not more than 0.74 (45). Nonetheless, prognostic tumor biomarkers also play a vital role in the diagnosis and treatment of tumors. For example, LOH for chromosomes 16q and 1p was demonstrated as a specific marker for increased relapse risk in favorable-histology WT (46). Several retrospective studies of heterogeneously-treated patients suggested an association between a gain of chromosomes 1q and tumor recurrence (47,48). Our previous study also indicated that astrocyte elevated gene-1 overexpression in histologically favorable WT was associated with poor patient prognosis (49).    (50). Very few studies currently focus on the development of prognostic models that integrate the clinicopathologic features and biomarkers based on comprehensive sequencing analysis for WT patient prognosis.
In this study, we performed multiple analyses to identify the prognostic lncRNA signature, including filtering of lowexpression lncRNAs, identification of lncRNAs associated with OS through univariate Cox analysis and LASSO method, and identification of independently predictive lncRNAs through multivariate Cox analysis. As a result, the potential minimum number of robust lncRNAs to predict WT patients' prognosis was obtained. Accordingly, a risk-score formula was constructed, and the subsequent log-rank tests and Kaplan-Meier analyses further confirmed the predictive value of the risk score for OS and EFS in the low-and high-risk groups of WT patients. Among the 10 prognostic lncRNAs in the present study, the low expression of EMX2OS (also known as EMX2-AS1 or NCRNA00045) was correlated with a poor prognosis in participants with WT [β=-0.55528, hazard ratio (HR) =0.573914; Table 2]. Interestingly, Gu et al. also reported that downregulation of EMX2OS was an independent prognosis factor for shorter recurrence-free survival of classical papillary thyroid cancer (51). Nomograms also called nomographs or alignment charts, which are widely used as prognostic devices in the field of medicine and oncology, can generate the individual probability of clinical events by integrating multiple prognostic variables (52). In this study, a genomicclinicopathologic predictive nomogram was developed through incorporating multiple prognostic factors including gender, cooperative group protocol, stage, and risk score to predict 1-, 3-, and 5-year OS of children with WT, which has shown stronger prediction with better AUCs (1-, 3-, and 5-year OS in the combined model of 0.9272, 0.9428, and 0.9259, respectively) (45,50). The combined model and the risk score model demonstrated significantly higher predictive accuracy, which was further verified by the time-dependent ROC curve and calibration curve analyses. Importantly, there was no significant difference in AUC between the combined and the risk score models, which indicated that the risk score played a critical role in the genomic-clinicopathologic nomogram. Hence, as a powerful and independent prognostic factor, the risk score achieved extremely high accuracy without combining clinicopathologic factors.
Nonetheless, several limitations need to be addressed in the present study. First, WT's sample size was relatively small,