Method for whole genome sequencing of picogram quantities of DNA

Abstract

The present invention relates to methods of whole genome sequencing of a single cell or a cell population for identifying single nucleotide variants in the genome of a single cell or a cell population, determining chromosomal structural variations in the genome of a single cell or a cell population or determining phasing information in the genome of a single cell or a cell population. Also described are methods of preparing indexed DNA libraries for sequencing of nucleic acid molecules; methods of preparing indexed DNA libraries for whole genome sequencing of a single cell or a cell population for identifying single nucleotide variants in the genome of a single cell or a cell population, determining chromosomal structural variations in the genome of a single cell or a cell population or determining phasing information in the genome of a single cell or a cell population; and methods of whole genome sequencing of a single cell or a cell population providing data for identifying single nucleotide variants (SNVs) in the genome of a single cell or a cell population, determining chromosomal structural variations in the genome of a single cell or a cell population or determining phasing information in the genome of a single cell or a cell population.

Description

Technical Field

[0001] This invention relates to methods for preparing indexed DNA libraries for sequencing, such as whole-genome sequencing of single cells or cell populations for identifying single nucleotide variants (SNVs), determining chromosomal structural variations, or determining phasing information in the genome of a single cell or cell group. Background Technology

[0002] Next-generation sequencing has revolutionized our understanding of the genetic evolution of human cells in health and disease. In high-volume cancer genome sequencing, the prevalence of variants—the proportion of cells containing a variant—can be inferred to calculate the clonal composition of a tumor. In turn, understanding the clonal composition allows for the construction of phylogenetic trees to tell the story of how a particular tumor has evolved over time (1-3). Analyzing common mutations within individual clones can be used to infer mutational processes that may have played a role in tumor evolution. Understanding which mutational processes occur within a tumor and what mechanisms drive them is highly desirable, as it provides opportunities for therapeutic interventions or predicting the evolutionary trajectory of a tumor. However, limitations in sequencing depth mean that the use of standard high-volume whole-genome sequencing (WGS) methods ( Figure 1A and Figure 4 Only very common mutations that appear early in tumor development can be detected. Therefore, the ability to model evolutionary events is limited to early events already identified in tumor evolution, rather than recent or current processes (1, 4). This limits the practical application of understanding mutational processes. Studying current or recent evolutionary events requires the confident identification of mutations with very low prevalence (1, 4).

[0003] Sequencing of spatially correlated single cells or small populations of cells offers hope for addressing this problem by detecting cell-specific or clone-specific mutations. Figure 4 This gives a reading of the mutational processes currently observed in cells. Figure 1AHowever, accurate sequencing of small (picometer) amounts of DNA obtained from single cells or spatially relevant cells is extremely challenging. When dealing with small amounts of DNA, the inevitable DNA damage, either through oxidation or spontaneous deamination, is particularly problematic (5). These sources of damage result in a disproportionate number of artificial C>A and C>T mutations, respectively (5–7). Identifying these artificial mutations as variants during variant calling leads to a large number of false positive (FP) variant calls. Therefore, whole-genome amplification can introduce serious errors in single nucleotide variant (SNV) identification, hindering accurate estimation of the mutational load (6). An important concern is that such mutations can also be attributed to biological processes such as the accumulation of oxidative DNA damage with age or the overactivity of members of the APOBEC deaminase family (8–10).

[0004] Previously, it was impossible to distinguish between biologically driven C>A and C>T mutations when processing picogram-sized amounts of DNA and artificial mutations that occurred during library preparation. Furthermore, the usual whole-genome amplification (WGA) step before sequencing increases the number of artificial mutations and amplifies errors caused by DNA damage (5). Several methods have been proposed to reduce DNA damage during library preparation or to filter out false positives during analysis (5, 11–13). However, to date, such techniques still result in the retention of thousands of false-positive mutations, thus requiring extensive validation before definitive biological conclusions can be drawn (5, 11, 12). Since extensive validation is often not feasible (5), a reliable method is needed to eliminate false-positive variants in whole-genome amplification sequencing data.

[0005] Complete Genomics previously published a long-fragment read (LFR) method for whole-genome sequencing and haplotyping from 10 to 20 human cells (Peters BA, et al. Accurate whole-genome sequencing and haplotyping from 10 to 20 human cells. Nature. 2012 Jul 11; 487(7406):190-5. doi:10.1038 / nature11236. PubMed PMID:22785314; PubMed Central PMCID:PMC3397394). However, this method is very complex, prone to bias due to index cross-contamination, and produces a large number of false positives. Summary of the Invention

[0006] Therefore, the object of the present invention is to provide an improved method for preparing DNA libraries for sequencing, SNV analysis, identification of chromosomal structural variations, or identification of phasing information.

[0007] According to a first aspect of the present invention, a method for whole-genome sequencing of a single cell or cell population is provided to identify single nucleotide variants (SNVs) in the genome of a single cell or cell population, determine chromosomal structural variations in the genome of a single cell or cell population, or determine phasing information in the genome of a single cell or cell population, the method comprising:

[0008] i) Provides a porous array plate including multiple rows and columns of reaction wells;

[0009] ii) Provide genomic DNA from a single cell or cell population, wherein the genomic DNA is distributed in multiple reaction wells on a multi-well array plate such that no more than one single-stranded genomic DNA molecule is present at any given site in each reaction well.

[0010] iii) Perform whole-genome amplification (WGA) on each genomic DNA molecule to provide multiple copies of the genomic DNA molecule in each reaction well;

[0011] iv) Fragment the DNA molecule in each reaction well and ligate a pair of circular adapters or label them with transposase delivery adapters at each end to form a suitable DNA fragment, wherein the circular adapters or transposase delivery adapters include a column index (Ci) sequence or a row index (Ri) sequence, wherein the Ci sequence is universal for each circular adapter or transposase delivery adapter for each reaction well in a column of the multi-well array plate, or wherein each Ri sequence is universal for each circular adapter or transposase delivery adapter for each reaction well in a row of the multi-well array plate;

[0012] vi) Providing an indexed DNA library by indexing PCR of a compatible DNA fragment, wherein the compatible DNA fragment is amplified using forward and reverse indexing primers to form an indexed PCR product, wherein a row index (Ri) sequence or a column index (Ci) sequence is introduced by each forward and reverse indexing primer to each end of the compatible DNA fragment, such that the resulting indexed PCR product includes both a pair of universal flanking column index (Ci) sequences per well for a column and a pair of universal flanking row index (Ri) sequences per well for a row; and

[0013] vii) Sequencing the indexed DNA library to provide data for identifying any single nucleotide variant in the genome of a single cell or cell population, identifying chromosomal structural variations in the genome of a single cell or cell population, or identifying phasing information in the genome of a single cell or cell population.

[0014] Advantageously, the present invention provides an indexed DNA library for single DNA molecule sequencing methods to obtain high-quality and data-rich sequencing results from picogram-sized amounts of DNA obtained from clinical samples (referred to as DigiPico; for use in sequencing). Pique DNA number (Sequencing). This invention also provides an advantageous indexing strategy to virtually eliminate cross-contamination and improve ligation efficiency. First, a set of indexes is introduced into the stem-loop of a universal adapter. In the ligation step, all wells in each column of the plate will receive different indexed loop adapters or transposase delivery adapters, thus a total of 24 different oligonucleotides are sufficient to index all columns of the plate with the first set of indexes (column indexing). After the ligation step, all wells in each row can be merged into a single tube, forming 16 different pools. These 16 different pools can be readily purified for use in the next indexing step. In the next step, the purified products in each pool can be indexed by PCR using only the 16 different index primers (row indexing). This enables unprecedented accuracy in single-cell sequencing, a significant improvement over known methods. This invention can be used to identify private mutations and potential neoantigens in single cells or very few cells that could serve as therapeutic targets. This invention can also be used to determine chromosomal structural variations, such as numerical or structural aberrations, or to determine phasing information. Figure 18 It is clearly demonstrated that the method of the present invention can greatly improve the accuracy of identifying true nucleotide variants, eliminating many false positives compared to the Complete Genomics LFR method, and more accurately distinguishing samples presenting different numbers of mutations.

[0015] Cells and cell populations

[0016] The cell or cell population may include eukaryotic cells, such as mammalian cells. In one embodiment, the cell is human. In one embodiment, the cell is at least a diploid cell. The cell may be a cancer cell or a pre-cancerous cell. The cell may include tumor islands. In one embodiment, the cell may be derived from a tissue biopsy of a subject.

[0017] Cells such as tumor islands can be microdissected cells captured by laser.

[0018] In the case of sequencing DNA from multiple cells, these cells may be spatially related. The cells may coexist in a tumor or tumor region. In another embodiment, these cells may or may not be neighbors.

[0019] The SNV to be identified can contain single nucleotide mutations. This method can be used to identify multiple distinct SNVs in genomic DNA.

[0020] Provides nucleic acid molecules and pore distribution

[0021] Nucleic acid can be purified or partially purified. In another embodiment, nucleic acid can be provided in cell lysates. Genomic DNA can be provided as purified DNA. In another embodiment, genomic DNA can be provided from resuspended cell nuclei or whole cells, such as laser-captured microdissection cells.

[0022] Genomic DNA can include the DNA of a single cell or a group of cells (cell population), such as spatially related cells. Genomic DNA can contain the DNA of about 1 to 30 cells. Genomic DNA can contain the DNA of about 1 to 100 cells. In another embodiment, genomic DNA can contain the DNA of about 1 to 80 cells. In another embodiment, genomic DNA can contain the DNA of about 1 to 50 cells. In another embodiment, genomic DNA can contain the DNA of about 1 to 40 cells. In another embodiment, genomic DNA can contain the DNA of about 10 to 30 cells. In another embodiment, genomic DNA can contain the DNA of about 20 to 30 cells. In another embodiment, genomic DNA can contain the DNA of about 20 to 40 cells. In another embodiment, genomic DNA can contain the DNA of about 10 to 40 cells.

[0023] In cases where nucleic acids, such as DNA, are double-stranded, they can be denatured before being dispensed into wells. Denaturation can be achieved by heating and / or using a denaturing buffer. In one embodiment, a denaturing buffer, such as D2 buffer from the Repli-g Single Cell Kit (Qiagen), can be used to denature nucleic acids, such as genomic DNA, or cell nuclei or cells containing genomic DNA.

[0024] Nucleic acids, such as DNA, can be distributed into the wells such that no more than one single-stranded genomic DNA molecule is present at any given site in each reaction well. The distribution of nucleic acids can be facilitated by dilution of the nucleic acids. Therefore, in one embodiment, the nucleic acid solution can be diluted. Those skilled in the art will readily determine the necessary dilution level and solution volume to achieve no more than one single-stranded genomic DNA molecule at any given site in each reaction well. Those skilled in the art will recognize that the necessary dilution level can be determined mathematically such that no more than one single-stranded genomic DNA molecule at any given site in each reaction well is statistically highly probable. For example, when the cell number is known, the Poisson distribution can be used for this calculation.

[0025] In one embodiment, the DNA contents of a single cell can be distributed in a single row or column of wells. Therefore, a multi-well array plate can be used to analyze a variety of different single cells, such as one per row or column. At least one well can be used to add cells and extract DNA contents. In another embodiment, the DNA contents of a cell or cell population are distributed in the wells of rows and columns of a single multi-well array plate.

[0026] Those skilled in the art will recognize that any standard multiwell plate can be used in the methods of this invention. Preferably, the multiwell plate is compatible with any PCR and / or sequencing instrument that can be used. The multiwell plate may include a 384-well plate, such as a 24x16-well plate. In another embodiment, the multiwell plate may include a 1536-well plate. Those skilled in the art will understand that a greater number of Ri and / or Ci sequences may be required to index a larger plate.

[0027] Using a 384-well multi-well array plate can advantageously provide enough wells to distribute a diluted genomic DNA strand of about 20-30 cells, thus providing a single DNA molecule for each well.

[0028] Amplification

[0029] In embodiments where the nucleic acid is genomic DNA, the amplification of the genomic DNA molecule may include whole genome amplification (WGA). WGA may include the step of adding amplification reagents for DNA amplification to the genomic DNA. The amplification reagents for DNA amplification may also be referred to as an "amplification mixture." Those skilled in the art will understand that an amplification mixture may contain all the reagents required to amplify DNA (i.e., to produce multiple copies of DNA). Such components may include reaction buffers, polymerases, and dNTPs. DNA polymerization reporter molecules, such as DNA-binding dyes (e.g., Evagreen), may also be used. TMThis can be provided in the amplification mixture, for example, to allow for monitoring of the amplification reaction using real-time PCR. The DNA-binding dye can consist of two monomeric DNA-binding dyes linked by a flexible spacer. In the absence of DNA, the dimer dye can present a circular conformation that is inactive in DNA binding. When DNA is available, the circular conformation can be equilibriumly converted to a random conformation that allows it to bind to DNA and emit fluorescence.

[0030] Amplification reagents can be provided in each well before or after adding DNA to the well.

[0031] Those skilled in the art will be able to provide suitable conditions for the amplification reaction to occur, including appropriate temperature and incubation time. For example, the plate can be incubated at about 30°C for at least about 1 hour, followed by heat inactivation, for example, at about 65°C for at least 5 minutes.

[0032] Fragmentation and ligation of circular or transposase delivery adapters

[0033] In one embodiment, a circular adapter is provided, such that the method includes the step of fragmenting the DNA molecule in each reaction well and a subsequent ligation reaction to ligate the circular adapter to the fragmented DNA. The fragmented DNA may undergo end repair prior to ligation. In an alternative embodiment, a transposase delivery adapter may be provided, such that the method includes fragmenting the DNA molecule via a tagging process. Tagging may include providing a transposase carrying an oligonucleotide, such as Tn5, referred herein as a transposase delivery adapter. Conventional techniques and reagents for performing tagging to form a suitable DNA molecule will be familiar to those skilled in the art.

[0034] Fragmenting the DNA molecule from each reaction well into multiple dsDNA fragments can include direct fragmentation, such as enzymatic or mechanical fragmentation. In one embodiment, DNA fragmentation includes enzymatic fragmentation.

[0035] DNA fragmentation or tagging can be provided by adding fragmentation or tagging reagents to the DNA in each well. This can be achieved, for example, by using a multi-well dispenser, such as an I-DOT (Dispendix, Germany) dispenser or similar, to simultaneously add the fragmentation or tagging reagents to each well. The fragmentation or tagging reaction can be timed to provide fragments of the desired size. Those skilled in the art will understand that the timing of the fragmentation or tagging reaction can depend on the method used, such as the type and level of enzyme provided for the reaction. Therefore, those skilled in the art can follow standard protocols for timing a given reaction, such as the timing of a reaction kit.

[0036] Fragmentation reagents may include restriction enzymes or nicking enzymes, such as DNase I (deoxyribonuclease I). When a nicking enzyme is provided, a single-strand-specific enzyme that recognizes the nick site and then cleaves the second strand can be provided. In one embodiment, a library preparation kit, such as the Lotus DNA Library Preparation Kit (IDT, USA), can be used.

[0037] After DNA fragmentation to form dsDNA fragments, end repair and dA-tailing can be performed on these fragments so that they can be ligated to other DNA molecules, such as circular linkers. Enzymes used for end repair and / or dA-tailing can include DNA polymerases, such as T4 DNA polymerase, and polynucleotide kinases (PNKs), such as T4 polynucleotide kinase. T4 DNA polymerase (in the presence of dNTPs) can fill in the 5' overhang and trim the 3' overhang down to the dsDNA interface to generate a blunt end. T4 PNK can then phosphorylate the 5' terminal nucleotide. DNA polymerases with terminal transferase activity that leave the 3' terminal adenine, such as Taq DNA polymerase, can be provided for A-tailing.

