Claude Cursor Skill

hic-aggregation

Build comprehensive chromatin contact maps by aggregating Hi-C loop calls (BEDPE) across multiple ENCODE experiments, donors, and labs. Use when the user wants to answer "what regions are in 3D contact in my tissue?" by creating a union catalog of chromatin loops. Handles resolut

LLM Mart · 0 points · 0 views 0 listing impressions 0 install-command copies
Virus-scanned Reviewed automatically before listing.

Full trust report

Download ammawla-encode-toolkit-plugin_skills_hic-aggregation-36836c8.zip · 18 KB
Part of ammawla/encode-toolkit — 90 skills

Install

skills CLI npx skills add https://github.com/ammawla/encode-toolkit/tree/main/plugin/skills/hic-aggregation
Claude Code claude plugin marketplace add https://llmmart.ai/marketplace.json && claude plugin install ammawla-encode-toolkit@llmmart
Git git clone https://github.com/ammawla/encode-toolkit.git

The skills CLI installs just this skill, for any of its supported agents. Claude Code installs the whole ammawla/encode-toolkit collection as a plugin from our marketplace. Git is the plain clone.

Skill manifest

Aggregate Hi-C Chromatin Contacts Across Studies

When to Use

  • User wants to build a comprehensive catalog of chromatin loops from multiple Hi-C experiments
  • User asks "what regions are in 3D contact in my tissue?" or "aggregate loop calls across donors"
  • User needs a union catalog of BEDPE loops with resolution-aware anchor matching
  • User wants to identify high-confidence loops supported by multiple experiments
  • Example queries: "aggregate Hi-C loops for K562", "combine chromatin contacts across labs", "find consensus TAD boundaries in liver"

Build a comprehensive catalog of chromatin loops for a tissue/cell type by merging BEDPE loop calls from multiple ENCODE Hi-C experiments.

Scientific Rationale

The question: "What regions are in 3D physical contact in my tissue?"

Like histone marks and accessibility, chromatin loops are a detection question. If a loop between Region A and Region B is detected in one donor but not another, the contact is still real — individual variation, sequencing depth, and computational resolution explain absence. We want the union of all detected contacts.

Key Concepts

Hi-C data measures pairwise chromatin interactions genome-wide. After processing:

  • Contact matrix (.hic file): Genome-wide interaction frequencies at multiple resolutions
  • Loop calls (BEDPE): Statistically significant point interactions (loops) identified by algorithms like HICCUPS or Juicer
  • TAD boundaries: Topologically associating domain boundaries
  • Compartments: A/B compartment assignments

BEDPE format (Paired-End BED):

chr1  start1  end1  chr2  start2  end2  name  score  strand1  strand2

Each row represents a contact between two genomic anchor regions.

Literature Support

  • Loop Catalog (Reyna et al. 2025, Nucleic Acids Research): Created a union catalog of 4.19M unique loops across 1,089 Hi-C datasets. Demonstrated that union approach captures tissue-specific and constitutive loops. Used resolution-aware merging at 5kb, 10kb, and 25kb bins.
  • AQuA Tools (Chakraborty et al. 2025): Toolkit for BEDPE intersection, union, and annotation. Handles paired-region arithmetic.
  • mariner (Flores et al. 2024, Bioinformatics): R/Bioconductor package for BEDPE manipulation including merging loops across experiments with configurable anchor tolerance.
  • ENCODE Phase 3 (Gorkin et al. 2020, Nature, 301 citations): Integrated Hi-C data across tissues to define regulatory loops connecting enhancers to promoters.
  • ENCODE Blacklist (Amemiya et al. 2019, Scientific Reports, 1,372 citations): Problematic genomic regions to filter from loop anchors. DOI
  • Mustache (Roayaei Ardakany et al. 2020, Genome Biology, 165 citations): Multi-scale loop caller that recovers more validated loops than HICCUPS. Different callers produce discordant loop sets.
  • Wolff et al. 2022 (GigaScience): Benchmark showing loop callers intersect by ~50% at most — critical context for why union approach is necessary.

Step 1: Find All Available Hi-C Data

encode_search_experiments(
    assay_title="Hi-C",
    organ="pancreas",           # user's tissue of interest
    biosample_type="tissue",
    limit=100
)

Present a summary to the user:

  • Total Hi-C experiments
  • Labs represented
  • Unique donors/biosamples
  • Resolution(s) available (check experiment metadata)

Use encode_get_facets to check availability:

encode_get_facets(assay_title="Hi-C", organ="pancreas")

Note: Hi-C data is computationally expensive to produce, so there are typically fewer experiments per tissue than ChIP-seq or ATAC-seq. Even 2-3 experiments can be valuable for union catalogs.

Step 2: Quality-Gate Each Experiment

encode_get_experiment(accession="ENCSR...")

Hi-C Quality Checks

  • Audit status: no ERROR flags
  • Sequencing depth: 400M+ valid read pairs for loop calling (ENCODE standard)
  • Cis/trans ratio: >60% cis contacts expected (low cis suggests noisy library)
  • Hi-C-specific QC: Library complexity, PCR duplicate rate
  • Has loop calls (BEDPE output) — not all Hi-C experiments have called loops
  • Resolution: at least 5-10kb resolution for loop detection

Include if:

  • Has BEDPE loop calls at consistent resolution
  • Passes ENCODE audit (no ERROR flags)
  • Adequate sequencing depth for loop resolution

Exclude if:

  • ERROR audit flags
  • Only contact matrices without loop calls
  • Very low sequencing depth (<200M valid pairs — insufficient for loop calling)

Track all included experiments:

encode_track_experiment(accession="ENCSR...")

Step 3: Download Loop Call Files

For each experiment, get BEDPE loop calls:

# Search for loop/interaction files
encode_list_files(
    experiment_accession="ENCSR...",
    file_format="bedpe",
    assembly="GRCh38"
)

# Or ask for loop calls by output type
encode_list_files(
    experiment_accession="ENCSR...",
    output_type="loops",
    assembly="GRCh38"
)

# Or contact domains
encode_list_files(
    experiment_accession="ENCSR...",
    output_type="contact domains",
    assembly="GRCh38"
)

File selection priority:

  1. Chromatin interactions (loop calls from HICCUPS or similar)
  2. Contact domains (TADs — different analysis, handle separately)
  3. Replicated loops (if available)

Prefer preferred_default=True files when available.

encode_download_files(
    file_accessions=["ENCFF...", ...],
    download_dir="/path/to/data/hic_loops",
    organize_by="flat"
)

Validate the downloaded BEDPE files before filtering. The report gives the detected anchor resolution, which Step 4 needs; gzipped inputs are read directly.

python3 scripts/validate_loops.py sample.bedpe [--min-distance 20000] [--expected-resolution 10000]

Step 4: Understanding Hi-C Resolution and Anchors

Critical: Resolution-Aware Processing

Hi-C loop anchors are binned regions, not precise positions. The resolution determines anchor size:

Resolution Anchor Width Best For Typical Loop Count
5 kb 5,000 bp Fine-scale promoter-enhancer loops More loops
10 kb 10,000 bp Standard analysis Moderate
25 kb 25,000 bp Large-scale domain contacts Fewer loops

All loops being merged must be at the same resolution, or anchors must be harmonized to a common resolution.

Harmonizing Resolution

If experiments have loops called at different resolutions:

# Expand 5kb anchors to 10kb resolution
awk -v res=10000 'BEGIN{OFS="\t"} {
    # Bin anchor 1
    bin1_start = int($2/res) * res
    bin1_end = bin1_start + res
    # Bin anchor 2
    bin2_start = int($5/res) * res
    bin2_end = bin2_start + res
    print $1, bin1_start, bin1_end, $4, bin2_start, bin2_end, $7, $8, $9, $10
}' fine_res_loops.bedpe > harmonized_loops.bedpe

Step 5: Per-Sample Filtering

5a. ENCODE Blocklist Filtering (Amemiya et al. 2019)

Remove loops with anchors in artifact-prone regions (download from https://github.com/Boyle-Lab/Blacklist/blob/master/lists/hg38-blacklist.v2.bed.gz):

# Filter loops where EITHER anchor overlaps a blocklist region
gunzip -k hg38-blacklist.v2.bed.gz

# First, extract anchor 1 and anchor 2 as separate BED files
awk 'BEGIN{OFS="\t"} {print $1,$2,$3,NR}' sample.bedpe > anchors1.bed
awk 'BEGIN{OFS="\t"} {print $4,$5,$6,NR}' sample.bedpe > anchors2.bed

# Find anchor rows NOT in blocklist
bedtools intersect -a anchors1.bed -b hg38-blacklist.v2.bed -v | cut -f4 > clean_rows_1.txt
bedtools intersect -a anchors2.bed -b hg38-blacklist.v2.bed -v | cut -f4 > clean_rows_2.txt

# Keep only rows where BOTH anchors pass
comm -12 <(sort clean_rows_1.txt) <(sort clean_rows_2.txt) > clean_rows.txt
awk 'NR==FNR{a[$1];next} FNR in a' clean_rows.txt sample.bedpe > sample.filtered.bedpe

5b. Score Filtering

Filter by interaction score/significance:

# If BEDPE has a score column (col 8), filter to significant interactions
# Keep top 75% by score (true distribution quantile, not range-based)
TOTAL=$(wc -l < sample.filtered.bedpe)
LINE_25=$(echo "$TOTAL" | awk '{printf "%d", $1 * 0.25}')
THRESHOLD=$(sort -k8,8n sample.filtered.bedpe | awk -v line="$LINE_25" 'NR==line{print $8}')
awk -v t="$THRESHOLD" '$8 >= t' sample.filtered.bedpe > sample.qfiltered.bedpe

5c. Remove Self-Ligation Artifacts

Loops where both anchors are very close are likely artifacts:

# Remove loops where anchors are on same chromosome and < 20kb apart
awk '{
    if ($1 != $4) print $0;  # inter-chromosomal: keep (rare but real)
    else if (($5 - $3) >= 20000) print $0;  # > 20kb apart: keep
}' sample.qfiltered.bedpe > sample.clean.bedpe

