CN115369159A - 一种基于双端测序重叠片段和dna双链互补片段的超低频突变检测方法 - Google Patents

Abstract

本发明公开了一种基于双端测序重叠片段和DNA双链互补片段的超低频突变检测方法，包括选择不依赖于PCR扩增的方法构建DNA文库，结合使用TE缓冲液、缩短片段化时间进行基因组DNA片段化和使用特异性酶末端修复减少测序错误，缩短文库片段长度产生双端测序重叠片段结合互补链信息获得DNA模板的拷贝片段进行测序错误校正，选择内源性标签对拷贝片段进行标记，根据互补链来源的双端读段在基因组上比对位置相同、比对方向相反提取拷贝片段，进行序列信息比对，保留一致的碱基信息用于突变分析。本发明构建的测序方法有效降低NGS突变检测的错误率，针对哺乳动物全基因组实现高准确性、低测序成本的超低频突变检测。

Description

一种基于双端测序重叠片段和DNA双链互补片段的超低频突 变检测方法

技术领域

本发明属于分子生物学领域，具体涉及一种基于双端测序重叠片段和DNA双链互补片段的超低频突变检测方法。

背景技术

DNA测序是识别基因组序列信息、解决大量生命科学领域问题的关键手段。自2005年以来，第二代测序技术(Next-generation sequencing，NGS)以Illumina Solexa和HiSeq技术、Roche 454技术、ABI SOLiD技术以及Life Technologies Ion Torrent等技术为代表，逐步实现了自动化的高通量测序，使得低成本的大规模基因组变异检测成为可能。经过近十几年的不断发展，NGS以其高通量和低成本，已然成为目前生命科学研究的常规技术手段，广泛应用于肿瘤研究、宏基因组学分析、进化研究和产前诊断等诸多领域。

然而，NGS的随机测序错误率很高，约为10^-3-10^-2；也就是说，平均每测序10²-10³个碱基，就会有1个碱基被错误地识别为其他碱基。遗传于父母的胚系突变及生命早期的体细胞突变会存在于绝大多数体细胞中，NGS的高错误率并不会影响对此类突变的检测。在肿瘤组织等样本中，由于存在细胞的克隆性增殖，相同突变会在同一个克隆来源的细胞群中存在；当这些突变的比例(即，变异等位基因频率)大于1％时，随机测序错误的影响可以被忽略，通过提高NGS测序深度也能够获得比较准确的突变信息(真实的变异等位基因频率＞～1％＞随机测序错误产生的变异等位基因频率)。而在正常组织细胞中，自发突变的发生率很低。以人类基因组为例，细胞每分裂一次碱基发生突变的概率约为10^-9-10^-8。同时，绝大多数自发突变并不会引起细胞表型和功能的改变，携带有这些突变的细胞其扩增往往是有限的。因此，在正常的组织细胞中，突变发生率低并且其变异等位基因频率通常远低于1％，属于“超低频突变”。此时，在常规NGS 10^-3-10^-2的测序错误率的条件下，难以对上述超低频突变进行准确的检测。

目前，有几种基于NGS的策略或技术被应用于超低频突变的检测，包括：(1)单细胞测序技术(分离单个细胞，进行全基因组扩增并测序)；(2)体外培养单克隆来源的细胞利用NGS进行突变检测；(3)使用微量组织(含有较高比例的克隆性增殖的细胞群)或者利用显微技术根据组织结构分离单克隆来源的细胞群进行突变检测以及(4)分子一致性测序(Molecular consensus sequencing)策略。前面三种方法主要是利用单细胞来源的序列信息(扩增片段或者单克隆来源的细胞)进行超低频突变检测。而分子一致性测序策略则是通过直接降低NGS进行碱基识别的错误率，实现对基因组片段进行直接的突变检测；其更适用于群体细胞的突变检测，以了解整体水平上的突变情形。

分子一致性测序策略降低NGS进行碱基识别错误率的基本原理是：利用同一待测DNA模板多个拷贝的测序信息(即“对模板进行多次重复测序”)校正随机测序错误，从而降低NGS进行碱基识别的错误率(图1)。在分子一致性测序方法中，通常在文库制备阶段使用PCR扩增来获得待测DNA模板的多个拷贝片段，并对这些拷贝片段分别进行测序，读取序列信息；随后利用拷贝片段携带的从模板来源的相同的分子标签(Molecular barcodes)的信息识别和标记这些从同一模板来源的多个拷贝片段的测序信息。理论上，这些同一模板来源的拷贝片段是同质的，其通过测序得到的序列信息应当是完全一致的；由于NGS的测序错误的存在，其中会存在少量不一致的变异信息(图1)。NGS的测序错误是随机出现的，因而尽管其发生率较高，但是在多个独立测序片段的相同位点出现相同类型的测序错误的概率是很低的。通过对从同一模板来源的多个拷贝片段的测序信息进行比对，在大多数片段中一致的碱基信息被认为反映了模板链上真实的碱基信息；而在少数片段中存在的不一致的碱基信息被认为是随机测序错误从而去除。通过这样的错误校正策略，分子一致性测序可以去除大量因测序错误导致的变异位点，得到与模板信息高度一致的序列信息用于突变检测。

近年来，基于分子一致性测序策略已经产生了许多改良的高准确性的NGS突变检测方法，包括：Safe-Sequencing System(Safe-SeqS)、Circle sequencing、Bottlenecksequencing system(BotSeqS)、Duplex sequencing(DupSeq)、Hypothesis alignmentwith weak overlap(Hawk-Seq^TM)、Jade-Seq^TM、CarcSeq和Nanorate sequencing(NanoSeq)等。前期工作已进行类似方法(以下简称“PECC-Seq^*”)的开发，即：利用双端测序重叠片段及DNA双链互补片段进行随机测序错误的校正。这些方法能够将NGS进行DNA测序的错误率降低至10^-9-10^-4，其中一些方法的错误率已经接近哺乳动物基因组的自发突变水平10^-9-10^-7，能够初步实现正常组织细胞上的超低频突变的检测。

超低频突变检测的准确性是应用分子一致性测序改良方法时需要考虑的关键问题。由于哺乳动物基因组自发突变率极低，突变检测方法的准确性越高(即，进行碱基识别的错误率越低)，所能检测到的正常组织细胞基因组的突变频率和突变特征越准确。理论上，要实现对于正常组织细胞基因组自发突变水平上的超低频突变检测，所需的突变检测方法的错误率需要接近甚至低于自发突变的发生率，即10^-9-10^-7。目前，除NanoSeq方法的错误率接近10^-9-10^-8外，BotSeqS、DupSeq、Hawk-Seq^TM以及PECC-Seq^*的错误率在～10^-7水平，仍然是略高于基因组自发突变水平；而其他一些方法的错误率则要高于10^-7水平(表1)。因此，绝大多数的分子一致性测序改良方法的准确性仍然不足以真正实现准确的基因组自发突变检测，准确性仍有待进一步提高。

表1：PECC-Seq^*与代表性的分子一致性测序改良方法比较

此外，与常规的NGS检测相比，分子一致性测序改良方法其所需的测序深度和测序成本相对较高、测序深度(测序成本)-数据产出的效率(以下简称“测序效率”)较低。在常规NGS中，每一个模板理论上至多只会被检测一次。而在分子一致性测序分析中，同一个模板来源的多个拷贝需要被测序以对模板序列信息进行测序错误校正。因而，分子一致性测序分析所需的测序深度往往比常规NGS要高很多、其测序效率较低，由此导致其测序成本非常昂贵。对于微型基因组(例如，微生物基因组)和目标基因片段的突变检测而言，分子一致性测序分析的测序效率和产生的测序成本升高尚可接受；但对于哺乳动物全基因组(例如，人类基因组和大小鼠基因组)突变检测而言，分子一致性测序分析所需的测序深度及其产生的测序成本限制了这类方法的广泛应用。目前，仅有少数几种方法，包括：NanoSeq、BotSeqS、Hawk-Seq^TM和PECC-Seq^*应用于哺乳动物全基因组范围的突变检测(表1)。这几种方法中，BotSeqS方法平均每检测～8700个碱基才能获得1个校正后的碱基信息；PECC-Seq^*方法平均每检测～500个碱基能够获得1个校正后的碱基信息(IlluminaHiSeqXTen测序平台)；尽管NanoSeq测序效率较高，平均每检测～30个碱基就获得1个校正后的碱基信息，但其在文库制备时使用内切酶进行片段化，因此最多只能覆盖～30％的基因组位点，并不能完全反映全基因组上的突变信息。DupSeq作为一种已经成熟的商业化的突变检测方法，主要适用于较小的基因组片段(微型基因组或目标基因)的突变检测。

现有的分子一致性测序改良方法在准确性、测序成本和在哺乳动物全基因组突变检测的适用性上都各有优缺点(表1)，PECC-Seq^*方法在现有的类似方法中准确性(检测到的基因组背景突变频率越接近于基因组自发突变频率水平，表明方法的错误率越低、准确性越高)和测序效率都相对较高，仅逊于目前NanoSeq方法，且能用于大型基因组全基因的突变检测；但其错误率仍然要略高于正常组织细胞基因组的自发突变频率，同时在测序效率上也有进一步提高的空间，仍然有必要在此基础上进一步优化和开发，实现更准确和效率更高的全基因组超低频突变检测。

发明内容

基于此，本发明的主要目的是提供一种基于双端测序重叠片段和DNA双链互补片段的超低频突变检测方法，基于分子一致性测序策略得到，降低NGS进行突变检测的错误率，实现高准确性、低测序成本的哺乳动物全基因组超低频突变检测。

为实现上述目的，本发明采用如下技术方案：

本发明提供的一种基于双端测序重叠片段和DNA双链互补片段的超低频突变检测方法(以下简称PECC-Seq方法)，包括以下步骤：

(1)以不依赖于PCR扩增方法(PCR-free)构建DNA文库，同时优化DNA文库制备条件，包括结合使用TE缓冲液、缩短片段化时间进行基因组DNA片段化和使用特异性酶进行末端修复来减少文库制备过程引入大量的末端修复错误导致的测序错误；

