Взаимодействие полов и генотипов в ассоциациях по борьбе с заболеваниями: комплексный сетевой анализ в Биобанке Соединенного Королевства.

Background

Sexual dimorphism is an integral component to consider in our attempts to fulfill the promise of personalized medicine [ 1 ]. In addition to anatomy, physiology, and behavior, biological sex has been found to affect the prevalence, onset, and severity of nearly all human diseases [ 1 ]. Without a knowledge of the impact of sex on disease, the pursuit of transformed patient care through personalized disease prediction and treatment will remain incomplete.

Multiple factors contribute to such sex-specific effects on disease, including hormonal profiles, epigenetic and transcriptomic influences, differences in immune- and endocrine-related processes, and behavior [ 1 – 5 ]. Numerous studies have investigated the underlying genetic architecture of these sex-based differences—the interactions between genotype and sex on phenotypic outcomes that have been uncovered through these publications can be referred to as “genotype-by-sex” or “GxS” effects [ 1 , 6 – 8 ]. For instance, in 1998, O’Donnell et al. found a genetic association with hypertension biased toward men through a quantitative trait linkage analysis of sibling pairs in the Framingham Heart Study [ 9 ]. In 2009, Zhuang and Morris found genetic influences biased toward women through a genome-wide association study (GWAS) of rheumatoid arthritis [ 10 ]. Most recently, in 2021, Bernabeu et al. conducted a sex-specific phenome-wide association study (PheWAS) in the UK Biobank (UKBB) and found evidence of GxS effects for a variety of independent phenotypes including ankylosing spondylitis, gout, and hypothyroidism [ 11 ]. In particular, the large-scale analysis conducted by Bernabeu et al. stands out as an example of the impact that massive patient datasets can have in accelerating the identification of genetic associations across diseases. Nevertheless, a key issue remains across these publications – complex diseases rarely impact patients one-at-a-time.

Shared biological, environmental, and genetic factors can all contribute to onset of multiple phenotypes in a single patient [ 12 ]. These resulting cross-phenotype associations can impact patients with increased health burdens and risk of morbidity [ 12 ]. Thus, it is insufficient to evaluate the genetic architecture of individual diseases if we wish to gain a deeper understanding of overall patient health [ 12 , 13 ]. We must consider the genetics underlying the entire landscape of the human phenome and genome, known as the human “diseasome [ 14 ],” with a particular focus on the effects of polygenicity and genetic pleiotropy.

Given that GxS interactions are prominent in the pathogenesis of many individual complex diseases, we hypothesize that sex-specific differences also exist in the genetic associations across diseases. The field of network medicine offers an intuitive way of investigating the interactions between phenotypes [ 13 , 15 ]—in particular, the analysis of disease-disease networks (DDNs), graphs where nodes represent diseases and edges represent shared components between diseases, such as associated genetic variants, can reveal global and local connectivity across multiple phenotypes [ 14 , 16 – 20 ] as well as highlight the influence of pleiotropy [ 21 ] in disease onset. By analyzing a DDN, a researcher or clinician can evaluate how diseases are linked to one another, with immediate insight into potentially shared genetic architecture through the identification of pleiotropic single-nucleotide polymorphisms (SNPs) at specific genomic locations. In particular, a DDN generated from the results of a PheWAS applied to an electronic health record (EHR)-linked biobank is an ideal approach for investigating genetically-derived phenotypic interactions for the entire diseasome of a given population [ 22 ].

Patterns of genetic associations between diseases specific to different populations can be represented by their respective DDNs [ 19 ]. Thus, comparing two DDNs to one another can reveal disparities in cross-phenotype associations across populations. Several approaches exist in the field of network science to compare graphs to one another [ 23 ], including graphlet, graph kernels, and graph isomorphisms. With the addition of biological context for individual nodes, such methods can provide insight into differences in the connections between disease traits across networks. Furthermore, analysis performed by Elgart et al. in 2022 demonstrated the importance of evaluating differences in genetic correlations across gender from a network medicine comparison [ 24 ].

Thus, the objective of our study was to perform a biologically informed comparison of sex-specific DDNs to evaluate genetic interactions in cross-phenotype associations across sex. We construct a male- and female-specific DDN using sex-stratified PheWAS summary statistics from the UK Biobank (UKBB) [ 25 ] and compare these networks using network statistics, cluster analysis, edge set comparison, and node embeddings to identify potential GxS interactions. Our analysis identifies multiple examples of GxS effects in disease associations, revealing multiple phenotypes and SNPs of interest for future studies of the influence of sex on the genetic architecture of multimorbidities.

Methods

Data

