CN115485389A - 皮克量dna的全基因组测序方法 - Google Patents

Abstract

本发明涉及单细胞或细胞群的全基因组测序方法，用于识别单细胞或细胞群基因组中的单核苷酸变体、确定单细胞或细胞群基因组中的染色体结构变异或确定单细胞或细胞群基因组中的定相信息。还描述了用于核酸分子测序的索引的DNA文库的制备方法；用于单细胞或细胞群的全基因组测序的索引的DNA文库的制备方法，用于识别单细胞或细胞群的基因组中的单核苷酸变体、确定单细胞或细胞群的基因组中的染色体结构变异或确定单细胞或细胞群的基因组中的定相信息；以及单个细胞或细胞群的全基因组测序的方法，提供用于识别单细胞或细胞群基因组中的单核苷酸变体(SNV)、确定单细胞或细胞群的基因组中的染色体结构变异或确定单细胞或细胞群的基因组中的定相信息的数据。

Description

皮克量DNA的全基因组测序方法

本发明涉及制备用于测序的索引DNA文库的方法，例如用于识别单核苷酸变体(SNV)、确定染色体结构变异或确定单细胞或细胞群(cell-group)的基因组中的定相(phasing)信息的单细胞或细胞群的全基因组测序。

下一代测序彻底改变了我们对人类细胞在健康和疾病中的遗传进化的理解。在大批量癌症基因组测序中，推断变体的流行率，即含有变体的细胞的比例，可以计算肿瘤的克隆组成。反过来，对克隆组成的了解可以构建进化树，以讲述特定肿瘤如何随时间进化的故事(1-3)。分析单个克隆内的共有突变可用于推断可能在肿瘤进化过程中已经起作用的突变过程。理解肿瘤内发生了哪些突变过程以及是什么机制驱动它们是非常可取的，因为这可以为治疗性干预或预测肿瘤的进化轨迹提供机会。然而，测序深度的限制意味着使用标准的大批量全基因组测序(WGS)方法(图1A和图4)只能检测到在肿瘤发生早期出现的非常普遍的突变。因此，对进化事件建模的能力仅限于在肿瘤进化过程中已经确定的早期事件，而不是最近或当前的过程(1，4)。这限制了理解突变过程的实际应用。研究当前或最近的进化事件需要确信的识别出流行率非常低的突变(1，4)。

对单细胞或小群体的空间相关细胞的测序提供了通过检测细胞特异性或克隆特异性突变来解决这个问题的希望(图4)。这给出了当前观察到的细胞中发生的突变过程的读数(图1A)。然而，可以从单细胞或空间相关细胞获得的少(皮克)量DNA的精确测序是非常具有挑战性的。当处理少量DNA时，通过氧化或自发的脱氨基而不可避免的DNA损伤尤其麻烦(5)。这些损伤源分别导致不成比例的人工C>A和C>T突变数量(5-7)。在变体识别(variant calling)期间将这些人工突变识别为变体导致大量假阳性(FP)变体检出(call)。因此，全基因组扩增会在单核苷酸变体(SNV)识别中引发严重错误，从而妨碍准确估计突变负荷(6)。一个重要的问题是，此类突变也可归因于生物学过程，例如随着年龄的增长氧化DNA损伤的积累或APOBEC脱氨酶家族成员的过度活跃(8-10)。

以前无法区分当处理皮克量DNA时生物学驱动的C>A和C>T突变和在文库制备过程中出现的人工突变。此外，测序前通常的全基因组扩增(WGA)步骤会增加人工突变的数量并放大由DNA损伤引起的误差(5)。已经提出了几种方法来减少文库制备过程中的DNA损伤或过滤掉在分析过程中的假阳性结果(5,11-13)。然而，迄今为止，此类技术仍会导致数千个假阳性突变的保留，因此，在得出确切的生物学结论之前需要进行广泛的验证(5,11,12)。由于在大多数情况下无法进行广泛验证(5)，因此需要一种可靠的方法来消除全基因组扩增测序数据中的假阳性变体。

Complete Genomics公司先前发表了一种用于10到20个人类细胞的全基因组测序和的单倍型分析的长片段读取(LFR)方法(Peters BA,et al.Accurate whole-genomesequencing and haplotyping from 10to 20human cells.Nature.2012Jul 11；487(7406):190-5.doi:10.1038/nature11236.PubMed PMID:22785314；PubMed CentralPMCID:PMC3397394)。然而，这种方法非常复杂，容易因索引交叉污染而产生偏差，并产生大量假阳性。

因此，本发明的目的是提供一种改进的方法以制备用于测序、SNV分析、确定染色体结构变异或确定定相信息的DNA文库。

根据本发明的第一方面，提供了一种单细胞或细胞群的全基因组测序方法，以识别单细胞或细胞群的基因组中的单核苷酸变体(SNV)、确定单细胞或细胞群的基因组中的染色体结构变异或确定单细胞或细胞群的基因组中的定相信息，该方法包括：

i)提供包括多行和多列反应孔的多孔阵列板；

ii)提供单细胞或细胞群的基因组DNA，其中基因组DNA分布在多孔阵列板上的多个反应孔中，使得每个反应孔的任何给定位点的单链基因组DNA分子不超过一个，

iii)对每个基因组DNA分子进行全基因组扩增(WGA)，以在每个反应孔中提供基因组DNA分子的多个拷贝(copy)；

iv)将每个反应孔的DNA分子片段化并在每一端连接一对环状接头或使用转座酶递送接头进行标记以形成适配的DNA片段，其中环状接头或转座酶递送接头包括列索引(Column Index，Ci)序列或行索引(Row Index，Ri)序列，其中所述Ci序列对于所述多孔阵列板的一列中的每个反应孔的每个环状接头或转座酶递送接头是通用的，或者其中每个Ri序列对于所述多孔阵列板的一排中的每个反应孔的每个环状接头或转座酶递送接头是通用的；

vi)通过对适配的DNA片段进行索引(indexing)PCR来提供索引DNA文库，其中使用正向和反向索引引物扩增适配的DNA片段以形成索引的PCR产物，其中行索引(Ri)序列或列索引(Ci)序列由每个正向和反向索引引物引入到适配的DNA片段的每一端，使得所得的索引PCR产物包括一对对于一列的每个孔都是通用的侧翼列索引(Ci)序列和一对对于一行的每个孔都是通用的侧翼行索引(Ri)序列两者；以及

vii)对索引的DNA文库进行测序以提供用于确定单细胞或细胞群的基因组中任何单核苷酸变体、确定单细胞或细胞群的基因组中染色体结构变异或者确定单细胞或细胞群的基因组中的定相信息的数据。

有利地，本发明提供了用于单个DNA分子测序方法的索引DNA文库，以从从临床样本获得的皮克量的DNA中获得高质量和数据丰富的测序结果(称为DigiPico；用于皮克DNA的数字测序)。本发明还提供了一种有利的索引策略，以基本上(virtually)消除交叉污染并提高连接效率。首先一组索引被引入通用接头的茎环中。在连接步骤中，板的每列中的所有孔都将接收不同的索引环状接头或转座酶递送接头，因此总共24种不同的寡核苷酸足以用第一组索引对板的所有列进行索引(列索引)。在连接步骤之后，每行中的所有孔都可以合并到一个单独的管中，从而形成16个不同的合并物(pool)。这16个不同的池可以方便地纯化，以用于下一步的索引。在下一步中，每个池中的纯化产物可以仅使用16种不同的索引引物(行索引)通过PCR反应进行索引。这使得单细胞测序具有前所未有的准确性，这是该技术相对于已知方法的重大改进。本发明可用于鉴定单细胞或极少数细胞中的私有突变和潜在新抗原(neo-antigen)，这些突变或新抗原可作为治疗靶点。本发明还可用于确定染色体结构变异，例如数目或结构畸变，或用于确定定相信息。本文中的图18清楚地表明，本发明的方法可以大大提高确定真实核苷酸变体的准确性，相对于Complete Genomics公司的LFR方法去除了许多假阳性，并且更准确地区分呈现不同数量突变的样本。

细胞和细胞群

细胞或细胞群可以包括真核细胞，例如哺乳动物细胞。在一个实施方式中，该细胞是人类的。在一个实施方式中，该细胞至少是二倍体细胞。该细胞可以是癌细胞或癌前(pre-cancerous)细胞。该细胞可以包括肿瘤小岛。在一个实施方式中，该细胞可以来源于受试者的组织活检。

诸如肿瘤小岛的细胞可以是激光捕获的显微切割细胞。

在对来自多个细胞的DNA进行测序的情况下，这些细胞可能是空间相关的细胞。细胞可以共同位于肿瘤或肿瘤区域中。在另一实施方式中，这些细胞可以是或者可以不是近邻。

待确定的SNV可以包含单核苷酸突变。该方法可用于确定基因组DNA中的多个不同SNV。

提供核酸分子和孔分布

核酸可以是纯化的核酸或部分纯化的核酸。在另一实施方式中，可以在细胞裂解物中提供核酸。基因组DNA可以作为纯化的DNA提供。在另一实施方式中，基因组DNA可以由重悬的细胞核或全细胞提供，例如激光捕获的显微切割细胞。

基因组DNA可以包括单细胞或一群细胞(细胞群)的DNA，例如空间相关细胞。基因组DNA可包含约1至30个细胞的DNA。基因组DNA可包含约1至100个细胞的DNA。在另一实施方式中，基因组DNA可包含约1至80个细胞的DNA。在另一实施方式中，基因组DNA可包含约1至50个细胞的DNA。在另一实施方式中，基因组DNA可包含约1至40个细胞的DNA。在另一实施方式中，基因组DNA可包含约10至30个细胞的DNA。在另一实施方式中，基因组DNA可包含约20至30个细胞的DNA。在另一实施方式中，基因组DNA可包含约20至40个细胞的DNA。在另一实施方式中，基因组DNA可包含约10至40个细胞的DNA。

在核酸，例如DNA，是双链的情况下，核酸可以在分配到孔中之前进行变性。变性可以通过加热和/或变性缓冲液来实现。在一个实施方式中，可以使用变性缓冲液，例如来自Repli-g单细胞试剂盒(Qiagen)的D2缓冲液使核酸，例如基因组DNA或含有基因组DNA的细胞核或细胞变性。

核酸，例如DNA，可以分布到孔中，使得每个反应孔中任何给定位点的单链基因组DNA分子不超过一个。核酸的分布可以通过核酸的稀释来促进。因此，在一个实施方式中，可以稀释核酸溶液。本领域技术人员将容易地确定对实现每个反应孔中任何给定位点的单链基因组DNA分子不超过一个所必需的稀释水平和溶液体积。本领域技术人员将认识到，所必需的稀释水平可以用数学方法确定，使得每个反应孔中任何给定位点的单链基因组DNA分子不超过一个在统计学上具有高概率。例如，当细胞数量已知时，Poisson分布可用于该计算。

在一个实施方式中，单细胞的DNA内含物可以分布在单行或单列的孔中。因此，多孔阵列板可用于分析多种不同的单细胞，例如每行或每列一种。至少一个孔可用于添加细胞和提取DNA内含物。在另一实施方式中，细胞或细胞群的DNA内含物分布在单个多孔阵列板的行和列的孔中。

本领域的技术人员将意识到任何标准的多孔阵列板都可以用于本发明的方法中。优选地，多孔阵列板与可以使用的任何PCR和/或测序仪器兼容。多孔阵列板可以包括384孔板，例如24x16孔板。在另一实施方式中，多孔阵列板可以包括1536孔板。本领域的技术人员将理解的是可能需要更多数量的Ri和/或Ci序列以索引更大的阵列板。

使用384孔多孔阵列板可以有利地提供足够的孔来分布约20-30个细胞的稀释基因组DNA链，从而可以为孔提供单个DNA分子。

扩增

在核酸是基因组DNA的实施方式中，基因组DNA分子的扩增可以包括全基因组扩增(Whole Genome Amplification，WGA)。WGA可以包括将用于DNA扩增的扩增试剂添加到基因组DNA的步骤。用于DNA扩增的扩增试剂也可以称为“扩增混合物”。本领域的技术人员将理解的是扩增混合物可包含扩增DNA(即产生DNA的多个拷贝)所需的所有试剂。这样的组分可以包括反应缓冲液、聚合酶和dNTP。DNA聚合报告分子，例如DNA结合染料(例如Evagreen^TM)可以在扩增混合物中提供，例如以允许使用实时PCR监测扩增反应。DNA结合染料可由通过柔性间隔物连接的两种单体DNA结合染料构成。在没有DNA的情况下，二聚体染料可以呈现在DNA结合中是无活性的环状构象。当DNA可用时，环状构象可以通过平衡转变为能够与DNA结合以发射荧光的随机构象。

可以在向孔中添加DNA之前或之后在每个孔中提供扩增试剂。

本领域的技术人员将能够为发生扩增反应提供合适的条件，包括合适的温度和培养时间。例如，板可以在约30℃培养至少约1小时，然后进行热灭活，例如在约65℃持续至少5分钟。

环状接头或转座酶递送接头的片段化和连接

在一个实施方式中，提供了环状接头，使得该方法包括将每个反应孔的DNA分子片段化的步骤和后续的连接反应，以将环状接头连接到片段化的DNA上。片段化的DNA可以在连接前进行末端修复。在一替代实施方式中，可以提供转座酶递送接头，使得该方法包括通过标签化(tagmentation)过程将DNA分子片段化。标签化可以包括提供携带寡核苷酸的转座酶，例如Tn5，其在本文中被称为转座酶递送接头。本领域技术人员将熟悉用于进行标签化以形成适配的DNA分子的常规技术和试剂。

将每个反应孔的DNA分子片段化成多个dsDNA片段可以包括直接片段化，例如酶促或机械片段化。在一个实施方式中，DNA的片段化包括酶促片段化。

DNA的片段化或标签化可以通过向每个孔中的DNA添加片段化或标签化试剂来提供。例如通过使用多孔分配器，例如I-DOT(Dispendix，德国)分配器或类似物，同时将片段化或标签化试剂添加到每个孔中。可以对片段化或标签化反应计时，以提供所需大小的片段。本领域的技术人员将理解的是片段化或标签化反应的计时可以取决于所使用的方法，例如为反应提供的酶的类型和水平。因此，本领域技术人员可以遵循给定反应的标准协议计时，例如反应试剂盒的计时。

片段化试剂可以包括限制酶或切口酶，例如DNase I(脱氧核糖核酸酶I)。在提供切口酶的情况下，可以提供识别切口位点然后切割第二条链的单链特异性酶。在一个实施方式中，可以使用文库制备试剂盒，例如Lotus DNA文库制备试剂盒(IDT，美国)。

在DNA片段化以形成dsDNA片段后，可以对dsDNA片段进行末端修复和dA加尾，以便它们可以连接到其他DNA分子，例如环状接头。用于末端修复和/或dA加尾的酶可以包括DNA聚合酶，例如T4 DNA聚合酶，和多核苷酸激酶(PNK)，例如T4多核苷酸激酶。T4 DNA聚合酶(在dNTP存在的情况下)可以填补5'突出端，并将3'突出端向下修剪到dsDNA接口以生成平端。然后T4 PNK可以磷酸化5'末端核苷酸。可以提供具有末端转移酶活性并留下3'末端腺嘌呤的DNA聚合酶，例如Taq DNA聚合酶，用于A加尾。

在一个实施方式中，dsDNA的片段化、末端修复和dA加尾都在单个反应中进行。

在一个实施方式中，环状接头可以通过连接引入到片段化的DNA上。环状接头与dsDNA片段的连接可以包括添加环状接头和连接酶，例如T4 DNA连接酶。

环状接头可包含具有二级茎环结构的寡核苷酸，例如DNA。茎环结构可以由单个寡核苷酸分子提供，该寡核苷酸分子包含一对在环区侧翼的互补序列区，其中该对互补序列排列成相互杂交以形成环状接头的茎环结构。环状接头进一步编码茎区中的列索引(Ci)序列或行索引(Ri)序列。

列索引(Ci)序列或行索引(Ri)序列可以包含能够分别将DNA标记为来自行或列的预定序列。列索引(Ci)序列或行索引(Ri)序列的长度可以是至少三个核苷酸。

