0基础学习基因组三代测序组装（6）——基因组polish

racon的基本用法如下

1

racon [options ...] <sequences> <overlaps> <target sequences>

需要用到三种输入文件：sequences是指用来纠错的三代基因组测序数据（后面以原始数据称呼）；target sequences指需要校正的组装后的基因组数据（后面以组装基因组称呼）；overlaps指回比到组装基因组的原始数据文件，其中包含了所有的overlaps，其文件格式为MHAP/PAF/SAM三种之一。

因此在使用Racon之前需要使用其他比对工具将三代数据回贴到组装基因组上，在转录组分析笔记中有介绍过相关软件，我这里用minimap2进行比对，这是专门针对三代测序数据开发的比对工具，运行速度较快。

1

minimap2 -a -t 20 <target sequences> <sequences> > minimap_1.sam

-a表示结果为sam格式，<target sequences>处传入组装基因组的绝对路径，<sequences>处传入原始数据的绝对路径，比对结果的sam文件命名为minimap_1.sam

1

racon -t 50 <sequences> minimap_1.sam <target sequences> > racon_minimap_1.fasta

如上一次循环下来（大约3小时），得到的racon_minimap_1.fasta就是经过一次三代数据校正的组装基因组。

一般要用三代数据polish 2-4次，之后用二代数据继续校正4次左右，可以写脚本循环，需要注意racon因为要一次读入三代原始数据和回比的sam数据，内存需求量非常大，申请的内核数需要自己计算一下，否则会报内存溢出的错误（220Mb的基因组，100X测序深度，申请50个核才能跑动）。

脚本文件racon.sh如下：

1
2
3
4
5
6
7

#!/bin/bash

minimap2 -a -t 50 $1 $2 > minimap_1.sam
racon -t 50 $2 minimap_1.sam $1 > racon_minimap_1.fasta

minimap2 -a -t 50 racon_minimap_1.fasta $2 > minimap_2.sam
racon -t 50 $2 minimap_2.sam racon_minimap_1.fasta > racon_minimap_2.fasta

recon.slurm：

1
2
3
4
5
6
7
8
9

#!/bin/bash
#SBATCH -J recon
#SBATCH -N 1
#SBATCH -n 50
#SBATCH -t 7200

date
bash racon.sh /public/home/wlxie/NextDenovo/03_rundir/03.ctg_graph/nd.asm.fasta /public/home/wlxie/luobuma/luobuma/baima_rawdata/Third_generation_sequencing/clean_filter.fq
date

1

sbatch recon.slurm

可以从输出日志中看到，racon运行主要分两步，分别是校准overlap和生成共有序列（也就是去重），在生成共有序列（consensus sequence）之后再进行二代数据的纠错。

这一步的Racon检测出两个contig可能是嵌合体（chimeric），所谓嵌合contig，该contig的某段区域可能可以比对上不同的染色体，或者头尾部分可能分别属于不同的染色体。第一遍racon的时候没有检测到，第二遍racon才出现这个提示，我不确定这两个contig是否真的是嵌合体，最终还是需要Hi-C数据来验证。

我这里两次循环得到polish的结果文件racon_minimap_2.fasta，接下来用NextPolish软件继续用二代数据polish。

NextPolish是武汉那希望组开发的与NextDenovo配套的基因组polish软件，支持二代短读长、三代长读长和HiFi数据进行纠错。

和之前的NextDenovo操作方法类似，首先需要准备一个input文件，写入二代数据的绝对路径到sgs.fofn：

1

realpath ./1.fq ./2.fq > sgs.fofn

从doc文件夹中copy一份配置文件run.cfg，修改参数：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

[General]
job_type = slurm	# local, slurm, sge, pbs, lsf塔大学校集群选择slurm
job_prefix = nextPolish
task = best		# all, default, best, 1, 2, 5, 12, 1212…  1，2是两个不同的短读长算法模块，5是长读长算法模块，默认best
rewrite = yes
deltmp = yes	# 删除临时结果文件
rerun = 3		# 重复执行polish次数
parallel_jobs = 20	# 每个job线程数
multithread_jobs = 5	# job数
genome = /public/home/wlxie/baima_polish/racon_minimap_2.fasta
genome_size = auto
workdir = ./01_rundir
polish_options = -p {multithread_jobs}

[sgs_option] #optional
sgs_fofn = ./sgs.fofn	# 输入文件位置（一行一条）
sgs_options = -max_depth 100 -bwa	# 使用bwa进行比对

长reads和HiFi的两段配置信息删除，只留下短读长sgs_options。

这个软件的优点是速度快（100线程，4次polish，220Mb的基因组，72G的二代数据量仅仅用了8小时），而且只需要提供配置和输入文件就可以到polish结束出结果，经过4次Polish结果的迭代，最终结果如下:

Contig N50比一开始NextDenovo组装结果大，也就是组装效果更好。