(2)缩短文库片段长度获得双端测序重叠片段结合互补链信息获得DNA模板的拷贝片段进行测序错误校正，包括：缩短文库插入片段的长度至～150bp，并先通过Illumina测序的2×150bp双端测序模式获得单链模板自然来源的两个拷贝，即双端测序的重叠片段；再用双链DNA模板互补链的一组拷贝片段获得非PCR来源的同一双链DNA模板来源的四条拷贝片段进行测序错误校正；

(3)选择内源性标签对拷贝片段进行标记；

(4)数据处理和突变分析，包括：

根据双端读段在基因组的比对信息进行同一模板来源的拷贝片段提取：互补链来源的双端读段在基因组上的比对位置相同、比对方向相反；

对提取的拷贝片段进行序列信息比对：当同一位点存在不一致碱基时，提示存在测序错误，则该位点被剔除；当且仅当同一位点序列信息一致，保留一致的碱基信息用于突变分析。

作为优选，上述PECC-Seq方法适用于哺乳动物基因组的全基因组突变检测，包括人类基因组和小鼠基因组等，测序深度≤40×。

作为优选，步骤(1)中，所述TE缓冲液的组成为10mM Tris-HCl,1mM EDTA-Na₂，pH8.0。

作为优选，步骤(1)中，使用Covaris Focused-ultrasonicator S220非接触式超声波破碎仪进行片段化的条件为：最大发射功率175W，工作系数10％，循环数200，水位12，片段化时间140s及温度7℃。

作为优选，步骤(1)中，所述特异性酶选自单链特异性核酸酶、DNA损伤修复酶或DNA聚合酶；其中所述单链特异性核酸酶包括S1核酸酶，所述DNA损伤修复酶包括Fpg酶，所述末端修复酶包括T4 DNA聚合酶。

作为优选，步骤(4)中，剔除分布于插入片段末端10bp区域内的碱基信息不用于数据分析。

作为优选，步骤(4)中，还包括给出提取拷贝片段、一致性碱基和突变分析过程的步骤。

本发明前述PECC-Seq方法用于超低频突变检测的原理在于：PECC-Seq基于

DNA PCR-Free文库制备流程进行文库制备改良，通过调整片段化条件及文库片段筛选的条件获得目标长度(～150bp)的DNA双链模板片段，同时添加特异性核酸酶对模板片段进行末端修复减少单链损伤固定的背景测序错误，获得插入片段长度～150bp的PCR-free的测序文库(图10A)。文库片段利用Illumina2×150bp双端测序模式进行测序，插入片段的长度与测序的读长接近，双端测序产生的R1和R2的序列信息会出现很大的重叠(图10B)，这部分重叠序列均反映同一单链模板的信息，是模板的两个拷贝。按图10和图11所示的PECC-Seq数据分析流程对原始的测序数据进行序列比对和数据过滤后，根据双端读段的比对信息(比对位置和比对方向)进行互补链来源的双端读段的提取，互补链来源的双端读段在基因组上的比对位置相同而比对方向相反，比对提取的同一模板来源的拷贝片段(互补链来源的双端读段的重叠序列)，仅保留同一位点碱基信息一致并且分布于文库插入片段内部(距片段末端＞10bp的位点)的碱基信息用于突变分析。

本发明还提供所述基于双端测序重叠片段和DNA双链互补片段的超低频突变检测方法在哺乳动物全基因组超低频突变检测中的应用。

作为优选，所述哺乳动物基因组包括人类基因组和小鼠基因组。

与现有技术相比，本发明的有益效果在于：

(1)构建PCR-free的DNA文库。现有的近似方法均选择依赖于PCR扩增的文库制备方法，即以PCR扩增为基础获得同一模板的多个拷贝片段用于NGS的随机测序的校正。PCR扩增过程本身就会向模板序列中引入错误，导致测序结果中的变异(测序错误)，并且PCR扩增过程固有的不均匀性，使得依赖于PCR扩增的文库制备方法得到的测序结果对基因组的覆盖度要低于不依赖于PCR扩增(PCR-free)的文库制备方法。此外，由于PCR扩增过程的不均匀性，很难控制用于测序错误校正所使用的拷贝片段的数目；其往往远多于理论上所需要的片段数目。本发明选择PCR-free的文库制备方案构建的检测方法对基因组有更好的覆盖度、能够减少PCR扩增过程引入的测序错误，并且能更好的控制用于错误校正的拷贝片段数目。

(2)互补链结合双端测序重叠片段进行测序错误校正。在分子一致性测序分析中，由随机测序错误和模板单链损伤导致的背景测序错误发生率的理论值可以分别由公式1和公式2计算得到；即至少需要4条来源于互补链的拷贝片段(单链的2条拷贝片段×1对互补链)进行错误校正，可使得分子一致性分析中由于随机测序错误((10^-2÷3)⁴×3＝3.7×10^-10)和模板单链损伤(10^-4×((10^-2÷3)²×3)＝3.3×10^-9，单链损伤导致的单链一致性测序错误的发生率以10^-4计算)导致的背景测序错误的发生率降低至自发突变水平。在PCR-free的文库制备基础上，要获得上述足够的模板的拷贝片段(互补链来源的至少4条拷贝片段)进行错误校正，利用双端测序的重叠片段。在正常的测序文库中插入片段平均长度≥350bp，经过2×150 bp双端测序后，其产生的双端读段R1和R2并不会出现重叠的序列信息。本发明通过调整插入片段长度至～150 bp，获得双端读段R1和R2直接的重叠的序列信息，即模板的两个拷贝，结合互补链的信息得到符合要求的足够的拷贝片段用于测序错误的校正。

校正后的测序错误率＝(测序平台的随机测序错误率÷3)^{拷贝片段数目}×3

(公式1)

基于互补链的一致性测序中单链损伤导致的测序错误的发生率

＝单链损伤导致的单链一致性测序错误的发生率×((测序平台的随机测序错误率÷3)^{互补链的拷贝片段数目}×3)

(公式2)

(3)互补链结合排他性扩增重复片段进行单链一致性测序分析。利用构建得到的双端测序重叠片段结合发现的数据中排他性扩增重复片段获得从同一单链模板来源的4条拷贝片段，由此进行基于单链的一致性测序分析，从而获得文库制备及测序过程中可能引入的单链损伤及单链变异的信息。本发明结合PECC-Seq和基于单链的一致性分析方法(PEEA-Seq)对测序过程中的错误来源进行系统分析。

(4)给出相应的数据分析步骤，利用该流程完成拷贝片段的提取、一致性碱基的提取和突变分析的过程，准确高效。

附图说明

图1是分子一致性测序策略的技术原理图。

图2是不同分子一致性测序方法的技术原理图比较。

图3是PECC-Seq^*方法的技术原理图。

图4是本发明中PEEA-Seq方法的技术原理图。

图5是本发明中PECC-Seq^*方法中背景测序错误产生原因分析。

图6是不同片段化强度和缓冲液下PECC-Seq^*方法检测到的背景突变频率。

图7是不同片段化时间下PECC-Seq^*方法检测到的背景突变频率。

图8是不同超声作用平均强度下PECC-Seq^*方法检测到的背景突变频率。

图9是不同文库制备条件下PECC-Seq^*方法检测到的背景突变频率。

图10是本发明中PECC-Seq方法用于突变检测的技术原理图。

图11是本发明中PECC-Seq方法数据分析流程示意图。

图12是实施例中AAI处理的转基因小鼠肾脏组织不同类型碱基置换突变的突变频率。

图13是实施例中AAI处理的转基因小鼠肾脏组织的突变谱。

图14是实施例中PECC-Seq^*方法检测得到的背景突变在插入片段上的位置分布。

图15是实施例中PECC-Seq^*方法检测得到的TK6细胞的背景突变谱(去除末端7bp)。

图16是实施例中PECC-Seq^*方法检测得到的插入片段末端7bp序列上的突变谱。

图17是实施例中基于酶打断的文库制备方法结合PECC-Seq*方法检测得到的背景突变谱(去除末端1bp)。

图18是实施例中分布于相同Tiles的重复片段的分布情况。

图19是实施例中分布于不同Tiles的重复片段的分布情况。

