Identification of biomarkers related to CD8⁺ T cell infiltration with gene co-expression network in clear cell renal cell carcinoma

RNA expression data

The research strategy is presented in Figure 1.

We obtained RNA expression data for 265 ccRCC samples Gene Expression Omnibus (GEO) database. All data for tumor samples in the dataset were obtained, and 4411 genes with Coefficient of variation values greater than 0.1 were selected for additional analysis (Supplementary Table 1).

Figure 1. The workflow of the study.

Evaluation of tumor-infiltrating immune cells (TIICs)

CIBERSORT is an analytical algorithm that analyzes RNA expression data to assess the abundance of different cell subtypes for each sample. The fractions of 22 TIICs were calculated by using the R package “CIBERSORT”. Then, the fractions of seven subtypes of T cells in every sample were selected as trait data of WGCNA (Supplementary Table 2).

Gene co-expression network of ccRCC

The expression values of the 4411 genes were used to construct a co-expression network using the R package “WGCNA” We calculated average linkage and Pearson’s correlation values to cluster the samples of GSE73731 (Supplementary Figure 1). To build a scale-free network, we picked β = 3 (scale free R² = 0.8723676) as the soft-thresholding power (Figure 2A, 2B). A hierarchical clustering tree was constructed using dynamic hybrid cutting. Each leaf on the tree represents a single gene, and genes with similar expression data are close together and form a branch of the tree, representing a gene module. Nine modules were generated (Figure 2C).

Figure 2. Selection of the appropriate beta value to construct a hierarchical cluster number. (A) Analyze the scale-free fit index of the 1-20 soft threshold power (β). (B) Analyze the average connectivity of 1-20 soft threshold power. (C) Genes are grouped into various modules by hierarchical clustering, and different colors represent different modules.

Identification of hub modules and enrichment analysis

Among the nine module, the green module was highly correlated to T cells CD8 (CD8⁺ T cells) (R²=0.5, P=2e-18), T cells CD4 memory activated (R²=0.45, P=1e-14), and T cells gamma delta (R²=0.62, P =3e-29), and the yellow module showed higher correlation with activated T cells CD4 memory activated (R²=0.5, P=4e-18; Figure 3A). The correlation between other modules and T cells was less than 0.5. We were interested specifically in the CD8⁺ T cells, so focused on the green module that showed correlation with CD8⁺ T cells was identified as a hub module. Genes included in this module were next analyzed using the web tool “Matascape” for pathway and process enrichment analysis. The 20 highest enrichment terms were all immune-related terms (Figure 3B), and the three most highly enriched terms were Lymphocyte Activation, Adaptive Immune Response, and Cytokine-mediated Signaling Pathway (Supplementary Table 3).

Figure 3. Key modules and feature notes. (A) Heatmap shows correlations of module eigengenes with T-cell infiltration. (B) The first 20 enriched terms are shown as a bar chart on the left. The network diagram on the right is constructed with each enrichment term as a node and the similarity of the node as the edge. Nodes with the same cluster ID are the same color.

Identification and validation of hub genes

The highly connected genes of the module were investigated as potential key factors related to CD8⁺ T cell infiltration level. According to the cut-off standard (Module-Membership > 0.8 and Gene-Significance > 0.5), 30 genes were selected as candidate hub genes (Figure 4B). From the protein-protein interactions (PPI) network, cut-offs of reliability > 0.7 and connectivity > 15 (node/edge) were applied to identify 30 genes as central nodes, and we visualized these results using Cytoscape (Figure 4A). Ten genes were selected in both analyses designated as hub genes (LCK, CD2, CD3D, CD3G, IRF1, IFNG, CCR5, CD8A, CCL5 and CXCL9) (Figure 4C). To investigate the relationship between these hub genes and CD8⁺ T cells, we analyzed the expression data for these genes in the TIMER database. The results showed positive correlation of the expression values of the 10 genes with the infiltration levels of CD8⁺ T cells (correlation coefficient of at least 0.75 for all genes except CXCL9) (Figure 5A). As an example, we show a scatter plot of CCL5 expression and CD8⁺ T cell infiltration level in Figure 5B. We next queried the TISIDB database to obtain the Spearman correlation values between abundance of tumor-infiltrating lymphocytes and gene expression. The results show a positive correlation between hub genes and tumor-infiltrating lymphocytes. The correlation values were highest for Activated CD8⁺ T cell (Act CD8) and Effector memory CD8⁺ T cell (Tem CD8) (Figure 5C). These analyses verified the identified hub genes as strongly associated with the CD8⁺ T cell infiltration level and playing significant roles in the immune microenvironment.

