Chromosome-Specific Enrichment Analysis of Annotations Using Block Jackknife

Tests how strongly a set of significant variants overlaps a genomic annotation relative to the genome, with block-jackknife standard errors.

Overview¶

A variant set that overlaps an annotation more often than chance would predict is evidence that the annotation marks functional sequence. Counting the overlap is easy; putting an error bar on it is not, because variants are correlated along the genome and so cannot be resampled independently.

This module computes an odds ratio and an enrichment statistic for each annotation, then leaves out one chromosome at a time and recomputes, using the spread across those 22 leave-one-out estimates as a block-jackknife standard error. Blocking by chromosome keeps correlated variants together rather than splitting them across resamples.

When to run it. Run this module after chromosome-level enrichment inputs are available when a combined enrichment estimate and block-jackknife standard error are needed.

Definitions and test statistics¶

Odds Ratio (OR)¶

Strength of association between significant variants and an annotation:

OR = \frac{\left| AB \right| / \left| A \setminus B \right|}{\left| B \setminus A \right| / \left| \text{noA-noB} \right|}

(1)

where $A$ = SNPs in the annotation, $B$ = significant SNPs, $AB$ = their intersection, and noA-noB = SNPs in neither.

Enrichment¶

Whether an annotation contains a higher proportion of significant SNPs than expected by chance:

\text{Enrichment} = \frac{|AB|/|B|}{|A|/|\text{Target Set}|}

(2)

i.e. proportion of significant SNPs in the annotation divided by proportion of all SNPs in the annotation, where the Target Set is all SNPs in the genome or study region.

Standard Error (LOCO block jackknife)¶

Standard errors are estimated by removing one chromosome at a time and recomputing the statistic, capturing variability due to genomic structure. With $\theta_i$ the statistic ( $OR_i$ or $\text{Enrichment}_i$ ) excluding chromosome $i$ , $\bar{\theta}$ its mean, and $N$ the number of chromosomes (22 for autosomes):

\text{SE}(\theta) = \sqrt{\frac{\sum_{i=1}^{N} (\theta_i - \bar{\theta})^2}{N \cdot (N-1)}}

(3)

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import PathPatch, Circle, FancyBboxPatch
from matplotlib.path import Path

def make_lens_path(c1, c2, r, n=400):
    x1, y1 = c1; x2, y2 = c2
    d = np.hypot(x2 - x1, y2 - y1)
    if d >= 2 * r:
        return None
    a   = np.arccos(d / (2 * r))
    phi = np.arctan2(y2 - y1, x2 - x1)
    t1  = np.linspace(phi - a, phi + a, n)
    t2  = np.linspace(phi + np.pi + a, phi + np.pi - a, n)
    verts = (list(zip(x1 + r*np.cos(t1), y1 + r*np.sin(t1))) +
             list(zip(x2 + r*np.cos(t2), y2 + r*np.sin(t2))) +
             [(x1 + r*np.cos(phi - a), y1 + r*np.sin(phi - a))])
    codes = [Path.MOVETO] + [Path.LINETO]*(2*n - 1) + [Path.CLOSEPOLY]
    return Path(verts, codes)

fig, ax = plt.subplots(figsize=(14, 11))
ax.set_xlim(0, 14); ax.set_ylim(0, 11); ax.axis('off')
fig.patch.set_facecolor('white')

# Background box (C)
bg = FancyBboxPatch((0.2, 3.1), 13.6, 7.5, boxstyle='round,pad=0.15',
                    facecolor='#F0F3F7', edgecolor='#99A9BB',
                    linewidth=2.0, zorder=0)
ax.add_patch(bg)

ax.text(7.0, 10.38, r'$\mathbf{C}$ : Target Set (All SNPs)',
        ha='center', va='center', fontsize=17, fontweight='bold', color='#1C2833')
ax.text(0.55, 9.85, r'$|C|$' + '\nAll SNPs',
        ha='left', va='top', fontsize=12, color='#707B7C', linespacing=1.5)

# Circles
cx_a, cy, r = 4.8, 6.7, 2.75
cx_b = 9.2
circ_a = Circle((cx_a, cy), r, facecolor='#AED6F1', edgecolor='#2471A3',
                linewidth=3.0, alpha=0.70, zorder=1)
