The RNA transcription data and corresponding clinical data of the sarcoma sample were downloaded from the TCGA-SARC dataset. After excluding samples with incomplete RNA transcription data, the dataset contained 259 primary tumors and 2 normal tissue samples from 259 sarcoma cases. The pathological classification of the sarcoma samples included leiomyosarcomas (n = 104), dedifferentiated liposarcomas (n = 59), undifferentiated sarcomas (n = 34), fibromyxosarcomas (n = 25), malignant fibrous histiocytomas (n = 12), synovial sarcomas (n = 10), malignant peripheral nerve sheath tumor (n = 9), and other types (n = 6). Meanwhile, the RNA transcription data of 515 normal adipose tissue were downloaded from the GTEx database. Furthermore, the transcribed data and correlated information of the two GEO datasets, GSE68591 and GSE40021, were downloaded and analyzed by the R package GEOquery. The GSE68591 dataset included mRNA transcription data of 68 human sarcoma cell lines and 5 normal human cells. The GSE40021 dataset included mRNA transcription data from 53 primary synovial sarcoma samples after excluding samples with no follow-up period.

Oncomine (www.oncomine.org) is the largest comprehensive bioinformatics analysis platform comprising 86,733 tumor and non-tumor samples from 715 databases (Pan et al., 2019; Zou et al., 2021). Therefore, we analyzed the expression of GINS subunits in sarcoma and non-tumor samples via the Oncomine database. The criteria for filtering datasets included a p-value of 0.01, a fold change of 2, and genes ranked in the top 10% (Sun et al., 2019).

TIMER is a bioinformatics website that comprehensively evaluates RNA expression and multiple immune cell infiltrations in 10,897 samples of 32 cancers from the TCGA database (Li et al., 2017). In this study, the TIMER database was used to verify the significant difference in GINS subunits expression level between tumor and non-tumor tissues. Meanwhile, the TIMER database was used to investigate the relationship between the transcription levels of GINS genes and immune cell infiltration in SARC (Chen et al., 2020).

Kaplan-Meier plotter is an online tool that analyzes the survival effects of 54,675 genes from 10,461 cancer samples across 21 cancer types (Lanczky and Gyorffy, 2021). The Log-rank p-value, HRs, 95% CI, and survival curve were calculated and displayed on the web page (Li et al., 2020b). Moreover, we analyzed the effect of the mRNA transcription level of GINS genes on survival in SARC by Kaplan-Meier Plotter.

GEPIA is a bioinformatics tool that comprehensively analyzes RNA expression sequencing from 9,736 tumors and 8,587 non-tumor samples, from TCGA and GTEx databases (Tang et al., 2017). The GEPIA database was used for single-gene survival analysis to assess the prognostic value of GINS.

The cBioPortal for Cancer Genomics is an online bioinformatics site for multi-dimensional, visual analysis of cancer genomes from multiple cancer genome databases (Gao et al., 2013). Through the cBioPortal website, we obtained the mutation profiles, mRNA expression levels, and neighbor genes of GINS1/2/3/4 in sarcoma. The threshold value of mRNA expression z scores was set as ± 2.0 (Zhang et al., 2020a).

CCLE is an online bioinformatics database containing the genetic characteristics of 1072 cell lines for various cancers (Ghandi et al., 2019). This study analyzed the mRNA transcription of GINS genes in different cancer cell lines using CCLE data to further understand GINS genes expression in sarcoma.

The University of Alabama at Birmingham Cancer Data Analysis Portal (UALCAN) is a bioinformatics website for comprehensive analyzing RNA-seq data and clinical information for 31 cancer types in TCGA (Chandrashekar et al., 2017). To identify the co-expressed genes of GINS1/2/3/4 in SARC, we retrieved the co-expressed genes of GINS1/2/3/4 through UALCAN, and then intersected these results.

