CN111354418B - 基于参考基因组注释文件的高通量测序技术动物tRFs数据分析方法 - Google Patents

Abstract

本发明提出了一种基于参考基因组注释文件的高通量测序技术动物tRFs数据分析方法，包括文件准备步骤、tRNA数据库整理步骤、GO、KEGG背景文件整理步骤、下机数据质控步骤、参考基因组比对步骤、tRFs比对注释步骤、差异tRFs分析步骤、靶基因预测及富集分析步骤、网页版报告整理步骤。本发明结果全面，包含涉及到的tRFs分析内容以及其靶基因预测，GO、KEGG富集分析以及对应的可视化展示；自动整理所有分析结果，每一步分析完成之后自动对结果进行汇总统计，可视化，以及逻辑化归类整理，结果文件可直接用于生成网页版报告；所有操作步骤可以溯源，方便错误查询，如果分析报错，会有对应的报错日志信息。

Description

基于参考基因组注释文件的高通量测序技术动物tRFs数据分 析方法

技术领域

本发明属于高通量转录组测序领域，具体涉及一种基于参考基因组注释文件的高通量测序技术动物tRFs数据分析方法。

背景技术

tRFs(长度为15-30nt)是在某些特定情况下，由成熟tRNA或前体tRNA被特异性剪切产生的；tRFs在动物中行使与miRNA类似的功能，参与转录后基因表达调控，研究tRFs的表达差异，tRFs的靶基因、以及其作用机制，是转录组数据分析中的重要一部分，对于动物基因的表达调控机制有重要意义。

tRNA在绝大多数物种的参考基因组注释文件中都有注释，但多数物种缺少 tRFs的注释，相应的trfdb数据库也没有收录相关物种信息，但是基于物种的参考基因组注释文件，也可以进行tRFs的分析。目前还没有针对无参考tRFs 数据的tRFs的数据分析工具，特别是没有自动化的分析实现tRFs测序结果的流程化分析工具，包括tRFs的鉴定，表达量分析和差异分析，靶基因位点分析， GO、KEGG功能富集分析等各个步骤的自动化整合。

现有的动物tRFs高通量分析方法存在如下缺陷：

(1)适用性不强：只能分析现有tRFs数据库的物种；

(2)结果展示不完整：分析结果过于简单，数据挖掘的不深入，缺少数据对应的可视化展示内容。

发明内容

为了克服现有技术所存在的上述缺陷，本发明的目的在于提供一种基于参考基因组注释文件的高通量测序技术动物tRFs数据分析方法。

为了实现本发明的目的，所采用的技术方案是：

本发明提供了一种基于参考基因组注释文件的高通量测序技术动物tRFs数据分析方法，包括如下步骤：

(1)文件准备步骤：

准备config文件，读取后用于进行数据自动化质控以及后续数据分析；

(2)tRNA数据库整理步骤：

使用bedtools软件，利用NCBI或者Ensembl数据库中物种现有的参考基因组序列文件、参考基因组gff注释文件获得物种的tRNA序列信息包括转运的氨基酸，长度，反向密码子以及碱基信息；

(3)GO、KEGG背景文件整理步骤：

使用所述参考基因组序列文件、所述参考基因组gff注释文件提取物种的 mRNA序列；根据所述mRNA序列文件、KEGG数据库、SwissProt数据库，整理GO、 KEGG背景文件，包括mRNA注释文件，GO富集背景文件，KEGG富集背景文件， KEGG注释文件，KEGG pathway文件；

(4)下机数据处理步骤：

将下机得到的原始数据，通过Cutadapt、FastQC、Fastx-Toolkit、 NGS_QC_Toolkit软件去除接头，保留15-41nt长度的序列，然后过滤低质量序列，即：以5个碱基长度为窗口对原始序列进行搜索，当窗口中碱基的平均测序质量低于20时，将从窗口最前端开始的部分截断并舍弃。将过滤后的数据进行去重，获得无重复的序列，并标记所有序列的数量。同时对原始数据和过滤数据量进行统计，并以柱状图展示各个样本不同长度序列的数量分布特征。过滤序列用于后续分析；

(5)参考基因组比对步骤：

使用bowtie软件将上述去重后的序列与参考基因组序列做比对，筛选出碱基错配数小于2的结果，获得参考基因组的比对率；

