三代测序数据基因组组装流程（一）

说明

基因组组装一直是生物信息学中的一个重点以及难点，想要得到高质量的组装结果非常的困难，不仅需要花费大量的钱进行测序，还得消耗很多的精力进行一遍又一遍的组装和polish。在这里会根据我们课题组在2023年的一篇文章所提出的方法，给出一个完整的组装流程，以供借鉴，希望能够帮到大家。
文章的链接如下（https://www.nature.com/articles/s41467-022-35670-y），有兴趣可以自行阅读，这篇文章是我师兄在德国的成果，这老哥是个埃及人，所以代码写的非常非常捞，抽象的一比，谁看谁懵逼，绝对不建议直接使用他提供的源码（https://github.com/ganlab/GALA），当然你如果想挑战一下自己也是没问题。这里我得感谢另一位JohnUrban大哥不辞辛苦的把源码改了一遍，你可以直接把原本复杂的步骤自动化的完整运行，所以建议使用这个（https://github.com/JohnUrban/GALA）。
这个方法主要是针对三代测序数据Nanopore和Pacbio HiFi进行组装，别的我没试过，不知道效果怎么样，目前我已经尝试过水稻、桃金娘、葡萄、人类、猪的基因组，都是有用的，虽然不一定能够达到目前学术界所追求的所谓T2T的组装结果，但是能够有效提升组装质量。
整个方法的核心思想就是希望将原始的测序reads按照染色体进行分离，各条染色体独立进行组装，一方面能够减少组装时的数据复杂度，提升组装效率，另一方面可以避免来自不同染色体数据之间的相互干扰，提高组装的准确度。这个思路我觉得没毛病，但是用起来就很难评，不过也确实提供基因组组装一个新的解决思路。原本是希望能够直接从测序reads入手，直接根据reads之间的关系将其按染色体分离，但是埃及老哥没搞出来，所以目前只能退一步，先用已有的组装软件对测序数据进行预组装，将预组装得到的contigs作为输入数据进行后续的操作。至于直接分离reads我另一个师兄目前也已经实现，但文章尚未发表，有兴趣可以关注后续的更新。

usage: gala -h  [options] <draft_names & paths> <fa/fq> <reads> <platform>

GALA Gap-free Long-read Assembler

positional arguments:
  draft_names           Draft names and paths [required]
  input_file            input type (fq/fa) [required]
  reads                 raw/corrected reads [required]
  sequencing_platform   -pacbio-raw -pacbio-corrected -nanopore-raw -nanopore-
                        corrected [required]

第一步 初步组装

想要运行GALA，我们需要两类输入数据，第一是原始的测序数据（fq/fa），第二则是利用已有组装工具进行初步组装后得到的contigs（fa）。因此第一步我们需要借助已有的组装软件对测序数据进行初步组装，我们还需要根据测序数据的实际情况选择适合的组装软件，经过实际测试，我发现以下组装软件比较好用，Canu（https://github.com/marbl/canu）、Flye（https://github.com/fenderglass/Flye）、Miniasm（https://github.com/lh3/miniasm）、Necat（https://github.com/xiaochuanle/NECAT）、NextDenovo（https://github.com/Nextomics/NextDenovo）、SMARTdenovo（https://github.com/ruanjue/smartdenovo）、Wtdbg2（https://github.com/ruanjue/wtdbg2）、Hifiasm（https://github.com/chhylp123/hifiasm），我会以水稻和桃金娘的测序数据为例，展示这些组装软件的使用方法和组装结果，至于他们的详细信息可以点击链接进行查看。
这里水稻采用的是公开的数据集，水稻品种为Dom Sufid，数据类型为Nanopore，原始数据共有42.7GB，覆盖深度约为56x，共有12条染色体；桃金娘则是HiFi数据，共67GB，深度约为74x，共有11条染色体。
假设水稻测序数据文件为rice.fasytq，桃金娘测序数据文件为rt.fastq详细的命令如下。

##Canu
#rice
canu -p rice -d canu genomeSize=380 -nanopore ./rice.fastq 
#rt
canu -p rt -d canu genomeSize=480 -pacbio-hifi ./rt.fastq 

Canu在组装之前会进行纠错，生成纠错后的文件rice.trimmedReads.fasta.gz，这个纠错后的文件也可以作为其他组装软件的输入数据，提高组装质量，canu.contigs.fasta是组装完成的contigs序列。

