一种宏病毒组数据的分析方法

本发明涉及一种宏病毒组数据的分析方法，属于宏病毒组学。

背景技术：

1、病毒是地球上数量最多的一种类生命体，对地球生物化学循环具有重要的贡献。病毒宏基因组测序又称宏病毒组(virome)，是在宏基因组学理论的基础上，结合现有的病毒分子生物学检测技术而兴起的一个新的学科分支。

2、宏病毒组直接以样本中所有病毒的遗传物质(dna或者rna)为研究对象，先富集病毒颗粒再获取基因组序列信息后鉴定其中所有的病毒组成及相对丰度，是一种发现新病毒、病毒感染预警和控制的有力手段，在病毒的起源和进化模式、遗传多样性和地理分布、以及病毒和宿主的相互关系等研究领域都具有重要意义。

3、宏病毒组研究可应用于人体或动物血液、组织、粪便等样本，植物组织样本，以及水体、土壤等各种环境样本，用以分析其中的病毒群落。然而，由于病毒基因组普遍较小，病毒核酸在样本中的相对含量非常低，宿主的基因组序列干扰严重，且已知病毒的数量和基因组信息有限等问题，使得宏病毒组研究数据分析存在着一定的困难。

4、目前的宏病毒组数据分析方法不够准确，系统性不强，较为繁琐，分析流程时间长，这阻碍了对病毒群落的全面了解，不利于测序数据的深入挖掘。

技术实现思路

1、为了解决目前存在的宏病毒组数据分析流程繁琐、准确性差、不系统的问题，本发明提供了一种宏病毒组数据的分析方法，包括：

2、步骤一：序列质控；

3、使用trimmomatic软件，去除低质量序列，使用bwa-men去除可能存在真核基因组的序列；

4、步骤二：序列组装和聚类；

5、使用megahit软件对所述步骤一获得的高质量序列进行组装，获得contigs序列，使用cd-hit对所述contigs序列进行去冗余并进行聚类；

6、步骤三：病毒序列识别和检查；

7、针对聚类后得到的contigs序列，使用virsorter2进行病毒序列的识别，使用check-v检查所述病毒序列的完整性及污染程度；

8、步骤四：病毒物种注释和丰度计算；

9、使用hmmer v3.1b2比对参考蛋白的hmm模型，根据蛋白比对结果进行物种注释，使用bwa软件计算病毒序列的丰度；

10、步骤五：病毒生活方式判断：使用phatyp软件进行病毒生活方式判断；

11、步骤六：病毒宿主预测：用crisprone搜索细菌序列中的crispr间隔、trna匹配、基因组同源性匹配三种方式对病毒宿主进行预测；

12、步骤七：病毒辅助代谢基因分析：使用vibrant和dram-v判定病毒编码的辅助代谢基因。

13、可选的，所述步骤六包括：

14、使用crisprone搜索细菌序列中的crispr间隔，使用blastn比对crispr间隔序列与病毒contig，选择同源性大于等于95％和单核苷酸多态性小于等于2个的对比结果，对应crispr间隔所在的细菌序列认定为假定的病毒宿主；

15、利用trnascan-se软件预测病毒contig中的trna基因，使用blastn对比trna序列与细菌scaffold，保留序列一致性大于等于95％的对比结果认定为病毒宿主；

16、通过blastn将细菌序列与病毒contig直接进行比对，依据参数：bitscore≥50，e-value≥0.003，identity≥70％，matching length≥2500bp筛选病毒的预测宿主。

17、可选的，所述步骤一使用trimmomatic软件从reads中剪切adapter和其他illumina特定序列；执行滑动窗口修剪，窗口大小为4bp，从5’端第一个碱基位置开始移动，窗口中碱基平均质量≥q15，从第一个平均质量低于q15的窗口的3’端碱基处截断序列；经上述质量筛查后，去除序列长度小于36bp的序列，去除低质量的序列。

18、可选的，所述步骤一使用bwa-men软件建立宿主基因组索引，然后将宿主基因组与样品reads进行比对，最后使用samtools去除宿主序列。

19、可选的，所述步骤二中使用megahit软件，设置最小的kmer为35，最大的kmer为95，kmer间隔为20，输出的contig最小值为500。

20、可选的，所述步骤二中使用cd-hit的est模式，相似性为0.9，两两序列进行对比时的word size为5，采用16gb ram，线程数为4。

21、可选的，所述步骤三中使用virsorter2，选择最小序列长度为1500bp，线程数为4。

22、可选的，所述步骤三中使用check-v软件的end_to_end模式，线程数40。

23、可选的，所述步骤四使用bwa软件将病毒contig文件与clean read进行对比，运行bwa的mem模式，kmer长度为30，线程数为40，输出文件使用pileup.sh进行覆盖率计算。