Figure 4. Identification of hub genes. (A) PPI network of genes from the green module. The higher the number of connected nodes, the larger the size of the node. The green nodes represent a central node with more than 15 connections. (B) A scatter plot of the genes in the green module. Each green dot represents a gene, and dots within the red box indicate genes of Module Membership > 0.8 and Gene Significance > 0.5. (C) Hub genes were selected based on overlap between PPI and co-expression networks.

Figure 5. Validation of hub genes and PPI map construction. (A) Relationship between 10 hub genes expression and CD8⁺ T cell infiltration level; P < 0.05 is considered statistically significant. (B) Scatter plot of CCL5 expression and CD8⁺ T cell infiltration level. (C) The heat map shows the correlation between the ten identified hub genes and the TIICs from the TISIDB database. The redder color indicates a higher correlation, and greener color indicates a lower correlation. (D) Protein-protein interaction map of the ccRCC immune microenvironment.

Determination of immune and clinical characteristics

We searched for the Spearman correlations of expression of these 10 hub genes with the expression of immune factors in TISIDB database, including immune-inhibitory factors, immune-stimulatory factors, chemokines, and receptors. These correlations are presented as heat maps (Supplementary Figure 2). We identified 38 immune-related factors with average correlations with the 10 hub genes greater than 0.5. We constructed an immune infiltration interaction network based on the 10 hub genes and the 38 immune-related factors to explore the infiltration mechanism of CD8⁺ T cells using the STRING database. The results were imported into Cytoscape for visualization (Figure 5D). We next obtained the expression levels of the 10 genes of ccRCC from the The Cancer Genome Atlas (TCGA). The expression levels of these genes were higher in tumor tissues than in normal tissues (P < 0.05) based on Wilcoxon signed-rank test (Figure 6A–6J). Volcanic map also showed that the expression of 10 genes in tumor tissues was higher than that in normal tissues (Figure 6K), with corrected p values of all genes less than 0.05, indicating statistical significance. The fold changes of CCL5 and CXCL9 were more than 2.5-fold higher in tumor tissues than the levels in normal tissues (Supplementary Table 4). The relationships between hub genes and pathological stages are shown by boxplot (Figure 7A). Expression levels of all hub genes show significant differences in pathological stages (p < 0.05) and showed an upward trend with increased stages. Finally, we investigated the connection between tumor grades and hub genes (Figure 7B), in which LCK, CD2, CD3D, IFNG, CD8A, and CCL5 showed significant correlation with different tumor grades (p < 0.05), with grade increase corresponding to increased gene expression. Although no significant difference was detected for CD3G, IRF1, CCR5, and CXCL9, the gene expression level showed an upward trend with increased tumor grades.

Figure 6. Differential expression of the hub genes in transcriptional data of TCGA. (A) LCK, blue dots represent normal tissue and red dots represent tumor tissue. The y-axis shows the expression value of the gene. (B) CD2. (C) CD3D. (D) CD3G. (E) IRF1. (F) IFNG. (G) CCR5. (H) CD8A. (I) CCL5. (J) CXCL9. (K) The volcano plot of differentially expressed genes. Red dots indicate overexpression genes, green dots indicate low expression genes, and black circles represent hub genes.

