Overview of Association Testing with Known Structure

Consider the problem of testing for association of a genetic marker with a particular phenotype in a sample of n genotyped individuals. For simplicity, we assume that the marker to be tested is a SNP, with alleles labeled “0” and “1.” (The extension to multiallelic markers can be obtained as in previous work.¹⁶) Let Y = (Y₁, … Y_n)^T where Y_i = ½ × (the number of alleles of type 1 in individual i), so the value of Y_i is 0, ½, or 1. We treat the genotype data on the n individuals as random and the available phenotype information as fixed in the analysis, an approach that is appropriate, for example, with either random or phenotype-based ascertainment. Under the null hypothesis of no association and no linkage between marker and trait, we assume that the expected value of Y is E₀(Y) = p1, where 1 is a column vector of 1s of length n and p is a parameter representing the frequency of the type 1 allele. In models incorporating population structure, p would typically be interpreted as an “ancestral” allele frequency or some kind of average allele frequency across subpopulations. We denote by Var₀(Y) = Σ the n × n covariance matrix of Y under the null hypothesis of no association. It is often convenient to write Σ = σ²Ψ, where σ² is defined to be the variance of Y for an outbred individual in the absence of population structure, and Ψ accounts for relatedness, inbreeding, and population structure. We use the term “known structure” to refer to the case when the matrix Ψ is known. We always take σ² to be unknown and estimate it from the sample. Denote by $\widehat{\sigma }$² a suitable estimator of σ² (where two examples of suitable estimators are given in the next subsection). Then, in the case of known structure, we consider test statistics for association that have the rather general form

where V is a fixed, nonzero column vector of length n such that V^T1 = 0. Note that Var₀(V^TY) = σ²V^TΨV, so the denominator in Equation 1 can be viewed as an estimator of Var₀(V^TY). In a test for association, V would naturally include phenotype information and could also include pedigree information. One could include covariate information in V as well, although we do not treat that situation in the present work. There are a number of case-control association test statistics that have the general form in Equation 1, including the Pearson χ² test, the Armitage trend test,²⁴ the corrected χ² test,¹⁵ the W_QLS test,¹⁵ and the M_QLS test¹⁶ (details on how these tests can be written in the form of Equation 1 are given in subsection Examples of Association Tests with Known Structure). Under standard regularity conditions, the test statistic given in Equation 1 has an asymptotic χ₁² distribution under the null hypothesis of no association and no linkage.

Estimation of σ² when Structure Is Known

In the context of Equation 1, when structure is known, we have two general approaches for estimating σ² under the null hypothesis. If we assume that, for an outbred individual in the absence of population structure, HWE holds at the marker, then σ² = ½p(1 − p), where p is the frequency of allele 1 at the SNP being tested, and a reasonable estimator of σ² under this assumption is ${\widehat{\sigma }}_1^2=0.5\widehat{p}\left(1-\widehat{p}\right)$, where $\widehat{p}$ is a suitable estimator of p, the frequency of allele 1 at the SNP being tested. Examples of suitable estimators of p are (1) the sample frequency, $\overline{Y}$; (2) the best linear unbiased estimator (BLUE),²⁵ given by

and (3) a Bayesian estimator²⁶ such as $\left(n\overline{Y}+0.5\right)/\left(n+1\right)$.

Alternatively, an approach to estimation of σ² that does not assume σ² = 0.5p(1 − p) could be used. When Ψ is known, a reasonable estimator is

which is RSS/(n − 1) for generalized regression of Y on 1, where RSS is the residual sum of squares. Note that when Ψ = I, the n × n identity matrix, e.g., with unrelated individuals in the absence of population structure, then ${\widehat{\sigma }}_2^2$ is just the sample variance of Y.

Examples of Association Tests with Known Structure

Outline of ROADTRIPS Approach for Unknown Structure

