Estimating the genetic variance of Major Depressive Disorder (MDD) due to all SNPs

A limited number of meta-analyses combining multiple Genome Wide Association (GWA) studies of psychiatric disorders have yielded replicated associations between such disorders and Single Nucleotide Polymorphisms (SNPs) (e.g., schizophrenia, [1]). In contrast, analyses of Major Depressive Disorder (MDD) lack genome-wide significant results [2,3]. Previously, we conducted a genome-wide association (GWA) study in 1,738 MDD cases and 1,802 controls selected to be at low liability for MDD [2]. None of the SNPs reached genome-wide significance. Subsequent GWA and meta-analysis studies of MDD were equally unsuccessful in detecting significant SNPs [3]. The lack of breakthrough results has spurred a discussion concerning the utility of GWA studies of psychiatric phenotypes. Part of this discussion focuses on whether SNPs analyzed in GWA studies are the relevant polymorphisms, or whether psychiatric disorders traits are more likely influenced by rare variants.

Yang et al. [4] and So et al. [5] added to this discussion by proposing two different methods which focus on the estimation of phenotypic variance explained by currently genotyped SNPs rather than on the detection of specific SNPs that contribute to a trait or disorder. Using human height as an example, Yang et al. [4] demonstrated that SNPs that are currently included on genotyping platforms explain a substantial proportion of the heritability of the phenotype. So et al. [5, 6] reached similar conclusions. The two approaches differ substantially, with Yang et al.’s approach falling under the broad umbrella of variance decomposition methods, and So et al.’s approach being a density estimation method. Both methods focus on estimating variance due to additive genetic effects.

The rationale of Yang et al.’s method is to decompose the variance of the phenotype into a component that is due to the additive effects of all measured SNPs, and a residual component that is due to random noise, unmeasured environmental influences, or effects of unmeasured genetic variants that are not captured by the measured variants. Their 2-step method involves first obtaining the genetic similarity between all pairs of subjects, which is calculated as a correlation. For each pair of subjects, their minor allele counts on a given SNP are subtracted from the SNP mean, and the cross-product of the resulting deviance scores is divided by the SNP variance. Summing these terms over all SNPs gives a measure of genetic similarity between two subjects that is due to additive genetic effects. The second step of the method uses this measure as a random effect to predict the phenotype in a linear mixed model. The original approach for continuous traits has been extended to case/control studies [7,8].

The method proposed by So et al. is entirely different, and is applied subsequent to a GWAS. The basic idea is to compare the distribution of z-statistics of the regression coefficients of genome-wide SNPs in a GWAS to the theoretical Null distribution of z-statistics representing no effects. Explained variance will differ from zero if more z-statistics from the GWAS have larger values than expected under the Null. Specifically, the observed z-statistics, which contain error due to sampling fluctuation, are first corrected to obtain z-statistics representing “true” effect sizes. The estimate of variance explained for continuous phenotypes can then be computed by summing contributions of all non-null “true” z-statistics using sums of squares as in ANOVA. For case/control phenotypes So et al. [6] have described a method using the odds ratio, allele frequency, and disease prevalence.

Our aim was to estimate the joint effect of all SNPs in explaining the variance in MDD. In addition, we investigated four other phenotypes that are currently targets of large GWA projects, namely height, fasting glucose, smoking initiation, and current smoking. The estimates of explained variance of all five phenotypes using the two different methods were compared to heritability estimates obtained in twin and family studies that were carried out in the same population.

