二代测序基因组拼接

生信喵

一、案例介绍
1.1 案例介绍
      该文章中对 20 个细菌基因组进行测序，每个样本分别进行了 illumina，pacbio 以及 nanopore测序。比较三种数据的拼接结果。其中两株细菌已包含发表出来的全基因组序列。
《Comparison of long-read sequencing technologies in the hybrid assembly of complex bacterial genomes》
      1. Raw sequencing data: NCBI BioProject Accession PRJNA422511(https://www.ncbi.nlm.nih.gov/bioproject/PRJNA422511).
      2. Assemblies: FigShare doi https://doi.org/10.6084/m9.figshare.7649051 (https://doi.org/10.6084/m9.figshare.7649051).
      3. NCBI GenBank reference sequences:
      a. CFT073: NC_004431.1 (chromosome)
      b. MGH78578: NC_009648.1 (chromosome); NC_009649

1.2 获取 PRJNA422511 项目数据
      下载测序数据只要获得该数据在 SRA 数据库中对应的 SRA 号即可，一般会在文章中的 Data部分。如果存在多样本，则需要得到 PROJECT 号，在 PROJECT 号下面找对应的数据。
      https://www.ncbi.nlm.nih.gov/bioproject/PRJNA422511

1.3 下载数据
      这部分和上部分一样，大家可以往前找。

二、软件安装
      可以使用前一起安装好。

1. #使用 bioconda 进行安装
   
2. conda install -y mamba
   
3. mamba install -y spades
   
4. mamba install -y soapdenovo2
   
5. mamba install -y flye
   
6. mamba install -y canu
   
7. mamba install -y wgsim
   
8. mamba install -y megahit
   
9. mamba install -y gapfiller
   
10. mamba install -y jellyfish
    
11. mamba install -y gapfiller
    
12. mamba install -y soapdenovo2-gapcloser
    
13. mamba create -n medaka -y medaka
    
14. mamba create -n quast -y quast
    
15. mamba install -y genomescope2
    
16. mamba create -n unicycler -y unicycler

复制代码

三、数据处理
      这部分和上部分质控一样，大家可以往前找。

四、spades 拼接
4.1 软件介绍
      SPAdes 序列拼接软件是序列拼接软件中的后起之秀，拼接效果很不错，目前很多拼接软件已经都不在更新了，而 spades 却持续进行更新。该软件使用非常简单，并且支持多种数据格式，还支持混合拼接，使用起来非常方便。对于小基因组拼接有很好的拼接效果，不过经过更新，目前对于大基因组，甚至多倍体的基因组也有不错的效果。
      SPAdes 软件不仅支持 illumina 测序数据，而且还可以用于 Ion Torrent 测序数据，PacBio 测序数据、sanger 数据，Nanopore。并且可以加入其它序列拼接结果，作为辅助。使用简单，非常适合用于混合组装来改善拼接效果。软件使用起来也并不难。因为之前对于 Ion Torrent测序数据并没有官网软件，这里经过测试，发现 SPAdes 非常适合 Ion Torrent 数据的拼接。而且，SPAdes 还可以用于二倍体杂合基因组的拼接。

4.2 选项参数
      -o 为输出目录，需要提前建立一个目录，否则会报错
      -sc 指定输入文件为 singlecell 单细胞测测序数据
      --iontorrent 指定输入文件类型为 Ion Torrent 测序数据
      --test 会以软件自带的测试数据运行程序
      -h 或者--help 给出帮助信息
输入文件：
      --pe<#>-1，其中井号表示第几个文库，第一 pairend 文库就可以写成--pe1-1 后面接 reads1文件。
      --pe1-2 reads2 文件。以此类推，如果是大片段的 matepair 文库，就使用--mp1-1，--mp1-2。
      如果觉得 matepair 文库质量很高，文库的 SD 偏差小，可以使用 hqmp 选项。
      最新版本的 SPAdes，支持 illumina 最新的 NextSeq® Long Mate Pair libraries，所以可以使用--nxmate 选项。--sanger 指定 sanger 数据输入，
      --pacbio 指定 pacibo 数据输入，
      --nanopore 指定 nanopore 数据输入。
      --trusted-contigs 指定输入可信的其他软件拼接的 contig 结果
      --untrusted-contigs 指定输入不一定可信的其他软件拼接的 contig 结果
      SPAdes 是不支持这些数据单独拼接的，这些都是辅助上面短序列拼接的结果。
控制流程执行功能：
      --only-error-correction 只做数据纠错
      --only-assembler 只组装，不做数据纠错
      --careful 减少错误和插入缺失，添加此选项，会消耗更多的时间
高级选项：
      --dataset SPAdes 是可以将数据写到 YAML 格式的配置文件中的，适合一次很多数据。不过YAML 格式也不容易阅读，一般我们还是直接在命令行直接写
      -m 用于内存限制
      -t 控制线程数，默认 16 个
      -k kmer 数，一次可以输入多个，用逗号分隔，kmer 最大为 127，一般自动选择即可
      --cov-cutoff 覆盖率阈值，过滤掉费盖度过低的结果，默认不执行
      --phred-offset phred 质量体系，在数据纠错中会用到，现在 illumina 数据一般采用 phred 33

4.3 使用案例

1. spades.py -o spades_result -t 12 -1 ../data/clean.1.fq.gz -2 ../data/clean.2.fq.gz  1>spades.log 2>spades.err
   
2. spades.py --isolate -o spades_result -t 24 -k 45,95,12 -1 ../data/clean.1.fq.gz -2 ../data/clean.2.fq.gz  1>spades.log 2>spades.err

复制代码

五、SOAPdenovo 拼接
5.1 软件介绍
      SOAPdenovo 是由华大基因开发的 SOAP 软件包的一部分，SOAPdenovo 主要用于短序列reads 拼接，尤其是 illumina 测序数据。从小的细菌基因组到大的动植物基因组，人基因组都适用。已经成功应用于大熊猫基因组，黄瓜基因组等众多基因组的拼接中。
      SOAPdenovo 的一个优点是使用起来比较简单，但是却拥有很好的拼接效果，尤其在基因组构建 Scaffold 方面，效果很好。对于内存控制的也比较好。通常只要给软件输入测序的数据，即可拼接出很好的全基因组。

5.2 配置文件

1. git clone https://github.com/aquaskyline/SOAPdenovo2.git
   
2. cd SOAPdenovo2
   
3. make
   
4. /share/home/xiehs/biosoft/SOAPdenovo2 写入环境变量
   
5. #安装完毕

复制代码

1. #SOAPdenovo配置文件
   
2. max_rd_len=150
   
3. [LIB]
   
4. avg_ins=439
   
5. reverse_seq=0
   
6. asm_flags=3
   
7. rank=1
   
8. pair_num_cutoff=3
   
9. #q1=#reads1文件路径 (),这个需要替换，不要那么教条直接复制这段
   
10. #q2=#reads2文件路径 (),这个需要替换，不要那么教条直接复制这段
    
11. q1=/share/home/xiehs/05.assembly/data/clean.1.fq.gz
    
12. q2=/share/home/xiehs/05.assembly/data/clean.2.fq.gz

复制代码

5.3 选项参数。
      -s STR 配置文件
      -o STR 输出文件的文件名前缀
      -g STR 输入文件的文件名前缀，这个主要用在分布运行程序的时候。
      -K INT 输入的 K-mer 值大小，默认值 23，取值范围 13-63
      -p INT 程序运行时设定的线程数，默认值 8
      -R 利用 read 鉴别短的重复序列，默认值不进行此操作
      -d INT 去除频数不大于该值的 k-mer，默认值为 0
      -D INT 去除频数不大于该值的由 k-mer 连接的边，默认值为 1，即该边上每个点的频数都小于等于 1 时才去除
      -M INT 连接 contig 时合并相似序列的等级，默认值为 1，最大值 3。
      -F 利用 read 对 scaffold 中的 gap 进行填补，默认不执行
      -u 构建 scaffold 前不屏蔽高覆盖度的 contig，这里高频率覆盖度指平均 contig 覆盖深度的 2 倍。默认屏蔽
      -G INT 估计 gap 的大小和实际补 gap 的大小的差异，默认值为 50bp。
      -L 用于构建 scaffold 的 contig 的最短长度，默认为：Kmer 参数值 ×2
      -k INT map 步骤中 kmer 的大小，默认是和 K 一样的 kmer 大小
      -N INT 基因组大小
      -V 输出可视化的组装信息
5.4 使用案例：

1. #soapdenovo2运行
   
2. mkdir kmer45  
   
3. echo "SOAPdenovo-63mer all -s lib.list -K 45 -o kmer45/kmer45 -D 1 -d 1 -u 2 -p 24 >kmer45.log" >soapdenovo.sh
   
4. sh soapdenovo.sh
   
5. echo "SOAPdenovo-127mer all -s lib.list -K 125 -o kmer125/kmer125 -D 1 -d 1 -u 2 -p 24 >kmer45.log" >soapdenovo1.sh
   
6. echo "SOAPdenovo-127mer all -s lib.list -K 85 -o kmer85/kmer85 -D 1 -d 1 -u 2 -p 24 2>kmer85.log" >soapdenovo2.sh

复制代码

六、补洞
6.1 为什么存在“洞”区域
      基因组上的洞也叫做 GAP，是由Ｎ碱基构成的。这个 N 碱基是在利用 reads 之间的 pairend关系将 contig 构建成 scaffold 的过程中产生的。将 reads 之间的 insertsize 值减去两侧已有碱基数目，然后在中间填上 N 碱基，这就是 gap 的来源。因为这个 insertsize 是一个固定的值，所有 N 区域都用这个固定的值，所以我们要知道一点是，得到这个 Gap 数是不精确的，存在一定的偏差。也就是虽然基因组上这个洞有 300 个 N 碱基，但是实际上，并不是就缺了刚好 300bp 的碱基，真实的基因组序列与这个数目存在一定的偏差，可能要大于 300bp，也可能远小于这个数目。

生信喵

6.2 补洞原理
      在了解了以上内容之后，我们就可以开始来进行补洞操作了。
      补洞有多种策略，首先最理想的策略是加入测序长片段，例如 sanger 测序数据或者 pacbio测序数据。在洞的上下游各取一段非 N 碱基序列用来设计引物，然后从剩余样品中 PCR 扩增出这段区域，采用 sanger 测序法，将这段区域测序出来。如果能 PCR 出来就说明这个区域连接的没有问题，我们将 PCR 扩增出来的区域直接填到这个洞区域即可。但是如果洞区域过程，例如大于 2K，PCR 出来，还是无法测通，这就不好办了，这个时候需要测序更长的序列才行，所以，现在有采用加入 pacbio 测序数据的方式来补洞。
      
      还有一种策略是不加入额外长序列，直接用组装的 reads 来进行补洞。
      首先是读取所有 reads，将具有 pair-end 关系的 reads 转换后存储起来，然后是根据 N 碱基对 scaffolds 进行分割，接下来是将具有 pairend 关系的 reads 比对到 scaffold 上，得到落在gap 区域内的 reads。mapping 到 scaffold 的 pairend reads 有多种情况，例如图中所示的12345 五种情况，我们挑选那种 pairend 的 reads，其中一条 reads 比对到洞的两侧，另一条没有比对上 scaffold 上的 reads，我们根据 reads 的 pairend 关系，认为另一条落在了洞区域。其中只有第一种和第四种的 pairend reads 被利用到了。
      补洞会让基因组序列更加完整，但是我们同时必须意识到，这些区域本来就是比较难拼接的区域，所以，补洞区域序列的准确性要稍微低一些，这在变异检测的时候会有一些影响。所以，如果检测到的变异位点发生在补洞区域，一定要进行进一步的验证。

1. #补洞
   
2. GapCloser -a kmer45/kmer45.scafSeq -b lib.list -o kmer45.fill.fa -l 100 -p 25
   
3. -t >gapcloser.log
   

复制代码

七、常见问题
7.1 影响拼接的因素有那些？
      影响基因组拼接的因素很多，包括内在因素来自基因组本身的重复序列，多倍体杂合，还包括外在因素测序错误，测序饱和度等。
      1、重复序列是基因组拼接最大的影响因素。测序数据无法跨过“重复序列”区域，遇到重复区则“断开”；
      2、多倍体杂合：多倍体需要测序更多的数据，杂合造成更多的“气泡”；
      3、测序错误：测序错误导致 kmer 之间无法连接，出现更多的“道路”；
      4、GC 含量：高低 GC 导致测序无法测序到或者测序覆盖不均匀；
      5、宿主污染：宿主污染导致测序数据中绝大部分数据为宿主序列，无法测序足够的目标区域数据；
      6、覆盖深度：测序深度不足或者不均匀，会影响软件判断重复序列区域；
      7、随机性：测序 reads 之间需要更多的 overlap 连接，随机性不好，或者 DNA 降解会造成过多的 duplication
      8、duplication：Duplication 会是数据在某一区域覆盖度增高，覆盖不均与，导致与重复区域相似，干扰软件判断重复序列。

7.2 如何改善拼接效果？
      知道了有哪些因素会影响到拼接效果之后，就可以来调整这些因素来改善序列拼接的效果。首先必须保证测序覆盖度。如果测序效果不好，可以设立更加严格的标准，重新过滤数据，利用 reads 与拼接好的数据进行短序列比对，重新估计更加精确的 insertsize 值，去除insertsize 偏差过大的成对reads。调整多个 kmer 值，以及其他阈值来改善组装结果的质量。对于 illumina 测序来说，增加大片段的文库会显著改善序列拼接的效果。有些软件还支持加入其他测序数据进行混合拼接，或者加入可信的近源参考序列，或者不同软件拼接出来的 contig 来进行辅助拼接。这些策略都可以尝试着去改善拼接结果。

7.3 不同测序平台之间测序数据混合拼接？
      将不同测序平台的数据进行混合拼接，是当前技术条件下改善序列拼接效果的一个很好的策略。我们可以使用 illumina 数据进行拼接，然后使用 sanger 测序来进行补洞。或者在细菌基因组测序中，利用 pacbio 或者 nanopore 测序长片段数据来确定 contig 之间的位置关系，或者用于补洞，这个在构建细菌完成图中效果比较显著。
      或者直接使用 pacbio 或者 nanopore 测序进行连接，使用二代测序进行纠错。

7.4 为什么不用短 reads 直接 overlap 拼接？
      这个主要是由于测序错误的影响，为了去除测序错误同时又想充分利用数据，所以就采取了将 reads 切割成 kmer 的方法。如果测序 reads 中没有测序错误，那么直接用 reads 进行overlap 就行了。

7.5 如何选择 kmer 的大小？
      假设测序数据中没有测序错误。那么 kmer 大小就取 reads 长度大小的奇数就行，其实无需取切割 reads 成 kmer 了。如果测序数据有很多错误，如果取大 kmer，包含测序错误的大 kmer，直接就丢掉了。那么取 kmer 值越小，可以利用的数据就越多。所以，测序错误率是影响 kmer大小非常重要的因素。测序准确率越高，选择大 kmer 越好，测序错误率越低，选择小 kmer越好。但实际上还需要综合考虑基因组本身特点，重复率，测序深度等因素，都会对 kmer取值造成一定影响。

7.6 为什么 kmer 只能是奇数？
      kmer 不能取偶数主要是由于回文序列的影响。因为在取 kmer 的时候，我们不仅要正常取一遍 kmer，还要从它的反向互补序列再取一次 kmer。如果是回文序列，在使用偶数的 kmer值的时候就会有问题。

7.7 如何判断基因组中是否有污染？
      首先看基因组大小，基因组有污染，则代表基因组是混合基因组，例如结果是有基因组 A和基因组 B 组成，假设 A 基因组大小为 5.4M，B 基因组大小为 4.6M。所以如果样品 A 拼接结果中混有样品 B，则基因组大小为 5.4M+4.6M。
      其次，由于两个样品测序时添加的样品数量不同，则测序深度不同。

7.8 有污染的样品如何处理？
      如果发现测序拼接得到的基因组中有污染，很难通过生物信息的方法完全进行拆分，建议重新取样测序。