一种用于宏基因组测序数据的微生物物种与功能组成分析方法[发明专利]

(19)中华人民共和国国家知识产权局
(12)发明专利申请
(10)申请公布号 (43)申请公布日 (21)申请号 202011592565.5
(22)申请日 2020.12.29
(71)申请人 上海派森诺生物科技股份有限公司
地址 200030 上海市徐汇区银都路218号2
号楼3、4层
(72)发明人 李鸿毅　曲昊淼　寇文伯　薛正晟　
孙子奎　
(74)专利代理机构 上海点威知识产权代理有限
公司 31326
代理人 胡志强
(51)Int.Cl.
G16B 30/00(2019.01)
G16B 50/00(2019.01)
(54)发明名称
一种用于宏基因组测序数据的微生物物种
与功能组成分析方法
(57)摘要
本发明公开了一种用于宏基因组测序数据
的微生物物种与功能组成分析方法，包括下列步
骤：1)切除原始数据中的接头序列片段、低质量
片段，滤除过短序列、含模糊碱基序列；2)使用上
述获得的数据进行物种注释，3)对剔除非目物种
后的序列进行拼接，4)对叠连群序列进行相似性
聚类，5)使用blastn算法，6)预测非冗余的叠连
群序列中的基因区域，7)将非冗余蛋白序列集与
各类蛋白注释数据库进行比对，8)计算基因序列
丰度，本发明满足了拼接产生的叠连群序列的进
一步分析，避免了结果假阳性高，依赖于专用数
据库的构建，适用面窄的问题；灵敏度高，提高了
准确性。

权利要求书2页 说明书8页 附图1页CN 112599198 A 2021.04.02
C N 112599198
A
1.一种用于宏基因组测序数据的微生物物种与功能组成分析方法，其特征在于，包括下列步骤：
1)切除原始数据中的接头序列片段、低质量片段，滤除过短序列、含模糊碱基序列；若已知宿主基因组，则将宿主序列剔除；
2)使用上述获得的数据进行物种注释，并统计物种序列数即为丰度，再基于注释结果剔除注释到非目的物种的序列；
3)对剔除非目物种后的序列进行拼接，获得叠连群序列；
4)对叠连群序列进行相似性聚类，并计算各样本非冗余的叠连群序列丰度，并去掉总丰度为零的序列；
5)使用blastn算法，对非冗余的叠连群序列进行核酸数据库比对，并采用共同祖先算法获取拼接序列的物种注释信息，再基于注释信息，将步骤2)中注释得到的物种分为已验证存在的物种与疑似存在物种；
6)预测非冗余的叠连群序列中的基因区域，获得基因序列及其翻译的蛋白序列，再对蛋白序列进行相似性聚类，获得非冗余蛋白序列集；
7)将非冗余蛋白序列集与各类蛋白注释数据库进行比对，获得基因序列和蛋白序列的功能注释信息；
8)计算基因序列丰度，再通过基因/蛋白对应功能信息，获得功能丰度表。

2.如权利要求1所述的一种用于宏基因组测序数据的微生物物种与功能组成分析方法，其特征在于：步骤1)中使用FastQC检查原始数据的测序质量情况；使用fastp或trimmomatic软件其中一种将原始数据中低质量片段切割，并滤除接头序列和过短序列，获得高质量序列；使用bowtie2或bmtagger其中一种软件，去除比对到宿主基因组上的序列。

3.如权利要求1所述的一种用于宏基因组测序数据的微生物物种与功能组成分析方法，其特征在于：步骤2)中进行物种注释为进行非拼接序列的物种注释，将质控后的序列基于物种注释数据库进行k‑mer检索或局部相似性比对，获取序列的物种注释信息以及物种丰度表一，其中使用kraken2及核酸序列数据库进行k‑mer检索或使用kaiju及蛋白序列数据库进行局部相似性比对；使用bracken计算物种组成丰度表；非目的序列默认为注释到后生动物及绿色植物的序列或者注释到自定义的物种的序列。

4.如权利要求1所述的一种用于宏基因组测序数据的微生物物种与功能组成分析方法，其特征在于：步骤3)中序列拼接为对剔除非目物种后的序列进行拼接，获得叠连群序列，并去冗余；所述拼接方式包括各样本单独拼接，再分别找出各样本无法比对上的序列，对无法比对上的序列进行混合拼接；或先对各样本按分组分别合并拼接，再分别找出各样本无法比对上的序列，对无法比对上的序列进行混合拼接；或对所有样本合并拼接；所述拼接用的拼接软件为megahit或metaspades；使用minimap2或bowtie2软件找出各样本无法比对上叠连群的序列。

