一种宏转录组测序数据处理方法及系统

摘要

本公开提供了一种宏转录组测序数据处理方法及系统，获取宏转录组测序数据；对获取的测序数据进行过滤；对过滤后的测序数据进行组装；根据组装后的测序数据和预设基于原核生物的基因结构特征的概率预测模型，得到测序数据中原核生物的基因预测结果；对得到的测序数据基因结构预测结果添加功能注释；根据基因预测结果的基因定量结果，对添加功能注释后的注释结果进行定量处理。本公开通过对测序数据的过滤、组装、基因预测和注释，实现了宏转录组测序数据更准确和快速的处理，实现了宏转录组基因的更精准和高效注释，提高了宏转录组测序数据的后续处理效率。

说明书

技术领域

本公开涉及生物基因测序数据自动处理技术领域，特别涉及一种宏转录组测序数据处理方法及系统。

背景技术

本部分的陈述仅仅是提供了与本公开相关的背景技术，并不必然构成现有技术。

随着宏基因组学的发展，使用宏基因组学与宏转录组学等多组学研究人体微生物群落与人体健康的相关性，为医疗团队的疾病治疗带来了新的启发与方向。

发明人发现：

(1)当下的宏转录组分析流程，通常仅支持第二代测序技术，然而，随着测序技术的发展，三代测序技术的应用越来越广泛，没有针对三代测序技术的完整的宏转录组分析流程；

(2)当下的宏转录组分析流程，在基因注释阶段使用与现有基因组数据库比对的方法获得其中的基因信息，然而对于宏转录组来讲，其中大多数微生物是未知的，因此使用与现有数据库比对的方法进行基因注释会遗漏掉其中的新的未曾发现的基因序列；

(3)当下宏转录组分析流程在基因表达分析阶段提供的功能注释较少，不能满足当下多组学分析对众多不能功能分析的需求，现有的宏转录组分析流程中，工具太过单一不能灵活选择分析工具。

发明内容

为了解决现有技术的不足，本公开提供了一种宏转录组测序数据处理方法及系统，通过对测序数据的过滤、组装、基因预测和注释，实现了宏转录组测序数据更准确和快速的处理，实现了宏转录组基因的更精准和高效注释，提高了宏转录组测序数据的后续处理效率。

为了实现上述目的，本公开采用如下技术方案：

本公开第一方面提供了一种宏转录组测序数据处理方法。

一种宏转录组测序数据处理方法，包括以下过程：

获取宏转录组测序数据；

对获取的测序数据进行过滤；

对过滤后的测序数据进行组装；

根据组装后的测序数据和预设的基于原核生物的基因结构特征的概率预测模型，得到测序数据中原核生物的基因预测结果；

对得到的测序数据基因结构预测结果添加功能注释；

根据基因预测结果的基因定量结果，对添加功能注释后的注释结果进行定量处理。

进一步的，对获取的测序数据进行过滤，包括：过滤掉宏转录组测序数据中的接头引物、低质量碱基或序列以及rRNA。

进一步的，对二代测序数据，在过滤之前进行质量评估，评估结果以html格式展示，评估结果至少包括：碱基质量箱线图、四种碱基的分布情况、序列GC含量的分布情况、模糊碱基N的分布以及序列长度分布；至少通过设定引物序列、质量阈值和长度阈值进行数据过滤。

进一步的，对三代测序数据，在过滤之前进行质量评估，评估结果至少包括：序列的长度分布、质量分布以及大于预设质量等级的序列数量；至少通过设定长度阈值、序列平均质量阈值、GC含量阈值、对序列头部与尾部要减掉的长度进行数据过滤。

进一步的，对原核生物的rRNA序列建立索引，根据索引对测序序列的rRNA过滤，得到过滤之后的mRNA转录本序列。

进一步的，若为多样本差异表达分析，将所有样本文件合并为一个文件，进行过滤与组装，对组装得到的序列进行基因预测；

对基因预测后的结果进行编码基因定量，首先对预测基因建立索引，然后输入索引文件与测序文件实现基因定量；

在结果文件中，给出每个基因的TPM值，并对差异表达基因进行功能富集。

进一步的，若为单样本分析，至少进行GO、KEGG_Pathway、CAZy、PHI、VFDB、CARD、TCDB和COG功能注释，其中GO、KEGG_Pathway、CAZy与COG功能注释使用emapper.py软件实现，VFDB、CARD与TCDB注释通过把数据库下载到本地，建立索引，进行比对后实现。