图20是实施例中涉及的COSMIC数据库中AAs引起的特征性突变谱COSMICSBS22(https://cancer.sanger.ac.uk/signatures/sbs/sbs22/)。

图21是实施例中Covaris基因剪切仪片段化后DNA片段的长度分布。

图22是实施例中文库片段长度筛选后DNA片段的长度分布。

图23是实施例中制备得到的文库片段的长度分布。

图24是实施例中“PEEA_PE_E”错误的分布。

图25是实施例中“PEEA_E”错误的分布。

图26是实施例中“PECC_PE_E”错误的分布。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例的附图，对本发明实施例的技术方案进行清楚、完整地描述。显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例。基于所描述的本发明的实施例，本领域普通技术人员在无需创造性劳动的前提下所获得的所有其它实施例，都属于本发明保护的范围。

以下实施例采用基于双端测序重叠片段和DNA双链互补片段的超低频突变检测方法，其技术路线如下：

(1)制备不依赖于PCR扩增(PCR-free)的文库。现有技术中利用PCR扩增过程来获得拷贝片段，使用PCR扩增很难控制每一模板DNA所获得的拷贝片段的数目。同时，为了保证大多数模板DNA(即，大多数基因组位点)获得足够的拷贝片段进行分子一致性测序分析，需要增加PCR扩增的循环数。基于PCR扩增过程获得拷贝片段的分子一致性测序方法难以保证最优的分子一致性测序分析拷贝片段的数目(即，两对从DNA互补双链来源的拷贝片段)。此外，依赖于PCR扩增过程的分子一致性测序方法还存在基因组覆盖度不均匀和可能引入额外PCR错误等问题，选择构建PCR-free的文库。

(2)缩短文库片段长度产生双端测序重叠片段结合互补链信息获得DNA模板的拷贝片段进行测序错误校正。对双链DNA模板的两条单链分别进行测序可获得自然来源的两条反向互补的序列信息，其反映相同的模板信息，是双链DNA模板的一对拷贝片段。在使用PCR-free的文库制备策略的前提下，要实现前述的最优设计(DNA互补双链来源的两对拷贝片段)，则意味着对于同一DNA互补链来源的两条单链模板需要分别获得非PCR拷贝来源的一对拷贝片段。在Illumina的双端测序模式中，会对单链的文库片段从两端分别进行测序，产生两条方向相反的读段Read1(R1，方向P5→P7)和Read 2(R2，方向P7→P5)。通常情况下，用于Illumina双端测序的DNA文库平均插入片段长度为～350bp或～550bp。因此，在Illumina常用的2×150bp的双端测序平台上，其产生的R1和R2不会出现重叠。为获得来源于同一单链模板的非PCR拷贝来源的一对拷贝序列进行测序错误校正，设计拟对标准的Illumina文库制备方法进行修改，用于获得较短的测序文库片段(目标平均插入片段长度约为150bp)。截短的文库片段其经双端测序产生的R1和R2会出现重叠；重叠序列反映了同一单链模板相应位置的相同的序列信息，是同一单链模板产生的一对拷贝序列。利用DNA互补双链产生的两条单链模板分别经双端测序产生的重叠片段，获得满足前述最优设计条件的一组同一模板来源的拷贝片段用于测序错误的校正。

(3)优化文库制备条件降低文库制备过程引入的测序错误。结合PEEA-Seq方法对测序过程中的错误来源进行系统分析；发现：基于超声的DNA片段化过程会引起以模板上G的氧化损伤为主的单链损伤，这些损伤经过末端修复后会被固定为DNA双链模板上的碱基变异(背景测序错误)，难以通过上述设计的拷贝片段进行测序错误校正以去除。针对此类测序错误产生的原因，优化文库制备条件，包括结合使用TE缓冲液、缩短片段化时间进行基因组DNA片段化和使用特异性酶进行末端修复来减少文库制备过程引入的大量测序错误，进一步提高突变检测的准确性。

(4)内源性标签“标记”拷贝序列。使用“分子标签”标记拷贝片段是实施分子一致性测序分析的关键步骤。模板“添加”分子标签标记后，同一模板来源的拷贝片段会携带与模板相同的分子标签；不同模板来源的片段携带不同的分子标签。根据分子标签的信息识别从同一模板来源的一组拷贝片段，从而进行后续的比对和测序错误的校正。分子标签包括两类：外源性标签和内源性标签。外源性标签通常是一段整合在接头序列中的人工合成的短序列，其会随着模板片段的PCR扩增过程标记每一条拷贝片段(图1)。内源性标签则是指文库片段的断裂位点信息。当细胞数量较少时，不同细胞单倍体来源的基因组DNA片段其断裂位点往往是不同的；根据文库片段的断裂位点信息(即，其比对到基因组上的位置信息)也可以标记识别同一模板来源的一组拷贝片段。外源性标签需要合成并添加特殊的接头序列，其在文库制备和数据分析中相对复杂。此外，由于分子标签在合成或者测序时可能出现错误，会导致不同模板来源的片段被错误标记为同一模板来源的片段，从而导致分析的偏倚。而使用内源性标签时，主要依靠测序片段在模板上的比对位置对拷贝片段进行识别和标记；能较好地克服前述外源性标签的不足。但内源性标签(基因组上的比对位置)的数目相对外源性标签较为有限。在使用内源性标签时，需要考虑“标签冲突”的情形：随着测序深度(基因组DNA模板)的增加，不同单倍体来源的基因组DNA片段具有相同断裂位点的概率会逐渐增加；由此也会导致不同模板来源的片段被标记为同一模板来源的现象。发生“标签冲突”时可能会引起结果假阴性的情况。针对于哺乳动物基因组(包括：人类基因组和小鼠基因组等)的全基因组突变检测，其应用的测序深度≤40×。在该深度下，内源性标签出现“标签冲突”的概率可以忽略不计，选择内源性标签作为拷贝片段的标记系统简化文库制备、测序及数据分析流程。

(5)数据分析。数据分析的关键是从大量的测序读段中识别和提取从同一模板来源的拷贝片段。选择内源性标签系统标记同一模板来源的拷贝片段，需要借助拷贝片段在基因组上的比对的位置信息进行拷贝片段的提取。根据拷贝片段产生的原理，同一模板来源的两组双端读段其在基因组上的比对位置应当相同。同时，DNA的两条单链在进行测序时其测序的方向相反，互补链产生的两组双端读段其在基因组上的比对方向相反，即比对方向信息分别为“F1R2”和“F2R1”。故根据双端读段在基因组上的比对信息(比对位置和比对方向)即可进行同一模板来源的拷贝片段的提取。提取拷贝片段后对其序列信息进行比对。只有当一致的碱基数目≥4个时，其发生错误的概率才会低于10^-9-10^-7。当4条拷贝片段的同一位点存在不一致的碱基时，认为其中存在测序错误，即将该位点剔除；当且仅当同一位点4个碱基信息均一致时，才保留该位点一致的碱基信息用于后续的突变分析。利用4条互补链来源的拷贝片段进行测序错误的校正即可达到很高的准确性，在测序数据分析过程中不进行额外的生物信息学的测序信息的校正。

按照如图3所示的技术路线构建PECC-Seq^*方法，基于

DNA PCR-Free文库制备流程进行设计，具体是：通过调整标准文库制备流程中的文库制备条件获得插入片段平均长度～150bp的DNA测序文库用于双端测序、测序错误校正和突变分析。插入片段平均长度～150bp的测序文库经过2×150bp双端测序后，其产生的R1和R2片段会出现较大的重叠；这部分重叠序列是模板片段相同位置序列的两个反向互补的拷贝。再结合互补链的一对拷贝信息，得到从一组互补双链DNA模板来源的4条拷贝片段对测序信息进行错误校正。根据片段在基因组上的比对位置信息：同一DNA互补双链来源的两组双端读段其在基因组上的比对位置相同且比对方向相反，从测序数据中提取满足设计的一组包含4条片段的同一模板来源的拷贝片段。对于一组拷贝片段，其同一位点上的4个碱基应当是完全一致的；当其中存在不一致的碱基时，则提示该位点存在测序错误。根据这一思路对拷贝片段信息进行比对，可以对测序错误进行校正、剔除存在测序错误的位点，从而得到可信度高的测序信息。

用PECC-Seq^*方法检测培养细胞TK6人淋巴瘤母细胞全基因组水平的背景突变情况。结果显示，在文库插入片段末端检测得到的突变体频率要远大于文库插入片段内部，提示有末端修复过程导致的“末端修复错误”存在(图14)。根据其分布特征去除末端容易发生测序错误的7bp的碱基信息后，PECC-Seq^*检测得到的背景突变频率为3.88×10^-7；即应用PECC-Seq^*方法可将NGS进行突变检测的错误率降低至＜5×10^-7水平。此外，从约40×全基因组测序数据共获得2.79×10⁸校正后的碱基信息，平均每检测500个拷贝可以获得1个高准确性的碱基信息(Illumina HiSeqXTen测序平台)。

PECC-Seq^*方法已经大幅降低NGS的测序错误率到～10^-7水平，但其仍然是略高于哺乳动物基因组自发突变水平(10^-9～10^-7)。同时，在对背景突变的突变谱分析中提示检测到的背景突变存在大量的异常增加的CG>GC突变(图15)，均提示PECC-Seq^*检测得到的背景突变中仍有较多的测序错误(尤其是CG>GC类型的突变)。此外，在测序效率上，PECC-Seq^*较部分适用于全基因组的分子一致性测序方法(如，BotSeqS)有明显提升，但仍然相对较低。其中，基因组DNA提取、测序文库制备、测序和后续的数据分析过程均可能产生碱基损伤或碱基的变异，表现为检测到的测序错误。为进一步提高方法的准确性，需要对其测序错误产生的原因进行识别，具体如下：

(1)背景测序错误与文库制备过程相关。通过对突变谱的比较分析发现PECC-Seq^*检测得到的背景突变其突变特征(图15)与剔除的末端7bp序列上的变异(即，末端修复错误)特征(图16)接近(余弦相似度＝0.811)。而当利用基于酶打断的文库制备流程(LotusDNA Library Prep Kit，Integrated Device Technology,CA,USA)制备截短的测序文库进行类似的PECC-Seq^*分析时，其得到的背景突变特征(图17)与前述背景突变特征(图15)则明显不同(余弦相似度＝0.164)。结果表明，PECC-Seq^*检测得到的背景突变中存在的大量测序错误主要可能与测序文库制备过程有关，具体来说，可能与

DNAPCR-Free文库制备流程中的片段化过程(即，基于超声的片段化过程)和末端修复过程有关。

(2)结合基于单链的一致性测序方法分析测序错误来源。PECC-Seq^*是一种基于互补链的一致性测序分析改良方法。在分析过程中，对单链上存在的损伤和随机测序错误并不会加以区分，而是根据其4条拷贝片段上存在不一致的碱基信息直接剔除单链损伤和随机测序错误。并且对于单链模板而言，仅有2条拷贝片段(即双端测序的重叠片段)用于随机测序错误的校正；其准确性并不能满足对单链模板上可能存在的损伤进行分析的要求。因而单纯通过PECC-Seq^*分析难以对单链上可能存在的与超声的片段化过程和末端修复过程有关的变异进行分析。需要进一步借助于基于单链的一致性测序分析方法才能较好的识别和表征可能的由于文库制备片段化过程和末端修复过程引入的测序错误(图2A)。

(a)排他性扩增重复片段。在对TK6细胞的测序数据进行拷贝片段提取的过程中发现，数据中存在大量的在基因组上比对位置相同且比对方向亦相同的“重复片段”。由于PECC-Seq^*是基于PCR-free的文库制备策略，因此这些重复片段并不是由于PCR过程产生的“PCR重复片段”。并且，这些重复片段的比例要远高于互补链来源的拷贝片段(比对位置相同且比对方向相反的片段)的比例(～14.87％vs.1.37％)，结合所应用的测序深度较低(～40×)，提示这些重复片段也并不太可能是由于不同单倍体基因组来源的模板片段碰巧具有相同的比对位置导致的。

进一步对这些片段在流动池表面的分布进行探究发现，绝大多数重复片段(89.3％)分布于流动池表面相同的Tiles；并且其倾向于在邻近位置出现，随着重复片段在流动池表面的距离的增加，重复片段出现的概率迅速降低(图18)。同时，对分布于不同Tiles的重复片段(10.7％)分析发现，其同样倾向于在空间位置上相邻的Tiles出现(图19)。简单来说，这些重复片段的分布并不是随机的，其倾向于在流动池表面邻近的位置出现，距离越近，出现的概率越高。结合Illumina平台的相关资料、文献和发现的重复片段的分布特征，提示这些重复片段是由于Illumina平台使用的“排他性扩增(ExclusionAmplification，ExAmp)”技术引入的排他性扩增重复片段(ExAmp重复片段)。Illumina在其HiSeq 3000、HiSeq 4000、HiSeqX、NovaSeq 5000和NovaSeq 6000等测序平台上，引入了模块化的流动池并结合排他性扩增技术用于“簇生成”及碱基识别，以改善流动池的使用效率、减少光学重复及提高测序通量。模块化的流动池表面被预先分割为固定距离排列的纳米孔，“簇生成”过程在每个纳米孔中进行。当开始测序时，变性的单链文库片段进入纳米孔中进行第一链合成。由于使用了排他性扩增技术，当合成的第一链与模板链解链后，会立刻开始后续的“桥式PCR扩增”进行“簇生成”。因为这个过程非常迅速，解链的模板链并不需要被从流动池表面洗脱。模板链是游离在流动池中的，当相邻的纳米孔存在空缺时，模板链会继续占据相邻的纳米孔再次形成新的簇。从其产生的原理可知，数据中观察到的重复片段是同一个单链模板在临近纳米孔中形成的两个独立的测序簇被独立测序产生的，可以看做是一个单链模板来源的两个拷贝。

(b)结合双端测序产生的重叠片段和ExAmp重复片段进行单链的一致性测序分析。当使用4条拷贝片段进行测序错误校正时，即可将随机测序错误导致的变异降低到～10^-10。因此，使用4条单链来源的拷贝即可以满足进行单链一致性测序分析的要求。根据从数据中发现的ExAmp重复片段，再利用数据中由双端测序产生的重叠片段得到满足条件的4条单链模板的拷贝进行单链一致性测序分析。设计类似PECC-Seq^*的单链一致性测序方法(图4)，将其命名为“基于双端测序重叠片段及ExAmp重复片段的一致性测序方法”(Paired-Endand ExAmp Consensus Sequencing，PEEA-Seq)。其文库制备方法与PECC-Seq^*一致；数据分析与PECC-Seq^*类似，唯一的不同是在提取拷贝片段时是根据比对位置相同和比对方向相同的条件进行单链来源的拷贝片段的提取。其从单链模板来源的一致性测序分析，因而能够较好的反映单链损伤情况以及由于末端修复导致的错配。

(c)结合PECC-Seq^*和单链一致性测序方法PEEA-Seq分析背景测序错误来源。双端测序是对同一个模板从两个方向进行测序，理论上双端测序重叠片段中相同位点的碱基信息是完全一致的。当存在不匹配时，主要反映了随机测序错误：包括碱基识别错误、簇生成过程中引入的PCR错误等。PECC-Seq^*背景测序错误主要可能与文库制备中的片段化过程和末端修复过程有关。因此忽略双端测序不匹配的碱基信息的影响，提取双端重叠片段中一致的碱基(后文简称“双端一致性碱基”)信息用于后续分析。双端一致性碱基在PECC-Seq^*和PEEA-Seq中可能出现的错配及不匹配的情况可以包括以下6种情形(表2和3)。

“PECC_E”即PECC-Seq^*检测到的背景突变，其主要特征为以CG碱基对上的变异为主，其中CG>GC的突变频率较高(图15)。在分布上，这些变异容易发生于插入片段末端7bp序列上(图14)。

“PEEA_PE_E”错误指的是仅在ExAmp重复片段中的一组双端序列中出现的错配。两组片段是从同一条单链模板来的，因此错配的信息反映了单链模板在簇生成及其后续测序过程中引入的随机错误，并可能包含导致碱基配对性质发生改变的单链损伤(如：8-oxoG损伤可以同时与正确的碱基C(碱基信息不发生改变)或错误的碱基A(引起G>T变异)发生配对)。分析发现，“PEEA_PE_E”错误最主要的错误类型为T>C和A>G，其次为G>T；这些错误在模板链上的分布均匀，插入片段末端的发生率未明显高于插入片段内部(图24)。

“PEEA_E”错误即单链一致性的测序错误，反映了单链模板上的变异，包括：单链损伤、末端修复过程固定的错配碱基等(图25)。在文库插入片段的内部(剔除末端7bp序列)，最主要的错误类型为G>T。G>T主要反映了超声片段化过程产生的G的氧化损伤，即8-oxoG，导致的错配。但在文库片段末端7bp序列中，除了G>T的变异，还观察到显著升高的C>G的变异。C>G变异在插入片段内部和末端不同的分布情况(插入片段末端C>G的发生率远高于插入片段内部)提示C>G可能与末端修复过程中的单链损伤的固定(即G发生损伤，末端修复过程在其互补链掺入错配的G)有关。

“PECC_PE_E”错误指的是仅在互补链来源的双端序列中的一组中出现的错配。两组片段是从互补链来源的，因此其既包含了单链损伤信息，也包含了单链模板在文库制备过程和测序过程中引入的随机错误(即，叠加了“PEEA_PE_E”错误和“PEEA_E”错误)。和前两种类型错误的结果类似，其主要的错误类型包括：T>C和A>G(“PEEA_PE_E”的主要错误类型)及G>T(“PEEA_E”的主要错误类型)(图26)。

“PEEA_PE_E’”错误和“PECC_PE_E’”错误分别指PEEA-Seq和PECC-Seq^*分析中同一位点的两组双端一致性碱基均出现错配的情况。“PEEA_PE_E’”错误表明可能存在单链损伤导致的不同类型的错配(损伤碱基N发生了N>X和N>Y的错配)；而“PECC_PE_E’”错误则可能提示单链损伤位点在末端修复条件下产生了不一致的错配(配对碱基N:N’发生了N>X和N’>Y的错配)。这两种错误的特征均提示了单链损伤发生于G上，其与A、T或G均可能产生配对。

综合上述对PECC-Seq^*和PEEA-Seq分析中不同类型双端一致性的碱基的测序错误的分析表明，单链模板在簇生成及其后续测序过程中引入的随机错误的特征以T>C和A>G为主，与检测到的背景测序错误的突变特征并不一致；进一步说明背景突变并不是由于随机测序错误引起的。超声片段化过程引起的单链损伤以G损伤为主，产生了大量的G>T的单链损伤，但经过互补链的校正后，最终背景中CG>AT的突变并不是最主要的。“PEEA_E”错误末端显著增加的C>G的变异、“PEEA_PE_E’”错误和“PECC_PE_E’”错误提示可能存在其他未知类型的G损伤，其在末端修复过程中存在与G配对的情形，从而产生前述背景突变中高比例的CG>GC的突变特征。

前述分析结果可以简单概括出PECC-Seq^*检测得到的背景测序错误其产生的原因(图5)。在进行文库制备时，需要对基因组DNA进行片段化。基于超声的片段化过程会引起以G的氧化损伤(以下以“G^*”表示)为主的DNA单链损伤。正常情况下，当G^*发生于DNA互补双链的一条单链上时，利用互补链上正确的C碱基信息可以校正单链损伤产生的影响。然而，当这类损伤发生于双链模板上的单链区域时(包括：末端的单链突出和片段内部的单链缺口)，后续的末端修复过程会补平或填补这些单链区域。由于G^*引起碱基配对性质的改变，末端修复过程会向这些单链损伤位置的互补链上掺入错配碱基。此时，单链损伤被“固定”为模板链上的变异，无法利用互补链信息对由单链损伤引起的测序错误的校正；其最终成为背景中可以观察到的测序错误。

表2：双端一致性的碱基在PECC-Seq^＊和PEEA-Seq分析中可能出现的错配情

表3：双端一致性的碱基在PECC-Seq^＊和PEEA-Seq分析中出现的错配及其分布

*Main souroe of the mismatches

**Number of mismatches in total(in distal region)

方法优化：降低背景测序错误率

背景中测序错误的形成包括两个因素：超声引入单链损伤和末端修复过程固定单链损伤，从以上两个因素进行优化，具体如下：

(1)超声条件对背景测序错误的影响

以CovarisS220设备为例探究不同的超声片段化条件对背景错误率的影响(表4)。利用CovarisS220进行片段化时，其设备参数包括：最大发射功率(Peak incident power)、工作系数(Dutyfactor)、循环数(Cyclesperburst)、水位(Waterlevel)、处理时间及温度等，其中主要可能影响结果的超声参数包括：超声缓冲液成分、片段化时间、片段化的瞬时功率等。

表4

(a)超声缓冲液成分对背景错误率的影响

进行片段化时常用的缓冲液成分包括Tris-EDTA缓冲液(CovarisS220操作手册推荐)和RSB缓冲液(Resuspension Buffer)(

DNA PCR-Free文库制备流程推荐)。固定其他的超声参数，选用不同的缓冲液成分(TE缓冲液或RSB缓冲液)进行基因组DNA的片段化并进行PECC-Seq^*分析。结果显示，在不同的超声条件下，使用TE缓冲液替代RSB缓冲液均可显著降低背景中的CG>AT测序错误(p＜0.05)，背景中CG>AT错误的发生率在TE缓冲液制备的文库中下降了约45-57％(图6)。TE缓冲液的使用能够降低背景中CG>GC错误，但其影响较小，差异并不显著。上述结果表明，使用TE缓冲液代替RSB缓冲液制备片段化的DNA，能显著减少超声过程引入的损伤，尤其是8-oxoG损伤导致的CG>AT测序错误。

(b)超声片段化作用时间对背景错误率的影响

固定超声作用的平均强度，使用不同的超声作用时间对基因组DNA进行片段化。随着超声作用时间的延长，PECC-Seq^*检测到的CG碱基对上的背景测序错误有明显增加，以CG>AT和CG>TA的测序错误的增加为主，CG>GC略有增加；而TA碱基对上的背景测序错误率变化不明显(图7)。结果表明，进行片段化时，应尽可能缩短超声作用时间以减少超声损伤的影响。

(c)超声片段化平均作用强度对背景错误率的影响

固定作用于基因组DNA上的总的超声强度(即，超声作用的平均强度×超声作用时间)，使用不同的超声作用的平均强度(即，瞬时的最大发射功率)和超声作用时间的组合进行基因组DNA的片段化。结果显示，在相同的超声作用总强度下，降低超声作用平均强度、延长超声作用时间制备的测序文库的背景测序错误率要略高于超声作用平均强度高、超声作用时间短的片段化条件下制备的测序文库的背景测序错误率(图8)。提示超声作用时间延长对单链损伤的影响要大于超声作用平均强度升高对单链损伤的影响。

结果表明，使用TE缓冲液能够较RSB缓冲液明显减少由于超声片段化过程引入的背景测序错误。在超声作用条件中，超声作用时间对背景测序错误的影响最大。应在进行片段化过程中尽可能缩短超声作用于基因组DNA的时间。

(2)DNA损伤修复酶对背景测序错误的影响

超声片段化过程向DNA模板中引入了大量单链的氧化损伤，其中以8-oxoG为主。文献表明，在文库制备过程中对受损碱基进行修复，可以减少测序结果中～70％由于8-oxoG引起的CG>AT的测序错误。甲酰胺嘧啶-DNA糖基化酶(Formamidopyrimidine-DNAGlycosylase，Fpg)是一种常用的DNA修复酶，可以修复多种类型的DNA(氧化)损伤，包括8-oxoG等。在末端修复步骤前选择Fpg酶对片段化的DNA进行氧化损伤碱基的修复。结果显示，相同条件下，Fpg预处理可以减少约80％的CG>AT类型的测序错误；而对其他类型的测序错误的发生率影响不大(图9，TK6P1T4vs.TK6_fpg)。

(3)末端修复条件对背景测序错误的影响

(a)高保真T4 DNA聚合酶进行末端修复。文库制备中的末端修复过程指的是利用3’→5’核酸外切酶对3’突出进行切除以及利用DNA聚合酶对5’突出进行补平的过程，该过程将DNA片段的黏末端补齐为平末端，以进行后续的接头连接(图5)。当文库片段内部存在单链缺口时，其也会在这个过程中被DNA聚合酶填补修复(图5)。在实际应用中，通常联用T4DNA聚合酶和大肠杆菌DNA聚合酶I Klenow片段进行末端的补平：两者均同时具有3’→5’核酸外切酶和5’→3’DNA聚合酶活性；其中T4 DNA聚合酶主要发挥3’→5’核酸外切酶作用，用于切除3’突出；而大肠杆菌Klenow片段主要发挥5’→3’DNA聚合酶活性，补平5’突出。T4DNA聚合酶是目前保真性最好的聚合酶之一，其错误率为10^-7-10^-6；而Klenow片段的保真性稍差，错误率为～10^-5水平。因此，当大肠杆菌Klenow片段在末端修复过程中发挥主要的DNA聚合酶活性修复5’突出时，可能更容易引入错配。

在PECC-Seq^*的文库制备过程中，末端修复使用Illumina TruSeq DNA PCR-freeLibrary Prep Kit中配套的末端修复酶的混合液，其实际成分不可知。但根据Illumina其他的文库制备流程，推测其包含了T4 DNA聚合酶和大肠杆菌Klenow片段；因此，在可能主要是Klenow片段进行片段5’突出的补齐以及内部单链缺口的修复。

T4 DNA聚合酶其保真性好，并且由于其同时具有3’→5’核酸外切酶和5’→3’DNA聚合酶活性，故可以单独发挥末端修复的功能。应用T4 DNA聚合酶代替流程中使用的Illumina TruSeq DNA PCR-free Library Prep Kit中配套的末端修复酶的混合液对片段化DNA进行末端修复和文库制备。结果显示，与标准的文库制备方法相比，使用T4 DNA聚合酶进行末端修复能够减少约一半CG碱基对上的测序错误(CG>AT，CG>GC和CG>TA)(图9，TK6P1T4vs.TK6_erp)。

(b)单链特异性核酸酶。DNA聚合酶在末端修复过程中对存在单链损伤的单链片段进行补齐是导致背景中末端修复错误的根本原因。使用高保真T4 DNA聚合酶可以减少补齐过程中的错配。而如果能够剔除DNA聚合酶的补齐过程，理论上能够直接从根本上减少末端修复错误的产生。需要在末端修复过程中使用一种或几种酶分别去除末端3’突出、5’突出以及够切开并去除双链DNA内部的单链缺口，同时所使用的酶需要对双链DNA的降解作用较弱。通过查阅相关资料并结合相关报道，单链特异性的核酸内切酶能够可以特异性地降解DNA或者RNA中的单链位点，实现切开双链DNA内部的单链缺口以及降解末端的3’突出和5’突出的功能。

单链特异性核酸酶主要包括：S1核酸酶、绿豆核酸酶和BAL31三种。BAL31降解单链DNA的核酸内切酶活性较弱。绿豆核酸酶尽管其补平末端的效果较好，但是其需要在双链DNA内部单链缺口扩大到一定大小时才能发挥核酸内切酶的作用。而S1核酸酶在单链缺口较小时即可发挥内切酶的功能。S1核酸酶是比较好的实现本方法文库片段末端修复过程的酶，可以同时切除文库片段3’和5’突出，以及片段内部的单链缺口。

在前期的条件摸索中，仅使用S1核酸酶处理片段化DNA，其随后的接头的连接效率较低，提示可能是S1核酸酶处理后其末端突出未能被充分消化补平。故选择在S1核酸酶处理片段化DNA后再进行标准的末端修复进一步补平文库片段末端。结果显示，增加S1核酸酶处理后，得到的背景测序错误率可以降低到5×10^-8以下。这一错误率较之前未优化的方法降低了约1个数量级，并且比较接近基因组的自发突变水平。

最终选择TE缓冲液、缩短超声作用时间和S1核酸酶进行末端修复的组合对原有的文库制备方法进行优化，进一步减少背景中的测序错误。经过优化后，背景测序错误率可以降低到5×10^-8以下。

方法优化：提高测序效率

在PECC-Seq^*方法中，一致性测序分析的测序效率主要受互补链的回收效率(即互补链来源的拷贝片段数目占全部测序读段的比例)的影响；提取得到的互补链来源的拷贝片段数目越多，其测序效率越高。双链模板在测序时需要被变性为单链，再分别进行测序。通常，测序簇的数量(实际被测序的单链模板的数量)要远低于投入的文库片段的数量。因此，双链模板的两条链并不是总能同时被测序，从而导致了互补链回收效率和测序效率的降低。

理论上，互补链的回收效率主要受投入模板片段的数量与生成的测序簇的数量(实际被测序的单链模板的数量)的比值的影响：比值越低，互补链的理论回收效率越高，测序效率越高。这一比值在各个测序平台和不同的流动池中是不同的。计算部分常用的Illumina高通量测序平台上使用的不同的流动池和不同试剂盒时推荐的上样单链文库片段数量和流动池表面簇数量的比值(表5)。结果表明NovaSeq 6000测序平台其所需的绝对上样文库片段数量与流动池表面簇数量的比值最低，提示NovaSeq 6000测序平台的互补链的回收效率和测序效率会优于前述PECC-Seq^*分析中使用的HiSeqXTen测序平台。

使用NovaSeq 6000测序平台进行PECC-Seq^*分析。与前述HiSeqX Ten平台数据相比，NovaSeq平台的互补链回收效率提高了约4倍，平均对一个模板碱基进行错误校正平均需要检测116个拷贝，利用NovaSeq 6000测序平台较HiSeqX Ten平台能够显著提高PECC-Seq^*的测序效率。

表5：部分Illumina测序平台上样文库量和流动池簇数量的比值

具体方法步骤

①主要的特殊试剂、耗材及设备

·Illumina TruSeq DNA PCR-Free Low Throughput Library Prep Kit(Illumina,CA,USA)；

·Illumina TruSeq DNA Single Indexes Set A(Illumina,CA,USA)；

·QIAGEN Blood&Cell Culture DNA Mini Kit(Qiagen,Germany)或其它同等的DNA提取试剂盒；

·KAPA Library Quantification Kit(Roche,Basel,Switzerland)；

·S1核酸酶(Thermo Fisher,MA,USA)；

·Tris-EDTA缓冲液(TE buffer)(pH 8.0)；

·Agencourt AMPure XP beads(Beckman Coulter,CA,USA)(可替代Illumina文库制备试剂盒中配套的SPB纯化磁珠使用)；

·Nanodrop^TM 8000分光光度计(Thermo Fisher,MA,USA)；

·Covaris S220或其他同等设备；及Covaris microTube(部件号：520045)或其他配套超声打断管(Covaris，MA，USA)。

·Invitrogen Qubit 3.0或其他同等设备及配套Qubit^TM Assay Tubes和Qubit^TMdsDNA HS Assay Kit(Thermo Fisher,MA,USA)；

·磁力架或其他同等设备(Invitrogen DynamagTM-96Side Magnet，ThermoFisher，MA，USA)；

·Agilent Technology 2100Bioanalyzer及配套Agilent High SensitivityDNA Kit(Agilent Technology,CA,USA)；或PerkinElmer LabChip GX Touch 24及配套24DNA Extended Range LabChip for use with GX/GXII 24、HT High SensitivityReagent Kit(PerkinElmer,MA,USA)；或其它同等设备；

·IlluminaNovaSeq 6000测序仪及配套试剂(Illumina,CA,USA)。

②主要使用的数据分析工具

Trimmomatic、SAMtools、R、Integrative Genomics Viewer

③测序文库制备

a)基因组DNA提取