5.如权利要求1所述的一种用于宏基因组测序数据的微生物物种与功能组成分析方法，其特征在于：所述叠连群核酸序列行相似性聚类采用的聚类软件为MMseqs2或cd‑hit，所述MMseqs2采用easy‑linclust模式；计算各样本非冗余的叠连群序列丰度采用基于比对的htseq‑count或不基于比对的salmon。

6.如权利要求1所述的一种用于宏基因组测序数据的微生物物种与功能组成分析方
法，其特征在于：所述步骤5)中比对用软件为minimap2或bowtie2。

7.如权利要求1所述的一种用于宏基因组测序数据的微生物物种与功能组成分析方法，其特征在于：所述步骤6)中基因预测软件为metagenemark、prodigal或FragGeneScan其中一种；蛋白序列聚类软件为cd‑hit或MMseqs2，所述MMseqs2的聚类模式为easy‑linclust 或easy‑cluster模式。

8.如权利要求1所述的一种用于宏基因组测序数据的微生物物种与功能组成分析方法，其特征在于：步骤7)中功能注释的蛋白或功能数据库包括eggNOG、KEGG、CAZy、GO、NCBI‑nr、swissprot；所述非冗余蛋白序列集与功能数据库的参考蛋白序列集进行比对使用使用MMseqs2或diamond软件。

9.如权利要求1所述的一种用于宏基因组测序数据的微生物物种与功能组成分析方法，其特征在于：所述步骤8)中计算基因序列丰度采用基于比对的htseq‑count或不基于比对的salmon，所述比对用软件为minimap2或bowtie2。

一种用于宏基因组测序数据的微生物物种与功能组成分析
方法
技术领域
[0001]本发明涉及微生物基因分析领域，尤其涉及一种用于宏基因组测序数据的微生物物种与功能组成分析方法。

背景技术
[0002]随着下一代高通量测序技术(Next Generation Sequencing，NGS)的不断发展，人们对于微生物群落方面的研究也越来越全面和深入。

区别于常见的靶向微生物核糖体RNA 基因的扩增子测序技术，宏基因组学是将整个群落系统中全体微生物的基因组作为研究对象，基于鸟枪法测序技术，全面展现整个群落的物种组成和功能潜能组成，进而阐明微生物群落的作用机制。

然而，由于测序样本类型的多样，样本量规模、测序深度的多变，以及宿主、污染基因组的多少等因素，以及分析本身的复杂性，研究者们发开出了数量庞大的各类软件，以及配套的更为复杂的各种分析模式、参数和数据库。

目前，只有少数几个流程软件或分析网站提供了自动化的宏基因组分析方法或服务。