Figure 7. Analysis of hub genes and clinical indicators in the TCGA dataset. (A) Box plot of the hub genes for different pathological stages. (B) Box plot of the hub genes for different tumor grades.

Identification of prognostic biomarker

We analyzed 10 hub genes by Kaplan-Meier analysis, and only the results of CCL5 and IFNG were statistically significant (p < 0.05). For these two genes, the survival prognosis of patients with high expression was poor (Figure 8A, 8B). To validate differential expression of 10 genes in tumor and normal tissues, we used Oncomine to perform a meta-analysis of five analyses using four data sets, all of which included both tumor tissues and normal controls for 10 genes (Figure 8C). Data for IFNG were not included in these data sets, so we obtained meta-analytical results for the other nine genes. LCK, CD2, CD3D, CCR5, CCL5, and CXCL9 have Median Rank values less than 1000, CD3G is 3007.0, IRF1 is 1445.0, and CD8A is 1577.0. The results show these genes exhibit significant overexpression in tumor tissues which was consistent with the TCGA datasets analyses. Through Kaplan-Meier and Oncomine meta-analysis, we selected CCL5 as a prognostic biomarker for further analysis.

Figure 8. Kaplan-Meier analysis and Oncomine meta-analysis. (A) The overall survival analysis of CCL5. (B) The overall survival analysis of IFNG. (C) A meta-analysis of gene expression from Oncomine datasets. Colored squares represent the median of genes (relative to normal tissue) in five analyses. Red represents overexpression, blue represents low expression. This P value gives the average rank analysis.

Gene set enrichment analysis of CCL5

According to CCL5 expression median value, ccRCC samples from TCGA were divided into high expression group and low expression group for pathway gene set enrichment analysis. The enrichment results showed that immune-related pathways were enriched in the high expression group, with a total of 23 pathways statistically significantly enriched (p-value < 0.05, q-value < 0.05). The three most enriched pathways were “Antigen procession and presentation”, “Cell adhesion molecules cams”, and “Autoimmune thyroid disease” (Figure 9A, 9B, Supplementary Table 5). There were no significantly enriched pathways for the low expression group.

Figure 9. GSEA and experiment of CCL5. (A) The above picture shows the enrichment fractional broken line of the three pathways, and the lines in the lower figure correspond to the genes of each pathway. (B) The circle diagram shows three enrichment pathways and the core genes that play a role in the enrichment process. The larger the circle corresponding to each gene, the larger the rank metric score value. (C) The relative mRNA levels of normal control and si-CCL5 of 769-P cells. (D) The results of CCK-8 assay showed decreased proliferation ability of 769-P cells when treated with si-CCL5. (E) The invasion assay results showed decreased invasion ability of 769-P cells treated with si-CCL5.

Knockdown of CCL5 significantly inhibited cell proliferation and invasion

Because CCL5 is overexpressed in ccRCC and is related to poor prognosis, we next conducted functional experiments to explore the potential biological function of CCL5. First of all, we used small interfering RNAs to lower expression levels of CCL5 in the renal cell carcinoma cell line 769-P (Figure 9C) and evaluated the cell proliferation ability by Cell Counting Kit-8 (CCK-8). The results showed decreased cell line proliferation ability after CCL5 knockdown (Figure 9D). Finally, the invasion ability of the cell line decreased significantly after CCL5 knockdown (Figure 9E).

Gene expression data and processing

We downloaded ccRCC RNA expression data from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/), which contains data related to 265 samples [33]. The dataset of GSE73731 was obtained using the platform Affymetrix Human Genome U133 Plus 2.0 Array (HG U133 Plus 2.0). We used the R package “limma” [34] to normalize the RNA-sequencing data. Small variation of gene expression data often represents noise, so we used Coefficient of Variation values to select the most variant genes, which were then used to construct the network.

Evaluation of tumor-infiltrating immune cells