参考QIAGEN Blood&Cell Culture DNA Mini Kit(Qiagen,Germany)试剂盒或其它同等的DNA提取试剂盒的操作手册进行基因组DNA的提取。使用Nanodrop^TM 8000分光光度计鉴定提取的gDNA纯度，其A260/280比值应在1.8～2.0。使用荧光定量方法(InvitrogenQubit 3.0或其他同等设备，及配套的Qubit^TM dsDNA HS Assay Kit)定量gDNA浓度。

b)PECC-Seq测序文库制备

i.基因组DNA片段化

(1)1-2μg gDNA样本稀释于130μL TE buffer(pH 8.0)，并转移至Covaris超声打断管(Covaris microTube，部件号：520045)，280×g离心5s；

(2)在Covaris S220设备上按以下设定对gDNA样本进行片段化：最大发射功率175W，工作系数10％，循环数200，水位12，温度7℃，处理时间140s；(在其他同等设备上进行片段化所设定的条件需根据方法设计进行相应的摸索及调整)

(3)样品管280×g离心5s，样品转移至新的EP管中。

ii.片段化DNA纯化

(4)使用Invitrogen Qubit 3.0(或其他同等设备，配套QubitTMdsDNA HS AssayKit)定量片段化浓度；

(5)取400ng片段化DNA转移至新的200μL PCR管；