其中，以MG‑RAST(https://www.mg‑/)和IMG/M(https:///cgi‑bin/m/main.cgi)等(宏)基因组数据库网站为代表，它们在提供数据上传和存储服务的同时，也提供了有限种类的分析项目和计算资源。

而另一类软件，如HUMAnN2(/humann2)以及发明专利“一种宏基因组数据分析方法及系统”(CN 108334750 B)等，则是基于非拼接序列数据的数据库比对注释，也就是说，其物种或基因的鉴定，主要依赖于数据库对目标样本的微生物基因组的覆盖度。

然而，对于环境样本，如土壤，底泥，极端环境，我们对于微生物群落的认识是极其有限的，因而此类分析流程，适用的样本类型也是有限，主要以人体相关样本为宜。

[0003]这些现有的宏基因组分析流程具有如下缺陷：
[0004](1)未包含拼接步骤，无法满足基于拼接产生的叠连群序列的进一步分析。

[0005](2)物种分类学注释依靠非拼接序列数据的比对分析，未包含序列拼接并基于拼接后的数据进行相应鉴定，因而结果假阳性高，依赖于专用数据库的构建，适用面有限。

[0006](3)对于样本中基因的鉴定，主要采用比对方式，依赖于专用数据库的构建，假阳性高，灵敏度低，适用范围有限。

[0007](4)未包含宿主和/或污染基因组序列去除步骤，或去除步骤依赖于宿主和/或污染基因组的存在。

发明内容
[0008]本发明的提供一种用于宏基因组测序数据的微生物物种与功能组成分析方法。

[0009]本发明的方案是：
[0010]一种用于宏基因组测序数据的微生物物种与功能组成分析方法，包括下列步骤：[0011]1)切除原始数据中的接头序列片段、低质量片段，滤除过短序列、含模糊碱基序
列；若已知宿主基因组，则将宿主序列剔除；
[0012]2)使用上述获得的数据进行物种注释，并统计物种序列数即为丰度，再基于注释结果剔除注释到非目的物种的序列；
[0013]3)对剔除非目物种后的序列进行拼接，获得叠连群序列；
[0014]4)对叠连群序列进行相似性聚类，并计算各样本非冗余的叠连群序列丰度，并去掉总丰度为零的序列；
[0015]5)使用blastn算法，对非冗余的叠连群序列进行核酸数据库比对，并采用共同祖先算法获取拼接序列的物种注释信息，再基于注释信息，将步骤2)中注释得到的物种分为已验证存在的物种与疑似存在物种；
[0016]6)预测非冗余的叠连群序列中的基因区域，获得基因序列及其翻译的蛋白序列，再对蛋白序列进行相似性聚类，获得非冗余蛋白序列集；
[0017]7)将非冗余蛋白序列集与各类蛋白注释数据库进行比对，获得基因序列和蛋白序列的功能注释信息；
[0018]8)计算基因序列丰度，再通过基因/蛋白对应功能信息，获得功能丰度表。

[0019]作为优选的技术方案，步骤1)中使用FastQC检查原始数据的测序质量情况；使用fastp或trimmomatic软件其中一种将原始数据中低质量片段切割，并滤除接头序列和过短序列，获得高质量序列；使用bowtie2或bmtagger其中一种软件，去除比对到宿主基因组上的序列。

[0020]序列质控：检查原始宏基因组数据的测序质量；使用移动划窗法将原始数据中低质量片段切除，并弃去过短序列；若已知宿主基因组，则将比对到宿主基因组上的序列剔除。

[0021]作为优选的技术方案，使用fastp软件将原始数据中低质量片段切割，并滤除接头序列和过短序列，获得高质量序列；
[0022]作为优选的技术方案，使用bmtagger软件去除比对到宿主基因组上的序列。

[0023]作为优选的技术方案，步骤2)中进行物种注释为进行非拼接序列的物种注释，将质控后的序列基于物种注释数据库进行k‑mer检索或局部相似性比对，获取序列的物种注释信息以及物种丰度表一，其中使用kraken2及核酸序列数据库进行k‑mer检索或使用kaiju及蛋白序列数据库进行局部相似性比对，对序列数据进行物种注释；使用bracken计算物种组成丰度表；非目的序列默认为注释到后生动物及绿色植物的序列或者注释到自定义的物种的序列。

[0024]非拼接序列物种注释：进行非拼接序列的物种注释，将质控后的序列基于物种注释数据库进行k‑mer检索或局部相似性比对，获取序列的物种注释信息以及物种丰度表一。

[0025]作为优选的技术方案，使用kraken2及核酸序列数据库，对质控后的序列数据进行物种注释，并依据注释结果剔除非目的序列，使用bracken计算物种组成丰度表。

[0026]作为优选的技术方案，步骤3)中序列拼接为对剔除非目物种后的序列进行拼接，获得叠连群序列，并去冗余；所述拼接方式包括各样本单独拼接，再分别找出各样本无法比对上的序列，对无法比对上的序列进行混合拼接；或先对各样本按分组分别合并拼接，再分别找出各样本无法比对上的序列，对无法比对上的序列进行混合拼接；或对所有样本合并拼接；所述拼接用的拼接软件为megahit或metaspades；使用minimap2或bowtie2软件找出
各样本无法比对上叠连群的序列。

[0027]序列拼接：对剔除非目物种后的序列进行拼接，获得叠连群序列，并去冗余。

[0028]作为优选的技术方案，所述叠连群核酸序列行相似性聚类采用的聚类软件为MMseqs2或cd‑hit，所述MMseqs2采用easy‑linclust模式；计算各样本非冗余的叠连群序列丰度采用基于比对的htseq‑count或不基于比对的salmon。

