Increases in the efficiency of high-throughput genotyping have made it feasible for human geneticists to survey large numbers of genetic markers at a relatively low cost. A genome-wide association study (GWAS) is based on genotyping each individual in a sample on 100,000 or more common single nucleotide polymorphisms (SNPs) and then investigating the association of these SNPs with phenotypes of interest (McCarthy et al., Reference McCarthy, Abecasis, Cardon, Goldstein, Little, Ioannidis and Hirschhorn2008). The initial round of GWAS has been successful in identifying hundreds of genetic variants associated with a wide range of complex phenotypes (Hindorff et al., Reference Hindorff, Sethupathy, Junkins, Ramos, Mehta, Collins and Manolio2009; Visscher et al., Reference Visscher, Brown, McCarthy and Yang2012). The SNP variants that have been discovered have, however, characteristically had only small phenotypic effects. Consequently, very large samples are needed to ensure adequate statistical power, especially given the multiple testing burden (Manolio et al., Reference Manolio, Collins, Cox, Goldstein, Hindorff, Hunter and Visscher2009). As a consequence, consortia involving pooled GWAS samples, in some cases with sample sizes (N) exceeding 100,000, have been organized for several complex human phenotypes, including lipids (Kathiresan et al., Reference Kathiresan, Willer, Peloso, Demissie, Musunuru, Schadt and Cupples2009), height (Allen et al., Reference Allen, Estrada, Lettre, Berndt, Weedon, Rivadeneira and Hirschhorn2010), body mass index (Speliotes et al., Reference Speliotes, Willer, Berndt, Monda, Thorleifsson, Jackson and Loos2010), and age at menarche (Elks et al., Reference Elks, Perry, Sulem, Chasman, Franceschini, He and Murray2010).
GWAS on behavioral phenotypes suggests that effect sizes are likely to be similarly small for intellectual ability (Davies et al., Reference Davies, Tenesa, Payton, Yang, Harris, Liewald and Deary2011), personality (De Moor et al., Reference De Moor, Costa, Terracciano, Krueger, de Geus, Toshiko and Boomsma2010), and psychopathology (Sullivan et al., Reference Sullivan, Daly and O'Donovan2012). Pooled GWAS samples are clearly also needed in the behavioral domain to achieve adequate statistical power to identify the specific genetic variants implied to exist by twin, adoption, and family studies. Consortia are already organizing around several key behavioral traits, including educational attainment, intellectual achievement, personality, and a range of mental health and substance use disorders (Cichon et al., Reference Cichon, Craddock, Daly, Faraone, Gejman and Kelsoe2009). The purpose of this paper was to describe a relatively large (N > 8,000) sample that has been deeply phenotyped longitudinally at a single center and which has recently been genotyped on a GWAS array. The genotyping and initial analysis of the data were supported through the US National Institute on Drug Abuse's Genes, Environment and Development Initiative (GEDI). GEDI's major goal is to characterize the nature of gene–environment interplay in the development of substance use disorders. The sample provides numerous opportunities for collaborative research.
The Minnesota Center for Twin and Family Research
Minnesota Center for Twin and Family Research Longitudinal Samples
The GWAS sample is drawn from participants in one of three longitudinal studies undertaken under the auspices of the Minnesota Center for Twin and Family Research (MCTFR): (1) the Minnesota Twin Family Study (MTFS; Iacono et al., Reference Iacono, Carlson, Taylor, Elkins and McGue1999); (2) the Sibling Interaction and Behavior Study (SIBS; McGue et al., Reference McGue, Keyes, Sharma, Elkins, Legrand, Johnson and Iacono2007); and (3) the Enrichment Study (ES; Keyes et al., Reference Keyes, Malone, Elkins, Legrand, McGue and Iacono2009). These studies utilized similar assessment protocols and a common sampling unit, a four-member family consisting of a pair of siblings and their rearing parents. The offspring in all three samples were initially assessed in adolescence and followed into at least early adulthood. In total, 9,827 individuals (5,001 offspring and 4,826 parents) have completed an MCTFR intake assessment. Table 1 provides a descriptive overview of the MCTFR samples.
TABLE 1 Descriptive Characteristics of the Minnesota Center for Twin and Family Research Longitudinal Samples
MTFS = Minnesota Twin Family Study; ES = Enrichment Study; SIBS = Sibling Interaction and Behavior Study. Assessments currently in progress are designated as such; follow-up assessments that have not yet been undertaken are designated as NA.
The MTFS sample consists of 1,197 monozygotic (MZ) and 684 like-sex dizygotic (DZ) twin pairs (including five additional members from triplet sets). All pairs were ascertained from Minnesota state birth records between 1971–1985 and 1988–1994. The sample includes two cohorts: one initially assessed at age 11 years (the younger cohort) and the other initially assessed at age 17 years (the older cohort). Intake and follow-up assessments of the twins were scheduled to coincide with major transitions in the lives of these adolescents and young adults. Target ages (rate of follow-up participation) are 11, 14 (92.9%), 17 (87.3%), 20 (88.6%), 24 (89.1%), and 29 (in progress) for the younger cohort and 17, 20 (88.7%), 24 (93.8%), and 29 (94.2%) for the older cohort. At intake, approximately 18% of the recruited MCTFR families refused our invitation to participate and, based on a brief survey with non-participants, we found that non-participating families differed minimally in parental education, parental occupational status, or parental mental health from participating families (Iacono et al., Reference Iacono, Carlson, Taylor, Elkins and McGue1999). The sample is, thus, broadly representative of the population of the state of Minnesota.
The SIBS sample includes 409 adoptive and 208 non-adoptive families, each consisting of a pair of adolescent siblings and their rearing parents. Adoptive families were recruited from records of the three largest adoption agencies in Minnesota, and non-adoptive families were recruited from Minnesota birth records. At least one offspring in every adoptive family was not genetically related to other family members, but many adoptive families had a biological offspring of one or both of the rearing parents. All offspring in non-adoptive families are full biological siblings. Among those families that were eligible to participate, 63% of adoptive and 57% of non-adoptive families completed an intake assessment. There are minimal differences in socio-economic status and in offspring mental health between participating and non-participating but eligible families (McGue et al., Reference McGue, Keyes, Sharma, Elkins, Legrand, Johnson and Iacono2007). Unlike the MTFS, it was not logistically possible to link SIBS assessments with specific targeted ages. Rather, on average, offspring in SIBS families were in mid-adolescence at intake, late adolescence at their first follow-up and early adulthood at their second follow-up. The first follow-up of the SIBS sample is complete (with a 94.2% participation rate) and the second follow-up is in progress.
The ES was designed to extend the MTFS by oversampling 11-year-old twins likely to develop substance use disorders by virtue of being high on an index of externalizing behavior. Like the MTFS twins, ES twins were located through Minnesota state birth records (1988–1994). However, roughly half of the ES sample was selected after the mothers in those families completed a screening interview that established that at least one member of the twin pair was high on the externalizing dimension. The final ES sample included 300 MZ and 199 like-sex DZ twin pairs and their parents. The scheduling of intake and follow-up assessments of the ES twins parallels that for the younger MTFS cohort, with assessments targeted at ages 11, 14 (93.2% participation), 17 (in progress), and 20 (in progress). Among families eligible for participation in ES, 75.4% completed an intake assessment and participation rates did not vary by screening status. Differences in socio-economic and demographic factors between participating and eligible but non-participating families were generally minimal (Keyes et al., Reference Keyes, Malone, Elkins, Legrand, McGue and Iacono2009). In addition, 48 pairs of MZ twins aged 14–16 at initial assessment were recruited from the same population as the ES families and assessed 1 year apart. The twins in these families were part of a pilot study on adolescent brain development (AdBrain) and completed a clinical assessment similar to that used in ES along with both structural and functional MRIs. Twins in this study were also included in the GWAS sample.
MCTFR Assessments
Although not identical, the assessment protocols in the three MCTFR studies overlap extensively. For all three studies, the intake assessment involved a daylong visit to our labs at the University of Minnesota. Assessments common to all three studies include (1) interview assessment of common child and adult mental and substance use disorders based on the Diagnostic and Statistical Manual of Mental Disorders (DSM-III-R and DSM-IV); (2) quantity and frequency assessment of licit (e.g., tobacco and alcohol) and illicit drug use; (3) self-report assessment of personality; (4) individually administered Wechsler IQ tests and reading scales from the Wide Range Achievement Test; (5) teacher ratings of behavioral problems, personality, and academic progress and grades; (6) anthropometric measures; and (7) self-report of environmental risk and protective factors, including peer group characteristics, family functioning, and socio-economic status. In addition, the MTFS and ES include a half-day psychophysiological assessment that includes electroencephalography, event-related potentials, and autonomic nervous system activity. The SIBS, which does not include a psychophysiological assessment, includes videotaped family and sibling interaction sessions.
DNA Samples
Although collection of DNA samples has been a routine component of the MCTFR assessments for many years, to meet GEDI requirements and supply a fresh DNA sample to the Rutgers University Cell and DNA Repository (RUCDR), it was necessary to obtain new DNA samples. Among the 9,515 MCTFR participants who were eligible to provide a DNA sample (e.g., were alive and had not withdrawn from the study), 7,278 (76.5%) provided a blood sample and an additional 567 (6.0%) provided a saliva sample. Most of the individuals who did not provide a DNA sample could either not be contacted within the time frame of the GEDI project or had privacy or general concerns over providing a DNA sample for a repository. The RUCDR followed standard DNA extraction and storage procedures in processing the MCTFR samples (Sahota et al., Reference Sahota, Brooks, Tischfield, Weiner, Gabriel and Stephens2007).
Genotyping Method and Quality Control Filters
Genotyping
Genome-wide genotyping was carried out using the Illumina Human660W-Quad array (Illumina, Inc., San Diego, CA) according to the manufacturer's protocol, as described in Li et al. (Reference Li, Sheu, Ye, de Andrade, Wang, Chang and Yang2010b). This Infinium HD Beadchip required 200ng DNA per sample and contains 657,366 variants, including 95,876 intensity-only markers for calling copy number variants that are not considered here. For quality control (QC) purposes, each 96-well plate included DNA from two members of a single three-member Centre d'Etude du Polymorphisme Humain (CEPH) family, with the two members rotating across plates, as well as a duplicate of another randomly selected sample on that plate. In addition to standard QC filters, we used GenCall scores, metrics of genotype reliability generated by the BeadStudio software (Illumina Corporation, San Diego, California), to assess sample quality (Cunningham et al., Reference Cunningham, Sellers, Schildkraut, Fredericksen, Vierkant, Kelemen and Goode2008). Specifically, both the GenCall_10, representing the 10th percentile, and GenCall_50, 50th percentile, scores were used to assess sample quality following standard guidelines. In addition, each sample was genotyped on a custom 96-plex panel using Illumina VeraCode chemistry (Lin et al., Reference Lin, Yeakley, McDaniel and Shen2009). This panel contains SNPs present on the Human660W-Quad and served as an additional check on quality control.
QC of Samples
Genotyping was attempted with a total of 7,681 samples, including 83 samples that were included as within-plate duplicates, 160 CEPH control samples, and 60 samples that were repeats of samples that had failed an initial genotyping attempt due to low call rate. Of the 60 repeated samples, 47 (78%) were successfully genotyped on the second pass, including 16 (76%) of 21 retested samples that had failed completely on initial pass with no genotype calls. The 7,438 non-control samples (i.e., including those that failed the first attempt at genotyping and were re-genotyped) were subjected to five separate QC screens. These screens, along with the number (%) of samples failing each are (1) non-calls in more than 5,000 markers (N = 130, 1.7%); (2) GenCall_50 score < 0.9009 (N = 8, 0.1%); (3) GenCall_10 score <0.75 (N = 6, <0.1%); (4) extreme heterozygosity or homozygosity (N = 1, <0.1%); (5) sample mix-ups or failure to confirm genetic relationship with other genotyped relatives (N = 15, 0.2%). In total, 160 (2.2%) of the non-control samples were eliminated, leaving a cleaned sample of 7,278.
QC of SNP Markers
There are a total of 561,490 SNP markers (including sex chromosome and mitochondrial markers) on the Illumina Human660W-Quad array, but a set of 1,508 SNPs could not be called in any batch. The remaining 559,982 markers were subjected to nine separate QC filters. These QC filters along with both the number and percentage of markers failing each are listed here: (1) identified by Illumina as a bad marker on the array (N = 23, <0.1%); (2) more than one mismatch in duplicated samples (N = 70, <0.1%); (3) call rate < 99% (N = 3,924, 0.7%); (4) minor allele frequency less than 1% (N = 19,999, 3.6%); (5) more than two Mendelian inconsistencies across families (N = 9,117, 1.6%); (6) significant deviation from Hardy–Weinberg genotype frequencies in the White sample at p < 10−7 (N = 1,200, 0.2%); (7) autosomal or X marker associated with participant sex at p < 10−7 (N = 16, <0.1%); (8) significantly associated with batch at p < 10−7 (N = 17, <0.1%); and (9) markers with more than two heterozygous calls if they were on the X chromosome in male samples or from mitochondrial DNA in the total sample (N = 160, <0.1%). A total of 32,153 markers, 5.7% of the markers attempted, failed one or more of these QC filters, leaving 527,829 markers that passed all QC filters.
Final GWAS Sample
Among the 7,278 genotyped samples that passed all QC filters were 1,127 samples from individuals with an MZ twin who had not been genotyped. The genotypes of the non-genotyped twins were set equal to the genotypes of their genotyped MZ co-twins, resulting in a final GWAS sample of 8,405 individuals. The genotyped sample includes 3,924 (47%) males and 4,481 (53%) females; 4,434 (53%) of the sample are in the offspring generation and 3,971 (47%) are in the parent generation. In total, 2,390 families are represented in the analysis. For purposes of analysis (see later), four-member families were divided among MZ twin families (N = 1156), DZ twin families (N = 665), SIBS families with two adopted offspring (N = 269), SIBS families with two biological offspring (N = 191), and SIBS families with mixed adopted and biological offspring (N = 109). In addition, there were 101 genotyped individuals, usually step-parents, who did not fit into any one of these five four-member family types. Table 2 provides an overview of the number of participants in each category. The family structure of the sample produces several unique features. In particular, it allows us to confirm the genetic relationships of a large proportion of the sample, providing an additional level of QC screen not existing in most GWAS samples. For 6,919 (82.3%) of the genotyped participants (not counting the 1,127 non-genotyped MZ co-twins), we were able to confirm their relationship with at least one genetic relative in the sample (e.g., DZ twin, parent–offspring). Even though there are 2,390 families represented in the GWAS sample, because of spouse pairs and adoptive families the maximum number of genetically unrelated individuals in the sample is 4,706. There are 1,631 genotyped spouse pairs and 1,404 families where both genetic parents and at least one offspring have been genotyped.
TABLE 2 Number of Families and Individuals in Final GWAS Sample
All families consist of two parents and two offspring, but some members of some families are missing genotype data. The terms ‘Bio’ and ‘Adopt’ are used with the SIBS families (see text) where offspring raised as siblings were either the biological offspring of the parents (Bio) or were adopted (Adopt). Adopt/Adopt means that the parents adopted both of their children, Bio/Bio means that both children were biological offspring of the parents, and Adopt/Bio means that one child was adopted and the other was the biological offspring of the parents.
Genetic Ancestry
In the MCTFR, ancestry was based initially on the ethnicity specified on a birth certificate, adoption records, or by self-report. Of the 8,405 GWAS samples, 7,599 (90.4%) self-reported as having primarily European ancestry (i.e., ‘White’), 382 (4.5%) as Asian, 83 (1%) as African American (i.e., ‘Black’), and 127 (1.5%) reported mixed ancestry. All other ethnicities had a self-reported frequency of <1% and there were approximately 1% with a missing self-report. Because the genetics of complex phenotypes can vary across different ancestral groups (Bamshad, Reference Bamshad2005), the primary sample used in our GWAS analysis will be comprised of individuals of European ancestry. However, to improve on the accuracy of the self-report data and to deal with missing self-report data, we ran EIGENSTRAT, http://genepath.med.harvard.edu/~reich/Software.htm (Price et al., Reference Price, Patterson, Plenge, Weinblatt, Shadick and Reich2006) analysis, extracting the first 10 principal components (PCs), to aid in the identification of a cluster of individuals with European ancestry.
The principal component analysis implemented in EIGENSTRAT is sensitive to pairs of close relatives, so one member of every such pair was excluded in initial computations, but the genotypes of relatives were then projected onto the components constructed from the set of unrelated subjects (this was accomplished using the ‘-w poplist’ option to smartpca.perl). We identified relative pairs, even a few between-family pairs, by running a PLINK ‘–genome’ analysis to produce a matrix of ‘Z1’ identity-by-descent (IBD) statistics for all pairs of individuals. The Z1 statistic is the estimated probability that two individuals share exactly one allele from a common ancestor, averaged across the autosomal genome. We chose to classify pairs of individuals as unrelated if their Z1 was less than 0.10, which would be roughly midway between first-cousin once-removed and second-cousin except that the ethnic structure of our sample tends to amplify distant relationships in the PLINK analysis (inflating Z1 somewhat, but inflating Z2 and Pi-Hat much more), especially for members of minority groups. The matrix of Z1 statistics was used to select the largest possible subgroup of unrelated individuals. This is known as the maximum clique problem in computer science and mathematics. To find the guaranteed maximum is extremely computationally demanding (it is an NP complete problem), so we used a greedy algorithm to achieve a reasonable answer in only a few seconds. This solution provides a subset of subjects such that no pairs have Z1 greater than 0.10. It is conceivable that a larger subset could be produced with much more computational effort, but the gain is likely to be minimal.
To identify a core sample of individuals with European ancestry for genetic analysis, we computed the Mahalanobis D 2 for each genotyped individual from the centroid of self-reported Whites using all 10 PCs. Figure 1 displays the distribution of the resulting D 2 values. The cluster of data points centered at a D 2 value of approximately 3,400 represents individuals who self-reported as Asian. We defined the European cluster as including all those individuals with D 2 < 84. We also confirmed that full siblings and DZ twins were never in separate clusters. After adding 1,087 MZ co-twins, the final sample of 7,702 individuals with European ancestry (91.6% of the total sample) included 101 who had missing self-reported ethnicity data and 46 whose ethnicity had been recorded as other than White (several of those cases were found to have been caused by data entry errors). Figure 2 provides a plot of the first two PCs from the EIGENSTRAT analysis. The first PC is a dimension anchored at one end by individuals of self-reported East Asian ancestry and at the other end by individuals of self-reported European ancestry. The second PC is a dimension that differentiates individuals of self-reported African ancestry from those of European ancestry. Additional PCs tended to distinguish other groups such as the Native Americans or those from India.
FIGURE 1 Two histograms of Mahalanobis distances (D 2) to the centroid of the European American (‘White’) group. Both histograms display the same data, but with different axis ranges and bin widths. Histogram A uses bins of width 42 so that the first two bars on the left include all 6,615 White subjects and no others. The subjects with D 2 near 1000 are mostly mixed-race, Hispanic, and Native American. The peak between 3,000 and 4,000 was caused by 378 Asian subjects, mostly Korean adoptees, and the subjects with D 2 exceeding 4,000 are Black. Histogram B used bin widths of 1 so that all visible bars (values less than 84 on the abscissa) represent counts of White subjects. All subjects with D 2 less than 84 were used in the core European American group for all initial GWAS analyses.
FIGURE 2 (Colour online) Scatter plot of the first two principal components from an EIGENSTRAT analysis of 10,000 markers with color coding to show reported ethnicity. Some of the gray dots of ‘Unknown’ ancestry had been reported as White, but they were sufficiently far from the White centroid that they have been recoded as Unknown. Note that ‘Mixed’ was a vague category that seems to include individuals with some European ancestry mixed with either African, Asian, or Native American ancestry. Several subjects have been recoded as White because of Mahalanobis D 2 less than 84, indicating close proximity to the centroid of the White group in principal component space (see Figure 1). These included 2 of 73 reportedly Black subjects, 5 of 60 Hispanic, 17 of 77 Mixed, 9 of 44 Native American, and 59 of 66 of Unknown ancestry. An additional 41 of 6,564 subjects reported to be White were recoded as Unknown based on the same analysis. The graph reflects the final ethnic classifications that followed the Mahalanobis distance analysis described in the text and in Figure 1.
Imputation
Untyped SNP genotypes were imputed using HapMap2 as the reference panel. Samples were first phased using Beagle (Browning & Browning, Reference Browning and Browning2009), which takes into account the family structure of the data. Genotypes were then imputed from the phased data using Minimac, which is a computationally efficient version of the MACH program (Li et al., Reference Li, Willer, Ding, Scheet and Abecasis2010a). HapMap2 provided 2,543,887 autosomal SNPs in the r22 reference panel and 64,621 X-chromosome SNPs in the r21 reference panel. Of those SNPs, 501,912 autosomal markers and 11,685 X-chromosome markers had been genotyped for our samples on the Human660W-Quad platform, thus leaving 2,041,975 autosomal and 52,936 X SNPs to be imputed. Of those 2,094,911 imputed SNPs, 96.4% had r 2 > 0.5, 90.9% had r 2 > 0.8, 85.2% had r 2 > 0.9, and 77.2% had r 2 > 0.95 in the European American sample.
Analytical Strategy
Most GWAS involve independently sampled individuals or families of fixed structure (e.g., parent–offspring trios), which can be analyzed using freely available software, such as PLINK, http://pngu.mgh.harvard.edu/~purcell/plink/ (Purcell et al., Reference Purcell, Neale, Todd-Brown, Thomas, Ferreira and Sham2007). The five different types of four-member families that comprise the MCTFR sample, however, present analytical challenges that cannot be handled efficiently within PLINK or other existing software for GWAS analysis. Consequently, we developed an efficient analytical approach for the MCTFR GWAS data based on a Rapid Feasible Generalized Least Squares (RFGLS) algorithm (Li et al., Reference Li, Basu, Miller, Iacono and McGue2011). Briefly, the analysis involves regressing a phenotype on a set of pre-specified covariates and a single SNP genotype (coded 0-1-2), repeated for each SNP in the GWAS. Because of the family structure, the residuals from the regression are non-independently distributed with a variance–covariance structure that depends on family type (e.g., it is not the same for MZ and DZ twin families). In RFGLS, rather than estimate the residual variance–covariance matrix for each family type in each SNP regression, we estimate it only once for each family type based on the model with covariates but no SNP effect. The residual variance–covariance matrix in the analysis of individual SNP effects is then fixed at the estimates from the model without any SNP effects. As we show in Li et al. (Reference Li, Basu, Miller, Iacono and McGue2011), because any SNP effect is likely to be very small, fixing the variance–covariance matrix in this way has no effect on the test statistics but greatly increases the speed of the analysis.
Future Plans and Summary
The MCTFR GWAS sample is large, >8,400, deeply phenotyped and longitudinally informative. We recently demonstrated the developmental utility of the MCTFR GWAS sample using height as a proof of principle (Vrieze et al., Reference Vrieze, McGue, Miller, Legrand, Schork and Iacono2011). Specifically, we used 176 SNP height variants identified in the GIANT consortium (Allen et al., Reference Allen, Estrada, Lettre, Berndt, Weedon, Rivadeneira and Hirschhorn2010) to create a polygenic score that accounted for 9.2% of adult height variance in our sample (comparable to the 10.5% accounted for in the original study). When the polygenic score was investigated longitudinally, however, we were able to show that almost all of its effect was on pre-pubertal height; that is, it was not significantly predictive of pubertal growth spurt. As consortia publish findings from pooled GWAS for other phenotypes, the rich longitudinal assessments in the MCTFR will allow us to similarly characterize these effects developmentally. The MCTFR also includes a rich assessment of environmental risk and protective factors (Hicks et al., Reference Hicks, South, DiRago, Iacono and McGue2009), which will provide an opportunity to explore genotype–environment interaction effects as findings emerge from GWAS consortia.
The MCTFR sample will also soon be genotyped on Illumina's HumanExome BeadChip, which includes more than 240,000 coding sequence variants. These data, when combined with the original GWAS data, will provide extensive coverage of both common and relatively rare genetic variants in key genetic regions. We are also at the initial stages of planning for whole genome sequencing a large number of MCTFR participants, as well as obtaining additional relevant phenotype information. In this era of large-scale consortia, the extensive phenotypic, genetic, and environmental assessments in the MCTFR will provide numerous opportunities for collaboration.
Acknowledgment
This research was supported in part by USPHS Grants from the National Institute on Alcohol Abuse and Alcoholism (AA09367 and AA11886), the National Institute on Drug Abuse (DA05147, DA13240, and DA024417), and the National Institute on Mental Health (MH066140).
