Permutation Calibration For Cis-QTL Mapping

Why Calibrate a Cis-QTL Scan by Permutation?

A cis-QTL scan tests many nearby variants for association with one molecular phenotype. Nearby genotypes are often correlated through linkage disequilibrium (LD), so the variant-level test statistics are correlated too. A significance calculation that treats those tests as independent does not represent the actual scan. Permutation calibration offers a way to measure how large a scan-wide statistic can be under a specified null model while retaining the observed dependence among genotypes.

The calibration is only as reliable as its permutation scheme. Shuffling phenotype values is justified when they are exchangeable under the null being tested. Covariates, family relationships, repeated measurements, or other sample structure can break that condition. The practical task is therefore twofold: define the scan statistic and construct null datasets that respect the model and design.

Definition: For one molecular phenotype and a specified cis window, let \(T_j\) be the test statistic for association with variant \(j\). The cis-window maximum statistic is \(T_{\max}=\max_{j=1,\ldots,m}|T_j|\) for two-sided tests. It records the strongest association statistic anywhere in the window. The empirical permutation p-value compares the observed \(T_{\max}\) with its distribution under null datasets generated by the chosen permutation scheme.

A cis window is the prespecified genomic region around the molecular trait's gene, as defined for the analysis. The genotype data are held fixed. In each permutation, the phenotype is reassigned or reconstructed under the null, and the same association model is refit for every variant in the window. Using the same permuted phenotype for all variants is essential: it retains the observed LD pattern among the genotypes and therefore the correlation among their tests.

Empirical Significance From Permuted Maxima

Suppose \(B\) null datasets are generated. Let \(T_{\max}^{(b)}\) be the maximum absolute test statistic in permutation \(b\), for \(b=1,\ldots,B\). The usual finite-simulation estimate is

\hat{p}_{\mathrm{emp}}= \frac{1+\sum_{b=1}^{B} I\!\left(T_{\max}^{(b)}\geq T_{\max}^{\mathrm{obs}}\right)} {B+1}.

Here \(I(\cdot)\) equals 1 when the stated condition holds and 0 otherwise. The added 1 in numerator and denominator prevents a reported estimate of zero when no simulated maximum exceeds the observed one. With \(B\) permutations, the smallest reportable value is \(1/(B+1)\). More permutations are needed to estimate very small tail probabilities with useful precision.

This calculation asks: under the specified null and permutation scheme, how often would a scan across this cis window produce a maximum statistic at least as large as the observed maximum? It is a scan-level calibration for that phenotype and window. It is not the probability that the null is true, nor does it identify a causal variant.

Worked Example: A Cis Window With Correlated Variants

Consider a hypothetical expression phenotype tested at three variants in a prespecified cis window. The dosage correlations are 0.82 between variants 1 and 2, 0.76 between variants 2 and 3, and 0.61 between variants 1 and 3. The observed two-sided \(t\)-statistics from the fitted association models are 2.10, \(-3.40\), and 2.75. The scan statistic is therefore

T_{\max}^{\mathrm{obs}}=\max(|2.10|,|-3.40|,|2.75|)=3.40.

For illustration, 99 phenotype permutations are performed. In every permutation, one null phenotype vector is used for all three variants; each model is refit and the maximum absolute \(t\)-statistic is saved. Suppose 4 of the 99 permuted maxima are at least 3.40. Then

\hat{p}_{\mathrm{emp}}=\frac{1+4}{99+1}=\frac{5}{100}=0.0500.

Under the null represented by this permutation scheme, an entire three-variant scan reaches a maximum statistic of 3.40 or larger about 5% of the time by this empirical estimate. Permuting separately for each variant would not preserve the joint pattern of correlated tests and would not estimate this same scan-level tail probability.

Choosing a Valid Permutation Scheme

In the simplest setting, the null model says the phenotype is unrelated to all genotypes in the cis window, and the phenotype observations are exchangeable across donors. Exchangeability means that under the null, rearranging the donor labels of the phenotype does not change the joint distribution of the data. If this condition holds, a permutation of the phenotype vector provides a null dataset while the genotypes remain fixed.

Exchangeability is not the same as independence of the variant tests. The variants may be highly correlated; that is precisely why their genotypes are kept fixed and the same permuted phenotype is tested against all of them. The relevant question is whether the phenotype assignments are exchangeable across the units being permuted.

If the association model includes covariates, such as age, sex, or technical factors, raw phenotype values generally cannot be shuffled without consequence. A shuffle can destroy the phenotype-covariate relationship even when the covariates are not associated with genotype. One common approach is the Freedman–Lane residual permutation:

Fit the reduced model.
Fit the null model containing the covariates but no cis genotype term. Save its fitted values and residuals.

Permute residuals.
Rearrange the reduced-model residuals across exchangeable donors. If the design requires restricted permutations, rearrange only within the valid exchangeability groups.

Construct a null phenotype.
Add the permuted residuals to the original reduced-model fitted values, preserving the fitted covariate pattern.

Refit the cis scan.
For every variant, fit the same covariate-adjusted association model and save the maximum statistic across the window.

This method is not a universal guarantee of exact calibration. Its validity relies on residuals being exchangeable under the null, and finite-sample behavior can depend on the design and error distribution. Heteroskedasticity, strong leverage differences, related donors, or repeated observations may call for a design-specific method, restricted permutations, or a model-based null simulation. The analysis should state which null is being tested and why the generated null datasets represent it.

