Variant Harmonization for Colocalization

Why Harmonization Comes Before Colocalization

Colocalization asks whether two traits have evidence for a shared causal variant in a region. That comparison only makes sense if a variant in one study is matched to the same biological variant in the other. Two records can have similar identifiers but refer to different alleles or coordinates; conversely, the same variant can have different identifiers or allele labels in different files.

Harmonization is the process of making variant records comparable across studies before combining their evidence. It includes checking genome build, chromosome and position, reference and alternate alleles, effect-allele coding, identifiers, and the region and variants actually observed. Harmonization does not itself establish that two traits share a cause. It makes the variant-by-variant comparison interpretable.

Definition: A harmonized variant is a variant identified consistently across data sources, with compatible genome build, genomic location, allele representation, and effect-allele orientation. A practical variant key is build, chromosome, position, reference allele, and alternate allele, for example, GRCh38:7:100100:A:G. An rsID can be stored as an additional identifier, but should not replace checks of the underlying location and alleles.

The effect allele is the allele for which an association estimate is reported. If a study reports a regression coefficient \(\beta\), the coefficient is the estimated change in the trait per additional copy of that effect allele under the study’s genotype coding and model. The other allele is the alternative at that biallelic site. When two studies use opposite effect alleles, their coefficients have opposite signs even when they describe the same underlying association.

Allele Orientation and Strand Checks

For the same biallelic variant, changing the counted allele reverses the dosage: if the original dosage is \(G\), the opposite-allele dosage is \(2-G\). In a regression with an intercept, replacing \(G\) by \(2-G\) reverses the coefficient’s sign, while leaving its standard error and two-sided association test unchanged. The frequency of the newly counted allele is one minus the original frequency.

\beta_{\text{new}}=-\beta_{\text{old}},\qquad SE_{\text{new}}=SE_{\text{old}},\qquad f_{\text{new}}=1-f_{\text{old}}.

The frequency relation assumes the two alleles are exhaustive at a biallelic site and that frequencies refer to the same sample and variant. It is a useful check, not a substitute for correct identification. A frequency discrepancy can arise from different ancestry composition, imputation or quality filters, or a variant mismatch, as well as from a mistaken allele orientation.

Strand complements require care. A and G on one strand correspond to T and C on the opposite strand; their labels may therefore differ while the biological variant is the same. A/T and C/G variants are palindromic: complementing their alleles does not distinguish the two strand orientations. Allele frequency may help resolve a palindromic variant when the frequency is clearly away from 0.5 and populations and frequency definitions are comparable. If orientation remains ambiguous, exclude the variant from allele-sensitive comparisons rather than guessing. In many colocalization implementations, effect direction is not the primary evidence: association strength and variant identity drive the comparison. Even so, unresolved allele orientation is a warning that records may have been paired incorrectly, and it invalidates directional interpretation.

Worked Example: Reorienting Effect Alleles

Two hypothetical studies report the same confirmed biallelic variant on the same genome build. The GWAS reports effect allele A, other allele G, \(\beta=0.12\), \(SE=0.04\), and effect-allele frequency 0.30. The eQTL study reports effect allele G, other allele A, \(\beta=0.08\), \(SE=0.03\), and G frequency 0.68. Put both estimates on the A-allele scale and check the frequencies.

The GWAS is already on the A-allele scale, so its coefficient remains \(0.12\), its standard error remains \(0.04\), and its A frequency is 0.30. The eQTL coefficient is reported per G allele. Changing to the A-allele scale reverses its sign:

\beta_{\text{eQTL, per A}}=-0.08,\qquad SE_{\text{eQTL}}=0.03.

The eQTL study reports a G frequency of 0.68, so the corresponding A frequency is \(1-0.68=0.32\). This is reasonably close to the GWAS A frequency of 0.30, which is a compatible check in this hypothetical example; it does not prove the variant match. Before comparing, the analyst should still verify the build, position, and alleles. The GWAS \(z\)-statistic is \(0.12/0.04=3.00\); the reoriented eQTL \(z\)-statistic is \(-0.08/0.03=-2.67\). Reorientation changes the sign, not the magnitude of the eQTL association statistic.

