🧰 甲基化下游统计分析：独立工具推荐与使用指南
您不需要 Nextflow 封装，以下是优秀的独立工具 (standalone tools)，可直接用于 modkit pileup 结果的统计分析，功能对标 TRANSIT：
📊 工具总览：按 TRANSIT 功能映射
| TRANSIT 功能 | 甲基化对应工具 | 类型 | 安装方式 | 核心优势 |
|---|---|---|---|---|
| TTR 标准化 | DSS::normalizeCoverage() |
R 包 | conda install bioconductor-dss |
贝叶斯收缩，小样本稳健 |
| ANOVA 多组比较 | DSS::DMLtest() / methylKit |
R 包 | conda install bioconductor-methylkit |
支持任意分组设计 + FDR 校正 |
| Tn5Gaps (非插入运行) | metilene / DMRfinder |
C++ / Python | conda install -c bioconda metilene |
快速检测连续低/高甲基化区域 |
| 基因长度标准化 | bedtools + 自定义脚本 |
CLI 工具 | conda install -c bioconda bedtools |
灵活适配细菌基因密集特点 |
| 置换检验 + FDR | DSS (内置) / methylKit |
R 包 | 同上 | 自动多重校正，输出可直接发表 |
| Circos 可视化 | pyGenomeTracks / deepTools |
Python | pip install pygenometracks |
出版级基因组轨道图，易定制 |
🔧 工具详解：安装 + 使用示例
1️⃣ DSS (Dispersion Shrinkage for Sequencing) — 首选推荐 ⭐
🎯 对标 TRANSIT ANOVA：位点水平差异甲基化检测，支持多样本、多组比较
🔹 安装
# 推荐：通过 conda (自动解决 R 依赖)
conda create -n methyl_stats -c conda-forge -c bioconda r-base r-dss
conda activate methyl_stats
# 或手动在 R 中安装
R
> if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager")
> BiocManager::install("DSS")🔹 输入准备 (modkit bed → DSS 格式)
# 提取核心列: chrom, start, end, mod_count, coverage
# modkit 输出: $4=修饰码, $5=覆盖度, $10=修饰计数, $11=概率(0-100)
awk '$4 == "a" {printf "%s\t%d\t%d\t%d\t%d\n", $1, $2, $3, $10, $5}' An6_6mA_raw.bed > An6.dss
awk '$4 == "a" {printf "%s\t%d\t%d\t%d\t%d\n", $1, $2, $3, $10, $5}' BG5_6mA_raw.bed > BG5.dss🔹 R 分析脚本 (run_dss_6mA.R)
#!/usr/bin/env Rscript
library(DSS)
library(qvalue) # 用于更灵活的 FDR 校正
# 1. 读取数据
read_dss <- function(f) {
df <- read.table(f, header=F, col.names=c("chr","start","end","X","N"))
df$pos <- df$start + 1 # BED 0-based → 1-based
df[, c("chr","pos","N","X")] # DSS 要求格式
}
An6 <- read_dss("An6.dss")
BG5 <- read_dss("BG5.dss")
# 2. 构建 BSseq 对象
BSobj <- makeBSseqData(list(An6, BG5), sampleNames=c("An6", "BG5"))
# 3. 差异检验 (贝塔 - 二项模型 + 平滑)
dmlTest <- DMLtest(BSobj, group1=1, group2=2, smoothing=TRUE)
# 4. 提取显著位点 (FDR < 0.05, |diff| > 20%)
sigDML <- getSig(dmlTest, pval=0.05, delta=0.2)
# 5. 输出结果
write.table(sigDML, "An6_vs_BG5_6mA_DMLs.tsv",
sep="\t", quote=F, row.names=F)
# 6. (可选) 按基因注释聚合
# 需准备 genome.gff，用 ChIPseeker 或自定义脚本
cat("✅ 完成！显著差异位点数:", nrow(sigDML), "\n")🔹 运行
Rscript run_dss_6mA.R🔹 输出解读 (对标 TRANSIT ANOVA)
| DSS 输出列 | 含义 | TRANSIT 对应列 |
|---|---|---|
chr, pos |
基因组坐标 | Orf, TA_pos |
coverage |
总覆盖深度 | TAs |
group1, group2 |
两组修饰计数 | Mean_[cond] |
pval, p.adj |
原始/FDR 校正 p 值 | Pval, Padj |
diff |
甲基化差异 (组1 – 组2) | LFC_[cond] |
abs.diff |
绝对差异值 | – |
2️⃣ metilene — 快速 DMR 检测 (对标 Tn5Gaps)
🎯 对标 TRANSIT Tn5Gaps：检测连续低/高甲基化区域，基于二元分割 + 置换检验
🔹 安装
conda install -c bioconda metilene
# 或下载预编译二进制: https://github.com/metilene/metilene/releases🔹 输入准备
# metilene 需要 bedGraph 格式: chrom, start, end, value
# 用 modkit 概率列 ($11) 作为 value
awk '$4 == "a" {printf "%s\t%d\t%d\t%.2f\n", $1, $2, $3, $11}' An6_6mA_raw.bed > An6.bedgraph
awk '$4 == "a" {printf "%s\t%d\t%d\t%.2f\n", $1, $2, $3, $11}' BG5_6mA_raw.bed > BG5.bedgraph
# 准备基因组大小文件 (chrom, size)
samtools faidx genome.fa | cut -f1,2 > genome.sizes🔹 运行 DMR 检测
metilene \
-i An6.bedgraph BG5.bedgraph \
-g genome.sizes \
-o dmr_results.tsv \
-t 8 \ # 线程数
-p 0.05 \ # FDR 阈值
-d 0.2 # 最小甲基化差异 (20%)🔹 输出解读
# dmr_results.tsv 示例:
chr1 12345 12890 0.85 0.32 0.53 1.2e-05 0.003 hyper
# │ │ │ │ │ │ │ │ │
# │ │ │ │ │ │ │ │ └─ 方向: hyper/hypo
# │ │ │ │ │ │ │ └─ FDR 校正 p 值
# │ │ │ │ │ │ └─ 原始 p 值
# │ │ │ │ │ └─ 甲基化差异 (An6 - BG5)
# │ │ │ │ └─ BG5 平均甲基化
# │ │ │ └─ An6 平均甲基化
# │ │ └─ 区域结束
# │ └─ 区域起始
# └─ 染色体3️⃣ methylKit — 基因/区域水平聚合分析
🎯 对标基因水平分析：将位点聚合到基因/操纵子，进行组间比较
🔹 安装
conda install -c bioconda bioconductor-methylkit🔹 R 脚本示例 (methylkit_gene_level.R)
#!/usr/bin/env Rscript
library(methylKit)
library(GenomicFeatures)
# 1. 读取 modkit bed (需先转为 methylKit 格式)
# 格式: chr, start, end, strand, coverage, freq
read_modkit <- function(f, sample_id) {
df <- read.table(f, header=F)
df <- df[df$V4 == "a" & df$V11 >= 70 & df$V5 >= 10, ] # 6mA, prob≥70%, cov≥10
data.frame(
chr = df$V1,
start = df$V2 + 1, # 0-based → 1-based
end = df$V3,
strand = ifelse(df$V6 == "+", 1, -1),
coverage = df$V5,
freq = df$V11 # 概率 0-100
)
}
# 2. 加载样本
An6 <- read_modkit("An6_6mA_raw.bed", "An6")
BG5 <- read_modkit("BG5_6mA_raw.bed", "BG5")
# 3. 创建 methylKit 对象
myobj <- methRead(
list("An6.tmp", "BG5.tmp"), # 需先保存为临时文件
sample.id = c("An6", "BG5"),
assembly = "An6_genome",
treatment = c(1, 2), # 实验分组
context = "CpG" # 细菌可设为"any"
)
# 4. 基因水平聚合 (需 GFF 注释)
txdb <- makeTxDbFromGFF("genome.gff", format="gff3")
gene_meth <- tileMethylCounts(myobj, win.size=1000, step.size=500)
# 5. 差异分析
diff_meth <- calculateMethylDiff(gene_meth, difference=20, qvalue=0.01)
# 6. 输出
write.table(diff_meth, "gene_level_6mA_diff.tsv", sep="\t", quote=F)4️⃣ 可视化工具：pyGenomeTracks (对标 Circos)
🎯 对标 Circos：生成线性基因组轨道图，支持多条件叠加
🔹 安装
pip install pygenometracks🔹 准备配置文件 (tracks.ini)
[x-axis]
where = top
# 6mA 信号轨道 (An6)
[An6_6mA]
file = An6_6mA.bw
title = An6 6mA
height = 2
color = #1f78b4
min_value = 0
max_value = 100
# 6mA 信号轨道 (BG5)
[BG5_6mA]
file = BG5_6mA.bw
title = BG5 6mA
height = 2
color = #e31a1c
min_value = 0
max_value = 100
# 差异区域标记
[DMRs]
file = dmr_results.bed
title = Significant DMRs
height = 1
color = #ff7f00🔹 生成 bigWig (modkit bed → bigWig)
# 1. bed → bedGraph
awk '$4 == "a" {printf "%s\t%d\t%d\t%.2f\n", $1, $2, $3, $11}' An6_6mA_raw.bed > An6_6mA.bg
# 2. bedGraph → bigWig (需 genome.sizes)
bedGraphToBigWig An6_6mA.bg genome.sizes An6_6mA.bw🔹 运行可视化
pygenometracks \
--tracks tracks.ini \
--region chr1:100000-200000 \ # 可选：指定区域
--width 10 \
--height 5 \
--output methylation_tracks.png \
--dpi 300🎯 针对您项目的最小可行流程
# ✅ 步骤 1: 格式转换 (modkit bed → 各工具输入)
# DSS 输入
awk '$4 == "a" {printf "%s\t%d\t%d\t%d\t%d\n", $1, $2, $3, $10, $5}' An6_6mA_raw.bed > An6.dss
# metilene 输入
awk '$4 == "a" {printf "%s\t%d\t%d\t%.2f\n", $1, $2, $3, $11}' An6_6mA_raw.bed > An6.bedgraph
# bigWig 输入 (可视化)
awk '$4 == "a" {printf "%s\t%d\t%d\t%.2f\n", $1, $2, $3, $11}' An6_6mA_raw.bed | \
bedGraphToBigWig stdin genome.sizes An6_6mA.bw
# ✅ 步骤 2: 运行统计分析
Rscript run_dss_6mA.R # 位点水平差异
metilene -i An6.bedgraph BG5.bedgraph -g genome.sizes -o dmr.tsv # 区域水平
# ✅ 步骤 3: 可视化
pygenometracks --tracks tracks.ini --output final_tracks.png
# ✅ 步骤 4: 整合结果 (可选)
# 将 DSS 显著位点 + metilene DMR 区域 + 基因注释合并为最终表格📦 一键安装脚本 (conda 环境)
#!/bin/bash
# install_methyl_tools.sh
conda create -n methyl_analysis -c conda-forge -c bioconda \
r-base=4.3 \
r-dss \
r-methylkit \
r-chipseeker \
r-qvalue \
metilene \
bedtools \
samtools \
ucsc-bedgraphtobigwig \
-y
conda activate methyl_analysis
# 额外安装 Python 可视化工具
pip install pygenometracks deepTools
echo "✅ 所有工具安装完成！激活环境: conda activate methyl_analysis"💡 工具选择决策树
您的分析目标？
│
├─ 🔍 "哪些位点在条件间甲基化水平不同？"
│ └─→ ✅ DSS (最严谨，输出可直接用于发表)
│
├─ 🧬 "哪些连续区域甲基化状态改变？" (对标 Tn5Gaps)
│ └─→ ✅ metilene (最快) 或 DMRfinder (最灵活)
│
├─ 📊 "基因/操纵子水平聚合比较？"
│ └─→ ✅ methylKit + ChIPseeker 注释
│
├─ 🎨 "生成出版级可视化？"
│ └─→ ✅ pyGenomeTracks (线性) 或 Circos (环形，需额外配置)
│
└─ 🔎 "交互式探索关键位点？"
└─→ ✅ IGV + bigWig (最直观)🎯 一句话总结
🔑 核心推荐组合：
- 差异分析 →
DSS(R 包，对标 TRANSIT ANOVA)- 区域检测 →
metilene(C++ 二进制，对标 Tn5Gaps)- 可视化 →
pyGenomeTracks(Python，对标 Circos)🚀 安装最简单：全部可通过
conda或pip一键安装，无需编译
📈 您的 2,234 个高置信 6mA 位点 + 上述工具 = 完整、可发表的甲基化分析流程！
如需具体脚本调试或结果解读帮助，请随时告知！🔧🚀
🔍 nf-core/methylong 内置 DMR 选项：能替代哪些独立工具？
✅ 简短回答
| 您之前需要的功能 | methylong 内置选项 | 能否替代独立工具？ |
|---|---|---|
| 位点水平差异分析 (对标 TRANSIT ANOVA) | --population_dmrer dss + --dmr_a/--dmr_b |
✅ 完全可以替代独立 DSS |
| 区域水平聚合比较 | --population_dmrer dss (自动区域聚合) |
✅ 部分替代 methylKit |
| 连续低/高甲基化区域检测 (对标 Tn5Gaps) | ❌ 无专门功能 | ❌ 仍需 metilene 或自定义脚本 |
| 基因/操纵子注释关联 | ❌ 不直接整合 GFF | ❌ 仍需 bedtools/ChIPseeker |
| 出版级可视化 | ❌ 只输出 bed/bigWig | ❌ 仍需 pyGenomeTracks/IGV |
🎯 核心结论：methylong 的
--population_dmrer dss可以完全替代您独立运行 DSS 进行组间差异甲基化分析，无需额外工具！
📋 methylong DMR 选项详解 (对标 TRANSIT)
🔹 参数功能映射表
| methylong 参数 | 作用 | TRANSIT 对应功能 | 统计方法 |
|---|---|---|---|
--skip_snvs |
跳过 SNP calling + phasing | – (细菌无需) | – |
--haplotype_dmrer [dss/modkit] |
单倍型水平 DMR 工具 | – (细菌无单倍型) | 贝叶斯层次模型 (DSS) |
--population_dmrer [dss/modkit] |
⭐ 群体水平 DMR 工具 | ✅ ANOVA 多组比较 | 贝塔 – 二项检验 + FDR |
--dmr_population_scale |
启用多样本组间比较 | ✅ 条件间差异分析 | 同上 |
--dmr_a, --dmr_b |
定义两个比较组的样本名 | ✅ 指定对比条件 | – |
🚀 用 methylong 内置 DSS 替代独立运行的完整命令
场景：细菌 6mA，An6 vs BG5 差异分析
nextflow run nf-core/methylong \
-r 2.0.0 \
-profile docker \
--input samplesheet_6mA.csv \ # 包含 An6 和 BG5 两行
--outdir methylome_out_6mA \
--ont_aligner minimap2 \
--m6a \
--skip_snvs \ # 细菌无需 phasing
--dmr_population_scale \ # ⭐ 启用群体水平比较
--population_dmrer dss \ # ⭐ 使用 DSS 进行统计 (默认值，可省略)
--dmr_a "An6" \ # ⭐ 定义组 A (samplesheet 中的 group 列)
--dmr_b "BG5" \ # ⭐ 定义组 B
-resume \
-work-dir methylome_out_6mA/work🔹 samplesheet_6mA.csv 格式要求
group,sample,path,ref,method
An6,An6,/path/An6_6mA_mapped.mod.bam,/path/genome.fa,ont
BG5,BG5,/path/BG5_6mA_mapped.mod.bam,/path/genome.fa,ont⚠️ 关键：
group列的值 (An6,BG5) 必须与--dmr_a/--dmr_b完全一致！
📁 输出文件：methylong 内置 DMR 结果位置
methylome_out_6mA/
└── results/
└── dmr/
├── population_scale/
│ ├── An6_vs_BG5.DSS_DMLs.tsv # ⭐ 位点水平差异结果 (对标 TRANSIT anova_out.xls)
│ ├── An6_vs_BG5.DSS_DMRs.tsv # 区域水平差异结果
│ └── DSS_summary.html # 简要统计摘要
└── qc/
└── dmr_multiqc_report.html # DMR 分析质控🔹 DSS_DMLs.tsv 列解读 (对标 TRANSIT ANOVA 输出)
| 列名 | 含义 | TRANSIT 对应列 |
|---|---|---|
chr, pos |
基因组坐标 | Orf, TA_pos |
coverage_A, coverage_B |
两组覆盖深度 | TAs |
mod_count_A, mod_count_B |
两组修饰计数 | Mean_[condition] |
methylation_A, methylation_B |
两组甲基化概率 (0-100) | – |
diff |
甲基化差异 (A – B) | LFC_[condition] |
pval, p.adj |
原始/FDR 校正 p 值 | Pval, Padj |
significant |
是否显著 (FDR20%) | status |
⚠️ methylong 内置 DMR 的局限性 (仍需外部工具的场景)
❌ 无法替代的功能 1：连续区域检测 (对标 Tn5Gaps)
# methylong 输出的是"离散位点"差异，不是"连续非修饰区域"
# 如果您需要检测类似 Tn5Gaps 的"长片段无 6mA 区域"，仍需：
metilene \
-i An6.bedgraph BG5.bedgraph \
-g genome.sizes \
-o continuous_dmr.tsv \
-p 0.05 -d 0.2❌ 无法替代的功能 2：基因/操纵子水平注释
# methylong 不直接整合 GFF 注释
# 如需"哪些基因的 6mA 水平显著改变"，仍需：
bedtools intersect -a An6_vs_BG5.DSS_DMLs.tsv -b genome.gff -wa -wb > genes_with_diff_6mA.tsv❌ 无法替代的功能 3：高级可视化
# methylong 只输出 bed/bigWig，不生成出版级图
# 如需 Circos 式基因组视图，仍需：
pygenometracks --tracks tracks.ini --output final_figure.png🎯 推荐工作流：methylong 内置 + 最小外部工具
graph LR
A[modkit pileup] --> B[methylong with --population_dmrer dss]
B --> C[An6_vs_BG5.DSS_DMLs.tsv]
C --> D{需要连续区域？}
D -->|是 | E[metilene for continuous DMRs]
D -->|否 | F[直接进入注释]
E & F --> G[bedtools + genome.gff for gene annotation]
G --> H[pyGenomeTracks for visualization]✅ 最小命令集：
# 1️⃣ 主分析：methylong 内置 DSS (替代独立 DSS 运行)
nextflow run nf-core/methylong ... --dmr_population_scale --dmr_a An6 --dmr_b BG5 ...
# 2️⃣ (可选) 连续区域检测：仍需 metilene
awk '$4 == "a" {printf "%s\t%d\t%d\t%.2f\n", $1, $2, $3, $11}' An6_6mA_raw.bed > An6.bg
metilene -i An6.bg BG5.bg -g genome.sizes -o continuous_dmr.tsv
# 3️⃣ (可选) 基因注释：仍需 bedtools
bedtools intersect -a methylome_out_6mA/results/dmr/population_scale/An6_vs_BG5.DSS_DMLs.tsv \
-b genome.gff -wa -wb > annotated_diff_6mA.tsv💡 一句话总结
🔑 methylong 内置 DMR 选项可以：
- ✅ 完全替代独立 DSS 进行位点水平差异分析 (
--population_dmrer dss)- ✅ 自动完成：标准化 + 贝塔 – 二项检验 + FDR 校正 + 结果输出
🔑 但仍需外部工具：
- ❌ 连续区域检测 (→
metilene)- ❌ 基因注释关联 (→
bedtools/ChIPseeker)- ❌ 出版级可视化 (→
pyGenomeTracks)🎯 您的最优路径：先用 methylong 内置 DSS 拿到差异位点列表，再按需添加 1-2 个外部工具完成注释和可视化！
如需我帮您编写 samplesheet 或调试 methylong DMR 参数，请随时告知！🚀
📊 modkit pileup 输出文件列详解 (bedMethyl 格式)
🎯 您的结果统计
| 过滤条件 | 位点数量 | 说明 |
|---|---|---|
$11 >= 70 && $5 >= 10 |
2,234 | ✅ 高置信度 6mA 位点 (推荐用于发表) |
$11 >= 50 && $5 >= 5 |
3,422 | 🔍 探索性分析 (包含中等置信度位点) |
💡 建议：对于细菌 6mA 分析，2,234 个高置信度位点是合理的数量级（取决于菌株甲基转移酶活性）。
📋 bedMethyl 格式列详解 (18 列)
您的输出行示例：
5 4 5 a 16 - 4 5 255,0,0 16 70.00 0 16 0 0 0 0 0
│ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18| 列号 | 列名 | 示例值 | 含义 | 分析建议 |
|---|---|---|---|---|
| 1 | chrom |
5 |
染色体/contig 名称 | 用于按基因组区域分组 |
| 2 | start |
4 |
起始位置 (0-based) | ⚠️ BED 格式起始为 0，转换坐标时 +1 |
| 3 | end |
5 |
结束位置 (1-based, exclusive) | 单碱基修饰时 end = start + 1 |
| 4 | name |
a |
⭐ 修饰码: a=6mA, m=5mC, c=4mC |
过滤时用 $4 == "a" 提取 6mA |
| 5 | score |
16 |
⭐ 覆盖深度 (该位点总 reads 数) | 过滤低深度: $5 >= 10 |
| 6 | strand |
- |
链方向: + / - |
6mA 分析通常合并双链 |
| 7 | thickStart |
4 |
可视化用 (通常 = start) | 可忽略 |
| 8 | thickEnd |
5 |
可视化用 (通常 = end) | 可忽略 |
| 9 | itemRgb |
255,0,0 |
颜色值 (RGB) | 可忽略 |
| 10 | blockCount |
16 |
检测到修饰的 reads 数 | 可计算修饰比例: $10/$5 |
| 11 | blockSizes |
70.00 |
⭐⭐ 修饰概率 (0-100 百分比!) | 核心过滤列: $11 >= 70 |
| 12 | starts |
|
内部偏移 (通常为 0) | 可忽略 |
| 13-18 | 扩展列 | 16 0 0 0 0 0 |
modkit 附加统计信息 | 见下方详解 |
🔍 扩展列 (13-18) 详解
| 列号 | 含义 | 示例 | 用途 | |
|---|---|---|---|---|
| 13 | mod_count |
16 |
该位点检测到修饰的 reads 数 | 与列10相同，冗余备份 |
| 14 | unmod_count |
|
该位点未修饰的 reads 数 | 可验证: $13 + $14 ≈ $5 |
| 15 | del_count |
|
该位点有 deletion 的 reads 数 | 质控: 高值可能表示比对问题 |
| 16 | no_call_count |
|
无法判断修饰状态的 reads 数 | 质控: 高值表示信号模糊 |
| 17 | coverage_plus |
|
正链覆盖度 | 链特异性分析时用 |
| 18 | coverage_minus |
|
负链覆盖度 | 链特异性分析时用 |
🔧 实用计算示例：
# 计算每个位点的"原始修饰比例" (未过滤前)
awk '{printf "%s\t%s\t%.2f%%\n", $1, $2, ($10/$5)*100}' An6_6mA_raw.bed | head -5
# 验证覆盖度一致性 (正负链之和 ≈ 总覆盖)
awk '{if($17+$18 != $5) print "Mismatch:", $0}' An6_6mA_raw.bed | wc -l
# 预期: 0 或极少 (浮点误差)🎯 常用分析命令模板
🔹 1. 提取核心信息 (简化表格)
# 输出: chrom, position(1-based), strand, probability, coverage
awk '$4 == "a" && $11 >= 70 && $5 >= 10 {
printf "%s\t%d\t%s\t%.2f\t%d\n", $1, $2+1, $6, $11, $5
}' An6_6mA_raw.bed > An6_6mA_simple.tsv
# 查看前5行
head -5 An6_6mA_simple.tsv
# 输出示例:
# 5 5 - 99.54 16
# 5 7 + 98.21 23🔹 2. 按染色体统计 6mA 密度
# 每百万碱基的 6mA 位点数 (需知道各 contig 长度)
# 先获取基因组大小
samtools faidx bacass_out/Medaka/An6-unicycler-medaka_polished_genome.fa | cut -f1,2 > genome.sizes
# 计算密度
awk '$4 == "a" && $11 >= 70 {count[$1]++}
END {for(c in count) print c, count[c]}' An6_6mA_raw.bed | \
join - genome.sizes | \
awk '{printf "%s\t%d sites\t%.2f sites/Mb\n", $1, $2, $2/($3/1000000)}' | \
sort -k3 -rn🔹 3. 提取 6mA 位点上下游序列 (用于 motif 分析)
# 提取 ±20bp 序列
bedtools slop -i <(awk '$4 == "a" && $11 >= 70 {print $1"\t"$2"\t"$3}' An6_6mA_raw.bed) \
-g genome.sizes -b 20 | \
bedtools getfasta -fi bacass_out/Medaka/An6-unicycler-medaka_polished_genome.fa -bed - \
-fo An6_6mA_motif_seqs.fa
# 运行 MEME 找富集 motif
meme An6_6mA_motif_seqs.fa -oc meme_6mA -dna -mod zoops -nmotifs 3 -minw 4 -maxw 12🔹 4. 与基因注释关联 (需 GFF/GTF 文件)
# 假设有 genome.gff 文件
bedtools intersect -a <(awk '$4 == "a" && $11 >= 70 {print $1"\t"$2"\t"$3"\t"$11}' An6_6mA_raw.bed) \
-b genome.gff -wa -wb > An6_6mA_with_genes.tsv
# 统计 6mA 在基因区域的分布
cut -f8 An6_6mA_with_genes.tsv | cut -f3 | sort | uniq -c | sort -rn
# 输出示例:
# 1234 CDS
# 567 promoter
# 234 intergenic🔹 5. 生成 IGV 可视化的 bigWig 文件
# 1. 先转为 bedGraph (chrom, start, end, score=probability)
awk '$4 == "a" {printf "%s\t%d\t%d\t%.2f\n", $1, $2, $3, $11}' An6_6mA_raw.bed > An6_6mA.bedgraph
# 2. 转为 bigWig (需 genome.sizes 文件)
bedGraphToBigWig An6_6mA.bedgraph genome.sizes An6_6mA.bw
# 3. 在 IGV 中加载:
# - File → Load Genome → 加载 genome.fa
# - File → Load from File → 加载 An6_6mA.bw⚠️ 重要注意事项
1️⃣ 坐标系统：0-based vs 1-based
BED 格式 (modkit 输出): start=4, end=5 → 表示第 5 个碱基 (0-based)
生物学常用: position = start + 1 = 5 (1-based)
# 转换命令:
awk '{print $1"\t"($2+1)"\t"$3}' An6_6mA_final.bed > An6_6mA_1based.bed2️⃣ 概率列是百分比 (0-100)，不是小数 (0-1)
# ❌ 错误: 用小数阈值
awk '$11 >= 0.7' An6_6mA_raw.bed # 几乎无结果!
# ✅ 正确: 用百分比阈值
awk '$11 >= 70' An6_6mA_raw.bed # 正确过滤3️⃣ 修饰码大小写敏感
| 修饰类型 | modkit 0.5.0 输出码 | 过滤命令 |
|---|---|---|
| 6mA | a (小写) |
$4 == "a" |
| 5mC | m (小写) |
$4 == "m" |
| 4mC | c (小写) |
$4 == "c" |
📊 您的结果质量评估
# ✅ 检查高置信度位点的概率分布
awk '$4 == "a" && $11 >= 70 && $5 >= 10 {print $11}' An6_6mA_raw.bed | \
sort -n | uniq -c | tail -10
# 预期输出示例:
# 234 95.00
# 456 98.21
# 892 99.54
# 652 100.00 ← 大量 100% 置信度，信号极佳! ✅
# ✅ 检查覆盖度分布
awk '$4 == "a" && $11 >= 70 && $5 >= 10 {print $5}' An6_6mA_raw.bed | \
sort -n | uniq -c | head -10
# 预期: 大部分位点覆盖度 10-50×，说明测序深度充足 ✅💡 一句话总结
🔑 核心三列：
$4= 修饰码 (a=6mA)$5= 覆盖深度 (过滤低质量:>=10)$11= 修饰概率 (0-100 百分比! 过滤高置信:>=70)✅ 您的 2,234 个高置信 6mA 位点是高质量的起点，可用于：
- Motif 富集分析 (验证已知限制修饰系统)
- 基因组分布分析 (启动子/基因体/间隔区)
- 与转录组数据整合 (如有)
如需帮助进行下游 motif 分析或可视化，随时告诉我！🚀
🔬 从 Tn-seq (TRANSIT) 到甲基化分析：下游统计工具对比指南
您提出的问题非常精准！正如 Tn-seq 用 TRANSIT 进行插入位点的统计分析（ANOVA、Tn5Gaps），甲基化测序（如 modkit pileup 输出）也有对应的下游统计工具。以下是详细对比和推荐：
📊 Tn-seq vs 甲基化分析：核心概念映射
| Tn-seq 分析目标 | 甲基化分析对应目标 | 统计本质 |
|---|---|---|
| 插入位点计数 (template counts) | 修饰概率/计数 (mod_count/coverage) | 计数型数据 |
| 基因必需性 (essential genes) | 差异甲基化区域 (DMRs) | 区域水平差异检测 |
| 条件间比较 (ANOVA in TRANSIT) | 条件间甲基化差异 (DSS, methylKit) | 多组比较 + 多重校正 |
| 非插入区域检测 (Tn5Gaps) | 低甲基化/高甲基化区域 (DMRfinder) | 连续区域模式识别 |
| 密度标准化 (TTR normalization) | 覆盖度/文库大小标准化 (TMM, RLE) | 文库效应校正 |
💡 关键区别：Tn-seq 关注”有没有插入”（二元事件），甲基化关注”修饰比例是多少”（连续概率 0-100%）。
🧰 甲基化下游统计工具推荐（按功能分类）
🔹 1. 差异甲基化分析 (对应 TRANSIT ANOVA)
| 工具 | 适用场景 | 输入格式 | 核心方法 | 细菌 6mA 兼容性 |
|---|---|---|---|---|
| DSS (Dispersion Shrinkage for Sequencing) | ✅ 多样本组间比较 | bed/bedGraph + 样本分组 | 贝叶斯层次模型 + 算术-贝塔分布 | ✅ 优秀 (支持任意 motif) |
| methylKit | ✅ 多条件/时间序列 | bedGraph/tabix | 逻辑回归 + FDR 校正 | ✅ 良好 (需格式转换) |
| metilene | ✅ 快速 DMR 检测 | bedGraph | 二元分割 + Kolmogorov-Smirnov 检验 | ⚠️ 需确认 motif 支持 |
| DMRseq | ✅ 区域水平 + 协变量调整 | bed + 样本信息 | 平滑 + 置换检验 | ✅ 良好 |
🚀 DSS 快速入门 (推荐首选)
# R 代码示例：细菌 6mA 差异分析
library(DSS)
# 1. 读取 modkit 输出 (假设已过滤高置信位点)
# 文件格式: chrom, start, end, motif, coverage, mod_count, probability
sample1 <- read.table("An6_6mA_final.bed", header=F)
sample2 <- read.table("BG5_6mA_final.bed", header=F)
# 2. 构建 DSS 对象 (使用原始计数: mod_count / coverage)
BSobj <- makeBSseqData(
list(sample1[, c(1,2,3,6,5)], # chrom, start, end, mod_count, coverage
sample2[, c(1,2,3,6,5)]),
sampleNames = c("An6", "BG5")
)
# 3. 差异检验 (贝塔-二项模型)
dmlTest <- DMLtest(BSobj, group1=1, group2=2, smoothing=TRUE)
# 4. 提取显著差异位点 (FDR < 0.05, |diff| > 0.2)
sigDML <- getSig(dmlTest, pval=0.05, delta=0.2, meth.diff=0.2)
# 5. 导出结果
write.table(sigDML, "An6_vs_BG5_6mA_DMLs.tsv", sep="\t", quote=F, row.names=F)📋 DSS 输出列解读 (类比 TRANSIT)
| DSS 输出列 | 含义 | TRANSIT 对应列 |
|---|---|---|
chr, pos |
基因组坐标 | Orf, TA_pos |
coverage |
总覆盖深度 | TAs (TA 位点数) |
X1, X2 |
两组的修饰计数 | Mean_[condition] |
pval, p.adj |
原始/FDR 校正 p 值 | Pval, Padj |
diff |
甲基化差异 (组1 – 组2) | LFC_[condition] (logFC) |
abs.diff |
绝对差异值 | – |
🔹 2. 区域水平分析 (对应 Tn5Gaps: 连续非修饰区域检测)
| 工具 | 功能 | 适用场景 |
|---|---|---|
| DMRfinder | 识别连续低/高甲基化区域 | 类似 Tn5Gaps 的”非插入运行”检测 |
| methylSig | 区域聚合 + 差异检验 | 基因/启动子水平的甲基化变化 |
| custom R/Python | 自定义滑动窗口分析 | 细菌操纵子/基因簇水平分析 |
🚀 滑动窗口分析示例 (Python)
# 检测连续低甲基化区域 (类似 Tn5Gaps 的 non-insertion runs)
import pandas as pd
import numpy as np
# 读取 modkit 输出
df = pd.read_csv("An6_6mA_final.bed", sep="\t", header=None,
names=["chr","start","end","motif","cov","strand",
"thickS","thickE","rgb","mod_cnt","prob",
"start2","cnt1","cnt0","del_cnt","no_call","cov+","cov-"])
# 按染色体分组，检测连续低甲基化区域 (prob < 30%)
def find_hypo_regions(group, min_len=5, max_gap=100, prob_thresh=30):
low_meth = group[group["prob"] < prob_thresh].sort_values("start")
if len(low_meth) < min_len:
return None
# 简单连续区域检测 (可扩展为更复杂算法)
regions = []
curr_start, curr_end = low_meth.iloc[0]["start"], low_meth.iloc[0]["end"]
for _, row in low_meth.iloc[1:].iterrows():
if row["start"] - curr_end <= max_gap: # 允许小间隙
curr_end = row["end"]
else:
if curr_end - curr_start >= min_len * 10: # 最小区域长度
regions.append((curr_start, curr_end))
curr_start, curr_end = row["start"], row["end"]
return regions
# 应用并输出
hypo_regions = df.groupby("chr").apply(find_hypo_regions).dropna()
hypo_regions.to_csv("An6_hypomethylated_regions.bed", sep="\t", header=False, index=False)🔹 3. 功能注释与可视化 (对应 Circos 基因组视图)
| 工具 | 功能 | 输出示例 |
|---|---|---|
| ChIPseeker (R) | 甲基化位点与基因/启动子关联 | 饼图: CDS/promoter/intergenic 分布 |
| deepTools | 甲基化信号热图/轮廓图 | 基因体 ±2kb 的 6mA 密度曲线 |
| pyGenomeTracks | 多轨道基因组浏览器式可视化 | 类似 Circos 的线性基因组视图 |
| IGV + bigWig | 交互式位点查看 | 手动验证关键基因修饰模式 |
🚀 基因区域关联分析 (R + ChIPseeker)
library(ChIPseeker)
library(TxDb.Ecoli.NCBI.NC_008800.1.refseq) # 替换为您的细菌注释
# 读取高置信 6mA 位点 (转换为 GRanges)
library(GenomicRanges)
gr <- GRanges(
seqnames = df$chr,
ranges = IRanges(start = df$start + 1, end = df$end), # BED 0-based → 1-based
score = df$prob
)
# 注释到基因区域
anno <- annotatePeak(gr, tssRegion=c(-300, 300),
TxDb = TxDb.Ecoli.NCBI.NC_008800.1.refseq,
annoDb = "org.EcK12.eg.db") # 或用 custom GFF
# 可视化
plotAnnoPie(anno) # 6mA 在基因组区域的分布
plotAvgProf(anno) # TSS 附近的 6mA 密度曲线🎯 针对您细菌 6mA 项目的推荐流程
graph LR
A[modkit pileup 输出] --> B[质量控制 + 过滤]
B --> C[差异甲基化分析]
B --> D[区域模式检测]
C --> C1[DSS: 位点水平差异]
C --> C2[methylKit: 基因水平聚合]
D --> D1[滑动窗口: 连续低/高甲基化区]
D --> D2[与操纵子/基因簇关联]
C1 & C2 & D1 & D2 --> E[功能注释]
E --> F[可视化 + 生物学解释]✅ 具体命令/代码整合
# 🔹 步骤 1: 准备输入文件 (modkit → DSS 格式)
# 提取核心列: chrom, start, end, mod_count, coverage
awk '$4 == "a" && $11 >= 70 && $5 >= 10 {
printf "%s\t%d\t%d\t%d\t%d\n", $1, $2, $3, $10, $5
}' An6_6mA_raw.bed > An6_6mA_dss_input.tsv
# 🔹 步骤 2: 运行 DSS (R 脚本)
Rscript run_dss_6mA.R \
--input An6_6mA_dss_input.tsv,BG5_6mA_dss_input.tsv \
--groups An6,An6,An6,BG5,BG5,BG5 \
--output An6_vs_BG5_6mA_DSS_results.tsv
# 🔹 步骤 3: 提取显著差异位点 + 注释
# (见上方 R 代码示例)
# 🔹 步骤 4: 可视化
# - IGV: 加载 bigWig 轨道查看关键基因
# - pyGenomeTracks: 生成出版级多轨道图⚠️ 细菌甲基化分析的特殊注意事项
| 问题 | 解决方案 |
|---|---|
| 基因组小、基因密集 | 使用基因/操纵子为单位聚合，而非固定窗口 |
| 6mA 富集于特定 motif (如 GATC) | 分析时按 motif 分层，避免背景噪声 |
| 无”基因长度”标准化需求 | 细菌基因长度变异小，可直接用计数/概率 |
| 限制修饰系统影响 | 将已知 RM 系统基因作为阳性对照验证流程 |
| 样本量小 (常为 2-3 重复) | 优先用 DSS (贝叶斯收缩) 而非纯频率比较 |
📋 工具选择决策树
您的分析目标？
│
├─ 🔍 "哪些位点在条件 A vs B 甲基化水平不同？"
│ └─→ ✅ DSS (位点水平) 或 methylKit (基因水平)
│
├─ 🧬 "哪些基因/区域整体甲基化状态改变？"
│ └─→ ✅ methylKit + ChIPseeker 注释
│
├─ 🎯 "哪些区域连续低/高甲基化 (类似 Tn5Gaps)？"
│ └─→ ✅ 自定义滑动窗口 + DMRfinder
│
├─ 📊 "全局甲基化模式可视化？"
│ └─→ ✅ pyGenomeTracks / IGV + bigWig
│
└─ 🧪 "验证已知限制修饰系统活性？"
└─→ ✅ 提取 RM 系统基因位点，手动检查概率分布💡 一句话总结
🔑 核心推荐：
- 差异分析 →
DSS(最接近 TRANSIT ANOVA 的统计严谨性)- 区域模式 →
自定义滑动窗口(最灵活适配细菌基因组)- 功能注释 →
ChIPseeker + pyGenomeTracks(最直观可视化)🎯 您的 2,234 个高置信 6mA 位点是优质输入，结合菌株分组信息，即可开展与 Tn-seq 同等严谨的统计推断！
如需具体 R/Python 脚本模板或帮助解读 DSS 输出，请随时告知！🚀