CN112599198A - 一种用于宏基因组测序数据的微生物物种与功能组成分析方法 - Google Patents

一种用于宏基因组测序数据的微生物物种与功能组成分析方法 Download PDF

Info

Publication number
CN112599198A
CN112599198A CN202011592565.5A CN202011592565A CN112599198A CN 112599198 A CN112599198 A CN 112599198A CN 202011592565 A CN202011592565 A CN 202011592565A CN 112599198 A CN112599198 A CN 112599198A
Authority
CN
China
Prior art keywords
sequences
species
sequence
splicing
abundance
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.)
Pending
Application number
CN202011592565.5A
Other languages
English (en)
Inventor
李鸿毅
曲昊淼
寇文伯
薛正晟
孙子奎
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shanghai Personal Biotechnology Co ltd
Original Assignee
Shanghai Personal Biotechnology Co ltd
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Shanghai Personal Biotechnology Co ltd filed Critical Shanghai Personal Biotechnology Co ltd
Priority to CN202011592565.5A priority Critical patent/CN112599198A/zh
Publication of CN112599198A publication Critical patent/CN112599198A/zh
Pending legal-status Critical Current

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
    • G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding
    • G16B40/30—Unsupervised data analysis
    • 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
    • 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
    • G16B50/00—ICT programming tools or database systems specially adapted for bioinformatics
    • 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
    • G16B50/00—ICT programming tools or database systems specially adapted for bioinformatics
    • G16B50/10—Ontologies; Annotations

Landscapes

  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • General Health & Medical Sciences (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Biophysics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biotechnology (AREA)
  • Evolutionary Biology (AREA)
  • Theoretical Computer Science (AREA)
  • Bioethics (AREA)
  • Databases & Information Systems (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Data Mining & Analysis (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Artificial Intelligence (AREA)
  • Epidemiology (AREA)
  • Evolutionary Computation (AREA)
  • Public Health (AREA)
  • Software Systems (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
  • Information Retrieval, Db Structures And Fs Structures Therefor (AREA)

Abstract

本发明公开了一种用于宏基因组测序数据的微生物物种与功能组成分析方法，包括下列步骤：1)切除原始数据中的接头序列片段、低质量片段，滤除过短序列、含模糊碱基序列；2)使用上述获得的数据进行物种注释，3)对剔除非目物种后的序列进行拼接，4)对叠连群序列进行相似性聚类，5)使用blastn算法，6)预测非冗余的叠连群序列中的基因区域，7)将非冗余蛋白序列集与各类蛋白注释数据库进行比对，8)计算基因序列丰度，本发明满足了拼接产生的叠连群序列的进一步分析，避免了结果假阳性高，依赖于专用数据库的构建，适用面窄的问题；灵敏度高，提高了准确性。

Description