We used sex-stratified PheWAS summary statistics from Bernabeu et al. [ 11 ] as the source dataset for the construction of our DDNs. These results were generated by applying DISSECT [ 26 ] for multiple phenotype GWAS using data from the July 2017 release of the UKBB. Corrections were made for age, array batch, UK Biobank Assessment Center, and the first 20 genomic principal components. Individuals whose self-reported sex did not correspond with their biological chromosomal content were excluded, as were individuals identified by the UKBB as outliers based on genotyping missingness rate or heterogeneity [ 27 ]. Furthermore, individuals with a first or second genetic principal component more than 5 standard deviations from the mean of self-reported white Europeans were excluded. These filtered data consisted of 452,264 British individuals of European ancestry, including 245,494 females and 206,770 males, corresponding to 54.28% and 45.72% of the full population respectively. Participants were genotyped with either the Affymetrix UK BiLEVE Axiom array or the Affymetrix UK Biobank Axiom array [ 28 ]. Genotypes were later augmented by imputation of roughly 90 million genetic variants from the Haplotype Reference Consortium, the 1000 Genomes project, and the UK10K project. Biallelic variants that did not pass the UK Biobank QC procedures, had a p-value less than 1e-50 for departure from Hardy–Weinberg equilibrium, and had a minor allele frequency (MAF) less than 1e-4 were excluded from sex-stratified PheWAS. Finally, a minor allele frequency threshold of 10% was applied to select for common variants in the data and mitigate effects of low sample size, yielding a final set of 4,229,346 autosomal (244,743 genotyped and 3,984,603 imputed) and 7,227 genotyped X-chromosome genetic variants [ 11 ]. Additional details regarding data generation and quality control can be found in the original publication from Bernabeu et al. [ 11 ].

The PheWAS performed by Bernabeu et al. used individual ICD-10 [ 29 ] encodings to identify phenotypes. GWASs were also run for “blocks” of ICD-10 codes, corresponding to ontological groupings of diseases [ 29 ]. We considered only these blocked groupings of diseases for our DDNs to ensure the interpretability of our analysis. Our final dataset for network generation included 104 phenotypes across 14 disease categories (Table  1 ).

Table 1 Counts of disease categories for the phenotypes used in analysis Disease Category Number of Phenotypes Circulatory System 9 Dermatologic 7 Digestive 10 Ear Disorders 4 Endocrine/Metabolic 6 Eye Disorders 10 Genitourinary 6 Hematopoietic 5 Infectious Diseases 6 Mental and Behavioral Disorders 5 Musculoskeletal 13 Neoplasms 7 Nervous System 8 Respiratory 8 Total 104

Counts of disease categories for the phenotypes used in analysis

Construction of DDNs

Curated sex-stratified PheWAS summary data for 104 binary traits were used to generate the male- and female-specific DDNs. The methodology described by Verma et al. [ 19 ] was applied to create each DDN. An edge in the set E = e ij V × | V | was established between each pair of binary phenotypes v i and v j if the two diseases shared associations with at least one common SNP at a genome-wide significance threshold of 1e-4. e ij represents the number of shared SNPs between two phenotypes. For all analysis that required edge weights to be constrained to values between 0 and 1, edge weight values were represented by the cosine similarities between the vectors of individual nodes in the adjacency matrix A of each network. The final DDN is an undirected, weighted graph G = ( V , E ) , where node set V represents the set of binary phenotypes and edge set E represents all connections between phenotypes. To generate, visualize, and analyze both graphs, we made use of Gephi 0.90 [ 30 ] and sigma.js [ 31 ], open-source network visualization software packages, as well as NETMAGE [ 18 ], a web-based tool that allows users to upload PheWAS summary statistics and generate corresponding interactive disease-disease networks. Further analysis and visualization of DDN network statistics were performed using R 4.4.1 [ 32 ] and Python 3.11.3 [ 33 ], including the following Python and R packages: pandas (v2.2.2), seaborn (v0.13.2), matplotlib (v3.9.3), numpy (v1.26.4), network (v3.4.2), scikit-learn (v1.5.2), community (v0.16), node2vec, ClustAssess (v0.3.0), patchwork (v.1.3.0), magrittr (v2.0.3), dplyr (v1.1.4), tidyr (v.1.3.1), tibble (v3.2.1), and ggplot2 (v3.5.1).

Network comparison

Our resulting male- and female-specific DDNs were compared to one another with respect to four attributes: network topology, clustering behavior, edge sets, and node embeddings (Fig.  1 ).

Overview of network construction and comparison. ( 1 ) Disease-variant associations are taken from PheWAS summary statistics to generate ( 2 ) corresponding sex-specific disease-disease networks. These networks are then ( 3 ) compared to one another through A network topology, B clustering behavior, C common edge sets, and D embedding of node features through graph representation learning