[0038] In one implementation, dsDNA fragmentation, end repair, and dA tailing are all performed in a single reaction.

[0039] In one implementation, the circular adapter can be introduced onto fragmented DNA via ligation. Ligation of the circular adapter to the dsDNA fragment may include adding the circular adapter and a ligase, such as T4 DNA ligase.

[0040] Circular adapters may contain oligonucleotides, such as DNA, with a secondary stem-loop structure. The stem-loop structure can be provided by a single oligonucleotide molecule containing a pair of complementary sequence regions flanking the loop region, wherein the pair of complementary sequences are aligned to hybridize to form the stem-loop structure of the circular adapter. The circular adapter further encodes a column index (Ci) sequence or a row index (Ri) sequence in the stem region.

[0041] The column index (Ci) sequence or row index (Ri) sequence can contain predetermined sequences that can label DNA as originating from a row or column, respectively. The length of the column index (Ci) sequence or row index (Ri) sequence can be at least three nucleotides.

[0042] In one implementation, the ends of the adapted DNA fragments can be symmetrical. Specifically, the circular adapters or transposase delivery adapters attached to each end of the dsDNA fragments are identical, such that each dsDNA fragment receives a pair of identical flanking circular adapters or transposase delivery adapters. The Ci sequence pairs on the same adapted DNA fragments can be identical. Alternatively, if Ri sequences are provided, the pair of Ri sequences on the same adapted DNA fragments may be identical.

[0043] Advantageously, providing two identical Ci or Ri sequences on the adapted DNA fragment provides a marker to avoid analyzing indexed DNA library sequences that may be contaminated by cross-contamination between different columns or rows. Specifically, any indexed DNA library sequences that do not match a Ci sequence at either end can be discarded from the data analysis. Alternatively, if a Ri sequence is provided, any indexed DNA library sequences that do not match a Ri sequence at either end can be discarded from the analysis. This provides a first level of redundancy to eliminate index cross-contamination in subsequent data, which is an important problem in index library preparation and its subsequent analysis.

[0044] Circular adapters can provide 3' or 5' overhangs to aid in ligation to dsDNA fragments. A 3' or 5' overhang can be provided when the stem regions of the circular adapter hybridize together (i.e., the circular adapter is in a secondary / stem-loop structure). The 3' or 5' overhang can correspond to complementary overhangs on the dsDNA fragment, which has been end-repaired and prepared for ligation. The overhang may contain a single thymine.

[0045] The ring connector sequence may include the sequence of SEQ ID NO: 1 or a functional variant thereof.

[0046] After attaching the circular adapter, the single-stranded region of the circular DNA can be cleaved. This single-stranded region of the circular DNA can be enzymatically cleaved, for example by a USER (Uracil-Specific Excision Reagent) enzyme, which creates a single nucleotide gap at the location of uracil present in the loop. Therefore, in one embodiment, the circular adapter can contain uracil within the circular region.

[0047] Merge one row of holes

[0048] If the Ci sequence is provided in the adapted DNA fragment, the method may additionally include a step of merging the adapted DNA fragments from each reaction well into a single row prior to indexing PCR. Alternatively, if the Ri sequence is provided in the adapted DNA fragment, the method may additionally include a step of merging the adapted DNA fragments from each reaction well into a single column prior to indexing PCR. The merged adapted DNA fragments can then be used for PCR indexing of each row or column in a single pooled reaction, depending on what is being merged. In an alternative embodiment, columns or rows are not merged prior to indexing PCR.

[0049] Advantageously, merging rows or columns before indexing PCR significantly improves the efficiency of library preparation. For example, for a 16x24 (384) well plate, if 16 rows are merged between the introduction of circular adapters or transposase delivery adapters and the indexing PCR step, only 16 separate indexing PCR reactions are needed instead of 384 if they are not merged.

[0050] Size selection and indexing PCR

[0051] Prior to indexing PCR, a suitable DNA fragment can be selected based on size, and self-ligating adapters can be removed from the reaction, if applicable. An example of a desired size could be approximately 300-400 bp in length. Size selection can be provided by isolating or purifying a suitable DNA fragment of the desired length, for example, using gels or beads. SPRI beads (Solid Phase Reversible Immobilization beads) can be used for size selection. SPRI beads can comprise magnetic particles coated with carboxyl groups (in the form of succinic acid), which can nonspecifically and reversibly bind DNA.

[0052] Indexed PCR may include the step of mixing a fitted DNA fragment with a set of forward and reverse indexed PCR primers and PCR reagents. The forward and reverse indexed PCR primers may contain sequences arranged to hybridize with the fitted DNA fragment sequence to initiate polymerization. The sequences arranged to hybridize with the fitted DNA fragment sequence to initiate polymerization may be complementary sequences. The sequences used to initiate polymerization from the forward and reverse indexed PCR primers may be provided by circular adapters or transposase delivery adapters. The sequences used to initiate polymerization from the forward and reverse indexed PCR primers may be flanked by the Ci or Ri sequence of the fitted DNA fragment, such that the Ci or Ri sequence binds to the indexed PCR product.

[0053] The complementary sequences for hybridization provided by the forward and reverse primers can each be between approximately 15 and 30 nucleotides in length, for example, approximately 26 nucleotides in length.

[0054] In embodiments where the adapted DNA fragments contain the Ci sequence, the forward and reverse indexing PCR primers may each include the Ri sequence to provide a pair of Ri sequences in the indexed PCR product. In the case of merged rows, the Ri sequence is added to each adapted DNA fragment in the pool (from all wells in one row). Alternatively, in the case of unmerged rows, the same Ri sequence can be provided for each well in one row.

[0055] In an alternative implementation, the adapted DNA fragment contains the Ri sequence, and the forward and reverse indexing PCR primers may each contain a Ci sequence to provide a pair of Ci sequences in the indexed PCR product. In the case of merged columns, the Ci sequence is added to each adapted DNA fragment in the pool (from all wells in one row). Alternatively, in the case of unmerged rows, the same Ci sequence can be provided for each well in a column.

[0056] The row index (Ri) sequence, which can be provided by the forward and reverse primers, can be the same for each adapted DNA fragment in or from a row. Alternatively, the column index (Ci) sequence, which can be provided by the forward and reverse primers, can be the same for each adapted DNA fragment in or from a column.

[0057] The row index (Ri) sequence or column index (Ci) sequence provided by the forward and reverse primers can each be at least 3 nucleotides long, for example, about 8 nucleotides long.

[0058] The resulting ends of the indexed PCR product can be symmetrical. For example, the flanking sequences of the original DNA fragment sequence can be symmetrical. The indexed PCR product can contain a DNA fragment sequence flanked by a pair of identical Ci sequences (i.e., inner flanking sequences) and also flanked by a pair of identical Ri sequences (i.e., outer flanking sequences). In an alternative embodiment, the indexed PCR product can contain a DNA fragment sequence flanked by a pair of identical Ri sequences (i.e., inner flanking sequences) and also flanked by a pair of identical Ci sequences (i.e., outer flanking sequences).

[0059] Advantageously, providing two identical Ci or Ri sequence pairs on the index DNA fragment, in addition to the previously provided Ri or Ci sequences (provided by circular adapters or transposase delivery adapters), also provides a marker to avoid analyzing index DNA library sequences that may be contaminated by cross-contamination between different columns or rows. Specifically, any index DNA library sequences that do not have a matching Ci sequence at each end can be discarded from the data analysis. Alternatively, if a Ri sequence is provided, any index DNA library sequences that do not have a matching Ri sequence at each end can be discarded from the analysis. Providing matching Ci and Ri sequence pairs on the index DNA fragment provides both first- and second-level redundancy to eliminate index cross-contamination in subsequent data, which is an important issue in index library preparation and its subsequent analysis.

[0060] Forward and reverse indexed PCR primers can also contain sequencing adapter sequences, allowing the sequencing adapter to bind to the indexed PCR product. The sequencing adapter sequence on the primer can be 5'.

[0061] Sequencing adapters can be the ends of indexed PCR products. When sequencing adapter sequences are provided, the resulting ends of the indexed PCR product may not be symmetrical. For example, one end of the indexed PCR product can be modified with a sequencing adapter that has a different sequencer at the other end. Those skilled in the art will understand the sequencing adapters required for a given sequencing technology. For example, in the case of dye sequencing (e.g., Illumina dye sequencing), the sequencing adapter can be a P5 and P7 sequencing adapter (i.e., P5 at one end of the indexed PCR product, and P7 at the other end). Index primers providing the P5 sequence can contain the sequence of SEQ ID NO: 2. Index primers providing the P7 sequence can contain the sequence of SEQ ID NO: 3.

[0062] Once formed, the indexed PCR product can be called an "indexed DNA library sequence" or an "indexed DNA fragment". The combined indexed PCR product, indexed DNA sequence, or indexed DNA fragment can be called an "indexed DNA library".

[0063] Indexed DNA library

[0064] Indexed DNA libraries can be filtered by the size of the indexed DNA fragments, ensuring that only DNA fragments with the desired or appropriate index length are available for sequencing. After indexing PCR, the indexed PCR fragments can be purified / isolated, for example, using beads (e.g., SPRI beads). Purification removes unwanted short fragments, primer dimers, or other PCR artifacts or reagents.

[0065] You can check if the size distribution of the indexed library is appropriate. We typically check the library size distribution, such as on Tapestry or a bioanalytical instrument (Agilent) or similar. The size of the indexed DNA library can be adjusted by dilution after sequencing preparation, for example, to approximately 4 nM.

[0066] Indexed DNA libraries can be saved for later use, such as sequencing. For example, DNA libraries can be frozen or refrigerated.

[0067] Sequencing of indexed DNA libraries

[0068] The indexed DNA library can be sequenced or adapted for sequencing. Sequencing can be next-generation sequencing (NGS). Sequencing can be dye sequencing (e.g., Illumina dye sequencing), nanopore sequencing, or ion-fluidic sequencing. Those skilled in the art will be familiar with the various sequencing technologies / methods available and the required sequencing adapters.

[0069] Sequencing can be multiplex sequencing, in which DNA libraries with multiple indexes are sequenced simultaneously.

[0070] Determination and data analysis of mutations / nucleotide variations

[0071] This method may include identifying any true SNVs in the genome of a single cell or cell population by determining whether DNA library sequences from substantially all indices originating from a single well contain the same SNV, or whether only a portion of the indices of the DNA library sequences contain the same SNV. SNVs appearing in DNA library sequences from substantially all indices originating from a single well can be identified as true SNVs in the genomic DNA. Alternatively, SNVs found only in a portion of the indices of the DNA library sequences originating from a single well can be identified as false positive (FP) SNVs. False positive SNVs may be errors caused by damage or replication errors.

[0072] This method may further include pairing a DNA library sequence indexed from a single well representing one strand of genomic DNA with a DNA library sequence indexed from another well representing the complementary strand of genomic DNA. SNVs that are substantially present in the DNA library sequences indexed on both complementary strands of the genomic DNA may be identified as true SNVs. SNVs that are substantially absent in the DNA library sequences indexed on both complementary strands of the genomic DNA may be identified as false positives (i.e., not true SNVs).

[0073] The step of determining whether DNA library sequences from virtually all indices of a single well contain the same SNV or whether only a subset of indices of the DNA library sequences contain the same SNV can be performed in a computer (in silico), for example, using BAM file data. Alternatively, the step of pairing DNA library sequences from indices of a single well representing one strand of genomic DNA with DNA library sequences from indices of another well representing the complementary strand of genomic DNA can be performed in a biological computer, for example, using BAM file data.

[0074] In one embodiment, sequencing data from tumor cells, suspected tumor cells, or precancerous cells can be compared with sequencing data obtained from normal cells (i.e., non-cancer cells) taken from normal tissue (i.e., non-cancerous tissue) (e.g., as a control). Therefore, in one embodiment, the method includes preparing indexed DNA libraries from tumor cells, suspected tumor cells, or precancerous cells, and normal (i.e., non-cancer) cells. Indexed DNA libraries can be prepared in parallel for each cell type, for example, in different wells of the same multi-well plate, or separately. Sequencing of indexed DNA libraries from different cell types can be performed in the same sequencing run. Different types of cells (e.g., cancer cells or normal cells) can be from the same subject.

[0075] The probability fraction of a specific nucleotide variant being a true SNV or a false positive can be calculated in a computer, thereby determining whether a given variant nucleotide has a statistically significant probability of being a true SNV or a false positive.

[0076] In one implementation, sequencing a DNA library to determine SNVs in the library includes generating multiplex sequencing data from multiple wells and analyzing the SNV data.

[0077] In one implementation, analyzing SNV data includes de-multiplexing sequencing data, such that data from each well is assigned to a single well group. Furthermore, DNA libraries with different indices can be sequenced within the same sequencing run; therefore, the method can further include demultiplexing the sequencing data to identify / group DNA libraries with different indices.

[0078] The provided sequence data can be in the form of a paired-read FastQ file. Sequence data can be trimmed to remove adapter sequences, for example, in a paired-read FastQ file. Sequence data can also be trimmed for quality reasons, for example, in a paired-read FastQ file. Those skilled in the art will be able to easily adjust the desired threshold for the quality fraction of each base read in the sequence, for example, using a program such as TrimGalore. The resulting data can be referred to as "trimmed data."

[0079] In one implementation, analyzing SNV data involves mapping sequence data to a reference genome, such as the human hg19 reference genome, to generate aligned sequencing data in the form of a sequence alignment map (SAM) or a binary file version thereof (e.g., a BAM file). Trimmed read data can be mapped to the reference genome. SAM or BAM file data can be used to determine the SNVs present in each well.

[0080] Mapping can be done using a program such as Bowtie2, where the ignore-quals parameter is activated and duplicate readings are flagged, for example using the Picard tool.

[0081] Joint variant identification can be performed on all individual BAM files as well as merged BAM files from all holes, for example using a variant identification program such as the Platypus variant identification program (caller).

[0082] Low-quality (i.e., low-confidence) variants can be filtered out from the data. For example, low-quality (i.e., low-confidence) variants can be removed from the data by applying a quality filter. Example quality filters in the Platypus identification program may include QUAL>60, FR>0.1, HP≤4, QD>10, and SbPval≤0.95. Those skilled in the art will recognize that filtering out low-confidence variants is a routine procedure, and that each variant identification program may assign different confidence scores to each variant based on its algorithm, and that these confidence scores can be used to filter low-confidence (quality) variants. Therefore, the specific parameters depend on the variant identification program used.

[0083] The total number of wells covering each site (Tw) and the number of wells supporting each variant (Vw) can be determined. Well count filters, such as Tw>5, Vw>2, and Vw / Tw>0.1, can be used to retain only high-confidence sites for analysis.

[0084] Genomic regions with poor mappability (i.e., known regions more likely to be affected by biased readings) can be removed from the analysis, for example, using VCFtools.

[0085] Variant identification (genotyping) can then be performed on WGS data (e.g., from blood and bulk tumors) using a list of high-confidence variants already identified in the data, for example, using Platypus. The PlatypusminPosterior parameter can be set to 0, and the minMapQual parameter can be set to 5. Any variant that is completely unsupported in standard WGS data can be extracted as a UTD (unique to DigiPico) variant. Any variant that is also fully present in bulk sequencing data of blood samples (based on GATK analysis) can be extracted as a TP (true positive) variant.

[0086] Using Artificial Neural Networks (ANN)

[0087] The computer can determine or pair indexed DNA sequences and / or calculate probability scores according to the methods and calculations described herein. In one embodiment, the computer determination or pairing of indexed DNA sequences and / or the calculation of probability scores can be performed using an artificial neural network (ANN) model, such as a multilayer perceptron.