circ_b = Circle((cx_b, cy), r, facecolor='#FADBD8', edgecolor='#C0392B',
                linewidth=3.0, alpha=0.70, zorder=1)
ax.add_patch(circ_a); ax.add_patch(circ_b)

lens = PathPatch(make_lens_path((cx_a, cy), (cx_b, cy), r),
                 facecolor='#F0B429', edgecolor='none', alpha=0.90, zorder=2)
ax.add_patch(lens)

ax.text(cx_a - 1.0, cy + r*0.78, r'$\mathbf{A}$ : Annotated SNPs',
        ha='center', va='center', fontsize=14, fontweight='bold', color='#1A5276')
ax.text(cx_b + 1.0, cy + r*0.78, r'$\mathbf{B}$ : Significant SNPs',
        ha='center', va='center', fontsize=14, fontweight='bold', color='#7B241C')

mid_x = (cx_a + cx_b) / 2
ax.text(cx_a - 1.55, cy, r'$|A \setminus B|$',
        ha='center', va='center', fontsize=13, color='#1A5276', style='italic')
ax.text(mid_x, cy, r'$|A \cap B|$',
        ha='center', va='center', fontsize=14, fontweight='bold', color='#6E2F05', style='italic')
ax.text(cx_b + 1.55, cy, r'$|B \setminus A|$',
        ha='center', va='center', fontsize=13, color='#7B241C', style='italic')

ax.text(13.05, 9.1, r'$|\neg A \cap \neg B|$',
        ha='right', va='center', fontsize=13, color='#4A4A4A', style='italic')

# Formula box
ax.add_patch(FancyBboxPatch((0.3, 0.08), 13.4, 2.75,
                             boxstyle='round,pad=0.12',
                             facecolor='#FDFEFE', edgecolor='#BDC3C7',
                             linewidth=1.4, zorder=0))

# OR: formula centered, explanation centered below
ax.text(7.0, 2.38,
        r'$\mathrm{Odds\ Ratio} = \dfrac{|A \cap B|\;/\;|A \setminus B|}'
        r'{|B \setminus A|\;/\;|\neg A \cap \neg B|}'
        r'\ =\ \dfrac{|A \cap B| \times |\neg A \cap \neg B|}'
        r'{|A \setminus B| \times |B \setminus A|}$',
        ha='center', va='center', fontsize=12.5, color='#1C2833')
ax.text(7.0, 1.72,
        'odds of overlap in annotated SNPs  /  odds of overlap in non-annotated SNPs',
        ha='center', va='center', fontsize=11, color='#7F8C8D', style='italic')

# Divider
ax.plot([0.6, 13.4], [1.28, 1.28], color='#D5D8DC', linewidth=0.8)

# Enrichment: formula centered, explanation centered below
ax.text(7.0, 0.92,
        r'$\mathrm{Enrichment} = \dfrac{|A \cap B|\;/\;|B|}{|A|\;/\;|C|}$',
        ha='center', va='center', fontsize=12.5, color='#1C2833')
ax.text(7.0, 0.35,
        'fraction of significant SNPs in annotation  /  background rate of annotation in all SNPs',
        ha='center', va='center', fontsize=11, color='#7F8C8D', style='italic')

plt.tight_layout(pad=0.2)
plt.show()

When to run it. After a significant variant set exists and an annotation is available to test it against.

Input¶

--significant-variants-path: the variant set being tested, two columns chr and pos with a header. Required. Example tests/fixtures/eoo_enrichment/protocol_example.eoo_significant_variants.tsv.gz:
```
chr  pos
1    100015157
1    1000291
```
--baseline-anno-path: the annotation table covering the genome, one row per variant, with the annotation columns starting at --annotations-start. Required. Example input/enrichment/protocol_example.eoo_baseline_annotation.tsv:
```
chr  BP        ref  alt  rsid    maf     ANNO_baseline
1    10000654  A    C    rs1_0   0.3734  0
1    10100228  C    T    rs1_1   0.0602  0
```
--annotations-start: the 1-based column index at which the annotation columns begin, 7 by default. Everything left of it is treated as variant identity and frequency.
--trait and --annotation-name: labels, both required. They carry no analytic meaning; the output file is named {trait}.{annotation-name}.enrichment_results.rds.
--cwd: the directory outputs are written to.
--numThreads: threads for the jackknife, 8 by default.