在一个实施方式中，适配的DNA片段的末端可以是对称的。具体而言，连接到dsDNA片段的每一端的环状接头或转座酶递送接头是相同的，使得每个dsDNA片段接收一对相同的侧翼的环状接头或转座酶递送接头。同一适配的DNA片段上的Ci序列对可以是相同的。替代地，如果提供了Ri序列，则同一适配的DNA片段上的一对Ri序列可能是相同的。

有利地，在适配的DNA片段上提供两个相同的Ci或Ri序列提供了标记物，以避免分析分别由不同列或不同行之间交叉污染而导致的索引DNA文库序列。具体而言，可以从数据分析中丢弃在每一端都没有匹配Ci序列的任何索引DNA文库序列。替代地，如果提供了Ri序列，则可以从分析中丢弃在每一端都没有匹配Ri序列的任何索引DNA文库序列。这可以提供第一级冗余以消除后续数据中的索引交叉污染，该索引交叉污染在索引库的制备及其后续分析中是一个重要问题。

环状接头可提供3'或5'突出端(overhang)以帮助连接到dsDNA片段。当环状接头的茎区杂交在一起时(即环状接头处于二级/茎环结构)，可以提供3'或5'突出端。3'或5'突出端可以对应于dsDNA片段上的互补突出端，这些片段已被末端修复并制备用于连接。突出端可包含单个胸腺嘧啶。

环状接头序列可以包含SEQ ID NO：1的序列或其功能变体。

连接环状接头后，可以切割环状DNA的单链区域。环状DNA的单链区域可以被酶促切割，例如通过USER(尿嘧啶特异性切除试剂(Uracil-Specific Excision Reagent))酶，其在环中存在的尿嘧啶位置处产生单核苷酸缺口。因此，在一个实施方式中，环状接头可以在环状区域中包含尿嘧啶。

合并一行孔

如果在适配的DNA片段中提供Ci序列，则该方法可以另外包括在索引PCR之前将每个反应孔的适配的DNA片段合并成一行的步骤。替代地，如果在适配的DNA片段中提供Ri序列，则该方法可以另外包括在索引PCR之前将每个反应孔的适配的DNA片段合并成一列的步骤。然后，合并的适配的DNA片段可用于在单个合并(pooled)反应中对每行或每列的PCR索引，其具体取决于合并的内容。在一替代实施方式中，在进行索引PCR之前，不会合并列或行。

有利地，在索引PCR之前合并行或列大大提高了文库制备的效率。例如，对于16x24(384)孔板，如果在引入环状接头或转座酶递送接头和索引PCR步骤之间合并16行，则只需要16个单独的索引PCR反应，而不是384个索引PCR反应(如果他们没有合并)。

尺寸选择和索引PCR

在索引PCR之前，可以根据尺寸选择适配的DNA片段，例如当可适用时，也可以从反应中去除自连接接头。所需尺寸的一个示例可以是长度大约为300-400bp。尺寸选择可以通过分离或纯化所需长度的适配的DNA片段来提供，例如使用凝胶或珠。SPRI珠(固相可逆固定珠(Solid Phase Reversible Immobilisation bead))可用于尺寸选择。SPRI珠可以包括涂覆有羧基基团(以琥珀酸的形式)的磁性颗粒，其可以非特异性和可逆地结合DNA。

索引PCR可以包括将适配的DNA片段与一组正向和反向索引PCR引物和PCR试剂混合的步骤。正向和反向索引PCR引物可以包含被布置成与适配的DNA片段序列杂交以引发聚合的序列。布置成与适配的DNA片段序列杂交以引发聚合的序列可以是互补序列。用于从正向和反向索引PCR引物引发聚合的序列可以由环状接头或转座酶递送接头提供。用于从正向和反向索引PCR引物引发聚合的序列可以位于适配的DNA片段的Ci或Ri序列的侧翼，使得Ci或Ri序列结合到索引的PCR产物中。

由正向和反向引物提供的用于杂交的互补序列各自的长度可以在大约15和30个核苷酸之间，例如长度为大约26个核苷酸。

在适配的DNA片段包含Ci序列的实施方式中，正向和反向索引PCR引物可以各自包括Ri序列，用于在索引的PCR产物中提供一对Ri序列。在合并行的情况下，将Ri序列添加到池中的每个适配的DNA片段中(来自一行的所有孔)。替代地，在未合并行的情况下，可以为一行中的每个孔提供相同的Ri序列。

在一替代实施方式中，其中适配的DNA片段包含Ri序列，正向和反向索引PCR引物可以各自包括Ci序列，用于在索引的PCR产物中提供一对Ci序列。在合并列的情况下，将Ci序列添加到池中的每个适配的DNA片段中(来自一行的所有孔)。替代地，在未合并行的情况下，可以为一列中的每个孔提供相同的Ci序列。

可以由正向和反向引物提供的行索引(Ri)序列对于一行的或来自一行的每个适配的DNA片段可以是相同的。替代地，可以由正向和反向引物提供的列索引(Ci)序列对于一列的或来自一列的每个适配的DNA片段可以是相同的。

由正向和反向引物提供的行索引(Ri)序列或列索引(Ci)序列各自的长度可以至少为3个核苷酸，例如长度为大约8个核苷酸。

索引PCR产物的所得末端可以是对称的。例如，原始DNA片段序列侧翼的序列可以是对称的。索引PCR产物可以包含侧接一对相同的Ci序列(即内侧翼)和还侧接一对相同的Ri序列(即外侧翼)的DNA片段序列。在一替代实施方式中，索引PCR产物可以包含侧接一对相同的Ri序列(即内侧翼)和还侧接一对相同的Ci序列(即外侧翼)的DNA片段序列。

有利地，除了先前分别提供的Ri或Ci序列(由环状接头或转座酶递送接头提供的)之外，在索引DNA片段上提供两个相同的Ci或Ri序列对也提供了标记物，以避免分析分别由不同列或不同行之间交叉污染而导致的索引DNA文库序列。具体而言，可以从数据分析中丢弃在每一端都没有匹配Ci序列的任何索引DNA文库序列。替代地，如果提供了Ri序列，则可以从分析中丢弃在每一端都没有匹配Ri序列的任何索引DNA文库序列。在索引DNA片段上提供匹配的Ci和Ri序列对两者可以提供第一级和第二级冗余，以消除后续数据中的索引交叉污染，索引交叉污染在索引库的制备及其后续分析中是一个重要问题。

正向和反向索引PCR引物还可以包含测序接头序列，使得测序接头结合到索引PCR产物中。引物上的测序接头序列可以是5'。

测序接头可以是索引的PCR产物的末端。在提供测序接头序列的情况下，索引的PCR产物所得的末端可能不是对称的。例如，可以用不同于另一端的测序接头的测序接头来调整索引的PCR产物的一端。本领域的技术人员将了解给定测序技术需要的测序接头。例如，在染料测序(例如Illumina染料测序)的情况下，测序接头可以是P5和P7测序接头(即P5在索引的PCR产物的一端，而P7在另一端)。提供P5序列的索引引物可以包含SEQ ID NO：2的序列。提供P7序列的索引引物可以包含SEQ ID NO：3的序列。

一旦形成，则索引的PCR产物可以称为“索引的DNA文库序列”或“索引的DNA片段”。合并的索引的PCR产物、索引的DNA序列或索引的DNA片段可以称为“索引的DNA文库”。

索引的DNA文库

可以过滤索引的DNA文库的索引的DNA片段大小，使得只有具有所需或合适长度的索引的DNA片段可用于测序。在索引PCR之后，可以纯化/分离索引的PCR片段，例如通过珠(例如SPRI珠)。纯化可以去除不需要的短片段、引物二聚体或其他PCR人工制品或试剂。

可以检查索引的文库的尺寸分布是否合适。我们通常检查文库尺寸分布，如在Tapestation或生物分析仪仪器(Agilent)或类似仪器上检查。索引的DNA文库的大小可以在准备测序后通过稀释来调整，例如调整到约4nM。

可以保存索引的DNA文库以供以后使用，例如测序。例如，DNA文库可以冷冻或冷藏保存。

索引的DNA文库的测序

索引的DNA文库可以进行测序，或被调整以进行测序。测序可以是下一代测序(NGS)。测序可以是染料测序(例如Illumina染料测序)、纳米孔测序或离子激流测序。本领域技术人员将熟悉可使用的多种不同的测序技术/方法，以及所需的测序接头。

测序可以是多重测序，其中同时对多个索引的DNA文库进行测序。

突变/核苷酸变异的测定和数据分析

该方法可以包括通过确定源自单孔的基本上所有索引的DNA文库序列是否包含相同的SNV，或者是否只有一部分索引的DNA文库序列包含相同的SNV，来确定单细胞或细胞群的基因组中的任何真实SNV。在源自单孔的基本上所有索引的DNA文库序列中显示的SNV可被确定为基因组DNA中的真实SNV。另外或替代地，仅在源自单孔的索引的DNA文库序列的一部分中发现的SNV可被确定为假阳性(FP)SNV。假阳性SNV可能是损伤引起的错误或复制错误。

该方法可以进一步包括将源自代表基因组DNA的一条链的单孔的索引的DNA文库序列与源自代表基因组DNA的互补链的另一孔的索引的DNA文库序列配对。基本上在基因组DNA的两条互补链的所有索引的DNA文库序列中存在的SNV，可能被确定为真实SNV。基本上在基因组DNA的两条互补链的所有索引的DNA文库序列中不存在的SNV，可能被确定为假阳性(即不是真实SNV)。

确定是否源自单个孔的基本上所有索引的DNA文库序列包含相同的SNV或者是否只有一部分索引的DNA文库序列包含相同的SNV的步骤可以在计算机中(in silico)进行，例如使用BAM文件数据。另外或替代地，将源自代表基因组DNA的一条链的单个孔的索引的DNA文库序列与源自代表基因组DNA的互补链的另一孔的索引的DNA文库序列配对的步骤在生物计算机中进行，例如使用BAM文件数据。

在一个实施方式中，来自肿瘤细胞、疑似肿瘤细胞或癌前细胞的测序数据可以与从取自正常组织(即非癌组织)的正常细胞(即非癌细胞)获得的测序数据(例如作为对照)进行比较。因此，在一个实施方式中，该方法包括从肿瘤细胞、疑似肿瘤细胞或癌前细胞以及正常(即非癌)细胞制备索引的DNA文库。可以为每种细胞类型并行制备索引的DNA文库，例如在同一多孔板的不同孔中，或分别制备索引的DNA文库。来自不同类型细胞的索引的DNA文库的测序可以在同一测序运行中运行。不同类型的细胞(例如癌细胞或正常细胞)可以来自同一受试者。

可以在计算机中计算特定核苷酸变体为真实SNV或假阳性的概率分数，从而确定给定的变体核苷酸具有成为真实SNV或假阳性的统计学显著概率。

在一个实施方式中，对DNA文库进行测序以确定文库中的SNV包括从多个孔生成多重测序数据并分析SNV的数据。

在一个实施方式中，分析SNV的数据包括解复用(de-multiplexing)测序数据，使得来自每个孔的数据被分配给单个孔组。此外，可以在同一测序运行中对单独的索引的DNA文库进行测序，因此，该方法可以进一步包括对测序数据进行解复用，从而将不同的索引的DNA文库进行识别/分组。

提供的序列数据可以是成对读取的FastQ文件的形式。可以为了去除接头序列而修剪序列数据，例如成对读取的(Paired-read)FastQ文件。也可以为了质量而修剪序列数据，例如成对读取的FastQ文件。本领域的技术人员将能够容易地调整序列中每个碱基读数的质量分数的所需阈值，例如使用诸如TrimGalore的程序。得到的数据可以称为“修剪数据(trimmed data)”。

在一个实施方式中，分析SNV的数据包括将序列数据映射到参考基因组，例如人类hg19参考基因组，以生成序列比对图(SAM)或其二进制文件版本(例如BAM文件)形式的比对的测序数据。可以使用修剪后的读取数据映射到参考基因组。SAM或BAM文件数据可用于确定存在于每个孔的SNV。

映射可以使用诸如Bowtie2之类的程序，其中激活了ignore-quals参数并标记了重复读数，例如使用Picard工具。

可以对所有单个BAM文件以及来自所有孔的合并BAM文件执行联合变体识别，例如使用变体识别程序，例如Platypus变体识别程序(caller)。

可以从数据中过滤掉低质量(即低置信度)的变体。例如，可以通过应用质量过滤器从数据中删除低质量(即低置信度)的变体。Platypus识别程序中的示例质量过滤器可包括QUAL>60、FR>0.1、HP≤4、QD>10和SbPval≤0.95。本领域技术人员将认识到过滤掉低置信度变体是常规程序，并且每个变体识别程序可以根据它们的算法对每个变体具有不同的置信度分数，并且置信度分数可以用来过滤低置信度(质量)的变体。因此，具体参数取决于所使用的变体识别程序。

可以确定覆盖每个位点的孔的总数(Tw)和支持每个变体的孔的数量(Vw)。孔计数过滤器，例如Tw>5、Vw>2和Vw/Tw>0.1，可用于仅保留高置信度位点进行分析。

可以从分析中删除具有不良可映射性的基因组区域(即更可能对偏读数的已知区域)，例如使用VCFtools。

然后可以使用已在数据中识别的高置信度变体的结果列表对WGS数据(例如来自血液和大批量肿瘤)执行变体重识别(基因分型)，例如使用Platypus。PlatypusminPosterior参数可以设置为0，minMapQual参数可以设置为5。可将在标准WGS数据中均完全不受支持的任何变体提取为UTD(DigiPico独有的)变体。可将在血液样本的大批量测序数据中也完全存在的任何变体(基于GATK分析)提取为TP(真阳性)变体。

使用人工神经网络(ANN)

可根据本文所述的方法和计算进行索引的DNA序列的计算机确定或配对，和/或概率分数的计算。在一个实施方式中，可以通过人工神经网络(ANN)模型，例如通过多层感知器，执行索引的DNA序列的计算机确定或配对，和/或概率分数的计算。

多层感知器可能有由N个神经元组成的输入层(例如N＝41)，其中N是每个实验中使用的特征数。ANN模型可以包括至少两个具有ReLU(修正线性单元(Rectified LinearUnit))激活的隐藏层。ANN的最后一层可能是具有sigmoid激活的单个输出神经元。损失函数可以是二元交叉熵。

例如，可以使用Keras在Python3中对ANN进行编程。本领域技术人员将认识到Keras是用于开发和评估深度学习模型的免费开源Python库。然而，可以使用其他库。

可以使用一个或多个数据集对ANN进行预训练。例如，ANN可以使用包含已知核苷酸变体的数据集进行训练。

其他方面

根据本发明的另一个方面，提供了一种制备用于核酸分子测序的索引的DNA文库的方法，该方法包括：

i)提供包括多行和多列反应孔的多孔阵列板；

ii)提供核酸分子，其中核酸分子分布在多孔阵列板上的多个反应孔中，使得每个反应孔的任何给定位点的单链核酸分子不超过一个，

iii)对核酸分子进行扩增，以在每个反应孔中提供核酸分子的多个DNA拷贝；

iv)将每个反应孔的DNA分子片段化并在每一端连接一对环状接头或使用转座酶递送接头进行标签化以形成适配的DNA片段，其中环状接头或转座酶递送接头包括列索引(Ci)序列或行索引(Ri)序列，其中所述Ci序列对于所述多孔阵列板的一列中的每个反应孔的每个环状接头或转座酶递送接头是通用的，或者其中每个Ri序列对于所述多孔阵列板的一行中的每个反应孔的环状接头或转座酶递送接头是通用的；

vi)通过对适配的DNA片段进行索引PCR来提供索引的DNA文库，其中使用正向和反向索引引物扩增适配的DNA片段以形成索引的PCR产物，其中行索引(Ri)序列或列索引(Ci)序列通过每个正向和反向索引引物引入到适配的DNA片段的每一端，使得所得的索引的PCR产物包括一对对于一列的每个孔都是通用的侧翼列索引(Ci)序列和一对对于一行的每个孔都是通用的侧翼行索引(Ri)序列两者；以及

