噬菌体作为细菌的病毒，具有感染并杀伤宿主菌的能力。目前，抗生素滥用情况日益严重，导致耐药菌不断增多。为了遏制耐药菌的传播，噬菌体疗法正逐步受到中外科学家的关注和应用^[1]。随着噬菌体疗法的发展，针对不同宿主菌的新噬菌体分离和其基因组测序工作势在必行。在基因组水平上，噬菌体分为4类：双链DNA (dsDNA)噬菌体、单链DNA (ssDNA)噬菌体、双链RNA (dsRNA)噬菌体和单链RNA (ssRNA)噬菌体。在分离出的噬菌体中，以dsDNA噬菌体最为常见，也应用最广。dsDNA噬菌体的基因组分为环状和线状2种，其中线状基因组的末端类型多种多样^[2]。

dsDNA噬菌体基因组的末端类型和其包装机制相关。噬菌体的包装机制相对保守。在dsDNA噬菌体感染周期接近尾声时，其基因组通常会形成串联体，在包装过程中被终止酶切割，形成成熟的染色体。噬菌体使用4种主要机制来识别它们自己的DNA并启动和终止其包装^[3]：(1) 终止酶识别交错切割的特定位点(cos位点)，从而产生具有黏性末端的固定DNA末端，该末端具有5′或3′突出部分(例如Lambda、HK97)。(2) 从噬菌体DNA上识别出一个固定位置，通过延伸合成将在此产生直接末端重复序列(Direct Terminal Repeats，DTR)。这些DTR的大小可以从一百多个碱基(例如T3、T7)到超过一万个碱基(例如T5、Spo1)不等。(3) 终止酶可以在特定的包装位点(pac位点)启动对噬菌体串联体的包装。当噬菌体头部被填满后，随后的切割将在不同位置进行。这导致衣壳含有冗余末端且循环排列的基因组，其注射到宿主细胞后可通过重组使噬菌体基因组成环。(4) 类T4噬菌体使用这种头部包装机制的变体。在这种机制中不识别pac位点，包装是随机启动的。这些噬菌体通常会降解宿主DNA，确保只包装病毒DNA。

除此之外还至少有3种不常见的噬菌体基因组末端包装策略^[3]：(1) 噬菌体P2携带cos位点，但包装底物是环形dsDNA。(2) 噬菌体Mu通过转座在宿主基因组中复制，并携带宿主DNA片段作为其末端。(3) 噬菌体phi29在其DNA末端携带共价结合蛋白。当考虑到只有少量的噬菌体精确研究了其基因组末端包装机制，很可能在自然界中还存在其他包装机制。因此判断出噬菌体的末端类型有助于辅助新分离噬菌体的分类。

随着近几年高通量测序技术的不断发展以及价格的不断降低，大部分噬菌体基因组测序已经应用高通量测序技术来完成。在高通量测序技术大规模应用之前，dsDNA噬菌体的末端序列主要利用大引物(Mega-Primer) PCR的方法进行检测^[4]。现在针对dsDNA噬菌体基因组的高通量测序结果，中外科学家已经建立了相对标准的流程来验证dsDNA噬菌体的基因组末端^{[3, 5-6]}。其中Zhang等提出了利用噬菌体基因组高通量测序的高频序列判定噬菌体基因组末端的方法，这一方法指出，如噬菌体的基因组有末端，那么匹配到此末端序列的Reads就会表现出高频的特点^[7]。在这一思想的指导下，为了判定噬菌体的末端序列，研究者首先应统计高通量测序结果中的高频Reads序列。

Zhang等给出了统计高频Reads序列的方法以及在Linux系统下的命令^[7]。但鉴于大部分的病毒学或微生物学工作者对Linux操作系统不熟，使得这一方法和命令的应用受到了限制。在这一背景下，我们对Zhang等提出的统计高频Reads序列的方法进行了优化，设计出了在Windows操作系统下统计高通量测序数据中高频Reads的软件PhageGT。

1 材料与方法 1.1 所需材料与开发环境

1.1.1 材料

选取我们课题组以前发表的2个病毒(SRT6^[8]和SRT7^[9])和在NCBI数据库随机下载的8个病毒(SRR13926753、SRR13783533、SRR10882832、SRR10882831、SRR10882571、SRR10760098、SRR14121700、SRR10210268)，利用其原始的测序文件(20个，分别为2个病毒组双末端测序所得的正向和反向2个文件)进行软件测试及验证。其中文件中数据的细节如图 1所示。

1.1.2 开发环境

Windows 10、Visual Studio 2019。

1.2 软件界面设计