Network topology

Network analysis was performed on each graph using Gephi [ 30 ]. Calculated attributes included degree, weighted degree, network diameter, graph density, number of connected components, modularity, average clustering coefficient, and average path length. For a full description of the meaning of each network metric, refer to Supplemental Table  2 . These features provide an overview of the general structure of each network and serve as an effective baseline to evaluate differences in population-level phenotypic interactions. Furthermore, the degree, weighted degree, and betweenness centrality measurements were evaluated per node to identify the most central diseases in each network. For each of the three metrics, the top 10 nodes were identified. The resulting overlap of these node sets is presented in Table  3 , and the full centrality results are presented in Supplemental Table  3 .

Table 3 Most central diseases in sex-stratified DDNs based on degree, weighted degree, and betweenness centrality Female only Male only Both DDNs G35 to G37: Demyelinating diseases of the central nervous system I20 to I25: Ischemic heart diseases I10 to I15: Hypertensive diseases M15 to M19: Arthritis C43 to C44: Melanoma and other malignant neoplasms of skin J40 to J47: Chronic lower respiratory diseases E70 to E90: Metabolic disorders K90 to K93: Other diseases of the digestive system L40 to L45: Papulosquamous disorders E10 to E14: Diabetes mellitus E00 to E07: Disorders of thyroid gland

Most central diseases in sex-stratified DDNs based on degree, weighted degree, and betweenness centrality

Cluster behavior

An assessment of clustering behavior in DDNs can reveal patterns of community structure and provide evidence of cross-phenotypic interactions across disease categories. In particular, the Louvain method [ 34 ] is a popular approach for hierarchical clustering within graphs based upon modularity optimization. We selected this method over other clustering approaches due to its ability to optimize modularity within networks, its scalability for a large number of nodes and edges, its efficiency in performing hierarchical clustering of communities and subcommunities in a graph, and its flexibility for both weighted and unweighted graphs, as well as for directed and undirected graphs. [ 34 ] The Louvain method was applied to each graph to identify community cluster assignments for each disease. Cluster categorizations for each node were then compared across graphs to evaluate similarity in groupings. The level of correspondence between the clustering behavior for the two graphs was quantified using the element-centric similarity (ECS) measure, an unbiased approach for evaluating cluster resemblances [ 35 ]. Element-centric similarity was chosen over other metrics due to its ability to evaluate clustering quality at the level of individual nodes instead of the cluster level. Unlike measures such as Adjusted Rand Index or Normalized Mutual Information, ECS directly assesses whether relationships between elements are preserved rather than on global properties of clusters. It is also sensitive to over- and under-clustering, treats all cluster assignments equally, and provides a similarity score in the range of [0, 1], making it simpler to interpret.

Edge sets

A direct comparison of edges was performed across the networks to evaluate the extent of overlap in disease associations across the male- and female-specific data. Jaccard similarity was calculated to quantify the intersections across datasets. A higher Jaccard similarity closer to 1 indicates concordance in the identified genetic links between diseases across populations. Furthermore, edge weights for each graph were normalized to a value between 0 and 1 based on the cosine similarity metric, and correlations were calculated for the weights of the edges shared across the two graphs to determine how close to one another the distributions of edge weights were.

Node embeddings

Graph representation learning encodes the high-dimensional structure of a network into node-specific vectors, allowing for the identification of node-level contributions to disease associations and the comparison of nodes to one another through their generated features. Graph representation learning was applied using node2vec [ 36 ] to evaluate the contributions of nodes across graphs. Node2vec was chosen to perform node embeddings due to its ease of use and computational efficiency, its ability to balance both homophily and structural equivalence in evaluating node characteristics, and its use of biased random walks to conduct more nuanced exploration of the graph compared to other methods. Node2vec takes as input the structure of the network represented as an adjacency matrix and returns a vector of values for each node in the graph, corresponding to each individual disease’s embedding. Its application allows for the characterization of each node by the role it plays in both its local neighborhood and the entire diseasome. The return parameter (p) and the in–out parameter (q) of node2vec were both set to 1 to ensure equal likelihood in exploring local and global structure of the networks, akin to a DeepWalk [ 37 ] approach. Based upon the structures of both graphs and in accordance with other default parameters for node2vec, the dimension of the embedding representation was set to 16, the walk length set to 20, the number of walks set to 5000, the window size set to 10, and the number of epochs set to 10. Alternate values for these parameters were explored but were not found to significantly affect embedding comparison analysis. The embedding results from node2vec are unique to each graph, and thus do not provide insight into cross-network comparisons. However, by evaluating embedding distances between pairs of nodes within a single graph, we can gain a deeper understanding of each node’s influence on the topology and structure of the overall network. Thus, based on the embedding output from node2vec, we calculated the average distance from each node to all other nodes in its network, offering us a view of its centrality in the diseasome. These average distances were then compared per node across networks to determine how similar each disease’s connections to other phenotypes was across the sexes.

