mcpField
H
homer-motif-analysis

homer-motif-analysis

github.com/jaechang-hits数据与数据库38467 次安装6 天前更新
代码仓库

关于 homer-motif-analysis

homer-motif-analysis 是使用 HOMER 进行基因组调控元件基序分析和峰注释的技能,适用于 ChIP-seq、ATAC-seq 等数据分析。它可发现新基序、富集已知基序,并将峰关联到附近基因和基因组特征。

主要功能

  1. 01从峰集合中发现新基序
  2. 02分析已知基序富集
  3. 03将峰注释到附近基因
  4. 04计算峰到 TSS 的距离
  5. 05按基因组特征分类峰

适用场景

  1. 01从 ChIP-seq 峰识别候选转录因子
  2. 02在 ATAC-seq 区域寻找新基序
  3. 03比较不同条件下的基序富集
  4. 04为峰关联附近基因以供后续分析

安装

复制下面的命令,在终端中运行。

$npx skills add https://github.com/jaechang-hits/sciagent-skills --skill homer-motif-analysis

安装后在 Claude Code 中输入/homer-motif-analysis确认已加载。

详细介绍

HOMER——基序分析和峰注释

概述

HOMER(基序富集的超几何优化)是一套用于分析基因组调控元件的 Perl/C++ 工具。它的两个主要命令是 findMotifsGenome.pl 和 annotatePeaks.pl:前者利用 JASPAR/HOMER 数据库进行从头基序发现和已知基序富集分析;后者将每个峰映射到最近的基因,并标注其与 TSS 的距离及基因组特征类别(启动子、内含子、基因间区、重复序列)。HOMER 接受来自 MACS3 或类似峰调用工具的 BED 格式峰文件以及参考基因组组装版本作为输入,并输出 HTML/文本报告,按照 p 值和相对于匹配背景的富集倍数对富集基序进行排名。

适用场景

  • 从 JASPAR 或 HOMER 基序库中富集已知基序,以识别 ChIP-seq 峰集中的结合转录因子
  • 在 ATAC-seq 开放染色质区域中从头发现新的序列基序,无需预先了解结合 TF
  • 将一组峰作为目标集、另一组作为背景集运行 HOMER,从而比较两种条件(例如处理组与未处理组)之间的基序分布
  • 为下游功能分析或与 DESeq2 结果整合,注释基因组峰的最近基因及其与 TSS 的距离
  • 验证 ChIP-seq 实验质量:成功的免疫沉淀应使目标 TF 的典型基序成为排名首位的结果
  • 先使用 macs3-peak-calling 生成 BED 格式峰文件,作为 HOMER 的输入
  • 使用 jaspar-database 将 HOMER 发现的基序与 JASPAR ID 及其他 TF 元数据交叉比对
  • 需要概率性更强的 ZOOPS/TCM 模型或 MEME Suite 生态系统时,使用 MEME-CHIP(网页版或本地版)
  • 将 AME(MEME Suite 的一部分)用作已知基序扫描的更快替代方案,无需从头发现

前置条件

  • 软件:HOMER(Perl + 编译后的二进制文件),使用 conda 安装或手动安装
  • 基因组:安装 HOMER 后,必须通过 installGenome.pl 下载基因组序列
  • 输入:峰的 BED 文件(至少包含 chr、start、end 列);理想情况下使用以峰顶为中心的 MACS3 峰
  • Python 包(用于解析/可视化):pandas、matplotlib、seaborn

安装前请检查:当前环境中可能已提供该工具(例如在 pixi / conda 环境中)。先运行 command -v findMotifsGenome.pl,如果返回路径,则跳过下面的安装命令。在 pixi 项目中运行时,请通过 pixi run findMotifsGenome.pl 调用该工具,而不是直接运行 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 输出准备以峰顶为中心的 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:从头发现基序

运行 findMotifsGenome.pl,同时进行从头基序发现和已知基序富集分析。

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:仅进行已知基序富集

扫描峰中的特定已知基序,或跳过从头发现以提高速度。

# 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 和从头发现的基序文件读入 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")

关键参数

参数默认值范围 / 选项作用
-size20050–1000基序搜索时以峰中点为中心的窗口大小(bp);TF ChIP-seq 使用 200,ATAC-seq 使用 150
-mask关闭标志屏蔽重复元件(N-mask);强烈建议使用,以减少假阳性
-p11–64并行 CPU 线程数;典型数据集建议使用 4–8
-S251–100要查找的从头基序数量;通常 25 个就足够
-nomotif关闭标志跳过从头基序发现;仅执行已知基序富集(速度快 10 倍)
-bg自动生成BED 文件路径自定义背景区域文件;若省略,则自动生成的背景区域会与 GC 含量匹配
-find—.motif 文件路径扫描峰中指定基序矩阵的出现情况;输出逐峰注释
-mis20–4匹配已知基序时允许的最大错配数
-len8,10,12逗号分隔的整数从头搜索时尝试的基序长度;添加 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/ 中的从头基序矩阵

加载从头基序 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.txtTSV测试的所有已知基序:名称、p 值、q 值、目标序列百分比、背景序列百分比;主要结果文件
motif_output/knownResults.htmlHTML交互式 HTML 报告,包含基序标识图、统计信息和链接
motif_output/homerResults/目录从头发现的基序:motifN.motif(PWM 矩阵)、motifN.logo.png、similar.motifs.txt
motif_output/homerResults.htmlHTML从头发现基序的交互式 HTML 报告
annotated_peaks.txtTSV每个峰一行:染色体、起始位置、结束位置、注释、与 TSS 的距离、最近基因、RefSeq ID
peaks_with_motif1.txtTSV逐峰基序出现评分(来自 -find 模式)

故障排除

问题原因解决方案
ERROR: Genome not found未通过 installGenome.pl 安装基因组运行 installGenome.pl hg38;确认 ~/.homer/data/genomes/hg38/ 存在
command not found: findMotifsGenome.plHOMER 二进制文件不在 $PATH 中将 HOMER bin 目录添加到 PATH:export PATH=~/.homer/bin:$PATH(或激活 conda 环境)
从头预测基序质量差(出现富含 GC 的伪影)重复元件未屏蔽始终使用 -mask;也可尝试将 -size 增大到 300,或减小到 150
没有显著的已知基序(所有 p 值均 > 0.01)峰数量过少或峰大小不合适确保至少有 500 个峰;检查 -size 是否符合预期的结合足迹(转录因子为 150–200)
Conda 中的 HOMER 冲突(Perl 模块错误)HOMER 的 Perl 脚本与 conda 环境中的 Perl 冲突在专用 conda 环境中安装 HOMER:conda create -n homer_only -c bioconda homer
annotatePeaks.pl 的结果全部为“Intergenic”未安装基因组注释HOMER 注释需要运行 installGenome.pl 下载 GTF;检查 ~/.homer/data/genomes/hg38/ 中是否包含 *.ann 文件
运行时间过长(>2 小时)单线程运行或峰集过大添加 -p 8;将 -S 减少到 10 以提高速度;仅分析已知基序时使用 -nomotif
-bg 自定义背景产生的基序质量差背景区域的 GC 含量差异很大使用 GC 匹配的背景;运行带有 -noOverlapping 的 bedtools shuffle,生成大小相同的随机区域

参考资料

评论

支持 Markdown