(6)tRFs比对注释步骤：

使用bowtie软件将所述过滤后的数据合并去重之后与所述tRNA序列做比对，通过比对位置筛选出碱基错配数为0的结果，获得各样本的tRFs的鉴定注释结果，同时计算tRFs的表达量，并进行表达模式聚类分析；

(7)差异tRFs分析步骤：

根据所述注释到的tRFs信息以及表达量结果，使用DESeq或DESeq2进行差异表达分析，筛选同时满足差异倍数(差异倍数>2)和显著性(P值<0.05)的差异表达tRFs，统计并展示可视化结果；

(8)靶基因预测、富集分析步骤：

将所述差异tRFs与所述物种mRNA序列使用miRanda软件进行靶标预测，统计靶标结合位点信息，绘制结合位点示意图；

对上一步预测到的差异tRFs靶基因，利用上述整理好的GO、KEGG背景文件使用超几何分布检验计算方法进行GO功能和KEGG通路的富集分析，计算GO、 KEGG条目在差异tRFs靶基因中是否显著富集的P值，再对P值经Benjamini& Hochberg多重检验纠正后得到FDR；针对富集结果做柱状图和气泡图统计，获得差异tRFs可能参与的功能和代谢通路；

(9)网页版报告整理步骤：

根据所述结果一键化生成tRFs分析的网页版报告，网页版报告对整个分析结果做了汇总，并对每个分析步骤做了描述和对应的图表展示以及弹窗式帮助文档，网页报告设置了内部快捷链接和分析方法介绍/外部网站的链接，方便网页版内部快速跳转以及快速查阅网上资料。

本发明中，第(4)部分是数据自动化质控，(5)-(9)是后续数据分析。

在本发明的一个优选实施方式中，所述文件准备步骤当中config文件中包括：下机数据位置以及对应的样本分析名和分组名、用于差异分析的分组信息、差异倍数参数，生物学重复参数，参考基因组信息等。

在本发明的另一个优选实施方式中，所述tRFs比对注释步骤当中，所述tRFs 比对注释步骤当中tRFs的命名方式为“tRNA-转运氨基酸简称-反向密码子-tRNA 长度-5’或3’端”名称的命名方式，这种命名方式可以直观了解到tRFs来源的tRNA转运的氨基酸信息、长度信息，以及tRFs的来源。

在本发明的又一个优选实施方式中，所述tRFs差异分析步骤中，所述绘制图像包括采用R语言的ggplot2软件包绘制差异表达tRFs上下调统计柱状图、火山图；采用Pheatmap包对差异表达tRFs的表达量绘制热图。

在本发明的又一个优选实施方式中，所述靶基因预测、富集分析步骤中，使用python对靶标结合分值前10的关系对绘制结合位点示例图。

本发明的主要创新点在于：

针对无参考tRFs数据的tRFs测序数据的分析方法。

结果全面，包含涉及到的tRFs分析内容以及其靶基因预测，GO、KEGG富集分析以及对应的可视化展示。

自动整理所有分析结果，每一步分析完成之后自动对结果进行汇总统计，可视化，以及逻辑化归类整理，结果文件可直接用于生成网页版报告。

所有操作步骤可以溯源，方便错误查询，如果分析报错，会有对应的报错日志信息。

为了达到上述技术效果，本发明克服了物种没有对应tRFs数据库的情况下，使用公用数据库中该物种现有的参考基因组、基因组gff文件，结合tRFs的形成机制、序列长度等特性针对性的分析高通量测序中该物种tRFs的表达及差异表达情况,预测tRFs的靶基因，对预测到的靶基因做GO、KEGG富集分析，进而了解tRFs可能影响的生物学功能。

本发明在使用序列比对tRNA数据库时，根据tRFs的形成机制，即成熟的tRNA 序列在一定条件下5’端和3’端会被切割成特定长度(15-30nt)的tRFs，最大限度的挖掘高通量测序结果中存在tRFs信息；这个技术特征对于实现上文所述的技术效果是至关重要的。

附图说明

图1为本发明的流程示意图。

图2为准备文件示意图。

图3位本发明的PCA主成分分析图。

图4为本发明的差异tRFs汇总示意图。

图5为本发明的火山图示意图。

图6为本发明的热图示意图。

图7为本发明的靶标结合位点示例图。

图8为本发明的GO富集分析柱状图。