[0029]叠连群序列丰度计算：计算叠连群丰度，并去除丰度为0的序列，对余下叠连群序列进行相似性聚类，获得非冗余的叠连群序列的丰度表。

[0030]作为优选的技术方案，叠连群核酸序列聚类软件为MMseqs2。

[0031]作为优选的技术方案，计算叠连群丰度采用salmon。

[0032]作为优选的技术方案，所述步骤5)中比对用软件为minimap2或bowtie2。

[0033]叠连群序列物种注释：使用blastn算法，将非冗余的叠连群序列与NCBI‑nt数据库进行比对，采用最近祖先算法获取叠连群序列的物种信息；再结合非冗余的叠连群序列的丰度表，获得物种组成丰度表二；同时基于物种信息，将非拼接序列物种注释模块中注释得到的物种分为已验证存在的物种和疑似存在物种。

[0034]作为优选的技术方案，所述步骤6)中基因预测软件为metagenemark、prodigal或FragGeneScan其中一种；蛋白序列聚类软件为cd‑hit或MMseqs2，所述MMseqs2的聚类模式为easy‑linclust或easy‑cluster模式。

[0035]基因预测：预测非冗余的叠连群序列中的基因区域，获得基因序列及其翻译的蛋白序列，再对蛋白序列进行相似性聚类，获得非冗余蛋白序列集。

[0036]作为优选的技术方案，蛋白序列聚类软件为MMseqs2。

[0037]作为优选的技术方案，MMseqs2的聚类模式为easy‑cluster模式。

[0038]作为优选的技术方案，步骤7)中功能注释的蛋白或功能数据库包括eggNOG、KEGG、CAZy、GO、NCBI‑nr、swissprot；所述非冗余蛋白序列集与功能数据库的参考蛋白序列集进行比对使用使用MMseqs2或diamond软件。

[0039]功能注释：将非冗余蛋白序列集与蛋白或功能数据库进行比对，获得蛋白序列的功能注释信息，也即是基因序列的功能注释信息。

[0040]作为优选的技术方案，使用MMseqs2将非冗余蛋白序列集与功能数据库的参考蛋白序列集进行比对。

[0041]作为优选的技术方案，所述步骤8)中计算基因序列丰度采用基于比对的htseq‑count或不基于比对的salmon，所述比对用软件为minimap2或bowtie2。

[0042]基因丰度计算：计算基因序列丰度，也即是冗余蛋白序列的丰度，进而通过蛋白序列聚类的对应表，获得非冗余蛋白序列的丰度表，再结合这些蛋白的功能注释信息，最终获得功能单元的丰度表。

[0043]作为优选的技术方案，使用salmon计算基因序列丰度。

[0044]作为优选的技术方案，基于KEGG数据注释获得的KO丰度表，使用MinPath推断存在的KEGG代谢通路并计算其丰度表。

[0045]本发明提供了一种用于宏基因组测序数据的微生物物种与功能组成分析方法，包括下列步骤：1)切除原始数据中的接头序列片段、低质量片段，滤除过短序列、含模糊碱基序列；若已知宿主基因组，则将宿主序列剔除；2)使用上述获得的数据进行物种注释，并统
计物种序列数即为丰度，再基于注释结果剔除注释到非目的物种的序列；3)对剔除非目物种后的序列进行拼接，获得叠连群序列；4)对叠连群序列进行相似性聚类，并计算各样本非冗余的叠连群序列丰度，并去掉总丰度为零的序列；5)使用blastn算法，对非冗余的叠连群序列进行核酸数据库比对，并采用共同祖先算法获取拼接序列的物种注释信息，再基于注释信息，将步骤2)中注释得到的物种分为已验证存在的物种与疑似存在物种；6)预测非冗余的叠连群序列中的基因区域，获得基因序列及其翻译的蛋白序列，再对蛋白序列进行相似性聚类，获得非冗余蛋白序列集；7)将非冗余蛋白序列集与各类蛋白注释数据库进行比对，获得基因序列和蛋白序列的功能注释信息；8)计算基因序列丰度，再通过基因/蛋白对应功能信息，获得功能丰度表。

[0046]本发明的优点：
[0047]满足了拼接产生的叠连群序列的进一步分析，避免了结果假阳性高，依赖于专用数据库的构建，适用面窄的问题；灵敏度高，提高了准确性；
[0048]将复杂的宏基因组分析流程拆解总结为八大模块，再基于这些模块而非软件搭建流程，允许在特定模块中，选择调用不同的分析策略、软件或参数以应对各类情形。

