ADuLT: An efficient and robust time-to-event GWAS

2.1 Model

The ADuLT model is an extension of the classical LTM[13, 10], and is the model underlying our previously proposed LT-FH++ method[1]. To estimate an individual’s genetic liability, ADuLT utilises birth year, sex, phenotype-specific age-of-onset for cases and current age for controls, as well as population-based cumulative incidences (i.e. the probability of having developed the disease at a given age). In contrast to LT-FH++, the ADuLT model does not incorporate family history as presented here. Instead, we focus on comparing ADuLT to standard time-to-event GWAS methods. ADuLT can account for cohort effects (changes in disease incidence by birth year), as well as differences by sex. This however requires population-based estimates to be available by age, sex and birth year for each phenotype of interest.

The ADuLT model extends the classical LTM by allowing the threshold used to determine case-control status to depend on sex, birth year, and (if available) age-of-onset for an individual. The LTM assumes that each individual has a liability ℓ that follows a standard normal distribution in the population. When this liability is larger than a given threshold, ℓ ≥ T, where P (ℓ ≥ T) = K and K is the trait’s lifetime prevalence, then the individual is a case (z = 1), otherwise it is a control (z = 0). This model does not account for time-to-event. Under the ADuLT model, the trait prevalence K becomes the available population-based cumulative incidence stratified by sex and birth year, if this information is available. In Figure S12, an example of those CIPs can be seen for depression. Additionally, we assume that the liability can be decomposed into two independent components, a genetic component, ℓ_g, and an environmental component, ℓ_e, such that ℓ = ℓ_g + ℓ_e. The genetic liability ℓ_g is normally distributed with mean 0 and variance h², where h² denotes the trait heritability on the liability scale. The environmental component is normally distributed with mean 0 and variance 1 − h² and independent of ℓ_g.

ADuLT aims to estimate an expected genetic liability. We do this by expressing the liability as a 2-dimensional normal distribution given by:

The mean of the genetic component is given by where the information we condition on, namely the case-control status and CIPs, result in an interval of (full) liabilities to integrate over. The CIPs set the threshold and the case status determines if the integration is above (a case) or below (control) the threshold.

2.2 Simulation Details

The default simulation setup uses two generative models, namely the Cox proportional hazards model and the LTM. We simulate under both generative models in order not to favour one method over the other.

Initially, genotypes are simulated for N = 1, 000, 000 individuals and M = 20, 000 independent SNPs. The genotypes are sampled from a binomial distribution Binom(2, AF) with the probability parameter set to the allele frequency (AF) of a given SNP. The AFs are sampled from a uniform distribution on the interval (0.01, 0.49). SNPs are standardised using the true AF, and for the scaled SNPs, the effect sizes of causal SNPs were drawn from the normal distribution N (0, h²/C), where C denotes the number of causal SNPs and h² denotes the liability-scale heritability. In the simulations, we used h² = 0.5 and either C = 250 or C = 1000 causal SNPs. With the simulated genotypes and causal effect sizes, we then assigned synthetic phenotypes using the two generative models.

For the proportional hazards model, we opted for a simulation setup as similar as possible to the one used in SPACox[5]. First, we simulated the censoring times, c_i, for each individual i from an exponential distribution with a scale parameter of 0.15. Next, we simulated onset times[4], , using a Weibull distribution[2] as follows where λ is the event rate, U_i ∼ Unif(0, 1), , with ϵ_i ∼ N (0, 1 − h²), and are the scaled genotypes multiplied by effect sizes, corresponding to the genetic liability ℓ_g in the LTM. The case-control status z_i is then 1 if , and 0 otherwise. The event time is the observed time. The event rate λ was chosen such that the lifetime prevalence is fixed at e.g. 1% or 5%. The simulation of onset times depends on all causal SNPs, which deviates from the simulations of onset times in the SPACox paper, where the onset times depended on a single causal SNP only. This change was made in order for the full genetic load of an individual to influence the onset times, instead of just a single SNP. Next, we calculated the CIP of the simulated event times, i.e. the fraction of cases observed before a given point in time. then the proportions were converted to the ages-of-onset (in years) using the logistic function given by Equation (1) with median age-of-onset x₀ = 50 and growth rate k = 0.2. Both age and age-of-onset were used to calculate the cumulative incidence proportions, which in turn defines the thresholds under the ADuLT model. For instance, with a lifetime prevalence of 1%, 90% of all individuals had an age or age-of-onset between 17 and 57 years.