一种用于宏基因组测序数据的微生物物种与功能组成分析 方法
技术领域
本发明涉及微生物基因分析领域，尤其涉及一种用于宏基因组测序数据的微生物物种与功能组成分析方法。
背景技术
随着下一代高通量测序技术(Next Generation Sequencing，NGS)的不断发展，人们对于微生物群落方面的研究也越来越全面和深入。区别于常见的靶向微生物核糖体RNA基因的扩增子测序技术，宏基因组学是将整个群落系统中全体微生物的基因组作为研究对象，基于鸟枪法测序技术，全面展现整个群落的物种组成和功能潜能组成，进而阐明微生物群落的作用机制。然而，由于测序样本类型的多样，样本量规模、测序深度的多变，以及宿主、污染基因组的多少等因素，以及分析本身的复杂性，研究者们发开出了数量庞大的各类软件，以及配套的更为复杂的各种分析模式、参数和数据库。目前，只有少数几个流程软件或分析网站提供了自动化的宏基因组分析方法或服务。其中，以MG-RAST(https://www.mg-rast.org/)和IMG/M(https://img.jgi.doe.gov/cgi-bin/m/main.cgi)等(宏)基因组数据库网站为代表，它们在提供数据上传和存储服务的同时，也提供了有限种类的分析项目和计算资源。而另一类软件，如HUMAnN2(http://www.huttenhower.org/humann2)以及发明专利“一种宏基因组数据分析方法及系统”(CN 108334750 B)等，则是基于非拼接序列数据的数据库比对注释，也就是说，其物种或基因的鉴定，主要依赖于数据库对目标样本的微生物基因组的覆盖度。然而，对于环境样本，如土壤，底泥，极端环境，我们对于微生物群落的认识是极其有限的，因而此类分析流程，适用的样本类型也是有限，主要以人体相关样本为宜。
这些现有的宏基因组分析流程具有如下缺陷：
(1)未包含拼接步骤，无法满足基于拼接产生的叠连群序列的进一步分析。
(2)物种分类学注释依靠非拼接序列数据的比对分析，未包含序列拼接并基于拼接后的数据进行相应鉴定，因而结果假阳性高，依赖于专用数据库的构建，适用面有限。
(3)对于样本中基因的鉴定，主要采用比对方式，依赖于专用数据库的构建，假阳性高，灵敏度低，适用范围有限。
(4)未包含宿主和/或污染基因组序列去除步骤，或去除步骤依赖于宿主和/或污染基因组的存在。
发明内容
本发明的提供一种用于宏基因组测序数据的微生物物种与功能组成分析方法。
本发明的方案是：
一种用于宏基因组测序数据的微生物物种与功能组成分析方法，包括下列步骤：
1)切除原始数据中的接头序列片段、低质量片段，滤除过短序列、含模糊碱基序列；若已知宿主基因组，则将宿主序列剔除；
2)使用上述获得的数据进行物种注释，并统计物种序列数即为丰度，再基于注释结果剔除注释到非目的物种的序列；
3)对剔除非目物种后的序列进行拼接，获得叠连群序列；
4)对叠连群序列进行相似性聚类，并计算各样本非冗余的叠连群序列丰度，并去掉总丰度为零的序列；
5)使用blastn算法，对非冗余的叠连群序列进行核酸数据库比对，并采用共同祖先算法获取拼接序列的物种注释信息，再基于注释信息，将步骤2)中注释得到的物种分为已验证存在的物种与疑似存在物种；
6)预测非冗余的叠连群序列中的基因区域，获得基因序列及其翻译的蛋白序列，再对蛋白序列进行相似性聚类，获得非冗余蛋白序列集；
7)将非冗余蛋白序列集与各类蛋白注释数据库进行比对，获得基因序列和蛋白序列的功能注释信息；
8)计算基因序列丰度，再通过基因/蛋白对应功能信息，获得功能丰度表。
作为优选的技术方案，步骤1)中使用FastQC检查原始数据的测序质量情况；使用fastp或trimmomatic软件其中一种将原始数据中低质量片段切割，并滤除接头序列和过短序列，获得高质量序列；使用bowtie2或bmtagger其中一种软件，去除比对到宿主基因组上的序列。
序列质控：检查原始宏基因组数据的测序质量；使用移动划窗法将原始数据中低质量片段切除，并弃去过短序列；若已知宿主基因组，则将比对到宿主基因组上的序列剔除。
作为优选的技术方案，使用fastp软件将原始数据中低质量片段切割，并滤除接头序列和过短序列，获得高质量序列；
作为优选的技术方案，使用bmtagger软件去除比对到宿主基因组上的序列。
作为优选的技术方案，步骤2)中进行物种注释为进行非拼接序列的物种注释，将质控后的序列基于物种注释数据库进行k-mer检索或局部相似性比对，获取序列的物种注释信息以及物种丰度表一，其中使用kraken2及核酸序列数据库进行k-mer检索或使用kaiju及蛋白序列数据库进行局部相似性比对，对序列数据进行物种注释；使用bracken计算物种组成丰度表；非目的序列默认为注释到后生动物及绿色植物的序列或者注释到自定义的物种的序列。
非拼接序列物种注释：进行非拼接序列的物种注释，将质控后的序列基于物种注释数据库进行k-mer检索或局部相似性比对，获取序列的物种注释信息以及物种丰度表一。
作为优选的技术方案，使用kraken2及核酸序列数据库，对质控后的序列数据进行物种注释，并依据注释结果剔除非目的序列，使用bracken计算物种组成丰度表。
作为优选的技术方案，步骤3)中序列拼接为对剔除非目物种后的序列进行拼接，获得叠连群序列，并去冗余；所述拼接方式包括各样本单独拼接，再分别找出各样本无法比对上的序列，对无法比对上的序列进行混合拼接；或先对各样本按分组分别合并拼接，再分别找出各样本无法比对上的序列，对无法比对上的序列进行混合拼接；或对所有样本合并拼接；所述拼接用的拼接软件为megahit或metaspades；使用minimap2或bowtie2软件找出各样本无法比对上叠连群的序列。
序列拼接：对剔除非目物种后的序列进行拼接，获得叠连群序列，并去冗余。
作为优选的技术方案，所述叠连群核酸序列行相似性聚类采用的聚类软件为MMseqs2或cd-hit，所述MMseqs2采用easy-linclust模式；计算各样本非冗余的叠连群序列丰度采用基于比对的htseq-count或不基于比对的salmon。
叠连群序列丰度计算：计算叠连群丰度，并去除丰度为0的序列，对余下叠连群序列进行相似性聚类，获得非冗余的叠连群序列的丰度表。
作为优选的技术方案，叠连群核酸序列聚类软件为MMseqs2。
作为优选的技术方案，计算叠连群丰度采用salmon。
作为优选的技术方案，所述步骤5)中比对用软件为minimap2或bowtie2。
叠连群序列物种注释：使用blastn算法，将非冗余的叠连群序列与NCBI-nt数据库进行比对，采用最近祖先算法获取叠连群序列的物种信息；再结合非冗余的叠连群序列的丰度表，获得物种组成丰度表二；同时基于物种信息，将非拼接序列物种注释模块中注释得到的物种分为已验证存在的物种和疑似存在物种。
作为优选的技术方案，所述步骤6)中基因预测软件为metagenemark、prodigal或FragGeneScan其中一种；蛋白序列聚类软件为cd-hit或MMseqs2，所述MMseqs2的聚类模式为easy-linclust或easy-cluster模式。
基因预测：预测非冗余的叠连群序列中的基因区域，获得基因序列及其翻译的蛋白序列，再对蛋白序列进行相似性聚类，获得非冗余蛋白序列集。
作为优选的技术方案，蛋白序列聚类软件为MMseqs2。
作为优选的技术方案，MMseqs2的聚类模式为easy-cluster模式。
作为优选的技术方案，步骤7)中功能注释的蛋白或功能数据库包括eggNOG、KEGG、CAZy、GO、NCBI-nr、swissprot；所述非冗余蛋白序列集与功能数据库的参考蛋白序列集进行比对使用使用MMseqs2或diamond软件。
功能注释：将非冗余蛋白序列集与蛋白或功能数据库进行比对，获得蛋白序列的功能注释信息，也即是基因序列的功能注释信息。
作为优选的技术方案，使用MMseqs2将非冗余蛋白序列集与功能数据库的参考蛋白序列集进行比对。
作为优选的技术方案，所述步骤8)中计算基因序列丰度采用基于比对的htseq-count或不基于比对的salmon，所述比对用软件为minimap2或bowtie2。
基因丰度计算：计算基因序列丰度，也即是冗余蛋白序列的丰度，进而通过蛋白序列聚类的对应表，获得非冗余蛋白序列的丰度表，再结合这些蛋白的功能注释信息，最终获得功能单元的丰度表。
作为优选的技术方案，使用salmon计算基因序列丰度。
作为优选的技术方案，基于KEGG数据注释获得的KO丰度表，使用MinPath推断存在的KEGG代谢通路并计算其丰度表。
本发明提供了一种用于宏基因组测序数据的微生物物种与功能组成分析方法，包括下列步骤：1)切除原始数据中的接头序列片段、低质量片段，滤除过短序列、含模糊碱基序列；若已知宿主基因组，则将宿主序列剔除；2)使用上述获得的数据进行物种注释，并统计物种序列数即为丰度，再基于注释结果剔除注释到非目的物种的序列；3)对剔除非目物种后的序列进行拼接，获得叠连群序列；4)对叠连群序列进行相似性聚类，并计算各样本非冗余的叠连群序列丰度，并去掉总丰度为零的序列；5)使用blastn算法，对非冗余的叠连群序列进行核酸数据库比对，并采用共同祖先算法获取拼接序列的物种注释信息，再基于注释信息，将步骤2)中注释得到的物种分为已验证存在的物种与疑似存在物种；6)预测非冗余的叠连群序列中的基因区域，获得基因序列及其翻译的蛋白序列，再对蛋白序列进行相似性聚类，获得非冗余蛋白序列集；7)将非冗余蛋白序列集与各类蛋白注释数据库进行比对，获得基因序列和蛋白序列的功能注释信息；8)计算基因序列丰度，再通过基因/蛋白对应功能信息，获得功能丰度表。
本发明的优点：
满足了拼接产生的叠连群序列的进一步分析，避免了结果假阳性高，依赖于专用数据库的构建，适用面窄的问题；灵敏度高，提高了准确性；
将复杂的宏基因组分析流程拆解总结为八大模块，再基于这些模块而非软件搭建流程，允许在特定模块中，选择调用不同的分析策略、软件或参数以应对各类情形。
对于物种的分类学鉴定和丰度定量，同时包含了：基于非拼接序列的比对或k-mer检索，用于快速鉴定物种并准确定量；基于拼接叠连群序列的blastn比对，用于严格判别丰度较高的物种或关注的物种的有无，降低物种鉴定假阳性。
对于功能(潜能)的鉴定和丰度定量，同时包含了：基于非拼接序列的比对或k-mer检索，用于快速鉴定物种并准确定量；基于拼接叠连群序列的blastn比对，用于严格判别丰度较高的物种或关注的物种的有无，降低物种鉴定假阳性。
输出的分析结果文件全面，使用者可依据不同需求进行相应选择，为后续个性化的统计分析提供可能。
附图说明
图1为本发明实施例2的流程示意图。
具体实施方式
为了弥补以上不足，本发明提供了一种用于宏基因组测序数据的微生物物种与功能组成分析方法以解决上述背景技术中的问题。
一种用于宏基因组测序数据的微生物物种与功能组成分析方法，包括下列步骤：
1)切除原始数据中的接头序列片段、低质量片段，滤除过短序列、含模糊碱基序列；若已知宿主基因组，则将宿主序列剔除；
2)使用上述获得的数据进行物种注释，并统计物种序列数即为丰度，再基于注释结果剔除注释到非目的物种的序列；
3)对剔除非目物种后的序列进行拼接，获得叠连群序列；
4)对叠连群序列进行相似性聚类，并计算各样本非冗余的叠连群序列丰度，并去掉总丰度为零的序列；
5)使用blastn算法，对非冗余的叠连群序列进行核酸数据库比对，并采用共同祖先算法获取拼接序列的物种注释信息，再基于注释信息，将步骤2)中注释得到的物种分为已验证存在的物种与疑似存在物种；
6)预测非冗余的叠连群序列中的基因区域，获得基因序列及其翻译的蛋白序列，再对蛋白序列进行相似性聚类，获得非冗余蛋白序列集；
7)将非冗余蛋白序列集与各类蛋白注释数据库进行比对，获得基因序列和蛋白序列的功能注释信息；
8)计算基因序列丰度，再通过基因/蛋白对应功能信息，获得功能丰度表。
步骤1)中使用FastQC检查原始数据的测序质量情况；使用fastp或trimmomatic软件其中一种将原始数据中低质量片段切割，并滤除接头序列和过短序列，获得高质量序列；使用bowtie2或bmtagger其中一种软件，去除比对到宿主基因组上的序列。
步骤2)中进行物种注释为进行非拼接序列的物种注释，将质控后的序列基于物种注释数据库进行k-mer检索或局部相似性比对，获取序列的物种注释信息以及物种丰度表一，其中使用kraken2及核酸序列数据库进行k-mer检索或使用kaiju及蛋白序列数据库进行局部相似性比对；使用bracken计算物种组成丰度表；非目的序列默认为注释到后生动物及绿色植物的序列或者注释到自定义的物种的序列。
步骤3)中序列拼接为对剔除非目物种后的序列进行拼接，获得叠连群序列，并去冗余；所述拼接方式包括各样本单独拼接，再分别找出各样本无法比对上的序列，对无法比对上的序列进行混合拼接；或先对各样本按分组分别合并拼接，再分别找出各样本无法比对上的序列，对无法比对上的序列进行混合拼接；或对所有样本合并拼接；所述拼接用的拼接软件为megahit或metaspades；使用minimap2或bowtie2软件找出各样本无法比对上叠连群的序列。
所述叠连群核酸序列行相似性聚类采用的聚类软件为MMseqs2或cd-hit，所述MMseqs2采用easy-linclust模式；计算各样本非冗余的叠连群序列丰度采用基于比对的htseq-count或不基于比对的salmon。
所述步骤5)中比对用软件为minimap2或bowtie2。
所述步骤6)中基因预测软件为metagenemark、prodigal或FragGeneScan其中一种；蛋白序列聚类软件为cd-hit或MMseqs2，所述MMseqs2的聚类模式为easy-linclust或easy-cluster模式。
步骤7)中功能注释的蛋白或功能数据库包括eggNOG、KEGG、CAZy、GO、NCBI-nr、swissprot；所述非冗余蛋白序列集与功能数据库的参考蛋白序列集进行比对使用使用MMseqs2或diamond软件。
所述步骤8)中计算基因序列丰度采用基于比对的htseq-count或不基于比对的salmon，所述比对用软件为minimap2或bowtie2。
为了使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合具体实施例，进一步阐述本发明。
实施例1：
一种用于宏基因组测序数据的微生物物种与功能组成分析方法，包括下列步骤：
1)切除原始数据中的接头序列片段、低质量片段，滤除过短序列、含模糊碱基序列；若已知宿主基因组，则将宿主序列剔除；
2)使用上述获得的数据进行物种注释，并统计物种序列数即为丰度，再基于注释结果剔除注释到非目的物种的序列；
3)对剔除非目物种后的序列进行拼接，获得叠连群序列；
4)对叠连群序列进行相似性聚类，并计算各样本非冗余的叠连群序列丰度，并去掉总丰度为零的序列；
5)使用blastn算法，对非冗余的叠连群序列进行核酸数据库比对，并采用共同祖先算法获取拼接序列的物种注释信息，再基于注释信息，将步骤2)中注释得到的物种分为已验证存在的物种与疑似存在物种；
6)预测非冗余的叠连群序列中的基因区域，获得基因序列及其翻译的蛋白序列，再对蛋白序列进行相似性聚类，获得非冗余蛋白序列集；
7)将非冗余蛋白序列集与各类蛋白注释数据库进行比对，获得基因序列和蛋白序列的功能注释信息；
8)计算基因序列丰度，再通过基因/蛋白对应功能信息，获得功能丰度表。
步骤1)中使用FastQC检查原始数据的测序质量情况；使用fastp或将原始数据中低质量片段切割，并滤除接头序列和过短序列，获得高质量序列；使用bowtie2软件，去除比对到宿主基因组上的序列。
步骤2)中进行物种注释为进行非拼接序列的物种注释，将质控后的序列基于物种注释数据库进行k-mer检索或局部相似性比对，获取序列的物种注释信息以及物种丰度表一，其中使用kraken2及核酸序列数据库进行k-mer检索；使用bracken计算物种组成丰度表；非目的序列默认为注释到后生动物及绿色植物的序列。
步骤3)中序列拼接为对剔除非目物种后的序列进行拼接，获得叠连群序列，并去冗余；所述拼接方式包括各样本单独拼接，再分别找出各样本无法比对上的序列，对无法比对上的序列进行混合拼接；所述拼接用的拼接软件为megahit；使用minimap软件找出各样本无法比对上叠连群的序列。
所述叠连群核酸序列行相似性聚类采用的聚类软件为MMseqs2，所述MMseqs2采用easy-linclust模式；计算各样本非冗余的叠连群序列丰度采用不基于比对的salmon。
所述步骤5)中比对用软件为minimap2。
所述步骤6)中基因预测软件为metagenemark；蛋白序列聚类软件为MMseqs2，所述MMseqs2的聚类模式为easy-linclust模式。
步骤7)中功能注释的蛋白或功能数据库包括eggNOG、KEGG、CAZy、GO、NCBI-nr、swissprot；所述非冗余蛋白序列集与功能数据库的参考蛋白序列集进行比对使用使用MMseqs2软件。
所述步骤8)中计算基因序列丰度采用不基于比对的salmon，所述比对用软件为minimap2。
实施例2：
在本发明实施例中，使用模拟宏基因组数据进行分析，模拟数据物种组成见附表1，其中已知宿主基因组为人基因组。
在步骤S101中，首先使用FastQC检查原始数据的测序质量情况；使用Cutadapt识别3’端潜在的接头序列，并在识别的接头处截断。这里。要求与接头序列(R1:AGATCGGAAGAGCACACGTCTGAACTCCAGTCA；R2:AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT)的匹配长度至少达到3bp，且允许至多20％的碱基错配率。再使用fastp软件将其中的低质量片段切割，具体为采用滑动窗口法对序列进行质量筛查：窗口大小为5bp，从5'端第一个碱基位置开始移动，要求窗口中碱基平均质量大于等于Q20(即碱基平均测序准确率大于99％)，从第一个平均质量值低于Q20的窗口的3'端碱基处截断序列；经上述质量筛查后，去除序列长度小于50bp的序列，以及含有模糊碱基的序列。使用bmtagger软件去除比对到宿主——人基因组上的序列。
在步骤S102中，进行非拼接序列的物种注释，将步骤S101获得的有效序列基于物种注释数据库进行搜索或比对，获取序列的物种注释信息。这里使用kraken2及核酸序列数据库，对质控后的序列数据进行物种注释，并依据注释结果剔除非目的序列，使用bracken计算物种组成丰度表。非目的序列设置为注释到后生动物及绿色植物的序列。使用Kraken2进行物种注释。Kraken2通过使用精准k-mer匹配来实现高精度且快速的物种注释。数据库使用NCBI的RefSeq基因组数据库中的细菌、古菌、真菌、原生动物、病毒基因组，以及后生动物和绿色植物基因组数据作为参考构建数据库，调整置信度为0.5，进行注释分析。基于注释结果，再剔除注释到非目的物种的序列，这里的非目的物种为Metazoa(后生动物)和Viridiplantae(绿色植物)。最后，使用Bracken计算物种组成丰度表
在步骤S103中，先采用megahit的meta-sensitive模式，先对各样本单独拼接，再采用minimap2软件分别找出各样本无法比对上的序列，对这些序列再采用megahit的meta-sensitive模式进行混合拼接。
在步骤S104中，先使用软件MMseqs2的easy-linclust模式将(合并的)叠连群序列集按相似度95％，对齐区域覆盖度90％(占短序列的比例)进行聚类去冗余，获得非冗余的叠连群集合。丰度计算采用salmon软件，采用k-mer检索算法，计算冗余叠连群序列丰度，并依据聚类的对应关系，获得非冗余的叠连群序列的丰度表。
在步骤S105中，使用blastn，将非冗余的叠连群序列与ncbi-nt数据库比对，E值设置为0.00001，我们保留比对结果中E值最小的5个结果，采用Blast2lca软件中的“最近共同祖先(Lowest Common Ancestor,LCA)”算法，将参考序列分化为不同物种分枝前的最后一级共同分类，作为目标序列的物种分类注释信息。并进一步结合非冗余的叠连群序列的丰度表，获得物种组成丰度表二。基于该信息，再将步骤二中注释得到的物种分为已验证存在的物种和疑似存在物种。
在步骤S106中，进行基因预测。预测非冗余的叠连群序列中的基因区域，获得基因序列及其翻译的蛋白序列，再对蛋白序列进行相似性聚类，获得非冗余蛋白序列集。对非冗余的叠连群序列采用MetaGeneMark软件识别其中的开放阅读框，并预测其中的编码区域，从而获得对应的基因序列文件，蛋白序列文件，GTF(gene transfer format)文件以及GFF(general feature format)文件。再使用软件MMseqs2的easy-cluster模式将(合并的)蛋白序列集按相似度95％，对齐区域覆盖度90％(占短序列的比例)进行去冗余，获得非冗余的蛋白序列集。
在步骤S107中，采用MMseqs2的search模式，将非冗余蛋白序列集与eggnog数据库蛋白参考序列进行比对，获得蛋白序列的eggNOG Orthologous Groups功能注释信息。基于eggNOG注释信息，使用eggnog-mapper，推定蛋白序列的GO的词条，获得其GO功能注释信息。采用MMseqs2的search模式，将非冗余蛋白序列集与KOBAS软件的蛋白数据库进行比对，获得蛋白序列对应的KO信息。采用MMseqs2的search模式，将非冗余蛋白序列集与KOBAS软件的蛋白数据库进行比对，获得蛋白序列对应的KO信息。采用MMseqs2的search模式，将非冗余蛋白序列集与CAZy数据库的蛋白序列进行比对，获得蛋白序列的碳水化合物活性酶注释信息。采用MMseqs2的search模式，将非冗余蛋白序列集与NCBI-nr和swissprot数据库进行比对，获得蛋白序列的多种功能注释信息。
在步骤S108中，使用salmon计算基因序列丰度，也即是冗余的蛋白序列丰度；再通过聚类的对应表，获得非冗余蛋白序列的丰度表，结合S107得到的蛋白的功能信息，获得功能丰度表。基于各样本KO的丰度，使用MinPath推断各样本存在的KEGG代谢通路并计算其丰度。
自此，获得了基于分类水平注释的物种丰度表，物种分为已验证存在的物种和疑似存在物种。同时，也获得了eggNOG、KEGG、CAZy、GO等功能数据库注释的功能丰度表。
附表1如下：
Figure BDA0002869065890000111
以上显示和描述了本发明的基本原理、主要特征及本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (9)

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。
CN202011592565.5A 2020-12-29 2020-12-29 一种用于宏基因组测序数据的微生物物种与功能组成分析方法 Pending CN112599198A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011592565.5A CN112599198A (zh) 2020-12-29 2020-12-29 一种用于宏基因组测序数据的微生物物种与功能组成分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011592565.5A CN112599198A (zh) 2020-12-29 2020-12-29 一种用于宏基因组测序数据的微生物物种与功能组成分析方法

