一种基于二代测序技术的consensus序列统计分析、可视化方法

技术领域

本发明属于二代测序技术领域，具体涉及一种基于二代测序技术的consensus序列统计分析、可视化方法。

背景技术

长期以来，测序技术一直是分子生物学相关研究中最常用的技术手段之一。基于测序技术的人类基因组计划、转录组分析、微生物基因组重测序，单核苷酸多态性等方面的分析也同时促进了生物学其他领域研究的进步和发展。DNA测序技术的发展为人类探索自身和其他生命的奥秘提供了可能，同时，基因组学时代的来临对DNA测序技术也提出了更高的要求，推动了DNA测序技术的不断进步。目前科学研究已经进入高通量测序时代，从单一、局部的基因或基因片段的研究转变成了对整个基因组的研究，在基因组从头测序和转录组测序中应用较广，继第一代测序技术之后，随着第二代、第三代测序技术的兴起和发展，测序技术逐步朝着高通量、低成本、高安全性和商业化的方向发展。近几年来，尤其是以单分子实时测序为代表的第三代测序技术开始进入人们的视野，该测序技术跨越了第一代、第二代较短读长而直接对DNA单个分子进行测序的新一代测序平台应用日益广泛。

通过结合二代测序技术及下游生信分析来获得微生物尤其是病毒的全长基因组序列已成为新物种发现，进化分析，溯源研究，流行病学调查等领域的常用方法。但是在低含量样品中的微生物尤其是病毒全长基因组获取时，常常伴随有少量或者没有二代测序reads覆盖的区域(gaps)产生，造成无法直接通过二代测序及下游生信分析获得全长基因组，需要结合一代Sanger测序的方法来验证和补充这些区域。当进行批量生物样品的第二代测序分析时，会产生大量的consensus序列，同时伴随着大量的gaps和简并碱基，给后续一代测序的引物设计步骤带来繁重的工作。

发明内容

针对上述现有技术中存在的问题，本发明的目的是提供一种基于二代测序技术的consensus序列统计分析、可视化方法，解决当前进行批量生物样品的第二代测序分析时，会产生大量的consensus序列，同时伴随着大量的gaps和简并碱基，给后续一代测序的引物设计步骤带来繁重的工作等问题。

本发明的一个目的在于提供一种基于二代测序技术的consensus序列统计分析方法。

所述统计分析方法，包括如下步骤：

S1、获取二代测序数据中的consensus序列；

S2、对所述consensus序列进行文件序列数的统计，之后提取每条序列中每个位置的碱基并进行分析判断，找到每条序列中gaps和/或简并碱基；

S3、对每条序列中所述gaps和/或简并碱基的类型、个数和位置进行统计，并获取每条序列的长度；

S4、分别计算所述每条序列的序列覆盖度和所述gaps和/或简并碱基位置的标准差，输出结果，完成数据统计分析。

通过对二代测序分析过程中产生的consensus序列中的gaps及简并碱基的位置、数量、大小和分散程度进行自动化统计分析，确定每条序列中的gaps及简并碱基的具体信息，从而能更快、更高效的剔除这些gaps和简并碱基的片段，完成整个基因组的测序工作。

进一步地，步骤S1中，所述consensus序列的文本格式为fasta格式。

进一步地，步骤S2具体操作方法如下：

S21、对所述consensus序列进行文件序列数的统计，并将统计的数值作为总的序列分析次数循环变量；

S22、设置两个gaps和简并碱基位置连续性判定变量，并对序列所含总碱基数进行统计，作为序列碱基分析次数循环变量；

S23、通过for函数和if函数对每条序列中每个位置的碱基进行提取并进行分析判断；

S24、基于所述序列分析次数循环变量和所述序列碱基分析次数循环变量分别对每条序列中每个位置的碱基进行提取和判定，找到每条序列中gaps和/或简并碱基。

本发明的另外一个目的在于提供一种基于二代测序技术的consensus序列可视化方法。

所述可视化方法，包括如下步骤：

