Familial ASD studies

Individuals included in our analyses were participants in one of three studies of high familial likelihood for ASD. These studies were designed to investigate the early life factors involved in ASD and neurodevelopmental outcomes by enrolling mothers who had already given birth to a previous child with a clinical ASD diagnosis. The Early Autism Risk Longitudinal Investigation (EARLI) Study [21] and Markers of Autism Risk in Babies-Learning Early Signs (MARBLES) [22] both enrolled mothers who already had a child with ASD to participate when they became pregnant with another child and followed them throughout pregnancy until 36 months of age. We also included participants from the Infant Brain Imaging Study (IBIS), [23,24,25] which enrolled these high likelihood mothers and infants at 6 months and were then re-evaluated at 12 and 24 months of age.

Population-based study to explore early development phase I

A total of 3,769 children were enrolled in the Study to Explore Early Development, Phase 1 (SEED 1), a multisite cohort initiative designed to obtain a representative sample of ASD and typically developing preschool-aged children in the US. Children between 2 and 5 years old, born between September 1, 2003, and August 31, 2005, and living in one of six study site vicinities (San Francisco Bay Area, Philadelphia metropolitan area, northeast Maryland, central North Carolina, and the Atlanta metropolitan area) were ascertained through a variety of methods, including diagnostic clinics, organizations providing evaluation or services for children with developmental problems, educational departments, and population vital records. Detailed recruitment procedures are described elsewhere [26].

Genetic dataCleaning and imputation

Samples from familial ASD cohorts were genotyped using the MEGA 1-Million (M)(MARBLES) and 5 M Illumina (EARLI and IBIS) chips. The Multi-Ethnic Global Array (MEGA) is a high-density array consisting of more than 1.7 million single nucleotide polymorphisms (SNPs) and is designed to represent diverse ancestries. The 5 M array is more comprehensive, consisting of around 5-million high-density SNPs, and has high overlap with the 1 M array. The SEED samples were genotyped on the more recently available Global Screening Array (GSA), which contains 640 K variants and represents diverse ancestry. For all studies, whole blood or buccal tissue was collected and samples were processed and stored. Genotyping was performed at the Johns Hopkins Genetic Resources Core Facility (GRCF).

The resulting genotypes were then subject to quality control filters according to standard criteria [27]. Briefly, using PLINK v1.9 [28], individuals failing QC checks for sex, relatedness, sample call rate of 0.03 (using –missing), and divergent ancestry were removed. SNPs were also removed based on the following criteria: MAF ≥ 0.05 among European samples, missing call rates exceeding 0.05, and Hardy–Weinberg equilibrium exact test p-value below 0.00001 among European samples.

Following cleaning, studies were imputed on the Michigan Imputation Server using the Minimac4 pipeline provided by the University of Michigan. We specified the 1000 genomes project (1000G) Phase 3 v5 reference panel, hg19 array build, Eagle v2.4 phasing [29], and the quality control and imputation mode. We imputed each target study separately and included 2504 1000G samples along with the target samples surviving QC. Post-imputation, we required imputed SNPs to correlate with the true unobserved genotypes at an r-squared (R2) > 0.80.

Genetic ancestry classification

After applying both the variant and sample level filters, measured genotypes were used to compute genetic ancestry variables from principal components analysis in Eigensoft [30] according to a recommended procedure [27], which includes pooling the target samples with the 2504 1000G samples from diverse ancestral groups [31].

Classification to a defined ancestral group was carried out using K-means (R v4.0.3 with “kmeans” function) based on the first 2 principal components (PC1 and PC2) resultant from this procedure. Only the 661 African, 504 East Asian, and 503 European samples from 1000G were used as anchors to define three ancestral clusters. The minimum, maximum and standard deviation of PC1 and PC2 for each of the three 1000G ancestry groups were computed. Target sample principal components were then compared to these values to classify into an ancestral group. European ancestry classification required that the target PC1 and PC2 value fell within 1.96 standard deviations of the minimum and maximum values for the 1000G European corresponding principal components. All other ancestries were classified as non-European.

Top ranking ASD discovery SNPs

This analysis uses the top 88 discovery variants by p-value that were implicated in the ASD discovery GWAS [12]. The discovery GWAS combined samples from the Danish iPSYCH population as well as samples from the Psychiatric Genomics Consortium (PGC), totaling 18,381 cases of ASD and 27,969 controls reflecting European ancestry. Following initial identification of 88 top loci, a replication analysis was conducted with another 2,119 cases of ASD and 142,379 controls pooled across five different populations of Northern European ancestry. The three identified loci were found to be significant, while two additional loci became significant when meta-analyzed with the discovery sample. Although the majority of the single variant tests did not achieve statistical significance, a test to replicate the direction of effects was significant. For each identified top variant, the minor allele frequency, p-value, and odds ratios were provided in Supplementary Table 1.

Identifying correlated SNPs