任选地，其中正向和反向索引引物进一步将各自的5'和3'测序接头提供到适合用于测序反应的索引的PCR产物上。

核酸可以是DNA或RNA。在一个实施方式中，核酸是基因组DNA。在另一实施方式中，核酸是可以mRNA。

根据本发明的另一方面，提供了一种制备用于单细胞或细胞群的全基因组测序的索引的DNA文库的方法，以识别单细胞或细胞群的基因组中的单核苷酸变体、确定单细胞或细胞群的基因组中的染色体结构变异或确定单细胞或细胞群的基因组中的定相信息，该方法包括：

i)提供包括多行和多列反应孔的多孔阵列板；

ii)提供单细胞或细胞群的基因组DNA，其中基因组DNA分布在多孔阵列板上的多个反应孔中，使得每个反应孔的任何给定位点的单链基因组DNA分子不超过一个，

iii)对每个基因组DNA分子进行全基因组扩增(WGA)，以在每个反应孔中提供基因组DNA分子的多个拷贝；

iv)将每个反应孔的DNA分子片段化并在每一端连接一对环状接头或使用转座酶递送接头进行标签化以形成适配的DNA片段，其中环状接头或转座酶递送接头包括列索引(Ci)序列或行索引(Ri)序列，其中所述Ci序列对于所述多孔阵列板的一列中的每个反应孔的每个环状接头或转座酶递送接头是通用的，或者其中每个Ri序列对于所述多孔阵列板的一行中的每个反应孔的环状接头或转座酶递送接头是通用的；

vi)通过对适配的DNA片段进行索引PCR来提供索引的DNA文库，其中使用正向和反向索引引物扩增适配的DNA片段以形成索引的PCR产物，其中行索引(Ri)序列或列索引(Ci)序列通过每个正向和反向索引引物引入到适配的DNA片段的每一端，使得所得的索引的PCR产物包括一对对于一列的每个孔都是通用的侧翼列索引(Ci)序列和一对对于一行的每个孔都是通用的侧翼行索引(Ri)序列两者；以及

任选地，其中正向和反向索引引物进一步将各自的5'和3'测序接头提供到适合用于测序反应的索引的PCR产物上。

可以对索引的核酸进行测序，例如如本文所述。

因此，根据本发明的另一方面，提供了一种单细胞或细胞群的全基因组测序方法，以为识别单细胞或细胞群的基因组中的单核苷酸变体(SNV)提供数据，该方法包括：

i)通过实施根据本发明的方法制备索引的DNA文库，或提供根据本发明的方法制备的索引的DNA文库；

ii)对索引的DNA文库进行测序以提供用于确定单细胞或细胞群的基因组中的任何单核苷酸变体(SNV)的数据。

测序数据可用于确定SNV，例如如本文所述。另外或替代地，测序数据可用于确定与染色体结构变异有关的遗传变化。染色体畸变可以包括数目畸变和/或结构畸变。

另外或替代地，测序数据可用于确定细胞或细胞群中的定相信息。

如说明书和/或附图中公开的，本发明还可以包括单独或组合的一个或多个特征。

定义

术语“空间相关细胞”被理解为是指彼此直接相邻的细胞。

术语“假阳性(FP)突变”或“假阳性(FP)SNV”被理解为是指在从完整细胞中提取DNA之前基因组中不存在的变体核苷酸，例如，假突变可以是损伤引起的错误或复制错误。

术语“真突变/SNV”或“真阳性突变/SNV”可以互换使用，并且被理解为是指在DNA提取之前存在于活细胞基因组DNA中的变体核苷酸。

术语“单核苷酸变体”(SNV)可包括单核苷酸多态性(SNP)或序列中任何其他变异，例如突变。突变或变异可包括给定序列中的核苷酸取代、增加或缺失。

“染色体畸变”被理解为染色体DNA的缺失的、额外的或不规则的部分。它可以来自一个或多个染色体中典型数量的染色体或结构的异常。它们包括各种畸变，例如缺失、重复和插入。会发生平衡畸变，例如倒位以及染色体间和染色体内易位。此外，可能发生可移动元件插入、片段重复、多等位基因染色体数目畸变。上述的最终多重组合可以产生复杂的重排。

“定相(phasing)”被理解为将等位基因(A、C、T和G)分配给父和母系染色体的任务或过程。定相有助于确定是在父系方还是在母系方匹配，是在双方匹配还是在双方都不匹配。定相也有助于染色体映射的过程——将片段分配给特定的祖先。

本领域技术人员将理解，在适当的情况下，本发明的一个实施方式或方面的任选特征可适用于本发明的其他实施方式或方面。

现在将参考附图，仅以示例的方式，更详细地描述本发明的实施方式。

图1.DigiPico测序原理、工作流程和性能。(A)WGS方法只能识别肿瘤(红色和蓝色)中显性扩增克隆中的早期突变过程(EM)。当前活跃的突变过程(CM)导致各种各样具有不同克隆特异性突变的亚克隆。这种多样性决定了肿瘤的进化轨迹。(B)在WGA之前进行模板分区，以便每个隔间(compartment)从每个位点接收不超过一个DNA分子，以识别人工突变。由于损伤引起的错误(红色)和复制错误(蓝绿色)在复制过程中随机发生，因此人工突变导致双等位基因隔间。注意，真突变始终存在于同一隔间内的所有产物DNA链中。(C)DigiPico测序工作流程。LCM：激光捕获显微切割。(D)使用来自EvaGreen标记的DNA的端点相对荧光单位(RFU)确保模板和WGA过程在整个板上的均匀分布。将RFU值标准化以在每次运行中达到中值1。(E)使用Illumina接头引物(P5和P7)的每个孔qPCR测量每个孔中接头连接产物的相对数量。将Ct值归一化以在每次运行中达到中值0。(F)简化Digipico库制备过程需要微小化的WGA，该WGA可以专门地和敏感地扩增每个孔中的亚皮克量的DNA。数值代表9次重复的平均RFU值。误差线代表SD。(G、H和I)如图所示，对每次运行中各个孔的DigiPico测序数据的初步分析证实了测序的高质量和映射率、覆盖深度和覆盖广度的均匀性。(J)DigiPico(UTD)变体独有的定义。从相应的DigiPico数据中减去在标准WGS数据中可识别的SNV导致UTD变体。这些将主要由人工突变以及一些克隆特异性突变组成。由于运行D1110中的模板实际上是用于标准WGS的模板的子集，因此预计DigiPico运行D1110中的所有真实变体也将出现在标准WGS数据中。相比之下，由于深度限制，标准WGS数据中可能不存在运行D1111中的克隆特异性变体，即使支持此类变体的DNA分子可能以非常低的频率存在于大批量(bulk)DNA样本中。在所有箱线图中，水平线代表中位数。方框代表四分位距(在第25百分位和第75百分位之间)。晶须(whisker)代表不包括异常值的范围。异常值定义为高于或低于四分位距1.5倍的数据点。

图2.MutLX算法、设计和结果。(A)比较运行D1110中支持各种突变类型的孔数证实，正如假设的那样，大多数UTD仅存在于少数孔中。水平线代表中位数。方框代表四分位距。晶须示出不包括异常值的范围，这些异常值被定义为在四分位间距范围的1.5倍之外。(B)类似地，与真实变体相比，UTD的双等位基因间隔率似乎明显更高。该值是通过将同时存在变体和参考等位基因的孔数除以具有变体等位基因证据的孔总数来计算的。(C)显示了使用ANN分析DigiPico数据的主要挑战的图表。每个圆圈/星号表示一个变体。红线显示了分类模型的行为。通过模型预测线上方和/或左侧的所有变体都是真实的变体。对没有克隆特异性变体的样本进行分析将导致真实突变和人工突变之间的精确分离。相反，对具有真实克隆特异性突变的样本进行分析会导致未达到最佳标准的模型，这可能导致对真实UTD的过度拟合。这将强制执行模型，该模型以丢失几乎所有克隆特异性变体为代价删除所有FP识别。(D)显示了MutLX中的两步训练过程的图表。第一次训练步骤识别UTD中一些被错误标记的真实突变(灰色圆圈)。从第二次训练中的分析中暂时删除所有可能被错误标记的数据点(颜色是黑色)，以便获得更好的模型来为所有突变分配概率分数。最后，将从模型中获得的概率分数与这些概率分数的不确定性估计(如E中所述)相结合，可以有效消除FP识别，同时保持对真实克隆特异性变体的优异的敏感性。(E)显示了用于计算概率分数的不确定性估计的测试时间缺失(drop-out)分析的图表。黑色神经元表示在缺失分析期间已关闭的神经元。仅接受具有高概率分数和低不确定性分数的变体应该允许消除FP变体识别。(F)显示了运行D1110、D1111、DE011和GM12885的MutLX分析输出的ROC曲线。圆圈代表由MutLX确定的默认截止值。(G)表示SCcaller、Platypus和MutLX的输出中通过的UTD数量的条形图。由于在运行D1110、DE011和GM12885中预计不存在真实UTD，因此这些运行中的UTD数量代表每种分析方法的FP率。Platypus的值基于应用MutLX之前DigiPico特定的过滤标准。

图3.使用DigiPico/MutLX识别活跃的突变过程。(A)HGSOC患者#11152中肿瘤进化的示意图。各种肿瘤样本的标准大批量WGS在所有位点中鉴定出～11000个共有体细胞突变。紫色虚线表示所研究的肿瘤样本的最近共同祖先发生分叉的点。大批量测序(bulksequencing)还分别鉴定了近5000、3000和2000个亚克隆突变，这些突变分别对化疗前网膜肿块、PT2R复发和PALNR复发具有特异性。然而，这些突变可能在这些克隆扩增期间的任何时候发生，并且偏向于较旧的突变。这是由于识别低流行率体细胞突变的局限性。然而，对五个化疗前肿瘤小岛、PT2R和PALNR复发位点进行DigiPico测序，识别出在这些样本中的每一个样本中不同数量的最近出现的克隆特异性突变(以红色数字表示)。PT2R中克隆特异性变异的数量显著更高表明存在活跃的突变过程。(B)在D1111运行中，17号染色体上存在强烈的克隆特异性kataegis事件突出了这种活跃的突变过程。Y轴表示对数尺度下连续体细胞突变的成对距离。仅显示了来自17号染色体的突变。与亚克隆Kataegis事件有关的突变在方框中突出显示，几乎所有突变都是以链特异性C>T或C>G突变的形式。这表明APOBEC酶参与了这种超诱变过程。(C)kataegis中涉及的一些突变的代表性示例。基因组正向链上所有突变的存在进一步证实了超诱变事件的参与(图13)。

图4.识别近期突变的挑战。虽然可以很容易地从肿瘤的大量测序数据中研究旧突变，但由于所涉及突变的变异等位基因分数(VAF)较低，因此从这些数据中研究近期突变受到了阻碍。因此，启发式过滤标准不足以识别近期突变。对近期突变的可靠研究需要研究通过激光捕获显微切割(LCM)分离的单个癌细胞或肿瘤小岛(图1)。然而，此类样本中模板数量有限的WGA导致大量假阳性变体识别，从而阻碍小岛特异性变体的鉴定。我们的分析流水线MutLX可以通过从DigiPico测序数据中消除FP变体识别来克服这个问题。

图5.DigiPico数据的分析工作流程。(1)来自正常组织、大批量肿瘤和DigiPico文库中的下一代测序读数首先映射到人类基因组以生成bam文件。DigiPico读数被分为384个FastQ文件，一个FastQ文件用于384孔板的每个孔。(2)来自DigiPico的384个单独的bam文件合并为一个bam文件。(3)使用Platypus变体识别程序对384个单独的bam文件以及合并的bam文件执行从头(De novo)联合变体识别。添加合并的bam文件确保在变体识别期间不会错过每孔覆盖率低的变体。(4)然后将得到的重新DigiPico变体用作参考用于来自正常组织和大批量肿瘤的标准WGS数据中的变体重识别(re-calling)。(5)然后，通过消除在标准WGS数据中具有支持读数的任何变体，变体重识别数据可用于提取对DigiPico(UTD)所特有的变体。(6)标准WGS数据也用于使用GATK进行变体识别，以获得高置信度的种系SNP列表。(7)该列表将用作从DigiPico数据中提取TP变体识别的指南。为此，在大量血液样本中使用GATK鉴定的任何变体，以及使用Platypus在DigiPico数据中鉴定的任何变体都被假定为真实的。(8)然后通过MutLX使用UTD和种系SNP来训练二元分类模型，以从UTD中提取克隆特异性变体(图6)。

图6.MutLX分析算法。(1)通过从DigiPico数据中减去WGS数据来鉴定UTD变体。(2)UTD和SNP作为训练集来训练初级二元分类模型。(3)该模型用于训练集的初步分析，其允许(4)生成改进的训练集。(5)改进后的训练集然后用于生成分类模型，其(6)可用于UTD的分析。(7)“概率分数(probability score)”表明突变为真实的可能性，以及(8)使用模型计算每个突变的“不确定性分数”，作为计算的概率分数的不可靠性的量度。(9)真正的变体由具有高确定性的高概率分数(低不确定性分数)来鉴定。

图7.运行D1110、D1111、DE011和GM12885的概率分数值。0.2的截止值从UTD中删除了大多数FP变体检出，同时保留了所有样本中的几乎所有种系SNP。

图8.数据模拟证实AUC与真实UTD变体的数量负相关。在运行D1110和DE111中，将各种数量的体细胞突变人为地错误标记为UTD(UTD*)以实现1％、5％和10％的UTD*/UTD比率。由于这两次运行都是对来自大批量肿瘤样本的200pg纯化DNA进行的，因此预计它们都不会具有真正的UTD变体。使用MutLX对每个模拟数据集的独立分析证实了数据集中真实UTD的数量与AUC的数量之间的负相关。值得注意的是，低至1％的UTD*/UTD比率(运行D1110和DE11中分别有16个和36个变体)似乎可以降低AUC，这表明即使存在一小部分真实克隆特异性变体也会扰乱ROC曲线。

图9.合成DigiPico数据集的分析。在运行(A)D1110和(B)DE111中，将各种数量的高置信度体细胞突变人为地错误标记为UTD(UTD*)，从而以各种UTD*/UTD比率综合放大真实UTD的数量。由于这两次运行都是对来自大批量肿瘤样本的200pg纯化DNA进行的，因此预计它们都不会具有真正的UTD变体。结果表明，在大多数人工UTD中存在真正的UTD变体似乎不会损害MutLX生成的分类模型的完整性。每个箱线图显示了用于分析的10个不同UTD*子集的结果。为了在运行之间实现可比较的FP率，分析的截止值用于在运行D1110(A)和DE111(B)的所有合成数据集的种系变异中分别实现90％和95％的TPR。箱线图显示中位数、四分位距和不包括异常值的范围。异常值定义为高于或低于四分位距的1.5倍。

图10.运行D1111中鉴定的一些克隆特异性变体的靶向测序。目标位点的扩增子测序在MiSeq平台上进行。14个目标中的3个似乎在血液样本(以橙色突出显示)中具有高噪声水平，因此其被认为是不确定的。在剩余的11个突变中，只有1个似乎在PT2R肿瘤(以蓝色突出显示)的大批量DNA样本中没有任何证据。VAWF：变异等位基因孔分数(Variant AlleleWell Fraction)。

图11.运行DE111中鉴定的一些人工变体的靶向测序。

图12.DigiPico数据中MutLX识别的FP检出中各种突变类型的频率。绿点(在左侧)是从来自患者#11152、#11513和OP1036的标准WGS数据中的体细胞变体的分析获得的。红点(在右侧)从来自于通过MutLX从同一患者的DigiPico数据中识别的所有人工突变。黑线代表每组值的中位数。MutLX消除的突变中C>A突变的比例较高，这与先前的研究一致，表明文库制备过程中的氧化DNA损伤导致人工C>A突变的形成。

图13.在PT2R样本的亚克隆kataegis中识别的SNV的IGV图像。注意，几乎所有突变都是以基因组正向链上的C>T或C>G突变的形式。

