高通量测序的数据处理与分析指北(二)--宏基因组篇
宏基因组篇
前言
之前的一篇文章已经从生物实验的角度讲述了高通量测序的原理，这篇文章旨在介绍宏基因组二代测序数据的处理方式及其原理。在正文开始之前，我们先来认识一下什么是宏基因组。以我的理解，宏基因组就是某环境中所有生物的基因组的合集，这个环境可以是下水道，河流等自然环境，也可以是人体内肠道，口腔等体环境。而宏基因组中的生物往往指的是微生物，如真菌，细菌，病毒，古细菌。
我们这里主要以肠道微生物为例，也就是人体内肠道的宏基因组。肠道菌群的测序样本往往是粪便样本，现在主流的测序方式有两种：一种是16sRNA测序，一种是WGS(Whole Genome Sequencing) 全基因组测序。WGS测序数据量更大，所包含的信息更多，能注释出物种-样本的丰度矩阵，也能注释出基因-样本的丰度矩阵。而16sRNA测序测的是细菌核糖体RNA中的小亚基，这个小亚基的沉降系数是 16s，故被称为 16s RNA，这个16s RNA有一段非常保守的序列和一段变异序列，可以根据16s RNA 的变异度来进行物种分类，所以16s RNA数据往往只能注释出物种-样本的丰度矩阵。
原理介绍
之前文章中也提到了，由于测序技术的限制，目前二代测序只能测较短的碱基片段，所以需要对基因进行碎片化，我们要思考的问题就是这些碎片化的基因如何重新拼回到完整的基因组或者这些碎片化的基因如何确定其属于什么物种从而得到物种的丰度矩阵。
目前对宏基因组原始数据如何注释到物种的方法有两类主流方法，一类是基于bin进行物种注释的方法，一类是不基于bin进行物种注释的方法
基于bin的物种注释
基于bin的物种注释的代表软件有 metawrap，metabat2等。
在宏基因组的原始数据也就是fastq数据中，含有大量的read序列，首先是将read按照序列拼接成contigs，如图所示，上面的的read按照序列重合程度拼接成下面的contigs。
-metagenome/image-20230415213552290.png)
然后把相类似的contigs归为一个bin，而具体如何归bin的方法各种软件所用的原理都有些区别，这里介绍两种方法，也是这个视频中提到的两种分类的方法，第一种是依据四碱基频率来进行区分，所谓四碱基频率就是ATGC四个碱基为一组，共256种碱基组合，同一种物种的这256碱基组合的频率是相似的，并且物种亲缘关系越远则四碱基频率差距越大，故这一个256维的向量可以进行PCA降维，然后用聚类方法将类似的contigs聚到一起作为一个bin。
-metagenome/image-20230416123704317.png)
第二种方法是基于测序深度的，他的基本原理是由于不同的物种基因组大小不同，而同一种物种的基因组大小是类似的，因此可以根据contigs的深度来判断其是否为同一个物种，物种的基因组越大，在随机打碎DNA时产生的碎片越多，read数越高，最后通过read拼接而成的contigs的深度越大。
总而言之，bin就是一堆亲缘关系较近的contigs的合集，也可以视为一个物种基因组的草图。
得到高质量bin后就是对他进行基于数据库的注释，将能注释出来的bin注释出来。而bin的丰度，也就是物种的丰度的计算方式就是bin上每个碱基的深度除以bin序列长度。这个计算方式不太确定，推测的，暂时没找到资料
不基于bin的物种注释
基因bin的物种注释更加准确，但是耗时更长，这里再介绍一类直接从read比对数据得到物种丰度的宏基因组数据处理的方式，代表的软件有 kraken，metaphlan等。
-metagenome/image-20230416130355219.png)
这里主要以kraken的原理为例，它实际上就是将read 分成了多个 k-mers。 这个k-mers的意思就是是k bp长度的子序列，只不过这个子序列覆盖了read所有碱基，如150bp的read 能拆分出 150-31+1 个 31-mer, k-mers中的k长度是自定义的，默认是31，然后将这些 k-mers去跟数据库比对，k-mers对上最多的分支就作为这个read的物种分类，如上图，这个序列就是被认为是来自与4号物种的序列。同样的，将每个物种比对上的read数量除以其基因组长度就得到了其丰度。