[0088] A multilayer perceptron may have an input layer consisting of N neurons (e.g., N = 41), where N is the number of features used in each experiment. An ANN model may include at least two hidden layers with ReLU (Rectified Linear Unit) activation. The final layer of the ANN may be a single output neuron with sigmoid activation. The loss function may be binary cross-entropy.

[0089] For example, ANNs can be programmed in Python 3 using Keras. Those skilled in the art will recognize that Keras is a free and open-source Python library for developing and evaluating deep learning models. However, other libraries can be used.

[0090] An ANN can be pre-trained using one or more datasets. For example, an ANN can be trained using a dataset containing known nucleotide variants.

[0091] Other aspects

[0092] According to another aspect of the present invention, a method for preparing a DNA library indexed for nucleic acid molecular sequencing is provided, the method comprising:

[0093] i) Provides a porous array plate including multiple rows and columns of reaction wells;

[0094] ii) Provide nucleic acid molecules distributed in multiple reaction wells of a multi-well array plate such that no more than one single-stranded nucleic acid molecule is present at any given site in each reaction well.

[0095] iii) Amplify the nucleic acid molecules to provide multiple DNA copies of the nucleic acid molecules in each reaction well;

[0096] iv) Fragment the DNA molecule in each reaction well and ligate a pair of circular adapters or tag it with a transposase delivery adapter to form a suitable DNA fragment, wherein the circular adapter or transposase delivery adapter includes a column index (Ci) sequence or a row index (Ri) sequence, wherein the Ci sequence is universal for each circular adapter or transposase delivery adapter for each reaction well in a column of the multi-well array plate, or wherein each Ri sequence is universal for each circular adapter or transposase delivery adapter for each reaction well in a row of the multi-well array plate;

[0097] vi) Providing an indexed DNA library by indexing PCR of adapted DNA fragments, wherein adapted DNA fragments are amplified using forward and reverse indexing primers to form indexed PCR products, wherein row index (Ri) sequences or column index (Ci) sequences are introduced into each end of the adapted DNA fragment through each forward and reverse indexing primer, such that the resulting indexed PCR product includes both a pair of universal flanking column index (Ci) sequences per well for a column and a pair of universal flanking row index (Ri) sequences per well for a row; and

[0098] Optionally, the forward and reverse index primers further provide their respective 5' and 3' sequencing adapters to the PCR product of the index, which is suitable for use in the sequencing reaction.

[0099] Nucleic acids can be DNA or RNA. In one embodiment, the nucleic acid is genomic DNA. In another embodiment, the nucleic acid can be mRNA.

[0100] According to another aspect of the present invention, a method is provided for preparing a DNA library indexed for whole-genome sequencing of a single cell or cell population, to identify single nucleotide variants in the genome of a single cell or cell population, to determine chromosomal structural variations in the genome of a single cell or cell population, or to determine phasing information in the genome of a single cell or cell population, the method comprising:

[0101] i) Provides a porous array plate including multiple rows and columns of reaction wells;

[0102] ii) Provide genomic DNA from a single cell or cell population, wherein the genomic DNA is distributed in multiple reaction wells on a multi-well array plate such that no more than one single-stranded genomic DNA molecule is present at any given site in each reaction well.

[0103] iii) Perform whole-genome amplification (WGA) on each genomic DNA molecule to provide multiple copies of the genomic DNA molecule in each reaction well;

[0104] iv) Fragment the DNA molecule in each reaction well and ligate a pair of circular adapters or tag it with a transposase delivery adapter to form a suitable DNA fragment, wherein the circular adapter or transposase delivery adapter includes a column index (Ci) sequence or a row index (Ri) sequence, wherein the Ci sequence is universal for each circular adapter or transposase delivery adapter for each reaction well in a column of the multi-well array plate, or wherein each Ri sequence is universal for each circular adapter or transposase delivery adapter for each reaction well in a row of the multi-well array plate;

[0105] vi) Providing an indexed DNA library by indexing PCR of adapted DNA fragments, wherein adapted DNA fragments are amplified using forward and reverse indexing primers to form indexed PCR products, wherein row index (Ri) sequences or column index (Ci) sequences are introduced into each end of the adapted DNA fragment through each forward and reverse indexing primer, such that the resulting indexed PCR product includes both a pair of universal flanking column index (Ci) sequences per well for a column and a pair of universal flanking row index (Ri) sequences per well for a row; and

[0106] Optionally, the forward and reverse index primers further provide their respective 5' and 3' sequencing adapters to the PCR product of the index, which is suitable for use in the sequencing reaction.

[0107] The indexed nucleic acids can be sequenced, for example as described in this article.

[0108] Therefore, according to another aspect of the present invention, a method for whole-genome sequencing of a single cell or cell population is provided to provide data for identifying single nucleotide variants (SNVs) in the genome of a single cell or cell population, the method comprising:

[0109] i) Prepare an indexed DNA library by implementing the method according to the invention, or provide an indexed DNA library prepared by the method according to the invention;

[0110] ii) Sequencing the indexed DNA library to provide data for identifying any single nucleotide variants (SNVs) in the genome of a single cell or cell population.

[0111] Sequencing data can be used to identify SNVs, as described herein, for example. Alternatively or concurrently, sequencing data can be used to identify genetic changes associated with chromosomal structural variations. Chromosomal aberrations can include numerical and / or structural aberrations.

[0112] Alternatively, sequencing data can be used to determine phasing information in cells or cell populations.

[0113] As disclosed in the specification and / or drawings, the invention may also include one or more features, individually or in combination.

[0114] definition

[0115] The term "spatially related cells" is understood to refer to cells that are directly adjacent to each other.

[0116] The term “false positive (FP) mutation” or “false positive (FP) SNV” is understood to refer to a variant nucleotide that is not present in the genome before DNA is extracted from an intact cell. For example, a false mutation could be an error caused by damage or a replication error.

[0117] The terms “true mutation / SNV” or “true positive mutation / SNV” are used interchangeably and are understood to refer to variant nucleotides present in the genomic DNA of living cells prior to DNA extraction.

[0118] The term "single nucleotide variant" (SNV) can include single nucleotide polymorphisms (SNPs) or any other variation in a sequence, such as a mutation. A mutation or variation can include a nucleotide substitution, addition, or deletion in a given sequence.

[0119] Chromosomal aberrations are understood as missing, extra, or irregular portions of chromosomal DNA. They can originate from a typical number of abnormalities in the chromosomes or their structures within one or more chromosomes. These include various aberrations such as deletions, duplications, and insertions. Balanced aberrations can occur, such as inversions and interchromosomal and intrachromosomal translocations. In addition, insertions of mobile elements, segment duplications, and multiple allele aberrations can occur. The final multiple combinations of the above can produce complex rearrangements.

[0120] "Phasing" is understood as the task or process of assigning alleles (A, C, T, and G) to paternal and maternal chromosomes. Phasing helps determine whether a match occurs on the paternal or maternal side, on both sides, or on neither side. Phasing also aids in the process of chromosome mapping—assigning segments to specific ancestors.

[0121] Those skilled in the art will understand that, where appropriate, optional features of one embodiment or aspect of the invention may be applied to other embodiments or aspects of the invention. Attached Figure Description

[0122] Embodiments of the invention will now be described in more detail by way of example only, with reference to the accompanying drawings.

[0123] Figure 1. DigiPico sequencing principles, workflow, and performance. (A) The WGS method can only identify early mutational processes (EM) in dominant amplified clones within tumors (red and blue). Currently active mutational processes (CM) result in a wide variety of subclones with different clone-specific mutations. This diversity determines the evolutionary trajectory of the tumor. (B) Template partitioning is performed before WGA so that each compartment receives no more than one DNA molecule from each site to identify artificial mutations. Artificial mutations result in biallelic compartments because errors caused by damage (red) and replication errors (blue-green) occur randomly during replication. Note that true mutations are always present in all product DNA strands within the same compartment. (C) DigiPico sequencing workflow. LCM: Laser Capture Microcutting. (D) End-point relative fluorescence units (RFUs) from EvaGreen-tagged DNA are used to ensure uniform distribution of the template and WGA process across the entire plate. RFU values ​​are normalized to achieve a median of 1 in each run. (E) The relative amount of adapter ligation products in each well was measured by qPCR using Illumina adapter primers (P5 and P7). Ct values ​​were normalized to reach a median of 0 in each run. (F) Simplifying the Digipico library preparation process requires a miniaturized WGA that can specifically and sensitively amplify subpicogram amounts of DNA in each well. The values ​​represent the average RFU value of 9 replicates. Error bars represent SD. (G, H, and I) Preliminary analysis of the DigiPico sequencing data from each well in each run, as shown in the figure, confirms the high quality of sequencing and the uniformity of mapping rate, coverage depth, and coverage width. (J) Definitions unique to DigiPico (UTD) variants. Subtracting SNVs identifiable in standard WGS data from the corresponding DigiPico data results in UTD variants. These will consist primarily of artificial mutations as well as some cloning-specific mutations. Since the template in run D1110 is actually a subset of the template used in standard WGS, all true variants in DigiPico run D1110 are expected to also appear in the standard WGS data. In contrast, due to depth limitations, clonal-specific variants may not be present in standard WGS data running D1111, even though DNA molecules supporting such variants may exist at very low frequencies in bulk DNA samples. In all boxplots, the horizontal line represents the median. The boxes represent the interquartile range (between the 25th and 75th percentiles). Whiskers represent the range excluding outliers. Outliers are defined as data points that are 1.5 times higher or lower than the interquartile range.

[0124] Figure 2. MutLX Algorithm, Design, and Results. (A) Comparison of the number of wells supporting various mutation types in run D1110 confirms that, as assumed, most UTDs exist only in a minority of wells. The horizontal line represents the median. The boxes represent interquartile ranges. Whiskers show the range excluding outliers, which are defined as being more than 1.5 times the interquartile range. (B) Similarly, the biallelic interval rate of UTDs appears to be significantly higher compared to true variants. This value is calculated by dividing the number of wells with both variant and reference alleles by the total number of wells with evidence of variant alleles. (C) A graph showing the main challenges of analyzing DigiPico data using ANN. Each circle / asterisk represents a variant. The red line shows the behavior of the classification model. All variants above and / or to the left of the model prediction line are true variants. Analyzing samples without clonal-specific variants will result in accurate separation between true and artificial mutations. Conversely, analyzing samples with true clonal-specific mutations will result in a suboptimal model, which may lead to overfitting to true UTDs. This forces the model to remove all FP identifications at the cost of losing almost all clone-specific variants. (D) shows a graph of the two-step training process in MutLX. The first training step identifies some mislabeled true mutations in the UTD (grey circles). All potentially mislabeled data points (black) are temporarily removed from the analysis in the second training step to obtain a better model that assigns probability scores to all mutations. Finally, the probability scores obtained from the model are combined with uncertainty estimates of these probability scores (as described in E) to effectively eliminate FP identifications while maintaining excellent sensitivity to true clone-specific variants. (E) shows a graph of the test-time drop-out analysis used to compute uncertainty estimates of the probability scores. Black neurons represent neurons that were turned off during the drop-out analysis. Accepting only variants with high probability scores and low uncertainty scores should allow for the elimination of FP variant identifications. (F) shows the ROC curves of the MutLX analysis outputs running D1110, D1111, DE011, and GM12885. Circles represent the default cutoff values ​​determined by MutLX. (G) represents a bar chart showing the number of UTDs passed through the outputs of SCcaller, Platypus, and MutLX. Since no true UTDs are expected in runs D1110, DE011, and GM12885, the number of UTDs in these runs represents the FP rate for each analysis method. The Platypus value is based on DigiPico-specific filtering criteria prior to applying MutLX.

[0125] Figure 3. Identification of active mutational processes using DigiPico / MutLX. (A) Schematic diagram of tumor evolution in HGSOC patient #11152. Standard bulk WGS of various tumor samples identified ~11,000 common somatic mutations at all sites. The purple dashed lines indicate the points where the most recent common ancestor of the tumor samples studied bifurcated. Bulk sequencing also identified nearly 5,000, 3,000, and 2,000 subclonal mutations, respectively, which were specific for prechemoradiotherapy omental masses, PT2R recurrence, and PALNR recurrence. However, these mutations can occur at any time during these clonal expansions and are biased towards older mutations. This is due to the limitations of identifying low-prevalence somatic mutations. However, DigiPico sequencing of five prechemoradiotherapy tumor islands, PT2R, and PALNR recurrence sites identified varying numbers of recently occurring clonal-specific mutations (indicated by red numbers) in each of these samples. The significantly higher number of clonal-specific variants in PT2R indicates the presence of an active mutational process. (B) The strong clonal-specific kataegis event on chromosome 17 during the D1111 run highlights this active mutagenesis process. The Y-axis represents the pairwise distance of continuous somatic mutations on a logarithmic scale. Only mutations from chromosome 17 are shown. Mutations associated with the subclonal kataegis event are highlighted in boxes, and almost all of these mutations are in the form of strand-specific C>T or C>G mutations. This indicates that the APOBEC enzyme is involved in this hypermutation process. (C) Representative examples of some mutations involved in kataegis. The presence of all mutations on the forward strand of the genome further confirms the involvement of the hypermutation event. Figure 13 ).

[0126] Figure 4 The challenge of identifying recent mutations. While older mutations can be readily studied from large amounts of tumor sequencing data, the study of recent mutations from these data is hampered by the low variant allele fraction (VAF) of the mutations involved. Therefore, heuristic filtering criteria are insufficient to identify recent mutations. Reliable studies of recent mutations require examination of individual cancer cells or tumor islands isolated by laser capture microdissection (LCM) (Figure 1). However, WGA with limited template quantity in such samples leads to a large number of false-positive variant identifications, thus hindering the identification of island-specific variants. Our analysis pipeline, MutLX, overcomes this problem by eliminating FP variant identification from DigiPico sequencing data.

[0127] Figure 5.Analysis workflow for DigiPico data. (1) Next-generation sequencing reads from normal tissues, bulk tumors, and DigiPico libraries are first mapped to the human genome to generate bam files. DigiPico reads are split into 384 FastQ files, one FastQ file for each well of a 384-well plate. (2) The 384 individual bam files from DigiPico are merged into one bam file. (3) De novo joint variant identification is performed on the 384 individual bam files and the merged bam file using the Platypus variant identification program. The merged bam file is added to ensure that variants with low coverage per well are not missed during variant identification. (4) The resulting re-DigiPico variants are then used as a reference for variant re-calling in standard WGS data from normal tissues and bulk tumors. (5) Variant re-identification data can then be used to extract variants specific to DigiPico (UTD) by eliminating any variants with supporting reads in the standard WGS data. (6) Standard WGS data was also used for variant identification using GATK to obtain a high-confidence list of phylogenetic SNPs. (7) This list will serve as a guide for extracting TP variant identification from the DigiPico data. For this purpose, any variant identified using GATK in a large number of blood samples, as well as any variant identified using Platypus in the DigiPico data, is assumed to be true. (8) A binary classification model was then trained using MutLX with UTD and phylogenetic SNPs to extract clone-specific variants from the UTD. Figure 6 ).

[0128] Figure 6 MutLX analysis algorithm. (1) UTD variants are identified by subtracting WGS data from DigiPico data. (2) UTD and SNP are used as training sets to train a primary binary classification model. (3) This model is used for preliminary analysis of the training set, which allows (4) the generation of an improved training set. (5) The improved training set is then used to generate a classification model, which (6) can be used for UTD analysis. (7) A "probability score" indicates the likelihood of a mutation becoming true, and (8) an "uncertainty score" is calculated for each mutation using the model as a measure of the unreliability of the calculated probability score. (9) True variants are identified by a high probability score (low uncertainty score) with high certainty.

