0基础学习基因组三代测序组装（5）——三代数据组装

前段时间比较忙，现在继续整理基因组测序组装系列的学习笔记。第四篇笔记写的二代测序基因组组装，主要是演示二代测序数据组装的主流工具SOAPdenovo 2.0是如何应用的。我这里有了二代和三代的测序数据，后续组装还是以三代数据为主，这里就继续记录下几款三代测序数据组装的主流工具和用法。

现在主流的三代测序公司是Pacbio和Nanopore，两家测序公司测序原理不同，产生的数据类型也有区别。

1. 主流三代测序平台

1.1 Pacbio测序平台

Pacbio测序平台是单分子实时测序（single molecule real time sequencing，SMRT），原理是当DNA与聚合酶形成的复合物被ZMW（零模波导孔）捕获后，4种不同荧光标记的dNTP随机进入检测区域并与聚合酶结合，与模板匹配的碱基生成化学键并激发荧光，生成化学键激发的荧光存在的时间远远长于其他碱基被激发荧光的时间，从而实现单碱基的实时检测。

现在Pacbio有两种测序模式，一种是CLR测序模式（超长测序模式），产生数据基于单循环测序结果；一种是HiFi测序模式，也就是高保真测序模式，产生数据基于滚环一致序列（Circular Consensus Sequencing ，CCS）。具体原理就不说了，两种测序模式中HiFi数据相对来说准确率会高一些，所以不同软件对两种测序模式的数据也会有不同的处理。

1.2 Nanopore测序平台

Nanopore测序平台前面第一篇博客介绍过了，可以点击这里。需要了解ONT测序平台测序产生的原始数据是电信号，经过basecalling之后才可以转化成我们要的测序数据。

2. 三代基因组测序组装软件

2.1 NextDenovo

NextDenovo是武汉希望组开发的集校正、比对和组装一体的，基于字符串图（string graph-based）的三代测序基因组组装软件。它的实现过程和另一款经典的三代基因组组装软件Canu类似，经过长读长数据的纠错校正后再进行组装。

官网上介绍原来可以对CLR、HiFi和ONT数据都可以组装，HiFi数据可以跳过数据的自我纠错过程，如今HiFi数据被划掉了，也许已经不再适用，但是对Pabio的CLR和Nanopore的ONT测序数据仍有较好的组装效果，其介绍是组装的准确率有98%-99.8%。

NextDenovo主要有两个核心模块 NextCorrect和NextGraph。NextCorrect用于原始数据纠错，NextGraph用于纠错后数据的组装。据作者介绍，与其他工具相比，NextDenovo在装配一致性和单碱基装配精度方面表现出较高的水平。我个人用起来是觉得这个软件运行时间相对Canu较短，需要的算力资源较小，可以很快地组装出结果（后面可以进行比较）。

安装并测试通过后，我们就可以开始使用这个工具了。

首先准备input文件，将前面质控后的三代测序数据的绝对路径写在input.fofn文件里：

1	/public/home/wlxie/luobuma/luobuma/baima_rawdata/Third_generation_sequencing/clean_filter.fq

接下来也是最重要的，修改配置文件：

[General]
job_type = slurm # local, slurm, sge, pbs, lsf塔大学校集群选择slurm
job_prefix = nextDenovo
task = all # all, correct, assemble选择只进行correct还是只进行assemble，或者两者都进行，基因组小的话可以直接all
rewrite = yes # yes/no覆写
deltmp = yes 
parallel_jobs = 20 # number of tasks used to run in parallel线程数，咱学校的集群20勉强够
input_type = raw # raw, corrected输入的数据情况
read_type = ont # clr, ont, hifi数据类型
input_fofn = input.fofn		# 输入数据的位置信息
workdir = 03_rundir		# 输出的文件夹名字

[correct_option]
read_cutoff = 1k	# 进行correct的时候截取的最小read
genome_size = 230m # estimated genome size预估的基因组大小
sort_options = -m 20g -t 15
minimap2_options_raw = -t 8
pa_correction = 3 # number of corrected tasks used to run in parallel, each corrected task requires ~TOTAL_INPUT_BASES/4 bytes of memory usage.
correction_options = -p 15 -dbuf	# 非常重要！-dbuf让每一步作业释放内存，防止节点卡死！

[assemble_option]
minimap2_options_cns = -t 8 
nextgraph_options = -a 1

以上是我的配置文件信息，中文标注的地方都很重要，根据情况修改。其他参数可以用默认。其中预估的基因大小也是很有必要的，前面在做基因组Survey的时候预测过，这里就直接写预测的基因组大小。

