组装结果纠错
本帖最后由 生信喵 于 2022-6-20 08:27 编辑1 组装结果优化原理
为什么需要对组装结果进行矫正（polishing)?
由于三代 nanopore 测序质量比较低，原始数据中存在大量测序错误，即使拼接前进行了纠错，组装结果中仍会存在错误，用长读长或短读长的数据对组装结果进行矫正可以，提高准确率，减少 Miscalls，Indels，改善由错装（mis-assemblies)导致的低比对区域。因此，序列拼接完需要对拼接结果进行优化，根据文献报道，经过 polish 之后，拼接结果与真实基因组（其他测序数据拼接结果）的一致性可以达到 99.99%以上。即使组装工具带有纠错更能，仍建议再次进行一轮或多轮的矫正。
为什么需要对组装结果进行多轮矫正（polishing)?
这是因为 nanopore 数据主要的错误来自于插入与缺失，每次将测序数据与拼接基因组比对能够发现一些错误。下一轮数据与纠错后的序列重新比对，可以发现新的错误，这样经过多轮之后，就可以逐渐减少错误。
常用纠错工具：medaka，pilon，racon，nanopolish，nextpolish 等，可以利用三代测序进行纠错，也可以加入二代数据进行纠错。
2 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比较纠错前后差别。
3 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剩余代码流程和前述一致，修改下文件名即可。
4 racon 组装结果纠错
Racon 是一个基于 minimap 和 miniasm 的，构建一致性序列（consensus）的一款软件，也可以用于纠错。既可以用于三代数据也可以用于二代数据的纠错。输入数据需要三个，首先是 contig，然后是测序的 reads，以及前面二者比对的结果，这个比对结果可以是 MHAP，PAF，SAM 等三种格式当中的一种即可。数据结果为纠错后的 contig 序列。一般 racon 纠错也可以进行多轮，一般3轮纠错。
案例：利用 racon 进行纠错
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
5如何对一个物种做全基因组鉴定或者对植物做全基因组测序？
第一步背景调研：查资料该物种是否测过序，若测过，技术上有无突破；
第二步基因组大小：查资料、近源参考序列等；（2G）
第三步测序方案：至少要测（2x30倍=60G或者200倍=400G）；3代测序+2代纠错，nanopore+illumina；
第四步质控过滤；
第五步flye，canu，wtdbg2，smartdenovo等拼接基因组；
第六步调整选项参数，保证拼接结果最优；
第七步纠错；
第八步HiC数据（植物常用，定位染色体）。
就可以得到一个基因组了。
页:
[1]