Analyze ENCODE Functional Genomics Screens

When to Use

• User wants to find or analyze CRISPR screen, MPRA, or STARR-seq data from ENCODE
• User asks about "functional screens", "CRISPR perturbation", "reporter assay", or "enhancer validation"
• User needs to identify functionally validated regulatory elements from screen results
• User wants to integrate screen results with epigenomic annotations (ChIP-seq, ATAC-seq peaks)
• Example queries: "find CRISPR screen data in ENCODE", "analyze MPRA results for enhancer activity", "which regulatory elements have functional validation?"

Discover and interpret functional validation data from CRISPR screens, MPRA (Massively Parallel Reporter Assays), and STARR-seq experiments in the ENCODE catalog. These assays directly test whether candidate regulatory elements have functional activity, complementing the correlative evidence from ChIP-seq, ATAC-seq, and Hi-C.

Scientific Rationale

The question: "Which of the candidate regulatory elements identified by ENCODE actually have functional activity, and what genes do they regulate?"

The central challenge in regulatory genomics is that biochemical signatures (histone marks, chromatin accessibility, TF binding) are correlative — they identify candidate regulatory elements but cannot prove function. ENCODE Phase 4 addressed this gap by investing heavily in functional characterization: large-scale CRISPR perturbation screens, MPRA experiments testing thousands of candidate elements in parallel, and STARR-seq for genome-wide enhancer activity mapping.

The Validation Gap

ENCODE catalogs 926,535 human candidate cis-regulatory elements (cCREs). But how many of these are truly functional?

• CRISPR screens (Gasperini et al. 2019): Of 5,920 candidate enhancers tested by CRISPRi, only ~12% showed significant effects on nearby gene expression
• MPRA (Inoue et al. 2017; Tewhey et al. 2016): Reporter assays confirm activity for 40-60% of predicted enhancers, depending on cell type and element class
• STARR-seq (Arnold et al. 2013): Genome-wide enhancer assays identify thousands of active elements, but episomal context differs from chromosomal

These functional assays provide the strongest evidence (short of genetic studies in humans) that a regulatory element has biological activity. ENCODE4 has scaled these approaches: the Functional Characterization Centers (Yao et al. 2024) performed 108 CRISPRi screens with >540,000 perturbations, targeting 3.27 million ENCODE SCREEN cCREs.

Assay Comparison

Assay	Tests	Context	Scale	Confidence	Key Limitation
CRISPR screen (CRISPRi/CRISPRa)	Endogenous perturbation	Native chromatin	5,000–500,000 elements	Highest	Limited to cell lines; delivery constraints
MPRA	Reporter activity	Episomal (plasmid)	10,000–100,000 variants	High for activity	Removed from chromatin context
STARR-seq	Self-transcription	Episomal (plasmid)	Genome-wide library	High for activity	Episomal; position effects

Literature Support

