一种基于二代和三代测序技术的宏病毒组分析方法

摘要

本发明公开了一种基于二代和三代测序技术的宏病毒组分析方法，包括对二代测序原始数据进行质量控制和过滤，得到clean reads，对质控后的数据进行单样本和多样本组装得到contigs序列；对三代测序原始数据进行校正和质量控制，得到clean长序列，对质控后的三代数据进行组装得到contigs序列；对二三代测序的质控后数据进行混合组装得到contigs序列；将所有contig合并构建非冗余的contigs集合；最后进行病毒的识别和判定、病毒的物种注释和功能注释。本发明提出了可靠的基于二代和三代测序技术的宏病毒组的分析方法，其实现方法简单，应用范围广泛。

说明书

技术领域

本发明涉及二代和三代测序技术领域以及宏病毒组学领域，具体为一种基于二代和三代测序技术的宏病毒组分析方法。

背景技术

病毒是微生物的重要组成部分，其组成变化与众多疾病相关，是地球生物圈最丰富的生物种类。可以侵袭细菌且导致宿主裂解的病毒又被称作噬菌体，其作为一种可移动遗传元件，可将遗传物质在细菌之间传播。近年来，由于抗生素的广泛使用，超级细菌层出不穷，病原菌对抗生素耐药性的威胁越来越大，抗生素治疗效果受到严峻的挑战，由于噬菌体可以入侵并改变(通常是降低)宿主菌的适应性，噬菌体治疗又重新成为热点。因此，病毒在调节细菌的生物量、维持生物多样性、基因水平转移以及整个生物圈内的生物化学循环都有着重要作用，还可作抗生素的替代品，具有很强的实用价值。

宏病毒组学是在宏基因组学理论的基础上，结合现有的病毒分子生物学检测技术而兴起的一个新的学科分支，是某类样本中所有病毒或病毒类似物及其所携带遗传信息的总称。宏病毒组直接以环境中所有病毒的遗传物质为研究对象，能够快速准确的鉴定出环境中所有的病毒组成，在病毒发现、病毒溯源、微生物预警等研究方面具有重要作用。宏病毒研究可应用于人或动物肠道或者血液样本、海洋、土壤等的研究，用以挖掘潜在的对人类和环境的危害。测序技术的不断进步和分析方法的改进使得人们得以一窥未知的科学领域。此外，病毒的很多特性是基于整个群落环境及个体间的相互影响的。因此，宏病毒学相比于做单个个体的研究更能发现其特性。

由于测序技术原理的有限性和病毒序列的复杂性，二代测序原始序列读长较短(通常只有50～250bp)，难以获得低丰度的病毒，且组装片段准确率受多因素制约；三代测序读长较长(平均10-15kb，最长可达80kb)，但因为成本问题，一直不能广泛应用。如何兼顾测序成本和测序通量，有效处理原始测序数据进而得到完整准确的病毒序列是生物信息领域一直关注的热点问题。随着宏病毒组学的发展，针对病毒测序数据的分析工具层出不穷，但能同时结合二代和三代测序技术的优点完成病毒序列的识别鉴定、物种分析和功能分析的工具并不多见，因此整合现有的这些工具，完成上述分析模块以及使分析尽量高效准确已经成为迫切需求。

发明内容

本发明的目的在于提供一种基于二代和三代测序技术的宏病毒组分析方法，本发明通过从二代和三代测序宏病毒组数据分别处理到二三代混合组装、病毒判定和识别与物种、功能注释的全部流程，为研究人员提供准确的分析结果，也让非生物信息专业的科研人员独立完成病毒高通量数据分析，达到优化科研人员的工作效率，降低科研成本的目的。本发明提出了可靠的基于二代和三代测序技术的宏病毒组的分析方法，其实现方法简单，应用范围广泛；以解决现有技术中存在的兼顾测序成本和测序通量的技术问题。

