Worked QTL Analysis From Phenotype Processing To Association Report

From Molecular Measurements to a Defensible QTL Report

A molecular QTL result is not just the smallest p-value in a table. It depends on how the molecular phenotype was processed, which covariates were included, what variants counted as cis, and how many hypotheses were tested. A concise report should make those choices visible and connect them to an interpretable effect estimate.

The worked analysis below is entirely hypothetical. It follows one gene measured in 120 unrelated donors and tests variants in a prespecified cis window around that gene. The aim is to estimate and test an association, not to establish that a variant causes a change in expression.

Process and Inspect the Phenotype

Start by defining the phenotype and its scale. For an expression analysis, a typical workflow filters genes with insufficient counts, normalizes for library size and other assay features, and transforms the normalized measurements if a Gaussian linear model is planned. The transformation, normalization method, and any filtering thresholds should be fixed and recorded rather than chosen to make a particular association appear stronger.

For illustration, suppose a hypothetical sample has 240 reads assigned to a gene and a library size of 20 million reads. Its count per million (CPM) is \(240/(20{,}000{,}000)\times 1{,}000{,}000=12\). With a prespecified pseudocount of 0.5, the transformed value is \(\log_2(12+0.5)\), approximately 3.644. The pseudocount keeps the logarithm defined at zero; it also means the transformed scale is not an exact log of unadjusted abundance, especially for very small values.

Before fitting a model, check that the phenotype has enough usable observations and that the processing did not introduce obvious artifacts. A transformation does not guarantee normally distributed residuals, and no transformation can repair sample swaps, a confounded batch design, or a poorly defined phenotype. Earlier material on diagnosing molecular QTL results and technical artifacts is relevant to these checks.

Specify the Association Model Before Testing

For the hypothetical analysis, let \(Y_i\) be the transformed expression measurement for donor \(i\), and let \(G_i\in\{0,1,2\}\) be the number of copies of the tested effect allele. The model uses an additive genotype dosage and adjusts for prespecified age, sex, assay batch, and three ancestry principal components. Age and the principal components are centered; categorical variables are represented by indicator variables. These choices must be supported by the study design and data, not selected only because they produce a small p-value.

Here \(\beta_G\) is the conditional mean difference in transformed expression per additional effect allele, holding the included covariates fixed. The usual linear-model test and confidence interval rely on a correctly specified conditional mean, independent observations, no perfect collinearity, and errors with constant variance; exact small-sample t-based inference additionally assumes normally distributed errors conditional on the predictors. Genotype dosage should be coded consistently, and the tested variants and phenotype must be aligned to the same samples. If donors are related or measurements are repeated, the independence assumption may fail and an appropriate mixed model is needed.

Test a Defined Family and Adjust for Multiplicity

A cis analysis typically tests many variants for a molecular trait. A small unadjusted p-value must therefore be interpreted in relation to the full family of hypotheses. Define that family before examining results. In this example, the family is 100 prespecified variants in the gene’s cis window, all tested for association with this one expression phenotype in the same 120 donors. If a study tests many genes as well, that broader family needs to be addressed separately or included in the stated correction.

For the Benjamini–Hochberg procedure, sort the \(m\) p-values from smallest to largest, \(p_{(1)}\le\cdots\le p_{(m)}\). At target false discovery rate \(q\), find the largest rank \(i\) satisfying \(p_{(i)}\le(i/m)q\); discoveries through that rank are declared. The adjusted p-value for rank \(i\) is the minimum, over ranks \(j\ge i\), of \(m p_{(j)}/j\), capped at 1. The procedure controls the expected false discovery proportion under independence and certain positive-dependence conditions; arbitrary dependence requires additional care.

Check Diagnostics Before Writing the Result

A significant test does not establish that a model fits adequately. Examine residuals against fitted values and covariates for nonlinearity, unequal spread, or systematic patterns; inspect a residual quantile plot for substantial departures from normality when using t-based inference; and check for influential observations. Also verify genotype counts, missingness, allele coding, and the covariate design. Diagnostics can identify concerns, but a clean-looking plot cannot prove that every assumption is true.

In the hypothetical analysis, genotype counts are 48 donors with dosage 0, 54 with dosage 1, and 18 with dosage 2. Suppose residual plots show no marked curvature or funnel shape, the quantile plot has only mild tail departure, and no single donor is highly influential. Suppose also that the maximum variance inflation factor among the model predictors is 1.8 and the design matrix is full rank. These checks support proceeding with the stated model, but the modest number of dosage-2 donors and the observational design remain relevant limitations. A serious failure of assumptions should prompt a justified alternative model or a qualified report, not simply a different test chosen to obtain significance.

A Short Justification for the Test and Interval

The genotype coefficient is estimated after accounting for the other model predictors. In ordinary least squares, the same coefficient can be understood by first removing the linear contributions of the covariates from both the phenotype and genotype, then regressing the residualized phenotype on residualized genotype. This makes clear why covariate specification affects the genotype estimate: the test uses the variation in genotype not explained by those covariates.

Under the stated linear-model assumptions, the estimated coefficient has standard error \(SE(\hat{\beta}_G)\). Under \(H_0\), the standardized estimate follows a t distribution with residual degrees of freedom, giving the statistic and p-value used above. Inverting the same two-sided test yields the confidence interval:

Thus the interval and test are consistent: a two-sided test at level \(\alpha=0.05\) rejects a zero coefficient exactly when the corresponding 95% interval excludes zero, subject to using the same model and standard-error method.

Report the Analysis, Not Just the Smallest P-Value

A concise association report should identify the phenotype and its scale, sample size, tested region and variant, genotype coding, covariates, effect estimate with units, standard error or interval, raw and adjusted p-values, and key diagnostics or limitations. State the correction family so readers can tell what the adjusted result means. Report the genome build and variant identifier, effect and other allele, and cis-window definition when those details are available.

Common Pitfalls

• Reporting only the smallest raw p-value. This omits how many hypotheses were tested. Define the family and report the relevant adjusted result.
• Treating transformation as a cosmetic choice. The coefficient’s units depend on the phenotype scale. Name the normalization and transformation so the effect is interpretable.
• Calling a lead variant causal. The lead is the strongest observed association in the tested set; linkage disequilibrium and untested variants can prevent identification of the causal variant.
• Assuming adjustment removes every confounder. Covariates address specified sources of variation under the model. Unmeasured confounding, misspecification, and technical artifacts may remain.
• Reading an adjusted p-value as a probability that the finding is false. It is a multiple-testing-adjusted measure used for decision-making, not a posterior probability about an individual result.