Under the LTM, we set the trait status z_i equal to 1 if the liability exceeds the threshold, i.e. if ℓ_i > T, and 0 otherwise, where . The threshold T is determined by the lifetime prevalence K. For instance, a lifetime prevalence of 5% and 10% results in thresholds T = 1.64 and T = 1.28, respectively. The relationship between the age-of-onset and the liabilities above the threshold T, is given by the logistic function where K denotes the maximal attainable value, k denotes the growth rate, and x₀ denotes the median age-of-onset. Using the age of controls, we know how long they have lived without being diagnosed. This information allows us to exclude liabilities, i.e. the period of risk lived through so far. For both cases and controls, the personalised thresholds are calculated as T_i = Φ(1 − CIP_i), where T_i is the personalised threshold and CIP_i is the CIP for individual i. The liabilities below the personalised threshold are considered for controls and the liabilities above the threshold are considered for cases. If the population-representative CIPs are stratified by birth year and sex, the full liaiblity for cases can be fixed at T_i. Ages for controls are sampled from a uniform distribution between 10 and 90. This resulted in 90% of individuals having an age between 14 and 86.

2.3 The ADuLT survival model

As we showed previously[1], the age-dependent liability threshold model can be considered a survival model. More specifically, consider the survival function S_i(age) = P (age_i > age), where age_i represents age-on-set for cases or censoring time for the i^th individual. The probability that an individual has not become a case for a given age is equal to the probability that the individual’s liability is larger than the (individualised) liability threshold T_i(age), which is a shorthand notation for the age-dependent threshold given by

. Here sex_i and birth year_i are the i^th individual’s sex and birth year, respectively. If we assume that the individual liability consists of a genetic and an environmental component, , where and are Gaussian distributed with mean 0 and variance h² and 1 − h², respectively, then we can write the survival function as follows where Φ is the standard Gaussian cumulative distribution function and we assume that the genetic liability contribution is known. In the last equality, we standardise the environmental contribution with the known genetic contribution and the variance. From this we can derive the event density, and the hazard function for the i^th individual as

We note that this survival model is unusual in a couple of ways. First, each individual has a slightly different parameterisation of the model, which comes through the individualised liability threshold T_i(age). Second, the genetic effects affect the hazard rate by shifting the individual liability. Third, T_i(age) does not have to approach negative infinity as age approaches positive infinity, but may instead simply become fixed for all values T_i(age) above some threshold, e.g. if every individual in a cohort has died and no new event are possible. This is not necessarily a problem for the interpretation as T_i(age) may still be piece-wise differentiable, and the hazard rate for all values t above this threshold then becomes 0.

2.4 GWAS in iPSYCH

With the second wave of genotyped individuals, the iPSYCH case-cohort reached ∼143,000 individuals, up from ∼80,000.[6] Both waves have been imputed with the RICOPILI imputation pipeline[26], and were then combined into a single dataset. We restricted the analysis to SNPs that passed RICOPILI quality controls for both waves, resulting in a total of 8,785,478 SNPs for the GWAS. The analysis was restricted to a group of individuals with European ancestry, which were identified by calculating a robust Mahalanobis distance based on the first 20 PCs and restricting to a log-distance below 4.5[43]. We filtered for relatedness by removing individuals (the second one in each pair) with a KING-relatedness above 0.088. Since the iPSYCH case-cohort has a population representative subcohort and oversampled cases for six major psychiatric disorders (here we focus on ADHD, autism, depression and schizophrenia), we restricted each analysis to the individuals in the subcohort (which is a random sample of the entire population) and the cases for the phenotype being analysed, i.e. oversampled cases from the other psychiatric disorders were not used. The final number of individuals used for the GWAS of each phenotype is presented in Table S2. The linear regression GWAS was performed using the bigsnpr package[41] for R and SPACox GWAS was performed using the original implementation in the SPACox package for R. We used 20 PCs, sex, and imputation wave as covariates for all analyses. We included age as a covariate when analysing case-control status. Age was not included as a covariate when using the ADuLT phenotype or SPACox. We chose not to use a mixed model approach for GWAS with case-control status or ADuLT phenotypes, as SPACox did not have a similar option for random effects.

