甲基化 + RNA 双组学： Motif 富集分析（基于 HOMER）
模块简介
本模块基于 HOMER (Hypergeometric Optimization of Motif EnRichment) 开发，用于在基因组区域（如差异甲基化区域 DMRs、ChIP-seq peaks 等）中进行 DNA 序列 motif 分析。
HOMER 是一款成熟且广泛使用的 motif 发现工具，能够对给定的基因组区域开展系统性的序列分析，主要实现以下分析目标：
- De novo motif finding（从头 motif 发现）：在目标序列中寻找新的、未知的 DNA 序列模式，不依赖已知的 motif 数据库，适用于发现新的转录因子结合位点或调控序列。
- Known motif enrichment（已知 motif 富集分析）：评估目标序列中已知转录因子 motif 的富集情况，并与 HOMER 内置的 motif 数据库进行比对，用于识别可能参与调控的转录因子。
本分析主要适用于人源和小鼠样本，同时也支持其他具备参考基因组 FASTA 文件的物种。
输入文件准备
本分析需提供一个 BED 格式的输入文件，用于描述待分析的基因组区域（如差异甲基化区域 DMRs、ChIP-seq peaks 等）的坐标信息。
输入文件格式
BED 文件应为制表符分隔的文本文件，且至少包含以下三列信息：
| 列号 | 列名 | 说明 | 示例 |
|---|---|---|---|
| 1 | 染色体名称 | 染色体名称（如 chr1, chr2, chr11 等） | chr2 |
| 2 | 起始位置 | 起始位置（0-based，即从 0 开始计数） | 196161361 |
| 3 | 结束位置 | 结束位置（1-based，即不包含该位置） | 196164361 |
文件说明：
- 坐标系统采用 0-based 编码方式（起始位置从 0 开始）。
- 需明确指定所使用的基因组版本（如 hg38、hg19、mm10、mm9 等）。
- 每个区域代表一个差异甲基化区域 (DMR)。
文件结构示例
chr1 116753470 116770470
chr2 196161361 196164361
chr2 203704361 203712361
...HOMER 分析流程
关键参数设置
在分析过程中，以下参数对分析结果的可靠性及计算效率具有重要影响：
| 参数名称 | 中文释义 | 默认值/建议值 | 详细说明 |
|---|---|---|---|
-len | Motif 长度 | 6,7,8,9,10,11,12 | Motif 发现参数。 指定要分析的 motif 长度，可以指定多个长度（用逗号分隔）。 • 较小长度 (6-8bp)：适合发现较短的转录因子结合位点 • 中等长度 (9-12bp)：适合发现中等长度的调控序列 • 推荐设置： 6,7,8,9,10,11,12（参考 Nature 2024）注意：默认值为 8,10,12。 |
-size | 序列提取长度 | 200 | 序列提取参数。 从每个 DMR 中心提取的序列长度 (bp)。 • 较小值 (100-150bp)：更聚焦于中心区域，适合紧密的转录因子结合位点 • 较大值 (200-300bp)：包含更多上下文序列，适合较分散的调控元件 注意：默认值为 200bp。 |
-mask | 屏蔽重复序列 | 启用 | 质量控制参数。 屏蔽重复序列区域，避免重复序列干扰 motif 发现。 建议：始终使用此选项，除非特别需要分析重复序列区域。 |
-p | 并行线程数 | 8 | 性能参数。 并行计算使用的线程数。 • 建议值：根据服务器配置设置，通常为 4-16 • 注意：线程数过多可能导致内存占用过高 |
最佳实践建议：
关于
-len参数的选择：建议使用6,7,8,9,10,11,12的组合，以覆盖不同长度的 motif 模式；若分析时间有限，可仅使用默认推荐值8,10,12。关于
-size参数的选择：多数情况下使用默认值 200bp 即可。如果 DMR 区域长度较小（< 500bp），可适当减小该参数；若希望捕获更大范围的潜在调控序列，可将该值增大至 300bp。关于环境配置：需确保 HOMER 的
bin目录已加入PATH环境变量（或在命令中使用完整路径），并确保参考基因组FASTA文件与输入BED文件所使用的基因组版本一致。
HOMER 分析执行
HOMER 分析的基本流程包括：参数配置、环境设置、命令构建和执行。
## homer_bin：HOMER findMotifsGenome.pl 脚本的完整路径
homer_bin <- "/jp_envs/envs/common/bin/findMotifsGenome.pl"
## dmr_bed_file：输入 DMR BED 文件路径
dmr_bed_file <- "./DMRs/DMRs_significant_2.bed"
## genome_fa：参考基因组 fasta 文件路径
# 注意：必须与输入 BED 文件使用的参考基因组版本一致（如 hg38）
genome_fa <- "./genome.fa"
## output_dir：结果输出目录
# HOMER 会自动创建此目录
output_dir <- "DMRs_HOMER/"
# 检查 BED 文件中的 region 数量
# 如果超过 500 个 region，则只使用前 500 个
bed_lines <- readLines(dmr_bed_file, warn = FALSE)
# 过滤空行和注释行（以 # 开头的行）
bed_lines <- bed_lines[bed_lines != "" & !grepl("^#", bed_lines)]
num_regions <- length(bed_lines)
# 用于存储实际使用的 BED 文件路径（可能是原文件或临时文件）
actual_bed_file <- dmr_bed_file
temp_bed_file <- NULL
if (num_regions > 500) {
cat("检测到 BED 文件包含", num_regions, "个 region，超过 500 个\n")
cat("将只使用前 500 个 region 进行分析\n")
# 创建临时文件，只包含前 500 行
temp_bed_file <- tempfile(fileext = ".bed")
writeLines(bed_lines[1:500], temp_bed_file)
actual_bed_file <- temp_bed_file
cat("已创建临时 BED 文件:", temp_bed_file, "\n")
} else {
cat("BED 文件包含", num_regions, "个 region，将全部使用\n")
}
# 设置 PATH 环境变量，确保 HOMER 的依赖脚本可以被找到
# HOMER 的 bin 目录包含所有 perl 脚本（bed2pos.pl, homerTools 等）
homer_bin_dir <- dirname(homer_bin)
current_path <- Sys.getenv("PATH")
Sys.setenv(PATH = paste(homer_bin_dir, current_path, sep = ":"))
# 设置 LD_LIBRARY_PATH 环境变量，确保可以找到动态链接库（如 libcrypt.so.2）
# 这是解决 "error while loading shared libraries" 问题的关键
homer_lib_dir <- file.path(dirname(homer_bin_dir), "lib")
current_ld_path <- Sys.getenv("LD_LIBRARY_PATH")
if (current_ld_path == "") {
Sys.setenv(LD_LIBRARY_PATH = homer_lib_dir)
} else {
# 确保 homer_lib_dir 在 LD_LIBRARY_PATH 的最前面，优先查找
Sys.setenv(LD_LIBRARY_PATH = paste(homer_lib_dir, current_ld_path, sep = ":"))
}
cat("LD_LIBRARY_PATH:", Sys.getenv("LD_LIBRARY_PATH"), "\n")
# findMotifsGenome.pl 的基本命令格式：
# findMotifsGenome.pl <DMR bed 文件> <参考基因组> <输出目录> [选项]
args <- c(
actual_bed_file, # DMR BED 文件（可能是原文件或临时文件，只包含前 500 个 region）
genome_fa, # 参考基因组 fasta 文件
output_dir, # 输出目录
"-len", "6,7,8,9,10,11,12", # Motif 长度：分析 6-12bp 的 motif
"-mask", # 屏蔽重复序列
"-p", "8" # 使用 8 个线程并行计算
)
# 在 Jupyter 中使用 system2，捕获输出以便调试
cat("开始执行 HOMER 分析...\n")
cat("HOMER bin 目录:", homer_bin_dir, "\n")
cat("命令:", homer_bin, paste(args, collapse = " "), "\n\n")
flush.console() # 强制刷新输出缓冲区
# 使用 stdout=TRUE 和 stderr=TRUE 捕获输出
result <- system2(
command = homer_bin,
args = args,
wait = TRUE,
stdout = TRUE,
stderr = TRUE
)
# 如果使用了临时 BED 文件，删除它
if (!is.null(temp_bed_file) && file.exists(temp_bed_file)) {
unlink(temp_bed_file)
cat("已清理临时 BED 文件\n")
}将只使用前 500 个 region 进行分析
已创建临时 BED 文件: /tmp/RtmpKd4Kmx/file4274c462b79.bed
LD_LIBRARY_PATH: /jp_envs/envs/common/lib
开始执行 HOMER 分析...n HOMER bin 目录: /jp_envs/envs/common/bin
命令: /jp_envs/envs/common/bin/findMotifsGenome.pl /tmp/RtmpKd4Kmx/file4274c462b79.bed ./genome.fa DMRs_HOMER/ -len 6,7,8,9,10,11,12 -mask -p 8
结果解读
HOMER 会在输出目录中生成多类结果文件，主要分为 De novo motif finding 结果和 Known motif enrichment 结果两大类。
De novo Motif Finding 结果
homerResults.html和homerResults/目录：De novo motif 发现的格式化网页报告，是主要查看文件。它包含每个 motif 的可视化 logo 图、统计信息（p-value、富集倍数等）、序列模式以及详细的motif<length>.motif文件，便于在浏览器中快速浏览和筛选候选 motif。homerMotifs.motifs<length>：按 motif 长度分类的结果文件（如homerMotifs.motifs6对应 6bp 长度）。包含该长度下所有发现 motif 的序列模式及统计信息，用于查看特定长度的分析结果。homerMotifs.all.motifs：所有长度 motif 的合并文件。它汇总了所有homerMotifs.motifs<length>的内容，用于查看本次分析发现的完整 motif 集合。
Known Motif Enrichment 结果
knownResults.txt：已知 motif 富集分析的文本结果（TSV 格式，可用 Excel 打开），是主要查看文件。该文件详细列出了Motif Name（名称）、Consensus（共有序列）、富集统计量（P-value、Log P-value、q-value）以及 Motif 在目标和背景序列中的分布情况（Target/Background Sequences with Motif的数量#和百分比%），用于识别富集的转录因子。knownResults.html和knownResults/目录：已知 motif 富集分析的格式化网页报告 (HTML) 及配套目录。内容包含富集 motif 的可视化 logo、统计信息表格、富集倍数以及详细的known<length>.motif文件，提供直观且交互友好的结果浏览方式。