The GSE55235 (13), GSE1919 (14), GSE48780 (15), GSE55457 (13), and GSE55584 (13) datasets were retrieved and downloaded from the Gene Expression Omnibus(GEO, https://www.ncbi.nlm.nih.gov/geo/). Limma package was used to normalize the data and convert between probe ID and gene symbol through platform information, removing probes without gene symbol and taking the average expression value of multiple probes under the same symbol (16). We used variation coefficients to select the most significantly varied genes for analysis.

The proportion of 22 immune cell types and the infiltration of immune cells in samples from GSE55235, GSE1919, GSE48780, GSE55457, and GSE55584 datasets was calculated using the CIBERSORT algorithm (17).

“WGCNA” was used to construct a gene co-expression network for genes in the GSE55235 dataset with a variation coefficient greater than 0.1. The top 25% of highly expressed variants were analyzed (18). The reliability of the constructed scale-free network is ensured by removing abnormal samples. For that, they were used to approximate appropriate soft threshold rates and obtain adjacency values between genes whose variances were more significant before applying power functions. The adjacency values were then converted into topological overlap matrices (TOM) and derived the dissimilarity (1-TOM) values. The dynamic tree-cutting method was finally used to identify modules by hierarchical clustering of genes.

Each module underwent component analysis according to its features, and the correlation between module characteristics and T-cell subtypes was then analyzed using the Pearson test. We considered the module significantly correlated with T cell subtypes when P<0.05. The central module was defined as that with the highest correlation coefficient with CD8+T cells.

Based on the module connectivity and clinical characteristic relationships of each gene in the central module, candidate central genes were chosen. Pearson’s correlation between genes’ absolute values is used to describe module connectivity. The overall Pearson’s correlation between each gene and each characteristic was used to identify the relationship between clinical traits. For candidate central genes, we set gene significance >0.75. Verify whether there is a significant correlation between the candidate hub genes and CD8+T, to determine the reliable hub genes. GSE1919, GSE48780, GSE55457, and GSE55584 datasets were used to verify the Spearman correlation between CD8+T cell and central gene expression. Hub genes were defined as candidate hub genes that were significantly associated with CD8+T cells in the GSE55235, GSE1919, GSE48780, GSE55457, and GSE55584. In addition, CIBERSORT, MCPCOUNTER, QUANTISEQ, TIMER, and XCELL algorithms were used to verify the correlation between genes and CD8+T cells in different tumors.

Critical genes for T cell CD8+ were identified using the LASSO regression and random forest algorithms. Least Absolute Shrinkage and Selection Operator (LASSO) regression is performed using the “glmnet” package. RF is implemented using the “randomForest” package (19). A Venn diagram was used to visualize the results of the two algorithms and obtain the intersection genes.

We examined genes using Gene Ontology (GO) analysis (20), which comprised molecular function (MF), cellular component (CC), and biological process (BP), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (21).

Genes GSE55235, GSE1919, GSE48780, GSE55457, and GSE55584 were correlated, and the analysis results were then visualized using “ggplots2”. We obtained high-throughput expression data of cancer and normal tissues from the TCGA (22) and GETx databases (http://commonfund.nih.gov/GTEx), respectively, and evaluated the correlation between genes in various tumors and various normal tissues to further illustrate the correlation between central genes.

Transcription factors are involved in various complex biological processes by regulating the transcription process through specific DNA sequence recognition. RNA-seq data were obtained from TCGA database and we used ChEA3 (https://maayanlab.cloud/chea3/) database to identify the transcription factors of hub gene and demonstrated them using Cytoscape software.

A gene set’s significance between two biological states can be determined using gene set enrichment analysis (GSEA) (23). A single gene GSEA analysis was carried out to better investigate the potential mechanism through which Hub genes affect RA. At the same time, h.all.v7.5.1.symbols.gmt in the Molecular Signatures Database (MSigDB: https://www.gsea-msigdb.org/gsea/index.jsp) was selected as the reference gene set, and Spearman rank correlation coefficient was obtained by using the package “corrplot”, and adjustment P value <0.05 was used as the screening criterion.

GeneMANIA (http://genemania.org/), a website for creating protein-protein interaction (PPI) networks, may be used to make predictions for the functions of genes and find genes with similar effects. Among them, providing physical interaction, co-expression, co-localization, gene enrichment analysis, genetic interaction, and locus prediction are some of the bioinformatics methods used by network integration algorithms.

Microarray expression data of rheumatoid arthritis (GSE55235) were downloaded from the GEO database, and genes with a variation coefficient greater than 0.1 were selected for subsequent analysis. We next used the R package “CIBERSORT” to analyze the corresponding expression data in the dataset GSE55235 to ascertain the percentage of various immune cell subtypes in each sample from the dataset. The proportion of each block’s seven T cell subtypes was then chosen as the trait data for the WGCNA.

Based on genes with a variation coefficient greater than 0.1 from the GSE55235 dataset, the R package “WGCNA” was used to construct a gene co-expression network. The average association coefficient and Pearson correlation coefficient were then calculated, and cluster analysis on all samples in the dataset was performed. A scale-first network was built using β=7 as the soft threshold rate ( Figures 1A, B ). To create hierarchical clustering trees, the dynamic hybrid cutting technique was employed. Each branch is a module that combines all genes with similar expression levels, and each leaf represents a gene ( Figures 1C, D ). Functionally equivalent modules were then combined into one large module, resulting in 37 modules ( Figures 1E, F , 2A ).

Among the 37 modules, the blue module had the most significant relationship with CD8+T cell (R^{2 =} 0.69, P=2E-05). We designate the blue module as the hub module as a result. In this module, potentially important genes that are most closely related to CD8+ T cells were shown to have a stronger correlation. Gene with gene significance >0.75 were selected as candidate central genes, and a total of 10 genes (NKG7, CD8A, DHRS9, CCL5, IL2RG, TNS3, PARP12, GZMA, PRF1, CYTH4) were selected ( Figure 2B ).

The genes in the blue module were then subjected to GO/KEGG analysis. The findings demonstrated that genes for immune response-activating cell surface receptor signaling pathway, antigen receptor-mediated signaling pathway, and T cell activation in BP were highly enriched in the blue module. The external side of the plasma membrane, MHC protein complex, and secretory membrane were significantly enriched in CC. While cytokine receptor activity, MHC protein complex binding, and SH2 domain binding were significantly enriched in MF ( Figure 2C ). Module genes were significantly enriched in Epstein-Barr virus infection, Chemokine signaling pathway, Antigen processing, and presentation in KEGG ( Figure 2D ).

We performed LASSO and random forest analysis on the above-mentioned 10 genes, and a total of 6 genes (NKG7, CD8A, TNS3, PARP12, GZMA, PRF1) were obtained by LASSO analysis ( Figure 3A ). Genes of Importance >0.6 were selected from RF analysis results for analysis (CD8A, DHRS9, CCL5, IL2RG, TNS3, GZMA, PRF1)( Figure 3B ). Next, a Venn diagram was used to visualize the intersection genes of LASSO and RF, and a total of 4 genes (CD8A, TNS3, GZMA, PRF1) were found ( Figure 3C ). GSE55235 ( Figures 4A–C ), GSE55457 ( Figures 4D–F ), GSE55584 ( Figures 4G–I ), GSE1919 and GSE48780 ( Figures 4J–L ) were used to identify and validate the correlation between the levels of CD8A, TNS3, GZMA, PRF1, and CD8+ T-cell infiltration. CD8A, GZMA, and PRF1 were identified as reliable Hub genes. The results of the analysis revealed that three genes were significantly positively correlated with the degree of CD8+ T-cell infiltration in these data sets. At the same time, we selected different cancers for analysis ( Figure 4 ). CIBERSORT, McP-counter, QUANTISEQ, TIMER, and XCELL were employed to explore the expression values of genes related to CD8+T cells in different cancers and utilized the “ggPlot” R package for visualization. It was discovered that CD8A ( Figure 5 ), GZMA ( Supplementary Figure 1 ), and PRF1 ( Supplementary Figure 2 ) were positively correlated with CD8+T cell infiltration in various malignancies. These analyses verified that the identified hub genes play a significant role in the tumor immune microenvironment and are highly correlated with the degree of CD8+ T-cell infiltration.

We verified the correlation between hub genes in different data sets and found that CD8A, GZMA, and PRF1 were significantly positively correlated with each other. The relationship between CD8A, GZMA, and PRF1 in various malignancies was then examined ( Figure 6 ). Except for THYM, CD8A, and GZMA being positively correlated in all tumors, CD8A and PRF1 were positively correlated in 33 tumors, and GZMA and PRF1 were positively correlated in 33 tumors ( Figures 7A–C ). The results of normal tissue analysis showed that CD8A and GZMA were positively correlated in various normal tissues, GZMA and PRF1 were positively correlated in various normal tissues, CD8A and PRF1 were negatively correlated in Bone Marrow, and the rest were positively correlated ( Figures 7D–F ).

We identified 8 common transcription factors (EOMES, TBX21, STAT4, ZNF80, GFI1, SCML4, ZNF831, RUNX3) associated with CD8A, GZMA, and PRF1 from the CHEA3 database ( Figure 8 ).

To explore the function of the hub gene, we carried out a single gene GSEA analysis. It was found that CD8A was mainly enriched in Rheumatoid arthritis, Asthma, Cholesterol metabolism, etc ( Figure 9A ). GZEA is mainly enriched in cholesterol metabolism, steroid biosynthesis, mineral absorption, etc ( Figure 9B ). PRF1 was mainly enriched in Coronavirus disease COVID-19, the Pentose phosphate pathway, and Linoleic acid metabolism ( Figure 9C ). We then used h.all.v7.5.1.symbols. gmt as a reference gene set to analyze the correlation between genes and genomes. CD8A, GZMA, and PRF1 were found to be positively correlated with interferon gamma response, interferon alpha response, inflammatory response, il6 jak stat3 signaling, complement, and allograft rejection, and negatively with uv response dn, pancreas beta cells, myogensis, cholesterol homeostasis, androgen response, and adipogenesis.

GeneMANIA database was employed to create PPI networks for Hub genes. To further investigate the function of hub genes, we constructed a 20-gene interaction network ( Figure 10A ). At the same time, 20 genes were investigated for functional enrichment. Results showed that, in BP, Genes were mainly enriched in T cell-mediated immunity, modification of morphology, or physiology of other organisms. Genes were mainly enriched in immunological synapses and the external side of the plasma membrane under CC. In MF, genes were mainly enriched in phospholipase activator activity and lipase activator activity ( Figure 10B ). Genes in the KEGG pathway may regulate Primary immunodeficiency, Viral protein interaction with cytokine, and cytokine receptor. These findings lead us to hypothesize those hub genes are crucial for the immune system ( Figure 10C ).