In this study, we used the R package “CIBERSORT” to estimate the fraction of immune cells of GSE73731 samples. Specifically, the CIBERSORT algorithm was used to calculate the fractions of the 22 types of TIICs [20]. CIBERSORT is considered better to previous deconvolution methods for analysis of unknown mixture content and noise. This algorithm can be used to statistically estimate the relative proportions of cell subpopulations from complex tissue expression profiles, making it a useful tool to estimate the abundances of special cells in mixed tissue.

Co-expression network construction

Expression values of 4411 genes were used to construct a weight co-expression network using the R package “WGCNA” [17]. First, based on the Pearson’s correlation value between paired genes, the expression levels of individual transcripts were converted into a similarity matrix. Next, the similarity matrix was transformed to an adjacency matrix, as calculated by amn = |cmn| β (cmn = Pearson’s correlation between paired genes; amn = adjacency between paired genes). Parameter β can improve strong correlations and decrease weak correlations between genes. The adjacency matrix was then converted into a topological overlap matrix when the power of β = 3. To categorize genes with similar expression patterns into different modules, we applied a dynamic hybrid cutting method, using a bottom-up algorithm with a module minimum size cutoff of 30.

Construct module trait relationships

Module eigengenes were used to perform component analysis of each module. We calculated the correlation between module eigengenes and the infiltration level of T cells to determine the significance of modules by Pearson test. An individual module was considered significantly correlated with T cells when p < 0.05. We selected the interest T cell subtype and module with the highest correlation coefficient and defined that as a hub module.

Pathway and process enrichment analysis