The idea behind ROADTRIPS is to extend tests of the form given in Equation 1 to the situation when there could be unknown population structure and/or cryptic relatedness in the sample. To do this, we use genome-screen data to form an appropriate estimator $\widehat{\mathbf{\Psi }}$ of Ψ and consider various tests of the form

where we can allow V to take into account any known pedigree information, in addition to phenotype information, while simultaneously accounting for pedigree errors and additional unknown structure through $\widehat{\mathbf{\Psi }}$. This approach allows us to easily adapt to different patterns of missing genotypes at different tested markers, by including only the rows and columns of $\widehat{\mathbf{\Psi }}$ (and the entries of V and Y) that correspond to the individuals genotyped at the particular marker being tested. In what follows, we first describe the population genetic modeling assumptions that underlie our estimation and testing procedures, then we describe the estimators $\widehat{\mathbf{\Psi }}$ and ${\widehat{\sigma }}^2$.

Population Genetic Modeling Assumptions

The modeling assumptions we make are weak and are satisfied by commonly used models of population structure and commonly used models for related individuals. We consider S SNPs in a genome screen, and we generalize the notation of the previous subsections to a set of S SNPs by letting Y^s be the genotype vector corresponding to SNP s, namely, Y^s = (Y₁^s, …, Y_n^s)^T, s = 1, …, S, where Y^s_i = ½ × (the number of alleles of type 1 at SNP s in individual i). Our modeling assumption on the null mean, generalized from the preceding subsections, can be stated as

We make the following assumption regarding the null covariance matrix:

where Ψ is an arbitrary, positive semidefinite matrix, and σ_s² > 0 for all 1 ≤ s ≤ S. Here, the key point is that the correlation structure, captured by Ψ, is assumed to be the same across SNPs, whereas the scalar multiplier σ_s² is allowed to vary across SNPs. (Of course, this presumes that the same individuals are genotyped at all SNPs. When some individuals have missing genotypes at SNP s, the entries of Y^s and the rows and columns of Ψ that correspond to individuals with missing genotypes would be deleted.)

Note that Ψ and σ_s² are defined only up to a constant multiple, in the sense that cΨ and c⁻¹σ_s² would give the same value of Σ_s. By convention, σ_s² is usually chosen to be the variance of an outbred individual in the absence of population structure. We now give two examples of population genetic models that satisfy our assumptions.

Example 1: Related Individuals without Additional Population Structure

An example of a simple model that satisfies the assumptions of the previous subsection is the model for related individuals in an unstructured population. In this model, individuals in the sample can be related by pedigrees, where the pedigree founders are assumed to be independently drawn from a population that is in Hardy-Weinberg equilibrium (HWE). Mendelian inheritance is assumed in the pedigrees. In this case, it has previously been shown¹⁵ that Equations 9 and 10 hold, where p_s is interpreted as the allele frequency of SNP s in the population from which the founders are drawn, ${\sigma }_s^2=0.5p_s\left(1-p_s\right)$, and Ψ = Φ, where Φ is the kinship matrix given in Equation 4.

If the pedigrees are fully known, then the structure matrix Ψ is known, but if some genealogical information is missing, then Ψ might be partially or completely unknown.

Example 2: Balding-Nichols Model with Admixture

In the Balding-Nichols model^6,28,29 with admixture, we let p_s denote the “ancestral” allele frequency at SNP s and let q_k^s denote the allele frequency of SNP s in subpopulation k, 1 ≤ k ≤ K. We assume that the q_k^s are random variables that are independent across both k and s, with q_k^s ~Beta(p_s(1 − f_k)/f_k, (1 − p_s)(1 − f_k)/f_k), where f_k ≥ 0 can be viewed as Wright's standardized measure of variation³⁰ for subpopulation k. For SNP s, let q^s = (q₁^s, …, q_K^s)^T denote the vector of subpopulation-specific allele frequencies. Individual i is assumed to have admixture vector a_i = (a_i1, …, a_iK)^T, where a_ik ≥ 0 for all i and k, and ${\sum }_{k=1}^K{a}_{ik}=1$ for all i. Conditional on the random variable q^s, the two alleles of individual i at SNP s are assumed to be independent, identically distributed (i.i.d.) Bernoulli(a_i^Tq^s) random variables. In the Balding-Nichols model with admixture, Equations 9 and 10 hold, where p_s is interpreted as the “ancestral” allele frequency, ${\sigma }_s^2=0.5p_s\left(1-p_s\right)$, and the entries of Ψ are given by ${\mathrm{\Psi }}_{ii}=1+{\sum }_{k=1}^K{a}_{ik}^2{f}_k$ and ${\mathrm{\Psi }}_{ij}=2{\sum }_{k=1}^K{a}_{ik}{a}_{jk}{f}_k$ if i ≠ j. In this context, the population structure captured by Ψ would typically be unknown.

