Clustering by phenotype and genome-wide association study in autism

Autism spectrum disorder (ASD) has phenotypically and genetically heterogeneous characteristics. A simulation study demonstrated that attempts to categorize patients with a complex disease into more homogeneous subgroups could have more power to elucidate hidden heritability. We conducted cluster analyses using the k-means algorithm with a cluster number of 15 based on phenotypic variables from the Simons Simplex Collection (SSC). As a preliminary study, we conducted a conventional genome-wide association study (GWAS) with a data set of 597 ASD cases and 370 controls. In the second step, we divided cases based on the clustering results and conducted GWAS in each of the subgroups vs controls (cluster-based GWAS). We also conducted cluster-based GWAS on another SSC data set of 712 probands and 354 controls in the replication stage. In the preliminary study, which was conducted in conventional GWAS design, we observed no significant associations. In the second step of cluster-based GWASs, we identified 65 chromosomal loci, which included 30 intragenic loci located in 21 genes and 35 intergenic loci that satisfied the threshold of P < 5.0 × 10−8. Some of these loci were located within or near previously reported candidate genes for ASD: CDH5, CNTN5, CNTNAP5, DNAH17, DPP10, DSCAM, FOXK1, GABBR2, GRIN2A5, ITPR1, NTM, SDK1, SNCA, and SRRM4. Of these 65 significant chromosomal loci, rs11064685 located within the SRRM4 gene had a significantly different distribution in the cases vs controls in the replication cohort. These findings suggest that clustering may successfully identify subgroups with relatively homogeneous disease etiologies. Further cluster validation and replication studies are warranted in larger cohorts.


Introduction
Autism spectrum disorder (ASD) has heterogeneous characteristics in terms of both phenotypic features and genetics. ASD is mainly characterized by difficulties in communication and repetitive behaviors, but ASD also shows many other symptoms 1 . Regarding genetics, previous studies have not consistently identified genetic variants that are associated with an increased risk of ASD 2 , although several lines of evidence suggest that genetic factors strongly contribute to the increased risk of ASD. Monozygotic twins have higher concordance rates of ASD (92%) than dizygotic twins (10%) 3 . The recurrence risk ratio is 22 for ASD among siblings 4 . The Human Gene module of the Simons Foundation Autism Research Initiative (SFARI) gene provides a comprehensive reference for suggested human ASD-related genes in an up-to-date manner 5 and currently demonstrates~1000 genes that may have links to ASD, potentially indicating the heterogeneity of ASD. In addition to phenotype and genotype heterogeneities, ASD shows heterogeneous responses to interventions. Several kinds of pharmacological treatments are suggested, but the effects of these treatments are controversial 6 .
If the heterogeneous phenotypes and responses to treatment in some way correspond to differences in genotype, grouping persons with ASD according to phenotype and responses to treatment variables may increase the chances of identifying genetic susceptibility factors. Traylor and colleagues 7 demonstrated that attempts to categorize patients with a complex disease into more homogeneous subgroups could have more power to elucidate the hidden heritability in a simulation study. Several studies on Alzheimer's disease, neuroticism, or asthma indicated that items or symptoms were to some degree more useful for identifying high-impact genetic factors than broadly defined diagnoses [8][9][10] , although a study of ASD demonstrated modest effects of two-way stratification by individual symptoms 11 . In addition, medical researchers have begun to use machine learning methods, which is an artificial intelligence technique that can reveal masked patterns of data sets. In view of the abovementioned circumstances, clustering algorithms of machine learning and subsequent genome-wide association studies (GWASs) could be hypothesized to reveal novel and more genetically homogeneous clusters, but a combinatorial approach of cluster analysis and GWASs, to the best of our knowledge, has not been applied to any diseases including ASD.
We therefore explored whether grouping persons with ASD using a clustering algorithm with phenotype and responses to treatment variables can be used to discriminate more genetically homogeneous persons with ASD. In the present study, we conducted cluster-based GWASs (named cluster-based GWASs) using real data based on the concept of a previous simulation study 7 adopting a machine learning k-means 12 algorithm for cluster analysis.

Subjects and methods
We conducted the present study in accordance with the guidelines of the Declaration of Helsinki 13 and all other applicable guidelines. The protocol was reviewed and approved by the institutional review board of Tohoku University Graduate School of Medicine, and written informed consent was obtained from all participants over the age of 18 by the SFARI 14 . For participants under the age of 18, informed consent was obtained from a parent and/or legal guardian. In addition, for participants 10-17 years of age, informed assent was obtained from the individuals.