(6)向样品管中加入2倍体积的SPB磁珠(Sample purification beads，Illumina文库制备试剂盒配备；需提前30min回复室温)，吹吸混匀；

(7)室温孵育5min；

(8)样品管280×g，1min离心；

(9)将样品管置于磁力架上，直至溶液澄清(～8min)；

(10)保持样品管在磁力架上，弃去上清液；

(11)保持样品管在磁力架上，清洗磁珠两次：

a)向样品管中加入200μL新鲜配制的80％EtOH；

b)室温孵育30s；

c)弃去上清液；

(12)用20μL枪头吸去磁珠表面残留的EtOH；

(13)保持样品管在磁力架上，室温干燥5min；

(14)向样品管中加入21.5μL RSB缓冲液(Resuspension buffer，Illumina文库制备试剂盒配备；需提前30min回复室温)，将样品管从磁力架上取下，吹吸混匀；

(15)室温孵育2min；

(16)样品管280×g，1min离心；

(17)将样品管置于磁力架上，直至溶液澄清(2～5min)；

(18)将20μL上清液转移至新的200μL PCR管中。

iii.片段化DNA末端修复

(19)向样品管中加入以下试剂，吹吸混匀：

3.7μLRSB缓冲液、6μL 5×反应缓冲液(S1核酸酶试剂配备)、0.3μLS1核酸酶溶(100000units/mL)

