homer-motif-analysis

HOMER — モチーフ解析とピークアノテーション

概要

HOMER（Hypergeometric Optimization of Motif EnRichment）は、ゲノムの制御エレメントを解析するための Perl/C++ ツール群です。主なコマンドは 2 つあります。findMotifsGenome.pl は de novo モチーフ探索と、JASPAR/HOMER データベースに対する既知モチーフのエンリッチメント解析を実行します。annotatePeaks.pl は各ピークを最も近い遺伝子、TSS までの距離、ゲノム上の特徴クラス（プロモーター、イントロン、遺伝子間領域、リピート）に対応付けます。HOMER は、MACS3 などのピークコーラーが出力する BED 形式のピークファイルと参照ゲノムアセンブリを入力として受け取り、マッチしたバックグラウンドと比較して、エンリッチしたモチーフを p 値とエンリッチメント倍率の順に並べた HTML/テキストレポートを出力します。

使用する場面

• JASPAR または HOMER モチーフライブラリの既知モチーフをエンリッチメント解析し、ChIP-seq のピークセットに結合している転写因子を特定する場合
• 結合する TF の事前知識がない状態で、ATAC-seq のオープンクロマチン領域から de novo で新規配列モチーフを発見する場合
• 一方の条件のピークセットをターゲット、もう一方をバックグラウンドとして HOMER を実行し、2 つの条件（例：処理群と未処理群）のモチーフ構成を比較する場合
• 下流の機能解析や DESeq2 の結果との統合に向けて、ゲノムピークに最も近い遺伝子と TSS までの距離を付与する場合
• ChIP-seq 実験の品質を検証する場合：プルダウンが成功していれば、ターゲット TF の標準的なモチーフが最上位のヒットとして示されるはずです
• HOMER の入力となるピーク BED ファイルを生成するには、まず macs3-peak-calling を使用してください
• HOMER が発見したモチーフを JASPAR ID や追加の TF メタデータと照合するには、jaspar-database を使用してください
• より確率的な ZOOPS/TCM モデルや MEME Suite エコシステムが必要な場合は、MEME-CHIP（ウェブまたはローカル）を使用してください
• de novo 探索を行わず、既知モチーフのスキャンをより高速に実行する代替手段として、AME（MEME Suite の一部）を使用してください

前提条件

• ソフトウェア：HOMER（Perl とコンパイル済みバイナリ）、conda または手動インストール
• ゲノム：HOMER のインストール後、installGenome.pl でゲノム配列をダウンロードする必要があります
• 入力：ピークの BED ファイル（最低限、chr、start、end の各列が必要）；MACS3 で summit を中心にしたピークを生成するのが理想的です
• Python パッケージ（解析・可視化用）：pandas、matplotlib、seaborn

インストール前に確認：ツールが現在の環境（例：pixi / conda 環境内）ですでに利用可能な場合があります。まず command -v findMotifsGenome.pl を実行し、パスが返された場合は以下のインストールコマンドを省略してください。pixi プロジェクト内で実行する場合は、findMotifsGenome.pl をそのまま実行するのではなく、pixi run findMotifsGenome.pl でツールを呼び出してください。

# 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)"

D1 は質問するのではなく、既存の情報から導き出します。ピーク座標は、それが呼び出されたアセンブリ以外では意味を持たないため、回答は上流のステージで使用したアセンブリです。D3 は、十分に検討されないままデフォルトのままにされることが最も多い判断事項です。ピークと GC 含量が一致しないバックグラウンドを使うと、生物学的な特徴ではなく GC 含量がエンリッチメントされます。

クイックスタート

# 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"

ワークフロー

ステップ 1：インストールとゲノムのセットアップ

HOMER をインストールし、モチーフ解析に必要な参照ゲノム配列をダウンロードします。

# 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  ...

ステップ 2：入力ピークファイルの準備

最適なモチーフ解像度を得るため、MACS3 の出力から summit を中心にした BED ファイルを準備します。

# 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

ステップ 3：de novo モチーフ探索

findMotifsGenome.pl を実行して、de novo モチーフ探索と既知モチーフのエンリッチメント解析を同時に行います。

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

