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