One short derivation explains why a simple allele swap does not change association strength. Let \(G\) be the dosage of the original effect allele and \(Y=a+\beta G+\epsilon\). The opposite-allele dosage is \(G'=2-G\), so \(G=2-G'\). Substitution gives \(Y=(a+2\beta)-\beta G'+\epsilon\): the intercept changes, the coefficient changes sign, and the residual variation is unchanged. Under the same regression and standard-error calculation, the standard error is therefore unchanged as well.

Identifiers, Coordinates, and Allele Representation

A variant identifier alone is not a sufficient join key. rsIDs can be missing, updated, associated with multiple representations, or used inconsistently in files. Coordinates are also meaningful only relative to a named genome assembly. A position from one build should not be compared directly with a position from another build as if both refer to the same base. If a study’s build is unknown, resolve it from metadata or recover the original reference before proceeding.

For single-nucleotide variants, a useful audit compares build, chromosome, position, reference allele, alternate allele, effect allele, and other allele. For insertions and deletions, representations can differ even when they describe the same sequence change. Normalize indels against the chosen reference genome, including consistent trimming and left alignment, before comparing records. Multiallelic sites need additional care: each alternate allele can define a distinct association, and a record for one alternate allele must not be silently matched to another.

Worked Example: A Misleading Identifier Join

Suppose two hypothetical data sets each contain records labelled “marker-X.” In study 1, marker-X is GRCh37:7:100100:A:G; in study 2, the same label is attached to GRCh38:7:100250:C:T. The two position numbers differ, but that alone does not establish a mismatch because the genome builds differ. The analyst must map coordinates to a common build and compare the resulting alleles and location. If the mapped records identify different biological variants, the shared label must not be used to pair their association evidence.

To see how a mistaken pairing can affect the result, consider a simplified two-variant colocalization calculation. Assume one causal variant per trait, prior probabilities per variant \(p_1=p_2=0.01\) for association with either trait, and \(p_{12}=0.001\) for association with both. Let the trait-specific Bayes factors be 10 and 2 for trait 1 at variants A and B, and 2 and 10 for trait 2 at those same variants. The evidence contributions are

H_1=p_1\sum_i BF_{1i},\quad H_2=p_2\sum_i BF_{2i},\quad H_3=p_1p_2\sum_{i\ne j}BF_{1i}BF_{2j},\quad H_4=p_{12}\sum_i BF_{1i}BF_{2i}.

These are unnormalized hypothesis weights; \(H_0=1\). With correct matching, \(H_1=0.01(10+2)=0.12\) and \(H_2=0.01(2+10)=0.12\). The distinct-variant evidence is \(H_3=0.0001[(10)(10)+(2)(2)]=0.0104\), while the shared-variant evidence is \(H_4=0.001[(10)(2)+(2)(10)]=0.0400\). The total weight is \(1+0.12+0.12+0.0104+0.0400=1.2904\), giving \(P(H_4\mid\text{data})=0.0400/1.2904=0.0310\) in this toy calculation.

Now suppose an erroneous join swaps the trait 2 Bayes factors between A and B, making them appear to be 10 and 2. Then \(H_3=0.0001[(10)(2)+(2)(10)]=0.00400\), and \(H_4=0.001[(10)(10)+(2)(2)]=0.104\). The total weight becomes 1.348, and the resulting \(H_4\) posterior is \(0.104/1.348=0.0772\). The apparent support for sharing has more than doubled, solely because evidence was attached to the wrong variant. The numbers are hypothetical and illustrate the direction of the error, not a general expected change.

Regional Coverage and the Variant Set

Colocalization compares evidence over a specified set of variants. Before interpreting its result, check that both studies cover the same genomic region, that the boundaries include the signal in each study, and that the analysis uses a clearly defined set of variants available and adequately measured in both. Compare the number of variants reported, the number retained after quality control, and the number in the shared analysis set. A variant absent from a file is not a measured variant with no association.

Variant coverage can be incomplete because of different region boundaries, genotyping or imputation filters, allele-frequency thresholds, or differences in how variants are represented. If the strongest variant for one trait is missing from the shared set, the analysis cannot directly compare its evidence with the other trait’s evidence. A restricted set may therefore weaken support for a true shared signal, or alter the balance of evidence among hypotheses. Record what was removed and why; do not interpret a result as if both studies tested every relevant variant.

Worked Example: Missing the Strongest Shared Variant

Consider a hypothetical region with two harmonized variants, A and B. For each trait, suppose the Bayes factors are 20 at A and 1 at B. Use \(p_1=p_2=0.01\), \(p_{12}=0.001\), and \(H_0=1\), with the same single-causal-variant hypothesis weights defined above.

With both variants included, \(H_1=H_2=0.01(20+1)=0.21\). The distinct-variant weight is \(H_3=0.0001[(20)(1)+(1)(20)]=0.004\), and the shared-variant weight is \(H_4=0.001[(20)(20)+(1)(1)]=0.401\). The total is \(1+0.21+0.21+0.004+0.401=1.825\), so \(P(H_4\mid\text{data})=0.401/1.825=0.220\).

Now suppose a region boundary or filtering step removes A from the shared variant set, leaving only B. In the restricted analysis, \(H_1=0.01\), \(H_2=0.01\), \(H_3=0\), and \(H_4=0.001\). The \(H_3\) sum is zero because there is no pair of distinct variants in a one-variant set. The total is 1.021, and \(P(H_4\mid\text{data})=0.001/1.021=0.000979\). This sharp reduction reflects loss of the informative variant, not evidence that the missing variant had no association. The example shows why variant coverage and the analyzed set must accompany a colocalization result.

A Practical Harmonization Audit

A reproducible audit separates record identity, allele orientation, and coverage checks. Each issue should be resolved or documented before evidence is combined.

Confirm the reference build.
Record the genome assembly for each file. Convert coordinates to a common build using an appropriate reference-based mapping procedure before comparing positions.

Match variant identity.
Compare chromosome, position, and reference and alternate alleles. Normalize indels and distinguish alternate alleles at multiallelic sites. Use rsIDs as supporting metadata, not as the only match criterion.

Align effect alleles.
State which allele each coefficient and frequency refer to. Swap coefficient signs and replace allele frequencies by one minus the original frequency only when changing to the opposite allele at the same confirmed biallelic variant. Resolve strand orientation; exclude unresolved palindromic variants when necessary.

Inspect regional coverage.
Compare region boundaries, reported variant counts, quality-control exclusions, and the shared variant set. Identify any missing lead or otherwise important variants rather than treating absence as a null result.

Recheck before interpreting.
Confirm that each merged row refers to the same variant in both studies, that allele-frequency checks are plausible, and that the analysis reports the variant set on which its colocalization evidence is based.

Common Mistake and Key Takeaway

A common mistake is to assume that a matching rsID, a matching position number, or a matching effect-allele letter is enough to establish a match. For example, “both files say A” does not confirm that the A allele is on the same strand, at the same build-specific location, or paired with the same other allele. The correction is to verify the build and variant identity first, then orient the effect allele and check frequencies. If identity or orientation remains uncertain, do not force a match.

A related mistake is to think that flipping a coefficient automatically changes the evidence for colocalization. For a correctly identified biallelic variant, changing the counted allele reverses the coefficient’s sign but preserves its association strength under the same model. By contrast, pairing a coefficient or Bayes factor with the wrong variant, or omitting an important variant from the analyzed region, changes which evidence is compared and can alter the relative support for shared and distinct signals.

Key takeaway: Before interpreting colocalization, verify that both studies use a common, accurately defined variant set. Check build, coordinates, alleles, effect-allele orientation, identifiers, and regional coverage; document unresolved variants instead of treating them as reliable matches.

Check Your Understanding

Answer these questions using the definitions and checks above.

A study reports a T-allele frequency of 0.27 at a confirmed biallelic variant. What is the frequency of the other allele, and what additional information would you check before comparing it with another study?
One study reports a coefficient per C allele and another per G allele at a confirmed C/G variant. What happens to the coefficient sign, standard error, and allele frequency when the second estimate is expressed per C allele?
Why is a matching rsID insufficient by itself to establish that two summary-statistic rows refer to the same variant?
A variant is absent from one study’s file because it lies outside that study’s reported region. Is it valid to treat that variant as having no association in that study? Explain.
How can a mistaken swap of variant-level evidence change the balance between \(H_3\) and \(H_4\)?