• Gasperini et al. 2019 (Cell, ~800 citations): CRISPRi screen of 5,920 candidate enhancers with single-cell RNA-seq readout in K562 cells. Identified 664 enhancer-gene pairs. Established the "crisprQTL" framework linking perturbation effects to gene expression at single-cell resolution. DOI
• Fulco et al. 2019 (Nature Genetics, ~1,200 citations): CRISPRi tiling screen combined with the Activity-By-Contact (ABC) model. Demonstrated that ABC predictions outperform distance-based enhancer-gene assignment. Quantitative relationship between enhancer activity, contact frequency, and gene regulation. DOI
• Shalem et al. 2014 (Science, ~5,500 citations): Genome-scale CRISPR-Cas9 knockout screening. The foundational paper for CRISPR loss-of-function screens. Established library design principles and analytical frameworks. DOI
• Arnold et al. 2013 (Science, ~1,100 citations): STARR-seq — Self-Transcribing Active Regulatory Region sequencing. First genome-wide quantitative enhancer activity assay. Tests millions of fragments simultaneously in Drosophila; adapted to human. DOI
• Inoue et al. 2017 (Genome Research, ~400 citations): MPRA for systematic variant effect prediction. Tested thousands of regulatory variants for allele-specific enhancer activity. Demonstrated that GWAS risk alleles frequently alter enhancer function. DOI
• Li et al. 2014 (Genome Biology, ~2,800 citations): MAGeCK — Model-based Analysis of Genome-wide CRISPR-Cas9 Knockout. The standard computational tool for CRISPR screen analysis. Robust negative binomial model for guide RNA count data. DOI
• Gordon et al. 2020 (Nature Protocols, ~200 citations): MPRAflow — standardized computational pipeline for MPRA data analysis. Reproducible barcode counting, normalization, and activity scoring. DOI
• Tewhey et al. 2016 (Cell, ~600 citations): High-throughput identification of regulatory variants using MPRA. Tested >30,000 allelic pairs across 3,642 GWAS loci. Identified hundreds of variants with allele-specific regulatory activity. DOI
• Nasser et al. 2021 (Nature, ~700 citations): ABC model for enhancer-gene prediction using ENCODE data. Linked 5,036 GWAS signals to 2,249 genes across 131 cell types. Validated predictions against CRISPRi perturbation data. DOI
• Klein et al. 2020 (Nature Genetics, ~350 citations): CRISPRi screen design principles for non-coding regulatory elements. Established that guide RNA positioning relative to regulatory element boundaries is critical. Quantified the "shadow" of CRISPRi repression (~1–2 kb). DOI
• Yao et al. 2024 (Nature Methods, ~26 citations): ENCODE4 Functional Characterization Centers — 108 CRISPRi screens, >540,000 perturbations across multiple cell types. Pre-designed sgRNA library targeting 3.27M ENCODE SCREEN cCREs. Establishes the largest functional characterization dataset for regulatory elements. DOI
• Lee et al. 2020 (Genome Biology, ~150 citations): STARRPeaker — peak caller designed specifically for STARR-seq data. Handles input library normalization and identifies significant enhancer peaks from STARR-seq enrichment. DOI
• Kim & Hart 2021 (Genome Medicine, ~200 citations): BAGEL2 — Bayesian Analysis of Gene Essentiality. Updated framework for identifying essential genes and functional elements from CRISPR screen data. DOI

Finding ENCODE Screen Data

CRISPR Screens

encode_search_experiments(assay_title="CRISPR screen")

ENCODE contains CRISPRi (inhibition) and CRISPRa (activation) screens targeting regulatory elements:

Screen Type	Mechanism	Effect on Target	Use Case
CRISPRi (dCas9-KRAB)	Transcriptional repression	Silences enhancer/promoter	Loss-of-function; identifies required elements
CRISPRa (dCas9-VP64/p65)	Transcriptional activation	Activates latent elements	Gain-of-function; identifies sufficient elements
CRISPR knockout	Cas9 nuclease	Deletes element	Irreversible loss-of-function

Typical ENCODE CRISPR screen outputs:

• Guide RNA quantifications (sgRNA counts per condition)
• Element quantifications (aggregated guide effects per target element)
• Differential expression results

# List available files for a CRISPR screen experiment
encode_list_files(
    experiment_accession="ENCSR...",
    file_format="tsv",
    assembly="GRCh38"
)

MPRA Experiments

encode_search_experiments(assay_title="MPRA")

ENCODE MPRA experiments test candidate cis-regulatory elements for enhancer/promoter activity:

Typical MPRA outputs:

• Barcode count matrices (RNA and DNA)
• Activity scores per tested element
• Differential activity between alleles (for variant testing)

STARR-seq Experiments

encode_search_experiments(assay_title="STARR-seq")

STARR-seq tests enhancer activity genome-wide using self-transcribing reporter constructs:

Typical STARR-seq outputs:

• Aligned reads (BAM)
• Signal tracks (bigWig) — enrichment over input
• Peak files (BED) — identified enhancer elements

Combined Discovery

To find ALL functional characterization data for a tissue or cell type:

# All perturbation experiments
encode_search_experiments(perturbed=True, organ="...")

# Functional characterization in a specific cell line
encode_search_experiments(assay_title="CRISPR screen", biosample_term_name="K562")
encode_search_experiments(assay_title="MPRA", biosample_term_name="K562")
encode_search_experiments(assay_title="STARR-seq", biosample_term_name="K562")

CRISPR Screen Analysis

Data Types and File Formats

Data Type	Format	Description
Guide RNA counts	TSV	Raw or normalized sgRNA counts per sample
Element quantifications	TSV	Aggregated effect sizes per target element
Differential expression	TSV	Genes with significant expression changes

Analytical Workflow

Step 1: Guide counts     → Quality filter low-representation guides
Step 2: Normalization     → Median ratio or total count normalization
Step 3: Statistical test  → MAGeCK, BAGEL2, or custom model
Step 4: Hit calling       → FDR correction, effect size thresholds
Step 5: Integration       → Overlay on ENCODE cCREs and epigenomic marks

Key Analysis Tools

MAGeCK (Li et al. 2014)

