CN108133122B - 基因聚类方法和基于该方法的宏基因组组装方法和装置 - Google Patents

Abstract

本发明公开了一种基因聚类方法和基于该方法的宏基因组组装方法和装置，所述基因聚类方法包括：根据每个基因序列或小基因簇在多个样品中的丰度谱，以及每个大基因簇在多个样品中的丰度谱，通过多元线性回归，计算每个大基因簇与每个基因序列或小基因簇的相关性；找到与每个大基因簇的相关性大于设定值的基因序列或小基因簇，并将所找到的基因序列或小基因簇合并至相关联的大基因簇中，组合成新的大基因簇。能够找到多个基因序列或小基因簇的多个源头大基因簇，弥补现有方法的不足，从而增加基因簇的大小，并且获得更好的基因组草图。

Description

基因聚类方法和基于该方法的宏基因组组装方法和装置

技术领域

本发明涉及基因序列信息分析技术领域，尤其涉及一种基因聚类方法和基于该方法的宏基因组组装方法和装置。

背景技术

宏基因组指来源于同一环境中所有微生物基因组的集合，如人肠道宏基因组表示人肠道中存在的所有微生物的基因组。快速发展的第二代DNA测序技术可以有效获取宏基因组数据，通过生物信息技术的进一步分析，解读微生物群落的组成以及解析微生物群落的功能。然而，由于宏基因组由多个基因组混合组成，从宏基因组中解析出单个基因组是研究的一个重点和难点。

目前，已有研究通过基因序列在多样品之间的丰度谱关联信息，将基因序列聚集成基因簇，进而对较大的基因簇进行独立的基因组组装，从而得到单个基因组草图(Nielsen H.B.等人,Nature Biotechnology,2014)。该技术具体步骤如下：(1)通过宏基因组组装、基因预测、基因序列去冗余获取多样品特有的非冗余基因序列集合；(2)将测序读长序列(Reads)比对到上一步获取的基因序列集合，计算每个基因在每个样品的丰度，得到基因在多样品中的丰度谱；(3)根据基因在多样品的丰度谱，计算基因与基因之间的相关系数，根据相关系数对基因进行聚类，将相关系数满足一定阈值的基因聚集成基因簇；(4)将测序Reads比对到较大的基因簇，获取单个基因簇特有的测序Reads，进行单独组装，获得单个基因组草图。

现有技术的关键假设是来源于同一DNA分子的基因序列，在多样品中的丰度谱表现出强的相关性。然而，有一些基因序列同时存在多个基因组，它们的丰度谱由多个基因组决定，并不与其中任何一个基因组有强相关性。现有技术无法将此类型的基因序列聚集到对应的基因簇，而且这些基因具有重要的生物功能，如持家基因和核糖体基因等，是基因组的重要组成。

发明内容

本发明提出一种基因聚类方法和基于该方法的宏基因组组装方法，本发明的基因聚类方法能够找到多个基因序列或小基因簇的多个源头大基因簇，弥补现有方法的不足，从而增加基因簇的大小，并且获得更好的基因组草图。

根据本发明的第一方面，本发明提供一种基因聚类方法，用于将多个基因序列或小基因簇聚类到多个大基因簇中，其中上述基因簇是指一组有相同来源的基因集合，上述小基因簇是指基因数量小于预设值的基因集合，上述大基因簇是指基因数量大于上述预设值的基因集合；上述方法包括：

根据每个上述基因序列或小基因簇在多个样品中的丰度谱，以及每个上述大基因簇在上述多个样品中的丰度谱，通过多元线性回归，计算每个上述大基因簇与每个上述基因序列或小基因簇的相关性；

找到与每个上述大基因簇的相关性大于设定值的基因序列或小基因簇，并将所找到的基因序列或小基因簇合并至相关联的大基因簇中，组合成新的大基因簇。

进一步地，上述方法中的多元线性回归依据如下方程式进行：

其中，E_ij表示基因序列或小基因簇j在样品i中的相对丰度，a_ik表示大基因簇k在样品i中的相对丰度，e_kj表示基因序列或小基因簇j与大基因簇k的相关性。

进一步地，上述方法还包括使用系数收缩技术计算出上述相关性e_kj。