Estimation of the Matrix Ψ

The matrix Ψ is a function of the genealogy of the sampled individuals, where genealogy is broadly interpreted as including both population structure and the pedigree relationships of close relatives. The matrix Ψ will be unknown when there is hidden population structure and/or cryptic relatedness in the sample. We allow a completely general form for Ψ, assuming only that it is positive semidefinite (psd). When genome-screen data are available on the sampled individuals, this information can be used to estimate Ψ. For any pair of individuals i and j, let 𝒮_ij be the set of markers for which both i and j have nonmissing genotype data. Then if the allele frequencies p_s were known, and if σ_s² = 0.5p_s(1 − p_s) as in Examples 1 and 2, an unbiased estimator of Ψ_ij would be

where |𝒮_ij| is the number of elements of 𝒮_ij. If one assumed, for example, that genotypes at different SNPs were independent with $\left|{\mathcal{S}}_{ij}\right|\to \mathrm{\infty }$ and p_s known and that σ_s² = 0.5p_s(1 − p_s) held at all but a finite number of SNPs, then Equation 11 would provide a consistent estimator of Ψ_ij. However, p_s will generally not be known, so we propose to further restrict 𝒮_ij to those markers that are polymorphic in the sample; let ${\widehat{p}}_s={\overline{Y}}^s$, the observed proportion of type 1 alleles in the sample at marker s; and use estimator

The estimator $\widehat{\mathbf{\Psi }}$ of Equation 12 is essentially the same as the estimated covariance matrix used in EIGENSTRAT.⁶ An alternative estimator could be obtained by using the sample variance of Y^s in the denominator instead of $0.5{\widehat{p}}_s\left(1-{\widehat{p}}_s\right)$.

If every sampled individual were genotyped at the same markers, with no missing genotypes, then $\widehat{\mathbf{\Psi }}$ would be psd and singular, with $\widehat{\mathbf{\Psi }}\mathbf{1}=0$, i.e., 1 would be in the null space of $\widehat{\mathbf{\Psi }}$. With missing genotypes, it is possible for $\widehat{\mathbf{\Psi }}$ to be nonsingular and non-psd. The fact that $\widehat{\mathbf{\Psi }}$ might be non-psd is not, in itself, particularly problematic from a practical point of view, provided ${\mathbf{V}}^T\widehat{\mathbf{\Psi }}\mathbf{V}>0$ for the chosen V in Equation 8 and assuming this provides a sufficiently accurate estimator of Var₀(V^TY^s)/σ²_s. The fact that $\widehat{\mathbf{\Psi }}$ might be singular (e.g., in the case of no missing data) or close to singular, with $\widehat{\mathbf{\Psi }}$ orthogonal or approximately orthogonal to the vector 1, means that in those cases, one would not be able to directly plug $\widehat{\mathbf{\Psi }}$ into formulae such as Equations 2, 3, 6, and 7. This is discussed further in the next subsection. With substantially different amounts of missing data at different markers, as in the RA and COGA data sets analyzed in the Results section, the matrix $\widehat{\mathbf{\Psi }}$ might be nonsingular and so could be directly used in Equations 2, 3, 6, and 7.

