CN117238376B

CN117238376B - 基于二代测序技术的病毒载体序列分析系统和方法

Info

Publication number: CN117238376B
Application number: CN202311261744.4A
Authority: CN
Inventors: 陈靓靓; 李源; 王佳伟; 林婕; 潘轶
Original assignee: Shanghai Xuzhenda Biotechnology Co ltd
Current assignee: Shanghai Xuzhenda Biotechnology Co ltd
Priority date: 2023-09-27
Filing date: 2023-09-27
Publication date: 2024-04-30
Anticipated expiration: 2043-09-27
Also published as: CN117238376A

Abstract

本发明提供一种基于二代测序技术的病毒载体序列分析系统和方法。该系统包含载体序列组装模块、载体序列和理论序列差异分析模块。本发明的目的在于为病毒载体提供基因组层面的有效性和安全性评估。该系统和方法完全基于二代测序技术，不需要一代测序和三代测序补充，具有成本低、通量高、周期短和准确性高的优点。在病毒载体序列组装方面，该系统和方法可以组装到完整基因组的精度，优于已有的传统分析方法。

Description

基于二代测序技术的病毒载体序列分析系统和方法

技术领域

本发明涉及病毒载体序列分析系统和方法，尤其涉及基于二代测序技术的病毒载体序列分析系统和方法。

背景技术

基因治疗制品通常由含有工程化基因构建体的载体或递送系统组成，其活性成分可为DNA、RNA、基因改造的病毒、细菌或细胞，通过将外源基因导入靶细胞或组织，替代、补偿、阻断、修正特定基因，以达到治疗疾病的目的。

用于基因治疗制品的常见的载体系统是病毒载体和质粒DNA载体，载体设计与构建方案基于临床有效性和安全性考虑，通常基于基因治疗制品的作用机制，如通过编码功能性蛋白质的转基因表达，或采用RNA干扰、小RNA或基因编辑等方式，采用基因沉默、外显子跳跃、基因调控或基因敲除等方式修复、添加或删除特定的基因序列，进行载体的设计与构建。

为了保证最终治疗产品的有效性及安全性，《人用基因治疗制品总论》等监管文件提出，应对基因治疗制品载体的完整序列进行分析，确认序列与理论预期相符。

病毒载体是基因治疗制品中最常用到的载体系统。病毒载体的序列分析关系到临床申报和批次放行，有一定的实效性要求，需要准确、经济、快速的方案。二代测序相较一代测序和三代测序存在于速度快、成本低、覆盖度高、准确度高的优点，但是二代测序产出的序列为短序列(reads)，需要较好的完整序列分析方法来匹配病毒载体序列分析需求。除了基因治疗制品外，病毒载体还在基因编辑研究、细胞治疗、疫苗制备等多个生物医药领域有广泛的应用，好的病毒载体序列分析方案对各个领域的病毒载体应用都有所帮助。

病毒载体序列分析与传统的病毒序列分析有一定到相似性，但是序列结构与自然界中的病毒有较大差异。病毒载体序列是人工设计的，其中包含的功能元件来自自然界多个不同生物物种，根据载体设计方法、应用场景、生产工艺不同存在多种多样的差异，导致病毒载体之间不像天然病毒同物种之间那样相似，没有通用的参考基因组。

病毒载体分析的目的也与传统病毒分析存在差异。病毒载体分析首先需要验证载体序列全长是否完整且正确，因而需要通过序列拼接能得到完整基因组的精度。而传统病毒分析主要研究病毒的分类和功能特性，通常仅要求序列拼接能得到长于分类单元或功能单元的序列，通常拼接到contig(也称为重叠群)或scaffold的精度。在基因组组装领域，contig、scaffold和完整基因组都是指基因组组装的不同阶段或结果。contig可以反映出基因组中不同区域的顺序和相对位置，但是由于拼接算法以及基因组中的重复序列等因素，存在区间不连续、缺失、重复的情况。scaffold是指利用其他信息(如配对末端序列、物理图谱等)对contig进行排序和连接，形成更长的序列，但也经常存在缺失、分段和结构错误的情况。完整基因组通常是指所有的区间都被覆盖，实现了对基因组的完整重建。

其次，病毒载体样本中包含多个拷贝的载体序列，其中部分拷贝的突变也可能对载体功能造成不良影响，分析需要关注数据中的杂合低频突变，而传统病毒分析中的突变分析通常在拼接后进行突变分析，损失了突变频率和大部分杂合突变信息。最后，载体生产过程中用到了动物、细菌等其他宿主细胞，需要关注载体中是否存在与宿主基因组重组的情况，这部分分析现有的病毒分析产品均没有涉及。

目前尚未有其他成熟的病毒载体完整序列分析方法公开，下面将列举一些已公开的其他病毒分析方法，并解释其在病毒载体序列应用场景的局限性(后文简称为局限性，皆特指本应用场景的局限性)。

专利“一种基于二代和三代测序技术的宏病毒组分析方法”(公开号CN115691679A)公开了一种基于二代和三代测序的病毒序列分析方法，主要用于病毒物种鉴定，能够结合二代和三代数据拼接得到相对完整的病毒基因组信息。局限性在于：需要用到较为昂贵的三代测序技术，不利于周期和成本控制；此外出于鉴定的目的，该发明仅将病毒序列拼接到contig的精度，没有拼接病毒完整基因组；并且不能进行突变和重组的分析。

专利“病毒基因组鉴定和拼接的方法及应用”(公开号CN116072222A)公开了病毒鉴定和基因组组装的方法，主要用于病毒物种鉴定和新病毒发现，能够组装得到较为完整的病毒基因组信息。局限性在于：该发明首先将病毒序列组装为contig，再基于NCBI NR库中基因最相似的参考基因组拼接为病毒基因组全长，不适用于病毒载体这种包含大量人工元件的情况；并且不能进行突变和重组的分析。