图14.DigiPico和DigiPico2工作流程的比较。(A)DigiPico工作流程耗时近12小时，由7个步骤组成，其中5个步骤发生在384孔板中。(A)DigiPico2工作流程仅耗时4.5小时，由5个步骤组成，其中只有3个步骤发生在384孔板中。蓝色反应发生在384孔板格式中，绿色反应发生在16个孔中，橙色反应发生在1个管中。

图15.DigiPico和DigiPico2索引策略的比较。(A)2个退火索引寡核苷酸的不对称连接引入了一个i5索引和一个i7索引。这些索引的组合可以在DigiPico中产生384个不同的索引。注意，每条链都接收一个i5索引和一个i7索引，因此不存在冗余。(B)在DigiPico2中，最初的列索引(Ci)通过环状接头的有效连接引入。接下来，使用行索引引物(Ri)进行索引PCR。注意，每条链都会收到两次Ci和Ri索引，这会引入去除索引交叉污染物所需的冗余。

图16.DigiPico和DigiPico2结果的比较。(A)两种WGA产物似乎在整个板中都是相对均匀的。数值代表来自测量EvaGreen(RFU)的相对荧光。(B)在DigiPico2中，板中每个孔的索引频率似乎与WGA产物的RFU值具有更好的相关性。(C)这一事实可以使用相关的图来量化。(D)DigiPico2数据的MutLX分析似乎可以更好地区分真实和人工突变。注意图的右上部分中较低的突变存在。该区域标记了在分析中获得高概率分数的人工突变。

图17.单细胞DigiPico(ScDigiPico)测序方法的评估。为了评估ScDigiPico方法在识别活跃的突变过程中的有效性，我们通过使用N-乙基-N-亚硝基脲(ENU)诱变培养的Kuramochi细胞来模拟这样的过程。ENU是一种烷化剂和非常有效的诱变剂，优先致使T>C、T>A和C>T突变。在这种设置下，预计每个细胞都会获得一组不同的突变，但由于诱变的潜在机制是相同的，因此预计突变属于同一类型。(A)将培养的Kuramochi细胞暴露于0.1g/LENU，持续48小时。大多数细胞在ENU的存在下死亡，但那些存活下来的细胞会积累大量的突变。然后将单个突变的Kuramochi细胞分选到384孔板中第一列的每个孔中。在ScDigiPico期间，来自每个细胞的DNA含量在WGA和文库制备之前被均匀分布到该行的所有孔中。(B)分析突变细胞的ScDigiPico结果，正如预期的那样，新突变的频率与ENU的频率明显匹配。(C)在紫外线照射后分析细胞中的仅有一个细胞中，我们在7号染色体上发现了一个明显的Kataegis事件。

图18.DigiPico/MutLX消除了全基因组扩增DNA的假阳性。点代表来自癌症患者或全基因组扩增血液DNA的测序(从皮克血液DNA开始)的单个样本(约20个癌细胞)。与标准DNA测序相比，种系(血液)DNA不应有数千个独特的变体。然而，现有的方法发现了数以万计的这种假阳性突变。相比之下，DigiPico/MutLX消除了这些假阳性。

实施例1-使用DigiPico/MutLX以前所未有的准确度揭示肿瘤中活跃的突变过程

概要

大批量(bulk)全基因组测序(WGS)可以分析肿瘤进化，但由于深度限制，只能识别旧的突变事件。发现当前用于预测肿瘤进化轨迹的突变过程需要对单个克隆或单细胞进行密集测序。然而，此类研究本身就存在问题，这是因为当对皮克量的DNA进行测序时，会发现过多的假阳性突变。数据池化以增加对已发现突变的置信度，将该发现追溯到过去的共同祖先。在这里，我们报告了一个稳固的全基因组测序和分析流水线(DigiPico/MutLX)，其几乎消除了所有假阳性结果，同时保留了极好比例的真阳性。使用我们的方法，我们首次在复发性卵巢癌的约30个癌细胞组中鉴定了超突变(kataegis)事件。这是无法从大批量WGS数据中鉴定出来的。总体而言，我们提出DigiPico/MutLX方法作为一个强大的框架，以前所未有的准确度鉴定克隆特异性变体。

引言

在这项工作中，我们开发了一种单DNA分子WGA和测序方法，以从在临床样本中获得的皮克数量的DNA中获得高质量和数据丰富的测序结果(我们称为DigiPico；用于皮克DNA的数字测序)。此外，我们使用基于人工神经网络(ANN)的算法(MutLX，用于突变学习)对DigiPico数据实施了补充分析工作流程，以消除假阳性结果，同时在全基因组规模上保持对真阳性突变的优异灵敏度。我们使用来自单个患者的广泛测序肿瘤的数据验证我们的方法，其中从三个不同时间点对DNA的45次全基因组测序运行中获得的累积深度约为4200倍。我们通过对来自另外4名癌症患者和淋巴母细胞系的样本进行测序，展示了这些方法的多功能性。

材料与方法

患者样本和同意书

患者#11152、#11502和#11513提供了参与前瞻性生物标志物验证研究妇科肿瘤靶向治疗研究(Gynaecological Oncology Targeted Therapy Study)01(GO-Target-01)的书面同意书，研究伦理批准号为11/SC/0014。患者OP1036参与了前瞻性牛津卵巢癌预测化疗反应试验(Oxford Ovarian Cancer Predict Chemotherapy Response Trial)(OXO-PCR-01)，研究伦理批准号为12/SC/0404。酌情获得了研究参与者的必要知情同意书。手术当天采集血样。在腹腔镜检查或减积手术期间对肿瘤样本进行活检，并立即在干冰上冷冻。所有样本都储存在-80℃冰箱中清晰标记的冷冻管中。

细胞系

GM12885类淋巴母细胞系(RRID：CVCL_5F01)获自Coriell研究所，细胞按照提供者的建议进行培养。

切片和LCM

将冷冻肿瘤样本包埋在OCT(NEG-50，Richard-Allan Scientific)中，并在CryoStar低温恒温切片机(ThermoFisher Scientific)中使用MB DynaSharp切片机刀片(ThermoFisher Scientific)获取10-15μm切片。然后将肿瘤切片转移到PEN膜玻璃载玻片(Zeiss)上并立即在冰上染色(在70％乙醇中2分钟，在50％乙醇的1％甲酚紫(Sigma-Aldrich)中染色2分钟，然后用100％乙醇冲洗)。使用PALM激光显微切割系统(LaserMicrodissection System)(Zeiss)将单个肿瘤小岛弹射到200μl不透明AdhesiveCap(Zeiss)中。

标准WGS和数据分析

使用DNeasy血液和组织试剂盒(Qiagen)提取DNA。使用Covaris S220聚焦超声仪仪器将多达1μg DNA稀释在50μl水中进行片段化，以获得250-300bp片段。然后按照制造商的方案，使用NEBNext Ultra II文库制备试剂盒(NEB)将所得的DNA片段用于文库制备。生成的文库在Illumina NextSeq或HiSeq平台上进行测序，深度为人类基因组的30-40倍。FastQ格式的测序读数最初使用TrimGalore(14)进行修剪，然后使用Bowtie2映射到人类hg19基因组(15)。使用GATK的HaplotypeCaller进行种系变体识别(16)。使用变体等位基因分数截断值为0.2(17)的Strelka2检出体细胞变体。

DigiPico测序

首先使用来自Repli-g单细胞试剂盒(Qiagen)的5μl D2缓冲液，将200pg的纯化DNA、20-30个重悬的细胞核或激光捕获显微切割的肿瘤小岛进行变性。在室温下培养5分钟后，将95μl水添加到样本中，然后使用Mosquito HTS液体处理器(TTP Labtech)将200nl变性模板添加到384孔反应板的每个孔，每个孔中已经包含800nl的WGA混合物(0.58μl Sc反应缓冲液、0.04μl Sc聚合酶(REPLI-g单细胞试剂盒，Qiagen)、0.075μl 1mM dUTP(Invitrogen)、0.04μl EvaGreen 20x(Biotium)和0.065μl水)。将板在30℃培养2小时，然后在65℃热灭活15分钟。如果需要，在反应中添加EvaGreen，允许使用实时PCR机器监测WGA反应(18)。接下来，在没有任何纯化步骤的情况下，对全基因组扩增的DNA依次进行受控酶促片段化(19)反应步骤。简言之，(A)加入1200nl UDG混合物(在1.8倍的NEBuffer 3中，0.08U/μl rSAP(NEB)、0.2U/μl UDG(NEB)、0.4U/μl EndoIV(NEB))，在37℃培养2小时，并在65℃热灭活15分钟。(B)加入1200nl PolI混合物(0.4U/μl DNA聚合酶I(NEB)、0.25mMdNTP、8mM MgCl2和0.8mM DTT)，在37℃培养1.5小时，并在70℃热灭活20分钟。(C)加入1200nl Klenow混合物(0.5U/μl Klenow exo-(NEB)、0.5mM dNTP、8mM MgCl2和0.8mMDTT)，在37℃培养45分钟，并在70℃热灭活20分钟。(D)将400nl的20μM的具有孔特异性索引的全长llumina接头寡核苷酸(表S1)添加到每个孔中，然后添加1100nl的连接混合物(40U/μl T4 DNA连接酶(NEB)、5mM ATP、11.5％PEG 8000(Qiagen)和6.8mM MgCl2)，在20℃培养30分钟，并在65℃热灭活15分钟。

然后合并所得的产物并使用等体积的异丙醇析出DNA。然后将DNA重悬在水中，并使用Agencourt AMPure XP SPRI磁珠(Beckmann coulter)，以0.45倍的珠比用于左侧选择和额外的0.32倍珠比用于右侧选择，对产物进行双重尺寸选择。然后将纯化的DNA重悬在水中，并立即用于使用P5和P7引物混合物进行有限循环PCR扩增(表S1)。PCR进行12个循环，在55℃退火10秒，在72℃延伸45秒。最终产物以0.9倍的比例进行磁珠纯化。然后在Illumina测序平台上以2x150双端测序模式对生成的文库进行测序，以在人类基因组上实现30-40倍的覆盖深度。目前，DigiPico文库制备所需的额外处理步骤使总试剂成本增加了近250英镑。

DigiPico测序数据分析

DigiPico测序数据的分析流程见补充图5。简而言之，在对Illumina序数据进行多路分解后，获得了每个孔的384个配对读取FastQ文件。FastQ文件对接头序列和质量进行了修剪(14)。每个读数的前12个核苷酸也被去除。使用激活了ignore-quals参数的Bowtie2(15)将修剪后的读数映射到人类hg19参考基因组，并使用Picard工具(20)标记重复读数。然后使用Platypus变体识别程序对所有384个单独的bam文件以及来自所有孔的合并bam文件执行联合变体识别(21)。接下来，通过应用质量过滤器(QUAL>60、FR>0.1、HP≤4、QD>10和SbPval≤0.95)去除所有低质量变体。此外，确定了覆盖每个位点的孔总数(Tw)和支持每个变体的孔数(Vw)，并应用孔计数过滤器(Tw>5、Vw>2和Vw/Tw>0.1)以只保留高置信度位点用于分析。最后，使用VCFtools(23)从分析中删除了所有具有不良映射能力的基因组区域(22)。然后，使用Platypus(minPosterior参数设置为0，minMapQual参数设置为5)，使用所得到的高置信度从头DigiPico变体列表对来自血液和大批量肿瘤的WGS数据执行变体重新识别(基因分型)。将在标准WGS数据中均肯定(confidently)不受支持的任何变体提取为UTD(DigiPico独有的)变体。将在血液样本的大批量测序数据(基于GATK分析)中也肯定存在的任何变体提取为TP(真阳性)变体(图5)。

MutLX算法

MutLX分析流水线总结在图6中。

人工神经网络架构

本研究中使用的神经网络模型是多层感知器，其输入层由N个神经元(N＝41)组成，其中N是每个实验中使用的特征数，并使用Keras(24)在Python3中实现。该模型具有两个带有ReLU激活的隐藏层。我们改变了这些数目，但在使用大量神经元时没有看到任何显著改善。最后一层是具有S形(sigmoid)激活的单个输出神经元。损失函数是二元交叉熵。对于训练，我们应用了具有动量的随机梯度下降优化(Adam(25))，学习率为0.001，批大小为8，10个时期(epoch)。在10个时期之后，我们没有观察到任何额外的性能改进。

用于训练的特征

从DigiPico数据的Platypus输出中提取的以下特征用作神经网络模型的输入：

Platypus质量参数：QUAL、BRF、FR、HP、HapScore、MGOF、MMLQ、MQ、QD、SbPval、NF、NR、TCF和TCR(21)。

序列环境复杂度：F₂₀[1]、F₂₀[2]、F₂₀[3]。其中F₂₀[i]是在变体位置两侧10bp序列中i个最丰富的核苷酸的频率之和。

读取分布数据：R_合并[ref+var]、R_合并[var]、W[R[ref]>0s和R[var]＝0]、W[R[ref]>0][0/0]、W[R[ref]>0][0/1]、W[R[var]>0]、W[R[var]>0][1/1]、W[R[var]>0][0/1]、W[R[var]>0和R[ref]＝0][0/1]、W[R[ref]>0和R[var]>0]、W[R[ref]＝0和R[var]>0]、W[R[ref]＝0和R[var]>1]、W[R[ref]＝0和R[var]>2]、W[R[ref]＝0和R[var]>3]、W[R[ref]＝0和R[var]>4]、W[R[ref]＝0和R[var]>5]、R_max[1][var]、R_max[2][var]、R_max[3][var]、R_max[1][ref+var]、R_max[2][ref+var]、R_max[3][ref+var]、Max_c+Max_r、W[R[var]>0]-(Max_c+Max_r)。其中R_合并[x]表示合并的bam文件中支持等位基因x的总读数(ref表示参考等位基因，var表示变体等位基因)。其中示出的W[i][j]表示与报告的基因型j匹配标准i的孔数。在标准i中，R[x]表示支持等位基因x的特定孔中的读数。R_max[y][x]显示支持等位基因x的第y个最高读数的孔中的读数。最后，Max_c是支持变体等位基因的孔数最多的列中的变体支持孔数，Max_r是支持变体等位基因的孔数最多的行中的变体支持孔数。

使用MutLX进行训练

对于每次DigiPico运行，我们将完整的训练集视为所有UTD变体(标记为0)和杂合种系SNP(标记为1)的集合。该集中，UTD变体的数量远小于杂合种系SNP，使得该集不平衡。因此，为了避免在训练中偏向特定标签，我们为每次DigiPico运行创建25个不同的平衡训练子集。这样做是为了使每个训练子集由所有UTD变体和随机选择的杂合种系SNP子集(其大小等于UTD变体的数量)组成。如前所述，大多数UTD变体是FP变体检出，其中真实的克隆特异性变体比例未知，因此使0标签噪声。为了考虑这些噪声标签执行两步训练，我们采用以下策略。在每个平衡的训练子集上训练初始模型后，将得到的模型应用于整个训练集中的突变，以获得每个突变的初始概率值。这些概率值表示属于标签1类别的突变的预测概率。因此，任何获得接近1的预测概率值的0标记突变都可能是错误标记的突变。因此，为了降低训练集中错误标记数据的水平，所有概率值大于0.7的UTD变体和所有概率值小于0.3的种系SNP都被视为错误标记的，并从训练集中删除。此步骤中的截止值是通过对各种模拟数据集的分析凭经验确定的。最后，遵循与初始训练类似的子采样策略，在训练集的剩余突变上训练新模型。然后将该模型用于分析所有UTD变体。

概率和不确定性分数的计算

如前所述，在MutLX中，训练过程使用不同随机选择的种系SNP子集重复25次，每次都会产生不同的模型，因此每个突变有25个不同的预测概率值。因此，我们将每个突变的“概率分数”定义为所有其预测概率值的平均值：

其中P_i是从第i个训练子集获得的概率值，n表示子集的数量(n＝25)。

此外，为了获得每个概率值的不确定性估计，我们进行了测试时间缺失分析(26)。将训练好的模型应用于每个突变进行100次迭代，在此期间，神经网络的第一和第二隐藏层分别以0.8和0.7的速率缺失不同的神经元。这个过程为每个突变产生了100个概率值。基于这些值，我们将每个突变的“不确定性分数”定义为来自25个不同子集的缺失(dropout)方差的平均值：