Statistical evaluation against distributions of null network models

As a method of determining how statistically significant the differences were between metrics for the female and male DDNs, 1000 random networks were generated using the “networkx” package in Python according to the estimated degree distributions of each of the original DDNs. Unique pairings of each of these 2000 networks (resulting in 2000*1999/2 = 1,999,000 comparisons) were then generated to produce distributions of expected differences for each of the metrics being evaluated. One-sample t-tests were conducted to produce p-values indicating the statistical significance of the difference between the female and male-specific DDN compared to the expected difference for each metric.

For the comparison of edges using Jaccard similarity, a separate random network generation method was used, focused purely on the generation of edges—we determined all possible edges from the male and female DDNs by looking at all combinations of nodes, and then for each DDN, we determined the proportion of edges present in the full set of possible edges. These proportions dictated the graph density for each network. We then generated 2000 new edge sets—1000 of these edge sets used the graph density from the female DDN, while the other 1000 used the graph density from the male DDN. For each edge in each edge set, based on the currently chosen graph density, we randomly generated a number between 0 and 1 to determine if the edge was present or absent. This operation generated 2000 different edge sets, each with different numbers of edges that closely reflected the graph density of the original male and female DDNs. Unique pairings of each of these 2000 networks, again resulting in 1,999,000 comparisons across edge sets, were then evaluated for Jaccard similarity to develop a distribution of expected Jaccard indices for pairs of randomly generated networks. Lastly, a one-sample t-test was conducted to identify the statistical significance of the Jaccard index between the female- and male-specific DDN compared to this expected distribution.

Results

Network topology

We generated our female- and male-stratified DDN from the source UKBB PheWAS summary statistics (Fig.  2 ). Both networks included 103 diseases in their main connected components. The full set of phenotypes that overlapped between the two graphs consisted of 102 common diseases (Supplemental Table  1 ). The female DDN included 676 edges, while the male DDN included 598 edges. Table 2 includes a summary of network statistics for both graphs. See Supplemental Table  2 for a description of each network statistic. The general structure and edge compositions of these networks can also be explored through an interactive network visualization tool at https://hdpm.biomedinfolab.com/ddn/ukbb-female and. https://hdpm.biomedinfolab.com/ddn/ukbb-male [ 18 ].

A visualization of the A female- and B male-specific DDNs. Nodes are colored by disease category and sized by degree, while edge widths are determined by their weight, corresponding to the number of SNPs shared between each pair of diseases. The zoomed-in subgraphs of both networks include the top 10 nodes from each graph according to PageRank

Table 2 Network statistics for the two DDNs Female value Male value Difference Expected difference (null model) Standard deviation (null model) P -value Average degree 13.13 11.61 1.51 0.637 0.890 < 1e−4 Average weighted degree 1391.36 904.45 486.91 226.82 234.29 < 1e−4 Average path length 2.161 2.35 − 1.89e−01 − 4.64e−02 6.69e−02 < 1e−4 Network diameter 4 6 − 2 − 1.96e−01 7.82e−01 < 1e−4 Graph density 0.129 0.114 1.48e−02 6.24e−03 7.94e−03 < 1e−4 Modularity (Louvain method) 0.131 0.239 − 1.08e−01 4.99e−02 6.41e−02 < 1e−4 Average clustering coefficient 0.481 0.398 8.33e−02 2.26e−02 4.45e−02 < 1e−4

Network statistics for the two DDNs

A direct comparison of these values indicates general similarity in the overall structure of the male- and female-specific DDNs, with relatively close diameters, densities, and degrees for nodes. However, after comparisons to expected distributions for the calculated metrics, we find that the differences across all metrics are significantly larger in magnitude than would be expected at random, suggesting the presence of a notable GxS effect on the presence of disease multimorbidities. Through a visual comparison of the two graphs, a stronger interconnected module is evident at the center of the male DDN, represented by thicker edges in the center and thinner edges along the periphery. This behavior in the graph indicates higher levels of genetic association between phenotypes for this subset of diseases. Conversely, the female-specific DDN exhibits a much higher average weighted degree across nodes (1377.981 compared to 895.75), represented by more consistent edge thickness throughout the entirety of the network. This imbalance reflects the higher number of disease-SNP associations within the original female PheWAS data.

Central diseases