进一步地，上述小基因簇是相似丰度谱基因簇，上述大基因簇是宏基因组物种基因组，在本发明的一个实施例中，上述系数收缩技术是LASSO收缩方法，上述相关性的设定值为0.1。

根据本发明的第二方面，本发明提供一种基因聚类装置，用于将多个基因序列或小基因簇聚类到多个大基因簇中，其中上述基因簇是指一组有相同来源的基因集合，上述小基因簇是指基因数量小于预设值的基因集合，上述大基因簇是指基因数量大于上述预设值的基因集合；上述装置包括：

相关性计算单元，用于根据每个上述基因序列或小基因簇在多个样品中的丰度谱，以及每个上述大基因簇在上述多个样品中的丰度谱，通过多元线性回归，计算每个上述大基因簇与每个上述基因序列或小基因簇的相关性；

基因簇合并单元，用于找到与每个上述大基因簇的相关性大于设定值的基因序列或小基因簇，并将所找到的基因序列或小基因簇合并至相关联的大基因簇中，组合成新的大基因簇。

根据本发明的第三方面，本发明提供一种宏基因组组装方法，用于将来自宏基因组的多个基因序列或小基因簇聚类到多个大基因簇中，其中上述基因簇是指一组有相同来源的基因集合，上述小基因簇是指基因数量小于预设值的基因集合，上述大基因簇是指基因数量大于上述预设值的基因集合；上述方法包括：

根据每个上述基因序列或小基因簇在多个样品中的丰度谱，以及每个上述大基因簇在上述多个样品中的丰度谱，通过多元线性回归，计算每个上述大基因簇与每个上述基因序列或小基因簇的相关性；

找到与每个上述大基因簇的相关性大于设定值的基因序列或小基因簇，并将所找到的基因序列或小基因簇合并至相关联的大基因簇中，组合成新的大基因簇；

将测序读长序列比对至上述新的大基因簇，获取单个基因簇特有的测序读长序列，进行单独组装，获得单个基因组草图。

进一步地，上述方法中的多元线性回归依据如下方程式进行：

其中，E_ij表示基因序列或小基因簇j在样品i中的相对丰度，a_ik表示大基因簇k在样品i中的相对丰度，e_kj表示基因序列或小基因簇j与大基因簇k的相关性。

进一步地，上述方法还包括使用系数收缩技术计算出上述相关性e_kj。

进一步地，上述小基因簇是相似丰度谱基因簇，上述大基因簇是宏基因组物种基因组，在本发明的一个实施例中，上述系数收缩技术是LASSO收缩方法，上述相关性的设定值为0.1。

根据本发明的第四方面，本发明提供一种宏基因组组装装置，用于将来自宏基因组的多个基因序列或小基因簇聚类到多个大基因簇中，其中上述基因簇是指一组有相同来源的基因集合，上述小基因簇是指基因数量小于预设值的基因集合，上述大基因簇是指基因数量大于上述预设值的基因集合；上述装置包括：

相关性计算单元，用于根据每个上述基因序列或小基因簇在多个样品中的丰度谱，以及每个上述大基因簇在上述多个样品中的丰度谱，通过多元线性回归，计算每个上述大基因簇与每个上述基因序列或小基因簇的相关性；

基因簇合并单元，用于找到与每个上述大基因簇的相关性大于设定值的基因序列或小基因簇，并将所找到的基因序列或小基因簇合并至相关联的大基因簇中，组合成新的大基因簇；

基因组组装单元，用于将测序读长序列比对至上述新的大基因簇，获取单个基因簇特有的测序读长序列，进行单独组装，获得单个基因组草图。

本发明的基因聚类方法基于基因丰度谱的关联信息，利用多元线性回归方法，寻找与基因序列或小基因簇关联的大基因蔟，提升大基因簇集合，从而提升基因组的组装结果。本发明的基因聚类方法能够有效改进基因聚类效果，得到更有代表性的基因簇；且通过进一步组装获取更好的基因组草图，为深入的生物信息挖掘提供基础。

附图说明

图1为本发明实施例的基因聚类方法流程示意图；

图2为本发明实施例的基因聚类装置结构框图；

图3为本发明实施例的宏基因组组装方法流程示意图；

图4为本发明实施例的宏基因组组装装置结构框图；

