Homer Motif Analysis

# Install HOMER via conda (recommended — handles Perl dependencies)
conda install -c bioconda homer

# Verify installation
findMotifsGenome.pl 2>&1 | head -3
# Usage: findMotifsGenome.pl <peak/BED file> <genome> <output directory> [options]

annotatePeaks.pl 2>&1 | head -3
# Usage: annotatePeaks.pl <peak/BED file> <genome> [options]

# Install reference genomes (downloads 2-way masker + sequence; ~3–10 GB each)
installGenome.pl hg38
installGenome.pl mm10

# Install Python parsing dependencies
pip install pandas matplotlib seaborn

decisions:
  - id: D1
    param: genomeAssembly
    kind: derived
    source: upstream
    ask: "Which assembly are the peak coordinates on?"
    default: "carried from the peak-calling stage"

  - id: D2
    param: searchWindow
    kind: required
    source: user
    ask: "How wide a window around each peak centre should be searched?"
    default: "200 bp for TF ChIP, 150 bp for ATAC"

  - id: D3
    param: backgroundRegions
    kind: required
    source: user
    ask: "Compare peaks against GC-matched random genomic regions, or against a specific control set?"
    default: "auto-generated GC-matched background"

  - id: D4
    param: repeatMasking
    kind: required
    source: user
    ask: "Mask repetitive sequence, which otherwise dominates the enrichment as false positives?"
    default: "masked"

  - id: D5
    param: denovoDiscovery
    kind: required
    source: user
    ask: "Search for previously undescribed motifs, or only test the known-motif library?"
    default: "known motifs only - de novo discovery is roughly ten times slower"

  - id: D6
    param: motifLengths
    kind: optional
    source: user
    depends_on: [D5]
    ask: "Which motif widths should the de novo search try?"
    default: "8, 10, 12"
    skip_if: "de novo discovery disabled"

  - id: D7
    param: numDenovoMotifs
    kind: optional
    source: user
    depends_on: [D5]
    ask: "How many de novo motifs should be reported?"
    default: 25
    skip_if: "de novo discovery disabled"

  - id: D8
    param: knownMotifMismatches
    kind: optional
    source: user
    ask: "How loosely may a sequence match a known motif and still count?"
    default: 2

  - id: D9
    param: threads
    kind: never_ask
    source: data
    reason: "Affects runtime only, not the enrichment"
    default: "min(8, available_cores)"

# Run de novo + known motif enrichment on TF ChIP-seq peaks (hg38, 200 bp window)
findMotifsGenome.pl peaks/tf_chip_summits.bed hg38 motif_output/ \
    -size 200 -mask -p 4

# Annotate peaks with nearest genes and genomic features
annotatePeaks.pl peaks/tf_chip_peaks.narrowPeak hg38 > annotated_peaks.txt

echo "Top known motif:"
head -2 motif_output/knownResults.txt | tail -1 | cut -f1-4

echo "Annotated peaks: $(wc -l < annotated_peaks.txt) lines"

# Activate conda environment (or use existing env)
conda create -n homer_env -c bioconda homer python=3.10 -y
conda activate homer_env

# List available genomes
installGenome.pl list

# Install human (hg38) and mouse (mm10) genomes
# Downloads masked genome sequence and annotation files
installGenome.pl hg38
# Output: Installing hg38... Done. (3-5 min, ~3 GB)

installGenome.pl mm10
# Output: Installing mm10... Done. (3-5 min, ~2.5 GB)

# Verify genome is installed
ls ~/.homer/data/genomes/hg38/
# genome.fa  chrom.sizes  ...

# Check HOMER motif database
ls ~/.homer/data/knownTFs/
# vertebrates.motifs  jaspar.motifs  ...

# Option A: Use MACS3 summit file directly (already 1 bp summit positions)
# Expand summits to ±100 bp (200 bp total) centered on summit
awk 'BEGIN{OFS="\t"} {
    start = ($2 - 100 < 0) ? 0 : $2 - 100;
    print $1, start, $2 + 100, $4, $5
}' peaks/tf_chip_summits.bed > peaks/tf_chip_200bp.bed

echo "Summit-centered peaks: $(wc -l < peaks/tf_chip_200bp.bed)"
# Summit-centered peaks: 12453