[0129] Figure 7 Run the probability scores for D1110, D1111, DE011, and GM12885. A cutoff of 0.2 removed most FP variant detections from the UTD while preserving almost all germline SNPs across all samples.

[0130] Figure 8 Data simulations confirmed a negative correlation between AUC and the number of true UTD variants. In runs D1110 and DE111, various numbers of somatic mutations were artificially mislabeled as UTDs (UTD*) to achieve UTD* / UTD ratios of 1%, 5%, and 10%. Since both runs were performed on 200 pg of purified DNA from large batches of tumor samples, it was expected that neither would contain true UTD variants. Independent analysis of each simulated dataset using MutLX confirmed a negative correlation between the number of true UTDs and the number of AUCs in the datasets. Notably, UTD* / UTD ratios as low as 1% (16 and 36 variants in runs D1110 and DE11, respectively) appeared to lower the AUC, suggesting that even a small fraction of true clone-specific variants can disrupt the ROC curve.

[0131] Figure 9. Analysis of the synthetic DigiPico dataset. In runs (A) D1110 and (B) DE111, various numbers of high-confidence somatic mutations were artificially mislabeled as UTDs (UTD*), thus amplifying the number of true UTDs by various UTD* / UTD ratios. Since both runs were performed on 200 pg of purified DNA from large batches of tumor samples, it was expected that neither would contain true UTD variants. The results indicate that the presence of true UTD variants in most artificial UTDs does not appear to compromise the integrity of the MutLX-generated classification models. Each boxplot shows the results for 10 different UTD* subsets used in the analysis. To achieve comparable FP rates between runs, cutoff values ​​were used to achieve 90% and 95% TPRs, respectively, across all synthetic datasets in runs D1110(A) and DE111(B). Boxplots show the median, interquartile range, and range excluding outliers. Outliers were defined as being 1.5 times higher or lower than the interquartile range.

[0132] Figure 10 Targeted sequencing of some clone-specific variants identified in D1111 was performed. Amplicon sequencing of the target sites was performed on the MiSeq platform. Three of the 14 targets appeared to have high noise levels in the blood sample (highlighted in orange) and were therefore considered indeterminate. Of the remaining 11 mutations, only one appeared to have no evidence in large batches of DNA samples from PT2R tumors (highlighted in blue). VAWF: Variant Allele Well Fraction.

[0133] Figure 11 Targeted sequencing of some artificial variants identified in DE111 was performed.

[0134] Figure 12 Frequency of various mutation types in FP detection identified by MutLX in DigiPico data. Green dots (on the left) are obtained from the analysis of somatic variants from standard WGS data from patients #11152, #11513, and OP1036. Red dots (on the right) are from all artificial mutations identified from DigiPico data from the same patients via MutLX. The black line represents the median value for each group. The higher proportion of C>A mutations among those eliminated by MutLX is consistent with previous studies, suggesting that oxidative DNA damage during library preparation leads to the formation of artificial C>A mutations.

[0135] Figure 13 Image of IGV for SNVs identified in the subclone *kataegis* of the PT2R sample. Note that almost all mutations are in the form of C>T or C>G mutations on the forward strand of the genome.

[0136] Figure 14. Comparison of DigiPico and DigiPico2 workflows. (A) The DigiPico workflow took nearly 12 hours and consisted of 7 steps, 5 of which occurred in a 384-well plate. (B) The DigiPico2 workflow took only 4.5 hours and consisted of 5 steps, 3 of which occurred in a 384-well plate. The blue reaction occurred in the 384-well plate format, the green reaction occurred in 16 wells, and the orange reaction occurred in one tube.

[0137] Figure 15. Comparison of indexing strategies for DigiPico and DigiPico2. (A) Asymmetric ligation of two annealed indexed oligonucleotides introduces an i5 index and an i7 index. Combinations of these indexes can generate 384 different indexes in DigiPico. Note that each strand receives one i5 index and one i7 index, so there is no redundancy. (B) In DigiPico2, the initial column index (Ci) is introduced through efficient ligation with a circular adapter. Next, indexing PCR is performed using row index primers (Ri). Note that each strand receives both the Ci and Ri indexes twice, which introduces redundancy required to remove index cross-contaminants.

[0138] Figure 16. Comparison of results from DigiPico and DigiPico2. (A) Both WGA products appear to be relatively homogeneous throughout the plate. The values ​​represent relative fluorescence from measurements of EvaGreen (RFU). (B) In DigiPico2, the index frequency of each well in the plate appears to correlate better with the RFU value of the WGA product. (C) This fact can be quantified using a correlation plot. (D) MutLX analysis of the DigiPico2 data appears to better distinguish between real and artificial mutations. Note the presence of lower mutations in the upper right portion of the plot. This area marks artificial mutations that obtained high probability scores in the analysis.

[0139] Figure 17. Evaluation of the single-cell DigiPico (ScDigiPico) sequencing method. To evaluate the effectiveness of the ScDigiPico method in identifying active mutational processes, we simulated such a process using Kuramochi cells cultured with N-ethyl-N-nitrosourea (ENU) mutagenesis. ENU is an alkylating agent and a very potent mutagen that preferentially induces T>C, T>A, and C>T mutations. In this setting, each cell was expected to acquire a different set of mutations, but since the underlying mechanism of mutagenesis was the same, the mutations were expected to belong to the same type. (A) Cultured Kuramochi cells were exposed to 0.1 g / LENU for 48 hours. Most cells died in the presence of ENU, but those that survived accumulated a large number of mutations. Single-mutant Kuramochi cells were then sorted into each well in the first column of a 384-well plate. During ScDigiPico, the DNA content from each cell was evenly distributed to all wells in that row prior to WGA and library preparation. (B) Analysis of the ScDigiPico results of the mutant cells showed, as expected, that the frequency of the new mutations clearly matched the frequency of ENU. (C) In a single cell analyzed after UV irradiation, we found a distinct Kataegis event on chromosome 7.

[0140] Figure 18 DigiPico / MutLX eliminates false positives in whole-genome amplified DNA. A dot represents a single sample (approximately 20 cancer cells) from sequencing whole-genome amplified blood DNA (starting from a pint of blood DNA) from a cancer patient. Germline (blood) DNA should not have thousands of unique variants compared to standard DNA sequencing. However, existing methods have detected tens of thousands of such false positive mutations. In contrast, DigiPico / MutLX eliminates these false positives. Detailed Implementation

[0141] Example 1 - Using DigiPico / MutLX to reveal active mutational processes in tumors with unprecedented accuracy

[0142] summary

[0143] Bulk whole-genome sequencing (WGS) can analyze tumor evolution, but due to depth limitations, it can only identify older mutational events. Discovering mutational processes currently used to predict tumor evolutionary trajectories requires intensive sequencing of individual clones or single cells. However, such studies are inherently problematic because too many false-positive mutations are found when sequencing picogram-sized amounts of DNA. Data pooling increases the confidence of discovered mutations, tracing the discovery back to a common ancestor. Here, we report a robust whole-genome sequencing and analysis pipeline (DigiPico / MutLX) that virtually eliminates all false-positive results while retaining an excellent proportion of true positives. Using our method, we identified, for the first time, a hypermutation (kataegis) event in approximately 30 cancer cell groups in recurrent ovarian cancer. This was impossible to identify from bulk WGS data. Overall, we present the DigiPico / MutLX approach as a powerful framework for identifying clone-specific variants with unprecedented accuracy.

[0144] introduction

[0145] In this work, we developed a single-DNA molecule WGA and sequencing method to obtain high-quality and data-rich sequencing results from picogram-sized amounts of DNA obtained from clinical samples (which we call DigiPico; for use in...). Pique DNA number Sequencing). Furthermore, we used an algorithm based on artificial neural networks (ANN) (MutLX, for...) Mutation learning A supplemental analysis workflow was implemented on the DigiPico data to eliminate false positives while maintaining excellent sensitivity to true positive mutations at the whole-genome scale. We validated our methods using data from extensively sequenced tumors from a single patient, where the cumulative depth obtained from 45 whole-genome sequencing runs of DNA at three different time points was approximately 4200-fold. We demonstrated the versatility of these methods by sequencing samples from four additional cancer patients and lymphoblastic cell lines.

[0146] Materials and Methods

[0147] Patient samples and consent forms

[0148] Patients #11152, #11502, and #11513 provided written consent to participate in the prospective biomarker validation study, the Gynaecological Oncology Targeted Therapy Study 01 (GO-Target-01), with research ethics approval number 11 / SC / 0014. Patient OP1036 participated in the prospective Oxford Ovarian Cancer Predict Chemotherapy Response Trial (OXO-PCR-01), with research ethics approval number 12 / SC / 0404. Necessary informed consent forms were obtained from study participants as appropriate. Blood samples were collected on the day of surgery. Tumor samples were biopsied during laparoscopy or debulking surgery and immediately frozen on dry ice. All samples were stored in clearly labeled cryovials at -80°C.

[0149] cell lines

[0150] The GM12885 lymphoblastoid cell line (RRID: CVCL_5F01) was obtained from the Coriell Institute, and the cells were cultured according to the provider's recommendations.

[0151] Slices and LCM

[0152] Frozen tumor samples were embedded in an OCT (NEG-50, Richard-Allan Scientific) and 10–15 μm sections were obtained using an MB DynaSharp microtome blade (ThermoFisher Scientific) on a CryoStar cryostat microtome (ThermoFisher Scientific). The tumor sections were then transferred to PEN-lined glass slides (Zeiss) and immediately stained on ice (2 minutes in 70% ethanol, followed by 2 minutes in 1% cresol purple (Sigma-Aldrich) in 50% ethanol, then rinsed with 100% ethanol). Individual tumor islands were ejected into 200 μl opaque adhesive caps (Zeiss) using a PALM Laser Microdissection System (Zeiss).

[0153] Standard WGS and Data Analysis

[0154] DNA was extracted using the DNeasy Blood and Tissue Kit (Qiagen). Up to 1 μg of DNA was fragmented in 50 μl of water using a Covaris S220 focused ultrasound instrument to obtain fragments of 250–300 bp. The resulting DNA fragments were then used for library preparation using the NEBNext Ultra II Library Preparation Kit (NEB) according to the manufacturer’s protocol. The generated libraries were sequenced on an Illumina NextSeq or HiSeq platform to a depth of 30–40 times that of the human genome. Sequencing reads in FastQ format were initially pruned using TrimGalore (14) and then mapped to the human hg19 genome using Bowtie2 (15). Germ variant identification was performed using GATK’s HaplotypeCaller (16). Somatic variants were detected using Strelka2 with a variant allele fraction cutoff of 0.2 (17).

[0155] DigiPico sequencing

[0156] First, denature 200 pg of purified DNA, 20–30 resuspended cell nuclei, or laser-captured microdissectioned tumor islands using 5 μl D2 buffer from the Repli-g Single Cell Kit (Qiagen). After incubating at room temperature for 5 minutes, add 95 μl of water to the sample, and then use a Mosquito HTS liquid processor (TTP Labtech) to add 200 nmol of denatured template to each well of a 384-well reaction plate. Each well already contains 800 nmol of WGA mixture (0.58 μl Sc reaction buffer, 0.04 μl Sc polymerase (REPLI-g Single Cell Kit, Qiagen), 0.075 μl 1 mM dUTP (Invitrogen), 0.04 μl EvaGreen 20x (Biotium), and 0.065 μl water). Incubate the plate at 30 °C for 2 hours, then heat-inactivate at 65 °C for 15 minutes. If necessary, EvaGreen was added to the reaction to allow monitoring of the WGA reaction using a real-time PCR machine (18). Next, the whole-genome amplified DNA was subjected to controlled enzymatic fragmentation (19) reaction steps without any purification steps. Briefly, (A) 1200 nl of UDG mixture (0.08 U / μl rSAP (NEB), 0.2 U / μl UDG (NEB), 0.4 U / μl EndoIV (NEB) in 1.8 times NEBuffer 3) was added, incubated at 37 °C for 2 h, and heat-inactivated at 65 °C for 15 min. (B) 1200 nl of PolI mixture (0.4 U / μl DNA polymerase I (NEB), 0.25 mM dNTP, 8 mM MgCl2 and 0.8 mM DTT) was added, incubated at 37 °C for 1.5 h, and heat-inactivated at 70 °C for 20 min. (C) Add 1200 nmol of Klenow mixture (0.5 U / μl Klenow exo-(NEB), 0.5 mM dNTP, 8 mM MgCl2 and 0.8 mM MDT), incubate at 37 °C for 45 min, and heat-inactivate at 70 °C for 20 min. (D) Add 400 nmol of 20 μM full-length Llumina linker oligonucleotides with well-specific indexes (Table S1) to each well, then add 1100 nmol of ligation mixture (40 U / μl T4 DNA ligase (NEB), 5 mM ATP, 11.5% PEG 8000 (Qiagen) and 6.8 mM MgCl2), incubate at 20 °C for 30 min, and heat-inactivate at 65 °C for 15 min.

[0157] The resulting products were then combined and the DNA was precipitated using an equal volume of isopropanol. The DNA was then resuspended in water and the product was subjected to dual size selection using Agencourt AMPure XP SPRI magnetic beads (Beckmann coulter) at a bead ratio of 0.45 for left-side selection and an additional 0.32 for right-side selection. The purified DNA was then resuspended in water and immediately used for limited-cycle PCR amplification using a mixture of P5 and P7 primers (Table S1). PCR was performed for 12 cycles, annealing at 55°C for 10 seconds and extending at 72°C for 45 seconds. The final product was purified by magnetic beads at a ratio of 0.9. The resulting library was then sequenced on an Illumina sequencing platform in 2x150 paired-end sequencing mode to achieve 30–40 times coverage depth on the human genome. Currently, the additional processing steps required for DigiPico library preparation increase the total reagent cost by nearly £250.

[0158] DigiPico sequencing data analysis

[0159] The analysis workflow for DigiPico sequencing data is provided in the supplement. Figure 5In summary, after multiplexing the Illumina sequence data, 384 paired reads FastQ files were obtained for each well. The FastQ files were pruned for adapter sequences and quality (14). The first 12 nucleotides of each read were also removed. The pruned reads were mapped to the human hg19 reference genome using Bowtie2 (15) with the ignore-quals parameter activated, and duplicate reads were marked using the Picard tool (20). Joint variant identification was then performed on all 384 individual bam files and the merged bam files from all wells using the Platypus variant identification program (21). Next, all low-quality variants were removed by applying quality filters (QUAL>60, FR>0.1, HP≤4, QD>10, and SbPval≤0.95). In addition, the total number of wells covering each site (Tw) and the number of wells supporting each variant (Vw) were determined, and well count filters (Tw>5, Vw>2, and Vw / Tw>0.1) were applied to retain only high-confidence sites for analysis. Finally, all genomic regions with poor mapping ability were removed from the analysis using VCFtools (23) (22). Then, variant re-identification (genotyping) was performed on WGS data from blood and bulk tumors using Platypus (minPosterior parameter set to 0, minMapQual parameter set to 5) with the obtained high-confidence de novo DigiPico variant list. Any variants that were confidently unsupported in standard WGS data were extracted as UTD (DigiPico-specific) variants. Any variants that were also confidently present in bulk sequencing data (based on GATK analysis) of blood samples were extracted as TP (true positive) variants. Figure 5 ).

[0160] MutLX algorithm

[0161] MutLX analysis pipeline summary Figure 6 middle.

[0162] Artificial Neural Network Architecture