2.5 Cumulative Incidence Proportions

The CIPs can be interpreted as the proportion of individuals diagnosed with a certain disorder before a given age. As a result, the CIPs are population and disorder specific and can be stratified by sex and birth year. The CIPs used here were stratified by sex and birth year to account for differences in incidences between sexes and for different birth years (cohorts). The CIPs were estimated from Danish population-based registers. The Danish Civil Registration System[39] was used to identify individuals and contains all 9,251,071 individuals that lived in Denmark at some point between April 2, 1968 and December 31, 2016. The Danish Civil Registration System has continually recorded information since its launch in 1968, and includes information about sex, date of birth, date of death, and date of emigration, or immigration. Each individual has a unique identifier that can be used to link information of several registers. Information on psychiatric disorders was obtained from the Danish Psychiatric Central Research Register[33]. It contains all admissions to psychiatric inpatient facilities since 1969 and visits to outpatient psychiatric departments and emergency departments since 1995. From 1969 to 1993, the International Classification of Diseases, eighth revision (ICD-8) was used as the diagnostic system. From 1994 onwards, the tenth revision (ICD-10) was used. The four disorders of interest were identified by the following ICD-8 and ICD-10 codes: ADHD (308.01 and F90.0), autism (299.00, 299.01, 299.02, 299.03 and F84.0, F84.1, F84.5, F84.8, F84.9), depression (296.09, 296.29, 298.09, 300.49 and F32, F33), and schizophrenia (295.x9 excluding 295.79 and F20). The age-of-onset was defined as the age of an individual at first contact with the psychiatric care system, either inpatient, outpatient, or emergency visits. In the analyses, each individual was followed from birth, immigration, or January 1, 1969 (whichever happened last) until death, emigration, or December 31, 2016 (whichever happened first). The cumulative incidence function was estimated separately for each sex and birth year, and the Aalen-Johansen approach was used with death and emigration as competing events[15].

3.1 Overview of method

The age-dependent liability threshold model presented here was first introduced in our previous paper extending the LT-FH method to account for family history as well as age-of-onset, sex, and cohort effects among all individuals, including the family members[21, 1]. In this paper, we focus on the ADuLT model as an alternative to commonly used time-to-event or linear regression GWAS methods, without considering any family history.

The ADuLT model modifies the LTM by assuming that the threshold used to determine an individual’s case-control status corresponds to the CIP at the age of diagnosis. In Figure 1, we present the CIPs for ADHD for individuals born in Denmark in the year 2000. The CIPs increase as the population gets older, which in turn leads to a decreased threshold. If additional information, such as sex and birth year, is available, the population CIPs should be stratified according to this additional information (as seen in Figure 1), as this improves estimation of the genetic liability[1]. In the first step, a personalised threshold is assigned to each individual based on their current age or the age-of-onset, as well as sex and birth year. In the second step, the ADuLT model uses the liability-scale heritability to estimate a genetic liability for each individual. The third step uses the ADuLT phenotype as a continuous outcome in a GWAS. There are no restrictions on the choice of GWAS method as long as it accepts continuous outcomes, allowing researchers to benefit from current and future advances in GWAS methods. Note that Figure 1 illustrates the use of CIP for cases. If an individual is a control, the area of possible liabilities will instead be from negative infinity to the threshold identified from the CIPs.

• Download figure
• Open in new tab

Figure 1: Overview of ADuLT and illustration of the information it can account for.

