一种宏病毒组流程的优化分析方法
阅读说明：本技术 一种宏病毒组流程的优化分析方法 (Optimized analysis method for macrovirome process ) 是由 薛正晟 李鸿毅 寇文伯 孙子奎 于 2020-12-29 设计创作，主要内容包括：本发明公开了一种宏病毒组流程的优化分析方法,1)序列质控,2)序列拼接,3)病毒序列预测,4)病毒序列聚类处理,5)物种注释,本发明分析数据更为准确,能拼接出长度较长的contig,从而提高比对结果的准确性和注释结果的可靠性。(The invention discloses an optimized analysis method of a macrovirus group process, which comprises the steps of 1) sequence quality control, 2) sequence splicing, 3) virus sequence prediction, 4) virus sequence clustering processing and 5) species annotation.)
技术领域
本发明涉及基因检测技术领域，尤其涉及一种宏病毒组流程的优化分析方法。
背景技术
宏病毒组测序是指对样本种所有病毒基因组进行高通量测序，通过分析测序序列的构成分析特定环境中病毒群体的构成情况或基因的组成以及功能。借助不同环境下病毒的构成差异分析我们可以分析病毒与宿主之间的关系，寻找标志性病毒或特定功能的基因。为了实现着一些目的。当前运用于宏病毒组分析的最常用方法，是基于reads分析方法。该方法的原理是根据reads与病毒数据库进行比对，将比对上的序列注释为该病毒。但是，由于reads比对的可靠性不高，准确性较差，所得到的结果往往会受到广大的质疑，因此，需要引入基于拼接的数据分析方法。
而现有的基于reads分析方法具有如下缺陷：
(1)Reads太短，往往只有150bp，由此去比对，往往得到的结果为假阳性；
(2)基于reads的分析方法，其结果往往包含了大量的病毒结果，而大部分结果都不可信。
发明内容
本发明的提供一种宏病毒组流程的优化分析方法。
本发明的方案是：
一种宏病毒组流程的优化分析方法，包括下列步骤：
1)序列质控，使用fastp软件，将原始序列进行质量过滤，获得高质量序列；
2)序列拼接，使用megahit软件包对上一步得到的高质量序列进行拼接，得到拼接好的contig文件；
3)病毒序列预测，将上一步得到的contig文件使用virsorter软件进行病毒序列的预测，得到预测为病毒的contig序列；
4)病毒序列聚类处理，将步骤3)中挑选预测为病毒的contig序列，使用vcontact2软件，对序列进行聚类分析，获得聚类的结果文件；
5)物种注释，将上一步得到的聚类结果，使用lca的算法，计算出病毒序列的注释结果。
作为优选的技术方案，所述步骤1)中的原始数据是采用Cutadapt识别3'端潜在的接头序列，并在识别的接头序列处截断；与接头序列的匹配长度≥3bp，允许≤20％的碱基错配率。
作为优选的技术方案，所述步骤1)中fast软件p采用滑动窗口法对序列进行质量筛查：窗口大小为5bp，从5'端第一个碱基位置开始移动，要求窗口中碱基平均质量≥Q20，从第一个平均质量值低于Q20的窗口的3'端碱基处截断序列；经上述质量筛查后，去除序列长度小于50bp的序列；经上述质量筛查后，去除序列中含有模糊碱基的序列。
作为优选的技术方案，所述步骤3)中所述contig长度大于1kb；Virsorter取出其中的category1，category2，category4，category5；余下大于5kb的contigs提取这部分contigs的蛋白序列做vHMM。
作为优选的技术方案，所述contigs的蛋白序列做vHMM如下：
一个contigs至少要求五个蛋白有HMM的hit，占10％以上，有KO注释的基因小于20％，PFAM注释的基因小于40％；
有HMM hit的基因数大于等于PFAM hit的基因数,；
HMM hit大于等于60％的总基因数。
作为优选的技术方案，所述步骤5)中使用vContact2软件对预测为病毒的contig序列进行物种注释。
本发明提供了一种宏病毒组流程的优化分析方法，1)序列质控，使用fastp软件，将原始序列进行质量过滤，获得高质量序列；2)序列拼接，使用megahit软件包对上一步得到的高质量序列进行拼接，得到拼接好的contig文件；3)病毒序列预测，将上一步得到的contig文件使用virsorter软件进行病毒序列的预测，得到预测为病毒的contig序列；4)病毒序列聚类处理，将步骤3)中挑选预测为病毒的contig序列，使用vcontact2软件，对序列进行聚类分析，获得聚类的结果文件；5)物种注释，将上一步得到的聚类结果，使用lca的算法，计算出病毒序列的注释结果。
本发明的优点：分析数据相较基于reads的方法更为准确,能拼接出长度较长的contig，从而提高比对结果的准确性和注释结果的可靠性；
以vcontact2软件将病毒contig进行聚类，基于lca算法对病毒contig进行注释，结果可靠，避免了假阳性问题，数据分析更精准。
附图说明
图1为本发明实施例的流程示意图。
具体实施方式
为了弥补以上不足，本发明提供了一种宏病毒组流程的优化分析方法以解决上述背景技术中的问题。
一种宏病毒组流程的优化分析方法，包括下列步骤：
1)序列质控，使用fastp软件，将原始序列进行质量过滤，获得高质量序列；
2)序列拼接，使用megahit软件包对上一步得到的高质量序列进行拼接，得到拼接好的contig文件；
3)病毒序列预测，将上一步得到的contig文件使用virsorter软件进行病毒序列的预测，得到预测为病毒的contig序列；
4)病毒序列聚类处理，将步骤3)中挑选预测为病毒的contig序列，使用vcontact2软件，对序列进行聚类分析，获得聚类的结果文件；
5)物种注释，将上一步得到的聚类结果，使用lca的算法，计算出病毒序列的注释结果。
所述步骤1)中的原始数据是采用Cutadapt识别3'端潜在的接头序列，并在识别的接头序列处截断；与接头序列的匹配长度≥3bp，允许≤20％的碱基错配率。
所述步骤1)中fast软件p采用滑动窗口法对序列进行质量筛查：窗口大小为5bp，从5'端第一个碱基位置开始移动，要求窗口中碱基平均质量≥Q20，从第一个平均质量值低于Q20的窗口的3'端碱基处截断序列；经上述质量筛查后，去除序列长度小于50bp的序列；经上述质量筛查后，去除序列中含有模糊碱基的序列。
所述步骤3)中所述contig长度大于1kb；Virsorter取出其中的category1，category2，category4，category5；余下大于5kb的contigs提取这部分contigs的蛋白序列做vHMM。
所述contigs的蛋白序列做vHMM如下：
一个contigs至少要求五个蛋白有HMM的hit，占10％以上，有KO注释的基因小于20％，PFAM注释的基因小于40％；
有HMM hit的基因数大于等于PFAM hit的基因数,；
HMM hit大于等于60％的总基因数。
所述步骤5)中使用vContact2软件对预测为病毒的contig序列进行物种注释。
为了使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合具体实施例，进一步阐述本发明。
实施例：
在步骤S101中，进行序列质控，将现有的原始序列进行质量控制处理，得到高质量序列。
在本发明实施例中，采用Cutadapt识别3'端潜在的接头序列(极少数测通的情况会出现)，并在识别的接头序列处截断。要求与接头序列(R1:AGATCGGAAGAGCACACGTCTGAACTCCAGTCA；R2:AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT)的匹配长度至少达到3bp，且允许至多20％的碱基错配率。
去除3'端接头序列后，使用fastp采用滑动窗口法对序列进行质量筛查：窗口大小为5bp，从5'端第一个碱基位置开始移动，要求窗口中碱基平均质量大于等于Q20(即碱基平均测序准确率大于99％)，从第一个平均质量值低于Q20的窗口的3'端碱基处截断序列。
经上述质量筛查后，去除序列长度小于50bp的序列。
经上述质量筛查后，去除序列中含有模糊碱基的序列。
在步骤S102中，使用megahit对上一步得到的高质量蓄力进行拼接，得到拼接好的contig文件。
在本发明实施例中，首先以K-mer～[27,127]的参数设置，调用megahit(https://hku-bal.github.io/megabox/)，对每个样本的双端序列分别进行从头(De novo)组装拼接，通过De Bruijn图构建Contigs序列，并对生成的Contigs序列进行组装拼接效果评价。
在步骤S103中，将上一步得到的contig使用virsorter软件进行病毒序列的预测，得到预测为病毒的contig序列。
在本发明实施例中，先使用Virsorter预测contigs(只取长度大于1kb；Virsorter更新PFAM数据库为最新)，取出其中的category1，2，4，5(依据文献“Uncovering Earth’svirome”等)；余下的contigs(要求大于5kb)，提取这部分contigs的蛋白序列做vHMM(依据文献Uncovering Earth’s virome，下载HMM模型)，要求如下：
即：1.一个contigs至少要求五个蛋白有HMM的hit，且占10％以上，有KO注释的基因小于20％，PFAM注释的基因小于40％。
2.有HMM hit的基因数大于等于PFAM hit的基因数。(与以上并列)
3.HMM hit大于等于60％的总基因数。(与以上并列)。
在步骤S104中，挑选预测为病毒的contig序列，使用vcontact2软件，对序列进行聚类分析，获得聚类的结果文件。
在步骤S105中，将上一步得到的聚类结果，使用lca的算法，计算出病毒序列的注释结果。
在本发明实施例中，使用vContact2软件对上述预测为病毒的contig序列进行物种注释，由于每一条目标序列可能和多多条参考序列聚类在一起，而这些匹配的参考序列又分属不同的分类单元，为使分析严谨可靠，同时又不丢失生物学意义，我们采取“最低共同祖先(Lowest Common Ancestor，LCA)”算法，将参考序列分化为不同物种分枝前的最后一级共同分类，作为目标序列的物种分类注释信息。
以上显示和描述了本发明的基本原理、主要特征及本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。