其他参数的设定和使用可以参考这篇博客使用NextDenovo组装Nanopore数据，以及官方的参数说明手册NextDenovo Parameter Reference — NextDenovo latest documentation

需要强调一点，correction_options = -p 15 -dbuf这项参数是我在华农的集群平台手册上看到的，之前确实一直会卡死在某一步直到24h后台强杀这个进程，目前未知原因，加上之后运行正常。以我的数据来看，一个200多Mb的植物基因组，测序深度100X左右，一次组装运行结束需要12小时左右，已经非常快了。

最后是运行程序，我写了一个run.slurm文件：

#!/bin/bash
#SBATCH -n 40

./nextDenovo run.cfg

1	sbatch run.slurm

基因组比较大的话，建议分步运行，先correct，再assemble。

因为是在集群中运行，所有输出都会在slurm-xxxx.out的文件夹中显示，打开以后可以看到每个时间节点完成了什么任务，当有任务卡住几个小时都没动的时候，就要检查是否是配置文件是否正确。

最后组装的基因组在03.ctg_graph目录下，文件名称nd.asm.fasta。最底下输出了组装结果概况，contig N50为9Mb，总共组装出225Mb的基因组序列，contig总数为59，组装结果还算不错。

2.2 Canu

Canu是三代测序数据组装的经典工具，也是主要用于Pacbio和Nanopore公司的测序结果组装。

这个软件在安装过程中有点曲折，从官网下数据包，最后一步编译的过程会报错。目前我暂时没办法解决，但是可以用conda安装（虽然官网不建议这么做）。

1	conda install -c conda-forge -c bioconda -c defaults canu

如果报错动态库出问题，可以参考第三篇博客中的方法，寻找根目录下的动态库中是否有对应的版本文件，如果有，直接修改软链接到对应的动态库下。

Canu运行分为三个步骤：纠错（Correct）、修剪（Trim）和组装（Assemble）。可以每一个步骤分开跑，比如纠错修剪之后的数据可以放到别的软件中组装，或者用别的软件纠错之后作为输入到Canu中组装。考虑到塔大集群24小时自动杀程序，保险起见还是三个步骤分开跑比较安全。

下面是官方的帮助文档，写的非常详细：

usage:   canu [-version] [-citation] \
              [-haplotype | -correct | -trim | -assemble | -trim-assemble] \
              [-s <assembly-specifications-file>] \
               -p <assembly-prefix> \	# 输出文件前缀
               -d <assembly-directory> \	# 输出目录
               genomeSize=<number>[g|m|k] \		# 预测基因组大小
              [other-options] \
              [-haplotype{NAME} illumina.fastq.gz] \
              [-corrected] \
              [-trimmed] \
              [-pacbio |
               -nanopore |
               -pacbio-hifi] file1 file2 ...