The enrichment analysis was performed, and clusters were visualized by the R package software of clusterProfiler (v 4.2.0) (Yu et al., 2012). The present study conducted Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of GINS subunits and their related genes in RStudio. GO enrichment analysis predicted gene function from biological process (BP), cellular component (CC), and molecular function (MF). A p-value < 0.05 was statistically significant (Gong et al., 2021).

The starBase and mirDIP websites were used to predict miRNAs targeting GINS1/2/3/4 genes, and lncRNAs binding the target miRNA were predicted by starBase (Li et al., 2014; Tokar et al., 2018). The expression correlations between target miRNA and GINS gene, as well as lncRNA and target miRNA, were verified by co-expression analysis using the R package limma (v 3.50.0) based on the TCGA-SARC dataset. Our co-expression analysis set the threshold value of the p-value as < 0.005.

The R (v 4.1.0) and online analysis tools were used to perform all statistical analyses. We used R package ggpubr (v 0.4.0), ggplot2 (v 3.3.5), corrplot (v 0.92), clusterProfiler (v 4.2.0) and limma (v 3.50.0) to visualize data. The Wilcoxon rank-sum test was conducted to analyze the differential expression of GINS genes between sarcoma and non-tumor samples. The log-rank test was performed to compare survival time between high and low gene expression groups of SARC patients. The Spearman test was applied for correlation expression analysis between genes. The threshold of statistically significant was set as p-value < 0.05.

It is well-established that the GINS complex is composed of GINS1/2/3/4 in mammalian cells. First, the Oncomine database was used to assess different mRNA levels of the GINS subunit in SARC patients. We found high GINS1/2/3/4 mRNA expression levels in sarcoma compared to non-tumor tissue (Figure 2). Then, mRNA levels of GINS1/2/3/4 in SARC and non-tumor tissues were detected in TCGA and GTEx databases (Figures 3A–E). The transcriptional levels of GINS1/2/3/4 were significantly higher in SARC than in non-tumor tissues. Meanwhile, the differential expression of GINS genes was verified in the GEO dataset (Figure 3F). From the GSE68591 dataset, SARC samples exhibited significantly high transcriptional levels of GINS1/2/3/4 compared to non-tumor samples. Moreover, the subcellular localization of GINS subunits was explored in the HPA database (Figure 3G). The immunofluorescence imaging results exhibited GINS1/2/3/4 proteins (green) mainly localized to nucleoplasm (blue) in the osteosarcoma cell line (U2OS). Subsequently, the cBioPortal database was used to assess the frequency of gene changes of the GINS subunits in sarcomas. According to these data, GINS1/2/3/4 were altered in 8, 7, 3, and 12 percent of the SARC samples (Figure 3H). Enhanced mRNA transcription was the most frequent GINS subunits alteration in sarcomas (Figure 3I).

The CCLE database analysis showed that GINS1/2/3/4 were relatively highly expressed in sarcoma compared to most cancer cell lines. Meanwhile, it was found that at the RNAseq level, GINS1/2/3/4 expression in sarcoma ranked eighth, fourth, sixth and third among various cancer cell lines (Figures 4A–D).