24、可选的，所述步骤四使用hmmer v3.1b2软件进行病毒物种注释，包括：

25、首先下载vogdb数据库及物种对应分类表，然后使用hummer下的程序hmmscan对比病毒contig与数据库的hmm模型，设置e-value为0.0001，cpu数目为40；

26、要求每条contig至少20％的基因有注释，然后在这20％的基因里至少60％来自同一个属，不符合条件的则认为是无法注释的。

27、本发明有益效果是：

28、本发明的病毒组数据分析方法，通过序列质控，使用trimmomatic软件，去除低质量序列，使用bwa-men去除可能存在真核基因组的序列，提高了病毒序列质量，减少了宿主基因组的干扰；使用megahit软件对所述步骤一获得的高质量序列进行组装，获得contigs序列，使用cd-hit对所述contigs序列进行去冗余并进行聚类；针对聚类后得到的contigs序列，使用virsorter2进行病毒序列的识别，使用check-v检查所述病毒序列的完整性及污染程度；使用hmmer v3.1b2比对参考蛋白的hmm模型，根据蛋白比对结果进行物种注释，使用bwa软件计算病毒序列的丰度；使用phatyp软件进行病毒生活方式判断；使用crisprone搜索细菌序列中的crispr间隔、trna匹配、基因组同源性匹配三种方式对病毒宿主进行预测；使用vibrant和dram-v判定病毒编码的辅助代谢基因。

29、本发明方法能够准确地对病毒物种进行鉴定和注释，全面系统的对宏病毒组数据进行深入解析和挖掘，步骤简洁明了，分析时间短，大幅优化了宏病毒鉴定研究效果。

技术特征：

1.一种宏病毒组数据的分析方法，其特征在于，所述方法包括：

2.根据权利要求1所述的宏病毒组数据的分析方法，其特征在于，所述步骤六包括：

3.根据权利要求1所述的宏病毒组数据的分析方法，其特征在于，所述步骤一使用trimmomatic软件从reads中剪切adapter和其他illumina特定序列；执行滑动窗口修剪，窗口大小为4bp，从5’端第一个碱基位置开始移动，窗口中碱基平均质量≥q15，从第一个平均质量低于q15的窗口的3’端碱基处截断序列；经上述质量筛查后，去除序列长度小于36bp的序列，去除低质量的序列。

4.根据权利要求1所述的宏病毒组数据的分析方法，其特征在于，所述步骤一使用bwa-men软件建立宿主基因组索引，然后将宿主基因组与样品reads进行比对，最后使用samtools去除宿主序列。

5.根据权利要求1所述的宏病毒组数据的分析方法，其特征在于，所述步骤二中使用megahit软件，设置最小的kmer为35，最大的kmer为95，kmer间隔为20，输出的contig最小值为500。

6.根据权利要求1所述的宏病毒组数据的分析方法，其特征在于，所述步骤二中使用cd-hit的est模式，相似性为0.9，两两序列进行对比时的word size为5，采用16gb ram，线程数为4。

7.根据权利要求1所述的宏病毒组数据的分析方法，其特征在于，所述步骤三中使用virsorter2，选择最小序列长度为1500bp，线程数为4。

8.根据权利要求1所述的宏病毒组数据的分析方法，其特征在于，所述步骤三中使用check-v软件的end_to_end模式，线程数40。

9.根据权利要求1所述的宏病毒组数据的分析方法，其特征在于，所述步骤四使用bwa软件将病毒contig文件与clean read进行对比，运行bwa的mem模式，kmer长度为30，线程数为40，输出文件使用pileup.sh进行覆盖率计算。

10.根据权利要求1所述的宏病毒组数据的分析方法，其特征在于，所述步骤四使用hmmer v3.1b2软件进行病毒物种注释，包括：

技术总结
本发明公开了一种宏病毒组数据的分析方法，属于宏病毒组学技术领域。所述方法包括：序列质控、序列组装、序列聚类、病毒序列识别、病毒序列检查、病毒丰度计算、物种注释、病毒生活方式判断、病毒宿主预测和病毒辅助代谢基因分析。本发明利用Trimmomatic软件、BWA‑MEN、Megahit软件、CD‑HIT等工具执行宏病毒组数据的分析流程，实践证明本发明能够准确地对病毒物种进行鉴定和注释，全面系统地对宏病毒组数据进行深入解析和挖掘，步骤简洁明了，分析时间短，大幅提高宏病毒鉴定研究效果。

技术研发人员：徐岩,张华东,张红霞
受保护的技术使用者：江南大学
技术研发日：
技术公布日：2024/1/15