To determine the function of genes in the identified hub module, we used the web tool “Metascape” (http://metascape.org) for pathway and process enrichment analysis [35]. The tool displays the first 20 enriched terms as a bar graph. In order to further explore the relationship between terms, terms with similarity greater than 0.3 are connected by edges and presented as a network graph.

Identification of hub genes

We selected candidate hub genes based on the modular connectivity and clinical traits relationship of each gene in the hub module. Module connectivity is defined as the absolute value of the Pearson’s correlation between genes (Module Membership). Clinical trait relationship is defined as the absolute value of the Pearson’s correlation between each gene and the trait (Gene Significance). We set the Module-Membership > 0.8 and the Gene-Significance > 0.5 for candidate hub genes. Meanwhile we selected all genes in the hub module and used the Search Tool for the Retrieval of Interacting Genes (STRING; https://string-db.org/) database to construct PPI network and looked for central nodes [36]. Genes with node connectivity > 15 were considered central nodes. We used Cytoscape to present the network (https://cytoscape.org/) [37]. We did Venn analysis to compare candidate hub genes and central nodes in the PPI network using the online tool (http://bioinformatics.psb.ugent.be/webtools/Venn/).

Validation of hub genes

We used two immune-associated databases that are based on TCGA to validate these hub genes, as described below. First, we obtained the content of CD8⁺ T cells in each sample of ccRCC based on data in the Tumor Immune Estimation Resource (TIMER; https://cistrome.shinyapps.io/timer/) [38]. Spearman correlations between the infiltration level of CD8⁺ T cells and the expression of hub genes were calculated, and the results were compared using the R package “ggstatsplot”. Second, the Tumor Immune System Interactions Database (TISIDB; http://cis.hku.hk/TISIDB) [39] was then searched to determine Spearman correlations between the hub genes and TIICs. We present these results as a heat map constructed using the R package “pheatmap”.

Immune and clinical characteristic identification

Spearman correlations between hub genes and different immune factors were obtained from the TISIDB database, which includes immune-inhibitory and immune-stimulatory factors, chemokines, and receptors. We then constructed a heat map by using the R package “heatmap”. Immune factors related to hub genes with average correlation coefficients greater than 0.5 were picked to construct a network using STRING and Cytoscape [37]. To explore the clinical characteristic of hub genes, RNA expression data and clinical data of ccRCC were obtained from The Cancer Genome Atlas (TCGA; https://cancergenome.nih.gov/). The R packages “limma” and “beeswarm” were used to construct a scatter differential diagram of the data. The statistical significance of differences in expression between normal and tumor samples was analyzed using the Wilcoxon signed-rank test. The differences of all coding genes were analyzed by R package “limma", and the volcano map was drawn using the R package “ggplot2”. Finally, a boxplot was constructed to display the relationship between genes and clinical features, and the statistical significances were analyzed by the Kruskal-Wallis test.

Kaplan-Meier analysis and Oncomine mata-analysis

We found best separation through the R package “survminer” to divide the patients into high and low expression groups, such grouping minimizes the p value of the survival curve. Then took Kaplan-Meier analysis for the groups by R package “survival”. To validate the expression patterns of the hub genes, four independent microarray datasets were used from the Oncomine Cancer Microarray database (Oncomine, https://www.oncomine.org) to perform meta-analysis [40–43].

Gene set enrichment analysis (GSEA)

GSEA is a computational method used to determine whether a set of basically defined gene sets exhibits statistically significant differences between two biological states [44]. According to the median value of gene expression, the samples were divided into two groups, and “c2.cp.kegg.v7.0.symbols” gene set enrichment analysis was carried out, with p-value < 0.05 and q-value < 0.05 as indicative of statistical significance. The enrichment pathway was visualized using the R packages “ggplot2” and “clusterProfiler”. In the circle diagram, we only show the core genes of the enrichment process, and the size of the dots represent the Rank metric score, which indicates the amount of sequencing value of the gene.

Cell culture and siRNA-PTEN construction

The 769-P cells were cultivated in DMEM with 10% of fetal calf serum, 100 U/ml penicillin and 100 μg/ml streptomycin. The cells were cultured at 37 °C, with 5% CO₂. The CCL5-310 siRNA sequences were as follows: 5’-GCUGAACAAGGGAAGCUUTT-3’ 5’-AAGCUU GCCCUUGUUCAGCTT-3’. The CCL5-286 siRNA sequences were: 5’-GCAGGAUUUCCUGUAUG ACTT-3’ and 5’-GUCAUACAGGAAAUCCUGCTT-3’. The CCL5-240 siRNA sequences were: 5’-UC GUCCACAGGUCAAGGAUTT-3’ 5’-AUCCUUGAC CUGUGGACGATT-3’

RNA extraction and quantitative RT-PCR

RNA samples were extracted from cells with the FastPure Cell/Tissue Total RNA Isolation kit (Vazyme Biotech) and they were reverse-transcribed by HiScript III RT SuperMix for qPCR (Vazyme Biotech). To determine the relative transcript level, PCR was quantified in real-time using a LifeECO PCR machine (BIOER Technology Co, Ltd). SYBR Green was used as the fluorophore. The CCL5 primers sequences were as follows: 5′-CTC ATTGCTACTACTGCCCTCTGCGCTCCTGC-3′ and 5′GCTCATCTCCAAAGAGTTGATGTACTC -3′. The PCR parameters were: 95 °C 5min, followed by 50 cycles of 30 s and 1 min at 60 °C. Each sample was measured in three independent reactions. The threshold values of each sample/primer were determined, and the average error and standard error were calculated. Melting curve analysis was performed, and the mRNA expression levels were normalized against that of β-actin.

Cell proliferation analysis and invasion experiment

769-P cells were cultivated in five 96-well cell culture plates (1500 cells/well) respectively for 5 days. A volume of 10 μl of CCK-8 solution was added to each well in a plate at an interval of 24 h. The cells were cultures for 3 h in 37 °C, and then a photometry test was performed at wavelength of 450 nm. The 769-P cells were then transferred into the top layer holes in a basement membrane plate with wells containing the medium without serum. The lower cavity was filled with 12% fetal calf serum medium to allow chemical attraction. The medium was removed after 12-hour incubation and the implant in the upper cavity was discarded. Cells invading the lower cavity were fixed with 700 ml of 4% PFA and stained with crystal for 1 h to visualize. The cells in the lower cavity were counted and normalized to the number in control conditions to measure the relative invasion capacity.