🧰 甲基化下游统计分析：独立工具推荐与使用指南

您不需要 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.bed

2️⃣ 概率列是百分比 (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 输出，请随时告知！🚀

🧰 甲基化下游统计分析：独立工具推荐与使用指南

📊 工具总览：按 TRANSIT 功能映射

🔧 工具详解：安装 + 使用示例

1️⃣ DSS (Dispersion Shrinkage for Sequencing) — 首选推荐 ⭐

🔹 安装

🔹 输入准备 (modkit bed → DSS 格式)

🔹 R 分析脚本 (run_dss_6mA.R)

🔹 运行

🔹 输出解读 (对标 TRANSIT ANOVA)

2️⃣ metilene — 快速 DMR 检测 (对标 Tn5Gaps)

🔹 安装

🔹 输入准备

🔹 运行 DMR 检测

🔹 输出解读

3️⃣ methylKit — 基因/区域水平聚合分析

🔹 安装

🔹 R 脚本示例 (methylkit_gene_level.R)

4️⃣ 可视化工具：pyGenomeTracks (对标 Circos)

🔹 安装

🔹 准备配置文件 (tracks.ini)

🔹 生成 bigWig (modkit bed → bigWig)

🔹 运行可视化

🎯 针对您项目的最小可行流程

📦 一键安装脚本 (conda 环境)

💡 工具选择决策树

🎯 一句话总结

🔍 nf-core/methylong 内置 DMR 选项：能替代哪些独立工具？

✅ 简短回答

📋 methylong DMR 选项详解 (对标 TRANSIT)

🔹 参数功能映射表

🚀 用 methylong 内置 DSS 替代独立运行的完整命令

场景：细菌 6mA，An6 vs BG5 差异分析

🔹 samplesheet_6mA.csv 格式要求

📁 输出文件：methylong 内置 DMR 结果位置

🔹 DSS_DMLs.tsv 列解读 (对标 TRANSIT ANOVA 输出)

⚠️ methylong 内置 DMR 的局限性 (仍需外部工具的场景)

❌ 无法替代的功能 1：连续区域检测 (对标 Tn5Gaps)

❌ 无法替代的功能 2：基因/操纵子水平注释

❌ 无法替代的功能 3：高级可视化

🎯 推荐工作流：methylong 内置 + 最小外部工具

✅ 最小命令集：

💡 一句话总结

📊 modkit pileup 输出文件列详解 (bedMethyl 格式)

🎯 您的结果统计

📋 bedMethyl 格式列详解 (18 列)

🔍 扩展列 (13-18) 详解

🔧 实用计算示例：

🎯 常用分析命令模板

🔹 1. 提取核心信息 (简化表格)

🔹 2. 按染色体统计 6mA 密度

🔹 3. 提取 6mA 位点上下游序列 (用于 motif 分析)

🔹 4. 与基因注释关联 (需 GFF/GTF 文件)

🔹 5. 生成 IGV 可视化的 bigWig 文件

⚠️ 重要注意事项

1️⃣ 坐标系统：0-based vs 1-based

2️⃣ 概率列是百分比 (0-100)，不是小数 (0-1)

3️⃣ 修饰码大小写敏感

📊 您的结果质量评估

💡 一句话总结

🔬 从 Tn-seq (TRANSIT) 到甲基化分析：下游统计工具对比指南

📊 Tn-seq vs 甲基化分析：核心概念映射

🧰 甲基化下游统计工具推荐（按功能分类）

🔹 1. 差异甲基化分析 (对应 TRANSIT ANOVA)

🚀 DSS 快速入门 (推荐首选)

📋 DSS 输出列解读 (类比 TRANSIT)

🔹 2. 区域水平分析 (对应 Tn5Gaps: 连续非修饰区域检测)

🚀 滑动窗口分析示例 (Python)

🔹 3. 功能注释与可视化 (对应 Circos 基因组视图)

🚀 基因区域关联分析 (R + ChIPseeker)

🎯 针对您细菌 6mA 项目的推荐流程

✅ 具体命令/代码整合

⚠️ 细菌甲基化分析的特殊注意事项

📋 工具选择决策树

💡 一句话总结

🔹 R 分析脚本 (`run_dss_6mA.R`)

🔹 R 脚本示例 (`methylkit_gene_level.R`)

🔹 准备配置文件 (`tracks.ini`)