其中σ_i ²是从第i个训练子集的缺失分析中获得的100个概率值的方差，n表示子集的数量(n＝25)。

概率分数高于0.2(图7)的所有变体的不确定性分数用于生成推定的接受者操作特征(ROC)曲线。该曲线是通过考虑不确定性分数的0.0到0.25之间的截止值范围而生成的。在每个截止值，根据相应的UTD数量绘制具有不确定性分数低于截止值的种系SNP的比率。然后在将UTD的数量归一化为0到1之后计算曲线下面积(AUC)。请注意，在真实的克隆特异性变体不是期望的情况下(所有UTD都是FP检出)，假设是一个完美的模型，则该图表示ROC曲线，并且该图的AUC应该接近1。相反，当AUC显著降低时，表明样本中存在真正的克隆特异性变体。使用模拟数据集验证了真实UTD数量和AUC之间的这种负相关性(图8)。基于这些观察，对于ROC曲线表明存在真实的克隆特异性变体(AUC<0.9)的样本，MutLX使用“不确定性分数”截止值，该值导致95％的TPR以提高克隆特异性变体的恢复率(recovery rate)。对于AUC≥0.9的数据集，用于过滤数据的截止值是根据阈值曲线和ROC曲线的交点确定的。

模拟DigiPico数据集的生成和分析

模拟数据用于：(a)验证真实UTD数量与AUC间的负相关性(图8)和(b)确保不会发生对潜在的真实克隆特异性变体的过度拟合(图9)。

为了生成模拟数据集，我们首先使用Strelka2体细胞变体识别程序在来自患者#11152的肿瘤样本PT2R的大批量WGS数据中鉴定体细胞突变。然后在运行D1110的从头变体识别数据中鉴定出这些体细胞变体，且将Tw>6且Vw/Tw>0.45的任何体细胞变体选为高置信度的体细胞变体。接下来，各种数量的随机选择的高置信度体细胞变体被人为地错误标记为UTD(UTD*)，以达到0.01、0.02、0.03、0.04、0.05、0.06、0.07、0.08、0.09和0.1UTD*/UTD比率。然后将生成的合成变体列表独立地用于MutLX分析，并计算每次运行通过MutLX过滤的UTD和UTD*的数量。为了确保稳健的分析，对于每个比率，分析10个不同的体细胞变体的子集。对DigiPico数据DE111也进行了类似的分析，该数据DE111是从患者#11513的腹水样本的大批量DNA提取物中获得的。

MutLX算法的验证

来自患者#11152的肿瘤样本PT2R用于验证MutLX算法。从冷冻标本中宏观解剖一小块肿瘤，并包埋在OCT介质中用于切片。将肿瘤的第一部分(15μm)收集在单独的试管中，并将细胞核重悬在50μl无菌PBS溶液中。测量悬浮液中的细胞核总数，并用来自Repli-gmini WGA试剂盒(Qiagen)的等量D2缓冲液将含有30个细胞核的体积直接变性。产生的粗裂解物直接用于运行D1111的DigiPico文库制备。然后使用DNeasy血液和组织试剂盒(Qiagen)将剩余的肿瘤样本用于大批量DNA提取物。得到的200pg DNA直接用于制备DigiPico文库D1110。使用NEBNext Ultra DNA文库制备试剂盒(NEB)，将1μg DNA用于标准文库制备。在此设置中，仅运行D1111预计会有真实克隆特异性变体。由于运行D1110的模板是批量WGS分析中使用的模板的子集，运行D1110中的几乎所有真实变体也将以相似的频率出现在WGS数据中，因此其不会被鉴定为UTD。类似的逻辑也适用于DigiPico运行DE011和GM12885的结果。由于这些DigiPico运行都是在来自大批量DNA提取物的200pg DNA上进行的，因此预计这些样本中不会出现真实UTD变体。然而，还值得注意的是，由于数据的数字化性质，具有非常低频率(<0.05％)的变体将在运行D1110、DE011和GM12885中显示出膨胀的变体等位基因分数，因为此类变体不太可能出现在超过一个孔中，它们将从基于Vw过滤器的数据中消除。因此，可以安全地假设这些运行中的几乎所有UTD变体都是FP检出。

SCcaller在DigiPico数据上的应用

SCcaller最初是为分析多重置换扩增的单细胞测序数据而开发的(11)。由于DigiPico文库制备还需要对有限数量的模板DNA进行多重置换扩增，因此所得数据与SCcaller的自然输入基本相似。因此，我们使用DigiPico数据的合并bam文件作为SCcaller的输入。对于分析，杂合SNP列表是使用GATK HaplotypeCaller从各自的大批量WGS数据中获得的，并且截止值用于α＝0.01。接下来，所有过滤的SNV都用于对各自的标准WGS数据进行变体重新识别，并且将WGS数据肯定不支持的所有变体都提取为UTD变体。

突变验证

通过与来自独立的测序平台的大批量肿瘤的深度测序数据进行比较，验证通过MutLX分析的变体。通过与从在Complete Genomics测序平台上测序的相同肿瘤块获得的39个深度测序数据集进行比较，来验证来自患者#11152的所有DigiPico数据(27)。这包括三个Complete Genomics批量测序和36个LFR(长片段读取)测序数据。由于网膜肿块的独立测序数据不是从与用于DigiPico测序的肿瘤块完全相同的肿瘤块中获得的，因此通过此类比较对这些运行的验证率预计不会很高。

对于靶向验证，使用primer3工具设计引物以获得包含变体的扩增子(表S1)。扩增子是通过使用

High-Fidelity PCR Master Mix和GC缓冲液在1ng模板上进行16个循环执行2步PCR获得的。在使用NEBNext Ultra II试剂盒进行接头连接和索引之前，然后将来自每个样本的所有扩增子汇集和纯化。将所得的文库在MiSeq平台上进行测序。使用Bowtie2将测序结果映射到人类hg19基因组，并使用Platypus变体识别程序计算支持每个变体的读取数目。

局部超突变(kataegis)分析

为了生成降雨图(rainfall plot)，使用R中的自定义脚本，将17号染色体上连续体细胞突变对之间的距离与它们在每对中第二个突变的基因组位置作图。局部突变簇的存在表明kataegis事件。在这些图中，每个点基于相对于hg19人类参考基因组在该对中第二个突变的突变类型进行着色。

结果

DigiPico测序方法的实施

有或没有先前DNA损伤的扩增错误的关键特征是它们在扩增过程中被随机引入(6,7,28)。因此，我们假设在对单个DNA分子进行扩增和测序时，人工突变将仅存在于从原始单个DNA分子测序所产生的一部分读取中。相比之下，预计真正的变体将出现在所有此类读取中。在WGA之前将模板DNA分区为单独的隔间，这样每个隔间从每个位点接收不超过一个DNA分子，将产生这样的单个DNA分子测序数据(图1B)。由于人工突变会导致具有支持多个等位基因的读取的隔间，因此这种方法能够鉴定这些人工制品并允许消除FP变体识别。此外，这种分区方法还导致每个基因位点的WGA反应独立进行。因此，为WGA过程提供了多个内部复制数据。虽然预计真正的变体会定期出现在复制中，但由于其随机性，人工突变可能会在少数隔间中有限存在。因此，考虑到这两点，与真实突变相比，这种WGA和测序方法可以在隔间中产生明显的人工突变分布模式。由于人工神经网络已证明能够从高维输入中提取复杂模式，因此它们是识别和消除此类数据中的假阳性突变的良好候选者。虽然已经描述了以前的分区和测序方法来获得单倍型信息，但没有这样用于区分真正的突变和人工突变的方法(19,29,30)。

为了充分受益于用于临床样本的准确基因组研究的分区和测序方法的数据丰富性，我们开发了DigiPico测序(图1C)。为了进行DigiPico测序，首先，我们将近200pg DNA(从20-30个人体细胞中获得)均匀分布到384孔板的各个孔中。这确保了来自同一位点的两种不同DNA分子在同一孔中共存的可能性小于10％(19)。在WGA之后，每个孔都被独立处理成索引文库，在池化和排序之前，每个文库都接收唯一的条形码序列(图1C)。因为，在我们的方法中，人工突变的关键区别因素在于其独特的分布模式，因此DNA分子的均匀分布和扩增以及跨孔测序覆盖的一致深度将是至关重要的。实现这种均匀性(homogeneity)确保了真实突变和人工突变分布模式的差异最大化。为确保达到所需的均匀性，在每次DigiPico文库制备过程中，我们都监测WGA反应的进程并量化所有孔的最终结果。前者是通过在WGA反应中添加EvaGreen染料并每5分钟实时监测荧光强度来实现的。EvaGreen是一种嵌入染料，其与DNA的小凹槽结合，因此不会干扰等温WGA反应(18)。对于后者，我们引入了每孔qPCR步骤，以在汇集之前使用接头特异性引物测量每个孔中接头连接片段的相对数量。只有通过这两种均匀性测试的文库才能用于测序(图1D和1E)。重要的是，我们还将WGA反应体积缩小到1μl。这只有在鉴定出兼容的多重置换扩增(MDA)方法后才能实现。比较六种不同的MDA策略，REPLI-g单细胞扩增是唯一满足我们目的所需的灵敏度和选择性的方法(图1F)。反应小型化使我们能够在单个384孔板中简化文库制备过程，而无需使用现成的自动移液仪器进行中间纯化步骤。最后，我们旨在优化用于冷冻临床样本的DigiPico文库制备过程。这是通过直接对通过LCM(激光捕获显微切割)分离的小群相邻细胞(肿瘤小岛)的粗裂解物进行WGA反应来实现的。这种策略确保了基因组材料的最小损失，同时最小化操作时间，从而减少了模板氧化的机会。

DigiPico测序平台从有限的临床样本中生成高质量的文库

