Detecting genes contributing to longevity using twin data

Abstract

Searching for genes contributing to longevity is a typical task in association analysis. A number of methods can be used for finding this association -- from the simplest method based on the technique of contingency tables to more complex algorithms involving demographic data, which allow us to estimate the genotype-specific hazard functions. The independence of individuals is the common assumption in all these methods. At the same time, data on related individuals such as twins are often used in genetic studies. This paper proposes an extension of the relative risk model to encompass twin data. We estimate the power and also discuss what happens if we treat the twin data using the univariate model.

Introduction

Most common diseases and traits have a complex structure, for which the phenotype is determined by interactions between genetic and environmental factors. As any individual genetic variant can have a relatively modest effect on a disease or trait, linkage analysis has less power than association analysis. Classical association studies in their simplest form compare the frequency of alleles or genotypes for candidate genes between cases and controls. These candidate genes are usually chosen on the basis of biological hypotheses or from previous linkage analyses.

To identify genes associated with longevity, information on genotype frequencies for two or more age groups is needed. A significant trend of genotype frequencies being associated with age can indicate a gene-longevity association. In the basic 'gene frequency method', only the genotype frequencies in different age groups are compared [1–3] Some extensions of this method involve the use of demographic information about the population under study and allow the estimation of initial frequencies, relative risks and the age trajectories of mortality for candidate genes. These methods are known in the literature as the 'parametric method', the 'semi-parametric method', the 'non-parametric method' and the 'relative risk method' [4]. The use of these methods, however, has two limitations. First, the initial gene frequencies in all cohorts represented in the study must be the same. Secondly, the mortalities for genotypes do not depend on the birth year of the cohort. In two recently published papers,[5, 6] the authors exclude the first limitation, assuming a time trend in the genetic frequencies of subsequent birth cohorts. In principle, the time and the cohort covariates influencing mortality can be incorporated into the models too. The flexible parameterisation in the extended relative risk model [6] also allows detection of the antagonistic pleiotropic effect.

The methods mentioned above have been developed for datasets consisting of independent individuals. In this paper, we propose a method for detecting longevity genes for the dataset consisting of twin pairs. This method retains all the advantages of the relative risk model for univariate data described previously [6].

Materials and methods

To analyse the gene-longevity association, two datasets are needed: the genotype data and the univariate survival data for the individuals involved in the study. To improve the accuracy and power of the study, the longevity data for twins can additionally be analysed. Denoting longevity and non-longevity alleles at an autosomal locus by a and A, respectively, assume that the frequencies P _g of genotypes AA, aA or Aa, and aa at the moment of birth are P _AA , P _Aa and P _aa , respectively. If the Hardy-Weinberg equilibrium holds, then P _AA = (1-P _a )², P _Aa = 2P _a (1-P _a ) and P_{aa}=P_a^2, where P _a is the frequency of the allele a at the moment of birth. We parameterise P _a as follows:

P_a=1-1/\left(1+e^{v+\delta x+R\phi \left(x,x_0\right)}\right),x=T-t.

(1)

In accordance with (1), the logit of P _a is a linear function of unknown parameters R, ν and δ with domain of definition R³. This parameterisation includes the sudden change in the allele frequency by the value Rφ(x, x₀) in the cohort T - x₀ and the slow linear cohort effect ν + δx of the allele frequency. Here, T stands for the year of data collection, x for the age, and t for the cohort. We assume that the value of x₀ is known. The step function φ(x, x₀) is defined by the interval equations φ(x, x₀) = 1 for 0 ≤ x ≤ x₀ and φ(x, x₀) = 0 for x >x₀.