专利“一种病毒测序序列的自动化分析方法及系统”(公开号CN112863599A)公开了一种包含病毒组装、鉴定和变异检测等步骤的分析方法和系统，主要用于病毒鉴定、变异检测、系统发育分析和功能注释。局限性在于：该发明仅将序列拼接到contig的精度，没有拼接完整基因组；该发明检测病毒变异的方法为将组装后的contig比对到参考基因组，检测contig和参考基因组之间的差异，损失了突变频率和大部分杂合突变信息；并且不能进行重组分析。

专利“一种适用于矮缩病毒科和双生病毒科病毒基因组拼接的转录组测序方法”(公开号CN107475449A)公开了一种病毒基因组组装方法，首先通过比对找出目标病毒的序列，再通过拼接获得重叠群(contig和scaffold)，中间使用Genebank中的已知参考基因组比对筛选掉不匹配的序列，最后通过基因克隆技术和测序(一代测序)对空白区域进行填充，获得病毒完整基因组。局限性：该发明通过基因克隆和一代测序技术来对空白区域进行填充，基因克隆需进行引物设计、PCR扩增、电泳回收、克隆菌转化和筛选等复杂的实验步骤，实验周期较长成本较高；并且该发明不能进行突变和重组的分析。

发明内容

病毒载体广泛应用于细胞基因治疗和疫苗研发领域，有较高的安全性检测需求。本发明的目的在于为病毒载体提供基因组层面的有效性和安全性评估。因此本发明公开了一种基于二代测序技术的病毒载体序列分析系统和方法。具体技术方案如下：

第一方面，本发明提供一种基于二代测序技术的病毒载体序列分析系统，包含：