Output¶

{cwd}/enrichment/{trait}.{annotation-name}.enrichment_results.rds - a list of 10 holding the point estimates, their block-jackknife standard errors and the per-chromosome estimates they derive from. Example output/eoo_enrichment/enrichment/protocol_example.baseline.enrichment_results.rds:
```
List of 10
 $ summary                 :'data.frame':  2 obs. of  13 variables:
  ..$ Annotation           : chr [1:2] "maf" "ANNO_baseline"
  ..$ OR                   : num [1:2] 1 1
  ..$ OR_SE                : num [1:2] 0 0
  ..$ Enrichment           : num [1:2] 1 1
  ..$ Enrichment_SE        : num [1:2] 0 0
 $ OR_blockJacknife        : num [1:22, 1:2] 1 1 1 1 1 1 1 1 1 1 ...
 $ Enrichment_blockJacknife: num [1:22, 1:2] 1 1 1 1 1 1 1 1 1 1 ...
 $ OR                      : Named num [1:2] 0 0
 $ Enrichment              : Named num [1:2] 1 1
```
summary carries one row per annotation column, here maf and ANNO_baseline. The two blockJacknife matrices are 22 by 2: one row per chromosome left out, one column per annotation, and the spread down each column is what the standard errors are computed from. On this toy annotation every variant scores 0, so the odds ratios come back at 1 with zero standard error; that is the expected result for an annotation with no signal, not a failed run.

Minimal Working Example¶

One workflow, enrichment, run once per annotation. --trait and --annotation-name are labels only; together they name the output file.

Timing: TBD (on toy dataset)

sos run pipeline/eoo_enrichment.ipynb enrichment \
    --significant_variants_path tests/fixtures/eoo_enrichment/protocol_example.eoo_significant_variants.tsv.gz \
    --baseline_anno_path tests/fixtures/eoo_enrichment/protocol_example.eoo_baseline_annotation.tsv.gz \
    --trait protocol_example \
    --annotation-name baseline \
    --cwd output/eoo_enrichment

Command Interface¶

sos run pipeline/eoo_enrichment.ipynb -h

[91mERROR[0m: [91mNotebook JSON is invalid: %s[0m
usage: sos run code/SoS/enrichment/eoo_enrichment.ipynb
               [workflow_name | -t targets] [options] [workflow_options]
  workflow_name:        Single or combined workflows defined in this script
  targets:              One or more targets to generate
  options:              Single-hyphen sos parameters (see "sos run -h" for details)
  workflow_options:     Double-hyphen workflow-specific parameters

Workflows:
  enrichment

Global Workflow Options:
  --cwd output (as path)
                        Path to the work directory of the analysis.
  --modular-script-dir code/script (as path)
  --significant-variants-path VAL (as path, required)
  --baseline-anno-path VAL (as path, required)
  --numThreads 8 (as int)
                        Number of threads
  --trait VAL (as str, required)
                        For cluster jobs, number commands to run per job
  --annotation-name VAL (as str, required)
  --job-size 1 (as int)
  --walltime 12h
  --mem 16G

Sections
  enrichment:
    Workflow Options:
      --annotations-start 7 (as int)

Workflow implementation¶

[global]
# Path to the work directory of the analysis.
parameter: cwd = path('output')
parameter: modular_script_dir = path('code/script')  # override with --modular-script-dir

parameter: significant_variants_path = path
parameter: baseline_anno_path = path
# Number of threads
parameter: numThreads = 8
# For cluster jobs, number commands to run per job
parameter: trait = str
parameter: annotation_name = str
name = f"{trait}.{annotation_name}"
parameter: job_size = 1
parameter: walltime = '12h'
parameter: mem = '16G'

[enrichment]
parameter: annotations_start = 7
output: enrichment = f'{cwd:a}/{step_name}/{name}.enrichment_results.rds'

task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]:bnn}'
bash: expand = '${ }', stderr = f'{_output[0]}.stderr', stdout = f'{_output[0]}.stdout'
    Rscript ${modular_script_dir}/enrichment/eoo_enrichment.R \
        --significant-variants-path ${significant_variants_path} \
        --baseline-anno-path ${baseline_anno_path} \
        --annotations-start ${annotations_start} \
        --output ${_output['enrichment']}