|Home | About | Journals | Submit | Contact Us | Français|
We read with interest Garriock and colleagues’ article “A Genomewide Association Study of Citalopram Response in Major Depressive Disorder” (1). This work is an important step in analyzing the Sequenced Treatment Alternatives to Relieve Depression (STAR*D) genomewide SNP data. Given STAR*D’s complexity, many additional analyses will undoubtedly be required to realize the full potential of this valuable resource. In this spirit, we quantitatively model response as measured by the Quick Inventory of Depressive Symptoms (QIDS). The gains to quantitatively modeling response are evident in our results—unlike the initial analysis, we identify multiple genomewide significant SNPs.
We apply three well-established psychometric insights to improve treatment response measurement. First, dichotomizing a quantitative variable sharply reduces statistical power (2, 3). For instance, in regression involving two normally distributed variables, dichotomizing the predictor at its median reduces variance explained by 38%, with further reductions as the dichotomization point moves away from the median (3). Thus, dichotomizing citalopram response likely substantially decreases power to detect SNP effects. Second, when measuring change (i.e., treatment response), increasing number of assessments per subject monotonically increases estimate precision, with dramatic gains when assessments per subject are numerous, as in STAR*D (4, 5). Thus, modeling all QIDS assessments (~7.2 per subject) increases response estimate precision more than four-fold compared to using only first and last assessments (5). Third, as questionnaire items generally vary in how well they describe underlying constructs (6), factor analyses can improve measurement through estimating the association of each QIDS item to the depressed affect construct. Cumulatively, these techniques help maximize power to find genetic effects, which is crucial in pharmacogenomics given tremendously expensive clinical trials and relatively small sample sizes.
Details of STAR*D and GWAS methods have been previously described (1, 7, 8). Briefly, we control for ancestry using five genomewide MDS dimensions as GWAS covariates and exclude subjects with ambiguous sex assignment. Otherwise, we employ QC criteria used by Garriock and colleagues.
We begin by calculating two quantitative treatment measures as GWAS outcomes. For the first measure, we model all QIDS-Self-report (QIDS-S) repeated assessments to estimate each individuals’ response. Specifically, using a mixed model-based method described previously (9, 10), we first determine the optimal model of overtime drug response by fitting a series of models specifying linear symptom change for a given number of days on drug and flat thereafter. This series begins with a model specifying maximal response achieved on day one. Each subsequent model specifies incrementally longer duration until maximal response, with the final model assuming no plateau in drug effect (i.e., constant linear change). The best-fitting model of the series indicates number of days until maximal response. Consistent with recent meta-analyses indicating most SSRI benefit is evident at 2–3 weeks with additional benefits tapering thereafter (11, 12), the optimal plateau is estimated at 21 days for the QID-S. Next, 55 covariates (including design, socio-demographic and clinical variables) are screened to identify those improving treatment effect precision, with only trial design characteristics and age at first MDD diagnosis improving precision. Finally, based on this modeling, treatment effects quantifying each subject’s citalopram-induced symptom change are output as best linear unbiased predictors (BLUPs) (9, 13).
For outcome two, we model the QIDS’ factor structure. Thus, as opposed to simply summing all questionnaire items, the association of each symptom to the underlying depressed affect construct is estimated and used to weight the item’s contribution to the efficacy measure. Additionally, we jointly analyze the 16 QIDS-Clinician-rated (QIDS-C) items with the 16 QIDS-S items to further refine measurement. Structural equation modeling of the QIDS items indicates the best-fitting model (CFI = 0.92, RMSEA = 0.05) consists of 3 factors—sleep-related, weight/appetite, and depressed affect, and includes correlations between identical QIDS-S and QIDS-C items (6). After generating factor scores for the depressed affect factor (for brevity we consider only the depressed affect factor, as it explains most item variance and is conceptually most central), they are analyzed using the procedure described above to generate treatment effects. The optimal plateau is estimated at day 30. Only covariates describing trial design and age at first MDD diagnosis improve precision.
Q-Q plots for the two GWAS communicate three important points (Figure 1). First, the points’ general linearity suggest no inflation in p-value distributions. This is confirmed by genomic control statistics (λ=1.007, 1.004), indicating population stratification is well-controlled. Second, the divergence of linearity in the plots’ upper right corners suggests associations stronger than expected by chance, which is confirmed using either FDR (q<.1) (14, 15) or the initial paper’s Bonferroni (p<1.16×10−7) genomewide significance criterion (Table 1). The improvement of these results over the initial analysis is further evident in comparing the number of SNPs with p<1×10−5—whereas the former analysis yielded 7 SNPs, our analysis yields 15. Finally, comparing the Q-Q plots in Figure 1 demonstrates the advantage of modeling the QIDS’ factor structure—the same SNPs appear as top findings for both outcomes, but p-values are superior for the factor score.
The most significant SNP (rs11143665) and neighboring significant markers (rs11143678, rs11143679 and rs11143683) are located on 9q21.13 (Table 1). These markers center on UniGene cluster Hs.243837, which is expressed in brain, with one hit, rs11143683, located directly in the expressed sequence. A second group of findings (rs9952222, rs11873656 and rs7229673) are located in a gene-rich region on 18q11.2, including histamine receptor H4 (HRH4 <270kb centromeric) and pseudogene RAC1P1 (~50–60kb). Other findings include rs10512091 on 9q21.31, ~15kb from hypothetical gene LOC347119 and ~85kb from transducin-like enhancer split 4 (TLE4), which has high expression in brain and a role in neuronal development (16). Finally, rs7836900, at 8q12.1, is intergenic, and rs12673398 is in a tensin 3 (TNS3) intron on 7p12.3. Two of the five regions indicated here (9q21.31 and 18q11.2) correspond, with stronger significance, to top findings from Garriock and colleagues initial paper.
In conclusion, we have built upon Garriock and colleagues’ work to show the value of quantitatively modeling treatment response. We hope this both introduces promising candidate loci and promotes methods for increasing power in pharmacogenomic studies.
Daniel E Adkins gratefully acknowledges support from NIMH (1F31MH079653).
Drs Adkins, Åberg, McClay, Hettema, Bukszár and van den Oord report no biomedical financial interests or potential conflicts of interest. Dr Kornstein has received research support from Pfizer, Bristol-Myers Squibb, Lilly, Forest Laboratories, Wyeth, Novartis Pharmaceuticals, Boehringer-Ingelheim, and Takeda. She has served on Advisory Boards for Pfizer, Wyeth, Lilly, Bristol-Myers Squibb, Forest Laboratories, Endo, and Takeda; and has received book royalties from Guilford Press.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.