Step 6: Union Merge of Loops

The Paired-Region Matching Problem

Unlike peaks (single regions), loops are pairs of regions. Two loops match if both anchors overlap:

Loop 1:  [anchor1A]--------[anchor1B]
Loop 2:    [anchor2A]------[anchor2B]

These should merge if anchor1A overlaps anchor2A AND anchor1B overlaps anchor2B.

Method A: bedtools pairToPair (Recommended for simple union)

# Concatenate all filtered loops
cat sample1.clean.bedpe sample2.clean.bedpe ... > all_loops.bedpe

# Sort by anchor 1 coordinates
sort -k1,1 -k2,2n -k4,4 -k5,5n all_loops.bedpe > all_loops.sorted.bedpe

# Use a custom merge approach:
# 1. Bin anchors to resolution, creating a loop ID
# 2. Group by loop ID
# 3. Count support

awk -v res=10000 'BEGIN{OFS="\t"} {
    # Create binned anchor coordinates as loop identifier
    a1_bin = $1 ":" int($2/res)*res
    a2_bin = $4 ":" int($5/res)*res
    # Canonical order (smaller coordinate first) to handle orientation
    if (a1_bin < a2_bin) loop_id = a1_bin "-" a2_bin
    else loop_id = a2_bin "-" a1_bin
    print loop_id, $0
}' all_loops.sorted.bedpe | \
sort -k1,1 | \
awk 'BEGIN{OFS="\t"} {
    if ($1 != prev_id) {
        if (NR > 1) print chr1, start1, end1, chr2, start2, end2, count, max_score
        prev_id = $1
        chr1=$2; start1=$3; end1=$4; chr2=$5; start2=$6; end2=$7
        count = 1; max_score = $9
    } else {
        count++
        if ($9 > max_score) max_score = $9
        # Expand anchors to encompass all overlapping calls
        if ($3 < start1) start1 = $3
        if ($4 > end1) end1 = $4
        if ($6 < start2) start2 = $6
        if ($7 > end2) end2 = $7
    }
} END {
    print chr1, start1, end1, chr2, start2, end2, count, max_score
}' > union_loops.bedpe

Method B: Resolution-Binned Approach (Loop Catalog method)

Following the Loop Catalog (Reyna et al. 2025) approach:

# Bin all loop anchors to a fixed resolution
awk -v res=10000 'BEGIN{OFS="\t"} {
    a1_start = int($2/res) * res
    a1_end = a1_start + res
    a2_start = int($5/res) * res
    a2_end = a2_start + res
    # Canonical ordering
    if ($1 < $4 || ($1 == $4 && a1_start <= a2_start))
        print $1, a1_start, a1_end, $4, a2_start, a2_end
    else
        print $4, a2_start, a2_end, $1, a1_start, a1_end
}' all_loops.sorted.bedpe | \
sort -u | \
sort -k1,1 -k2,2n -k4,4 -k5,5n | \
uniq -c | \
awk 'BEGIN{OFS="\t"} {print $2,$3,$4,$5,$6,$7,$1}' > union_loops_binned.bedpe
# Columns: chr1, start1, end1, chr2, start2, end2, n_supporting_samples

Method C: Using Specialized Tools

mariner (R/Bioconductor):

library(mariner)
# Read BEDPE files as GInteractions
loops <- lapply(bedpe_files, read.table)
# Convert to GInteractions and merge
gi <- as_ginteractions(loops)
merged <- mergePairs(gi, radius = 10000)  # 10kb tolerance

AQuA Tools (Python):

# BEDPE union with anchor overlap tolerance
aqua bedpe-union -i sample1.bedpe sample2.bedpe -o union.bedpe --slop 5000

Step 7: Confidence Annotation

Given N total experiments:

Confidence Criteria Interpretation
High Detected in >=50% of samples Constitutive loop, present across individuals
Supported Detected in 2+ samples Likely real, some variation
Singleton Detected in 1 sample only May be individual-specific or depth-dependent
awk -v N=4 '{
    if ($7 >= N*0.5) conf="HIGH";
    else if ($7 >= 2) conf="SUPPORTED";
    else conf="SINGLETON";
    print $0"\t"conf"\t"$7"/"N
}' union_loops_binned.bedpe > union_loops.annotated.bedpe

Context for singletons: Hi-C loop detection is very sensitive to sequencing depth. Many singletons may simply be under-powered in other samples rather than biologically absent. The Loop Catalog found that a union approach captures ~3x more loops than any individual experiment.

Step 8: Separate Analysis for TADs and Compartments

TAD boundaries and A/B compartments require different aggregation than loops:

TAD Boundaries

TAD boundaries are single genomic positions. Aggregate like narrow peaks:

# Extract TAD boundary BED from contact domain files
# Each boundary is a narrow region
cat tad_boundaries_sample*.bed | \
bedtools sort -i - | \
bedtools merge -i - -d 40000 -c 1 -o count > union_tad_boundaries.bed
# 40kb gap tolerance because TAD boundaries are resolution-dependent

A/B Compartments

Compartment calls (eigenvector sign at each bin) should be aggregated by majority vote:

# For each resolution bin, assign A or B based on majority of samples
# This is more complex and typically done in R/Python

Step 9: Log Provenance

encode_log_derived_file(
    file_path="/path/to/union_loops.annotated.bedpe",
    source_accessions=["ENCSR...", "ENCSR...", ...],
    description="Union chromatin loops across N pancreas Hi-C experiments",
    file_type="aggregated_loops",
    tool_used="bedtools + custom merge at 10kb resolution",
    parameters="blocklist filtered, score >= 25th pctl, self-ligation >= 20kb removed, 10kb resolution binning"
)

Step 10: Summary Statistics

Report to the user:

  • Total input experiments: N
  • Experiments passing QC: M
  • Resolution used: Xkb
  • Total loops before merge: X
  • Union loops after merge: Y
  • High-confidence loops: Z (≥50% support)
  • Supported loops: W (2+ support)
  • Singleton loops: V (1 sample only)
  • Distance distribution: median and range of loop sizes (anchor-to-anchor)
  • Inter-chromosomal loops: count (expect very few)

Pitfalls Specific to Hi-C Data

  1. Resolution mismatch: Loop calls at 5kb vs 25kb resolution will have very different anchor sizes. Always harmonize to a common resolution before merging.

  2. Sequencing depth sensitivity: Loop calling requires deep sequencing (400M+ valid pairs). Shallowly sequenced experiments will call far fewer loops — this is under-detection, not absence.

  3. Algorithm differences are LARGE: Wolff et al. 2022 (GigaScience) found that HICCUPS, Mustache, Fit-Hi-C, and HiCExplorer loop callers intersect by ~50% at most. Mustache tends to recover more validated loops (Roayaei Ardakany et al. 2020). If mixing callers, note this in provenance — and this discordance is itself a reason to prefer the union approach.

  4. Orientation matters: BEDPE anchors should be canonically ordered (anchor1 < anchor2 by genomic coordinate) before merging to avoid duplicate counting.

  5. Inter-chromosomal contacts: These are rare but real. Handle separately — they cannot be distance-filtered.

  6. Distance distribution: Most loops are 100kb-2Mb. Very short-range contacts (<20kb) are often noise from undigested chromatin. Very long-range (>10Mb) are rare.

  7. Do NOT mix assemblies: All files must be GRCh38 or all hg19. Hi-C resolution binning makes liftOver of loops particularly error-prone.

  8. TADs vs loops: These are different features. TADs are domains (regions), loops are point contacts (pairs). Do not mix them in the same union.

  9. Micro-C as complement: Micro-C achieves higher resolution than Hi-C and can detect sub-TAD loops. Treat Micro-C loops as compatible with Hi-C loops in a union (Mustache works on both).

Walkthrough: Building a Cross-Tissue Loop Catalog for the MYC Locus

Goal: Aggregate Hi-C chromatin loops across tissues to identify conserved and tissue-specific 3D contacts at the MYC gene locus. Context: Cancer research — MYC is regulated by distal enhancers via chromatin looping.

Step 1: Find Hi-C experiments across tissues

encode_search_experiments(assay_title="Hi-C", organism="Homo sapiens", limit=50)

Expected output:

{
  "results": [
    {"accession": "ENCSR000AKA", "assay_title": "Hi-C", "biosample_summary": "GM12878", "status": "released"},
    {"accession": "ENCSR489OCU", "assay_title": "Hi-C", "biosample_summary": "K562", "status": "released"},
    {"accession": "ENCSR382RFU", "assay_title": "Hi-C", "biosample_summary": "liver", "status": "released"}
  ],
  "total": 89,
  "limit": 50,
  "offset": 0,
  "has_more": true,
  "next_offset": 50
}

Interpretation: 89 Hi-C experiments available. Select 5–10 spanning diverse tissue types for cross-tissue comparison.

Step 2: List loop files for each experiment

encode_list_files(experiment_accession="ENCSR000AKA", file_format="bedpe", assembly="GRCh38")

Expected output (a JSON array of file records; fields abridged):

[
  {"accession": "ENCFF001ABC", "output_type": "contact domains", "file_format": "bedpe", "file_size_human": "2.4 MB"},
  {"accession": "ENCFF002DEF", "output_type": "loops", "file_format": "bedpe", "file_size_human": "1.8 MB"}
]

Interpretation: Use "loops" files for loop aggregation. Contact domains are TADs, not loops.

Step 3: Download loop files

encode_download_files(file_accessions=["ENCFF002DEF", "ENCFF003GHI", "ENCFF004JKL"], download_dir="/data/hic_loops")

Expected output (one of the three downloaded entries shown):