We used the Kaplan–Meier Plotter Database, GEPIA database, and the R package survival to determine the effect of the GINS gene on the survival of SARC patients. Kaplan–Meier Plotter Database analysis showed that higher expression levels of GINS1/2/3/4 mRNA correlated with lower OS in sarcoma (Figures 5A–D). SARC patients with high GINS1 expression exhibited shorter median OS than those with low GINS1 expression (51.2 vs. 85.83 months; HR = 1.89, p = 0.0017). High GINS2 expression in SARC patients was significantly poorer median OS in contrast with low GINS2 expression (49.27 vs. 82.13 months; HR = 1.81, p = 0.0031). High GINS3 expression was associated with a shorter upper quartile OS than low GINS3 expression (22.93 vs. 38.07 months; HR = 1.74, p = 0.011). Moreover, high GINS4 expression in SARC patients was associated with a poor upper quartile OS in contrast with low GINS4 expressed patients (19.97 vs. 54.23 months; HR = 2.88, p = 2.6e-05). In terms of Recurrence-free survival (RFS), high expression of GINS1/2/3/4 predicted a poorer RFS in sarcoma (Figures 5E–H). High GINS1 expression in SARC patients exhibited poorer median RFS compared to low GINS1 expression (18.6 vs. 88.63 months; HR = 2.61, p = 9e-05). Moreover, high GINS2 expression in SARC patients correlated with a significantly shorter median RFS in contrast with low GINS2 expression (17.2 vs. 88.63 months; HR = 2.79, p = 2.4e-05). High GINS3 expression was associated with a poorer upper quartile RFS compared to low GINS3 expression (11.63 vs. 21.03 months; HR = 2, p = 0.023). High GINS4 expression in SARC patients correlated with a shorter upper quartile RFS than low GINS4 expression (10.03 vs. 17.77 months; HR = 2.24, p = 0.0018).

Similar results were obtained during the GEPIA database. We found that increased mRNA expression of GINS1/3 in sarcoma was closely related to poorer OS (Figures 6A,C). In this respect, high GINS1 expression in SARC patients correlated with shorter OS than low GINS1 expression (HR = 1.8, p = 0.0062). SARC patients with high GINS3 expression had a poorer OS than those with low GINS3 (HR = 1.5, p = 0.0035). GINS2 and GINS4 expression in sarcoma also affected the OS, although there was no statistical significance (Figures 6B,D). In terms of Disease-free survival (DFS), high expression of GINS1/2/3 also predicted shorter DFS in sarcoma (Figures 6E–G). We found that high GINS1 expression was associated with poorer DFS than low GINS1 expression (HR = 1.7, p = 0.004). The DFS of SARC patients with high GINS2 expression was shorter than those with low GINS2 expression (HR = 1.9, p = 0.00026). High GINS3 expression correlated with a poorer DFS compared with low GINS3 expression (HR = 1.4, p = 0.045). GINS4 expression was also associated with a relatively poor DFS in sarcoma; however, there was no statistical significance (Figure 6H).

Finally, the prognostic value of GINS genes was verified in the GEO dataset. ROC analysis of the GSE40021 dataset showed that upregulated GINS1/2/3 expression was closely related to sarcoma metastasis (Figures 7A–C). The AUC values of GINS1 for predicting survival at 1, 3, and 5 years were 0.77, 0.72, and 0.7, while the corresponding AUC values for GINS2 were 0.89, 0.78, and 0.82, respectively. ROC analysis also demonstrated a good predictive value of GINS3 for survival at 1, 3, and 5 years with AUC values of 0.85, 0.74, and 0.76, respectively. The corresponding AUC values for GINS4 were 0.68, 0.55, and 0.65, respectively (Figure 7D). The above results indicate that GINS1/2/3 have prognostic significance in predicting sarcoma metastasis. Furthermore, Kaplan-Meier analysis of metastasis-free survival (MFS) of GINS subunits was performed in the GSE40021 dataset (Figures 7E–H). SARC patients with high GINS1 expression exhibited shorter mean MFS than patients with low GINS1 expression (28.95 vs. 104.17 months; p = 2.00e-05). The mean MFS of SARC patients with high GINS2 expression was poorer than those with low GINS2 expression (42.85 vs. 116.59 months; p = 3.03e-05). Moreover, patients with high GINS3 expression had a shorter mean MFS than those with low GINS3 expression (48.71 vs. 102.85 months; p = 0.00025). Finally, SARC patients with high GINS4 expression had a poor mean MFS compared with patients with low GINS4 expression (34.49 vs. 97.30 months; p = 0.0018).