(20)样品管简单离心后，置于PCR仪上，按以下设定运行“S1程序”：

a)预热盖板：40℃

b)反应体系：30μL

c)程序：

30℃，30min，4℃保持

(21)重复步骤(6)-(13)，纯化S1核酸酶处理后的片段化DNA；

(22)向样品管中加入52.5μL RSB缓冲液，将样品管从磁力架上取下，吹吸混匀；

(23)室温孵育2min；

(24)样品管280×g，1min离心；

(25)将样品管置于磁力架上，直至溶液澄清(2～5min)；

(26)将50μL上清液转移至新的200μL PCR管中；

(27)向样品管中加入以下试剂，吹吸混匀：

10μL RSB缓冲液；40μL ERP2或ERP3(Illumina文库制备试剂盒配备)

(28)样品管简单离心后，置于PCR仪上，按以下设定运行“ERP程序”：

a)预热盖板：100℃

b)反应体系：100μL

c)程序：

30℃，30min，4℃保持

iv.目标文库片段筛选

(29)向样品管中加入90μL SPB磁珠，吹吸混匀；

(30)室温孵育5min；

(31)280×g，1min离心；

(32)将样品管置于磁力架上，直至溶液澄清(～5min)；

(33)将185μL上清液转移至新的200μL PCR管中；

(34)向样品管中加入125μL SPB磁珠，吹吸混匀；

(35)室温孵育5min；

(36)280×g，1min离心；

(37)将样品管置于磁力架上，直至溶液澄清(～5min)；

(38)保持样品管在磁力架上，弃去上清液；

(39)保持样品管在磁力架上，清洗磁珠两次：

a)向样品管中加入200μL新鲜配置的80％EtOH；

b)室温孵育30s；

c)弃去上清液；

(40)用20μL枪头吸去磁珠表面残留的EtOH；

(41)保持样品管在磁力架上，室温干燥5min；

(42)向样品管中加入17.5μL RSB缓冲液，将样品管从磁力架上取下，吹吸混匀；

(43)室温孵育2min；

(44)280×g，1min离心；

(45)将样品管置于磁力架上，直至溶液澄清(2～5min)；

(46)将15μL上清液转移至新的200μL PCR管中。

v.片段化DNA 3’腺苷酸化

(47)向样品管中加入以下试剂，吹吸混匀：

2.5μL RSB缓冲液、12.5μL ATL或ATL2(Illumina文库制备试剂盒配备)

(48)样品管简单离心后，置于PCR仪上，按以下设定运行“ATAIL70程序”：

a)预热盖板：100℃

b)反应体系：30μL

c)程序：

37℃，30min；70℃，5min；4℃保持

vi.接头连接

(49)向样品管中依次加入以下试剂，吹吸混匀：

2.5μLRSB缓冲液

2.5μLLIG2(Illumina文库制备试剂盒配备)

2.5μLDNA接头序列(Illumina TruSeq DNA Single Indexes Set A)

(50)样品管280×g，1min离心后，置于PCR仪上，按以下设定运行“LIG程序”：

a)预热盖板：100℃

b)反应体系：37.5μL

c)程序：

30℃，10min；4℃保持

(51)向样品管中加入5μL STL，吹吸混匀；

(52)280×g，1min离心。

vii.文库片段纯化

(53)向样品管中加入42.5μL SPB磁珠，吹吸混匀；

(54)室温孵育5min；

(55)280×g，1min离心；

(56)将样品管置于磁力架上，直至溶液澄清(2～5min)；

(57)保持样品管在磁力架上，弃去上清液；

(58)保持样品管在磁力架上，两次清洗磁珠：

a)向样品管中加入200μL新鲜配置的80％EtOH；

b)室温孵育30s；

c)弃去上清液；

(59)用20μL枪头吸去磁珠表面残留的EtOH；

(60)保持样品管在磁力架上，室温干燥5min；

(61)向样品管中加入52.5μL RSB缓冲液，将样品管从磁力架上取下，吹吸混匀；

(62)室温孵育2min；

(63)280×g，1min离心；

(64)将样品管置于磁力架上，直至溶液澄清(2～5min)；

(65)将50μL上清液转移至新的200μL PCR管中；

(66)向样品管中加入50μL SPB磁珠，吹吸混匀；

(67)重复步骤(54)-(60)；

(68)向样品管中加入22.5μL RSB缓冲液，将样品管从磁力架上取下，吹吸混匀；

(69)室温孵育2min；

(70)280×g，1min离心；

(71)将样品管置于磁力架上，直至溶液澄清(2～5min)；

(72)将20μL上清液转移至新的200μL PCR管中。(制备好的文库片段，经过文库质量鉴定后，可直接用于测序)。

c)测序文库质量鉴定

在Aglient Technology 2100Bioanalyzer上利用Agilent High SensitivityDNA Kit(或在LabChip GX Touch 24上利用24DNA Extended Range LabChip for usewith GX/GXII 24和HT High Sensitivity Reagent Kit)鉴定测序文库片段大小，并使用KAPA Library Quantification Kit绝对定量文库浓度。

④Illumina测序

在Illumina NovaSeq 6000平台按照相应的操作手册进行2×150bp双端测序。

⑤测序数据分析

a)数据前处理

i.Trimmomatic去接头

使用的去接头的参数为：

ILLUMINACLIP:/path/to/adapters/TruSeq3-PE-2.fa:2:30:10:1:true。

ii.序列比对

利用BWA软件将去接头后的原始读段比对到参考基因组。使用SAMtools软件对比对文件进行压缩及排序获得BAM文件。

iii.比对文件过滤

按以下条件对比对文件进行过滤：1.剔除存在多重比对的读段；2.剔除存在嵌合比对的读段；3.剔除读段5’端含有软切除序列的读段；4.剔除未恰当比对的双端读段；5.剔除比对质量值低的读段；6.去除经上述过滤后仅含单端数据的读段。该步骤目的是去除低质量的比对序列、保留比对质量高的双端读段。

b)同一模板来源的拷贝片段提取

根据双端读段在基因组上的比对位置及比对方向进行同一模板来源的拷贝片段的提取。比对位置指的是互补链来源的一组双端读段其在基因组上两侧的比对坐标应当是相同的。比对方向指的是互补链来源的一组双端读段其R1和R2在基因组上的比对方向是相反的，两者应分别为F1R2和F2R1，即4条读段的FLAG值应为(99,147)和(83,163)的组合。

在Linux终端利用提取拷贝片段的代码得到按序排列的拷贝片段序列信息(test_cpl_dup.txt，每4行为一组拷贝片段信息)及用于一致性碱基提取的文件(test_crg.txt)。

c)一致性碱基信息提取

该步骤首先依次提取各组拷贝片段各位点上的碱基信息；并结合SAMtools软件生成的pileup文件进行过滤：剔除存在超过三个双端一致性碱基不匹配的读段，剔除分布于超过约2倍平均测序深度的位点及去除存在不一致碱基的位点。最终得到包含各位点碱基信息、在片段分布位置及样本该位点其他读段碱基信息的文件，用于后续突变分析。

d)突变分析

对提取的一致性碱基进行过滤，剔除插入片段末端10bp范围内的位点，筛选候选的突变位点。考虑到基因组单核苷酸多态性位点以及累积突变的存在，设定以下条件过滤，减少SNPs和累积突变的影响：1.突变仅在一组拷贝序列中被观察到，而未在同一测序数据的其他读段中观察到；2.突变未在其他组别的测序数据(即，同一背景基因组的测序数据)中观察到。所有候选的突变位点经基因组浏览器(IGV)检查核实。

本发明实现基于NGS技术的高准确性、低成本的哺乳动物全基因组超低频突变检测，以下以动物实验(应用PECC-Seq方法检测小鼠全基因组突变)为例，证明PECC-Seq用于哺乳动物全基因组超低频突变(点突变)检测所能达到的效果，包括准确性、灵敏度和测序效率(测序成本)。

实施例1动物实验