本公开第二方面提供了一种宏转录组测序数据处理系统。

一种宏转录组测序数据处理系统，包括：

数据获取模块，被配置为：获取宏转录组测序数据；

数据过滤模块，被配置为：对获取的测序数据进行过滤；

数据组装模块，被配置为：对过滤后的测序数据进行组装；

基因预测模块，被配置为：根据组装后的测序数据和预设的基于原核生物的基因结构特征的概率预测模型，得到测序数据中原核生物的基因预测结果；

注释添加模块，被配置为：对得到的测序数据的基因结构预测结果添加功能注释；

定量分析模块，被配置为：根据基因预测结果的基因定量结果，对添加功能注释后的注释结果进行定量处理。

本公开第三方面提供了一种计算机可读存储介质，其上存储有程序，该程序被处理器执行时实现如本公开第一方面所述的宏转录组测序数据处理方法中的步骤。

本公开第四方面提供了一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的程序，所述处理器执行所述程序时实现如本公开第一方面所述的宏转录组测序数据处理方法中的步骤。

与现有技术相比，本公开的有益效果是：

1、本公开所述的方法、系统、介质或电子设备，通过对测序数据的过滤、组装、基因预测和注释，实现了宏转录组测序数据更准确和快速的处理，实现了宏转录组基因的更精准和高效注释，提高了宏转录组测序数据的后续处理效率。

2、本公开所述的方法、系统、介质或电子设备，同时支持二代、三代测序技术的宏转录组数据的处理和分析，提供了可选择的多种功能注释，并对处理结果实现了可视化。

附图说明

构成本公开的一部分的说明书附图用来提供对本公开的进一步理解，本公开的示意性实施例及其说明用于解释本公开，并不构成对本公开的不当限定。

图1为本公开实施例1提供的宏转录组测序数据处理方法的流程示意图。

图2为本公开实施例5提供的宏转录组测序数据处理方法的流程示意图。

具体实施方式

下面结合附图与实施例对本公开作进一步说明。

应该指出，以下详细说明都是例示性的，旨在对本公开提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本公开所属技术领域的普通技术人员通常理解的相同含义。

需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本公开的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。

在不冲突的情况下，本公开中的实施例及实施例中的特征可以相互组合。

实施例1：

如图1所示，本公开实施例1提供了一种宏转录组测序数据处理方法，包括以下过程：

获取宏转录组测序数据；

对获取的测序数据进行过滤；

对过滤后的测序数据进行组装；

根据组装后的测序数据和预设的基于原核生物的基因结构特征的概率预测模型，得到测序数据中原核生物的基因预测结果；

对得到的测序数据基因结构预测结果添加功能注释；

根据基因预测结果的基因定量结果，对添加功能注释后的注释结果进行定量处理。

对获取的测序数据进行过滤，包括：过滤掉宏转录组测序数据中的接头引物、低质量碱基或序列以及rRNA。

对二代测序数据，在过滤之前进行质量评估，评估结果以html格式展示，评估结果至少包括：碱基质量箱线图、四种碱基的分布情况、序列GC含量的分布情况、模糊碱基N的分布以及序列长度分布；至少通过设定引物序列、质量阈值和长度阈值进行数据过滤。

对三代测序数据，在过滤之前进行质量评估，评估结果至少包括：序列的长度分布、质量分布以及大于预设质量等级的序列数量；至少通过设定长度阈值、序列平均质量阈值、GC含量阈值、对序列头部与尾部要减掉的长度进行数据过滤。

对原核生物的rRNA序列建立索引，根据索引对测序序列的rRNA过滤，得到过滤之后的mRNA转录本序列。

若为多样本差异表达分析，将所有样本文件合并为一个文件，进行过滤与组装，对组装得到的序列进行基因预测；

对基因预测后的结果进行编码基因定量，首先对预测基因建立索引，然后输入索引文件与测序文件实现基因定量；

在结果文件中，给出每个基因的TPM值，并对差异表达基因进行功能富集。

若为单样本分析，至少进行GO、KEGG_Pathway、CAZy、PHI、VFDB、CARD、TCDB和COG功能注释，其中GO、KEGG_Pathway、CAZy与COG功能注释使用emapper.py软件实现，VFDB、CARD与TCDB注释通过把数据库下载到本地，建立索引，进行比对后实现。

