基于二代全基因组测序的快速基因组组装技术流程
Technical Process of Rapid Genome Assembly Based on NGS Whole-genome Sequencing   

研究背景

近年来，二代 (高通量) 测序技术的快速发展以及测序成本的降低，为基因组数据的获 取提供了便利，也使得它目前已经广泛应用于解决系统发生关系、生物进化等领域的科学问题。然而，由于二代测序技术本身的一些局限性，例如测序序列长度远远小于分子可能的原有长度，且难以区分基因组中的重复序列；测序时覆盖度不均匀，影响组装质量；测序时存在一定比例的错误，影响组装结果；测序得到的大量数据也增加了计算的复杂性，故基于二代测序的基因组组装具有一定的难度，也会占用大量计算资源与时间。
        本文提出了一套基于二代全基因测序的快速基因组组装流程，它整合了一系列快速、高效的生物信息工具，可以从低深度测序的原始数据中快速组装基因组。使用该流程与技术，在普通配置的台式电脑上组装完成单个物种基因组 (大小一般在0.1-2 Gbp) 仅需2-24 h。该方法已在弹尾纲、昆虫纲、蛛形纲等类群中得到应用，组装的 基因组数据也已成功应用于六足类 (如虱目) 的分子系统学研究中。

仪器设备

1. 普通电脑 (最低4核8线程，16 G内存)

软件

本教程可以自行在服务器或个人电脑上安装好相关软件进行。相关软件的安装版本及使用说明如下。