在优化了DigiPico文库制备过程的所有必要方面后，我们决定评估从临床样本中获得的DigiPico文库的质量。为此，我们从高度浆液性卵巢癌患者(#11152)获得的冷冻复发肿瘤样本(PT2R)中制备了DigiPico文库D1110和D1111。在本实验中，虽然D1110文库是从PT2R样本的大批量DNA提取物中提取的200pg模板制备的，但D1111文库制备是直接在该肿瘤样本剩余一小块冰冻部分(包含近30个癌细胞)上进行的。每个文库都在IlluminaNextSeq平台上进行测序，以获得150x2双端形式的近400,000,000个读取。对获得的测序数据的初步评估表明，D1110和D1111文库均产生了高质量的测序数据，对人类hg19基因组的总体映射率分别为91.35％和94.27％(图1G)。分析跨板分布的均匀性表明，在运行D1110和D1111中，平均每个孔覆盖近4.4％和4.6％的基因组，每个孔的平均深度分别为1.7x和2.1x，具有出色的跨板的均匀性(图1H和1I)。这累计产生了92.1％和91.1％的覆盖广度，每次运行的深度分别为30x和43x。这些结果证实，DigiPico测序可用于从有限数量的冷冻临床样本中产生高质量的测序数据，并具有出色的覆盖度。

最后，我们评估了我们关于不同突变类型的独特分布模式的最初假设在实际的DigiPico数据集中是否成立。为此，我们假设DigiPico数据集和同一肿瘤样本的标准批量测序数据之间共有的任何变体都必须是真实的变体。这些应主要由种系SNP和克隆体细胞变体组成。因此，根据定义，所有FP变体检出和大多数克隆特异性突变(如果它们存在于所研究的样本中)将属于仅存在于DigiPico数据中而不存在于大批量WGS数据中的变体。为简单起见，这些变体在下文中称为UTD(DigiPico独有的)。因此，考虑到PT2R样本的标准大批量测序数据已经从用于D1110文库制备的相同DNA提取物中获得，D1110 DigiPico运行中的几乎所有UTD变体都应该是人工制品(图1J)。相反，运行D1111中的UTD可能包含一些克隆特异性突变以及人工突变(图1J)。因此，我们使用运行D1110中的UTD变体作为我们分析中人工突变的代表。在运行D1110中比较同一基因座的两个等位基因共存的孔的频率(图2A)以及支持每个变体的孔的数量(图2B)表明，与任何其他类别的突变相比，UTD具有明显更高比例的前者和较低数量的后者(对于两种分析，p<2e-16，单向方差分析，然后是Tukey HSD测试)。这清楚地支持了如假设的DigiPico数据集中人工突变的不同分布模式。

DigiPico数据的MutLX分析流水线

使用DigiPico测序获得高质量数据后，我们决定实施分析流水线，以根据突变的分布模式消除FP变体检出。如前所述，ANN算法理想地适合处理这种复杂模式的问题。给定一组代表性的正确标记的示例(训练集)，ANN可以学习对突变进行分类，而无需任何专用于类的信息。然而，对于从测序数据中消除FP突变的问题，实施ANN算法有两个主要问题；(a)难以获得可推广的模型和(b)无法获得具有代表性的准确标记的训练集。首先，不可能生成可概括用于分析每个DigiPico数据集的模型，这是因为突变的分布模式取决于各种专用于运行的初始条件，这些初始条件不容易解释(例如基因组的拷贝数状态)。因此，需要为每个DigiPico运行量身定制的专用于运行的模型。这意味着需要选择运行特异性突变的子集作为每次DigiPico运行的训练集。其次，虽然可以很容易地从基因组中的已知SNP中提取正确标记的真实突变示例，但不可能确定一组具有代表性且准确的人工突变示例。为了解决这个问题，我们认为UTD是一组有代表性的人工突变的合理近似值，假设UTD主要由此类突变组成。然而，这种假设可能会导致关键挑战。根据定义，UTD由人工突变和真实的克隆特异性突变组成。虽然预计在所有DigiPico运行中都会大量存在人工突变，但根据样本，真实的克隆特异性突变可能以不同的频率出现(图1J)。因此，当UTD都被视为人工突变的示例时，具有更多克隆特异性变体的样本将具有更嘈杂的训练集。如果不考虑这一点，它可能会使具有真实克隆特定变体的样本处于分析劣势，因为更嘈杂的训练集会导致更差的分类模型(图2C)。具体来说，在具有真实克隆特异性变体的样本中，在人工突变的示例中存在真正的突变可能导致模型对这些变体的过度拟合(图2C)。这可能会削弱模型准确鉴定此类样本中真实的克隆特异性变体的能力。鉴于DigiPico测序的主要目标是鉴定克隆特异性变体，因此必须确保在分析DigiPico数据集时不会发生这种过度拟合。

考虑到所有上述限制和问题，我们设计并实施了一个基于ANN的二元分类器MutLX，用于分析DigiPico数据集。DigiPico分析流水线的重点是有效消除FP检出和准确鉴定来自UTD的真实克隆特异性变体。为了解决用不完善训练集训练的问题，我们在训练MutLX时采用了以下方法。最初，我们将所有UTD变体视为人工突变(标记为0)的示例，并将相似数量的随机选择的杂合种系SNP视为真实变体(标记为1)的示例。由于大多数UTD变体是FP检出，真实的克隆特异性变体比例未知，因此在此阶段0标签被认为是“噪声”。换句话说，虽然真实的克隆特异性变体必须被标记为1，但由于它们在此阶段的匿名性，它们仍被标记为0等等。为了适应训练数据集中的这种类型的噪声，我们采用了两步训练过程(图2D)。第一步，在给定所有标记示例数据的情况下训练模型。然后使用这个初始模型来计算每个突变属于其标签类别的概率。基于这些模型预测，任何似乎被错误标记的突变都会从数据集中暂时消除(裁剪)。这里的根本假设是，即使在噪声样本上训练的模型可能不如在假设的干净数据集上训练的模型那么稳健，但它仍然会偏向于对正确示例的更好预测，因为它们与数据集中的噪声相比具有更高的比例。因此，针对其原始标签的模型预测的示例很可能一开始就被错误标记。第二步，在这个裁剪后的训练集上训练新的分类模型。由于第二个训练集可能包含较少错误标记的数据点，因此无论所有UTD变体是否确实是人工制品，最终模型预计更有效地鉴定真实的突变(31，32)。我们接下来使用这个模型为每个假定的突变分配“概率分数”。该分数表示某个突变属于真实变体类别的可能性。虽然这个两步训练过程有望显著改善分类模型，但由于训练集的不完善，最终模型仍然容易出错。因此，我们添加了另一层次的分析，以进一步提高我们流水线的准确性。这是通过将不确定性估计分配给每个突变的“概率分数”来实现的。这种不确定性估计是基于这样一个假设，即ANN隐藏层中的大多数激活的神经元支持稳健的预测。因此，这些神经元的任何子集也会始终产生相似的概率分数，因此，从不同神经元子集获得的各种“概率分数”之间的差异将很小(图2E)。相比之下，人工突变的看似高的“概率分数”最可能仅由ANN隐藏层中的一些神经元支持。因此，不同的神经元子集将导致不同的“概率分数”，这将导致从不同神经元子集获得的人工突变的分数差异很大(图2E)。因此，“不确定性分数”可以计算为从MutLX中经过训练的ANN的多个不同随机选择的神经元子集获得的“概率分数”的方差(variance)(26)。因此，每个突变的“概率分数”和“不确定性分数”的组合应该使我们能够准确地确定被检出的变体是真实的突变还是模板中人工改变的结果(图6)。

MutLX算法的验证

为了验证我们的策略，我们选择在运行D1110和D1111上测试MutLX分析流水线。这是因为这些DigiPico运行是从HGSOC中获得的，该HGSOC之前使用从48个独立全基因组测序数据集获得的数据进行了广泛测序，这些数据集跨越三个不同时间点(患者#11152)，总深度约为4200x，来自两个独立测序平台(33)。据我们所知，这包括迄今为止最广泛的全基因组测序肿瘤。这个格外大的数据集允许对该肿瘤中的突变进行可靠的交叉验证。为此，我们使用MutLX算法分析来自运行D1110和D1111的测序数据。如前所述，当使用来自PT2R位点的大批量测序数据与这些DigiPico数据集进行比较时，预计真正的UTD变体(克隆特异性变体)仅出现在运行D1111中，而运行D1110中几乎所有UTD预计是人工制品(图1J)。此外，我们还分析了从血液样本(运行DE011)的纯化DNA以及培养的GM12885类淋巴母细胞制备的DigiPico测序数据，预计这两种细胞都没有真实的UTD突变。对这些DigiPico运行进行从头变体识别，然后根据孔数进行初始过滤，导致在每个样本中鉴定出数千个UTD变体，几乎所有这些变体都被认为是FP检出。然而，在运行D1110、DE011和GM12885的UTD变体上应用MutLX算法，导致在仅分别对于4、7和3个全基因组范围的FP突变的这些运行中，有效消除超过99％的FP变体检出，同时保持约85％的灵敏度以检测真实的突变(图2F和2G)。相比之下，SCcaller(11)对相同数据的分析分别导致713、712和13,280个FP变体检出(图2G)。另一方面，MutLX在运行D1111中鉴定出264个假定的克隆特异性变体，其中238个(90％)克隆特异性变体通过与该肿瘤样本的独立高深度数据集进行比较得到了验证(图2G)。此外，通过对肿瘤的大批量DNA进行靶向测序，进一步验证了这些观察结果。因此，在11个含有来自运行D1111的克隆特异性变体的分析扩增子之外，发现10个最终以低频率存在于PT2R样本的大批量DNA中(图S7)。此外，对来自运行DE111的37个看似高质量的UTD变体进行扩增子测序(这些变体被MutLX算法标记为人工制品)，表明没有证据表明它们存在于大批量DNA样本中(图11)。这些结果清楚地证实，MutLX可以学习将人工突变与真实变体区分开来的准确分类模型，并能够有效地鉴定DigiPico数据中真正的克隆特异性变体。

此外，我们研究了存在真正的克隆特异性突变是否会由于过度拟合而损害模型的灵敏度。为此，我们人为地将运行D1110和DE111中不同数量的体细胞突变错误标记为人工UTD变体(UTD*)，以生成具有不同比例的真实UTD的合成数据集。然后通过MutLX单独分析这些合成数据集，并在所有合成数据集中检查不同UTD*/UTD比率下的FP率以及UTD*的恢复率。结果表明，高达10％的UTD*/UTD比率不会显著影响UTD*变体的恢复率，表明MutLX中不会发生过度拟合(图9)。

DigiPico/MutLX测序和分析方法的通用性

最后，为了确保我们提出的方法的通用性，对来自四个不同HGSOC患者的各种模板DNA来源进行了DigiPico测序，并使用MutLX算法分析得到的UTD。结果清楚地表明，MutLX可以可靠地鉴定和消除来自各种DigiPico文库的人工变体检出(表1)。这有力地表明DigiPico/MutLX可以有效地研究实体瘤中近期获得的突变。重要的是，分析这些数据中不同突变类型的频率表明，在已鉴定的人工突变中更高水平的C>A突变的存在，这与这种FP检出是模板DNA氧化损伤的结果的观点一致(图12)。

使用DigiPico/MutLX研究活跃突变过程

我们接下来测试了研究HGSOC患者的突变过程的可行性(#11152)。对于该患者，可获得来自化疗前网膜肿块的各种测序数据(以30x标准批量测序，以及5次肿瘤小岛DigiPico运行)。患者随后复发，并从骨盆(盆腔复发性肿瘤；PT2R)和主动脉旁淋巴结(PALNR)收集肿瘤样本，用于标准批量测序以及肿瘤小岛的DigiPico测序。对大量化疗前测序数据的分析鉴定了13,721个体细胞突变。来自DigiPico数据的这些突变中有84.6％存在于至少三个肿瘤小岛中，其中91.4％还存在于先前公布的LFR数据中的至少三个额外的小岛中(33)。突变的高发生率表明它们是在肿瘤中固定的早期突变。对来自肿瘤小岛的DigiPico数据的分析表明，存在数量有限的大体积肿瘤(bulk tumor)中缺少的克隆特异性突变。与其他小岛相比，五个化疗前小岛中的每一个都具有许多真实独特的突变(2、6、8、8和36)，表明它们是最近发生的(图3A)。PT2R复发的大批量WGS数据表明出现了3,009个新的体细胞突变，这些突变在化疗前的大量测序数据、DigiPico数据或LFR数据中不存在。这些突变可能发生在任何时候，因为网膜肿块和PT2R复发的共同祖先彼此偏离(图3A)。对来自患者#11152的复发样本中的肿瘤小岛的分析表明，与主动脉旁淋巴结复发(PALNR)或化疗前肿瘤相比，盆腔复发肿瘤(PT2R)具有高负荷的克隆特异性突变。这一观察结果表明，在该患者中，在SNV诱变下的分子机制最近可能已被激活(图3A)。此外，使用降雨图分析PT2R样本的克隆特异性突变揭示了17号染色体上存在强烈的亚克隆局部超诱变(kataegis)事件(8)(图3B、3C和11)。将构成此kataegis事件的突变与化疗前网膜肿块的批量测序数据、DigiPico数据和LFR数据进行比较，发现仅在DigiPico PT2R数据中发现它们，表明它们是真实的克隆特异性突变。

结论

在这项工作中，我们将DigiPico/MutLX作为一个集成平台，用于在全基因组规模上以前所未有的准确度鉴定来自小细胞群的突变。我们相信，这项工作为发现癌症和正常组织中发生的当前或最近的体细胞突变过程提供了重要的垫脚石。了解当前的突变过程是预测肿瘤进化轨迹的关键，并且可能是在治疗上干扰这些轨迹的关键。在肿瘤的大批量测序中鉴定的突变一定发生在肿瘤从发生到呈现的延长历史中的某个时间点。相反，细胞特异性突变一定发生在该细胞有限的生命周期内。类似地，源自单细胞的小克隆中的突变也是近期的。这种突变的年龄不能超过克隆的年龄，克隆的年龄由生成该克隆所需的细胞分裂次数定义。研究细胞特异性或小克隆特异性突变的模式可以鉴定近期或当前的突变过程(1)。定义这样的过程是非常可取的，因为它们可能与生物或化学现象有因果关系，因此可以产生重要的机理性见解。识别这些机制具有重要的实际意义，因为它们可能适合治疗干预或预测未来的肿瘤行为。现有技术不允许直接准确鉴定来自来自肿瘤的单个细胞或单个小克隆的突变。DigiPico/MutLX首次实现了这一目标。

为了克服主要与发现假阳性突变相关的重大技术缺陷，当前的单细胞WGS分析方法要么需要广泛的验证研究(11)，要么依赖于组合来自多个细胞的数据以获得细胞之间共有的可靠突变(12，34)。然后将这些细胞分组为源自共同祖先的克隆。虽然与大批量测序相比，这些技术适用于更近的共同祖先，但它们仍然不理想，因为从这些方法获得的数据不能反映现有细胞中发生的突变过程。此外，降低每个细胞的测序深度以实现对大量细胞的测序降低覆盖广度，而覆盖广度已经因制备步骤期间遗传物质的损失而受到影响。这增加了需要分析的细胞数量，以推断和鉴定克隆，进而将祖先进一步追溯到过去。此外，在单细胞分析方法中缺乏有关物理(physical)相关性的信息，导致失去对可能来自单个克隆的细胞进行分组的机会。这增加了推断克隆的祖先与现在之间的差距，使得难以定义肿瘤中当前存在的细胞内活跃的过程。

DigiPico/MutLX具有能够保存空间信息的独特优势。分析空间相关细胞，保留物理相关性，并假设物理相关细胞属于单个克隆(9)。还建议定义可能源自组织驻留干细胞的不同结构，以识别和分析克隆。例如，来自单个小肠隐窝或单个子宫内膜腺体的细胞可以合理地预期来自单个组织驻留干细胞(35，36)。在这些情况下，每个解剖单元定义一个克隆，该克隆可能具有或不具有能够与突变驱动因素相关的克隆特异性突变。此外，来自克隆的测序数据可在计算方面用于推断亚克隆且预测可能在克隆中出现的更近期事件。这类似于大批量测序和分析所实现的，但在由有限数量的细胞组成的单个克隆水平上。由于空间转录组学技术的最新发展，保存空间信息也特别有趣(37)。可以想象，将高度准确的DNA测序与空间转录组学相结合，可以剖析组织中的遗传和非遗传异质性。简而言之，目前用于分析小克隆的技术会产生大量假阳性结果，其使得在没有详尽验证的情况下，不能够获得在基因组规模上直接准确的克隆特异性信息。组合来自多个克隆的数据是一种常见的解决方案，但会将祖先进一步追溯到过去。我们之前曾使用这种方法分析少量肿瘤细胞(肿瘤小岛)(33)。由于与来自单个小岛的突变检出相关的不确定性，有必要仅识别所有肿瘤小岛之间共有的突变并仅有效鉴定主要的突变。然后使用靶向测序对大约700个突变进行单独验证。虽然这仍然产生了重要的生物学见解，但我们无法研究小岛特异性突变。DigiPico/MutLX现在使得能够研究这种突变。我们证明了约30个癌细胞的DNA的直接分析如何导致亚克隆kataegis事件的成功鉴定。

总的来说，我们在这里展示了DigiPico和MutLX可以从从临床样本中获得的有限数量的细胞中超准确地识别体细胞突变，这是对现有方法的重要改进。此外，与其他依赖基因组二倍体区域来计算扩增偏差的计算方法不同，我们的方法还与经历广泛拷贝数改变的基因组兼容，例如在HGSOC中。我们相信DigiPico/MutLX方法的通用性能够研究肿瘤和正常组织中的活跃突变过程。

可获性

MutLX的源代码可以在Github上获得(https://github.com/mmdknr/DigiPico)。

登记号

本研究中使用的所有测序数据可在EGA(EGAD0001005118)上获得。

参考文献

1.Turajlic,S.,Sottoriva,A.,Graham,T.and Swanton,C.(2019)Resolvinggenetic heterogeneity in cancer.Nat.Rev.Genet.,10.1038/s41576-019-0114-6.

2.Zhang,J.,

S.S.,Marjani,S.L.,Zhang,W.and Pan,X.(2018)Characterization of cancer genomic heterogeneity by next-generationsequencing advances precision medicine in cancer treatment.Precis.Clin.Med.,1,29–48.

3.Gerstung,M.,Jolly,C.,Leshchiner,I.,Dentro,S.C.,Gonzalez,S.,Mitchell,T.J.,Rubanova,Y.,Anur,P.,Rosebrock,D.,Yu,K.,et al.(2017)Theevolutionary history of 2,658cancers.bioRxiv,10.1101/161562.

4.Barber,L.J.,Davies,M.N.and Gerlinger,M.(2015)Dissecting cancerevolution at the macro-heterogeneity and micro-heterogeneityscale.Curr.Opin.Genet.Dev.,30,1–6.

5.Bohrson,C.L.,Barton,A.R.,Lodato,M.A.,Rodin,R.E.,Luquette,L.J.,Viswanadham,V.V,Gulhan,D.C.,Cortes-Ciriano,I.,Sherman,M.A.,Kwon,M.,et al.(2019)Linked-read analysis identifies mutations in single-cell DNA-sequencingdata.Nat.Genet.,10.1038/s41588-019-0366-2.

6.Chen,L.,Liu,P.,Evans,T.C.J.and Ettwiller,L.M.(2017)DNA damage is apervasive cause of sequencing errors,directly confounding variantidentification.Science,355,752–756.

7.Costello,M.,Pugh,T.J.,Fennell,T.J.,Stewart,C.,Lichtenstein,L.,Meldrim,J.C.,Fostel,J.L.,Friedrich,D.C.,Perrin,D.,Dionne,D.,et al.(2013)Discovery and characterization of artifactual mutations in deep coveragetargeted capture sequencing data due to oxidative DNA damage during samplepreparation.Nucleic Acids Res.,41,e67.

8.Nik-Zainal,S.,Alexandrov,L.B.,Wedge,D.C.,Van Loo,P.,Greenman,C.D.,Raine,K.,Jones,D.,Hinton,J.,Marshall,J.,Stebbings,L.A.,et al.(2012)Mutationalprocesses molding the genomes of 21breast cancers.Cell,149,979–993.

9.Martincorena,I.,Fowler,J.C.,Wabik,A.,Lawson,A.R.J.,Abascal,F.,Hall,M.W.J.,Cagan,A.,Murai,K.,Mahbubani,K.,Stratton,M.R.,et al.(2018)Somaticmutant clones colonize the human esophagus with age.Science,362,911–917.

10.Tubbs,A.and Nussenzweig,A.(2017)Endogenous DNA Damage as a Sourceof Genomic Instability in Cancer.Cell,168,644–656.

11.Dong,X.,Zhang,L.,Milholland,B.,Lee,M.,Maslov,A.Y.,Wang,T.and Vijg,J.(2017)Accurate identification of single-nucleotide variants in whole-genome-amplified single cells.Nat.Methods,14,491–493.

12.Zafar,H.,Wang,Y.,Nakhleh,L.,Navin,N.and Chen,K.(2016)Monovar:single-nucleotide variant detection in single cells.Nat.Methods,13,505–507.

13.Chen,C.,Xing,D.,Tan,L.,Li,H.,Zhou,G.,Huang,L.and Xie,X.S.(2017)Single-cell whole-genome analyses by Linear Amplification via TransposonInsertion(LIANTI).Science,356,189–194.

14.Krueger F.(2016)Trim Galore！

15.Langmead,B.and Salzberg,S.L.(2012)Fast gapped-read alignment withBowtie 2.Nat Meth,9,357–359.

16.McKenna,A.,Hanna,M.,Banks,E.,Sivachenko,A.,Cibulskis,K.,Kernytsky,A.,Garimella,K.,Altshuler,D.,Gabriel,S.,Daly,M.,et al.(2010)The GenomeAnalysis Toolkit:a MapReduce framework for analyzing next-generation DNAsequencing data.Genome Res.,20,1297–1303.

17.Kim,S.,Scheffler,K.,Halpern,A.L.,Bekritsky,M.A.,Noh,E.,

M.,Chen,X.,Kim,Y.,Beyter,D.,Krusche,P.,et al.(2018)Strelka2:fast and accuratecalling of germline and somatic variants.Nat.Methods,15,591–594.

18.Hosokawa,M.,Nishikawa,Y.,Kogawa,M.and Takeyama,H.(2017)Massivelyparallel whole genome amplification for single-cell sequencing using dropletmicrofluidics.Sci.Rep.,7,5199.

19.Peters,B.A.,Kermani,B.G.,Sparks,A.B.,Alferov,O.,Hong,P.,Alexeev,A.,Jiang,Y.,Dahl,F.,Tang,Y.T.,Haas,J.,et al.(2012)Accurate whole-genomesequencing and haplotyping from 10 to 20 human cells.Nature,487,190–195.

20.Picard Tools(2018).

21.Rimmer,A.,Phan,H.,Mathieson,I.,Iqbal,Z.,Twigg,S.R.F.,Wilkie,A.O.M.,McVean,G.and Lunter,G.(2014)Integrating mapping-,assembly-andhaplotype-based approaches for calling variants in clinical sequencingapplications.Nat.Genet.,46,912–918.

22.Derrien,T.,Estellé,J.,Marco Sola,S.,Knowles,D.G.,Raineri,E.,Guigó,R.and Ribeca,P.(2012)Fast computation and applications of genomemappability.PLoS One,7,e30377–e30377.

23.Danecek,P.,Auton,A.,Abecasis,G.,Albers,C.A.,Banks,E.,DePristo,M.A.,Handsaker,R.E.,Lunter,G.,Marth,G.T.,Sherry,S.T.,et al.(2011)The variantcall format and VCFtools.Bioinformatics,27,2156–2158.

24.Chollet,F.and others(2015)Keras.

25.Kingma,D.P.and Ba,J.(2014)Adam:A Method for StochasticOptimization.CoRR,abs/1412.6.

26.Gal,Y.and Ghahramani,Z.(2015)Dropout as a Bayesian Approximation:Representing Model Uncertainty in Deep Learning.arXiv e-prints.

27.Drmanac,R.,Sparks,A.B.,Callow,M.J.,Halpern,A.L.,Burns,N.L.,Kermani,B.G.,Carnevali,P.,Nazarenko,I.,Nilsen,G.B.,Yeung,G.,et al.(2010)HumanGenome Sequencing Using Unchained Base Reads on Self-Assembling DNANanoarrays.Science(80-.).,327.

28.Arbeithuber,B.,Makova,K.D.and Tiemann-Boege,I.(2016)Artifactualmutations resulting from DNA lesions limit detection levels in ultrasensitivesequencing applications.DNA Res.,23,547–559.

29.Amini,S.,Pushkarev,D.,Christiansen,L.,Kostem,E.,Royce,T.,Turk,C.,Pignatelli,N.,Adey,A.,Kitzman,J.O.,Vijayan,K.,et al.(2014)Haplotype-resolvedwhole-genome sequencing by contiguity-preserving transposition andcombinatorial indexing.Nat.Genet.,46,1343–1349.

30.Zheng,G.X.Y.,Lau,B.T.,Schnall-Levin,M.,Jarosz,M.,Bell,J.M.,Hindson,C.M.,Kyriazopoulou-Panagiotopoulou,S.,Masquelier,D.A.,Merrill,L.,Terry,J.M.,et al.(2016)Haplotyping germline and cancer genomes with high-throughput linked-read sequencing.Nat.Biotechnol.,34,303–311.

31.Northcutt,C.G.,Wu,T.and Chuang,I.L.(2017)Learning with ConfidentExamples:Rank Pruning for Robust Classification with Noisy Labels.InProceedings of the Thirty-Third Conference on Uncertainty in ArtificialIntelligence,UAI’17.AUAI Press.

32.Natarajan,N.,Dhillon,I.S.,Ravikumar,P.K.and Tewari,A.(2013)Learning with noisy labels.In Advances in neural information processingsystems.pp.1196–1204.

33.Hellner,K.,Miranda,F.,Fotso Chedom,D.,Herrero-Gonzalez,S.,Hayden,D.M.,Tearle,R.,Artibani,M.,KaramiNejadRanjbar,M.,Williams,R.,Gaitskell,K.,etal.(2016)Premalignant SOX2 overexpression in the fallopian tubes of ovariancancer patients:Discovery and validation studies.EBioMedicine,10,137–149.

34.Laks,E.,Zahn,H.,Lai,D.,McPherson,A.,Steif,A.,Brimhall,J.,Biele,J.,Wang,B.,Masud,T.,Grewal,D.,et al.(2018)Resource:Scalable whole genomesequencing of 40,000single cells identifies stochastic aneuploidies,genomereplication states and clonal repertoires.bioRxiv,10.1101/411058.

35.Moore,L.,Leongamornlert,D.,Coorens,T.H.H.,Sanders,M.A.,Ellis,P.,Dawson,K.,Maura,F.,Nangalia,J.,Tarpey,P.S.,Brunner,S.F.,et al.(2018)Themutational landscape of normal human endometrial epithelium.bioRxiv,10.1101/505685.

36.Lee-Six,H.,Ellis,P.,Osborne,R.J.,Sanders,M.A.,Moore,L.,Georgakopoulos,N.,Torrente,F.,Noorani,A.,Goddard,M.,Robinson,P.,et al.(2018)The landscape of somatic mutation in normal colorectal epithelialcells.bioRxiv,10.1101/416800.

37.Burgess,D.J.(2019)Spatial transcriptomics coming ofage.Nat.Rev.Genet.,20,317.

本文所述的所有参考文献均可通过引用并入。

DigiPico2

实施例2-DigiPico2，一种以前所未有的准确度对皮克量DNA进行全基因组测序的新方法

引言

先前我们将DigiPico文库制备流水线和MutLX分析平台描述为一种从有限数量的临床材料中准确识别单核苷酸变体(SNV)的方法。这是一个重要的方法学进步，主要是因为从临床样本中获得的有限数量的遗传物质必须在测序之前进行全基因组扩增(WGA)。然而，WGA过程在扩增的DNA中引入了多达100,000个人工突变，从而使最终分析结果充满假阳性变体检出，从而阻碍了对原始样本的任何有意义的遗传解释。在DigiPico/MutLX策略中，我们通过在WGA步骤之前将单个DNA分子分离成独立的隔间并在该过程之后对它们进行索引来克服这一障碍。通过这样做，我们将对于真实突变的信息数字化，这意味着每个隔间将携带或不携带突变的等位基因。然而，由于它们在WGA过程中产生的方式，人工突变将导致包含突变的和参考等位基因信息的隔间(图1B)。基于该信息，我们随后开发了一种基于人工神经网络(ANN)的算法MutLX，以有效地识别和消除我们数据中的这些人工突变(图2)。我们在模拟数据和患者样本上广泛测试了我们的策略，并表明我们的方法在消除FP变体检出方面确实有效。

然而，在产生高质量数据的同时，DigiPico文库制备方法几乎没有技术限制。首先，从前述的方法中借用的文库制备(CoREF)的片段化步骤非常复杂且耗时。此外，CoREF要求在WGA过程中使用dUTP。由于dUTP是一种非天然核苷酸，因此它可能会在最终产物中引入更多人工突变。接下来，我们发现在DigiPico中的接头连接效率非常低，这有时会影响文库质量。最后，由于索引数量多且索引信息缺乏冗余，存在索引交叉污染的机会，其可能会对最终结果产生不利影响。因此，我们开发了DigiPico2文库制备方法来解决所有这些问题。

结果

改进DigiPico文库制备工作流程

如前所述，在DigiPico方法中使用dUTP是因为需要CoREF片段过程，这是一种非常复杂的片段策略(图14A)。因此，使用替代的片段方法将同时解决复杂性和dUTP问题。为此，我们决定使用Lotus DNA文库制备试剂盒(美国IDT)提供的现有片段化和末端修复策略。在Lotus DNA文库制备试剂盒中，酶混合物用于以时间依赖性方式将大的DNA分子片段化成更小的段，并准备片段的末端用于连接步骤。为了使这一新策略与DigiPico兼容，我们最初使用I-DOT(Dispendix，德国)分配器来确保所有隔间或多或少同时接收酶混合物。接下来，我们优化了反应条件以达到DigiPico测序所需的片段长度。这样做，我们能够将文库制备时间从原来的12小时减少到4.5小时，并且无需在WGA反应中使用dUTP(图14B)。

接下来，我们旨在解决DigiPico中低连接效率的问题。最初，我们的接头连接和索引依赖于使用不对称连接方法。在这种方法中，具有短互补区域的长索引寡核苷酸用于连接，其效率极低(图15A)。更有效的方法将需要连接到片段两端的环状公共接头。但由于这些接头在其标准形式中不包含任何索引，因此需要单独纯化产物，然后使用索引引物通过PCR反应引入索引。然而，这增加了另一个挑战，因为纯化384种单独产物将非常复杂、耗时，并且会增加交叉污染的重大风险。为了克服这些问题，我们为DigiPico2设计了一种新颖的索引策略。在DigiPico2中，首先在通用接头的茎环中引入一组索引(图15B)。在连接步骤中，板的每列中的所有孔都将接收不同的索引环状共用接头，因此总共24种不同的寡核苷酸足以用第一组索引(列索引)对板的所有列进行索引。在连接步骤之后，每行中的所有孔都将合并到一个单独的管中，从而形成16个不同的合并物。这16个不同的池可以方便地纯化，以用于下一步的索引。在下一步中，来自每个池的纯化产物将使用16种不同的索引引物(行索引)通过PCR反应进行索引。在此阶段结束时，板的每个孔都将收到不同的列-行索引组合。使用这种索引策略不仅会显著提高连接效率，而且还引入了2组冗余，可用于消除数据中的索引交叉污染。在第一组中，列索引将结合到每个片段的两端。因此，任何交叉污染最可能会导致在其末端具有不同索引的片段，因此可以很容易地从数据中删除。在第二组中，每行的索引寡核苷酸可以被双重索引，以便标准索引1(i7)和标准索引2(i5)序列都可以唯一标识特定行。结合这些冗余组，索引交叉污染率可以降低至少2个数量级。

DigiPico2工作流程显着提高了文库质量

为了测试这些修改对数据最终质量的影响，对来自患者11152血液的120pg DNA进行DigiPico2测序。使用这个样本的原因在于我们之前已经对该患者的肿瘤和正常细胞进行了广泛的测序。正如预期的那样，WGA与之前的版本相似，导致产物分布非常均匀(图16A)。然而，在文库制备和测序之后，很明显与DigiPico不同，在DigiPico2中，最终文库中每个孔的表现似乎与WGA产物的数量强相关(图16A-C)。这可能是连接效率提高的直接结果。这种提高的相关性还可以允许引入QC措施，仅基于WGA产物的均匀性，这在以前是不可能的。此外，分析索引冗余信息，我们发现使用索引交叉污染过滤器消除了近5％的读取。在没有这种新过滤器情况下，这些污染物可能会对分析结果产生负面影响。最后，我们使用MutLX算法分析来自DigiPico2的数据。最终结果显示人工突变和真实突变之间有更清晰的区别，表明当使用MutLX算法进行分析时，DigiPico2的性能至少与DigiPico方法一样好(图16D)。

将DigiPico2工作流程扩展到单细胞全基因组测序

建立了DigiPico2工作流程后，我们测试了它是否可以应用于单细胞全基因组测序。这很重要，因为活跃的突变过程很可能在单个细胞内开始。因此，我们通过将单个细胞的DNA分区到384孔板的整行中，引入了单细胞DigiPico(ScDigiPico)测序的工作流程(图17A)。为了评估ScDigiPico方法在识别活跃的突变过程中的效力，我们通过使用N-乙基-N-亚硝基脲(ENU)诱变培养的Kuramochi细胞来模拟这样的过程。ENU是一种烷化剂和非常有效的诱变剂，优先致使T>C、T>A和C>T突变。在这种设置下，预计每个细胞获得一组不同的突变，但由于诱变的潜在机制是相同的，因此预计突变属于同一类型。ScDigiPico能够识别上述突变类型的富集(图17B)。我们还通过在单细胞分选和ScDigiPico文库制备之前用紫外线照射Kuramochi细胞进一步研究了ScDigiPico的潜力(图17C)，并在其中一个细胞中鉴定出Kataegis事件。这些结果累积表明，ScDigiPico是准确鉴定由单个癌细胞中的活跃突变过程引起的真正突变的有效策略。

DigiPico2方案

首先使用来自Repli-g单细胞试剂盒(Qiagen)的5μl D2缓冲液，将200pg的纯化DNA、20-30个重悬的细胞核或激光捕获显微切割的肿瘤小岛进行变性。在室温下培养5分钟后，将95μl水添加到样本中，然后使用Mosquito HTS液体处理器(TTP Labtech)将200nl变性模板添加到384孔反应板的每个孔中，每个孔中已经包含800nl的WGA混合物(0.58μl Sc反应缓冲液、0.04μl Sc聚合酶(REPLI-g单细胞试剂盒，Qiagen)、0.04μl EvaGreen 20x(Biotium)和0.065μl水)。将板在30℃培养1.5小时，然后在65℃热灭活15分钟。如果需要，在反应中添加EvaGreen允许使用实时PCR机器监测WGA反应。接下来，将250nl的WGA反应物转移到一个新板中，并使用I-DOT分配器(Dispendix)将1.1μl的NEBNext Ultra II FS反应混合物(753nl水、270nl超II FS反应缓冲液和77nl超II FS酶混合物)添加到每个孔中。将板在37℃培养6分钟，然后在65℃热灭活30分钟。接下来将150nl带有列索引的DigiPico索引环状接头添加到所有孔中。请注意，同一列中的所有孔在此阶段都会收到相同的索引寡核苷酸。接下来，使用Mosquito液体处理器将1.2μl超II连接混合物(1150nl超II连接主混合物、38nl连接增强剂和12nl水)添加到每个孔中，混合5个循环，并将板在20℃培养15分钟，然后在65℃热灭活10分钟。接下来，使用Mosquito液体处理器将同一行内的所有孔合并在一起。然后将1.5μl的USER酶(NEB)添加到每行的20μl合并物中，并将反应物在37℃培养15分钟。USER酶在尿嘧啶位置切割环状接头。接下来，使用SPRI珠对产物进行尺寸选择，以达到300-400bp的尺寸范围。然后使用行索引引物将每行的产物扩增4个循环。最终产品合并在一起，并使用SPRI珠纯化最终文库。

ScDigiPico方案

将单个细胞分选到384孔板第一列的孔中。每个孔含有4.5μl的MyPK缓冲液()。将板在55℃培养30分钟。然后向每个孔中加入900nl终止溶液(Stop Solution)，并将板在95℃保持5分钟以灭活蛋白酶K。然后使用Mosquito液体处理器将裂解的细胞跨行分配，每孔200nl。接下来将800nl WGA反应混合物添加到每个孔中，然后类似于DigiPico2进行WGA和文库制备。

索引环状接头-列索引

(文库制备说明手册中的部分序列(

Multiplex Oligos for

(用于

的

多重寡核苷酸)(索引引物组1))

https://international.neb.com/-/media/nebus/files/manuals/ manuale7335.pdf？rev＝4bf1622b342b4d73a2b01443068ed2c5&hash＝B049D91A18CDB471 AB388DC6E67E06B79263E5C5)

P-[索引’]

AGATCGGAAGAGCACACGTCTGAACTUCCCTACACGACGCTCTTCCGATCT

[索引]*T(SEQ ID NO:1)

其中P是5'磷酸基团，*表示硫代磷酸酯键。索引＝列索引(Ci)或行索引(Ri)序列，分别作为每一列或每一行的唯一条形码。

行索引引物(寡核苷酸序列

2007-2013Illumina,Inc.保留所有权利。)

P5:AATGATACGGCGACCACCGAGATCTACAC[r-索引]ACACTCTTTCCCTACACGACGCTCTTCCGATC*T(SEQ ID:NO:2)

P7:CAAGCAGAAGACGGCATACGAGAT[r-索引]GTGACTGGAGTTCAGACGTGTGCTCTTCCGATC*T(SEQ ID:NO:3)

*表示硫代磷酸酯键。

序列表

<110> 牛津大学创新有限公司

<120> 皮克量DNA的全基因组测序方法

<130> JDM104471P.WOP

<150> GB1918043.9

<151> 2019-12-09

<160> 22

<170> PatentIn 3.5版本

<210> 1

<211> 52

<212> DNA

<213> 人工序列

<220>

<223> 索引环状接头

<400> 1

agatcggaag agcacacgtc tgaactuccc tacacgacgc tcttccgatc tt 52

<210> 2

<211> 62

<212> DNA

<213> 人工序列

<220>

<223> 行索引引物

<400> 2

aatgatacgg cgaccaccga gatctacaca cactctttcc ctacacgacg ctcttccgat 60

ct 62

<210> 3

<211> 58

<212> DNA

<213> 人工序列

<220>

<223> 行索引引物

<400> 3

caagcagaag acggcatacg agatgtgact ggagttcaga cgtgtgctct tccgatct 58

<210> 4

<211> 23

<212> DNA

<213> 智人（Homo sapiens）

<400> 4

tgatcgcttc cgagcaataa gaa 23

<210> 5

<211> 23

<212> DNA

<213> 智人

<400> 5

ccttatttct gatgctctta gat 23

<210> 6

<211> 23

<212> DNA

<213> 智人

<400> 6

gtatcagtca gccagaaaaa agg 23

<210> 7

<211> 42

<212> DNA

<213> 智人

<400> 7

tttattgaag tttgttttcc tctttgatcc taccactttt tt 42

<210> 8

<211> 42

<212> DNA

<213> 智人

<400> 8

aagccaatgt attgatcgct tccgagcaat aagaatagtg at 42

<210> 9

<211> 42

<212> DNA

<213> 智人

<400> 9

gttgttgttt gccaagctaa tctgcctggt tttatttata tc 42

<210> 10

<211> 42

<212> DNA

<213> 智人

<400> 10

aagtctacat taaacaatga tcacatctaa agctttatct tt 42

<210> 11

<211> 42

<212> DNA

<213> 智人

<400> 11

ctcatatata aagccttatt tctgatgctc ttagatttct ga 42

<210> 12

<211> 42

<212> DNA

<213> 智人

<400> 12

gggactacag atgtgtgcca tcacacccag ctagtttttt gt 42

<210> 13

<211> 42

<212> DNA

<213> 智人

<400> 13

tgcataggta taggtatcag tcagccagaa aaaaggactt tg 42

<210> 14

<211> 42

<212> DNA

<213> 智人

<400> 14

gtatatatac aaatactttg tccatttaaa aaattaggtt at 42

<210> 15

<211> 42

<212> DNA

<213> 智人

<400> 15

cctatataga ctaacatgga tctaactttt tgactatctt cc 42

<210> 16

<211> 42

<212> DNA

<213> 智人

<400> 16

gaaatgcttt gtgaaatatg tcaacatact ggttgcaaat gc 42

<210> 17

<211> 42

<212> DNA

<213> 智人

<400> 17

gagtatggct atctatacct gccttttaag tttgaaacta ac 42

<210> 18

<211> 42

<212> DNA

<213> 智人

<400> 18

gtctttcctc tctctgtcct tccccgaaag tctactcggg tg 42

<210> 19

<211> 42

<212> DNA

<213> 智人

<400> 19

ggcatgatca ctgcagcctc tctgcttccc agattcaagt ga 42

<210> 20

<211> 42

<212> DNA

<213> 智人

<400> 20

cattaggggc tggacactca tcgagatgac ctgcctacaa at 42

<210> 21

<211> 42

<212> DNA

<213> 智人

<400> 21

gattgaaact gtccatttaa tctccttcct cccattatca at 42

<210> 22

<211> 42

<212> DNA

<213> 智人

<400> 22

tacctattta tctatatatt tcaacttata aaactttctt tc 42

Claims (29)

1.一种单细胞或细胞群的全基因组测序方法，用于识别所述单细胞或细胞群的基因组中的单核苷酸变体、确定所述单细胞或细胞群的基因组中的染色体结构变异或确定所述单细胞或细胞群的基因组中的定相信息，所述方法包括：

i)提供包括多行和多列反应孔的多孔阵列板；

ii)提供单细胞或细胞群的基因组DNA，其中所述基因组DNA分布在所述多孔阵列板上的多个反应孔中，使得每个反应孔的任何给定位点的单链基因组DNA分子不超过一个，

iii)对每个基因组DNA分子进行全基因组扩增(WGA)，以在每个反应孔中提供所述基因组DNA分子的多个拷贝；

iv)将每个反应孔的DNA分子片段化并在每一端连接一对环状接头或使用转座酶递送接头进行标签化以形成适配的DNA片段，其中所述环状接头或转座酶递送接头包括列索引(Ci)序列或行索引(Ri)序列，其中所述Ci序列对于所述多孔阵列板的一列中的每个反应孔的每个环状接头或转座酶递送接头是通用的，或者其中每个Ri序列对于所述多孔阵列板的一行中的每个反应孔的每个环状接头或转座酶递送接头是通用的；

vi)通过对所述适配的DNA片段进行索引PCR来提供索引的DNA文库，其中使用正向索引引物和反向索引引物扩增所述适配的DNA片段以形成索引的PCR产物，其中行索引(Ri)序列或列索引(Ci)序列通过每个正向索引引物和反向索引引物引入到所述适配的DNA片段的每一端，使得所得的索引的PCR产物包括一对对于一列的每个孔都是通用的侧翼列索引(Ci)序列和一对对于一行的每个孔都是通用的侧翼行索引(Ri)序列两者；以及