[0163] The neural network model used in this study is a multilayer perceptron with an input layer consisting of N neurons (N=41), where N is the number of features used in each experiment, and implemented in Python 3 using Keras (24). The model has two hidden layers with ReLU activation. We varied these numbers but did not see any significant improvement when using a large number of neurons. The final layer is a single output neuron with sigmoid activation. The loss function is binary cross-entropy. For training, we applied stochastic gradient descent optimization with momentum (Adam (25)) with a learning rate of 0.001, a batch size of 8, and 10 epochs. After 10 epochs, we did not observe any additional performance improvement.

[0164] Features used for training

[0165] The following features extracted from the Platypus output of the DigiPico data are used as input to the neural network model:

[0166] Platypus quality parameters: QUAL, BRF, FR, HP, HapScore, MGOF, MMLQ, MQ, QD, SbPval, NF, NR, TCF and TCR (21).

[0167] Sequence environment complexity: F 20 [1], F 20 [2], F 20 [3]. Among them, F 20 [i] is the sum of the frequencies of the i most abundant nucleotides in the 10bp sequence flanking the variant position.

[0168] Reading distributed data: R 合并 [ref+var]、R 合并 [var], W[R[ref]>0s and 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 and R[ref]=0][0 / 1], W[R[ref ]>0 and R[var]>0], W[R[ref]=0 and R[var]>0], W[R[ref]=0 and R[var]>1], W[R[ref]=0 and R[var]>2], W[R[ref]=0 and R[var]>3], W[R[ref]=0 and R[var]>4], W[R[ref]=0 and R[var]>5], R max [1][var]、R max [2][var]、Rmax [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 ). Among them, R 合并 [x] represents the total readings supporting allele x in the merged BAM file (ref represents the reference allele, var represents the variant allele). W[i][j] shown here represents the number of wells matching standard i with the reported genotype j. In standard i, R[x] represents the readings in a specific well supporting allele x. R max [y][x] shows the readings in the well that support the y-th highest reading for allele x. Finally, Max c Max is the number of variant-supporting wells in the column with the most wells supporting variant alleles. r It is the number of variant-supporting wells in the row with the most wells supporting variant alleles.

[0169] Training using MutLX

[0170] For each DigiPico run, we treat the entire training set as a collection of all UTD variants (labeled 0) and heterozygous SNPs (labeled 1). In this set, the number of UTD variants is much smaller than the number of heterozygous SNPs, making the set imbalanced. Therefore, to avoid bias towards specific labels during training, we create 25 distinct balanced training subsets for each DigiPico run. This is done so that each training subset consists of all UTD variants and a randomly selected subset of heterozygous SNPs (of equal size to the number of UTD variants). As mentioned earlier, most UTD variants are FP variant detections, where the proportion of true clone-specific variants is unknown, thus making the 0-label noisy. To account for these noisy labels and perform two-step training, we employ the following strategy: After training an initial model on each balanced training subset, the resulting model is applied to mutations across the entire training set to obtain initial probability values ​​for each mutation. These probability values ​​represent the predicted probability of a mutation belonging to the label 1 category. Therefore, any 0-labeled mutation that obtains a predicted probability value close to 1 is likely a mislabeled mutation. Therefore, to reduce the level of mislabeled data in the training set, all UTD variants with a probability greater than 0.7 and all phylogenetic SNPs with a probability less than 0.3 were considered mislabeled and removed from the training set. The cutoff values ​​in this step were determined empirically through analysis of various simulated datasets. Finally, following a subsampling strategy similar to the initial training, a new model was trained on the remaining mutations in the training set. This model was then used to analyze all UTD variants.

[0171] Calculation of probability and uncertainty fractions

[0172] As mentioned earlier, in MutLX, the training process is repeated 25 times using different randomly selected subsets of phylogenetic SNPs, each time producing a different model, thus each mutation has 25 different predicted probability values. Therefore, we define the "probability score" of each mutation as the average of all its predicted probability values:

[0173]

[0174] Where P i It is the probability value obtained from the i-th training subset, where n represents the number of subsets (n = 25).

[0175] Furthermore, to obtain an uncertainty estimate for each probability value, we performed a test-time dropout analysis (26). The trained model was applied to each mutation for 100 iterations, during which different neurons in the first and second hidden layers of the neural network were dropped at rates of 0.8 and 0.7, respectively. This process generated 100 probability values ​​for each mutation. Based on these values, we defined the “uncertainty score” for each mutation as the average of the dropout variances from 25 distinct subsets:

[0176]

[0177] Where σ i 2 It is the variance of 100 probability values ​​obtained from the missing data analysis of the i-th training subset, where n represents the number of subsets (n = 25).

[0178] The probability score is higher than 0.2. Figure 7 Uncertainty scores for all variants of the sample were used to generate a putative acceptor operating characteristic (ROC) curve. This curve was generated by considering a cutoff range of uncertainty scores between 0.0 and 0.25. At each cutoff, the ratio of phylogenetic SNPs with uncertainty scores below the corresponding UTD number was plotted. The area under the curve (AUC) was then calculated after normalizing the number of UTDs to 0 to 1. Note that in the case where the true clonal-specific variant is not expected (all UTDs are detected by FP), assuming a perfect model, this plot represents the ROC curve, and the AUC of this plot should be close to 1. Conversely, a significant decrease in AUC indicates the presence of a true clonal-specific variant in the sample. This negative correlation between the true UTD number and AUC was validated using a simulated dataset. Figure 8 Based on these observations, for samples where the ROC curve indicates the presence of true clone-specific variants (AUC < 0.9), MutLX uses an "uncertainty score" cutoff value that results in a 95% TPR to improve the recovery rate of clone-specific variants. For datasets with AUC ≥ 0.9, the cutoff value used to filter the data is determined based on the intersection of the threshold curve and the ROC curve.

[0179] Simulate the generation and analysis of the DigiPico dataset

[0180] The simulated data was used to: (a) verify the negative correlation between the actual number of UTDs and AUC. Figure 8 (a) and (b) ensure that overfitting to potentially real clone-specific variants does not occur (Figure 9).

[0181] To generate the simulated dataset, we first identified somatic mutations in a large batch of WGS data from the PT2R tumor sample of patient #11152 using the Strelka2 somatic variant identification procedure. These somatic variants were then identified in the de novo variant identification data run D1110, with any somatic variant with Tw > 6 and Vw / Tw > 0.45 selected as high-confidence somatic variants. Next, various numbers of randomly selected high-confidence somatic variants were artificially mislabeled as UTDs (UTDs*) to achieve UTD* / UTD ratios of 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, and 0.1. The resulting list of synthetic variants was then used independently for MutLX analysis, and the number of UTDs and UTD* filtered by MutLX was calculated for each run. To ensure robust analysis, a subset of 10 distinct somatic variants was analyzed for each ratio. A similar analysis was performed on DigiPico data DE111, which was obtained from a large-volume DNA extract of ascites sample from patient #11513.

[0182] Verification of the MutLX algorithm

[0183] Tumor samples from patient #11152 were used for PT2R validation of the MutLX algorithm. A small tumor fragment was macroscopically dissected from the frozen specimen and embedded in OCT medium for sectioning. The first portion (15 μm) of the tumor was collected in a separate tube, and the cell nuclei were resuspended in 50 μl of sterile PBS solution. The total number of cell nuclei in the suspension was measured, and a volume containing 30 cell nuclei was directly denatured using an equal volume of D2 buffer from the Repli-gmini WGA kit (Qiagen). The resulting crude lysate was used directly for DigiPico library preparation running D1111. The remaining tumor sample was then used for bulk DNA extraction using the DNeasy Blood and Tissue Kit (Qiagen). The resulting 200 pg of DNA was used directly for DigiPico library D1110 preparation. 1 μg of DNA was used for standard library preparation using the NEBNext Ultra DNA Library Preparation Kit (NEB). In this setup, running D1111 alone is expected to yield true clonal-specific variants. Since the template used in run D1110 is a subset of the templates used in bulk WGS analysis, almost all true variants in run D1110 will also appear in the WGS data at similar frequencies and therefore will not be identified as UTDs. Similar logic applies to the results of DigiPico runs DE011 and GM12885. Since these DigiPico runs were performed on 200 pg of DNA from bulk DNA extracts, it is expected that no true UTD variants will appear in these samples. However, it is also worth noting that due to the digitized nature of the data, variants with very low frequencies (<0.05%) will show inflated variant allele fractions in runs D1110, DE011, and GM12885, as such variants are unlikely to appear in more than one well and will be eliminated from the Vw-based filter data. Therefore, it is safe to assume that almost all UTD variants in these runs are FP detected.

[0184] Application of SCcaller on DigiPico data

[0185] SCcaller was originally developed for analyzing single-cell sequencing data with multiple substitution amplification (11). Since DigiPico library preparation also requires multiple substitution amplification of a limited amount of template DNA, the resulting data is largely similar to the natural input of SCcaller. Therefore, we used the merged bam file of the DigiPico data as input to SCcaller. For analysis, the list of heterozygous SNPs was obtained from their respective bulk WGS data using GATK HaplotypeCaller, with a cutoff value of α = 0.01. Subsequently, all filtered SNVs were used for variant re-identification against their respective standard WGS data, and all variants that were definitely not supported by the WGS data were extracted as UTD variants.

[0186] Mutation verification

[0187] Variants analyzed by MutLX were validated by comparison with deep sequencing data from large batches of tumors from independent sequencing platforms. All DigiPico data from patient #11152 were validated (27) by comparison with 39 deep sequencing datasets obtained from the same tumor blocks sequenced on the Complete Genomics sequencing platform. This included three Complete Genomics batch sequencing datasets and 36 LFR (long read) sequencing datasets. Since the independent sequencing data for the omental tumors were not obtained from the exact same tumor blocks used for DigiPico sequencing, the validation rate for these runs through such comparisons is not expected to be high.

[0188] For targeted validation, primers were designed using the Primer3 tool to obtain amplicones containing variants (Table S1). The amplicon was obtained by using... High-Fidelity PCR Master Mix and GC buffer were used to perform 16 cycles of two-step PCR on 1 ng template. All amplicons from each sample were then pooled and purified before adapter ligation and indexing using the NEBNext Ultra II kit. The resulting libraries were sequenced on the MiSeq platform. The sequencing results were mapped to the human hg19 genome using Bowtie2, and the number of reads supporting each variant was calculated using the Platypus variant identification program.

[0189] Local hypermutation (kataegis) analysis

[0190] To generate the rainfall plot, a custom script in R was used to plot the distances between contiguous somatic mutation pairs on chromosome 17 against the genomic location of the second mutation in each pair. The presence of local mutation clusters indicates a kataegis event. In these plots, each point is colored based on the mutation type of the second mutation in that pair relative to the hg19 human reference genome.

[0191] result

[0192] Implementation of the DigiPico sequencing method

[0193] A key characteristic of amplification errors, with or without prior DNA damage, is that they are randomly introduced during amplification (6,7,28). Therefore, we hypothesize that when amplifying and sequencing a single DNA molecule, artificial mutations will only be present in a subset of reads generated from sequencing the original single DNA molecule. In contrast, true variants are expected to appear in all such reads. Partitioning the template DNA into individual compartments before WGA, such that each compartment receives no more than one DNA molecule from each site, will produce such single-DNA molecule sequencing data (…). Figure 1B Since artificial mutations result in compartments with reads supporting multiple alleles, this approach is able to identify these artifacts and allow for the elimination of FP variant identification. Furthermore, this partitioning approach results in independent WGA responses at each gene locus. Therefore, multiple internal replication data are provided for the WGA process. While true variants are expected to appear regularly in replication, artificial mutations may have limited presence in a few compartments due to their randomness. Therefore, considering both of these factors, this WGA and sequencing approach can produce a distinct distribution pattern of artificial mutations within compartments compared to true mutations. Since artificial neural networks have demonstrated the ability to extract complex patterns from high-dimensional inputs, they are good candidates for identifying and eliminating false positive mutations in such data. While previous partitioning and sequencing methods for obtaining haplotype information have been described, none have been used to differentiate between true and artificial mutations (19,29,30).

