🔍 核心原则：BAM 文件已在 Dorado 阶段“锁定”修饰类型
您提供的 Dorado 命令非常关键。--modified-bases 参数决定了 BAM 文件中 MM/ML 标签包含哪些修饰信号。Nextflow 管道无法凭空检测未写入 BAM 的修饰，它只能解析、过滤和定量 BAM 里已有的标签。
| Dorado Basecalling 参数 | BAM 内含有的修饰信号 | Nextflow 能调出的结果 |
|---|---|---|
--modified-bases 6mA |
仅 A (6mA) |
只能输出 6mA |
--modified-bases 4mC_5mC |
仅 C (4mC) 和 CG (5mC) |
只能输出 4mC + 5mC |
⚠️ 重要结论：因为您用 Dorado 分两次生成了不同的 BAM，所以必须分两次运行 Nextflow。管道无法在一个 BAM 里“补全”另一种修饰。
🚀 最佳 Nextflow 命令推荐
1️⃣ 处理 6mA BAM（仅含 6mA 信号）
nextflow run nf-core/methylong -r 2.0.0 -profile docker \
--input samplesheet_6mA.csv \
--outdir methylome_out_6mA \
--ont_aligner minimap2 \
--m6a \
--skip_snvs \
-resume \
-work-dir methylome_out_6mA/work2️⃣ 处理 4mC_5mC BAM（含 4mC + 5mC 信号）
nextflow run nf-core/methylong -r 2.0.0 -profile docker \
--input samplesheet_c4mC_5mC.csv \
--outdir methylome_out_4mC_5mC \
--ont_aligner minimap2 \
--all_contexts \
--skip_snvs \
-resume \
-work-dir methylome_out_4mC_5mC/work🧩 为什么这样配置是最佳实践？
| 参数 | 在您当前场景中的真实作用 | 为什么必须这样设 |
|---|---|---|
--m6a |
强制 modkit 启用 --motif A 0 提取逻辑 |
Dorado 的 6mA BAM 里只有 A 标签。若不写此参数，modkit 默认只输出 CG，会导致 An6.bed.gz 为空或仅含背景噪声。 |
--all_contexts |
告诉 modkit：不要过滤，输出 BAM 中所有存在的 motif | Dorado 的 4mC_5mC 模型会同时写入 C (4mC) 和 CG (5mC) 标签。使用 --all_contexts 可一次性完整保留两者，无需手动拆分参数。 |
--skip_snvs |
跳过 SNP calling + WhatsHap phasing | 细菌为单倍体，无杂合位点。跳过此步骤可节省 30%~50% 运行时间，并彻底避免后续 DSS_HAPLOTYPE_LEVEL 因非回文 motif 报错。 |
独立 --outdir & -work-dir |
隔离两次运行的缓存与结果 | 防止 Nextflow 因输入路径不同但目录相同而混淆中间文件，也避免结果被覆盖。 |
🛑 关键避坑指南
❌ 常见误区
“我加了
--all_contexts，是不是就能在一个命令里同时跑出 6mA + 4mC + 5mC？”
✅ 事实
不能。 修饰类型在 Dorado basecalling 阶段已物理固化。Nextflow 的 pileup 步骤只是“读取器”：
- 如果 BAM 里没有
A的MM/ML标签，--m6a或--all_contexts都不会凭空变出 6mA。 - 如果 BAM 里只有
4mC_5mC标签，加了--m6a也只会输出空的A行（或忽略）。
💡 如果未来想“一次运行检测所有修饰”？
必须在 Dorado 阶段生成全修饰 BAM：
# Dorado 一步生成包含 6mA + 4mC + 5mC 的 BAM
dorado basecaller \
--modified-bases 6mA_4mC_5mC \ # 或根据 Dorado 版本使用对应模型
--emit-moves --device cuda:0 \
--reference "${REF}" \
dna_r10.4.1_e8.2_400bps_sup@v5.0.0 \
"${POD5}" | samtools view -b - > "An6_all_mods_mapped.mod.bam"然后 Nextflow 只需：
nextflow run nf-core/methylong ... --input samplesheet_all.csv --all_contexts --skip_snvs ...✅ 运行后验证命令（强烈建议执行）
# 1. 验证 6mA 结果是否成功提取
echo "=== 6mA 结果 ==="
zcat methylome_out_6mA/results/methylation/An6.bed.gz | awk '$4=="A"' | wc -l
zcat methylome_out_6mA/results/methylation/An6.bed.gz | head -3
# 2. 验证 4mC_5mC 结果是否同时包含 CG 和 C
echo "=== 4mC_5mC 结果 ==="
zcat methylome_out_4mC_5mC/results/methylation/An6.bed.gz | cut -f4 | sort | uniq -c
# 预期输出示例: 12345 CG (5mC)
# 67890 C (4mC)📊 总结对照表
| 目标修饰 | Dorado --modified-bases |
Nextflow 关键参数 | 输出 BED 第4列内容 |
|---|---|---|---|
| 6mA | 6mA |
--m6a |
A |
| 4mC + 5mC | 4mC_5mC |
--all_contexts |
CG (5mC) + C (4mC) |
| 全修饰 (未来) | 6mA_4mC_5mC |
--all_contexts |
A + CG + C + 其他 |
🎯 最终建议：您当前分两次生成 BAM 的做法非常规范。只需按上述两条命令分别运行，即可获得干净、无冲突、完整覆盖目标修饰的甲基化图谱。
如需后续 4mC/5mC 分离提取、motif 富集分析或与 6mA 结果整合对比的代码，随时告诉我！🚀
📊 比较：--all_contexts vs --m6a 用于不同甲基化类型
🔍 选项概述
| 选项 | 类型 | 对 modkit pileup 的影响 |
默认 Motif | 额外 Motif |
|---|---|---|---|---|
| 无选项 (默认) | – | 仅 CG 上下文 | CG (5mC/5hmC) |
– |
--m6a |
布尔值 | 添加腺嘌呤 Motif | CG + A |
A 位点的 6mA |
--all_contexts |
布尔值 | BAM 中所有存在的上下文 | 所有检测到的 | 5mC, 4mC, 6mA, 5hmC 等 |
🎯 根据甲基化类型推荐选项
| 甲基化类型 | 化学名称 | 常见于 | 推荐选项 | 理由 |
|---|---|---|---|---|
| 5mC | 5-甲基胞嘧啶 | 真核生物 (CpG) | --all_contexts 或默认 |
默认已捕获 CG；--all_contexts 用于 CHG/CHH |
| 5hmC | 5-羟甲基胞嘧啶 | 真核生物 (脑组织) | --all_contexts |
常与 5mC 在 CG 上下文中共同检测 |
| 6mA | N6-甲基腺嘌呤 | 细菌、藻类、低等真核生物 | ✅ --m6a |
需要明确的 A Motif |
| 4mC | N4-甲基胞嘧啶 | 细菌 | ✅ --all_contexts |
自动捕获非 CG 的 C Motif |
| 组合 (如 4mC + 5mC) | 多种修饰 | 具有复杂限制修饰系统的细菌 | ✅ --all_contexts |
一条命令捕获所有存在的修饰 |
🧪 具体应用场景
场景 1：仅 6mA (您当前的情况)
nextflow run nf-core/methylong \
--input samplesheet_6mA.csv \
--outdir methylome_out_6mA \
--m6a \ # ⭐ 6mA 专用
--skip_snvs \
-profile docker结果： An6.bed.gz 包含 motif = "A" 的行
场景 2：4mC + 5mC 组合 (您的第二个情况)
nextflow run nf-core/methylong \
--input samplesheet_c4mC_5mC.csv \
--outdir methylome_out_4mC_5mC \
--all_contexts \ # ⭐ 捕获 BAM 中 ALL 修饰
--skip_snvs \
-profile docker结果： An6.bed.gz 包含以下行：
motif = "CG"(5mC)motif = "C"或"CCWGG"(4mC)- 其他可能的细菌 Motif
场景 3：仅 5mC/5hmC (真核生物)
nextflow run nf-core/methylong \
--input samplesheet_5mC.csv \
--outdir methylome_out_5mC \
# 无需特殊选项 (默认 = CG) \
--skip_snvs \ # 用于单倍体生物
-profile docker结果： An6.bed.gz 主要包含 motif = "CG"
📋 比较表：输出内容
| 设置选项 | .bed.gz 中包含的 Motif |
第4列 (motif) 示例 |
|---|---|---|
| 默认 (无选项) | 仅 CG 上下文 | CG |
--m6a |
CG + 腺嘌呤 | CG, A |
--all_contexts |
BAM 中所有 | CG, C, CCWGG, A, GATC 等 |
⚠️ 重要提示
1️⃣ --m6a 不排除其他上下文
# 使用 --m6a 您会获得 BOTH (两者):
--m6a → CG (默认) + A (6mA)2️⃣ --all_contexts 更全面
# 使用 --all_contexts 您会获得 ALL (全部):
--all_contexts → CG + C + A + CHG + CHH + ...3️⃣ 前提条件：Basecalling 需包含正确的修饰
管道只能检测 modBAM 中已有的内容！
| Dorado Basecalling 命令 | BAM 中可用的修饰 | 管道选项 |
|---|---|---|
--modified-bases 5mC_5hmC |
仅 CG | 默认足够 |
--modified-bases 6mA |
仅 A | 必需 --m6a |
--modified-bases 5mC_5hmC 6mA |
CG + A | --m6a 或 --all_contexts |
--modified-bases 5mC_5hmC 6mA 4mC |
CG + A + C | ✅ 推荐 --all_contexts |
🎯 决策流程图
您期望哪些修饰？
│
├── 仅 5mC/5hmC (CG)
│ └─→ 默认 (无选项)
│
├── 仅 6mA (A)
│ └─→ --m6a
│
├── 6mA + 5mC 组合
│ └─→ --m6a (简单) 或 --all_contexts (全面)
│
└── 4mC + 5mC + 可能其他
└─→ ✅ --all_contexts (您情况的最佳选择！)💡 您的两个 Samplesheet 总结
| Samplesheet | 包含的修饰 | 推荐选项 | 命令 |
|---|---|---|---|
samplesheet_6mA.csv |
仅 6mA | --m6a |
已执行 ✅ |
samplesheet_c4mC_5mC.csv |
4mC + 5mC | ✅ --all_contexts |
见场景 2 |
🎯 核心区别对比
| 特性 | --m6a |
--all_contexts |
|---|---|---|
| 主要用途 | 专门检测 6mA | 检测所有存在的修饰 |
| 输出 Motif | CG + A | CG + C + A + 其他所有 |
| 适用场景 | 仅关注 6mA | 多种修饰共存 (如 4mC+5mC) |
| 文件大小 | 较小 (仅 CG+A) | 较大 (所有上下文) |
| 处理速度 | 较快 | 稍慢 (需处理更多位点) |
| 灵活性 | 针对性强 | 最全面，无需预先知道修饰类型 |
💡 一句话总结
🎯 对于您的
4mC_5mC数据，--all_contexts是最安全和最全面的选择，因为它会自动捕获 BAM 中所有存在的修饰，无需您单独指定每个 Motif！🎯 对于您的
6mA数据，--m6a是明确且高效的选择，专门针对腺嘌呤修饰优化！
如有其他问题，欢迎继续讨论！🚀
✅ 是的，您必须为 4mC_5mC 样本运行一条独立的命令。
📝 推荐命令
nextflow run nf-core/methylong -r 2.0.0 -profile docker \
--input samplesheet_c4mC_5mC.csv \
--outdir methylome_out_4mC_5mC \
--ont_aligner minimap2 \
--all_contexts \ # ⭐ 关键：同时捕获 5mC 和 4mC
--skip_snvs \
-resume \
-work-dir methylome_out_4mC_5mC/work🔍 为什么必须分开运行？
| 原因 | 说明 |
|---|---|
| 1. 修饰类型与 Motif 不同 | --m6a 会让底层 modkit 强制添加 --motif A 0（腺嘌呤修饰）。而 4mC/5mC 基于胞嘧啶修饰，默认使用 --motif CG 0（5mC）和 --motif C 0（4mC）。同一流程无法同时动态切换 motif 组合。 |
| 2. 输出文件会互相覆盖 | 管道默认将结果写入 --outdir/results/methylation/An6.bed.gz 等固定路径。若共用 --outdir，第二次运行会覆盖第一次的 6mA 结果。 |
| 3. 管道设计逻辑 | nf-core/methylong 按单次运行 = 单一修饰上下文设计。不同修饰需独立堆叠（pileup）与质控流程。 |
⚙️ 参数调整详解
| 参数 | 6mA 命令 | 4mC/5mC 命令 | 原因 |
|---|---|---|---|
--input |
samplesheet_6mA.csv |
samplesheet_c4mC_5mC.csv |
指向不同的 modBAM 文件 |
--outdir |
methylome_out_6mA |
methylome_out_4mC_5mC |
隔离结果，防止覆盖 |
--m6a |
✅ 添加 | ❌ 移除 | 6mA 专属开关，不适用于胞嘧啶修饰 |
--all_contexts |
❌ 可选 | ✅ 添加 | 告诉 modkit 输出 BAM 中所有检测到的修饰（含 5mC 和 4mC） |
--skip_snvs |
✅ 添加 | ✅ 添加 | 细菌单倍体，无需 SNP/Phasing，两次运行都建议保留 |
-work-dir |
.../6mA/work |
.../4mC_5mC/work |
建议分开，避免缓存哈希冲突或清理困难 |
💡 关于 --all_contexts 的重要说明
- 默认行为：不指定该参数时，管道仅输出
CG上下文（即 5mC/5hmC）。 - 添加
--all_contexts后：modkit会解析 BAM 中的MM/ML标签，自动输出所有修饰碱基（包括5mC、4mC、甚至6mA如果存在）。 -
结果文件结构：
methylome_out_4mC_5mC/results/methylation/ ├── An6.bed.gz # 包含所有 motif 行，第4列显示: CG, C, CCWGG 等 └── BG5.bed.gz您可通过以下命令快速分离 4mC 和 5mC：
# 提取 5mC (CG 位点) zcat methylome_out_4mC_5mC/results/methylation/An6.bed.gz | \ awk '$4=="CG"' > An6_5mC.bed # 提取 4mC (非CG的C位点，通常第4列为 "C") zcat methylome_out_4mC_5mC/results/methylation/An6.bed.gz | \ awk '$4=="C"' > An6_4mC.bed
🛡️ 运行注意事项
-resume会智能复用缓存：由于输入 BAM 路径不同，Nextflow 会自动重新计算4mC/5mC流程，不会与6mA缓存冲突。- 磁盘空间：两次运行会产生独立的
work/目录，总占用约是单次的 1.5~2 倍。运行完成后可安全删除work/：rm -rf methylome_out_6mA/work methylome_out_4mC_5mC/work - 验证输出：运行结束后，确认两个输出目录均包含：
ls methylome_out_6mA/results/methylation/*.bed.gz ls methylome_out_4mC_5mC/results/methylation/*.bed.gz
📌 一句话总结
不同修饰类型 = 不同 motif = 必须独立运行 + 独立输出目录。
用--all_contexts替代--m6a，切换 samplesheet 和 outdir，即可安全获取 4mC/5mC 结果。
如需帮助解析 An6.bed.gz 中的 motif 分布或做下游 4mC/5mC 基因注释，随时告知！🚀
📁 单样本无重复实验：nf-core/methylong 输出文件详解 (中文)
🎯 您的实验设计总结
✅ 单样本 (Single sample): An6
✅ 细菌基因组 (Haploid): 无父源/母源染色体
✅ 无生物学重复 (No replicates): 无法做统计差异分析
✅ 无组间比较 (No comparisons): 不计算DMR
✅ 目标: 6mA甲基化谱描述性分析📌 核心结论：您将获得该样本的6mA甲基化全景图谱，用于位点鉴定、motif分析、基因组分布等描述性研究。
📂 预期输出文件结构
methylome_out_6mA/
│
├── 📊 results/ # 🎯 主要结果目录 (您最关心的!)
│ │
│ ├── 🔬 methylation/ # ⭐ 6mA甲基化核心结果
│ │ ├── An6.bed.gz # 🏆 主文件: 6mA位点及修饰概率 (BED格式)
│ │ ├── An6.bed.gz.tbi # BED索引文件 (用于快速检索)
│ │ ├── An6.bigWig # 🎨 可视化文件 (用于IGV/JBrowse)
│ │ └── An6.bedGraph.gz # 备选格式 (bedGraph)
│ │
│ ├── 🗂️ alignment/ # 比对结果
│ │ ├── An6.bam # 比对后的BAM文件 (可能较大)
│ │ ├── An6.bam.bai # BAM索引
│ │ ├── An6.flagstat # 📈 比对统计 (总reads, 比对率等)
│ │ └── An6.stats # 详细比对指标
│ │
│ ├── 🔀 preprocessing/ # 预处理中间结果
│ │ ├── An6.repaired.bam # modkit repair后的BAM
│ │ └── An6.trimmed.fastq.gz # 去接头后的FASTQ (如有trimming)
│ │
│ ├── 🧬 variants/ # ⚠️ 因--skip_snvs可能为空或极简
│ │ └── (可能无文件或仅占位)
│ │
│ └── 📋 qc/ # 质量控制报告
│ ├── multiqc_report.html # 🎯 综合QC报告 (用浏览器打开)
│ ├── fastqc/ # 原始FastQC结果
│ └── modkit_stats/ # modkit运行统计
│
├── 🔧 work/ # Nextflow工作目录 (中间文件, 可删除)
│ └── (大量临时文件, 调试时用)
│
├── 📝 pipeline_info/ # 管道执行信息
│ ├── execution_report.html # Nextflow执行报告
│ ├── execution_timeline.html # 时间线可视化
│ ├── execution_trace.txt # 各步骤资源消耗
│ └── software_versions.yml # 所有工具版本记录 (用于复现)
│
└── 📄 *.log # 日志文件
├── .nextflow.log # 主日志 (排查错误时用)
└── nextflow.log # 备用日志🏆 核心文件详解：您最需要的3个文件
1️⃣ An6.bed.gz — 6mA甲基化位点主文件 ⭐⭐⭐
# 查看文件内容 (前5行)
zcat methylome_out_6mA/results/methylation/An6.bed.gz | head -5📋 BED格式列说明 (6mA专用):
| 列号 | 列名 | 示例值 | 含义 |
|---|---|---|---|
| 1 | chrom | chr1 |
染色体/contig名称 |
| 2 | start | 12345 |
起始位置 (0-based) |
| 3 | end | 12346 |
结束位置 (通常为单碱基) |
| 4 | motif | A |
修饰motif类型 (6mA = “A”) |
| 5 | methylation_prob | 0.87 |
⭐ 修饰概率 (0-1, 越高越可信) |
| 6 | strand | + |
链方向 (+/-) |
| 7 | coverage | 45 |
该位点总覆盖深度 |
| 8 | mod_count | 39 |
检测到修饰的reads数 |
| 9 | mod_prob_low | 0.82 |
修饰概率95%置信区间下限 |
| 10 | mod_prob_high | 0.92 |
修饰概率95%置信区间上限 |
| 11+ | … | … | 其他辅助信息 |
🔍 常用分析命令：
# ✅ 统计高置信度6mA位点 (概率≥0.8, 覆盖度≥10)
zcat methylome_out_6mA/results/methylation/An6.bed.gz | \
awk '$4=="A" && $5>=0.8 && $7>=10' | wc -l
# ✅ 提取所有6mA位点坐标 (用于后续motif分析)
zcat methylome_out_6mA/results/methylation/An6.bed.gz | \
awk '$4=="A" && $5>=0.8 {print $1"\t"$2"\t"$3"\t"$6}' > high_conf_6mA.bed
# ✅ 按染色体统计6mA密度
zcat methylome_out_6mA/results/methylation/An6.bed.gz | \
awk '$4=="A" && $5>=0.8 {count[$1]++} END {for(c in count) print c, count[c]}' | sort2️⃣ An6.bigWig — 可视化文件 ⭐⭐
# 用IGV桌面版打开:
# 1. 加载参考基因组: File → Load Genome from File → genome.fa
# 2. 加载甲基化轨道: File → Load from File → An6.bigWig
# 3. 调整轨道显示: 右键轨道 → Change Track Settings → 设置颜色/高度🎨 可视化效果示例：
基因组位置: ████████████████████
12000 13000
6mA信号: ▁▂▃▅▇█▇▅▃▂▁▁▂▃▅▆▇███▇▆▅▃▂
↑ 高峰 = 高置信度6mA位点🔧 如需自定义bigWig (可选):
# 如果管道未生成bigWig, 可手动转换:
bedGraphToBigWig \
methylome_out_6mA/results/methylation/An6.bedGraph \
genome.fa.sizes \
methylome_out_6mA/results/methylation/An6_custom.bw3️⃣ multiqc_report.html — 综合质量报告 ⭐
# 用浏览器打开 (本地或服务器):
firefox methylome_out_6mA/results/qc/multiqc_report.html
# 或 (无图形界面服务器):
python -m http.server 8000 # 然后通过本地浏览器访问 http://服务器IP:8000📊 报告包含的关键信息：
| 模块 | 内容 | 您需要关注的指标 |
|---|---|---|
| 🔹 FastQC | 原始reads质量 | Per base sequence quality > Q20 ✅ |
| 🔹 Porechop | 接头去除效率 | Adapter trimmed % < 10% ✅ |
| 🔹 ModKit repair | 修饰标签修复 | Successfully repaired % > 95% ✅ |
| 🔹 Minimap2 | 比对统计 | Mapping rate > 90% ✅ |
| 🔹 ModKit pileup | 甲基化调用 | Valid sites called, motif distribution |
| 🔹 General stats | 汇总 | Total reads, coverage, 6mA sites count |
🗂️ 其他辅助文件说明
🔹 An6.flagstat — 比对统计摘要
cat methylome_out_6mA/results/alignment/An6.flagstat# 示例输出:
1000000 + 0 in total (QC-passed reads + QC-failed reads)
950000 + 0 primary
950000 + 0 mapped (95.00% : N/A)
# → 95%比对率，质量良好 ✅🔹 software_versions.yml — 可复现性保障
# 记录所有工具版本，便于论文方法部分或复现:
dorado: 2.0.0+20e87c8b
minimap2: 2.26-r1175
modkit: 0.5.0
samtools: 1.18
# ...🔹 execution_trace.txt — 资源消耗监控
# 用于评估是否需要更多计算资源:
task_name, duration, memory, cpu
ONT_MINIMAP2_ALIGN, 45min, 8GB, 8 cores
MODKIT_PILEUP, 120min, 16GB, 12 cores
# → 如果某步骤内存不足, 下次可增加--max_memory参数❌ 您不会看到哪些文件？(因实验设计)
| 文件类型 | 原因 | 是否正常 |
|---|---|---|
*.vcf.gz |
--skip_snvs 跳过变异检测 |
✅ 正常 |
*_haplotagged.bam |
无phasing步骤 | ✅ 正常 |
*_DMR*.bed |
单样本无法计算差异 | ✅ 正常 |
*_population*.tsv |
无多样本/组间比较 | ✅ 正常 |
*_phasing*.stats |
无WhatsHap执行 | ✅ 正常 |
🎯 这些”缺失”的文件恰恰说明管道正确识别了您的实验设计！
🧪 快速验证：您的结果是否完整可用？
# ✅ 1. 检查核心甲基化文件是否存在且非空
ls -lh methylome_out_6mA/results/methylation/An6.bed.gz
# 应看到: -rw-r--r-- ... An6.bed.gz (大小 > 0)
# ✅ 2. 验证6mA位点数量合理 (细菌基因组通常10^3-10^5个6mA位点)
zcat methylome_out_6mA/results/methylation/An6.bed.gz | \
awk '$4=="A"' | wc -l
# 预期: 1,000 ~ 500,000 (取决于菌株和测序深度)
# ✅ 3. 检查修饰概率分布 (应呈双峰: 低概率背景 + 高概率真实位点)
zcat methylome_out_6mA/results/methylation/An6.bed.gz | \
awk '$4=="A" {print $5}' | \
datamash mean 1 sstdev 1 min 1 max 1 2>/dev/null || \
awk '$4=="A" {sum+=$5; sum2+=$5*$5; n++; if(n==1||$5<min)min=$5; if($5>max)max=$5} END {print "Mean:",sum/n, "Std:",sqrt((sum2-sum*sum/n)/n), "Min:",min, "Max:",max}'
# 预期: Mean ~0.3-0.6, Max >0.9 (有高质量位点)
# ✅ 4. 打开MultiQC确认无严重警告
grep -i "error\|fail" methylome_out_6mA/results/qc/multiqc_report.html | head -3
# 应无关键错误🎯 下一步分析建议 (单样本6mA)
🔬 描述性分析流程：
graph LR
A[An6.bed.gz] --> B[6mA位点过滤]
A --> C[基因组分布分析]
A --> D[Motif富集分析]
B --> B1[高置信度位点: prob≥0.8]
B --> B2[覆盖度过滤: depth≥10]
C --> C1[基因区域注释: 启动子/基因体/间隔区]
C --> C2[与已知功能元件重叠]
D --> D1[提取6mA位点±20bp序列]
D --> D2[MEME/HOMER找motif]
D --> D3[验证是否为已知限制修饰系统motif]💻 实用代码片段：
# 🔹 1. 提取高置信6mA位点
zcat methylome_out_6mA/results/methylation/An6.bed.gz | \
awk '$4=="A" && $5>=0.8 && $7>=10' > An6_6mA_highConf.bed
# 🔹 2. 关联基因注释 (需有genome.gff)
bedtools intersect -a An6_6mA_highConf.bed -b genome.gff -wa -wb > An6_6mA_with_genes.tsv
# 🔹 3. 提取motif序列 (用于MEME分析)
bedtools getfasta -fi genome.fa -bed An6_6mA_highConf.bed -fo An6_6mA_seqs.fa
# 🔹 4. 用MEME找富集motif (需安装MEME Suite)
meme An6_6mA_seqs.fa -oc meme_out -dna -mod zoops -nmotifs 3 -minw 4 -maxw 12
# 🔹 5. 统计6mA在基因区域的分布
awk '{print $NF}' An6_6mA_with_genes.tsv | cut -f3 | sort | uniq -c | sort -rn
# 输出示例: 1250 CDS, 340 promoter, 89 intergenic📊 可视化建议：
| 工具 | 用途 | 命令/链接 |
|---|---|---|
| IGV | 基因组浏览器查看6mA峰 | 加载 An6.bigWig + genome.fa |
| pyGenomeTracks | 出版级多轨道图 | pip install pyGenomeTracks |
| R/ggplot2 | 6mA密度分布图 | geom_density(aes(x=position)) |
| WebLogo | Motif序列标识图 | https://weblogo.berkeley.edu |
💡 一句话总结
🎯 单样本细菌6mA分析的核心产出 =
An6.bed.gz(位点列表) +An6.bigWig(可视化) +multiqc_report.html(质控)
这三个文件已足够开展：位点鉴定、motif发现、基因组分布、功能关联等描述性研究！🇩🇪 Für Einzelproben-Analyse von 6mA in Bakterien sind drei Dateien ausreichend:
An6.bed.gz(Modifikationsstellen),An6.bigWig(Visualisierung),multiqc_report.html(Qualitätskontrolle).
如有具体下游分析需求 (如motif发现、基因关联)，欢迎继续讨论！🚀
📋 nf-core/methylong DMR 选项详解 (中文)
🔍 什么是 DMR？
DMR = Differentially Methylated Region (差异甲基化区域)
指在不同样本组/条件之间，甲基化水平存在统计学显著差异的基因组区域。
示例：
样本组A (野生型): ████████░░░░░░░░ (甲基化比例 ~60%)
样本组B (突变型): ░░░░░░░░████████ (甲基化比例 ~20%)
↓
这是一个 DMR! 🎯🧩 DMR 选项逐项解释
1️⃣ --skip_snvs [boolean]
| 属性 | 说明 |
|---|---|
| 类型 | 布尔值 (开关参数，无需赋值) |
| 作用 | 跳过 SNP/变异检测 (SNV calling) 和 单倍型分型 (phasing) 步骤 |
| 默认 | false (即：默认会执行SNV+phasing) |
| 使用方式 | 添加该参数即启用：--skip_snvs |
🔬 技术流程影响：
默认流程 (不加 --skip_snvs):
原始BAM → Clair3(SNV calling) → WhatsHap(phasing) → haplotagged BAM → haplotype-level DMR
添加 --skip_snvs 后:
原始BAM → 直接进行 reference-level pileup → 跳过 phasing 相关步骤🎯 适用场景：
| 场景 | 是否添加 --skip_snvs |
原因 |
|---|---|---|
| 🦠 细菌/单倍体基因组 | ✅ 推荐添加 | 无杂合位点，phasing无生物学意义 |
| 🔬 仅需参考水平甲基化谱 | ✅ 添加 | 简化流程，加快运行 |
| 👥 人类/二倍体 + 研究等位基因特异性甲基化 | ❌ 不添加 | 需要单倍型信息 |
| 💻 计算资源有限 | ✅ 添加 | SNV+phasing 耗时较长 |
💡 细菌6mA分析示例：
# 细菌是单倍体，直接跳过无意义的SNV/phasing步骤
nextflow run nf-core/methylong \
--input samplesheet.csv \
--outdir results \
--m6a \
--skip_snvs \ # ⭐ 关键参数
-profile docker2️⃣ --haplotype_dmrer [string]
| 属性 | 说明 |
|---|---|
| 类型 | 字符串 (需赋值) |
| 作用 | 选择单倍型水平DMR分析使用的统计工具 |
| 可选值 | dss (默认) 或 modkit |
| 使用方式 | --haplotype_dmrer dss 或 --haplotype_dmrer modkit |
🔧 两种工具对比：
| 工具 | 原理 | 优点 | 缺点 | 适用场景 |
|---|---|---|---|---|
| dss (default) | 基于贝叶斯分层模型的统计检验 | ✅ 成熟稳定，广泛使用 ✅ 支持小样本 ✅ 输出详细统计量 |
❌ 仅支持CpG/CG上下文 ❌ 对非回文motif(如6mA)支持有限 |
5mC/5hmC分析，真核生物 |
| modkit | 基于modkit pileup结果的差异比较 | ✅ 支持多种修饰类型 (5mC/6mA/4mC等) ✅ 与ONT数据原生兼容 |
❌ 较新，文档较少 ❌ 统计方法相对简单 |
6mA/4mC等原核修饰，或需要统一工具链 |
⚠️ 重要注意事项：
# 如果您分析的是6mA (motif "A"，非回文):
# 即使使用 --haplotype_dmrer modkit，
# 底层modkit pileup仍可能使用 --combine-strands 导致报错!
# 解决方案：配合 --skip_snvs 使用，完全跳过haplotype-level DMR
nextflow run nf-core/methylong \
--m6a \
--skip_snvs \ # ⭐ 跳过整个haplotype流程，避免motif冲突
--haplotype_dmrer dss # 即使写了也不会执行，但语法上需要3️⃣ --dmr_population_scale [boolean]
| 属性 | 说明 |
|---|---|
| 类型 | 布尔值 (开关参数) |
| 作用 | 启用群体水平的DMR分析 (多样本/多组比较) |
| 默认 | false (即：默认只做单样本或两两比较) |
| 使用方式 | 添加该参数即启用：--dmr_population_scale |
🔬 分析层级对比：
📊 单样本水平 (默认):
输出: 每个样本的甲基化谱 (bed文件)
用途: 描述性分析，位点鉴定
📊 两两比较 (基础DMR):
输入: 2个样本/组
输出: 差异甲基化区域列表
用途: 简单对照实验 (如: 野生型 vs 突变型)
📊 群体水平 (--dmr_population_scale):
输入: 多个样本，按组分类 (如: 组A有5个样本，组B有6个样本)
输出: 考虑组内变异的统计学显著DMR
用途: 复杂实验设计，需要统计功效的场景📋 群体水平分析所需参数：
启用 --dmr_population_scale 后，必须同时指定：
--dmr_a: 对照组/第一组的样本列表或标识--dmr_b: 实验组/第二组的样本列表或标识
4️⃣ --population_dmrer [string]
| 属性 | 说明 |
|---|---|
| 类型 | 字符串 (需赋值) |
| 作用 | 选择群体水平DMR分析使用的统计工具 |
| 可选值 | dss (默认) 或 modkit |
| 使用方式 | --population_dmrer dss 或 --population_dmrer modkit |
🔧 与 --haplotype_dmrer 的区别：
| 参数 | 分析层级 | 输入数据 | 统计模型 |
|---|---|---|---|
--haplotype_dmrer |
单倍型水平 (H1 vs H2) | 单个样本的haplotagged BAM | 单样本内等位基因比较 |
--population_dmrer |
群体水平 (组A vs 组B) | 多样本的pileup结果 | 多样本间的组间比较 |
💡 实际使用建议：
# 大多数情况下，两个参数设为相同工具即可:
--haplotype_dmrer dss --population_dmrer dss
# 如果只做参考水平+群体比较 (细菌常见场景):
--skip_snvs --dmr_population_scale --population_dmrer modkit5️⃣ --dmr_a & 6️⃣ --dmr_b [string]
| 属性 | 说明 |
|---|---|
| 类型 | 字符串 (需赋值) |
| 作用 | 定义群体水平DMR分析中的两个比较组 |
| 依赖 | 仅在 --dmr_population_scale 启用时生效 |
| 格式 | 通常为样本名称列表，用逗号分隔或引用文件 |
📋 格式示例：
方式1: 逗号分隔的样本名
--dmr_a "sample1,sample2,sample3" \
--dmr_b "sample4,sample5,sample6"方式2: 引用samplesheet中的group列
# samplesheet.csv
group,sample,path,ref,method
WT,sample1,...,...,ont
WT,sample2,...,...,ont
Mut,sample3,...,...,ont
Mut,sample4,...,...,ont--dmr_a "WT" \
--dmr_b "Mut"方式3: 外部文件列表 (如管道支持)
--dmr_a "group_A_samples.txt" \
--dmr_b "group_B_samples.txt"🎯 细菌多菌株比较示例：
# 比较野生型菌株 (3个重复) vs 突变型菌株 (3个重复)
nextflow run nf-core/methylong \
--input samplesheet_6mA.csv \
--outdir results_6mA \
--m6a \
--skip_snvs \
--dmr_population_scale \
--population_dmrer modkit \
--dmr_a "WT_rep1,WT_rep2,WT_rep3" \
--dmr_b "Mut_rep1,Mut_rep2,Mut_rep3" \
-profile docker🧭 决策流程图：如何选择DMR选项？
graph TD
A[开始: 确定分析目标] --> B{研究对象?}
B -->|🦠 细菌/单倍体 | C[添加 --skip_snvs]
B -->|👥 人类/二倍体 | D{需要等位基因特异性分析?}
D -->|是 | E[保留默认: 执行SNV+phasing]
D -->|否 | C
C --> F{需要组间比较?}
E --> F
F -->|否，只需甲基化谱 | G[✅ 完成: 输出 reference-level bed]
F -->|是，两两比较 | H[✅ 完成: 默认执行基础DMR]
F -->|是，多样本群体比较 | I[添加 --dmr_population_scale + --dmr_a/b]
I --> J{选择统计工具}
J -->|5mC/真核 | K[--population_dmrer dss]
J -->|6mA/原核 | L[--population_dmrer modkit]
style C fill:#90EE90
style G fill:#90EE90
style K fill:#E6F3FF
style L fill:#E6F3FF🎯 针对您当前项目 (细菌6mA) 的最终配置
# 场景: 单个细菌菌株的6mA甲基化谱分析
nextflow run nf-core/methylong \
-r 2.0.0 \
-profile docker \
--input samplesheet_6mA.csv \
--outdir methylome_out_6mA \
--ont_aligner minimap2 \ # 避免Docker认证问题
--m6a \ # 启用6mA motif
--skip_snvs \ # ⭐ 细菌无需SNV/phasing
-resume \
-work-dir methylome_out_6mA/work# 场景: 多个细菌菌株的6mA差异甲基化分析
nextflow run nf-core/methylong \
-r 2.0.0 \
-profile docker \
--input samplesheet_6mA.csv \
--outdir methylome_out_6mA \
--ont_aligner minimap2 \
--m6a \
--skip_snvs \
--dmr_population_scale \ # ⭐ 启用群体比较
--population_dmrer modkit \ # ⭐ modkit更好支持6mA
--dmr_a "strain_WT1,strain_WT2" \ # ⭐ 定义组A
--dmr_b "strain_Mut1,strain_Mut2" \# ⭐ 定义组B
-resume \
-work-dir methylome_out_6mA/work💡 一句话总结
| 参数 | 核心作用 | 细菌6mA分析建议 |
|---|---|---|
--skip_snvs |
跳过SNV+phasing | ✅ 必加 (单倍体无需) |
--haplotype_dmrer |
单倍型水平DMR工具 | ⚪ 可选 (配合–skip_snvs时无效) |
--dmr_population_scale |
启用多样本组间比较 | 🔶 按需 (有多样本时添加) |
--population_dmrer |
群体水平DMR工具 | 🔶 如用群体分析，选modkit更兼容6mA |
--dmr_a / --dmr_b |
定义比较的两组样本 | 🔶 配合--dmr_population_scale使用 |
🎯 细菌6mA分析黄金法则:
--m6a + --skip_snvs + --ont_aligner minimap2= ✅ 简洁、高效、无报错！
如有其他参数疑问，欢迎继续讨论.
✅ 您完全正确！抱歉我之前的表述有误
🔧 参数正确写法
--haplotype_dmrer 是一个独立参数，接受 dss 或 modkit 作为值：
# ✅ 正确语法（参数 + 值，空格分隔）
--haplotype_dmrer dss
--haplotype_dmrer modkit
# ❌ 错误语法（连在一起）
--haplotype_dmrer dss # 这不是一个选项名，而是 参数+值从您提供的 help 输出可见：
DMR options
--skip_snvs [boolean] Indicate if to skip snvcall and phase
--haplotype_dmrer [string] DMRer option in DMR analysis for haplotype level, default is dss, specify modkit to switch (accepted: dss, modkit) [default: dss] 📌 默认值已经是
dss，所以如果您不指定该参数，管道会自动使用dss。
🎯 针对细菌6mA分析的最终推荐命令
由于细菌是单倍体 (haploid)，不需要单倍型水平分析，最简洁有效的方案是：
nextflow run nf-core/methylong \
-r 2.0.0 \
-profile docker \
--input samplesheet_6mA.csv \
--outdir methylome_out_6mA \
--ont_aligner minimap2 \
--m6a \
--skip_snvs \ # ⭐ 关键：跳过SNP calling + phasing + haplotype DMR
-resume \
-work-dir methylome_out_6mA/work参数说明：
| 参数 | 作用 | 为什么对细菌分析重要 |
|---|---|---|
--ont_aligner minimap2 |
使用minimap2替代dorado aligner | 避免 quay.io/nanoporetech/dorado 的Docker认证问题 |
--m6a |
启用6mA motif的pileup | 您的研究目标是6mA修饰 |
--skip_snvs |
跳过 SNP calling + WhatsHap phasing + haplotype-level DMR | 🟢 细菌无杂合位点，无需单倍型分析；同时避免 --combine-strands 与 非回文motif “A” 的冲突 |
🔁 如果您想保留 SNP calling（例如用于菌株变异分析）
nextflow run nf-core/methylong \
-r 2.0.0 \
-profile docker \
--input samplesheet_6mA.csv \
--outdir methylome_out_6mA \
--ont_aligner minimap2 \
--m6a \
--haplotype_dmrer dss \ # 保持默认值，显式写出便于阅读
-resume \
-work-dir methylome_out_6mA/work⚠️ 但请注意：即使运行了 phasing 步骤，由于细菌基因组是单倍体，不会产生有意义的 haplotype-level DMR 结果。且如果
modkit pileup在 haplotype 模式下仍使用--combine-strands，仍可能因 “A” 非回文而报错。
🛠️ 如果管道已在中途失败，如何恢复？
# 1. 确认之前的运行已生成主要结果
ls -la methylome_out_6mA/results/methylation/
# 应看到: An6.bed.gz (这就是您需要的6mA参考水平pileup结果！)
# 2. 用 --skip_snvs 跳过失败步骤，完成剩余流程
nextflow run nf-core/methylong \
-r 2.0.0 \
-profile docker \
--input samplesheet_6mA.csv \
--outdir methylome_out_6mA \
--ont_aligner minimap2 \
--m6a \
--skip_snvs \
-resume \
-work-dir methylome_out_6mA/work-resume 会自动复用已完成的步骤，只运行剩余部分 ✅
📋 快速验证：您的6mA结果是否已生成？
# 检查核心输出文件
zcat methylome_out_6mA/results/methylation/An6.bed.gz | head -3
# 示例输出 (6mA bed格式):
# chr1 12345 12346 A 0.85 + 12 10 0.83 0.91 ...
# 统计高置信度6mA位点 (修饰概率 ≥ 0.8)
zcat methylome_out_6mA/results/methylation/An6.bed.gz | \
awk '$4=="A" && $5>=0.8' | wc -l
# 检查MultiQC报告
firefox methylome_out_6mA/results/qc/multiqc_report.html💡 总结 (Zusammenfassung)
🇨🇳 细菌是单倍体 → 无父源/母源染色体 → 无需单倍型分析
✅ 使用 --skip_snvs 跳过 phasing 相关步骤
✅ 使用 --ont_aligner minimap2 避免 Docker 认证问题
✅ 核心结果 An6.bed.gz 已包含您需要的6mA甲基化信息
🇩🇪 Bakterien sind haploid → keine väterlichen/mütterlichen Chromosomen → keine Haplotyp-Analyse nötig
✅ --skip_snvs überspringt Phasing-Schritte
✅ --ont_aligner minimap2 vermeidet Docker-Auth-Probleme
✅ An6.bed.gz enthält bereits alle benötigten 6mA-Methylierungsdaten再次感谢您的指正！🙏 如有其他问题，随时讨论！🚀
✅ 您的理解完全正确！细菌基因组分析不需要单倍型水平分析
🧬 细菌的遗传特性总结
| 特性 | 细菌 (原核生物) | 人类/动植物 (真核生物) |
|---|---|---|
| 染色体倍性 | 🟢 单倍体 (Haploid) | 🔴 二倍体/多倍体 (Diploid/Polyploid) |
| 繁殖方式 | 🟢 无性繁殖 (二分裂) | 🔴 有性繁殖 (减数分裂+受精) |
| 父源/母源染色体 | ❌ 不存在 | ✅ 存在 (väterlich / mütterlich) |
| 杂合SNP (Heterozygous) | ❌ 通常不存在 | ✅ 常见 |
| Phasing/单倍型分型 | ❌ 无生物学意义 | ✅ 必要步骤 |
| Haplotype-level甲基化分析 | ❌ 不适用 | ✅ 用于印记基因/ASM研究 |
🔍 为什么细菌没有 “väterlich / mütterlich” 染色体？
细菌繁殖流程 (二分裂 Binary Fission)：
母细胞 (1个单倍体基因组)
↓
[DNA复制]
↓
2个相同的子代基因组 (克隆)
↓
[细胞分裂]
↓
子细胞A + 子细胞B (遗传物质完全相同)📌 关键点：细菌没有减数分裂 (Meiosis) 和受精过程，因此不存在”父源”和”母源”染色体的概念。
SNP在细菌中的含义：
| 场景 | SNP类型 | 生物学解释 |
|---|---|---|
| 同一菌株内 | 几乎全是纯合 (Homozygous) | 单倍体基因组，每个位点只有一个等位基因 |
| 不同菌株间 | 存在差异 | 用于菌株分型、进化分析，但不是杂合 |
| 混合感染/污染 | 看似”杂合” | 实际是多个菌株的混合，需去卷积而非phasing |
⚠️ 细菌甲基化分析的特殊注意事项
1️⃣ 6mA在细菌中的生物学意义
细菌6mA主要功能：
✅ 限制修饰系统 (Restriction-Modification): 区分"自我"vs"外来"DNA
✅ DNA复制起始调控
✅ 基因表达调控 (如相变 Phase variation)
✅ 错配修复 (MutHLS系统)
→ 这些功能都是"全基因组水平"或"位点特异性"的，
不需要区分父源/母源！2️⃣ modkit参数设置 (针对细菌6mA)
# ✅ 正确命令：细菌6mA分析
modkit pileup sample \
An6_6mA_mapped.mod.bam \
An6_6mA.bed \
--ref An6-unicycler-medaka_polished_genome.fa \
--motif A 0 \ # 6mA motif
--no-combine-strands \ # ⭐ A不是回文，不能合并链
--threads 12
# ❌ 错误命令：会导致报错
modkit pileup ... --motif A 0 --combine-strands ...
# Error! cannot combine strands with a motif that is not a palindrome3️⃣ Nextflow管道参数优化 (细菌场景)
nextflow run nf-core/methylong \
-r 2.0.0 \
-profile docker \
--input samplesheet_6mA.csv \
--outdir methylome_out_6mA \
--ont_aligner minimap2 \
--m6a \
--skip_snvs \ # ⭐ 跳过无意义的SNP phasing步骤
--haplotype_dmrer dss \ # 即使保留，也不会产生有意义的haplotype DMR
-resume \
-work-dir methylome_out_6mA/work🧪 细菌甲基化分析推荐流程
graph TD
A[原始Pod5/FAST5] --> B[Dorado basecalling + 6mA mod-call]
B --> C[生成 modBAM]
C --> D[比对到参考基因组 minimap2]
D --> E[modkit pileup --motif A 0 --no-combine-strands]
E --> F[生成 bed/bigWig]
F --> G[下游分析]
G --> G1[甲基化位点鉴定]
G --> G2[motif富集分析]
G --> G3[与基因注释关联]
G --> G4[菌株间比较甲基化]
style E fill:#90EE90,stroke:#333
style G fill:#E6F3FF,stroke:#333推荐下游分析工具 (细菌专用)：
| 分析目标 | 推荐工具 | 说明 |
|---|---|---|
| 6mA位点鉴定 | modkit, nanom6A |
基于ONT信号或modBAM |
| Motif发现 | MEME, HOMER |
寻找6mA富集的序列motif |
| 基因关联 | bedtools intersect |
关联6mA位点与基因/启动子 |
| 菌株比较 | DSS, methylKit (适配后) |
注意：用reference-level模式 |
| 可视化 | IGV, pyGenomeTracks |
加载bigWig轨道 |
🎯 针对您当前项目的最终建议
1. 立即修复管道错误：
# 方案A: 跳过单倍型分析 (最简单)
nextflow run nf-core/methylong \
--input samplesheet_6mA.csv \
--outdir methylome_out_6mA \
--ont_aligner minimap2 \
--m6a \
--skip_snvs \ # ⭐ 关键：细菌不需要phasing
-resume
# 方案B: 如果管道已生成主要结果，手动完成最后一步
cd methylome_out_6mA/results/methylation/
# An6.bed.gz 已经是您需要的6mA参考水平结果！✅2. 验证结果质量：
# 检查6mA位点数量和质量
zcat methylome_out_6mA/results/methylation/An6.bed.gz | \
awk '$4=="A" && $10>=0.8' | wc -l # 高置信度6mA位点
# 检查motif富集 (示例：GATC是常见6mA motif)
zcat methylome_out_6mA/results/methylation/An6.bed.gz | \
awk '$4=="A"' | cut -f1-3 | \
bedtools getfasta -fi genome.fa -bed - | \
grep -o "GATC" | wc -l3. 后续分析方向：
🔹 6mA位点是否富集在启动子区域？
🔹 是否集中在特定motif (如GATC, GANTC)?
🔹 与基因表达数据关联 (如果有RNA-seq)?
🔹 与其他菌株比较甲基化差异?💡 一句话总结 (Zusammenfassung)
🇨🇳 细菌是单倍体，无父源/母源染色体之分，因此单倍型水平甲基化分析无生物学意义。请使用参考基因组水平分析 +
--no-combine-strands参数进行6mA研究。🇩🇪 Bakterien sind haploid und haben keine väterlichen/mütterlichen Chromosomen. Daher ist eine Haplotyp-Level-Methylierungsanalyse biologisch nicht sinnvoll. Verwenden Sie Reference-Level-Analyse mit
--no-combine-strandsfür 6mA-Studien.
Ihre Verständnis ist absolut korrekt! 🎯 如有其他问题，随时讨论！🚀
🧬 单倍型水平分析 (Haplotype-level) vs 参考基因组水平分析 (Reference-level) 详解
📊 核心概念对比
| 维度 | 参考基因组水平 (Reference-level) | 单倍型水平 (Haplotype-level) |
|---|---|---|
| 基本单位 | 整个样本的所有reads合并统计 | 按父源/母源染色体分开统计 |
| 输入数据 | 原始align的BAM文件 | 经过phasing的haplotagged BAM |
| 甲基化值计算 | 所有reads在该位点的平均修饰概率 | 分别计算H1/H2两个单倍型的修饰概率 |
| 输出结果 | 每个位点1个甲基化值 | 每个位点2个甲基化值 (H1 + H2) |
| 生物学意义 | 样本整体甲基化状态 | 等位基因特异性甲基化 (ASM) |
🔍 直观示例：同一个CpG位点
假设基因组位置 chr1:1000 是一个CpG位点：
📋 Reference-level 分析结果
# 所有reads合并统计
位置 修饰类型 修饰次数 总覆盖度 甲基化比例
chr1:1000 5mC 45 100 45%表示：该样本在这个位置平均有45%的DNA分子被甲基化
📋 Haplotype-level 分析结果
# 按单倍型分开统计
位置 单倍型 修饰类型 修饰次数 总覆盖度 甲基化比例
chr1:1000 H1 5mC 40 50 80% ← 父源染色体
chr1:1000 H2 5mC 5 50 10% ← 母源染色体表示：父源染色体在该位点高度甲基化(80%)，母源染色体几乎未甲基化(10%)
→ 这是典型的等位基因特异性甲基化 (Allele-Specific Methylation, ASM)
🔄 技术流程对比
graph LR
A[原始reads] --> B[比对到参考基因组]
B --> C[Reference-level pileup]
B --> D[SNP calling]
D --> E[Phasing/单倍型分型]
E --> F[给reads添加HP标签]
F --> G[Haplotype-level pileup]
C --> H[输出: 每个位点1个值]
G --> I[输出: 每个位点2个值 H1+H2]关键步骤说明：
- SNP Calling (Clair3): 找出样本中的杂合位点
- Phasing (WhatsHap): 根据SNP连锁关系，将reads分配到父源/母源
- HP tagging: 在BAM文件中给每条read添加
HP:i:1或HP:i:2标签 - Strand处理:
- 参考水平: 可用
--combine-strands合并正负链 - 单倍型水平: 非回文motif(如6mA的”A”)不能合并链 ⚠️
- 参考水平: 可用
⚠️ 为什么6mA分析时会出现 “cannot combine strands” 错误？
回文 vs 非回文 Motif
| 修饰类型 | Motif | 是否回文 | 正负链关系 | 能否combine-strands |
|---|---|---|---|---|
| 5mC | CG |
✅ 是 | CG(正链) ↔ CG(负链) | ✅ 可以 |
| 6mA | A |
❌ 否 | A(正链) ↔ T(负链) | ❌ 不可以 |
原理解释：
- 5mC在CG位点: 正链的
C和负链的C是同一个生物学事件的两面，可以合并统计 - 6mA在A位点: 正链的
A甲基化和负链的T(对应正链A的互补)是两个独立的生物学事件，不能合并！
# ❌ 错误命令 (6mA + combine-strands)
modkit pileup --motif A 0 --combine-strands ...
# Error! cannot combine strands with a motif that is not a palindrome
# ✅ 正确命令 (6mA + 分开链)
modkit pileup --motif A 0 --no-combine-strands ...
# 或分别处理: --motif A+ 0 (正链) 和 --motif A- 0 (负链)🎯 应用场景：什么时候用哪种分析？
✅ 使用 Reference-level 分析，如果：
- 初步探索性分析，了解整体甲基化模式
- 样本是单倍体生物（如细菌、线粒体）
- 关注的是区域水平的甲基化变化（如启动子、增强子）
- 计算资源有限，需要快速出结果
- 分析6mA且不需要等位基因特异性信息
✅ 使用 Haplotype-level 分析，如果：
- 研究印记基因 (Imprinted genes) 的甲基化调控
- 分析杂合位点附近的等位基因特异性甲基化 (ASM)
- 探索顺式调控变异 (cis-regulatory variants) 对甲基化的影响
- 癌症研究中分析等位基因失衡的表观遗传改变
- 需要区分父源/母源染色体的表观遗传状态
📈 结果解读示例
案例：印记基因 IGF2 的启动子区域
# Reference-level 结果
区域: IGF2_promoter
平均甲基化: 52%
→ 结论：中等甲基化，无明显特征
# Haplotype-level 结果
区域: IGF2_promoter
H1(父源): 95% 甲基化
H2(母源): 8% 甲基化
→ 结论：典型的父源特异性甲基化，符合已知印记模式！✅案例：6mA在细菌基因组
# 细菌通常是单倍体，无杂合位点
→ Haplotype-level 分析无意义
→ 直接使用 Reference-level + --motif A 0 --no-combine-strands🔧 针对您当前6mA分析的实操建议
方案1：跳过单倍型分析（推荐 ✅）
# 添加 --skip_snvs 参数，跳过phasing相关步骤
nextflow run nf-core/methylong \
--input samplesheet_6mA.csv \
--outdir methylome_out_6mA \
--ont_aligner minimap2 \
--m6a \
--skip_snvs \ # ⭐ 关键：跳过单倍型分析
-resume方案2：保留单倍型BAM，但手动运行modkit
# 1. 让管道运行到生成 haplotagged BAM 后暂停
# 2. 手动执行正确参数的modkit:
modkit pileup sample \
An6_haplotagged.bam \
An6_6mA_haplo.bed \
--ref genome.fa \
--motif A 0 \
--partition-tag HP \
--no-combine-strands \ # ⭐ 关键：不合并链
--threads 12
# 3. 用R/Python后续分析DSS等方案3：如果确实需要单倍型6mA分析
# 分别计算正负链，后期在R中合并
# 正链6mA (A+)
modkit pileup sample ... --motif A+ 0 --partition-tag HP ...
# 负链6mA (实际是T位的互补信息，通常不关注)
# 一般6mA分析只关注正链的A即可📦 输出文件对比
| 文件 | Reference-level | Haplotype-level |
|---|---|---|
*.bed.gz |
✅ 每个位点1行 | ✅ 每个位点2行 (H1+H2) |
*.bigWig |
✅ 单轨可视化 | ✅ 双轨可视化 (H1轨 + H2轨) |
| DMR结果 | 样本间差异 | 单倍型间差异 + 样本间差异 |
| 计算时间 | ⚡ 较快 | 🐌 较慢 (需phasing+双倍pileup) |
💡 总结一句话
参考水平分析回答：”这个位点在整个样本中甲基化程度如何？”
单倍型水平分析回答：”这个位点在父源和母源染色体上的甲基化是否不同？”
对于6mA在原核生物或无需等位基因分辨的研究，参考水平分析 + 正确设置 --no-combine-strands 通常就足够了！🎯
如有具体数据分析需求，欢迎继续讨论！🚀
🔧 modkit Error: “cannot combine strands with a motif that is not a palindrome”
📋 Fehleranalyse
✅ Gute Nachrichten: Die Pipeline ist jetzt fast durchgelaufen!
- Alignment mit minimap2 ✔️
- Modkit pileup (Referenz-Level) ✔️
- Clair3 SNP-Calling ✔️
- WhatsHap Phasing ✔️
❌ Fehler im letzten Schritt:
Error! cannot combine strands with a motif that is not a palindrome🔍 Ursache
Der Prozess MODKIT_PILEUP_HAPLOTYPE_LEVEL versucht, Stränge zu kombinieren (--combine-strands) für das 6mA-Motiv “A”:
modkit pileup --motif CG 0 --motif A 0 --combine-strands --partition-tag HP ...| Motiv | Palindrom? | Kombinierbar? |
|---|---|---|
CG (5mC) |
✅ Ja (CG ↔ CG) | ✅ Ja |
A (6mA) |
❌ Nein (A ↔ T) | ❌ Nein! |
Für nicht-palindromische Motive wie 6mA (
A) ist die Modifikations-Kontext-Information strang-spezifisch. Forward- und Reverse-Strand können nicht sinnvoll kombiniert werden [[modkit docs]].
✅ Lösungen (absteigend nach Einfachheit)
🔹 Lösung 1: Haplotype-Level DMR überspringen (empfohlen)
Da Sie 6mA analysieren, ist die Referenz-Level-Pileup-Datei (An6.bed.gz) bereits das wichtigste Ergebnis. Die haplotype-level Analyse ist optional.
Option A: Mit --skip_snvs (überspringt SNP-Calling + Phasing + haplotype DMR):
nextflow run nf-core/methylong \
-r 2.0.0 \
-profile docker \
--input samplesheet_6mA.csv \
--outdir methylome_out_6mA \
--ont_aligner minimap2 \
--m6a \
--skip_snvs \ # ⭐ Haplotype-Level Schritte überspringen
-resume \
-work-dir methylome_out_6mA/workOption B: Nur DMR-Schritt deaktivieren via Custom-Config
Erstellen Sie skip_haplo_dmr.config:
process {
withName: '*DSS_HAPLOTYPE_LEVEL*' {
enabled = false
}
}Dann starten mit:
nextflow run nf-core/methylong ... -c skip_haplo_dmr.config -resume🔹 Lösung 2: --haplotype_dmrer modkit verwenden (experimentell)
Der alternative DMRer modkit verwendet möglicherweise andere Parameter:
nextflow run nf-core/methylong \
-r 2.0.0 \
-profile docker \
--input samplesheet_6mA.csv \
--outdir methylome_out_6mA \
--ont_aligner minimap2 \
--m6a \
--haplotype_dmrer modkit \ # ⭐ Alternative zu 'dss' (Standard)
-resume \
-work-dir methylome_out_6mA/work⚠️ Falls auch
modkitals DMRer--combine-strandsverwendet, tritt derselbe Fehler auf. Dann Lösung 1 verwenden.
🔹 Lösung 3: Custom modkit-Parameter via Config (fortgeschritten)
Wenn Sie unbedingt haplotype-level pileup für 6mA benötigen, entfernen Sie --combine-strands für das “A”-Motif:
-
Custom-Config
modkit_6mA.config:process { withName: '*MODKIT_PILEUP_HAPLOTYPE_LEVEL*' { // Override command via script customization // Hinweis: Dies erfordert ggf. Pipeline-Fork oder pre-run script hook } } -
Alternativ: Manueller modkit-Lauf nach Pipeline-Ende:
# Nach erfolgreichem Pipeline-Lauf (bis vor DSS_HAPLOTYPE_LEVEL) cd methylome_out_6mA/
Modkit pileup für 6mA OHNE –combine-strands
modkit pileup sample \ work/…/An6_haplotagged.bam \ An6_6mA_haplo.bed \ –ref bacass_out/Medaka/An6-unicycler-medaka_polished_genome.fa \ –motif A 0 \ –partition-tag HP \ –threads 12 \ –no-combine-strands # ⭐ Explizit deaktivieren
Ergebnis weiterverwenden (z.B. mit DSS in R)
---
## 📁 Wo finde ich meine Ergebnisse jetzt schon?
Auch ohne haplotype-level DMR sind die **wichtigsten Outputs** bereits vorhanden:
```bash
methylome_out_6mA/
├── results/
│ ├── methylation/
│ │ ├── An6.bed.gz # ⭐ Hauptergebnis: 6mA-Pileup (Referenz-Level)
│ │ ├── An6.bed.gz.tbi
│ │ └── An6.bigWig # Für IGV-Visualisierung
│ ├── variants/
│ │ ├── An6.vcf.gz # SNP-Calling-Ergebnisse
│ │ └── An6_haplotagged.bam # Phased BAM für Haplotyp-Analyse
│ └── qc/
│ └── multiqc_report.html # Qualitätsübersicht
└── work/.../ # Zwischenergebnisse (inkl. haplotagged BAM)🔎 Schnell-Check der 6mA-Ergebnisse:
# 1. Pileup-Datei prüfen
zcat methylome_out_6mA/results/methylation/An6.bed.gz | head -5
# 2. 6mA-Sites zählen (Motif "A")
zcat methylome_out_6mA/results/methylation/An6.bed.gz | \
awk '$4 == "A"' | wc -l
# 3. MultiQC-Report öffnen
firefox methylome_out_6mA/results/qc/multiqc_report.html🎯 Empfohlener nächster Befehl
# Haplotype-level DMR überspringen, Rest fortsetzen:
nextflow run nf-core/methylong \
-r 2.0.0 \
-profile docker \
--input samplesheet_6mA.csv \
--outdir methylome_out_6mA \
--ont_aligner minimap2 \
--m6a \
--skip_snvs \
-resume \
-work-dir methylome_out_6mA/work💡 Tipp: Die Referenz-Level-Pileup-Datei (
An6.bed.gz) enthält bereits alle 6mA-Modifikationswahrscheinlichkeiten pro Position – für die meisten Analysen völlig ausreichend!
Bei weiteren Fragen helfe ich gerne! 🚀
在您的命令中，--m6a 是一个布尔型开关参数，核心作用是：
🔹 指示管道在甲基化堆叠（pileup）步骤中，专门识别、提取并统计 6mA（N6-甲基腺嘌呤）修饰位点。
🔍 详细机制说明
| 维度 | 说明 |
|---|---|
| 默认行为 | modkit（管道内部调用的甲基化定量工具）默认主要针对真核生物常见的 胞嘧啶修饰（如 CG 位点的 5mC/5hmC）进行统计。 |
添加 --m6a 后 |
管道会在底层自动为 modkit pileup 命令追加 --motif A 0 参数，强制其扫描并输出所有 腺嘌呤（A） 位点的修饰概率。 |
| 底层命令映射 | bash<br># 不加 --m6a<br>modkit pileup ... --motif CG 0 ...<br><br># 加上 --m6a 后<br>modkit pileup ... --motif CG 0 --motif A 0 ...<br> |
🦠 为什么对您的细菌分析必不可少？
- 细菌的主要修饰是 6mA：原核生物基因组中，6mA 通常由限制修饰系统（如 Dam, CcrM 甲基转移酶）催化，富集在
GATC、GANTC等 motif 中。 - 不加的后果：如果不写
--m6a，管道只会输出默认的 5mC/CG 结果。您的核心输出文件An6.bed.gz中将几乎没有或完全没有 6mA 数据，后续 motif 分析和功能注释将全部落空。 - 结果验证：运行后，您可以用以下命令确认 6mA 信号已成功捕获：
zcat methylome_out_6mA/results/methylation/An6.bed.gz | \ awk '$4=="A" && $5>=0.8' | wc -l # 预期输出：数百至数万行（取决于测序深度和菌株甲基化活性）
⚠️ 使用注意事项
- ✅ 正确写法：直接写
--m6a即可（布尔开关，不需要写--m6a true或--m6a=1，否则可能触发参数校验警告）。 - 🔄 与其他参数配合：通常与
--ont_aligner minimap2和--skip_snvs一起使用，构成细菌 6mA 分析的最小完整命令集。
💡 一句话总结
--m6a= “请专门为我提取细菌基因组中的 6mA 修饰信号，不要只跑默认的 5mC/CG。”
对于原核表观遗传研究，这是必加参数，直接决定您能否拿到目标数据！🎯