Estimation of σ² when Structure Is Unknown

In this subsection, we drop the subscript s and use notations Y, p, and σ² for the SNP being tested; e.g., we assume E₀(Y) = p1 and Var₀(Y) = Σ = σ²Ψ. As we did for the case of known structure, we consider two general approaches to estimation of σ² when structure is unknown. The first approach is to take estimators of the form ${\widehat{\sigma }}_1^2=0.5\widehat{p}\left(1-\widehat{p}\right)$, where $\widehat{p}$ is a suitable estimator of p. When $\widehat{\mathbf{\Psi }}$ is orthogonal or approximately orthogonal to the vector 1, we cannot plug it into Equation 2 to obtain the BLUE of p. (Note that use of the Moore-Penrose generalized inverse ${\widehat{\mathbf{\Psi }}}^{-}$ in place of ${\widehat{\mathbf{\Psi }}}^{-1}$ also does not work, because ${\widehat{\mathbf{\Psi }}}^{-}$ is also orthogonal or approximately orthogonal to 1, so plugging into Equation 2 would result in both numerator and denominator being exactly or approximately zero.) Instead, we use the more stable estimator $\widehat{p}=\overline{Y}$, the sample allele frequency. Thus, our first estimator becomes

As we did in the case of known structure, we also consider an estimator of σ² that does not assume σ² = 0.5p(1 − p) at the SNP being tested. When $\widehat{\mathbf{\Psi }}$ is orthogonal or approximately orthogonal to the vector 1, we replace Equation 3 by

where ${\widehat{\mathbf{\Psi }}}^{-}$ is the Moore-Penrose generalized inverse of $\widehat{\mathbf{\Psi }}$.

Association Tests when Structure Is Partially or Completely Unknown

We apply Equation 8 to extend association tests developed for the situation of known structure (as described in subsection Examples of Association Tests with Known Structure) to association tests that are appropriate for situations of partially or completely unknown structure. We call the tests based on Equation 8 the ROADTRIPS versions of the corresponding tests given by Equation 1, and we now give several examples. Table 1 gives the weight vector for each statistic defined below.

GAW 15 Rheumatoid Arthritis Data

We apply ROADTRIPS to perform association analysis of RA data provided for GAW 15 by a UK group led by Jane Worthington and Sally John (these data are described in detail elsewhere).²² Data are available on 157 nuclear families, where 156 of these have at least two affected individuals. Individuals were diagnosed as affected according to the American College of Rheumatology (ACR) criteria. There are 550 individuals with available genotype data. After exclusion of 2 duplicate individuals and 4 outlier individuals who have estimated inbreeding coefficients more than 3 standard deviations (SDs) above the average (where the estimated inbreeding coefficient of individual i is taken to be ${\widehat{\mathrm{\Psi }}}_{ii}-1$), there are 339 affected individuals, 198 unaffected controls, and 7 controls of unknown phenotype in the analysis. The data set includes 10,156 autosomal SNPs that passed quality control filters. We exclude 285 SNPs that are not polymorphic (minor allele frequency less than 0.01). The remaining 9871 SNPs are tested for association.

Relationship Configurations

Both relationship configurations 1 and 2 include 100 unrelated affected individuals, 400 unrelated unaffected individuals, and individuals sampled from 120 outbred, three-generation pedigrees, where each pedigree has a total of 16 individuals, the phenotypes of all individuals in the pedigrees are observed, and genotypes are observed for only a subset of individuals in each pedigree. We sample six types of these pedigrees; the types are described in Table 2. Relationship configuration 1 has 40 pedigrees of type 1, 40 of type 2, and 40 of type 3, as well as 100 unrelated affected and 400 unrelated unaffected individuals. Relationship configuration 2 contains all six types of pedigrees in Table 2, with 20 of each type, and it also contains 100 unrelated affected and 400 unrelated unaffected individuals.

Population Structure Settings

