组装结果纠错

一、组装结果优化原理

1.1为什么需要对组装结果进行矫正（polishing)?

由于三代 nanopore 测序质量比较低，原始数据中存在大量测序错误，即使拼接前进行了纠错，组装结果中仍会存在错误，用长读长或短读长的数据对组装结果进行矫正可以，提高准确率，减少 Miscalls，Indels，改善由错装（mis-assemblies)导致的低比对区域。因此，序列拼接完需要对拼接结果进行优化，根据文献报道，经过 polish 之后，拼接结果与真实基因组（其他测序数据拼接结果）的一致性可以达到 99.99%以上。即使组装工具带有纠错功能，仍建议再次进行一轮或多轮的矫正。

1.2为什么需要对组装结果进行多轮矫正（polishing)?

这是因为 nanopore 数据主要的错误来自于插入与缺失，每次将测序数据与拼接基因组比对能够发现一些错误。下一轮数据与纠错后的序列重新比对，可以发现新的错误，这样经过多轮之后，就可以逐渐减少错误。

常用纠错工具：medaka，pilon，racon，nanopolish，nextpolish 等，可以利用三代测序进行纠错，也可以加入二代数据进行纠错。

二、 medaka 组装结果纠错

Medaka 是一个基于叠加序列的一致性序列修正工具，高度推荐使用以获得最佳的一致性准确性。Medaka 现可以用于变体识别（variant calling）。使用纳米孔 R9.4.1 版芯片和最佳的工具，现在你可以进行 SNPs 识别，获得 99%准确率。例如，使用当前的 R9.4.1 版纳米孔，利用 Flip-flop 碱基识别器和 Medaka, 测序金黄色葡萄球菌（S.aureus）基因组，准确性达到 Q44，其它的小型基因组准确率约 Q40。

软件特色：

✓ 由 Oxford Nanopore 开发的开源软件

✓ 仅需使用.fasta 或.fastq 数据

✓ 速度比 Nanopolish 快 50 倍，支持 CPU 和 GPU

✓ 通常在 Pomoxis 组装后使用

✓ 用 FASTQ 文档和组装结果作为输入文件

✓ 50X5Mbase 基因组用时 20 分钟

✓ 在 Racon 基础上，进一步提升了数据准确率

✓ 可针对不同数据进行个性化纠错方法训练

✓ 兼容 Linux 和 MacOS 系统

官网：https://github.com/nanoporetech/medaka

软件安装运行

#创建文件夹
mkdir 41.polish;cd 41.polish
cp /share/home/xiehs/05.assembly/40.nanopore/flye/assembly.fasta .

conda activate medaka
#运行软件
READ=/share/home/xiehs/05.assembly/data/nanopore.fastq.gz
medaka_consensus -i $READ -d assembly.fasta -o medaka_result -m r941_min_high_g360 -v medaka.vcf -t 24 >medaka.log

-i 输入测序 reads

-d 需要纠错的拼接结果

-o 输出结果文件

-m 芯片类型,需要选好。

-t 并行计算

seqkit stat assembly.fasta ./medaka_result/consensus.fasta

比较纠错前后差别。

三、 pilon 组装结果纠错

pilon 是由 broadinstitute 研究所开发的纠错工具，输入原始拼接结果以及原始测序数据比对到拼接结果的 bam 文件即可。pilon 通过比对后的 bam 文件，可以识别拼接中非一致性的序列，包括单碱基的不同，小的 indel，大的 indel，后者空位 gap，以及错误拼接的区域。

输入的 bam 可以来自于二代测序数据的比对，也可以来自于三代测序数据比对得到的 bam，bam 文件需要排序并建立索引。

pilon 纠错流程图

echo "java -jar ~/biosoft/pilon-1.24.jar" >pilon
chmod u+x pilon
./pilon
mv pilon ~/bin/
which pilon
#对拼接结果建立索引
mkdir pilon
ln -s ../medaka/medaka_result/consensus.fasta medaka.fasta
bwa index medaka.fasta

#illumina比对排序建索引
READ1=/ifs1/TestDatas/nanopore7/data/MGH78578/illumina.sra_1.fastq.gz
READ2=/ifs1/TestDatas/nanopore7/data/MGH78578/illumina.sra_2.fastq.gz
bwa mem -t 4 -R '@RG\tID:foo\tSM:bar:\tPL:ILLUMINA' medaka.fasta $READ1 $READ2 >illumina.sam
samtools sort -@ 4 -O bam -o illumina.sorted.bam illumina.sam
samtools index illumina.sorted.bam
#利用pilon进行纠错
java -Xmx32G -jar /ifs1/Software/biosoft/pilon/pilon-1.23.jar --genome medaka.fasta --fix all --changes --frags illumina.sorted.bam --output pilon --outdir pilon_result --threads 24 --vcf 2> pilon.log

还可以继续第二次纠错。

bwa index pilon.fasta

剩余代码流程和前述一致，修改下文件名即可。

四、racon 组装结果纠错

Racon 是一个基于 minimap 和 miniasm 的，构建一致性序列（consensus）的一款软件，也可以用于纠错。既可以用于三代数据也可以用于二代数据的纠错。输入数据需要三个，首先是 contig，然后是测序的 reads，以及前面二者比对的结果，这个比对结果可以是 MHAP，PAF，SAM 等三种格式当中的一种即可。数据结果为纠错后的 contig 序列。一般 racon 纠错也可以进行多轮，一般3轮纠错。

mkdir racon
#连接原始拼接结果
DRAFT=../pilon/pilon_result/pilon.fasta
READ=/ifs1/TestDatas/nanopore7/data/MGH78578/clean.filtlong.fq.gz

#minimap2比对
minimap2 -t 4 ${DRAFT} ${READ} > round_1.paf
#racon进行纠错
racon -t 4 ${READ} round_1.paf ${DRAFT} > racon_round1.fasta

#第二轮纠错    
minimap2 -t 4 racon_round1.fasta ${READ} > round_2.paf
racon -t 4 ${READ} round_2.paf racon_round1.fasta> racon_round2.fasta

#第三轮纠错
minimap2 -t 4 racon_round2.fasta ${READ}  > round_3.paf
racon -t 4 ${READ}  round_3.paf racon_round2.fasta> racon_round3.fasta

#将最终结果修改为样品名
mv racon_round3.fasta MGH78578.fasta

五、如何对一个物种做全基因组鉴定或者对植物做全基因组测序？

第一步背景调研：查资料该物种是否测过序，若测过，技术上有无突破；

第二步基因组大小：查资料、近源参考序列等；（2G）

第三步测序方案：至少要测（2x30倍=60G或者200倍=400G）；3代测序+2代纠错，nanopore+illumina；

第四步质控过滤；

第五步flye，canu，wtdbg2，smartdenovo等拼接基因组；

第六步调整选项参数，保证拼接结果最优；

第七步纠错；

第八步HiC数据（植物常用，定位染色体）。

就可以得到一个基因组了。

写在最后：有时间我们会努力更新的。大家互动交流可以前去论坛，地址在下面，复制去浏览器即可访问，弥补下公众号没有留言功能的缺憾。原地址暂未启用（bioinfoer.com）。

sx.voiceclouds.cn

有些板块也可以预设为大家日常趣事的分享等，欢迎大家来提建议。