data sets
We used phenotypic variables, history of treatment, and genotypic data from the Simons Simplex Collection (SSC) 14 . The SSC establishes a repository of phenotypic data and genetic data/samples from mainly simplex families.
The SSC data were publicly released in October 2007 and are directly available from the SFARI. From the SSC data set, we used data from 614 affected white male probands who had no missing information regarding Autism Diagnostic Interview-Revised (ADI-R) scores 15 and vitamin treatment 16,17 and 391 unaffected brothers for whom genotype data, generated by the Illumina Human Omni2.5 (Omni2.5) array, were available for subsequent clustering and genetic analyses. We excluded participants whose ancestries were estimated to be different from the other participants using principal component analyses (PCAs) performed by EIGENSOFT version 7.2.1 18 for the genotype data. Based on the PCAs, we excluded data beyond four standard deviations of principal components 1 or 2 ( Supplementary Fig. 1). Therefore, we used data from 597 probands and 370 unaffected brothers.
In the replication study, we used another SSC data set genotyped using the Illumina 1Mv3 (1Mv3) array. In the data set, data from 735 affected male probands with no missing information regarding ADI-R scores or vitamin treatment and 387 unaffected brothers were available. After conducting PCA, we excluded data beyond four standard deviations of principal components 1 or 2 as outliers. In this way, we used data from 712 probands and 354 unaffected brothers in the replication study.

Clustering
We conducted cluster analyses using phenotypic variables of ADI-R scores and history of vitamin treatment 16,17 . We chose these variables because the ADI-R is one of the most reliable estimates of ASD and has the ability to evaluate substructure domains of ASD 15 . Among the ADI-R scores, "the total score for the Verbal Communication Domain of the ADI-R minus the total score for the Nonverbal Communication Domain of the ADI-R", "the total score for the Nonverbal Communication Domain of the ADI-R", "the total score for the Restricted, Repetitive, and Stereotyped Patterns of Behavior Domain of the ADI-R", and "the total score for the Reciprocal Social Interaction Domain of the ADI-R" were included in the preprocessed data set.
Among the treatments, we selected the variable of history of vitamin treatment because we recently found that a cluster of persons with ASD is associated with potential responsiveness to vitamin B6 treatment 16,17 . The history of treatment is not always compatible with responsiveness, but we considered that continuous treatment indicates responsiveness to some degree. The SSC data set includes history of treatment but not variables of responsiveness.
We applied the machine learning k-means 12 algorithm to conduct a cluster analysis to divide the data set obtained from ASD persons into subgroups using phenotypic variables and history of treatment. The k-means algorithm requires a cluster number (k) determined by researchers. We set a priori k of 5, 10, 15, and 20 under the hypothesis that ASD consists of hundreds of subgroups 5,14 and considering statistical power by sample size calculations 19 . We performed the analyses using the scikit-learn toolkit in Python 2.7 (Supplementary Information 1).
Clustering is an exploratory data analysis technique, and the validity of the clustering results may be judged by external knowledge, such as the purpose of the segmentation 20 . Several methods have proposed to prespecify a cluster number of k, such as visual examination of the data, and likelihood and error-based approaches; however, these methods do not necessarily provide results that are consistent with each other 21 . Although there are measures for evaluating the quality of the clusters 22 , the number of clusters should still be determined according to the research purposes. We regarded the inflation factor (λ) of quantile-quantile (Q-Q) plots of the logarithm of the P value to base 10 (−log 10 P) as one of the indicators of successful clustering in the present study. We calculated λ for each cluster number.
When conducting clustering, we combined the two data sets of male probands, one genotyped using the Omni2.5 array and the other genotyped using the 1Mv3 array. After clustering, we redivided the new data set according to the SNP arrays used. In the discovery stage, we used the Omni2.5 data set and the 1Mv3 data set in the replication stage.

Genotype data and quality control
We used the SSC data set, in which probands and unaffected brothers had already been genotyped in other previous studies 14,23 . In the discovery stage, we used the data set genotyped by the Omni2.5 array, which has 2,383,385 probes. We excluded SNPs with a minor allele frequency < 0.01, call rate < 0.95, and Hardy-Weinberg equilibrium test P < 0.000001.
In the replication study, where we used the data set genotyped using the 1Mv3 array, we applied the same cutoff values for quality control as those used in the discovery stage. The 1Mv3 array includes 1,147,689 SNPs. The Omni2.5 array and the 1Mv3 array shared 675,923 SNPs.

