技术领域
[0001] 本发明涉及生物信息领域，具体地，本发明涉及检测DNA分子序列突变的方法及装置，更具体地，本发明涉及检测DNA分子序列突变的方法及装置、计算机程序产品、计算设备以及计算机可读存储介质。
相关背景技术
[0002] DNA测序技术是遗传学和生物医学研究领域重要的分析工具之一。随着现代高通量测序技术的发展，基因组信息能够被更深入的解析，并通过比较个体之间的DNA序列差异来鉴定疾病或进行新药开发等。
[0003] 然而，在进行分子测序时存在许多误差来源，例如PCR扩增、DNA片段化和测序仪器本身的误差等。这些误差可能导致无法准确地检测出个体之间的细微差异，进而影响疾病诊断和治疗方案的制定或误判新药研发靶点。
[0004] 为了解决这些问题，科学家们已经提出了一系列的分子检测方法，以便更加准确地检测DNA序列变异。如利用一代测序技术等。但目前还检测变异准确率较低，仍有较大发展空间。
具体实施方式
[0065] 下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。
[0066] 定义和说明
[0067] 在本文中，除非另有说明，单数形式“一种”、“一个”等包括复数指代物(一个以上)；“一组”或者“多个”指两个或两个以上。
[0068] 在本文中，除非另有说明，术语“包含”或“包括”为开放式表达，即包括本发明所指明的内容，但并不排除其他方面的内容。
[0069] 在本文中，所称的“测序”为序列测定，同“核酸测序”或“基因测序”，指核酸序列中碱基次序的测定。本发明实施例采用的是Nanopore测序，利用一个纳米孔将DNA或RNA单链引导通过，同时监测电流信号变化，从而实现对DNA/RNA序列的测序。相比于传统的基于荧光信号的测序技术，Nanopore测序具有优势，如速度更快、成本更低、操作更简便等。此外，它还可以直接测序RNA，避免了转录过程中可能出现的错误。
[0070] 在本文中，所述ONT接头是指用于连接DNA序列到Nanopore技术测序芯片的短DNA分子。
[0071] 需要说明的是，本发明实施例中所述命令的参数可根据序列的需求进行更改，例如，“canu‑pVP64‑P‑d VP64‑P genomeSize＝11k useGrid＝false maxThreads＝8maxMemory＝10G‑nanopore VP64‑P.1k.fq”，所述genomeSize、useGrid、maxThreads、maxMemory等参数以及文件的命名(VP64‑P.1k.fq)均可做适应性修改。
[0072] 检测DNA分子序列突变的方法
[0073] 在本发明的一方面，本发明提出了一种检测DNA分子序列突变的方法，所述方法包括：获取所述DNA分子的测序结果，所述测序结果由多个测序读段构成；将所述测序读段进行校正处理，以便获得校正测序读段；和对所述校正测序读段进行序列突变分析处理，以便获得所述DNA分子序列突变信息。
[0074] 具体的，为了方便理解，下面对本申请的技术方案(以质粒Nanopore文库构建至序列组装为例)进行详细解释和说明，所述技术方案具体步骤如下：
[0075] 1)文库构建：将环状结构的质粒用TN5酶进行酶切打断，并加标签；Pooling及纯化；加Nanopore接头；
[0076] 2)Nanopore测序：测序步骤不局限于特定的模式和平台，例如，可选则性使用Fast、HAC和Super三种测序模式；
[0077] 3)去除接头序列：使用可以去除Nanopore测序数据中的接头序列的软件或程序执行。例如，可以使用Porechop程序(https://github.com/rrwick/Porechop)来去除Nanopore测序数据中的接头序列；
[0078] 4)测序读段筛选处理(可选)：选取平均质量值最高的50倍及以上质粒长度的测序读段数据，用于后续校正处理和组装处理。被选取的数据可以是随机选取的，也可以平均质量值最高的；但选取的测序读段的平均质量值不能低于10倍以下；
[0079] 5)测序读段校正处理。用canu软件对步骤4)选取数据进行校正，去除其中的错误和噪音，然后再进行组装，以提高组装结果的准确性和连续性；
[0080] 6)测序读段切割处理(可选)：将步骤5)得到的校正测序读段按一定阈值切割为多个子读段，阈值范围为300bp～30000bp。本步骤可选择性采取无交集切割及有交集切割(交集的阈值一般选取200bp以下)两种方法实现；
[0081] 7)序列组装(可选)：使用Newbler软件对步骤6)得到的子读段进行组装，获得组装结果；
[0082] 8)目标序列获取(可选)：比对软件可使用blastn，本步骤需要知道引物序列或者载体序列，如无则可不执行此步骤。本步骤主要用于解决样本是否存在污染，若存在污染，则步骤7)就会组装出多条contig序列，通过引物序列或载体序列和组装结果比对，即可确定目标contig序列。若样品不存在污染，组装结果一般只有1条contig。
[0083] 9)参考序列建库；
[0084] 10)测序读段与参考序列库比对：所述测序读段为校正后测序读段；
[0085] 11)比对结果格式转换并进行排序；
[0086] 12)测序读段序列突变分析。
[0087] 最终可获得如图5所示的信息文件，包括conitg ID(第一列)，碱基位置(第二列)，具体碱基(第三列)，位点覆盖深度(第四列)，每条序列比对信息(第五列)，质量值(第六列)。
[0088] 需要说明的是，上述完整突变检测方法适于无参考序列分析，如果有参考序列，可对上述步骤进行适应性删减。
[0089] 有益效果
[0090] 基于上述方法进行Nanopore文库构建、文库测序、测序数据组装、突变信息检测等处理，能够准确组装，同时通过比对和变异检测，可得到每个碱基的覆盖深度和变异信息，提高检测结果的准确率和检测范围，用时较短，且花费成本较低。
[0091] 下面将更详细地描述本发明的实施例，所述实施例的示例在附图中示出。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。实施例中未注明具体技术或条件的，按照本领域内的文献所描述的技术或条件或者按照产品说明书进行。所用试剂或仪器未注明生产厂商者，均为可以通过市购获得的常规产品。
[0092] 实施例1
[0093] 本实施例展示文库构建实验流程：
[0094] 1)Tagmentation(标签化)
[0095] 按照表1所示进行反应体系配制，充分混匀后，置于PCR仪上反应(反应条件如表2所示)。
[0096] 表1：反应体系
[0097]试剂 用量(μL)
50ng DNA模版 9
快速条码(Rapid Barcodes,RB01‑96,one for each sample) 1
总体积 10
[0098] 表2：反应条件
[0099]温度(℃) 时间(min)
30 2
80 2
4 ∞(无穷)
[0100] 2)Pooling(混合样本)
[0101] 将包含不同barcode样本的反应液汇集到新的1.5mL离心管中，形成一个样本。
[0102] 3)磁珠纯化
[0103] ①将SPRI beads解冻并进行重悬混匀；
[0104] ②向pooling后的样本中加入等体积重悬好的SPRI beads，充分混匀，置于常温混匀机上5min；
[0105] ③将离心管短暂离心并置于磁力架上，待溶液澄清后，小心移除上清；
[0106] ④保持离心管始终处于磁力架上，加入1.5mL新鲜配置的80％乙醇漂洗磁珠，室温孵育30秒后小心移除上清；
[0107] ⑤重复上一步骤一次，总计漂洗两次；
[0108] ⑥吸除任何残留的80％乙醇，开盖短暂晾干，将离心管从磁力架上取出，加入15μL Elution Buffer(洗脱缓冲液)轻轻吹打混匀，常温放置10min；
[0109] ⑦将离心管短暂离心并置于磁力架上，待溶液澄清后，小心吸取全部上清至新的1.5mL EP管中；
[0110] ⑧取1μL纯化产物进行Qubit检测。
[0111] 4)加ONT(Oxford Nanopore Technologies)接头
[0112] 取11μL纯化产物到新的1.5mL离心管中，加入1μL快速接头F(Rapid Adapter F，RAP F)，吹打混匀后离心，室温孵育5min。
[0113] 对比例1
[0114] 本实施例利用常用方法对上述Nanopore混样测序结果进行组装(canu软件)。以VP64‑P、VP64‑C、KRAB‑P质粒为示例，其中，质粒1(VP64‑P)、质粒2(VP64‑C)、质粒3(KRAB‑P)的参考序列已知，长度分别为12,133bp、11,307bp、11,397bp。具体步骤如下：
[0115] 1)质粒1(VP64‑P)、质粒2(VP64‑C)、质粒3(KRAB‑P)按照实施例1的建库方法完成文库构建，三个样品分别对应barcode1、barcode2和barcode3；
[0116] 2)选择Nanopore的Super测序模式，用型号PromethION型号的测序仪对文库进行测序；
[0117] 3)测序数据产量到达100MB后停止测序。获得3个样本的数据量如表3所示；
[0118] 表3：样本数据产量统计
[0119]样本名称 数据量(MB)
VP64‑P 39MB
VP64‑C 31MB
KRAB‑P 30MB
[0120] 4)使用Porechop程序来去除Nanopore测序数据中的接头序列。命令为“porechop‑iraw.fq‑o clean.fq”，其中，raw.fq表示原始下机数据，clean.fq表示去除掉接头后的数据，为了方便后续描述，三个样品去完接头后的文件分别命名为“VP64‑P.clean.fq”、“VP64‑C.clean.fq”、“KRAB‑P.clean.fq”；
[0121] 5)分别对“VP64‑P.clean.fq”、“VP64‑C.clean.fq”、“KRAB‑P.clean.fq”每条序列的平均质量值从高到低排序，并分别取平均质量值最高的1000条序列用于后续的校正和组装，分别命名为“VP64‑P.1k.fq”、“VP64‑C.1k.fq”、“KRAB‑P.1k.fq”；
[0122] 6)用canu软件分别对“VP64‑P.1k.fq”、“VP64‑C.1k.fq”、“KRAB‑P.1k.fq”进行校正及组装，分别获得组装结果“VP64‑P.contigs.fasta”、“VP64‑C.contigs.fasta”和“KRAB‑P.contigs.fasta”；
[0123] VP64‑P的校正和组装命令为“canu‑p VP64‑P‑d VP64‑P genomeSize＝11k useGrid＝false maxThreads＝8maxMemory＝10G‑nanopore VP64‑P.1k.fq”；
[0124] VP64‑C的校正和组装命令为“canu‑p VP64‑C‑d VP64‑C genomeSize＝11k maxThreads＝8useGrid＝false maxMemory＝10G‑nanopore VP64‑C.1k.fq”；
[0125] KRAB‑P的校正和组装命令为“canu‑p KRAB‑P‑d KRAB‑P genomeSize＝11k useGrid＝false maxThreads＝8maxMemory＝10G‑nanopore KRAB‑P.1k.fq”；
[0126] 组装后结果如表4所示，可以明显看出，组装结果长度于实际长度差异较大。
[0127] 表4：三个质粒实际长度与组装长度统计
[0128]
[0129]
[0130] 实施例2
[0131] 本实施例基于本发明技术方案对相同样本(样本通对比例1)进行组装。其中，质粒1(VP64‑P)、质粒2(VP64‑C)、质粒3(KRAB‑P)的参考序列已知，长度分别为12,133bp、11,
307bp、11,397bp。本实施例步骤1)到5)和对比例1的建库和测序方法相同，此处不再赘述，后面的实施步骤如下所示：
[0132] 6)用canu软件分别对“VP64‑P.1k.fq”、“VP64‑C.1k.fq”、“KRAB‑P.1k.fq”进行校正。VP64‑P的校正和组装命令为“canu‑correct‑p VP64‑P‑d VP64‑P genomeSize＝11k useGrid＝false maxThreads＝8maxMemory＝10G‑nanopore VP64‑P.1k.fq”，获得校正后的序列“VP64‑P.correctedReads.fasta.gz”；
[0133] VP64‑C的校正和组装命令为“canu‑correct‑p VP64‑C‑d VP64‑C genomeSize＝11k maxThreads＝8useGrid＝false maxMemory＝10G‑nanopore VP64‑C.1k.fq”。获得校正后的序列“VP64‑C.correctedReads.fasta.gz”；
[0134] KRAB‑P的校正和组装命令为“canu‑correct‑p KRAB‑P‑d KRAB‑P genomeSize＝11kuseGrid＝false maxThreads＝8maxMemory＝10G‑nanopore KRAB‑P.1k.fq”，获得校正后的序列“KRAB‑P.correctedReads.fasta.gz”；
[0135] 7)以1000bp为单位对校正后的序列进行无交集切割为示例(同样也可进行有交集切割)，分别得到“VP64‑P.split.fa”、“VP64‑C.split.fa”和“KRAB‑P.split.fa”三个文件；
[0136] 8)用newbler软件分别对“VP64‑P.split.fa”、“VP64‑C.split.fa”和“KRAB‑P.split.fa”三个文件进行组装，分别得到文件名为“454AllContigs.fna”的组装结果，为方便后续描述，我们将三个组装结果文件分别命名为“VP64‑P.asm.fa”、“VP64‑C.asm.fa”和“KRAB‑P.asm.fa”。这三个组装结果均只有1条contig；
[0137] 9)分别将三个质粒的引物序列用blastn软件比对到组装结果。结果是三个引物序列均能比对上组装结果；
[0138] 10)进一步的，为了验证本发明方法的有效性，分别将3个组装结果(VP64‑P.asm.fa、VP64‑C.asm.fa和KRAB‑P.asm.fa)与参考序列(VP64‑P.fa、VP64‑C.fa、KRAB‑P.fa)比对，检验两者的相似性。VP64‑P比对命令“blasr VP64‑P.asm.fa VP64‑P.fa‑out VP64‑P.m4‑m 4”，并将VP64‑P.asm.fa未比对上的区域提取出来，再与VP64‑P.fa比对，之后统计比对结果。其余两个样本进行类似操作。比对结果如表5所示，相似度均大于99.9％。
[0139] 表5：组装结果与参考序列相似度统计
[0140]
[0141]
[0142] 其中，上述三个质粒实际长度与本实施例检测的长度对比统计如表6所示，可以看出，仅相差几bp。
[0143] 表6：三个质粒实际长度与本实施例的组装长度统计
[0144]样本名 实际长度 组装长度
VP64‑P 12,133bp 12,130bp
VP64‑C 11,307bp 11,303bp
KRAB‑P 11,397bp 11,393bp
[0145] 综上所述，使用本发明实施例从实验建库、测序及信息分析，总花费时间小于6小时。相比Sanger测序花费的时间能大幅的缩减，准确率也有明显的提高。
[0146] 实施例3
[0147] 本实施例基于实施例2步骤6)中获得的校正后序列进行突变信息检测。利用如下命令得到样本VP64‑P的变异信息统计文件：(1)bwa index VP64‑P.asm.fa；(2)bwa mem‑x ont2d VP64‑P.asm.fa VP64‑P.correctedReads.fasta.gz>VP64‑P.bwa.sam；(3)samtools view‑bS VP64‑P.bwa.sam>VP64‑P.bwa.bam；(4)samtools sort VP64‑P.bwa.bam‑o VP64‑P.sort.bam；(5)samtoolsmpileup‑fVP64‑P.asm.fa VP64‑P.sort.bam>VP64‑P.pileup.txt。
[0148] 样本VP64‑C和KRAB‑P进行类似操作，得到变异信息统计文件VP64‑C.pileup.txt和KRAB‑P.pileup.txt。
[0149] 如图6所示，利用样本KRAB‑P的校正序列对参考序列做变异检测，结果发现，参考序列第737碱基位点G，有43条校正后的序列覆盖，其中，有22条序列支持这个位点是C，21条支持是G，表明该位点为突变位点。而第746碱基位点A，有41条校正后的序列覆盖，其中，有41条序列支持这个位点是G，表明该位点也是一个突变位点。
[0150] 如图7所示，利用样本KRAB‑P的校正序列对参考序列做变异检测，结果发现，参考序列第10197碱基位点A，有44条校正后的序列覆盖，其中有22条序列支持这个位点是A，22条序列支持该位点插入CGCG四位碱基。
[0151] 如图8所示，利用样本KRAB‑P的校正序列对参考序列做变异检测，结果发现，参考序列第754碱基位点G，有41条校正后的序列覆盖，有41条序列支持缺失“GCT”3个碱基。
[0152] 应当理解，本发明的各部分可以用硬件、软件、固件或它们的组合来实现。在上述实施方式中，多个步骤或方法可以用存储在存储器中且由合适的指令执行系统执行的软件或固件来实现。例如，如果用硬件来实现，可用本领域公知的下列技术中的任一项或他们的组合来实现：具有用于对数据信号实现逻辑功能的逻辑门电路的离散逻辑电路，具有合适的组合逻辑门电路的专用集成电路，可编程门阵列(PGA)，现场可编程门阵列(FPGA)等。
[0153] 本技术领域的普通技术人员可以理解实现上述实施例方法携带的全部或部分步骤是可以通过程序来指令相关的硬件完成，所述的程序可以存储于一种计算机可读存储介质中，该程序在执行时，包括方法实施例的步骤之一或其组合。
[0154] 此外，在本发明各个实施例中的各功能单元可以集成在一个处理模块中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个模块中。上述集成的模块既可以采用硬件的形式实现，也可以采用软件功能模块的形式实现。所述集成的模块如果以软件功能模块的形式实现并作为独立的产品销售或使用时，也可以存储在一个计算机可读取存储介质中。
[0155] 在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。
[0156] 尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在不脱离本发明的原理和宗旨的情况下在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。
