Three HCP utilities

Posted on 01.August.2016 by A. M. Winkler

If you are working with data from the Human Connectome Project (HCP), perhaps these three small Octave/MATLAB utilities may be of some use:

hcp2blocks.m: Takes the restricted file with information about kinship and zygosity and produces a multi-level exchangeability blocks file that can be used with PALM for permutation inference. It is fully described here.
hcp2solar.m: Takes restricted and unrestricted files to produce a pedigree file that can be used with SOLAR for heritability and genome-wide association analyses.
picktraits.m: Takes either restricted or unrestricted files, a list of traits and a list of subject IDs to produce tables with selected traits for the selected subjects. These can be used to, e.g., produce design matrices for subsequent analysis.

These functions need to parse relatively large CSV files, which is somewhat inefficient in MATLAB and Octave. Still, since these commands usually have to be executed only once for a particular analysis, a 1-2 minute wait seems acceptable.

If downloaded directly from the above links, remember also to download the prerequisites: strcsvread.m and strcsvwrite.m. Alternatively, clone the full repository from GitHub. The link is this. Other tools may be added in the future.

A fourth utility

For the HCP-S1200 release (March/2017), zygosity information is provided in the fields ZygositySR (self-reported zygosity) and ZygosityGT (zygosity determined by genetic methods for select subjects). If needed, these two fields can be merged into a new field named simply Zygosity. To do so, use a fourth utility, command mergezyg.

Permutation tests in the Human Connectome Project

Posted on 07.December.2015 by A. M. Winkler

Permutation tests are known to be superior to parametric tests: they are based on only few assumptions, essentially that the data are exchangeable, and allow the correction for the multiplicity of tests and the use of various non-standard statistics. The exchangeability assumption allows data to be permuted whenever their joint distribution remains unaltered. Usually this means that each observation needs to be independent from the others.

In many studies, however, there are repeated measurements on the same subjects, which violates exchangeability: clearly, the various measurements obtained from a given subject are not independent from each other. In the Human Connectome Project (HCP) (Van Essen et al, 2012; 2013; see references at the end), subjects are sampled along with their siblings (most of them are twins), such that independence cannot be guaranteed either.

In Winkler et al. (2014), certain structured types of non-independence in brain imaging were addressed through the definition of exchangeability blocks (EBs). Observations within EB can be shuffled freely or, alternatively, the EBs themselves can be shuffled as a whole. This allows various designs that otherwise could not be assessed through permutations.

The same idea can be generalised for blocks that are nested within other blocks, in a multi-level fashion. In the paper Multi-level Block Permutation (Winkler et al., 2015) we presented a method that allows blocks to be shuffled a whole, and inside them, sub-blocks are further allowed to be shuffled, in a recursive process. The method is flexible enough to accommodate permutations, sign-flippings (sometimes also called “wild bootstrap”), and permutations together with sign-flippings.

In particular, this permutation scheme allows the data of the HCP to be analysed via permutations: subjects are allowed to be shuffled with their siblings while keeping the joint distribution intra-sibship maintained. Then each sibship is allowed to be shuffled with others of the same type.

In the paper we show that the error type I is controlled at the nominal level, and the power is just marginally smaller than that would be obtained by permuting freely if free permutation were allowed. The more complex the block structure, the larger the reductions in power, although with large sample sizes, the difference is barely noticeable.

Importantly, simply ignoring family structure in designs as this causes the error rates not to be controlled, with excess false positives, and invalid results. We show in the paper examples of false positives that can arise, even after correction for multiple testing, when testing associations between cortical thickness, cortical area, and measures of body size as height, weight, and body-mass index, all of them highly heritable. Such false positives can be avoided with permutation tests that respect the family structure.

The figure at the top shows how the subjects of the HCP (terminal dots, shown in white colour) can be shuffled or not, while respecting the family structure. Blue dots indicate branches that can be permuted, whereas red dots indicate branches that cannot (see the main paper for details). This diagram includes 232 subjects of an early public release of HCP data. The tree on the left considers dizygotic twins as a category on their own, i.e., that cannot be shuffled with ordinary siblings, whereas the tree on the right considers dizygotic twins as ordinary siblings.

The first applied study using our strategy has just appeared. The method is implemented in the freely available package PALM — Permutation Analysis of Linear Models, and a set of practical steps to use it with actual HCP data is available here.

References