图5为本发明实施例的CAG578以及9个Faecalibacterium MGS在396个样品中的丰度谱图；横坐标表示样品，样品顺序按照CAG578的丰度由高到低排列；纵坐标表示基因丰度；

图6为本发明实施例的宏基因组组装结果统计图，灰色柱表示原有研究结果，黑色柱表示运用本发明方法后的结果；纯色填充表示与参考基因组匹配的序列长度，非纯色填充表示不能与参考基因组匹配的序列长度。

具体实施方式

下面通过具体实施方式结合附图对本发明作进一步详细说明。

如图1所示，本发明实施例的基因聚类方法，用于将多个基因序列或小基因簇聚类到多个大基因簇中，其中基因簇是指一组有相同来源的基因集合，小基因簇是指基因数量小于预设值的基因集合，大基因簇是指基因数量大于预设值的基因集合。

基因序列是指单个基因序列或片段，基因采用其广泛性含义，包括任何能够表达蛋白质的核苷酸序列和任何其他序列片段。小基因簇是指两个以上基因序列的集合。基因簇是指一组有相同来源的基因集合，例如来源于同一细菌基因组的基因集合。每个基因簇中的基因一般具有相似的丰度谱。小基因簇与大基因簇的划分标准是基因集合中基因的数量多少，一般基因数量小于预设值的基因集合称为小基因簇，而基因数量大于预设值的基因集合称为大基因簇。预设值是人为设定的数值，根据具体需要和应用场景可以设定适当的数值，例如设置为500、600、700或800等。在本发明的一个实施例中，预设值为700，即基因的数量大于700的基因集合称为大基因簇，基因的数量小于700的基因集合称为小基因簇。

本发明实施例的基因聚类方法包括：

S110：根据每个基因序列或小基因簇在多个样品中的丰度谱，以及每个大基因簇在多个样品中的丰度谱，通过多元线性回归，计算每个大基因簇与每个基因序列或小基因簇的相关性。

基因序列、小基因簇或大基因簇在多个样品中的丰度谱，可以按照现有技术的方法得到，例如依据Nielsen H.B.等人,Nature Biotechnology,2014中公开的方法。具体通过如下步骤实现：(1)通过宏基因组组装、基因预测、基因序列去冗余获取多样品特有的非冗余基因序列集合；(2)将测序读长序列(Reads)比对到上一步获取的基因序列集合，计算每个基因在每个样品的丰度，得到基因在多样品中的丰度谱；(3)根据基因在多样品的丰度谱，计算基因与基因之间的相关系数，根据相关系数对基因进行聚类，将相关系数满足一定阈值的基因聚集成基因簇。得到的基因簇根据预设值(例如700)可以划分为小基因簇和大基因簇，同时根据基因簇中每个基因在每个样品的丰度就可以得到小基因簇或大基因簇在多个样品中的丰度谱。

本发明进行多元线性回归的理论依据是基因序列的丰度是其源头基因组丰度的线性加和。相应的，基因序列或小基因簇的丰度是其源头大基因簇丰度的线性加和。在一个典型但非限定性的实施例中，多元线性回归依据如下方程式进行：

其中，因变量E_ij表示基因序列或小基因簇j在样品i中的相对丰度，自变量a_ik表示大基因簇k在样品i中的相对丰度，相关系数e_kj表示基因序列或小基因簇j与大基因簇k的相关性，即表示因变量和自变量的相关性强弱。

在一个典型但非限定性的实施例中，小基因簇是相似丰度谱基因簇(Co-abundance Gene Group，CAG)，大基因簇是宏基因组物种(Metagenomics Species，MGS)基因组。

值得说明的是，当前的生物研究中，很多情况下都是样品数量小于未知系数数量，例如在本发明中有可能样品i数量小于未知系数e_kj数量。在这种情况下，需要引入LASSO(Least Absolute Shrinkage and Selection Operator，最小绝对收缩和选择算子)、子集回归或岭回归等技术进行系数收缩，计算出方程的解。

S120：找到与每个大基因簇的相关性大于设定值的基因序列或小基因簇，并将所找到的基因序列或小基因簇合并至相关联的大基因簇中，组合成新的大基因簇。

