bio-atac-seq-enhancer-gene-linking

What this skill does

## Version Compatibility

Reference examples tested with: ABC-Enhancer-Gene-Prediction 0.2.2+ (Engreitz lab), ENCODE-rE2G v1.0+ (2024 release), Cicero 1.20+, GenomicInteractions 1.36+, FitHiChIP 9.1+, hicpro 3.1+, FAN-C 0.9+, MACS3 3.0+, samtools 1.19+, bedtools 2.31+.

Verify before use:
- CLI: `<tool> --version` then `<tool> --help` to confirm flags
- R: `packageVersion('<pkg>')` then `?function_name` to verify parameters
- Python: `pip show <package>` then `help(module.function)` to check signatures

If code throws unexpected errors, introspect the installed package and adapt rather than retrying.

# Enhancer-Gene Linking

**"Which gene does this distal accessible region regulate?"** -> Predict the enhancer's target gene using a model that combines accessibility activity, 3D contact frequency, and (optionally) sequence-based chromatin predictions. Output is a per-(enhancer, gene) score that can be thresholded for high-confidence calls.

- CLI: ABC pipeline (`run.neighborhoods.py`, `predict.py` from Engreitz lab)
- CLI: ENCODE-rE2G (Snakemake-based; ENCODE 4 standard 2024)
- R: Cicero (ATAC-only; covered in atac-seq/co-accessibility)
- CLI: FitHiChIP / hichipper for HiChIP H3K27ac loops
- Database: EpiMap (Boix 2021), GeneHancer, FANTOM5 (pre-computed reference)

ABC and ENCODE-rE2G are the canonical predictors when Hi-C/Micro-C data is available. Cicero is the ATAC-only fallback. CRISPRi-FlowFISH (Fulco 2019) is the gold-standard experimental validation.

## Algorithmic Taxonomy

| Method | Inputs | Mathematics | Strength | Fails when |
|--------|--------|-------------|----------|------------|
| ABC (Fulco 2019, Nasser 2021) | ATAC + H3K27ac + Hi-C/Micro-C | ABC = (Activity_E x Contact_E,G) / sum_e(Activity_e x Contact_e,G); threshold typically >= 0.02 | Mechanistically grounded; published gold-standard for human cell lines | Requires matched Hi-C / Micro-C; cell-type-specific; default contact uses average across 5 ENCODE cell types if Hi-C not available |
| ENCODE-rE2G (ENCODE 2024) | ATAC + H3K27ac + (Hi-C optional) | Logistic regression trained on CRISPRi-FlowFISH ground truth; uses ABC features + sequence features + distance | ENCODE 4 standard; pre-trained models for many cell types | Pre-trained models only available for ENCODE cell types; retraining requires CRISPRi data |
| Cicero (Pliner 2018) | scATAC peak-cell matrix | Graphical lasso on metacell co-accessibility | ATAC-only; works without Hi-C | Less concordant with Hi-C than ABC; cis-distance-limited; alpha-sensitive |
| HiChIP H3K27ac + FitHiChIP | H3K27ac HiChIP | Statistically significant loops at FDR < 0.05 | Direct experimental loop measurement; cell-type-specific; orthogonal to ATAC | Requires HiChIP wet-lab; only captures loops within HiChIP resolution (~10 kb) |
| Hi-C + HiCCUPS | Bulk Hi-C | Fold-enrichment loop calling | Most-validated 3D contact method | Resolution typically 5-25 kb; misses sub-loop fine structure |
| Capture Hi-C / PCHi-C (CHiCAGO) | Promoter Capture Hi-C | Asymptotic CHiCAGO score | High-resolution promoter-anchored | Wet-lab cost; promoter capture only |
| EpiMap (Boix 2021) reference | None (pre-computed lookup) | Bulk-derived enhancer-gene predictions in 833 epigenomes | Fast, comprehensive | Cell-type-agnostic for tissues outside the reference set |
| GeneHancer / FANTOM5 (legacy) | None (pre-computed lookup) | Pre-computed; varied methods per database | Comprehensive lookup; widely cited | Older; less reliable than ABC for cell-type-specific |

Methodology evolves; verify against current Engreitz lab releases (ABC), ENCODE 4 publications (ENCODE-rE2G), and Mumbach 2017 (HiChIP) before locking pipelines.

## ABC Mathematics

For each candidate (enhancer E, gene G) pair within the cis window (default 5 Mb):