Publications (1)

Publication Number Publication Date
CN112599198A true CN112599198A (zh) 2021-04-02

Family

ID=75203360

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011592565.5A Pending CN112599198A (zh) 2020-12-29 2020-12-29 一种用于宏基因组测序数据的微生物物种与功能组成分析方法

Country Status (1)

Country Link
CN (1) CN112599198A (zh)

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113539369A (zh) * 2021-07-14 2021-10-22 江苏先声医学诊断有限公司 一种优化的kraken2算法及其在二代测序中的应用
CN113611359A (zh) * 2021-08-13 2021-11-05 江苏先声医学诊断有限公司 一种提高宏基因组纳米孔测序数据菌种组装效率的方法
CN114783519A (zh) * 2022-05-30 2022-07-22 南京农业大学 一种利用宏基因组分析土壤生物复合污染的方法
CN115044689A (zh) * 2022-06-17 2022-09-13 北京大学第三医院（北京大学第三临床医学院） 一种呼吸机相关性肺炎下呼吸道抽吸物的样本处理方法
WO2023098152A1 (zh) * 2021-11-30 2023-06-08 深圳零一生命科技有限责任公司 一种微生物基因数据库的构建方法及系统
CN116312798A (zh) * 2023-02-22 2023-06-23 江苏先声医学诊断有限公司 一种宏基因组测序数据物种验证的方法及应用
WO2024066461A1 (zh) * 2022-09-26 2024-04-04 华东理工大学 基于宏基因组和宏转录组识别油藏驱油功能微生物的方法
CN118212987A (zh) * 2024-05-21 2024-06-18 中国医学科学院北京协和医院 一种基因数据处理方法、装置、存储介质及电子设备
WO2024138840A1 (zh) * 2022-12-30 2024-07-04 中国农业科学院深圳农业基因组研究所 用于宏基因组物种成分及丰度定量的数据处理方法、装置及存储介质
WO2025050635A1 (zh) * 2023-09-08 2025-03-13 上海派森诺生物科技股份有限公司 一种宏基因组功能注释校正模型的构建方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110033826A (zh) * 2018-12-10 2019-07-19 上海派森诺生物科技股份有限公司 一种应用于宏病毒组高通量测序数据的分析方法
CN110349629A (zh) * 2019-06-20 2019-10-18 广州赛哲生物科技股份有限公司 一种利用宏基因组或宏转录组检测微生物的分析方法
US20200160936A1 (en) * 2017-06-28 2020-05-21 Icahn School Of Medicine At Mount Sinai Methods for high-resolution microbiome analysis
CN112071366A (zh) * 2020-10-13 2020-12-11 南开大学 一种基于二代测序技术的宏基因组数据分析方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20200160936A1 (en) * 2017-06-28 2020-05-21 Icahn School Of Medicine At Mount Sinai Methods for high-resolution microbiome analysis
CN110033826A (zh) * 2018-12-10 2019-07-19 上海派森诺生物科技股份有限公司 一种应用于宏病毒组高通量测序数据的分析方法
CN110349629A (zh) * 2019-06-20 2019-10-18 广州赛哲生物科技股份有限公司 一种利用宏基因组或宏转录组检测微生物的分析方法
CN112071366A (zh) * 2020-10-13 2020-12-11 南开大学 一种基于二代测序技术的宏基因组数据分析方法