根据本发明的目的，提供了一种基于二代和三代测序技术的宏病毒组分析方法，包括以下步骤：

(1)将二代基因组测序所得的原始高通量测序数据进行过滤，得到过滤后的数据；将该过滤后的数据进行组装，得到二代组装contigs；

(2)将三代基因组测序所得的原始高通量测序数据进行过滤，得到过滤后的数据；将该过滤后的数据进行组装，得到三代组装contigs；该步骤中的组装根据样本测序数据量的大小进行，单样本测序数据量超过100G时，使用以下方法1或方法2的单一组装方法；若小于等于100G，则同时使用以下两种组装方法；其中：

方法1：通过reads长度、测序深度和基因组大小运用de Bruijn算法计算出性能最佳的K-mer，将其作为固定K-mer值检测重复序列并组装成contigs序列；

方法2：采用少数服从多数原则对reads进行纠错修正，根据重叠度确定碱基一致的区域为高质量区域，修整不能完全重叠的低质量区域；重新计算重叠区的序列一致性，序列完全重叠区域为unitigs，重叠区域包含部分重叠的片段为contigs；然后对unitigs序列文件使用k-medoids聚类算法进行迭代binning，得到bins；再进行序列判定：即当contig的所有unitigs都属于同一bin，那么保留该contig，反之，则保留所有unitigs；

(3)二代和三代宏病毒组测序数据进行混合组装：同时采用以下两种组装方法，得到二代和三代混合组装的contigs；其中：

方法1：构建步骤(1)中过滤后的数据的拼接集，得到二代拼接集，利用贝叶斯模型方法将步骤(2)中过滤后的数据来填补二代拼接集中间不能相连的片段，得到contigs；

方法2：使用HYBRIDSPADES算法对步骤(1)中过滤后的数据和步骤(2)中过滤后的数据进行组装得到contigs；

(4)构建非冗余的contigs集合：将步骤(1)得到的二代组装contigs、步骤(2)得到的三代组装contigs以及步骤(3)得到的二代和三代混合组装的contigs合并，根据序列的相似度将大于95％的序列聚类成一个cluster，选取最长的序列作为该cluster的代表序列，所有的代表序列即为获得的非冗余contigs集合；

(5)病毒序列的识别和判定：对步骤(4)中得到的非冗余的congtigs集合采用不少于三种方法进行病毒识别和判定，至少被两种方法识别为病毒的contigs被保留；将保留的contigs与病毒组数据库进行比对来判定完整性，过滤掉大于10％的错误匹配，最终保留的contigs被判定为病毒序列；病毒识别的方法包括：

方法1：应用机器学习使用基因组特征进行病毒识别；所述基因组特征包括结构、功能、分类注释和病毒标志基因；

方法2：基于病毒序列K-mers词频和机器学习识别病毒序列；

方法3：基于深度学习算法，设计出双通道卷积神经网络的网络结构，有效地对DNA片段进行病毒识别。

优选地，步骤(1)中将二代基因组测序所得的原始高通量测序数据进行过滤具体为：对二代测序宏病毒组的原始高通量数据进行过滤，去除adapter序列，去除N的比例大于10％的reads；数据经过过滤后通过分析碱基的组成及质量分布，进行碱基质量分析，去除质量值Q≤30的低质量reads；质量控制后的reads基于序列比对算法，去除碱基数超过85％比对到细菌污染物基因组的reads，最终得到高质量的clean reads。

优选地，步骤(1)中的组装具体为：基于K-mer迭代的de Bruijn图的组装算法得到组装结果contig，过滤掉长度在1.5kb以下的短序列；根据样本数和测序量，选择不同的数据组装方法；其中：

方法1：当单样本测序量小于等于20G且样本数少于100，必须同时进行单样本组装和多样本组装；单样本组装即单一样本组装，多样本组装即多个样本reads混合组装；

方法2：当单样本测序量大于20G，选择单样本组装方法。