8周龄的雄性gpt delta转基因小鼠(C57BL/6J背景)分别给予溶剂对照(生理盐水)，0.125mg/kg/d、0.25mg/kg/d、0.5mg/kg/d、1mg/kg/d、2mg/kg/d和4mg/kg/d马兜铃酸I(Aristolochic acidI，AAI)(溶解于生理盐水)灌胃处理(对照组n＝5/组,剂量组n＝6/组)，每日给药1次，连续给药28天。期间，小鼠可以自由摄取食物和水。给药结束3天后处死小鼠，取小鼠肝脏和肾脏组织，迅速置于液氮中冷冻后转移至-80℃冰箱长期保存，用于后续的突变检测。

1)PECC-Seq应用于全基因组超低频突变检测的效果

①方法的准确性

理论上，通过检测标准参比DNA(即序列已知、不存在突变的基因组DNA)可以获得PECC-Seq用于突变检测的错误(率)信息。然而真实情况下，完全不存在突变的基因组是很难获得的，直接评价方法的错误率比较困难。PECC-Seq检测基因组得到的变异信息包括了真实的变异和方法产生的测序错误。故，PECC-Seq检测得到的突变频率＝真实的突变频率+方法的错误率；方法的错误率≤PECC-Seq检测得到的突变频率。利用PECC-Seq检测对照动物得到的背景突变频率可以近似估算方法的错误率。

应用PECC-Seq检测对照组gpt delta转基因小鼠肝脏和肾脏组织得到的背景突变频率如表6所示。从小鼠肝脏和肾脏组织得到的背景突变频率范围在(3.36±0.73)×10^-8。故PECC-Seq检测基因组超低频突变的错误率低于此水平，即低于～3.36×10^-8。

表6：PECC-Seq检测得到的小鼠肝脏和肾脏组织背景突变频率

与目前几种代表性的分子一致性测序改良方法相比，PECC-Seq检测得到的背景突变频率较DupSeq和BotSeqS低约一个数量级，与NanoSeq接近(表7)。表明PECC-Seq的准确性要好于前两种方法，与NanoSeq接近。

表7：PECC-Seq与代表性的分子一致性测序改良方法比较

②基因组覆盖度

PECC-Seq是基于PCR-free的文库制备策略进行文库制备。因此，其不会受PCR扩增的影响，理论上较其他所有依赖于PCR扩增所制备的文库有更好的基因组覆盖度。对从各样本得到的一致性碱基的位置信息进行统计。单次检测获得约1×，共3.19×10⁹个一致性碱基，其覆盖了基因组62.5％(1.70×10⁹)的位点。提示PECC-Seq文库对基因组的覆盖比较均匀，有较高的全基因组覆盖度。现有的几种代表性的分子一致性测序改良方法均基于PCR扩增进行文库制备；其中NanoSeq方法由于还使用限制性内切酶进行片段化，其仅能覆盖约30％的基因组位点。结果表明，PECC-Seq可以对全基因组水平的突变进行检测，其有较高的基因组覆盖度。

③检测突变的灵敏度

基于转基因动物模型的基因突变试验是最经典也是较为灵敏的检测化学物诱发突变情况的方法，已被列入经济合作与发展组织(OECD)化学物毒性评价指南中。当转基因动物基因组上转入的报告基因片段发生突变时，其相应的蛋白功能会发生变化(即表型改变)。通过检测报告基因表型的变化推断报告基因的突变情况，并以报告基因的突变情况反映全基因组的突变情况。基因突变试验反映出基因组碱基水平的突变情况，其检测基因组突变情况的准确性和灵敏度很高。常用的转基因动物模型包括：Big Blue大鼠/小鼠,MutaMouse小鼠和gpt delta大鼠/小鼠等。本实验中选择了gpt delta转基因小鼠作为动物模型。gpt delta转基因小鼠其报告基因gpt基因的突变频率较低，因此相对其他转基因动物模型在检测化学物诱发突变上更为灵敏。选择基于gpt delta转基因小鼠模型的gpt基因突变试验作为参照能较好地反映PECC-Seq方法的灵敏度。

本实验给予gpt delta转基因小鼠不同剂量的致突变剂AAI处理，分别用gpt基因突变试验和PECC-Seq方法检测小鼠肾脏组织(AAI作用的靶器官)基因组突变频率的变化，比较两种方法得到的突变频率的变化情况以估计PECC-Seq方法相较于经典方法的灵敏度。

a)应用经典的gpt基因突变试验检测AAI诱导的突变

应用gpt基因突变试验检测得到的AAI处理的gptdelta转基因小鼠肾脏组织的gpt基因突变频率如表8所示。≥0.25mg/kg/d的AAI处理28天可引起小鼠肾脏组织gpt基因突变频率的显著升高(p＜0.05)，具有剂量-效应关系。本实验条件下，AAI处理引起突变频率升高的LOAEL值为0.25mg/kg/d。

表8：gpt基因突变试验检测得到的小鼠肾脏组织gpt基因突变频率

b)应用PECC-Seq方法检测AAI诱导的突变

应用PECC-Seq方法检测得到的AAI处理的小鼠肾脏组织的突变频率如表9和图12所示。结果显示，AAI处理28天可引起小鼠肾脏组织突变频率的显著升高，具有剂量-效应关系。在本实验中，最低剂量的AAI处理即可引起突变的明显升高(p＜0.05)。表明本实验条件下AAI处理引起突变频率升高的LOAEL值为≤0.125mg/kg/d。

表9：PECC-Seq方法检测得到的小鼠肾脏组织突变频率

研究表明，AAI引起的突变以TA>AT为主，并且好发于5’-CpApG-3’三核苷酸序列上。COSMIC数据库将AAs暴露引起的前述突变特征命名为突变标签SBS22(图20)。如图12所示，PECC-Seq检测得到的AAI处理引起的突变以TA>TA为主，并且其突变频率随着AAI剂量增加而明显升高。此外，仅在高剂量下观察到CG>AT和CG>TA突变的增加。对高剂量检测到的突变绘制三核苷酸突变谱(图13)，结果显示其在5’-CpTpG-3’(互补链上为5’-CpApG-3’)三核苷酸序列上有明显的偏向性。高剂量组(1、2和4mg/kg/d剂量组)得到的突变谱与AAs特征性的突变标签COSMIC SBS22非常接近，其余弦相似度分别为0.847、0.942和0.973。

c)PECC-Seq方法与gpt基因突变试验比较

PECC-Seq检测得到的突变频率为碱基发生突变的概率；而gpt基因突变试验的结果为gpt基因片段发生突变的概率。两者概念上并不一致，不能直接将两种方法得到的突变频率结果进行比较。因此，通过比较两种突变分析方法得到的突变频率的变化倍数来初步指征方法检测的灵敏度。与gpt基因突变试验的结果相比，PECC-Seq检测得到的突变频率升高的LOAEL值≤0.125mg/kg/d，低于gpt基因突变试验检测到突变频率变化的LOAEL值0.25mg/kg/d。同时，在低剂量下，相同剂量组PECC-Seq检测得到的突变频率升高倍数要高于gpt基因突变试验。结果表明PECC-Seq方法的灵敏度优于目前经典的gpt基因突变试验。

④测序效率和测序成本

如表7所示，PECC-Seq其平均每检测～132个碱基能够获得1个校正后的碱基信息,。其测序效率要明显优于BotSeqS方法。与NanoSeq相比，其测序效率略低。从测序成本而言，按背景突变频率为～10⁸、每样本需要检测10⁸碱基计算，平均每样本需要检测约13G数据，按￥50/G测序费用计算，每样本的测序成本约为￥650。因此，尽管其测序效率较NanoSeq略低，但其测序成本已能满足日常的应用需求。

⑤小结

PECC-Seq应用于小鼠全基因组突变检测的数据表明，PECC-Seq可以实现对哺乳动物全基因组的突变检测。其检测的错误率≤3.36×10^-8，非常接近基因组的自发突变水平。同时，其在进行基因组自发突变和诱发突变的检测时，其检测的灵敏度与经典的gpt基因突变试验相当。与传统的gpt基因突变试验相比，PECC-Seq还可以直接获得全基因组上的突变特征信息。结果表明，PECC-Seq可以用于哺乳动物全基因组水平的自发突变和诱发突变的检测。

实施例2

以哺乳动物的全基因组为检测对象，并不适合小型基因组或特定基因片段的突变检测，以下以小鼠肝脏组织全基因组突变检测为例进行说明。

(1)肝脏组织gDNA提取。选择新鲜或者冷冻的肝脏组织约15mg-25mg，按照QIAGENBlood&Cell Culture DNA Mini Kit(Qiagen,Germany，Cat.No.13323)或者QIAGEN DNeasyBlood&Tissue Kit(Qiagen,Germany，Cat.No.69506)试剂盒的操作手册进行基因组DNA提取。使用NanodropTM 8000分光光度计(Thermo Fisher，MA，USA)鉴定提取的gDNA纯度，gDNA的A260/280值应在1.8～2.0。

(2)测序文库制备(根据前述“具体方法步骤”部分，进行测序文库制备)

(a)gDNA定量及片段化

使用Invitrogen Qubit 3.0(Thermo Fisher，MA，USA)定量gDNA浓度。将1～3μgDNA样本稀释于130μL TE buffer(pH 8.0)后，转移至Covaris超声打断管(CovarismicroTube，Part Number 520045，Covaris，MA，USA)。在Covaris S220上按以下设定对gDNA样本进行超声打断：最大发射功率175W，工作系数10％，循环数200，水位12，处理时间140s及温度7℃。

经过Covaris S220基因剪切仪片段化后的DNA片段典型的长度分布见图21。片段长度分布于100-500bp之间。由于随着超声时间的延长，超声导致的氧化损伤的程度会增加，有可能引起背景中随机测序错误的增加。因此，在该步骤中并未按照Covaris S220推荐的条件直接将片段打断至目标长度，而是尽可能减少超声打断的时间；故所得到的片段长度要大于目标长度。

在其他同等设备上进行片段化所设定的条件需根据方法设计进行相应的摸索及调整。在获得目标片段的前提下，尽可能缩短超声片段化的时间。

(b)片段化DNA纯化(见“具体方法步骤”部分)

使用Invitrogen Qubit 3.0定量gDNA浓度。取～400ng片段化DNA进行纯化。

(c)片段化DNA末端修复

