Potential common mechanisms between primary Sjögren's syndrome and Hashimoto's thyroiditis: a public databases-based study

Abstract

Objective: Primary Sjögren's syndrome (pSS) and Hashimoto thyroiditis (HT) can occur in the same patient population, but the mechanism of co-occurrence remains unknown. This study aims to explore the underlying mechanism.

Methods: We screened differentially expressed genes (DEGs) in the pSS and HT-related transcriptomic microarrays. Based on KEGG, PID, Reactome, and BioCarta enrichment analysis, pathway annotations were performed. A PPI network was developed using STRING. Betweenness, BottleNeck, MNC, Radiality EPC, and Stress topological analyses were performed to identify hub genes. Then, we used two more datasets to validate the key genes. Immune infiltration landscape of pSS and HT were profiled based on CIBERSORT, Xcell, MCPCounter, and EPIC. Correlation between T/B cells and key genes was performed. Single gene GSEA analysis was performed to further explore enriched pathways of key genes. Finally, we predicted the drugs of key genes and the cross-talk genes targeted in the protein domain.

Results: A total of 93 cross-talk genes were found. These genes were mainly related to the immune system. STAT1, CD8A, and PTPRC were identified as hub genes using six topological methods. STAT1 and PTPRC are considered key genes after in silico validation. STAT1 and PTPRC were linked to CD8⁺ Tcm and other immune cells in the pSS and HT dataset. GSEA analysis showed that STAT1 and PTPRC may play a role in pSS and HT through several pathways, including IFNγ response, IFNα response, allograft rejection, E2F targets, complement, G2M checkpoint, IL6-JAK-STAT3 signaling, KRAS signaling up, IL2-STAT5 signaling, IL6-JAK-STAT3-signaling, and inflammatory response. Guttiferone K and picoplatin may be the candidate drugs for the treatment of pSS and HT. Cross-talk genes were mainly enriched in IGc1, MHCIIα and SCY.

Conclusion: We analysed databases and gene expression data for pSS and HT. We identified two genes (STAT1, PTPRC) as potential biomarkers of disease activity in pSS and HT. We also gained new insights into the cellular and molecular mechanisms associated with pSS and HT. Based on the key genes and cross-talk genes, we predicted potential drugs and protein domains for pSS and HT.

Keywords: Hashimoto’s thyroiditis; common mechanisms; cross-talk genes; hub genes; primary Sjögren’s syndrome.

FIGURE 1

Study flowchart.

FIGURE 2

Identification of cross-talk genes in pSS and HT. (A) Volcano plot shows up- (red) and down- (blue) regulated DEGs in GSE40611. (B) Heatmap plot GSE40611 shows TOP20 up- and downregulated DEGs. (C) Volcano plot shows up- (red) and down- (blue) regulated DEGs in GSE138198. (D) Heatmap plot GSE138198 shows TOP20 up- and downregulated DEGs. (E) Cross-talk up- and downregulated genes in pSS and HT are plotted using Venn after the intersection of GSE40611 and GSE138198. In volcano plots, FDR <0.05 and |Fold Change| >1.5. The samples in heatmap plots are grouped and clustered according to disease and control using the Euclidean method with FDR <0.05.

FIGURE 3

Functional annotation of cross-talk genes. (A) KEGG enrichment. (B) PID enrichment. (C) Reactome enrichment. (D) BioCatar enrichment. The X-axis represents gene ratio while the Y-axis represents TOP5 enriched pathways. The size of the dots represents the number of genes, and the color of the dots represents the magnitude of the FDR. FDR <0.05.

FIGURE 4

Hub genes screening after protein-protein interaction analysis. (A) TOP10 hub genes screened by the Betweenness method. (B) TOP10 hub genes screened by the BottleNeck method. (C) TOP10 hub genes screened by the MNC method. (D) TOP10 hub genes screened by the Radiality method. (E) TOP10 hub genes screened by the EPC method. (F) TOP10 hub genes screened by the Stress method. (G) Three hub genes (STAT1, CD8A, PTPRC) after intersection of six screening methods using Venn plot.

FIGURE 5

Hub genes in silico validation. (A) STAT1, CD8A, and PTPRC expressions are validated in another pSS dataset of GSE127952. (B) ROC curve of three hub genes in GSE127952. (C) STAT1, CD8A, and PTPRC expressions are validated in another HT dataset of GSE29315. (D) ROC curve of three hub genes in GSE29315. ROC curves reflect the relationship between sensitivity and specificity. The closer the AUC is to 1, the better the diagnostic effect of the variable in predicting the outcome.

FIGURE 6

Immune infiltration landscape and correlation between T/B cells and key genes. (A) Immune infiltration of pSS using the CIBERSORT method. (B) Immune infiltration of HT using the CIBERSORT method. (C) Correlation between T cells and STAT1/PTPRC in pSS. (D) Correlation between T cells and STAT1/PTPRC in HT. (E) Correlation between B cells and STAT1/PTPRC in pSS. (F) Correlation between B cells and STAT1/PTPRC in HT. Correlation analysis is based on Pearson’s method. A red correlation coefficient indicates a positive relationship between T/B cells and STAT1/PTPRC; a blue one, a negative one. The threshold criteria are p < 0.05 and -log10 (p-value) < 1.3.

FIGURE 7

Single-gene GSEA of key genes. (A) Single-gene GSEA of STAT1 in GSE40611 shows two entries. (B) Single-gene GSEA of PTPRC in GSE40611 shows eight entries. (C) Single-gene GSEA of STAT1 in GSE138198 seven entries. (D) Single-gene GSEA of PTPRC in GSE1138198 three entries. (E) Intersection of GSEA enriched entries for STAT1 in pSS and HT shows one entry. (F) Intersection of GSEA enriched entries for PTPRC in pSS and HT shows three entries. The cut-off criterion is p < 0.05.

FIGURE 8

Drug prediction of key genes and prediction of cross-talk genes targeted protein domain. (A) Drug prediction of STAT1 and PTPRC shows the flow path of gene, drug, regulatory approval, indication, and interaction score using a Sankey plot. (B) Protein domain prediction of cross-talk genes of pSS and HT indicates IGc1 (6.5% genes), MHCIIα (3.9% genes) and SCY (6.5% genes), using a bar plot. The cut-off criteria are p < 0.05.