S1、获取二代测序数据中的consensus序列；

S2、对所述consensus序列进行文件序列数的统计，之后提取每条序列中每个位置的碱基并进行分析判断，找到每条序列中gaps和/或简并碱基；

S3、设置两个初始变量，用于对所述gaps和/或简并碱基的位置以及所述每条序列的位置进行标定；

S4、将所述标定的gaps和/或简并碱基的位置信息以及每条序列的位置信息载入ggplot2画图包，进行可视化分析结果的展示。

进一步地，步骤S1中，所述consensus序列的文本格式为fasta格式。

进一步地，步骤S2具体操作方法如下：

S21、对所述consensus序列进行文件序列数的统计，并将统计的数值作为总的序列分析次数循环变量；

S22、设置两个gaps和简并碱基位置连续性判定变量，并对序列所含总碱基数进行统计，作为序列碱基分析次数循环变量；

S23、通过for函数和if函数对每条序列中每个位置的碱基进行提取并进行分析判断；

S24、基于所述序列分析次数循环变量和所述序列碱基分析次数循环变量分别对每条序列中每个位置的碱基进行提取和判定，找到每条序列中gaps和/或简并碱基。

与现有技术相比，本发明具有如下优点：

1)通过对二代测序分析过程中产生的consensus序列中的gaps及简并碱基的位置、数量、大小和分散程度进行自动化统计分析，确定每条序列中的gaps及简并碱基的具体信息，从而能更快、更高效的剔除这些gaps和简并碱基的片段，完成整个基因组的测序工作；

2)将得到的consensus序列中的gaps及简并碱基的位置、数量、大小进行可视化展示，便于后续补gaps时，节省大量人力和时间来进行引物设计工作。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单的介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明consensus序列统计分析、可视化方法原理图；

图2为本发明实施例1中consensus序列统计分析结果图；

图3为本发明实施例1中consensus序列进一步统计分析结果图；

图4为本发明实施例2中consensus序列可视化结果图。

具体实施方式

下面将结合本发明实施例中的附图对本发明实施例中的技术方案进行清楚、完整地描述。显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例，基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明提供了一种基于二代测序技术的consensus序列统计分析、可视化方法，具体原理图见图1。

实施例1一种基于二代测序技术的consensus序列统计分析方法

首先设计一个gaps_degenerate_bases_statistics.sh脚本，之后按照如下步骤进行consensus序列统计分析：

S1、准备由二代测序数据分析产生的fasta格式的consensus序列文件；

S2、gaps_degenerate_bases_statistics.sh脚本对所述consensus序列进行文件序列数的统计，之后提取每条序列中每个位置的碱基并进行分析判断，找到每条序列中gaps和/或简并碱基；

进一步地，步骤S2具体操作方法如下：

S21、对所述consensus序列进行文件序列数的统计，并将统计的数值作为总的序列分析次数循环变量；

S22、设置两个gaps和简并碱基位置连续性判定变量，都赋值为“0”作为初始值或者不连续判定值，随后脚本通过将序列名字和序列内容分别进行提取，并对序列所含总碱基数进行统计，作为序列碱基分析次数循环变量；

S23、通过for函数和if函数对每条序列中每个位置的碱基进行提取并进行分析判断，对于是gaps和简并碱基的位置，进行位置统计和碱基提取，同时通过对预先设置的两个gaps和简并碱基位置连续性判定变量是否为“0”进行判断，来确定gaps和简并碱基的位置是否连续，如果为“0”则，将其修改为“1”，表示从该位置起出现了gaps或者简并碱基，如果不为“0”，则不对变量的值进行修改。对于不是gaps和简并碱基的位置同样会用if函数进行位置连续性判定变量是否为“0”的判断，如果为“0”，则不对变量值进行修改，如果不为零，则将其赋值为“0”，表示位置连续性在此中断；

S24、基于所述序列分析次数循环变量和所述序列碱基分析次数循环变量分别对每条序列中每个位置的碱基进行提取和判定，找到每条序列中gaps和/或简并碱基；

