Unicycler+medaka_consensus+Polypolish/Pilon manually

Leave a reply

在 Oxford Nanopore (ONT) 的测序技术和 Medaka 模型命名中，e82 实际上是 e8.2 的简写。

它代表的是纳米孔蛋白的工程化迭代版本（Pore Engineering Version）。

为了让你完全看懂这串像密码一样的模型名称（例如 r1041_e82_400bps_sup），我们可以把它拆解成 4 个核心部分：

1. ONT 命名公式拆解

缩写 完整含义 通俗解释（代表什么？）
r1041 R10.4.1 纳米孔的“物理结构”。代表第 10 代孔蛋白的 4.1 亚型。它的孔径比早期的 R9 更长，能同时容纳更多碱基，从而大幅提高了读取准确度。
e82 e8.2 纳米孔的“化学/工程微调”。代表孔蛋白的第 8 大代、第 2 次工程化微调版本（Engineering version）。
400bps 400 bases per second 马达蛋白的“测序速度”。DNA 以每秒 400 个碱基的速度穿过纳米孔（早期版本是 450bps，故意降到 400bps 是为了让电信号采样更密集、更清晰）。
sup Super Accuracy Basecalling 神经网络的“精度等级”。SUP 是目前最高精度的标准模型（比 HAC 更准，但计算更慢）。

2. e8.2 具体解决了什么问题？

ONT 的科学家在研发 R10.4.1 孔蛋白时，发现仅仅改变物理孔径还不够。为了让电信号特征更明显，他们会对孔蛋白内部的氨基酸序列或化学修饰进行微调（这就是 engineering 的含义）。

  • e8.1 是早期的工程版本。
  • e8.2 是目前 R10.4.1 芯片上最成熟、最广泛使用的标准版本。它专门优化了均聚物（homopolymers，比如连续的 AAAAA 或 CCCCC） 和相似碱基的电信号分辨率，极大减少了 Indel（插入/缺失）错误。 因此，只要诺禾致源的报告上写了 R10.4.1 芯片，它在化学层面上对应的就是 e8.2 版本。

3. 为什么 Medaka 里写的是 e82 而不是 e8.2？

这只是 Medaka 开发者为了避免文件名和命令行中出现多个小数点（.）导致系统解析错误，而人为去掉了小数点，把 e8.2 简写成了 e82。 同理，R10.4.1 被简写成了 r1041。

4. 结合你的项目（诺禾致源报告）的最终确认

根据你提供的报告信息：

  • 芯片：FLO-PRO114M (R10.4.1) 👉 对应 r1041 和 e82
  • 速度：400 bp/s 👉 对应 400bps
  • 精度：SUP (超高精度) 👉 对应 sup
  • 版本：Medaka 默认推荐最新版的 v4.3.0 👉 对应 v4.3.0

结论： 你的数据完美匹配 r1041_e82_400bps_sup_v4.3.0（或者你列表里的 r1041_e82_400bps_sup_g632，两者在很多 Medaka 版本中是等价的，推荐用带 v4.3.0 的最新版）。

⚠️ 再次强调关于“甲基化 (Methylation)”的避坑指南

正如你提供的文本中所提醒的，这是一个极其重要的实操细节：

如果你的文件夹叫 Methylation，说明诺禾致源在 Basecalling（碱基识别）时用了甲基化模型（比如 ...sup@v4.2.0_5mCG_5hmC）。这个模型输出的文件里包含了 5mC/5hmC 的修饰标签。

但是，Medaka 的 consensus polishing 神经网络是“瞎子”，它只认识标准的 A/T/C/G，不认识甲基化标签。

  • 如果你把带有甲基化标签的数据喂给 Medaka：它可能会直接报错，或者把甲基化信号误认为是测序噪音，导致纠错失败甚至引入错误。
  • 正确做法：在跑 Medaka 之前，确保你输入的 .fastq 文件是用标准 DNA 模型（不带 5mCG 后缀）Basecalling 出来的。如果你只有 .pod5 原始电信号文件，你需要用 Dorado 或 Guppy 的标准 SUP 模型重新 Basecall 一遍，生成纯粹的 A/T/C/G fastq 文件，然后再丢给 Medaka 去抛光。

Don’t recommend running Pilon + Polypolish back-to-back.** Here’s why, and what to do instead.

Why chaining them is a bad idea

Pilon and Polypolish are both short-read polishers solving the same problem. They’re redundant, not complementary:

  • Pilon aggressively “fixes” bases, indels, and even gaps — but it can introduce misassemblies in repeat regions (it maps reads into repeats and “corrects” sequence incorrectly).
  • Polypolish (by the Unicycler author, Ryan Wick) was built specifically to be more conservative and safer around repeats — it only changes positions supported by reads and doesn’t touch structure.

If you chain them, you get diminishing returns and risk them fighting each other: e.g., Pilon introduces a bad change in a repeat, and conservative Polypolish can’t revert it (it needs read support). Or Pilon runs after Polypolish and re-breaks the repeat regions Polypolish left safe.

Recommended approach: pick ONE (Polypolish is the safer choice here)

Given your case — a bacterial genome with IS elements/repeats, destined for SNP analysis + NCBI submission — use Polypolish only:

cd adeAB_unicycler_conservative

# Index the draft
bwa index assembly.fasta

# Map short reads with short k-mer (lets reads span repeats)
bwa mem -t 32 -k 15 assembly.fasta \
    ../raw_data/19606adeAB_R1.fastq.gz \
    ../raw_data/19606adeAB_R2.fastq.gz | samtools sort -o map.bam -

# Filter + polish
polypolish filter --in map.bam --out-fwd fwd.sam --out-bwd bwd.sam
polypolish polish \
    --draft assembly.fasta \
    --in-fwd fwd.sam --in-bwd bwd.sam \
    --out polypolished.fasta