ステップ 4：既知モチーフのエンリッチメントのみを実行

ピーク内に特定の既知モチーフが出現するかをスキャンするか、処理を高速化するため de novo 探索を省略します。

# 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

ステップ 5：ピークのアノテーション

annotatePeaks.pl を使用して、各ピークをゲノム上の特徴と最も近い遺伝子に対応付けます。

# 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

ステップ 6：Python で HOMER の結果を解析

knownResults.txt と de novo モチーフファイルを読み込み、pandas で下流解析を行います。

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
#  ...

ステップ 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")

主なパラメーター

パラメーター	デフォルト	範囲／オプション	効果
-size	200	50–1000	モチーフ検索の対象とする、ピークの中点を中心としたウィンドウサイズ（bp）。TF ChIP-seq では 200、ATAC-seq では 150 を使用します
-mask	オフ	フラグ	反復配列要素をマスクします（N マスク）。偽陽性を減らすため、強く推奨します
-p	1	1–64	並列処理する CPU スレッド数。一般的なデータセットでは 4～8 を使用します
-S	25	1–100	検出する de novo モチーフの数。通常は 25 で十分です
-nomotif	オフ	フラグ	de novo モチーフ検出をスキップし、既知モチーフのエンリッチメントのみを実行します（10 倍高速）
-bg	自動生成	BED ファイルのパス	カスタム背景領域ファイル。省略した場合、自動生成される背景は GC 含量を一致させます
-find	—	.motif ファイルのパス	特定のモチーフ行列の出現をピーク内でスキャンし、ピークごとの注釈を出力します
-mis	2	0–4	既知モチーフとのマッチングで許容される最大ミスマッチ数
-len	8,10,12	カンマ区切りの整数	de novo 検索で試すモチーフ長。6 を追加すると、より短いモチーフを検出します
genome	必須	hg38、mm10、hg19、dm6、ce11 など	参照ゲノムアセンブリ。installGenome.pl でインストールしておく必要があります

よく使われるレシピ

レシピ 1：findMotifsGenome.pl 用 Python subprocess ラッパー

Python スクリプトから HOMER を実行し、完了ステータスを取得します。

import subprocess
import sys
from pathlib import Path

def run_homer_motifs(
    peaks_bed: str,
    genome: str,
    output_dir: str,
    size: int = 200,
    mask: bool = True,
    cpus: int = 4,
    n_motifs: int = 25,
    denovo: bool = True,
    bg_bed: str = None
) -> int:
    """Run findMotifsGenome.pl and return exit code."""
    Path(output_dir).mkdir(parents=True, exist_ok=True)
    cmd = [
        "findMotifsGenome.pl", peaks_bed, genome, output_dir,
        "-size", str(size),
        "-p", str(cpus),
        "-S", str(n_motifs),
    ]
    if mask:
        cmd.append("-mask")
    if not denovo:
        cmd.append("-nomotif")
    if bg_bed:
        cmd.extend(["-bg", bg_bed])

    print(f"Running: {' '.join(cmd)}")
    result = subprocess.run(cmd, capture_output=True, text=True)
    if result.returncode != 0:
        print(f"HOMER stderr:\n{result.stderr[-2000:]}", file=sys.stderr)
        return result.returncode
    print(f"Done. Output in {output_dir}/")
    return 0

# Usage
exit_code = run_homer_motifs(
    peaks_bed="peaks/tf_chip_200bp.bed",
    genome="hg38",
    output_dir="motif_output/",
    size=200, mask=True, cpus=4, n_motifs=25
)
print(f"Exit code: {exit_code}")
# Running: findMotifsGenome.pl peaks/tf_chip_200bp.bed hg38 motif_output/ -size 200 -p 4 -S 25 -mask
# Done. Output in motif_output/
# Exit code: 0

レシピ 2：複数のピークセットに対するモチーフ解析のバッチ処理

複数の ChIP-seq サンプルに対して、ループ内で HOMER を実行します。

#!/bin/bash
# Batch motif analysis for several ChIP-seq experiments (same genome)
GENOME="hg38"
PEAK_DIR="peaks"
MOTIF_DIR="motif_results"
mkdir -p "$MOTIF_DIR"