[0194] To fully benefit from the data richness of partitioning and sequencing methods for accurate genomic studies of clinical samples, we developed DigiPico sequencing (… Figure 1C To perform DigiPico sequencing, we first evenly distributed nearly 200 pg of DNA (obtained from 20–30 human cells) into each well of a 384-well plate. This ensured that the probability of two different DNA molecules from the same site coexisting in the same well was less than 10% (19). After WGA, each well was independently processed into an index library, and each library received a unique barcode sequence before pooling and sorting. Figure 1C Because the key differentiating factor of artificial mutations in our method lies in their unique distribution patterns, the uniform distribution and amplification of DNA molecules, as well as the consistent depth of cross-well sequencing coverage, are crucial. Achieving this homogeneity ensures that the difference between the distribution patterns of real and artificial mutations is maximized. To ensure the desired homogeneity, we monitored the progress of the WGA reaction and quantified the final results for all wells during each DigiPico library preparation. The former was achieved by adding EvaGreen dye to the WGA reaction and monitoring the fluorescence intensity in real time every 5 minutes. EvaGreen is an intercalating dye that binds to the small grooves of DNA and therefore does not interfere with the isothermal WGA reaction (18). For the latter, we introduced a per-well qPCR step to measure the relative number of adapter-ligated fragments in each well using adapter-specific primers before pooling. Only libraries that passed both homogeneity tests were used for sequencing. Figure 1D and 1E Importantly, we also reduced the WGA reaction volume to 1 μl. This was only possible after identifying a compatible multiple substitution amplification (MDA) method. Comparing six different MDA strategies, REPLI-g single-cell amplification was the only method that met the sensitivity and selectivity required for our purposes. Figure 1F Miniaturization of the reaction allows us to streamline the library preparation process in a single 384-well plate without the need for intermediate purification steps using readily available automated pipetting instruments. Finally, we aimed to optimize the DigiPico library preparation process for frozen clinical samples. This was achieved by directly performing a WGA reaction on crude lysates of small clusters of adjacent cells (tumor islands) separated by LCM (laser capture microdissection). This strategy ensures minimal loss of genomic material while minimizing operation time, thereby reducing the chance of template oxidation.

[0195] The DigiPico sequencing platform generates high-quality libraries from limited clinical samples.

[0196] After optimizing all necessary aspects of the DigiPico library preparation process, we decided to evaluate the quality of DigiPico libraries obtained from clinical samples. To this end, we prepared DigiPico libraries D1110 and D1111 from frozen recurrent tumor samples (PT2R) obtained from a patient with high-grade serous ovarian cancer (#11152). In this experiment, while library D1110 was prepared from 200 pg of template extracted from a large batch of DNA extracts from the PT2R sample, library D1111 was prepared directly from a small remaining frozen fragment of the tumor sample (containing nearly 30 cancer cells). Each library was sequenced on the Illumina NextSeq platform to obtain nearly 400,000,000 reads in 150x2 paired-end format. Preliminary evaluation of the obtained sequencing data indicated that both libraries D1110 and D1111 produced high-quality sequencing data with overall mapping rates of 91.35% and 94.27% of the human hg19 genome, respectively. Figure 1G Analysis of the uniformity of the cross-plate distribution showed that, in runs D1110 and D1111, each well covered an average of nearly 4.4% and 4.6% of the genome, respectively, with average depths of 1.7x and 2.1x per well, demonstrating excellent cross-plate uniformity. Figure 1H and 1I This resulted in cumulative coverage of 92.1% and 91.1%, respectively, with depths of 30x and 43x per run. These results confirm that DigiPico sequencing can be used to generate high-quality sequencing data with excellent coverage from a limited number of frozen clinical samples.

[0197] Finally, we evaluated whether our initial hypothesis regarding the unique distribution patterns of different mutation types held true in the actual DigiPico dataset. To this end, we hypothesized that any variants shared between the DigiPico dataset and standard batch sequencing data of the same tumor sample must be authentic variants. These should consist primarily of germline SNPs and clonal somatic variants. Therefore, by definition, all FP variant detections and most clone-specific mutations (if present in the studied samples) would belong to variants present only in the DigiPico data and not in the bulk WGS data. For simplicity, these variants are referred to as UTDs (DigiPico-specific) below. Therefore, given that the standard bulk sequencing data for PT2R samples were obtained from the same DNA extract used for D1110 library preparation, almost all UTD variants in the D1110 DigiPico run should be artifacts (…). Figure 1J Conversely, the UTD in running D1111 may contain some clone-specific mutations as well as artificial mutations (). Figure 1JTherefore, we used the UTD variants in run D1110 as representatives of the artificial mutations in our analysis. The frequency of wells with coexisting alleles at the same locus was compared in run D1110. Figure 2A ) and the number of holes supporting each variant ( Figure 2B The results showed that UTD had a significantly higher proportion of the former and a lower number of the latter compared to any other category of mutations (for both analyses, p < 2e-16, one-way ANOVA, and then the Tukey HSD test). This clearly supports the hypothetical different distribution patterns of artificial mutations in the DigiPico dataset.

[0198] DigiPico's MutLX Analysis Pipeline

[0199] After obtaining high-quality data using DigiPico sequencing, we decided to implement an analysis pipeline to eliminate FP variant detection based on mutation distribution patterns. As mentioned earlier, the ANN algorithm is ideally suited to handle such complex patterns. Given a representative set of correctly labeled examples (training set), the ANN can learn to classify mutations without any class-specific information. However, implementing the ANN algorithm for eliminating FP mutations from sequencing data presents two main challenges: (a) the difficulty in obtaining a generalizable model and (b) the inability to obtain a training set with representative and accurate labels. First, it is impossible to generate a model that can generalize for analyzing each DigiPico dataset because mutation distribution patterns depend on various run-specific initial conditions that are not easily interpreted (e.g., copy number state of the genome). Therefore, a run-specific model tailored to each DigiPico run is required. This means selecting a subset of run-specific mutations as the training set for each DigiPico run. Second, while it is easy to extract correctly labeled real mutation examples from known SNPs in the genome, it is impossible to determine a representative and accurate set of artificially labeled mutation examples. To address this issue, we consider the UTD to be a reasonable approximation of a representative set of artificial mutations, assuming that the UTD is primarily composed of such mutations. However, this assumption can lead to key challenges. By definition, the UTD consists of both artificial mutations and real clone-specific mutations. While artificial mutations are expected to be abundant in all DigiPico runs, real clone-specific mutations may appear at varying frequencies depending on the sample. Figure 1J Therefore, when all UTDs are treated as examples of artificial mutations, samples with more clone-specific variants will have a noisier training set. If this is not considered, it could put samples with true clone-specific variants at an analytical disadvantage, as a noisier training set leads to a worse classification model. Figure 2CSpecifically, in samples with true clone-specific variants, the presence of real mutations in artificially mutated examples may lead to overfitting of the model to these variants. Figure 2C This could weaken the model's ability to accurately identify true clone-specific variants in such samples. Given that the primary goal of DigiPico sequencing is to identify clone-specific variants, it is essential to ensure that such overfitting does not occur when analyzing the DigiPico dataset.

[0200] Considering all the aforementioned limitations and problems, we designed and implemented an ANN-based binary classifier, MutLX, for analyzing the DigiPico dataset. The focus of the DigiPico analysis pipeline is to effectively eliminate FP detections and accurately identify true clone-specific variants from UTDs. To address the problem of training with an imperfect training set, we adopted the following approach when training MutLX. Initially, we treated all UTD variants as examples of artificial mutations (labeled 0) and a similar number of randomly selected heterozygous SNPs as examples of true variants (labeled 1). Since most UTD variants are FP detections and the proportion of true clone-specific variants is unknown, the 0 label is considered "noise" at this stage. In other words, although true clone-specific variants should be labeled 1, they are still labeled 0 due to their anonymity at this stage, and so on. To accommodate this type of noise in the training dataset, we adopted a two-step training process ( Figure 2DThe first step involves training a model given all labeled example data. This initial model is then used to calculate the probability that each mutation belongs to its labeled class. Based on these model predictions, any mutations that appear to be mislabeled are temporarily removed from the dataset (pruning). The fundamental assumption here is that even though a model trained on noisy samples may not be as robust as one trained on a hypothetically clean dataset, it will still be biased towards better predictions of correct examples because they have a higher proportion of noise in the dataset. Therefore, examples predicted by the model for its original labels are likely to be mislabeled from the outset. The second step involves training a new classification model on this pruned training set. Since the second training set likely contains fewer mislabeled data points, the final model is expected to identify true mutations more effectively, regardless of whether all UTD variants are indeed artifacts (31, 32). We then use this model to assign a “probability score” to each hypothetical mutation. This score represents the likelihood that a given mutation belongs to the true variant class. While this two-step training process promises to significantly improve the classification model, the final model remains error-prone due to the imperfections of the training set. Therefore, we add another layer of analysis to further improve the accuracy of our pipeline. This is achieved by assigning an uncertainty estimate to a “probability score” for each mutation. This uncertainty estimate is based on the assumption that the majority of activated neurons in the ANN hidden layers support robust predictions. Therefore, any subset of these neurons will always produce similar probability scores, and thus the differences between the various “probability scores” obtained from different subsets of neurons will be small. Figure 2E In contrast, the seemingly high "probability score" of artificial mutations is most likely supported by only a few neurons in the ANN hidden layers. Therefore, different subsets of neurons will result in different "probability scores," leading to significant differences in the scores of artificial mutations obtained from different subsets of neurons. Figure 2E Therefore, the "uncertainty score" can be calculated as the variance of the "probability scores" obtained from multiple randomly selected subsets of neurons in a MutLX-trained ANN (26). Thus, the combination of the "probability score" and "uncertainty score" for each mutation should enable us to accurately determine whether the detected variant is a genuine mutation or the result of artificial alteration in the template. Figure 6 ).

[0201] Verification of the MutLX algorithm

[0202] To validate our strategy, we chose to test the MutLX analysis pipeline on runs D1110 and D1111. This is because these DigiPico runs were obtained from HGSOC, which had previously performed extensive sequencing using data from 48 independent whole-genome sequencing datasets spanning three different time points (patient #11152), with a total depth of approximately 4200x, from two independent sequencing platforms (33). To our knowledge, this includes the most extensive whole-genome sequencing of tumors to date. This exceptionally large dataset allows for reliable cross-validation of mutations in this tumor. For this purpose, we used the MutLX algorithm to analyze sequencing data from runs D1110 and D1111. As previously mentioned, when comparing these DigiPico datasets with large-volume sequencing data from PT2R sites, true UTD variants (clone-specific variants) were expected to appear only in run D1111, while almost all UTDs in run D1110 were expected to be artifacts (clone-specific variants). Figure 1J Furthermore, we analyzed purified DNA from blood samples (run DE011) and DigiPico sequencing data prepared from cultured GM12885 lymphoblastoid cells, neither of which were expected to contain true UTD mutations. De novo variant identification was performed on these DigiPico runs, followed by initial filtering based on well number, resulting in the identification of thousands of UTD variants in each sample, almost all of which were considered FP detections. However, applying the MutLX algorithm to the UTD variants in runs D1110, DE011, and GM12885 resulted in efficient elimination of over 99% of FP variant detections in these runs, targeting only 4, 7, and 3 genome-wide FP mutations, respectively, while maintaining approximately 85% sensitivity for detecting true mutations. Figure 2F and 2G In contrast, SCcaller(11) analysis of the same data resulted in the detection of 713, 712, and 13,280 FP variants, respectively. Figure 2G On the other hand, MutLX identified 264 putative clonal-specific variants in the D1111 run, of which 238 (90%) were validated by comparison with independent high-depth datasets of the tumor sample. Figure 2G Furthermore, these observations were further validated by targeted sequencing of large-volume DNA from the tumor. Thus, in addition to the 11 analytical amplicon containing clone-specific variants from running D1111, 10 amplicon strains were found that ultimately existed at low frequencies in large-volume DNA from the PT2R sample. Figure 10Furthermore, amplicon sequencing of 37 seemingly high-quality UTD variants from DE111 (these variants were flagged as artifacts by the MutLX algorithm) showed no evidence of their presence in large-volume DNA samples. Figure 11 These results clearly demonstrate that MutLX can learn an accurate classification model that distinguishes artificial mutations from real variants and can effectively identify true clone-specific variants in DigiPico data.

[0203] Furthermore, we investigated whether the presence of true clone-specific mutations could impair model sensitivity due to overfitting. To this end, we artificially labeled different numbers of somatic mutation mislabels in runs D1110 and DE111 as artificial UTD variants (UTD*) to generate synthetic datasets with different proportions of true UTDs. These synthetic datasets were then analyzed individually using MutLX, and the FP rate and UTD* recovery rate were examined at different UTD* / UTD ratios across all synthetic datasets. The results showed that a UTD* / UTD ratio up to 10% did not significantly affect the recovery rate of UTD* variants, indicating that overfitting does not occur in MutLX (Figure 9).

[0204] Universality of DigiPico / MutLX sequencing and analysis methods

[0205] Finally, to ensure the generality of our proposed method, DigiPico sequencing was performed on various template DNA sources from four different HGSOC patients, and the resulting UTDs were analyzed using the MutLX algorithm. The results clearly demonstrate that MutLX can reliably identify and eliminate artificial variant detections from various DigiPico libraries (Table 1). This strongly suggests that DigiPico / MutLX can effectively investigate recently acquired mutations in solid tumors. Importantly, analysis of the frequency of different mutation types in these data indicates the presence of higher levels of C>A mutations among the identified artificial variants, consistent with the view that this FP detection is a result of oxidative damage to the template DNA. Figure 12 ).

[0206] Using DigiPico / MutLX to study active mutation processes

[0207] We next tested the feasibility of studying the mutational process in a patient with HGSOC (#11152). For this patient, various sequencing data from the prechemotherapy omental mass were available (30x standard batch sequencing and 5 DigiPico runs on tumor islands). The patient subsequently relapsed, and tumor samples were collected from the pelvis (pelvic recurrent tumor; PT2R) and para-aortic lymph nodes (PALNR) for standard batch sequencing and DigiPico sequencing of tumor islands. Analysis of the large prechemotherapy sequencing data identified 13,721 somatic mutations. Of these mutations from the DigiPico data, 84.6% were present in at least three tumor islands, and 91.4% were also present in at least three additional islands in previously published LFR data (33). The high incidence of mutations suggests that they are early mutations fixed in the tumor. Analysis of the DigiPico data from the tumor islands revealed the presence of a limited number of clonal-specific mutations that are absent in bulk tumors. Compared to the other islands, each of the five pre-chemotherapy islands had many truly unique mutations (2, 6, 8, 8, and 36), indicating that they were recent occurrences. Figure 3A Large-scale WGS data on PT2R recurrence revealed 3,009 novel somatic mutations not present in extensive pre-chemotherapy sequencing, DigiPico, or LFR data. These mutations could occur at any time, as the common ancestor of omental masses and PT2R recurrence has diverged from each other. Figure 3A Analysis of tumor islands in a recurrent sample from patient #11152 revealed a high burden of clonal-specific mutations in pelvic recurrent tumors (PT2R) compared to para-aortic lymph node recurrence (PALNR) or pre-chemotherapy tumors. This observation suggests that molecular mechanisms induced by SNV mutagenesis may have been recently activated in this patient. Figure 3A Furthermore, analysis of clonal-specific mutations in the PT2R samples using rainfall mapping revealed a strong subclonal local hypermutation (kataegis) event on chromosome 17 (8). Figure 3B , 3C (and 11). Comparison of the mutations constituting this kataegis event with bulk sequencing data of prechemotherapy omental masses, DigiPico data, and LFR data revealed that they were found only in DigiPico PT2R data, indicating that they are true clonal-specific mutations.

[0208] in conclusion

[0209] In this work, we use DigiPico / MutLX as an integrated platform to identify mutations from small cell populations with unprecedented accuracy at the whole-genome scale. We believe this work provides an important stepping stone for discovering current or recent somatic mutational processes occurring in cancer and normal tissues. Understanding current mutational processes is key to predicting the evolutionary trajectory of tumors and may be key to interfering with these trajectories in therapy. Mutations identified in high-volume sequencing of tumors must have occurred at some point in the extended history of the tumor from its occurrence to its presentation. In contrast, cell-specific mutations must have occurred within the finite lifespan of that cell. Similarly, mutations in small clones originating from single cells are also recent. The age of such mutations cannot exceed the age of the clone, which is defined by the number of cell divisions required to generate that clone. Studying patterns of cell-specific or small clone-specific mutations can identify recent or current mutational processes (1). Defining such processes is highly desirable because they may be causally related to biological or chemical phenomena and thus can yield important mechanistic insights. Identifying these mechanisms has important practical implications because they may be suitable for therapeutic interventions or predicting future tumor behavior. Current techniques do not allow for the direct and accurate identification of mutations from single cells or single small clones originating from tumors. DigiPico / MutLX has achieved this goal for the first time.

[0210] To overcome major technical limitations primarily associated with the discovery of false-positive mutations, current single-cell WGS analysis methods either require extensive validation studies (11) or rely on combining data from multiple cells to obtain reliable mutations shared between cells (12, 34). These cells are then grouped into clones derived from a common ancestor. While these techniques are applicable to more recent common ancestors compared to bulk sequencing, they remain unsatisfactory because the data obtained from these methods do not reflect the mutational processes occurring in existing cells. Furthermore, reducing the sequencing depth per cell to achieve sequencing of a large number of cells reduces coverage breadth, which is already compromised by the loss of genetic material during preparation steps. This increases the number of cells that need to be analyzed to infer and identify clones, further tracing ancestry back to the past. In addition, the lack of information on physical relevance in single-cell analysis methods leads to the loss of opportunities to group cells that may originate from a single clone. This increases the gap between the inferred ancestor of the clone and the present, making it difficult to define currently active intracellular processes in the tumor.

[0211] DigiPico / MutLX has a unique advantage in preserving spatial information. Analyzing spatially related cells preserves physical relevance and assumes that physically related cells belong to a single clone (9). It is also suggested to define different structures that may originate from tissue-resident stem cells to identify and analyze clones. For example, cells from a single small intestinal crypt or a single endometrial gland can reasonably be expected to originate from a single tissue-resident stem cell (35, 36). In these cases, each anatomical unit defines a clone that may or may not have clone-specific mutations that can be associated with mutation drivers. Furthermore, sequencing data from clones can be used computationally to infer subclones and predict more recent events that may occur within clones. This is similar to what is achieved with bulk sequencing and analysis, but at the level of a single clone composed of a limited number of cells. Preserving spatial information is also particularly interesting due to recent developments in spatial transcriptomics technologies (37). It is conceivable that combining highly accurate DNA sequencing with spatial transcriptomics could dissect genetic and non-genetic heterogeneity in tissues. In short, current techniques for analyzing small clones produce a large number of false positives, making it impossible to obtain direct and accurate clone-specific information at the genome scale without thorough validation. Combining data from multiple clones is a common solution, but it pushes ancestry further back. We previously used this approach to analyze a small number of tumor cells (tumor islands) (33). Due to the uncertainty associated with mutation detection from individual islands, it was necessary to identify only mutations common to all tumor islands and to effectively identify only major mutations. Approximately 700 mutations were then individually validated using targeted sequencing. While this still yielded important biological insights, we were unable to investigate island-specific mutations. DigiPico / MutLX now makes it possible to investigate such mutations. We demonstrate how direct analysis of the DNA from approximately 30 cancer cells led to the successful identification of a subclonal kataegis event.

[0212] In summary, we demonstrate that DigiPico and MutLX can identify somatic mutations with ultraaccurate precision from a limited number of cells obtained from clinical samples, representing a significant improvement over existing methods. Furthermore, unlike other computational methods that rely on diploid regions of the genome to calculate amplification bias, our method is also compatible with genomes undergoing extensive copy number alterations, such as in HGSOC. We believe the versatility of the DigiPico / MutLX method enables the study of active mutational processes in both tumor and normal tissues.

[0213] Availability

[0214] The source code for MutLX is available on GitHub. https: / / github.com / mmdknr / DigiPico ).

[0215] Registration number

[0216] All sequencing data used in this study were available on EGA (EGAD0001005118).

[0217] References

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

[0219] 2. Zhang, J., SS, Marjani, SL, Zhang, W. and Pan, X. (2018) Characterization of cancer genomic heterogeneity by next-generation sequencing advances precision medicine in cancer treatment. Precis. Clin. Med., 1, 29-48.

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

[0221] 4. Barber, LJ, Davies, MN and Gerlinger, M. (2015) Dissecting cancerevolution at the macro-heterogeneity and micro-heterogeneityscale. Curr. Opin. Genet. Dev., 30, 1-6.

[0222] 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.

[0223] 6.Chen,L.,Liu,P.,Evans,T.C.J.and Ettwiller,L.M.(2017)DNA damage isapervasive cause of sequencing errors,directly confounding variantidentification.Science,355,752-756.

[0224] 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.

[0225] 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.

[0226] 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.

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

[0228] 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-amplifiedsingle cells.Nat.Methods,14,491-493.

[0229] 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.

[0230] 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.

[0231] 14.Krueger F.(2016)Trim Galore!

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

[0233] 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 GenomeAnalysisToolkit:a MapReduce framework for analyzing next-generation DNAsequencing data.Genome Res.,20,1297-1303.

[0234] 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 ofgermline and somatic variants.Nat.Methods,15,591-594.

[0235] 18.Hosokawa,M.,Nishikawa,Y.,Kogawa,M.and Takeyama,H.(2017)Massivelyparallel whole genome amplification for single-cell sequencing usingdroplet microfluidics.Sci.Rep.,7,5199.

[0236] 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-genomesequencingand haplotyping from 10 to 20 human cells.Nature,487,190-195.

[0237] 20.Picard Tools(2018).

[0238] 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-basedapproaches for calling variants in clinical sequencingapplications.Nat.Genet.,46,912-918.

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

[0240] 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 formatand VCFtools.Bioinformatics,27,2156-2158.

[0241] 24.Chollet,F.and others(2015)Keras.

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

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

[0244] 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.

[0245] 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.

[0246] 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.

[0247] 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.

[0248] 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.AUAIPress.

[0249] 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.

[0250] 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.

[0251] 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 of40,000single cells identifies stochastic aneuploidies,genomereplication states and clonal repertoires.bioRxiv,10.1101 / 411058.

[0252] 35. Moore, L., Leongamornlert, D., Coorens, THH, Sanders, MA, Ellis, P., Dawson, K., Maura, F., Nangalia, J., Tarpey, PS, Brunner, SF, et al. (2018) Themutational landscape of normal human endometrial epithelium. bioRxiv, 10.1101 / 505685.

[0253] 36. Lee-Six, H., Ellis, P., Osborne, RJ, Sanders, MA, 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.

[0254] 37. Burgess, DJ (2019) Spatial transcriptomics coming of age. Nat. Rev. Genet., 20, 317.

[0255] All references mentioned in this article may be incorporated by way of citation.

[0256] DigiPico2

[0257] Example 2 - DigiPico2, a new method for whole-genome sequencing of picogram-sized DNA samples with unprecedented accuracy.

[0258] introduction

[0259] Previously, we described the DigiPico library preparation pipeline and the MutLX analysis platform as a method for accurately identifying single nucleotide variants (SNVs) from a limited amount of clinical material. This is a significant methodological advance, primarily because the limited amount of genetic material obtained from clinical samples must undergo whole-genome amplification (WGA) prior to sequencing. However, the WGA process introduces up to 100,000 artificial mutations into the amplified DNA, resulting in a final analysis riddled with false-positive variant detections, thus hindering any meaningful genetic interpretation of the original sample. In the DigiPico / MutLX strategy, we overcome this obstacle by isolating individual DNA molecules into independent compartments prior to the WGA step and indexing them afterward. By doing so, we digitize information about true mutations, meaning each compartment will carry or not carry the mutated allele. However, due to the way they are generated during the WGA process, artificial mutations will result in compartments containing both mutated and reference allele information (…). Figure 1B Based on this information, we subsequently developed an artificial neural network (ANN) based algorithm, MutLX, to effectively identify and eliminate these artificial mutations in our data (Figure 2). We extensively tested our strategy on simulated data and patient samples, demonstrating that our method is indeed effective in eliminating the detection of FP variants.

[0260] However, while generating high-quality data, the DigiPico library preparation method has few technical limitations. First, the fragmentation step of the library preparation (CoREF) borrowed from the aforementioned methods is very complex and time-consuming. Furthermore, CoREF requires the use of dUTP during the WGA process. Since dUTP is a non-natural nucleotide, it may introduce more artificial mutations into the final product. Next, we found that the adapter ligation efficiency in DigiPico is very low, which can sometimes affect library quality. Finally, due to the large number of indexes and the lack of redundancy in index information, there is a chance of index cross-contamination, which may adversely affect the final results. Therefore, we developed the DigiPico2 library preparation method to address all these issues.

[0261] result

[0262] Improved DigiPico library preparation workflow

[0263] As mentioned earlier, dUTP is used in the DigiPico method because the CoREF fragmentation process is required, which is a very complex fragmentation strategy. Figure 14ATherefore, using an alternative fragmentation method would address both the complexity and dUTP issues. To this end, we decided to use the existing fragmentation and end-repair strategy provided by the Lotus DNA Library Preparation Kit (IDT, USA). In the Lotus DNA Library Preparation Kit, the enzyme mixture is used to fragment large DNA molecules into smaller fragments in a time-dependent manner and prepare the ends of the fragments for the ligation step. To make this new strategy compatible with DigiPico, we initially used an I-DOT (Dispendix, Germany) dispenser to ensure that all compartments received the enzyme mixture more or less simultaneously. Next, we optimized the reaction conditions to achieve the fragment lengths required for DigiPico sequencing. By doing so, we were able to reduce the library preparation time from 12 hours to 4.5 hours and eliminate the need for dUTP in the WGA reaction. Figure 14B ).

[0264] Next, we aim to address the issue of low ligation efficiency in DigiPico. Initially, our linker ligation and indexing relied on an asymmetric ligation method. In this method, long indexed oligonucleotides with short complementary regions are used for ligation, which is extremely inefficient. Figure 15A A more efficient approach would require circular common adapters attached to both ends of the fragment. However, since these adapters do not contain any indexes in their standard form, the products need to be purified separately, and then the indexes introduced via PCR using index primers. This adds another challenge, as purifying 384 individual products would be extremely complex, time-consuming, and carries a significant risk of cross-contamination. To overcome these issues, we devised a novel indexing strategy for DigiPico2. In DigiPico2, a set of indexes is first introduced into the stem-loop of the universal adapter (…). Figure 15BIn the ligation step, all wells in each column of the plate receive a different indexed circular common adapter, so a total of 24 different oligonucleotides are sufficient to index all columns of the plate with the first set of indexes (column indexes). After the ligation step, all wells in each row are merged into a single tube, resulting in 16 different pools. These 16 different pools can be easily purified for use in the next indexing step. In the next step, the purified products from each pool are indexed by PCR using 16 different index primers (row indexes). At the end of this stage, each well of the plate will receive a different column-row index combination. Using this indexing strategy not only significantly improves ligation efficiency but also introduces two sets of redundancy that can be used to eliminate index cross-contamination in the data. In the first set, the column index is bound to both ends of each fragment. Therefore, any cross-contamination is most likely to result in fragments with different indices at their ends, which can then be easily removed from the data. In the second set, the indexed oligonucleotides in each row can be doubly indexed so that both standard index 1 (i7) and standard index 2 (i5) sequences can uniquely identify a particular row. By combining these redundant groups, the index cross-contamination rate can be reduced by at least two orders of magnitude.

[0265] DigiPico2 workflow significantly improves library quality.

[0266] To test the impact of these modifications on the final data quality, 120 pg of DNA from patient 11152's blood sample was sequenced using DigiPico2. This sample was used because we had previously performed extensive sequencing on both the patient's tumor and normal cells. As expected, the WGA was similar to previous versions, resulting in a very uniform product distribution. Figure 16A However, after library preparation and sequencing, it is clear that, unlike DigiPico, in DigiPico2, the performance of each well in the final library appears to be strongly correlated with the amount of WGA product. Figure 16A -C). This is likely a direct result of improved link efficiency. This improved correlation also allows for the introduction of QC measures based solely on the uniformity of WGA products, which was previously impossible. Furthermore, analyzing index redundancy information, we found that using an index cross-contamination filter eliminated nearly 5% of reads. Without this new filter, these contaminants could negatively impact the analysis results. Finally, we analyzed the data from DigiPico2 using the MutLX algorithm. The final results show a clearer distinction between artificial and real mutations, indicating that DigiPico2 performs at least as well as the DigiPico method when analyzed using the MutLX algorithm. Figure 16D ).

