Constructing and Reading GWAS QQ Plots

Why Plot Observed Against Expected Null P-Values?

A genome-wide association study can produce hundreds of thousands or millions of p-values. A QQ plot compresses their overall pattern into one graph: it compares the ordered observed p-values with reference values expected under a specified null model. The plot can reveal whether the distribution broadly follows the reference or departs from it, especially among the smallest p-values.

A QQ plot is a diagnostic summary, not a list of discoveries and not a map of the genome. Sorting the p-values removes their variant identities. To find which genomic regions contribute to a departure, the ranks must be connected back to variant-level results and genomic coordinates.

Definition: Let $p_1,\ldots,p_m$ be the p-values for the $m$ variants in a defined testing family, and write their ordered values as $p_{(1)}\leq p_{(2)}\leq\cdots\leq p_{(m)}$. In a conventional GWAS QQ plot, the horizontal coordinate for rank $i$ is the transformed reference plotting position $-\log_{10}\{i/(m+1)\}$, and the vertical coordinate is the observed transformed p-value $-\log_{10}\{p_{(i)}\}$. Smaller p-values therefore appear higher on the graph.

The reference plotting position $i/(m+1)$ is the expected value of the $i$th ordered p-value under independent, continuous Uniform(0,1) null p-values. Its negative logarithm is obtained by transforming that expected p-value. It is not, in general, the expected value of the negative-log-transformed order statistic. This distinction matters because a nonlinear transformation does not commute with taking an expectation.

Constructing the Plot and Reading Its Reference

First define the testing family: decide which variant-level p-values belong in the plot, and use the same family size $m$ throughout. Then sort those p-values from smallest to largest. For each rank $i$, calculate the reference plotting position and the observed transformed p-value. Plot the reference quantity horizontally and the observed quantity vertically. With this convention, points near the line $y=x$ follow the reference pattern; points above it have smaller observed p-values than the reference plotting position.

Under the ideal global null, every tested variant is null and its p-value is Uniform(0,1). If the p-values are also independent, the ordered p-values have a standard order-statistic distribution. Real GWAS p-values may be only approximately uniform under null hypotheses, and variants in linkage disequilibrium are correlated. Thus, the line is a useful reference, but a confidence envelope based on independent p-values may not be well calibrated for a real genome-wide analysis.

Worked Example: Calculating Two QQ-Plot Coordinates

Suppose a hypothetical analysis tests $m=1{,}000$ variants. The p-value at rank $i=500$ is $p_{(500)}=0.45$. Calculate its conventional coordinates.

The reference plotting position is $500/(1{,}000+1)=0.4995$. Its transformed reference plotting coordinate is

x_{500}=-\log_{10}(0.4995)=0.3015.

The observed coordinate is

y_{500}=-\log_{10}(0.45)=0.3468.

Thus, this point is at approximately $(0.3015,0.3468)$, slightly above the diagonal. The observed p-value at this rank is somewhat smaller than the reference plotting position. This one point does not establish a systematic departure; the overall pattern and a suitable envelope provide context.

What Is Expected on the Negative-Log Scale?

The conventional plotting position is convenient and widely used, but it should not be described as the exact expected coordinate on the negative-log scale. Under independent Uniform(0,1) null p-values, $P_{(i)}$, the $i$th ordered p-value, has a Beta distribution with shape parameters $i$ and $m+1-i$. Its mean is $i/(m+1)$. The expected transformed value instead depends on the full distribution of $P_{(i)}$.

Result: If the p-values are independent Uniform(0,1) variables, then $P_{(i)}\sim\operatorname{Beta}(i,m+1-i)$, and $$ E[-\log_{10}\{P_{(i)}\}]=\frac{\psi(m+1)-\psi(i)}{\ln(10)} =\frac{H_m-H_{i-1}}{\ln(10)}. $$ Here $\psi$ is the digamma function and $H_j=\sum_{k=1}^{j}1/k$, with $H_0=0$. This exact expected transformed value generally differs from $-\log_{10}\{i/(m+1)\}$, the conventional transformed reference plotting position.

