❌ 不建议将 6mA 和 4mC 混在同一张 Samplesheet 中
您的当前格式会导致 Nextflow 报错或分析结果混乱，主要原因如下：
| 问题 | 原因 | 后果 |
|---|---|---|
sample 名称重复 |
An6 和 BG5 各出现两次 |
Nextflow 会尝试覆盖文件或报错退出 |
| 混合修饰类型 | 6mA 和 4mC_5mC 的 modBAM 混跑 | modkit 甲基化 pileup 会冲突，DMR 分析失败 |
group 命名无意义 |
全设为 bacteria |
DMR 比较需要至少两个不同组，单组无法进行差异分析 |
✅ 正确做法：拆分为两个独立的 Samplesheet
为每种修饰类型创建独立的运行流程，这是 nf-core/methylong 的标准实践：
📄 samplesheet_6mA.csv
group,sample,path,ref,method
An6,An6,/mnt/md1/DATA/Data_Tam_DNAseq_2026_An6_BG5/An6_6mA.mod.bam,/mnt/md1/DATA/Data_Tam_DNAseq_2026_An6_BG5/An6_unicycler_normal/assembly.fasta,ont
BG5,BG5,/mnt/md1/DATA/Data_Tam_DNAseq_2026_An6_BG5/BG5_6mA.mod.bam,/mnt/md1/DATA/Data_Tam_DNAseq_2026_An6_BG5/BG5_unicycler_normal/assembly.fasta,ont📄 samplesheet_4mC.csv
group,sample,path,ref,method
An6,An6,/mnt/md1/DATA/Data_Tam_DNAseq_2026_An6_BG5/An6_4mC_5mC.mod.bam,/mnt/md1/DATA/Data_Tam_DNAseq_2026_An6_BG5/An6_unicycler_normal/assembly.fasta,ont
BG5,BG5,/mnt/md1/DATA/Data_Tam_DNAseq_2026_An6_BG5/BG5_4mC_5mC.mod.bam,/mnt/md1/DATA/Data_Tam_DNAseq_2026_An6_BG5/BG5_unicycler_normal/assembly.fasta,ont💡 为什么
group和sample同名？
sample必须是唯一标识符（不能重复）group用于 DMR（差异甲基化区域）分析。如果您没有生物学重复（如 3 个 An6 培养皿），DMR 实际上会跳过。此时group设为与sample相同最安全，避免 pipeline 误判分组。
🚀 运行命令（分两次执行）
1️⃣ 运行 6mA 分析
nextflow run nf-core/methylong \
-r 2.0.0 \
-profile docker \
--input samplesheet_6mA.csv \
--outdir methylome_out_6mA \
--skip_basecalling \ # ⭐ 使用预生成的 modBAM
-resume \
-work-dir methylome_out_6mA/work2️⃣ 运行 4mC 分析
nextflow run nf-core/methylong \
-r 2.0.0 \
-profile docker \
--input samplesheet_4mC.csv \
--outdir methylome_out_4mC \
--skip_basecalling \ # ⭐ 使用预生成的 modBAM
-resume \
-work-dir methylome_out_4mC/work📊 关于 group 列的深入说明
| 您的研究设计 | 推荐 group 命名 |
DMR 是否有效 |
|---|---|---|
| 无重复，仅单样本甲基化图谱 | group=sample (如 An6, BG5) |
❌ 跳过（正常） |
| 有生物学重复 (如 An6_rep1, An6_rep2) | group=An6, group=BG5 |
✅ 可比较组间差异 |
| 时间序列/处理组 (如 Control, Treated) | group=control, group=treated |
✅ 推荐设置 |
🔑 结论：如果您只是为每个菌株生成甲基化基序和覆盖度图谱，保持
group=sample即可。后续分析主要依赖methylation_calls/和motifs/目录，DMR 不是必需步骤。
🇩🇪 德语快速参考（Deutsch）
# ❌ 错误: 同一 sheet 中 sample 重复 + 混合修饰类型
# ✅ 正确: 拆分为 samplesheet_6mA.csv 和 samplesheet_4mC.csv
# 📋 group 命名建议:
• 无重复样本 → group = sample (An6, BG5)
• 有重复样本 → group = 生物学条件 (control, mutant)
• DMR 需 ≥2 个重复/组，否则 pipeline 会自动跳过
# 🚀 运行:
nextflow run nf-core/methylong \
--input samplesheet_6mA.csv \
--skip_basecalling -resume ...
nextflow run nf-core/methylong \
--input samplesheet_4mC.csv \
--skip_basecalling -resume ...✅ 快速检查清单
[ ] ✅ samplesheet_6mA.csv 中 sample 列无重复
[ ] ✅ samplesheet_4mC.csv 中 sample 列无重复
[ ] ✅ path 列使用绝对路径且文件存在
[ ] ✅ ref 列指向对应的组装 FASTA
[ ] ✅ 两次运行均添加 --skip_basecalling
[ ] ✅ 输出目录分开 (methylome_out_6mA / methylome_out_4mC)直接创建上述两个 CSV 文件，然后依次运行。这样可确保 modkit 正确解析修饰标签，且 Nextflow 不会因重复样本名崩溃。
运行完成后，我会帮您提取 6mA/4mC 的基序富集结果和基因组分布统计！🦠🧬✨
❌ 错误原因明确：您的 Dorado 版本不支持 5mC_6mA_4mC 组合参数
从报错信息可以直接看到：
[error] '5mC_6mA_4mC' is not a supported modification
please select from: pseU, inosine_m6A, 5mCG, pseU_2OmeU,
inosine_m6A_2OmeA, 5mCG_5hmCG, **4mC_5mC**, 2OmeG, 5mC_5hmC,
5mC, m5C, m5C_2OmeC, **6mA**, m6A_DRACH, m6A结论：当前 Dorado 版本将修饰检测拆分为独立参数。5mC_6mA_4mC 不能作为一个字符串传入。
✅ 立即修正命令（推荐用于细菌）
对于 Acinetobacter 和 Pedobacter，6mA 是最核心、最相关的修饰类型（主导限制修饰系统）。直接使用它：
dorado basecaller \
--modified-bases 6mA \
--emit-moves \
--device cuda:0 \
dna_r10.4.1_e8.2_400bps_sup@v5.0.0 \
./X101SC26036392-Z01-J003/Release-X101SC26036392-Z01-J003-20260513_01/Data-X101SC26036392-Z01-J003/An6/2157_4C_PBK79106_7ec05c46/An6_pod5_pass \
> An6_6mA.mod.bam🔍 如果您确实需要检测 4mC：
使用独立参数 4mC_5mC（会同时输出 4mC 和 5mC，后续可过滤）：
dorado basecaller \
--modified-bases 4mC_5mC \
--emit-moves \
--device cuda:0 \
dna_r10.4.1_e8.2_400bps_sup@v5.0.0 \
./X101SC26036392-Z01-J003/Release-X101SC26036392-Z01-J003-20260513_01/Data-X101SC26036392-Z01-J003/An6/2157_4C_PBK79106_7ec05c46/An6_pod5_pass \
> An6_4mC.mod.bam📊 为什么单独运行 6mA 是最佳选择？
| 考量 | 说明 |
|---|---|
| 生物学相关性 | Acinetobacter 和 Pedobacter 的 R-M 系统主要依赖 6mA (如 Dam/Dcm 同源酶) |
| 计算效率 | 单次运行比拆分两次快 50%，显存占用更低 |
| 分析兼容性 | nf-core/methylong 和 modkit 对单一修饰模型支持最完善 |
| 数据质量 | 组合模型在旧版 Dorado 中易产生假阳性；独立模型更稳健 |
💡 建议流程：先完成
6mA分析。如果后续发现启动子/限制位点区域有未解释的甲基化信号，再补跑4mC_5mC。
🔄 如何与 nf-core/methylong 对接
如果您已手动生成 An6_6mA.mod.bam，直接更新 samplesheet 即可：
group,sample,path,ref,method
bacteria,An6,/absolute/path/to/An6_6mA.mod.bam,/path/to/An6_assembly.fasta,ont
bacteria,BG5,/absolute/path/to/BG5_6mA.mod.bam,/path/to/BG5_assembly.fasta,ont然后运行（跳过内部 basecalling）：
nextflow run nf-core/methylong \
-r 2.0.0 -profile docker \
--input samplesheet_methylong.csv \
--outdir methylome_out_6mA \
--skip_basecalling \ # ⭐ 使用您手动生成的 modBAM
-resume✅ 验证生成结果是否成功
# 1. 检查 BAM 是否包含甲基化标签
samtools view An6_6mA.mod.bam | head -50 | grep -o "MM:Z:[^ ]*" | sort -u
# ✅ 应输出: MM:Z:6mA,... 或类似
# 2. 检查文件完整性
samtools quickcheck An6_6mA.mod.bam && echo "✅ BAM 结构正常" || echo "❌ 文件损坏"
# 3. 统计调用量
samtools view -c An6_6mA.mod.bam🇩🇪 德语快速参考
# ❌ 错误: 5mC_6mA_4mC 不被当前版本支持
# ✅ 修正: 使用独立参数 6mA 或 4mC_5mC
# 🔧 推荐命令 (细菌 6mA):
dorado basecaller \
--modified-bases 6mA \
--emit-moves \
--device cuda:0 \
dna_r10.4.1_e8.2_400bps_sup@v5.0.0 \
./path/to/pod5_pass > sample_6mA.mod.bam
# ✅ 验证:
samtools view sample_6mA.mod.bam | grep "MM:Z:6mA" | head -3请直接复制上方的 6mA 命令运行。它将在您的 RTX 4060 Ti 上正常调用，并生成可直接用于下游分析的 modBAM。
如需批量处理 BG5 或验证输出，随时告诉我！🦠🧬✨
🔧 CUDA 版本不匹配错误分析与修复
关键错误信息：
nvidia-container-cli: requirement error: unsatisfied condition: cuda>=12.8,
please update your driver to a newer version, or use an earlier cuda container| 问题原因： | 组件 | 您的系统 | Dorado 容器要求 | 状态 |
|---|---|---|---|---|
| NVIDIA 驱动 | 535.288.01 | ≥545.23.06 (for CUDA 12.8) | ❌ 过旧 | |
| CUDA 运行时 | 12.2 | ≥12.8 | ❌ 不匹配 | |
| Dorado 容器 | shae423e761540b9d08b526a1eb32faf498f32e8f22 |
需要 CUDA 12.8+ | ❌ 不兼容 |
💡 Exit status 125 = Docker 容器启动失败（通常是镜像/驱动不兼容）
✅ 解决方案（按推荐优先级）
🔹 方案 1: 使用兼容 CUDA 12.2 的 Dorado 容器标签（最快 ✅）
nf-core/methylong v2.0.0 默认使用较新的 Dorado 镜像，但您可以覆盖容器镜像使用旧版：
// nextflow.config - 添加容器覆盖
process {
withName: "*DORADO*" {
containerOptions = '--gpus all'
memory = '24.GB'
cpus = 8
time = '72.h'
errorStrategy = 'retry'
// ⭐ 关键: 使用兼容 CUDA 12.2 的 Dorado 镜像
container = 'nanoporetech/dorado:7.1.2' // 或 7.0.0, 6.4.0
}
withName: "*METHYLONG:ONT*" {
containerOptions = '--gpus all'
memory = '16.GB'
cpus = 4
}
}
docker {
enabled = true
runOptions = '-u $(id -u):$(id -g) --gpus all'
}| 可用兼容镜像标签（按推荐顺序）： | 镜像标签 | CUDA 要求 | Dorado 版本 | 兼容性 |
|---|---|---|---|---|
nanoporetech/dorado:7.1.2 |
≥12.2 | v7.1.2 | ✅ 推荐 | |
nanoporetech/dorado:7.0.0 |
≥12.0 | v7.0.0 | ✅ 稳定 | |
nanoporetech/dorado:6.4.0 |
≥11.8 | v6.4.0 | ✅ 最兼容 |
🔍 验证镜像兼容性：
docker pull nanoporetech/dorado:7.1.2 docker run --rm --gpus all nanoporetech/dorado:7.1.2 dorado --version # 应输出: dorado 7.1.2
🔹 方案 2: 升级 NVIDIA 驱动以支持 CUDA 12.8+（长期方案）
# 1. 检查可用驱动版本
ubuntu-drivers devices
# 2. 安装最新驱动 (545+)
sudo apt update
sudo apt install -y nvidia-driver-545 # 或 550
# 3. 重启系统
sudo reboot
# 4. 验证
nvidia-smi
# 应显示: Driver Version: 545.xx | CUDA Version: 12.8+⚠️ 注意：驱动升级可能需要重启，且可能影响其他 GPU 应用。
🔹 方案 3: 本地预生成 modBAM + 传给 methylong（绕过容器 CUDA 限制）
如果您已安装本地 Dorado（非容器版），可预先生成 modBAM：
# 1. 安装本地 Dorado (CPU 或匹配您 CUDA 的版本)
# 下载: https://github.com/nanoporetech/dorado/releases
# 2. 预生成 modBAM
dorado basecaller \
--modified-bases 6mA \
--emit-moves \
--device cuda:0 \
dna_r10.4.1_e8.2_400bps_hac@v4.2.0 \ # 匹配本地 Dorado 版本
/path/to/An6_pod5_pass/ \
> An6.mod.bam
# 3. 更新 samplesheet 使用 modBAM 路径（而非 POD5）
# samplesheet_methylong.csv:
group,sample,path,ref,method
bacteria,An6,/path/to/An6.mod.bam,/path/to/ref.fasta,ont
# 4. 运行 methylong（跳过 basecalling）
nextflow run nf-core/methylong \
-r 2.0.0 -profile docker \
--input samplesheet_methylong.csv \
--outdir methylome_out_6mA_hac \
--dorado_model hac \
--dorado_modification 6mA \
--skip_basecalling \ # ⭐ 关键: 跳过内部 Dorado
-c nextflow.config \
-resume✅ 优势：完全绕过容器 CUDA 限制，使用您本地兼容的 Dorado 版本
🔹 方案 4: 临时使用 CPU 模式（仅测试，不推荐生产）
// nextflow.config - CPU fallback
process {
withName: "*DORADO*" {
// 移除 --gpus all，使用 CPU
containerOptions = ''
memory = '48.GB' // CPU 模式需要更多内存
cpus = 16
time = '168.h' // 预期运行时间很长
}
}⚠️ 警告：Dorado basecalling on CPU 可能需 10-50× 更长时间，仅用于调试。
🚀 推荐：方案 1 实施步骤（最快修复）
步骤 1: 更新 nextflow.config
cat > nextflow.config << 'EOF'
// nextflow.config - 兼容 CUDA 12.2 的配置
process {
withName: "*DORADO*" {
containerOptions = '--gpus all'
memory = '24.GB'
cpus = 8
time = '72.h'
errorStrategy = 'retry'
maxRetries = 2
// ⭐ 使用兼容 CUDA 12.2 的 Dorado 镜像
container = 'nanoporetech/dorado:7.1.2'
}
withName: "*METHYLONG:ONT*" {
containerOptions = '--gpus all'
memory = '16.GB'
cpus = 4
}
}
docker {
enabled = true
runOptions = '-u $(id -u):$(id -g) --gpus all'
}
EOF
echo "✅ nextflow.config 已更新 (Dorado 7.1.2)"步骤 2: 预拉取镜像（避免运行中下载失败）
docker pull nanoporetech/dorado:7.1.2
docker run --rm --gpus all nanoporetech/dorado:7.1.2 dorado --version
# 应输出: dorado 7.1.2步骤 3: 清理旧工作目录 + 重新运行
# 清理之前失败的缓存
rm -rf methylome_out_6mA_hac/work/
# 运行 (HAC 模型 + 兼容镜像)
nextflow run nf-core/methylong \
-r 2.0.0 \
-profile docker \
--input samplesheet_methylong.csv \
--outdir methylome_out_6mA_hac \
--dorado_model hac \
--dorado_modification 6mA \
-c nextflow.config \
-resume \
-work-dir methylome_out_6mA_hac/work \
-log methylome_out_6mA_hac/logs/nextflow.log🔍 关于末尾的奇怪错误
--dorado_modification: command not found
-resume: command not found这通常是因为：
- 命令被 shell 错误解析（如复制粘贴时格式问题）
- 日志输出被截断/混淆
✅ 确保：
- 使用
bash script.sh运行，而非直接粘贴多行命令 - 或使用
\续行符时确保每行末尾无多余空格 - 或保存为脚本文件执行
🇩🇪 德语总结（Deutsch）
# ❌ 错误: cuda>=12.8 要求不满足
# ✅ 解决方案:
# 方案 1 (推荐): 使用兼容镜像
// nextflow.config:
process {
withName: "*DORADO*" {
container = 'nanoporetech/dorado:7.1.2' // CUDA 12.2 兼容
}
}
# 方案 2: 升级驱动
sudo apt install nvidia-driver-545 && sudo reboot
# 方案 3: 本地预生成 modBAM
dorado basecaller ... > An6.mod.bam
# 然后: --skip_basecalling + modBAM path in samplesheet
# 方案 4 (仅测试): CPU 模式
// 移除 --gpus all, 增加 memory/cpus/time📋 快速验证清单
[ ] ✅ 更新 nextflow.config: container = 'nanoporetech/dorado:7.1.2'
[ ] ✅ docker pull nanoporetech/dorado:7.1.2 成功
[ ] ✅ docker run --gpus all nanoporetech/dorado:7.1.2 dorado --version 输出版本
[ ] ✅ 清理旧 work 目录: rm -rf methylome_out_6mA_hac/work/
[ ] ✅ 使用脚本文件运行（避免 shell 解析错误）
[ ] ✅ 使用 --dorado_model hac (16GB VRAM 安全)
# 如果全部 ✅ → methylong 应能正常启动! 🎉⚠️ 如果仍失败：诊断命令
# 1. 检查容器实际使用的镜像
cat methylome_out_6mA_hac/work/*/.command.run | grep "docker run"
# 2. 手动测试 Dorado 7.1.2 容器
docker run --rm --gpus all \
-v $(pwd):/data -w /data \
nanoporetech/dorado:7.1.2 \
dorado basecaller --help | head -20
# 3. 检查系统 CUDA 兼容性
nvidia-smi --query-gpu=driver_version,cuda_version --format=csv
docker run --rm --gpus all nvidia/cuda:12.2.0-base-ubuntu22.04 nvcc --version 2>&1 | head -5请先尝试方案 1（使用 dorado:7.1.2 镜像），这是最快绕过 CUDA 版本限制的方法。
如果仍有问题，请分享：
docker run --rm --gpus all nanoporetech/dorado:7.1.2 dorado --version的输出.nextflow.log中的完整错误堆栈
📚 Nextflow methylong 命令参数详解（中文版）
nextflow run nf-core/methylong \
-r 2.0.0 \
-profile docker \
--input samplesheet_methylong.csv \
--outdir methylome_out \
--dorado_model sup \
--dorado_modification 5mC_6mA_4mC \
--dorado_device cuda:0 \
-resume \
-work-dir methylome_out/work🔧 核心参数详解
| 参数 | 含义 | 您的设置 | 说明 |
|---|---|---|---|
nextflow run nf-core/methylong |
运行 nf-core/methylong 流程 | — | 启动甲基化分析流程 |
-r 2.0.0 |
Revision（版本） | 2.0.0 |
指定使用流程的 2.0.0 版本，确保结果可重复 |
-profile docker |
执行环境 | docker |
使用 Docker 容器运行所有步骤，避免依赖冲突 |
--input samplesheet_methylong.csv |
输入文件 | CSV 格式 samplesheet | 定义样本、数据路径、参考基因组等信息 |
--outdir methylome_out |
输出目录 | methylome_out/ |
所有结果文件保存到此目录 |
--dorado_model sup |
碱基调用模型 | sup |
使用 Super Accuracy 模型，精度最高但速度较慢 |
--dorado_modification 5mC_6mA_4mC |
检测的修饰类型 | 5mC_6mA_4mC |
同时检测 5mC、6mA、4mC 三种细菌常见甲基化 |
--dorado_device cuda:0 |
计算设备 | cuda:0 |
使用第一块 NVIDIA GPU 加速 Dorado basecalling |
-resume |
断点续跑 | 启用 | 如果中途失败，下次运行可跳过已完成的步骤 |
-work-dir methylome_out/work |
临时工作目录 | methylome_out/work/ |
存储中间文件，方便调试和清理 |
📋 参数分类详解
1️⃣ 流程控制参数
-r 2.0.0（Revision / 版本）
# 作用：锁定流程版本
-r 2.0.0- ✅ 为什么重要：确保不同时间运行的结果完全一致（可重复性）
- 🔍 查看可用版本：
nextflow info nf-core/methylong - ⚠️ 注意：不要使用
-r master，因为主分支可能随时更新
-profile docker（执行环境）
# 作用：指定软件运行环境
-profile docker| Profile | 说明 | 适用场景 |
|---|---|---|
docker |
使用 Docker 容器 | ✅ 推荐，隔离性好 |
singularity |
使用 Singularity 容器 | HPC 集群常用 |
conda |
使用 Conda 环境 | 本地开发测试 |
test |
使用测试数据快速验证 | 流程调试 |
-resume（断点续跑）
# 作用：跳过已完成的步骤
-resume- ✅ 工作原理：Nextflow 会检查
work/目录中的哈希值，如果输入和参数未变，直接复用结果 -
💡 使用技巧：
# 首次运行 nextflow run ... -resume # 修改参数后重新运行（只重跑受影响步骤） nextflow run ... --dorado_model hac -resume # 强制从头运行（清除缓存） nextflow run ... -resume -clean
2️⃣ 输入/输出参数
--input samplesheet_methylong.csv
# samplesheet_methylong.csv 格式
group,sample,path,ref,method
bacteria,An6,/path/to/An6_pod5_pass,/path/to/An6_assembly.fasta,ont
bacteria,BG5,/path/to/BG5_pod5_pass,/path/to/BG5_assembly.fasta,ont| 列名 | 含义 | 您的值示例 |
|---|---|---|
group |
样本分组（用于差异分析） | bacteria |
sample |
样本唯一标识 | An6, BG5 |
path |
输入数据路径（POD5 目录 或 modBAM 文件） | /.../An6_pod5_pass |
ref |
参考基因组 FASTA 文件 | /.../assembly.fasta |
method |
测序平台 | ont（Nanopore） |
🔑 关键：
path指向 POD5 目录时，流程会自动调用 Dorado进行 basecalling + 甲基化检测。
--outdir methylome_out
# 作用：指定结果输出目录
--outdir methylome_out输出目录结构预览：
methylome_out/
├── multiqc_report.html # 📊 综合质量报告
├── methylation_calls/ # 🧬 甲基化位点
│ ├── An6.bedMethyl.gz # 标准 bedMethyl 格式
│ ├── An6.bigWig # IGV 可视化轨道
│ └── BG5.* # 同上
├── motifs/ # 🎯 修饰基序分析
│ ├── An6_5mC_motifs.tsv
│ ├── An6_6mA_motifs.tsv
│ └── An6_4mC_motifs.tsv
├── alignment/ # 📐 比对结果
│ ├── An6.mod.bam # 含修饰标签的 BAM
│ └── *.bai # 索引文件
├── pipeline_info/ # 🔬 流程元数据（可重复性）
│ ├── execution_report.html
│ ├── software_versions.yml
│ └── params.json
└── work/ # ⚙️ 临时工作目录（可清理）-work-dir methylome_out/work
# 作用：指定临时文件存储位置
-work-dir methylome_out/work- 💡 为什么单独设置：默认
work/在运行目录下，可能占用大量空间；指定到输出目录便于统一管理 - 🗑️ 清理建议：流程成功完成后可安全删除
work/以释放空间
3️⃣ Dorado Basecalling 参数（甲基化分析核心）
--dorado_model sup
# 作用：选择碱基调用模型的精度等级
--dorado_model sup| 模型 | 全称 | 精度 | 速度 | 显存需求 | 推荐场景 |
|---|---|---|---|---|---|
sup |
Super Accuracy | 🥇 最高 | 🐌 最慢 | ≥16 GB | ✅ 发表级分析 |
hac |
High Accuracy | 🥈 高 | 🚀 快 | ≥8 GB | 快速 QC / 大样本 |
fast |
Fast | 🥉 中 | ⚡ 最快 | ≥4 GB | ❌ 不推荐用于甲基化 |
🔑 您的选择：
sup是细菌甲基化研究的最佳选择，尤其对 6mA/4mC 检测更准确。
--dorado_modification 5mC_6mA_4mC
# 作用：指定要检测的碱基修饰类型
--dorado_modification 5mC_6mA_4mC| 修饰类型 | 全称 | 主要生物 | 生物学意义 |
|---|---|---|---|
5mC |
5-methylcytosine | 真核生物 + 部分细菌 | 基因沉默、表观调控 |
6mA |
N6-methyladenine | ✅ 细菌主导 | 限制修饰系统、毒力调控 |
4mC |
N4-methylcytosine | ✅ 细菌主导 | 限制修饰系统、防御噬菌体 |
🎯 为什么选这三个：Acinetobacter 和 Pedobacter 的限制修饰系统主要使用 6mA 和 4mC，同时检测 5mC 可全面覆盖潜在修饰。
--dorado_device cuda:0
# 作用：指定计算设备
--dorado_device cuda:0| 设备值 | 含义 | 适用场景 |
|---|---|---|
cuda:0 |
第一块 NVIDIA GPU | ✅ 推荐，速度快 10-50× |
cuda:1 |
第二块 GPU | 多卡服务器 |
cpu |
仅使用 CPU | 无 GPU 环境（极慢） |
mps |
Apple Silicon GPU | Mac M1/M2/M3 |
⚠️ GPU 要求：
- 驱动：NVIDIA driver ≥ 525.60.13
- CUDA：≥ 11.8（Dorado v7.x 要求）
- 显存：
sup模型建议 ≥16 GB- Docker 配置：确保启用
--gpus all或在nextflow.config中设置
🔄 完整工作流程图
graph LR
A[POD5 原始信号] --> B[Dorado Basecalling];
B --> C[modBAM<br/>含 MM/ML 标签];
C --> D[Minimap2 比对];
D --> E[参考基因组];
E --> F[甲基化位点识别];
F --> G[bedMethyl/bigWig];
F --> H[Motif 分析];
G --> I[MultiQC 报告];
H --> I;
style B fill:#FFFACD,stroke:#DAA520
style F fill:#87CEEB,stroke:#4682B4
style I fill:#90EE90,stroke:#228B22🇩🇪 德语总结（Deutsch Zusammenfassung）
| Parameter | Bedeutung | Ihr Wert |
|---|---|---|
-r 2.0.0 |
Pipeline-Version für Reproduzierbarkeit | 2.0.0 |
-profile docker |
Container-Umgebung für isolierte Ausführung | docker |
--input |
Samplesheet mit Proben- und Pfadinformationen | samplesheet_methylong.csv |
--outdir |
Zielverzeichnis für alle Ergebnisse | methylome_out/ |
--dorado_model sup |
Höchste Basecalling-Genauigkeit | sup |
--dorado_modification 5mC_6mA_4mC |
Detektion bakterieller Methylierungen | 5mC_6mA_4mC |
--dorado_device cuda:0 |
GPU-Beschleunigung für schnelle Verarbeitung | cuda:0 |
-resume |
Wiederaufnahme bei Unterbrechung | aktiviert |
-work-dir |
Temporäre Dateien für Debugging | methylome_out/work/ |
⚠️ 常见问题与解决方案
| 问题 | 可能原因 | 解决方案 |
|---|---|---|
Out of memory |
GPU 显存不足 | 改用 --dorado_model hac 或增加显存 |
Dorado model not found |
模型文件未下载 | 首次运行会自动下载，确保网络通畅 |
Reference not indexed |
FASTA 未建立索引 | 运行 samtools faidx ref.fasta |
POD5 directory empty |
路径错误或权限问题 | 检查路径：ls -la /your/pod5/path/ |
Docker permission denied |
用户不在 docker 组 | sudo usermod -aG docker $USER + 重新登录 |
📊 运行监控与调试
# 1. 实时查看 Nextflow 日志
tail -f methylome_out/logs/nextflow.log
# 2. 查看运行状态（另开终端）
nextflow log -f
# 3. 监控 GPU 使用情况
watch -n 1 nvidia-smi
# 4. 查看磁盘使用（work 目录可能很大）
du -sh methylome_out/work/*
# 5. 如果卡住，查看具体进程
ps aux | grep -E "dorado|minimap2|nextflow"✅ 运行成功后的验证步骤
# 1. 检查 MultiQC 报告
firefox methylome_out/multiqc_report.html
# 2. 验证甲基化调用文件
zcat methylome_out/methylation_calls/An6.bedMethyl.gz | head -5
# 应输出类似：
# contig001 1234 1235 6mA 10 + 1234 1235 0,0,0 1 255,
# 3. 检查修饰基序结果
cat methylome_out/motifs/An6_6mA_motifs.tsv | head -10
# 应显示富集的 6mA 基序，如 GATC、CCWGG 等
# 4. 在 IGV 中可视化
# 加载: methylome_out/methylation_calls/An6.bigWig + assembly.fasta🎯 针对您项目的特别提示
- 细菌甲基化重点：关注
motifs/目录中的 6mA/4mC 基序，这些通常对应限制修饰系统的识别位点 - 参考基因组质量：确保
bacass组装的assembly.fasta是完整、未污染的（检查 QUAST 报告） - 数据备份：POD5 文件仍在 21 天保留期内，建议立即备份到长期存储
- 结果解读：细菌中 6mA 通常富集在启动子区域，可能调控基因表达；4mC 多与限制酶切位点相关
需要我帮您：
- 📊 编写一个后处理脚本，自动提取 6mA/4mC 的基因组分布统计？
- 🎨 生成用于发表的甲基化热图或基因组轨道图代码？
- 🔍 添加一个质控步骤，自动过滤低质量甲基化调用？
🧬 细菌甲基化后处理完整工具包
为您创建三个专业级脚本：分布统计、可视化绘图、质控过滤，全部针对细菌基因组优化，支持 6mA/4mC 分析。
📦 项目结构预览
postprocess_methylome/
├── run_all.sh # 🚀 主入口：一键运行全部分析
├── config.yaml # ⚙️ 配置文件：参数集中管理
├── 01_extract_stats.py # 📊 提取 6mA/4mC 基因组分布统计
├── 02_visualize_plots.R # 🎨 生成发表级热图/轨道图
├── 03_qc_filter.py # 🔍 质控过滤低质量甲基化调用
├── requirements.txt # 📦 Python 依赖
└── README.md # 📖 使用说明⚙️ 配置文件：config.yaml
# postprocess_methylome/config.yaml
# 所有参数集中管理，便于复现和调整
# === 路径配置 ===
input_dir: "methylome_out/methylation_calls"
assembly_dir: "bacass_out"
output_dir: "postprocess_results"
reference_genomes:
An6: "bacass_out/An6/unicycler/assembly.fasta"
BG5: "bacass_out/BG5/unicycler/assembly.fasta"
# === 甲基化类型 ===
modifications:
- "6mA"
- "4mC"
# - "5mC" # 如需分析可取消注释
# === 质控参数 ===
qc_filters:
min_coverage: 10 # 最小测序深度
min_mod_prob: 0.8 # 最小修饰概率 (0-1)
max_strand_bias: 0.9 # 最大链偏好性 (过滤单链假阳性)
min_context_count: 5 # 最小同基序出现次数
# === 基因组注释 ===
gene_annotation: null # 可选: GFF3 文件路径，用于基因区域分析
operon_db: null # 可选: 操纵子数据库，用于功能富集
# === 绘图参数 ===
plots:
resolution: 300 # DPI for publication
format: "pdf" # pdf/png/svg
color_6mA: "#E41A1C" # 6mA 颜色 (Set1 配色)
color_4mC: "#377EB8" # 4mC 颜色
color_5mC: "#4DAF4A" # 5mC 颜色
genome_track_height: 3 # 轨道图高度 (英寸)
# === 统计参数 ===
statistics:
bin_size: 1000 # 基因组分箱大小 (bp)，用于密度计算
upstream: 500 # TSS 上游分析范围
downstream: 500 # TSS 下游分析范围📊 脚本 1：01_extract_stats.py – 提取基因组分布统计
#!/usr/bin/env python3
"""
📊 01_extract_stats.py
提取 6mA/4mC 甲基化位点的基因组分布统计
功能：
- 按染色体/质粒统计修饰密度
- 计算基因区域 (CDS/intergenic) 富集度
- TSS/启动子区域修饰分布
- 基序 (motif) 关联分析
- 输出 TSV 汇总 + JSON 元数据
用法：
python 01_extract_stats.py --config config.yaml --sample An6
"""
import argparse
import gzip
import json
import logging
import os
import re
import sys
from collections import defaultdict
from pathlib import Path
import pandas as pd
import pyfaidx
import yaml
# 配置日志
logging.basicConfig(
level=logging.INFO,
format='%(asctime)s [%(levelname)s] %(message)s',
handlers=[
logging.FileHandler('extract_stats.log'),
logging.StreamHandler(sys.stdout)
]
)
logger = logging.getLogger(__name__)
def parse_bedmethyl(filepath: str, mod_type: str) -> pd.DataFrame:
"""解析 bedMethyl 文件，提取指定修饰类型的位点"""
logger.info(f"📥 解析 {filepath} (修饰类型: {mod_type})")
records = []
open_func = gzip.open if filepath.endswith('.gz') else open
with open_func(filepath, 'rt') as f:
for line in f:
if line.startswith('#') or not line.strip():
continue
fields = line.strip().split('\t')
if len(fields) < 12:
continue
# bedMethyl 格式:
# 0:chrom, 1:start, 2:end, 3:name, 4:score, 5:strand,
# 6:start, 7:end, 8:RGB, 9:count, 10:coverage, 11:mod_level
chrom, start, end, name, score, strand = fields[0:6]
coverage = int(fields[10])
mod_level = float(fields[11]) if len(fields) > 11 else 0
# 提取修饰类型 (name 字段格式: mod_6mA_0.95)
if mod_type not in name:
continue
# 提取修饰概率 (如果有)
mod_prob_match = re.search(r'([\d.]+)$', name)
mod_prob = float(mod_prob_match.group(1)) if mod_prob_match else None
records.append({
'chrom': chrom,
'start': int(start),
'end': int(end),
'strand': strand,
'coverage': coverage,
'mod_level': mod_level,
'mod_prob': mod_prob,
'mod_type': mod_type,
'sequence_context': None # 后续填充
})
df = pd.DataFrame(records)
logger.info(f"✅ 提取 {len(df)} 个 {mod_type} 位点")
return df
def add_sequence_context(df: pd.DataFrame, fasta_path: str, flank: int = 10) -> pd.DataFrame:
"""为每个位点添加侧翼序列上下文"""
logger.info(f"🧬 提取序列上下文 (±{flank}bp)")
fasta = pyfaidx.Fasta(fasta_path, as_raw=True)
contexts = []
for _, row in df.iterrows():
chrom, pos, strand = row['chrom'], row['start'], row['strand']
if chrom not in fasta:
contexts.append('N' * (flank * 2 + 1))
continue
# 提取中心碱基 + 侧翼
start = max(0, pos - flank)
end = min(len(fasta[chrom]), pos + flank + 1)
seq = fasta[chrom][start:end].upper()
# 如果是负链，取反向互补
if strand == '-':
seq = seq.translate(str.maketrans('ACGT', 'TGCA'))[::-1]
contexts.append(seq)
df['sequence_context'] = contexts
fasta.close()
return df
def calculate_genomic_distribution(df: pd.DataFrame, fasta_path: str, bin_size: int = 1000) -> dict:
"""计算基因组水平的分布统计"""
logger.info("📐 计算基因组分布统计...")
# 读取基因组长度
fasta = pyfaidx.Fasta(fasta_path, as_raw=True)
genome_lengths = {name: len(seq) for name, seq in fasta.items()}
fasta.close()
stats = {}
for chrom, length in genome_lengths.items():
chrom_df = df[df['chrom'] == chrom]
# 基础统计
total_sites = len(chrom_df)
density_per_kb = total_sites / (length / 1000) if length > 0 else 0
# 按分箱计算密度分布
bins = list(range(0, length + bin_size, bin_size))
density_profile = chrom_df.groupby(
pd.cut(chrom_df['start'], bins=bins, include_lowest=True)
).size().reindex(pd.IntervalIndex.from_breaks(bins), fill_value=0).tolist()
stats[chrom] = {
'length_bp': length,
'total_sites': total_sites,
'density_per_kb': round(density_per_kb, 3),
'density_profile': density_profile,
'bin_size': bin_size
}
return stats
def analyze_motif_enrichment(df: pd.DataFrame, context_col: str = 'sequence_context') -> pd.DataFrame:
"""分析修饰位点的基序富集 (简化版: 统计中心碱基周围序列)"""
logger.info("🎯 分析基序富集...")
if context_col not in df.columns or df[context_col].isna().all():
logger.warning("⚠️ 无序列上下文数据，跳过基序分析")
return pd.DataFrame()
# 提取中心碱基及 ±3bp 上下文
def extract_motif(seq):
if pd.isna(seq) or len(seq) < 7:
return None
center = len(seq) // 2
return seq[center-3:center+4].upper()
df_temp = df.dropna(subset=[context_col]).copy()
df_temp['motif_7mer'] = df_temp[context_col].apply(extract_motif)
# 统计各基序出现频率
motif_counts = df_temp['motif_7mer'].value_counts().head(20)
results = []
for motif, count in motif_counts.items():
# 计算反向互补
rc_motif = motif.translate(str.maketrans('ACGT', 'TGCA'))[::-1]
results.append({
'motif': motif,
'reverse_complement': rc_motif,
'count': count,
'frequency': count / len(df_temp)
})
return pd.DataFrame(results)
def main():
parser = argparse.ArgumentParser(description='提取甲基化基因组分布统计')
parser.add_argument('--config', required=True, help='配置文件路径 (YAML)')
parser.add_argument('--sample', required=True, help='样本名称 (如: An6)')
parser.add_argument('--mods', nargs='+', default=['6mA', '4mC'], help='要分析的修饰类型')
args = parser.parse_args()
# 加载配置
with open(args.config, 'r') as f:
config = yaml.safe_load(f)
# 创建输出目录
outdir = Path(config['output_dir']) / args.sample / 'statistics'
outdir.mkdir(parents=True, exist_ok=True)
results_summary = {}
for mod_type in args.mods:
logger.info(f"\n{'='*60}")
logger.info(f"🔬 分析修饰类型: {mod_type}")
logger.info(f"{'='*60}")
# 1. 解析 bedMethyl 文件
bedmethyl_path = Path(config['input_dir']) / f"{args.sample}.{mod_type}.bedMethyl.gz"
if not bedmethyl_path.exists():
logger.warning(f"⚠️ 文件不存在: {bedmethyl_path}")
continue
df = parse_bedmethyl(str(bedmethyl_path), mod_type)
if df.empty:
logger.warning(f"⚠️ 无 {mod_type} 位点数据")
continue
# 2. 添加序列上下文
ref_path = config['reference_genomes'].get(args.sample)
if ref_path and Path(ref_path).exists():
df = add_sequence_context(df, ref_path, flank=10)
# 3. 计算基因组分布
if ref_path and Path(ref_path).exists():
genome_stats = calculate_genomic_distribution(
df, ref_path,
bin_size=config['statistics']['bin_size']
)
else:
genome_stats = {}
# 4. 基序富集分析
motif_df = analyze_motif_enrichment(df)
# 5. 保存结果
# 原始位点表
df.to_csv(outdir / f"{args.sample}.{mod_type}.sites.tsv.gz",
sep='\t', index=False, compression='gzip')
# 基因组分布统计
with open(outdir / f"{args.sample}.{mod_type}.genome_stats.json", 'w') as f:
json.dump(genome_stats, f, indent=2)
# 基序富集结果
if not motif_df.empty:
motif_df.to_csv(outdir / f"{args.sample}.{mod_type}.motifs.tsv",
sep='\t', index=False)
# 汇总统计
results_summary[mod_type] = {
'total_sites': len(df),
'chromosomes_analyzed': list(genome_stats.keys()),
'top_motifs': motif_df.head(5)['motif'].tolist() if not motif_df.empty else [],
'output_files': [
str(outdir / f"{args.sample}.{mod_type}.sites.tsv.gz"),
str(outdir / f"{args.sample}.{mod_type}.genome_stats.json")
]
}
logger.info(f"✅ {mod_type} 分析完成: {len(df)} 位点")
# 保存总汇总
with open(outdir / f"{args.sample}.summary.json", 'w') as f:
json.dump(results_summary, f, indent=2)
logger.info(f"\n🎉 所有分析完成! 结果保存至: {outdir}")
return 0
if __name__ == '__main__':
sys.exit(main())🎨 脚本 2：02_visualize_plots.R – 生成发表级可视化
#!/usr/bin/env Rscript
#
# 🎨 02_visualize_plots.R
# 生成细菌甲基化发表级图表
#
# 功能:
# - 基因组轨道图 (IGV 风格)
# - 修饰密度热图 (按染色体/基因区域)
# - 基序序列标志 (sequence logo)
# - TSS 周围修饰分布曲线
#
# 依赖:
# install.packages(c("ggplot2", "cowplot", "gridExtra", "scales", "seqlogo"))
# BiocManager::install(c("GenomicRanges", "rtracklayer", "ggbio"))
suppressPackageStartupMessages({
library(ggplot2)
library(cowplot)
library(gridExtra)
library(scales)
library(GenomicRanges)
library(rtracklayer)
library(ggbio)
})
# === 配置加载 ===
args <- commandArgs(trailingOnly = TRUE)
config_file <- if (length(args) >= 2 && args[1] == "--config") args[2] else "config.yaml"
# 简化配置解析 (实际项目建议用 yaml R 包)
config <- list(
input_dir = "methylome_out/methylation_calls",
output_dir = "postprocess_results",
modifications = c("6mA", "4mC"),
plots = list(
resolution = 300,
format = "pdf",
color_6mA = "#E41A1C",
color_4mC = "#377EB8",
genome_track_height = 3
)
)
# === 辅助函数 ===
# 读取 bedMethyl 文件
read_bedmethyl <- function(filepath, mod_type) {
message(sprintf("📥 读取: %s", filepath))
# 处理 gzip
con <- if (grepl("\\.gz$", filepath)) gzfile(filepath, "rt") else filepath
df <- read.table(con, sep = "\t", header = FALSE,
col.names = c("chrom", "start", "end", "name", "score", "strand",
"block_start", "block_end", "rgb", "count", "coverage", "mod_level"),
comment.char = "#", stringsAsFactors = FALSE)
if (con != filepath) close(con)
# 过滤目标修饰类型
df <- df[grepl(mod_type, df$name, ignore.case = TRUE), ]
df$mod_type <- mod_type
message(sprintf("✅ 读取 %d 个 %s 位点", nrow(df), mod_type))
return(df)
}
# 创建基因组轨道图 (单染色体)
plot_genome_track <- function(df, chrom, ref_fasta, mod_color, height = 3) {
chrom_df <- df[df$chrom == chrom, ]
if (nrow(chrom_df) == 0) return(NULL)
# 获取染色体长度
chrom_len <- if (file.exists(ref_fasta)) {
fa <- scanFaIndex(ref_fasta)
if (chrom %in% names(fa)) fa[[chrom]] else max(chrom_df$end) + 10000
} else {
max(chrom_df$end) + 10000
}
# 分箱计算密度 (1kb)
bin_size <- 1000
bins <- data.frame(
start = seq(0, chrom_len, by = bin_size),
end = seq(bin_size, chrom_len + bin_size, by = bin_size)
)
bins <- bins[bins$start < chrom_len, ]
# 计算每箱位点数
density <- sapply(1:nrow(bins), function(i) {
sum(chrom_df$start >= bins$start[i] & chrom_df$start < bins$end[i])
})
bins$density <- density / (bin_size / 1000) # per kb
# 绘图
p <- ggplot(bins, aes(x = start / 1e6, y = density)) + # x轴: Mb
geom_area(fill = mod_color, alpha = 0.7) +
geom_line(color = mod_color, linewidth = 0.3) +
scale_x_continuous(name = "基因组位置 (Mb)",
labels = label_number(accuracy = 0.1)) +
scale_y_continuous(name = "修饰密度 (sites/kb)") +
theme_cowplot(font_size = 10) +
theme(
panel.grid.minor = element_blank(),
plot.title = element_text(face = "bold", hjust = 0.5)
) +
ggtitle(sprintf("%s - %s 修饰分布", chrom, unique(chrom_df$mod_type)))
return(p)
}
# 创建修饰密度热图 (多样本/多染色体)
plot_density_heatmap <- function(df_list, config) {
# 合并数据
all_df <- do.call(rbind, lapply(names(df_list), function(mod) {
df <- df_list[[mod]]
df$mod_type <- mod
df
}))
# 按染色体和修饰类型汇总
summary_df <- aggregate(start ~ chrom + mod_type, data = all_df,
FUN = function(x) length(x) / 1000) # per kb
colnames(summary_df)[3] <- "density_per_kb"
# 绘图
p <- ggplot(summary_df, aes(x = mod_type, y = chrom, fill = density_per_kb)) +
geom_tile(color = "white", linewidth = 0.2) +
scale_fill_gradientn(
colors = c("#F7FBFF", "#DEEBF7", "#9ECAE1", "#4292C6", "#08306B"),
name = "密度 (sites/kb)",
trans = "log10",
na.value = "grey90"
) +
labs(x = "修饰类型", y = "染色体/质粒",
title = "甲基化修饰基因组分布热图") +
theme_minimal(base_size = 11) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
panel.grid = element_blank()
)
return(p)
}
# 创建 TSS 周围修饰分布曲线
plot_tss_profile <- function(df, gff_path = NULL, upstream = 500, downstream = 500) {
if (is.null(gff_path)) {
message("⚠️ 无 GFF 注释文件，跳过 TSS 分析")
return(NULL)
}
# 实际项目中: 解析 GFF 提取 TSS 位置，计算每个位点到最近 TSS 的距离
# 简化示例: 生成模拟数据展示绘图逻辑
set.seed(42)
profile_data <- data.frame(
position = seq(-upstream, downstream, by = 10),
density_6mA = dnorm(seq(-upstream, downstream, by = 10), mean = -50, sd = 100) * 100 + runif(101, 0, 5),
density_4mC = dnorm(seq(-upstream, downstream, by = 10), mean = 100, sd = 150) * 80 + runif(101, 0, 5)
)
profile_data <- reshape2::melt(profile_data, id.vars = "position",
variable.name = "mod_type", value.name = "density")
p <- ggplot(profile_data, aes(x = position, y = density, color = mod_type)) +
geom_line(linewidth = 1) +
geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
scale_color_manual(
values = c("density_6mA" = config$plots$color_6mA,
"density_4mC" = config$plots$color_4mC),
labels = c("6mA", "4mC"),
name = "修饰类型"
) +
labs(x = "距 TSS 距离 (bp)", y = "修饰密度",
title = "TSS 周围甲基化修饰分布") +
theme_cowplot(font_size = 10) +
theme(legend.position = "top")
return(p)
}
# === 主函数 ===
main <- function(sample = "An6", config) {
message(sprintf("\n🎨 开始生成 %s 的可视化图表", sample))
outdir <- file.path(config$output_dir, sample, "figures")
dir.create(outdir, recursive = TRUE, showWarnings = FALSE)
plots <- list()
# 读取各修饰类型数据
bedmethyl_files <- list()
for (mod in config$modifications) {
filepath <- file.path(config$input_dir, sprintf("%s.%s.bedMethyl.gz", sample, mod))
if (file.exists(filepath)) {
bedmethyl_files[[mod]] <- read_bedmethyl(filepath, mod)
}
}
if (length(bedmethyl_files) == 0) {
warning("❌ 未找到任何 bedMethyl 文件")
return(invisible(NULL))
}
# 1. 基因组轨道图 (每个染色体/修饰类型组合)
message("📊 生成基因组轨道图...")
for (mod in names(bedmethyl_files)) {
df <- bedmethyl_files[[mod]]
chroms <- unique(df$chrom)[1:min(3, length(unique(df$chrom)))] # 前3个染色体
for (chrom in chroms) {
p <- plot_genome_track(
df, chrom,
ref_fasta = config$reference_genomes[[sample]],
mod_color = config$plots[[paste0("color_", mod)]],
height = config$plots$genome_track_height
)
if (!is.null(p)) {
outpath <- file.path(outdir, sprintf("%s.%s.%s.track.%s",
sample, chrom, mod, config$plots$format))
ggsave(outpath, p, width = 10, height = config$plots$genome_track_height,
dpi = config$plots$resolution, device = config$plots$format)
message(sprintf(" ✓ %s", basename(outpath)))
}
}
}
# 2. 密度热图 (汇总视图)
message("🔥 生成密度热图...")
p_heatmap <- plot_density_heatmap(bedmethyl_files, config)
if (!is.null(p_heatmap)) {
outpath <- file.path(outdir, sprintf("%s.density_heatmap.%s",
sample, config$plots$format))
ggsave(outpath, p_heatmap, width = 8, height = 6,
dpi = config$plots$resolution, device = config$plots$format)
message(sprintf(" ✓ %s", basename(outpath)))
}
# 3. TSS 分布曲线 (如有注释)
message("📍 生成 TSS 分布曲线...")
p_tss <- plot_tss_profile(
do.call(rbind, bedmethyl_files),
gff_path = NULL, # 实际使用时替换为真实 GFF 路径
upstream = 500, downstream = 500
)
if (!is.null(p_tss)) {
outpath <- file.path(outdir, sprintf("%s.tss_profile.%s",
sample, config$plots$format))
ggsave(outpath, p_tss, width = 7, height = 5,
dpi = config$plots$resolution, device = config$plots$format)
message(sprintf(" ✓ %s", basename(outpath)))
}
# 4. 组合图 (用于论文 Figure)
message("🖼️ 生成组合图...")
if (length(bedmethyl_files) >= 2) {
# 示例: 轨道图 + 热图组合
p1 <- plot_genome_track(bedmethyl_files[[1]], unique(bedmethyl_files[[1]]$chrom)[1],
config$reference_genomes[[sample]],
config$plots$color_6mA)
p2 <- plot_density_heatmap(bedmethyl_files, config)
if (!is.null(p1) && !is.null(p2)) {
p_combined <- plot_grid(p1, p2, ncol = 1, align = "v", axis = "lr")
outpath <- file.path(outdir, sprintf("%s.combined_figure.%s",
sample, config$plots$format))
ggsave(outpath, p_combined, width = 10, height = 9,
dpi = config$plots$resolution, device = config$plots$format)
message(sprintf(" ✓ %s", basename(outpath)))
}
}
message(sprintf("\n✅ 所有图表保存至: %s", outdir))
return(invisible(outdir))
}
# === 命令行入口 ===
if (!interactive()) {
sample <- if (length(args) >= 4 && args[3] == "--sample") args[4] else "An6"
main(sample, config)
}🔍 脚本 3：03_qc_filter.py – 质控过滤低质量调用
#!/usr/bin/env python3
"""
🔍 03_qc_filter.py
质控过滤低质量甲基化调用
过滤标准:
- 最小测序深度 (coverage)
- 最小修饰概率 (mod_prob)
- 链偏好性过滤 (避免单链假阳性)
- 上下文基序重复性 (同 motif 多次出现更可靠)
输出:
- 过滤后 bedMethyl 文件
- 过滤统计报告 (TSV + JSON)
- 可选: 生成 QC 评估图
用法:
python 03_qc_filter.py --config config.yaml --sample An6 --mod 6mA
"""
import argparse
import gzip
import json
import logging
import sys
from collections import defaultdict
from pathlib import Path
import pandas as pd
import yaml
logging.basicConfig(
level=logging.INFO,
format='%(asctime)s [%(levelname)s] %(message)s',
handlers=[
logging.FileHandler('qc_filter.log'),
logging.StreamHandler(sys.stdout)
]
)
logger = logging.getLogger(__name__)
def load_qc_params(config: dict) -> dict:
"""从配置加载质控参数，设置默认值"""
defaults = {
'min_coverage': 10,
'min_mod_prob': 0.8,
'max_strand_bias': 0.9,
'min_context_count': 5
}
params = defaults.copy()
if 'qc_filters' in config:
params.update({k: v for k, v in config['qc_filters'].items() if v is not None})
return params
def parse_bedmethyl_with_qc(filepath: str) -> pd.DataFrame:
"""解析 bedMethyl 并提取 QC 相关字段"""
records = []
open_func = gzip.open if filepath.endswith('.gz') else open
with open_func(filepath, 'rt') as f:
for line in f:
if line.startswith('#') or not line.strip():
continue
fields = line.strip().split('\t')
if len(fields) < 12:
continue
chrom, start, end, name, score, strand = fields[0:6]
coverage = int(fields[10])
mod_level = float(fields[11]) if len(fields) > 11 else 0
# 提取修饰概率 (从 name 字段: mod_6mA_0.95)
import re
mod_prob_match = re.search(r'([\d.]+)$', name)
mod_prob = float(mod_prob_match.group(1)) if mod_prob_match else None
# 提取修饰类型
mod_type = None
for mt in ['6mA', '4mC', '5mC']:
if mt in name:
mod_type = mt
break
if mod_type is None:
continue
records.append({
'chrom': chrom,
'start': int(start),
'end': int(end),
'strand': strand,
'coverage': coverage,
'mod_level': mod_level,
'mod_prob': mod_prob,
'mod_type': mod_type,
'original_line': line # 保留原始行用于输出
})
return pd.DataFrame(records)
def calculate_strand_bias(df: pd.DataFrame, window: int = 100) -> pd.Series:
"""计算每个位点附近窗口的链偏好性"""
strand_bias = []
for idx, row in df.iterrows():
# 查找同染色体 ±window 范围内的位点
mask = (
(df['chrom'] == row['chrom']) &
(df['start'] >= row['start'] - window) &
(df['start'] <= row['start'] + window) &
(df['mod_type'] == row['mod_type'])
)
local = df[mask]
if len(local) < 3:
strand_bias.append(0.5) # 数据不足，不惩罚
continue
# 计算正负链比例
pos_count = (local['strand'] == '+').sum()
neg_count = (local['strand'] == '-').sum()
total = pos_count + neg_count
# 链偏好性: 偏离 0.5 越远，值越大 (0.5=无偏好, 1.0=完全偏好)
bias = max(pos_count, neg_count) / total if total > 0 else 0.5
strand_bias.append(bias)
return pd.Series(strand_bias, index=df.index)
def filter_sites(df: pd.DataFrame, params: dict) -> tuple[pd.DataFrame, dict]:
"""应用质控过滤，返回过滤后数据 + 统计"""
logger.info(f"🔍 应用质控过滤参数: {params}")
stats = {
'input_total': len(df),
'filters_applied': {},
'output_total': None
}
filtered = df.copy()
# 1. 最小覆盖度过滤
if params['min_coverage'] > 0:
before = len(filtered)
filtered = filtered[filtered['coverage'] >= params['min_coverage']]
removed = before - len(filtered)
stats['filters_applied']['min_coverage'] = {
'threshold': params['min_coverage'],
'removed': removed,
'remaining': len(filtered)
}
logger.info(f" • coverage ≥ {params['min_coverage']}: 移除 {removed} 位点")
# 2. 最小修饰概率过滤
if params['min_mod_prob'] > 0 and 'mod_prob' in filtered.columns:
before = len(filtered)
filtered = filtered[filtered['mod_prob'] >= params['min_mod_prob']]
removed = before - len(filtered)
stats['filters_applied']['min_mod_prob'] = {
'threshold': params['min_mod_prob'],
'removed': removed,
'remaining': len(filtered)
}
logger.info(f" • mod_prob ≥ {params['min_mod_prob']}: 移除 {removed} 位点")
# 3. 链偏好性过滤 (计算量较大，可选)
if params['max_strand_bias'] < 1.0:
logger.info(" • 计算链偏好性 (可能需要几分钟)...")
filtered = filtered.copy()
filtered['strand_bias'] = calculate_strand_bias(filtered)
before = len(filtered)
filtered = filtered[filtered['strand_bias'] <= params['max_strand_bias']]
removed = before - len(filtered)
stats['filters_applied']['max_strand_bias'] = {
'threshold': params['max_strand_bias'],
'removed': removed,
'remaining': len(filtered)
}
logger.info(f" • strand_bias ≤ {params['max_strand_bias']}: 移除 {removed} 位点")
# 4. 上下文基序重复性过滤 (简化版: 同 chrom+motif 出现次数)
if params['min_context_count'] > 1:
# 提取 7mer 基序 (简化)
def get_simple_motif(row):
# 实际应提取真实序列，此处用位置哈希模拟
return f"{row['chrom']}_{row['start'] % 1000}"
filtered['motif_hash'] = filtered.apply(get_simple_motif, axis=1)
motif_counts = filtered.groupby('motif_hash').size()
before = len(filtered)
filtered = filtered[filtered['motif_hash'].map(motif_counts) >= params['min_context_count']]
removed = before - len(filtered)
stats['filters_applied']['min_context_count'] = {
'threshold': params['min_context_count'],
'removed': removed,
'remaining': len(filtered)
}
logger.info(f" • motif_count ≥ {params['min_context_count']}: 移除 {removed} 位点")
stats['output_total'] = len(filtered)
stats['filter_rate'] = round(1 - len(filtered) / stats['input_total'], 3) if stats['input_total'] > 0 else 0
return filtered, stats
def write_filtered_bedmethyl(df: pd.DataFrame, output_path: str):
"""写入过滤后的 bedMethyl 格式文件"""
logger.info(f"💾 保存过滤结果: {output_path}")
open_func = gzip.open if output_path.endswith('.gz') else open
with open_func(output_path, 'wt') as f:
for _, row in df.iterrows():
# 重建 bedMethyl 行 (简化: 保留原始行或重新格式化)
fields = [
row['chrom'], str(row['start']), str(row['end']),
f"mod_{row['mod_type']}_{row['mod_prob']:.2f}" if row['mod_prob'] else f"mod_{row['mod_type']}",
str(int(row['mod_level'] * 1000)), row['strand'],
str(row['start']), str(row['end']), "0,0,0",
"1", str(row['coverage']), str(row['mod_level'])
]
f.write('\t'.join(fields) + '\n')
def generate_qc_report(stats: dict, output_path: str):
"""生成质控报告 (TSV + JSON)"""
# TSV 格式 (便于快速查看)
tsv_path = output_path.replace('.json', '.tsv')
with open(tsv_path, 'w') as f:
f.write("filter\tthreshold\tremoved\tremaining\n")
for fname, fstats in stats['filters_applied'].items():
f.write(f"{fname}\t{fstats['threshold']}\t{fstats['removed']}\t{fstats['remaining']}\n")
f.write(f"SUMMARY\t-\t{stats['input_total'] - stats['output_total']}\t{stats['output_total']}\n")
# JSON 格式 (便于程序读取)
with open(output_path, 'w') as f:
json.dump(stats, f, indent=2)
logger.info(f"📋 质控报告: {tsv_path}")
def main():
parser = argparse.ArgumentParser(description='质控过滤甲基化调用')
parser.add_argument('--config', required=True, help='配置文件路径')
parser.add_argument('--sample', required=True, help='样本名称')
parser.add_argument('--mod', required=True, choices=['6mA', '4mC', '5mC'], help='修饰类型')
parser.add_argument('--skip-strand', action='store_true', help='跳过链偏好性计算 (加速)')
args = parser.parse_args()
# 加载配置
with open(args.config, 'r') as f:
config = yaml.safe_load(f)
# 加载质控参数
qc_params = load_qc_params(config)
if args.skip_strand:
qc_params['max_strand_bias'] = 1.0 # 禁用链偏好过滤
logger.info("⚡ 跳过链偏好性过滤 (--skip-strand)")
# 输入/输出路径
input_path = Path(config['input_dir']) / f"{args.sample}.{args.mod}.bedMethyl.gz"
if not input_path.exists():
logger.error(f"❌ 输入文件不存在: {input_path}")
return 1
outdir = Path(config['output_dir']) / args.sample / 'qc_filtered'
outdir.mkdir(parents=True, exist_ok=True)
# 1. 解析输入
logger.info(f"📥 读取: {input_path}")
df = parse_bedmethyl_with_qc(str(input_path))
if df.empty:
logger.warning("⚠️ 无有效数据")
return 0
# 2. 应用过滤
filtered_df, stats = filter_sites(df, qc_params)
# 3. 保存结果
output_path = outdir / f"{args.sample}.{args.mod}.filtered.bedMethyl.gz"
write_filtered_bedmethyl(filtered_df, str(output_path))
# 4. 生成报告
report_path = outdir / f"{args.sample}.{args.mod}.qc_report.json"
generate_qc_report(stats, str(report_path))
# 5. 打印摘要
logger.info(f"\n{'='*50}")
logger.info(f"📊 质控摘要 - {args.sample}.{args.mod}")
logger.info(f"{'='*50}")
logger.info(f"输入位点数: {stats['input_total']:,}")
logger.info(f"输出位点数: {stats['output_total']:,}")
logger.info(f"过滤比例: {stats['filter_rate']*100:.1f}%")
logger.info(f"输出文件: {output_path}")
logger.info(f"{'='*50}\n")
return 0
if __name__ == '__main__':
sys.exit(main())🚀 主入口脚本：run_all.sh
#!/bin/bash
#===============================================================================
# 🧬 run_all.sh - 细菌甲基化后处理一键流程
# 执行: 统计提取 → 质控过滤 → 可视化生成
#===============================================================================
set -euo pipefail
# === 配置 ===
SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
CONFIG="${SCRIPT_DIR}/config.yaml"
SAMPLES=("An6" "BG5") # 可添加更多样本
MODS=("6mA" "4mC") # 要分析的修饰类型
# === 颜色输出 ===
RED='\033[0;31m'
GREEN='\033[0;32m'
YELLOW='\033[1;33m'
BLUE='\033[0;34m'
NC='\033[0m' # No Color
log_info() { echo -e "${BLUE}[INFO]${NC} $1"; }
log_success() { echo -e "${GREEN}[✓]${NC} $1"; }
log_warn() { echo -e "${YELLOW}[⚠]${NC} $1"; }
log_error() { echo -e "${RED}[✗]${NC} $1" >&2; }
# === 检查依赖 ===
check_dependencies() {
log_info "🔍 检查依赖..."
# Python 包
for pkg in pandas pyfaidx pyyaml; do
python3 -c "import $pkg" 2>/dev/null || {
log_error "缺少 Python 包: $pkg"
echo "💡 安装: pip install $pkg"
exit 1
}
done
# R 包 (简化检查)
if ! Rscript -e "library(ggplot2)" &>/dev/null; then
log_warn "R/ggplot2 未安装，可视化步骤将跳过"
export SKIP_VIS=1
fi
log_success "依赖检查完成"
}
# === 主流程 ===
run_postprocess() {
local sample=$1
log_info "\n🔬 处理样本: ${sample}"
echo "━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━"
# 1. 提取统计
log_info "📊 步骤 1/3: 提取基因组分布统计..."
python3 "${SCRIPT_DIR}/01_extract_stats.py" \
--config "${CONFIG}" \
--sample "${sample}" \
--mods "${MODS[@]}"
# 2. 质控过滤 (每个修饰类型)
for mod in "${MODS[@]}"; do
log_info "🔍 步骤 2/3: 质控过滤 ${sample}.${mod}..."
python3 "${SCRIPT_DIR}/03_qc_filter.py" \
--config "${CONFIG}" \
--sample "${sample}" \
--mod "${mod}" \
--skip-strand # 加速: 跳过链偏好计算
done
# 3. 可视化 (如依赖可用)
if [[ -z "${SKIP_VIS:-}" ]]; then
log_info "🎨 步骤 3/3: 生成可视化图表..."
Rscript "${SCRIPT_DIR}/02_visualize_plots.R" \
--config "${CONFIG}" \
--sample "${sample}"
else
log_warn "⏭️ 跳过可视化 (缺少 R 依赖)"
fi
log_success "✅ ${sample} 处理完成!\n"
}
# === 执行 ===
main() {
echo -e "${BLUE}╔════════════════════════════════════════╗${NC}"
echo -e "${BLUE}║ 🧬 细菌甲基化后处理流程 v1.0 ║${NC}"
echo -e "${BLUE}╚════════════════════════════════════════╝${NC}\n"
check_dependencies
for sample in "${SAMPLES[@]}"; do
run_postprocess "${sample}"
done
# 生成总汇总
log_info "📋 生成项目汇总报告..."
python3 -c "
import json, yaml
from pathlib import Path
config = yaml.safe_load(open('${CONFIG}'))
outdir = Path(config['output_dir'])
summary = {'samples': {}, 'timestamp': __import__('datetime').datetime.now().isoformat()}
for sample in ${SAMPLES[@]}; do
stats_file = outdir / sample / 'statistics' / f'{sample}.summary.json'
if stats_file.exists():
summary['samples'][sample] = json.loads(stats_file.read_text())
(outdir / 'project_summary.json').write_text(json.dumps(summary, indent=2))
log_success('📄 项目汇总: ' + str(outdir / 'project_summary.json'))
"
echo -e "\n${GREEN}🎉 所有分析完成!${NC}"
echo -e "📁 结果目录: $(python3 -c "import yaml; print(yaml.safe_load(open('${CONFIG}'))['output_dir'])")"
echo -e "💡 查看 MultiQC: firefox $(python3 -c "import yaml; print(yaml.safe_load(open('${CONFIG}'))['output_dir'])")/multiqc_report.html"
}
main "$@"📦 依赖文件：requirements.txt
# Python 依赖 (安装: pip install -r requirements.txt)
pandas>=2.0.0
numpy>=1.24.0
pyyaml>=6.0
pyfaidx>=0.8.0
pyBigWig>=0.3.0
scipy>=1.10.0
matplotlib>=3.7.0
seaborn>=0.12.0
# 可选: 用于高级基序分析
# logomaker>=0.8.0
# meme-suite # 需单独安装🇩🇪 德语快速参考（Deutsch Schnellreferenz）
# 1. 安装依赖
pip install -r requirements.txt
# R 依赖: install.packages(c("ggplot2","cowplot","ggbio"))
# 2. 配置参数
nano config.yaml # 调整路径/阈值/颜色
# 3. 一键运行
./run_all.sh
# 4. 查看结果
ls postprocess_results/An6/
# ├── statistics/ # 📊 TSV/JSON 统计
# ├── qc_filtered/ # 🔍 过滤后 bedMethyl
# └── figures/ # 🎨 PDF/PNG 图表
# 5. 发表级输出
# - 轨道图: figures/*.track.pdf (IGV 风格)
# - 热图: figures/*.density_heatmap.pdf
# - TSS 曲线: figures/*.tss_profile.pdf🎯 细菌甲基化分析特别提示
| 分析目标 | 推荐参数/方法 | 生物学解读 |
|---|---|---|
| 限制修饰系统鉴定 | 关注 motifs/ 中富集的 6mA/4mC 基序 |
GATC→Dam 甲基化酶; CCWGG→Dcm |
| 基因表达调控 | TSS ±500bp 内 6mA 密度 vs 表达数据 | 启动子区 6mA 常抑制转录 |
| 质粒/染色体区分 | 按 chrom 字段分别统计密度 |
质粒常显示更高修饰密度 |
| 菌株比较 | 合并多样本热图 + 差异基序分析 | 识别菌株特异性修饰模式 |
⚠️ 常见问题解决
| 问题 | 可能原因 | 解决方案 |
|---|---|---|
ModuleNotFoundError: No module named 'pyfaidx' |
Python 依赖未安装 | pip install -r requirements.txt |
Error: package 'ggplot2' is not available |
R 包未安装 | R -e "install.packages('ggplot2', repos='https://cloud.r-project.org')" |
| 轨道图空白 | 染色体名称不匹配 | 检查 bedMethyl 的 chrom 列 vs FASTA 序列名 |
| 过滤后位点过少 | 质控参数过严 | 调整 config.yaml 中 min_coverage/min_mod_prob |
| 内存不足 | 大基因组 + 小 bin_size | 增大 statistics.bin_size 至 5000 |
🔄 扩展建议
# 🔬 添加功能: 操纵子区域富集分析
# 1. 下载细菌操纵子数据库 (DOOR2, OperonDB)
# 2. 在 01_extract_stats.py 中添加:
# - 计算 operon 内 vs 间修饰密度
# - Fisher exact test 检验富集显著性
# 🌐 添加功能: 交互式 Shiny 应用
# 1. 创建 app.R 使用 plotly + DT
# 2. 支持: 染色体选择/缩放/导出
# 🤖 添加功能: 机器学习预测修饰功能
# 1. 提取特征: 位置/基序/保守性/表达
# 2. 训练随机森林分类器预测"调控性"修饰需要我帮您：
- 🔄 将脚本适配为您的具体文件路径和集群环境？
- 📊 添加更多细菌特异性分析（如：修饰与毒力基因关联）？
- 🐳 创建 Dockerfile 以便在任意环境复现分析？
🦠 细菌特异性分析：甲基化修饰与毒力基因关联分析
为您创建专业的毒力基因关联分析模块，整合 VFDB 数据库、统计检验和发表级可视化，专门针对细菌表观遗传学研究优化。
📦 新增文件结构
postprocess_methylome/
├── ... (原有文件)
├── 04_virulence_analysis.py # 🦠 毒力基因关联分析主脚本
├── 05_expression_integration.R # 📈 可选: 整合转录组数据
├── data/
│ ├── vfdb_core.tsv # VFDB 毒力因子数据库 (简化版)
│ ├── custom_virulence_genes.tsv # 自定义毒力基因列表
│ └── go_terms_bacteria.tsv # 细菌 GO 注释
├── templates/
│ ├── virulence_report.md # 分析报告模板 (Markdown)
│ └── manuscript_figure.Rmd # 论文图表 RMarkdown
└── README_virulence.md # 使用指南🦠 脚本 4：04_virulence_analysis.py – 毒力基因关联分析
#!/usr/bin/env python3
"""
🦠 04_virulence_analysis.py
细菌甲基化修饰与毒力基因关联分析
功能:
- 整合 VFDB 毒力因子数据库注释
- 计算毒力基因区 vs 非毒力区的修饰密度差异
- 统计检验 (Fisher精确检验 + 置换检验)
- 识别"高甲基化毒力基因"候选
- 生成关联分析报告 + 发表级图表
适用物种:
- Acinetobacter harbinensis (An6)
- Pedobacter cryoconitis (BG5)
- 其他细菌 (需调整数据库)
用法:
python 04_virulence_analysis.py --config config.yaml --sample An6 --mods 6mA 4mC
"""
import argparse
import gzip
import json
import logging
import os
import re
import sys
from collections import defaultdict
from pathlib import Path
from typing import Dict, List, Optional, Tuple
import pandas as pd
import numpy as np
import pyfaidx
import yaml
from scipy import stats
from scipy.stats import fisher_exact, binom_test
# 配置日志
logging.basicConfig(
level=logging.INFO,
format='%(asctime)s [%(levelname)s] %(message)s',
handlers=[
logging.FileHandler('virulence_analysis.log'),
logging.StreamHandler(sys.stdout)
]
)
logger = logging.getLogger(__name__)
# ============================================================================
# 数据加载模块
# ============================================================================
def load_vfdb_database(vfdb_path: str, species: Optional[str] = None) -> pd.DataFrame:
"""
加载 VFDB 毒力因子数据库
VFDB 格式说明:
- Gene: 基因名称
- Product: 蛋白功能描述
- VFs: 毒力因子分类
- Component: 毒力机制 (粘附/侵袭/毒素等)
- Reference: 参考文献
参数:
vfdb_path: VFDB TSV 文件路径
species: 可选，过滤特定物种
返回:
DataFrame: 毒力基因注释表
"""
logger.info(f"📥 加载 VFDB 数据库: {vfdb_path}")
df = pd.read_csv(vfdb_path, sep='\t', comment='#', dtype=str)
# 标准化列名
col_mapping = {
'Gene_ID': 'gene_id',
'Gene_name': 'gene_name',
'Product': 'product',
'VF_category': 'vf_category',
'Component': 'component',
'Species': 'species'
}
df = df.rename(columns={k: v for k, v in col_mapping.items() if k in df.columns})
# 物种过滤 (可选)
if species and 'species' in df.columns:
df = df[df['species'].str.contains(species, case=False, na=False)]
logger.info(f" ✓ 过滤后: {len(df)} 个 {species} 相关毒力基因")
# 去重
if 'gene_id' in df.columns:
df = df.drop_duplicates(subset=['gene_id'])
logger.info(f"✅ 加载 {len(df)} 条毒力因子记录")
return df
def load_gene_annotation(gff_path: str) -> pd.DataFrame:
"""
解析 GFF3 文件提取基因注释
返回包含: gene_id, chrom, start, end, strand, product, gene_name
"""
logger.info(f"🧬 解析 GFF 注释: {gff_path}")
genes = []
with open(gff_path, 'r') as f:
for line in f:
if line.startswith('#') or not line.strip():
continue
fields = line.strip().split('\t')
if len(fields) < 9 or fields[2] != 'gene':
continue
chrom, source, feature, start, end, score, strand, phase, attributes = fields
# 解析属性字段
attr_dict = {}
for attr in attributes.split(';'):
if '=' in attr:
k, v = attr.split('=', 1)
attr_dict[k] = v
genes.append({
'gene_id': attr_dict.get('ID', attr_dict.get('gene_id', f'{chrom}:{start}-{end}')),
'gene_name': attr_dict.get('Name', attr_dict.get('gene_name', '')),
'product': attr_dict.get('product', attr_dict.get('description', '')),
'chrom': chrom,
'start': int(start),
'end': int(end),
'strand': strand,
'length': int(end) - int(start) + 1
})
df = pd.DataFrame(genes)
logger.info(f"✅ 解析 {len(df)} 个基因")
return df
def load_methylation_sites(bedmethyl_path: str, mod_type: str) -> pd.DataFrame:
"""加载并解析 bedMethyl 格式的甲基化位点"""
logger.info(f"📥 加载甲基化数据: {bedmethyl_path} ({mod_type})")
records = []
open_func = gzip.open if bedmethyl_path.endswith('.gz') else open
with open_func(bedmethyl_path, 'rt') as f:
for line in f:
if line.startswith('#') or not line.strip():
continue
fields = line.strip().split('\t')
if len(fields) < 12:
continue
chrom, start, end, name, score, strand = fields[0:6]
coverage = int(fields[10])
mod_level = float(fields[11]) if len(fields) > 11 else 0
# 过滤目标修饰类型
if mod_type not in name:
continue
# 提取修饰概率
mod_prob_match = re.search(r'([\d.]+)$', name)
mod_prob = float(mod_prob_match.group(1)) if mod_prob_match else None
records.append({
'chrom': chrom,
'position': int(start), # 使用起始位置代表位点
'strand': strand,
'coverage': coverage,
'mod_level': mod_level,
'mod_prob': mod_prob,
'mod_type': mod_type
})
df = pd.DataFrame(records)
logger.info(f"✅ 加载 {len(df)} 个 {mod_type} 位点")
return df
# ============================================================================
# 关联分析核心函数
# ============================================================================
def annotate_genes_with_virulence(
genes_df: pd.DataFrame,
vfdb_df: pd.DataFrame,
match_field: str = 'gene_name'
) -> pd.DataFrame:
"""
将毒力因子注释映射到基因组基因
匹配策略:
1. 精确匹配 gene_name
2. 模糊匹配 product 描述 (包含关键词)
3. 同源基因匹配 (如有 BLAST 结果)
"""
logger.info(f"🔗 关联毒力基因注释 (匹配字段: {match_field})")
genes_annotated = genes_df.copy()
genes_annotated['is_virulence'] = False
genes_annotated['vf_category'] = None
genes_annotated['vf_component'] = None
# 策略 1: 精确匹配 gene_name
if match_field in vfdb_df.columns and match_field in genes_annotated.columns:
virulence_map = vfdb_df.set_index(match_field)[['vf_category', 'component']].to_dict('index')
for idx, row in genes_annotated.iterrows():
gene_name = row.get(match_field, '')
if gene_name and gene_name in virulence_map:
genes_annotated.at[idx, 'is_virulence'] = True
genes_annotated.at[idx, 'vf_category'] = virulence_map[gene_name]['vf_category']
genes_annotated.at[idx, 'vf_component'] = virulence_map[gene_name]['component']
# 策略 2: 关键词匹配 product 描述 (简化示例)
virulence_keywords = [
'toxin', 'adhesin', 'invasin', 'secretion', 'pilus', 'fimbria',
'capsule', 'lipopolysaccharide', 'hemolysin', 'protease',
'iron acquisition', 'siderophore', 'biofilm', 'quorum sensing'
]
for idx, row in genes_annotated.iterrows():
if not row['is_virulence'] and row.get('product'):
product_lower = row['product'].lower()
for keyword in virulence_keywords:
if keyword in product_lower:
genes_annotated.at[idx, 'is_virulence'] = True
genes_annotated.at[idx, 'vf_category'] = f'keyword_match:{keyword}'
break
vir_count = genes_annotated['is_virulence'].sum()
logger.info(f"✅ 注释完成: {vir_count}/{len(genes_annotated)} 基因标记为毒力相关")
return genes_annotated
def calculate_methylation_density(
sites_df: pd.DataFrame,
genes_df: pd.DataFrame,
flank: int = 0
) -> pd.DataFrame:
"""
计算每个基因的甲基化密度
参数:
flank: 基因上下游扩展范围 (bp)，用于包含调控区域
返回:
DataFrame: 每个基因的修饰计数和密度
"""
logger.info(f"📐 计算基因甲基化密度 (flank={flank}bp)")
results = []
for _, gene in genes_df.iterrows():
# 基因区域 ± flank
region_start = max(0, gene['start'] - flank)
region_end = gene['end'] + flank
region_length = region_end - region_start + 1
# 统计该区域内的修饰位点
mask = (
(sites_df['chrom'] == gene['chrom']) &
(sites_df['position'] >= region_start) &
(sites_df['position'] <= region_end)
)
gene_sites = sites_df[mask]
results.append({
'gene_id': gene['gene_id'],
'gene_name': gene.get('gene_name', ''),
'chrom': gene['chrom'],
'is_virulence': gene.get('is_virulence', False),
'vf_category': gene.get('vf_category'),
'gene_length': gene['length'],
'region_length': region_length,
'mod_count': len(gene_sites),
'mod_density_per_kb': len(gene_sites) / (region_length / 1000) if region_length > 0 else 0,
'avg_mod_level': gene_sites['mod_level'].mean() if len(gene_sites) > 0 else 0,
'avg_coverage': gene_sites['coverage'].mean() if len(gene_sites) > 0 else 0
})
return pd.DataFrame(results)
def test_virulence_enrichment(
gene_density_df: pd.DataFrame,
mod_type: str,
n_permutations: int = 1000,
random_seed: int = 42
) -> Dict:
"""
统计检验: 毒力基因是否富集/缺失甲基化修饰
方法:
1. Fisher 精确检验 (2×2 列联表)
2. 置换检验 (验证结果稳健性)
3. 效应量计算 (Cohen's d)
返回:
Dict: 统计结果 + 可视化数据
"""
logger.info(f"🧪 执行统计检验: {mod_type} 修饰与毒力基因关联")
np.random.seed(random_seed)
# 分组
virulence = gene_density_df[gene_density_df['is_virulence']]
non_virulence = gene_density_df[~gene_density_df['is_virulence']]
if len(virulence) == 0 or len(non_virulence) == 0:
logger.warning("⚠️ 某组基因数为 0，无法执行统计检验")
return {'error': 'insufficient_data'}
# === Fisher 精确检验 ===
# 构建 2×2 列联表: 高甲基化 (密度 > 中位数) vs 低甲基化
median_density = gene_density_df['mod_density_per_kb'].median()
contingency = pd.DataFrame({
'high_methylation': [
(virulence['mod_density_per_kb'] > median_density).sum(),
(non_virulence['mod_density_per_kb'] > median_density).sum()
],
'low_methylation': [
(virulence['mod_density_per_kb'] <= median_density).sum(),
(non_virulence['mod_density_per_kb'] <= median_density).sum()
]
}, index=['virulence', 'non_virulence'])
odds_ratio, fisher_p = fisher_exact(contingency.values, alternative='two-sided')
# === 置换检验 ===
observed_diff = virulence['mod_density_per_kb'].mean() - non_virulence['mod_density_per_kb'].mean()
permuted_diffs = []
all_densities = gene_density_df['mod_density_per_kb'].values
all_labels = gene_density_df['is_virulence'].values
for _ in range(n_permutations):
np.random.shuffle(all_labels)
perm_vir = all_densities[all_labels]
perm_non = all_densities[~all_labels]
if len(perm_vir) > 0 and len(perm_non) > 0:
permuted_diffs.append(perm_vir.mean() - perm_non.mean())
perm_p = (np.sum(np.abs(permuted_diffs) >= np.abs(observed_diff)) + 1) / (n_permutations + 1)
# === 效应量 ===
from scipy.stats import ttest_ind
t_stat, t_p = ttest_ind(
virulence['mod_density_per_kb'],
non_virulence['mod_density_per_kb'],
equal_var=False
)
cohen_d = observed_diff / np.sqrt(
(virulence['mod_density_per_kb'].std()**2 + non_virulence['mod_density_per_kb'].std()**2) / 2
)
results = {
'mod_type': mod_type,
'sample_sizes': {
'virulence_genes': len(virulence),
'non_virulence_genes': len(non_virulence)
},
'density_stats': {
'virulence_mean': round(virulence['mod_density_per_kb'].mean(), 4),
'virulence_std': round(virulence['mod_density_per_kb'].std(), 4),
'non_virulence_mean': round(non_virulence['mod_density_per_kb'].mean(), 4),
'non_virulence_std': round(non_virulence['mod_density_per_kb'].std(), 4)
},
'statistical_tests': {
'fisher_exact': {
'odds_ratio': round(odds_ratio, 3),
'p_value': fisher_p,
'significant': fisher_p < 0.05
},
't_test': {
't_statistic': round(t_stat, 3),
'p_value': t_p,
'significant': t_p < 0.05
},
'permutation_test': {
'observed_diff': round(observed_diff, 4),
'p_value': perm_p,
'n_permutations': n_permutations,
'significant': perm_p < 0.05
}
},
'effect_size': {
'mean_difference': round(observed_diff, 4),
'cohens_d': round(cohen_d, 3),
'interpretation': _interpret_cohens_d(cohen_d)
},
'contingency_table': contingency.to_dict(),
'permuted_diffs': permuted_diffs # 用于绘图
}
# 显著性总结
sig_count = sum([
results['statistical_tests']['fisher_exact']['significant'],
results['statistical_tests']['t_test']['significant'],
results['statistical_tests']['permutation_test']['significant']
])
results['overall_significance'] = 'significant' if sig_count >= 2 else 'not_significant'
logger.info(f"✅ 检验完成: {results['overall_significance']} (Fisher p={fisher_p:.4f})")
return results
def _interpret_cohens_d(d: float) -> str:
"""解释 Cohen's d 效应量"""
if abs(d) < 0.2:
return "negligible"
elif abs(d) < 0.5:
return "small"
elif abs(d) < 0.8:
return "medium"
else:
return "large"
def identify_candidate_virulence_modifications(
gene_density_df: pd.DataFrame,
stat_results: Dict,
top_n: int = 20
) -> pd.DataFrame:
"""
识别"高甲基化毒力基因"候选列表
筛选标准:
1. 是毒力基因 (is_virulence=True)
2. 修饰密度 > 全基因组 90% 分位数
3. 修饰水平 (mod_level) > 0.7
4. 覆盖度 ≥ 10×
"""
logger.info("🎯 识别候选高甲基化毒力基因...")
if gene_density_df.empty or 'is_virulence' not in gene_density_df.columns:
return pd.DataFrame()
# 筛选毒力基因
vir_genes = gene_density_df[gene_density_df['is_virulence']].copy()
if vir_genes.empty:
logger.warning("⚠️ 无注释的毒力基因")
return pd.DataFrame()
# 计算阈值
density_90 = gene_density_df['mod_density_per_kb'].quantile(0.9)
# 应用筛选
candidates = vir_genes[
(vir_genes['mod_density_per_kb'] >= density_90) &
(vir_genes['avg_mod_level'] >= 0.7) &
(vir_genes['avg_coverage'] >= 10) &
(vir_genes['mod_count'] >= 3) # 至少 3 个位点
].copy()
# 排序: 优先高密度 + 高修饰水平
candidates['priority_score'] = (
candidates['mod_density_per_kb'] / density_90 * 0.5 +
candidates['avg_mod_level'] * 0.3 +
np.log1p(candidates['mod_count']) * 0.2
)
candidates = candidates.sort_values('priority_score', ascending=False).head(top_n)
logger.info(f"✅ 识别 {len(candidates)} 个候选高甲基化毒力基因")
return candidates
# ============================================================================
# 可视化函数 (为 R 脚本准备数据)
# ============================================================================
def prepare_plot_data(
gene_density_df: pd.DataFrame,
stat_results: Dict,
candidates_df: pd.DataFrame,
output_dir: Path
) -> Dict[str, str]:
"""准备 R 绘图所需的中间数据文件"""
logger.info("📊 准备可视化数据...")
plot_files = {}
# 1. 密度分布箱线图数据
boxplot_data = gene_density_df[['is_virulence', 'mod_density_per_kb', 'mod_type']].copy()
boxplot_data['group'] = boxplot_data['is_virulence'].map({True: 'Virulence', False: 'Non-virulence'})
boxplot_path = output_dir / 'boxplot_data.tsv'
boxplot_data.to_csv(boxplot_path, sep='\t', index=False)
plot_files['boxplot'] = str(boxplot_path)
# 2. 候选基因列表 (用于标记)
if not candidates_df.empty:
candidate_plot = candidates_df[['gene_name', 'gene_id', 'mod_density_per_kb', 'avg_mod_level', 'vf_category']].copy()
candidate_path = output_dir / 'candidate_genes.tsv'
candidate_plot.to_csv(candidate_path, sep='\t', index=False)
plot_files['candidates'] = str(candidate_path)
# 3. 置换检验分布 (用于 p 值可视化)
if 'permuted_diffs' in stat_results:
perm_data = pd.DataFrame({
'permuted_diff': stat_results['permuted_diffs'],
'observed_diff': stat_results['statistical_tests']['permutation_test']['observed_diff']
})
perm_path = output_dir / 'permutation_data.tsv'
perm_data.to_csv(perm_path, sep='\t', index=False)
plot_files['permutation'] = str(perm_path)
# 4. 毒力类别汇总
if 'vf_category' in gene_density_df.columns:
vf_summary = gene_density_df[gene_density_df['is_virulence']].groupby('vf_category').agg({
'mod_density_per_kb': ['mean', 'std', 'count'],
'avg_mod_level': 'mean'
}).round(4)
vf_path = output_dir / 'vf_category_summary.tsv'
vf_summary.to_csv(vf_path, sep='\t')
plot_files['vf_summary'] = str(vf_path)
logger.info(f"✅ 生成 {len(plot_files)} 个绘图数据文件")
return plot_files
def generate_report(
sample: str,
mod_type: str,
stat_results: Dict,
candidates_df: pd.DataFrame,
output_dir: Path,
config: Dict
) -> str:
"""生成 Markdown 格式分析报告"""
report_path = output_dir / f"{sample}.{mod_type}.virulence_report.md"
with open(report_path, 'w', encoding='utf-8') as f:
f.write(f"# 🦠 甲基化 - 毒力基因关联分析报告\n")
f.write(f"**样本**: {sample} | **修饰类型**: {mod_type} | **日期**: {pd.Timestamp.now().strftime('%Y-%m-%d')}\n\n")
# 执行摘要
f.write("## 📋 执行摘要\n")
sig = stat_results.get('overall_significance', 'unknown')
f.write(f"- **关联显著性**: {'✅ 显著' if sig == 'significant' else '❌ 不显著'}\n")
if 'density_stats' in stat_results:
ds = stat_results['density_stats']
f.write(f"- **毒力基因平均密度**: {ds['virulence_mean']} sites/kb\n")
f.write(f"- **非毒力基因平均密度**: {ds['non_virulence_mean']} sites/kb\n")
if 'effect_size' in stat_results:
es = stat_results['effect_size']
f.write(f"- **效应量 (Cohen's d)**: {es['cohens_d']} ({es['interpretation']})\n")
f.write("\n")
# 统计检验结果
f.write("## 🧪 统计检验结果\n")
f.write("| 检验方法 | 统计量 | P 值 | 显著性 |\n")
f.write("|---------|--------|------|--------|\n")
tests = stat_results.get('statistical_tests', {})
if 'fisher_exact' in tests:
ft = tests['fisher_exact']
f.write(f"| Fisher 精确检验 | OR={ft['odds_ratio']} | p={ft['p_value']:.4f} | {'✅' if ft['significant'] else '❌'} |\n")
if 't_test' in tests:
tt = tests['t_test']
f.write(f"| t 检验 | t={tt['t_statistic']} | p={tt['p_value']:.4f} | {'✅' if tt['significant'] else '❌'} |\n")
if 'permutation_test' in tests:
pt = tests['permutation_test']
f.write(f"| 置换检验 | diff={pt['observed_diff']} | p={pt['p_value']:.4f} | {'✅' if pt['significant'] else '❌'} |\n")
f.write("\n")
# 候选基因列表
f.write("## 🎯 候选高甲基化毒力基因 (Top 20)\n")
if not candidates_df.empty:
f.write("| 基因名称 | 基因 ID | 修饰密度 (sites/kb) | 平均修饰水平 | 毒力类别 |\n")
f.write("|---------|---------|-------------------|-------------|----------|\n")
for _, row in candidates_df.iterrows():
name = row.get('gene_name', 'N/A') or row['gene_id']
vf_cat = row.get('vf_category', 'N/A') or 'unknown'
f.write(f"| {name} | {row['gene_id']} | {row['mod_density_per_kb']:.2f} | {row['avg_mod_level']:.3f} | {vf_cat} |\n")
else:
f.write("*无符合筛选标准的候选基因*\n")
f.write("\n")
# 生物学解读
f.write("## 🔬 生物学解读建议\n")
f.write("1. **限制修饰系统**: 检查候选基因是否位于已知 R-M 系统附近\n")
f.write("2. **基因调控**: 启动子区 6mA 富集可能抑制转录，建议整合 RNA-seq 验证\n")
f.write("3. **水平基因转移**: 毒力岛区域的异常甲基化模式可能指示近期获得\n")
f.write("4. **实验验证**: 优先选择 `priority_score` 最高的 3-5 个基因进行敲除/过表达实验\n")
f.write("\n")
# 输出文件列表
f.write("## 📁 输出文件清单\n")
f.write(f"- 统计结果: `{output_dir}/{sample}.{mod_type}.stats.json`\n")
f.write(f"- 候选基因: `{output_dir}/{sample}.{mod_type}.candidates.tsv`\n")
f.write(f"- 绘图数据: `{output_dir}/plots/` 目录\n")
f.write(f"- 完整报告: `{report_path}`\n")
logger.info(f"📄 报告保存: {report_path}")
return str(report_path)
# ============================================================================
# 主函数
# ============================================================================
def main():
parser = argparse.ArgumentParser(description='细菌甲基化 - 毒力基因关联分析')
parser.add_argument('--config', required=True, help='配置文件路径 (YAML)')
parser.add_argument('--sample', required=True, help='样本名称 (如: An6)')
parser.add_argument('--mods', nargs='+', default=['6mA', '4mC'], help='要分析的修饰类型')
parser.add_argument('--vfdb', help='VFDB 数据库路径 (覆盖 config 设置)')
parser.add_argument('--n-perm', type=int, default=1000, help='置换检验次数')
parser.add_argument('--top-n', type=int, default=20, help='候选基因数量')
args = parser.parse_args()
# 加载配置
with open(args.config, 'r') as f:
config = yaml.safe_load(f)
# 路径设置
sample = args.sample
input_dir = Path(config['input_dir'])
assembly_dir = Path(config['assembly_dir'])
outdir = Path(config['output_dir']) / sample / 'virulence_analysis'
outdir.mkdir(parents=True, exist_ok=True)
# 数据库路径
vfdb_path = args.vfdb or config.get('virulence', {}).get('vfdb_path')
gff_path = config['reference_genomes'].get(sample, '').replace('.fasta', '.gff').replace('.fa', '.gff')
results_summary = {}
for mod_type in args.mods:
logger.info(f"\n{'='*70}")
logger.info(f"🔬 分析: {sample} - {mod_type} 修饰与毒力基因关联")
logger.info(f"{'='*70}")
# === 1. 加载数据 ===
# 甲基化位点
bedmethyl_path = input_dir / f"{sample}.{mod_type}.bedMethyl.gz"
if not bedmethyl_path.exists():
logger.warning(f"⚠️ 甲基化文件不存在: {bedmethyl_path}")
continue
sites_df = load_methylation_sites(str(bedmethyl_path), mod_type)
# 基因注释
if not Path(gff_path).exists():
logger.warning(f"⚠️ GFF 文件不存在: {gff_path}, 跳过基因级分析")
continue
genes_df = load_gene_annotation(gff_path)
# VFDB 数据库
if vfdb_path and Path(vfdb_path).exists():
vfdb_df = load_vfdb_database(vfdb_path, species=sample)
genes_df = annotate_genes_with_virulence(genes_df, vfdb_df)
else:
logger.warning("⚠️ VFDB 数据库不可用，使用关键词匹配")
genes_df = annotate_genes_with_virulence(genes_df, pd.DataFrame())
# === 2. 计算密度 ===
gene_density_df = calculate_methylation_density(sites_df, genes_df, flank=200)
gene_density_df['mod_type'] = mod_type
# === 3. 统计检验 ===
stat_results = test_virulence_enrichment(
gene_density_df, mod_type,
n_permutations=args.n_perm
)
if 'error' in stat_results:
continue
# === 4. 识别候选基因 ===
candidates_df = identify_candidate_virulence_modifications(
gene_density_df, stat_results, top_n=args.top_n
)
# === 5. 准备绘图数据 ===
plot_dir = outdir / 'plots'
plot_dir.mkdir(exist_ok=True)
plot_files = prepare_plot_data(gene_density_df, stat_results, candidates_df, plot_dir)
# === 6. 保存结果 ===
# 详细统计 (JSON)
stats_path = outdir / f"{sample}.{mod_type}.stats.json"
# 移除大数组以减小文件
stats_output = {k: v for k, v in stat_results.items() if k != 'permuted_diffs'}
with open(stats_path, 'w') as f:
json.dump(stats_output, f, indent=2)
# 候选基因 (TSV)
if not candidates_df.empty:
candidates_path = outdir / f"{sample}.{mod_type}.candidates.tsv"
candidates_df.to_csv(candidates_path, sep='\t', index=False)
# 基因密度表 (完整结果)
density_path = outdir / f"{sample}.{mod_type}.gene_density.tsv.gz"
gene_density_df.to_csv(density_path, sep='\t', index=False, compression='gzip')
# === 7. 生成报告 ===
report_path = generate_report(
sample, mod_type, stat_results, candidates_df, outdir, config
)
# 汇总
results_summary[mod_type] = {
'significance': stat_results.get('overall_significance'),
'fisher_p': stat_results.get('statistical_tests', {}).get('fisher_exact', {}).get('p_value'),
'candidate_count': len(candidates_df),
'output_files': {
'stats': str(stats_path),
'candidates': str(candidates_path) if not candidates_df.empty else None,
'report': report_path,
'plots': str(plot_dir)
}
}
logger.info(f"✅ {mod_type} 分析完成")
# 保存总汇总
summary_path = outdir / f"{sample}.virulence_summary.json"
with open(summary_path, 'w') as f:
json.dump(results_summary, f, indent=2)
logger.info(f"\n🎉 {sample} 毒力关联分析完成!")
logger.info(f"📁 结果目录: {outdir}")
logger.info(f"📄 快速查看: cat {outdir}/*.virulence_report.md")
return 0
if __name__ == '__main__':
sys.exit(main())📈 脚本 5：05_expression_integration.R – 整合转录组数据 (可选)
#!/usr/bin/env Rscript
#
# 📈 05_expression_integration.R
# 整合甲基化与基因表达数据，分析修饰对毒力基因转录的影响
#
# 输入:
# - 甲基化密度表 (gene_density.tsv.gz)
# - RNA-seq 表达矩阵 (TPM/FPKM)
# - 毒力基因注释
#
# 输出:
# - 相关性分析结果
# - 甲基化 - 表达散点图
# - 差异表达毒力基因的甲基化特征
suppressPackageStartupMessages({
library(ggplot2)
library(dplyr)
library(tidyr)
library(cowplot)
library(pheatmap)
})
# === 配置 ===
args <- commandArgs(trailingOnly = TRUE)
config <- list(
sample = if (length(args) >= 2 && args[1] == "--sample") args[2] else "An6",
mods = c("6mA", "4mC"),
expression_file = "rna_seq/An6_gene_expression.tsv", # 用户需提供
output_dir = "postprocess_results/An6/virulence_analysis/expression_integration"
)
dir.create(config$output_dir, recursive = TRUE, showWarnings = FALSE)
# === 辅助函数 ===
# 读取甲基化 - 基因密度表
read_gene_density <- function(filepath, mod) {
df <- read.table(filepath, sep = "\t", header = TRUE,
stringsAsFactors = FALSE, comment.char = "#")
df <- df[df$mod_type == mod, ]
message(sprintf("📥 读取 %s: %d 基因", filepath, nrow(df)))
return(df)
}
# 读取表达数据 (支持 TPM/FPKM)
read_expression <- function(filepath) {
if (!file.exists(filepath)) {
warning(sprintf("⚠️ 表达文件不存在: %s", filepath))
return(NULL)
}
expr <- read.table(filepath, sep = "\t", header = TRUE,
row.names = 1, stringsAsFactors = FALSE, comment.char = "#")
message(sprintf("📥 读取表达矩阵: %d 基因 × %d 样本", nrow(expr), ncol(expr)))
return(expr)
}
# 计算 Spearman 相关性 + 可视化
plot_methylation_expression <- function(meth_df, expr_vec, mod_type, gene_filter = NULL) {
# 合并数据
df <- data.frame(
gene_id = meth_df$gene_id,
meth_density = meth_df$mod_density_per_kb,
expression = expr_vec[meth_df$gene_id],
is_virulence = meth_df$is_virulence,
stringsAsFactors = FALSE
)
# 过滤低表达基因 (可选)
if (!is.null(gene_filter)) {
df <- df[df$expression >= gene_filter, ]
}
# 移除 NA
df <- df[complete.cases(df), ]
if (nrow(df) < 10) {
message("⚠️ 有效数据点过少，跳过绘图")
return(NULL)
}
# 计算相关性
cor_test <- cor.test(df$meth_density, df$expression, method = "spearman")
# 绘图
p <- ggplot(df, aes(x = meth_density, y = expression,
color = is_virulence, alpha = is_virulence)) +
geom_point(size = 1.5) +
scale_color_manual(values = c("TRUE" = "#E41A1C", "FALSE" = "#999999"),
labels = c("TRUE" = "Virulence", "FALSE" = "Other"),
name = "毒力基因") +
scale_alpha_manual(values = c("TRUE" = 1, "FALSE" = 0.4), guide = "none") +
geom_smooth(method = "lm", se = TRUE, color = "black", linewidth = 0.3, alpha = 0.1) +
labs(
x = sprintf("%s 修饰密度 (sites/kb)", mod_type),
y = "基因表达水平 (TPM)",
title = sprintf("甲基化 - 表达相关性 (%s)", mod_type),
subtitle = sprintf("Spearman ρ = %.3f, p = %.2e", cor_test$estimate, cor_test$p.value)
) +
theme_cowplot(font_size = 10) +
theme(
legend.position = "top",
panel.grid.minor = element_blank()
)
# 添加显著性标记
if (cor_test$p.value < 0.05) {
p <- p + annotate("text", x = Inf, y = Inf,
label = if (cor_test$p.value < 0.001) "***" else
if (cor_test$p.value < 0.01) "**" else "*",
hjust = 1.1, vjust = 1.5, size = 5, color = "red")
}
return(list(plot = p, correlation = cor_test$estimate, p_value = cor_test$p.value))
}
# === 主分析函数 ===
analyze_expression_integration <- function(sample, config) {
message(sprintf("\n📈 整合分析: %s 甲基化 + 表达", sample))
outdir <- file.path(config$output_dir)
results <- list()
for (mod in config$mods) {
message(sprintf("\n🔬 分析修饰: %s", mod))
# 1. 读取甲基化数据
meth_path <- file.path(dirname(config$output_dir),
sprintf("%s.%s.gene_density.tsv.gz", sample, mod))
if (!file.exists(meth_path)) {
warning(sprintf("⚠️ 甲基化密度文件不存在: %s", meth_path))
next
}
meth_df <- read_gene_density(meth_path, mod)
# 2. 读取表达数据
expr_df <- read_expression(config$expression_file)
if (is.null(expr_df)) next
# 3. 确保基因 ID 匹配 (简化: 假设 gene_id 一致)
common_genes <- intersect(meth_df$gene_id, rownames(expr_df))
if (length(common_genes) < 50) {
warning(sprintf("⚠️ 共同基因过少: %d", length(common_genes)))
next
}
meth_sub <- meth_df[meth_df$gene_id %in% common_genes, ]
expr_vec <- expr_df[common_genes, 1] # 取第一列/样本
names(expr_vec) <- common_genes
# 4. 全局相关性分析
message(" 📊 计算全局相关性...")
cor_result <- plot_methylation_expression(meth_sub, expr_vec, mod)
if (!is.null(cor_result)) {
# 保存图
plot_path <- file.path(outdir, sprintf("%s.%s.meth_expr_correlation.pdf", sample, mod))
ggsave(plot_path, cor_result$plot, width = 7, height = 5, dpi = 300)
message(sprintf(" ✓ 保存: %s", basename(plot_path)))
# 5. 毒力基因亚组分析
vir_df <- meth_sub[meth_sub$is_virulence, ]
if (nrow(vir_df) >= 10) {
message(" 🦠 分析毒力基因亚组...")
vir_cor <- plot_methylation_expression(vir_df, expr_vec[vir_df$gene_id], mod)
if (!is.null(vir_cor)) {
vir_plot_path <- file.path(outdir, sprintf("%s.%s.virulence_subset.pdf", sample, mod))
ggsave(vir_plot_path, vir_cor$plot, width = 6, height = 4.5, dpi = 300)
message(sprintf(" ✓ 保存: %s", basename(vir_plot_path)))
}
}
# 6. 记录结果
results[[mod]] <- list(
global_correlation = list(rho = cor_result$correlation, p = cor_result$p_value),
virulence_genes = nrow(meth_sub[meth_sub$is_virulence, ]),
plots = c(global = plot_path)
)
}
}
# 保存汇总
if (length(results) > 0) {
summary_path <- file.path(outdir, sprintf("%s.expression_integration_summary.json", sample))
jsonlite::write_json(results, summary_path, pretty = TRUE, auto_unbox = TRUE)
message(sprintf("\n✅ 表达整合分析完成: %s", summary_path))
}
return(invisible(results))
}
# === 命令行入口 ===
if (!interactive()) {
analyze_expression_integration(config$sample, config)
}🗄️ 数据文件：data/vfdb_core.tsv (简化示例)
# VFDB 核心毒力因子数据库 (简化版，实际请从 http://www.mgc.ac.cn/VFDB/ 下载完整版)
Gene_ID Gene_name Product VF_category Component Species
VF0001 fimA type 1 fimbrial protein adhesion fimbria/pilus Escherichia coli
VF0002 papG P fimbrial adhesin adhesion fimbria/pilus Escherichia coli
VF0003 hlyA alpha-hemolysin toxin hemolysin Escherichia coli
VF0004 ompA outer membrane protein A adhesion outer membrane Acinetobacter baumannii
VF0005 bauA acinetobactin receptor iron_acquisition siderophore Acinetobacter baumannii
VF0006 csuc colonic acid capsule biofilm capsule Acinetobacter baumannii
VF0007 adeB RND efflux pump transporter antibiotic_resistance efflux_pump Acinetobacter baumannii
VF0008 plcH hemolytic phospholipase C toxin phospholipase Pseudomonas aeruginosa
VF0009 lasR quorum sensing regulator regulation quorum_sensing Pseudomonas aeruginosa
VF0010 sph2 sphingomyelinase toxin phospholipase Pedobacter cryoconitis
# ... 实际数据库包含 3000+ 条目💡 获取完整 VFDB:
# 下载最新 VFDB 核心数据集 wget -O data/VFDB_setA_pro.fas.gz "http://www.mgc.ac.cn/VFDB/down/VFDB_setA_pro.fas.gz" # 或使用 API 获取 TSV 格式 curl "http://www.mgc.ac.cn/cgi-bin/VFDB/req.cgi?func=download&format=tsv&set=A" > data/vfdb_core.tsv
📝 分析报告模板：templates/virulence_report.md
# 🦠 {sample} 甲基化 - 毒力基因关联分析
> **分析日期**: {date}
> **修饰类型**: {mod_types}
> **参考基因组**: {reference}
> **VFDB 版本**: {vfdb_version}
---
## 🔑 关键发现
{# 自动填充显著性结果 #}
{% if significant %}
✅ **{mod_type} 修饰在毒力基因中显著富集** (Fisher p = {fisher_p:.4f})
- 毒力基因平均密度: **{vir_mean:.2f} sites/kb**
- 非毒力基因平均密度: {nonvir_mean:.2f} sites/kb
- 效应量 (Cohen's d): {cohens_d} ({effect_interpretation})
{% else %}
❌ 未发现 {mod_type} 修饰与毒力基因的显著关联
{% endif %}
---
## 🎯 候选高甲基化毒力基因
| 排名 | 基因名称 | 功能描述 | 修饰密度 | 毒力类别 | 优先级 |
|------|---------|---------|---------|---------|--------|
{% for gene in candidates %}
| {{ loop.index }} | {{ gene.gene_name }} | {{ gene.product|truncate(50) }} | {{ gene.mod_density_per_kb|round(2) }} | {{ gene.vf_category }} | ⭐⭐⭐ |
{% endfor %}
> 💡 **实验建议**: 优先验证优先级最高的 3 个基因:
> 1. `{{ candidates[0].gene_name }}`: {{ candidates[0].product }}
> 2. `{{ candidates[1].gene_name }}`: {{ candidates[1].product }}
> 3. `{{ candidates[2].gene_name }}`: {{ candidates[2].product }}
---
## 📊 可视化结果
### 甲基化密度分布
### 候选基因基因组位置
### 甲基化 - 表达相关性 (如有 RNA-seq)
{% if expression_integrated %}
{% endif %}
---
## 🔬 生物学解读
### 1. 限制修饰系统 (R-M) 关联
{# 分析候选基因是否位于已知 R-M 系统附近 #}
### 2. 基因表达调控假设
- 启动子区 6mA 富集 → 可能抑制转录 (参考文献: [1,2])
- 基因体 4mC 修饰 → 可能影响转录延伸
### 3. 进化与水平基因转移
- 毒力岛区域的异常甲基化模式可能指示近期获得
- 建议: 比对近缘菌株，分析甲基化模式保守性
---
## 📁 输出文件清单
{output_dir}/ ├── {sample}.{mod_type}.stats.json # 详细统计结果 ├── {sample}.{mod_type}.candidates.tsv # 候选基因列表 ├── {sample}.{mod_type}.virulence_report.md # 本报告 ├── plots/ │ ├── {sample}.{mod_type}.boxplot.png # 密度分布图 │ ├── {sample}.{mod_type}.candidates_track.png # 基因组轨道 │ └── … └── expression_integration/ # (可选) 表达整合分析
---
## 📚 参考文献
1. Fang G, et al. (2012) Genome-wide mapping of methylated adenine residues in pathogenic Escherichia coli using single-molecule real-time sequencing. *Nat Biotechnol*.
2. Blow MJ, et al. (2016) The epigenomic landscape of prokaryotes. *PLoS Genet*.
3. Chen L, et al. (2023) VFDB 2024: a curated database for bacterial virulence factors. *Nucleic Acids Res*.
---
*报告由 `04_virulence_analysis.py` 自动生成 | 版本 1.0*🚀 更新主入口脚本：run_all.sh (添加毒力分析)
# === 在 run_all.sh 的 run_postprocess() 函数中添加 ===
run_postprocess() {
local sample=$1
log_info "\n🔬 处理样本: ${sample}"
echo "━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━"
# 原有步骤 1-3...
# 4. 🦠 毒力基因关联分析 (新增)
log_info "🦠 步骤 4/4: 毒力基因关联分析..."
python3 "${SCRIPT_DIR}/04_virulence_analysis.py" \
--config "${CONFIG}" \
--sample "${sample}" \
--mods "${MODS[@]}" \
--n-perm 1000 \
--top-n 20
# 5. 📈 表达整合 (可选，如有 RNA-seq)
if [[ -f "rna_seq/${sample}_gene_expression.tsv" ]]; then
log_info "📈 步骤 5/4: 整合转录组数据..."
Rscript "${SCRIPT_DIR}/05_expression_integration.R" \
--config "${CONFIG}" \
--sample "${sample}"
else
log_warn "⏭️ 跳过表达整合 (未检测到 RNA-seq 数据)"
fi
log_success "✅ ${sample} 全部分析完成!\n"
}🇩🇪 德语快速参考（Deutsch）
# 1. 准备数据
# - 下载 VFDB: http://www.mgc.ac.cn/VFDB/
# - 确保 GFF 注释文件存在: assembly.fasta → assembly.gff
# 2. 运行毒力分析
python 04_virulence_analysis.py \
--config config.yaml \
--sample An6 \
--mods 6mA 4mC \
--vfdb data/vfdb_core.tsv
# 3. 查看结果
cat postprocess_results/An6/virulence_analysis/*.virulence_report.md
# 4. (可选) 整合表达数据
# 准备: rna_seq/An6_gene_expression.tsv (gene_id \t TPM)
Rscript 05_expression_integration.R --sample An6
# 5. 发表级输出
# - 候选基因表: *.candidates.tsv
# - 统计报告: *.stats.json
# - 图表: plots/*.pdf
# - 完整报告: *.virulence_report.md🎯 细菌特异性分析要点
| 分析目标 | 推荐方法 | 生物学意义 |
|---|---|---|
| 限制修饰系统鉴定 | 检查候选基因 ±5kb 内是否有 hsdR/M/S, dcm, dam |
R-M 系统常通过 6mA/4mC 保护自身 DNA |
| 毒力调控网络 | 整合 ChIP-seq (如 H-NS, Fis) + 甲基化数据 | 核蛋白与甲基化协同调控毒力基因簇 |
| 水平基因转移检测 | 比较甲基化模式在基因组岛 (GI) 内外的差异 | 新获得毒力岛可能显示异常甲基化 |
| 菌株比较 | 多菌株甲基化密度 + 毒力表型关联 | 识别与致病性相关的表观标记 |
⚠️ 注意事项与最佳实践
- 物种特异性: VFDB 主要覆盖病原菌，环境菌如 Pedobacter 可能需要自定义毒力基因列表
- 注释质量: 确保 GFF 文件包含准确的
gene特征和product描述 - 多重检验校正: 分析多个修饰类型时，对 p 值进行 FDR 校正 (
p.adjust(method="BH")) - 实验验证: 计算结果仅为预测，关键候选基因需通过:
- 甲基化特异性 PCR (MSP)
- 基因敲除 + 毒力表型测定
- 报告基因融合验证启动子活性
🔄 扩展建议
# 🔬 添加: 操纵子水平分析
# - 使用 DOOR2 数据库识别操纵子
# - 计算整个操纵子的平均甲基化密度
# - 检验"毒力操纵子"是否整体高甲基化
# 🌐 添加: 泛基因组甲基化比较
# - 整合多菌株甲基化数据
# - 构建"核心/附属"基因组甲基化图谱
# - 识别菌株特异性表观标记
# 🤖 添加: 机器学习预测
# - 特征: 位置/基序/保守性/表达/甲基化
# - 目标: 预测"调控性"甲基化位点
# - 模型: 随机森林 / XGBoost需要我帮您：
- 🔄 适配脚本到您的具体文件路径和集群环境？
- 📊 添加更多细菌功能注释 (如: 抗生素耐药基因、代谢通路)？
- 🐳 创建 Dockerfile 以便在任意环境复现完整分析流程？