Each simulation setting specifies a particular relationship configuration combined with a particular setting of population structure. Each setting of population structure is a special case of the Balding-Nichols model with admixture described in Example 2, in which we take f_k = 0.01 for every subpopulation. Population structure 1 has individuals sampled from two subpopulations, with 60% of the pedigrees and affected unrelated individuals sampled from subpopulation 1 and the remaining 40% of the pedigrees and affected unrelated individuals sampled from subpopulation 2. Among the unrelated unaffecteds, 40% are sampled from subpopulation 1 and 60% from subpopulation 2. Population structure 2 is similar to population structure 1, except that the proportions 60% and 40% are replaced by 80% and 20%, respectively. Population structure 3 is similar to population structure 1, except that there are three subpopulations, with all of the unrelated unaffecteds sampled from subpopulation 3. Population structure 4 has individuals sampled from an admixed population, formed from two subpopulations. Individuals in the admixed population are assumed to have i.i.d. admixture vectors of the form (a, 1 − a), where a is a Uniform(0,1) random variable.

Within a given setting of population structure, to simulate the unrelated individuals needed in relationship configurations 1 and 2, we first simulate genotypes according to the chosen setting of population structure, simulate phenotypes conditional on genotypes, and then randomly ascertain 100 affected and 400 unaffected individuals. To sample particular pedigree types within a given setting of population structure, we first simulate genotypes for pedigree founders according to the chosen setting of population structure, drop alleles down the pedigree, simulate phenotypes conditional on genotypes, and then do rejection sampling to obtain 100 replicates of each of the pedigree configuration types. Then, in the simulations, pedigrees of each type are sampled with replacement from the 100 previously obtained replicates of that pedigree type.

Random and Causal SNPs

We use a trait model that has two unlinked causal SNPs (which we call “SNP 1” and “SNP 2”) with epistasis between them.¹⁶ The ancestral frequencies of the type 1 alleles at SNPs 1 and 2 are taken to be 0.1 and 0.5, respectively. Individuals with at least one copy of allele 1 at SNP 1 and at least one copy of allele 1 at SNP 2 have a penetrance of 0.3. All other individuals have a penetrance of 0.05. In the power studies, association is tested with SNP 2. In contrast, “random” SNPs are assumed to be unlinked and unassociated with the trait, and their ancestral allele frequencies are obtained as i.i.d. draws from a uniform (0.1, 0.9) distribution.

Assessment of Type 1 Error

For each of the eight combinations of population structure and relationship configuration, we generate genotype data for 100,000 random SNPs that are neither linked nor associated with the trait, and we test each of them for association at the 0.0001 level, using various test statistics. In Table 3, for each test statistic, we report the empirical type 1 error, which we calculate as the proportion of simulations in which the test statistic exceeds the χ₁² quantile corresponding to nominal type 1 error level 0.0001. The statistics compared are the four ROADTRIPS statistics, Rχ₂, RM₂, RM_NI2, and RW_NI2, where the subscript “2” denotes use of the estimator ${\widehat{\sigma }}_2^2$ of Equation 14; GC; EIGENSTRAT; the M_QLS; the corrected Armitage χ²; and the uncorrected Armitage trend test. Note that because there are only two types of controls, the Rχ₂ and RM_NI2 tests are identical. The method of Rakovski and Stram¹³ corresponds to the ROADTRIPS statistic RW_NI2, which is included in the comparison. The correct type 1 error of FBAT has been established previously.³¹ Using an exact binomial calculation, we determine that empirical type 1 error rates falling in the range of 0.00004–0.00016 are not significantly different from the nominal 0.0001 level.

