二代测序宏基因组拼接

生信喵

一、宏基因组拼接原理
       基因组拼接一直是整个基因组数据分析中最重要和最核心的工作，因为基因组包含了一个物种全部的遗传信息。得到的基因组越完整，包含的基因组信息也越多，对于后续对整个基因组的功能分析，变异检测都有非常大的帮助。由于基因组本身具有的高度重复序列，多倍体杂合位点，低复杂度区域以及测序错误等诸多条件的影响，基因组拼接一直是一项非常复杂且困难的工作。尤其是基因组重复序列的影响，一直是二代短读长测序最难解决的问题，尽管后来基于二代测序数据开发除了一些辅助拼接方案，例如大片段文库，Optical mapping光学图谱，三位基因组等辅助方案，都无法彻底解决基因组拼接难题。而利用 nanopore 长度长测序，将革命性的解决重复序列对于基因组拼接的影响。
       nanopore 测序数据不适合利用 NGS 数据组装的 de Bruijn 图法进行组装，主要存在两方面的原因。第一，de Bruijn 图法等方法依赖测序 reads 拆分的 k-mer 测序准确，而高错误率的MinION 测序 reads 不能保证这一点；第二，de Bruijn 图的结构不适用长 reads。
       nanopore 测序数据的长 reads 更适合 Sanger 测序时期基于有 overlap 的共有（consensus）序列组装的方法。这种算法被称为 Overlap Layout Consensus，简写为 OLC 算法，主要分为三步：
       （1）Overlap：，对所有 reads 进行两两比对，找到片段间的重叠信息；
       （2）Layout：根据得到的重叠信息将存在的重叠片段建立一种组合关系，形成重叠群，即Contig；
       （3）根据构成 Contig 的片段的原始质量数据，在重叠群中寻找一条质量最重的序列路径，并获得与路径对应的序列，即 Consensus。
       宏基因组拼接 Metagenomically-Assembled Genomes (MAGs)，与单个物种拼接里类似，但由于宏基因组属于混合样品，并且各个样品之间丰度不同，因此拼接难度更大。通过宏基因组拼接，可以得到更加完整的基因组序列和基因序列。纳米孔测序的宏基因组拼接，由于测序长度更长，可以直接拼接出一些细菌完整的基因组序列，而这些细菌往往无法通过传统纯培养的方法获得，这为获得无法纯培养样品得到完整基因组序列提供了新思路。
二、为什么宏基因组拼接比较难？
       基因组拼接一直是生物信息分析当中最重要，也是最困难的过程。尽管纳米孔的长读长测序显著改善了基因组拼接的效果，但基因组拼接依然有很多问题当前技术无法解决。在现阶段，拼接单个物种基因组依然存在较大的技术难度，而宏基因组属于“混合”样品，里面包含了多个物种，因此难度更大。
       影响基因组拼接的因素很多，包括内在因素来自基因组本身的重复序列，多倍体杂合，还包括外在因素测序错误，测序饱和度等。而对于宏基因来说还存在各个样品测序深度不均匀的影响，这些都会影响到宏基因组的拼接效果。
       在以上所有因素当中，重复序列是基因组拼接最大的影响因素。而在宏基因组拼接中的影响比普通的动植物基因组更大。这是因为对于单个物种来说，所有的重复序列都只属于物种内重复，例如一个基因的多个拷贝或者是一些非编码的串联重复等等。而对于宏基因组来说，由于属于“混合样品”，群落中存在大量微生物，微生物之间的基因交流（如水平基因转移），促使物种间大量的相似甚至相同序列片段出现，从而形成物种间重复。因此，宏基因组实际上是受物种内重复（Intragenomic Repeats）和物种间重复（Intergenomic Repeats）双重因素影响，拼接难度更大。
      
       物种内重复（上）与物种间重复（下）示意图
       另外，由于群落中的样品来自同一种或者同一属或更高层级，序列之间具有相似性，例如同样来自肠杆菌属，会将来自不同物种的序列拼接在一起，形成嵌合体，造成拼接错误。
