Skip to content

mecfs_bio.build_system.task.protein_coloc.coloc_abf

Colocalization of two traits at one locus by approximate Bayes factors (coloc.abf).

Follows Giambartolomei et al. 2014 and the coloc R package (approx.bf.estimates, combine.abf, sdY.est). The method assumes at most one causal variant per trait in the locus and needs no LD.

Glossary: - M: number of variants in the locus (both traits measured on the same M variants, same allele orientation). - z = beta / se, V = se^2. - W: prior standard deviation of the true effect. W = 0.15 * sdY for a quantitative trait (sdY is the trait's phenotypic standard deviation) and W = 0.2 for a case-control trait (log odds scale). - r = W^2 / (W^2 + V). - lABF: log approximate Bayes factor of a variant being causal versus no association, lABF = 0.5 * (log(1 - r) + r * z^2).

Hypotheses and their unnormalized log posteriors, with L1 = logsumexp(lABF1), L2 = logsumexp(lABF2) and L12 = logsumexp(lABF1 + lABF2): - H0, no association with either trait: 0. - H1, trait 1 only: log p1 + L1. - H2, trait 2 only: log p2 + L2. - H3, both traits, distinct causal variants: log p1 + log p2 + log(exp(L1 + L2) - exp(L12)). - H4, both traits, one shared causal variant: log p12 + L12. p1, p2 and p12 are the prior probabilities that a given variant is causal for trait 1 only, trait 2 only, or both.

When sdY of a quantitative trait is unknown it is estimated from the data: for a variant with minor allele frequency f, sample size n and standard error se, se^2 is approximately sdY^2 / (2 * n * f * (1 - f)). Rearranging, 2 * n * f * (1 - f) is approximately sdY^2 * (1 / se^2). So sdY^2 is the slope of the no-intercept regression

y_j = c * x_j + e_j,  with y_j = 2 * n_j * f_j * (1 - f_j) and x_j = 1 / se_j^2,

over variants j, where e_j is the regression error (the deviation of y_j from c * x_j). The least squares solution is c = sum_j(x_j * y_j) / sum_j(x_j^2), and sdY = sqrt(c).

Classes:

Functions:

  • combine_abf –
  • estimate_sd_y –

    Estimate sdY as sqrt(c), c = sum_j(x_j * y_j) / sum_j(x_j^2), x = 1 / se^2, y = 2 * n * f * (1 - f).

  • log_abf –

    Per-variant lABF = 0.5 * (log(1 - r) + r * z^2), r = W^2 / (W^2 + V).

Attributes:

TraitModel module-attribute 

TraitModel = QuantitativeTraitModel | CaseControlTraitModel

CaseControlTraitModel

A case-control trait with effects on the log odds scale.

Attributes:

case_fraction instance-attribute 

case_fraction: float

ColocAbf

coloc.abf: single causal variant per trait, no LD.

Methods:

Attributes:

name class-attribute instance-attribute 

name: str = 'coloc_abf'

priors class-attribute instance-attribute 

priors: ColocPriors = ColocPriors()

run 

run(data: ColocInput) -> ColocPosteriors
Source code in mecfs_bio/build_system/task/protein_coloc/coloc_abf.py 
227
228
229
230
231
232
233
234
def run(self, data: ColocInput) -> ColocPosteriors:
    w1 = _prior_sd(data.model1, data.se1, data.maf1, data.n1)
    w2 = _prior_sd(data.model2, data.se2, data.maf2, data.n2)
    return combine_abf(
        labf1=log_abf(data.beta1, data.se1, w1),
        labf2=log_abf(data.beta2, data.se2, w2),
        priors=self.priors,
    )

ColocInput

Summary statistics of two traits on the same M variants, in the same allele orientation.

maf and n of a trait are needed only when its model is quantitative with sd_y None.

Attributes:

beta1 instance-attribute 

beta1: ndarray

beta2 instance-attribute 

beta2: ndarray

maf1 instance-attribute 

maf1: ndarray | None

maf2 instance-attribute 

maf2: ndarray | None

model1 instance-attribute 

model1: TraitModel

model2 instance-attribute 

model2: TraitModel

n1 instance-attribute 

n1: ndarray | None

n2 instance-attribute 

n2: ndarray | None

se1 instance-attribute 

se1: ndarray

se2 instance-attribute 

se2: ndarray

ColocMethod

Bases: Protocol

Methods:

Attributes:

name property 

name: str

priors property 

priors: ColocPriors

run 

run(data: ColocInput) -> ColocPosteriors
Source code in mecfs_bio/build_system/task/protein_coloc/coloc_abf.py 
154
def run(self, data: ColocInput) -> ColocPosteriors: ...

ColocPosteriors

Attributes:

n_variants instance-attribute 

n_variants: int

pp_h0 instance-attribute 

pp_h0: float

pp_h1 instance-attribute 

pp_h1: float

pp_h2 instance-attribute 

pp_h2: float

pp_h3 instance-attribute 

pp_h3: float

pp_h4 instance-attribute 

pp_h4: float

ColocPriors

Attributes:

p1 class-attribute instance-attribute 

p1: float = 0.0001

p12 class-attribute instance-attribute 

p12: float = 1e-05

p2 class-attribute instance-attribute 

p2: float = 0.0001

QuantitativeTraitModel

A quantitative trait; sd_y None means estimate sdY from allele frequency and N.

Attributes:

sd_y instance-attribute 

sd_y: float | None

combine_abf 

combine_abf(
    labf1: ndarray, labf2: ndarray, priors: ColocPriors
) -> ColocPosteriors
Source code in mecfs_bio/build_system/task/protein_coloc/coloc_abf.py 
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
def combine_abf(
    labf1: np.ndarray, labf2: np.ndarray, priors: ColocPriors
) -> ColocPosteriors:
    l1 = float(logsumexp(labf1))
    l2 = float(logsumexp(labf2))
    l12 = float(logsumexp(labf1 + labf2))
    log_h = np.array(
        [
            0.0,
            np.log(priors.p1) + l1,
            np.log(priors.p2) + l2,
            np.log(priors.p1) + np.log(priors.p2) + _log_diff_exp(l1 + l2, l12),
            np.log(priors.p12) + l12,
        ]
    )
    pp = np.exp(log_h - logsumexp(log_h))
    return ColocPosteriors(
        pp_h0=float(pp[0]),
        pp_h1=float(pp[1]),
        pp_h2=float(pp[2]),
        pp_h3=float(pp[3]),
        pp_h4=float(pp[4]),
        n_variants=len(labf1),
    )

estimate_sd_y 

estimate_sd_y(
    se: ndarray, maf: ndarray, n: ndarray
) -> float

Estimate sdY as sqrt(c), c = sum_j(x_j * y_j) / sum_j(x_j^2), x = 1 / se^2, y = 2 * n * f * (1 - f).

Source code in mecfs_bio/build_system/task/protein_coloc/coloc_abf.py 
157
158
159
160
161
162
163
164
def estimate_sd_y(se: np.ndarray, maf: np.ndarray, n: np.ndarray) -> float:
    """Estimate sdY as sqrt(c), c = sum_j(x_j * y_j) / sum_j(x_j^2), x = 1 / se^2, y = 2 * n * f * (1 - f)."""
    x = 1 / se**2
    y = 2 * n * maf * (1 - maf)
    keep = np.isfinite(x) & np.isfinite(y)
    assert keep.any(), "no finite variants to estimate sd_y from"
    coefficient = float(np.dot(x[keep], y[keep]) / np.dot(x[keep], x[keep]))
    return float(np.sqrt(coefficient))

log_abf 

log_abf(
    beta: ndarray, se: ndarray, prior_sd: float
) -> np.ndarray

Per-variant lABF = 0.5 * (log(1 - r) + r * z^2), r = W^2 / (W^2 + V).

Source code in mecfs_bio/build_system/task/protein_coloc/coloc_abf.py 
179
180
181
182
183
184
def log_abf(beta: np.ndarray, se: np.ndarray, prior_sd: float) -> np.ndarray:
    """Per-variant lABF = 0.5 * (log(1 - r) + r * z^2), r = W^2 / (W^2 + V)."""
    v = se**2
    r = prior_sd**2 / (prior_sd**2 + v)
    z = beta / se
    return 0.5 * (np.log1p(-r) + r * z**2)