在该步骤中，相关性的设定值是人为设定的经验值，在不同应用场景下该值可能不同。在一个典型但非限定性的实施例中，在人肠道中存在的所有微生物的基因组作为人肠道宏基因组，相似丰度谱基因簇(Co-abundance Gene Group，CAG)代表小基因簇，宏基因组物种(Metagenomics Species，MGS)基因组代表大基因簇的情况下，相关性的设定值为0.1，即大基因簇MGS与小基因簇CAG之间的相关性大于0.1表示二者之间有强相关性，二者可以合并在一起；大基因簇MGS与小基因簇CAG之间的相关性小于0.1表示二者之间没有强相关性，二者不可以合并在一起。

本发明实施例的基因聚类方法，基于基因丰度谱的关联信息，利用多元线性回归方法，寻找与基因序列或小基因簇关联的大基因蔟，提升大基因簇集合，从而提升基因组的组装结果。此外，本发明实施例的基因聚类方法，在多个领域可以有应用，如RNA测序中基因表达量的计算，宏转录组分析也是本发明的一个潜在应用领域。可以运用本发明的方法解析基因表达量的组成。

对应于本发明实施例的基因聚类方法，本发明还提供一种基因聚类装置，如图2所示，包括：相关性计算单元210，用于根据每个基因序列或小基因簇在多个样品中的丰度谱，以及每个大基因簇在多个样品中的丰度谱，通过多元线性回归，计算每个大基因簇与每个基因序列或小基因簇的相关性；基因簇合并单元220，用于找到与每个大基因簇的相关性大于设定值的基因序列或小基因簇，并将所找到的基因序列或小基因簇合并至相关联的大基因簇中，组合成新的大基因簇。

如图3所示，本发明实施例的宏基因组组装方法，用于将来自宏基因组的多个基因序列或小基因簇聚类到多个大基因簇中，其中基因簇是指一组有相同来源的基因集合，小基因簇是指基因数量小于预设值的基因集合，大基因簇是指基因数量大于预设值的基因集合。

本发明实施例的宏基因组组装方法包括：

S310：根据每个基因序列或小基因簇在多个样品中的丰度谱，以及每个大基因簇在多个样品中的丰度谱，通过多元线性回归，计算每个大基因簇与每个基因序列或小基因簇的相关性。

S320：找到与每个大基因簇的相关性大于设定值的基因序列或小基因簇，并将所找到的基因序列或小基因簇合并至相关联的大基因簇中，组合成新的大基因簇。

S330：将测序读长序列比对至新的大基因簇，获取单个基因簇特有的测序读长序列，进行单独组装，获得单个基因组草图。

对应于本发明实施例的宏基因组组装方法，本发明还提供一种宏基因组组装装置，如图4所示，包括：相关性计算单元410，用于根据每个基因序列或小基因簇在多个样品中的丰度谱，以及每个大基因簇在多个样品中的丰度谱，通过多元线性回归，计算每个大基因簇与每个基因序列或小基因簇的相关性；基因簇合并单元420，用于找到与每个大基因簇的相关性大于设定值的基因序列或小基因簇，并将所找到的基因序列或小基因簇合并至相关联的大基因簇中，组合成新的大基因簇；基因组组装单元430，用于将测序读长序列比对至新的大基因簇，获取单个基因簇特有的测序读长序列，进行单独组装，获得单个基因组草图。

本发明的宏基因组组装方法能够有效改进基因聚类效果，得到更有代表性的基因簇；且通过进一步组装获取更好的基因组草图，为深入的生物信息挖掘提供基础。

本领域技术人员可以理解，上述实施方式中各种方法的全部或部分步骤可以通过程序来指令相关硬件完成，该程序可以存储于计算机可读存储介质中，存储介质可以包括：只读存储器、随机存储器、磁盘或光盘等。

以下通过实施例详细说明本发明的技术方案和效果，应当理解，实施例仅是示例性的，不能理解为对本发明保护范围的限制。

实施例

本实施例选用文章Nielsen H.B.等人,Nature Biotechnology,2014中的数据，该数据可以从公开网站http://www.cbs.dtu.dk/databases/CAG/获取。网站上可以下载基因簇的丰度谱文件(名称为CAG abundance profiles)，该文件包含基因簇在396个样品中的丰度谱信息。其中，CAG指Co-abundance Gene Group(相似丰度谱基因簇)，表示多个具有相似丰度谱的基因聚集成的基因簇，并且基因数量大于700的基因簇被定义为MGS，指Metagenomics Species(宏基因组物种)，表示宏基因组物种，具体方法在文章中有详细描述。这个数据集合包含6,640个CAG和741个MGS基因簇。

