Deriving A TWAS Test Statistic

From SNP Associations to a Gene-Level Score

A transcriptome-wide association study (TWAS) tests whether a gene’s genetically predicted expression is associated with a trait. The prediction is built from SNP genotypes and expression-prediction weights. If SNPs are correlated through linkage disequilibrium, their association statistics are correlated too. Therefore, a weighted sum of SNP-level statistics needs a variance calculation that accounts for their covariance.

This tutorial derives that score from a linear-model perspective and then expresses it using summary statistics. The target is the test statistic for one gene and one trait. The statistic measures evidence of association; by itself, it does not establish that changing the gene’s expression would cause a change in the trait.

Model and Assumptions

For the derivation, suppose each SNP genotype has been centered and scaled to variance one in the association sample. Write \(X\) for the \(n\times m\) genotype matrix, with \(n\) individuals and \(m\) SNPs, so that \(X^\mathsf{T}X/n\) estimates \(R\). Let \(y\) be a quantitative trait after centering and, if needed, adjusting for covariates. The expression predictor for each individual is the weighted genotype sum \(Xw\). The weights are treated as fixed for this test; for example, they may have been estimated in an independent expression reference sample.

For a simple derivation, take the null model to be \(y=\varepsilon\), where the residuals have mean zero, variance \(\sigma^2\), and are independent across individuals conditional on the genotypes. The linear-model score calculation also applies after suitable covariate adjustment, using genotypes and outcomes residualized against those covariates. The normal approximation requires enough information for the SNP association statistics to be approximately jointly normal. The association statistics and LD matrix must refer to harmonized SNPs, alleles, and a sufficiently compatible population. In practice, differences in sample composition or allele coding can violate that alignment.

Define the SNP-level score component \(z_j\) as \(X_j^\mathsf{T}y/(\sigma\sqrt{n})\), where \(X_j\) is column \(j\) of \(X\). In applied analyses, this component is represented by a SNP association \(z\)-statistic, typically an estimated effect divided by its standard error. Under the null, the components have variance approximately one, and their covariance is determined by the genotype correlations.

Why the Denominator Contains the SNP Covariance

The unstandardized gene-level score is \(w^\mathsf{T}z\): each SNP association is multiplied by its expression-prediction weight and the contributions are summed. To turn this into a test statistic, calculate the variance of that sum under the null.

The quadratic form expands to the sum of each SNP’s variance contribution and every pairwise covariance contribution:

Since the SNP statistics have unit variance under the stated null approximation, their covariance is approximately \(R_{jk}\). Dividing the weighted sum by its null standard deviation therefore gives \(Z_{\mathrm{TWAS}}\). This is a standardization, not a second weighting step: the numerator aggregates associations, and the denominator scales that aggregate according to its variability.

The same result follows directly from the individual-level regression score. The score for testing association of \(y\) with \(Xw\) at a zero gene-level effect is \((Xw)^\mathsf{T}y\). Under the null, its variance is \(\sigma^2(Xw)^\mathsf{T}(Xw)\). Since \((Xw)^\mathsf{T}y=w^\mathsf{T}X^\mathsf{T}y\) and \((Xw)^\mathsf{T}(Xw)\approx n w^\mathsf{T}Rw\), dividing the score by its standard deviation yields the same expression as the summary-statistic formula.

Using Effect Estimates and Standard Errors

The formula can also be written with SNP effect estimates, but the units and covariance need to be handled consistently. Let \(\hat b_j\) be the SNP-trait effect estimate and \(s_j\) its standard error. If \(z_j=\hat b_j/s_j\), collect the standard errors in the diagonal matrix \(D=\operatorname{diag}(s_1,\ldots,s_m)\), so that \(\hat b=Dz\). When the null covariance of \(z\) is \(R\), the covariance of the effect estimates is approximately \(DRD\). Consequently, the weighted effect sum has variance approximately \(w^\mathsf{T}DRDw\), and its standardized statistic is $$ \frac{w^\mathsf{T}\hat b}{\sqrt{w^\mathsf{T}DRDw}}. $$ When all standard errors are the same value \(s\), this reduces to \(w^\mathsf{T}\hat b/(s\sqrt{w^\mathsf{T}Rw})\). A formula that mixes effect estimates with a denominator intended for \(z\)-statistics will generally have the wrong scale.

Covariance Can Increase or Decrease the Score’s Scale

The off-diagonal terms matter in relation to both the SNP correlations and the signs of the weights. Positively correlated SNPs with weights of the same sign add positive covariance terms, increasing the variance of the sum. If the weights have opposite signs, the same positive correlation can reduce that variance. The numerator alone does not reveal either effect.

Assumptions, Checks, and Alternatives

Before interpreting a calculated score, check that the weight and association vectors refer to the same variants and allele orientation. If an allele is reversed for a SNP, its effect estimate or \(z\)-statistic must be reoriented consistently; otherwise, a contribution can acquire the wrong sign. Also check that \(R\) is a correlation matrix for the relevant SNPs in a population compatible with the trait-association sample. A mismatched or poorly estimated covariance matrix can distort the denominator, even if the numerator is calculated correctly.

The quadratic form must be nonnegative because it is a variance. A negative computed value indicates a problem such as an invalid covariance matrix, inconsistent SNP ordering, or numerical error. If \(w^\mathsf{T}Rw=0\), the predicted genotype score has zero variance in the modeled data, so the standardized statistic is undefined. Very small values can also make the calculation sensitive to small errors in weights or covariance estimates. These are reasons to inspect the inputs, not to force a p-value from an invalid denominator.

This covariance-adjusted score differs from simply summing SNP \(z\)-statistics or applying an independence-based denominator. Those shortcuts are valid only when their covariance assumptions hold. An individual-level regression on predicted expression is a related alternative: when it uses the same individuals, genotypes, weights, and covariate adjustment, its score calculation gives the same standardized result. It may be unavailable when only summary statistics can be analyzed. With summary data, the matrix expression makes the required covariance explicit, but its quality depends on having appropriate LD information.

The null normal approximation is not a guarantee that every practical analysis is calibrated. It relies on suitable association statistics and covariance estimates, and can be affected by ancestry mismatch, sample structure, or dependence between the expression-weight training data and the trait sample. The score also tests the weighted predicted-expression component represented by the chosen \(w\); it is not a direct test of every possible expression effect at the gene. In particular, association from this test alone does not demonstrate a causal expression mechanism.