·S1核酸酶处理

纯化后的DNA片段根据起始的gDNA量确定S1核酸酶的使用量。按照30units S1核酸酶/400ngDNA配制反应体系，即：20μL DNA(～400ng起始gDNA)，6μL 5×ReactionBuffer，0.3μL S1核酸酶(100000units/mL即30units S1核酸酶)，3.7μL ddH₂O配制30μL反应体系；置于PCR热循环仪25℃孵育30min。

·处理后DNA片段纯化(见“具体方法步骤”部分)

·末端修复(见“具体方法步骤”部分)

理论上，S1核酸酶能够特异性降解插入片段末端的单链位点，实现末端补平的功能。但文库制备的结果表明，如果不对S1核酸酶处理后的片段进行末端修复，其接头连接效率很低，实际操作中对于S1核酸酶处理后的片段仍需再进行末端修复。

(d)目标文库片段筛选

末端修复后的文库片段根据前述“具体方法步骤”部分的条件进行目标文库片段的筛选。经过文库片段长度筛选后DNA片段典型的长度分布见附图9。通过上述文库片段筛选条件，可以从超声片段化后获得的较长的片段中去除长片段获得目标片段。筛选后的片段长度分布于100bp-300bp之间；平均长度为～170bp；接近于目标插入片段长度。这些文库片段经过后续的2×150bp的双端测序，能够获得较多的双端测序重叠片段的序列信息。

·片段化DNA 3’腺苷酸化(见“具体方法步骤”部分)

·接头连接(见“具体方法步骤”部分)

·文库片段纯化(见“具体方法步骤”部分)

·测序文库质量控制

在LabChip GX Touch 24(PerkinElmer，MA，USA)利用24DNA Extended RangeLabChip for use with GX/GXII 24(PerkinElmer，MA，USA)和HT HighSensitivityReagent Kit(PerkinElmer，MA，USA)或在鉴定测序文库片段大小，并使用KAPA LibraryQuantification Kit(Roche，Basel，Switzerland)绝对定量文库浓度。

制备的文库片段的长度分布见图23；峰值片段长度在690bp附近。由于制备的文库片段末端为不互补的“Y”型接头，其会影响文库片段在泳道的迁移，从而导致检测的片段长度大于实际长度(58bp 5’TruSeq通用接头序列+插入序列(～150bp)+63bp 3’接头序列(含6bp标签序列)，即～270bp)。因此，仪器检测得到的峰值片段长度较实际长度偏大，但满足方法设计要求的。使用不同仪器设备和不同试剂盒检测得到的峰值片段长度的偏移是不同的，需根据片段筛选的结果或测序数据比对的结果确定。

(3)DNA测序

选择NovaSeq 6000(Illumina，CA，USA)测序平台进行2×150bp的双端测序。平均测序数据量约为10G。根据突变数量可以相应调整测序数据量。当突变频率低时，需要相应增加数据量。

(4)测序数据分析

按照“具体方法步骤”部分分析流程进行数据分析。

(a)对测序得到的数据进行前处理

·Trimmomatic过滤

PECC-Seq使用的文库插入片段平均长度约为150bp。相较于标准的350bp或550bp的插入片段，其在2×150bp的读长下很容易出现“测通”的现象，导致读段末端出现接头污染因此，需对原始的读段使用Trimmomatic软件中的“回文模式”(Palindrome mode)去接头。

·BWA比对后，过滤异常比对读段

PECC-Seq是根据内源性标签(即双端序列比对位置和比对方向)进行同一模板来源的拷贝序列的提取。因此，测序读段的准确比对对后续的数据分析至关重要。由于序列比对到基因组位置的不准确而产生的比对错误是短读长测序平台常见的测序错误来源。同时，这些比对错误也可能会错误地将不同模板来源的片段匹配为同一模板互补链来源的拷贝序列，从而导致数据分析时出现假阴性错误)。因此，需尽可能地减少比对错误，使用比对质量较好及比对准确性较高的读段进行后续的分析，按照PECC-Seq设定条件进行比对读段过滤，83.67％读段可以通过过滤；表明所设定的过滤条件并不会导致测序信息的大量丢失。

(b)提取拷贝片段，进行碱基信息校正

根据互补链来源的片段的特征，利用双端读段在基因组上的比对位置和比对方向信息进行拷贝片段的提取。由于存在ExAmp重复片段，因此数据中存在一定的比例的含有相同比对位置、包含不同比对方向，但双端读段超过3组的情况。为了简化分析，在提取拷贝片段时，仅保留一组读段中刚好含有2对相同比对坐标的双端序列的情况；再从中提取比对方向相反的互补链来源的拷贝片段。此外，在后续的分析中，继续剔除了部分拷贝片段：包括存在≥3个互补链来源的两个双端一致性碱基的不匹配情况的拷贝片段，以及比对到深度超过两倍平均深度的位点的片段；以进一步减少由于可能的重复序列上的比对不准确导致的偏倚。

从拷贝片段中提取一致性碱基时，所有存在不一致碱基的位点均被去除，仅保留4个碱基信息一致的位点。当4个位点信息一致时，其由于随机测序错误及单链损伤导致的背景测序错误的发生率可以降低到10^-9水平，满足超低频突变检测的需要。同时，根据末端修复错误的分布，进一步将分布于末端10bp范围内的序列剔除以减少残留的末端修复错误的干扰。

(c)提取突变信息(过滤SNPs和germline突变位点)

根据方案的设计，检测到的变异是测序错误的概率约为10^-9。因此，检测到的变异的可信度很高，不需要借助额外的生物信息学分析方法对检测到的突变进行校正；其可以直接用于后续的突变分析。由于超低频突变的发生率较低，其VAF很低；当测序深度比较低时，其不太可能在其他的片段中被观察到。因此，设置了一下条件进行过滤，以减少SNPs、germline突变和生命早期发生的突变累积的影响：1.突变仅在一组拷贝序列中被观察到，而未在同一WGS数据的其他读段中观察到；2.突变未在其他组别的WGS数据(即，同一背景基因组的测序数据)中观察到。

结果分析

(1)片段长度分布

将测序数据进行比对后，利用SAMtools统计比对信息。比对片段插入片段的平均长度～180bp。考虑存在异常比对的情况，实际的插入片段平均长度会略低于这个数值。利用比对数据得到的信息与图22所示的检测长度接近。即利用该方法的设计，得到插入片段长度接近150bp的文库片段。

(2)突变分析

如表10所示，应用PECC-Seq对小鼠肝脏进行突变检测。其末端检测得到的变异频率要远高于插入片段内部，高约两个数量级。对末端10bp范围内碱基进行剔除后，得到的背景突变频率在(2.78±0.53)×10^-8接近基因组的自发突变水平。

表10：PECC-Seq检测得到的小鼠肝脏组织背景突变频率

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。

Claims (10)

1.一种基于双端测序重叠片段和DNA双链互补片段的超低频突变检测方法，其特征在于，包括以下步骤：

(1)以不依赖于PCR扩增方法构建DNA文库，同时优化DNA文库制备条件，包括结合使用TE缓冲液、缩短片段化时间进行基因组DNA片段化和使用特异性酶进行末端修复；

(2)缩短文库片段长度产生双端测序重叠片段结合互补链信息，获得DNA模板的拷贝片段并进行测序错误校正，包括：缩短文库插入片段的长度至～150bp，通过Illumina测序的2×150bp双端测序模式获得双端测序重叠片段，用双链DNA模板互补链的一组拷贝片段获得非PCR来源的同一双链DNA模板来源的四条拷贝片段进行测序错误校正；

(3)选择内源性标签对拷贝片段进行标记；

(4)数据处理和突变分析，包括：

根据双端读段在基因组的比对信息进行同一模板来源的拷贝片段提取：互补链来源的双端读段在基因组上的比对位置相同、比对方向相反；

对提取的拷贝片段进行序列信息比对：当同一位点存在不一致碱基时，提示存在测序错误，则该位点被剔除；当且仅当同一位点序列信息一致，保留一致的碱基信息用于突变分析。

2.根据权利要求1所述基于双端测序重叠片段和DNA双链互补片段的超低频突变检测方法，其特征在于，适用于包括人类基因组和小鼠基因组的哺乳动物全基因组，测序深度≤40×。

3.根据权利要求1所述基于双端测序重叠片段和DNA双链互补片段的超低频突变检测方法，其特征在于，步骤(1)中，所述TE缓冲液的组成为10mM Tris-HCl,1mM EDTA-Na₂，pH8.0。

4.根据权利要求1所述基于双端测序重叠片段和DNA双链互补片段的超低频突变检测方法，其特征在于，步骤(1)中，用Covaris Focused-ultrasonicator S220非接触式超声波破碎仪片段化的条件为：最大发射功率为175W，工作系数为10％，循环数为200，水位为12，片段化时间为140s，温度为7℃。

5.根据权利要求1所述基于双端测序重叠片段和DNA双链互补片段的超低频突变检测方法，其特征在于，步骤(1)中，所述特异性酶选自单链特异性核酸酶、DNA损伤修复酶或DNA聚合酶，所述单链特异性核酸酶包括S1核酸酶，所述DNA损伤修复酶包括Fpg酶，所述末端修复酶包括T4 DNA聚合酶。

6.根据权利要求1所述基于双端测序重叠片段和DNA双链互补片段的超低频突变检测方法，其特征在于，步骤(4)中，剔除分布于插入片段末端10bp区域内的碱基信息。

7.根据权利要求1所述基于双端测序重叠片段和DNA双链互补片段的超低频突变检测方法，其特征在于，步骤(4)中，还包括给出提取拷贝片段、一致性碱基和突变分析过程的步骤。

8.根据权利要求1所述基于双端测序重叠片段和DNA双链互补片段的超低频突变检测方法，其特征在于，步骤(4)中，还包括结合基于单链的一致性分析方法对测序过程中的错误来源进行系统分析的步骤。

9.权利要求1-8任一项所述基于双端测序重叠片段和DNA双链互补片段的超低频突变检测方法在哺乳动物全基因组超低频突变检测中的应用。

10.根据权利要求9所述的应用，其特征在于，所述哺乳动物基因组包括人类基因组和小鼠基因组。

Images