vii)对所述索引的DNA文库进行测序以提供用于确定所述单细胞或细胞群的基因组中的任何单核苷酸变体、确定所述单细胞或细胞群的基因组中的染色体结构变异或者确定所述单细胞或细胞群的基因组中的定相信息的数据。

2.根据权利要求1所述的单细胞或细胞群的全基因组测序方法，其中，所述细胞或细胞群来自受试者的组织活检。

3.根据权利要求1或2所述的单细胞或细胞群的全基因组测序方法，其中，所述细胞或细胞群包括癌细胞、癌前细胞或疑似癌细胞、或其细胞的组合。

4.根据任一项前述权利要求所述的单细胞或细胞群的全基因组测序方法，其中，所述基因组DNA包括约1至30个细胞的DNA。

5.根据任一项前述权利要求所述的单细胞或细胞群的全基因组测序方法，其中，单细胞的DNA内含物分布在单行的孔中；或者

单细胞或细胞群的DNA内含物分布在单个多孔阵列板的行和列两者的孔中。

6.根据任一项前述权利要求所述的单细胞或细胞群的全基因组测序方法，其中，所述多孔阵列板包括384孔板。

7.根据任一项前述权利要求所述的单细胞或细胞群的全基因组测序方法，其中，在扩增混合物中提供DNA聚合报告分子。