Within each DDN, central nodes represent diseases with a high level of genetic associations to other phenotypes. Using three different definitions of network centrality, degree, weighted degree, and betweenness centrality, we identified the most central diseases in both graphs (Table  3 , Supplemental Table  3 ). These nodes correspond to phenotypes that exhibit the highest amount of genetic pleiotropy with other diseases in each network.

Common central diseases for both graphs included hypertensive disorders, bronchitis and asthma, diabetes mellitus, thyroid disorders, and celiac diseases. Central diseases in only the female-specific DDN included demyelinating diseases of the nervous system such as multiple sclerosis and arthritis. Central diseases in only the male-specific DDN included ischemic heart diseases such as anginas and myocardial infarction, skin cancers, metabolic diseases, and papulosquamous diseases such as psoriasis.

Clustering behavior

The Louvain method segregated both the female- and male-specific DDNs into four distinct modules. However, the element-centric similarity score between these two assignments came out to only 0.304, indicating close to random agreement in the clusters assigned to each disease across the two networks. Examining a breakdown of disease categories (Fig.  3 ) in the clusters identified by Louvain clustering, we see a general level of agreement in the largest, most central cluster. Cluster ID 1 includes diseases from the circulatory, digestive, eye disorder, and musculoskeletal categories (Supplemental Table  4 )—however, individual diseases and the category breakdowns for the remaining clusters differ more broadly. Thus, it appears that this core set of phenotypes shares strong genetic associations with one another in both sexes, and that the clustering behavior for the remaining phenotypes outside of this cluster varies more significantly across the two graphs. This variability reflects a potential “core-periphery” structure in both networks [ 38 ], with a core sub-graph of highly interconnected diseases shared across the sexes and varied peripheral sets of nodes loosely connected to each core.

Bar charts of cluster assignments after Louvain clustering for both DDNs, colored by disease category

Edge sets

Edge sets for both DDNs were compared to one another to evaluate consistency of cross-phenotype associations across both populations. Out of the 676 edges in the female-specific DDN and the 598 edges in the male-specific DDN, the two graphs shared 266 edges in common (Supplemental Fig.  1 ).

The union of these two edge sets, representing the full list of unique cross-phenotype associations identified in the creation of our DDNs, was 1008, yielding a Jaccard similarity coefficient of 26.39%. To account for the scale-free nature of the network and reduce the impact of weaker weights on the calculation of Jaccard similarity, we applied an edge-cut threshold [ 39 ], filtering out edges from each graph that had weights lower than the median edge weight of the original networks (5 for the female-specific DDN and 8 for the male-specific DDN). After dropping edges from both networks with weights lower than these values, we found that the Jaccard similarity of our two graphs was 30.21%. To evaluate the significance of this overlap, we developed a null model distribution of randomly-generated edges for each graph to produce simulated edge weights. Table 4 presents the results of this analysis.

Table 4 Jaccard similarity of edge sets across the two DDNs Original graphs Null model mean Null model standard deviation Null model p -value Jaccard similarity 30.21% 42.54% 1.96% < 1e−4

Jaccard similarity of edge sets across the two DDNs

Based on our results, we can see regardless of the model that our edge sets differ significantly from similarities expected from a null model. Our Jaccard similarity indicates that the genetic associations identified for the two sexes exhibit a low level of similarity compared to the expected agreement of edge sets.

To evaluate the edges common across our networks, we calculated the correlations between edge weights after normalization to a range of [0,1] based on cosine similarity. The resulting Pearson correlation coefficient (r = 0.586) indicates that even though a low number of edges are shared between the two graphs, there is a moderate similarity in the weights of the edges that are shared.

Within the context of disease categories, we find that seemingly similar proportions of cross-category phenotype associations exist across the two networks, with 89.35% of associations connecting diseases of different categories for females and 88.13% of associations connecting diseases of different categories for males. However, comparing the difference in proportions to the average difference in percentage of cross-category edges between pairs of networks from our randomly-generated null distribution, we see that our two original networks exhibit statistically different proportions of cross-category edges from one another (Table  5 ).

Table 5 Cross-category edge proportions Female value Male value Difference Expected difference (null model) Standard deviation (null model) P -value Percent of cross-category edges 89.35% 88.13% 1.22% − 0.17% 1.20% < 1e−4

Cross-category edge proportions

Different patterns of disease interactions across the sexes become more apparent when considering the number of edges that connect nodes across individual disease categories (Fig.  4 ). For the female-specific DDN, numerous interactions exist across the musculoskeletal, digestive, and endocrine/metabolic disease categories, reflecting strong genetic associations present in immune disorders such as arthritis and endocrine-related phenotypes such as hyper- and hypothyroidism [ 40 ]. For the male-specific DDN, associations between diseases are primarily concentrated within circulatory system phenotypes. This discrepancy may reflect known male-specific genetic associations with cardiovascular traits [ 41 , 42 ].