Statistical analysis
As a preliminary study, we conducted a conventional GWAS in the whole Omni2.5 data set, with a total of 597 male probands and 370 unaffected brothers. Here, we used the brothers of the cases as controls, in contrast to many previous studies in which genetically unrelated controls were used. We thus adopted the sib transmission disequilibrium test (sib-TDT) 24 , a family-based association test, to take into account familial relationships among the participants. In the second step, in the discovery stage, we conducted cluster-based GWAS in each subgroup of the cases, which had been divided using the k-means 12 algorithm, and the controls. As mentioned above, the controls were the brothers of the cases, and we then excluded the unaffected brothers of the cases belonging to the subgroup being analyzed. Details of the study design are shown in Supplementary Fig. 2. We applied the Cochran-Armitage trend test 25 , which examines the risk of disease in those who do not have the allele of interest, those who have a single copy, and those who are homozygous.
We further tested the significantly associated loci found in the discovery studies in the replication stage. The level of significance for association was set as P < 0.05 in the replication studies.
Association analyses were performed with the PLINK software package 26 . The detected SNPs were subsequently annotated using ANNOVAR 27 . Manhattan plots and Q-Q plots were generated using the 'qqman' package in R version 3.0.2.

Cluster-based GWAS
As a preliminary study, we conducted a conventional GWAS with the Omni2.5 data set using the sib-TDT. We observed no significant associations (Fig. 1). Although we adopted the sib-TDT here because we used the brothers of the cases as controls, we also used the Cochran-Armitage trend test and found that the −log 10 P values were distributed downward compared with the expected values, as shown in Supplementary Fig. 3.
We also applied the sib-TDT to cluster 1, which was obtained by dividing all the cases using k-means with k of 15, and all the controls and found that the observed −logP values were lower than expected, as shown in Supplementary Fig. 3. As the sib-TDT may efficiently work in a population consisting of a substantial number of sibs, a limited number of brothers of the probands among all the controls probably contributed to a substantial loss of power. Thus, we excluded the brothers of the probands in each subset from the controls so that each subset of probands has no genetic relations with the rest of the controls and conducted the Cochran-Armitage trend test, as in many other studies. In the present study, therefore, we applied the sib-TDT to the GWAS of the whole data set, whereas in the cluster-based GWAS, we excluded in turn the unaffected brothers of the cases belonging to the subgroup being analyzed and used the Cochran-Armitage trend test to account for the relationships between participants.
The average inflation factor λ for the cluster-based GWAS with k of 5, 10, 15, and 20 were 1.021, 1.024, 1.038, and 1.053, respectively. Several lines of evidence suggest that regarding an appropriate threshold of λ, empirically, a value <1.050 is deemed safe for avoiding false positives 28 . Under the hypothesis that ASD consists of hundreds of subgroups 14 , we compared λ values giving larger numbers of clusters as priority. We therefore considered the cluster-based GWAS using k-means cluster analysis with k of 15 to be the most appropriate approach to the present data set. The characteristics of each cluster are presented in Table 1.

Gene interpretation
We observed 65 chromosomal loci that satisfied the threshold of P < 5.0 × 10 −8 (Fig. 2); 30 out of the 65 loci were located within 21 genes, and the remaining 35 loci were intergenic ( Table 2). Among them, eight loci were located within or near the genes associated with the Human Gene module of the SFARI Gene scoring system 5  The SFARI Gene scoring system ranges from "Category 1", which indicates "high confidence", through "Category 6", which denotes "evidence does not support a role". Genes of a syndromic disorder (e.g., fragile X syndrome) related to ASD are categorized in a different category. Rare single-gene variants, disruptions/mutations, and submicroscopic deletions/duplications related to ASD are categorized as "Rare Single Gene Mutation".

Replication study
We conducted replication studies with another independent data set that included a total of 712 male Table 1 Characteristics of each of 15 k-means clusters in the Omni2.5 data set.
Cluster no.  probands and 354 unaffected brothers and had been genotyped using the 1Mv3 array. As mentioned before, we had previously carried out cluster analyses in the combined data set genotyped with either Omni2.5 or 1Mv3 and then redivided it according to the SNP arrays used. The characteristics of each of the 15 clusters in the 1Mv3 data set are presented in Supplementary Table 1.
Among the 65 genome-wide significant chromosomal loci found in the discovery study, seven chromosomal loci were included in the 1Mv3 array. Of these loci, rs11064685, within SRRM4 in Cluster 13, had a significantly different distribution (p = 0.03) in cases vs controls in the replication cohort (Table 3).