example: canu -d run1 -p godzilla genomeSize=1g -nanopore-raw reads/*.fasta.gz 

  To restrict canu to only a specific stage, use:	# 描述canu要执行的主程序
    -haplotype     - generate haplotype-specific reads
    -correct       - generate corrected reads
    -trim          - generate trimmed reads
    -assemble      - generate an assembly
    -trim-assemble - generate trimmed reads and then assemble them

  Reads can be either FASTA or FASTQ format, uncompressed, or compressed with gz, bz2 or xz.

  Reads are specified by the technology they were generated with, and any processing performed.

  [processing]		# 描述reads状态
    -corrected
    -trimmed

  [technology]		# 描述测序平台（数据类型）
    -pacbio      <files>
    -nanopore    <files>
    -pacbio-hifi <files>
    
  Some common options:
    useGrid=string
      - Run under grid control (true), locally (false), or set up for grid control
        but don't submit any jobs (remote)
    rawErrorRate=fraction-error		# 降低这个参数会提高第一步的速度
      - The allowed difference in an overlap between two raw uncorrected reads.  For lower
        quality reads, use a higher number.  The defaults are 0.300 for PacBio reads and
        0.500 for Nanopore reads.  
    correctedErrorRate=fraction-error	# 降低这个参数可以提高组装效率
      - The allowed difference in an overlap between two corrected reads.  Assemblies of
        low coverage or data with biological differences will benefit from a slight increase
        in this.  Defaults are 0.045 for PacBio reads and 0.144 for Nanopore reads.
    gridOptions=string
      - Pass string to the command used to submit jobs to the grid.  Can be used to set
        maximum run time limits.  Should NOT be used to set memory limits; Canu will do
        that for you.
    minReadLength=number
      - Ignore reads shorter than 'number' bases long.  Default: 1000.
    minOverlapLength=number
      - Ignore read-to-read overlaps shorter than 'number' bases long.  Default: 500.
  A full list of options can be printed with '-options'.  All options can be supplied in
  an optional sepc file with the -s option.	# 可以用-s来提供自己修改的参数文件

Complete documentation at http://canu.readthedocs.org/en/latest/

也可以点击这里，进官方手册看原文。下面用我自己的三代数据跑一个案例。

2.2.1 纠错（Correct）

因为三代测序数据错误率较高，纠错的步骤是通过序列之间的一致性比较获得高可信的碱基。

创建一个canu的空文件夹，写入以下内容到correct.slurm:

#!/bin/bash
#SBATCH -N 1
#SBATCH -n 40

canu -correct -p baima -d baima_nanopore genomeSize=230m -nanopore-raw /public/home/wlxie/luobuma/luobuma/baima_rawdata/Third_generation_sequencing/pass.fq

1	sbatch correct.slurm

主要就是注意下参数，-p是输出文件的前缀，-d是输出文件的目录名，需要声明这个数据是什么平台测的，以及数据是什么状态。虽然我这里只申请了40个核，但是canu会自动提交作业直到你能申请的核数上限……在塔大集群我的权限是200个核，通过scontrol show job 可以查到，我这边一次性提交了136个作业，排队100多个任务，占用192个核…..

纠错、修整和组装每一个步骤会依次进行以下各个阶段，需要的内存和核数挺高的，所以推荐在集群中运行。

如果运行时间很长，建议到设置的输出文件夹目录下查看canu.out文件，详细记载了正在执行哪一步以及花了多少时间。如果不确定程序是否卡死，直接通过scontrol show job命令查看状态为RUNNING的作业，进入作业的输出目录，如果文件夹中的内容一直到最近的时间点都有更新，则可以放心地继续运行。

仅仅这一个步骤花了30个小时。并且这一步会将100X的测序数据量降到40X（默认，可以调整，见官方readSampleingCoverage参数介绍）。最后生成文件baima.correctedReads.fasta.gz，我为了方便复制到了前一个文件夹，修改文件权限为0755。

2.2.2 修整（trim）

修整是在上一步纠错的基础上，再对reads进行修剪，删去可疑区域。

这一步经历的步骤与上一步是一样的，虽然上一步纠错已经将数据量降了一大半，对于我的基因组，这一步依然要跑32小时.

同样在canu文件夹写入如下内容到trim.slurm：

#!/bin/bash
#SBATCH -n 50
#SBATCH -t 7200

canu -trim -p baima -d baima_trim genomeSize=230m -corrected -nanopore /public/home/wlxie/canu/baima.correctedReads.fasta.gz

1	sbatch trim.slurm

注意修改参数以及reads所处的状态。

在canu.out输出文件中，可以找到trim步骤用了什么参数，处理了哪些类型的reads：

生成的结果文件在baima_trim文件夹中，名称为baima.trimmedReads.fasta.gz，同样修改文件权限，移到前一个文件夹方便操作。

2.2.3 组装（Assemble）

经过前两步的数据纠错和修整，这一步才是正式组装基因组。

在canu文件夹写入如下内容到assemble.slurm：

#!/bin/bash
#SBATCH -n 50
#SBATCH -t 10800

canu -assemble -p baima -d baima_assemble genomeSize=230m correctedErrorRate=0.144 -trimmed -corrected -nanopore /public/home/wlxie/canu/baima.trimmedReads.fasta.gz

1	sbatch assemble.slurm

在官方的介绍中，correctedErrorRate这个参数可以根据前面纠错和修整的reads质量做修改的，默认是Pacbio数据0.045，Nanopore数据0.144，降低这个参数值可以加快组装的效率，但是存在遗漏overlap和组装片段断裂的风险。低于30X测序深度以下的数据，官方建议可以略微提高这个值，对于60X测序深度以上的数据可以略微降低这个值，每次改变1%左右比较合适。

我这里就用默认参数了，组装时间可能会比较长，就将作业的时间调整为7天。

实际运行时间为30小时，最终组装结果如下：

这个工具组装出的结果比预期大很多，前面基因组survey预测的基因组大小为230Mbp，实际通过canu组装出有276Mbp，且contig数明显比NextDenovo多，导致contig N50指标低。我觉得可能是correctedErrorRate这个值比较高，可以适当调低一些，过于严格的纠错标准可能导致组装的contig比较碎。

因为前面的NextDenovo组装的效果已经比较理想，因此这一步我也就不再细调参数了。以NextDenovo组装出的基因组继续后面的分析。以目前我的植物基因组来看，用NextDenovo组装三代基因组的效率和质量都比Canu高。