Heatmaps of edge sets across disease categories for female- and male-specific DDNs. Darker colors indicate more edges shared between disease categories

Looking closer at the edges that are shared between the two networks, we see a concentration of connections between the endocrine/metabolic, circulatory, and musculoskeletal disease categories (Fig.  5 ), in keeping with the central roles we saw such phenotypes play across both sexes in the results of our centrality analysis.

Heatmap of edge sets across disease categories shared by both the female- and male-specific DDNs. Darker colors indicate more edges shared between disease categories

Node embeddings

Graph representation learning using node2vec was applied to the nodes of each graph to represent disease associations within distinct embedding spaces. To identify diseases that were centrally located within clusters of other phenotypes, we pinpointed nodes with the shortest average embedding distances from other nodes. While most phenotypes exhibited no sort of correlation in embedding centrality across the two graphs, we found that eye disorders occupied the central position in clusters for both sexes. These phenotypes included visual disturbances and blindness (H53-H54) for both sexes, disorders of conjunctiva (H15-H22) for males, and disorders of lens (H25-H28) for females. This consistent pattern suggests that eye disorders tend to cluster together across both sexes. In visualizing the distances between centrally embedded diseases common to both networks (Fig.  6 ), we found that, on average, the embedding distances are more condensed within the female-specific DDN than the male-specific DDN.

Heatmaps of embedding distances for most central diseases for female- and male-specific DDNs

These shorter distances in the female-specific DDN align with the higher density of the network and the increased number of disease-variant associations in the original female PheWAS data.

Effect of sexually-dimorphic SNPs on cross-phenotype associations with thyroid disorders

Within the original source data paper, Bernabeu et al. identified two sexually-dimorphic SNPs (sdSNPs) which were found to have significantly different statistical associations with E00-E07, disorders of the thyroid gland: rs9357120 and rs3130552. We evaluated the effect of these sdSNPs on sex-stratified DDN topology to identify additional instances of GxS interactions in cross-phenotype associations. rs9357120 is located on chromosome 6 at position 31,262,940. Despite its association with E00-E07 for both strata of the PheWAS data, based upon the p-value thresholding required for inclusion in our network analysis, rs9357120 was not included in the male-specific DDN. Within the female-specific DDN, this SNP contributes to an association between E00-E07 and G35-G37, corresponding to demyelinating diseases of the central nervous system. This association reflects the significance of autoimmune disorders for women [ 40 , 43 ], linking both thyroid disorders and phenotypes including multiple sclerosis. rs3130552 contributes to links between E00-E07 and K90-K93 (Other diseases of the digestive system) for both sexes. Furthermore, this SNP generates an association with L40-L45 (Papulosquamous disorders) only for the female-specific DDN, as well as a link with E10-E14 (Diabetes mellitus) only for the male-specific DDN. These disparate associations reflect the relevance for metabolic disorders for males and immune disorders for females [ 40 , 43 ]. We also scraped the Genotype-Tissue Expression (GTEx) Portal to identify relevant expression quantitative trait loci (eQTLs) [ 44 ] associated with these SNPs. Both SNPs map to multiple genes, including HLA-C , PSORS1C2 , PSORS1C3 , POU5F1 , HCG22 , HCG27 , and CCHCR1 . In particular, CCHCR1 , corresponding to the transcript ENSG00000204536.13, is a sex-biased eQTL expressed in thyroid tissue, with a p-value of 1.5e-14 and a normalized effect size of − 0.31 (Fig.  7 ). This expression pattern suggests a potential genetic contributor to the sex-specific cross-phenotype associations involving E00-E07.

Bulk tissue gene expression for CCHCR1 from the GTEx Portal. Distributions are sorted by median log value of gene expression. Thyroid tissue exhibits the highest gene expression for a sex-agnostic tissue

Discussion

Given that individual diseases are often affected by GxS interactions, we hypothesized that cross-phenotype associations may also be influenced by GxS effects. Sex-stratified GWAS and PheWAS offer a significant opportunity to identify sex-stratified genetic associations with disease onset. The identification of these possible instances of genetic pleiotropy also remains unbiased because PheWASs are both disease- and variant-agnostic [ 21 , 22 , 45 ]. Furthermore, the results of a PheWAS can lend themselves to the network-based analysis of associations between phenotypes [ 19 ].

In this study, we generated and compared DDNs of genetic associations between binary phenotypes using significant SNPs from sex-stratified PheWAS summary data. The similar structures of our male- and female-specific DDNs suggest consistency of the behavior of disease interactions across the sexes. However, the significantly lower Jaccard similarity score of the networks’ edge sets relative to a null distribution as well as the inconsistent clustering behavior of diseases across these two networks indicate the presence of GxS interactions in the onset of disease multimorbidities.