Van Essen DC, Ugurbil K, Auerbach E, Barch D, Behrens TE, Bucholz R, Chang A, Chen L, Corbetta M, Curtiss SW, Della Penna S, Feinberg D, Glasser MF, Harel N, Heath AC, Larson-Prior L, Marcus D, Michalareas G, Moeller S, Oostenveld R, Petersen SE, Prior F, Schlaggar BL, Smith SM, Snyder AZ, Xu J, Yacoub E; WU-Minn HCP Consortium. The Human Connectome Project: a data acquisition perspective. Neuroimage. 2012;62(4):2222-31.
Van Essen DC, Smith SM, Barch DM, Behrens TEJ, Yacoub E, Ugurbil K. The WU-Minn Human Connectome Project: An overview. Neuroimage. 2013;80:62-79.
Winkler AM, Ridgway GR, Webster MA, Smith SM, Nichols TE. Permutation inference for the general linear model. Neuroimage. 2014;92:381-97. (Open Access)
Winkler AM, Webster MA, Vidaurre D, Nichols TE, Smith SM. Multi-level block permutation. Neuroimage. 2015;123:253-68. (Open Access)

Variance components in genetic analyses

Posted on 02.August.2015 by A. M. Winkler

Pedigree-based analyses allow investigation of genetic and environmental influences on anatomy, physiology, and behaviour.

Methods based on components of variance have been used extensively to assess genetic influences and identify loci associated with various traits quantifying aspects of anatomy, physiology, and behaviour, in both normal and pathological conditions. In an earlier post, indices of genetic resemblance between relatives were presented, and in the last post, the kinship matrix was defined. In this post, these topics are used to present a basic model that allows partitioning of the phenotypic variance into sources of variation that can be ascribed to genetic, environmental, and other factors.

A simple model

Consider the model:

$\mathbf{Y} = \mathbf{X}\mathbf{B} + \boldsymbol{\Upsilon}$

where, for $S$ subjects, $T$ traits, $P$ covariates and $K$ variance components, $\mathbf{Y}_{S \times T}$ are the observed trait values for each subject, $\mathbf{X}_{S \times P}$ is a matrix of covariates, $\mathbf{B}_{P \times T}$ is a matrix of unknown covariates’ weights, and $\boldsymbol{\Upsilon}_{S \times T}$ are the residuals after the covariates have been taken into account.

The elements of each column $t$ of $\boldsymbol{\Upsilon}$ are assumed to follow a multivariate normal distribution $\mathcal{N}\left(0;\mathbf{S}\right)$ , where $\mathbf{S}$ is the between-subject covariance matrix. The elements of each row $s$ of $\boldsymbol{\Upsilon}$ are assumed to follow a normal distribution $\mathcal{N}\left(0;\mathbf{R}\right)$ , where $\mathbf{R}$ is the between-trait covariance matrix. Both $\mathbf{R}$ and $\mathbf{S}$ are seen as the sum of $K$ variance components, i.e. $\mathbf{R} = \sum_{k} \mathbf{R}_{k}$ and $\mathbf{S} = \sum_{k} \mathbf{S}_{k}$ . For a discussion on these equalities, see Eisenhart (1947) [see references at the end].

An equivalent model

The same model can be written in an alternative way. Let $\mathbf{y}_{S \cdot T \times 1}$ be the stacked vector of traits, $\mathbf{\tilde{X}}_{S \cdot T \times P \cdot T} = \mathbf{X} \otimes \mathbf{I}_{T \times T}$ is the matrix of covariates, $\boldsymbol{\beta}_{P \cdot T \times 1}$ the vector with the covariates’ weights, $\boldsymbol{\upsilon}_{S \cdot T \times 1}$ the residuals after the covariates have been taken into account, and $\otimes$ represent the Kronecker product. The model can then be written as:

$\mathbf{y} = \mathbf{\tilde{X}}\boldsymbol{\beta} + \boldsymbol{\upsilon}$

The stacked residuals $\boldsymbol{\upsilon}$ is assumed to follow a multivariate normal distribution $\mathcal{N}\left(0;\boldsymbol{\Omega}\right)$ , where $\boldsymbol{\Omega}$ can be seen as the sum of $K$ variance components:

$\boldsymbol{\Omega} = \sum_{k} \mathbf{R}_k \otimes \mathbf{S}_k$

The $\boldsymbol{\Omega}$ here is the same as in Almasy and Blangero (1998). $\mathbf{S}_k$ can be modelled as correlation matrices. The associated scalars are absorbed into the (to be estimated) $\mathbf{R}_k$ . $\mathbf{R}$ is the phenotypic covariance matrix between the residuals of the traits:

$\mathbf{R} = \left[ \begin{array}{ccc} \mathsf{Var}(\upsilon_1) & \cdots & \mathsf{Cov}(\upsilon_1,\upsilon_T) \\ \vdots & \ddots & \vdots \\ \mathsf{Cov}(\upsilon_T,\upsilon_1) & \cdots & \mathsf{Var}(\upsilon_T) \end{array}\right]$

whereas each $\mathbf{R}_k$ are the share of these covariances attributable to the $k$ -th component:

$\mathbf{R}_k = \left[ \begin{array}{ccccc} \mathsf{Var}_k(\upsilon_1) & \cdots & \mathsf{Cov}_k(\upsilon_1,\upsilon_T) \\ \vdots & \ddots & \vdots \\ \mathsf{Cov}_k(\upsilon_T,\upsilon_1) & \cdots & \mathsf{Var}_k(\upsilon_T) \end{array}\right]$

$\mathbf{R}$ can be converted to a between-trait phenotypic correlation matrix $\mathbf{\mathring{R}}$ as:

$\mathbf{\mathring{R}} = \left[ \begin{array}{ccc} \frac{\mathsf{Var}(\upsilon_1)}{\mathsf{Var}(\upsilon_1)} & \cdots & \frac{\mathsf{Cov}(\upsilon_1,\upsilon_T)}{\left(\mathsf{Var}(\upsilon_1)\mathsf{Var}(\upsilon_T)\right)^{1/2}} \\ \vdots & \ddots & \vdots \\ \frac{\mathsf{Cov}(\upsilon_1,\upsilon_T)}{\left(\mathsf{Var}(\upsilon_1)\mathsf{Var}(\upsilon_T)\right)^{1/2}} & \cdots & \frac{\mathsf{Var}(\upsilon_T)}{\mathsf{Var}(\upsilon_T)} \end{array}\right]$

The above phenotypic correlation matrix has unit diagonal and can still be fractioned into their $K$ components:

$\mathbf{\mathring{R}}_k = \left[ \begin{array}{ccc} \frac{\mathsf{Var}_k(\upsilon_1)}{\mathsf{Var}(\upsilon_1)} & \cdots & \frac{\mathsf{Cov}_k(\upsilon_1,\upsilon_T)}{\left(\mathsf{Var}(\upsilon_1)\mathsf{Var}(\upsilon_T)\right)^{1/2}} \\ \vdots & \ddots & \vdots \\ \frac{\mathsf{Cov}_k(\upsilon_T,\upsilon_1)}{\left(\mathsf{Var}(\upsilon_T)\mathsf{Var}(\upsilon_1)\right)^{1/2}} & \cdots & \frac{\mathsf{Var}_k(\upsilon_T)}{\mathsf{Var}(\upsilon_T)} \end{array}\right]$

The relationship $\mathbf{\mathring{R}} = \sum_k \mathbf{\mathring{R}}_k$ holds. The diagonal elements of $\mathbf{\mathring{R}}_k$ may receive particular names, e.g., heritability, environmentability, dominance effects, shared enviromental effects, etc., depending on what is modelled in the corresponding $\mathbf{S}_k$ . However, the off-diagonal elements of $\mathbf{\mathring{R}}_k$ are not the $\rho_k$ that correspond, e.g. to the genetic or environmental correlation. These off-diagonal elements are instead the signed $\text{\textsc{erv}}$ when $\mathbf{S}_k=2\cdot\boldsymbol{\Phi}$ , or their $\text{\textsc{erv}}_k$ -equivalent for other variance components (see below). In this particular case, they can also be called “bivariate heritabilities” (Falconer and MacKay, 1996). A matrix $\mathbf{\breve{R}}_k$ that contains these correlations $\rho_k$ , which are the fraction of the variance attributable to the $k$ -th component that is shared between pairs of traits is given by:

$\mathbf{\breve{R}}_k = \left[ \begin{array}{ccc} \frac{\mathsf{Var}_k(\upsilon_1)}{\mathsf{Var}_k(\upsilon_1)} & \cdots & \frac{\mathsf{Cov}_k(\upsilon_1,\upsilon_T)}{\left(\mathsf{Var}_k(\upsilon_1)\mathsf{Var}_k(\upsilon_T)\right)^{1/2}} \\ \vdots & \ddots & \vdots \\ \frac{\mathsf{Cov}_k(\upsilon_T,\upsilon_1)}{\left(\mathsf{Var}_k(\upsilon_T)\mathsf{Var}_k(\upsilon_1)\right)^{1/2}} & \cdots & \frac{\mathsf{Var}_k(\upsilon_T)}{\mathsf{Var}_k(\upsilon_T)} \end{array}\right]$