S3、对每条序列中所述gaps和/或简并碱基的类型、个数和位置进行统计，并获取每条序列的长度，结果见图2；

图2为gaps_degenerate_bases_statistics_outcome.txt结果文件，从左到右每列依次代表“所分析的consensus序列名”，“gap或者简并碱基位置”，“gap或者简并碱基类型”。

S4、分别计算所述每条序列的序列覆盖度和所述gaps和/或简并碱基位置的标准差，输出结果，完成数据统计分析，结果见图3。

图3为gaps_degenerate_bases_statistics_outcome.tsv结果文件(可用excel软件打开)，每行从左到右依次表示的是“序列名”，“consensus序列长度”，“未确定的碱基个数”，“consensus序列的覆盖度”，“gaps位置个数”，“简并碱基位置个数”，“未确定的碱基的标准差(用于表示gaps和简并碱基的分散程度，越小越集中，越容易设计引物补gaps)”。

实施例2一种基于二代测序技术的consensus序列可视化方法

设计gaps_degenerate_bases_display.sh，gaps_degenerate_bases_display.R脚本，之后按照如下步骤进行consensus序列可视化展示：

S1、准备由二代测序数据分析产生的fasta格式的consensus序列文件；

S2、gaps_degenerate_bases_display.sh脚本对所述consensus序列进行文件序列数的统计，之后提取每条序列中每个位置的碱基并进行分析判断，找到每条序列中gaps和/或简并碱基；

进一步地，步骤S2具体操作方法如下：

S21、对所述consensus序列进行文件序列数的统计，并将统计的数值作为总的序列分析次数循环变量；

S22、设置两个gaps和简并碱基位置连续性判定变量，并都赋值为“0”作为初始值或者不连续判定值，随后脚本对序列名字和序列内容分别进行提取，并对序列所含总碱基数进行统计，作为序列碱基分析次数循环变量；

S23、通过for函数和if函数分别对序列中每个位置的碱基进行提取并进行是否为gaps或者简并碱基的判断，对于判断为“是”的碱基，进行位置统计和碱基提取，同时通过对gaps和简并碱基位置连续判性定变量是否为“0”进行判断，来确定gaps和简并碱基的位置是否连续，如果为“0”则，将其修改为“1”，表示从该位置起出现了gaps或者简并碱基，如果不为“0”，则不对变量的值进行修改；对于不是gaps和简并碱基的位置，同样会用if函数进行gaps和简并碱基位置连续判性定变量是否为“0”的判断，如果为“0”，则不对变量值进行修改，如果不为零，则将其赋值为“0”，表示位置连续性在此中断；

S24、基于所述序列分析次数循环变量和所述序列碱基分析次数循环变量分别对每条序列中每个位置的碱基进行提取和判定，找到每条序列中gaps和/或简并碱基。

S3、设置两个初始变量，用于对所述gaps和/或简并碱基的位置以及所述每条序列的位置进行标定；

S4、利用gaps_degenerate_bases_display.R脚本，将所述标定的gaps和/或简并碱基的位置信息以及每条序列的位置信息载入ggplot2画图包，进行可视化分析结果的展示，结果见图4。

图4为gaps_degenerate_bases_display.R脚本运行后的可视化结果，每一列表示一条consensus序列，黑色部分表示gaps区域，浅灰色部分表示简并碱基区域，图中最下面的黑线表示的是序列位置坐标；另外，每个gaps和简并碱基的位置和数量也表示在了黑色或浅灰色区域的中心，可用pdf编辑器打开放大查看。

以上实施例，仅为本发明的具体实施方式，用以说明本发明的技术方案，而非对其限制，本发明的保护范围并不局限于此，尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，其依然可以对前述实施例所记载的技术方案进行修改或可轻易想到变化，或者对其中部分技术特征进行等同替换；而这些修改、变化或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围，都应涵盖在本发明的保护范围之内。