图9为本发明的KEGG富集分析气泡图。

具体实施方式

结合以下具体实施例和附图，对发明作进一步的详细说明。实施本发明的过程、条件、实验方法等，除以下专门提及的内容之外，均为本领域的普遍知识和公知常识，本发明没有特别限制内容。

本发明提供了一种基于参考基因组注释文件的高通量测序技术动物tRFs数据分析方法，包括如下步骤：

(1)文件准备步骤：

准备config文件，读取后用于进行数据自动化质控以及后续数据分析；

(2)tRNA数据库整理步骤：

使用bedtools软件，利用NCBI或者Ensembl数据库中物种现有的参考基因组序列文件、参考基因组gff注释文件获得物种的tRNA序列信息包括转运的氨基酸，长度，反向密码子以及碱基信息；

(3)GO、KEGG背景文件整理步骤：

使用参考基因组序列文件、参考基因组gff注释文件提取物种的mRNA序列；根据mRNA序列文件、KEGG数据库、SwissProt数据库，整理GO、KEGG背景文件，包括mRNA注释文件，GO富集背景文件，KEGG富集背景文件，KEGG注释文件， KEGG pathway文件；

(4)下机数据处理步骤：

将下机得到的原始数据，通过Cutadapt、FastQC、Fastx-Toolkit、 NGS_QC_Toolkit软件去除接头，保留15-41nt长度的序列，然后过滤低质量序列，即：以5个碱基长度为窗口对原始序列进行搜索，当窗口中碱基的平均测序质量低于20时，将从窗口最前端开始的部分截断并舍弃。将过滤后的数据进行去重，获得无重复的序列，并标记所有序列的数量。同时对原始数据和过滤数据量进行统计，并以柱状图展示各个样本不同长度序列的数量分布特征。过滤序列用于后续分析；

(5)参考基因组比对步骤：

使用bowtie软件将上述去重后的序列与参考基因组序列做比对，筛选出碱基错配数小于2的结果，获得参考基因组的比对率；

(6)tRFs比对注释步骤：

使用bowtie软件将过滤后的数据合并去重之后与tRNA序列做比对，通过比对位置筛选出碱基错配数为0的结果，获得各样本的tRFs的鉴定注释结果，同时计算tRFs的表达量，并进行表达模式聚类分析；

(7)差异tRFs分析步骤：

根据注释到的tRFs信息以及表达量结果，使用DESeq或DESeq2进行差异表达分析，筛选同时满足差异倍数(差异倍数>2)和显著性(P值<0.05)的差异表达 tRFs，统计并展示可视化结果；

(8)靶基因预测、富集分析步骤：

将差异tRFs与物种mRNA序列使用miRanda软件进行靶标预测，统计靶标结合位点信息，绘制结合位点示意图；

对上一步预测到的差异tRFs靶基因，利用上述整理好的GO、KEGG背景文件使用超几何分布检验计算方法进行GO功能和KEGG通路的富集分析，计算GO、 KEGG条目在差异tRFs靶基因中是否显著富集的P值，再对P值经Benjamini& Hochberg多重检验纠正后得到FDR；针对富集结果做柱状图和气泡图统计，获得差异tRFs可能参与的功能和代谢通路；

(9)网页版报告整理步骤：

根据结果一键化生成tRFs分析的网页版报告，网页版报告对整个分析结果做了汇总，并对每个分析步骤做了描述和对应的图表展示以及弹窗式帮助文档，网页报告设置了内部快捷链接和分析方法介绍/外部网站的链接，方便网页版内部快速跳转以及快速查阅网上资料。

实施例1

接受用户的小RNA测序数据，以及相关的数据库信息，同时整理该物种的 tRFs数据库以及GO、KEGG富集的背景文件，然后对所有的数据进行相关的分析，得到每个样本中tRFs的注释信息，并对tRFs进行定量分析，以及样本/组间差异表达分析，靶基因预测和GO、KEGG富集分析，文件准备如图2。

首先是对下机数据进行过滤和数量统计。对下机数据进行去除接头序列，并通过5bp的滑动窗口，对原始序列进行搜索，当窗口中碱基的平均测序质量低于 20时，将从窗口最前端开始的部分截断并舍弃过滤低质量序列。然后过滤掉长度小于15或者大于41nt的序列。然后对高质量数据的重复序列做去重处理，得到无冗余序列。并对原始数据和高质量数据进行数量统计。