优选地，步骤(2)中将三代基因组测序所得的原始高通量测序数据进行过滤具体为：对三代测序宏病毒组的原始高通量数据进行过滤，去除adapter序列，去除少于2次环绕次数的subreads；基于环状一致性测序CCS模式，来源于同一条模板链的subreads经过一致性校正，得到高精度单分子一致性的reads；质量控制后的reads基于序列比对算法，去除碱基数超过85％比对到细菌污染物基因组的reads，最终得到高质量的clean reads。

优选地，步骤(5)之后还包括病毒的物种注释，采用步骤(5)得到的最终保留的contigs，计算平均核苷酸一致性ANI，保留ANI>95％的病毒组数据集，即genus水平；通过将氨基酸水平的基因与TrEMBL数据库中的病毒子集进行比对，得到family水平的注释；基于K-mer特征通过监督机器学习方法对病毒基因组进行genus水平的物种注释；依据病毒的family和genus水平的注释结果，从已知的taxonomy库补充病毒序列的其他分类水平注释。

优选地，步骤(5)之后还包括病毒的功能注释，通过卷积神经网络检测已知的温和和毒性病毒序列的特异性特征，识别步骤(5)得到的最终保留的contigs是否含有这些特征，最终将病毒序列识别为温和和毒性病毒。

优选地，将保留的contigs与病毒组数据库进行比对来判定完整性，0～5％错误匹配的被认为高质量contig，5～10％错误匹配的是中等质量contig，大于10％的错误匹配是低质量contig需要被过滤。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，主要具备以下的技术优点：

(1)本发明提供了从二代和三代测序宏病毒组数据分别处理到二三代混合组装、病毒判定和识别与物种、功能注释的全部流程，为研究人员提供准确的分析结果，也让非生物信息专业的科研人员独立完成病毒高通量数据分析，达到优化科研人员的工作效率，降低科研成本的目的。本发明提出了可靠的基于二代和三代测序技术的宏病毒组的分析方法，其实现方法简单，应用范围广泛。

(2)本发明分析流程更加合理、全面，分析结果更加准确，解决了目前多代测序技术相结合的病毒组组装和鉴定流程不规范的问题，为研究人员提供便利和技术支撑。

附图说明

图1是本发明流程图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

本发明的一些实例的具体过程如下：

二代测序数据过滤：对二代测序宏病毒组的原始高通量数据进行过滤，去除adapter序列，去除N的比例大于10％的reads；基于碱基的组成及质量分布进行碱基质量分析，去除质量值Q≤30的低质量reads；质量控制后的reads基于短序列比对算法，去除碱基数超过85％比对到细菌等污染物基因组的reads，最终得到高质量的clean reads；

具体地，使用fastqc首先对原始测序数据进行质量评估，根据评估结果使用软件Trimmomatic和软件bowtie2进行质量过滤和DNA过滤。

二代数据组装：基于K-mer迭代的de Bruijn图的组装算法得到组装结果contig，过滤掉长度在1.5kb以下的短序列；根据样本数和测序量，选择不同的数据组装策略。其中：

策略1：当单样本测序量小于20G且样本数少于100，必须同时进行单样本组装和多样本组装。单样本组装即单一样本组装，多样本组装即多个样本reads混合组装；

策略2：当单样本测序量大于20G，可以选择单样本组装策略；

具体地，使用软件IDBA_UD完成单样本组装，使用软件megahit完成多样本组装。

三代测序数据过滤：对三代测序宏病毒组的原始高通量数据进行过滤，去除adapter序列，去除少于2次环绕次数的subreads；基于环状一致性测序CCS模式，来源于同一条模板链的subreads经过一致性校正，得到高精度单分子一致性的reads；质量控制后的reads基于序列比对算法，去除碱基数超过85％比对到细菌等污染物基因组的reads，最终得到高质量的clean reads；

具体地，使用CCS软件完成质量过滤，使用软件bowtie2进行DNA过滤。