SAMPLES=(CTCF H3K4me3 FOXA2 RUNX1)

for sample in "${SAMPLES[@]}"; do
    echo "=== Processing $sample ==="
    peak_file="${PEAK_DIR}/${sample}_summits.bed"

    if [ ! -f "$peak_file" ]; then
        echo "  Skipping $sample: $peak_file not found"
        continue
    fi

    # Center on summit ±100 bp
    awk 'BEGIN{OFS="\t"} {s=$2-100<0?0:$2-100; print $1,s,$2+100,$4,$5}' \
        "$peak_file" > "${PEAK_DIR}/${sample}_200bp.bed"

    findMotifsGenome.pl "${PEAK_DIR}/${sample}_200bp.bed" "$GENOME" \
        "${MOTIF_DIR}/${sample}/" \
        -size 200 -mask -p 4 -nomotif \
        2> "${MOTIF_DIR}/${sample}.log"

    n=$(wc -l < "${MOTIF_DIR}/${sample}/knownResults.txt")
    echo "  $sample: $n motifs tested. Top hit: $(sed -n '2p' "${MOTIF_DIR}/${sample}/knownResults.txt" | cut -f1)"
done
# === Processing CTCF ===
#   CTCF: 392 motifs tested. Top hit: CTCF(Zf)/GM12878-CTCF-ChIP-Seq(GSE32465)/Homer
# === Processing FOXA2 ===
#   FOXA2: 392 motifs tested. Top hit: Foxa2(Forkhead)/Liver-Foxa2-ChIP-Seq(GSE25694)/Homer

レシピ 3：homerResults/ から de novo モチーフ行列を解析する

後続の比較やプロットに使用するため、de novo モチーフの PWM 行列を読み込みます。

import os
import re
import pandas as pd

def parse_homer_motif(motif_file: str) -> dict:
    """Parse a HOMER .motif file into name, log_odds_threshold, and PWM."""
    with open(motif_file) as f:
        header = f.readline().strip()  # >motif_name\tlog_odds\tlog_p-value\t0\tnucs
        rows = []
        for line in f:
            line = line.strip()
            if line:
                rows.append([float(x) for x in line.split("\t")])
    parts = header.lstrip(">").split("\t")
    name = parts[0]
    log_odds = float(parts[1]) if len(parts) > 1 else 0.0
    log_p = float(parts[2]) if len(parts) > 2 else 0.0
    pwm = pd.DataFrame(rows, columns=["A", "C", "G", "T"])
    return {"name": name, "log_odds": log_odds, "log_p": log_p, "pwm": pwm}

# Load all de novo motifs
motif_dir = "motif_output/homerResults/"
motifs = []
for fn in sorted(os.listdir(motif_dir)):
    if fn.endswith(".motif"):
        motif = parse_homer_motif(os.path.join(motif_dir, fn))
        motifs.append(motif)
        print(f"{fn}: {motif['name']} ({len(motif['pwm'])} positions, log_p={motif['log_p']:.1f})")

print(f"\nLoaded {len(motifs)} de novo motifs")
# motif1.motif: CTCF-motif (19 positions, log_p=-8234.1)
# motif2.motif: CTCFL-motif (17 positions, log_p=-3421.7)
# ...
# Loaded 25 de novo motifs

レシピ 4：ピークに注釈を付け、差次的発現解析の結果と結合する

ピーク周辺の制御対象遺伝子を特定するため、ピークの注釈を RNA-seq の DE 結果と結合します。

import pandas as pd

# Load annotated peaks
annot = pd.read_csv("annotated_peaks.txt", sep="\t", header=0, low_memory=False)
annot.columns = annot.columns.str.strip()

# Key columns from HOMER annotation
peak_genes = annot[["PeakID (cmd=annotatePeaks.pl peaks.bed hg38)",
                     "Chr", "Start", "End",
                     "Annotation", "Distance to TSS",
                     "Nearest RefSeq", "Gene Name"]].copy()
peak_genes.columns = ["peak_id", "chr", "start", "end",
                      "annotation", "tss_dist", "refseq", "gene_name"]

