Multi-PheWAS intersection approach to identify sex differences across comorbidities in 59 140 pediatric patients with autism spectrum disorder

Abstract Objective To identify differences related to sex and define autism spectrum disorder (ASD) comorbidities female-enriched through a comprehensive multi-PheWAS intersection approach on big, real-world data. Although sex difference is a consistent and recognized feature of ASD, additional clinical correlates could help to identify potential disease subgroups, based on sex and age. Materials and Methods We performed a systematic comorbidity analysis on 1860 groups of comorbidities exploring all spectrum of known disease, in 59 140 individuals (11 440 females) with ASD from 4 age groups. We explored ASD sex differences in 2 independent real-world datasets, across all potential comorbidities by comparing (1) females with ASD vs males with ASD and (2) females with ASD vs females without ASD. Results We identified 27 different comorbidities that appeared significantly more frequently in females with ASD. The comorbidities were mostly neurological (eg, epilepsy, odds ratio [OR] > 1.8, 3-18 years of age), congenital (eg, chromosomal anomalies, OR > 2, 3-18 years of age), and mental disorders (eg, intellectual disability, OR > 1.7, 6-18 years of age). Novel comorbidities included endocrine metabolic diseases (eg, failure to thrive, OR = 2.5, ages 0-2), digestive disorders (gastroesophageal reflux disease: OR = 1.7, 6-11 years of age; and constipation: OR > 1.6, 3-11 years of age), and sense organs (strabismus: OR > 1.8, 3-18 years of age). Discussion A multi-PheWAS intersection approach on real-world data as presented in this study uniquely contributes to the growing body of research regarding sex-based comorbidity analysis in ASD population. Conclusions Our findings provide insights into female-enriched ASD comorbidities that are potentially important in diagnosis, as well as the identification of distinct comorbidity patterns influencing anticipatory treatment or referrals. The code is publicly available (https://github.com/hms-dbmi/sexDifferenceInASD).


INTRODUCTION
One of the most consistent features reported for autism spectrum disorders (ASDs) has been the imbalance of sex ratio. [1][2][3] The current sex ratio is estimated to be 3:1 male-to-female, 4,5 varying according to the age and population under study. The reasons for this disproportionate sex imbalance, however, remain unclear. Various biological theories have been proposed, including those related to genes, hormones, and brain structure. 2,6 The female protective effect theory 7,8 suggests that females require a higher number of genetic mutations to develop ASD, but then present greater ASD severity. 9,10 By contrast, males with ASD present a higher genetic variability and decreased ASD severity. 11 Other studies have identified sex chromosomes as potential risk factors for neurodevelopmental disorders, such as ASD. Some have argued that possession of a second X chromosome might be protective for females. 12 An additional hypothesis postulates that sexually dimorphic hormonal exposures, 13,14 such as high testosterone levels during pregnancy, could account for the preponderance of ASD in males. 15 While all these hypotheses point to potential mechanisms associated with ASD sex differences, biology is not the only potential explanation for the sexual dimorphism. The development of instruments to detect ASD relied mainly on male samples, 7,8 contributing to a possible misdiagnosis of females, late diagnosis when compared with males, or a lack of diagnosis overall. [9][10][11][12] The working hypothesis is that some key features for identifying males with ASD are not as prominent in females, [13][14][15][16][17] and these cues are often overlooked during diagnosis. Therefore, applying the same ASD diagnostic criteria to both sexes 12,18,19 might have inflated the maleto-female sex ratio.
To study these and other causes further, researchers have conducted analyses of comorbidity patterns. One widely used approach is the phenome-wide association analysis (PheWAS). This approach consists of comparing 2 groups, cases and controls, to identify statistically significant phenotypes more likely to be present in cases than in controls. Cases and controls are defined by a unique variable (eg, having or not a disease). These are useful and arguably critical to unlocking some of the complexities of the disorder, as much of the diversity seen in ASD relates to co-occurring conditions. In addition, the phenotypes that co-occur with ASD may help distinguish distinct ASD subgroups, which, upon thorough investigation, could enhance understanding of the disease. Evidence shows that the comorbidity burden, defined as the coexistence of 2 or more diseases in the same person, 20 is higher in individuals diagnosed with ASD, as compared with the general population, 21,22 including gastrointestinal conditions, [23][24][25] epilepsy, 30 and psychiatric comorbidities. 21,25,32 Yet, very few analyses have examined whether these comorbidity patterns in patients with ASD are distinct for each sex [20][21][22][23][24] and they have been focused on a predefined list of comorbidities. In addition, there is a dearth of comorbidity research in some age groups, including infants with ASD. [20][21][22][23][24] Further, investigators have been hindered by insufficient sample size. 13,20,21,25 Larger, diverse, and more representative samples are required to draw meaningful conclusions about differences between ASD subgroups. We sought to obtain such a sample through the use of realworld data, defined as data relating to patient health status and/or the delivery of health care routinely collected from a variety of sources. We drew on 2 such independent real-world datasets to yield a sizable enough population of individuals with ASD to perform studies that would yield clinically meaningful results.
The primary objective of this study is to develop an approach, the multi-PheWAS intersection, that enable investigators to explore sex differences in real-world datasets, across all diagnostic codes and potential comorbidities. We apply this method to identify female ASD characteristic phenotypes by comparing (1) females with ASD vs males with ASD and (2) females with ASD vs females without ASD. We use another dataset collected from a pediatric university hospital to further validate those results .
The biomedical informatics contributions of this study include (1) the design, implementation of a multi-PheWAS intersection pipeline; (2) reducing the potential bias in the identification of mental disorder patients when using real-world data by increasing the specificity using a threshold of 3 diagnostic codes per patient; and (3) making the full code available in GitHub to enable the reproduction of these results to anyone with access to the datasets analyzed in this study and that can also be easily adapted to other datasets and diseases changing the value of the input parameters.

Data
We analyzed 2 independent real-world datasets: Aetna, 26 a nationwide claims database of 67 million patients, collected between January 1, 2008, and December 1, 2019, and an independent set of data from 1.7 million patients at Boston Children's Hospital, collected between January 1, 2008, and September 19, 2016.
To be eligible for this study ( Figure 1A), individuals 0 to 18 years of age had to be diagnosed upon at least 3 different occasions, based on an International Classification of Diseases-Ninth Revision-Clinical Modification (ICD-9-CM) 27 code for ASD (299.00, 299.01, 299.80, 299.81, 299.90, 299.91). Eligible children also had at least 1 phenotypic category in addition to ASD and their records collected for more than 12 months ( Figure 1B). As a control, we studied all 7 million children without ASD (3 763 418 females and 3 858 927 males) presented in the claims database. To avoid overlap between the 2 datasets, patients from Boston Children's Hospital were not eligible for this study if their insurance was the same as the claim database.

Phenotype comorbidity analysis
To identify statistically significant phenotypic codes that were more likely present in females with ASD, we deployed a multi-PheWAS intersection approach, analyzing the intersection of phenotypes that are significantly different in distinct phenome-wide comparisons ( Figure 1C). We first selected and classified 2 groups of patients: cases (females with ASD) vs controls (males with ASD and females without ASD) (Figure 2A). We then systematically examined each phenotypic category to find associations between either sex (females with ASD vs males with ASD) or disease (females with ASD vs females without ASD) ( Figure 2B). To identify statistically significant phenotypic codes that were more likely present in females without ASD, we deployed the same phenotype comorbidity approach, comparing females without ASD with males without ASD ( Figure  2C). This allows us to identify comorbidities more likely in females in the general population. We then excluded these phenotypic codes from our results because they are not associated with the fact of having ASD.
Diagnoses were initially presented as ICD-9-CM codes (10 279 in the claims and 8585 in the hospital datasets). We aggregated them using PheWAS 28 codes (1860 in the claims and 1709 in hospital datasets). We then determined the counts for each patient and phenotypic category. To increase the reliability of the reuse of diag- Figure 1. General workflow for the phenotype comorbidity analysis. For our method, we began with our selection of patient data. (A) Autism spectrum disorder (ASD) patient inclusion criteria: 0-18 years of age, assigned an ASD International Classification of Diseases-Ninth Revision-Clinical Modification (ICD-9-CM) code at least 3 times, and records collected for at least 12 months. (B) Patients without ASD inclusion criteria: 0-18 years of age, never assigned an ASD ICD-9-CM code, and records collected for at least 12 months. We then proceeded through our (C) phenotype comorbidity analysis. For each of the phenome-wide association analysis (PheWAS) codes present per set, we ran a Fisher exact test, comparing patients with or without ASD, males vs females. nostic codes as a proxy for clinical diagnoses, we restricted our cohort to patients that received ASD diagnostic codes at least 3 distinct times. 28 The same approach was used for the comorbid diagnosis, only those recorded on at least 3 distinct occasions were included in our analysis.
We analyzed all phenotypic codes presented by patients as follows: phenotypic codes diagnosed in patients 0 to 2 years of age (infants), 3 to 5 years of age (prepubertal), 6 to 11 years of age (peripubertal), and 12 to 18 years of age (postpubertal). 29 Each patient can contribute across several age groups with different phenotypes during their development. For each phenotypic code and age group, we ran a Fisher exact test, based on diagnosis or lack of diagnosis with ASD, and separately, sex. We estimated the confidence interval (CI), odds ratio (OR), and P value. The conservative Bonferroni correction was applied to correct for multiple testing. To be considered more likely to occur in females with ASD, a phenotypic code had to be statistically significant with a corrected P value <.01/total. Phenotype codes at each age group (detailed information in Supplementary Material 1), as measured by comorbidity testing, and have an OR >1.5 (meaning 1.5Â more likely to be present in one group as compared with other groups). This cutoff of at least 1.5 was applied to get closer to results that could potentially be clinically significant. We explored how many statistically significant results were obtained from multiple OR cutoffs from 1.1 to 2.0. For visualization purposes, the heatmaps have been organized based on the phenotype category (eg, "Epilepsy," "Partial epilepsy," and "Convulsions" are plotted in consecutive rows as they all belong to the category "Neurological"). These higher phenotype categories have been extracted from the phecode definition file available at the PheWAS catalog site (https://phewascatalog.org/phecodes_v1_1). The analysis was performed using SQL and R version 3.4.1 (R Foundation for Statistical Computing, Vienna, Austria). Open-source code is publicly available at GitHub (https://github.com/hms-dbmi/sexDifferen-ceInASD).
The study was approved by the central Institutional Review Board at the National Human Genome Research Institute (IRB-P00026337).

RESULTS
After the inclusion criteria, using claims data, we analyzed 59 140 individuals with ASD (11 440 females and 47 700 males) that met our inclusion criteria. We further validated our findings using an independent pediatric hospital dataset, composed of 3728 patients with ASD (762 females and 2966 males) and 28 740 randomly selected individuals without ASD (13 988 females and 14 752 males) after applying the inclusion criteria.

Identification of phenotypes in the claims dataset
We discovered 27 distinct phenotype codes that were more likely present in females with ASD, as compared with males with ASD and females without ASD (Table 1, Figure 3, and Supplementary Tables S1 and S2, which contains the raw numbers of females and males with and without ASD for each specific phenotype category). Demonstrating specificity to ASD, we showed that these phenotype codes were not more likely present in females without ASD, as compared with males without ASD (Supplementary Table S3). Supplementary  Table S3 does not contain all the phenotypes that were more likely in females without ASD vs males without ASD, it only contains those phenotypes that were also significantly more prevalent in females with ASD in the other 2 compared groups.
We found 3 phenotypic codes shared in the infant group, 9 in prepubertal, 17 in peripubertal, and 18 in adolescence ( Table 1). The 3 phenotypic infant codes were other specified congenital anomalies of the nervous system, failure to thrive, and symptoms related to nutrition, metabolism, and development ( For individuals diagnosed in peripubertal, we discovered 17 phenotypic codes that more likely appeared in females with ASD. As with prepubertal, most phenotypic codes in the peripubertal group were neurological and congenital. Digestive issues (gastroesophageal reflux disease and constipation) as well as sense organs (strabismus) were also found in this age group as more likely to occur in females with ASD. When comparing females and males with ASD in the same age group, the phenotypic code with the highest OR of 20.1 (95% CI, 8.734-53.75) was "other cerebral degenerations." By contrast, when comparing females with and without ASD, the phenotypic code that demonstrated the highest OR of 131 (95% CI, 111.943-153.279) was intellectual disability.
For diagnosis in adolescence, we found 18 phenotypic codes, most of them of neurological and psychiatric. Of these, 10 were also reported in other age groups. Digestive issues were not more likely to occur in females with ASD in this age group, while lack of normal physiological development unspecified (ASD female and ASD male: When comparing phenotypic codes that were nonspecific to age (Figure 3), we found 6 phenotypic codes shared by children with ASD between 3 and 18 years of age (strabismus, other cerebral degenerations, chromosomal anomalies, and epilepsy-related disorders). Specifically, strabismus (mostly esotropia, exotropia not otherwise specified, and accommodative component in esotropia) was more likely to occur in females with ASD when compared with females without ASD (in patients 3-5 years of age: OR, 4 Phenotypic codes that were statistically significant (corrected P value < .01) and more likely to appear in females with ASD, as compared with males with ASD and females without ASD (OR >1.5). For each phenotypic code, we presented the OR and the 95% CI (based on a comparison between females and males with ASD and females with or without ASD). Phenotypic codes present in this table were not statistically significant in females without ASD in the general population. The table is sorted according to female vs male OR in decreasing order.

Identification of phenotypes in the pediatric hospital dataset
Using the pediatric hospital dataset, we detected 7 phenotypic codes that were more likely to be present in females with ASD (Table 2, Supplementary Figure S2, Tables S4 and S5). Of those, 6 were identical to the ones found in the claims data, specifically chromosomal anomalies, developmental delays and disorders, other cerebral degenerations, epilepsy, recurrent seizures and convulsions, and epi-lepsy and convulsions. These 6 phenotypic codes were found as statistically significant only in the 6-to 11-years-old group. Only one phenotypic code milestone was specific to the pediatric hospital dataset, and it was only associated with sex dimorphism in patients 3 to 5 years of age.

DISCUSSION
In this study, we applied the multi-PheWAS intersection approach in a large-scale, real-world claims dataset and we identified 27 phenotypic codes more likely present in females with ASD, 10 not previously reported. We further validated our findings using an independent pediatric hospital dataset. With this approach, we over- Figure 3. Comorbidities more likely to appear in females with autism spectrum disorder (ASD) in the claims database by age groups: (A) novel comorbidities and (B) previously reported comorbidities. The y-axis on the left represents the phenotype categories (added for visualization purposes to sort the phecodes on the right in a meaningful way), the y-axis on the right represent the specific phenotype codes (phecodes). The x-axis represents the 4 different age groups. The cell color is related to the odds ratio when comparing females with ASD vs males with ASD. GERD: gastroesophageal reflux disease.
came one of the main limitations experienced in the field of ASD research: population sample size for females. 20,21 The strength of this approach was further enhanced by the performance of a systematic and large-scale analysis of comorbidity patterns. We were able to analyze all possible ASD comorbidities at the PheWAS level, without having to reduce the analysis to a specific set. Our results support prior findings. We showed that neurological disorders, such as epilepsy, were more likely to occur in females with ASD, as compared with females without ASD and with males with ASD, in all children 3 to 18 years of age. Epilepsy was also significantly prevalent in the pediatric hospital dataset for patients 6 to 11 years of age. Supekar et al 23 and Amiet et al 30 also reported a higher prevalence of epilepsy in females due to the increased prevalence of cognitive and motor deficits. These results align with our analysis: we found the prevalence of epilepsy and intellectual disabilities to be twice as high in females with ASD compared with males with ASD for groups 3 to 18 years of age (Supplementary Table S6). Notably, females and males without ASD exhibited the same prevalence of the 2 comorbidities, between 0.03% and 0.04%. Our study also reported muscle weakness to be more prevalent in females without ASD (OR, 1.7 [95% CI, 1.37-2.10]), confirming Ben-Itzchak's study, 33 in which musculoskeletal deficits were significantly higher overall in females (73.8%) than in males (57.1%).
Two other comorbidities have been previously reported as being more prevalent in females with ASD compared with males: urinary incontinence (OR, 1.627 [95% CI, 1.376-1.918] in ASD females and ASD males; and OR, 24.3 [95% CI, 21.014-28.085] in ASD females and non-ASD females, 12-18 years of age) and strabismus (OR >1.8 in ASD females and ASD males 3-18 years of age; OR, 4.9 [95% CI, 4.718-5.948] in ASD females and non-ASD females 3-5 years of age; OR, 5.3 [95% CI, 4.718-5.948] in ASD females and non-ASD females 6-11 years of age; OR, 8.5 [95% CI, 7.122-10.042] in ASD females and non-ASD females 12-18 years of age). However, the underlying reason for the discrepancy was the use of risperidone for urinary incontinence 31 and the use of valproate for strabismus. In the claims dataset, none of the patients with strabismus were treated with valproate, and for the individuals with urinary incontinence, a higher percentage of males (17.49%) than females (14.08%) with ASD were treated with risperidone. Other underlying causes could explain the higher number of females with ASD that show these 2 comorbidities.
Davignon et al 32 analyzed comorbidities within a dataset of 4123 individuals with ASD and more than 1000 distinct ICD-9 codes, grouped into 175 conditions. However, the study did not focus on a pediatric population. Other ASD comorbidity studies had a small sample size and did not factor in sex differences. They were also restricted to a predefined comorbidity list, derived mainly from the most well-known ASD comorbidities. 20,23,25 Finally, other studies have focused only on one comorbidity: suicidal ideation 34 or epilepsy. 23,30 By contrast, our study simultaneously targeted all existing comorbidities without relying on a predefined list of ASD comorbidities.
Beyond supporting prior findings, our comorbidity analysis further uncovered evidence of 10 novel phenotypes-gastroesophageal reflux disease; constipation; lack of normal physiological development; failure to thrive; dysthymic disorder; myopia; symptoms concerning nutrition, metabolism, and development; fever of unknown origin; dehydration; and abdominal pain. While known to be associated with ASD, these were not previously associated with sex dimorphism in ASD. These results could be explained by selection bias: females are only diagnosed when they have more severe symptoms, that is, epilepsy, owing to underlying pathology in females. Still, the signs and symptoms that we identified likely signal underlying common pathways and the need for further study.
When we compared the 2 datasets, claims and pediatric hospital, we found that for the age group 6 to 11 years old, they shared 6 of the 17 comorbidities first identified in the claim data. Out of the 7 comorbidities found in the pediatric hospital, 6 were identified in the claim dataset. The one present in the pediatric hospital but not in the claims dataset, "Delayed milestones" in the age group 3 to 5 years old, might be explained by the distinct type of claims data we used. Boston Children's Hospital treats a disproportionate number of patients with complex and rare conditions. By contrast, the nationwide claim database contains patients from a more representative and national population. The high concentration of complexity in the overall pediatric hospital dataset patient population could have easily contributed to the differences appearing between datasets in this study.
In summary, sex an important specifier that could help in defining ASD subgroups, identifying mechanisms for disease and potential shared druggable pathways. This systematic, large-scale analysis of sex differences between comorbidities presented in individuals Phenotypic codes that were statistically significant (corrected P value < .01) and more likely to appear in females with ASD, as compared with males with ASD, and males and females without ASD (OR > 1.5). For each phenotypic code, we presented the OR and the 95% CI (based on the comparison between females and males with ASD and females with or without ASD). Phenotypic codes present in this table were not statistically significant in females without ASD in the general population. The table is sorted according to female vs male OR in decreasing order.
with ASD addresses the skewed sex ratio historically reported in ASD. Results support the notion that the apparent sex difference can be explained in part by differences in the nature of the phenotypic codes that manifest in males vs females with ASD. By understanding those phenotypes more likely present in females, ASD diagnosis in females may be improved. Such phenotypic categorization could help clinicians in guidance and specialist referrals appropriate for the specific anticipated comorbidities. Moreover, as clinicians account for sex as an important variable, it might aid in the identification of distinct comorbidity patterns that could define new ASD subtypes and the best course of treatment. These additional clinical correlates could also help to develop clinical decisive support models and to identify potential ASD subgroups.

Limitations
The sources of data used in this study may have introduced bias. Specifically, ICD-9-CM terminology was designed for billing and administrative purposes and therefore may include coding bias. Even if the sample size would have been bigger if applying less stringent thresholds (109 305 patients with 1 ASD code, 80 374 patients with 2 ASD codes), we preferred to reduce this potential bias by mapping to a greater level of aggregation, using PheWAS codes, and limiting the eligibility of patients to those diagnosed with the same code at least 3 independent times. Using datasets such as claims data allows you to be more restrictive and still have a big sample size . We were interested to show that this is possible. We increased the specificity by decreasing the number of patients. This is important for psychiatric conditions such as autism. Taking as cases only those patients diagnosed on 3 different occasions with the phenotype, we reduce the possibility that, for example, the code was used related to ASD testing, and not necessarily to ASD diagnosis. We would not use the same approach if our target disorder would be an acute one (eg, myocardial infarction). In this case, one code approach would have been robust enough. Patients only diagnosed 1 or 2 times with the phenotype were considered neither cases nor controls, and were excluded from the analysis. Because we searched for associations without a hypothesis, there is a chance that some findings may involve noise. There also may be systematic bias and potential misclassification in the underlying database. For instance, in our results, dementia appeared in the claims data as a statistically significant phenotypic category of patients younger than 18 years of age (0.001%).
We also assert that claims data and electronic health records have been recognized as excellent sources of information for comorbidity clinical research. 35 In fact, data repurposing provides an unprecedented opportunity to use healthcare data for additional research purposes. Finally, the alignment of our results with those from previous reports supports the robustness of our method and the validity of our results.
In this study, we have focused on sex differences, understanding sex as the biological sex, not the social construct gender. Recent studies have shown that transgender and gender-diverse individuals have higher rates of autism diagnosis and traits as well as related neurodevelopmental and psychiatric conditions. 36 Future studies should include gender as a covariate. Although gender is not automatically recorded in claims and electronic health records, gender dysphoria is. In our dataset we did not find statistically significant differences in terms of gender dysphoria when comparing females with ASD and males with ASD, but it is statistically significantly higher in females with ASD than in females without ASD for the age group 12 to 18 years of age. Further studies would be required and in this specific case without aggregating at a higher level. Using the ICD-10 codification, F64-related codes could be used as a starting point when using this medical data. This will broaden the way we think about autism and other disorders, going from the binary idea of just "women" and "men".

CONCLUSIONS
The present study uniquely contributes to the growing body of research regarding sex differences in ASD by overcoming one of the main limitations of previous studies: small sample size. Further, we did not focus our analysis a priori on a selected group of specific disorders. Instead, we analyzed all clinical phenotypic codes applying the multi-PheWAS intersection approach on big real-world data, detecting comorbidity patterns in all patients with ASD to uncover potential sex differences. Overall, when looking for ASD subgroups, sex is an important specifier that can ultimately help in defining ASD subgroups and in the identification of novel mechanisms for disease and so advance our understanding of the disorder.

FUNDING
This work has been supported by the National Institutes of Health BD2K grant U54HG007963. JMZ received grants from Stichting de Drie Lichten and Stichting Sophia Kinderziekenhuis Fonds for a research internship at Harvard Medical School.

AUTHOR CONTRIBUTIONS
AG-S conceptualized and designed the study, collected the data, carried out the analyses, drafted the initial manuscript, and reviewed and revised the manuscript. CS and CDN provided feedback in the design of the study and statistical methods, and critically reviewed the manuscript for important intellectual content. NJ collected data, carried out the validation analyses, and reviewed and revised the manuscript. RK and TND supervised the data extraction methodology and reviewed and revised the manuscript. KPF, JMZ, NP, and IK provided critical feedback on the methods and results and critically reviewed the final draft for important intellectual content. PA conceptualized and designed the study, coordinated and supervised data collection, and critically reviewed the manuscript for important intellectual content. All authors approved the final manuscript as submitted and agree to be accountable for all aspects of the work.

SUPPLEMENTARY MATERIAL
Supplementary material is available at Journal of the American Medical Informatics Association online .

CONFLICT OF INTEREST
The autors have no biomedical financial interest or potential conflict of interest to disclose.

DATA AVAILABILITY STATEMENT
Individual participant data will not be available.