##Flye
#rice
flye --nano-raw ./rice.fastq --genome-size 380m -o flye_raw --threads 10
flye --nano-corr ./canu/rice.trimmedReads.fasta.gz --genome-size 380m -o flye_corr --threads 10
#rt
flye --pacbio-hifi ./rt.fastq --genome-size 480m -o flye --threads 10

Flye对于ONT测序数据有着非常好的表现，在使用中也可以增加参数--iteration n让他进行校准，n表示进行n次循环，能够有效提升组装的连续性。同时我们还用Flye组装经过Canu纠错的数据，因为ONT数据的错误率较高，会导致组装结果准确度下降，若是先进行纠错然后再组装，会提高组装结果的质量。

##Miniasm
#rice
name=miniasm
mkdir -p $name
cd $name
reads=./canu/rice.trimmedReads.fasta.gz
minimap2 -x ava-ont -t10 $reads $reads | gzip -1 > name.paf.gz
miniasm -f $reads $name.paf.gz > $name.gfa
awk ‘/^S/{print “>” $2 “\n” $3}’ $name.gfa > $name.fasta
name1=miniasm_first_round
genome=$name.fasta
minimap2 -x map-ont $genome $reads > name1.paf
racon $reads $name1.paf $genome > $name1.fasta

Miniasm的核心算法是利用reads之间的重叠关系不断地延长，因此在组装超长的ONT数据时表现较好，因为HiFi测序reads长度较短，并不适用。在Miniasm运行完成后，使用Racon进行纠错能够有效改善组装质量，在结果较差时也可以选择多次Racon进行纠错。这里Miniasm组装的是Canu纠错后生成的FASTA文件，若是组装原始数据，只需进行替换即可。

##NextDenovo
#rice
nextDenovo ./rice.cfg
#rt
nextDenovo ./rt.cfg

NextDenovo 需要设置配置文件，以水稻测序数据为例，配置文件的内容如下：

job_type = local
job_prefix = rice
task = all
rewrite = yes
deltmp = yes
rerun = 3
parallel_jobs = 10
input_type = raw
input_fofn = input.fofn
workdir = rice
[correct_option]
read_cutoff = 2k
seed_cutoff = 0
genome_size = 380000000
seed_depth = 40
pa_correction = 2
sort_options = -m 1g -t 2
minimap2_options_raw = -t 8
correction_options = -p 15
[assemble_option]
minimap2_options_cns = -t 8
nextgraph_options = -a 1

input.fofn中是rice.fastq文件的绝对路径，对于其他数据，一般需要调整genome_size、seed_cutoff、seed_depth等参数，各个参数的详细说明，请点击NextDenovo的原链接。

##Wtdbg2
#rice
wtdbg2.pl -i ./rice.fastq -t 10 -x ont -g 380m -o rice_raw/wtdbg2
wtdbg2.pl -i ./canu/rice.trimmedReads.fasta.gz -t 10 -x ont -g 380m -o rice_corr/wtdbg2
#rt
wtdbg2.pl -i ./rt.fastq -t 10 -x sq -g 480m -o rt/wtdbg2

Wtdbg2在使用时需要注意设置数据的类型，-x中有以下几种参数对应不同的测序平台：rs指PacBio RSII、sq指PacBio Sequel、ccs指PacBio CCS、ont指Oxford Nanopore。

##Necat
#rice
necat.pl correct rice_config.txt
necat.pl assemble rice_config.txt
necat.pl bridge rice_config.txt

运行Necat也需要建立配置文件，Necat是专门用来组装ONT测序数据的，若是想要组装PacBio数据，需要用到另一个软件Mecat（https://github.com/xiaochuanle/MECAT），配置文件可以用命令来生成：

necat.pl config config.txt

以水稻为例，配置文件如下：