实施例2：

本公开实施例2提供了一种宏转录组测序数据处理系统，包括：

数据获取模块，被配置为：获取宏转录组测序数据；

数据过滤模块，被配置为：对获取的测序数据进行过滤；

数据组装模块，被配置为：对过滤后的测序数据进行组装；

基因预测模块，被配置为：根据组装后的测序数据和预设的基于原核生物的基因结构特征的概率预测模型，得到测序数据中原核生物的基因预测结果；

注释添加模块，被配置为：对得到的测序数据基因结构预测结果添加功能注释；

定量分析模块，被配置为：根据基因预测结果的基因定量结果，对添加功能注释后的注释结果进行定量处理。

所述系统的工作方法与实施例1提供的宏转录组测序数据处理方法相同，这里不再赘述。

实施例3：

本公开实施例3提供了一种计算机可读存储介质，其上存储有程序，该程序被处理器执行时实现如本公开实施例1所述的宏转录组测序数据处理方法中的步骤，所述步骤为：

获取宏转录组测序数据；

对获取的测序数据进行过滤；

对过滤后的测序数据进行组装；

根据组装后的测序数据和预设的基于原核生物的基因结构特征的概率预测模型，得到测序数据中原核生物的基因预测结果；

对得到的测序数据基因结构预测结果添加功能注释；

根据基因预测结果的基因定量结果，对添加功能注释后的注释结果进行定量处理。

详细步骤与实施例1提供的宏转录组测序数据处理方法相同，这里不再赘述。

实施例4：

本公开实施例4提供了一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的程序，所述处理器执行所述程序时实现如本公开实施例1所述的宏转录组测序数据处理方法中的步骤，所述步骤为：

获取宏转录组测序数据；

对获取的测序数据进行过滤；

对过滤后的测序数据进行组装；

根据组装后的测序数据和预设的基于原核生物的基因结构特征的概率预测模型，得到测序数据中原核生物的基因预测结果；

对得到的测序数据基因结构预测结果添加功能注释；

根据基因预测结果的基因定量结果，对添加功能注释后的注释结果进行定量处理。

详细步骤与实施例1提供的宏转录组测序数据处理方法相同，这里不再赘述。

实施例5：

如图2所示，本公开实施例5提供了一种宏转录组测序数据处理方法，包括以下过程：

对宏转录组测序数据进行质量过滤，过滤掉其中的引物与低质量碱基，除此之外，由于微生物遗传物质中存在大量的rRNA序列，尽管在测序之前会对其中的rRNA进行去除，但是仍然有残留的可能，因此对于宏转录组测序数据需要对其中的rRNA进行过滤。

经过上述处理，接下来本实施例对宏转录组基因表达情况进行分析，首先对序列进行组装得到包含信息更多的较长的序列片段，之后通过对原核生物的基因结构特征进行训练来预测其中的基因的方法进行基因预测。在完成测序序列的基因预测之后，如果为单一样本分析，接下来为基因提供了众多的功能注释。如果为多样本，则可以分析不同分组间的差异表达基因，并对差异表达基因进行功能富集。

具体的，包括如下步骤：质量评估、质量控制、去除rRNA、序列组装、基因预测、基因表达分析。

S1：质量控制

由于测序技术的有限性，测序得到的原始序列通常需要过滤，过滤掉其中的接头引物、低质量碱基或序列以及宿主DNA或rRNA等，在过滤之前，本实施例提供了FastQC对二代测序测序序列的质量进行评估，评估结果以html展示，包括碱基质量箱线图、四种碱基的分布情况、序列GC含量的分布情况、模糊碱基N的分布以及序列长度分布等信息。

接下来，本实施例提供了Trimmomatic、Cutadapt与fastp三种针对二代测序数据的质量控制软件，用户可自行选择使用哪款软件，以及通过设置指定软件的引物序列、质量阈值、长度阈值等参数实现序列的质量控制。

三代测序数据的质量评估使用Nanoplot质量评估软件，评估结果包括序列的长度分布、质量分布、大于Q5、Q7、Q10、Q15等质量的序列数量等。接下来，本实施例使用NanoFilt与Filtlong实现三代测序数据的质量控制，Nanofilt通过设置长度阈值、序列平均质量阈值、GC含量阈值、对序列头部与尾部要减掉的长度等对序列进行过滤。