Cited By (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113539369A (zh) * 2021-07-14 2021-10-22 江苏先声医学诊断有限公司 一种优化的kraken2算法及其在二代测序中的应用
CN113611359A (zh) * 2021-08-13 2021-11-05 江苏先声医学诊断有限公司 一种提高宏基因组纳米孔测序数据菌种组装效率的方法
CN113611359B (zh) * 2021-08-13 2022-08-05 江苏先声医学诊断有限公司 一种提高宏基因组纳米孔测序数据菌种组装效率的方法
WO2023098152A1 (zh) * 2021-11-30 2023-06-08 深圳零一生命科技有限责任公司 一种微生物基因数据库的构建方法及系统
CN114783519A (zh) * 2022-05-30 2022-07-22 南京农业大学 一种利用宏基因组分析土壤生物复合污染的方法
CN114783519B (zh) * 2022-05-30 2025-03-25 南京农业大学 一种利用宏基因组分析土壤生物复合污染的方法
CN115044689A (zh) * 2022-06-17 2022-09-13 北京大学第三医院（北京大学第三临床医学院） 一种呼吸机相关性肺炎下呼吸道抽吸物的样本处理方法
WO2024066461A1 (zh) * 2022-09-26 2024-04-04 华东理工大学 基于宏基因组和宏转录组识别油藏驱油功能微生物的方法
WO2024138840A1 (zh) * 2022-12-30 2024-07-04 中国农业科学院深圳农业基因组研究所 用于宏基因组物种成分及丰度定量的数据处理方法、装置及存储介质
CN116312798A (zh) * 2023-02-22 2023-06-23 江苏先声医学诊断有限公司 一种宏基因组测序数据物种验证的方法及应用
CN116312798B (zh) * 2023-02-22 2023-11-10 江苏先声医学诊断有限公司 一种宏基因组测序数据物种验证的方法及应用
WO2025050635A1 (zh) * 2023-09-08 2025-03-13 上海派森诺生物科技股份有限公司 一种宏基因组功能注释校正模型的构建方法
CN118212987A (zh) * 2024-05-21 2024-06-18 中国医学科学院北京协和医院 一种基因数据处理方法、装置、存储介质及电子设备

Similar Documents

Publication Publication Date Title
CN112599198A (zh) 一种用于宏基因组测序数据的微生物物种与功能组成分析方法
Feng et al. Metagenome assembly of high-fidelity long reads with hifiasm-meta
Srivastava et al. Alevin efficiently estimates accurate gene abundances from dscRNA-seq data
Siegwald et al. Assessment of common and emerging bioinformatics pipelines for targeted metagenomics
Saheb Kashaf et al. Recovering prokaryotic genomes from host-associated, short-read shotgun metagenomic sequencing data
US9189594B2 (en) Method and systems for processing polymeric sequence data and related information
US20250182850A1 (en) Creation or use of anchor-based data structures for sample-derived characteristic determination
CN111009286A (zh) 对宿主样本进行微生物分析的方法和装置
CN102521528A (zh) 一种基因序列数据的筛选方法
CN115631789B (zh) 一种基于泛基因组的群体联合变异检测方法
CN112687344B (zh) 一种基于宏基因组的人腺病毒分子分型和溯源方法及系统
CN115064215B (zh) 一种通过相似度进行菌株溯源及属性鉴定的方法
CN111192630B (zh) 一种宏基因组数据挖掘方法
Chung et al. FADU: a quantification tool for prokaryotic transcriptomic analyses
Hesse K-Mer-based genome size estimation in theory and practice
Wu et al. DeepRetention: a deep learning approach for intron retention detection
Chan et al. The artefactual branch effect and phylogenetic conflict: species delimitation with gene flow in mangrove pit vipers (Trimeresurus purpureomaculatus-erythrurus complex)
Darvish et al. Needle: a fast and space-efficient prefilter for estimating the quantification of very large collections of expression experiments
CN118212987B (zh) 一种基因数据处理方法、装置、存储介质及电子设备
CN110111847A (zh) 基于its2鉴定植物物种的方法及设备
CN110970093B (zh) 一种筛选引物设计模板的方法、装置及应用
JP7649822B2 (ja) ナノポア及びハイスループット配列決定データに基づくファージ同定及び宿主予測方法
Fu et al. VIGA: a one-stop tool for eukaryotic virus identification and genome assembly from next-generation-sequencing data
WO2024187428A1 (zh) 基于stLFR宏基因组测序数据构建高质量微生物基因组的组装流程
CN115862740B (zh) 一种面向大规模病毒基因组数据的快速分布式多序列比对方法

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
RJ01 Rejection of invention patent application after publication
RJ01 Rejection of invention patent application after publication

Application publication date: 20210402