PROJECT=rice ## 输出结果文件的文件名
ONT_READ_LIST=read_list.txt ## nanopore数据路径及文件名
GENOME_SIZE=380000000  ## 基因组大小
THREADS=10 ## 线程数量
MIN_READ_LENGTH=5000  ## 最短length
PREP_OUTPUT_COVERAGE=40  ## 设定corrected reads的覆盖度
OVLP_FAST_OPTIONS=-n 500 -z 20 -b 2000 -e 0.5 -j 0 -u 1 -a 1000
OVLP_SENSITIVE_OPTIONS=-n 500 -z 10 -e 0.5 -j 0 -u 1 -a 1000
CNS_FAST_OPTIONS=-a 2000 -x 4 -y 12 -l 1000 -e 0.5 -p 0.8 -u 0
CNS_SENSITIVE_OPTIONS=-a 2000 -x 4 -y 12 -l 1000 -e 0.5 -p 0.8 -u 0
TRIM_OVLP_OPTIONS=-n 100 -z 10 -b 2000 -e 0.5 -j 1 -u 1 -a 400
ASM_OVLP_OPTIONS=-n 100 -z 10 -b 2000 -e 0.5 -j 1 -u 0 -a 400
NUM_ITER=2
CNS_OUTPUT_COVERAGE=30   ## 最长的30X数据用于后续组装
CLEANUP=1
USE_GRID=false
GRID_NODE=0
GRID_OPTIONS=
SMALL_MEMORY=0
FSA_OL_FILTER_OPTIONS=
FSA_ASSEMBLE_OPTIONS=
FSA_CTG_BRIDGE_OPTIONS=
POLISH_CONTIGS=true  ## 设置最后bridge结果是否polish

SMARTDenovo和HiFiasm的运行比较简单，而且不需要调整参数，建议以后的开发者都向他们看齐。

##SMARTDenovo
#rt
mkdir smartdenovo
cd smartdenovo
smartdenovo.pl -p smartdenovo -c 1 ./rt.fastq > smartdenovo.mk
make -f smartdenovo.mak

Hifiasm是专门用来组装HiFi数据的软件，也是比较好用的。

##Hifiasm
#rt
mkdir hifiasm
cd hifiasm
hifiasm -o rt -t 10 ./rt.fastq
awk ‘/^S/{print “>” $2;print $3}’ rt.bp.p_ctg.gfa > hifiasm.fasta

初步组装结果比较

等待所有软件运行完成后，我们需要根据组装结果，选择最合适的作为输入数据，然后再运行GALA。
水稻ONT测序数据组装结果：

Statistics	Canu	Flye	Wtdbg2	Miniasm	Necat	NextDenovo
Contigs	81	49	1128	62	66	23
Genome Size	383,971,102	380,844,286	357,519,400	382,869,794	382,439,828	380,767,888
N50	26,421,520	26,556,431	4,039,823	20,917,597	19,897,488	24,838,933
L50	7	6	26	8	8	6

经Canu纠错后的水稻基因组组装结果：

Statistics	Flye	Wtdbg2	Miniasm
Contigs	73	2166	74
Genome Size	376,571,620	417,867,573	380,573,187
N50	1,343,061	1,680,652	17,003,590
L50	12	60	9

桃金娘HiFi测序数据组装结果：

Statistics	Canu	Flye	SMARTDenovo	NextDenovo	Hifiasm
Contigs	2984	2118	225	125	413
Genome Size	756,136,501	567,348,138	470,579,731	468,746,210	506,040,249
N50	7,433,737	785,927	4,922,020	10,942,526	32,554,573
L50	22	131	29	12	7

目前来说，我们主要根据N50和L50判断组装结果的连续性，连续性越高，后续聚类结果的质量也会越好。从表中可知，对于ONT数据来说，NextDenovo的组装结果连续性最好，而Hifiasm最适合组装HiFi数据。然后Flye、Wtdbg2、Miniasm组装纠错后的数据在表中看不出什么优势，甚至单从连续性来看会变差，但是经过我们后续的测试，发现纠错后的数据在聚类的过程中不容易出错，会极大的改善染色体聚类结果，因此还是建议对ONT数据先纠错然后再组装。由HiFi数据的组装结果可知，所有的组装结果中contigs的数量都远远超过染色体的实际数目，因此组装结果中包含大量的零散判断，这在后续的步骤中需要额外关注。

总结

以上便是第一步预组装的全部内容，这里的步骤仅供参考。进行预组装并非简单的公式化的跑流程，而且需要了解自己想要组装的数据和用到的组装软件，根据物种以及测序数据的实际特征，灵活的选择最适合的组装软件，才能得到高质量的组装结果。