[0267] Extending the DigiPico2 workflow to single-cell whole-genome sequencing

[0268] After establishing the DigiPico2 workflow, we tested its applicability to single-cell whole-genome sequencing. This is important because active mutational processes are likely to begin within a single cell. Therefore, we introduced a single-cell DigiPico (ScDigiPico) sequencing workflow by partitioning the DNA of a single cell into an entire row of a 384-well plate. Figure 17A To evaluate the potency of the ScDigiPico method in identifying active mutational processes, we simulated such processes using Kuramochi cells cultured with N-ethyl-N-nitrosourea (ENU). ENU is an alkylating agent and a highly potent mutagen, preferentially causing T>C, T>A, and C>T mutations. In this setting, each cell was expected to acquire a different set of mutations, but since the underlying mechanism of mutagenesis was the same, the mutations were expected to belong to the same type. ScDigiPico was able to identify enrichments of the aforementioned mutation types. Figure 17B We further investigated the potential of ScDigiPico by irradiating Kuramochchi cells with ultraviolet light prior to single-cell sorting and ScDigiPico library preparation. Figure 17C The study identified the Kataegis event in one of the cells. These cumulative results suggest that ScDigiPico is an effective strategy for accurately identifying true mutations caused by active mutational processes within a single cancer cell.

[0269] DigiPico2 scheme

[0270] First, denature 200 pg of purified DNA, 20–30 resuspended cell nuclei, or laser-captured microdissectioned tumor islands using 5 μl D2 buffer from the Repli-g Single Cell Kit (Qiagen). After incubating at room temperature for 5 minutes, add 95 μl of water to the sample, and then use a Mosquito HTS liquid processor (TTP Labtech) to add 200 nmol of denatured template to each well of a 384-well reaction plate. Each well already contains 800 nmol of WGA mixture (0.58 μl Sc reaction buffer, 0.04 μl Sc polymerase (REPLI-g Single Cell Kit, Qiagen), 0.04 μl EvaGreen 20x (Biotium), and 0.065 μl water). Incubate the plate at 30 °C for 1.5 h, then heat-inactivate at 65 °C for 15 min. If desired, adding EvaGreen to the reaction allows for monitoring of the WGA reaction using a real-time PCR machine. Next, transfer 250 nmol of WGA reaction mixture to a new plate and add 1.1 nmol of NEBNext Ultra II FS reaction mixture (753 nmol water, 270 nmol Ultra II FS reaction buffer, and 77 nmol Ultra II FS enzyme mixture) to each well using an I-DOT dispenser (Dispendix). Incubate the plate at 37 °C for 6 min, then heat-inactivate at 65 °C for 30 min. Next, add 150 nmol of DigiPico indexed loop adapters with column indexes to all wells. Note that all wells in the same column will receive the same indexed oligonucleotides at this stage. Next, add 1.2 nmol of Ultra II ligation mixture (1150 nmol Ultra II ligation master mixture, 38 nmol ligation enhancer, and 12 nmol water) to each well using a Mosquito liquid processor, mix for 5 cycles, and incubate the plate at 20 °C for 15 min, then heat-inactivate at 65 °C for 10 min. Finally, merge all wells in the same row using a Mosquito liquid processor. Then, 1.5 μl of USER enzyme (NEB) was added to each 20 μl merged product, and the reaction mixture was incubated at 37 °C for 15 min. USER enzyme cleaves the circular linker at the uracil site. Next, SPRI beads were used to select the product size to achieve a range of 300–400 bp. Each product was then amplified for 4 cycles using row index primers. The final products were merged together, and the final library was purified using SPRI beads.

[0271] ScDigiPico solution

[0272] Single cells were sorted into the wells of the first column of a 384-well plate. Each well contained 4.5 μl of MyPK buffer. The plate was incubated at 55 °C for 30 min. Then, 900 nmol of stop solution was added to each well, and the plate was incubated at 95 °C for 5 min to inactivate proteinase K. The lysed cells were then distributed across rows using a Mosquito liquid processor, 200 nmol per well. Next, 800 nmol of WGA reaction mixture was added to each well, and WGA and library preparation were performed similarly to DigiPico2.

[0273] Indexed ring connector - column index

[0274] (Partial sequences from the library preparation manual) Multiplex Oligos for (for) of Multiple oligonucleotides (index primer set 1)

[0275] https: / / international.neb.com / - / media / nebus / files / manuals / manuale7335.pdf? rev=4bf1622b342b4d73a2b01443068ed2c5&hash=B049D91A18CDB471 AB388DC6E67E06B79263E5C5 )