{
  "downloaded": [
    {
      "accession": "ENCFF002DEF",
      "file_path": "/data/hic_loops/ENCFF002DEF.bedpe",
      "file_size": 1887436,
      "file_size_human": "1.8 MB",
      "success": true,
      "error": "",
      "md5_verified": true
    }
  ],
  "errors": [],
  "summary": {
    "total_requested": 3,
    "successful": 3,
    "failed": 0,
    "total_size": 5872025,
    "total_size_human": "5.6 MB"
  }
}

Step 4: Aggregate loops with resolution-aware anchor matching

Apply union merge across tissues:

  • Expand loop anchors by ±resolution (e.g., ±5kb for 5kb resolution data)
  • Merge overlapping anchors using bedtools pairToPair
  • Assign tissue support counts to each union loop
  • Filter: require ≥2 tissue support for conserved loops

Step 5: Filter to MYC locus

# MYC locus: chr8:127,700,000-128,000,000
awk '$1=="chr8" && $2>=127700000 && $3<=128000000' union_loops.bedpe > myc_loops.bedpe

Interpretation: Loops anchored at the MYC promoter connecting to distal enhancers. Conserved loops (≥3 tissues) likely represent fundamental regulatory architecture; tissue-specific loops may drive context-dependent MYC activation.

Integration with downstream skills

  • Feed loop anchors into → peak-annotation for gene assignment at anchor regions
  • Overlay with → histone-aggregation H3K27ac peaks to identify active enhancer-promoter loops
  • Cross-reference loop-disrupting variants via → variant-annotation
  • Visualize in → ucsc-browser as interaction tracks

Code Examples

1. Survey available Hi-C data by tissue

encode_get_facets(assay_title="Hi-C", organism="Homo sapiens")

Expected output (facet field names are the top-level keys):

{
  "biosample_ontology.organ_slims": [
    {"term": "brain", "count": 24},
    {"term": "blood", "count": 15}
  ]
}

2. Get details for a specific Hi-C experiment

encode_get_experiment(accession="ENCSR000AKA")

Expected output (fields abridged):

{
  "accession": "ENCSR000AKA",
  "assay_title": "Hi-C",
  "biosample_summary": "GM12878",
  "assembly": ["GRCh38"],
  "bio_replicate_count": 2,
  "status": "released",
  "lab": "Erez Lieberman Aiden, Baylor",
  "audit_error_count": 0,
  "audit_not_compliant_count": 0,
  "audit_warning_count": 1,
  "audit_internal_action_count": 0
}

3. Compare loop sets between two cell types

Both experiments must already be tracked; otherwise the tool returns {"error": "Experiment ... not tracked. Track it first."}.

encode_compare_experiments(accession1="ENCSR000AKA", accession2="ENCSR489OCU")

Expected output:

{
  "experiment_1": {"accession": "ENCSR000AKA", "assay": "Hi-C", "biosample": "GM12878"},
  "experiment_2": {"accession": "ENCSR489OCU", "assay": "Hi-C", "biosample": "K562"},
  "verdict": "COMPATIBLE_WITH_CAVEATS",
  "recommendation": "These experiments can be compared, but the warnings should be addressed in your analysis.",
  "compatible_aspects": [
    "Same organism: Homo sapiens",
    "Same assembly: GRCh38",
    "Same assay: Hi-C",
    "Same biosample type: cell line"
  ],
  "issues": [],
  "warnings": [
    "Different labs: Erez Lieberman Aiden, Baylor vs Job Dekker, UMass. Batch effects possible."
  ]
}

Integration

This skill produces... Feed into... Purpose
Union loop catalog (BEDPE) peak-annotation Assign genes to loop anchors
Conserved loop coordinates histone-aggregation Overlay H3K27ac at anchors to find active enhancer-promoter loops
Tissue-specific loops accessibility-aggregation Check if loop anchors overlap open chromatin
Loop anchor BED intervals variant-annotation Find GWAS/clinical variants disrupting loop anchors
Loop anchor coordinates liftover-coordinates Convert hg19 loops to GRCh38
Aggregated loop statistics visualization-workflow Generate loop frequency heatmaps
Loop-gene assignments disease-research Connect loop disruptions to disease phenotypes

Related Skills

  • histone-aggregation: Loop anchors often overlap with H3K27ac/H3K4me1 peaks — integrate with histone union sets to annotate loop function
  • accessibility-aggregation: Loop anchors frequently coincide with accessible chromatin — validate loops by requiring anchor accessibility
  • regulatory-elements: Use loops to connect distal enhancers (H3K27ac) to target promoters (H3K4me3)
  • epigenome-profiling: Loops add 3D context to 1D chromatin state maps
  • pipeline-hic: Process raw Hi-C data through the full ENCODE-aligned pipeline
  • batch-analysis: Batch processing workflows for systematic Hi-C loop aggregation
  • publication-trust: Verify literature claims backing analytical decisions

Presenting Results

  • Present aggregated loops as: chr | anchor1_start | anchor1_end | anchor2_start | anchor2_end | sample_count | resolution. Show loop statistics. Suggest: "Would you like to check if any GWAS variants overlap loop anchors?"

For the request: "$ARGUMENTS"