To estimate the genotype frequencies for twin pairs, we need to calculate the bivariate survival functions. One possible approach to doing so is to use the correlated gamma-frailty model, which provides simple analytical expressions for the bivariate survival functions [7]. Assume that that individual's instantaneous risk of death μ for genotype g∈{aa, Aa, AA} at age x, as measured by the hazard of mortality, is μ(x, Z, g) = Zμ_{0, g}(x), where Z is the gamma distributed frailty (non-observed risk of mortality) with mean 1 and variance σ², and μ_{0, g}(x) is the baseline hazard. The univariate survival function S_g\left(x\right)=E{e}^{-Z{H}_g\left(x\right)}={\left(1+{\sigma }^2{H}_g\left(x\right)\right)}^{-1/{\sigma }^2} is the Laplace transform for the gamma probability density function at the point H_g\left(x\right)={\int }_0^x{\mu }_{0,g}\left(t\right)dt (cumulative hazard function). For related individuals, we assume that life spans T₁ and T₂ are conditionally independent, given frailties Z₁, Z₂ and genotypes g₁, g₂. In general, frailties Z₁ and Z₂ have unequal variances. Below, we shall assume, for simplicity, that Z₁ and Z₂ are identically distributed. If Corr(Z₁, Z₂) = ρ, E(Z₁) = E(Z₂) = 1 and Var(Z₁) = Var(Z₂) = σ², then:

\begin{array}{c}P\left\{T_1>x_1,T_2>x_2\right\}=S_{g1,g2}\left(x_1,x_2\right)\hfill \\ =\frac{S_{g1}{\left(x_1\right)}^{1-\rho }{S}_{g2}{\left(x_2\right)}^{1-\rho }}{{\left(S_{g1}{\left(x_1\right)}^{-{\sigma }^2}+S_{g2}{\left(x_2\right)}^{-{\sigma }^2}-1\right)}^{\rho /{\sigma }^2}}\hfill \end{array}

(2)

Here, S_{g1, g2}(x₁, x₂) is the bivariate survival function at ages x₁ and x₂ for twins with genotypes g₁ and g₂, respectively. We relate cumulative hazard functions with some unknown function H₀(x) as follows:

H_g\left(x\right)=c_{g}x+a_g{H}_0{\left(x\right)}^{b_g}

(3)