[0049]对于物种的分类学鉴定和丰度定量，同时包含了：基于非拼接序列的比对或k‑mer 检索，用于快速鉴定物种并准确定量；基于拼接叠连群序列的blastn比对，用于严格判别丰度较高的物种或关注的物种的有无，降低物种鉴定假阳性。

[0050]对于功能(潜能)的鉴定和丰度定量，同时包含了：基于非拼接序列的比对或k‑mer 检索，用于快速鉴定物种并准确定量；基于拼接叠连群序列的blastn比对，用于严格判别丰度较高的物种或关注的物种的有无，降低物种鉴定假阳性。

[0051]输出的分析结果文件全面，使用者可依据不同需求进行相应选择，为后续个性化的统计分析提供可能。

附图说明
[0052]图1为本发明实施例2的流程示意图。

具体实施方式
[0053]为了弥补以上不足，本发明提供了一种用于宏基因组测序数据的微生物物种与功能组成分析方法以解决上述背景技术中的问题。

[0054]一种用于宏基因组测序数据的微生物物种与功能组成分析方法，包括下列步骤：[0055]1)切除原始数据中的接头序列片段、低质量片段，滤除过短序列、含模糊碱基序列；若已知宿主基因组，则将宿主序列剔除；
[0056]2)使用上述获得的数据进行物种注释，并统计物种序列数即为丰度，再基于注释结果剔除注释到非目的物种的序列；
[0057]3)对剔除非目物种后的序列进行拼接，获得叠连群序列；
[0058]4)对叠连群序列进行相似性聚类，并计算各样本非冗余的叠连群序列丰度，并去掉总丰度为零的序列；
[0059]5)使用blastn算法，对非冗余的叠连群序列进行核酸数据库比对，并采用共同祖先算法获取拼接序列的物种注释信息，再基于注释信息，将步骤2)中注释得到的物种分为
已验证存在的物种与疑似存在物种；
[0060]6)预测非冗余的叠连群序列中的基因区域，获得基因序列及其翻译的蛋白序列，再对蛋白序列进行相似性聚类，获得非冗余蛋白序列集；
[0061]7)将非冗余蛋白序列集与各类蛋白注释数据库进行比对，获得基因序列和蛋白序列的功能注释信息；
[0062]8)计算基因序列丰度，再通过基因/蛋白对应功能信息，获得功能丰度表。

[0063]步骤1)中使用FastQC检查原始数据的测序质量情况；使用fastp或trimmomatic软件其中一种将原始数据中低质量片段切割，并滤除接头序列和过短序列，获得高质量序列；使用bowtie2或bmtagger其中一种软件，去除比对到宿主基因组上的序列。

[0064]步骤2)中进行物种注释为进行非拼接序列的物种注释，将质控后的序列基于物种注释数据库进行k‑mer检索或局部相似性比对，获取序列的物种注释信息以及物种丰度表一，其中使用kraken2及核酸序列数据库进行k‑mer检索或使用kaiju及蛋白序列数据库进行局部相似性比对；使用bracken计算物种组成丰度表；非目的序列默认为注释到后生动物及绿色植物的序列或者注释到自定义的物种的序列。

[0065]步骤3)中序列拼接为对剔除非目物种后的序列进行拼接，获得叠连群序列，并去冗余；所述拼接方式包括各样本单独拼接，再分别找出各样本无法比对上的序列，对无法比对上的序列进行混合拼接；或先对各样本按分组分别合并拼接，再分别找出各样本无法比对上的序列，对无法比对上的序列进行混合拼接；或对所有样本合并拼接；所述拼接用的拼接软件为megahit或metaspades；使用minimap2或bowtie2软件找出各样本无法比对上叠连群的序列。

[0066]所述叠连群核酸序列行相似性聚类采用的聚类软件为MMseqs2或cd‑hit，所述MMseqs2采用easy‑linclust模式；计算各样本非冗余的叠连群序列丰度采用基于比对的htseq‑count或不基于比对的salmon。

[0067]所述步骤5)中比对用软件为minimap2或bowtie2。

[0068]所述步骤6)中基因预测软件为metagenemark、prodigal或FragGeneScan其中一种；蛋白序列聚类软件为cd‑hit或MMseqs2，所述MMseqs2的聚类模式为easy‑linclust或easy‑cluster模式。