Overview of the information used, and the different steps needed to perform a GWAS based on the ADuLT phenotype. The cumulative incidence proportions (CIPs) stratified by sex and birth year (here ADHD for individuals born in Denmark in 2000) are converted to a threshold for the age-dependent liability threshold model. Females are represented by the red line, while males are represented by the blue line. The CIPs has been marked at the age of 10 and 15 for both sexes (dotted lines). Finally, a genetic liability is estimated for each individual, and this ADuLT phenotype can be used as the outcome in a GWAS.

3.2 Simulation Results

We used two generative models for the simulations, namely the LTM and the proportional hazards model (see Methods). The performance of a simple case-control GWAS, SPACox, and the ADuLT phenotype used as the outcome in a linear regression-based GWAS was assessed under both generative models. Sex and age or age-of-onset were simulated for 1 million individuals, each with 20,000 independent SNPs. To examine the effect of ascertainment of cases, which is common in GWAS data, similar analyses were performed where the total number of individuals was randomly downsampling from 1 million to 20,000 individuals, leaving 10,000 controls and 10,000 cases in each downsampled dataset.

Figure 2 displays the power for each method under both generative models with 250 causal SNPs. A similar plot showing the power of the same generative models but with 1000 causal SNPs can be found in Figure S1. Without downsampling, the power of all three methods is similar under both generative models (Figure 2A). In Figure 2B, which is based on a downsampled data set simulating case ascertainment, the power of all three methods decreased due to a reduced sample size, but the power of SPACox was disproportionately affected by the downsampling. For simulated traits that have been downsampled and have a lifetime prevalence of 5% or below, SPACox performs worse than linear regression for both the case-control status and the ADuLT phenotype by more than a factor of 10 in the worst case, and approximately 25% worse in the best case. Under the proportional hazards model and a lifetime prevalence of 20%, and with downsampling, SPACox has an average power on par with ADuLT.

• Download figure
• Open in new tab

Figure 2: Power simulation results with 250 causal SNPs under both generative models and varying prevalences.

The power is shown for different population prevalence, varying from 1% to 20%. A) The power, i.e. the fraction of causal SNPs detected for each method, without downsampling. B) The power with downsampling, i.e. the number of individuals is subsampled to 10k cases and 10k controls.

In Table S1, which is based on data simulated under the LTM and without downsampling, the relative power of all methods are within 3% of one another. ADuLT obtained the highest power in all cases, while the lowest power was observed in connection to SPACox. With downsampling and 1000 causal SNPs, the increase in power was 117% with ADuLT over SPACox across all prevalences considered, and it was 96% for case-control status over SPACox. With downsampling and 250 causal SNPs, we observed an average increase of 34% in power over SPACox with ADuLT, and a 29% increase in power with case-control status, showing that SPACox has a comparatively low power for low effect sizes.

In Figures S2 and S3, the average χ²-statistics for the null SNPs is reported. Plots were achieved for 250 and 1000 causal SNPs, respectively, and each plot contains results for four different lifetime prevalences, with and without downsampling, and for both generative models. All models were well calibrated, since no inflation of the null statistics is observed. Figures S4 and S5 show the relative power with SPACox as the baseline method. As before each plot holds the relative power for four lifetime prevalences, with and without downsampling and both generative models. In addition, different plots were achieved for 250 and 1000 causal SNPs. For 250 causal SNPs and no downsampling, performance of all methods were similar. However, with downsampling, SPACox only identified a few causal SNPs, which resulted in large relative power gains for ADuLT and linear regression (see Figure 2B). In Table S1, simulation results for all parameter setups are available, including the power, relative power compared to SPACox, and mean chi-square statistic of null SNPs.

3.2.1 Computation Times

The computational time for estimating the ADuLT phenotype depends solely on the number of individuals. The running time for the GWAS step depends heavily on the implementation of the GWAS method used. In Figure 3, the combined running times of estimating the ADuLT phenotype and performing a GWAS using the bigsnpr package[41] are reported. We used 4 CPU cores for both steps, which is a conservative number of cores. The SPACox implementation does not support parallelization, which is why SPACox was run sequentially. We find that ADuLT together with a linear regression is faster than SPACox, even with only modest parallelization. Logistic regression of a binary phenotype is slower than linear regression of the same phenotype[41], which means ADuLT together with a linear regression may be faster and have higher power to detect causal SNPs.