Files (encode-toolkit)
  • references
    • literature.md 7.2 KB
      # Hi-C Aggregation — Literature References
      
      **Last updated:** 2026-03-07
      **Purpose:** Reference catalog for the hic-aggregation skill — papers supporting the union-based approach for aggregating Hi-C chromatin loop calls across experiments, resolution-aware anchor matching, and the critical finding that loop callers produce highly discordant results.
      
      ---
      
      ## Hi-C Loop Catalogs
      
      ---
      
      ### Reyna et al. 2025 — Loop Catalog: comprehensive chromatin loop resource
      
      - **Citation:** Reyna J, et al. Loop Catalog: a comprehensive, cell type-specific resource of chromatin loops. Nucleic Acids Research, 2025.
      - **Citations:** ~20
      - **Key findings:** Created the largest union catalog of chromatin loops, containing 4.19 million unique loops across 1,089 Hi-C datasets from human cell types and tissues. Used resolution-aware merging at 5kb, 10kb, and 25kb bins, matching loop anchors within one bin-width tolerance to account for resolution-dependent positional uncertainty. Demonstrated that the union approach captures both constitutive loops (shared across >50% of cell types) and cell-type-specific loops, providing the most comprehensive 3D genome map available. This catalog validates the union-based aggregation strategy used in this skill.
      
      ---
      
      ### Wolff et al. 2022 — Loop caller benchmarking
      
      - **Citation:** Wolff J, Bhardwaj V, Nothjunge S, Richard G, Renschler G, Gilsbach R, Manke T, Backofen R, Ramírez F, Grüning BA. Galaxy HiCExplorer 3: a web server for reproducible Hi-C, capture Hi-C and single-cell Hi-C data analysis, quality control and visualization. Nucleic Acids Research, 50(W1):W697-W705, 2022.
      - **DOI:** [10.1093/nar/gkac407](https://doi.org/10.1093/nar/gkac407)
      - **PMID:** 35639517
      - **Citations:** ~100
      - **Key findings:** Benchmark showing that different Hi-C loop callers (HiCCUPS, HICCUPS2, Mustache, HiCExplorer, Fit-Hi-C2, chromosight) intersect by approximately 50% at most when applied to the same dataset. This ~50% concordance rate is the critical context for why union-based aggregation is necessary: no single caller captures all real loops, and combining calls from multiple callers and experiments provides the most complete contact map. This finding parallels the principle from histone ChIP-seq that absence of a peak does not mean absence of a binding site.
      
      ---
      
      ### Roayaei Ardakany et al. 2020 — Mustache: multi-scale loop detection
      
      - **Citation:** Roayaei Ardakany A, Gezer HT, Lonardi S, Ay F. Mustache: multi-scale detection of chromatin loops from Hi-C and Micro-C maps using scale-space representation. Genome Biology, 21:256, 2020.
      - **DOI:** [10.1186/s13059-020-02167-0](https://doi.org/10.1186/s13059-020-02167-0)
      - **PMID:** 33023656 | **PMC:** PMC7537270
      - **Citations:** ~165
      - **Key findings:** Introduced Mustache, a multi-scale loop caller that uses scale-space theory to detect loops at multiple resolutions simultaneously. Recovers more experimentally validated loops than HiCCUPS, particularly at lower sequencing depths, because multi-scale analysis captures both sharp and diffuse contact enrichments. Relevant to aggregation because combining Mustache and HiCCUPS calls from the same dataset substantially increases the loop catalog — further supporting the union approach.
      
      ---
      
      ## Hi-C Data Processing
      
      ---
      
      ### Rao et al. 2014 — 3D map of the human genome
      
      - **Citation:** Rao SSP, Huntley MH, Durand NC, Stamenova EK, et al. A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping. Cell, 159(7):1665-1680, 2014.
      - **DOI:** [10.1016/j.cell.2014.11.021](https://doi.org/10.1016/j.cell.2014.11.021)
      - **PMID:** 25497547 | **PMC:** PMC5635824
      - **Citations:** ~5,000
      - **Key findings:** Introduced HiCCUPS (Hi-C Computational Unbiased Peak Search) for systematic loop detection, the primary loop caller used in ENCODE Hi-C experiments. Identified ~10,000 loops in GM12878 cells with convergent CTCF motifs at anchors. Loop calls from HiCCUPS are the primary input for aggregation in this skill, stored in BEDPE format.
      
      ---
      
      ### Durand et al. 2016 — Juicer tools
      
      - **Citation:** Durand NC, Shamim MS, Machol I, et al. Juicer provides a one-click system for analyzing loop-resolution Hi-C experiments. Cell Systems, 3(1):95-98, 2016.
      - **DOI:** [10.1016/j.cels.2016.07.002](https://doi.org/10.1016/j.cels.2016.07.002)
      - **PMID:** 27467249 | **PMC:** PMC5846465
      - **Citations:** ~2,000
      - **Key findings:** Juicer pipeline and tools including HiCCUPS for loop calling and Arrowhead for TAD annotation. ENCODE Hi-C experiments are processed through Juicer, and HiCCUPS loop calls in BEDPE format are the primary input for aggregation. The .hic file format stores contact matrices at multiple resolutions enabling resolution-aware analysis.
      
      ---
      
      ## BEDPE Manipulation Tools
      
      ---
      
      ### Chakraborty et al. 2025 — AQuA Tools
      
      - **Citation:** Chakraborty A, et al. AQuA Tools: toolkit for paired-region analysis of Hi-C and related data. 2025.
      - **Key findings:** Computational toolkit for BEDPE intersection, union, and annotation, specifically designed for Hi-C loop analysis. Handles the paired-region arithmetic required for loop aggregation: two loops are considered "the same" only if both anchors overlap (not just one). Provides functions for resolution-aware anchor matching where anchor overlap is defined within one bin-width tolerance (e.g., anchors within 10kb of each other at 10kb resolution are considered matching).
      
      ---
      
      ### Flores et al. 2024 — mariner: R/Bioconductor for BEDPE operations
      
      - **Citation:** Flores EK, et al. mariner: an R/Bioconductor package for exploring Hi-C data. Bioinformatics, 2024.
      - **Key findings:** R/Bioconductor package for BEDPE manipulation including merging loops across experiments with configurable anchor tolerance. Provides functions for loop intersection, union, and annotation within the Bioconductor framework, integrating with InteractionSet and GenomicRanges. Useful for users who prefer R-based analysis of aggregated Hi-C loop catalogs.
      
      ---
      
      ## Quality Framework
      
      ---
      
      ### Yardimci et al. 2019 — Hi-C data quality and reproducibility
      
      - **Citation:** Yardimci GG, Ozadam H, Sauria MEG, et al. Measuring the reproducibility and quality of Hi-C data. Genome Biology, 20:57, 2019.
      - **DOI:** [10.1186/s13059-019-1658-7](https://doi.org/10.1186/s13059-019-1658-7)
      - **PMID:** 30890173 | **PMC:** PMC6425651
      - **Citations:** ~250
      - **Key findings:** Established that cis/trans ratio and long-range cis fraction are the most informative quality metrics for Hi-C data. Experiments with poor quality metrics should be excluded from loop aggregation. Loop calling requires substantially deeper sequencing than compartment or TAD analysis — minimum 500M+ valid contacts for loops at 5-10kb resolution. These quality thresholds gate which experiments are included in the aggregation workflow.
      
      ---
      
      ### Amemiya et al. 2019 — ENCODE Blacklist
      
      - **DOI:** [10.1038/s41598-019-45839-z](https://doi.org/10.1038/s41598-019-45839-z) | **PMID:** 31249361 | **Citations:** ~1,372
      - **Hi-C aggregation role:** Loop anchors overlapping blacklisted regions should be filtered before aggregation because these regions produce artifactual high-contact signals due to multi-mapping reads, which can create false loops that would appear as high-confidence contacts in the aggregated catalog.
      
    • loop-caller-comparison.md 6.7 KB
      # Hi-C Loop Caller Comparison
      
      Reference guide covering concordance between Hi-C loop callers, based on Wolff et al. 2022 benchmarking and individual tool characteristics.
      
      ## The Concordance Problem
      
      Wolff et al. (2022, GigaScience) performed the most comprehensive benchmark of Hi-C loop callers to date. The critical finding:
      
      > Loop callers **intersect by approximately 50% at most**. Different algorithms applied to the same Hi-C data produce substantially different loop sets, even when tuned to similar sensitivity levels.
      
      This discordance is the primary motivation for the union approach in Hi-C aggregation. Any single caller misses loops that another caller recovers.
      
      ## Major Loop Callers
      
      ### HiCCUPS (Rao et al. 2014)
      
      - **Developer**: Aiden Lab (Baylor College of Medicine)
      - **Method**: Identifies enriched pixels in the contact matrix relative to local background using a donut-shaped kernel
      - **Resolution**: Works best at 5kb and 10kb
      - **Strengths**: Gold standard for point-source loop detection; GPU-accelerated
      - **Weaknesses**: Conservative (high specificity, lower sensitivity); requires deep sequencing (>1B contacts for 5kb resolution); misses weaker loops
      - **ENCODE usage**: Primary loop caller in ENCODE Hi-C pipeline
      
      ### Mustache (Roayaei Ardakany et al. 2020, 165 citations)
      
      - **Developer**: Ay Lab (La Jolla Institute)
      - **Method**: Scale-space theory with blob detection across multiple resolutions
      - **Resolution**: Multi-scale (5kb to 100kb simultaneously)
      - **Strengths**: Recovers more validated loops than HiCCUPS; works at multiple resolutions simultaneously; no GPU required; works with both Hi-C and Micro-C
      - **Weaknesses**: More liberal calling (higher sensitivity, more potential false positives)
      - **Key advantage**: Multi-scale approach avoids resolution-specific artifacts
      
      ### Fit-Hi-C / FitHiC2 (Ay et al. 2014; Kaul et al. 2020)
      
      - **Developer**: Ay Lab
      - **Method**: Statistical model of contact frequency as function of genomic distance; identifies outlier contacts
      - **Resolution**: Flexible (1kb to 1Mb)
      - **Strengths**: Principled statistical framework; handles distance decay properly; provides p-values and q-values
      - **Weaknesses**: Calls many short-range interactions that may not be structural loops; less specific for point interactions
      
      ### HiCExplorer (Wolff et al. 2020)
      
      - **Developer**: Manke Lab (Max Planck Institute)
      - **Method**: Z-score based enrichment relative to expected contact frequency
      - **Resolution**: Configurable, typically 10-50kb
      - **Strengths**: Full Hi-C analysis suite; integrates with other HiCExplorer tools
      - **Weaknesses**: Z-score approach sensitive to normalization method
      
      ### Chromosight (Matthey-Doret et al. 2020)
      
      - **Developer**: Koszul Lab (Institut Pasteur)
      - **Method**: Template matching using convolution with loop kernel
      - **Resolution**: Flexible
      - **Strengths**: Fast; pattern-based approach detects loops, borders, and stripes
      - **Weaknesses**: Less commonly used in ENCODE context
      
      ## Wolff et al. 2022 Benchmark Results
      
      ### Pairwise Overlap Between Callers
      
      | Caller A | Caller B | Overlap (A in B) | Overlap (B in A) |
      |----------|----------|-----------------|-----------------|
      | HiCCUPS | Mustache | ~45% | ~40% |
      | HiCCUPS | FitHiC2 | ~50% | ~30% |
      | HiCCUPS | HiCExplorer | ~40% | ~35% |
      | Mustache | FitHiC2 | ~45% | ~35% |
      
      Note: Exact percentages vary by dataset and parameters. The consistent finding is that **no two callers agree on more than half their loops**.
      
      ### Factors Affecting Concordance
      
      1. **Sequencing depth**: Deeper datasets show better concordance (more power to detect weak loops)
      2. **Resolution**: 10kb resolution shows better concordance than 5kb (larger bins reduce stochastic variation)
      3. **Cell type**: Cell types with strong compartmentalization (e.g., GM12878) show better concordance than those with weaker structure
      4. **Normalization**: KR vs ICE vs VC normalization affects each caller differently
      
      ## Implications for Union Aggregation
      
      ### Why Union Is Appropriate
      
      1. **Different callers have different strengths**: HiCCUPS excels at strong point interactions; Mustache finds weaker but validated loops; FitHiC2 provides robust statistics
      2. **False negatives are the bigger problem**: For a catalog of "what loops exist," missing real loops is worse than including a few false positives
      3. **Loop Catalog precedent**: Reyna et al. (2025) created a union catalog of 4.19M loops across 1,089 datasets using this rationale
      
      ### Handling Mixed Caller Output
      
      When merging loops from different experiments that may have used different callers:
      
      ```bash
      # All loops go into the same union, regardless of which caller produced them
      # The sample/experiment ID tracks provenance, not the caller
      # Resolution harmonization is the critical step, not caller harmonization
      ```
      
      ### Confidence Annotation Accounts for Caller Variation
      
      The confidence system (HIGH / SUPPORTED / SINGLETON) naturally handles caller discordance:
      - A loop detected by multiple callers in the same experiment counts once per experiment
      - A loop detected across multiple experiments (regardless of caller) gets higher confidence
      - Caller-specific loops that do not replicate across experiments remain singletons
      
      ## Resolution Effects on Loop Detection
      
      | Resolution | Loops Detected | Concordance | Best For |
      |-----------|---------------|-------------|---------|
      | 1 kb | Most (finest scale) | Lowest | Micro-C data, promoter-enhancer |
      | 5 kb | Many | Moderate | Standard analysis, ENCODE |
      | 10 kb | Moderate | Highest | Cross-study comparison |
      | 25 kb | Fewest | High | Large-scale domain contacts |
      
      **Recommendation**: When merging loops across studies, harmonize to 10kb resolution for maximum concordance, unless all datasets support 5kb resolution.
      
      ## Practical Recommendations
      
      1. **Accept that callers disagree** -- this is a feature of loop detection, not a bug
      2. **Use union approach** for catalogs ("where are loops?")
      3. **Use intersection** only for high-confidence subsets needed for specific analyses
      4. **Record the caller** in provenance metadata for each experiment
      5. **Prefer Mustache** if re-calling loops from contact matrices (more validated loops)
      6. **Prefer HiCCUPS** if using ENCODE pre-called loops (already the pipeline standard)
      
      ## References
      
      - Wolff et al. 2022, GigaScience -- comprehensive loop caller benchmark (~50% concordance finding)
      - Roayaei Ardakany et al. 2020, Genome Biology -- Mustache multi-scale loop caller (165 citations)
      - Rao et al. 2014, Cell -- HiCCUPS algorithm and original loop catalog (4,900+ citations)
      - Kaul et al. 2020, Nature Methods -- FitHiC2 with improved statistical model
      - Reyna et al. 2025, Nucleic Acids Research -- Loop Catalog with 4.19M union loops
      - Matthey-Doret et al. 2020, Nature Communications -- Chromosight pattern-based detection
      
  • scripts
    • validate_loops.py 11.7 KB
      #!/usr/bin/env python3
      """Validate BEDPE loop files from ENCODE Hi-C aggregation.
      
      Checks BEDPE format, anchor validity, resolution consistency, canonical ordering,
      cis/trans classification, and reports summary statistics.
      
      Usage:
          python validate_loops.py input.bedpe [--min-distance 20000] [--expected-resolution 10000]
          python validate_loops.py input.bedpe --expected-resolution 5000
      
      Plain and gzipped (.gz) inputs are both accepted.
      """
      
      import argparse
      import gzip
      import statistics
      import sys
      from collections import Counter
      from pathlib import Path
      
      VALID_CHROMS = {f"chr{i}" for i in range(1, 23)} | {"chrX", "chrY", "chrM"}
      MIN_BEDPE_COLS = 6
      
      MAX_COLUMN_ERRORS = 5
      
      
      def parse_args():
          parser = argparse.ArgumentParser(
              description="Validate BEDPE loop files from ENCODE Hi-C aggregation.",
              formatter_class=argparse.RawDescriptionHelpFormatter,
              epilog=(
                  "Examples:\n"
                  "  python validate_loops.py loops.bedpe\n"
                  "  python validate_loops.py loops.bedpe --expected-resolution 10000\n"
                  "  python validate_loops.py loops.bedpe --min-distance 20000\n"
              ),
          )
          parser.add_argument("input", type=Path, help="Input BEDPE loop file")
          parser.add_argument(
              "--min-distance",
              type=int,
              default=20000,
              help="Minimum anchor-to-anchor distance for cis loops (bp). Default: 20000",
          )
          parser.add_argument(
              "--expected-resolution",
              type=int,
              default=None,
              help="Expected resolution in bp (e.g., 5000, 10000, 25000). Checks anchor size consistency.",
          )
          return parser.parse_args()
      
      
      def open_text(path):
          """Open a plain or gzipped text file for reading."""
          if str(path).endswith(".gz"):
              return gzip.open(path, "rt")
          return open(path)
      
      
      def fmt_distance(bp):
          """Render a genomic distance, in bp below 1kb so sub-kb values do not read as 0kb."""
          return f"{bp}bp" if bp < 1000 else f"{bp // 1000}kb"
      
      
      def quartiles(values):
          """25th percentile, median and 75th percentile (interpolated)."""
          median = statistics.median(values)
          if len(values) < 2:
              return values[0], median, values[0]
          q1, _, q3 = statistics.quantiles(values, n=4, method="inclusive")
          return q1, median, q3
      
      
      def validate_loops(input_path, min_distance, expected_resolution):
          errors = []
          warnings = []
          chrom_counts = Counter()
          anchor1_sizes = []
          anchor2_sizes = []
          loop_distances = []
          total_lines = 0
          # Every line excluded from the statistics below counts as malformed, so
          # total_lines == valid_loops + bad_lines always holds.
          bad_lines = 0
          column_errors = 0
          cis_loops = 0
          trans_loops = 0
          short_range = 0
          non_canonical = 0
      
          if not input_path.exists():
              print(f"ERROR: Input file not found: {input_path}", file=sys.stderr)
              sys.exit(1)
      
          with open_text(input_path) as f:
              for line_num, line in enumerate(f, 1):
                  if line.startswith("#") or line.startswith("track") or line.startswith("browser"):
                      continue
                  line = line.strip()
                  if not line:
                      continue
      
                  total_lines += 1
                  fields = line.split("\t")
      
                  # BEDPE requires at least 6 columns
                  if len(fields) < MIN_BEDPE_COLS:
                      bad_lines += 1
                      column_errors += 1
                      if column_errors <= MAX_COLUMN_ERRORS:
                          errors.append(f"Line {line_num}: expected >= {MIN_BEDPE_COLS} columns (BEDPE), got {len(fields)}")
                      elif column_errors == MAX_COLUMN_ERRORS + 1:
                          errors.append("... suppressing further column-count errors")
                      continue
      
                  chr1 = fields[0]
                  chr2 = fields[3]
      
                  # Chromosome validation for both anchors: a row with either anchor on an
                  # unrecognized chromosome is one malformed line, and is skipped entirely.
                  invalid_anchors = [
                      (label, value)
                      for label, value in (("anchor1", chr1), ("anchor2", chr2))
                      if value not in VALID_CHROMS and not value.startswith("chr")
                  ]
                  if invalid_anchors:
                      for label, value in invalid_anchors:
                          errors.append(f"Line {line_num}: invalid {label} chromosome '{value}'")
                      bad_lines += 1
                      continue
      
                  # Coordinate validation
                  try:
                      start1 = int(fields[1])
                      end1 = int(fields[2])
                      start2 = int(fields[4])
                      end2 = int(fields[5])
                  except ValueError:
                      errors.append(f"Line {line_num}: non-integer coordinates in anchor fields")
                      bad_lines += 1
                      continue
      
                  # Non-negative coordinates, and start < end for each anchor. A row that breaks either
                  # rule is malformed: count it once and keep it out of the statistics.
                  impossible_anchor = False
                  for label, val in [("start1", start1), ("end1", end1), ("start2", start2), ("end2", end2)]:
                      if val < 0:
                          errors.append(f"Line {line_num}: negative {label} ({val})")
                          impossible_anchor = True
                  if start1 >= end1:
                      errors.append(f"Line {line_num}: anchor1 start ({start1}) >= end ({end1})")
                      impossible_anchor = True
                  if start2 >= end2:
                      errors.append(f"Line {line_num}: anchor2 start ({start2}) >= end ({end2})")
                      impossible_anchor = True
                  if impossible_anchor:
                      bad_lines += 1
                      continue
      
                  a1_size = end1 - start1
                  a2_size = end2 - start2
                  anchor1_sizes.append(a1_size)
                  anchor2_sizes.append(a2_size)
      
                  # Cis vs trans classification
                  if chr1 == chr2:
                      cis_loops += 1
                      chrom_counts[chr1] += 1
      
                      # Distance check (midpoint-to-midpoint)
                      mid1 = (start1 + end1) // 2
                      mid2 = (start2 + end2) // 2
                      distance = abs(mid2 - mid1)
                      loop_distances.append(distance)
      
                      if distance < min_distance:
                          short_range += 1
      
                      # Canonical ordering: anchor1.start < anchor2.start
                      if start1 > start2:
                          non_canonical += 1
                  else:
                      trans_loops += 1
      
                  # Expected resolution check
                  if expected_resolution is not None:
                      if a1_size != expected_resolution:
                          if len(warnings) < 10:
                              warnings.append(f"Line {line_num}: anchor1 size {a1_size} != expected {expected_resolution}")
                      if a2_size != expected_resolution:
                          if len(warnings) < 10:
                              warnings.append(f"Line {line_num}: anchor2 size {a2_size} != expected {expected_resolution}")
      
          if total_lines == 0:
              print(f"ERROR: no data rows in {input_path} (only comments, headers or blank lines)", file=sys.stderr)
              sys.exit(1)
      
          valid_loops = total_lines - bad_lines
      
          # --- Detect resolution from anchor sizes ---
          detected_resolution = None
          if anchor1_sizes:
              all_sizes = anchor1_sizes + anchor2_sizes
              size_counts = Counter(all_sizes)
              most_common_size, most_common_count = size_counts.most_common(1)[0]
              size_pct = 100 * most_common_count / len(all_sizes)
              if size_pct > 50:
                  detected_resolution = most_common_size
      
          # --- Report Statistics ---
          print("=== Hi-C BEDPE Loop Validation Report ===")
          print(f"File: {input_path}")
          print()
      
          print("--- Summary ---")
          print(f"Data lines: {total_lines:,}")
          print(f"Valid loops: {valid_loops:,}")
          print(f"Malformed lines: {bad_lines}")
          print(f"Cis loops (same chromosome): {cis_loops:,} ({100 * cis_loops / max(valid_loops, 1):.1f}%)")
          print(f"Trans loops (inter-chromosomal): {trans_loops:,} ({100 * trans_loops / max(valid_loops, 1):.1f}%)")
          print(f"Non-canonical ordering (anchor1 > anchor2): {non_canonical:,}")
          print(f"Short-range loops (<{fmt_distance(min_distance)}): {short_range:,}")
          print()
      
          if detected_resolution:
              print("--- Resolution ---")
              print(f"Detected resolution: {detected_resolution:,} bp ({fmt_distance(detected_resolution)})")
              print(f"Anchor size consistency: {size_pct:.1f}% of anchors match detected resolution")
          elif anchor1_sizes:
              print("--- Resolution ---")
              print("WARNING: No consistent resolution detected. Anchor sizes vary.")
              top3 = size_counts.most_common(3)
              for sz, cnt in top3:
                  print(f"  {sz:>8,} bp: {cnt:,} anchors ({100 * cnt / len(all_sizes):.1f}%)")
          print()
      
          if loop_distances:
              sorted_dist = sorted(loop_distances)
              n = len(sorted_dist)
              dist_q1, dist_median, dist_q3 = quartiles(sorted_dist)
              print("--- Loop Distance Distribution (cis only) ---")
              print(f"Min:    {sorted_dist[0]:>12,} bp ({fmt_distance(sorted_dist[0])})")
              print(f"25th:   {dist_q1:>12,.1f} bp ({fmt_distance(int(dist_q1))})")
              print(f"Median: {dist_median:>12,.1f} bp ({fmt_distance(int(dist_median))})")
              print(f"75th:   {dist_q3:>12,.1f} bp ({fmt_distance(int(dist_q3))})")
              print(f"Max:    {sorted_dist[-1]:>12,} bp ({fmt_distance(sorted_dist[-1])})")
      
              # Distance buckets
              buckets = [
                  ("< 100kb", 0, 100_000),
                  ("100kb - 500kb", 100_000, 500_000),
                  ("500kb - 1Mb", 500_000, 1_000_000),
                  ("1Mb - 5Mb", 1_000_000, 5_000_000),
                  (">= 5Mb", 5_000_000, float("inf")),
              ]
              print("\n  Distance buckets:")
              for label, lo, hi in buckets:
                  count = sum(1 for d in loop_distances if lo <= d < hi)
                  print(f"    {label:<15} {count:>8,}  ({100 * count / n:.1f}%)")
              print()
      
          if cis_loops > 0:
              print("--- Chromosome Distribution (cis loops) ---")
              for chrom in sorted(chrom_counts.keys(), key=lambda c: (len(c), c)):
                  count = chrom_counts[chrom]
                  pct = 100 * count / max(cis_loops, 1)
                  print(f"  {chrom:<6} {count:>8,}  ({pct:5.1f}%)")
              print()
      
          # --- Warnings ---
          if trans_loops > valid_loops * 0.05:
              msg = (
                  f"WARNING: {100 * trans_loops / max(valid_loops, 1):.1f}% of loops are trans "
                  f"(inter-chromosomal). Expect <5% for typical Hi-C data."
              )
              print(msg, file=sys.stderr)
      
          if short_range > 0:
              msg = (
                  f"WARNING: {short_range:,} loops have anchors <{fmt_distance(min_distance)} apart. "
                  f"These may be self-ligation artifacts."
              )
              print(msg, file=sys.stderr)
      
          if non_canonical > 0:
              msg = (
                  f"WARNING: {non_canonical:,} loops have non-canonical ordering "
                  f"(anchor1.start > anchor2.start). Canonicalize before merging."
              )
              print(msg, file=sys.stderr)
      
          for w in warnings[:20]:
              print(w, file=sys.stderr)
      
          # --- Errors ---
          if errors:
              print(f"\n--- Errors ({len(errors)}) ---", file=sys.stderr)
              for e in errors[:50]:
                  print(f"  {e}", file=sys.stderr)
              if len(errors) > 50:
                  print(f"  ... and {len(errors) - 50} more errors", file=sys.stderr)
      
          has_errors = len(errors) > 0
          if has_errors:
              print(f"\nRESULT: FAIL -- {len(errors)} error(s) found", file=sys.stderr)
          else:
              print("\nRESULT: PASS -- file is valid BEDPE")
      
          return 1 if has_errors else 0
      
      
      if __name__ == "__main__":
          args = parse_args()
          exit_code = validate_loops(args.input, args.min_distance, args.expected_resolution)
          sys.exit(exit_code)
      
  • SKILL.md 22.8 KB
    ---
    name: hic-aggregation
    description: Build comprehensive chromatin contact maps by aggregating Hi-C loop calls (BEDPE) across multiple ENCODE experiments, donors, and labs. Use when the user wants to answer "what regions are in 3D contact in my tissue?" by creating a union catalog of chromatin loops. Handles resolution-aware anchor matching, cross-lab variation, and Hi-C-specific quality metrics.
    ---
    
    # Aggregate Hi-C Chromatin Contacts Across Studies
    
    ## When to Use
    
    - User wants to build a comprehensive catalog of chromatin loops from multiple Hi-C experiments
    - User asks "what regions are in 3D contact in my tissue?" or "aggregate loop calls across donors"
    - User needs a union catalog of BEDPE loops with resolution-aware anchor matching
    - User wants to identify high-confidence loops supported by multiple experiments
    - Example queries: "aggregate Hi-C loops for K562", "combine chromatin contacts across labs", "find consensus TAD boundaries in liver"
    
    Build a comprehensive catalog of chromatin loops for a tissue/cell type by merging BEDPE loop calls from multiple ENCODE Hi-C experiments.
    
    ## Scientific Rationale
    
    **The question**: "What regions are in 3D physical contact in my tissue?"
    
    Like histone marks and accessibility, chromatin loops are a **detection question**. If a loop between Region A and Region B is detected in one donor but not another, the contact is still real — individual variation, sequencing depth, and computational resolution explain absence. We want the **union of all detected contacts**.
    
    ### Key Concepts
    
    **Hi-C data** measures pairwise chromatin interactions genome-wide. After processing:
    - **Contact matrix** (`.hic` file): Genome-wide interaction frequencies at multiple resolutions
    - **Loop calls** (BEDPE): Statistically significant point interactions (loops) identified by algorithms like HICCUPS or Juicer
    - **TAD boundaries**: Topologically associating domain boundaries
    - **Compartments**: A/B compartment assignments
    
    **BEDPE format** (Paired-End BED):
    ```
    chr1  start1  end1  chr2  start2  end2  name  score  strand1  strand2
    ```
    Each row represents a contact between two genomic anchor regions.
    
    ### Literature Support
    - **Loop Catalog** (Reyna et al. 2025, Nucleic Acids Research): Created a union catalog of 4.19M unique loops across 1,089 Hi-C datasets. Demonstrated that union approach captures tissue-specific and constitutive loops. Used resolution-aware merging at 5kb, 10kb, and 25kb bins.
    - **AQuA Tools** (Chakraborty et al. 2025): Toolkit for BEDPE intersection, union, and annotation. Handles paired-region arithmetic.
    - **mariner** (Flores et al. 2024, Bioinformatics): R/Bioconductor package for BEDPE manipulation including merging loops across experiments with configurable anchor tolerance.
    - **ENCODE Phase 3** (Gorkin et al. 2020, Nature, 301 citations): Integrated Hi-C data across tissues to define regulatory loops connecting enhancers to promoters.
    - **ENCODE Blacklist** (Amemiya et al. 2019, Scientific Reports, 1,372 citations): Problematic genomic regions to filter from loop anchors. [DOI](https://doi.org/10.1038/s41598-019-45839-z)
    - **Mustache** (Roayaei Ardakany et al. 2020, Genome Biology, 165 citations): Multi-scale loop caller that recovers more validated loops than HICCUPS. Different callers produce discordant loop sets.
    - **Wolff et al. 2022** (GigaScience): Benchmark showing loop callers intersect by **~50% at most** — critical context for why union approach is necessary.
    
    ## Step 1: Find All Available Hi-C Data
    
    ```
    encode_search_experiments(
        assay_title="Hi-C",
        organ="pancreas",           # user's tissue of interest
        biosample_type="tissue",
        limit=100
    )
    ```
    
    Present a summary to the user:
    - Total Hi-C experiments
    - Labs represented
    - Unique donors/biosamples
    - Resolution(s) available (check experiment metadata)
    
    Use `encode_get_facets` to check availability:
    ```
    encode_get_facets(assay_title="Hi-C", organ="pancreas")
    ```
    
    **Note**: Hi-C data is computationally expensive to produce, so there are typically fewer experiments per tissue than ChIP-seq or ATAC-seq. Even 2-3 experiments can be valuable for union catalogs.
    
    ## Step 2: Quality-Gate Each Experiment
    
    ```
    encode_get_experiment(accession="ENCSR...")
    ```
    
    ### Hi-C Quality Checks
    - Audit status: no ERROR flags
    - **Sequencing depth**: 400M+ valid read pairs for loop calling (ENCODE standard)
    - **Cis/trans ratio**: >60% cis contacts expected (low cis suggests noisy library)
    - **Hi-C-specific QC**: Library complexity, PCR duplicate rate
    - Has loop calls (BEDPE output) — not all Hi-C experiments have called loops
    - Resolution: at least 5-10kb resolution for loop detection
    
    ### Include if:
    - Has BEDPE loop calls at consistent resolution
    - Passes ENCODE audit (no ERROR flags)
    - Adequate sequencing depth for loop resolution
    
    ### Exclude if:
    - ERROR audit flags
    - Only contact matrices without loop calls
    - Very low sequencing depth (<200M valid pairs — insufficient for loop calling)
    
    Track all included experiments:
    ```
    encode_track_experiment(accession="ENCSR...")
    ```
    
    ## Step 3: Download Loop Call Files
    
    For each experiment, get BEDPE loop calls:
    
    ```
    # Search for loop/interaction files
    encode_list_files(
        experiment_accession="ENCSR...",
        file_format="bedpe",
        assembly="GRCh38"
    )
    
    # Or ask for loop calls by output type
    encode_list_files(
        experiment_accession="ENCSR...",
        output_type="loops",
        assembly="GRCh38"
    )
    
    # Or contact domains
    encode_list_files(
        experiment_accession="ENCSR...",
        output_type="contact domains",
        assembly="GRCh38"
    )
    ```
    
    **File selection priority:**
    1. **Chromatin interactions** (loop calls from HICCUPS or similar)
    2. **Contact domains** (TADs — different analysis, handle separately)
    3. **Replicated loops** (if available)
    
    Prefer `preferred_default=True` files when available.
    
    ```
    encode_download_files(
        file_accessions=["ENCFF...", ...],
        download_dir="/path/to/data/hic_loops",
        organize_by="flat"
    )
    ```
    
    Validate the downloaded BEDPE files before filtering. The report gives the detected
    anchor resolution, which Step 4 needs; gzipped inputs are read directly.
    
    ```bash
    python3 scripts/validate_loops.py sample.bedpe [--min-distance 20000] [--expected-resolution 10000]
    ```
    
    ## Step 4: Understanding Hi-C Resolution and Anchors
    
    ### Critical: Resolution-Aware Processing
    
    Hi-C loop anchors are **binned regions**, not precise positions. The resolution determines anchor size:
    
    | Resolution | Anchor Width | Best For | Typical Loop Count |
    |-----------|-------------|---------|-------------------|
    | 5 kb | 5,000 bp | Fine-scale promoter-enhancer loops | More loops |
    | 10 kb | 10,000 bp | Standard analysis | Moderate |
    | 25 kb | 25,000 bp | Large-scale domain contacts | Fewer loops |
    
    **All loops being merged must be at the same resolution**, or anchors must be harmonized to a common resolution.
    
    ### Harmonizing Resolution
    If experiments have loops called at different resolutions:
    ```bash
    # Expand 5kb anchors to 10kb resolution
    awk -v res=10000 'BEGIN{OFS="\t"} {
        # Bin anchor 1
        bin1_start = int($2/res) * res
        bin1_end = bin1_start + res
        # Bin anchor 2
        bin2_start = int($5/res) * res
        bin2_end = bin2_start + res
        print $1, bin1_start, bin1_end, $4, bin2_start, bin2_end, $7, $8, $9, $10
    }' fine_res_loops.bedpe > harmonized_loops.bedpe
    ```
    
    ## Step 5: Per-Sample Filtering
    
    ### 5a. ENCODE Blocklist Filtering (Amemiya et al. 2019)
    Remove loops with anchors in artifact-prone regions (download from https://github.com/Boyle-Lab/Blacklist/blob/master/lists/hg38-blacklist.v2.bed.gz):
    ```bash
    # Filter loops where EITHER anchor overlaps a blocklist region
    gunzip -k hg38-blacklist.v2.bed.gz
    
    # First, extract anchor 1 and anchor 2 as separate BED files
    awk 'BEGIN{OFS="\t"} {print $1,$2,$3,NR}' sample.bedpe > anchors1.bed
    awk 'BEGIN{OFS="\t"} {print $4,$5,$6,NR}' sample.bedpe > anchors2.bed
    
    # Find anchor rows NOT in blocklist
    bedtools intersect -a anchors1.bed -b hg38-blacklist.v2.bed -v | cut -f4 > clean_rows_1.txt
    bedtools intersect -a anchors2.bed -b hg38-blacklist.v2.bed -v | cut -f4 > clean_rows_2.txt
    
    # Keep only rows where BOTH anchors pass
    comm -12 <(sort clean_rows_1.txt) <(sort clean_rows_2.txt) > clean_rows.txt
    awk 'NR==FNR{a[$1];next} FNR in a' clean_rows.txt sample.bedpe > sample.filtered.bedpe
    ```
    
    ### 5b. Score Filtering
    Filter by interaction score/significance:
    ```bash
    # If BEDPE has a score column (col 8), filter to significant interactions
    # Keep top 75% by score (true distribution quantile, not range-based)
    TOTAL=$(wc -l < sample.filtered.bedpe)
    LINE_25=$(echo "$TOTAL" | awk '{printf "%d", $1 * 0.25}')
    THRESHOLD=$(sort -k8,8n sample.filtered.bedpe | awk -v line="$LINE_25" 'NR==line{print $8}')
    awk -v t="$THRESHOLD" '$8 >= t' sample.filtered.bedpe > sample.qfiltered.bedpe
    ```
    
    ### 5c. Remove Self-Ligation Artifacts
    Loops where both anchors are very close are likely artifacts:
    ```bash
    # Remove loops where anchors are on same chromosome and < 20kb apart
    awk '{
        if ($1 != $4) print $0;  # inter-chromosomal: keep (rare but real)
        else if (($5 - $3) >= 20000) print $0;  # > 20kb apart: keep
    }' sample.qfiltered.bedpe > sample.clean.bedpe
    ```
    
    ## Step 6: Union Merge of Loops
    
    ### The Paired-Region Matching Problem
    
    Unlike peaks (single regions), loops are **pairs of regions**. Two loops match if **both anchors overlap**:
    
    ```
    Loop 1:  [anchor1A]--------[anchor1B]
    Loop 2:    [anchor2A]------[anchor2B]
    ```
    
    These should merge if anchor1A overlaps anchor2A AND anchor1B overlaps anchor2B.
    
    ### Method A: bedtools pairToPair (Recommended for simple union)
    
    ```bash
    # Concatenate all filtered loops
    cat sample1.clean.bedpe sample2.clean.bedpe ... > all_loops.bedpe
    
    # Sort by anchor 1 coordinates
    sort -k1,1 -k2,2n -k4,4 -k5,5n all_loops.bedpe > all_loops.sorted.bedpe
    
    # Use a custom merge approach:
    # 1. Bin anchors to resolution, creating a loop ID
    # 2. Group by loop ID
    # 3. Count support
    
    awk -v res=10000 'BEGIN{OFS="\t"} {
        # Create binned anchor coordinates as loop identifier
        a1_bin = $1 ":" int($2/res)*res
        a2_bin = $4 ":" int($5/res)*res
        # Canonical order (smaller coordinate first) to handle orientation
        if (a1_bin < a2_bin) loop_id = a1_bin "-" a2_bin
        else loop_id = a2_bin "-" a1_bin
        print loop_id, $0
    }' all_loops.sorted.bedpe | \
    sort -k1,1 | \
    awk 'BEGIN{OFS="\t"} {
        if ($1 != prev_id) {
            if (NR > 1) print chr1, start1, end1, chr2, start2, end2, count, max_score
            prev_id = $1
            chr1=$2; start1=$3; end1=$4; chr2=$5; start2=$6; end2=$7
            count = 1; max_score = $9
        } else {
            count++
            if ($9 > max_score) max_score = $9
            # Expand anchors to encompass all overlapping calls
            if ($3 < start1) start1 = $3
            if ($4 > end1) end1 = $4
            if ($6 < start2) start2 = $6
            if ($7 > end2) end2 = $7
        }
    } END {
        print chr1, start1, end1, chr2, start2, end2, count, max_score
    }' > union_loops.bedpe
    ```
    
    ### Method B: Resolution-Binned Approach (Loop Catalog method)
    
    Following the Loop Catalog (Reyna et al. 2025) approach:
    
    ```bash
    # Bin all loop anchors to a fixed resolution
    awk -v res=10000 'BEGIN{OFS="\t"} {
        a1_start = int($2/res) * res
        a1_end = a1_start + res
        a2_start = int($5/res) * res
        a2_end = a2_start + res
        # Canonical ordering
        if ($1 < $4 || ($1 == $4 && a1_start <= a2_start))
            print $1, a1_start, a1_end, $4, a2_start, a2_end
        else
            print $4, a2_start, a2_end, $1, a1_start, a1_end
    }' all_loops.sorted.bedpe | \
    sort -u | \
    sort -k1,1 -k2,2n -k4,4 -k5,5n | \
    uniq -c | \
    awk 'BEGIN{OFS="\t"} {print $2,$3,$4,$5,$6,$7,$1}' > union_loops_binned.bedpe
    # Columns: chr1, start1, end1, chr2, start2, end2, n_supporting_samples
    ```
    
    ### Method C: Using Specialized Tools
    
    **mariner** (R/Bioconductor):
    ```r
    library(mariner)
    # Read BEDPE files as GInteractions
    loops <- lapply(bedpe_files, read.table)
    # Convert to GInteractions and merge
    gi <- as_ginteractions(loops)
    merged <- mergePairs(gi, radius = 10000)  # 10kb tolerance
    ```
    
    **AQuA Tools** (Python):
    ```python
    # BEDPE union with anchor overlap tolerance
    aqua bedpe-union -i sample1.bedpe sample2.bedpe -o union.bedpe --slop 5000
    ```
    
    ## Step 7: Confidence Annotation
    
    Given N total experiments:
    
    | Confidence | Criteria | Interpretation |
    |-----------|----------|----------------|
    | **High** | Detected in >=50% of samples | Constitutive loop, present across individuals |
    | **Supported** | Detected in 2+ samples | Likely real, some variation |
    | **Singleton** | Detected in 1 sample only | May be individual-specific or depth-dependent |
    
    ```bash
    awk -v N=4 '{
        if ($7 >= N*0.5) conf="HIGH";
        else if ($7 >= 2) conf="SUPPORTED";
        else conf="SINGLETON";
        print $0"\t"conf"\t"$7"/"N
    }' union_loops_binned.bedpe > union_loops.annotated.bedpe
    ```
    
    **Context for singletons**: Hi-C loop detection is very sensitive to sequencing depth. Many singletons may simply be under-powered in other samples rather than biologically absent. The Loop Catalog found that a union approach captures ~3x more loops than any individual experiment.
    
    ## Step 8: Separate Analysis for TADs and Compartments
    
    **TAD boundaries** and **A/B compartments** require different aggregation than loops:
    
    ### TAD Boundaries
    TAD boundaries are single genomic positions. Aggregate like narrow peaks:
    ```bash
    # Extract TAD boundary BED from contact domain files
    # Each boundary is a narrow region
    cat tad_boundaries_sample*.bed | \
    bedtools sort -i - | \
    bedtools merge -i - -d 40000 -c 1 -o count > union_tad_boundaries.bed
    # 40kb gap tolerance because TAD boundaries are resolution-dependent
    ```
    
    ### A/B Compartments
    Compartment calls (eigenvector sign at each bin) should be aggregated by majority vote:
    ```bash
    # For each resolution bin, assign A or B based on majority of samples
    # This is more complex and typically done in R/Python
    ```
    
    ## Step 9: Log Provenance
    
    ```
    encode_log_derived_file(
        file_path="/path/to/union_loops.annotated.bedpe",
        source_accessions=["ENCSR...", "ENCSR...", ...],
        description="Union chromatin loops across N pancreas Hi-C experiments",
        file_type="aggregated_loops",
        tool_used="bedtools + custom merge at 10kb resolution",
        parameters="blocklist filtered, score >= 25th pctl, self-ligation >= 20kb removed, 10kb resolution binning"
    )
    ```
    
    ## Step 10: Summary Statistics
    
    Report to the user:
    - Total input experiments: N
    - Experiments passing QC: M
    - Resolution used: Xkb
    - Total loops before merge: X
    - Union loops after merge: Y
    - High-confidence loops: Z (≥50% support)
    - Supported loops: W (2+ support)
    - Singleton loops: V (1 sample only)
    - Distance distribution: median and range of loop sizes (anchor-to-anchor)
    - Inter-chromosomal loops: count (expect very few)
    
    ## Pitfalls Specific to Hi-C Data
    
    1. **Resolution mismatch**: Loop calls at 5kb vs 25kb resolution will have very different anchor sizes. Always harmonize to a common resolution before merging.
    
    2. **Sequencing depth sensitivity**: Loop calling requires deep sequencing (400M+ valid pairs). Shallowly sequenced experiments will call far fewer loops — this is under-detection, not absence.
    
    3. **Algorithm differences are LARGE**: Wolff et al. 2022 (GigaScience) found that HICCUPS, Mustache, Fit-Hi-C, and HiCExplorer loop callers **intersect by ~50% at most**. Mustache tends to recover more validated loops (Roayaei Ardakany et al. 2020). If mixing callers, note this in provenance — and this discordance is itself a reason to prefer the union approach.
    
    4. **Orientation matters**: BEDPE anchors should be canonically ordered (anchor1 < anchor2 by genomic coordinate) before merging to avoid duplicate counting.
    
    5. **Inter-chromosomal contacts**: These are rare but real. Handle separately — they cannot be distance-filtered.
    
    6. **Distance distribution**: Most loops are 100kb-2Mb. Very short-range contacts (<20kb) are often noise from undigested chromatin. Very long-range (>10Mb) are rare.
    
    7. **Do NOT mix assemblies**: All files must be GRCh38 or all hg19. Hi-C resolution binning makes liftOver of loops particularly error-prone.
    
    8. **TADs vs loops**: These are different features. TADs are domains (regions), loops are point contacts (pairs). Do not mix them in the same union.
    
    9. **Micro-C as complement**: Micro-C achieves higher resolution than Hi-C and can detect sub-TAD loops. Treat Micro-C loops as compatible with Hi-C loops in a union (Mustache works on both).
    
    ## Walkthrough: Building a Cross-Tissue Loop Catalog for the MYC Locus
    
    **Goal**: Aggregate Hi-C chromatin loops across tissues to identify conserved and tissue-specific 3D contacts at the MYC gene locus.
    **Context**: Cancer research — MYC is regulated by distal enhancers via chromatin looping.
    
    ### Step 1: Find Hi-C experiments across tissues
    
    ```
    encode_search_experiments(assay_title="Hi-C", organism="Homo sapiens", limit=50)
    ```
    
    Expected output:
    ```json
    {
      "results": [
        {"accession": "ENCSR000AKA", "assay_title": "Hi-C", "biosample_summary": "GM12878", "status": "released"},
        {"accession": "ENCSR489OCU", "assay_title": "Hi-C", "biosample_summary": "K562", "status": "released"},
        {"accession": "ENCSR382RFU", "assay_title": "Hi-C", "biosample_summary": "liver", "status": "released"}
      ],
      "total": 89,
      "limit": 50,
      "offset": 0,
      "has_more": true,
      "next_offset": 50
    }
    ```
    
    **Interpretation**: 89 Hi-C experiments available. Select 5–10 spanning diverse tissue types for cross-tissue comparison.
    
    ### Step 2: List loop files for each experiment
    
    ```
    encode_list_files(experiment_accession="ENCSR000AKA", file_format="bedpe", assembly="GRCh38")
    ```
    
    Expected output (a JSON array of file records; fields abridged):
    ```json
    [
      {"accession": "ENCFF001ABC", "output_type": "contact domains", "file_format": "bedpe", "file_size_human": "2.4 MB"},
      {"accession": "ENCFF002DEF", "output_type": "loops", "file_format": "bedpe", "file_size_human": "1.8 MB"}
    ]
    ```
    
    **Interpretation**: Use "loops" files for loop aggregation. Contact domains are TADs, not loops.
    
    ### Step 3: Download loop files
    
    ```
    encode_download_files(file_accessions=["ENCFF002DEF", "ENCFF003GHI", "ENCFF004JKL"], download_dir="/data/hic_loops")
    ```
    
    Expected output (one of the three `downloaded` entries shown):
    ```json
    {
      "downloaded": [
        {
          "accession": "ENCFF002DEF",
          "file_path": "/data/hic_loops/ENCFF002DEF.bedpe",
          "file_size": 1887436,
          "file_size_human": "1.8 MB",
          "success": true,
          "error": "",
          "md5_verified": true
        }
      ],
      "errors": [],
      "summary": {
        "total_requested": 3,
        "successful": 3,
        "failed": 0,
        "total_size": 5872025,
        "total_size_human": "5.6 MB"
      }
    }
    ```
    
    ### Step 4: Aggregate loops with resolution-aware anchor matching
    
    Apply union merge across tissues:
    - Expand loop anchors by ±resolution (e.g., ±5kb for 5kb resolution data)
    - Merge overlapping anchors using bedtools pairToPair
    - Assign tissue support counts to each union loop
    - Filter: require ≥2 tissue support for conserved loops
    
    ### Step 5: Filter to MYC locus
    
    ```bash
    # MYC locus: chr8:127,700,000-128,000,000
    awk '$1=="chr8" && $2>=127700000 && $3<=128000000' union_loops.bedpe > myc_loops.bedpe
    ```
    
    **Interpretation**: Loops anchored at the MYC promoter connecting to distal enhancers. Conserved loops (≥3 tissues) likely represent fundamental regulatory architecture; tissue-specific loops may drive context-dependent MYC activation.
    
    ### Integration with downstream skills
    - Feed loop anchors into → **peak-annotation** for gene assignment at anchor regions
    - Overlay with → **histone-aggregation** H3K27ac peaks to identify active enhancer-promoter loops
    - Cross-reference loop-disrupting variants via → **variant-annotation**
    - Visualize in → **ucsc-browser** as interaction tracks
    
    ## Code Examples
    
    ### 1. Survey available Hi-C data by tissue
    ```
    encode_get_facets(assay_title="Hi-C", organism="Homo sapiens")
    ```
    
    Expected output (facet field names are the top-level keys):
    ```json
    {
      "biosample_ontology.organ_slims": [
        {"term": "brain", "count": 24},
        {"term": "blood", "count": 15}
      ]
    }
    ```
    
    ### 2. Get details for a specific Hi-C experiment
    ```
    encode_get_experiment(accession="ENCSR000AKA")
    ```
    
    Expected output (fields abridged):
    ```json
    {
      "accession": "ENCSR000AKA",
      "assay_title": "Hi-C",
      "biosample_summary": "GM12878",
      "assembly": ["GRCh38"],
      "bio_replicate_count": 2,
      "status": "released",
      "lab": "Erez Lieberman Aiden, Baylor",
      "audit_error_count": 0,
      "audit_not_compliant_count": 0,
      "audit_warning_count": 1,
      "audit_internal_action_count": 0
    }
    ```
    
    ### 3. Compare loop sets between two cell types
    
    Both experiments must already be tracked; otherwise the tool returns `{"error": "Experiment ... not tracked. Track it first."}`.
    ```
    encode_compare_experiments(accession1="ENCSR000AKA", accession2="ENCSR489OCU")
    ```
    
    Expected output:
    ```json
    {
      "experiment_1": {"accession": "ENCSR000AKA", "assay": "Hi-C", "biosample": "GM12878"},
      "experiment_2": {"accession": "ENCSR489OCU", "assay": "Hi-C", "biosample": "K562"},
      "verdict": "COMPATIBLE_WITH_CAVEATS",
      "recommendation": "These experiments can be compared, but the warnings should be addressed in your analysis.",
      "compatible_aspects": [
        "Same organism: Homo sapiens",
        "Same assembly: GRCh38",
        "Same assay: Hi-C",
        "Same biosample type: cell line"
      ],
      "issues": [],
      "warnings": [
        "Different labs: Erez Lieberman Aiden, Baylor vs Job Dekker, UMass. Batch effects possible."
      ]
    }
    ```
    
    ## Integration
    
    | This skill produces... | Feed into... | Purpose |
    |---|---|---|
    | Union loop catalog (BEDPE) | **peak-annotation** | Assign genes to loop anchors |
    | Conserved loop coordinates | **histone-aggregation** | Overlay H3K27ac at anchors to find active enhancer-promoter loops |
    | Tissue-specific loops | **accessibility-aggregation** | Check if loop anchors overlap open chromatin |
    | Loop anchor BED intervals | **variant-annotation** | Find GWAS/clinical variants disrupting loop anchors |
    | Loop anchor coordinates | **liftover-coordinates** | Convert hg19 loops to GRCh38 |
    | Aggregated loop statistics | **visualization-workflow** | Generate loop frequency heatmaps |
    | Loop-gene assignments | **disease-research** | Connect loop disruptions to disease phenotypes |
    
    ## Related Skills
    
    - **histone-aggregation**: Loop anchors often overlap with H3K27ac/H3K4me1 peaks — integrate with histone union sets to annotate loop function
    - **accessibility-aggregation**: Loop anchors frequently coincide with accessible chromatin — validate loops by requiring anchor accessibility
    - **regulatory-elements**: Use loops to connect distal enhancers (H3K27ac) to target promoters (H3K4me3)
    - **epigenome-profiling**: Loops add 3D context to 1D chromatin state maps
    - **pipeline-hic**: Process raw Hi-C data through the full ENCODE-aligned pipeline
    - **batch-analysis**: Batch processing workflows for systematic Hi-C loop aggregation
    - **publication-trust**: Verify literature claims backing analytical decisions
    
    ## Presenting Results
    
    - Present aggregated loops as: chr | anchor1_start | anchor1_end | anchor2_start | anchor2_end | sample_count | resolution. Show loop statistics. Suggest: "Would you like to check if any GWAS variants overlap loop anchors?"
    
    ## For the request: "$ARGUMENTS"
    

Comments (0)

Sign in to join the conversation.

No comments yet.

Reviews (0)

No reviews yet.

Related