Furthermore, although the proportion of cross-category links is quite consistent across the two networks, an analysis of the specific edges reveals a concentration of edges between circulatory diseases for males, suggesting the highest level of genetic association within this phenotype category. For the female-specific DDN, a wider spread of connections exists across musculoskeletal, digestive, and endocrine/metabolic diseases, indicating a broader swath of cross-phenotype associations for females. Thus, sex seems to be significantly associated with the underlying genetic architecture of cross-phenotype associations.

Our DDNs also reflect a stronger overall spread of disease-SNP associations for females compared to males. This is evinced through the larger number of edges in the female-specific DDN, a higher average weighted degree, and shorter average embedding distances after applying node2vec. This difference across the sexes may be a result of a higher sample size for females [ 46 ], or the result of factors outside of genetics playing a stronger role in disease onset for males [ 1 ].

Considering the placement of individual diseases within each DDN allows us to evaluate the prominence of their genetic associations with other phenotypes. The more central a node is in the DDN, the more influential its genetic underpinnings are with respect to its cross-phenotype associations. A phenotype that is more central in one DDN than the other indicates genotype-by-sex interactions that influence its onset and its comorbidities. For both of our sex-stratified DDNs, cardiometabolic diseases and lower respiratory diseases display high centrality. Our DDNs also reflect known prominent roles for metabolic disorders in males and autoimmune disorders in females [ 40 ]. More specifically, our male-specific DDN highlights the following diseases as central in the network: I20-I25 (Ischemic heart diseases), C43-C44 (Melanoma and other malignant neoplasms of skin), L40-L45 (Papulosquamous disorders), and E70-E90 (Metabolic disorders). Meanwhile, our female-specific DDN highlights G35-G37 (Demyelinating diseases of the central nervous system) and M15-M19 (Osteoarthritis) as central. These central diseases are all confirmed in the literature as exhibiting sex-specific prevalence and onset [ 42 , 47 ]. Our DDNs suggest that genetic components may factor not only in their sex-specific onset but also their interactions with other phenotypes.

We identified centrally embedded diseases for both networks using node2vec. Our embeddings reveal a high level of clustering among eye disorders including visual disturbances and blindness, disorders of conjunctiva, and disorders of lens, suggesting that these phenotypes exhibit lower levels of genetic association with other types of diseases.

Investigating sdSNPs that have been found to have differential associations with individual phenotypes across sexes [ 11 ] provides evidence that these SNPs also have differential correlations with cross-phenotype associations. We highlight two sdSNPs, rs9357120 and rs3130552, associated with disorders of the thyroid gland, and find their differential associations with metabolic disorders for males and autoimmune disorders for females. Furthermore, perusal of the GTEx Portal demonstrates that both of these SNPs are associated with a corresponding sexually-dimorphic eQTL: CCHCR1 . This result provides evidence of how a genetic component may differentially affect the cross-phenotype associations observed across the sexes.

Our comparison of sex-specific DDNs serves as a landmark foray into the study of genotype-by-sex interactions in cross-phenotype associations. In particular, our network-based representation of disease interactions offers an intuitive, methodical way of characterizing the role of individual diseases across the human diseasome. Identifying individual genetic associations involved in these disease links offers further opportunities to pinpoint genetic contributors to disease comorbidities and interactions.

There are a few limitations to consider in our study. First, factors external to genetics may also interact with sex – hormonal profiles, including boosted immune support from estrogen, epigenetic and transcriptomic influences, endocrine profiles, and the onset of menopause can all lead to sex-specific differences in disease onset [ 3 , 5 ]. Indeed, the consideration of age as a variable is crucial to improve our understanding of the onset of pathological phenotypes – future work should compare age-stratified DDNs to improve our ability to infer disease risk profiles. The fulfillment of gender-based roles and societal biases related to gender are also likely to influence the source data with which we work [ 3 , 4 ]. For instance, diagnoses from clinicians are prone to bias, such as overdiagnoses of mental disorders and underdiagnoses of pain-related phenotypes for women [ 3 , 4 ]. Furthermore, men have been found to be riskier in their daily activities and less likely to seek professional care when needed [ 3 , 4 ]. In terms of the specific phenotyping definitions within the source data, the ICD-10 coding system is known to be an imperfect method of capturing true occurrences of diseases in patient datasets [ 48 ]. Our use of the “blocked” ICD-10 phenotype groupings was applied to improve our definitions of diseases but is still not necessarily the most accurate representation of patient health. Similarly, our binary representations of disease are limited in the extent of information that can be represented regarding patient diseases compared to continuous phenotypes [ 49 ]. With respect to data quality control and processing, the minor allele frequency threshold of 10% in the source data may lead to missed phenotype-SNP associations. However, as clarified through response papers to Bernabeu et al., these genetic associations are still represented through other SNPs in LD with true causative or correlative SNPs [ 50 , 51 ]. Further, we used a strict p-value threshold of 1e-4 to establish a SNP as being significantly associated with a phenotype in our DDNs [ 17 ]. While this stringency may lead to our missing some genetic associations between phenotypes in our networks, we can also be more confident in the connections that are ultimately included [ 17 ]. We also note that our DDNs represent data only for the UKBB population, meaning that conclusions drawn from our analysis can only be interpreted from a British European perspective [ 25 ]. Thus, any conclusions drawn from our analysis need to bear these potential inaccuracies and biases in mind. Most significantly, it is unclear how much of an effect the difference in sample sizes between the male and female populations in the source UKBB dataset have affected the genetic associations to which we have access. Further research must be conducted externally to determine the extent to which sample size affects results in genome-wide association studies.