[0069]步骤7)中功能注释的蛋白或功能数据库包括eggNOG、KEGG、CAZy、GO、NCBI‑nr、swissprot；所述非冗余蛋白序列集与功能数据库的参考蛋白序列集进行比对使用使用MMseqs2或diamond软件。

[0070]所述步骤8)中计算基因序列丰度采用基于比对的htseq‑count或不基于比对的salmon，所述比对用软件为minimap2或bowtie2。

[0071]为了使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合具体实施例，进一步阐述本发明。

[0072]实施例1：
[0073]一种用于宏基因组测序数据的微生物物种与功能组成分析方法，包括下列步骤：[0074]1)切除原始数据中的接头序列片段、低质量片段，滤除过短序列、含模糊碱基序列；若已知宿主基因组，则将宿主序列剔除；
[0075]2)使用上述获得的数据进行物种注释，并统计物种序列数即为丰度，再基于注释
结果剔除注释到非目的物种的序列；
[0076]3)对剔除非目物种后的序列进行拼接，获得叠连群序列；
[0077]4)对叠连群序列进行相似性聚类，并计算各样本非冗余的叠连群序列丰度，并去掉总丰度为零的序列；
[0078]5)使用blastn算法，对非冗余的叠连群序列进行核酸数据库比对，并采用共同祖先算法获取拼接序列的物种注释信息，再基于注释信息，将步骤2)中注释得到的物种分为已验证存在的物种与疑似存在物种；
[0079]6)预测非冗余的叠连群序列中的基因区域，获得基因序列及其翻译的蛋白序列，再对蛋白序列进行相似性聚类，获得非冗余蛋白序列集；
[0080]7)将非冗余蛋白序列集与各类蛋白注释数据库进行比对，获得基因序列和蛋白序列的功能注释信息；
[0081]8)计算基因序列丰度，再通过基因/蛋白对应功能信息，获得功能丰度表。

[0082]步骤1)中使用FastQC检查原始数据的测序质量情况；使用fastp或将原始数据中低质量片段切割，并滤除接头序列和过短序列，获得高质量序列；使用bowtie2软件，去除比对到宿主基因组上的序列。

[0083]步骤2)中进行物种注释为进行非拼接序列的物种注释，将质控后的序列基于物种注释数据库进行k‑mer检索或局部相似性比对，获取序列的物种注释信息以及物种丰度表一，其中使用kraken2及核酸序列数据库进行k‑mer检索；使用bracken计算物种组成丰度表；非目的序列默认为注释到后生动物及绿色植物的序列。

[0084]步骤3)中序列拼接为对剔除非目物种后的序列进行拼接，获得叠连群序列，并去冗余；所述拼接方式包括各样本单独拼接，再分别找出各样本无法比对上的序列，对无法比对上的序列进行混合拼接；所述拼接用的拼接软件为megahit；使用minimap软件找出各样本无法比对上叠连群的序列。

[0085]所述叠连群核酸序列行相似性聚类采用的聚类软件为MMseqs2，所述MMseqs2采用easy‑linclust模式；计算各样本非冗余的叠连群序列丰度采用不基于比对的salmon。

[0086]所述步骤5)中比对用软件为minimap2。

[0087]所述步骤6)中基因预测软件为metagenemark；蛋白序列聚类软件为MMseqs2，所述MMseqs2的聚类模式为easy‑linclust模式。

[0088]步骤7)中功能注释的蛋白或功能数据库包括eggNOG、KEGG、CAZy、GO、NCBI‑nr、swissprot；所述非冗余蛋白序列集与功能数据库的参考蛋白序列集进行比对使用使用MMseqs2软件。

[0089]所述步骤8)中计算基因序列丰度采用不基于比对的salmon，所述比对用软件为minimap2。

[0090]实施例2：
[0091]在本发明实施例中，使用模拟宏基因组数据进行分析，模拟数据物种组成见附表1，其中已知宿主基因组为人基因组。

[0092]在步骤S101中，首先使用FastQC检查原始数据的测序质量情况；使用Cutadapt识别3’端潜在的接头序列，并在识别的接头处截断。

这里。

要求与接头序列(R1:AGATCGGAAGA GCACACGTCTGAACTCCAGTCA；R2:AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT)的匹配长度至少达。