For GC and all of the ROADTRIPS statistics, empirical type 1 error is not significantly different from the nominal level. In contrast, type 1 error is inflated for EIGENSTRAT, M_QLS, the corrected Armitage χ², and the Armitage trend test. This is to be expected, because these tests either (1) correct for related individuals but not for population structure (M_QLS, corrected Armitage χ²), (2) correct for population structure but not for related individuals (EIGENSTRAT), or (3) correct for neither (Armitage trend test). In particular, the top principal components in EIGENSTRAT are not able to capture the complicated covariance structure due to the related individuals in the samples. The results for EIGENSTRAT in Table 3 are obtained with the default setting of ten principal components and with outlier removal. The results without outlier removal and with different numbers of principal components are similar (results not shown). We also performed all the ROADTRIPS tests with variance estimator ${\widehat{\sigma }}_1^2$ of Equation 13 instead of ${\widehat{\sigma }}_2^2$ and obtained nearly identical empirical type 1 error rates (results not shown).

Type 1 Error with Missing Genotypes: Rχ₂ and GC

When all SNPs have the same pattern of missing genotypes, we expect Rχ₂ to perform similarly to GC, because in this case, both Rχ₂ and GC are equivalent to correcting all the Armitage χ² statistics across the genome by a common factor, though this factor differs between the two methods. However, when different SNPs have different rates of missing genotypes, we expect the Rχ₂ statistic to have better control of type 1 error than GC, because the Rχ₂ statistic allows different SNPs to have different correction factors appropriate to the level of genotype information available, whereas GC applies the same correction factor to all SNPs.

To assess the magnitude of this effect, we perform a simulation study under four different settings of population structure and relationship configuration, given in columns 1 and 2 of Table 4. For each combination of settings, genotype data are generated for 100,000 random SNPs that are neither linked nor associated with the phenotype. The proportions of individuals with missing genotype data at different SNPs are taken to be i.i.d. random variables drawn from a Beta(3, 12) distribution, which has a mean of 0.2 and a SD of 0.1. Given the proportion of missing genotypes at a marker, the individuals whose genotypes will be set to missing for that marker are chosen uniformly at random from the sample.

The empirical type 1 error rates for Rχ₂ and GC are given in Table 4. GC has inflated type 1 error for all of the simulation settings, because of undercorrection of test statistics from SNPs that have relatively low levels of missing genotypes, whereas the empirical type 1 error for Rχ₂ is not significantly different from the nominal 0.0001 level. The results illustrate that ROADTRIPS is not only robust to cryptic population and pedigree structure, but also to varying rates of randomly missing genotype data. This results from the fact that the entire empirical covariance matrix is estimated in ROADTRIPS, so when some individuals have missing genotype data at the SNP being tested, the corresponding rows and columns can be deleted from the empirical covariance matrix, allowing one to obtain a variance estimator that accounts for missing genotype data at the marker being tested.

Power Comparison

We assess power to detect association in the presence of population and pedigree structure only for those tests that maintain correct nominal type 1 error, namely the four ROADTRIPS tests, FBAT, and GC (although, as can be seen in Table 4, the type 1 error of GC might not be correct when the rate of missing genotype data varies across markers). The simulations are performed with relationship configuration 2 under each of the four different settings of population structure. For each setting, 5000 simulated replicates are performed, and SNP 2 of the trait model is tested for association with the trait by each method. Table 5 reports power for each statistic for settings 1, 2, and 4 of population structure. Here, power is calculated as the proportion of simulations for which the statistic exceeds the χ₁² quantile corresponding to nominal type 1 error level 0.0001. Power for population structure 3 is close to 0 for all of the statistics (data not shown). The ROADTRIPS statistics are all calculated with estimator ${\widehat{\sigma }}_2^2$ of Equation 14. The FBAT and RM methods are given the information of the correct pedigree structure, but not the population structure. The FBAT statistic is calculated with offset value set equal to the prevalence, and the RM statistic is calculated with k set equal to the prevalence. As expected, the RM test is the most powerful in all settings, because it uses the known pedigree information to incorporate phenotype information about relatives with missing genotype data, and it corrects for additional unknown population structure by means of the empirical covariance matrix. In contrast, the FBAT test, which was given all the same information as the RM test, performs very poorly in these simulations, because it is not able to incorporate the data on the 500 unrelated individuals, and it is also not able to incorporate the data from pedigree types 1, 2, and 3, because they do not meet the FBAT criteria for “informative families.” If we assume that no pedigree information is available, then the Rχ, RM_NI, and GC tests all give identical or almost identical power, assuming that all markers have comparable amounts of missing genotype data. When different SNPs have different amounts of missing genotype data, GC might not adequately control type 1 error, so Rχ and RM_NI would be preferred. The RW_NI test, which corresponds to the test of Rakovski and Stram,¹³ has lower power than Rχ and RM_NI, which is not surprising in light of the fact that the W_QLS test, of which the RW_NI is an extension, was shown¹⁶ to have generally lower power than the M_QLS and corrected χ² tests, of which the RM_NI and Rχ, respectively, are extensions.