载体序列组装模块，用于通过de Bruijn图算法将样本的二代测序数据reads拼接为contigs序列；还用于通过从addgene(https://www.addgene.org/)获取所有目标类别载体的骨架序列，使用分子遗传相关性计算工具mash计算各个目标类别载体的骨架序列和contigs序列之间的遗传相关性，保留和样本序列分子遗传特征最相似的n个载体的骨架序列作为载体参考序列库；以及还用于将载体参考序列库输入shiver程序进行初始化后与contigs序列进行比对融合为单条参考序列同时生成softclip数据，再将样本的二代测序数据reads比对到融合得到的单条参考序列上重构出初步全长序列，用softclip数据对初步全长序列进行矫正得到完整的载体全长序列；n为5～20，优选为8～12，更优选为10；

载体序列和理论序列差异分析模块，用于对载体序列组装模块得到的载体全长序列与病毒载体设计的理论序列的差异位点和差异片断进行分析得到突变信息。

在一些实施方案中，载体序列和理论序列差异分析模块利用基因组比对软件MUMmer4(https://mummer4.github.io/)的核酸序列比对程序nucmer对载体序列组装模块得到的载体全长序列和病毒载体设计的理论序列进行比对；利用MUMmer4的比对筛选程序delta-filter采用LIS(longest increasing subsequence)算法对nucmer比对结果进行过滤，将载体序列组装模块得到的载体全长序列映射到病毒载体设计的理论序列的最佳参考位置；利用MUMmer4的核酸序列差异分析程序dnadiff基于过滤后的比对文件分析差异位点和差异片断。

在一些实施方案中，还包含载体突变分析模块，用于利用二代数据比对软件bwa(https://github.com/lh3/bwa)将样本的二代测序数据reads比对到病毒载体设计的理论序列上，并利用比对数据处理软件samtools(https://github.com/samtools/samtools)将比对数据处理为排序过的bam文件；还用于利用突变分析软件VarDict(https://github.com/AstraZeneca-NGS/VarDict)的VarDict程序分析位点突变，并利用VarDict的teststrandbias.R和var2vcf_valid.pl过滤掉低可信度突变，获取vcf表格格式的突变结果；该突变结果包括突变位点与突变频率。

在一些实施方案中，还包含载体重组分析模块，用于利用微生物分析软件vsearch(https://github.com/torognes/vsearch)的fastq_mergepairs指令基于重叠序列将二代测序的两端reads合并为单端序列；还用于利用基因组重组比对程序ViReMa(https://sourceforge.net/projects/virema/)将样本的二代测序数据reads比对到病毒载体设计的理论序列上，对于reads上不能比对上理论序列的片断，进一步基于宿主细胞参考基因组索引将其与宿主序列进行比对，分析样本的载体序列与宿主基因组序列的重组情况。

在一些实施方案中，载体序列组装模块还用于过滤掉与样本序列mash遗传距离为1(即相关性为0)的不相关载体的骨架序列，进一步基于mash遗传距离对载体的骨架序列进行排序从而获得和样本序列分子遗传特征最相似的n个载体的骨架序列。

在一些实施方案中，载体序列组装模块还用于在进行de Bruijn图算法组装contigs序列时，使用多个不同的k值进行切割，从而得到不同长度的k-mers，从而提高contigs序列组装的准确性和连续性。

第二方面，本发明还提供一种基于二代测序技术的病毒载体序列分析方法，包括以下步骤：

S1、载体序列组装：

S1.1、de Bruijn图算法初步组装：使用多尺度de Bruijn图算法将样本二代测序数据reads初步拼接为contigs序列；

S1.2、构建载体参考序列库：a.从addgene(https://www.addgene.org/)获取所有目标类别载体的骨架序列；b.使用分子遗传相关性计算工具mash计算各个载体骨架序列和初步组装步骤得到的样本contigs序列之间的遗传相关性，保留和样本序列分子遗传特征最相似的n个载体的骨架序列作为载体参考序列库；

S1.3、基因组结构重构：a.将载体参考序列库输入shiver程序进行初始化；b.将步骤S1.1得到的contigs序列与载体参考序列库进行比对融合为单条参考序列同时生成softclip数据；c.将样本的二代测序数据reads比对到融合得到的单条参考序列上，重构出初步全长序列，并利用比对产生的softclip数据对初步全长序列进行矫正，获得完整的载体全长序列；n为5～20，优选为8～12，更优选为10；

S2、载体序列和理论序列差异分析：对步骤S1得到的载体全长序列与病毒载体设计的理论序列的差异位点和差异片断进行分析得到突变信息。

在一些实施方案中，包括步骤S2包括：

S2.1、序列全长比对：利用基因组比对软件MUMmer4(https://mummer4.github.io/)的核酸序列比对程序nucmer对步骤S1得到的载体全长序列和病毒载体设计的理论序列进行比对；

S2.2、序列比对过滤：利用基因组比对软件MUMmer4(https://mummer4.github.io/)的比对筛选程序delta-filter利用LIS(longest increasingsubsequence)算法对nucmer比对结果进行过滤，将步骤S1得到的载体全长序列映射到病毒载体设计的理论序列的最佳参考位置；

S2.3、利用基因组比对软件MUMmer4(https://mummer4.github.io/)的核酸序列差异分析程序dnadiff基于过滤后的比对文件分析差异位点和差异片断。

在一些实施方案中，还包括以下步骤：

S3、载体突变分析：

S3.1、二代测序reads比对到病毒载体设计的理论序列：利用二代数据比对软件bwa(https://github.com/lh3/bwa)将样本的二代测序数据reads比对到病毒载体设计的理论序列上，并利用比对数据处理软件samtools(https://github.com/samtools/samtools)将比对数据处理为排序过的bam文件；

S3.2、突变位点与突变频率分析：利用突变分析软件VarDict(https://github.com/AstraZeneca-NGS/VarDict)的VarDict程序分析位点突变，并利用VarDict的teststrandbias.R和var2vcf_valid.pl过滤掉低可信度突变，获取vcf表格格式的突变结果；该突变结果包括突变位点与突变频率。

在一些实施方案中，还包括以下步骤：

S4、载体重组分析：

S4.1、双端数据合并：利用微生物分析软件vsearch(https://github.com/torognes/vsearch)的fastq_mergepairs指令基于重叠序列将二代测序的两端reads合并为单端序列；

S4.2、序列重组位点分析：利用基因组重组比对程序ViReMa(https://sourceforge.net/projects/virema/)将样本的二代测序数据reads比对到病毒载体设计的理论序列上，对于reads上不能比对上理论序列的片断，进一步基于宿主细胞参考基因组索引将其与宿主序列进行比对，分析样本的载体序列与宿主基因组序列的重组情况。

在一些实施方案中，步骤S1还包括过滤掉与样本序列mash遗传距离为1(即相关性为0)的不相关载体的骨架序列，进一步基于mash遗传距离对载体的骨架序列进行排序从而获得和样本序列分子遗传特征最相似的n个载体的骨架序列。

在一些实施方案中，步骤S1还包括在进行de Bruijn图算法组装contigs序列时，使用多个不同的k值进行切割，从而得到不同长度的k-mers，从而提高contigs序列组装的准确性和连续性。

第三方面，本发明还提供一种电子设备，包括：处理器以及存储器，所述处理器与所述存储器连接；

所述存储器，用于存储所述处理器的计算机程序；

其中，所述处理器被配置为通过执行所述计算机程序来实现第二方面所述的方法。

第四方面，本发明还提供一种计算机可读存储介质，所述计算机可读存储介质上存储有计算机程序，所述计算机程序被处理器执行时实现第二方面所述的方法的步骤。

在一些实施方案中本发明可以一站式地完成载体序列组装、载体序列和理论序列差异分析、载体突变分析、载体重组分析等序列分析工作，从而对可确认病毒载体序列是否与原本设计相符。本发明所涉及的系统和方法完全基于二代测序技术，不需要一代测序和三代测序补充，具有成本低、通量高、周期短和准确性高的优点；具有组装基因组完整度高、可检测低频杂合突变和基因组重组的技术优点。在病毒载体序列组装方面，该系统和方法可以组装到完整基因组的精度，优于已有的传统分析方法。在病毒载体突变检测方面，该系统和方法可以检测病毒载体样本中存在的点突变和插入缺失突变，该系统和方法可以检测纯合突变和最低至5％的低频率杂合突变。在病毒载体重组检测方面，该系统和方法可以检测病毒载体序列与宿主细胞等序列的重组情况，可以检测纯合重组和低至5％的低频杂合重组。

以下将结合附图对本发明的构思、具体结构及产生的技术效果作进一步说明，以充分地了解本发明的目的、特征和效果。

附图说明

图1是本发明一实施例的整体流程图。

图2是病毒载体序列组装模块的流程图。

图3是病毒载体序列和理论序列差异分析模块的流程图。

图4是病毒载体突变分析模块的流程图。

图5是病毒载体重组分析模块的流程图。

具体实施方式

为了使发明实现的技术手段、创造特征、达成目的和功效易于明白了解，现结合具体图示，进一步阐述本发明。但本发明不仅限于以下实施的案例。

须知，本说明书所附图式所绘示的结构、比例、大小等，均仅用以配合说明书所揭示的内容，以供熟悉此技术的人士了解与阅读，并非用以限定本发明可实施的限定条件，故不具技术上的实质意义，任何结构的修饰、比例关系的改变或大小的调整，在不影响本发明所能产生的功效及所能达成的目的下，均应仍落在本发明所揭示的技术内容得能涵盖的范围内。

为了说明本发明的方法和效果，以下将以基于一例已知的慢病毒载体序列构建模拟二代测序数据，以之为例进行分析说明。慢病毒载体是病毒载体中常见的一类，慢病毒载体由人类免疫缺陷病毒改造而来，具有整合效率高、目标基因表达稳定等优点，在细胞治疗和基因治疗领域有广泛的应用。该慢病毒载体序列来源于NCBI，总长度4555bp，数据编号DL196061.1(https://www.ncbi.nlm.nih.gov/nuccore/DL196061.1)，为一例可表达黄病毒科蛋白的重组慢病毒载体。

本例数据的载体种类、应用领域、分析参数仅作为示例。本发明应用范围不限于该领域、该载体种类，本发明的实施方式包括但不限于以下示例。

测试数据构建步骤如下：

(1)载体序列获取：从NCBI目标序列页面(https://www.ncbi.nlm.nih.gov/nuccore/DL196061.1)下载fasta格式的载体全长序列文件，将文件命名为DL196061.1.fasta。

(2)突变和重组载体序列生成：为了对载体突变检测模块和载体重组检测模块的方法和效果进行说明，我们对载体序列DL196061.1进行编辑，生成了变异序列DL196061.1_variant1和DL196061.1_variant2，各序列导入的变异列表见表1，具体导入点突变、插入缺失突变信息见表2，导入基因组重组信息见表3。

表1变异序列导入变异列表

表2具体导入点突变、插入缺失突变信息表

表3具体导入基因组重组信息表

注：rec1指将120bp(hg19 chr12 25362728-25362848)从载体序列DL196061.1的1270处插入形成的重组基因组；

rec2指将140bp(hg19 chr17 7588908-7589048)从载体序列DL196061.1的2629处插入形成的重组基因组。

(3)生成模拟二代测序数据：利用模拟二代测序软件ART(https://www.niehs.nih.gov/research/resources/software/biostatistics/art/index.cfm)生成DL196061.1载体的二代测序数据，具体使用命令如下：

art_illumina-ss MSv3-i DL196061.1.fasta-l 250-f 50000-

o./DL196061.1-p-m 300-s 100

上述命令为使用软件“ART”生成模拟二代测序数据的具体方法，其中各个词段含义如下：

art_illumina：利用ART软件进行模拟illumina测序平台的二代

测序数据生成；

-ss MSv3：模拟的设备型号为illumina MiSeq；

-i DL196061.1.fasta：对NCBI获取的已知的载体序列DL196061.1.fasta进行模拟测序；

-l 250：生成序列读长为250bp；

-f 50000：生成平均基因组覆盖深度为50000x的测序数据；

-o./DL196061.1：模拟生成的fastq文件输出路径前缀

为./DL196061.1；

-p：模拟进行双端测序；

-m 300：生成平均插入片段长度为300的数据；

-s 100：生成插入片段长度标准差为100的数据。

以同样的方法生成DL196061.1_variant1和DL196061.1_variant2的模拟二代测序数据。

(4)生成低频突变和重组的样本数据：为了后续演示本发明对低频突变和重组的检测效果，我们将DL196061.1_variant1和DL196061.1_variant2的模拟二代测序数据以10％、5％、1％的频率分别混入DL196061.1的模拟二代测序数据中生成模拟低频变异样本数据，具体混入比例如表4。

表4低频突变和重组样本混样比例

(5)测试数据汇总：综上所述，我们使用了NCBI获取的载体序列DL196061.1和2种人工编辑变异序列DL196061.1_variant1和DL196061.1_variant2共3种载体序列。人工引入了4种点突变、2种插入突变、2种缺失突变，共8种突变。并引入了2种基因组重组。生成了3个纯合样本数据(样本命名与来源序列一致，分别称为DL196061.1、DL196061.1_variant1和DL196061.1_variant2)、6个包含低频杂合变异的样本(命名方式见表4)，共9个模拟二代测序数据样本。

后续将以3个纯合样本数据为例说明本发明载体序列组装模块和载体序列和理论序列差异分析模块的实施方式和效果。将以6个包含低频杂合变异的样本数据为例说明本发明载体突变分析模块和载体重组分析模块的实施方式和效果。将以实施例1-4作为优选实施例来详细具体阐述本发明在第一方面和第二方面的所涉及的发明构思，但本发明并不限于实施例1-4。图1示出了本发明一实施例的整体流程图。病毒载体二代测序数据reads经过载体序列组装后得到载体全长序列，再通过载体序列和理论序列差异分析，比较载体全长序列与病毒载体设计的理论序列进行差异分析；通过载体突变分析，比较病毒载体二代测序数据reads与病毒载体设计的理论序列，得出载体突变信息；通过载体重组分析，比对到病毒载体设计的理论序列上，对于reads上不能比对上理论序列的片断，进一步基于宿主细胞参考基因组索引将其与宿主序列进行比对，分析样本的载体序列与宿主基因组序列的重组情况。

实施例1

图2展示了本发明中载体序列组装模块的流程框架，具体流程以上文提到的3个纯合样本模拟二代测序数据为例演示如下：

(1)de Bruijn图算法初步组装

使用基于多尺度de Bruijn图算法的基因组组装工具spades程序(https://github.com/ablab/spades)的RNA病毒组装方案将样本二代短读长序列(reads)组装为长序列重叠群(contigs)。具体方法和参数以DL196061.1样本举例如下：

spades.py--rnaviral-1DL196061.1_R1.fastq.gz-2DL196061.1_R2.fastq.gz-oassembly-t 24

其中各个词段含义如下：

spades.py：基于多尺度de Bruijn图算法的基因组组装程序；

--rnaviral：使用RNA病毒组装方案；

-1DL196061.1_R1.fastq.gz：将DL196061.1样本的二代测序数据的R1端输入组装程序，基于Bruijn图原理进行de novo数据组装；

-2DL196061.1_R2.fastq.gz：将DL196061.1样本的二代测序数据的R2端输入组装程序，基于Bruijn图原理进行de novo数据组装；

-o assembly：输出组装中间文件及组装完成的contigs序列文件到assembly文件夹；

-t 24：为组装程序配置使用的CPU核数为24核。

用同样的方法初步拼接DL196061.1_variant1和DL196061.1_variant2样本的二代测序数据，初步组装得到的contigs数据如表5。可以看到初步组装后没有形成单条完整的基因组，而是形成了2-4个重叠群，组装后的序列覆盖了约98％的实际序列，还有接近2％的实际序列没有覆盖到，没有形成完整的全长。

表5初步组装数据

(2)构建“载体参考序列库”

本次目标载体类型为慢病毒载体(lentivirus vector)从全球性质粒分享平台addgene(https://www.addgene.org/viral-vectors/lentivirus/)的质粒序列库获取lentivirus标签下的所有转移质粒骨架序列。并使用基于k-mer相似性的分子遗传差异性计算工具mash(https://github.com/marbl/mash)计算各个载体骨架序列和初步组装步骤得到的样本contigs序列之间的遗传相关性，过滤掉与样本序列mash遗传距离为1(即相关性为0)的不相关载体的骨架序列，进一步基于mash遗传距离对载体的骨架序列进行排序，最终保留和样本序列分子遗传特征最相似的10个载体的骨架序列作为“载体参考序列库”。本实施例以10个载体的骨架序列作为“载体参考序列库”，在其他实施方案中还可以以5至20个载体的骨架序列作为“载体参考序列库”，优选为8至12个。

本次实施例构建的“载体参考序列库”包含的载体的骨架序列以及与样本contigs序列的遗传距离数据见表6。

表6“载体参考序列库”中的载体序列与样本contigs序列的遗传距离

(3)基因组结构重构

使用基因组重构工具shiver(https://github.com/ChrisHIV/shiver)基于“载体参考序列库”将初步拼接得到的contigs序列和二代测序数据reads重构为单条长慢病毒载体序列(后续称为初步全长序列)。使用shiver重构主要包含程序初始化、contigs比对和reads比对三步，具体方法描述如下：

a.将“载体参考序列库”输入shiver程序进行初始化：shiver程序初始化需要输入比对后的参考序列库，首先使用mafft程序(https://mafft.cbrc.jp/alignment/software/)的E-INS-i(Suitable for sequences with long unalignable regions,veryslow)算法对上文步骤获得的“载体参考序列库”进行多重序列比对，并使用shiver程序的shiver_init.sh命令执行初始化。

b.将样本contigs与“载体参考序列库”进行比对，使用shiver程序的shiver_align_contigs.sh命令将样本contigs与“载体参考序列库”进行比对和裁剪融合，将样本contigs与“载体参考序列库”融合为单条参考序列；

c.将样本二代测序数据reads比对到融合的单条参考序列上。使用shiver程序的shiver_map_reads.sh命令将样本二代测序数据reads比对到上一步获得的融合的单条参考序列上，基于融合的单条参考序列重构获得初步全长序列。并利用比对产生的softclip数据对初步全长序列进行矫正，获得完整的载体全长序列。

基因组结构重构得到的载体全长序列数据如表7。可以看到基因组结构重构后形成了单条完整的载体全长序列。

表7基因组结构重构数据

实施例2

图3展示了本发明中载体序列和理论序列差异分析模块的流程框架，具体流程以上文提到的3个纯合样本模拟二代测序数据为例演示如下：

(1)序列全长比对

利用基因组比对软件MUMmer4(https://mummer4.github.io/)的核酸序列比对程序nucmer对组装得到的载体全长序列和病毒载体设计的理论序列进行比对。具体方法和参数以DL196061.1样本举例如下：

nucmer--maxmatch-p DL196061.1-t 24DL196061.1.fasta DL196061.1_assembled.fasta

其中各个词段含义如下：

nucmer：MUMmer4的核酸序列比对程序；

--maxmatch：不限唯一比对；

-p DL196061.1：输出文件标题为DL196061.1；

-t 24：为比对序程序配置使用的CPU核数为24核；

DL196061.1.fasta DL196061.1_assembled.fasta：将理论载体序列DL196061.1.fasta和载体序列组装模块获得的载体全长序列DL196061.1_assembled.fasta输入程序进行比对得到DL196061.1.delta。

(2)序列比对过滤

利用基因组比对软件MUMmer4(https://mummer4.github.io/)的比对筛选程序delta-filter利用LIS(longest increasing subsequence)算法对nucmer比对结果进行过滤掉短比对，将载体全长序列映射到理论序列的最佳参考位置。具体方法和参数以DL196061.1样本举例如下：

delta-filter-q DL196061.1.delta>DL196061.1.filtered.delta

其中各个词段含义如下：

delta-filter：MUMmer4的比对过滤程序；

-q：将载体全长序列映射到理论序列的最佳参考位置；

DL196061.1.delta：将nucmer分析得到的比对数据DL196061.1.delta输入过滤程序；

>DL196061.1.filtered.delta：将过滤后的比对保存到DL196061.1.filtered.delta文件。

(3)序列全长差异分析

利用基因组比对软件MUMmer4(https://mummer4.github.io/)的核酸序列差异分析程序dnadiff基于过滤后的比对文件分析位点差异、差异片断等。具体方法和参数以DL196061.1样本举例如下：

dnadiff-d DL196061.1.filtered.delta-p DL196061.1.filtered

其中各个词段含义如下：

dnadiff：MUMmer4的核酸序列差异分析程序；

-d DL196061.1.filtered.delta：将过滤后的比对文件输入差异分析程序；

-p DL196061.1.filtered：输出文件标题为DL196061.1.filtered。

本次测试数据的病毒载体设计的理论序列均为DL196061.1.fasta，本步骤检测目标为人工引入的变异。以上用同样的方法分析DL196061.1_variant1和DL196061.1_variant2样本载体全长序列和病毒载体设计的理论序列的差异。

得到位点差异分析结果见表8。可以看到DL196061.1_variant1和DL196061.1_variant2样本结果与实际引入的差异位点相符。

得到差异片断分析结果见表9，可以看到DL196061.1_variant1和DL196061.1_variant2因为引入了宿主重组，与理论序列相比多出了差异片断，差异片断长度与引入的重组片断长度相符。

表8载体全长序列和理论序列的位点差异分析结果

表9载体全长序列和理论序列的差异片断分析结果

实施例3

图4展示了本发明中载体突变分析模块的流程框架，具体流程以上文提到的6个包含低频杂合变异的样本的模拟二代测序数据为例演示如下。

(1)二代测序reads比对到载体序列全长：

利用二代数据比对软件bwa(https://github.com/lh3/bwa)将二代测序reads比对到载体理论序列上，并利用比对数据处理软件samtools(https://github.com/samtools/samtools)将比对数据处理为排序过的bam文件。具体方法和参数以mix_variant1_0.1样本举例如下：

bwa mem\

-M-t 24\

-R

"@RG\tID:mix_variant1_0.1\tPL:ILLUMINA\tSM:mix_variant1_0.1"\

-K 10000000\

DL196061.1.fasta\

mix_variant1_0.1_R1.fastq.gz\

mix_variant1_0.1_R2.fastq.gz|\

samtools view-bS-@24|\

samtools sort-@24-o mix_variant1_0.1.bam

其中各个词段含义如下：

bwa mem：利用二代数据比对软件bwa

(https://github.com/lh3/bwa)的mem算法进行比对；

-M：将同一个reads比对到参考基因组上不同位置的比对标记为次要比对；

--t 24：为比对程序配置使用的CPU核数为24核；

-R

"@RG\tID:mix_variant1_0.1\tPL:ILLUMINA\tSM:mix_variant1_0.1"：reads标题行；

-K 10000000：每批次处理10000000个输入碱基；

DL196061.1.fasta：以载体理论序列DL196061.1.fasta为参考基因组；

mix_variant1_0.1_R1.fastq.gz mix_variant1_0.1_R2.fastq.gz：将mix_variant1_0.1样本的二代测序数据的R1和R2端输入比对程序；

samtools view-bS-@24：使用samtools程序用24个CPU核将比对数据转换为二进制比对数据；

samtools sort-@24-o mix_variant1_0.1.bam：使用samtools程对比对数据基于比对坐标进行排序，并输出为mix_variant1_0.1.bam文件。

(2)突变位点与突变频率分析：

利用突变分析软件VarDict(https://github.com/AstraZeneca-NGS/VarDict)的VarDict程序分析位点突变，并利用VarDict的teststrandbias.R和var2vcf_valid.pl过滤掉低可信度突变，获取vcf表格格式的突变结果。具体方法和参数以mix_variant1_0.1样本举例如下：

VarDict-f 0.001\

-c 1-S2-E 3-g 4-th 24\

-G DL196061.1.fasta-N mix_variant1_0.1-b mix_variant1_0.1.bam\

DL196061.1.fasta.bed|\

teststrandbias.R|\

var2vcf_valid.pl-N mix_variant1_0.1-E-f 0.001\

>mix_variant1_0.1.vcf

其中各个词段含义如下：

VarDict：利用突变分析软件VarDict的VarDict程序进行突变检测；

-f 0.001：过滤掉reads比例在千分之一以下的突变；

-c 1-S2-E 3-g 4：设置输入bed格式为标准格式；

-th 24：为变异分析程序配置使用的CPU核数为24核；

-G DL196061.1.fasta：以载体理论序列DL196061.1.fasta为参考基因组；

-N mix_variant1_0.1：设置样本名为mix_variant1_0.1；

-b mix_variant1_0.1.bam：输入比对程序输出的bam格式比对文件进行突变分析；

mix_variant1_0.1_assembled.fasta.bed：输入突变分析范围；

teststrandbias.R：使用VarDict的teststrandbias.R脚本计算突变的链偏好显著性；

var2vcf_valid.pl-N mix_variant1_0.1-E-f 0.001：使用VarDict的var2vcf_valid.pl脚本对突变数据进行过滤，去掉低可信度突变；

>mix_variant1_0.1.vcf：将过滤后的突变数据保留为vcf表格。

得到的突变和突变频率数据见表10，可以看到检出低频杂合变异与实际混样导入变异相符，检出杂合变异频率与实际混样比例正相关。受算法影响，检出杂合变异频率低于实际混样比例，但是不影响位点检出。

表10低频杂合变异分析结果

实施例4

图5展示了本发明中载体重组分析模块的流程框架，具体流程以上文提到的模拟二代测序数据为例演示如下：

(1)双端数据合并：

利用微生物分析软件vsearch

(https://github.com/torognes/vsearch)的--fastq_mergepairs指令基于重叠序列将二代测序的两端reads合并为单端序列。具体方法和参数以mix_variant1_0.1样本举例如下：

vsearch--fastq_mergepairs\

mix_variant1_0.1_R1.fastq.gz\

--reverse mix_variant1_0.1_R2.fastq.gz\

--fastqout mix_variant1_0.1.merged.fastq\

--threads 24\

--gzip_decompress&&\

gzip mix_variant1_0.1.merged.fastq

其中各个词段含义如下：

vsearch--fastq_mergepairs：利用微生物分析软件vsearch的--fastq_mergepairs指令基于重叠序列将二代测序的两端reads合并为单端序列；

mix_variant1_0.1_R1.fastq.gz：将mix_variant1_0.1样本的二代测序数据的R1输入序列合并程序；

--reverse mix_variant1_0.1_R2.fastq.gz：将mix_variant1_0.1样本的二代测序数据的R2端输入序列合并程序；

--fastqout mix_variant1_0.1.merged.fastq：将合并后的单端序列保存到mix_variant1_0.1.merged.fastq文件；

--threads 24：为合并程序配置使用的CPU核数为24核；

--gzip_decompress：使用gzip对输入文件进行解压缩；

gzip mix_variant1_0.1.merged.fastq：使用gzip对输出文件进行压缩。

(2)序列重组位点分析：

利用基因组重组比对程序ViReMa

(https://sourceforge.net/projects/virema/)将双端合并后的二代测序数据reads比对到载体理论序列上，对于reads上不能比对上载体理论序列的片断，进一步基于宿主细胞参考基因组将其与宿主序列进行比对，分析载体序列与宿主基因组序列的重组情况。具体方法和参数以mix_variant1_0.1样本举例如下：

ViReMa.py\

--Output_Dir mix_variant1_0.1\

--Host_Index ucsc.hg19.fasta\

--p 24\

-BED\

DL196061.1.fasta\

mix_variant1_0.1.merged.fastq.gz\

mix_variant1_0.1.virema.sam

其中各个词段含义如下：

ViReMa.py：利用基因组重组比对程序ViReMa

(https://sourceforge.net/projects/virema/)进行重组比对分析；

--Output_Dir mix_variant1_0.1：将结果输出到mix_variant1_0.1目录；

--Host_Index ucsc.hg19.fasta：将宿主细胞参考基因组输入基因组重组比对程序；

--p 24：为基因组重组比对程序程序配置使用的CPU核数为24核；

-BED：输出bed格式的结果文件；

DL196061.1.fasta：将载体理论序列输入输入基因组重组比对程序；

mix_variant1_0.1.merged.fastq.gz：将双端合并后的二代测序数据reads输入基因组重组比对程序；

mix_variant1_0.1.virema.sam：将重组比对数据保存到mix_variant1_0.1.virema.sam文件。

得到载体重组分析结果见表11。检测到的重组位点与实际引入重组位点相符，发现的重组reads数与实际混入的重组序列比例成正比。

表11载体重组分析结果

第三方面，本发明实施例还提供一种电子设备，包括：处理器以及存储器，处理器与存储器连接；

存储器，用于存储处理器的计算机程序；

其中，处理器被配置为通过执行计算机程序来实现第二方面的方法。

例如，存储器可以包括随机存储器、闪存、只读存储器、可编程只读存储器、非易失性存储器或寄存器等。处理器可以是中央处理器(Central Processing Unit，CPU)等。或者是图像处理器(Graphic Processing Unit，GPU)存储器可以存储可执行指令。处理器可以执行在存储器中存储的可执行指令，从而实现本文描述的各个过程。

可以理解，本实施例中的存储器可以是易失性存储器或非易失性存储器，或可包括易失性和非易失性存储器两者。其中，非易失性存储器可以是ROM(Read-OnlyMemory，只读存储器)、PROM(ProgrammableROM，可编程只读存储器)、EPROM(ErasablePROM，可擦除可编程只读存储器)、EEPROM(ElectricallyEPROM，电可擦除可编程只读存储器)或闪存。易失性存储器可以是RAM(RandomAccessMemory，随机存取存储器)，其用作外部高速缓存。通过示例性但不是限制性说明，许多形式的RAM可用，例如SRAM(StaticRAM，静态随机存取存储器)、DRAM(DynamicRAM，动态随机存取存储器)、SDRAM(SynchronousDRAM，同步动态随机存取存储器)、DDRSDRAM(DoubleDataRate SDRAM，双倍数据速率同步动态随机存取存储器)、ESDRAM(Enhanced SDRAM，增强型同步动态随机存取存储器)、SLDRAM(SynchlinkDRAM，同步连接动态随机存取存储器)和DRRAM(DirectRambusRAM，直接内存总线随机存取存储器)。本文描述的存储器42旨在包括但不限于这些和任意其它适合类型的存储器。

在一些实施方式中，存储器存储了如下的元素，升级包、可执行单元或者数据结构，或者他们的子集，或者他们的扩展集：操作系统和应用程序。

其中，操作系统，包含各种系统程序，例如框架层、核心库层、驱动层等，用于实现各种基础业务以及处理基于硬件的任务。应用程序，包含各种应用程序，用于实现各种应用业务。实现本发明实施例方法的程序可以包含在应用程序中。

在本发明实施例中，处理器通过调用存储器存储的程序或指令，具体的，可以是应用程序中存储的程序或指令，处理器用于执行第二方面所提供的方法步骤。

第四方面，本发明还提供一种计算机可读存储介质，计算机可读存储介质上存储有计算机程序，计算机程序被处理器执行时实现本发明第二方面的方法的步骤。

例如，机器可读存储介质可以包括但不限于各种已知和未知类型的非易失性存储器。

第五方面，本发明实施例还提供一种芯片，用于执行上述第二方面中的方法。具体地，该芯片包括：处理器，用于从存储器中调用并运行计算机程序，使得安装有该芯片的设备用于执行上述第二方面中的方法。

第六方面，本发明实施例还提供一种计算机程序产品，包括计算机程序指令，该计算机程序指令使得计算机执行上述第二方面中的方法。

本领域技术人员可以明白的是，结合本文中所公开的实施例描述的各示例的单元及算法步骤能够以电子硬件、或者软件和电子硬件的结合来实现。这些功能是以硬件还是软件方式来实现，取决于技术方案的特定应用和设计约束条件。本领域技术人员可以针对每个特定的应用，使用不同的方式来实现所描述的功能，但是这种实现并不应认为超出本申请的范围。

在本申请实施例中，所公开的系统、装置和方法可以通过其它方式来实现。例如，单元的划分仅仅为一种逻辑功能划分，在实际实现时还可以有另外的划分方式。例如，多个单元或组件可以进行组合或者可以集成到另一个系统中。另外，各个单元之间的耦合可以是直接耦合或间接耦合。另外，在本申请实施例中的各功能单元可以集成在一个处理单元中，也可以是单独的物理存在等等。

应理解，在本申请的各种实施例中，各过程的序号的大小并不意味着执行顺序的先后，各过程的执行顺序应以其功能和内在逻辑确定，而不应对本申请的实施例的实施过程构成任何限定。

所述功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在机器可读存储介质中。因此，本申请的技术方案可以以软件产品的形式来体现，该软件产品可以存储在机器可读存储介质中，其可以包括若干指令用以使得电子设备执行本申请实施例所描述的技术方案的全部或部分过程。上述存储介质可以包括ROM、RAM、可移动盘、硬盘、磁盘或者光盘等各种可以存储程序代码的介质。

以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术无需创造性劳动就可以根据本发明的构思作出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。

Claims

1.一种基于二代测序技术的病毒载体序列分析系统，其特征在于，包含：

载体序列组装模块，用于通过de Bruijn图算法将样本的二代测序数据reads拼接为contigs序列；还用于通过从addgene获取所有目标类别载体的骨架序列，使用分子遗传相关性计算工具mash计算各个目标类别载体的骨架序列和contigs序列之间的遗传相关性，保留和样本序列分子遗传特征最相似的n个载体的骨架序列作为载体参考序列库；以及还用于将载体参考序列库输入shiver程序进行初始化后与contigs序列进行比对融合为单条参考序列同时生成softclip数据，再将样本的二代测序数据reads比对到融合得到的单条参考序列上重构出初步全长序列，用softclip数据对初步全长序列进行矫正得到完整的载体全长序列；n为5～20；

2.根据权利要求1所述的基于二代测序技术的病毒载体序列分析系统，其特征在于，载体序列和理论序列差异分析模块利用基因组比对软件MUMmer4的核酸序列比对程序nucmer对载体序列组装模块得到的载体全长序列和病毒载体设计的理论序列进行比对；利用MUMmer4的比对筛选程序delta-filter采用LIS算法对nucmer比对结果进行过滤，将载体序列组装模块得到的载体全长序列映射到病毒载体设计的理论序列的最佳参考位置；利用MUMmer4的核酸序列差异分析程序dnadiff基于过滤后的比对文件分析差异位点和差异片断。

3.根据权利要求1所述的基于二代测序技术的病毒载体序列分析系统，其特征在于，还包含载体突变分析模块，用于利用二代数据比对软件bwa将样本的二代测序数据reads比对到病毒载体设计的理论序列上，并利用比对数据处理软件samtools将比对数据处理为排序过的bam文件；还用于利用突变分析软件VarDict的VarDict程序分析位点突变，并利用VarDict的teststrandbias.R和var2vcf_valid.pl过滤掉低可信度突变，获取vcf表格格式的突变结果；该突变结果包括突变位点与突变频率。

4.根据权利要求1所述的基于二代测序技术的病毒载体序列分析系统，其特征在于，还包含载体重组分析模块，用于利用微生物分析软件vsearch的fastq_mergepairs指令基于重叠序列将二代测序的两端reads合并为单端序列；还用于利用基因组重组比对程序ViReMa将样本的二代测序数据reads比对到病毒载体设计的理论序列上，对于reads上不能比对上理论序列的片断，进一步基于宿主细胞参考基因组索引将其与宿主序列进行比对，分析样本的载体序列与宿主基因组序列的重组情况。

5.一种基于二代测序技术的病毒载体序列分析方法，其特征在于，包括以下步骤：

S1、载体序列组装：

S1.2、构建载体参考序列库：a.从addgene获取所有目标类别载体的骨架序列；b.使用分子遗传相关性计算工具mash计算各个载体骨架序列和初步组装步骤得到的样本contigs序列之间的遗传相关性，保留和样本序列分子遗传特征最相似的n个载体的骨架序列作为载体参考序列库；

S1.3、基因组结构重构：a.将载体参考序列库输入shiver程序进行初始化；b.将步骤S1.1得到的contigs序列与载体参考序列库进行比对融合为单条参考序列同时生成softclip数据；c.将样本的二代测序数据reads比对到融合得到的单条参考序列上，重构出初步全长序列，并利用比对产生的softclip数据对初步全长序列进行矫正，获得完整的载体全长序列；n为5～20；

6.根据权利要求5所述基于二代测序技术的病毒载体序列分析方法，其特征在于，包括步骤S2包括：

S2.1、序列全长比对：利用基因组比对软件MUMmer4的核酸序列比对程序nucmer对步骤S1得到的载体全长序列和病毒载体设计的理论序列进行比对；

S2.2、序列比对过滤：利用基因组比对软件MUMmer4的比对筛选程序delta-filter利用LIS算法对nucmer比对结果进行过滤，将步骤S1得到的载体全长序列映射到病毒载体设计的理论序列的最佳参考位置；

S2.3、利用基因组比对软件MUMmer4的核酸序列差异分析程序dnadiff基于过滤后的比对文件分析差异位点和差异片断。

7.根据权利要求5所述基于二代测序技术的病毒载体序列分析方法，其特征在于，还包括以下步骤：

S3、载体突变分析：

S3.1、二代测序reads比对到病毒载体设计的理论序列：利用二代数据比对软件bwa将样本的二代测序数据reads比对到病毒载体设计的理论序列上，并利用比对数据处理软件samtools将比对数据处理为排序过的bam文件；

S3.2、突变位点与突变频率分析：利用突变分析软件VarDict的VarDict程序分析位点突变，并利用VarDict的teststrandbias.R和var2vcf_valid.pl过滤掉低可信度突变，获取vcf表格格式的突变结果；该突变结果包括突变位点与突变频率。

8.根据权利要求5所述基于二代测序技术的病毒载体序列分析方法，其特征在于，还包括以下步骤：

S4、载体重组分析：

S4.1、双端数据合并：利用微生物分析软件vsearch的fastq_mergepairs指令基于重叠序列将二代测序的两端reads合并为单端序列；

S4.2、序列重组位点分析：利用基因组重组比对程序ViReMa将样本的二代测序数据reads比对到病毒载体设计的理论序列上，对于reads上不能比对上理论序列的片断，进一步基于宿主细胞参考基因组索引将其与宿主序列进行比对，分析样本的载体序列与宿主基因组序列的重组情况。

9.电子设备，其特征在于，包括：处理器以及存储器，所述处理器与所述存储器连接；

所述存储器，用于存储所述处理器的计算机程序；

其中，所述处理器被配置为通过执行所述计算机程序来实现权利要求5至8任一项所述的方法。

10.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质上存储有计算机程序，所述计算机程序被处理器执行时实现如上权利要求5至8中任一项所述的方法的步骤。