The data analyzed in this study came from individuals who took part in the Netherlands Twin Register Biobank study [9] and the Netherlands Study of Depression and Anxiety [10]. Genotyping on these samples was performed on the Affymetrix 6.0 (N=298), Affymetrix 5.0/Perlegen (N=3697), Illumina 370 (N=290), Illumina 660 (N=1439) and Illumina Omni Express 1M (N=445) platforms. Quality control was carried out per platform. Thresholds for SNPs were MAF > 1%, HWE > 0.00001, missing >95% and 0.30 < Heterozygosity < 0.35. Samples were excluded from the data if their expected sex and IBD status did not match, or if the genotype missing rate > 10%. Subsequent to alignment using the Hapmap 2 Build 36 release 24 CEU reference set, and exclusion of SNPs if allele frequencies differed more than 15% with the reference set and/or the other platforms, the data were merged into a single dataset (N=5856). Imputation was done using IMPUTE v2, and SNPs were removed based on HWE < 0.00001, proper info < 0.40, MAF < 1%, allele frequency difference > 0.15 against reference and a missing rate > 5% assuming a 90% calling threshold. Individuals were removed with IBS > 0.025, resulting in a sample of 4605 individuals. A second round of QC was carried out on the imputed data closely following steps and criteria described in Anderson et al. [11] and Purcell et al. [12]. In addition, we removed SNPs that were significant in a Cochran-Mantel-Haenszel test of SNP-phenotype associations across platforms. For case-control phenotypes, we further removed SNPs that did not pass an adaptation of a two-locus test, which controls for differential case/control calling errors [13]. Quality control resulted in ~1,140,000 SNPs for the analyses. The exact number varied slightly for the different phenotypes since the two-locus and platform association tests are phenotype specific. The number of subjects varied between N=3245 and N=4240 depending on the phenotype. For MDD we used a stringent definition of case (i.e., diagnosis of MDD) and control (i.e., no reported psychiatric illness), resulting in 1620 cases and 1625 controls. The data include a large proportion (~2/3) of the original Genetic Association Information Network (GAIN) MDD data set, which was selected for MDD. Note that these subjects were genotyped on the same platform, with MDD cases and controls randomly assigned to plates [2]. Approximately 1/3 of the N=4240 sample that passed qc were not selected for MDD.

For both methods, we used software provided on the respective developers’ websites (http://gump.qimr.edu.au/gcta/, and http://sites.google.com/site/honcheongso/software/total-vg). Prior to applying the density estimation method we carried out LD pruning as suggested by So et al. [5], and used the kernel estimate as it is more stable in pruned data [5]. We repeated all analyses with different pruned sets and obtained very similar results, which demonstrates the robustness regarding the SNPs selected for the analysis.

Table 1 summarizes the results for all 5 phenotypes. For Yang et al.’s method likelihood ratio tests comparing models with and without the SNP effects were carried out, which were significant for all phenotypes. Additional permutation tests substantiated these results. Including 4 principal components to control for population stratification did not impact results. Note that Yang et al. [4], and also So et al. [5] pointed out that their estimates are conservative due to imperfect LD of typed with causal SNPs. Results in Table 1 are uncorrected for this underestimation, the extent of which might differ for the two methods. Although point estimates vary between the two methods, all 95% CI’s overlap.

The additive genetic variance of MDD due to all SNPs was estimated at 32% (95% CI: [15.0%, 48.8%]) and 28% (95% CI: [16.5%, 39.3%]), respectively. Even when considering the lower bounds of the confidence intervals, this is a substantial proportion of the heritability estimate of 36% established in a twin study this population [14], and that agrees with numerous other twin and family studies. The estimate concerning height using Yang et al.’s method is consistent with results obtained by Yang et al. [4]. A limitation regarding the two smoking phenotypes is that ~2/3 of the sample is selected by MDD status, and both smoking phenotypes in our sample have rank correlations of ~0.24 with MDD. Height and fasting glucose have essentially zero correlation with MDD.

The results concerning MDD are in stark contrast to GWA studies that aim at detecting specific SNPs [2,3]. In GWAS of MDD, SNPs explain less than 1% of the variance [15]. Our analyses show that SNPs typed on currently available platforms contain substantial information concerning the additive genetic variance of MDD and other complex phenotypes. The failure of MDD GWA studies is, among other reasons, likely due to small effect sizes of the involved SNPs, resulting in insufficient power to discriminate between signals and false positives [16]. In addition to improving power in genome-wide analyses by increasing the number of samples in meta-analyses, novel analytical approaches including data-mining procedures can be customized to efficiently extract the information that is present in genome-wide SNP data. Given appropriate power it should be possible to detect at least some of the locations in the genome associated with MDD, and to provide a basis for the understanding of the genetic underpinnings of this devastating disorder.