192_基因组组装polish¶

一句话概述¶

基因组组装polish是利用原始reads（短读长/长读长）对初始组装进行碱基级别错误校正的过程，通过Pilon（短读长）、Medaka/Pepper（Nanopore）和DeepVariant/gcpp（PacBio）等工具将组装的碱基准确率从~99%提升到>99.99%（QV40+）。

核心知识点表格¶

知识点	说明
Polish目的	修正组装中的碱基错误（SNP和InDel）
短读长Polish	Pilon, NextPolish, Polypolish, POLCA
Nanopore Polish	Medaka（官方推荐）, 已废弃的Nanopolish
PacBio Polish	gcpp (Arrow), DeepVariant
Polish轮次	通常1-3轮，多轮边际收益递减
评估指标	QV值（BUSCO完整率、MERQURY QV）
HiFi数据	PacBio HiFi数据准确率已很高（>Q30），通常不需要Polish
推荐策略	先用同类型reads polish，再用Illumina短reads精修

步骤详解¶

第一步：理解polish的必要性¶

白话解释：组装就像拼拼图，拼完了大致形状对了，但有些边角可能不够精确。Polish就是对着原始图片（测序reads）仔细修正每一个像素（碱基），让拼图更精确。

技术细节：不同测序平台的系统性错误不同。Nanopore主要是同聚物（homopolymer）区域的InDel错误。PacBio CLR有随机错误。Illumina有低比例的替换错误但几乎无InDel错误。Polish策略是先用同平台数据校正系统性错误，再用Illumina校正残余错误。

第二步：使用Medaka校正Nanopore组装¶

# 安装Medaka
pip install medaka

# Medaka一步polish
medaka_polish \
    -i nanopore_reads.fastq.gz \
    -d draft_assembly.fa \
    -o medaka_output \
    -t 16 \
    -m r1041_e82_400bps_sup_v4.3.0  # 选择匹配的模型

# 手动分步运行（更灵活）
# Step1: 比对reads到组装
mini_align \
    -i nanopore_reads.fastq.gz \
    -r draft_assembly.fa \
    -P \
    -t 16 \
    -m

# Step2: 运行Medaka consensus
medaka consensus \
    medaka_output/calls_to_draft.bam \
    medaka_output/consensus_probs.hdf \
    --model r1041_e82_400bps_sup_v4.3.0 \
    --threads 16 \
    --batch_size 100

# Step3: 生成polish后的序列
medaka stitch \
    medaka_output/consensus_probs.hdf \
    draft_assembly.fa \
    medaka_output/consensus.fasta

第三步：使用Pilon进行短读长polish¶

# Step1: 将Illumina reads比对到组装
bwa index draft_assembly.fa
bwa mem -t 16 draft_assembly.fa illumina_R1.fq.gz illumina_R2.fq.gz | \
    samtools sort -@ 8 -o aligned.bam
samtools index aligned.bam

# Step2: 运行Pilon
java -Xmx64G -jar pilon.jar \
    --genome draft_assembly.fa \
    --frags aligned.bam \
    --output pilon_polished \
    --outdir pilon_output \
    --changes \
    --vcf \
    --threads 16

# 多轮Polish
for round in 1 2 3; do
    echo "=== Round $round ==="
    ref="draft_assembly.fa"
    [[ $round -gt 1 ]] && ref="pilon_round$((round-1)).fa"

    bwa index $ref
    bwa mem -t 16 $ref R1.fq.gz R2.fq.gz | samtools sort -@ 8 -o round${round}.bam
    samtools index round${round}.bam

    java -Xmx64G -jar pilon.jar \
        --genome $ref \
        --frags round${round}.bam \
        --output pilon_round${round} \
        --changes

    # 检查改变数量
    wc -l pilon_round${round}.changes
done

第四步：使用NextPolish（推荐，速度更快）¶

# 安装NextPolish
conda install -c bioconda nextpolish

# 准备配置文件
cat << 'EOF' > run.cfg
[General]
job_type = local
job_prefix = nextPolish
task = best
rewrite = yes
rerun = 3
parallel_jobs = 4
multithread_jobs = 8
genome = draft_assembly.fa
genome_size = auto
workdir = ./nextpolish_output
polish_options = -p {multithread_jobs}

[sgs_option]
sgs_fofn = sgs.fofn
sgs_options = -max_depth 100 -bwa
EOF