为了让更多的科研工作者使用起来简单、便捷，此软件的界面设计较为简单，软件界面如图 2和图 3所示。软件全自动化执行，用户只需输入2个特定的初始值并点击“open a file”按钮选择需要操作的序列文件，然后等待软件自动运行即可。软件运行完毕可自动弹出并保存数据文件及展示软件运行情况。

1.3 主要功能

1.3.1 筛选序列

目前高通量测序已经成为研究病毒基因组的首选实验方案。随着测序价格的降低，现在高通量测序的数据量通常较为庞大，结果文件大多为Gb数量级。我们开发的PhageGT软件可以在庞大的序列文件中进行序列筛选，然后提取出文件中的所有序列，而且对重复序列进行计数。

1.3.2 计算频率

此软件可以计算每个序列出现的频率。

1.3.3 根据频率将序列从大到小进行排列

每个序列出现的频率大小不一，其中末端序列出现的次数最为频繁，其频率也比其他序列高。利用此种方法，可以筛选出高频序列，进而对噬菌体基因组的末端序列进行分析。

1.3.4 R值计算

R值代表病毒组测序原始序列中Reads的最高频率和Reads平均频率的比值。可根据R值简单地推出此病毒组是否存在末端。

1.4 主要算法

1.4.1 界面设计

在Visual Studio 2019中按照文件丨新建丨项目丨MFC应用的顺序新建基于对话框的桌面应用程序。然后在工具箱中添加编辑框及“open a file” “drop out”等按钮，并删除“确定” “取消”等按钮。新建模态类型的关于对话框并在“about”对话框中添加软件信息，当用户点击“about”按钮时将跳转到关于对话框。

在主对话框中提示用户应导入的文件类型，在打开文件对话框中将用户导入的文件格式强制为.fq丨.txt和所有文件2种方式，目的是防止用户导入错误的文件格式，但如果用户有其他类型的测序文件时还可以在所有文件选项中导入测序文件。

基于测序进度添加必要的文字提示，例如“Please import the virus group high-throughput sequencing file (**.fq)” “File is being read, please wait patiently...” “Sorting data...” “Calculating R value...” “File is exported to the root directory of Disk D”。

1.4.2 自定义所取碱基个数

对多组序列进行测试后发现对于不同病毒基因组的不同数据，所需碱基数不同结果会有一定的差异，所以在PhageGT软件中设定了自定义读取Reads中碱基个数的功能。建议碱基个数的设定数值为20，这也是Zhang等所建议的数值^[7]。发现选取的碱基个数越多其准确率越高，但随着选取碱基个数的增加，软件运行时间将会越长。当选取碱基个数接近单行序列总长度时，准确率会降低，因为测序时存在的微小误差会造成同一序列不同Reads的个别碱基不同。如果选取长度过长会将此类Reads归为不同类。所以选取合适碱基数量尤为重要。

获取用户输入的字符串，然后将其类型转换为int类型，目的是后面字符读取数量时避免因数字类型不同造成不必要的错误。

1.4.3 序列筛选

因每个序列前都有一个符号“/”，所以可将其作为标识符，遍历序列文件中的每个字符，当遍历到“/”符号时再向后遍历的下一行就是基因序列。新建一个结构体，结构体中私有成员为基因序列和其出现的次数。当第一次遍历到基因序列时将其存入已经新建的数组中，当第二次遍历到基因序列时将其与存入数组中的序列进行比较，若数组中已经存在此序列则继续进行遍历，否则将其存入数组。由此可以得到筛选出的不重复序列。

1.4.4 序列排序

定义comp()函数，此函数的功能为比较2个序列的出现次数大小。然后利用sort()函数进行排序，排序方法为comp()。

1.4.5 序列重复次数及频率计算

在将得到的序列与数组中的序列进行对比时若此序列在数组中已经存在则将其成员变量num加一，最后程序运行完毕时每个变量其成员num为其出现的次数。

设置一个变量记录遍历过的基因序列的总数，然后用每个序列其出现的次数除以总数即是对应基因序列的频率。

1.4.6 R值计算^[7]

R=文件中序列的最高频率/序列的平均频率。最高频率=最高频Reads的频率，平均频率=序列文件中Reads总数/(2×基因组长度)。R < 30，说明没有末端或者末端很少；30 < R < 60，说明具有部分末端；R > 60，说明末端存在。

1.4.7 运行时间

利用colcok()函数记录程序运行的初始时间和结束时间，然后两者相减即为程序运行时间。因其所得结果以秒为单位，因此要对所得结果进行判断若结果小于60则不需进行处理；若结果大于60小于3 600则结果应除以60转换为分钟计时方式；若结果大于3 600则结果应除以3 600转换为小时计时方式。

1.5 使用方法

1.5.1 打开软件

双击打开PhageGT.exe可执行文件。

1.5.2 自定义要读取的碱基个数

