The AD terms (category C20.111), including 68 diseases, were acquired from the Medical Subject Headings (MeSH, https://www.nlm.nih.gov/mesh/meshhome.html). After removing the complications of ADs, we extracted human disease-related genes from the Genetic Association Database (GAD, https://geneticassociationdb.nih.gov/) (Becker et al., 2004) and mapped these genes to the AD terms for integration. The disease gene sets consisted of 267 related genes of 26 ADs (Supplementary Table S1).

We downloaded the information on protein-protein interactions of human genes from the Human Protein Reference Database (HPRD, http://www.hprd.org/) (Peri et al., 2003) and used Cytoscape software (v3.8.2) (Shannon et al., 2003) to construct a human protein-protein interaction network. The topological properties of AD-related genes in this network were computed (Hidalgo et al., 2009; Chavali et al., 2010). The gene set of disease d was defined as G = {g ₁, g ₂, g ₃, … , g _i , … , g _k }. In order to assessed the network similarity between ADs, we first calculated the average topological properties of each AD in the network as follows: R = ∑ 1 ≤ i ≤ k r i k , C = ∑ 1 ≤ i ≤ k c i k , B = ∑ 1 ≤ i ≤ k b i k , S = ∑ 1 ≤ i ≤ k s i k , a n d H = ∑ 1 ≤ i ≤ k h i k (1-5) where k represents the number of genes in G, g _i is the ith gene of G, r _i is the degree of g _i , c _i is the clustering coefficient of g _i , b _i is the betweenness centrality of g _i , s _i is the average shortest path length of g _i , h _i is the neighborhood connectivity of g _i , R is the degree of d, C is the clustering coefficient of d, B is the betweenness centrality of d, S is the average shortest path length of d, and H is the neighborhood connectivity of d. As shown in Figure 1, d ₁ and d ₂ are two ADs from MeSH, and G ₁ and G ₂ are gene sets related to d ₁ and d ₂. The average topological properties of d ₁ and d ₂ were defined as vector T ₁ = {R ₁, C ₁, B ₁, S ₁, H ₁} and vector T ₂ = {R ₂, C ₂, B ₂, S ₂, H ₂}, respectively. We defined the network similarity (NetSim) score between d ₁ and d ₂ as the Pearson correlation coefficient (PCC) calculated with T ₁ and T ₂. The formula that was used as follows: ρ T 1 , T 2 = cov ( T 1 , T 2 ) σ T 1 , σ T 2 (6) where cov(T ₁,T ₂) is the covariance of variables T ₁ and T ₂, σT ₁ and σT ₂ are the standard deviations for T ₁ and T ₂.

The data on the functional interactions of genes was downloaded from HumanNet, which is a human gene functional interaction network based on Gene Ontology annotation (Lee et al., 2011). Each interaction in HumanNet has a log likelihood score (LLS) that measures the probability of a functional association between genes (Cheng et al., 2014).

We downloaded LLSs between human genes from HumanNet and normalized the LLSs as follows: L L S N ( g i , g j ) = L L S ( g i , g j ) − L L S min L L S max − L L S min (7) where g _i and g _j are the ith and jth gene respectively. LLS _N (g _i , g _j ) indicates LLS between g _i and g _j after normalization. LLS(g _i , g _j ) indicates LLS between g _i and g _j . LLS _min and LLS _max are the minimum LLS and the maximum LLS of HumanNet respectively.

The functional similarity (FunSim) score between a pair of genes was defined as follows: F u n S i m ( g i , g j ) = { 1 L L S N 0 i = j i ≠ j i ≠ j a n d a n d e ( i , j ) ∈ E ( H u m a n N e t ) e ( i , j ) ∉ E ( H u m a n N e t ) (8) where e(i, j) represents the interaction edge between g _i and g _j . E(HumanNet) is a set including all the edges of HumanNet.

Next, we defined the functional association between a gene g and a gene set G = {g ₁, g ₂, g ₃, … , g _i , … , g _k } as follows: F G ( g ) = max 1 ≤ i ≤ k ( F u n S i m ( g , g i ) ) , g i ∈ G (9) where k represents the number of genes in G, g _i is the ith gene of G.

G ₁ = {g ₁₁, g ₁₂, … , g _1i , … , g _1m } and G ₂ = {g ₂₁, g ₂₂, … , g _2j , … , g _2n } are gene sets related to d ₁ and d ₂ respectively. m is the number of genes in G ₁, and n is the number of genes in G ₂. We defined FunSim score of d ₁ and d ₂ as follows: F u n S i m ( d 1 , d 2 ) = ∑ 1 ≤ i ≤ m F G 2 ( g 1 i ) + ∑ 1 ≤ j ≤ n F G 1 ( g 2 j ) m + n , g 1 i ∈ G 1 , g 2 j ∈ G 2 (10)

Semantic similarity is a method to measure the closeness between two terms, according to a given ontology (Zhang and Lai, 2015; Zhang and Lai, 2016; Del Prete et al., 2018). The Resnik method was applied to our study. The human disease terms were obtained from the Human Disease Ontology (DO, http://www.disease-ontology.org) (Schriml et al., 2019). The DO includes the breadth of common and rare diseases, organized as a directed acyclic graph in which a term represents a DO term and an edge represents an “IS_A” relationship between diseases. The information content (IC) of each DO term could be calculated as follows: I C ( d ) = − log ( n N ) (11) where d is a disease term of DO, n is the number of genes related to d, and N is the total number of genes related to DO. As shown in Figure 1, d ₁ and d ₂ are two AD terms of DO, and d _MICA is the most informative common ancestor (MICA) of d ₁ and d ₂. The MICA means the ancestor that has the maximum IC among all the common ancestors between terms of ontology. And we defined the semantic similarity (SemSim) score of d ₁ and d ₂ as follows: S e m S i m ( d 1 , d 2 ) = I C ( d M I C A ) (12)

To identify more reliable pairs of similar ADs, we integrated above three kinds of similarity (NetSim, FunSim, and SemSim) scores to comprehensively determine the levels of similarity between ADs as follows: I n t e g r a t e d S i m ( d 1 , d 2 ) = N e t S i m ( d 1 , d 2 ) ⋅ F u n S i m ( d 1 , d 2 ) ⋅ S e m S i m ( d 1 , d 2 ) 3 (13) where d ₁ and d ₂ are two AD terms of MeSH.

Functional enrichment analysis of Gene Ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) for related genes of similar ADs was performed to infer potential biological processes and pathways using the DAVID Bioinformatics Tool (http://david.abcc.ncifcrf.gov/, version 6.7) (Huang Da et al., 2009). The p-values for the biological processes and pathways were adjusted for false discovery rate (FDR) by the Benjamini-Hochberg method. The biological processes and pathways with FDR less than 0.05 were considered statistically significant functional categories.

To determine whether multiple ADs had high genetic similarity, hierarchical clustering was performed on the integrated similarity scores between 26 ADs based on the Euclidean distance. The d ₁ and d ₂ are two diseases of 26 ADs. We defined the Euclidean distance between d ₁ and d ₂ as follows: E ( d 1 , d 2 ) = ∑ 1 ≤ i ≤ 26 ( x 1 i − x 2 i ) 2 (14) where x _1i is the integrated similarity score between d ₁ and ith disease of 26 ADs, x _2i is the integrated similarity score between d ₂ and ith disease of 26 ADs. The 26 ADs were clustered based on the Euclidean distances between diseases and there were shorter Euclidean distances between ADs belonging to the same disease cluster.

To identify pairs of similar diseases in genetics, we calculated similarity between 26 ADs based on three measurements. All of AD pairs were sorted in descending according to their scores of NetSim, FunSim, and SemSim. To enhance the reliability of the study, we selected ten AD pairs from intersection of top 50 AD pairs with the highest NetSim score, top 50 AD pairs with the highest FunSim score, and top 50 AD pairs with the highest SemSim score for further analysis (Figure 2 and Supplementary Tables S2–4), which were considered as potential pairs of similar ADs.

Further, a network was generated on the ten AD pairs and their related genes, including 14 ADs and 247 genes (Figure 3). The topological properties of the network were investigated. We extracted 35 genes each of whose degree was greater than or equal to three from the network as the possible AD relevant genes. The top five genes regarding degree were HLA-DRB1, HLA-DQB1, CTLA4, HLA-DQA1, and TNF, indicating that these genes were critical in multiple ADs. Previous researches have demonstrated that these genes are associated with many ADs and exert various functions in human autoimmune disorders. Besides, we ascertained that RA and multiple sclerosis (MS) involved more similarity relationships than other AD in this network, implying that the pathogenesis of RA and MS might exist in most of ADs from the network. What’s more, a pair of diseases with higher number of shared genes, suggesting they likely have higher genetic similarity. Thus, the ten potential pairs of similar ADs were ranked by their number of shared genes and shown in Table 1. And the top-rank disease pair is RA and SLE, followed by T1D and RA, and RA and MS.

To explore common genetic mechanisms of a variety of ADs, we performed functional enrichment analysis of GO and KEGG for the 35 AD relevant genes (Figure 4A). The cutoff criterion was a FDR less than 0.05. The top ten significant GO terms in the BP were mainly associated with immune response, antigen processing and presentation, and interferon gamma (IFNγ)-related functions (Figure 4B). Notably, IFNG was involved in the top three GO terms and was defined as a hub gene. IFNG can encode IFNγ that is a cytokine that is critical for innate and adaptive immunity against viral, bacterial, and protozoan infections. And aberrant IFNγ expression is associated with a number of ADs, such as RA and SLE (Hu and Ivashkiv, 2009; Barrat et al., 2019). The top ten significant KEGG pathways contained four AD-correlated pathways, such as “inflammatory bowel disease (IBD),” “autoimmune thyroid disease,” “type 1 diabetes mellitus,” and “rheumatoid arthritis,” which illustrated that the 35 genes might induce the initiation and development of multiple ADs (Figure 4C). HLA-DQB1, HLA-DRB1, HLA-DPB1, HLA-DQA1 were involved in all of the top ten KEGG pathways and were defined as hub genes. The four genes belonged to HLA class II alleles which were suggested to contribute to the susceptibility and resistance to ADs (Wieber et al., 2021).

To identify more reliable pairs of similar ADs, we integrated the three measurements to compute integrated similarity scores between 26 ADs (see Materials and Methods). All of AD pairs were sorted in descending according to their integrated similarity scores. And the top ten AD pairs were extracted for further analysis (Table 2). We found that the ten AD pairs contained three potential pairs of similar ADs consisting of RA and SLE, myasthenia gravis (MG) and autoimmune thyroiditis (AIT), and autoimmune polyendocrinopathies (AP) and uveomeningoencephalitic syndrome (Vogt-Koyanagi-Harada syndrome, VKH) which were defined as the significant pairs of similar ADs.

To reveal the underlying mechanisms shared by two similar ADs, we performed functional enrichment analysis of related genes of RA and SLE, MG and AIT, and AP and VKH. The cutoff criterion was a FDR less than 0.05. The related genes of RA and SLE were significantly enriched in GO terms mainly involved in immune response, inflammatory response, and IFNγ. The significant enriched pathways including RA, inflammatory bowel disease (IBD), tuberculosis, etc (Figures 5A,B). The related genes of MG and AIT were mainly related to immune response (GO), antigen processing and presentation (GO), autoimmune thyroid disease (AITD) (KEGG), and allograft rejection (KEGG) (Figures 5C,D). Moreover, the related genes of AP and VKH were mainly associated with the antigen processing and presentation in GO and some pathways such as viral myocarditis, Staphylococcus aureus infection, AITD, intestinal immune network for IgA production, etc (Figures 5E,F).

To determine whether there was high genetic similarity among multiple ADs, we applied hierarchical clustering to the integrated similarity matrix of 26 ADs. The disease clusters consisting of AD pairs with integrated similarity scores greater than 0.3 were considered to be significant. As shown in Figure 6A, three significant disease groups were identified from the 26 ADs, and the ADs of each disease group might have high genetic similarity. We found that the three significant pairs of similar ADs were located in three different clusters, respectively. In the cluster one, bullous pemphigoid might be similar to AP and VKH. The ADs of cluster one are involved in endocrine autoimmunity. For instance, bullous pemphigoid has been proved to be related to immunodysregulation polyendocrinopathy enteropathy X-linked syndrome (Mcginness et al., 2006). In the cluster two, pemphigus, Graves disease, Sjogren’s syndrome, Addison Disease, and autoimmune hepatitis might be similar to MG and AIT. The ADs of cluster two contain the main AITDs (AIT and Graves disease) and frequent ADs involved in PolyA (AIT, Graves disease, and Sjogren’s syndrome) (Amador-Patarroyo et al., 2012; Botello et al., 2020). In the cluster three, T1D and MS might be similar to RA and SLE. The ADs of cluster three are all chronic inflammatory ADs and share multiple genetic susceptibility loci (Richard-Miceli and Criswell, 2012). Then, we performed pathway enrichment analysis of the related genes of ADs of three significant clusters. The related genes of ADs of three disease clusters were significantly enriched in pathways of T1D and IBD (Figures 6B–D). Therefore, we infer that T1D and IBD can participate in PolyA in various AD patients.

Next, we detected the causal genes in common among the ADs which belonged to the same significant disease cluster, and these genes could be used for PolyA research. As shown in Table 3, the ADs of cluster one shared one gene (HLA-DQA1); the ADs of cluster two shared two genes (HLA-DQB1 and HLA-DRB1); the ADs of cluster three shared eight genes (TNF, HLA-DRB1, PDCD1, PTPN22, CCR5, IL6, HLA-DQB1, and CTLA4). Identification of these genes will contribute to the discovery of novel prognostic, diagnostic, and therapeutic markers and justification of drug repurposing for ADs.