The standard tool for CRISPR screen analysis:

# Count sgRNAs from FASTQ
mageck count -l library.tsv -n experiment \
    --sample-label "control,treatment" \
    --fastq control_R1.fastq.gz treatment_R1.fastq.gz

# Test for enrichment/depletion
mageck test -k experiment.count.txt \
    -t treatment -c control \
    -n results --remove-zero both

MAGeCK outputs:

• results.gene_summary.txt — Gene/element-level results (RRA and MLE)
• results.sgrna_summary.txt — Individual guide-level results
• Key columns: neg|score, neg|fdr, pos|score, pos|fdr

BAGEL2 (Kim & Hart 2021)

Bayesian framework for gene essentiality from CRISPR screens:

# Calculate Bayes Factors
BAGEL.py fc -i counts.txt -o foldchange.txt -c control_columns
BAGEL.py bf -i foldchange.txt -o bayes_factors.txt \
    -e essential_genes.txt -n nonessential_genes.txt
BAGEL.py pr -i bayes_factors.txt -o precision_recall.txt

Quality Control Metrics

QC Metric	Threshold	Description
Guide representation	>200 reads/guide	Minimum coverage for statistical power
Replicate correlation	Pearson r > 0.7	Between biological replicates
Positive control enrichment	p < 0.01	Known essential genes/elements should score
Negative control depletion	~50% at FDR 0.1	Random non-targeting guides show no effect
Gini index	<0.3	Measures guide count distribution evenness
Mapping rate	>70%	Reads mapping to library sequences

CRISPRi-Specific Considerations

• Repression window: CRISPRi silences a ~1–2 kb window around the guide target (Klein et al. 2020). Multiple guides within this window are not independent.
• TSS proximity bias: Elements near gene TSSs show stronger effects due to proximity to promoters. Control for distance-to-TSS in analysis.
• Multiple guide averaging: Require 2+ guides per element with concordant effects to call a hit.

Integration with ENCODE Epigenomic Data

After identifying screen hits:

# Overlay CRISPRi hits on ENCODE cCREs
encode_search_experiments(assay_title="Histone ChIP-seq", target="H3K27ac", biosample_term_name="K562")

# Check chromatin accessibility at hit locations
encode_search_experiments(assay_title="ATAC-seq", biosample_term_name="K562")

# Download peak files for intersection
encode_list_files(
    experiment_accession="ENCSR...",
    file_format="bed",
    output_type="IDR thresholded peaks",
    assembly="GRCh38",
    preferred_default=True
)

Intersection analysis:

# Intersect CRISPR hits with H3K27ac peaks
bedtools intersect -a crispr_hits.bed -b h3k27ac_peaks.bed -wa -wb > hits_in_enhancers.bed

# Calculate enrichment of hits in specific chromatin states
# Compare: (hits in enhancers / total hits) vs (all tested elements in enhancers / total tested)

MPRA Analysis

Data Types and File Formats

Data Type	Format	Description
Barcode counts (DNA)	TSV	Input library representation
Barcode counts (RNA)	TSV	Transcriptional output
Activity scores	TSV	RNA/DNA ratio per element

Analytical Workflow

Step 1: Barcode counting   → Count barcodes in DNA and RNA libraries
Step 2: DNA normalization  → Normalize RNA counts by DNA representation
Step 3: Activity scoring   → Calculate RNA/DNA ratio per element
Step 4: Statistical test   → Compare to negative controls
Step 5: Allelic comparison → Test for allele-specific activity (if applicable)

MPRAflow (Gordon et al. 2020)

Standardized computational pipeline:

# Run MPRAflow
nextflow run MPRAflow/MPRAflow.nf \
    --design design_file.txt \
    --fastq_insert insert_reads/ \
    --fastq_bc barcode_reads/ \
    --outdir results/

MPRAflow handles:

• Barcode-to-element association
• Count normalization
• Activity scoring with confidence intervals
• Allele-specific comparisons

Activity Score Interpretation

The core MPRA measurement is the activity ratio:

Activity = log2(RNA_counts / DNA_counts)

Activity Score	Interpretation
Activity >> 0 (e.g., >1.5)	Strong enhancer/promoter activity
Activity ~ 0	No regulatory activity (or activity equal to minimal promoter)
Activity << 0 (e.g., <-1.0)	Potential silencer activity (reduces transcription)

Statistical Framework