1. 删除重复、质控、降低丰度：BBTools v38.32 (https://sourceforge.net/projects/bbmap/)
2. 组装 (Contig水平组装)：Minia v3.2.4 (https://github.com/GATB/minia) 或SPA des v3.14.1 (http://bioinf.spbau.ru/spades)
3. 去杂合：Redundans v0.14c (https://github.com/lpryszcz/redundans) 或cd-hit (http://www.bioinformatics.org/cd-hit/)
4. 序列延长 (Scaffold水平组装)：BESST v2.2.8 (https://github.com/ksahlin/BESST)
5. 填补缺口：GapCloser v1.12 (http://soap.genomics.org.cn/)
   
6. 基因组完整性评估：BUSCO v3.0.2 (http://busco.ezlab.org/)
   

也可在Windows系统预先安装好虚拟机VMware或VirtualBox，导入我们准备的虚拟镜像 (.vmdk)。该镜像已安装CentOS 7.3系统和所需的生物信息学工具，建议磁盘空间超过100 Gb (至少50 Gb)。CentOS 7.3系统具体信息及下载链接详见https://1drv.ms/f/s!Ak7sQGBGqlguq1gdk1KKIR7uwAfj。

实验步骤

该流程主要指基于二代全基因组测序的原始数据进行快速基因组组装。可以分为三个过程：数据获取、基因组组装和基因组完整性评估，具体操作流程如图1所示。


图1. 基于全基因组测序的基因组组装及基因组评估流程图

一、数据获取 (见图1)

二代全基因组测序一般是指通过Illumina高通量和BGI测序平台获取的全基因组序列，国内几所较为大型的测序公司，例如：北京贝瑞禾康生物技术有限公司、上海生工生物工程股份有限公司、北京诺禾致源科技股份有限公司、深圳华大生物科技有限公司等都可以进行二代全基因组测序。
      二代全基因组测序原始数据还可以从NCBI数据库中 (https://www.ncbi.nlm.nih.g ov/) 使用SRA Toolkit v2.9.0软件 (SRA工具包开发团队，https://trace.ncbi.nlm.nih.g ov/Traces/sra/sra.cgi?view=software) 下载，并将其转换为gzip fastq格式，代码如下：
fastq-dump --split-3 SRR****.sra -O out_path
      本文测试的数据均下载自NCBI数据库已公开发表的二代全基因组测序原始数据，具体数据信息见Zhang et al., 2019。

二、基因组组装 (见图1)

在正式基因组组装前，需要根据测序结果进行基因组大小、重复序列含量和物种杂合 度等信息的评估。基于这些信息，确定后续的组装策略，从而降低计算资源，节约时 间。本文提出的快速基因组组装流程包括删除重复、质控、降低丰度、组装 (Contig 水平组装)、去杂合、序列延长 (Scaffold水平组装) 和填补缺口。

1. 删除重复：BBTools (Bushnell, 2014) 组件拥有对高通量测序数据进行质量控制、 修剪、纠错、比对等功能，因此，本文使用BBTools组件中的clumpify.sh脚本对数据进行重复数据的删除。代码如下：
   clumpify.sh in1=illumina.R1.fq.gz in2=illumina.R2.fq.gz out1=1.clumped.fq.gz out2=2.clumped.fq.gz pigz dedupe
   注意：clumpify.sh脚本详细参数和具体操作步骤参考软件使用说明。
   
2. 质控：为了降低数据偏差，需要在正式组装前进行质控，目的是在减少错误的情 况下提高整体质量。
   质控的软件有很多，例如BBTools，fastp (Chen et al., 2018)，Trimmomatic (Bolger et al., 2014) 等。本文使用BBTools组件中的bbduk.sh脚本对数据进行质控：利用Phred算法将数据两侧修剪至Q15，并过滤小于15 bp以及超过连续5个N的序列，数据两端至少修剪10 bp的Poly A或Poly T，并校正重叠的成对读取。代码如下：
   bbduk.sh in1=1.clumped.fq.gz in2=2.clumped.fq.gz out1=1.trim.fq.gz out2=2.trim.fq.gz pigz ordered qtrim=rl trimq=20 minlen=15 ecco=t maxns=5 trimpolya=10 trimpolyg=10 trimpolyc=10
   注意：bbduk.sh脚本详细参数和具体操作步骤参考软件使用说明。
   
3. 降低丰度 (可选)：为了降低组装数据的复杂性和加快基因组组装速度，我们使用bbnorm.sh脚本 (BBTools组件之一) 对测序数据进行归一化处理，对序列深度在10x以上的高覆盖度序列进行深度为10x的二次采样，对于序列深度低于10x的低覆盖度序列则全部数据保留。代码如下：
   bbnorm.sh in1=1.trim.fq.gz in2=2.trim.fq.gz out1=1.nor.fq.gz out2=2.nor.fq.gz target=10 min=2 histcol=2 khist=khist.txt peaks=peaks.txt
   注意：bbnorm.sh脚本详细参数和具体操作步骤参考软件使用说明。
   
4. 组装 (Contig水平组装)：Contig组装是基因组组装的第一步，也是非常重要的一步，在实际组装过程中可以根据具体的要求选择不同的组装软件进行基因组组装。Minia3 v3.2.4 (Chikhi and Rizk, 2013) 是内存资源消耗最低的基因组组装工具，组装速度极快，但是组装质量较低；SPAdes v3.14.1 (Bankevich et al., 2012) 组装速度较慢，计算资源需求多，但是组装质量高，全自动化流程。
   1)

   使用Minia3 v3.2.4软件基于GATB-Minia-Pipeline (详细参数和具体操作见https://github.com/GATB/gatb-minia-pipeline) 进行Contig组装。对于100 bp的读长选择21、41、61、81的K-mer值，对于150 bp的读长选择21、41、61、81、101、121的K-mer值。代码如下：
   for KMER in $(cat kmer.list)
       do
       minia -in reads.list -kmer-size $KMER -abundance-min 2 -out k$KMER -max-memory $MEMORY
       rm -rf dummy* *unitigs* *h5 trashme*
       cp ../reads.list reads.list
       for n in 1 2 3; do echo "k$KMER.contigs.fa" >> reads.list; done
       done
   注意：以上操作步骤会根据基因组大小调整K-mer值，Minia3 v3.2.4的详细参数和具体操作步骤参考软件使用说明。
   
   2)

   若想要获得较高质量的组装结果，可以考虑使用SPAdes v3.14.1软件进行组装。代码如下：
   spades.py -1 1.nor.fq.gz -2 2.nor.fq.gz -o spades -t $THREADS
   注意：SPAdes v3.14.1的详细参数和具体操作步骤参考软件使用说明。
   

5. 去杂合：同一基因组区域的多条单倍体有时由于差异分化过大，组装过程中可能被组装为多条序列，一般组装时仅保留1条。因此，需要对初步组装后的数据进行去杂合。Redundans v0.13c (Pryszcz et al., 2016) 软件速度慢，但是去杂合质量高；而cd-hit (Fu et al., 2012) 软件速度快，但是去杂合质量低。为了得到高质量的组装结果，本文使用Redundans v0.13c进行去冗余操作。代码如下：
   redundans.py -v -f contigs.fa -o out --noscaffolding --nogapclosing --identity 0.7
   注意：identity参数指序列相似度，一般设置在0.5-0.9之间；Redundans v0.13c的详细参数和具体操作步骤参考软件使用说明。

6. 序列延长 (Scaffold水平组装)：本文使用SOAPdenovo2 (Luo et al., 2012) 套件中的BESST v2.2.8 (Sahlin et al., 2014)，通过使用碱基N将Contig连接延长组装为Scaffolds。代码如下：
   runBESST -c contig.fa -f map.bam -o output -orientation fr --iter 10000
   注意：BESST v2.2.8的详细参数和具体操作步骤参考软件使用说明。
7. 填补缺口：使用SOAPdenovo2组件中的GapCloser v1.12 (Xu et al., 2020)，将上一步骤所生成的Scaffolds中的碱基N替换成碱基A\T\G\C。代码如下：
   GapCloser -a Scaffolds.fa -b gapcloser.config -o scaffolds.gapcloser.fa -l $READ_LENGTH -t $THREADS
   注意：GapCloser v1.12的详细参数和具体操作步骤参考软件使用说明。

三、基因组完整性评估 (见图1)

组装评估结果可查看scaffold和contig N50长度、scaffolds/contigs数量等。使用BUSCO v3.0.2 (Waterhouse et al., 2018) 进行完整性评估，通常认为基因组BUSCO 完整性> 90%的组装结果较好，完整性> 50%就可以进行系统发育分析。
查看组装指标：seqkit stat -a genome.fa或stats.sh in=genome.fa
查看多个基因组组装后的结果：statswrapper.sh in=genome1.fa,genome2.fa
完整性评估：run_BUSCO.py -i scaffolds.fa -c 8 -m geno -o $SPECIES -l $LINEAGE
注意：$LINEAGE为USCO参考数据集，需要根据物种选择适合的参考数据集 (详见：https://busco-data.ezlab.org/v5/data/lineages/)，结果在run_$SPECIES目录下 的short_summary_result.txt查看。BUSCO v3.0.2详细参数及具体操作参考软件使用说明。
上述流程可以使用Zhang et al. (2019) 的Script1进行，包括基因组快速组装的所有步骤，即删除重复 (BBTools组件之clumpify.sh)、数据质控 (BBTools组件之 bbduk.sh)、降低丰度 (BBTools组件之bbnorm.sh)，Contig水平组装 (Minia3)，去杂合 (Redundans)，Scaffold水平组装 (BESST)，填补缺口 (GapCloser) 和完整性评估 (BUSCO)。详细内容及脚本见https://github.com/xtmtd/PLWS#，可实现一键操作：
sh script1_Genome_assembly1.sh forward_reads_file reverse_reads_file
产生的主要结果文件如表1所示。

表1. 产生的主要结果文件及内容

注：实际操作过程中，需输入Script1和正反向原始序列所在路径。在执行该脚本时，根据实际需要选择或输入相关参数。

数据分析

本文使用Zhang et al. (2019) 中的数据，通过PLWS中script1_Genome_assembly1.sh成功地快速组装了基因组，并用于后续通用单拷贝核基因 (Universal Single-Copy Orthologs, USCOs) 和超保守元件 (Ultraconserved Elements, UCEs) 的提取。基于该方案组装的基因组数据已成功应用于六足类的分子系统学研究 (具体数据和实验效果详见Zhang et al., 2019)。

致谢

感谢国家自然科学基金项目 (31970434、31772491) 的资助。该方案摘本实验室已发 表的文章Zhang, F., Ding, Y., Zhu, C. D., Zhou, X., Orr, M. C., Scheu, S. and Luan, Y. X. (2019). Phylogenomics from low-coverage whole-genome sequencing. Methods Ecol Evol 00: 1-11.

竞争性利益声明

作者声明没有利益冲突。

参考文献

1. Bankevich, A., Nurk, S., Antipov, D., Gurevich, A. A., Dvorkin, M., Kulikov, A. S., Lesin, V. M., Nikolenko, S. I., Pham, S., Prjibelski, A. D., Pyshkin, A. V., Sirotkin, A. V., Vyahhi, N., Tesler, G., Alekseyev, M. A. and Pevzner, P. A. (2012). SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing. J Comput Biol 19(5): 455-477.
2. Bolger, A. M., Marc, L. and Bjoern, U. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30(15): 2114-2120.
3. Bushnell, B. (2014). BBtools. Retrieved from https://sourceforge.net/projects/bbmap/.
4. Chen, S., Zhou, Y., Chen, Y. and Gu, J. (2018). fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 34(17): 884-890.
5. Chikhi, R. and Rizk, G. (2013). Space-efficient and exact de Bruijn graph representation based on a Bloom filter. Algorithm Mol Biol 8(1): 22.
6. Fu, L., Niu, B., Zhu, Z., Wu, S. and Li, W. (2012). CD-HIT: accelerated for clustering the next-generation sequencing data. Bioinformatics 28(23): 3150-3152.
7. Luo, R., Liu, B., Xie, Y., Li, Z., Huang, W., Yuan, J., He, G., Chen, Y., Pan, Q., Liu, Y., Tang, J., Wu, G., Zhang, H., Shi, Y., Liu, Y., Yu, C., Wang, B., Lu, Y., Han, C., Cheung, D. W., Yiu, S. M., Peng, S., Zhu, X., Liu, G., Liao, X., Li, Y., Yang, H., Wang, J., Lam, T. W. and Wang, J. (2012). SOAPdenovo2: an empirically improved memory-efficient short-read de novo assembler. Gigascience 1(1): 18.
8. Pryszcz, P. L. and Gabaldón, T. (2016). Redundans: an assembly pipeline for highly heterozygous genomes. Nucleic Acids Res 44(12): e113.
9. Sahlin, K., Vezzi, F., Nystedt, B., Lundeberg, J. and Arvestad, L. (2014). BESST–efficient scaffolding of large fragmented assemblies. BMC Bioinformatics 15: 281.
10. Waterhouse, R. M., Seppey, M., Simão, F. A., Manni, M., Ioannidis, P., Klioutchnikov, G., Kriventseva, E. V. and Zdobnov, E. M. (2018). BUSCO applications from quality assessments to gene prediction and phylogenomics. Mol Biol Evol 35(3): 543-548.
11. Xu, M., Guo, L., Gu, S., Wang, O., Zhang, R., Peters, B. A., Fan, G., Liu, X., Xu, X., Deng, L. and Zhang, Y. (2020). TGS-GapCloser: a fast and accurate gap closer for large genomes with low coverage of error-prone long reads. GigaScience 9: 9.
12. Zhang, F., Ding, Y., Zhu, C. D., Zhou, X., Orr, M. C., Scheu, S. and Luan, Y. X. (2019). Phylogenomics from low-coverage whole-genome sequencing. Methods Ecol Evol 00: 1-11.