# Option B: Use narrowPeak file directly (HOMER accepts multi-column BED)
# HOMER uses columns 1-3 (chr, start, end) and centers internally with -size
cp peaks/tf_chip_peaks.narrowPeak peaks/input_peaks.bed

# Option C: Prepare a custom background region file (matched GC content)
# HOMER auto-generates background if not provided, but explicit background
# is recommended when comparing two peak sets
# Use the control peak set or random genomic regions as background:
bedtools shuffle -i peaks/tf_chip_peaks.narrowPeak \
    -g ~/.homer/data/genomes/hg38/chrom.sizes \
    -excl peaks/tf_chip_peaks.narrowPeak > peaks/background_regions.bed

echo "Background regions: $(wc -l < peaks/background_regions.bed)"
# Background regions: 12453

mkdir -p motif_output/

# Full run: de novo + known motif enrichment
# -size 200: use 200 bp window centered on peak midpoint
# -mask: mask repetitive elements (recommended for clean motifs)
# -p 4: use 4 CPU threads
# -S 25: find top 25 de novo motifs (default)
findMotifsGenome.pl peaks/tf_chip_200bp.bed hg38 motif_output/ \
    -size 200 \
    -mask \
    -p 4 \
    -S 25

# Check progress output:
# Reading genome sizes for hg38 ...
# Scanning for motifs...
# Optimizing 25 motifs...
# Done! Output in motif_output/

echo "Known results: $(wc -l < motif_output/knownResults.txt) motifs"
echo "De novo motifs: $(ls motif_output/homerResults/*.motif 2>/dev/null | wc -l) motifs"
# Known results: 392 motifs
# De novo motifs: 25 motifs

# For mouse peaks (mm10)
# findMotifsGenome.pl peaks/atac_peaks.bed mm10 motif_output_mm10/ \
#     -size 200 -mask -p 4

# Skip de novo discovery (faster when you only need known motifs)
findMotifsGenome.pl peaks/tf_chip_200bp.bed hg38 motif_output_known/ \
    -size 200 \
    -mask \
    -p 4 \
    -nomotif

echo "Known motif results: $(wc -l < motif_output_known/knownResults.txt)"
# Known motif results: 392

# Find occurrences of a specific motif across peaks (outputs peak-level annotation)
# Extract the motif matrix file for the TF of interest from homerResults/
findMotifsGenome.pl peaks/tf_chip_200bp.bed hg38 motif_scan_out/ \
    -size 200 \
    -mask \
    -find motif_output/homerResults/motif1.motif \
    > peaks_with_motif1.txt

echo "Peaks containing motif1: $(wc -l < peaks_with_motif1.txt)"
# Peaks containing motif1: 8941

# Custom background: compare treated vs. control peak sets
findMotifsGenome.pl peaks/treated_peaks.bed hg38 motif_treated_vs_ctrl/ \
    -size 200 \
    -mask \
    -p 4 \
    -bg peaks/control_peaks.bed

# Annotate peaks with nearest gene and TSS distance
# Outputs a tab-delimited file with genomic context for each peak
annotatePeaks.pl peaks/tf_chip_peaks.narrowPeak hg38 \
    > annotated_peaks.txt

echo "Annotated peaks: $(($(wc -l < annotated_peaks.txt) - 1)) peaks"
# Annotated peaks: 12453 peaks

# Preview column headers and first peak
head -2 annotated_peaks.txt | cut -f1-10

# Annotate ATAC-seq peaks (same command, different input)
annotatePeaks.pl peaks/atac_sample_peaks.narrowPeak hg38 \
    > annotated_atac.txt

# Generate TSS-distance histogram (for tag density plots)
# annotatePeaks.pl can compute read density around peaks with -d flag
# annotatePeaks.pl tss hg38 -size 4000 -hist 10 \
#     -d chip_tagdir/ > tss_histogram.txt

import pandas as pd
import subprocess
import io

# --- Parse known motif enrichment results ---
# knownResults.txt columns:
# Motif Name | Consensus | P-value | Log P-value | q-value | # Target Seqs w/ motif | % Target | # Bg Seqs w/ motif | % Bg
known_cols = [
    "motif_name", "consensus", "pvalue", "log_pvalue",
    "qvalue", "n_target_seqs", "pct_target",
    "n_bg_seqs", "pct_bg"
]
known = pd.read_csv(
    "motif_output/knownResults.txt",
    sep="\t", header=0, names=known_cols
)