• Null distribution: Negative control elements (scrambled sequences) define background
• Test: Compare each element's activity to the negative control distribution
• Multiple testing: Benjamini-Hochberg FDR across all tested elements
• Effect size: log2(activity ratio) relative to negative controls
• Significance threshold: FDR < 0.05 AND |log2FC| > 1.0 (typical, adjust per experiment)

Integration with ENCODE Annotations

# Correlate MPRA activity with ENCODE histone mark signals
# Elements with strong MPRA activity should show:
# - High H3K27ac signal (active enhancer mark)
# - Chromatin accessibility (ATAC-seq/DNase-seq signal)
# - TF binding (ChIP-seq signal at element)

# Elements with NO MPRA activity despite H3K27ac may be:
# - Context-dependent (active only in specific conditions)
# - False-positive biochemical marks
# - Silencer elements (negative activity in MPRA)

STARR-seq Analysis

Data Types and File Formats

Data Type	Format	Description
Input library	BAM/FASTQ	Cloned genomic fragments
STARR-seq output	BAM/FASTQ	Self-transcribed fragments (enriched)
Signal tracks	bigWig	Enrichment over input
Peaks	BED	Identified enhancer elements

Analytical Workflow

Step 1: Alignment         → Map input and STARR-seq reads to genome
Step 2: Enrichment        → Calculate STARR-seq / input ratio
Step 3: Peak calling      → Identify enriched regions (enhancers)
Step 4: Quantification    → Measure enhancer strength per peak
Step 5: Integration       → Compare with ENCODE cCRE predictions

STARRPeaker (Lee et al. 2020)

Purpose-built peak caller for STARR-seq data:

# Call peaks from STARR-seq
starrpeaker \
    --prefix output_prefix \
    --chromsize hg38.chrom.sizes \
    --bam input.bam starrseq.bam \
    --threshold 0.05

Alternatively, standard peak callers can be applied to the enrichment:

# Using MACS2 on STARR-seq enrichment
macs2 callpeak -t starrseq.bam -c input.bam \
    -f BAM -g hs --nomodel \
    -n starr_enhancers -q 0.05

Comparison with ENCODE cCRE Predictions

STARR-seq Result	ENCODE cCRE Status	Interpretation
Active in STARR-seq	dELS or pELS	Validated enhancer
Active in STARR-seq	PLS	Promoter with enhancer activity
Active in STARR-seq	Not in cCRE catalog	Novel enhancer (or context-dependent)
Inactive in STARR-seq	dELS or pELS	Possible false positive cCRE, or context-dependent
Inactive in STARR-seq	Not in cCRE catalog	Confirmed non-enhancer

STARR-seq Caveats

• Episomal context: STARR-seq tests elements outside their native chromatin. Elements requiring specific chromatin context may fail to score.
• Position effects: The reporter construct itself introduces bias (minimal promoter, polyA signal placement).
• Quantitative, not binary: STARR-seq enrichment is continuous. Use a statistical threshold (e.g., input-normalized enrichment > 2-fold, FDR < 0.05) rather than arbitrary cutoffs.
• Species considerations: STARR-seq was originally developed in Drosophila (Arnold et al. 2013). Human STARR-seq requires modifications (ori-free constructs) to reduce plasmid replication artifacts.

Integrating Screens with Epigenomic Data

Overlay on ChromHMM States

For each functionally validated element, determine its chromatin state:

# Retrieve ChromHMM annotations for the cell type
# (Pre-computed by Roadmap Epigenomics for 111 reference epigenomes)

# Enrichment analysis: are screen hits preferentially in specific chromatin states?
# Expected: hits enriched in active enhancer (Enh, EnhG) and active promoter (TssA) states
# Unexpected enrichment in quiescent or heterochromatin states suggests novel regulatory mechanisms

Correlate with ENCODE cCRE Classes

cCRE Class	Expected Screen Hit Rate	Rationale
PLS (Promoter-like)	30–50%	Promoters are consistently active
pELS (Proximal enhancer-like)	15–25%	Proximity to promoters increases detection
dELS (Distal enhancer-like)	5–15%	Distal enhancers are often cell-type-specific
CTCF-only	<5%	Insulators rarely show enhancer activity
No cCRE overlap	1–3%	Novel elements or context-dependent activity

ABC Model Integration (Nasser et al. 2021)

The Activity-By-Contact model predicts enhancer-gene links using ENCODE data. Cross-reference screen hits:

# For each CRISPR hit:
# 1. Check if the hit is predicted by ABC to regulate the observed target gene
# 2. ABC-predicted enhancers that are also CRISPR-validated are highest confidence
# 3. CRISPR hits NOT predicted by ABC may act through mechanisms ABC does not model

