Development and validation of novel prognostic models for zinc finger proteins-related genes in soft tissue sarcoma

Data processing

Raw data of 1555 zinc finger genes (Supplementary Table 1) were initially obtained from the HUGO Gene Nomenclature Committee database (https://www.genenames.org/). Differentially expressed genes (DEGs) were identified from expression profilings of STS and normal tissue samples obtained from GSE2719 using the “limma” package [16] of R 4.0.0 software. The following two criteria were adopted to select DEGs: adjusted p-values <0.05 and log2|fold change| values >1. The “ggplot2” package was utilized to create volcano plots while heatmaps of this study were constructed using the “pheatmap” package.

RNA-seq data and clinical information of patients diagnosed with STS were retrieved from TCGA (https://portal.gdc.cancer.gov/), GSE21050 [17], and GSE30929 [18]. And we subsequently downloaded data of genomic mutation in STS including both and copy number variation from TCGA database. Then these data regarding genomic mutation were analyzed using the R package ‘Rcircos’ to plot the copy number variation landscape that will enable us to better and more directly observe changes of ZNFs on human chromosomes.

Construction and validation of prognostic models

After excluding patients without sufficient follow-up time, we then randomly divided patients into training cohort or testing cohort. We adopted the survival R package to perform univariate Cox regression analysis to determine factors that were significantly associated with prognosis.

Then factors confirmed by the univariate Cox regression analysis to be significantly associated with prognosis were subsequently incorporated into the LASSO-penalized Cox regression analysis to establish prognostic models. The LASSO algorithm was adopted to select variables and for shrinkage using the “glmnet” R package. Additionally risk score of each individual patient in this study was calculated according to the following formula: Risk score = coef gene1 × Exp gene1 + coef gene2 × Exp gene2 + coef genei × Exp genei. Then each patient was assessed based on his or her median risk score calculated according to the aforementioned formula and divided into the high-risk group or low-risk group. Overall survival (OS) of patients in the high-risk group was compared with that of patients in the low-risk group by performing the log-rank test. And predictive capability of the prognostic model was evaluated by performing receiver operating characteristics (ROC) analysis using the “survivalROC” package. Testing cohort obtained from TCGA database and external cohort retrieved from GEO database (GSE21050) were used as the validation group. In a similar way, we established a predictive model for progression-free survival (PFS) as well. In this case, testing cohort obtained from TCGA database and external cohort retrieved from GEO database (GSE30929) were used as the validation group. Then prognostic significance of risk score and the aforementioned clinical variables was evaluated among patients in the TCGA training cohort and testing cohort by performing univariate and multivariate Cox regression analyses. Ultimately, we established nomograms using the “rms” package to help us more efficiently predict OS and PFS of patients with STS.

GO enrichment and gene set variation analysis (GSVA)

We conducted GO enrichment analysis utilizing the “clusterProfiler” package of R software. Through a hypergeometric distribution using P<0.05 as the significance threshold, we then performed functional enrichment analyses, results of which were made visualized using the “ggplot2” packages of R software.

Variations of key gene sets were estimated using the GSVA package of R software as an unsupervised, non-parametric method [19]. A gene expression matrix was used as the input for GSVA algorithm. GSVA scores were calculated in a non-parametrical way using the Kolmogorov-Smirnov (KS)-like random walk statistic and a negative value for a particular sample and gene set.

Consensus clustering

Consensus clustering was used to classify CAD cased into different subgroups [20]. Consensus clustering was accomplished using the K-means algorithm with the Spearman distance. Maximum number of clusters was set at 10. Final cluster number was identified by cluster consensus score and the consensus matrix.

Immune scores and tumor infiltrating immune cell (TIIC) analysis

The ESTIMATE algorithm was adopted to calculate immune and stromal scores of the samples from TCGA datasets. Based on the median immune score, patients were assigned into the high-score group or low-score group. OS of high-score group was compared with that of low-score group via log-rank test. TIICs in TCGA datasets were assessed using the CIBERSORT analytical tool (https://cibersort.stanford.edu/) [21]. 22 immune cells, abundance ratio matrix was identified at p values less than 0.05.

Cell culture and transfection

A673 and SW982 cell lines were purchased from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China). A673 and SW982 cells were cultured in DMEM (Biological Industries, Shanghai, China) or L-15 (Sigma) containing 10% fetal bovine serum (Biological Industries). ZNF141 siRNA and negative control siRNA oligonucleotides were designed and synthesized by Biosyntech (Suzhou, China). The sequences of si-1 and si-2 are shown in Supplementary Table 2. The siRNA transfections were performed using Lipofectamine RNAiMAX (Invitrogen, Carlsbad, CA, USA). Then we transfected ZNF141-overexpression-promoting plasmid (GV-ZNF141) and negative control plasmid (GV-Vector) that were designed and synthesized by GeneChem (Shanghai, China) into SW982 cells with Lipofectamine 3000 (Invitrogen, Carlsbad, CA, USA).

Real-time PCR and Western blotting

Following the manufacturer’s instructions, total RNA was retrieved with AG RNAex Pro Reagent (Accurate Biology, Changsha, China) and was reverse-transcribed into cDNA using the Evo M-MLV RT Premix Kit (Accurate Biology). Real-time PCR (RT-PCR) assays were performed using the SYBR Green Premix Pro Taq HS qPCR Kit (Accurate Biology) according to the manufacturer’s protocols. The primer sequences are listed in Supplementary Table 2.

Total protein of STS cells were extracted with RIPA lysis buffer (Beyotime Biotechnology, Shanghai, China), which were then separated with the method of SDS-PAGE. Then the proteins separated by SDS-PAGE were transferred onto PVDF membranes (Merck Millipore, Billerica, MA, USA). After being loaded with transferred proteins, these PVDF membranes were soaked in 5% BSA solution for at least one hour at room temperature, which was immediately followed by incubation with different primary antibodies at 4° C overnight. On the second day, the PVDF membranes were incubated with secondary antibodies for one hour at room temperature after washing three times. After washing with TBST solution, the PVDF membranes were visualized using a BeyoECL Plus Kit (Beyotime Biotechnology). Anti-GAPDH (10494-1-AP) were purchased from Proteintech (Wuhan, China) while anti-ZNF141 antibody (TA339080) was produced by OriGene Technologies (Wuxi, China).

Cell proliferation and plate clone formation

Cell proliferation was detected using the Cell Counting Kit-8 (CCK8, Beyotime Biotechnology) according to the manufacturer’s instructions. The cells were cultured in a 96-well plate (2000 cells/well) at 37° C for 0, 24, 48, and 72 h. For the plate clone formation assay, 1000 cells were seeded into 6-well plates and cultivated for two weeks. The number of colonies with >100 cells was counted under a light microscope, and the colonies were visualized.

Statistical analysis

R software was used for most bioinformatics and statistical analyses, including RNA-seq data normalization and transformation, DEG analysis, survival analyses, ROC analysis, CIBERSORT, ESTIMATE, GSVA, and Spearman rank correlation analysis. For the in vitro experiments, all quantitative data are presented as mean ± standard deviation of three independent experiments. Differences between the three groups were analyzed with one-way ANOVA using GraphPad Prism 8.0 (GraphPad, La Jolla, CA, USA). Statistical significance was set at p<0.05.

Data availability statement

All datasets presented in this study are included in the article/Supplementary Material.

Differentially expressed ZNFs in STS

The study design is illustrated in Figure 1. Differentially expressed ZNFs between normal tissues and tumor tissues were identified from GSE2719 using the “limma” package in R software. The heatmap demonstrating expression of ZNFs was presented in Figure 2A. A total of differentially expressed 110 ZNFs were identified with 82 ones up-regulated and 28 down-regulated (Figure 2B and Supplementary Table 3).

Construction and validation of prognostic models

Patients without sufficient follow-up time length were excluded from this study and remaining patients were randomly divided into the training cohort or the testing cohort. Then univariate Cox regression analysis was accomplished to determine variables significantly associated with OS of patients in the TCGA training cohort. A total of nine differentially expressed ZNFs were demonstrated to be significantly associated with OS (HLTF, ZNF292, ZNF141, LDB3, PHF14, ZNF322, PDLIM1, NR3C2, and LIMS2) (Figure 2C). LASSO regression analysis was accomplished using the “glmnet” package of R software to select these identified ZNFs. Then the predictive model was established using five-fold cross-validation and we demonstrated the confidence interval under each lambda in Figure 3A. Trajectory of the coefficient for each gene with a value of -in(lambda) is shown in Figure 3B. There was a total of nine genes identified as the model gene signature. The formula was as follows:

$$\begin{array}{l}\text{Risk score}= (0.225 \times \text{Exp} \text{HLTF})+( 0.495\\ \times \text{Exp ZNF}292)+( 0.971\\ \times \text{Exp} \text{ZNF}141)+(-0.136\\ \times \text{Exp LDB}3)+(0.717\\ \times \text{Exp PHF}14)+ (0.889\\ \times \text{Exp ZNF}322)+(-0.125\\ \times \text{Exp} \text{PDLIM}1)+(-0.809\\ \times \text{Exp} \text{NR}3\text{C}2)+ (-0.098\\ \times \text{Exp LIMS}2).\end{array}$$

Risk score of each patient was calculated and each patient from TCGA training cohort was divided into low-risk group or high-risk group. Then we performed survival analyses, results of which revealed that in comparison with patients in the low-risk group, those belonging to the high-risk group had significantly worse OS (p<0.001, Figure 3C). Additionally, time-dependent ROC analysis was performed, results of which revealed that areas under ROC curve (AUC) at one, three, and five years were 0.807, 0.792, and 0.799, respectively (Figure 3C). We then conducted model verification using data of the TCGA testing cohort (Figure 3D) and the GSE21050 cohort (Figure 3E) to further assess the validity and stability of the prognostic models, results of which revealed that compared with those in the low-risk group, patients of the high-risk group had significantly poorer prognosis. The AUC of the model in TCGA testing cohort at one, three, and five years was 0.789, 0.699, and 0.636, respectively (Figure 3D), and those of the GSE21050 cohort were 0.613, 0.607, and 0.577, respectively (Figure 3E).

In a similar way, we then established and validated a model predicting PFS of patients with STS. Firstly, datasets of STS from TCGA were randomly divided into the training cohort or the testing cohort. Univariate Cox regression analysis was accomplished to evaluate the impacts of expression levels of ZNFs on PFS of patients from the training cohort. A total of nine ZNFs including ZIC1, RNF128, ZNF141, ZHX2, TRIM9, ZNF281, NR3C2, ZNHIT2, and LIMS2 were demonstrated to be significantly associated with PFS (Figure 2D). Then a subsequent LASSO regression analysis was performed to screen ZNFs (Figure 4A, 4B), results of which revealed that seven genes were included in the predictive model for PFS. The following was the formula:

$$\begin{array}{l}Risk score= (0.077 \times Exp ZIC1) +( 0.195\\ \times Exp ZNF141)+( -0.273\\ \times Exp ZHX2)+( 0.041\\ \times Exp ZNF281)+ (-0.420\\ \times Exp ZNHIT2)+ (-0.398\\ \times Exp NR3C2)+ (-0.041\\ \times Exp LIMS2). \end{array}$$

For patients in the training cohort, patients with high risk had significantly shorter PFS than those with low risk (p<0.001, Figure 4C). AUC of the PFS model at one, three, and five years was 0.711, 0.725, and 0.716, respectively (Figure 4C). Then patients from TCGA testing cohort (Figure 4D) and the GSE30929 cohort (Figure 4E) were utilized to further validate this predictive model for PFS. Like in the training cohort, patients with high risk were demonstrated to have worse PFS than those with low risk. The AUC of the model in TCGA testing cohort was 0.681, 0.635, and 0.624 at one, three, and five years, (Figure 4D), and those in the GSE30929 cohort were 0.694, 0.679, and 0.668, respectively (Figure 4E).

Moreover, prognostic significance of risk score and the aforementioned clinical variables was evaluated among patients from TCGA training cohort and testing cohort by accomplishing Cox regression analysis, results of which revealed that for patients from both cohorts, risk score was an independent predictive factor for both OS (p<0.01, Supplementary Figure 1A, 1B) and PFS (p<0.05, Supplementary Figure 2A, 2B). Thus, through the aforementioned statistical analyses, we identified prognosis-associated ZNFs and established models that could reliably and efficiently predict long-term outcomes of patients diagnosed with STS.

Building predictive nomograms

A clinically useful nomogram that will help physicians more accurately predict OS (Figure 5A) or PFS (Figure 5B) of patients with STS was generated. According to the results of multivariate analysis among patients from the validation cohort, each variable was vested with a corresponding point based on the scale from the nomogram. A horizontal line was drawn to determine each variable^’s point. The total score of each patient was obtained by adding all the points together and based on the total score, we could more accurately estimate the one-, three-, and five-year survival rates of patients with STS.

GO enrichment analysis and GSVA analysis

In order to evaluate the possible functions of differentially expressed ZNFs in STS, we then performed GO enrichment analysis using the “clusterProfiler” package of R software. As shown in Figure 6A, we revealed that ZNFs were mainly located in PML bodies, stress fibers, contractile actin filament bundles, filamentous actin, and actin filament bundles. As for the specific molecular functions of these ZNFs, we found that ZNFs were mainly involved in DNA-binding transcriptional repressor and activator activity, ubiquitin-like protein transferase activity, ubiquitin-protein transferase activity and steroid hormone receptor activity.

In order to investigate the possible mechanisms through which ZNFs-based risk score model predicted prognosis of patients with STS, we then performed GSVA analysis. In the model predicting OS, the primary differences in the KEGG pathway between high-risk group and low-risk group were as follows: vascular smooth muscle contraction, phosphatidylinositol signaling system, inositol phosphate metabolism, adipocytokine signaling pathway and inositol phosphate metabolism (Figure 6B), which were significantly more activated among patients in the low-risk group and this might lead to better OS. In the model predicting PFS, the primary differences in the KEGG pathway between high-risk group and low-risk group were as follows: histidine, tyrosine, tryptophan, arachidonic acid, fatty acid, and drug metabolism cytochrome P450, and metabolism of xenobiotics by cytochrome P450 (Figure 6C), which were significantly more activated among patients in the low-risk group and this might lead to better PFS.

Molecular subtypes and immune infiltration of STS based on ZNFs

Thirteen genes (HLTF, LDB3, LIMS2, NR3C2, PDLIM1, PHF14, ZHX2, ZIC1, ZNF141, ZNF281, ZNF292, ZNF322, and ZNHIT2) were identified to be significantly associated with prognosis. Then we evaluated the prevalence of somatic mutations in 1,555 ZNFs in STS. Of the 237 samples, a total of 51 (21.52%) ones were demonstrated to be with genetic alterations, primarily missense mutations, nonsense mutations, and multiple hits in ZNFs. PCLO was the frequent mutation (Figure 7A). Then through further analysis of the 13 prognosis-related ZNFs, we demonstrated that copy number variation (CNV) mutations were prevalent. Specifically, PHF14, ZHX2, and ZNF141 showed widespread amplifications, while LDB3, PDLIM1, ZNF292, ZNF322, and ZNHIT2 showed prevalent deletions (Figure 7B). The locations of CNV alterations in the 13 prognosis-related ZNFs on the chromosomes are shown in Figure 7C.

Consensus ClusterPlus package was adopted to divide STS samples from TCGA into k (k=2-9) subtypes based on 13 ZNFs significantly associated with prognosis. Best division of the CDF curve based on consensus scores could be achieved when k was 4. Accordingly, we identified four distinct subtypes (Figure 8A), including 34, 67, 87, and 65 cases in cluster A, B, C, and D, respectively (Figure 8B). It was demonstrated through survival analysis that cluster B had prominent survival advantages over other clusters and cluster D was found to have the worst OS (p<0.001, Figure 8B) and PFS (p=0.001, Figure 8C). The heatmap of the 13 prognosis-related ZNFs for four different subtypes was presented in Figure 8D.

Subsequently, we further evaluated the immunological characteristics of the four subtypes. The immune and stromal scores of samples in the TCGA cohort were quantified using the ESTIMATE algorithm, results of which demonstrated that in comparison with patients with high scores, OS (p<0.01, Supplementary Figure 3A) and PFS (p<0.05, Supplementary Figure 3B) of those with low scores were significantly worse. Four subtypes were revealed to have significantly different immune scores. Immune scores of subtypes B and C were significantly higher than those of subtypes A and D (Figure 8E). Results of analysis of the immunological characteristics implied that unlike subtypes B and C, immune responses within subtypes A and D STS were significantly suppressed. Furthermore, we assessed the expression of several immune checkpoints in four different subtypes, results of which demonstrated that compared with those in subtypes A and D, expression levels of most immune checkpoints were significantly higher in subtypes B and C (Figure 8F), which was consistent with results of immune score analysis. Then the immune infiltration profile of the TCGA cohort was characterized with CIBERSORT. For STS samples from TCGA, the major constituents of TIICs were CD8 T cells, resting CD4 memory T cells, resting mast cells, M2 and M0 macrophages (Figure 8G). Additionally, it was also observed that infiltration proportions of five immune cells in four subtypes were significantly different (Figure 8G). Significantly much more resting CD4 memory T cells and mast cells were observed in subtype A STS than in other three subtypes. Moreover, it was demonstrated that CD8 T cells within subtype B STS were significantly more prominent while significantly more abundant infiltration of M2 within subtypes C and D STS and M0 within subtype D STS were detected. Therefore, given the fact that the dominant immune cells within subtype B were CD8+ T cells (the main killers of tumor cells, Philip and Schietinger, 2021) while the dominant immune cells within subtype C were M2 macrophages (the main tumor promoters, Najafi et al., 2019), we could explain why prognosis of subtype B STS was significantly better than that of subtype C STS. Similarly, the most prominent immune cells within subtype A STS were resting CD4 memory T cells and mast cells and the most prominent immune cells within subtype D STS were M0 and M2 macrophages, which might explain why patients with subtype A STS had better prognosis than those with subtype D STS (Supplementary Figure 3C). Therefore, considering all the results, we could draw the conclusion that proportion of immune cells within STS would significantly affect patients' prognosis.

ZNF141 promoted the proliferation and viability of STS cells

ZNF141 has been identified as an OS- and RFS-related ZNF, but its role in STS is still unclear. Therefore, we verified the effect of ZNF141 on STS cells in vitro. ZNF141 mRNA expression was upregulated in tumor tissues compared to that in normal tissues (Figure 9A). The results from TCGA dataset show that patients with high ZNF141 mRNA levels had poorer OS (Figure 9B) and RFS (Figure 9C) than those with low ZNF141 mRNA levels. Subsequently, siRNAs (si-nc, si-1, and si-2) were transfected to A673 cells. RT-PCR and western blotting assays showed that ZNF141 expression was significantly downregulated in A673 cells transfected with si-1 or si-2 compared with those transfected with si-nc (Figure 9D). The results of the CCK8 showed that the proliferation of A673 cells with knockdown ZNF141 (si-1 or si-2) was weaker than that of the control group (si-nc) (Figure 9E). Plate clone formation assays showed that the colony numbers of A673 cells with ZNF141 knockdown were significantly reduced (Figure 9F). Next, ZNF141 overexpression (GV- ZNF141) or control plasmid (GV-Vector) were transfected to SW982 cells. The successful overexpression of ZNF141 in SW982 cells was verified by both RT-PCR and western blotting assays (Figure 9G). CCK8 assays showed that the proliferation of SW982 cells with ZNF141 overexpression (OE) was stronger than that of the control group (NC) (Figure 9H). Plate clone formation assays showed that the colony numbers of SW982 cells with ZNF141 overexpression were significantly increased (Figure 9I). These results suggested that ZNF141 promoted the proliferation and viability of STS cells, and we could conclude that ZNF141 played critical roles in promoting progression of STS.