Discussion
One of the most important findings of our study was that reasonably decreasing the sample size could increase the statistical power. A plausible explanation is that our clustering may have successfully identified subgroups that are etiologically more homogeneous. At least two reasons could reduce the possibility of false positives of the present results of statistically significant SNPs in clusterbased GWAS. First, the present study validated the usefulness and feasibility of the concept of a previous simulation study 7 , which indicated that homogeneous case subgroups increase power in genetic association studies by Traylor and colleagues, using measurement data in the real world. Second, a substantial number of statistically significant SNPs in cluster-based GWAS observed in the present study were located within or near previously reported candidate genes for ASD 5,[30][31][32][33][34][35] .
We observed many statistically significant SNPs in cluster-based GWAS: CDH5, CNTN5, CNTNAP5, DNAH17, DPP10, DSCAM, FOXK1, GABBR2, GRIN2A5, ITPR1, NTM, SDK1, SNCA, and SRRM4. In particular, loci within the SRRM4 gene had significantly different distributions in the cases vs controls in the replication cohort. Previous studies indicate that SRRM4 is strongly associated with ASD, indicating that our results may be valid to some degree. The gene regulates neural microexons. In the brains of individuals with ASD, these microexons are frequently dysregulated 48 . In addition, nSR100/SRRM4 haploinsufficiency in mice induced autistic features such as sensory hypersensitivity and altered social behavior and impaired synaptic transmission and excitability 49 .
In addition to SRRM4, we observed several genes located within or near previously reported candidate genes for ASD. The relatively high correspondence between our results in part and the SFARI Gene scoring system 5 indicates that the statistically significant loci we found may be associated with ASD subgroups (Fig. 2). We also observed several important genes associated with ASD and other related disorders 29 from previous reports. Powers were calculated using the method based on the results in Nam's study 19 .
These findings suggest that the statistically significant SNPs might explain autistic symptoms because these diseases are suggested to have shared etiology, even in part, with ASD 29 . Associations at the remaining significant loci that were not in the SFARI module or described above have not been previously reported, and to the best of our knowledge, some of them might be novel findings. These results might suggest that novel genetic loci of ASD could be found by identifying better defined subgroups, although further confirmation is needed in future cohorts with larger sample sizes.
Previous studies regarding Alzheimer's disease, neuroticism, or asthma found that items or symptoms showed, to some degree, increased ORs between the case loci and control loci compared with those from previous studies using broadly defined disease diagnoses [8][9][10] . These findings may indicate that GWAS based on a symptom or an item could identify genetically more homogeneous subgroups and let us hypothesize that a relatively reasonable combination of symptoms or items could identify more genetically homogeneous subgroups. In contrast, Chaste and colleagues showed that stratifying children with ASD based on the phenotype only modestly increased power in GWAS 11 . The discrepancy between their findings and ours might be explained by usage of phenotype variables. Chaste and colleagues used one item or symptom alone with limited number of subgroups, whereas we used combinations of them with a machine learning method with a potentially sufficient number of clusters. DeMichele-Sweet and colleagues reported that subgrouping only by having psychosis could lead to the identification of limited loci that had small effects 50 , but Mukherjee and colleagues found a substantial number of suggestive loci that had extreme ORs after categorizing persons with Alzheimer's disease based on relative performance across cognitive domains by modern psychometric approaches 8 .
Validation of clusters is essential. In the present study, we selected the k-means algorithm, focused on ADI-R items and treatment as variables, and determined cluster numbers based on the λ of the Q-Q plots. Although we believe this approach is one of the relevant ways, selection of variables, selection of algorithms and selection of cluster numbers still remain to be considered in future mathematical and biological cluster validation studies because controversies surrounding evaluation of the quality of the clusters are important issues and are still ongoing and because validated clusters may lead to elucidate the genetic architectures of ASD 7 .
The present study has a limitation to be noted. Substantial differences in the two genotyping platforms may have affected the results of the replication study. The Omni2.5 array includes 2,383,385 autosomal SNPs, whereas the 1Mv3 array includes 1,147,689 SNPs, with 675,923 shared SNPs between the two. Of the 65 statistically significant chromosomal loci in the discovery data, only seven chromosomal loci were shared between the two arrays.
Our study demonstrated that if the data set consists of multiple heterogeneous subgroups, even a subgroup that includes a much smaller number of homogeneous individuals could detect high-impact genetic factors. Hypothetical examples of the concept of cluster-based GWAS are shown in Supplementary Fig. 4. As shown in Table 2, only 30 etiologically homogeneous probands and 300 controls can have a statistical power of~1.00, calculated using the method based on the results in Nam's study 19 . Although the integral model, which assumes many genetic variants have a small effect, may contribute to the formation of some subgroups of ASD, our results indicate that clustering by specific phenotypic variables may provide a candidate example for identifying etiologically similar cases of ASD.
Our data indicate the relevance of cluster-based GWAS as a means to identify more homogeneous subgroups of ASD than broadly defined subgroups. Future investigation of cluster validation and replication with a larger sample Table 3 Results of replication studies in the 1Mv3 data set for statistically significant chromosomal loci in the discovery studies. size is therefore warranted. Such studies will provide clues to elucidate the genetic structures and etiologies of ASD and facilitate the development of precision medicine for ASD.