同时使用bedtools根据参考基因组gff注释文件和参考基因组序列获得 tRNA序列，接着使用bowtie软件将所有样本高质量序列数据合并去重后与所述整理好的tRNA序列进行0碱基错配比对，reads能从tRNA序列的5’或者3’端开始完全匹配且长度在15-30范围内的认定为潜在的tRFs，对于比对上同一 tRNA的同一位置(5’或3’端)的reads，选取丰度最高的reads作为鉴定到的tRFs，根据获得的tRFs整理tRFs数据库。

接着根据参考mRNA序列文件整理GO、KEGG的背景文件，包括靶基因注释文件，GO背景文件，KEGG背景文件，KEGG注释文件，KEGG pathway文件。

接下来使用bowtie软件将样本无冗余序列与整理好的tRFs数据库比对，筛选出碱基错配数为0的结果汇总得到各样本的tRFs的种类和reads数，同时计算tRFs表达量，进行表达模式分析，使用R语言ggplot2绘制PCA主成分分析图(直观了解样本的聚类情况)，参见图3。

根据所述tRFs的表达量结果，使用DESeq或DESeq2进行差异表达分析，筛选同时满足差异倍数(FoldChange>2)和显著性(Pvalue<0.05)的差异表达tRFs，使用R语言的ggplot2软件包绘制差异表达tRFs的统计柱状图(直观了解各个差异比较组的差异数量)、火山图(直观了解差异tRFs的分布)，采用R语言的 Pheatmap包对差异表达tRFs的表达量绘制热图(直观了解差异tRFs的表达高低)，参见图4、5、6。

根据序列相似性及碱基互补配对，对筛选到的显著差异表达的tRFs进行靶基因预测。以本物种的mRNA序列为目标序列,使用miRanda软件对差异表达的 tRFs进行靶基因预测，使用python对靶标结合分值前10的关系对绘制结合位点示例图，参见图7。然后使用超几何检验计算靶基因富集到哪些GO功能和KEGG 代谢通路上，了解差异tRFs可能行使的功能，参见图8、9。

最终整理所有的分析结果，所有分析内容按类别排放在不用目录下。质控统计结果，序列长度分布图形放在质控目录；将参考基因组的比对注释结果放在参考基因组比对目录；将tRFs比对注释结果放在tRFs比对结果目录；将tRFs表达量以及PCA、样本聚类结果放在tRFs表达目录；将tRFs的差异表达相关的分析放在tRFs差异目录；将差异表达的tRFs对应的靶基因预测结果放在tRFs靶基因预测目录；将靶基因的GO富集分析结果放在GO富集目录；将靶基因的KEGG 通路富集分析结果放在KEGG富集目录。

根据上述目录使用python脚本一键化生成对应的网页化分析报告。

本发明的保护内容不局限于以上实施例。在不背离发明构思的精神和范围下，本领域技术人员能够想到的变化和优点都被包括在本发明中，并且以所附的权利要求书为保护范围。

SEQUENCE LISTING

<110> 基于参考基因组注释文件的高通量测序技术动物tRFs数据分析方法

<120> 上海欧易生物医学科技有限公司

<160> 3

<170> PatentIn version 3.3

<210> 1

<211> 12

<212> DNA

<213> 人工序列

<400> 1

cggggggctc tg 12

<210> 2

<211> 8

<212> DNA

<213> 人工序列

<400> 2

gcaggccc 8

<210> 3

<211> 21

<212> DNA

<213> 人工序列

<400> 3

cgcccccgag tctcgtccgg a 21

Claims (6)

1.一种基于参考基因组注释文件的高通量测序技术动物tRFs数据分析方法，其特征在于，包括如下步骤：

(1)文件准备步骤：

准备config文件，读取后用于进行数据自动化质控以及后续数据分析；

(2)tRNA数据库整理步骤：

利用NCBI或者Ensembl数据库中物种现有的参考基因组序列文件、参考基因组gff注释文件获得物种的tRNA序列信息包括转运的氨基酸，长度，反向密码子以及碱基信息；

(3)GO、KEGG背景文件整理步骤：