• Download figure
• Open in new tab

Figure 3: Running times of ADuLT combined with a linear regression GWAS compared to SPACox.

Each point represents the mean value of 10 replications, while the error bars are represented by the estimate ±1.96 standard errors. Run times were assessed for a varying number of individuals and SNPs (M). SPACox uses a single CPU core, as no parallelization is available. We used 4 CPU cores for estimating the ADuLT phenotype and performing the linear regression GWAS for this phenotype. The means and corresponding standard errors of the run times can be found in Table S1.

3.3 GWAS of psychiatric disorders in iPSYCH

The iPSYCH data has been linked to the Danish registers, which means that detailed information on age-of-onset, age, sex, and birth year can be assessed for all genotyped individuals that are part of the iPSYCH cohort[6] This supplementary information was used to analyse four psychiatric traits, namely ADHD, autism, depression, and schizophrenia. For each of these phenotypes, population-based CIPs were obtained by birth year and sex (see Figures S6, S9, S12 and S15 for plots of the CIPs used, and see Cumulative Incidence Porportions for details). The prevalences were used to tailor the thresholds to each individual under the ADuLT model (see Methods).

We performed GWASs for each of the four phenotypes and for each of the methods considered, i.e. using either the case-control status or the estimated genetic liability by ADuLT as the outcome in a linear regression-based GWAS or SPACox (see Methods for details). Figure 4 displays the Manhattan plots for ADHD for all methods, where the case-control GWAS included age as a covariate, while the ADuLT GWAS and SPACox did not. To report nearly independent findings, LD clumping was performed on the summary statistics with a r² threshold of 0.1 and a window size of 500kb, prioritising the SNPs with the lowest p-values. This was done for each combination of phenotype and method. The lowest p-value LD-clumped SNPs that are unique to ADuLT and ADHD can be found in Table S3 and the LD-clumped snps that are unique to case-control status and ADHD can be found in Table S4. For ADHD, we found 12 independent genome-wide significant associations when using the ADuLT phenotype as the outcome, while case-control status and SPACox found 11 and 5 associations, respectively. The ADuLT GWAS had two independent associations that were not identified by case-control associations, and case-control GWAS found one association that was not found by the ADuLT GWAS. One of the associations unique to ADuLT is rs4660756. The gene closest to this SNP is ST3GAL3, which has previously been associated with educational attainment[36] and ADHD[52]. SPACox also identified ST3GAL3, but through rs11810109 instead. The association unique to case-control GWAS is rs8085882 on chromosome 18. The closest gene is ZNF521, which has previously been associated with education attainment[29], ADHD[46], and smoking initiation[27]. The association with the lowest p-value that is shared among all methods is rs4916723 on chromosome 5 with LINC00461 as the closest gene. This gene has also been reported as being associated with educational attainment[27] and ADHD[9].

• Download figure
• Open in new tab

Figure 4: Manhattan plots from GWAS with the ADuLT phenotype, SPACox, and case-control status for ADHD.

Manhattan plots for ADHD for all three methods. Case-control GWAS uses the age of individuals as a covariate, whereas the ADuLT GWAS and SPACox do not. The orange dots indicate suggestive SNPs with a p-value threshold of 5 × 10⁻⁶. The red dots correspond to genome-wide significant SNPs with a p-value threshold of 5 × 10⁻⁸. The diamonds correspond to the lowest p-value LD clumped SNP in a 500k base pair window with an r² = 0.1 threshold.

Across the four psychiatric disorders, ADuLT found 20 independent genome-wide significant associations, while case-control status found 17 and SPACox found 8. The Manhattan plots for each of the methods, each of the remaining disorders (autism, depression, and schizophrenia), and with and without age as a covariate can be found in Figures S7, S8, S10, S11, S13, S14, S16 and S17. Notably, SPACox consistently identified fewer associations than the ADuLT and case-control status GWASs, and was the only method that did not identify any significant association for major depression and schizophrenia.