# Convert string percentages to floats
known["pct_target"] = known["pct_target"].str.rstrip("%").astype(float)
known["pct_bg"] = known["pct_bg"].str.rstrip("%").astype(float)
known["fold_enrichment"] = known["pct_target"] / known["pct_bg"].replace(0, 0.001)
known["-log10_pvalue"] = -known["log_pvalue"] / 2.303  # log10 from natural log

print(f"Total known motifs tested: {len(known)}")
print(f"Significant (p < 1e-5): {(known['pvalue'].astype(float) < 1e-5).sum()}")
print("\nTop 5 enriched motifs:")
print(
    known[["motif_name", "pct_target", "pct_bg", "fold_enrichment", "-log10_pvalue"]]
    .head(5)
    .to_string(index=False)
)
# Total known motifs tested: 391
# Significant (p < 1e-5): 47
# Top 5 enriched motifs:
#  motif_name  pct_target  pct_bg  fold_enrichment  -log10_pvalue
#  CTCF(Zf)/GM12878-CTCF-ChIP-Seq...   78.21    8.32           9.40         312.4
#  ZNF143(Zf)/...                        42.11    5.10           8.26         188.7
#  ...

import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import seaborn as sns

# --- Parse annotated peaks ---
# annotatePeaks.pl output header: PeakID chr start end strand score ...
# Key columns: "Annotation", "Distance to TSS", "Nearest RefSeq", "Gene Name"
annot = pd.read_csv("annotated_peaks.txt", sep="\t", header=0, low_memory=False)
annot.columns = annot.columns.str.strip()

# Simplify annotation categories
def simplify_annotation(ann):
    if pd.isna(ann):
        return "Other"
    ann = str(ann)
    if "promoter" in ann.lower():
        return "Promoter (<2 kb)"
    elif "exon" in ann.lower():
        return "Exon"
    elif "intron" in ann.lower():
        return "Intron"
    elif "tts" in ann.lower() or "downstream" in ann.lower():
        return "Downstream"
    elif "intergenic" in ann.lower():
        return "Intergenic"
    else:
        return "Other"

annot["simple_ann"] = annot["Annotation"].apply(simplify_annotation)
ann_counts = annot["simple_ann"].value_counts()

print("Peak annotation summary:")
print(ann_counts.to_string())
# Peak annotation summary:
# Intron             5832
# Intergenic         3102
# Promoter (<2 kb)   2218
# Exon                891
# Downstream          287
# Other               123

# --- Plot 1: Annotation pie chart ---
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

colors = ["#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00", "#A65628"]
axes[0].pie(
    ann_counts.values,
    labels=ann_counts.index,
    colors=colors[:len(ann_counts)],
    autopct="%1.1f%%",
    startangle=140,
    textprops={"fontsize": 9}
)
axes[0].set_title("Peak Genomic Annotation\n(n=12,453 peaks)", fontsize=11)

# --- Plot 2: Top known motif enrichments ---
known = pd.read_csv("motif_output/knownResults.txt", sep="\t", header=0)
known.columns = [
    "motif_name", "consensus", "pvalue", "log_pvalue",
    "qvalue", "n_target_seqs", "pct_target", "n_bg_seqs", "pct_bg"
]
known["pct_target"] = known["pct_target"].str.rstrip("%").astype(float)
known["pct_bg"] = known["pct_bg"].str.rstrip("%").astype(float)
known["-log10_p"] = (-known["log_pvalue"].astype(float)) / 2.303
known["short_name"] = known["motif_name"].str.split("/").str[0]

top_motifs = known.head(15).copy()
sns.barplot(
    data=top_motifs,
    x="-log10_p",
    y="short_name",
    palette="Blues_r",
    ax=axes[1]
)
axes[1].set_xlabel("-log10(p-value)", fontsize=10)
axes[1].set_ylabel("Motif", fontsize=10)
axes[1].set_title("Top 15 Enriched Known Motifs\n(HOMER/JASPAR database)", fontsize=11)
axes[1].axvline(x=5, color="red", linestyle="--", linewidth=0.8, label="p=1e-5")
axes[1].legend(fontsize=8)

plt.tight_layout()
plt.savefig("homer_summary.png", dpi=150, bbox_inches="tight")
plt.close()
print("Saved: homer_summary.png")