使用所述参考基因组序列文件、所述参考基因组gff注释文件提取物种的mRNA序列；根据所述mRNA序列文件、KEGG数据库、SwissProt数据库，整理GO、KEGG背景文件，包括mRNA注释文件，GO富集背景文件，KEGG富集背景文件，KEGG注释文件，KEGG pathway文件；

(4)下机数据质控步骤：

将下机得到的原始数据去除接头，保留15-41nt长度的序列，然后过滤低质量序列，即以5个碱基长度为窗口对原始序列进行搜索，当窗口中碱基的平均测序质量低于20时，将从窗口最前端开始的部分截断并舍弃；将过滤后的数据进行去重，获得无重复的序列，并标记所有序列的数量；同时对原始数据和过滤数据量进行统计，并以柱状图展示各个样本不同长度序列的数量分布特征；过滤序列用于后续分析；

同时根据参考基因组gff注释文件和参考基因组序列获得tRNA序列，接着将所有样本高质量序列数据合并去重后与整理好的tRNA序列进行0碱基错配比对，reads能从tRNA序列的5’或者3’端开始完全匹配且长度在15-30范围内的认定为潜在的tRFs，对于比对上同一tRNA的同一位置，即5’或3’端的reads，选取丰度最高的reads作为鉴定到的tRFs，根据获得的tRFs整理tRFs数据库；

(5)参考基因组比对步骤：

将上述去重后的序列与参考基因组序列做比对，筛选出碱基错配数小于2的结果，获得参考基因组的比对率；

(6)tRFs比对注释步骤：

将所述过滤后的数据合并去重之后与所述tRNA序列做比对，通过比对位置筛选出碱基错配数为0的结果，获得各样本的tRFs的鉴定注释结果，同时计算tRFs的表达量，并进行表达模式聚类分析；

(7)差异tRFs分析步骤：

根据所述注释到的tRFs信息以及表达量结果，使用DESeq或DESeq2进行差异表达分析，筛选同时满足差异倍数和显著性的差异表达tRFs，统计并展示可视化结果；

(8)靶基因预测、富集分析步骤：

将所述差异tRFs与所述物种mRNA序列进行靶标预测，统计靶标结合位点信息，绘制结合位点示意图；

对上一步预测到的差异tRFs靶基因，利用上述整理好的GO、KEGG背景文件使用超几何分布检验计算方法进行GO功能和KEGG通路的富集分析，计算GO、KEGG条目在差异tRFs靶基因中是否显著富集的P值，再对P值经Benjamini&Hochberg多重检验纠正后得到FDR；针对富集结果做柱状图和气泡图统计，获得差异tRFs可能参与的功能和代谢通路；

(9)网页版报告整理步骤：

根据结果一键化生成tRFs分析的网页版报告，网页版报告对整个分析结果做汇总，并对每个分析步骤做描述和对应的图表展示以及弹窗式帮助文档，网页报告设置内部快捷链接和分析方法介绍/外部网站的链接，实现网页版内部快速跳转以及快速查阅网上资料。

2.如权利要求1所述的基于参考基因组注释文件的高通量测序技术动物tRFs数据分析方法，其特征在于，所述文件准备步骤当中config文件中包括：下机数据位置以及对应的样本分析名和分组名、用于差异分析的分组信息、差异倍数参数，生物学重复参数，参考基因组信息。

3.如权利要求1所述的基于参考基因组注释文件的高通量测序技术动物tRFs数据分析方法，其特征在于，所述tRFs比对注释步骤当中tRFs的命名方式为“tRNA-转运氨基酸简称-反向密码子-tRNA长度-5’或3’端”名称的命名方式。

4.如权利要求1所述的基于参考基因组注释文件的高通量测序技术动物tRFs数据分析方法，其特征在于，所述差异tRFs分析步骤中，所述绘制图像包括采用R语言的ggplot2软件包绘制差异表达tRFs上下调统计柱状图、火山图；采用Pheatmap包对差异表达tRFs的表达量绘制热图。

5.如权利要求1所述的基于参考基因组注释文件的高通量测序技术动物tRFs数据分析方法，其特征在于，所述靶基因预测、富集分析步骤中，使用python对靶标结合分值前10的关系对绘制结合位点示例图。

6.如权利要求1所述的基于参考基因组注释文件的高通量测序技术动物tRFs数据分析方法，其特征在于，所述步骤(7)中，差异倍数>2，显著性P值<0.05。

Images