In addition to the original list of 88 top variants, nearby SNPs that were highly correlated with the discovery index variants were identified as proxies. We accomplished this by accessing the GRCH37 REST API database (http://grch37.rest.ensembl.org) via R software. The Ensemble database uses the 1000G phase 3 reference to perform searches for correlated SNPs in specific windows. Proxy searches were performed specifying a reference panel made up of samples with European ancestry. Proxies were kept if they were found to be in linkage disequilibrium (LD) with the index SNP at R2 >  = 0.80. When multiple proxies were available for a missing SNP, the proxy was selected based on highest R2 with the index SNP and closest physical distance.

Calculating single variant coverage among top ASD hits

Representation of the original index or proxy SNP was determined for cleaned and imputed target datasets by evaluating the overlap using chromosome, base pair, and variant identifier (rs number). An identical process was applied to obtain coverage for the Global Diversity Array-8 (GDA) and the Infinium PsychArray. Because we did not have target samples typed on these arrays, we were only able to evaluate coverage for variants on the pre-cleaned manifest files. The manifest files we downloaded for evaluation are publicly available on the Illumina website ( [1, 2]).

Literature search for published reports using an ASD polygenic score

To characterize the methods for deriving and describing ASD-PGS that are commonly employed by researchers to date, we conducted a literature search in PubMed to identify manuscripts published through October 2022 that reported on ASD-PGS associations, where the target sample was in children and the outcome was either ASD or an ASD-related trait. We searched on the terms “ASD” and “Polygenic Risk Score” to obtain 109 potential hits. We also supplemented our search with the same terms in Google, evaluated the suggested literature from each identified manuscript, and considered references returned by the PGS Catalog [32] when searching for “autism”. Two researchers evaluated each abstract to rule out those studies that did not report on an ASD-PGS, did not perform analyses in children, or did not report on child ASD status or child ASD trait. We did not consider randomized trials. After a review of abstracts, 36 potential manuscripts met our criteria, and of these, 24 were selected for extraction. Extraction included the method and software used to derive ASD-PGS, parameters specified for the method, and whether the number of SNPs informing the ASD-PGS was reported.

Creation of an information metric for polygenic scores

While there exist several methods to derive PGS, the majority of researchers employ the clumping and thresholding (C + T) method [33], where redundant SNPs due to linkage disequilibrium (LD) are removed (i.e. clumping), and only the SNPs that fall below an established discovery p-value threshold level are included in the score (i.e. thresholding). Scores are derived using PLINK software [28, 34], which can also be implemented via PRSice version 1 [35] or version 2 [36]. To derive a score, the researcher must specify: 1) a reference panel, 2) a target panel and 3) a target population. The reference panel will be used to determine the amount of LD between nearby SNPs, supervised by discovery effect sizes or p-values, to make decisions about which SNPs to prune in the “clump” procedure. The target panel will limit the choice of SNPs to those available after cleaning and/or imputing the target array data. Thus, together, the reference and target panels will incorporate the discovery information to select the list of SNPs that inform the PGS. Finally, using the clumped SNP list, the “score” command will sum and weight each genotype for each sample in the target population.

With this C + T method, if a SNP, or any proxy SNP, is not available to represent a genomic region, its failure to be incorporated into the PGS calculation is likely to result in information loss. Even in the case when a proxy SNP is selected from the target panel to represent the index discovery SNP, some information loss will occur. Our goal was to assess the impact of these missing representations of discovery loci due to lack of coverage when using a C + T method to derive ASD-PGS, even after high quality imputation recovered a number of SNPs in the target data that were not present on the cleaned, measured GWAS array.

To accomplish this, we made use of the 1000G phase 3 v5 reference panel [31], which contains a comprehensive set of genetic variants from whole genome sequencing across 26 different populations. We limited the reference panel to the 498 unrelated individuals of European ancestry (eur1kg). The eur1kg panel can serve as the reference panel when clumping, and can also act as a high coverage, comprehensive target panel, in addition to serving as the target population. Alternatively, the high quality post-imputed data from a target cohort can also serve as a target panel, and the samples can be used as the target population in whom scoring takes place. We compute two different PGS as outlined in Fig. 1. All PGS are informed by the Grove et al. discovery results, and the eur1kg panel is always specified as the reference panel of choice. The first score represents the score computed from a comprehensive panel and specifies the eur1kg as both the target panel and the target population (full-eur1kg PGS); the process to create full-eur1kg PGS is shown in Fig. 1 on the left. The second score represents the PGS computed using the cohort target panel, which is not as comprehensive as eur1kg. This second PGS is scored in the eur1kg sample to serve as an “apples to apples” comparison to the full-eur1kg PGS, as shown on the right of Fig. 1 (referred to as target-eur1kg PGS). Our coverage metric is the correlation between the full-eur1kg and target-eur1kg PGS, which will range from 0.0 to 1.0.

Fig. 1

PGS Computation Workflow. Legend: Our genome-wide information metric is simply the direct correlation between full-eur1kg and target-eur1kg PGS, which reflects information loss genome-wide. A correlation of 1.0 indicates that no information loss occurred, whereas a low correlation suggests a substantial loss of information. We derive each ASD-PGS using a C + T method, limiting to biallelic, high quality (info > 0.80) SNPs, for a suite of p-value discovery thresholds (5 × 10–8, 1 × 10–6, 10–4,10–3, .01, 0.05, 0.10, 0.20, 0.50 and 1.0) scored in PLINK software. Because target-eur1kg-ASD-PGS were clumped and scored for each cohort separately, a different SNP selection informed each of the target-eur1kg PGS, leading to cohort-specific correlations with the full-eur1kg PGS