在第一个编辑框中输入要选择的碱基个数，如图 4所示。

1.5.3 输入基因组的长度

在第二个编辑框中输入基因组的长度，如图 5所示，其中基因组的长度由CLC Genomics Workbench 20软件将一组正反序列拼接后获得。

1.5.4 浏览文件

点击打开文件进行文件浏览，如图 6所示。

1.5.5 选择文件类型

在文件浏览对话框中选择要导入的文件类型，如图 7所示。

1.5.6 选择文件

选中要导入的文件，然后点击打开导入文件，如图 8所示。

1.5.7 等待软件运行

根据软件进度提示：数据读取，数据排序，数据导出，数据导出完毕，如图 9、图 10、图 11所示，等待测序完成，需等待时间决定机器配置。若读取序列文件较大可能会出现软件未响应，这属于正常现象，软件实际仍在运行，耐心等待即可。

2 结果与讨论 2.1 基因组序列处理完毕后的输出文件

同一病毒正向和反向高通量测序的结果，经此软件处理后输出的前几个碱基序列应是相同的，并且其出现次数与频率也应大致相同，如表 1、表 2、表 3、表 4所示。

2.2 两对基因组的R值大小与处理时间

由表 5可得对于同一组序列其R值大小接近，都可以判断是否存在末端，但是随着序列文件大小的增大处理时间也随之升高。对于不同电脑配置其软件运行时间也是不同的，电脑配置越高其处理时间越短。

2.3 噬菌体基因组末端分析

由表 5可得2个噬菌体基因组均存在末端，由表 1、表 2、表 3、表 4可知噬菌体SRT6、噬菌体SRT7分别存在2个主要的高频序列5′-AACCACGAACCCACTAGGCA-3′和5′-AGCGA CAGCTTACACCCGTT-3′、5′-AGGGACTGAGGC ATAACCAA-3′和5′-TCTCAAAGATACAACTCC CA-3′。如表 5所示，末端序列与一般序列出现频率的比值分别为178.13 (噬菌体STR6的正向)、167.20 (噬菌体STR6的反向)、177.24 (噬菌体STR7的正向)、154.17 (噬菌体STR7的反向)。根据SRT6和SRT7的全基因组序列以及已报道的计算方法^[7]可计算出噬菌体基因组分别存在有1 165 bp (F19968，R21132)和175 bp (F1，R175)的直接重复末端。

2.4 软件通用性和准确性分析

将STR6、STR7和NCBI数据库中随机选取的8个病毒(20个序列文件)同时在Windows系统下使用此软件(碱基个数设定数值为20)和在Linux系统下使用命令“awk 'NR%4==2' XXX.fastq|sort|uniq- c -w20|sort -g -r -o XXX.Freq”进行分析。

根据表 6可得此软件的计算结果与在Linux系统下使用命令行界面计算的结果完全相同，所以可以判断出此软件具有较高的准确性。因为选取的序列是在NCBI数据库中随机选取的，所以可以判断出此软件具有较高的通用性。

在表 7中发现Windows系统下此软件的计算时间比在Linux系统下使用命令操作高了数倍，所以此软件在运行时间上不如Linux命令操作有优势。但此软件可以直接得出R值并判断是否具有末端，并且采用Windows系统下的图形化界面，使操作简洁、直观，方便了不熟悉Linux系统的工作人员。

3 结论

高通量测序技术的发展与使用，使科研人员需要处理的基因序列信息越来越大。针对噬菌体，已有科学家提出可利用基因组高通量测序的数据分析得到其基因组末端信息。但具体的分析过程需要相应的序列处理软件。与人工手动处理这些序列相比借助于计算机处理更节省时间，科研人员可以将更多的时间投入到更有价值的方向。将此前发布的Linux系统下噬菌体基因组高通量测序结果分析的方法^[7]，在Windows系统下设计相应软件，可以让不熟悉Linux系统的科研工作者省去了学习Linux系统的时间，使用熟悉的Windows系统的应用软件即可简单地进行数据处理。

当前挖掘各种生物数据之间的关系变得尤为重要，计算机技术的进步更是推动了这一生物信息学领域的发展，为相应科研人员提供了一个更好的平台，让他们可以借助于计算机快捷、简单地进行数据处理。根据我们测试的2个噬菌体基因组序列末端的验证结果可以看出此软件确实可行，可以较为简单地得到各种我们期望的结果。但是由于目前不同公司给出高通量测序的序列文件的细节可能有所不同，针对不同数据，软件也应相应的进行修改，需要开发与之相匹配的处理软件。在后续版本中我们将会对软件进行升级，以提供可以处理尽可能全部公司原始测序数据的软件。

所有数据集和软件均可从https://zenodo.org/record/4674231#.YHADb-gzZxc免费获得。