We also assess power to detect association when there is pedigree structure but no population structure. The simulation is carried out as in the preceding paragraph except that relationship configuration 2 with no population structure is simulated. In addition to the tests compared in Table 5, we also calculate power for the M_QLS and corrected Armitage χ² tests, both of which have correct type 1 error in this setting, provided that the pedigree structure is known. Estimated power for this setting is given in Table 6. The M_QLS and RM tests are the most powerful among the tests considered. The RM test is able to match the power of M_QLS even though RM uses the estimator $\widehat{\mathbf{\Psi }}$ in the variance calculation, whereas M_QLS uses the true Ψ. As in the previous power comparison, FBAT has almost no power in this setting, and the RW_NI test has lower power than all of the other statistics except FBAT. There is no significant difference in power among Rχ, RM_NI, GC, and the corrected Armitage χ² tests. The M_QLS and the corrected Armitage χ² tests have power identical to their corresponding ROADTRIPS extensions, RM and Rχ, which illustrates that power is not compromised by use of the empirical matrix $\widehat{\mathbf{\Psi }}$ in the variance correction for the ROADTRIPS statistics.

GAW 15 Rheumatoid Arthritis Data

We apply RM, Rχ, GC, and FBAT to the GAW 15 data to test for association of SNPs with RA. Using a previously reported²² prevalence of 0.8% for RA in people of European descent, we set both the offset value in FBAT and the prevalence value in RM to 0.008. The entries of the empirical covariance matrix $\widehat{\mathbf{\Psi }}$ and the correction factor for GC are calculated using SNP data across the autosomal chromosomes for the study individuals. Table 7 gives the results of all tests for those SNPs for which at least one of the tests has a nominal p value <5.0 × 10⁻⁵. For 8 of these 10 SNPs, the RM test has the smallest p value among the four tests used. After Bonferroni correction to adjust for four different tests of association at each of 9,871 SNPs, the RM test is significant at the 5% level for 2 SNPs: snp264363 on chromosome 11 (p = 5.2 × 10⁻⁷ uncorrected, 0.021 corrected) and snp152076 on chromosome 9 (p = 5.4 × 10⁻⁷ uncorrected, 0.021 corrected). The Rχ test is significant at the 5% level for an additional SNP, snp547632 on chromosome 11 (p = 1.2 × 10⁻⁶ uncorrected, 0.047 corrected).

A histogram of the estimated self-kinship values (where the estimated self-kinship value of individual i is taken to be $.5{\widehat{\mathrm{\Psi }}}_{i,i}$) of the individuals included in the analysis can be found in Figure 1. The histogram shows that the values are not centered around 0.5, which is the self-kinship value in the absence of population structure and inbreeding. The self-kinship mean is 0.512, and the majority of the kinship values (77%) are greater than 0.5. There are a few pairs of individuals that should have a kinship coefficient value equal to 0.25 based on the available pedigree information (i.e., parent-offspring pairs and sibling pairs), but have estimated kinship coefficient values close to 0 and thus appear to be unrelated. There are also a few pairs that are not members of the same pedigree but have kinship coefficient estimates that indicate that they are related. Both these phenomena could be caused by sample switches. An attractive feature of the ROADTRIPS methods is that they automatically correct for misspecified relationships in a sample, in addition to hidden population structure, while simultaneously allowing for different SNPs to have different rates of missing genotypes.