This is cleaner than Pilon (no Java dependency — remember your trycycler env was missing Pilon) and safer for repeats.

If you genuinely want extra polishing rounds

Run the same tool twice rather than mixing two tools — but with Polypolish, one round is usually enough. If you insist on maximum accuracy for SNP calling, the more valuable addition is a long-read polish (Medaka) before the short-read polish, because that’s genuinely complementary:

# Optional complementary step: long-read polish first
medaka_consensus -i long_reads.fastq.gz -d assembly.fasta \
    -m r1041_e82_400bps_sup_v5.0.0 -t 32 -o medaka_out
# THEN short-read polish the medaka output with Polypolish

Always verify after polishing

quast.py assembly.fasta polypolished.fasta -o polish_compare
# Confirm: genome size still ~3.9 Mb, circularity preserved, no new fragmentation

⚠️ And as before — since the ΔadeAB assembly had that ~129 kb gap / circularity issue, make sure you polish the corrected assembly, not the gapped one.

Bottom line: Use Polypolish alone (optionally preceded by Medaka for long reads). Don’t stack Pilon + Polypolish. Want me to write a single drop-in script that does Medaka → Polypolish → QUAST for all your strains in a loop?

下面给你一套可直接使用的 mamba 环境配置，覆盖 Unicycler（混合组装）→ Medaka（长读长抛光）→ Polypolish（短读长抛光）的完整流程，并附带比对、质控和可视化工具。

1. 创建 environment.yml

新建一个文件 environment.yml，粘贴以下内容：

name: hybrid_polish
channels:
  - conda-forge
  - bioconda
  - defaults
dependencies:
  # ===== 混合组装 =====
  - unicycler>=0.5.0

  # ===== 长读长抛光 (ONT) =====
  - medaka
  - racon

  # ===== 短读长抛光 (Illumina) =====
  - polypolish
  - pilon

  # ===== 比对工具 =====
  - minimap2>=2.24
  - bwa
  - samtools>=1.17
  - bcftools

  # ===== 质控 (QC) =====
  - fastqc
  - multiqc
  - fastp
  - nanofilt
  - nanostat

  # ===== 组装评估 / 可视化 =====
  - quast
  - bandage
  - ragtag

  # ===== 辅助工具 =====
  - seqtk
  - bedtools

💡 我没有把 circlator 放进去，因为它依赖较老、容易和 unicycler/medaka 产生冲突。如果你后续确实需要环化，可以单独建一个环境装它（见第 5 节）。

2. 创建并激活环境

# 创建环境（mamba 会自动解决依赖）
mamba env create -f environment.yml

# 激活环境
mamba activate hybrid_polish

如果 mamba env create 解析依赖很慢或失败，可以改用显式安装（更直观）：

mamba create -n hybrid_polish -c conda-forge -c bioconda \
  unicycler medaka racon polypolish pilon \
  minimap2 bwa samtools bcftools \
  fastqc multiqc fastp nanofilt nanostat \
  quast bandage ragtag seqtk bedtools

mamba activate hybrid_polish

3. 验证安装

unicycler --version
medaka --version          # 或: medaka_consensus --help
polypolish --version
pilon --version
minimap2 --version
samtools --version
racon --version

4. ⚠️ Medaka 模型注意事项（关键）

medaka 安装后不会自动带模型，模型需要匹配你的测序化学（flowcell + basecaller）。先列出可用模型：

medaka tools list_models

常见模型对照：

测序化学 推荐模型
R10.4.1 + SUP (v5) r1041_e82_400bps_sup_v5.0.0
R10.4.1 + HAC (v5) r1041_e82_400bps_hac_v5.0.0
R9.4.1 + SUP r941_prom_sup_g506
R9.4.1 + HAC r941_prom_hac_g507

如果 medaka_consensus 报模型缺失，它通常会在首次运行时自动下载；若处于离线环境，需手动 medaka tools download_models。

5. 完整工作流示例（Unicycler → Medaka → Polypolish）

# --- Step 1: 混合组装 (Unicycler) ---
unicycler \
  -1 R1.fq.gz -2 R2.fq.gz \
  -l long_reads.fq.gz \
  -o unicycler_out -t 32

# --- Step 2: 长读长抛光 (Medaka) ---
medaka_consensus \
  -i long_reads.fq.gz \
  -d unicycler_out/assembly.fasta \
  -m r1041_e82_400bps_sup_v5.0.0 \
  -t 32 \
  -o medaka_out

# --- Step 3: 短读长抛光 (Polypolish) ---
bwa mem -t 32 -k 15 medaka_out/consensus.fasta R1.fq.gz R2.fq.gz | samtools sort -o map.bam -
samtools index map.bam

polypolish filter --in map.bam --out-fwd fwd.sam --out-bwd bwd.sam
polypolish polish \
  --draft medaka_out/consensus.fasta \
  --in-fwd fwd.sam --in-bwd bwd.sam \
  --out final_polished.fasta

# --- Step 4: 质量评估 (QUAST) ---
quast.py final_polished.fasta -o quast_report

6. 如果遇到依赖冲突

  • medaka 冲突：medaka 对 Python/TensorFlow 版本敏感。如果和 unicycler 冲突，可单独建环境： 
    mamba create -n medaka_env -c conda-forge -c bioconda medaka
  • circlator 需要时： 
    mamba create -n circlator_env -c conda-forge -c bioconda circlator

需要我把这套流程写成一个可复用的 bash 脚本（自动判断样本名、循环处理多个菌株），或者帮你适配成 Nextflow/bacass 的 config 吗？

Leave a Reply

Your email address will not be published. Required fields are marked *