CN114121160B - 一种检测样本中宏病毒组的方法和系统 - Google Patents
一种检测样本中宏病毒组的方法和系统 Download PDFInfo
- Publication number
- CN114121160B CN114121160B CN202111407911.2A CN202111407911A CN114121160B CN 114121160 B CN114121160 B CN 114121160B CN 202111407911 A CN202111407911 A CN 202111407911A CN 114121160 B CN114121160 B CN 114121160B
- Authority
- CN
- China
- Prior art keywords
- virus
- contigs
- sequencing data
- sequences
- nanopore
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
- G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
- G16B30/10—Sequence alignment; Homology search
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
- G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
- G16B30/20—Sequence assembly
Abstract
本发明公开了一种检测样本中宏病毒组的方法，属于宏基因组分析技术领域，包括将待检测样本的二代测序数据和三代测序数据进行相互校正和混合组装，得到混合组装contigs，并将三代测序数据单独组装得到Nanopore contigs；进一步将所述混合组装contigs和Nanopore contigs分别进行比对注释，得到候选病毒contigs、非病毒contigs，和Nanopore病毒contigs；最后将3种数据集进行聚类分析，根据聚类结果进一步地综合分析，补充遗漏的病毒序列，并校正物种注释结果，得到更灵敏、准确且全面的病毒鉴定结果。
Description
技术领域
本发明属于宏基因组分析技术领域，具体地，涉及一种检测样本中宏病毒组的方法和系统。
背景技术
随着宏基因组学的发展，越来越多的研究证明，病毒在不同生态系统中发挥着关键作用，因此对宏基因组数据进行病毒的分析是非常有必要的。
近年来，以高通量测序为基础的病毒宏基因组学技术以其时效性、高通量的优势，使得人们能够对不同类型样本进行微生物测序，对大量的不可培养的病毒进行研究。目前宏基因组测序方法主要包括二代测序(Next-generation Sequencing)和三代测序技术。二代测序因通量高、准确性高而广泛应用于病毒宏基因组领域；而三代测序技术中的纳米孔单分子测序平台(Oxford Nanopore Technologies，ONT)因速度快、测序读长长等优势，也越来越多应用于病毒宏基因组领域。然而，现有的病毒宏基因组学分析方法，存在如下不足：
(1)现有的病毒序列数据库不全，大量的新病毒未被测序，造成病毒序列的鉴定敏感性和准确性均较低；
(2)基于二代测序技术产生的contigs较短，鉴定的病毒序列完整度低；
(3)基于Nanopore测序平台产生的reads准确性较低，鉴定的病毒序列质量低。
发明内容
为了解决上述技术问题中的至少一个，本发明采取的技术方案如下：
本发明第一方面提供一种检测样本中宏病毒组的方法，包括以下步骤：
S1，分别获得待检测样本二代测序数据和三代测序数据；
S2，将所述二代测序数据和所述三代测序数据进行相互校正和混合组装，得到混合组装contigs，将三代测序数据单独组装得到Nanopore contigs；
S3，病毒鉴定与物种注释：
将所述混合组装contigs进行比对注释，得到候选病毒contigs和非病毒contigs，并对鉴定的病毒contigs进行物种注释，得到候选病毒的contigs的注释结果，
将所述Nanopore contigs进行比对注释，得到Nanopore病毒contigs，同样进行物种注释得到Nanopore病毒contigs的注释结果；
S4，病毒鉴定的校正：
S41，将候选病毒contigs、非病毒contigs和Nanopore病毒contigs三个数据集进行聚类分析，
S42，筛选出代表序列来自Nanopore病毒contigs的聚类单元，该聚类单元中如果某个非病毒contig能够比对上所述代表序列，则该非病毒contig作为遗漏的序列补充进候选病毒contigs，得到校正后的候选病毒contigs，将该聚类单元里的校正后的候选病毒contigs作为一个类病毒株，并使用代表序列的物种注释结果作为该类病毒株的注释结果，
S43，将不在筛选出的聚类单元里的每条候选病毒contig均作为类病毒株，与经过校正得到的类病毒株的注释结果进行汇总，得到最终的病毒鉴定结果。
在本发明的一些实施方案中，二代测序数据利用二代测序平台测序得到的数据，所述第二代测序平台包括但不限于Illumina-Solexa(ATM，HiSeq2000TM等)、ABI-Solid和Roche-454(焦磷酸测序)；三代测序数据利用单分子测序平台测序得到的数据所述单分子测序平台包括但不限于Helicos公司的真实单分子测序技术(True Single MoleculeDNAsequencing)，Pacific Biosciences公司单分子实时测序(single molecule real-time(SMRTTM))，以及Oxford Nanopore Technologies公司的纳米孔测序技术等(Rusk,Nicole(2009-04-01).Cheap Third-Generation Sequencing.Nature Methods 6(4):2446(4)。
在本发明的一些实施方案中，用于二代测序和用于三代测序的核酸样本利用同样的方法抽提或者利用不同的方法抽提，在一个优选的实施方案中，均采用酚-氯仿抽提法提取核酸样本。
进一步地，获得的所述二代测序数据和所述三代测序数据为经过质控和去除宿主序列后的二代测序数据和三代测序数据。
在本发明的一些实施方案中，对于二代测序数据，进行质控的步骤如下：
(1)过滤含有2个以上N碱基的read，对保留下来的序列，在序列首尾两端使得滑动窗口检测碱基质量，切除质量值低于Q15的短片段序列，同时检测序列尾端的polyX序列并去除；
(2)过滤低复杂度序列、重复序列和接头序列；
(3)过滤长度小于15bp的序列，
在本发明中，对于二代测序数据进行质控，可以采取任意能够实现上述功能的软件、算法或程序，例如fastp v0.20，参数设置-n 2 -q 15 -x -5 -3 -y –dedup --dup_calc_accuracy 6。
在本发明的一些实施方案中，对于三代测序数据，进行质控的步骤为：只保留质量值大于7的序列。
在本发明的一些实施方案中，对于二代测序数据和三代测序数据去除宿主序列，均采用首先与宿主参考基因组进行比较，提取未比对上的序列，即得到去除宿主后的二代测序数据或三代测序数据。在本发明的一些优选的实施方案中，利用不同的软件、算法或程序分别将二代测序数据或三代测序数据与宿主参考基因组，并用同样的软件、算法或程序提取未比对上的序列，在本发明的一些具体实施方案中，利用bwa v0.7.17将二代测序数据与宿主参考基因组进行比对，利用minimap2 v2.2将三代测序数据与宿主参考基因组进行比对，之后均利用samtools v1.9提取未比对上的序列。
在本发明的一些具体实施方案中，利用OPERA-MS v0.9.0软件将所述二代测序数据和三代测序数据进行相互校正与混合组装，得到混合组装contigs。在本发明中，也可以利用任意具有相同功能的软件、算法或程序进行相互校正与混合组装。
在本发明的一些具体实施方案中，利用flye v2.8.3对所述三代测序数据进行组装，得到Nanopore contigs。在本发明中，也可以利用任意具有相同功能的软件、算法或程序对三代测序数据进行组装。
进一步地，根据权利要求1所述的一种检测样本中宏病毒组的方法，其特征在于，步骤S2中，混合组装contigs和Nanopore contigs只保留长度大于或等于500bp的序列。
在本发明的一些实施方案中，步骤S3中，所述比对注释是指将混合组装contigs或Nanopore contigs与病毒特异性蛋白序列数据库比对。
在本发明的一些具体实施方案中，利用CheckV v0.7软件对所述混合组装contigs和所述Nanopore contigs进行病毒序列的比对注释，得到病毒注释结果。在本发明中，也可以利用任意具有相同功能的软件、算法或程序对所述混合组装contigs和所述Nanoporecontigs进行病毒序列的比对注释。
在本发明的一些实施方案中，所述的CheckV软件配套的数据库来自https://portal.nersc.gov/CheckV/checkv-db-v1.0.tar.gz。
在本发明的一些实施方案中，步骤S3中，过滤掉比对结果中未鉴定到病毒基因的序列，并过滤掉注释上的宿主基因数量比病毒基因数量大于10倍的序列，得到候选病毒contigs、非病毒contigs和Nanopore病毒contigs。
在本发明的一些实施方案中，步骤S41中，利用PSI-CD-HIT脚本进行聚类分析，优选的参数为-c 0.95 -aS 0.85 -para 4 -blp 20。在本发明中，也可以利用任意具有相同功能的软件、算法或程序进行所述聚类分析。
在本发明的一些实施方案中，步骤S42中，筛选出的聚类单元中，还满足以下条件：
(1)序列数大于3；
(2)序列来自所述三个数据集中的至少2个。
在本发明的一些实施方案中，步骤S42中，所述能够比对上是指比对结果evlaue<1e-8，且输入序列的coverage>80％。
在本发明的一些实施方案中，步骤S42中，将聚类单元中的代表序列作为目标序列，校正后的病毒序列集作为输入序列，进行共线性分析，对序列进行排序，并记录位置坐标及序列方向信息，方便比较类病毒株内的序列间的关系。
本发明第二方面提供一种检测样本中宏病毒组的系统，能够实现本发明第一方面的方法，包括：
测序数据存储模块，用于获得待检测样本二代测序数据和三代测序数据并存储；
组装模块，与所述测序数据存储模块连接，用于将所述二代测序数据和所述三代测序数据进行相互校正和混合组装，得到混合组装contigs，并用于将三代测序数据单独组装得到Nanopore contigs；
病毒鉴定与注释模块，与所述组装模块连接，用于将所述混合组装contigs进行比对注释，得到候选病毒contigs和非病毒contigs，并对鉴定的病毒contigs进行物种注释，得到候选病毒的contigs的注释结果；并用于将所述Nanopore contigs进行比对注释，得到Nanopore病毒contigs，同样进行物种注释得到Nanopore病毒contigs的注释结果；
病毒鉴定校正模块，与所述病毒鉴定与注释模块连接，用于：
将候选病毒contigs、非病毒contigs和Nanopore病毒contigs三个数据集进行聚类分析，筛选出代表序列来自Nanopore病毒contigs的聚类单元，该聚类单元中如果某个非病毒contig能够比对上所述代表序列，则该非病毒contig作为遗漏的序列补充进候选病毒contigs，得到校正后的候选病毒contigs，将该聚类单元里的校正后的候选病毒contigs作为一个类病毒株，并使用代表序列的物种注释结果作为该类病毒株的注释结果；将不在筛选出的聚类单元里的每条候选病毒contig均作为类病毒株，与经过校正得到的类病毒株的注释结果进行汇总。
在本发明的一些实施方案中，进一步包括测序数据前处理模块，用于对二代测序数据或三代测序数据进行质控和和去除宿主序列。
本发明的有益效果
相对于现有技术，本发明具有以下有益效果：
(1)本发明的方法和系统基于现有的病毒数据库，结合二代测序技术和三代测序技术的优势，对样本进行二代宏基因组测序和Nanopore三代宏基因组测序，加上优异的分析与筛选算法，可有效提高病毒序列的鉴定敏感性和准确性，并可有效提高病毒序列完整度。
(2)本发明方法和系统适用范围广，可覆盖不同类型的样本，并同时结合了短读长测序和长读长测序技术各自的优势和系统，可鉴定样本内的高完整度病毒序列。
(3)本发明分析方法和系统，提高了病毒鉴定的灵敏性和准确性，并且比现有分析方法快3.5倍以上。
(4)本发明的方法和系统可从测序数据中灵敏且准确地分析自然界的新病毒种，为病毒学研究提供技术支持。
附图说明
图1示出了本发明实施例2检测样本中宏病毒组的系统的示意图。
图2示出了本发明实施例3样本鉴定的cluster示例的共线性图。
具体实施方式
除非另有说明、从上下文暗示或属于现有技术的惯例，否则本申请中所有的份数和百分比都基于重量，且所用的测试和表征方法都是与本申请的提交日期同步的。在适用的情况下，本申请中涉及的任何专利、专利申请或公开的内容全部结合于此作为参考，且其等价的同族专利也引入作为参考，特别这些文献所披露的关于本领域中的合成技术、产物和加工设计、聚合物、共聚单体、引发剂或催化剂等的定义。如果现有技术中披露的具体术语的定义与本申请中提供的任何定义不一致，则以本申请中提供的术语定义为准。
本申请中的数字范围是近似值，因此除非另有说明，否则其可包括范围以外的数值。数值范围包括以1个单位增加的从下限值到上限值的所有数值，条件是在任意较低值与任意较高值之间存在至少2个单位的间隔。例如，如果记载组分、物理或其它性质(如分子量，熔体指数等)是100至1000，意味着明确列举了所有的单个数值，例如100，101，102等，以及所有的子范围，例如100到166，155到170，198到200等。对于包含小于1的数值或者包含大于1的分数(例如1.1，1.5等)的范围，则适当地将1个单位看作0.0001，0.001，0.01或者0.1。对于包含小于10(例如1到5)的个位数的范围，通常将1个单位看作0.1。这些仅仅是想要表达的内容的具体示例，并且所列举的最低值与最高值之间的数值的所有可能的组合都被认为清楚记载在本申请中。
关于化学化合物使用时，除非明确地说明，否则单数包括所有的异构形式，反之亦然(例如，“己烷”单独地或共同地包括己烷的全部异构体)。另外，除非明确地说明，否则用“一个”，“一种”或“该”形容的名词也包括其复数形式。
术语“包含”，“包括”，“具有”以及它们的派生词不排除任何其它的组分、步骤或过程的存在，且与这些其它的组分、步骤或过程是否在本申请中披露无关。为消除任何疑问，除非明确说明，否则本申请中所有使用术语“包含”，“包括”，或“具有”的组合物可以包含任何附加的添加剂、辅料或化合物。相反，出来对操作性能所必要的那些，术语“基本上由……组成”将任何其他组分、步骤或过程排除在任何该术语下文叙述的范围之外。术语“由……组成”不包括未具体描述或列出的任何组分、步骤或过程。除非明确说明，否则术语“或”指列出的单独成员或其任何组合。
为了使本发明所解决的技术问题、技术方案及有益效果更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。
实施例
以下例子在此用于示范本发明的优选实施方案。本领域内的技术人员会明白，下述例子中披露的技术代表发明人发现的可以用于实施本发明的技术，因此可以视为实施本发明的优选方案。但是本领域内的技术人员根据本说明书应该明白，这里所公开的特定实施例可以做很多修改，仍然能得到相同的或者类似的结果，而非背离本发明的精神或范围。
除非另有定义，所有在此使用的技术和科学的术语，和本发明所属领域内的技术人员所通常理解的意思相同，在此公开引用及他们引用的材料都将以引用的方式被并入。
那些本领域内的技术人员将意识到或者通过常规试验就能了解许多这里所描述的发明的特定实施方案的许多等同技术。这些等同将被包含在权利要求书中。
下述实施例中未作具体说明的分子生物学实验方法，均按照《分子克隆实验指南》(第四版)(J.萨姆布鲁克、M.R.格林，2017)一书中所列的具体方法进行，或者按照试剂盒和产品说明书进行。其他实验方法，如无特殊说明，均为常规方法。下述实施例中所用的仪器设备，如无特殊说明，均为实验室常规仪器设备；下述实施例中所用的试验材料，如无特殊说明，均为自常规生化试剂商店购买得到的。
实施例1检测样本中宏病毒组的方法
一、样本处理
将样品分为2份，分别进行二代和三代测序平台的实验处理，具体如下：
1.二代测序
1)首先使用酚-氯仿抽提法提取总DNA。核酸提取后，使用Thermo NanoDrop 2000和Qubit两种仪器分别进行核酸浓度的测定和质量检测；
2)用磁珠法对DNA片段进行纯化和洗脱，并进行片段末端修复；
3)对DNA片段进行加接头和扩增，构建测序文库；
4)使用Illumina Novaseq平台进行PE150的宏基因组测序，测序下机后对原始bam文件进行数据格式转换和index拆分，得到样本的测序原始数据(fastq格式)。
2.三代测序：
1)首先使用酚-氯仿抽提法提取总DNA。核酸提取后，使用Thermo NanoDrop 2000和Qubit两种仪器分别进行核酸浓度测定和质量检测；
2)用磁珠法对DNA片段进行纯化和洗脱，并进行片段末端修复，最后加barcode，构建测序文库；
3)使用Nanopore MinION平台进行宏基因组测序，测序下机后对原始fast5文件使用Guppy v3.0.3进行basecalling和barcode拆分，得到样本的测序原始数据(fastq格式)。
二、测序数据分析
1.测序数据前处理
二代数据：使用fastp v0.20对测序原始数据进行质控，参数设置为-n 2-q 15-x-5-3-y--dedup。首先过滤含有2个以上N碱基的read，对保留下来的序列，在序列首尾两端使用滑动窗口检测碱基质量，切除质量值低于Q15的短片段序列，同时检测序列尾端的polyX序列并去除。此外，还过滤低复杂度序列、重复序列和接头序列。上述所有过滤步骤通过后，过滤长度小于15bp的序列，得到二代clean reads。
三代数据：只保留质量值大于7的序列，得到三代clean reads。
2.去除宿主序列
如样本有明确的宿主来源，则：
1)对于二代数据：使用bwa v0.7.17的默认参数将二代clean reads比对到宿主参考基因组，并用samtools v1.9提取未比对上的序列，得到二代非宿主序列；
2)对于三代数据：使用minimap2 v2.2的默认参数将三代clean reads比对到宿主参考基因组，并用samtools v1.9提取未比对上的序列，得到三代非宿主序列。
3.数据组装
1)混合组装：使用OPERA-MS v0.9.0软件将上述二代和三代数据进行相互校正与混合组装；
2)纯三代组装：使用flye v2.8.3对三代数据进行组装。
对两种组装结果均只保留序列长度大于或等于500bp以上的序列，分别得到混合组装contigs和Nanopore contigs。
4.病毒鉴定和物种注释：
1)混合组装contigs：
a)病毒注释：使用CheckV v0.7软件的默认参数对混合组装contigs进行病毒序列的比对注释，得到病毒注释结果；
其中CheckV共分为三个步骤来注释病毒序列，并对潜在的病毒序列进行质量评价：①使用Prodigal v2.6.3对序列进行基因预测；②使用hhmer v3.3将序列与病毒特异性HMM模型(隐马模型)数据和非病毒微生物特异性HMM模型数据进行比对，识别并去除序列边缘的非病毒区域，并统计宿主基因和病毒基因数量；③使用Diamond v2.0.6将上述序列比对到病毒特异性蛋白库(氨基酸比对方法)，通过计算平均氨基酸同一性(aai_id)和比对覆盖率(aai_af)来评估病毒基因组的完整性和病毒序列的置信度；④对于氨基酸比对方法鉴定的低置信度(比对错误率高于10％)病毒序列，进一步使用隐马模型比对的方法来评估病毒基因组的完整性和病毒序列的置信度。
b)病毒筛选：
i.过滤checkv_quality栏标记为Not-determined(未鉴定到病毒基因)的序列；
ii.过滤注释上的宿主基因数量比病毒基因数量大于10倍的序列，得到候选病毒contigs和非病毒contigs；
c)病毒物种注释：使用vpf-tools软件对候选病毒contigs进行病毒物种注释，获得分类学谱系信息，得到候选病毒contigs的注释结果。
2)纯三代组装：
按上述步骤鉴定Nanopore contigs里包含的病毒，并进行物种注释，记为Nanopore病毒contigs。
5.病毒鉴定的校正：
步骤4得到的结果存在病毒序列遗漏、同一株病毒的多个碎片序列不能有效归纳和物种注释不全等问题。
聚类分析：将步骤4得到的候选病毒contigs、非病毒contigs和Nanopore病毒contigs共3个数据集使用cd-hit v4.8.1的PSI-CD-HIT脚本(参数设置为-c 0.95 -aS0.85 -para 4 -blp 20)进行聚类分析。并按以下3个条件过滤聚类结果的聚类单元：
a)Nanopore病毒contig作为代表序列；
b)序列数大于等于3；
c)序列来自2个以上数据集，得到最终的聚类结果。
病毒鉴定结果优化：对聚类单元(cluster)逐个进行比对分析，
a)使用Megablast v2.2.26将非代表序列集比对至该cluster的代表序列(Nanopore病毒contig)；
b)只保留evlaue<1e-8，且输入序列的coverage>80％的比对记录；
c)若b)中包含非病毒contigs，则作为遗漏的病毒序列进行补充进候选病毒contigs，得到校正后的候选病毒contigs。需注意：Nanopore病毒contigs只作为代表序列，不补充进病毒序列集；
d)将cluster里的校正后的候选病毒contigs作为一个类病毒株，并使用代表序列的物种注释结果作为该类病毒株的分类学谱系信息；
e)对b)的结果进行共线性分析，对序列进行排序，并记录位置坐标及序列方向信息。
病毒鉴定结果汇总：将不在筛选的聚类结果里的每条候选病毒contig均作为类病毒株，与经过校正得到的类病毒株的注释结果汇总，并作为最终的病毒鉴定结果。
实施例2检测样本中宏病毒组的系统
本实施例提供实现实施例1检测样本中宏病毒组的方法的系统，如图1所示，包括：
测序数据存储模块，用于获得待检测样本二代测序数据和三代测序数据并存储；
组装模块，与测序数据存储模块连接，用于将二代测序数据和三代测序数据进行相互校正和混合组装，得到混合组装contigs，并用于将三代测序数据单独组装得到Nanopore contigs；
病毒鉴定与注释模块，与组装模块连接，用于将混合组装contigs进行比对注释，得到候选病毒contigs和非病毒contigs，并对鉴定的病毒contigs进行物种注释，得到候选病毒的contigs的注释结果；并用于将Nanopore contigs进行比对注释，得到Nanopore病毒contigs，同样进行物种注释得到Nanopore病毒contigs的注释结果；
病毒鉴定校正模块，与病毒鉴定与注释模块连接，用于：
将候选病毒contigs、非病毒contigs和Nanopore病毒contigs三个数据集进行聚类分析，筛选出代表序列来自Nanopore病毒contigs的聚类单元，该聚类单元中如果某个非病毒contig能够比对上代表序列，则该非病毒contig作为遗漏的序列补充进候选病毒contigs，得到校正后的候选病毒contigs，将该聚类单元里的校正后的候选病毒contigs作为一个类病毒株，并使用代表序列的物种注释结果作为该类病毒株的注释结果；将不在筛选出的聚类单元里的每条候选病毒contig均作为类病毒株，与经过校正得到的类病毒株的注释结果进行汇总。
实施例3人粪便样本进行宏基因组测序并分析
一、样本处理
按实施例1的方法将样本分为2份以后，分别进行二代宏基因组测序和三代宏基因组测序。
二、测序数据分析
按实施例1方法得到样本的二代测序数据量约为10Gb bases，三代测序数据量约为6Gb bases。去除人宿主序列后，二代数据剩余9.2Gb，三代数据剩余5.64Gb。对二代和三代质控数据进行混合组装得到23万个contigs，共512.8Mb bases；而对三代质控数据进行单独的组装得到2000多个contigs，共74.1Mb bases(表1)。
表1 样本数据组装统计表
说明：
Total number：组装contigs数量，长度均大于500bp；
Total len(Mb)：组装contigs总长度；
Max len：组装contigs最大长度；
N50 len：将组装contigs按长度从大到小排序，并对其长度进行累加，当累加长度达到全部contigs长度的一半时，该contigs的长度。
将组装的序列进行初步的病毒鉴定，混合组装有7497条contig注释为病毒序列，而三代组装有296条病毒序列(表2)。而初步的病毒物种注释统计结果见表3，可以明显看到序列长度对病毒物种注释的精确度有显著影响：每个分类水平下，三代组装序列的注释百分比均比混合组装要高(科：58.4％vs 20.7％)。
表2 样本初步鉴定的病毒统计表
说明：
Totalnumber：病毒contigs数量，长度均大于500bp；
Totallen(Mb)：病毒contigs总长度；
Maxlen：病毒contigs最大长度；
N50len：将病毒contigs按长度从大到小排序，并对其长度进行累加，当累加长度达到全部contigs长度的一半时，该contigs的长度。
表3 样本初步鉴定的病毒物种注释统计表
说明：
(1)Total：病毒contigs总数；
(2)Level_D：分类为域的病毒contigs数目；
(3)Level_P：分类为门的病毒contigs数目；
(4)Level_C：分类为纲的病毒contigs数目；
(5)Level_O：分类为目的病毒contigs数目；
(6)Level_F：分类为科的病毒contigs数目；
(7)Level_G：分类为属的病毒contigs数目；
(8)Unknown：未知分类的病毒contigs数目。
经过病毒鉴定结果的校正分析后，得到最终的病毒鉴定结果，统计结果见表4和表5：
表4 样本最终鉴定的病毒统计表
说明：
Cluster数：类病毒株数目；
Totalnumber：最终鉴定的病毒contigs数量，长度均大于500bp；
Totallen(Mb)：最终鉴定的病毒contigs总长度；
Maxlen：最终鉴定的病毒contigs最大长度；
N50len：将最终鉴定的病毒contigs按长度从大到小排序，并对其长度进行累加，当累加长度达到全部contigs长度的一半时，该contigs的长度。
表5 样本最终鉴定的病毒物种注释统计表
说明：
(1)Total：最终鉴定的病毒contigs总数；
(2)Level_D：最终鉴定的分类为域的病毒contigs数目；
(3)Level_P：最终鉴定的分类为门的病毒contigs数目；
(4)Level_C：最终鉴定的分类为纲的病毒contigs数目；
(5)Level_O：最终鉴定的分类为目的病毒contigs数目；
(6)Level_F：最终鉴定的分类为科的病毒contigs数目；
(7)Level_G：最终鉴定的分类为属的病毒contigs数目；
(8)Unknown：未知分类的病毒contigs数目。
可知，病毒数量提高到11547条(相比初步鉴定的7497条序列，提高了1.54倍)，且显著提高了序列的物种注释精确度(科水平：30.8％vs 20.7％；属水平：2.6％vs 1.2％)，并按cluster的关系，有效归纳了片段序列，共得到1675个类病毒株。而其中一个典型的类病毒株(cluster)的信息见表6，此类病毒株来自有尾噬菌体目，代表序列长363kb，包括候选病毒contigs集的3条病毒contig，以及非病毒contigs集的17条序列，大大提高了病毒鉴定的敏感性，通过共线性可视化分析对序列集进行了排序，并记录位置坐标及序列方向的信息(图2)。
表6 样本鉴定的cluster示例
实施例4方法的时效性、计算资源消耗检测
一、数据
为了测试本发明实施例1的运算时间和计算资源消耗，使用5套病毒宏基因组数据进行测试(Illumina PE150：碱基数约10Gb；Nanopore：碱基数约6Gb)，比较本发明与普通的宏病毒组分析方法(Guoyan Zhao,Guang Wub,Efrem S.Limb,et al.VirusSeeker,acomputational pipeline for virus discovery and virome compositionanalysis.virology.503(2017)21-30.doi:10.1016/j.virol.2017.01.005.，全文并入此处)在分析用时、计算资源消耗的差异。
二、结果
从表7可了解到，同样的CPU核数分析约10Gb数据量的样本，普通的病毒宏基因组分析流程比本发明方法耗时高约4倍，服务器内存资源消耗高约3倍。
表7 实施例1方法与普通的宏病毒组分析流程的时效性、计算资源消耗对比表
普通方法使用NCBI nt/nr全库，造成分析非常耗时且漫长，并且病毒的鉴定不灵敏，而本发明通过使用病毒特异性蛋白序列数据库、加上精心设计的分析策略和参数优化节省了分析时长、降低了内存消耗，达到快速且准确分析宏基因组测序数据中的病毒组成。
实施例5方法的准确性分析
一、数据
为了评价本发明实施例1的分析准确性，用RefSeq(ftp://ftp.ncbi.nlm.nih.gov/genomes)数据库制作模拟数据，步骤如下：
1.参考序列组成：
随机挑选130种病毒科，每种科随机挑选1株病毒参考基因组组成病毒序列库；随机挑选3种细菌参考基因组；纳入3种常见宿主(人、老鼠和猪)的参考基因组；
2.组装序列模拟：
i.模拟混合组装序列：采用期望值1kb的正态分布模型，对参考序列集进行随机打断(window size:0)，生成模拟数据集1；
ii.模拟Nanopore组装序列：采用期望值25kb的正态分布模型，对参考序列集进行随机打断(window size:0)，生成模拟数据集2。另外，两份模拟数据集里，3种宿主均只纳入100条打断的序列(随机挑选)。具体的数据情况见表8。
表8 模拟数据表
二、结果
模拟数据集按实施例1描述的方法进行宏病毒组分析，并同时也用实施例4使用的普通的宏病毒测序数据分析流程进行分析(只支持输入数据集1)，得到病毒鉴定结果的统计表9。
表9 实施例1方法与普通的宏病毒组分析流程的准确性对比
从结果可知，利用实施例1的方法初步的病毒鉴定结果(未进行病毒鉴定结果的校正)，其敏感性和精确度均优于普通方法(敏感性：82.5％vs 64.8％；精确度：89.6％vs71.5％)，其得益于本发明使用的CheckV软件配套的病毒特异性蛋白序列数据库和病毒特异性隐马模型数据库。
更为重要的是，利用实施例1的方法对病毒鉴定结果进行校正后的分析效果明显优于校正前和普通方法，敏感性高达97.8％，而精确度高达98.6％，其得益于本发明精心设计的分析策略，充分利用了Nanopore三代测序技术的优势，使得分析结果兼具高敏感度和高精确度。
实施例6方法的新病毒种鉴定
一、数据
为了评价实施例1的方法鉴定新病毒种的性能，采集1例土壤样本，按照实施例1的方法对该例样本进行病毒宏基因组测序。
二、测序数据分析
按照实施例1中的模式进行宏病毒组分析，得到新病毒种鉴定表10。
表10 样本鉴定的新病毒种信息表
| 注释栏目 | 注释信息 |
| 分类学谱系 | D__Viruses.P__Uroviricota.C__Caudoviricetes.O__Caudovirales |
| 病毒类型 | 噬菌体 |
| 序列长度(bp) | 36,378 |
| 基因预测与注释 | 基因总数：60；病毒基因：7；宿主基因：0 |
| 完整度 | 完整基因组 |
| 置信度 | 高 |
| 比对一致性(％) | 53.55 |
| 比对覆盖度(％) | 96.1 |
结果显示此病毒序列与目标序列的相似性只有53.5％，而比对覆盖度高达96.1％，且具备完整病毒基因组的结构，故提示是一个有尾噬菌体目下的一个高质量的新病毒种，表明本发明的分析方法具有优异的高质量新病毒种的鉴定性能。
在本发明提及的所有文献都在本申请中引用作为参考，就如同每一篇文献被单独引用作为参考那样。此外应理解，在阅读了本发明的上述讲授内容之后，本领域技术人员可以对本发明作各种改动或修改，这些等价形式同样落于本申请所附权利要求书所限定的范围。
Claims (10)
1.一种检测样本中宏病毒组的方法，其特征在于，包括以下步骤：
S1，分别获得待检测样本二代测序数据和三代测序数据；
S2，将所述二代测序数据和所述三代测序数据进行相互校正和混合组装，得到混合组装contigs，将三代测序数据单独组装得到Nanopore contigs；
S3，病毒鉴定与物种注释：
将所述混合组装contigs进行比对注释，得到候选病毒contigs和非病毒contigs，并对鉴定的病毒contigs进行物种注释，得到候选病毒的contigs的注释结果，
将所述Nanopore contigs进行比对注释，得到Nanopore病毒contigs，同样进行物种注释得到Nanopore病毒contigs的注释结果；
S4，病毒鉴定的校正：
S41，将候选病毒contigs、非病毒contigs和Nanopore病毒contigs三个数据集进行聚类分析，
S42，筛选出代表序列来自Nanopore病毒contigs的聚类单元，该聚类单元中如果某个非病毒contig能够比对上所述代表序列，则该非病毒contig作为遗漏的序列补充进候选病毒contigs，得到校正后的候选病毒contigs，将该聚类单元里的校正后的候选病毒contigs作为一个类病毒株，并使用代表序列的物种注释结果作为该类病毒株的注释结果，
S43，将不在筛选出的聚类单元里的每条候选病毒contig均作为类病毒株，与经过校正得到的类病毒株的注释结果进行汇总，得到最终的病毒鉴定结果。
2.根据权利要求1所述的一种检测样本中宏病毒组的方法，其特征在于，步骤S1中，获得的所述二代测序数据和所述三代测序数据为经过质控和去除宿主序列后的二代测序数据和三代测序数据。
3.根据权利要求2所述的一种检测样本中宏病毒组的方法，其特征在于，对于二代测序数据，进行质控的步骤如下：
(1)过滤含有2个以上N碱基的read，对保留下来的序列，在序列首尾两端使得滑动窗口检测碱基质量，切除质量值低于Q15的短片段序列，同时检测序列尾端的polyX序列并去除；
(2)过滤低复杂度序列、重复序列和接头序列；
(3)过滤长度小于15bp的序列，
对于三代测序数据，进行质控的步骤为：只保留质量值大于7的序列。
4.根据权利要求1所述的一种检测样本中宏病毒组的方法，其特征在于，步骤S2中，混合组装contigs和Nanopore contigs只保留长度大于或等于500bp的序列。
5.根据权利要求1所述的一种检测样本中宏病毒组的方法，其特征在于，步骤S3中，所述比对注释是指将混合组装contigs或Nanopore contigs与病毒特异性蛋白序列数据库比对。
6.根据权利要求5所述的一种检测样本中宏病毒组的方法，其特征在于，步骤S3中，过滤掉比对结果中未鉴定到病毒基因的序列，并过滤掉注释上的宿主基因数量比病毒基因数量大于10倍的序列，得到候选病毒contigs、非病毒contigs和Nanopore病毒contigs。
7.根据权利要求1所述的一种检测样本中宏病毒组的方法，其特征在于，步骤S42中，筛选出的聚类单元中，还满足以下条件：
(1)序列数大于3；
(2)序列来自所述三个数据集中的至少2个。
8.根据权利要求1所述的检测样本中宏病毒组的方法，其特征在于，步骤S42中，所述能够比对上是指比对结果evlaue<1e-8，且输入序列的coverage>80％。
9.一种检测样本中宏病毒组的系统，其特征在于，包括：
测序数据存储模块，用于获得待检测样本二代测序数据和三代测序数据并存储；
组装模块，与所述测序数据存储模块连接，用于将所述二代测序数据和所述三代测序数据进行相互校正和混合组装，得到混合组装contigs，并用于将三代测序数据单独组装得到Nanopore contigs；
病毒鉴定与注释模块，与所述组装模块连接，用于将所述混合组装contigs进行比对注释，得到候选病毒contigs和非病毒contigs，并对鉴定的病毒contigs进行物种注释，得到候选病毒的contigs的注释结果；并用于将所述Nanopore contigs进行比对注释，得到Nanopore病毒contigs，同样进行物种注释得到Nanopore病毒contigs的注释结果；
病毒鉴定校正模块，与所述病毒鉴定与注释模块连接，用于：
将候选病毒contigs、非病毒contigs和Nanopore病毒contigs三个数据集进行聚类分析，筛选出代表序列来自Nanopore病毒contigs的聚类单元，该聚类单元中如果某个非病毒contig能够比对上所述代表序列，则该非病毒contig作为遗漏的序列补充进候选病毒contigs，得到校正后的候选病毒contigs，将该聚类单元里的校正后的候选病毒contigs作为一个类病毒株，并使用代表序列的物种注释结果作为该类病毒株的注释结果；将不在筛选出的聚类单元里的每条候选病毒contig均作为类病毒株，与经过校正得到的类病毒株的注释结果进行汇总。
10.根据权利要求9所述的一种检测样本中宏病毒组的系统，其特征在于，进一步包括测序数据前处理模块，用于对二代测序数据或三代测序数据进行质控和和去除宿主序列。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111407911.2A CN114121160B (zh) | 2021-11-25 | 2021-11-25 | 一种检测样本中宏病毒组的方法和系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111407911.2A CN114121160B (zh) | 2021-11-25 | 2021-11-25 | 一种检测样本中宏病毒组的方法和系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114121160A CN114121160A (zh) | 2022-03-01 |
| CN114121160B true CN114121160B (zh) | 2022-06-21 |
Family
ID=80372716
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111407911.2A Active CN114121160B (zh) | 2021-11-25 | 2021-11-25 | 一种检测样本中宏病毒组的方法和系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114121160B (zh) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115101129B (zh) * | 2022-06-27 | 2023-03-24 | 青岛华大医学检验所有限公司 | 一种基于宏基因组测序数据组装病原微生物基因组的方法 |
| CN115198036B (zh) * | 2022-09-13 | 2022-12-30 | 江苏省环境工程技术有限公司 | 一种基于纳米孔和高通量测序数据的噬菌体鉴定和宿主预测方法 |
| CN116072222B (zh) * | 2023-02-16 | 2024-02-06 | 湖南大学 | 病毒基因组鉴定和拼接的方法及应用 |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106536756A (zh) * | 2014-06-26 | 2017-03-22 | 10X基因组学有限公司 | 核酸序列的分析 |
| CN108629156A (zh) * | 2017-03-21 | 2018-10-09 | 深圳华大基因科技服务有限公司 | 三代测序数据纠错的方法、装置和计算机可读存储介质 |
| CN110168105A (zh) * | 2016-12-09 | 2019-08-23 | 瑞泽恩制药公司 | 用于对t细胞受体进行测序的系统和方法及其用途 |
| CN110858503A (zh) * | 2018-08-11 | 2020-03-03 | 中国科学院昆明动物研究所 | 综合应用第三代超长测序读段和第二代链接式读段从头组装基因组的方法 |
| WO2020072858A1 (en) * | 2018-10-04 | 2020-04-09 | T2 Biosystems, Inc. | Methods and compositions for high sensitivity detection of drug resistance markers |
| CN113077842A (zh) * | 2021-03-25 | 2021-07-06 | 北京百迈客生物科技有限公司 | 一种三代全长转录组辅助基因预测的方法 |
| CN113611357A (zh) * | 2020-11-17 | 2021-11-05 | 上海美吉生物医药科技有限公司 | 基于宏基因组的抗性基因分析方法、装置、介质及终端 |
-
2021
- 2021-11-25 CN CN202111407911.2A patent/CN114121160B/zh active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106536756A (zh) * | 2014-06-26 | 2017-03-22 | 10X基因组学有限公司 | 核酸序列的分析 |
| CN110168105A (zh) * | 2016-12-09 | 2019-08-23 | 瑞泽恩制药公司 | 用于对t细胞受体进行测序的系统和方法及其用途 |
| CN108629156A (zh) * | 2017-03-21 | 2018-10-09 | 深圳华大基因科技服务有限公司 | 三代测序数据纠错的方法、装置和计算机可读存储介质 |
| CN110858503A (zh) * | 2018-08-11 | 2020-03-03 | 中国科学院昆明动物研究所 | 综合应用第三代超长测序读段和第二代链接式读段从头组装基因组的方法 |
| WO2020072858A1 (en) * | 2018-10-04 | 2020-04-09 | T2 Biosystems, Inc. | Methods and compositions for high sensitivity detection of drug resistance markers |
| CN113611357A (zh) * | 2020-11-17 | 2021-11-05 | 上海美吉生物医药科技有限公司 | 基于宏基因组的抗性基因分析方法、装置、介质及终端 |
| CN113077842A (zh) * | 2021-03-25 | 2021-07-06 | 北京百迈客生物科技有限公司 | 一种三代全长转录组辅助基因预测的方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114121160A (zh) | 2022-03-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN114121160B (zh) | 一种检测样本中宏病毒组的方法和系统 | |
| Almeida et al. | Bioinformatics tools to assess metagenomic data for applied microbiology | |
| US20200294628A1 (en) | Creation or use of anchor-based data structures for sample-derived characteristic determination | |
| García-López et al. | Fragmentation and coverage variation in viral metagenome assemblies, and their effect in diversity calculations | |
| CN113160882B (zh) | 一种基于三代测序的病原微生物宏基因组检测方法 | |
| CN110875082B (zh) | 一种基于靶向扩增测序的微生物检测方法和装置 | |
| CN107533589A (zh) | 生物信息学数据处理系统 | |
| CN114420212B (zh) | 一种大肠杆菌菌株鉴定方法和系统 | |
| CN112687344A (zh) | 一种基于宏基因组的人腺病毒分子分型和溯源方法及系统 | |
| CN113470743A (zh) | 一种基于bd单细胞转录组和蛋白组测序数据的差异基因分析方法 | |
| CN115662516A (zh) | 一种基于二代测序技术的高通量预测噬菌体宿主的分析方法 | |
| Schmitz et al. | Quality control and evaluation of plant epigenomics data | |
| Hebert et al. | Interrogating 1000 insect genomes for NUMTs: A risk assessment for estimates of species richness | |
| Kaiser et al. | Automated structural variant verification in human genomes using single-molecule electronic DNA mapping | |
| CN110751985B (zh) | 与大体重鸡只高度关联的肠道微生物标记物 | |
| CN115938491B (zh) | 一种用于临床病原诊断的高质量细菌基因组数据库构建方法及系统 | |
| CN115198036B (zh) | 一种基于纳米孔和高通量测序数据的噬菌体鉴定和宿主预测方法 | |
| CN116646010B (zh) | 人源性病毒检测方法及装置、设备、存储介质 | |
| Sengupta et al. | Classification and identification of fungal sequences using characteristic restriction endonuclease cut order | |
| CN113355438B (zh) | 一种血浆微生物物种多样性评估方法、装置和存储介质 | |
| CN110331210B (zh) | 一套用于获取馆藏甲虫标本DNA条形码的Mini-Barcoding引物及其应用 | |
| Teo et al. | A comparative study of metagenomics analysis pipelines at the species level | |
| Deng et al. | Bioinformatics Analysis for NGS Amplicon Sequencing | |
| Bester et al. | Validation of High-Throughput Sequencing (HTS) for Routine Detection of Citrus Viruses and Viroids | |
| CN115449545A (zh) | 利用麦类作物d基因组外显子捕获测序芯片鉴定节节麦衍生种质片段的方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |