Immune Cell-Related Genes in Juvenile Idiopathic Arthritis Identified Using Transcriptomic and Single-Cell Sequencing Data

Abstract

Juvenile idiopathic arthritis (JIA) is the most common chronic rheumatic disease in children. The heterogeneity of the disease can be investigated via single-cell RNA sequencing (scRNA-seq) for its gap in the literature. Firstly, five types of immune cells (plasma cells, naive CD4 T cells, memory-activated CD4 T cells, eosinophils, and neutrophils) were significantly different between normal control (NC) and JIA samples. WGCNA was performed to identify genes that exhibited the highest correlation to differential immune cells. Then, 168 differentially expressed immune cell-related genes (DE-ICRGs) were identified by overlapping 13,706 genes identified by WGCNA and 286 differentially expressed genes (DEGs) between JIA and NC specimens. Next, four key genes, namely SOCS3, JUN, CLEC4C, and NFKBIA, were identified by a protein-protein interaction (PPI) network and three machine learning algorithms. The results of functional enrichment revealed that SOCS3, JUN, and NFKBIA were all associated with hallmark TNF-α signaling via NF-κB. In addition, cells in JIA samples were clustered into four groups (B cell, monocyte, NK cell, and T cell groups) by single-cell data analysis. CLEC4C and JUN exhibited the highest level of expression in B cells; NFKBIA and SOCS3 exhibited the highest level of expression in monocytes. Finally, real-time quantitative PCR (RT-qPCR) revealed that the expression of three key genes was consistent with that determined by differential analysis. Our study revealed four key genes with prognostic value for JIA. Our findings could have potential implications for JIA treatment and investigation.

Keywords: WGCNA; juvenile idiopathic arthritis; machine learning analysis; protein–protein interaction; single-cell RNA sequencing.

Figure 1

PBMC immune microenvironment and WGCNA of JIA and NC samples. (A) Abundance histogram of the proportion of immune cells. (B) A box plot of comparison of immune cells between JIA and NC samples. (C) Sample clustering dendrogram. (D) Network topology analysis under various soft-threshold and various soft-thresholding powers. The red line represents signed R². (E) A dendrogram of module eigengenes. The red line represents MEDissThres value. (F) A gene dendrogram with different similarities and module colors. (G) A heat map of the correlations between 18 modules and JIA-related immune cells. p values are shown as * p < 0.05; ** p < 0.01; *** p < 0.001.

Figure 2

GO and KEGG enrichment analyses of DE-ICRGs in JIA. (A) A volcano plot of the DEGs. (B) A heat map of the DEGs. (C) A Venn diagram of genes identified by WGCNA and DEGs. (D) GO analysis of DE-ICRGs. (E) KEGG analysis of DE-ICRGs.

Figure 3

DE-ICRGs identified through PPI network and machine learning analysis. (A) PPI network of 168 shared genes, i.e., DE-ICRGs. (B) The top 20 genes exhibiting the highest degree of interaction in the PPI network. (C) LASSO regression model with gene coefficients, in which lines with different colors represent different genes. (D) LASSO regression model with 10-fold cross-validation. (E) ROC curves of the random forest model. (F) Feature importance scores of the random forest model. (G) Kappa of the RFE model. (H) Venn diagram of genes identified using the LASSO, random forest, and RFE models.

Figure 4

Validation of the hub genes and single-gene GSEA. (A) ROC curves of hub genes identified from the data sourced from GSE13501. (B) ROC curves of hub genes in identified from the data sourced from GSE112057. (C) ROC curves of logistic regression models for the key genes. (D–G) GSEA of CLEC4C, JUN, NFKBIA, and SOCS3. GSEA, gene set enrichment analysis; ROC, receiver operating characteristic.

Figure 5

ScRNA-seq analysis of JIA and NC samples. (A) Violin plots following quality control and normalization. (B) Distribution of cell counts and intracellular gene counts following standardization. (C) A scatter plot of 2000 highly mutative genes and the top 10 among them. (D) PCA without distinct separations of cells in samples. (E,F) Elbow and JackStraw plots of PCs after linear dimensionality reduction. (G) A t-SNE plot of cell distribution in eight samples following t-SNE dimensionality reduction.

Figure 6

Clustering of immune cells and key gene expression patterns. (A) Clustering analysis and the constitution of clusters in JIA and NC. (B) A heat map of marker genes within nine clusters. (C) A heat map of cell annotation. (D) t-SNE visualization following cell subset annotation. (E) Dot plot of key gene expression patterns.

Figure 7

Cellular distribution of the key genes and cell quasi-temporal trajectory analysis. (A) Distribution of key genes in cells. (B) Comparison of the distribution of key genes in cells between JIA and NC samples. (C,D) Pseudotime analysis of the overall temporal alterations of cell populations. (E–H) The expression of CLEC4C, JUN, NFKBIA, and SOCS3 as revealed by pseudotime analysis.

Figure 8

CeRNA network, IPA, and prediction of potential drugs targeting the key genes. (A) CeRNA network of key genes. Circles represent mRNA (red indicates upregulated expression, and green indicates downregulated expression); diamonds represent lncRNAs (purple indicates linkage = 1, and pink indicates linkage >1); and triangles represent miRNAs (orange indicates linkage = 2, and yellow indicates linkage >2). (B) The top 10 pathways as revealed by IPA. (C,D) The most significant upstream regulators of activation and inhibition and their pathways of action. (E,F) Identification of biological functions associated with the top 10 pathways as revealed by IPA. (G) The drug–gene interaction network.