基于CAG和MGS的丰度谱文件，进行如下分析：

1)依据公式

计算741个MGS对应于每个CAG的回归系数e_kj。使用R语言中的软件包，调用LASSO回归分析功能，并且选取截距为0、回归系数非负的设置，其中R语言中的软件包的来源如下：

https://cran.r-project.org/web/packages/penalized/index.html。

2)对每个CAG，只保留回归系数大于0.1的MGS；

3)统计每个MGS对应的所有关联的CAG；

4)将MGS和所有关联CAG合并在一起，组合成新的MGS基因簇；

5)将单个样品的测序Reads比对到新的MGS基因簇，之后利用这些Reads进行基因组组装。比对软件可以选用SOAPaligner(http://soap.genomics.org.cn/ soapaligner.html)。

6)将组装得到的重叠群(Contig)序列(表示测序Reads经过拼接后得到的较长的DNA序列)与参考基因组比较，统计比较结果。组装软件可以选用Velvet-sc(http:// bix.ucsd.edu/projects/singlecell/)。

为了展示方便，该实施例选取细菌Faecalibacterium(栖粪杆菌属)对应的MGS进行展示。栖粪杆菌属是人粪便中广泛存在的一类细菌，普遍存在于人的肠道中。采用的数据中，栖粪杆菌属包含9个MGS，编号分别为：MGS82、MGS97、MGS106、MGS125、MGS128、MGS133、MGS142、MGS195和MGS206，这些MGS可以理解成9个不同的栖粪杆菌属基因组。

经过上述分析，找到与这9个MGS关联的CAG，结果统计参见表一。在基因数量上，这9个MGS有29.18％～51.58％的提升。更有意思的是，这9个MGS都与CAG578关联；然而，这9个MGS与CAG578在丰度谱上的关联性并不强，如图5所示。这种弱相关性是已有研究方法所不能解决的难点，本发明可以有效解决这一问题，同时找到与CAG578相关联的多个MGS。CAG578包含的基因具有编码核糖体蛋白、代谢相关等重要功能，是基因组的核心组成。进一步的组装结果显示，这9个MGS的基因组在总长度上平均有15.81％的提升，并且与参考基因组匹配的序列长度上平均有13.50％的提升(如表二和图6)。表明本发明能有效提高基因组的组装结果。

表一：9个栖粪杆菌属的MGS以及关联的CAG统计

表二：9个MGS改进前后的组装结果统计

*参考序列对应和NCBI数据库中四个栖粪杆菌属的参考基因组。A2_165表示Faecalibacterium prausnitzii A2-165，KLE12555表示Faecalibacteriumcf.prausnitzii KLE1255，SL3_3表示Faecalibacterium prausnitzii SL3/3，L2_6表示Faecalibacterium prausnitzii L2-6。

以上内容是结合具体的实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。

Claims (8)

1.一种基因聚类方法，其特征在于，所述方法用于将多个基因序列或小基因簇聚类到多个大基因簇中，其中所述基因簇是指一组有相同来源的基因集合，所述小基因簇是指基因数量小于预设值的基因集合，所述大基因簇是指基因数量大于所述预设值的基因集合；所述方法包括：

根据每个所述基因序列或小基因簇在多个样品中的丰度谱，以及每个所述大基因簇在所述多个样品中的丰度谱，通过多元线性回归，计算每个所述大基因簇与每个所述基因序列或小基因簇的相关性；

找到与每个所述大基因簇的相关性大于设定值的基因序列或小基因簇，并将所找到的基因序列或小基因簇合并至相关联的大基因簇中，组合成新的大基因簇；

所述多元线性回归依据如下方程式进行：

其中，E_ij表示基因序列或小基因簇j在样品i中的相对丰度，a_ik表示大基因簇k在样品i中的相对丰度，e_kj表示基因序列或小基因簇j与大基因簇k的相关性。

2.根据权利要求1所述的基因聚类方法，其特征在于，所述方法还包括使用系数收缩技术计算出所述相关性e_kj。

3.根据权利要求2所述的基因聚类方法，其特征在于，所述小基因簇是相似丰度谱基因簇，所述大基因簇是宏基因组物种基因组，所述系数收缩技术是LASSO收缩方法，所述相关性的设定值为0.1。

4.一种基因聚类装置，其特征在于，所述装置用于将多个基因序列或小基因簇聚类到多个大基因簇中，其中所述基因簇是指一组有相同来源的基因集合，所述小基因簇是指基因数量小于预设值的基因集合，所述大基因簇是指基因数量大于所述预设值的基因集合；所述装置包括：

相关性计算单元，用于根据每个所述基因序列或小基因簇在多个样品中的丰度谱，以及每个所述大基因簇在所述多个样品中的丰度谱，通过多元线性回归，计算每个所述大基因簇与每个所述基因序列或小基因簇的相关性；

基因簇合并单元，用于找到与每个所述大基因簇的相关性大于设定值的基因序列或小基因簇，并将所找到的基因序列或小基因簇合并至相关联的大基因簇中，组合成新的大基因簇；

所述多元线性回归依据如下方程式进行：

其中，E_ij表示基因序列或小基因簇j在样品i中的相对丰度，a_ik表示大基因簇k在样品i中的相对丰度，e_kj表示基因序列或小基因簇j与大基因簇k的相关性。

5.一种宏基因组组装方法，其特征在于，所述方法用于将来自宏基因组的多个基因序列或小基因簇聚类到多个大基因簇中，其中所述基因簇是指一组有相同来源的基因集合，所述小基因簇是指基因数量小于预设值的基因集合，所述大基因簇是指基因数量大于所述预设值的基因集合；所述方法包括：

根据每个所述基因序列或小基因簇在多个样品中的丰度谱，以及每个所述大基因簇在所述多个样品中的丰度谱，通过多元线性回归，计算每个所述大基因簇与每个所述基因序列或小基因簇的相关性；

找到与每个所述大基因簇的相关性大于设定值的基因序列或小基因簇，并将所找到的基因序列或小基因簇合并至相关联的大基因簇中，组合成新的大基因簇；

将测序读长序列比对至所述新的大基因簇，获取单个基因簇特有的测序读长序列，进行单独组装，获得单个基因组草图；

所述多元线性回归依据如下方程式进行：

其中，E_ij表示基因序列或小基因簇j在样品i中的相对丰度，a_ik表示大基因簇k在样品i中的相对丰度，e_kj表示基因序列或小基因簇j与大基因簇k的相关性。

6.根据权利要求5所述的宏基因组组装方法，其特征在于，所述方法还包括使用系数收缩技术计算出所述相关性e_kj。

7.根据权利要求6所述的宏基因组组装方法，其特征在于，所述小基因簇是相似丰度谱基因簇，所述大基因簇是宏基因组物种基因组，所述系数收缩技术是LASSO收缩方法，所述相关性的设定值为0.1。

8.一种宏基因组组装装置，其特征在于，所述装置用于将来自宏基因组的多个基因序列或小基因簇聚类到多个大基因簇中，其中所述基因簇是指一组有相同来源的基因集合，所述小基因簇是指基因数量小于预设值的基因集合，所述大基因簇是指基因数量大于所述预设值的基因集合；所述装置包括：

相关性计算单元，用于根据每个所述基因序列或小基因簇在多个样品中的丰度谱，以及每个所述大基因簇在所述多个样品中的丰度谱，通过多元线性回归，计算每个所述大基因簇与每个所述基因序列或小基因簇的相关性；

基因簇合并单元，用于找到与每个所述大基因簇的相关性大于设定值的基因序列或小基因簇，并将所找到的基因序列或小基因簇合并至相关联的大基因簇中，组合成新的大基因簇；

基因组组装单元，用于将测序读长序列比对至所述新的大基因簇，获取单个基因簇特有的测序读长序列，进行单独组装，获得单个基因组草图；

所述多元线性回归依据如下方程式进行：

其中，E_ij表示基因序列或小基因簇j在样品i中的相对丰度，a_ik表示大基因簇k在样品i中的相对丰度，e_kj表示基因序列或小基因簇j与大基因簇k的相关性。

Images