[0276] P-[index']

[0277] AGATCGGAAGAGCACACGTCTGAACTUCCCTACACGACGCTCTTCCGATCT

[0278] [Index]*T(SEQ ID NO:1)

[0279] Where P is a 5' phosphate group, and * represents a thiophosphate bond. The index is either a column index (Ci) or a row index (Ri) sequence, serving as a unique barcode for each column or row, respectively.

[0280] Line index primers (oligonucleotide sequences) © 2007-2013 Illumina, Inc. All rights reserved.

[0281] P5:AATGATACGGCGACCACCGAGATCTACAC[r-index]

[0282] ACACTCTTTCCCTACACGACGCTCTTCCGATC*T(SEQ ID:NO:2)

[0283] P7:CAAGCAGAAGACGGCATACGAGAT[r-index]

[0284] GTGACTGGAGTTCAGACGTGTGCTCTTCCGATC*T(SEQ ID:NO:3)

[0285] * indicates a thiophosphate bond. sequence list <110> Oxford University Innovation Ltd. <120> Whole-genome sequencing methods for picogram-sized DNA <130> JDM104471P.WOP <150> GB1918043.9 <151> 2019-12-09 <160> twenty two <170> PatentIn version 3.5 <210> 1 <211> 52 <212> DNA <213> Artificial sequence <220> <223> Index ring connector <400> 1 agatcggaag agcacacgtc tgaactuccc tacacgacgc tcttccgatc tt 52 <210> 2 <211> 62 <212> DNA <213> Artificial sequence <220> <223> Line index primers <400> 2 aatgatacgg cgaccaccga gatctacaca cactctttcc ctacacgacg ctcttccgat 60 ct 62 <210> 3 <211> 58 <212> DNA <213> Artificial sequence <220> <223> Line index primers <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> Homo sapiens <400> 5 ccttatttct gatgctctta gat 23 <210> 6 <211> 23 <212> DNA <213> Homo sapiens <400> 6 gtatcagtca gccagaaaaa agg 23 <210> 7 <211> 42 <212> DNA <213> Homo sapiens <400> 7 tttattgaag tttgttttcc tctttgatcc taccactttt tt 42 <210> 8 <211> 42 <212> DNA <213> Homo sapiens <400> 8 aagccaatgt attgatcgct tccgagcaat aagaatagtg at 42 <210> 9 <​​​​​​​gttgttgttt gccaagctaa tctgcctggt tttatttata tc 42 <210> 10 <211> 42 <212> DNA <213> Homo sapiens <400> 10 aagtctacat taaacaatga tcacatctaa agctttatct tt 42 <210> 11 <211> 42 <212> DNA <213> Homo sapiens <400> 11 ctcatatata aagccttatt tctgatgctc ttagatttct ga 42 <210> 12 <211> 42 <212> DNA <213> Homo sapiens[[ID=�3]] <400> 12 gggactacag atgtgtgcca tcacacccag ctagtttttt gt 42 <210> 13 <211> 42 <212> DNA <213> Homo sapiens <400> 13 tgcataggta taggtatcag tcagccagaa aaaaggactt tg 42 <210> 14 <211> 42 <212> DNA <213> Homo sapiens <400> 14 gtatatatac aaatactttg tccatttaaa aaattaggtt at 42 <210> 15 <211> 42 <212> DNA <213> Homo sapiens <400> 15 cctatataga ctaacatgga tctaactttt tgactatctt cc 42 <210> 16 <211> 42 <212> DNA <213> Homo sapiens <400> 16 gaaatgcttt gtgaaatatg tcaacatact ggttgcaaat gc 42 <210> 17 <211> 42 <212> DNA <213> Homo sapiens <400> 17 gagtatggct atctatacct gccttttaag tttgaaacta ac 42 <210> 18 <211> 42 <212> DNA <213> Homo sapiens <400> 18 gtctttcctc tctctgtcct tccccgaaag tctactcggg tg 42 <210> 19 <211> 42 <212> DNA <213> Homo sapiens <400> 19 ggcatgatca ctgcagcctc tctgcttccc agattcaagt ga 42 <210> 20 <211> 42 <212> DNA <213> Homo sapiens <400> 20 cattaggggc tggacactca tcgagatgac ctgcctacaa at 42 <210> 21 <211> 42 <212> DNA <213> Homo sapiens <400> 21 gattgaaact gtccatttaa tctccttcct cccattatca at 42 <210> 22 <211> 42 <212> DNA <213> Homo sapiens <400> 22 tacctattta tctatatatt tcaacttata aaactttctt tc 42

Claims

1. A method of whole genome sequencing of a single cell or a cell population for identifying single nucleotide variants in the genome of the single cell or cell population, determining chromosomal structural variations in the genome of the single cell or cell population, or determining phasing information in the genome of the single cell or cell population, the method comprising: i) providing a multi-well array plate comprising a plurality of rows and a plurality of columns of reaction wells; ii) providing genomic DNA of a single cell or a cell population, wherein the genomic DNA is distributed in a plurality of reaction wells on the multi-well array plate such that there is no more than one single-stranded genomic DNA molecule of any given locus per reaction well, iii) performing whole genome amplification (WGA) on each genomic DNA molecule to provide multiple copies of the genomic DNA molecule in each reaction well; iv) fragmenting the DNA molecules of each reaction well and ligating a pair of circular adaptors or tag using transposase delivery adaptors to each end to form adapted DNA fragments, wherein the circular adaptors or transposase delivery adaptors comprise a column index (Ci) sequence or a row index (Ri) sequence, wherein the Ci sequence is common to each circular adaptor or transposase delivery adaptor for each reaction well in a column of the multi-well array plate, or wherein each Ri sequence is common to each circular adaptor or transposase delivery adaptor for each reaction well in a row of the multi-well array plate; v) providing an indexed DNA library by indexing PCR on the adapted DNA fragments, wherein the adapted DNA fragments are amplified using forward index primers and reverse index primers to form indexed PCR products, wherein a row index (Ri) sequence or a column index (Ci) sequence is introduced to each end of the adapted DNA fragments by each forward index primer and reverse index primer, such that the resulting indexed PCR products comprise both a pair of flanking column index (flanking Ci) sequences that are common to each well of a column and a pair of flanking row index (flanking Ri) sequences that are common to each well of a row; and vi) sequencing the indexed DNA library to provide data for determining any single nucleotide variants in the genome of the single cell or cell population, determining chromosomal structural variations in the genome of the single cell or cell population, or determining phasing information in the genome of the single cell or cell population.

2. The method of whole genome sequencing of a single cell or cell population of claim 1, wherein, The cell or cell population is from a tissue biopsy of a subject.

3. The method of whole genome sequencing of a single cell or cell population according to claim 1 or 2, wherein, The cell or cell population comprises cancer cells, precancerous cells, or suspected cancer cells, or a combination of cells thereof.

4. The method of whole genome sequencing of a single cell or cell population of claim 1, wherein, The genomic DNA comprises DNA of about 1 to 30 cells.

5. The method of whole genome sequencing of a single cell or cell population of claim 1, wherein, The DNA content of a single cell is distributed in wells of a single row; or The DNA content of a single cell or cell population is distributed in wells of both a row and a column of a single multi-well array plate.

6. The method of whole genome sequencing of a single cell or cell population of claim 1, wherein, The multi-well array plate comprises a 384-well plate.

7. The method of whole genome sequencing of a single cell or cell population of claim 1, wherein, A DNA polymerization reporter molecule is provided in the amplification mix.

8. The method of whole genome sequencing of a single cell or cell population of claim 1, wherein, The circular adaptors are provided such that the method comprises a step of fragmenting the DNA molecules of each reaction well and a subsequent ligation reaction to ligate the circular adaptors to the fragmented DNA; or wherein the transposase delivery adaptor is provided such that the method comprises fragmenting the DNA molecules by a tagmentation process.

9. The method of whole genome sequencing of a single cell or cell population of claim 1, wherein, Fragmenting the DNA molecules in each reaction well into a plurality of dsDNA fragments comprises direct fragmentation by an enzyme.

10. The method of whole genome sequencing of a single cell or cell population of claim 1, wherein, The fragmentation or tagmentation reagents are added to each well simultaneously.

11. The method of whole genome sequencing of a single cell or cell population of claim 1, wherein, After the DNA is fragmented to form DNA fragments, the DNA fragments are end-repaired and / or dA-tailed such that they can be ligated to the circular adaptor.

12. The method of whole genome sequencing of a single cell or cell population of claim 1, wherein, The circular adaptor comprises an oligonucleotide having a secondary stem-loop structure, and wherein the circular adaptor comprises a pair of complementary sequence regions flanking the loop region, wherein the pair of complementary sequence regions are arranged to hybridize to each other to form the stem-loop structure of the circular adaptor.

13. The method of whole genome sequencing of a single cell or cell population of claim 1, wherein, The ends of the adapted DNA fragments are symmetrical.

14. The method of whole genome sequencing of a single cell or cell population of claim 1, wherein, The circular adaptor comprises a uracil in the loop region, and after ligation of the circular adaptor, a single-stranded region of the circular DNA is cleaved at the uracil.

15. The method of whole genome sequencing of a single cell or cell population of claim 1, wherein, Ci sequences are provided in the adapted DNA fragments, and the method can additionally comprise the step of pooling the adapted DNA fragments from each reaction well in a row prior to the indexing PCR; or wherein Ri sequences are provided in the adapted DNA fragments, and the method additionally comprises the step of pooling the adapted DNA fragments from each reaction well in a column prior to the indexing PCR.

16. The method of whole genome sequencing of a single cell or cell population of claim 1, wherein, The adapted DNA fragments comprise Ci sequences, and the forward and reverse indexing PCR primers each comprise Ri sequences for providing a pair of Ri sequences in the indexed PCR products.

17. The method of whole genome sequencing of a single cell or cell population of claim 1, wherein, The forward and reverse indexing PCR primers further comprise sequencing adaptor sequences such that sequencing adaptors are bound to the indexed PCR products.

18. The method of whole genome sequencing of a single cell or cell population of claim 1, wherein, The indexed DNA fragments of the indexed DNA library are size filtered.

19. The method of whole genome sequencing of a single cell or cell population of claim 1, wherein, The method comprises determining any true SNVs in the genome of the single cell or cell group by determining whether substantially all indexed DNA library sequences originating from a single well comprise the same SNV, or whether only a fraction of the indexed DNA library sequences comprise the same SNV, wherein SNVs shown in substantially all indexed DNA library sequences originating from a single well are determined to be true SNVs in the genomic DNA, and SNVs found in only a fraction of the indexed DNA library sequences originating from a single well are determined to be false positive (FP) SNVs.

20. The method of whole genome sequencing of a single cell or cell population of claim 1, wherein, The method further comprises matching indexed DNA library sequences originating from a single well representing one strand of the genomic DNA to indexed DNA library sequences originating from another well representing the complementary strand of the genomic DNA, wherein SNVs present in substantially all indexed DNA library sequences of both complementary strands of the genomic DNA are determined to be true SNVs, and SNVs not present in substantially all indexed DNA library sequences of both complementary strands of the genomic DNA are determined to be false positives.

21. The method of whole genome sequencing of a single cell or cell population according to claim 19 or 20, wherein, The determination is performed in silico using BAM file data generated from mapping the sequence data to a reference genome.

22. The method of whole genome sequencing of a single cell or cell population of claim 21, wherein, The computer-determined or matched and / or the calculation of the probability score of the indexed DNA sequences is performed by an artificial neural network (ANN) model, optionally by a multilayer perceptron.

23. The method of whole genome sequencing of a single cell or cell population of claim 1, wherein, The method prepares an indexed DNA library from tumor cells, suspected tumor cells, or precancerous cells, and non-cancerous cells, and wherein sequencing data from the tumor cells, suspected tumor cells, or precancerous cells is compared to sequencing data obtained from the non-cancerous cells taken from non-cancerous tissue as a control.

24. The method of whole genome sequencing of a single cell or cell population of claim 1, wherein, The probability score of a particular nucleotide variant being a true SNV or false positive is calculated in a computer, thereby determining that a given variant nucleotide has a statistically significant probability of being a true SNV or false positive.

25. The method of whole genome sequencing of a single cell or cell-group according to claim 1, wherein, The sequence data is provided in the form of paired-read FastQ files.

26. The method of whole genome sequencing of a single cell or cell population of claim 25, wherein, The sequence data of the paired-read FastQ files is trimmed to remove adapter sequences and for quality, thereby providing trimmed data.

27. A method of preparing an indexed DNA library for sequencing of nucleic acid molecules, the method comprising: i) providing a multi-well array plate comprising a plurality of rows and a plurality of columns of reaction wells; ii) providing nucleic acid molecules, wherein the nucleic acid molecules are distributed among a plurality of reaction wells on the multi-well array plate such that there is no more than one single-stranded nucleic acid molecule of any given position per reaction well, iii) amplifying the nucleic acid molecules to provide a plurality of DNA copies of the nucleic acid molecules in each reaction well; iv) fragmenting the DNA molecules of each reaction well and ligating a pair of circular adapters or tagmentation with transposase to each end to form adapted DNA fragments, wherein the circular adapters or transposase-delivered adapters comprise a column index (Ci) sequence or a row index (Ri) sequence, wherein the Ci sequence is common to each circular adapter or transposase-delivered adapter of each reaction well in a column of the multi-well array plate, or wherein each Ri sequence is common to each circular adapter or transposase-delivered adapter of each reaction well in a row of the multi-well array plate; v) providing the indexed DNA library by indexing PCR of the adapted DNA fragments, wherein the adapted DNA fragments are amplified using forward index primers and reverse index primers to form indexed PCR products, wherein a row index (Ri) sequence or a column index (Ci) sequence is introduced to each end of the adapted DNA fragments by each forward index primer and reverse index primer, such that the resulting indexed PCR products comprise both a pair of flanking column index (flanking Ci) sequences that are common to each well of a column and a pair of flanking row index (flanking Ri) sequences that are common to each well of a row; and optionally, wherein the forward index primers and reverse index primers further provide respective 5’ and 3’ sequencing adapters onto the indexed PCR products suitable for a sequencing reaction.

28. A method of preparing an indexed DNA library for whole genome sequencing of a single cell or cell population for identifying single nucleotide variants in the genome of the single cell or cell population, determining chromosomal structural variants in the genome of the single cell or cell population, or determining phasing information in the genome of the single cell or cell population, the method comprising: i) providing a multi-well array plate comprising a plurality of rows and a plurality of columns of reaction wells; ii) providing genomic DNA of a single cell or cell population, wherein the genomic DNA is distributed among a plurality of reaction wells on the multi-well array plate such that there is no more than one single-stranded genomic DNA molecule of any given locus per reaction well, iii) performing whole genome amplification (WGA) on each genomic DNA molecule to provide multiple copies of the genomic DNA molecule in each reaction well; iv) fragmenting the DNA molecules of each reaction well and ligating a pair of circular adaptors or tag using transposase delivery adaptors to each end to form adapted DNA fragments, wherein the circular adaptors or transposase delivery adaptors comprise a column index (Ci) sequence or a row index (Ri) sequence, wherein the Ci sequence is universal for each circular adaptor or transposase delivery adaptor of each reaction well in a column of the multi-well array plate, or wherein each Ri sequence is universal for each circular adaptor or transposase delivery adaptor of each reaction well in a row of the multi-well array plate; v) providing an indexed DNA library by performing index PCR on the adapted DNA fragments, wherein the adapted DNA fragments are amplified using forward index primers and reverse index primers to form indexed PCR products, wherein a row index (Ri) sequence or a column index (Ci) sequence is introduced to each end of the adapted DNA fragments by each forward index primer and reverse index primer, such that the resulting indexed PCR products comprise both a pair of flanking column index (flanking Ci) sequences that are universal for each well of a column and a pair of flanking row index (flanking Ri) sequences that are universal for each well of a row; and optionally, wherein the forward index primers and reverse index primers further provide respective 5’ and 3’ sequencing adaptors onto the indexed PCR products suitable for use in a sequencing reaction.

29. A method of whole genome sequencing of a single cell or cell population for identifying single nucleotide variants (SNVs) in the genome of the single cell or cell population, determining chromosomal structural variants in the genome of the single cell or cell population, or determining phasing information in the genome of the single cell or cell population, the method comprising: i) preparing an indexed DNA library by performing the method of claim 27 or 28, or providing an indexed DNA library prepared according to the method of claim 27 or 28; ii) sequencing the indexed DNA library to provide data for determining any single nucleotide variants (SNVs) in the genome of the single cell or cell group, determining chromosomal structural variations in the genome of the single cell or cell group, or determining phasing information in the genome of the single cell or cell group. ii) sequencing the indexed DNA library to provide data for determining any single nucleotide variants (SNVs) in the genome of the single cell or cell group, determining chromosomal structural variations in the genome of the single cell or cell group, or determining phasing information in the genome of the single cell or cell group. ii) sequencing the indexed DNA library to provide data for determining any single nucleotide variants (SNVs) in the genome of the single cell or cell group, determining chrom

Images