三代测序组装：根据样本测序数据量的大小，单样本测序数据量超过100G时，可以使用单一组装策略；若不足100G，则同时使用两种组装策略；其中：

策略1：通过reads长度、测序深度和预估基因组大小计算出性能最佳的K-mer，将其作为固定K-mer值检测重复序列并组装成contigs序列，具体的使用软件flye实现；

策略2：采用少数服从多数原则对reads进行纠错修正，根据重叠度确定碱基一致的区域为高质量区域，修整不能完全重叠的低质量区域；重新计算重叠区的序列一致性，组装得到片段较短但碱基更准确的unitigs和片段更长但质量相对差一点的contigs，具体的使用软件CANU实现；对unitigs序列文件使用k-medoids聚类算法进行迭代binning，得到bins，具体的使用MetaBAT实现；序列判定：当contig的所有unitigs都属于同一bin，那么保留该contig，反之，则保留所有unitigs，具体的使用R、python或者shell语言实现。

二代和三代宏病毒组测序数据组装：采用两种组装策略。其中：

策略1：构建短读长(二代测序数据过滤后的reads)的拼接集，利用贝叶斯模型方法将低覆盖率的长读长(三代测序数据过滤后的reads)来改善二代拼接集中间不能相连的片段，得到contigs，具体的使用软件OPERA-MS实现；

策略2：使用汇编短读长和长读长的HYBRIDSPADES算法进行组装得到contigs，具体的使用软件metaSPAdes实现。

构建非冗余的contigs集合：将二代组装的contigs、三代组装的contigs和二三代混合组装的contigs合并，根据序列的相似度将大于95％的序列聚类成一个cluster，优选的，选取最长的序列作为该cluster的代表序列，所有的代表序列即为获得的非冗余contigs集合，具体的使用软件CD-hit实现。

病毒的识别和判定：对非冗余的congtigs采用不少于三种方法进行病毒识别和判定，至少被两种方法识别为病毒的contigs被保留，具体的使用R、python或者shell语言实现；将保留的contigs与公共完整病毒组数据库进行比对来估计完整性，0～5％错误匹配的被认为高质量contig，5～10％错误匹配的是中等质量contig，大于10％的错误匹配是低质量contig需要被过滤，最终保留的contigs被判定为病毒序列，具体的使用软件CheckV实现。该病毒集合并用于后续的物种注释和功能注释。所述的病毒识别方法包括：

策略1：应用机器学习使用基因组特征(包括结构/功能/分类注释和病毒标志基因)来进行病毒识别，具体的使用软件VirSorter实现；

策略2：基于病毒序列K-mers词频和机器学习识别病毒序列，具体的使用软件VirFinder实现；

策略3：基于深度学习算法，设计出"双通道卷积神经网络"的网络结构，有效地对DNA片段进行病毒识别，具体的使用软件PPP-META实现。

病毒的物种注释：全部病毒contigs计算平均核苷酸一致性(ANI)，保留ANI>95％的病毒组数据集，即genus水平，具体的使用软件FastANI实现；通过将氨基酸水平的基因与TrEMBL数据库中的病毒子集进行比对，得到family水平的注释，具体的使用软件demovir实现；基于K-mer特征通过监督机器学习方法对病毒基因组进行genus水平的物种注释，具体的使用软件VirusTaxo实现；依据病毒的family和genus水平的注释结果，从已知的taxonomy库完善病毒序列的其他分类水平注释，具体的使用R或者python语言实现。

病毒的功能注释：通过卷积神经网络检测病毒序列特征以获取有价值的局部特征，将病毒序列识别为温和和毒性病毒，具体的使用DeePhage实现。

综上所述，本发明构建了更加全面和更为准确的病毒鉴定识别方法。本发明开发同时基于二代测序、三代测序和二三代测序的病毒组组装方法，从而解决多代测序的宏病毒组组装鉴定流程不规范的问题，同时整合多种组装和鉴定工具，分析结果更为准确。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。