```
ABC(E -> G) = Activity_E * Contact_E,G / sum_{all e in window}(Activity_e * Contact_e,G)
```

- **Activity_E** = ATAC reads at E * H3K27ac reads at E (geometric mean of normalized signals; reflects "enhancer strength")
- **Contact_E,G** = Hi-C/Micro-C contact frequency from E to G's TSS (after distance-correction)
- **Window** = +/- 5 Mb cis (default; ENCODE-rE2G uses 1 Mb)

Threshold typical: ABC >= 0.02 for high-confidence; >= 0.01 for exploratory.

When Hi-C is unavailable, ABC uses an "average contact" derived from 5 ENCODE Hi-C cell types as proxy. Performance degrades but it's not zero (Fulco 2019).

## ENCODE-rE2G Differences from ABC

ENCODE-rE2G (2024) is a reformulation:

- **Logistic regression** trained on CRISPRi-FlowFISH ground truth (~10 cell types)
- **Features:** ABC score components + 3D contact + distance + activity ratios
- **Multiple feature configurations:** "abc-features", "no-hic-features" for cells without 3D data
- **Output:** Per-pair probability of regulatory connection
- **Pre-trained models** for ENCODE cell lines; logistic params vary by cell type

ENCODE-rE2G generally outperforms ABC at CRISPRi recall, especially at modest distances (50-500 kb). For ENCODE cell types, prefer ENCODE-rE2G; for novel cell types, ABC remains the default.

## Per-Tool Failure Modes

### ABC -- Wrong cell-type-matched Hi-C

**Trigger:** Using K562 Hi-C contact when actual cell type is GM12878.

**Mechanism:** Contact frequencies differ across cell types at compartment and TAD boundaries; using mismatched Hi-C produces wrong ABC scores.

**Symptom:** ABC predictions concentrate at known K562-specific loci even when ATAC data is from GM12878.

**Fix:** Use cell-type-matched Hi-C or Micro-C. If unavailable, ABC's "average HiC" (5-cell-type pooled) is the documented fallback with acknowledged degradation. Document the proxy in methods.

### ABC -- H3K27ac normalization

**Trigger:** H3K27ac ChIP-seq with different sequencing depth than ATAC.

**Mechanism:** ABC's "Activity" is the geometric mean of accessibility and H3K27ac signals; both must be normalized to the same scale.

**Symptom:** Activity scores skewed; some peaks have very high activity from H3K27ac alone, others from ATAC alone.

**Fix:** Normalize both signals to reads-per-million in peaks (RPM-IP) before combining. Use ABC's `--qnorm` flag with a quantile-normalization reference file (e.g. `--qnorm src/EnhancersQNormRef.K562.txt` from the ABC repo).

### ENCODE-rE2G -- Cell type not in pre-trained set

**Trigger:** Running pre-trained model on a primary cell type not in CRISPRi training.

**Mechanism:** Logistic regression coefficients learned from ENCODE cell types may not transfer to primary tissues.

**Fix:** Use the closest ENCODE cell type (myeloid lineage -> K562; lymphoid -> GM12878; hepatic -> HepG2). Document the proxy. For high-stakes use, custom retraining requires CRISPRi-FlowFISH data.

### Cicero -- No Hi-C concordance benchmark

**Trigger:** Reporting Cicero connections as enhancer-gene calls without external validation.

**Mechanism:** Cicero is statistical co-accessibility; correlation with Hi-C 3D contacts is ~30-50%. Many strong Cicero connections are NOT Hi-C-validated.

**Fix:** When Hi-C is available, cross-validate; report both. When only ATAC, use Cicero with the explicit caveat that connections are co-accessibility hypotheses, not contact predictions.

### HiChIP -- Loop calling threshold

**Trigger:** Default FitHiChIP at FDR < 0.05.

**Mechanism:** HiChIP loops are abundant (10k-100k per dataset); FDR alone produces a long tail of weak loops.

**Fix:** Threshold at FDR < 0.05 AND number of contacts per loop >= 5; or use the top N most significant where N = expected number of loops based on cell type.

### EpiMap / GeneHancer -- Cell-type-agnostic limitation

**Trigger:** Using EpiMap or GeneHancer pre-computed pairs for a specific cell type.

**Mechanism:** These references aggregate across many tissues / experiments; cell-type-specific connections are diluted.

**Fix:** Use as a baseline / sanity check, not as the primary call. ABC or ENCODE-rE2G in the actual cell type is preferred.

## Decision Tree by Available Data

| Available data | Recomm