This was an observational study using transcriptomic data of glioma in TCGA and CGGA. Based on our previous study of the function of SDC1 in glioma, our goal is to further explore the function of SDC1 in glioma. To study whether the expression of SDC1 in glioma is related to its immune infiltration, the samples in the data sets of glioma were evaluated through the comprehensive immune landscape of the TME through ESTIMATE algorithms. To ensure the reliability of the calculated immune landscape, the consistency of relevant clinical features with the immune landscape was evaluated. To find the specific immune cells that SDC1 may affect, the immune cell composition in the samples was also calculated by ssGSEA algorithms, the correlation of SDC1 between different immune cells was estimated. To find the genes that are not only related to SDC1, but also involved in glioma immune invasion, WGCNA was used to further data mining of differential genes of SDC1. After obtaining hub genes, GO and KEGG analysis was done to find the specific signal pathway. Finally, the reliability of the hub genes is evaluated by clinical prognosis analysis, and PPI predicts the possible mechanism of SDC1. Figure 1 illustrates the workflow chart of data preparation, processing, analysis, and validation.

The raw gene expression levels and clinical information of glioma were downloaded from the Chinese Glioma Genome Atlas (CGGA) (http://www.cgga.org.cn, downloaded on 07/21/2021), which includes two data sets, mRNAseq_693 (n = 693) and mRNAseq_325 (n = 693). TCGA RNAseq and clinical information of glioma are downloaded from the TCGA database (https://portal.gdc.cancer.gov, downloaded on 07/21/2021, n = 702 after removing six duplicate samples). The leucocyte gene signature matrix refer to Charoentong et al., 2017, download from https://www.cell.com/cms/10.1016/j.celrep.2016.12.019/attachment/f353dac9-4bf5-4a52-bb9a-775e74d5e968/mmc3.xlsx).

To visualize immune-related cell heterogeneity in the TME, two machine learning methods of gene expression profiles were used: the Estimation of STromal and Immune cells in Malignant Tumors using Expression data (ESTIMATE) (Yoshihara et al., 2013) and ssGSEA algorithms (Charoentong et al., 2017). First, the ESTIMATE algorithm was used to estimate the proportion of immune and stromal components in glioma samples. Spearman correlation analysis is used to assess the correlation between SDC1 expression and the degree of immune infiltration in glioma specimens. Second, to estimate the proportion of tumor-infiltrating immune cells (TICs), gene set variation analysis (GSVA) was used to estimate the variation of the immune cells through reference gene expression data set enrichment, that is, the leucocyte gene signature matrix (Charoentong et al., 2017); the proportions of 28 tumor-infiltrating lymphocytes were assessed in glioma tissues.

Based on the discovery that the expression of SDC1 increases with the increase of glioma grade, the tumor specimens with high and low expression of SDC1 were selected for differential analysis to study the possible function of SDC1 in glioma. After normalization of TCGA RNA-seq with the linear models for microarray data (LIMMA) package in the R language environment and log2 transformation, then the TCGA RNAseq, mRNAseq_693, and mRNAseq_325 were analyzed for the difference. The DEGs between the high and the low expression of SDC1 samples were identified (DEGs, p value < .05). The clustering heat map of differential genes is to assess the homogeneity of the samples. The Volcano maps are used to show genes that are significantly differentially expressed. Identification of unique genes within the DEG data sets in the three data sets, and the Venn diagram shows the overlapping results.

The WGCNA package is implemented in R (http://www.r-project.org/). WGCNA is an analysis method that analyzes gene expression patterns of multiple samples and is widely used in genome research. First, it is assumed that the gene network obeys a scale-free network, and the gene co-expression correlation matrix and the adjacency function formed by the gene network are defined, and then, the dissimilarity coefficients of different nodes are calculated, and a hierarchical clustering tree is constructed accordingly, which clusters genes with similar expression patterns. Different branches of the clustering tree represent different gene modules. The degree of gene co-expression in the module is high, and the degree of co-expression of genes belonging to different modules is low. Finally, the relationship between modules and specific traits or phenotypes was analyzed. Given the above algorithm, the expression of SDC1in glioma samples was selected as a clinical phenotype in our data mining. After removing outlier samples and those with incomplete clinical information, the adjacency matrix was constructed with 145 GBM samples. The soft thresholding power β was set to eight to ensure a correlation coefficient close to .9. A topological overlap matrix (TOM) plot was assessed as independent of each module. Module membership (MM) refers to the correlation coefficient between genes and module eigengenes and is used to describe the reliability of a gene belonging to a module. Finally, we calculated the correlation between the modules and the clinical phenotype, the modules that are not only related to SDC1, but also tumor immune infiltration, are selected to be the hub module.

The intro module connectivity reflects the importance of nodes; it represents the sum of the degree of correlation between this gene and other genes in that module. For each expression profile, gene significance (GS) was calculated to assess the correlation between each gene and each trait; MM represents the correlation between a single gene and its module, indicating whether it is consistent with the trend of the module. The top 50% of genes with the highest intro module connectivity, GS, and MM were selected as hub genes in the hub modules.

To understand the biological meaning of the hub genes, the cluster profile (R package) was used for performing GO enrichment and KEGG pathway analyses (p-value <.05 was regarded as significant). The genes related to immune cells and immune signaling pathways are selected as target genes.

To verify the importance of the expression of the key gene, the survival package was loaded in R and survminer was used for the survival analysis; 696 glioma tumor cases in TCGA had a detailed clinicopathological grading and survival time record, which were used for clinical correlation analysis. We evaluated the correlation between the target genes and clinical grade in the data set of glioma in the TCGA. The Kaplan–Meier method was used to plot the survival curve, and log-rank as the statistical significance test; p < .05 was considered significant.

GeneMANIA (http://www.genemania.org) is a flexible, user-friendly web interface, which is used to generate assumptions about gene function, analyze gene lists, and determine gene priority for functional analysis (Szklarczyk et al., 2011). STRING is a database of known and predicted interactions, including direct (physical) and indirect (functional) associations from various sources, including high-throughput experiments and genomic contexts (Warde-Farley et al., 2010). SDC1 and Hub Genes were submitted to the above two websites to analyzes the PPI network.

To explore the correlation of SDC1 expression with immune infiltration, we first used the ESTIMATE algorithm to assessed glioma purity as well as stromal and immune scores. The landscape of clinical and molecular characteristics is described in three data sets, respectively (Figure 2A). As our previous reports show, the expression of SDC1 increased with the increase of tumor grade. Consistent with recent studies, Stromal and immune scores significantly increased along with malignancy progression, whereas glioma purity decreased in higher grades (Zhang et al., 2017). The IDH mutation and 1p/19q codeletion take place more in the low-grade gliomas, but the MGMT promoter mutation did not display this rule. Besides this, the incidence rate of men is higher than that of women. Most important of all, we found that SDC1 expression is a positive correlation with immune scores (Figure 2B). The correlation coefficient between SDC1 and immune scores is .4, .35, and .28 in the mRNAseq_325, mRNAseq_693, and TCGA RNAseq data sets, respectively.

The ssGSEA was introduced to quantify the relative infiltration of 28 immune cell types in the TME in the three data sets. We calculate the correlation of SDC1 with each immune cell (Figure 3A); the correlation coefficient between SDC1 and activated CD4+T cells is .71, .57, and .63 in mRNAseq_325, mRNAseq_693, and TCGA RNAseq data sets, respectively. The correlation coefficient between SDC1 and activated CD8+T cell is .47, .43, and .38 in mRNAseq_325, mRNAseq_693, and TCGA RNAseq data sets, respectively (Figure 3B).

The DEG analysis between the high and low expression of SDC1 of glioma samples.

The DEGs in the three data sets were calculated though the LIMMA package. The best homogeneity between tumor samples in the high and low expression of SDC1 of glioma samples is shown in the TCGA data sets. The volcano maps show the multiple changes of different genes more clearly (Figure 4A). The black points represent the genes that abs (logFc) < 1, p-value < .05; the red and blue points represent the genes which 1 < abs (logFc) < 3, p-value < .05; and the green points represent the genes which 3 < abs (logFc), p-value < .05. To get more reliable differential genes within the DEG data sets, the unique genes sets were determined using Draw Venn Diagram (Figure 4B). The number of 8,887 genes was selected to WGCNA in TCGA-seq data sets in the next data analysis.

The GBM samples with clinic characteristics were selected to find the gene sets not only related to SDC1 but also tumor immune infiltration. After removing outlier samples, 145 samples were used to construct an adjacency matrix (Figure 5A). Mapping of clinical trait variables and aggregation trees display that WGCNA analysis can be performed. The dynamic tree display different modules were constructed (Figure 5B). The soft thresholding power β was set to eight to ensure a correlation coefficient close to .9 (Figure 5C). A total of 22 different color-coded co-expression modules were identified (Figure 5D). A TOM plot showed that each module in the network was independent of each other, further indicating that gene expression in each module was also relatively independent (Figure 5E). After combining the module conservatism test, the Z-scores >10 of the module were selected. No modules are excluded. Subsequently, we analyzed the module related to both SDC1 and immune invasion in glioma. The blue module (contained 685 genes) was significantly correlated with immune scores and correlates with SDC1; therefore, this module was selected as a significant module for further analyses. Besides this, we painted scatterplots of MM versus gene significance (GS) of SDC1 and immune scores, and a correlation was obtained between the blue module and SDC1 (Figures 5F,G). The genes of the hub modules with the top 50% intro module connectivity, MM, and GS were selected as hub genes.

The hub genes in the blue module were selected for GO analysis, the BP group genes were mainly associated with activation of lymphocytes, and MF group genes were mainly enriched in the binding of antigen and antibody. We selected the group genes that correlated with T cells as target genes, which include 55 genes (Figure 6A). To explore potential signaling pathways relating SDC1 to immune scores, KEGG pathways were employed to identify the function of the signaling pathway. We found eight genes correlated with antigen processing and presentation (Figure 6B). The eight genes are CTSL, HLA-DPA1, HLA-DPB1, HLA-DRB1, CD74, HLA-DOA, HLA-DMB, and CD4.

To determine the associations between target genes in human glioma and increases in tumor progression and poor prognosis, the LGG and GBM data from TCGA data sets were analyzed. The target gene expression levels (Figures 7A–H) were significantly correlated with more advanced tumor grades (Kruskal–Wallis nonparametric test, p < .0001 for both). Moreover, analysis of TCGA patient survival revealed that elevated target gene expression predicted prognosis (Figures 8A–H). Glioma patients with a high level of target genes above the median had dramatically decreased survival compared with those with target gene levels below the median (Kaplan–Meier survival analysis, p < .0001 in log-rank). These data confirm that target genes play a crucial role in the tumorigenesis of glioma.

To explore the protein interaction between SDC1 and the target genes as well as the interaction with other genes, the GeneMANIA (https://genemania.org) and String online database (https://string-db.org) was used to perform PPI analysis, and the results showed that SDC1 may interact with CD4, CTSL, and HPSE directly or be indirectly involved in the regulation of T cell activation (Figure 9).