# 准备reads列表
echo "illumina_R1.fq.gz" > sgs.fofn
echo "illumina_R2.fq.gz" >> sgs.fofn

# 运行
nextPolish run.cfg

第五步：评估polish效果¶

# 1. MERQURY评估QV值
meryl count k=21 output reads.meryl illumina_R1.fq.gz illumina_R2.fq.gz
merqury.sh reads.meryl polished_assembly.fa merqury_output
# QV值越高越好，>40表示每10000个碱基不到1个错误

# 2. BUSCO完整率
busco -i polished_assembly.fa -l eukaryota_odb10 -m genome -c 16 -o busco_output

# 3. 前后对比
echo "=== Before Polish ==="
busco -i draft_assembly.fa -l eukaryota_odb10 -m genome -c 16 -o busco_before
echo "=== After Polish ==="
busco -i polished_assembly.fa -l eukaryota_odb10 -m genome -c 16 -o busco_after

# 4. 统计polish改变的数量
grep -c "^" pilon_round1.changes
grep -c "^" pilon_round2.changes
grep -c "^" pilon_round3.changes
# 应该逐轮递减，趋近于零

实战命令速查¶

# Medaka一步polish
medaka_polish -i reads.fq.gz -d draft.fa -o output/ -t 16

# Pilon polish
bwa mem -t 16 draft.fa R1.fq.gz R2.fq.gz | samtools sort -o aligned.bam
java -Xmx64G -jar pilon.jar --genome draft.fa --frags aligned.bam --output polished

# Polypolish（多轮短读长polish）
bwa mem -t 16 -a draft.fa R1.fq.gz > R1.sam
bwa mem -t 16 -a draft.fa R2.fq.gz > R2.sam
polypolish filter --in1 R1.sam --in2 R2.sam --out1 R1_filtered.sam --out2 R2_filtered.sam
polypolish polish draft.fa R1_filtered.sam R2_filtered.sam > polished.fa

# MERQURY评估
meryl count k=21 output reads.meryl reads.fq.gz
merqury.sh reads.meryl assembly.fa output_prefix

面试常问点¶

Q1: 为什么PacBio HiFi组装通常不需要polish？ A: HiFi reads的碱基准确率已达99.9%（Q30+），使用CCS（Circular Consensus Sequencing）对同一分子多次测序并取共识。hifiasm等工具直接组装HiFi reads产生的consensus已经足够准确。额外polish的改善微乎其微，甚至可能引入错误。

Q2: Polish应该做几轮？ A: 通常2-3轮足够。每轮校正的错误数应该递减。当一轮的改变数量趋近于零（<10个改变）时应停止。过多轮次可能引入新错误（over-polishing）。建议用MERQURY QV值监控每轮效果。

Q3: Medaka和Nanopolish有什么区别？ A: Nanopolish使用HMM模型直接分析原始电信号，准确但极慢且需要fast5文件。Medaka使用神经网络模型，只需要basecalled FASTQ即可，速度快很多。2024年后推荐使用Medaka，Nanopolish已逐步废弃。

Q4: 短读长和长读长polish的推荐顺序？ A: 先用同类型长读长polish（如Medaka校正Nanopore系统性错误），再用Illumina短读长polish精修残余错误。这是因为短读长在重复区域比对不准确，先校正大部分错误后再用短读长处理效果更好。

Q5: 如何判断组装质量是否足够好？ A: (1)MERQURY QV>40（每万碱基<1个错误）；(2)BUSCO完整率>95%；(3)基因预测的完整性提升；(4)连续两轮polish改变数<100。

易错点¶

使用错误的Medaka模型：必须匹配实际使用的测序化学和basecaller版本
Pilon内存不足：大基因组需要64-128GB内存，使用-Xmx参数设置
Over-polishing：过多轮次可能在杂合区域引入错误
忘记索引BAM文件：Pilon需要排序并索引的BAM文件
HiFi数据不必要的polish：可能降低质量而非提升

补充知识¶

Polish工具对比¶

工具	数据类型	速度	准确率提升	推荐度
Medaka	Nanopore	快	高	强推荐
Pilon	Illumina	中	中	推荐
NextPolish	Illumina/ONT	快	高	推荐
Polypolish	Illumina	快	中	推荐
POLCA	Illumina	快	中	可选
gcpp(Arrow)	PacBio CLR	慢	高	推荐(CLR)