三、二代模拟数据宏基因组拼接
3.1 混合模拟数据
       人为添加 10 种微生物，其中包括 8 株细菌，两株真菌。分为两种模式，一种按比例平均分配，称为 Even 数据集，8 株细菌各占 8%，2 株真菌各占 4%。另一种按照对数进行分配，称为 Log 数据集。并且包括 illumina 与 nanopore 数据。

      
       mock 数据集

1. https://github.com/LomanLab/mockcommunity
   
2. https://lomanlab.github.io/mockcommunity/

复制代码

      选择 Enterococcus faecalis ERR2935850，Escherichia coli ERR2935852 ， Salmonella enterica ERR2935848 以及 Staphylococcus aureus ERR2935849 四种细菌进行混合。然后对混合数据进行拼接。

1. #下载数据
   
2. axel -n 100 https://sra-pub-run-odp.s3.amazonaws.com/sra/ERR2935850/ERR2935850
   
3. axel -n 100 https://sra-pub-run-odp.s3.amazonaws.com/sra/ERR2935852/ERR2935852
   
4. axel -n 100 https://sra-pub-run-odp.s3.amazonaws.com/sra/ERR2935848/ERR2935848
   
5. axel -n 100 https://sra-pub-run-odp.s3.amazonaws.com/sra/ERR2935849/ERR2935849

复制代码

      下载基因组序列

1. https://github.com/al-mcintyre/mCaller_analysis_scripts

复制代码

1. #下载单株基因组序列
   
2. git clone https://github.com/al-mcintyre/mCaller_analysis_scripts.git
   
3. #合并参考序列
   
4. cd /share/home/xiehs/18.mags/mCaller_analysis_scripts/assemblies
   
5. cat efaecalis_pb.fasta k12_pb.fasta senterica_pb.fasta saureus_pb.fasta >>ref.fasta
   
6. mv ref.fasta ../../
   
7. cd ../../

复制代码

1. #数据转换
   
2. mv ERR2935848 ERR2935848.sra
   
3. mv ERR2935849 ERR2935849.sra
   
4. mv ERR2935850 ERR2935850.sra
   
5. mv ERR2935852 ERR2935852.sra
   
6. 
   
7. fasterq-dump ERR2935848.sra
   
8. fasterq-dump ERR2935849.sra
   
9. fasterq-dump ERR2935850.sra
   
10. fasterq-dump ERR2935852.sra

复制代码

1. #随机抽样
   
2. #Enterococcus faecalis ERR2935850，
   
3. #Escherichia coli ERR2935852，
   
4. #Salmonella enterica ERR2935848
   
5. #Staphylococcus aureus ERR2935849
   
6. conda activate rnaseq
   
7. seqkit sample -p 0.03 ERR2935848_1.fastq >meta_1.fastq
   
8. seqkit sample -p 0.03 ERR2935848_2.fastq >meta_2.fastq
   
9. seqkit sample -p 0.02 ERR2935849_1.fastq >>meta_1.fastq
   
10. seqkit sample -p 0.02 ERR2935849_2.fastq >>meta_2.fastq
    
11. seqkit sample -p 0.01 ERR2935850_1.fastq >>meta_1.fastq
    
12. seqkit sample -p 0.01 ERR2935850_2.fastq >>meta_2.fastq
    
13. seqkit sample -p 0.05 ERR2935852_1.fastq >>meta_1.fastq
    
14. seqkit sample -p 0.05 ERR2935852_2.fastq >>meta_2.fastq

复制代码

3.2 基因组拼接

1. #基因组拼接
   
2. #写脚本 spades路径需要自己确定
   
3. echo "/share/home/xiehs/Software/miniconda3/envs/unicycler/bin/spades.py -o spades_assembly -t 24 --meta -1 meta_1.fastq -2 meta_2.fastq 1>spades.log 2>spades.err"> spades.sh
   
4. #后台运行
   
5. #bsub -q fat -n 24 -o %J.log -e %J.err sh spades.sh #可选集群作业
   
6. nohup sh spades.sh &

复制代码

3.3 拼接结果评估

1. #拼接结果评估
   
2. conda activate quast
   
3. nohup metaquast -o metaQUAST -t 12 ref.fasta spades_assembly/scaffolds.fasta &

复制代码

      
       拼接结果评估