组装结果纠错

生信喵

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
软件安装运行

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

复制代码

1. conda activate medaka
   
2. #运行软件
   
3. READ=/share/home/xiehs/05.assembly/data/nanopore.fastq.gz
   
4. 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 并行计算

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

复制代码

比较纠错前后差别。

3 pilon 组装结果纠错
      pilon 是由 broadinstitute 研究所开发的纠错工具，输入原始拼接结果以及原始测序数据比对到拼接结果的 bam 文件即可。pilon 通过比对后的 bam 文件，可以识别拼接中非一致性的序列，包括单碱基的不同，小的 indel，大的 indel，后者空位 gap，以及错误拼接的区域。
      输入的 bam 可以来自于二代测序数据的比对，也可以来自于三代测序数据比对得到的 bam，bam 文件需要排序并建立索引。
      
      pilon 纠错流程图

1. echo "java -jar ~/biosoft/pilon-1.24.jar" >pilon
   
2. chmod u+x pilon

复制代码

1. ./pilon
   
2. mv pilon ~/bin/
   
3. which pilon

复制代码

1. #对拼接结果建立索引
   
2. mkdir pilon
   
3. ln -s ../medaka/medaka_result/consensus.fasta medaka.fasta
   
4. bwa index medaka.fasta
   
5. 
   
6. #illumina比对排序建索引
   
7. READ1=/ifs1/TestDatas/nanopore7/data/MGH78578/illumina.sra_1.fastq.gz
   
8. READ2=/ifs1/TestDatas/nanopore7/data/MGH78578/illumina.sra_2.fastq.gz
   
9. bwa mem -t 4 -R '@RG\tID:foo\tSM:bar:\tPL:ILLUMINA' medaka.fasta $READ1 $READ2 >illumina.sam
   
10. samtools sort -@ 4 -O bam -o illumina.sorted.bam illumina.sam
    
11. samtools index illumina.sorted.bam
    
12. #利用pilon进行纠错
    
13. 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
    
14. ​

复制代码

还可以继续第二次纠错。

1. bwa index pilon.fasta

复制代码

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

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

1. mkdir racon
   
2. #连接原始拼接结果
   
3. DRAFT=../pilon/pilon_result/pilon.fasta
   
4. READ=/ifs1/TestDatas/nanopore7/data/MGH78578/clean.filtlong.fq.gz
   
5. 
   
6. #minimap2比对
   
7. minimap2 -t 4 ${DRAFT} ${READ} > round_1.paf
   
8. #racon进行纠错
   
9. racon -t 4 ${READ} round_1.paf ${DRAFT} > racon_round1.fasta
   
10. 
    
11. #第二轮纠错   
    
12. minimap2 -t 4 racon_round1.fasta ${READ} > round_2.paf
    
13. racon -t 4 ${READ} round_2.paf racon_round1.fasta> racon_round2.fasta
    
14. 
    
15. #第三轮纠错
    
16. minimap2 -t 4 racon_round2.fasta ${READ}  > round_3.paf
    
17. racon -t 4 ${READ}  round_3.paf racon_round2.fasta> racon_round3.fasta
    
18. 
    
19. #将最终结果修改为样品名
    
20. mv racon_round3.fasta MGH78578.fasta

复制代码

5如何对一个物种做全基因组鉴定或者对植物做全基因组测序？
      第一步背景调研：查资料该物种是否测过序，若测过，技术上有无突破；
      第二步基因组大小：查资料、近源参考序列等；（2G）
      第三步测序方案：至少要测（2x30倍=60G或者200倍=400G）；3代测序+2代纠错，nanopore+illumina；
      第四步质控过滤；
      第五步flye，canu，wtdbg2，smartdenovo等拼接基因组；
      第六步调整选项参数，保证拼接结果最优；
      第七步纠错；
      第八步HiC数据（植物常用，定位染色体）。
      就可以得到一个基因组了。