GWAS Variant Enrichment in Screen-Validated Elements

Screen-validated elements provide the strongest evidence for variant interpretation:

# Workflow:
# 1. Identify GWAS variants in LD (r2 > 0.8) from gwas-catalog skill
# 2. Intersect with functionally validated enhancers
# 3. Variants in CRISPR-validated enhancers are highest-priority causal candidates
# 4. Test for enrichment: are GWAS variants over-represented in screen hits?

bedtools intersect \
    -a gwas_variants_ld.bed \
    -b crispr_validated_enhancers.bed \
    -wa -wb > gwas_in_validated_enhancers.bed

Experimental Design Considerations

Library Complexity

Screen Type	Minimum Library Size	Recommended Coverage
CRISPR (gene-level)	4–6 guides per gene	500x per guide
CRISPR (element-tiling)	1 guide per ~100bp	200x per guide
MPRA	10–20 barcodes per element	100x per barcode
STARR-seq	Genome-wide fragmentation	10x genome coverage

Control Element Selection

Positive controls (expected to score):

• Known essential genes (for knockout screens)
• Validated enhancers from VISTA or previous CRISPR studies
• Strong constitutive promoters (e.g., CMV, EF1a for MPRA)

Negative controls (expected to show no effect):

• Non-targeting sgRNAs (CRISPR) — 100–1,000 random guides
• Scrambled sequences (MPRA) — sequence-matched but shuffled
• Gene deserts (STARR-seq) — genomic regions far from genes

Cell Type Matching

Critical: The functional screen must be performed in a cell type that is biologically relevant to the regulatory elements being tested. An enhancer active in hepatocytes may show no MPRA activity in HEK293 cells.

ENCODE Screen Cell Types	Tissue Relevance
K562 (CML)	Hematopoietic lineage, Tier 1
HepG2	Liver / hepatocyte
GM12878	B-lymphocyte, Tier 1
WTC-11 (iPSC-derived)	Multiple differentiated cell types
A549	Lung epithelial

When integrating screen data with ENCODE epigenomic data:

# Ensure cell type match
encode_search_experiments(assay_title="Histone ChIP-seq", target="H3K27ac", biosample_term_name="K562")
# Match the screen cell type exactly for meaningful correlation

Power Analysis

| Elements Tested | Guides/Barcodes Per | Replicates | Expected Power (FDR<0.05, |FC|>1.5) | |----------------|-------------------|-----------|--------------------------------------| | 5,000 | 4 guides | 3 bio | ~80% for strong effects | | 5,000 | 6 guides | 3 bio | ~90% for moderate effects | | 50,000 | 2 guides | 2 bio | ~60% for strong effects only |

Pitfalls & Edge Cases

• CRISPR screen vs MPRA measure different things: CRISPR screens test necessity (loss-of-function), while MPRA tests sufficiency (gain-of-function). An element can be MPRA-active but CRISPR-dispensable if redundant enhancers compensate.
• MPRA context effects: MPRA tests elements on episomal plasmids outside their native chromatin context. Elements that require specific chromatin states or 3D contacts may score negative in MPRA but be active in vivo.
• STARR-seq input library bias: STARR-seq results depend heavily on input library complexity. Low-complexity libraries produce false negatives. Always check input coverage before interpreting negative results.
• Cell-type specificity of screens: A CRISPR screen in K562 may miss enhancers active only in primary tissues. Screen results are specific to the cell type used — do not assume they generalize across tissues.
• Effect size vs statistical significance: In pooled CRISPR screens, large guide RNA libraries can produce statistically significant but biologically trivial effect sizes. Always set minimum effect size thresholds alongside p-value cutoffs.
• Guide RNA efficiency confounds: Poor guide RNA cutting efficiency creates false negatives in CRISPR screens. Use multiple guides per element and aggregate with MAGeCK or similar tools.

Provenance Integration

Log all screen analysis operations:

encode_track_experiment(
    accession="ENCSR...",
    notes="CRISPR screen analysis for [cell type] regulatory elements"
)

encode_log_derived_file(
    file_path="/path/to/screen_results.tsv",
    source_accessions=["ENCSR...", "ENCFF..."],
    description="CRISPR screen hit list: [N] significant elements at FDR<0.05 from [total] tested in [cell type]",
    file_type="screen_results",
    tool_used="MAGeCK v0.5.9.5",
    parameters="mageck test -t treatment -c control --remove-zero both; FDR<0.05, |LFC|>0.5"
)