The distinction follows directly from the Beta density. For $a=i$ and $b=m+1-i$, the density is proportional to $u^{a-1}(1-u)^{b-1}$. Differentiating the beta-function integral with respect to $a$ brings down a factor of $\ln(u)$, so

E[\ln\{P_{(i)}\}] =\frac{\partial}{\partial a}\ln B(a,b) =\psi(a)-\psi(a+b).

Negating and dividing by $\ln(10)$ gives the stated expected negative-log value. Since $\psi(j)=H_{j-1}-\gamma$ for positive integers $j$, where $\gamma$ is Euler's constant, the harmonic-number form follows. The conventional reference coordinate remains a transformed expected p-value, not this expectation of a transformed p-value.

Worked Example: The Smallest P-Value and the Two Reference Quantities

Consider $m=100$ independent, continuous Uniform(0,1) null p-values. For the smallest p-value, $i=1$, the expected p-value is $1/101=0.009901$. The conventional transformed reference plotting position is

-\log_{10}(1/101)=2.004.

But the expected transformed smallest p-value is $H_{100}/\ln(10)$. Using $H_{100}\approx5.1874$,

E[-\log_{10}\{P_{(1)}\}] =\frac{5.1874}{2.3026} \approx2.253.

The two values, 2.004 and 2.253, are not contradictory: one transforms the expected p-value, while the other averages the transformed p-value across its null distribution. A QQ plot using $i/(m+1)$ uses the first quantity as its plotting convention.

Pointwise Confidence Envelopes

A confidence envelope gives a range of null values against which to compare the observed order statistics. Under the independent Uniform(0,1) model, $P_{(i)}\sim\operatorname{Beta}(i,m+1-i)$. Let $q_{i,L}$ and $q_{i,U}$ be the lower and upper $\alpha/2$ quantiles of this Beta distribution. A central $1-\alpha$ pointwise interval for $P_{(i)}$ is $[q_{i,L},q_{i,U}]$. Because $-\log_{10}(p)$ decreases as $p$ increases, the corresponding interval on the vertical axis is reversed:

[-\log_{10}(q_{i,U}),\;-\log_{10}(q_{i,L})].

“Pointwise” means that the interval has the stated coverage for one specified rank $i$, under the model. It does not mean that the entire collection of plotted ranks will lie inside all the intervals with probability $1-\alpha$. A simultaneous envelope requires a procedure designed to control coverage across ranks, such as calibrating a maximum-deviation statistic under the null.

Worked Example: Checking the Minimum Against a Pointwise Envelope

Continue with $m=100$ independent Uniform(0,1) null p-values. For rank $i=1$, the distribution is $\operatorname{Beta}(1,100)$, whose cumulative distribution function is $F(u)=1-(1-u)^{100}$. Its quantile at probability $r$ is $1-(1-r)^{1/100}$. The central 95% interval for the minimum is therefore

q_{1,L}=1-0.975^{1/100}\approx0.0002532,\qquad q_{1,U}=1-0.025^{1/100}\approx0.03622.

Suppose the hypothetical observed minimum is $p_{(1)}=0.004$. Its transformed value is $-\log_{10}(0.004)=2.398$. The transformed pointwise interval is

[-\log_{10}(0.03622),\;-\log_{10}(0.0002531)] \approx[1.441,3.597].

Since $2.398$ lies inside this interval, the minimum is not outside this rank-specific 95% envelope. This does not show that the whole QQ plot is compatible with the global null, nor does it account for dependence between variants.

Reading Departures and Tracing Them to Regions

A QQ plot displays p-value ranks, not variant locations. A sustained departure above the diagonal means that observed p-values are smaller than the reference values across a range of ranks. Points below the diagonal indicate larger p-values than the reference at those ranks. A few extreme points can create a sharp upper-tail departure, while a broad departure over many ranks reflects a wider distributional difference. The plot alone does not determine why a departure occurs; possible explanations include genuine association, model misspecification, or imperfect null calibration.