Open in a separate window

Figure 1

Histogram of Estimated Self-Kinship Coefficient Values for the GAW 15 UK RA Data

The vertical line at 0.5 represents the self-kinship value in the absence of population structure and inbreeding.

GAW 14 COGA Data

We apply RM, Rχ, and GC to test for association with an alcoholism-related phenotype in the GAW 14 COGA data. A previous analysis¹⁶ of these data used the M_QLS (with k set to 0.05) and corrected χ² tests; these tests correct for known pedigree information, but do not make any correction for unknown structure. In our reanalysis, we compare the results of the tests with and without the correction for unknown structure. To make the results comparable, we set k = 0.05 in the RM test. Table 8 lists SNPs for which at least one of the RM, Rχ, and GC tests has a nominal p value <1.0 × 10⁻⁵. For 11 of these 12 SNPs, the RM test has the smallest p value among the three tests used. After Bonferroni correction to adjust for three different tests of association at each of 10,407 SNPs, the RM test is significant at the 5% level for 4 SNPs: tsc1750530 on chromosome 16 (p = 2.3 × 10⁻⁸ uncorrected, 0.0007 corrected), tsc0046696 on chromosome 18 (p = 1.4 × 10⁻⁷ uncorrected, 0.004 corrected), tsc1177811 on chromosome 1 (p = 1.5 × 10⁻⁷ uncorrected, 0.005 corrected), and tsc1637642 on chromosome 5 (p = 6.8 × 10⁻⁷ uncorrected, 0.02 corrected). The Rχ test is significant at the 5% level for an additional SNP, tsc0571038 on chromosome 11 (p = 8.0 × 10⁻⁷ uncorrected, 0.025 corrected). Of the 5 SNPs that were previously¹⁶ identified as genome-wide significant by M_QLS or $W_{{\chi }_{corr}^2}$, 4 were also identified by ROADTRIPS. The exception is tsc0057290 on chromosome 18, which was identified by M_QLS and is no longer significant after correcting for cryptic structure. ROADTRIPS identifies an additional SNP, tsc1637642 on chromosome 5, which was not identified in the previous analysis. GC did not identify any significant SNPs. There are some SNPs in Table 8 that are in or near genes of interest; the details have previously been reported.¹⁶ A previous analysis³² of these data with FBAT, in which a slightly larger sample of individuals was used, detected one SNP with nominal p value 6 × 10⁻⁵, which is not significant after Bonferroni correction for the number of SNPs tested.

Figure 2 gives a histogram of the estimated self-kinship coefficient values for the genotyped individuals who were included in the analysis. The center of the histogram is shifted from 0.5 (the self-kinship coefficient in the absence of population structure and inbreeding). Seventy-one percent of the values are greater than 0.5, and the mean self-kinship value is 0.506. Just as there were in the UK RA data, in the COGA data there are a few pairs of individuals who appear to have misspecified relationships or be cryptically related. As previously mentioned, ROADTRIPS adjusts for this in the variance correction.

Open in a separate window

Figure 2

Histogram of Estimated Self-Kinship Coefficient Values for the GAW 14 COGA Data

The vertical line at 0.5 represents the self-kinship value in the absence of population structure and inbreeding.

This study was supported in part by the University of California President's Postdoctoral Fellowship and the Lamond Family Foundation Postdoctoral Fellowship (to T.T.) and by National Institutes of Health grant R01 HG001645 (to M.S.M.). Data on alcohol dependence were provided by the Collaborative Study on the Genetics of Alcoholism (U10AA008401) through the Genetics Analysis Workshop (R01GM031575). Data on RA were provided by a UK group led by Jane Worthington and Sally John through the Genetics Analysis Workshop (R01GM031575).