Our comprehensive analysis identifed 675 genes and the top 50 most frequently genes co-expressed with GINS1/2/3/4 in sarcoma by UALCAN and cBioPortal database, respectively (Figure 8A). Finally, 45 co-expressed genes were simultaneously predicted by these two databases (Figure 8B). These data suggest that EXO1, UHRF1, DLGAP5, MCM10, ESPL1, TPX2, NCAPH, CENPI, ERCC6L, CCNA2, POLQ, TTK, AURKA, BUB1, KIF18B, CDCA8, KIF4A, KIF4B, NCAPG, CENPA, CENPE, CCNB1, KIF2C, FOXM1, CCNB2, HJURP, KIF18A, BUB1B, NEK2, DEPDC1, CENPF, PLK1, TROAP, KIFC1, ASPM, CDC20, TOP2A, NUF2, SKA1, FAM72D, MELK, UBE2C, AURKB, RAD54L, and GTSE1 were closely related to the biological functions and signaling pathways of differentially expressed GINS subunits in SARC. The protein-protein interaction results showed that GINS1/2/3/4 subunits and 45 co-expressed genes were closely related (Figure 8C). In addition, the cBioPortal and GEPIA databases were used to analyze the co-expression relationships between different GINS genes. Analysis in the cBioPortal database showed GINS1 was positively correlated with GINS2/3/4, with Pearson correlation coefficient (r) values of 0.56, 0.46, and 0.57, respectively (P_s < 0.05). GINS2 was significantly correlated with GINS3 and GINS4, with r values of 0.47 and 0.39, respectively (P_s < 0.05). GINS3 was significantly correlated with GINS4 with a correlation coefficient value of 0.29 (p < 0.05) (Figure 8D). Based on GEPIA database analysis, similar results were found. GINS1 was positively associated with GINS2/3/4, with r values of 0.67, 0.6, and 0.64, respectively (P_s < 0.05) (Figures 8E–G). GINS2 was significantly associated with GINS3 and GINS4, with r values of 0.55 and 0.47, respectively (P_s < 0.05) (Figures 8H,I). GINS3 was correlated with GINS4 with a correlation coefficient value of 0.4 (p < 0.05) (Figure 8J).

The GINS genes and their neighboring genes in SARC were analyzed by R package clusterProfiler for GO and KEGG enrichment analysis. These genes were found primarily involved in chromosome segregation and cell division, associated with the genesis and progression of malignant tumors. Consistent with the literature, in the BP category, these genes were significantly enriched in the mitosis, spindle, and microtubule cytoskeleton (Figure 9A). The spindle, chromosome centromeric region, chromosomal region, condensed chromosome, condensed chromosome centromeric region, kinetochore, microtubule, mitotic spindle, spindle microtubule, and spindle midzone were the top ten most significantly enriched items in the CC category (Figure 9B). In terms of MF, the GO terms were significantly enriched in processes, including microtubule binding, protein serine/threonine/tyrosine kinase activity, microtubule motor activity, cytoskeletal motor activity, ATP hydrolysis activity, and cyclin−dependent protein serine/threonine kinase regulator activity (Figure 9C). The KEGG pathway analysis showed that GINS genes and their neighboring genes were remarkably enriched in cell cycle, oocyte meiosis, progesterone-mediated oocyte maturation, human T-cell leukemia virus one infection, and cellular senescence (Figure 9D). More detailed information is shown in Table 1.