# Filter promoter-proximal peaks (within 2 kb of TSS)
promoter_peaks = peak_genes[peak_genes["tss_dist"].abs() < 2000].copy()
print(f"Promoter-proximal peaks (|TSS| < 2kb): {len(promoter_peaks)}")
# Promoter-proximal peaks (|TSS| < 2kb): 2218

# Load DESeq2 results (gene_name, log2FC, padj)
de_results = pd.read_csv("deseq2_results.csv")
de_sig = de_results[de_results["padj"] < 0.05].copy()

# Merge: find DE genes with a nearby peak
merged = promoter_peaks.merge(de_sig, on="gene_name", how="inner")
print(f"DE genes with promoter-proximal peak: {len(merged)}")
print(merged[["gene_name", "tss_dist", "annotation", "log2FC", "padj"]].head())
# DE genes with promoter-proximal peak: 347

想定される出力

出力	形式	説明
motif_output/knownResults.txt	TSV	検定したすべての既知モチーフ：名前、p 値、q 値、% target、% background。主な結果ファイルです
motif_output/knownResults.html	HTML	モチーフロゴ、統計情報、リンクを含むインタラクティブな HTML レポートです
motif_output/homerResults/	ディレクトリ	de novo モチーフ：motifN.motif（PWM 行列）、motifN.logo.png、similar.motifs.txt
motif_output/homerResults.html	HTML	de novo モチーフのインタラクティブな HTML レポートです
annotated_peaks.txt	TSV	ピークごとに 1 行：chr、start、end、注釈、TSS までの距離、最も近い遺伝子、RefSeq ID
peaks_with_motif1.txt	TSV	ピークごとのモチーフ出現スコア（-find モードの出力）

トラブルシューティング

問題	原因	解決策
ERROR: Genome not found	installGenome.pl でゲノムがインストールされていません	installGenome.pl hg38 を実行し、~/.homer/data/genomes/hg38/ が存在することを確認します
command not found: findMotifsGenome.pl	HOMER のバイナリが $PATH に含まれていません	HOMER の bin ディレクトリを PATH に追加します：export PATH=~/.homer/bin:$PATH（または conda 環境を有効化します）
de novo モチーフの品質が低い（GC 含量の高いアーティファクト）	反復配列がマスクされていません	必ず -mask を使用します。また、-size を 300 に増やすか、150 に減らしてみてください
有意な既知モチーフがない（すべて p > 0.01）	ピークが少なすぎるか、ピークサイズが不適切です	ピークが 500 以上あることを確認し、-size が予想される結合フットプリントに合っているか確認します（TF の場合は 150～200）
conda で HOMER を使用すると競合する（Perl モジュールエラー）	HOMER の Perl スクリプトが conda 環境の Perl と競合しています	専用の conda 環境に HOMER をインストールします：conda create -n homer_only -c bioconda homer
annotatePeaks.pl の結果がすべて "Intergenic" になる	ゲノムアノテーションがインストールされていません	HOMER のアノテーションには GTF をダウンロードする installGenome.pl が必要です。~/.homer/data/genomes/hg38/ に *.ann ファイルがあるか確認します
実行時間が非常に長い（2 時間超）	シングルスレッドで実行しているか、ピークセットが大きすぎます	-p 8 を追加します。速度を上げるには -S を 10 に減らし、既知モチーフのみを解析する場合は -nomotif を使用します
-bg で指定したカスタム背景から質の低いモチーフが得られる	背景領域の GC 含量が大きく異なります	GC 含量を一致させた背景を使用します。bedtools shuffle に -noOverlapping を指定して、同じサイズのランダムな領域を生成します

参考資料

• HOMER 公式ドキュメント — HOMER のすべてのコマンド、パラメーター、出力形式を網羅したガイド
• Heinz S et al. (2010) "Simple Combinations of Lineage-Determining Transcription Factors Prime cis-Regulatory Elements Required for Macrophage and B Cell Identities." Molecular Cell 38(4):576–589. DOI:10.1016/j.molcel.2010.05.004
• HOMER GitHub: samtools/homer — ソースコードと課題トラッカー
• JASPAR 2024 database — HOMER が既知モチーフのエンリッチメント解析に使用する、キュレーション済みの TF 結合モチーフデータベース