Future work that extends our analysis includes the comparison of sex-specific DDNs derived from different baseline populations such as the All of Us dataset [ 52 ] or the Penn Medicine Biobank [ 53 ], the application of different graph representation learning methods [ 54 ] and the comparison [ 55 ] of these embeddings, an incorporation of mendelian randomization [ 56 ] to identify potential causative genetic relationships, and a further investigation of individual SNPs that may contribute to sex-specific associations between diseases.

Overall, our sex-stratified network comparison confirms the presence of GxS interactions in cross-phenotype associations and helps to navigate the study of complex diseases through the lens of genetic pleiotropy and polygenicity. The results of our analysis can facilitate future explorations of precision medicine, providing further insight into sex-specific genetic targets for drug development [ 57 ] and improving our ability to conduct disease risk prediction [ 1 ].

Conclusions

Through the construction and comparison of disease-disease networks from sex-stratified PheWAS summary data, we find evidence of numerous genotype-by-sex interactions that influence cross-phenotype associations across the human diseasome. We highlight disease categories as well as individual phenotypes that exhibit evidence of GxS influence in their associations with other phenotypes. Our analysis confirms several known sex-biased traits and identifies new associations as well, suggesting genetic contributors to their connections with other phenotypes. Ultimately, our work provides a view into the interplay between sex and genetics on disease associations and comorbidities across the landscape of human disease.

Supplementary Information

Supplementary Material 1.

Acknowledgements

The Genotype-Tissue Expression (GTEx) Project was supported by the Common Fund of the Office of the Director of the National Institutes of Health, and by NCI, NHGRI, NHLBI, NIDA, NIMH, and NINDS. The data used for the analyses described in this manuscript were obtained from the GTEx Portal on 07/03/23. We would like to thank Manu Shivakumar for his help in developing the web-based visualizations of the female- and male-specific networks. All figures were created or edited using BioRender.com.

Author contributions

All authors were involved in designing and conceptualizing the study. YA and DK supervised the study. VS performed data curation and analysis. VS wrote the original manuscript. All authors revised and approved the final manuscript.

Funding

This work has been supported by the National Institute of General Medical Sciences (NIGMS) R01 GM138597 and S10OD023495.

Availability of data and materials

Summary statistics for binary disease phenotypes can be accessed at https://datashare.ed.ac.uk/handle/10283/3908 . Additionally, sex-specific disease connections can be explored on our web visualization tool at https://hdpm.biomedinfolab.com/ddn/ukbb-female and https://hdpm.biomedinfolab.com/ddn/ukbb-male . Lastly, all analysis materials can be found at https://github.com/vsriram24/sexSpecificNetworkComparison .

Declarations

Ethics approval and consent to participate

The UK Biobank project was approved by the National Research Ethics Service Committee North West-Haydock (REC reference: 11/NW/0382) and all participants agreed to their inclusion.

Consent for publication

All the authors agreed to publish this article.

Competing interests

The authors declare no competing interests.

Footnotes

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Contributor Information

Yong-Yeol Ahn, Email: yyahn@iu.edu.

Dokyoon Kim, Email: dokyoon.kim@pennmedicine.upenn.edu.

References

Associated Data

Supplementary Materials

Supplementary Material 1.

Data Availability Statement

Summary statistics for binary disease phenotypes can be accessed at https://datashare.ed.ac.uk/handle/10283/3908 . Additionally, sex-specific disease connections can be explored on our web visualization tool at https://hdpm.biomedinfolab.com/ddn/ukbb-female and https://hdpm.biomedinfolab.com/ddn/ukbb-male . Lastly, all analysis materials can be found at https://github.com/vsriram24/sexSpecificNetworkComparison .