As for the phenotypic correlation matrix, each $\mathbf{\breve{R}}_k$ has unit diagonal.

The most common case

A particular case is when $\mathbf{S}_1 = 2\cdot\boldsymbol{\Phi}$ , the coefficient of familial relationship between subjects, and $\mathbf{S}_2 = \mathbf{I}_{S \times S}$ . In this case, the $T$ diagonal elements of $\mathbf{\mathring{R}}_1$ represent the heritability ( $h_t^2$ ) for each trait $t$ . The diagonal of $\mathbf{\mathring{R}}_2$ contains $1-h_t^2$ , the environmentability. The off-diagonal elements of $\mathbf{\mathring{R}}_1$ contain the signed $\text{\textsc{erv}}$ (see below). The genetic correlations, $\rho_g$ are the off-diagonal elements of $\mathbf{\breve{R}}_1$ , whereas the off-diagonal elements of $\mathbf{\breve{R}}_2$ are $\rho_e$ , the environmental correlations between traits. In this particular case, the components of $\mathbf{R}$ , i.e., $\mathbf{R}_k$ are equivalent to $\mathbf{G}$ and $\mathbf{E}$ covariance matrices as in Almasy et al (1997).

Relationship with the ERV

To see how the off-diagonal elements of $\mathbf{\mathring{R}}_k$ are the signed Endophenotypic Ranking Values for each of the $k$ -th variance component, $\text{\textsc{erv}}_k$ (Glahn et al, 2011), note that for a pair of traits $i$ and $j$ :

$\mathring{R}_{kij} = \frac{\mathsf{Cov}_k(\upsilon_i,\upsilon_j)}{\left(\mathsf{Var}(\upsilon_i)\mathsf{Var}(\upsilon_j)\right)^{1/2}}$

Multiplying both numerator and denominator by $\left(\mathsf{Var}_k(\upsilon_i)\mathsf{Var}_k(\upsilon_j)\right)^{1/2}$ and rearranging the terms gives:

$\mathring{R}_{kij} = \frac{\mathsf{Cov}_k(\upsilon_i,\upsilon_j)}{\left(\mathsf{Var}_k(\upsilon_i)\mathsf{Var}_k(\upsilon_j)\right)^{1/2}} \left(\frac{\mathsf{Var}_k(\upsilon_i)}{\mathsf{Var}(\upsilon_i)}\frac{\mathsf{Var}_k(\upsilon_j)}{\mathsf{Var}(\upsilon_j)}\right)^{1/2}$

When $\mathbf{S}_k = 2\cdot\boldsymbol{\Phi}$ , the above reduces to $\mathring{R}_{kij} = \rho_k \sqrt{h^2_i h^2_j}$ , which is the signed version of $\text{\textsc{erv}}=\left|\rho_g\sqrt{h_i^2h_j^2}\right|$ when $k$ is the genetic component.

Positive-definiteness

$\mathbf{R}$ and $\mathbf{R}_k$ are covariance matrices and so, are positive-definite, whereas the correlation matrices $\mathbf{\mathring{R}}$ , $\mathbf{\mathring{R}}_k$ and $\mathbf{\breve{R}}_k$ are positive-semidefinite. A hybrid matrix that does not have to be positive-definite or semidefinite is:

$\mathbf{\check{R}}_k = \mathbf{I} \odot \mathbf{\mathring{R}}_k + \left(\mathbf{J}-\mathbf{I}\right) \odot \mathbf{\breve{R}}_k$

where $\mathbf{J}$ is a matrix of ones, $\mathbf{I}$ is the identity, both of size $T \times T$ , and $\odot$ is the Hadamard product. An example of such matrix of practical use is to show concisely the heritabilities for each trait in the diagonal and the genetic correlations in the off-diagonal.

Cauchy-Schwarz

Algorithmic advantages can be obtained from the positive-definiteness of $\mathbf{\mathring{R}}_k$ . The Cauchy–Schwarz theorem (Cauchy, 1821; Schwarz, 1888) states that:

$\mathring{R}_{kij} \leqslant \sqrt{\mathring{R}_{kii}\mathring{R}_{kjj}}$

Hence, the bounds for the off-diagonal elements can be known from the diagonal elements, which, by their turn, are estimated in a simpler, univariate model.