The GINS gene-related ceRNA networks were predicted and constructed using online tools and Rstudio based on TCGA-SARC dataset. Through the starBase and mirDIP websites, 230 and 1216 miRNAs were predicted to target GINS1, respectively, of which 192 miRNAs were simultaneously predicted by both websites (Figure 10A). Meanwhile, the co-expression analysis of these 192 miRNAs in TCGA-SARC dataset illustrate that hsa-miR-10a-5p, hsa-miR-26b-5p, hsa-miR-29a-3p, hsa-miR-34a-5p, hsa-miR-139-5p and hsa-miR-150-5p were negatively co-expressed with GINS1, with R values of −0.22, −0.18, −0.3, −0.21, −0.18 and −0.26, respectively (P_s < 0.005) (Figures 10E–J). Moreover, the starBase and mirDIP websites predicted the binding of 72 and 686 miRNAs to GINS2, respectively, of which 51 were predicted by both websites (Figure 10B). The co-expression analysis of these 51 miRNAs in TCGA-SARC dataset showed that only hsa-miR-29a-3p was negatively correlated to the transcriptional levels of GINS2, with the r values of −0.21 (p = 0.00057) (Figure 10K). In addition, we predicted 94 and 815 miRNAs targeting GINS3 using starBase and mirDIP websites, respectively, of which 77 miRNAs were predicted by both websites (Figure 10C). Then, a correlation analysis of these 77 miRNAs expression in TCGA-SARC dataset showed that only hsa-miR-145-5p was negatively co-expressed with GINS3, with the r values of −0.2 (p = 0.0013) (Figure 10L). Furthermore, starBase and mirDIP websites predicted 92 and 1587 miRNAs were bound to GINS4, respectively, of which 83 miRNAs were predicted by both websites (Figure 10D). The co-expression analysis of these 83 miRNAs in TCGA-SARC dataset showed that has-miR-26a-5p, has-miR-26b-5p, has-miR-127-3p, and has-miR-150-5p were negatively correlated with the transcriptional levels of GINS4, with the r values −0.31, −0.19, −0.19 and −0.2, respectively (P_s < 0.005) (Figures 10M–P). Subsequently, starBase predicted 77, 92, 92, 114, 125, 28, 75, 113 and 187 lncRNAs were bound to has-miR-10a-5p, has-miR-26a-5p, has-miR-26b-5p, has-miR-29a-3p, has-miR-34a-5p, has-miR-127-3p, has-miR-139-5p, has-miR-145-5p, and has-miR-150-5p, respectively. Co-expression analysis in TCGA-SARC dataset found that AC110769.2, AL356299.2, and AL121832.3 were found to be negatively co-expressed with has-miR-10a-5p (r = −0.22, −0.22 and −0.32) and positively correlated with the transcriptional levels of GINS1 (r = 0.3, 0.36 and 0.48). Moreover, AC016026.1, THUMPD3-AS1, LINC00205, ENTPD1-AS1, PSMD6-AS2, AC000120.1, NNT-AS1, MALAT1, AC105339.2, TUG1, AC023355.1, AL035425.3 and AP000974.1, were found to be negatively co-expressed with hsa-miR-26a-5p (r = −0.3, −0.39, −0.23, −0.21, −0.21, −0.22, −0.37, −0.32, −0.25, −0.34, −0.26, −0.22, and −0.33) and positively correlated with the transcriptional levels of GINS4 (r = 0.24, 0.31, 0.35, 0.26, 0.28, 0.26, 0.29, 0.23, 0.39, 0.39, 0.25, 0.21, and 0.3). Furthermore, NNT-AS1 and AP000974.1 were found to be negatively co-expressed with hsa-miR-26b-5p (r = −0.31 and −0.25) and positively correlated with the transcriptional levels of GINS1 (r = 0.3 and 0.28) and GINS4 (r = 0.29 and 0.3). In addition, THUMPD3-AS1, CRNDE, and MIR762HG were found to be negatively co-expressed with hsa-miR-29a-3p (r = −0.3, −0.28 and −0.26) and positively correlated with the transcriptional levels of GINS1 (r = 0.36, 0.28 and 0.3) and GINS2 (r = 0.23, 0.31 and 0.31). Subsequently, we found that LINC01521, SLC9A3-AS1, AC253536.3, LINC00665, AC073529.1, TMEM147-AS1, AC099811.1, CKMT2-AS1, AC120114.1, AC139887.2, TUG1, AC008147.2, AC093249.6, TERC, AC108704.2, AC104447.1 and SLFNL1-AS1 were negatively co-expressed with hsa-miR-34a-5p (r = −0.32, −0.27, −0.37, −0.25, −0.31, −0.24, −0.41, −0.41, −0.45, −0.27, −0.38, −0.23, −0.27, −0.21, −0.34, −0.26, and −0.45) and positively correlated with the transcriptional levels of GINS1 (r = 0.36, 0.22, 0.22, 0.34, 0.56, 0.26, 0.21, 0.22, 0.48, 0.22, 0.33, 0.23, 0.25, 0.21, 0.21, 0.36, and 0.32). Additionally, AP000253.1 and NFYC-AS1 were negatively co-expressed with hsa-miR-127-3p (r = −0.28 and −0.46) and positively correlated with the transcriptional levels of GINS4 (r = 0.23 and 0.3); AC125611.4, AC108704.2 and AC026401.3 were found to be negatively co-expressed with hsa-miR-139-5p (r = −0.29, −0.21 and −0.26) and positively correlated with the transcriptional levels of GINS1 (r = 0.4, 0.21 and 0.39). LINC00852, ATP2B1-AS1 and HELLPAR were found to be negatively co-expressed with hsa-miR-145-5p (r = −0.28, −0.39 and −0.24) and positively correlated with the transcriptional levels of GINS3 (r = 0.29, 0.24 and 0.21). Finally, ENTPD1-AS1, AC021078.1, FOXD2-AS1, PTPRG-AS1, OIP5-AS1, NNT-AS1, OTUD6B-AS1, AC125611.4, AC023355.1, AL035425.3, AC011447.3, AC022150.4, AC073957.3, and AC007191.1 were found to be negatively co-expressed with hsa-miR-150-5p (R = −0.24, −0.21, −0.23, −0.40, −0.32, −0.28, −0.24, −0.39, −0.22, −0.41, −0.31, −0.26, −0.25, and −0.37) and positively correlated with the transcriptional levels of GINS1 (r = 0.27, 0.21, 0.37, 0.35, 0.24, 0.30, 0.29, 0.40, 0.30, 0.27, 0.41, 0.42, 0.24, and 0.31) and GINS4 (r = 0.26, 0.27, 0.25, 0.25, 0.31, 0.20, 0.29, 0.25, 0.31, 0.25, 0.21, 0.36, 0.52, 0.35, and 0.38) (P_s < 0.005). The interactive relationship of the ceRNA network is shown in a Sankey diagram (Figure 11). Finally, THUMPD3-AS1 was found to be a lncRNA that regulates most mRNAs of GINS genes (GINS1/2/4) in sarcoma by binding to hsa-miR-26a-5p and hsa-miR-29a-3p simultaneously (Figures 12A,B). A sub-network of lncRNA THUMPD3-AS1, miRNA and GINS subunits in SARC was constucted by Cytoscape software (Figure 12C).

The TIMER database was utilized to detect the expression of GINS genes with tumor purity and immune cell infiltration and the effect of immune cell infiltration on survival in SARC. We found that GINS1 expression was positive correlated with tumor purity (r = 0.338, p = 5.94e-08), and negative correlated with the levels of tumor-infiltrating CD4+ T cell (r = −0.234, p = 2.56e-04) and macrophages (r = −0.233, p = 3.21e-04) (Figure 13A). GINS2 was positive correlated with tumor purity (r = 0.26, p = 3.78e-05), and negative correlated with the levels of tumor-infiltrating macrophages (r = −0.137, p = 3.57e-02) (Figure 13B). GINS3 was positive correlated with tumor purity (r = 0.102, p = 1.11e-01) (Figure 13C). Finally, GINS4 was positive correlated with tumor purity (r = 0.148, p = 2.06e-02), and negative correlated with the levels of tumor-infiltrating CD4+ T cell (r = −0.175, p = 6.57e-03) (Figure 13D). Meanwhile, high CD4+ T cell and neutrophil infiltration levels were associated with better survival outcomes in SARC (p < 0.05) (Figure 13E). These findings suggest that GINS genes may influence patient survival by interacting with the infiltrated immune cells in sarcomas.