8.根据任一项前述权利要求所述的单细胞或细胞群的全基因组测序方法，其中，提供所述环状接头，使得所述方法包括将每个反应孔的DNA分子片段化的步骤和后续的连接反应，以将环状接头连接到片段化的DNA上；或者

其中，提供所述转座酶递送接头，使得所述方法包括通过标签化过程将所述DNA分子片段化。

9.根据任一项前述权利要求所述的单细胞或细胞群的全基因组测序方法，其中，将每个反应孔的DNA分子片段化为多个dsDNA片段包括通过酶直接片段化。

10.根据任一项前述权利要求所述的单细胞或细胞群的全基因组测序方法，其中，将片段化或标签化试剂同时添加到每个孔中。

11.根据任一项前述权利要求所述的单细胞或细胞群的全基因组测序方法，其中，DNA片段化形成DNA片段后，将所述DNA片段进行末端修复和/或dA加尾，使得它们能够连接到所述环状接头。

12.根据任一项前述权利要求所述的单细胞或细胞群的全基因组测序方法，其中，所述环状接头包括具有二级茎环结构的寡核苷酸，以及，其中所述环状接头包括在环区侧翼的一对互补序列区，其中该对互补序列布置成相互杂交以形成所述环状接头的茎环结构。

13.根据任一项前述权利要求所述的单细胞或细胞群的全基因组测序方法，其中，所述适配的DNA片段的末端是对称的。

14.根据任一项前述权利要求所述的单细胞或细胞群的全基因组测序方法，其中，所述环状接头在环状区域中包含尿嘧啶，并且在所述环状接头连接后，环状DNA的单链区在所述尿嘧啶处被切割。

15.根据任一项前述权利要求所述的单细胞或细胞群的全基因组测序方法，其中，在所述适配的DNA片段中提供Ci序列，并且所述方法能够另外包括在索引PCR之前将一行中每个反应孔的适配的DNA片段合并的步骤；或者

其中，在所述适配的DNA片段中提供Ri序列，并且所述方法另外包括在索引PCR之前将一列中每个反应孔的适配的DNA片段合并的步骤。

16.根据任一项前述权利要求所述的单细胞或细胞群的全基因组测序方法，其中，所述适配的DNA片段包括Ci序列，并且所述正向索引PCR引物和反向索引PCR引物各自包括Ri序列，用于在所述索引的PCR产物中提供一对Ri序列。

17.根据任一项前述权利要求所述的单细胞或细胞群的全基因组测序方法，其中，所述正向索引PCR引物和反向索引PCR引物还包括测序接头序列，使得测序接头结合到所述索引的PCR产物中。

18.根据任一项前述权利要求所述的单细胞或细胞群的全基因组测序方法，其中，对所述索引的DNA文库的索引DNA片段尺寸进行过滤。

19.根据任一项前述权利要求所述的单细胞或细胞群的全基因组测序方法，其中，所述方法包括通过确定源自单孔的基本上所有索引的DNA文库序列是否包含相同的SNV，或者是否只有一部分的索引的DNA文库序列包含相同的SNV，来确定所述单细胞或细胞群的基因组中的任何真实SNV，

其中，在源自单孔的基本上所有索引的DNA文库序列中显示的SNV被确定为所述基因组DNA中的真实SNV，并且在源自单孔的仅一部分的索引的DNA文库序列中发现的SNV被确定为假阳性(FP)SNV。

20.根据任一项前述权利要求所述的单细胞或细胞群的全基因组测序方法，其中，所述方法进一步包括将源自代表所述基因组DNA的一条链的单孔的索引的DNA文库序列与源自代表基因组DNA的互补链的另一孔的索引的DNA文库序列匹配，

其中，基本上在所述基因组DNA的两条互补链的所有索引的DNA文库序列中存在的SNV被确定为真实的SNV，并且基本上在基因组DNA的两条互补链的所有索引的DNA文库序列中都不存在的SNV基本被确定为假阳性。

21.根据权利要求19或20所述的单细胞或细胞群的全基因组测序方法，其中，所述确定是使用从将序列数据映射到参考基因组中生成的BAM文件数据在计算机上进行的。

22.根据权利要求21所述的单细胞或细胞群的全基因组测序方法，其中，通过人工神经网络(ANN)模型，任选地通过多层感知器，执行索引的DNA序列的计算机确定或匹配和/或概率分数的计算。

23.根据任一项前述权利要求所述的单细胞或细胞群的全基因组测序方法，其中，所述方法从肿瘤细胞、疑似肿瘤细胞或癌前细胞以及正常(即非癌)细胞制备索引的DNA文库，并且，其中来自所述肿瘤细胞、疑似肿瘤细胞或癌前细胞的测序数据与作为对照的从取自正常组织的所述正常细胞(即非癌细胞)获得的测序数据进行比较。

24.根据任一项前述权利要求所述的单细胞或细胞群的全基因组测序方法，其中，在计算机中计算特定核苷酸变体为真实SNV或假阳性的概率分数，从而确定给定变体核苷酸具有为真实SNV或假阳性的统计上显著的概率。

25.根据任一项前述权利要求所述的单细胞或细胞群的全基因组测序方法，其中，以成对读取的FastQ文件的形式提供所述序列数据。

26.根据权利要求25所述的单细胞或细胞群的全基因组测序方法，其中，所述成对读取的FastQ文件的序列数据进行修剪以去除接头序列以及为了质量，从而提供修剪后的数据。

27.一种制备用于核酸分子测序的索引的DNA文库的方法，所述方法包括：

i)提供包括多行和多列反应孔的多孔阵列板；

ii)提供核酸分子，其中所述核酸分子分布在所述多孔阵列板上的多个反应孔中，使得每个反应孔的任何给定位点的单链核酸分子不超过一个，

iii)对所述核酸分子进行扩增，以在每个反应孔中提供所述核酸分子的多个DNA拷贝；

iv)将每个反应孔的DNA分子片段化并在每一端连接一对环状接头或使用转座酶递送接头进行标签化以形成适配的DNA片段，其中所述环状接头或转座酶递送接头包括列索引(Ci)序列或行索引(Ri)序列，其中所述Ci序列对于所述多孔阵列板的一列中的每个反应孔的每个环状接头或转座酶递送接头是通用的，或者其中每个Ri序列对于所述多孔阵列板的一行中的每个反应孔的每个环状接头或转座酶递送接头是通用的；

vi)通过对所述适配的DNA片段进行索引PCR来提供所述索引的DNA文库，其中使用正向索引引物和反向索引引物扩增所述适配的DNA片段以形成索引的PCR产物，其中行索引(Ri)序列或列索引(Ci)序列通过每个正向索引引物和反向索引引物引入所述适配的DNA片段的每一端，使得所得的索引的PCR产物包括一对对于一列的每个孔都是通用的侧翼列索引(Ci)序列和一对对于一行的每个孔都是通用的侧翼行索引(Ri)序列两者；以及

任选地，其中所述正向索引引物和反向索引引物进一步将各自的5'和3'测序接头提供到适合用于测序反应的所述索引的PCR产物上。

28.一种制备用于单细胞或细胞群的全基因组测序的索引的DNA文库的方法，用于识别所述单细胞或细胞群的基因组中的单核苷酸变体、确定所述单细胞或细胞群的基因组中的染色体结构变异或确定所述单细胞或细胞群的基因组中的定相信息，所述方法包括：

i)提供包括多行和多列反应孔的多孔阵列板；

ii)提供单细胞或细胞群的基因组DNA，其中所述基因组DNA分布在所述多孔阵列板上的多个反应孔中，使得每个反应孔的任何给定位点的单链基因组DNA分子不超过一个，

iii)对每个基因组DNA分子进行全基因组扩增(WGA)，以在每个反应孔中提供所述基因组DNA分子的多个拷贝；

iv)将每个反应孔的DNA分子片段化并在每一端连接一对环状接头或使用转座酶递送接头进行标签化以形成适配的DNA片段，其中所述环状接头或转座酶递送接头包括列索引(Ci)序列或行索引(Ri)序列，其中所述Ci序列对于所述多孔阵列板的一列中的每个反应孔的每个环状接头或转座酶递送接头是通用的，或者其中每个Ri序列对于所述多孔阵列板的一行中的每个反应孔的每个环状接头或转座酶递送接头是通用的；

vi)通过对所述适配的DNA片段进行索引PCR来提供索引的DNA文库，其中使用正向索引引物和反向索引引物扩增所述适配的DNA片段以形成索引的PCR产物，其中行索引(Ri)序列或列索引(Ci)序列通过每个正向索引引物和反向索引引物引入到所述适配的DNA片段的每一端，使得所得的索引的PCR产物包括一对对于一列的每个孔都是通用的侧翼列索引(Ci)序列和一对对于一行的每个孔都是通用的侧翼行索引(Ri)序列两者；以及

任选地，其中所述正向索引引物和反向索引引物进一步将各自的5'和3'测序接头提供到适合用于测序反应的所述索引的PCR产物上。

29.一种单细胞或细胞群的全基因组测序方法，用于识别所述单细胞或细胞群的基因组中的单核苷酸变体(SNV)、确定所述单细胞或细胞群的基因组中的染色体结构变异或确定所述单细胞或细胞群的基因组中的定相信息，所述方法包括：

i)通过实施根据权利要求27或28所述的方法制备索引的DNA文库，或提供根据权利要求27或28所述的方法制备的索引的DNA文库；

ii)对所述索引的DNA文库进行测序以提供用于确定所述单细胞或细胞群的基因组中的任何单核苷酸变体(SNV)、确定所述单细胞或细胞群的基因组中的染色体结构变异或者确定所述单细胞或细胞群的基因组中的定相信息的数据。

Images