The Cauchy-Schwarz inequality imposes limits on the off-diagonal values of the matrix that contains the genetic covariances (or bivariate heritabilities).

Parameter estimation

Under the multivariate normal assumption, the parameters can be estimated maximising the following loglikelihood function:

$\mathcal{L}\left(\mathbf{R}_k,\boldsymbol{\beta}\Big|\mathbf{y},\mathbf{\tilde{X}}\right) = -\frac{1}{2} \left(N \ln 2\pi + \ln \left|\boldsymbol{\Omega}\right| + \left(\mathbf{y}-\mathbf{\tilde{X}}\boldsymbol{\beta}\right)'\boldsymbol{\Omega}\left(\mathbf{y}-\mathbf{\tilde{X}}\boldsymbol{\beta}\right)\right)$

where $N=S \cdot T$ is the number of observations on the stacked vector $\mathbf{y}$ . Unbiased estimates for $\boldsymbol{\beta}$ , although inefficient and inappropriate for hypothesis testing, can be obtained with ordinary least squares (OLS).

Parametric inference

Hypothesis testing can be performed with the likelihood ratio test (LRT), i.e., the test statistic is produced by subtracting from the loglikelihood of the model in which all the parameters are free to vary ( $\mathcal{L}_1$ ), the loglikelihood of a model in which the parameters being tested are constrained to zero, the null model ( $\mathcal{L}_0$ ). The statistic is given by $\lambda = 2\left(\mathcal{L}_1-\mathcal{L}_0\right)$ (Wilks, 1938), which here is asymptotically distributed as a 50:50 mixture of a $\chi^2_0$ and $\chi^2_{\text{df}}$ distributions, where df is the number of parameters being tested and free to vary in the unconstrained model (Self and Liang, 1987). From this distribution the p-values can be obtained.

References

Almasy L, Blangero J. Multipoint quantitative-trait linkage analysis in general pedigrees. Am J Hum Genet. 1998 May;62(5):1198–211.
Almasy L, Dyer TD, Blangero J. Bivariate quantitative trait linkage analysis: pleiotropy versus co-incident linkages. Genet Epidemiol. 1997 Jan;14(6):953–8.
Cauchy, A-L. Œuvres Completes D’Augustin Cauchy (published in 1882). Vol. XV. Gauthier-Villars et Fils, Paris, 1821.
Eisenhart C. The assumptions underlying the analysis of variance. Biometrics. 1947;3(1):1–21.
Falconer DS, Mackay TFC. Introduction to Quantitative Genetics. Addison Wesley Longman, Harlow, Essex, UK, 1996.
Glahn DC, Curran JE, Winkler AM, Carless MA, Kent JW, Charlesworth JC, et al. High dimensional endophenotype ranking in the search for major depression risk genes. Biol Psychiatry. 2012 Jan 1;71(1):6–14.
Schwarz HA. Über ein Flächen kleinsten Flächeninhalts betreffendes Problem der Variationsrechnung. Acta Societatis Scientiarum Fennicæ XV, 319–332, 1888.
Self SG, Liang KY. Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. J Am Stat Assoc. 1987;82(398):605–10.
Wilks SS. The Large-Sample Distribution of the Likelihood Ratio for Testing Composite Hypotheses. Ann Math Stat. 1938;9(1):60–2.

The photograph at the top (elephants) is by Anja Osenberg and was generously released into public domain.

Understanding the kinship matrix

Posted on 29.July.2015 by A. M. Winkler

Coefficients to assess the genetic resemblance between individuals were presented in the last post. Among these, the coefficient of kinship, $\phi$ , is probably the most interesting. It gives a probabilistic estimate that a random gene from a given subject $i$ is identical by descent (ibd) to a gene in the same locus from a subject $j$ . For $N$ subjects, these probabilities can be assembled in a $N \times N$ matrix termed kinship matrix, usually represented as $\mathbf{\Phi}$ , that has elements $\phi_{ij}$ , and that can be used to model the covariance between individuals in quantitative genetics.

Consider the pedigree in the figure below, consisted of 14 subjects:

The corresponding kinship matrix, already multiplied by two to indicate expected covariances between subjects (i.e., $2\cdot\mathbf{\Phi}$ ), is:

Note that the diagonal elements can have values above unity, given the consanguineous mating in the family (between s09 and s12, indicated by the double line in the pedigree diagram).

In the next post, details on how the kinship matrix can be used investigate heritabilities, genetic correlations, and to perform association studies will be presented.