To identify which genomic regions contribute to an upper-tail departure, retain the variant identifiers and coordinates alongside the p-values. After plotting, map the ranks of interest back to those records, inspect the corresponding genomic positions, and use a genome-position display or regional association results to see where the small p-values cluster. A cluster can show that several top-ranked variants lie in one region, but the QQ plot itself cannot show that clustering because its ranks have been detached from genomic order.

Worked Example: Mapping an Upper-Tail Departure to Genomic Positions

A hypothetical study tests $m=10{,}000$ variants. Its three smallest p-values are $1.0\times10^{-9}$, $3.0\times10^{-8}$, and $2.0\times10^{-7}$. Suppose the variant records show that the first two are on chromosome 6 within one region and the third is on chromosome 2 elsewhere. Calculate the first two reference coordinates and compare them with the observed values.

For rank 1, the reference plotting position is $1/10{,}001\approx0.00009999$, so $x_1\approx4.000$. The observed coordinate is $-\log_{10}(10^{-9})=9.000$. For rank 2, the reference position is $2/10{,}001\approx0.00019998$, so $x_2\approx3.699$; its observed coordinate is $-\log_{10}(3.0\times10^{-8})\approx7.523$.

Both points are far above the diagonal, and the variant records show that these two extreme ranks come from the same chromosome 6 region. The third-ranked p-value has observed coordinate $-\log_{10}(2.0\times10^{-7})\approx6.699$, but its variant is in a different region. This illustrates how rank-level departures can be connected to genomic regions only by returning to the variant-level results. It does not establish that any region is causal or explain why its p-values are small.

For actual GWAS data, independent-null Beta envelopes may be inappropriate because nearby variants are correlated through linkage disequilibrium. One way to obtain a more suitable reference is to simulate or permute under a valid null procedure that preserves the relevant dependence and repeats the same analysis pipeline. The resulting simulated QQ plots can be used to assess rank-wise or simultaneous departures. Such calibration depends on the validity of the null procedure; an envelope cannot correct an invalid association model simply by being drawn around the observed points.

Common Mistakes

Calling the conventional reference coordinate an exact expected transformed value. The expression $-\log_{10}\{i/(m+1)\}$ transforms the expected p-value. In general it differs from $E[-\log_{10}\{P_{(i)}\}]$.
Reversing the vertical envelope bounds. Since negative logarithms decrease as p-values increase, the upper p-value quantile gives the lower transformed bound.
Treating pointwise limits as a simultaneous guarantee. Each pointwise interval addresses one rank under its model; it does not guarantee that every rank lies within its interval together.
Assuming the independent-null envelope automatically fits a GWAS. Linkage disequilibrium induces dependence among variant test results, changing the joint distribution of the ordered p-values.
Reading genomic location from a QQ plot. Sorting discards location. Use the variant identifiers and coordinates to find the regions contributing to selected ranks.
Interpreting departure as a specific cause. A QQ pattern is a diagnostic signal, not proof of association, confounding, or any other single explanation.

Key takeaway: A conventional GWAS QQ plot compares observed ordered p-values with transformed expected-null plotting positions. Keep that plotting convention distinct from the exact expected value on the transformed scale, interpret envelopes under their stated dependence assumptions, and use variant-level coordinates to locate genomic contributions.

Check Your Understanding

Answer each question using the plotting convention and assumptions stated in the tutorial.

For $m=499$, what is the conventional reference plotting position for rank $i=10$, and which quantity does it represent before transformation?
Why is $-\log_{10}\{i/(m+1)\}$ generally not equal to $E[-\log_{10}\{P_{(i)}\}]$?
Why are the transformed bounds from a p-value interval $[q_{i,L},q_{i,U}]$ written in reversed order?
What does a point above the QQ reference line indicate about its observed p-value relative to its reference plotting position?
Why might a pointwise Beta envelope based on independent null p-values be unreliable for correlated GWAS variants?
What additional information is needed to identify which genomic regions contribute to the smallest p-values?