with unknown a _g ≥ 0, b _g ≥ 0 and c _g ≥ 0. Such parameterisation, where cumulative hazards H _g (x) rather than survival functions S _g (x) for different genotypes are parametrically related (eg S_g\left(x\right)=S_0{\left(x\right)}^{b_g}, is more flexible and allows us to detect the antagonistic pleiotropic effect [6]. Without loss of generality, we can assume that a _AA = b _AA = 1.

For univariate and bivariate survival functions in the whole population, it holds that:

\begin{array}{c}S\left(x\right)=\sum _g{P}_g{S}_g\left(x\right),\hfill \\ S^{MZ}\left(x_1,x_2\right)=\sum _{g,g}{P}_{g,g}^{MZ}{S}_{g,g}^{MZ}\left(x_1,x_2\right)\hfill \\ S^{DZ}\left(x_1,x_2\right)=\sum _{g1,g2}{P}_{g1,g2}^{DZ}{S}_{g1,g2}^{DZ}\left(x_1,x_2\right)\hfill \end{array}

(4)

Here, P _g , P_{g,g}^{MZ} and P_{g1,g2}^{DZ} are the univariate and the bivariate genotype frequencies for monozygotic (MZ) and dizygotic (DZ) twin pairs, respectively, at the moment of birth. Since the frailty correlation ρ can be different for MZ and DZ twins, we use the upper index MZ or DZ in the notation for bivariate survival. Given univariate survival S(x) and parameters, we calculate the baseline cumulative hazard H₀(x) using the simple bisectional procedure [6]. For univariate genotype frequencies, we will use the values given above. To calculate the bivariate genotype frequencies, note that for MZ twin pairs, g₁ = g₂ = g and P_{g,g}^{MZ}=P_g. Assuming independent transmission of the maternal and paternal alleles to the offspring, we can use the standard formulae for DZ twin pairs:

\begin{array}{c}{P}_{aa,aa}^{DZ}=P_{aa}^2+\left(1/2\right)P_{aa}{P}_{Aa}+\left(1/16\right)P_{Aa}^2\hfill \\ P_{aa,Aa}^{DZ}=P_{Aa,aa}^{DZ}=\left(1/2\right)P_{aa}{P}_{Aa}+\left(1/8\right)P_{Aa}^2\hfill \\ P_{aa,AA}^{DZ}=P_{AA,aa}^{DZ}=\left(1/16\right)P_{Aa}^2\hfill \\ P_{Aa,AA}^{DZ}=P_{AA,Aa}^{DZ}=\left(1/2\right)P_{Aa}{P}_{AA}+\left(1/8\right)P_{Aa}^2\hfill \\ P_{Aa,Aa}^{DZ}=\left(1/2\right)P_{aa}{P}_{Aa}+2P_{aa}{P}_{AA}\hfill \\ +\left(1/2\right)P_{Aa}{P}_{AA}+\left(1/4\right)P_{Aa}^2\hfill \\ P_{AA,AA}^{DZ}=P_{AA}^2+\left(1/2\right)P_{AA}{P}_{Aa}+\left(1/16\right)P_{Aa}^2\hfill \end{array}

(5)

The frequencies {\pi }_{g,g}^{MZ}\left(x\right) and {\pi }_{g1,g2}^{DZ}\left(x\right) of the genotype (g, g) and (g₁, g₂) at any age x for MZ and DZ twin pairs can be calculated from the formulae:

{\pi }_{g,g}^{MZ}\left(x\right)=\frac{P_g{S}_{g,g}^{MZ}\left(x,x\right)}{S^{MZ}\left(x,x\right)}

(6)

{\pi }_{g1,g2}^{DZ}\left(x\right)=\frac{P_{g1,g2}^{DZ}{S}_{g1,g2}^{DZ}\left(x,x\right)}{S^{DZ}\left(x,x\right)}

(7)

Assuming that the variance σ² does not depend on genotype and zygosity, we have the following unknown vector parameter:

\begin{array}{c}\theta =(R,\delta ,\nu ,a_{aa},a_{aA+Aa},b_{aa},b_{aA+Aa},c_{aa},c_{aA+Aa},\hfill \\ c_{AA},{\sigma }^2,{\rho }_{MZ},{\rho }_{DZ}).\hfill \end{array}

Here, ρ _MZ and ρ _DZ are the frailty correlations for MZ and DZ twins. We estimate unknown vector parameter θ maximising the likelihood function:

Li{k}_g=\left(\prod _{i=1}^{N_g^{MZ}}{\pi }_{g_i,g_i}^{MZ}\left(x_i,\theta \right)\right)\left(\prod _{i=1}^{N_g^{DZ}}{\pi }_{g_{i1,i2}}^{DZ}\left(x_i,\theta \right)\right)

(8)

(the maximum likelihood estimates [MLE]), where x _i is the age of twin pair i at the moment of data collection, N_g^{MZ} and N_g^{DZ} are the observed numbers of MZ and DZ twin pairs in the genetic dataset (twin pairs with known genotypes and ages), respectively. To choose the optimal model, we can use the likelihood ratio test for nested models and either the Akaike information criterion (AIC) or the Bayesian information criterion (BIC) for non-nested models. Under the null hypothesis, we assume that a _aa = a _Aa = b _aa = b _Aa = 1 and c _aa = c _Aa = c _AA = 0. Significant deviation from this hypothesis can indicate a gene-longevity association.

If, in addition to genetic data, the data on the longevity of related individuals such as twins are also available, we can use this information to improve the accuracy of statistical estimates and to increase the power. Denote the life spans of the twin pair i in the demographic dataset by (x_i1, x_i2), where I=1,\dots ,N_d^{MZ} for MZ twin pairs and I=1,\dots ,N_d^{DZ} for DZ twin pairs. We assume that twin pairs in the sample are chosen at random and that all twins are deceased. Although the censored data are less informative than non-censored data, they can be also included in the analysis. The bivariate probability density function for a twin pair with longevities x_i1 and x_i2 can be calculated as follows:

\begin{array}{c}\frac{{\partial }^2{S}^j\left(x_{i1},x_{i2}\right)}{\partial x_{i1}\partial x_{i2}}=\sum _{g1,g2}{P}_{g1,g2}^j\frac{S_{g1,g2}^j\left(x_{i1},x_{i1}\right)\partial S_{g1}\left(x_1\right)}{S_{g1}\left(x_1\right)S_{g2}\left(x_2\right)\partial S\left(x_1\right)}\hfill \\ \times \frac{\partial S_{g2}\left(x_2\right)}{\partial S\left(x_2\right)}\frac{\partial S\left(x_1\right)}{\partial x_1}\frac{\partial S\left(x_2\right)}{\partial x_2}\hfill \\ \times \left(1-{\rho }_j^2\right.+\frac{{\rho }_j\left(1-{\rho }_j\right)\left(S_{g1}{\left(x_1\right)}^{-{\sigma }^2}+S_{g2}{\left(x_2\right)}^{-{\sigma }^2}\right)}{\left(S_{g1}{\left(x_1\right)}^{-{\sigma }^2}+S_{g2}{\left(x_2\right)}^{-{\sigma }^2}-1\right)}\hfill \\ \left.+\frac{{\rho }_j\left({\rho }_j+{\sigma }^2\right)S_{g1}{\left(x_1\right)}^{-{\sigma }^2}{S}_{g2}{\left(x_2\right)}^{-{\sigma }^2}}{{\left(S_{g1}{\left(x_1\right)}^{-{\sigma }^2}+S_{g2}{\left(x_2\right)}^{-{\sigma }^2}-1\right)}^2}\right)\hfill \end{array}

(9)

with j = MZ, DZ. We can write the likelihood function for the demographic dataset in the form:

\begin{array}{c}Li{k}_d=\left(\prod _{i=1}^{N_d^{MZ}}\frac{{\partial }^2{S}^{MZ}\left(x_{i1},x_{i2}\right)}{\partial x_{i1}\partial x_{i2}}\left(x_{i1},x_{i2},\theta \right)\right)\\ \times \left(\prod _{i=1}^{N_d^{DZ}}\frac{{\partial }^2{S}^{DZ}\left(x_{i1},x_{i2}\right)}{\partial x_{i1}\partial x_{i2}}\left(x_{i1},x_{i2},\theta \right)\right)\end{array}

(10)

Now, unknown parameters can be found through maximising the joint likelihood function Lik _g × Lik _d .

Results

To carry out the numerical experiments, we used simulated data. To generate datasets with a sample size of N_g^{MZ}=1000, N_g^{DZ}=2000 for genotype data and of N_d^{MZ}=150, N_d^{DZ}=300 for longevity data we assumed that:

• The action of the dominant allele a on longevity can be characterised by parameters a _AA = b _AA = 1, c _aa = c _AA = 0, a _aa = a _Aa = 0.8, b _aa = b _Aa = 1.2;

• The survival function for genotype AA has a form

  \begin{array}{c}\tilde{S}\left(x\right)={\left(1+s^2\tilde{H}\left(x\right)\right)}^{-1/s^2},\\ \tilde{H}\left(x\right)=\tilde{c}x+\mathit{ã}\left(e^{\tilde{b}x}-1\right)/\tilde{b}\end{array}

  (11)

with ã = 2.5·10^-5, \tilde{b}=0.1, \tilde{c}=0 and ln s² = -4.5;

• Individual frailty for twins are gamma-distributed, with mean 1 and variance σ² = 1. Frailty correlations ρ _MZ and ρ _DZ are equal to 0.5 and 0.25, respectively;

• The Hardy-Weinberg equilibrium at the moment of conception holds. There is no genotype selection before birth;

• The slow continuous component of the cohort effect has parameters ν = -2 and δ = 0.005. This corresponds to the frequency P _a ≈ 0.182 for individuals born in year T (the year of data collection) and decreases in the frequency by 0.4 per cent per year. The sudden jump of P _a with parameter R = 0.5 occurred in the cohort T-50;

• The birth dates of all twin pairs from the longevity dataset are uniformly distributed over the cohort interval [T-110, T-100]. The ages of the twins from the genetic dataset at the moment of data collection are uniformly distributed over the age interval [0-105] years.

Nearly one in every 100 deliveries is a twin birth, and the DZ/MZ ratio is approximately equal to 2. From this, it follows that in the stationary population consisting of 300,000 individuals with crude birth and death rates q₀ equal to 15 per 1,000, the life expectancy at birth e₀ is equal to 1,000/15 ≈ 66.7 years and we will find approximately (1/300) × (300,000 × q₀e₀) = 1,000 MZ and (2/300) × (300,000 × q₀e₀) = 2,000 DZ twin pairs. We will also find 150 MZ and 300 DZ newborn twin pairs over the ten-year cohort interval. Since the influence of a decrease in child mortality before the age of 11-13 years on the univariate survival and, therefore, selection is relatively small, we have not included this effect in the simulated data. In general, chosen simulation parameters produce a bivariate lifespan distribution which is similar to the true one.

The estimates of unknown parameters and of the power for 1,000 simulations are given in Table 1. The power was calculated at the 5 per cent significance level. We have used the bivariate and the univariate models applied to the joint bivariate genetic and longevity data or to the bivariate genetic data only. The age dynamics of the hazard functions for genotypes with/without allele a and the age dynamics of the frequencies for the longevity allele/genotypes with the longevity allele are shown in Figures 1 and 2. To establish how often the true bivariate model applied to the bivariate genetic data turns out to be optimal compared with the false univariate model, we used the likelihood ratio test. Significant differences between two these models at a significance level of p < 0.05 were observed in 100 per cent of cases.

Table 1 Parameter estimates (sample means) and their standard deviations (in brackets) for 1,000 simulations, calculated using the bivariate (univariate) model applied to the joint bivariate genetic and longevity data* (***) or to the bivariate genetic data **(****)

Figure 1

Hazard function for genotypes with/without allele a (solid line/dashed line).

Figure 2

Frequency of allele a /genotypes containing allele a (solid line/dashed line).

Discussion

The maximum likelihood method yields correct estimates if the model is correctly specified. In this case, the MLE of unknown parameters under certain regularity conditions are asymptotically unbiased, normal and efficient. If we treat the bivariate data in the same way as the univariate data, and the marginal model is correctly specified, then the robust Hubert-White 'sandwich' estimator of the covariance matrix of parameter estimates yields an asymptotically consistent covariance matrix [8–10]. As we see in Table 1, there is an increase in statistical power when using the more robust univariate model compared with the bivariate model. Nevertheless, the estimates of parameters a _aa and b _aa for the relative risk of the longevity genotype and the estimate of σ for the standard deviation of frailty are closer to their true values if we use the bivariate model. Including the information on longevity in the dataset, however, can substantially improve statistical estimates, increase the power and decrease the variance. It seems that implementation of the approach based on the more robust univariate model, compared with the bivariate model, is preferable for the sample sizes used in this study. Based on the correlation estimates in the MZ and DZ twins, we are able to estimate the contribution of the candidate gene to the heritability [11]. Under the null hypothesis (no heritability), we put ρ _MZ = ρ _DZ . The effect of antagonistic pleiotropy is clearly seen in Figure 1. The presence of allele a in an individual's genotype guarantees the lower hazard of mortality only up to the age of approximately 76 years. The hazard of individuals with genotype AA is then lower than that of individuals with allele a in the genotype. Similar to the univariate model, the bivariate model effectively identifies not only the slow cohort trend of P _g , including the antagonistic pleiotropic effect, but also the sudden change in this parameter. As expected, the frequencies of allele a and of the genotypes containing allele a increase continuously in the age intervals [0-50] and [50-80], fall abruptly at the age of 50 and decrease continuously after the age of 80 (see Figure 2). Univariate and bivariate (for twins) genotype frequencies at the longevity locus at the moment of conception depend on the genotype frequencies in the parental population and the transmission probabilities. In the model we have used, two assumptions were made relating to the longevity locus. First, that the Hardy-Weinberg equilibrium holds for the parental population. Secondly, that the segregation ratio does not deviate from 0.5 [12]. In principle, we can dispense with both of these assumptions and include them as null hypotheses in the study. Significant deviation from the null hypotheses can be tested using the likelihood ratio test. Rejection of the hypothesis about the Hardy-Weinberg equilibrium can indicate possible genotype selection during the gestation period. Significant deviation from Mendelian transmission can mean, for example, that longevity is not governed by the alleles at a single locus. Population admixture and stratification can lead to linkage disequilibrium between longevity and marker loci. In such situations, the study may reveal evidence for ('spurious') association with the marker, even if it is unlinked to the longevity locus. If the sub-population factors influencing the allele frequencies in the marker and longevity loci are identified (eg ethnicity, geographical origin, etc), they can be included in the study. Another solution for this problem is to partition the association effects into between- and within-family components [13, 14]. It was shown that admixture impacts the between-family component estimate, and that the within-family component estimate is independent of any 'spurious' effects when samples from a number of population strata are combined.