Worked Example: Preserving a Covariate Relationship

Suppose six hypothetical donors have a quantitative covariate \(x=(-2,-1,0,0,1,2)\) and phenotype \(y=(4,5,5,6,7,9)\). The reduced model contains an intercept and \(x\), but no genotype. Since \(\bar{x}=0\) and \(\bar{y}=6\), the fitted slope is

\hat{\beta}_x= \frac{\sum_i (x_i-\bar{x})(y_i-\bar{y})}{\sum_i (x_i-\bar{x})^2} =\frac{12}{10}=1.20.

The fitted values are \((3.6,4.8,6.0,6.0,7.2,8.4)\), and the residuals are \((0.4,0.2,-1.0,0,-0.2,0.6)\). For one illustrative permutation, rearrange these residuals as \((0,-1.0,0.6,0.4,0.2,-0.2)\). Adding them to the original fitted values gives the null phenotype

$$ y^*=(3.6,3.8,6.6,6.4,7.4,8.2). $$

This construction retains the fitted trend with \(x\) while rearranging residuals. The resulting phenotype is then tested against every cis variant using the same covariate-adjusted model. By contrast, directly shuffling the six raw \(y\) values would not preserve their systematic relationship with \(x\). This example illustrates the construction; with only six donors, it would not by itself justify a reliable tail estimate.

Justification: Why the Maximum Calibrates the Scan

Under the complete null, assume the phenotype assignments are exchangeable under the chosen scheme and the genotype data are treated as fixed. Conditional on those genotypes, each valid rearrangement represents a possible null assignment of phenotype to donors. For each assignment, calculating the maximum across the same \(m\) variants reproduces both the scan procedure and the dependence among its tests.

The permutation estimate is the fraction of these null scans whose maximum is at least as extreme as the observed maximum, with a finite-simulation correction. Its key calibration feature follows directly: it compares like with like, namely an observed scan maximum with null scan maxima generated by repeating the whole scan. The argument does not require independence among variant statistics. It does require that the null assignments be valid and that the same analysis and window be used for observed and permuted data.

A variant-specific adjusted comparison can also compare its observed statistic with the distribution of permuted window maxima. Such a comparison addresses whether that variant's result is unusually large relative to the scan-wide null threshold. Strong family-wise error control beyond the complete null requires additional conditions on the joint null distribution; it should not be assumed merely because permutations were used. Broader choices about error-rate control across variants and molecular traits require separate consideration.

Worked Example: Detecting a Sample-Structure Problem

In a hypothetical study, 40 donors come from two recruitment centers with different average expression levels and different allele frequencies. A naive unrestricted phenotype permutation gives 12 exceedances among 999 permutations for the observed window maximum. The finite-permutation estimate is

\hat{p}_{\mathrm{emp,naive}}=\frac{1+12}{999+1}=0.0130.

The number alone does not establish valid calibration: shuffling expression values across centers breaks the center-specific phenotype distribution. Suppose instead the analysis includes center in the reduced model and uses a scheme that permutes residuals only among donors for whom residual exchangeability is justified within center. In a hypothetical set of 999 such permutations, 39 maxima meet or exceed the observed value. The estimate becomes

\hat{p}_{\mathrm{emp,restricted}}=\frac{1+39}{999+1}=0.0400.

The difference indicates that the calibration depends on the null construction. It is not evidence that one p-value is automatically correct: the restricted approach is appropriate only if its within-center exchangeability assumption is defensible and the model addresses the structure relevant to the study. With family relatedness or repeated measures, simply permuting within center may still be inadequate; the permutation units or null simulation must reflect those dependencies.

Common Mistakes and Practical Checks

Permuting independently for each variant. This loses the joint structure of scan statistics. Use a shared null phenotype assignment across all variants in each permutation.
Shuffling raw phenotypes despite covariates. This can break covariate relationships and create an incorrect null. Use a justified covariate-aware scheme, such as residual permutation, and check its assumptions.
Assuming that the word “permutation” guarantees validity. The rearranged units must be exchangeable under the null. Consider ancestry, recruitment center, family, batch, and repeated observations when defining those units.
Reporting zero when no permutation exceeds the observed statistic. With the finite-permutation estimate, zero exceedances gives \(1/(B+1)\), not zero. Report the number of permutations because it determines tail resolution.
Interpreting a scan-level result as a causal-variant finding. A significant cis-window maximum indicates evidence against the specified no-association null for the scan. LD can make several variants show similar statistics, and the maximum does not resolve which variant is causal.

Key takeaway: A defensible cis-QTL permutation calibration repeats the entire scan on null phenotypes, uses the same phenotype assignment across correlated variants, and respects covariates and sample structure. The empirical p-value is meaningful only relative to the null and exchangeability assumptions used to generate those phenotypes.

Check Your Understanding

Answer each question by identifying both the calculation and the assumption it depends on.

A cis scan tests 80 variants. Why should the same permuted phenotype be used for all 80 variants in a given permutation?
In 4,999 permutations, 24 maximum statistics equal or exceed the observed maximum. Calculate the finite-permutation empirical p-value.
Why can shuffling raw phenotype values be inappropriate when the fitted model includes a covariate associated with phenotype?
What feature of the study design might require restricted permutations or model-based null simulation rather than unrestricted donor shuffling?
Does a significant empirical p-value for the cis-window maximum identify the causal variant? Explain what it does establish.