Unicycler + Medaka (HAC) + Polypolish manually

Leave a reply

Table 1: Structural Variants in Complete Long-Read Assemblies (vs. CP059040)

Strain Variant Type Size (bp) Reference Coordinates (CP059040) Classification / Biological Significance
19606WT Tandem Contraction 52,668 2,810,861 – 2,863,470 Background variation (inherent to the 19606 lab strain)
Tandem Contraction 198 3,124,916 – 3,125,037 Background variation
19606_adeAB Deletion 4,282 1,844,323 – 1,848,605 Targeted ΔadeAB knockout
Tandem Contraction 52,668 2,810,861 – 2,863,470 Background variation
Tandem Contraction 198 3,124,916 – 3,125,037 Background variation
19606_adeIJ Deletion 4,436 737,224 – 741,667 Targeted ΔadeIJ knockout
Tandem Contraction 52,668 2,810,861 – 2,863,470 Background variation
Tandem Contraction 198 3,124,916 – 3,125,037 Background variation

抓到了！真相大白。这是一个非常关键的发现！

从你提取出的 fastq header 信息中，我们看到了真实的 Basecalling 模型： 👉 dna_r10.4.1_e8.2_400bps_hac@v5.2.0_barcode12

⚠️ 关键纠正：是 HAC，不是 SUP！

之前我们根据诺禾致源的常规细菌完成图经验，推测他们使用的是超高精度的 SUP (Super Accuracy) 模型。但实际数据证明，他们这批数据使用的是 HAC (High Accuracy) 模型，版本正是 v5.2.0。

在 Medaka 中，模型必须与 Basecalling 模型严格匹配。 如果你用 SUP 模型去抛光 HAC 的数据，神经网络会因为错误特征不匹配而导致抛光效果变差，甚至引入错误。

🎯 最终 100% 匹配的 Medaka 模型

根据你之前 medaka tools list_models 的列表，完美对应的模型是： r1041_e82_400bps_hac_v5.2.0

🛠️ 最终正确的 Medaka 运行命令

现在你可以放心地使用以下命令对你的 Unicycler 组装结果（特别是需要修复的 $\Delta$adeAB）进行长读长抛光了：

# 确保你在 medaka 或 hybrid_polish 环境中
mamba activate medaka 

# 运行 Medaka (以 adeAB 为例)
medaka_consensus \
    -i merged_19606_adeAB_longreads.fastq.gz \
    -d 19606_adeAB-unicycler.scaffolds.fa \
    -m r1041_e82_400bps_hac_v5.2.0 \
    -t 32 \
    -o medaka_adeAB_hac_v520

💡 关于 HAC 模型的补充说明（不用担心精度）

你可能会担心 HAC 精度不如 SUP。实际上：

  1. 对于细菌基因组：R10.4.1 芯片的 HAC 模型精度已经极高（Q20+，准确率 >99%），对于组装和寻找大片段结构变异（SV）完全足够。
  2. 你有短读长兜底：Medaka (长读长) 跑完后，你还会用 Polypolish / Pilon (Illumina 短读长) 进行最后一步抛光。短读长会把 HAC 残留的极少量单碱基错误（特别是 homopolymer 区域）彻底修正。
  3. 最终结果：Unicycler + Medaka (HAC) + Polypolish 产出的基因组，质量与使用 SUP 模型几乎没有区别，完全满足 NCBI 提交和 SNP 分析的要求。

现在参数已经完全锁定，你可以开始跑 Medaka 了！跑完后记得用 quast.py 检查一下组装质量，然后再上 Polypolish。

不是的（No）。 这是一个非常容易混淆的概念，让我为您彻底理清 .pod5 文件和“带甲基化标签的序列”之间的区别。

简单来说：.pod5 是“原始原料”，而“带甲基化标签”是指“加工后的产品”。

1. 概念澄清

文件类型 本质是什么？ 包含甲基化信息吗？ 能给 Medaka 用吗？
.pod5 原始电信号文件 (Raw Signal)。相当于测序仪记录下的“未经翻译的电流波动图”。 天然包含。因为甲基化修饰会改变电流信号，这些物理信号都原封不动地保存在 .pod5 里。 ❌ 绝对不能。Medaka 不认识电信号，它只认识 A/T/C/G 字母。
.fastq / .fastq.gz 碱基序列文件。是 Basecalling 软件（如 Dorado）把 .pod5 的电信号“翻译”成的 A/T/C/G 字母。 取决于翻译时用的模型。如果用“标准模型”翻译，就只有 ATCG；如果用“甲基化模型”翻译，可能会带有修饰信息。 ✅ 可以。Medaka 需要的是标准的 .fastq.gz。
.bam (带 tags) 比对后的序列文件。诺禾致源做甲基化分析时，通常会交付带有 MM/ML 标签的 .bam 文件，专门用来标记哪个 C 碱基被甲基化了（5mC/5hmC）。 包含。这就是所谓的“带甲基化标签的数据”。 ❌ 不能。Medaka 处理不了这些复杂的修饰标签。

2. 诺禾致源（Novogene）交付数据的真实情况

在诺禾致源交付给您的文件夹中，通常是这样的结构：

  1. 01.RawData (或类似目录) 里的 .pod5 文件：
    • 这是最原始的电信号。无论您后续想做标准基因组组装，还是想做甲基化分析，源头都是这一套 .pod5 文件。
  2. 02.CleanData (或 Basecalling 目录) 里的 .fastq.gz 文件：
    • 这是诺禾致源用标准 DNA 模型（如 dna_r10.4.1_e8.2_400bps_sup）对 .pod5 进行 Basecalling 后生成的纯净序列。这套文件只包含 A/T/C/G，没有甲基化标签。
    • 👉 这正是您用来跑 Unicycler 组装和 Medaka 抛光的文件！
  3. 甲基化分析专属目录（如果有）里的 .bam 或特定报告：
    • 这是诺禾致源用甲基化模型（如 ..._5mCG_5hmC）重新处理 .pod5 后生成的，专门用于后续分析甲基化位点。

3. 您现在该怎么做？（实操指南）

您不需要去管 .pod5 文件（除非您想自己用 Dorado 重新跑一遍 Basecalling，这非常耗时且没必要）。

您只需要确认一件事： 您输入给 Medaka 的 .fastq.gz 文件，是不是诺禾致源交付的标准 Clean Data？

如果是，那它里面只有纯粹的 A/T/C/G，完全没有任何甲基化标签的干扰，您可以放心大胆地直接用它跑 Medaka！

验证您的 fastq 是否“纯净”的极简方法：

您可以随便打开一个诺禾致源交付的 .fastq.gz 文件的前几行看看：

zcat merged_19606_adeAB_longreads.fastq.gz | head -n 8

正常的输出应该是这样的（只有 ATCG 和.quality 分数）：

@read_id_1
ATCGATCGATCGATCG...
+
IIIIIIIIIIIIIIII...

只要序列行（第二行）里只有 A, T, C, G（偶尔有 N），没有乱七八糟的化学修饰符号，它就是标准的 DNA 序列，完美适配 Medaka。

总结

  • .pod5 = 包含一切信号的原始电报（不能直接给 Medaka）。
  • 诺禾交付的 .fastq.gz = 翻译好的标准 ATCG 字母（直接拿来跑 Medaka，非常安全）。
  • 带甲基化标签的数据 = 诺禾用特殊模型翻译并打上标记的 .bam 文件（用于甲基化分析，别喂给 Medaka）。

所以，直接用您跑 Unicycler 时用的那个 merged_...longreads.fastq.gz 去跑 Medaka 就可以了，不用担心甲基化标签的问题！

在 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 *