Filtlong提供的质量控制参数更多，包括长度阈值、要保留的序列的百分比、序列平均质量阈值、窗口长度与窗口平均质量阈值、对序列的头部与尾部要减掉的长度等，此外，Filtlong还可以使用二代测序数据作为参考序列进行质量控制。

在完成上述原始序列基本的质量控制之后，本实施例使用SortMeRNA软件的indexdb_rna工具对rfam-5.8-database-id98.fasta，rfam-5s-database-id98.fasta，silva-arc-16s-id95.fasta，silva-arc-23s-id98.fasta，silva-bac-16s-id90.fasta，silva-bac-23s-id98.fasta，silva-euk-18s-id95.fasta，silva-euk-28s-id98.fasta等原核生物的rRNA序列建立索引，通过使用sortmerna工具实现了对测序序列的rRNA过滤，并且通过--other参数的设置得到过滤之后的mRNA转录本序列。

S2：序列组装

本实施例提供了Trinity、IDBA-UD、MEGAHIT以及SPAdes四种组装软件用于二代测序数据的组装，Trinity是专门用于转录本组装的组装工具，而IDBA-UD、MEGAHIT以及SPAdes适用于组装测序非常不均匀的测序数据，因此适合于组装宏转录组数据。

本实施例提供了6种用于三代测序数据的组装工具，分别为Canu、Flye、Miniasm、SMARTdenovo、Wtdbg以及NextDenovo。由于三代测序数据错误率较高的特点，因此还需要对三代测序数据进行纠错，本实施例提供了racon与Pilon两种三代数据的纠错工具。用户可灵活选择上述组装工具，并对相应工具的参数进行设置完成组装，之后，本实施例使用Quast软件对组装结果进行评估，包括contig的数量(contig是经过组装的较长的序列)、N50、L50等参数。

S3：基因表达分析

基因表达分析分为两种，一种为单样本分析，一种为多样本差异表达分析，若为多样本差异表达分析，要求用户在最开始时将所有样本文件合并为一个文件，进行上述质量控制与组装等步骤。

本实施例使用Prodigal与MetaGeneMark软件对组装得到的序列进行基因预测。接下来，流程使用salmon软件实现编码基因定量，首先使用salmon index命令对预测基因建立索引，然后使用salmon quant命令输入索引文件与测序文件实现基因定量。在结果文件中，salmon定量给出了每个基因的TPM值。

若为多样本差异分析，本实施例将不同样本的定量结果合并，输入样本的分组信息，以及要对照的分组，实现edgeR差异分析，并绘制热图、火山图，并且使用GSEA软件实现差异基因的GO富集与KEGG Pathway富集。

若为单样本分析，本实施例提供了GO、KEGG_Pathway、CAZy、PHI、VFDB、CARD、TCDB、COG等功能注释。流程使用emapper.py软件实现GO、KEGG_Pathway、COG与CAZy注释，VFDB、CARD与TCDB数据库的注释是通过把上述数据库下载到本地，建立索引，然后使用比对软件比对实现，流程提供了Blastp与Diamond两种比对工具。

本实施例结合上述基因定量的结果对功能注释的结果实现了定量，并使用柱状图进行可视化。此外，本实施例也可以将上述单样本的功能注释结果进行合并使用LEfSe软件发现不同分组间的差异功能，并且也可以对合并结果实现热图可视化。

本实施例提出的针对宏转录组数据分析方法，提供了针对二代、三代测序技术的不同测序特点的工具，使分析结果更精准；提出的宏转录组数据分析方法，提供了可灵活选择的分析工具，增加了分析流程的灵活性；提出的宏转录组基因注释方法，通过使用预测方法与传统参考数据库比对方法相比，更容易挖掘新型基因，避免了基因遗漏；提出的宏转录组数据分析方法，提供了更全面的功能注释信息，帮助用户在多方面分析基因功能。

本领域内的技术人员应明白，本公开的实施例可提供为方法、系统、或计算机程序产品。因此，本公开可采用硬件实施例、软件实施例、或结合软件和硬件方面的实施例的形式。而且，本公开可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器和光学存储器等)上实施的计算机程序产品的形式。

本公开是参照根据本公开实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于一计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体(Read-Only Memory，ROM)或随机存储记忆体(RandomAccessMemory，RAM)等。

以上所述仅为本公开的优选实施例而已，并不用于限制本公开，对于本领域的技术人员来说，本公开可以有各种更改和变化。凡在本公开的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本公开的保护范围之内。