Patent Application
20200402613
This application claims priority from GB1803596.4 filed 6 Mar. 2018 and from GB1819134.6 filed 23 Nov. 2018, the contents and elements of which are herein incorporated by reference for all purposes.
The present invention relates in part to methods for detecting the presence of variant DNA, such as circulating tumour DNA (ctDNA) from, e.g., a cell-free DNA (cfDNA) source, such as blood plasma, or for detecting variant DNA in forensic applications, in pathogen identification, in agricultural and environmental monitoring of species contamination. In particular, the methods of the invention find use in the diagnosis, treatment and especially monitoring of cancer, including monitoring that is done following tumour resection.
The work leading to this invention has received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement no 337905.
Cell-free DNA (cfDNA), such as circulating tumour DNA (ctDNA) is increasingly being used as a non-invasive tool to monitor disease burden, response to treatment, and risk of relapse1,2. Following treatment, patients may have low ctDNA levels and even in advanced disease, concentrations can be lower than a few copies per sample volume3. In such cases, an individual sample may contain less than one detectable copy of a given mutation due to sampling statistics, resulting in undetected ctDNA even if its average concentration is non-zero: i.e. a false-negative underestimate of the ctDNA level1,3,4.
Next-generation sequencing (NGS) offers the possibility to analyse a large number of mutations in plasma in a single reaction. This has been demonstrated through amplicon-based5,6 and hybrid-capture methods for targeted sequencing7-9, using either standardised panels5,9 or bespoke panels covering regions that are specific to each patient5-7. These approaches have generally been applied to screen or monitor individual mutations. Despite targeting ˜20 patient-specific loci, a recent study detected ctDNA in <50% of patients with early-stage NSCLC and did not detect ctDNA immediately post-surgery in most patients who later relapsed6; this suggests that greater sensitivity would be required to effectively achieve this important clinical goal. The use of highly multiplexed capture panels that cover thousands of mutations has been suggested1,7, but this has not so far been demonstrated for analysis of ctDNA. These approaches for ctDNA analysis relied on identification of individual mutations across panels of variable sizes.
The detection of individual mutations is limited by both sampling error and sequencing background noise; when signals do not reach a pre-specified threshold for mutation calling, the information in these signals is lost.
Pécuchet et al., Clin. Chem., 2016, Vol. 62, No. 11, pp. 1492-1503, describes analysis of base-position error rate of next-generation sequencing to detect tumour mutations in circulating DNA. WO2016/009224 describes a method for detecting a genetic variant. WO2015/164432 describes a method for detecting mutations and ploidy in chromosomal segments. WO2013/138510 describes measurement of nucleic acid variants using highly-multiplexed error-suppressed deep sequencing. Ahn et al., Scientific Reports, 2017, 7:46678|DOI: 10.1038/srep46678 describes asymmetrical barcode adapter-assisted recovery of duplicate reads and error correction strategy to detect rare mutations in circulating tumor DNA. Kockan et al., Bioinformatics, 2017, Vol. 33, No. 1, pp. 26-34, describes ultrasensitive detection of single nucleotide variants and indels in circulating tumour DNA. WO2014/039556 describes systems and method to detect rare mutations and copy number variation. These references generally concern methods for reducing background noise rates of sequencing, including by use of Unique molecular identifiers (UMIs).
Newman et al., 2016 describe improvements to the CAPP-Seq method for detecting ctDNA in which integrated digital error suppression is employed (iDES CAPP-Seq)7. However, the iDES CAPP-Seq method involves the use of position-specific error rate for error correction. This requires determination of the error rate for each locus, which in turn requires a least 1/(position-specific error rate) molecules to be targeted at every locus to be interrogated. There remains an unmet need for methods of ctDNA detection that reduce the number of samples required to be analysed in order to carry out error suppression.
While detection of ctDNA shows promise in the field of cancer care, there remains an unmet need for methods and systems that maximise signal-to-noise ratio in the context of low ctDNA fraction. The present invention seeks to provide solutions to these needs and provides further related advantages.
The present inventors hypothesised that by integrating signals across a large number of mutation loci, it would be possible to mitigate the effects of sampling noise and obtain a more sensitive and accurate estimate of ctDNA levels, even when ctDNA was present at very low concentrations (
To use ctDNA information more efficiently, the present inventors circumvented the “calling” of individual mutations, and aimed to combine the information from mutant reads across multiple, e.g., all, the tumour-mutated loci. The present inventors found that by generating and combining a large number of sequencing reads from plasma DNA covering multiple loci that are mutated in a patient's tumour, it is possible to achieve detection that surpasses the sensitivity of previous methods. The inventors developed an algorithm, called INtegration of VAriant Reads (INVAR), that aggregates mutant signal across hundreds or thousands of mutation loci, to assess whether overall genome-wide signal is significantly above background, or non-distinguishable from background (
In a further optimisation of the INVAR approach, integration may be targeted to focus on the integration of residual disease signal. In particular, a focussed INVAR approach, described herein, aggregates minimal residual disease (MRD) ‘MRD-like signal’ by selecting signal from loci with up to 2 mutant molecules only. Secondly, only molecules with a mutation supported by the forward and reverse (F+R) read are considered for contributing to signal, which constitutes both an error-suppression and size-selection step. Thirdly, mutant reads per locus are weighted based on their tumour allele fraction to up-weight the mutations that are more prevalent in the tumour. Fourthly, the signal is then aggregated—in some cases by trinucleotide context. Fifthly, P-values are integrated using a suitable method (e.g. Fisher's method or Brown's method), but only across the top N classes in order to focus on MRD-like signal
The end result is a focused INVAR algorithm that is optimised for detection of residual disease.
Accordingly, in a first aspect the present invention provides a method (optionally a computer-implemented method) for detecting and/or quantifying cell-free DNA (cfDNA), such as circulating tumour DNA (ctDNA), in a DNA-containing sample obtained from a patient, the method comprising:
In certain cases, calculating mutant allele fraction may comprise calculating the weighted mean of the allele fractions at each of the patient-specific loci. In certain cases, calculating mutant allele fraction may comprise counting the number of mutant reads and comparing this to a pre-determined threshold. The pre-determined threshold may in some cases be a function of sequencing depth, but need not be a simple sum. In particular, a threshold model on the number of mutant reads may be applied.
Step (c) may be considered optional because its function is to reduce noise which may not be necessary in certain cases. In some embodiments a noise reduction step may be employed, which may for example comprise reads collapsing. In certain embodiments a noise reduction step may be omitted. In particular, where confidence arises from other mechanisms (such as replicates, use of classes, etc.), or as a result of improvements in sequencing quality, which may arise in the future. In particular, where step (c) is performed, reads collapsing may be as defined further herein. As used herein there terms “reads collapsing” and “read collapsing” are intended to be interchangeable.
Step (d) may be considered optional and/or may be carried out in different ways in different embodiments. In particular, in some embodiments step (d) is performed in order to compare the mutant fraction to background. This need not include the step of calculating the sum of mutant divided by the sum of total reads.
Some key statistical tests which do not use that calculation, but instead use only the number of mutant reads and the number of total reads, without dividing them. In some embodiments the method comprises interrogation of the individual at each mutant read and assess the background rate of that individual to compare with it.
In this embodiment it is not necessary to calculate the sum of the mutant reads across all loci.
In some embodiments, the method further comprises:
In some embodiments, the method comprises quantifying the concentration or amount of cfDNA (e.g., ctDNA) in the sample obtained from the patient, wherein quantifying the concentration or amount of cfDNA (e.g., ctDNA) comprises subtracting the background sequencing error rate from the mutant allele fraction calculated in step (d). In some embodiments, the calculation of Fisher's exact test may be independent of said step (d).
As described herein, differences were observed in the background sequencing error rate per class of mutation, i.e. the error rate of different single nucleotide substitutions differs (e.g. see
In some embodiments, the background sequencing error rate is or has been determined for each class of mutation (e.g., each class of base substitution) (“mutation class”) represented in said at least 2, 3, 4, 5, 6, 7, 8, 9, 10 or more patient-specific loci, and the mutant allele fraction calculation in step (d) is performed for each mutation class, taking into account the background sequencing error rate of that mutation class; the mutant allele fraction of each class is combined to provide a measure of the global mutant allele fraction of the sample. In particular, the global mutant allele fraction may be calculated as the mean of all of the individual per-class background-subtracted mutant allele fractions, weighted by the total number of read families observed in that class. In certain embodiments, particularly where the number of mutant and non-mutant reads is used to determine the presence of cfDNA without determining the mutant allele fraction, the calculation step (d) may be omitted.
In some embodiments, the method comprises making a determination of the statistical significance or otherwise of the calculated mutant allele fraction, taking into account the background sequencing error rate. In cases where the mutant allele fraction is calculated per mutation class and then combined into a global mutant allele fraction, the determination of statistical significance of the calculated global mutant allele fraction, may comprise determining the individual statistical significance of the mutant allele fraction of each mutation class and then combining the individual statistical significance determinations into a global statistical significance determination for the global mutant allele fraction. Various statistical methods may be suitable for the determination of statistical significance of the mutant allele fraction. In particular cases, the determination of the statistical significance of the mutant allele fraction may comprise carrying out a one-sided Fisher's exact test, given a contingency table comprising: the number of mutant reads from the sample, the total number of reads from the sample, and the number of mutant reads expected from the background sequencing error rate. In certain embodiments where mutant allele fraction is calculated on a per mutation class basis, the determination of mutant allele statistical significance may comprise carrying out multiple one-sided Fisher's exact tests to determine the statistical significance of the number of mutant reads observed given the background sequencing error rate for that mutation class, thereby generating a p-value for each mutation class, and combining the p-values using the Empirical Brown's method to provide a global measure of statistical significance for the mutant allele fraction of the sample.
When mutant allele fraction is calculated on a per mutation class basis, the number of mutation classes will generally be governed by the mutations found to be present in the at least 2, 3, 4, 5, 6, 7, 8, 9, 10, 100, 1000 or at least 5000 mutation-containing loci representative of a tumour of the patient (“patient-specific loci”). For many cases, the mutation classes may comprise at least 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 or all 12 of the following mutation classes: C>G, G>C, T>G, A>C, C>A, G>T, T>C, A>G, T>A, A>T, C>T and T>C. Preferably, the tumour-specific mutations at the patient-specific loci include mutations belonging to at least 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 or 12 different mutation classes. Further mutation classes are contemplated herein. For example, mutations may be split based on a greater number of sequence subsets such as by di-nucleotide context, tri-nucleotide context or by individual locus, which may further improve resolution of error rates.
As described herein (see Example 8 and
In some embodiments, the sequence data comprising sequence reads obtained in step (b) represent Tailored Panel Sequencing (TAPAS) sequence reads, focussed-exome sequence reads, whole-exome sequence reads or whole-genome sequence reads. The choice of sequence reads may reflect, inter alia, the mutation rate of the cancer being studied. Tumour-derived mutations can be identified using exome sequencing as demonstrated herein, but also across smaller focused panels or larger scales such as whole genome. In examples described herein where the patients had melanoma, exome sequencing was sufficient to identify hundreds to thousands of mutations per patient. Based on known mutation rates of cancer types, exome sequencing would also suffice for many cancer types with relatively high mutation rates, for example: lung, bladder, oesophageal, or colorectal cancers. For cancers with a mutation rate of ˜1 per megabase or less, whole-genome sequencing of tumours for mutation profiling may be desirable. For ovarian and brain cancers, this would result in thousands of mutations identified per patient. Moreover, the sequence data comprising sequence reads may cover a sufficient portion of the exome or genome of the sequence tumour to identify at least 2, 3, 4, 5, 6, 7, 8, 9, 10, 50, 100, 500, 1000, 2500 or at least 5000 mutation-containing loci. Additionally or alternatively, the sequence data comprising sequence reads may cover a sufficient portion of the exome or genome of the sequence tumour to ensure that the tumour-specific mutations at the patient-specific loci include mutations belonging to at least 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 or 12 different mutation classes. Additionally or alternatively, the sequence data comprising sequence reads may cover a sufficient portion of the exome or genome of the tumour to ensure that the tumour-specific mutations at the patient-specific loci include at least 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63 or at least 64 trinucleotide contexts, in particular, trinucleotide contexts selected from the group consisting of: CGC, GGC, TCG, ACG, GCG, TGC, CCG, GCA, CGA, GCC, CGG, CGT, AGC, GCT, TCA, TGA, AGT, ACC, CCC, CCA, CTT, GGG, CCT, GAG, CTG, AGG, CAG, CTC, AGA, TCC, GGT, TGG, CTA, ACA, TCT, TAG, AAG, TGT, ACT, GTC, GGA, TAC, TTG, CAA, TTC, TTA, ATC, ATG, TAA, TAT, CAT, GTT, ATT, ATA, GAA, GAC, GAT, CAC, GTG, TTT, GTA, AAT, AAA and AAC.
In some embodiments, the 2, 3, 4, 5, 6, 7, 8, 9, 10, 50, 100, 500, 1000, 2500 or at least 5000 mutation-containing loci representative of the tumour of the patient are obtained by sequencing DNA obtained directly from a tumour sample from the patient or sequencing DNA obtained from a plasma sample from the patient at a time of high tumour disease burden (e.g. prior to the start of therapeutic treatment or prior to surgical resection). In this way the determination of the tumour sequence, e.g., tumour exome or portion thereof or tumour genome or portion thereof, can be made using a relatively more abundant source of tumour-derived DNA and then the information about which loci contain tumour-specific mutations (step (a)) can be employed in the methods of the present invention practiced on sequence reads (step (b) obtained at a time when tumour-derived DNA is more scarce (e.g. after the patient has received at least one course of treatment and/or after surgical tumour resection)). For example, the method may be used to monitor recurrence of the tumour by detecting low levels of ctDNA. The determination of loci of interest comprising the 2, 3, 4, 5, 6, 7, 8, 9, 10, 50, 100, 500, 1000, 2500 or at least 5000 mutation-containing loci representative of the tumour of the patient will generally involve comparison with germline DNA sequencing of the patient so as to identify which loci contain mutations that are specific to the tumour relative to, or compared with, the patient's germline genomic sequence. For example, DNA extracted from buffy coat or any other suitable source of germline DNA (e.g. saliva, hair follicles, skin, cheek swab, white blood cells).
In some embodiments, the loci of interest are filtered by removing loci known to be single nucleotide polymorphisms (SNPs), e.g., by removing those positions found in common SNP databases (such as 1000 Genomes ALL or EUR). This filtering focuses on signal, i.e. tumour mutated loci, by excluding those loci that may be SNPs (see Example 10 herein).
In some embodiments, the sequence data comprising sequence reads provided in step (b) represent sequence reads of a plurality of DNA fragments from a substantially cell-free plasma sample from the patient. In some embodiments, the sequence data comprising sequence reads provided in step (b) represent sequence reads of a plurality of DNA fragments from any of the sample types as defined herein. The use of cell-free DNA (cfDNA) as a sample source provides a relatively non-invasive method for obtaining sample (so called “liquid biopsy”). The sequence reads obtained from cfDNA will comprise sequence reads of both that fraction of circulating DNA fragments that has its origins from the tumour or tumours of the patient (ctDNA fraction), if present, and that fraction of circulating DNA fragments that has its origins from non-tumour tissue or cells.
In some embodiments, the sequence data comprising sequence reads obtained in step (b) represent sequence reads of a plurality of polynucleotide fragments from a sample obtained from the patient after the patient has begun a course of treatment of the tumour and/or after the patient has had surgical resection of the tumour, and wherein the method is for monitoring the presence, growth, treatment response, or recurrence of the tumour. In particular embodiments, the method is for monitoring the presence and/or recurrence of minimal residual disease (MRD).
In accordance with this and other aspects of the present invention, the patient may be a patient who has, or has had, a cancer selected from melanoma, lung cancer, bladder cancer, oesophageal cancer, colorectal cancers, ovarian cancer brain cancer, and/or breast cancer. In particular, the patient may have been diagnosed as having melanoma, including advanced and/or invasive melanoma with or without metastases.
In some embodiments, the reads collapsing step (c) comprises the grouping of duplicate sequencing reads into read families based on fragment start and end positions and at least one molecular barcode, which uniquely label individual starting cfDNA molecules. As defined further herein, “barcode” or “molecular barcode” as used herein means a unique string of bases, generally of length <20, such as <10 bp, that may be ligated to DNA molecules as the first step during library preparation. As a result, read families may be uniquely identified, and thus linked to their starting molecule. This allows error-suppression via “reads collapsing”. Therefore, duplicate reads with the same start and end positions and molecular barcode can be identified computationally as having originated from the same starting cfDNA molecule, termed as a ‘read family’. In particular, a minimum 60%, 70%, 75%, 80%, 85%, 90% or even a 95% consensus (‘consensus threshold’) may be required between all family members for a read to be included in a read family. Thus, for example, where there are three reads in a read family and two of those reads show consensus while one shows, e.g., an alternative base a given nucleotide position, the read family would have a resulting consensus of ⅔ or 66%. In the case where there is a mutation, but the mutant base is not supported by a consensus greater than, or equal to, the consensus threshold for inclusion in a read family may be discarded (i.e. not used further in the analysis). In particular cases, a minimum family size of 2, 3, 4 or 5 reads may be required. In some cases, read families not satisfying this minimum family size may be disregarded in the analysis. The greater the family size, the greater the extent of error-suppression, because the consensus across the read family is supported by a larger number of independent reads. Therefore, in order to set a limit for the error-suppression step, it may be advantageous to specify a particular minimum family size threshold.
As described herein, the present inventors found that an in silico size selection, even with relaxed settings, was able to enrich for mutant signal (i.e. ctDNA) while minimising loss of rare mutant alleles. The enrichment was in some cases greater for lower initial allele fractions (see
Non-tumour DNA has been observed to peak at 166 bp, thus in some aspects size-selection windows may be adjusted to exclude or minimise non-tumour DNAs of length proximal to this maxima. Also contemplated herein are one or more narrower size windows for the size selection that would be expected to result in greater enrichment. For example, size ranges of 120-155 bp, 120-180 bp, 260-390 bp and/or 445-455 may be employed. Alternatively, the size selection may be less stringent with wider size selection windows such as 110-200 bp, 240-410 bp and/or 430-470 bp. In some embodiments, the in silico size selection may size select to one or more (e.g. 2 or 3) size windows that have been predetermined based on experimentally-determined size windows that enrich for ctDNA in the sample(s) in question. For example, the sequence reads from one or more samples may be combined, the size distribution of fragments determined, and the ratio between the proportion of mutant, and wild-type (i.e. germline sequence) reads determined. The size windows for the methods of the present invention may be those that display enrichment in the proportion of mutant reads relative to wild-type reads.
In certain embodiments, one or more filters are applied to the read families in order to focus on those families more likely to be tumour-derived. In some cases, the one or more filters may be minimal residual disease (MRD) filters, such as those described in Example 10 herein. In particular, a filter step may comprise excluding those loci with >2 mutant molecules. Alternatively or additionally, a filter step may comprise selecting (i.e. including) only those fragments which have been sequence in both forward (F) and reverse (R) direction. As described in Example 10, the requirement that mutant reads only be considered as contributing to signal at a locus if there is at least one F and at least one R read at the locus serves a dual purpose of suppressing sequencing artifacts, and selecting for mutant reads from short cfDNA fragments (supported by reads in both directions), which are slightly enriched in ctDNA (
In certain embodiments, a tumour allele fraction weighting is applied in order to increase the weighting (up-weight) the signal applied from mutations that are more prevalent in the tumour. As described in Example 11 herein, the present inventors found that the likelihood of observing a given mutation in cfDNA from plasma is proportional to the tumour allele fraction for the given mutation in the tumour (see
wherein:
AFcontext is the allele frequency of a given (e.g. trinucleotide) context; tumourAF is the allele frequency of the locus as determined by analysis of the tumour (e.g by sequencing DNA obtained directly from the tumour); and MRD-like loci are the mutation-containing loci determined from the tumour of the patient and which have been filtered to select for minimal residual disease signal. The effect of weighting by tumour allele fraction can be seen in Example 11, particularly comparing
In certain embodiments, the p-value for each trinucleotide context is determined by comparing samples against background error rates. The top (i.e. most significant) n p-values from trinucleotide contexts are then combined using a suitable technique, such as Fisher's method or Brown's method. In some cases n may be 2, 3, 4, 5, 6, 7, 8, 10, 12 or more. For example, where n=6 p-values from the top 6 trinucleotide contexts may be combined according to the formula:
In certain embodiments, the global allele fraction, AFglobal, is calculated based on all signal in all contexts, taking into account the background error, E. Preferably, AFglobal is determined according to the formula:
In a second aspect, the present invention provides a method for monitoring the presence of, growth of, prognosis of, regression of, treatment response of, or recurrence of a cancer in a patient, the method comprising:
In some cases, the sequencing step (i) may comprise Next-generation Sequencing (NGS), including Illumina® sequencing, or Sanger sequencing. NGS offers the speed and accuracy required to detect mutations, either through whole-genome sequencing (WGS) or by focusing on specific regions or genes using whole-exome sequencing (WES) or targeted gene sequencing. Examples of NGS techniques include methods employing sequencing by synthesis, sequencing by hybridisation, sequencing by ligation, pyrosequencing, nanopore sequencing, or electrochemical sequencing.
In some cases, the method of this aspect of the present invention further comprises a step, prior to sequencing, of preparing a DNA library from a sample (e.g. a plasma sample) obtained from the patient or from more than one patient. Optionally, the library may be barcoded.
In some cases, the method of this aspect of the present invention further comprises a step prior to sequencing of obtaining a sample from the patient. For example, a blood sample may be collected from a patient who has been diagnosed as having, or being likely to have, a cancer. The sample may be subjected to one or more extraction or purification steps, such as centrifugation, in order to obtain substantially cell-free DNA source (e.g. to obtain a plasma sample). The method may further comprise determining the cfDNA concentration of the sample. It is specifically contemplated that the sample may be transported and/or stored (optionally after freezing). The sample collection may take place at a location remote from the sequencing location and/or the computer-implemented method steps may take place at a location remote from the sample collection location and/or remote from the sequencing location (e.g. the computer-implemented method steps may be performed by means of a networked computer, such as by means of a “cloud” provider). Nevertheless, the entire method may in some cases be performed at single location, which may be advantageous for “on-site” determination or monitoring of cancer.
In some cases, the method of this aspect of the invention may further comprise obtaining tumour imaging data and/or measuring or detecting one or more tumour biomarkers to assist with the determination of presence of, growth of, treatment response of, or recurrence of the cancer. In particular, the tumour imaging data may comprise computed tomography (CT) data, e.g. to measure tumour volume. In particular, cases, the biomarker may comprise lactate dehydrogenase (LDH) concentration. Such additional means of tumour detection and/or quantification may corroborate a determination made by the methods of the present invention and/or may assist with resolving ambiguous determinations.
In some cases, the method of this aspect may further comprise a step of recommending or selecting the patient for anti-cancer treatment, including follow-on or continuing treatment(s). For example, where the sample is determined to contain ctDNA (e.g. where the mutant allele fraction is found to be greater, including statistically significantly greater, than the background sequencing error rate), the patient may be determined to have, or to have a recurrence of, a cancer which may benefit from anti-cancer treatment, including chemotherapy, immunotherapy, radiotherapy, surgery or a combination thereof. Likewise, where the sample is determined not to contain ctDNA or to have a ctDNA level below the limit of detection of the method of the present invention (e.g. where the mutant allele fraction is found to be not greater, or not statistically significantly greater, than the background sequencing error rate), the patient may be determined not to have, or to be in remission from, a cancer. The patient may therefore benefit from the avoidance of unnecessary anti-cancer treatment which could be associated with unwanted side effects.
In a third aspect, the present invention provides a method of treatment of a patient who has or has had a cancer, the method comprising:
In some cases, the anti-cancer treatment may be selected from chemotherapy, immunotherapy, radiotherapy and surgery. In particular, the anti-cancer treatment may comprise one or more of: vemurafenib, ipilimumab, pazopanib, dabrafenib and trametinib. In particular, where the patient has or has had melanoma, and the sample is determined to contain cfDNA (e.g., ctDNA), the aforementioned anti-cancer treatments may be suitable.
Without wishing to be bound by any particular theory, the present inventors believe that the methods of the present invention may find application beyond the field of cancer monitoring and cf DNA e.g., ctDNA detection. In particular, the INVAR algorithm may find use in forensic science (e.g. detecting trace amounts of a suspected perpetrator's (or victim's) DNA in a sample containing a larger fraction of another person's DNA, such as a suspected victim (or perpetrator, as the context directs), agriculture and food (e.g. to detect contamination), lineage tracing, clinical genetics, and transplant medicine. The ability of the INVAR method to improve signal-to-noise ratio by aggregating across many, e.g., all, mutant reads and optionally splitting (further analysing) by mutation class, makes this method attractive in applications where a sample is suspected of containing a minor fraction of a target DNA or other polynucleotide (e.g. RNA), including fragments thereof, that may differ in sequence at a number of loci from the DNA or other polynucleotide (e.g. RNA), including fragments thereof, making up the larger fraction of the sample.
Accordingly, in a fourth aspect the present invention provides a method for detecting a target polynucleotide in a sample in which the target polynucleotide is a minor fraction of the total polynucleotide in the sample, wherein the target polynucleotide and the non-target polynucleotide differ in sequence at a plurality of loci, the method comprising:
In some cases, the background sequencing error rate is or has been determined for each class of base substitution represented in said at least 2, 3, 4, 5, 6, 7, 8, 9, 10, 50, 100, 500, 1000, 2500 or 5000 loci, wherein the target polynucleotide fraction calculation in step (d) is performed for each base substitution class,
The target polynucleotide may be DNA or RNA.
In accordance with any aspect of the present invention, the patient is a mammal, preferably a human. The patient may have been diagnosed as having a cancer. In some cases, the patient may have had a course of treatment for the cancer and/or had surgery to excise a cancer.
In accordance with any aspect of the present invention, the method may comprise analysing a given sample in a plurality of (e.g. 2, 3, 4, 5, 6, or more) replicates, and using the signal in the replicates to improve confidence in the determination of presence or absence of cfDNA in the sample. In such cases, it is possible to relax other constraints of the methods of the present invention. For example, by using sample replicates it may be possible to omit the reads collapsing step. Nevertheless, the use of sample replicates and reads collapsing are not mutually exclusive, and so both sample replicates and reads collapsing may, in certain embodiments, both be employed in methods of the present invention.
In accordance with any aspect of the present invention, in some embodiments the analysis of the sample includes a size-selection step which separates out different fragment sizes of DNA.
In some embodiments the sample obtained from the patient is a limited volume sample comprising less than one tumour-derived haploid genome. In some embodiments the sequencing data from the sample represents sequencing coverage or depth of less than 1, 2, 3, 4, 10 or 12 haploid genomes.
In some embodiments the sample obtained from the patient is a limited volume sample selected from the group consisting of:
In some embodiments the sample obtained from the patient is:
In some embodiments the patient is healthy or has a disease (e.g. a cancer) and/or wherein the patient is human or a non-human animal (e.g. a rodent).
In some embodiments the animal has xenografted or xenotransplanted human tumour tissue.
In some embodiments the sample being analysed (e.g. the sample obtained from the patient) is subjected to a size-selection step in which genomic DNA (gDNA) fragments of >200 bp, >300 bp, >500 bp, >700 bp, >1000 bp, >1200 bp, >1500 bp or >2000 bp are filtered-out, depleted or removed from the sample prior to analysis, e.g. prior to DNA sequencing, to generate a size-selected sample.
In some embodiments the size selection step is carried out prior to sequencing library preparation or after sequencing library preparation.
In some embodiments the size selection step is a right-sided size selection employing bead-based capture of gDNA fragments.
In some embodiments the method comprises a step of suppressing outlier noise at a locus that is not consistent with the remaining distribution of patient-specific mutant signal in that sample (“outlier suppression”).
In some embodiments the likelihood of ctDNA presence in the sample is determined by the generalized likelihood ratio:
wherein the terms of the generalized likelihood ratio are as defined in the Supplementary Methods for Example 14.
In a fifth aspect the present invention provides a method for detecting variant cell-free DNA (cfDNA) in a sample obtained from a patient, where analysis of the sample includes a size-selection step which separates out different fragment sizes of DNA.
In some embodiments the sample obtained is a limited volume sample selected from the group consisting of:
In some embodiments said size-selection step comprises filtering-out, depleting or removing genomic DNA (gDNA) fragments of >200 bp, >300 bp, >500 bp, >700 bp, >1000 bp, >1200 bp, >1500 bp, or >2000 bp prior to analysis, e.g. prior to DNA sequencing.
In some embodiments the method comprises:
In some embodiments the sample obtained from the patient is:
In some embodiments the patient is healthy or has a disease (e.g. a cancer) and/or wherein the patient is human or a non-human animal (e.g. a rodent).
In some embodiments the animal model has xenografted or xenotransplanted human tumour tissue.
In some embodiments said analysis comprises next-generation sequencing (NGS) of the size-selected sample or a library generated from the size-selected sample.
In some embodiments said analysis comprises sequencing the size-selected sample or a library generated from the size-selected sample to generate sequence reads and further comprises analysis of the sequence reads selected from:
In some embodiments the t-MAD score is determined by trimming regions of genome that exhibit high copy number variability in whole genome datasets derived from healthy subjects and then calculating the median absolute deviation from log2R=0 of the non-trimmed regions of the genome.
In some embodiments the size selection step is carried out prior to, or after, a sequencing library preparation step.
In some embodiments the size selection step is a right-sided size selection employing bead-based capture of gDNA fragments.
In some embodiments the variant cell-free DNA is circulating tumour DNA (ctDNA).
In some embodiments the method is for early detection of cancer, monitoring of cancer treatment, detection of residual disease, is used to guide treatment decisions, to assess the status of a cancer in the patient or cancer progression or cancer response to treatment, or the need for or the type of further treatment for the patient.
In some embodiments the patient is a human or an animal model of a cancer (e.g. a rodent).
In some embodiments the variant cell-free DNA comprises:
In some embodiments the method is used to provide information to guide medical treatment, changes in diet, or physical training, or is used for forensic analysis or to identify individuals whose biological material is present in the sample or to identify organisms whose biological material is present in the sample.
In some embodiments the patient is a human child having or suspected of having a paediatric cancer. Paediatric cancers are often associated with difficulties in sample collection, e.g. due to the age of the patient, and samples may be of small volume and/or contain low levels of ctDNA. Paediatric cancers include: various brain tumours, lymphomas, leukaemia, neuroblastoma, Wilms tumour, Non-Hodgkin lymphoma, Childhood rhabdomyosarcoma, Retinoblastoma, Osteosarcoma, Ewing sarcoma, Germ cell tumors, Pleuropulmonary blastoma, Hepatoblastoma and hepatocellular carcinoma.
In a sixth aspect, the present invention provides a system for detecting target cell-free DNA (cfDNA) in a DNA-containing sample obtained from a patient, comprising:
In some embodiments the system is for use in a method of the present invention.
In a seventh aspect, the present invention provides a non-transitory computer readable medium for detecting target cell-free DNA (cfDNA) in a DNA-containing sample obtained from a patient, comprising instructions that, when executed by at least one processor, cause the at least one processor to perform operations comprising:
In some embodiments the medium is for use in a method of the present invention.
Embodiments of the present invention will now be described by way of examples and not limited thereby, with reference to the accompanying figures. However, various further aspects and embodiments of the present invention will be apparent to those skilled in the art in view of the present disclosure.
The present invention includes the combination of the aspects and preferred features described except where such a combination is clearly impermissible or is stated to be expressly avoided. These and further aspects and embodiments of the invention are described in further detail below and with reference to the accompanying examples and figures.
In describing the present invention, the following terms will be employed, and are intended to be defined as indicated below.
“and/or” where used herein is to be taken as specific disclosure of each of the two specified features or components with or without the other. For example “A and/or B” is to be taken as specific disclosure of each of (i) A, (ii) B and (iii) A and B, just as if each is set out individually herein.
“Computer-implemented method” where used herein is to be taken as meaning a method whose implementation involves the use of a computer, computer network or other programmable apparatus, wherein one or more features of the method are realised wholly or partly by means of a computer program.
“Patient” as used herein in accordance with any aspect of the present invention is intended to be equivalent to “subject” and specifically includes both healthy individuals and individuals having a disease or disorder (e.g. a proliferative disorder such as a cancer). The patient may be a human, a companion animal (e.g. a dog or cat), a laboratory animal (e.g. a mouse, rat, rabbit, pig or non-human primate), an animal having a xenografted or xenotransplanted tumour or tumour tissue (e.g. from a human tumour), a domestic or farm animal (e.g. a pig, cow, horse or sheep). Preferably, the patient is a human patient. In some cases, the patient is a human patient who has been diagnosed with, is suspected of having or has been classified as at risk of developing, a cancer.
A “sample” as used herein may be a biological sample, such as a cell-free DNA sample, a cell (including a circulating tumour cell) or tissue sample (e.g. a biopsy), a biological fluid, an extract (e.g. a protein or DNA extract obtained from the subject). In particular, the sample may be a tumour sample, a biological fluid sample containing DNA, a blood sample (including plasma or serum sample), a urine sample, a cervical smear, a cerebrospinal fluid (CSF) sample, or a non-tumour tissue sample. It has been found that urine and cervical smears contains cells, and so may provide a suitable sample for use in accordance with the present invention. Other sample types suitable for use in accordance with the present invention include fine needle aspirates, lymph nodes, surgical margins, bone marrow or other tissue from a tumour microenvironment, where traces of tumour DNA may be found or expected to be found. The sample may be one which has been freshly obtained from the subject (e.g. a blood draw) or may be one which has been processed and/or stored prior to making a determination (e.g. frozen, fixed or subjected to one or more purification, enrichment or extractions steps, including centrifugation). The sample may be derived from one or more of the above biological samples via a process of enrichment or amplification. For example, the sample may comprise a DNA library generated from the biological sample and may optionally be a barcoded or otherwise tagged DNA library. A plurality of samples may be taken from a single patient, e.g. serially during a course of treatment. Moreover, a plurality of samples may be taken from a plurality of patients. Sample preparation may be as described in the Materials and Methods section herein. Moreover, the methods of the present invention have been demonstrated to detect tumour-derived mutant DNA in urine samples (data not shown). Accordingly, use of blood or urine samples as a source of patient DNA potentially containing mutant tumour DNA to be detected is specifically contemplated herein. For forensic applications, the sample may be any fluid or tissue or item having or suspected of having a mixed DNA or RNA (e.g. target and background, such as perpetrator DNA or RNA and victim DNA or RNA). For the analysis of contamination, the sample may be any fluid, organism, item, foodstuff or plant having or suspect of having mixed DNA or RNA (e.g. target and background, such as contamination source (e.g. pathogen) DNA or RNA and non-contamination source DNA or RNA).
“Right-sided size selection” as used herein in some embodiments employ AMPure beads as described at https://research.fhcrc.org/content/dam/stripe/hahn/methods/mol_biol/SPRIselect %20User %20Guide.pdf (the entire contents of which is incorporate herein by reference). In particular, a 1× selection step as used in some embodiments implies a cut-off in-between the curves for 1.2× and 0.95×, so is estimated at around 200-300 bp.
“Blood spot” as used herein may in some embodiments be a dried blood spot sample. Typically, blood samples are blotted and dried on filter paper. Dried blood spot specimens may be collected by applying one or a few drops of blood (e.g. around 50 μl), drawn by lancet from the finger, heel or toe, onto specially manufactured absorbent filter paper. The blood may be allowed to thoroughly saturate the paper and may typically be air dried for several hours. Specimens may be stored in low gas-permeability plastic bag with desiccant added to reduce humidity, and may be kept at ambient temperature.
In accordance with some embodiments of the present invention, loci that carry mutations that are specific to the tumour of the patient may be identified. In some cases tumour DNA is sequenced so as to give an average 8 Gb of unique mapped reads per sample with an average of 80% of base pairs covered by >20 reads. In some cases, single nucleotide variants (SNVs) (versus a germline sequence, e.g., from a buffy coat sample) may be selected from the sequence data obtained from the tumour sample. In some cases, patient-specific loci are those which display SNVs with 1 mutant read and 10 total reads as determined from tumour sequencing. In some cases loci may be excluded if they show 1 forward (F) and 1 reverse (R) non-reference read (following read de-duplication) in the germline sequence (e.g. a buffy coat sample). Optionally, loci may be excluded if they are SNPs identified in common SNP databases, such as the 1000 Genomes database.
The sequence reads data may be provided or obtained directly, e.g., by sequencing the cfDNA sample or library or by obtaining or being provided with sequencing data that has already been generated, for example by retrieving sequence read data from a non-volatile or volatile computer memory, data store or network location. Where the sequence reads are obtained by sequencing a sample, the median mass of input DNA may in some cases be in the range 1-100 ng, e.g., 2-50 ng or 3-10 ng. The DNA may be amplified to obtain a library having, e.g. 100-1000 ng of DNA. The median sequencing depth of sequence reads (e.g. quality-filtered sequence reads) at each patient-specific loci may be in the range 500×-2000×, e.g., 750×-1500× or even 1200×-1400×. The sequence reads may be in a suitable data format, such as FASTQ.
The sequence read data, e.g., FASTQ files, may be subjected to one or more processing or clean-up steps prior to or as part of the step of reads collapsing into read families. For example, the sequence data files may be processed using one or more tools selected from as FastQC v0.11.5, a tool to remove adaptor sequences (e.g. cutadapt v1.9.1). The sequence reads (e.g. trimmed sequence reads) may be aligned to an appropriate reference genome, for example, the human genome hg19.
As used herein “read” or “sequencing read” may be taken to mean the sequence that has been read from one molecule and read once. Each molecule can be read any number of times, depending on the sequencing performed.
As used herein, “read family” may be taken to mean multiple sequencing reads arising from the same molecule (hence replicates). Because they are from the same starting molecule, each read will have the same start and end position in the human genome following alignment of the read. In addition, when molecular barcodes are ligated to starting molecules prior to PCR and sequencing, each read family will also have the same molecular barcode. The process of error suppression by molecular barcodes is described at the following URL: https://github.com/umich-brcf-bioinf/Connor/blob/master/doc/METHODS.rst (the contents of which as shown in 5 Mar. 2018 are expressly incorporated herein by reference).
As used herein “collapsing” or “reads collapsing” may be taken to mean that, given a read family (set of replicate reads), error-suppression for PCR and sequencing errors may be performed by generating a consensus sequence across the family for every base position. Thus, a family of N (number of) reads is ‘collapsed’ into a consensus sequence of one read, which consensus sequence can be expected to contain fewer errors.
Reads collapsing may be performed based on fragment start and end position and custom inline barcodes. A suitable tool is CONNOR described at https://github.com/umich-brcf-bioinf/Connor/blob/master/doc/METHODS.rst (the entire contents of which as shown in 5 Mar. 2018 are expressly incorporated herein by reference). CONNOR may be used with a consensus frequency threshold -f set to 8.8, 0.85, 0.9 or 0.95. CONNOR may be used with a minimum family size threshold -s set as 2, 3, 4, 5, 6, 7, 8, 9, or 10. Preferably, the consensus frequency threshold is 0.9 and the minimum family size threshold is 5.
Quality filters may be applied in the process of determining the number of mutant and wild-type reads/read families, as described in the Materials and Methods section herein.
In some cases, one or more MRD-filters are applied to focus on tumour-derived MRD read families. In particular, the MRD-filter step may comprise one or both of:
As used herein “barcode” or “molecular barcode” may be taken to mean a unique string of bases, generally but not necessarily of length <10 bp, e.g. a molecular barcode employed by the invention may be 6, 7, 8, 9, or 10 bp in length), that may be ligated to one or more DNA molecules as the first step during library preparation. As a result, read families (from above) may be uniquely identified, and thus linked to their starting molecule. This allows error-suppression via ‘reads collapsing’, as described above.
In some cases, a region either side of each patient-specific loci (e.g. 20, 15, 10 or 5 bp either side) may be employed to determine the error rate for each mutation class. In some cases, non-reference bases are only accepted if found to be present in both the forward F and reverse R read. In some cases if a locus displays mutant error-suppressed families in ≥3 separate libraries, it may be filtered out (“blacklisted”) on the basis of having a higher locus-specific error rate.
The sequencing error analysis may be carried out to determine the background error rates regardless of mutation class, and by separating data by mutation class. The error rate may be determined by taking the ratio between the sum of mutant reads in a class and the total number of reads in a class. In some cases, this ratio data may be resampled 100 times with replacement to obtain 95% confidence intervals of the error rate.
In accordance with some embodiments of the present invention a variant read for a particular patient-specific loci may be accepted only where the observed variant (e.g. SNV) matches the mutation determined in the tumour sequence at that loci. For example, if a C>T mutation was expected based on tumour sequencing/genotyping, but a C>A is observed in the mutant reads, then the mutant reads may ignored and may be excluded from the patient-specific signal. Alternatively or additionally, loci may be only considered as contributing to signal if there are at least ≥1 F and ≥1 R read family at that position. This has two advantages: to reduce single-stranded artefacts from sequencing, and to bias detection towards short fragments which have greater overlap between F and R reads in certain sequencing platforms, e.g. PE150 sequencing.
For each sample, the mutant allele fraction may be calculated across all patient-specific loci as follows:
In certain cases, the mutant allele fraction may be calculated by trinucleotide context. The mutant allele fraction by context may be based on tumour-weighted read families according to the formula:
wherein:
AFcontext is the allele frequency of a given (e.g. trinucleotide) context; tumourAF is the allele frequency of the locus as determined by sequencing DNA obtained directly from the tumour; and MRD-like loci are the mutation-containing loci determined from the tumour of the patient and which have been filtered to select for minimal residual disease signal.
The significance of the observed number of mutant reads may be determined using a one-sided Fisher's exact test, given a contingency table with the number of mutant reads and total reads for both the sample of interest, and from the background error rate.
In some embodiments of the present invention each sample may be split into a plurality of mutation classes (e.g. 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 or all 12 of the following SNV classes: C>G, G>C, T>G, A>C, C>A, G>T, T>C, A>G, T>A, A>T, C>T and T>C) based on the mutation class expected at that locus from tumour sequencing.
Variant reads may be integrated for each class, as above. Multiple one-sided Fisher's exact tests may be used to determine the significance of the number of mutant read families observed given the background error rate for that mutation class. This method will generate 12 P-values per sample, which may then be combined using the Empirical Brown's method. If a sample had no data in a class, that class may be treated as having zero mutant reads and thus a P-value of 1.
To improve specificity further, in some embodiments, the method of the present invention may require samples to have mutant reads in 2 separate classes; this ensures that detection is based on signal being present in multiple loci subject to different types of error processes.
The significance threshold for combined P-values obtained by INVAR may in some cases be determined using Receiver Operating Characteristic analysis on patient-specific (test) and non-patient-specific (control) samples. For example, the analysis may employ the OptimalCutpoints package in R, with the ‘MaxEfficiency’ method, which maximises classification accuracy.
In some cases, background error rates may be subtracted from the observed allele fraction. This can be performed with or without taking into account differences in error rate by class. If the observed mutant allele fraction is less than the background error rate, then the background-subtracted allele fraction may be set at zero. For background subtraction by mutation class for a sample, the error rate of each of the classes may be subtracted from the mutant allele fraction of that class. A mean allele fraction may then be calculated from each of the individual background-subtracted allele fractions, weighted by the total number of read families observed in that class.
The following is presented by way of example and is not to be construed as a limitation to the scope of the claims.
MelResist (REC number 11/NE/0312) is a translational study of response and resistance mechanisms to systemic therapies of melanoma, including BRAF targeted therapy and immunotherapy. For each patient in this cohort, fresh frozen metastatic tumour biopsies and plasma samples were collected before the initiation of treatment and plasma was collected at varying time points during treatment. Patients may have received multiple lines of treatment over time. Demographics and clinical outcomes are collected prospectively. The study was coordinated by the Cambridge Cancer Trials Unit-Cancer Theme.
Peripheral blood samples were collected longitudinally at each clinic visit in S-Monovette 9 mL EDTA tubes. For this study, up to 8 samples per patient were analysed from their serially-collected samples. One aliquot of whole blood at baseline was stored at −80° C. for germline DNA. For plasma collection, samples were centrifuged at 1600 g for 10 minutes within an hour of the blood draw, and then an additional centrifugation of 20,000 g for 10 minutes was carried out. Plasma aliquots were stored at −80° C.
Extraction of DNA from Fresh Frozen Tissue and Plasma
Up to 30 mg of each fresh frozen tissue biopsy sample was combined with 600 μL RLT buffer (QIAGEN), then placed in a Precellys CD14 tube (Bertin Technologies) and homogenised at 6,500 rpm for two bursts of 20 seconds separated by 5 seconds. DNA was then extracted using the AllPrep extraction kit (Qiagen) as per the manufacturer's protocol.
Genomic DNA was extracted from 10 mL whole blood using the Gentra Puregene Blood Kit (Qiagen) as per the manufacturer's protocol. Eluted DNA concentration was quantified using Qubit (ThermoFisher Scientific).
Plasma samples were extracted using the QIAsymphony instrument (Qiagen) using a 2 mL QIAamp protocol. For each QIAsymphony batch, 24 samples were extracted, which included a healthy individual control sample (Seralab). Plasma samples were eluted in 90 μl water and stored at −80° C.
CT imaging was acquired as part of the standard of care for each patient and were examined retrospectively. The slice thickness was 5 mm in all cases. All lesions with a largest diameter greater than ˜5 mm were outlined slice by slice on the CT images by an experienced operator, under the guidance of a radiologist, using custom software written in MATLAB (Mathworks, Natick, Mass.). The outlines were subsequently imported into the LIFEx software application25 in NifTI format for processing. Tumour volume was then reported by LIFEx as an output parameter from its texture based processing module.
To quantify the cfDNA concentration of each sample, digital PCR was carried out using a Biomark HD (Fluidigm) using Taq-man probes for the housekeeping gene RPP30 (Sigma Aldrich), and a unique XenT locus, labelled with ROX and FAM, respectively. 55 PCR cycles were used. The RPP30 assay was 65 bp in length. The estimated number of RPP30 DNA copies per μl of eluate was used to determine the cfDNA concentration in the original sample.
Tumour and buffy coat (germline) library preparation, sequencing and variant calling were performed as described by Varela et al.26, using the SureSelectXT Human All Exon 50 Mb (Agilent) bait set or a custom targeted sequencing bait set. Eight samples were multiplexed per pool and each pool loaded on to two lanes of a HiSeq 2000 (Illumina), giving an average 8 Gb of unique mapped reads per sample with an average of 80% of base pairs covered by >20 reads. Targeted sequencing was carried out using the Sanger CGP Cancer Genes V3 panel for 365 genes relevant to cancer, as described previously27. For this exploratory analysis, all mutation calls from tumour sequencing were included in the TAPAS panel design (see Results). Loci were excluded if they showed 1 forward (F) and 1 reverse (R) non-reference read (following read de-duplication) in the buffy coat sample.
TAPAS libraries from 10 patients were prepared in duplicate using the Rubicon ThruPLEX Tag-seq kit. The median input mass was 4.4 ng for plasma DNA libraries (IQR 3.2-10.0 ng). In order to compare error rates between molecularly barcoded libraries and non-molecularly barcoded libraries, additional plasma libraries were prepared using the Rubicon ThruPLEX Plasma-seq kit. cfDNA samples were vacuum concentrated at 30° C. using a SpeedVac (ThemoFisher) prior to library preparation where required.
Based on the starting concentration of DNA in each sample, the number of PCR amplification cycles during the ThruPLEX protocol was varied between 7-15 cycles as recommended by the manufacturer28. Following amplification and sample barcoding, libraries were purified using Ampure XT beads (Beckman Coulter) at a 1:1 ratio. Library concentration was determined using the Illumina/ROX low Library Quantification kit (Roche) with two sample dilutions, in triplicate. 1:10 diluted libraries were run on Bioanalyser HS chip (Agilent) to determine library fragment size.
333-750 ng of each library was captured using the Agilent SureSelectXT protocol, with the addition of i5 and i7 blocking oligos (IDT) as recommended by the manufacturer29. Libraries were pooled for capture between 1 to 3-plex pools up to a maximum capture input of 1000 ng. 13 cycles were used for post-capture amplification. Post-capture libraries were purified with Ampure XT beads at a 1:1.8 ratio, then were quantified and library fragment size was determined as before. A median of 9 TAPAS libraries were pooled per lane of HiSeq 4000.
FastQC v0.11.5 was run on all FASTQ files, then cutadapt v1.9.1 was used to remove known 5′ and 3′ adaptor sequences specified in a separate FASTA of adaptor sequences. Trimmed FASTQ files were aligned to the hg19 genome using BWA-mem v0.7.13 with a seed length of 19. Duplicates were marked using Picardtools v2.2.4 MarkDuplicates. BAM files were indexed using Samtools v1.3.1. Local realignment for known indels and base quality recalibration were carried out using GATK v3.7. Next, regions to be disregarded on the basis of having a high level of sequencing noise (also known as “blacklisted regions”) identified by the ENCODE consortia were removed from BAM files.
Error suppression was carried out on ThruPLEX Tag-seq library BAM files using Connor30, which generates a consensus sequence between replicate sequencing reads based on fragment start and end position, and custom inline molecular barcodes. The consensus frequency threshold -f was set as 0.9, and the minimum family size threshold -s was set as 5 following analysis of error rates vs. proportion of data retained; read families below these thresholds were discarded. ThruPLEX Plasma-seq libraries were also used as input for Connor with the same settings using a custom shell script. This script adds a false barcode and stem to the appropriate end of each read, and modifies the CIGAR string.
Samtools mpileup v1.3.1 was used at patient-specific loci to determine the number of mutant and wild-type reads/read-families for raw and error-suppressed data. The following settings were used: -d 10000 (maximum depth threshold) --ff UNMAP (exclude unmapped reads) -q 13 (minimum Phred mapping quality score) -Q 13 (minimum Phred base quality score) -x (ignore overlaps) -f ucsc.hgl9.fasta. VCF Parser31 v1.6--split was used to separate multi-allelic calls, and SnpSift extractFields was used to extract the columns of interest. For analysis of non-error-suppressed TAPAS data, a minimum of 5 reads were required at a locus; the threshold for error-suppressed data was a minimum of 1 read family (consisting of 5 members). Individual data points (i.e. a single locus in a single sample) were filtered if the mapping quality/strand bias (MQSB) at that locus, as determined by Samtools mpileup, was <0.01.
TAPAS was applied to patients' first plasma time point in order to call variants in either the genes of interest that were tiled across, or in the bait regions either side of the patient-specific variants, that may have been missed from tumour exome sequencing alone. Mutect2 (GATK) was used for initial mutation calling, and was given the hg19 COSMIC database VCF, the dbSNP database VCF, the baitset BED file (including resistance loci and genes of interest). The matched buffy coat exome BAM was used as the germline sample.
To learn the background error rates, off-target bases from TAPAS data were used. Sequencing data from patients was used for this purpose, since germline events could be removed based on exome sequencing of buffy coats, and known tumour loci could be excluded. Thus, 10 bp either side of each patient-specific variant was used to determine the error rate for each class of SNV. We specified that non-reference bases had to be present in both the F and R read. To avoid possible biological contamination of error-rates, loci were excluded if they had ≥1 overlapping mutations in COSMIC. In addition, following error suppression, each locus was assessed individually in all the samples belonging to the same patient, and if a locus showed mutant error-suppressed families in ≥3 separate libraries, it was disregarded from further analysis. Given a background error rate per read family of ˜6×10−5, the probability of observing mutant read families at a single locus in ≥3 samples (out of a median of 6 samples per patient) from the same individual by chance, with an average of 200 read families per locus, gives a binomial probability of ˜1×10−12.
This analysis was carried out to determine the background error rates regardless of mutation class, and by splitting data by mutation class. The error rate was determined by taking the ratio between the sum of mutant reads in a class and the total number of reads in a class. This data was resampled 100 times with replacement to obtain 95% confidence intervals of the error rate.
Detection of ctDNA was carried out for patient-specific loci only, i.e. if a C>T mutation was expected based on tumour genotyping, but a C>A was observed, then the mutant reads were ignored and did not contribute to the patient-specific signal. Furthermore, loci were only considered as contributing to signal if there was at least ≥1 F and ≥1 R read family at that position. This has two advantages: to reduce single-stranded artefacts from sequencing, and to bias detection towards short fragments which have greater overlap between F and R reads using PE150 sequencing.
For each sample, the mutant allele fraction was calculated across all patient-specific loci as follows:
The significance of the observed number of mutant reads was determined using a one-sided Fisher's exact test, given a contingency table with the number of mutant reads and total reads for both the sample of interest, and from the background error rate.
Since differences in error rates were observed between SNV classes, each sample was split into 12 based on the mutation class expected at that locus from tumour sequencing. Variant reads were integrated for each class, as above. Multiple one-sided Fisher's exact tests were used to determine the significance of the number of mutant read families observed given the background error rate for that mutation class. This generated 12 P-values per sample, which were then combined using the Empirical Brown's method, which is an extension of Fisher's method that can be used to combine dependent P-values16. If a sample had no data in a class, that class was treated as having zero mutant reads and thus a P-value of 1. To improve specificity of this approach further, we required samples to have mutant reads in 2 separate classes; this was in order to ensure that detection was based on signal being present in multiple loci subject to different types of error processes.
All patients were sequenced with the same sequencing panel, and since 99.9% of the variants were private to each patient (i.e. unique to that individual only), all other patients could be used to determine the false positive for ctDNA detection and thus set the P-value threshold for the cohort. This approach leverages the redundant sequencing performed, and utilises the multiple samples sequenced from each individual to exclude germline variants. As such, TAPAS data was split into patient-specific and non-patient-specific based on whether each locus was mutated in the patient's tumour. Non-patient-specific data was used for determining significance thresholds.
To use patients as controls, the technical noise should be separated from any true biological signal that may be detected in plasma but was missed in the tumour. Therefore, using error-suppressed non-patient-specific data, loci were disregarded from further analysis (“blacklisted”) if they contained mutant read families in 3 separate libraries from the same individual, which we calculated to be sufficiently unlikely to be observed to warrant ignoring these loci (P=1×10−12, see Determination of background error rates). As a result, 44 loci out of 12,558 (0.35%) were disregarded from further analysis (“blacklisted”). Although imperfect tumour and buffy coat genotyping of patients may result in residual biological signal in control samples, this was preferable to the cost of sequencing many control samples with the same panel and discarding non-patient-specific data.
The significance threshold for combined P-values obtained by INVAR was determined using Receiver Operating Characteristic analysis on patient-specific (test) and non-patient-specific (control) samples using the OptimalCutpoints package in R, with the ‘MaxEfficiency’ method, which maximises classification accuracy.
Plasma cfDNA was obtained from one healthy individual (Seralab), and mutant cfDNA was obtained from one patient at a high tumour burden time point (MR1004; 2,746 patient-specific mutations). The cfDNA concentrations of the eluates were equalised using water, then the patient's sample was serially diluted by healthy cfDNA in a 1:5 ratio to give a final 15,625× dilution of the original cfDNA eluate.
Library preparation was carried out in duplicate using the ThruPLEX Plasma-seq kit with 3.7 ng input for all libraries. 50 ng spike-in dilution experiment Equal masses of plasma cfDNA from 6 patients were pooled to produce a hypothetical patient with a total of 9,636 patient-specific variants. A pool of plasma cfDNA was generated from 11 healthy individuals (Seralab). The cfDNA concentrations of the patient sample and healthy pool were equalised using water, then the patient sample was serially diluted by healthy cfDNA in a 1:10 ratio to give a 100,000× dilution of the original 1× pooled sample. Library preparation was carried out with the ThruPLEX Tag-seq kit, in duplicate, with input amounts of up to 50 ng per library. For libraries with an expected allele fraction greater than the limit of detection of TAPAS without error suppression, we reduced the input material into library preparation to conserve patient plasma DNA where it was certain to be detected.
To test the limit of detection of INVAR-TAPAS with varying numbers of mutations, both the patient-specific mixture experiment and all non-patient-specific data were downsampled to between 50-5,000 mutations. BRAF was always included in each set of sampled mutations to simulate a panel design for BRAFout patients. Mutations were sampled iteratively 100 times, and detection of ctDNA was tested using INVAR.
Background-Subtraction for ctDNA Quantification
To accurately determine mutant allele fractions down to parts per million, background error rates were subtracted from the observed allele fraction. This can be performed with or without taking into account differences in error rate by class. If the observed mutant allele fraction was less than the background error rate, then the background-subtracted allele fraction was set at zero.
For background subtraction by mutation class for a sample, the error rate of each of the 12 classes was subtracted from the mutant allele fraction of that class. A mean allele fraction was then calculated from each of the individual background-subtracted allele fractions, weighted by the total number of read families observed in that class.
Variants removed by blacklisting (i.e. filtered out on the basis of having a higher locus-specific error rate as described above) were previously excluded on the basis of showing evidence for being biological signal. We attempted to call mutations from this blacklist for variants that were known mutations. Therefore, the data were intersected with the COSMIC database for known driver mutations (number of overlapping mutations 5). For each mutation locus, the background error rate for that locus was determined using non-patient-specific data (i.e. patients whose tumours had been genotyped as negative for that mutation). A one-tailed Fisher's exact test was used to test the significance of the number of mutant reads in the sample given the total depth in that sample, and the mutant reads and total depth in the background. The P-value threshold was set as 0.05 and corrected for multiple hypotheses by the Bonferroni method. Individual mutation calls were confirmed by aggregating mutant reads across multiple, temporally-separated samples.
To achieve high sequencing depth at defined loci that were mutated in the patients' tumours, a tailored hybrid-capture sequencing panel was designed based on single nucleotide variants (SNVs) identified in sequencing of tumour biopsies. SNVs with 1 mutant read and 10 total reads were selected from exome sequencing (9 patients) or targeted sequencing (1 patient) of a baseline metastatic biopsy. The median number of SNVs identified per patient was 673 (IQR 250-1,209;
The finalised bait set was applied to libraries generated in duplicate from serially collected plasma cfDNA samples, collected over two years (maximum 8 samples per patient). DNA was extracted from 2 mL of plasma, and the median input mass was 4.4 ng for plasma DNA libraries (IQR 3.2-10.0 ng). A median of 9 TAPAS libraries (IQR 8-12) were pooled per lane of HiSeq 4000 (PE150). For each of the patient-specific loci, the median depth of quality-filtered reads (Methods) was 1,367× for each sample (IQR 761-1,886×).
To identify additional mutations covered by the panel that may have been missed by tumour sequencing, an additional mutation calling step was performed on the first plasma time point, pre- or on starting drug treatment, when ctDNA levels were expected to be higher. Plasma mutation calling added a median of 19 SNVs mutations per patient (IQR 9-22; not shown) for subsequent analysis, giving a total of 12,558 patient-specific SNVs across the cohort. The observed rate of de novo identification of SNVs in our cohort was in line with the estimate of 14.4 coding mutations per Mb in melanoma (IQR=8.0-24.9) reported previously10. The BRAF V600E mutation was found in 9 out of 10 patients, and a further 18 mutations were shared between any two patients. Overall, 99.9% of the targeted mutated loci were unique to an individual patient.
We sought to learn the background error rates (i.e. the rate of observing a mutant base that was not expected) with and without error suppression in TAPAS sequencing data. Bases either side of patient-specific variants were studied as they would have a comparable sequencing depth to patient-specific variants, and would be subject to the same technical biases. In order to leverage this off-target sequencing of patient samples, germline events and potential biological signals were excluded if they occurred multiple times in samples from the same individual (Methods); these loci were set aside for subsequent de novo mutation calling.
Error suppression can be achieved by determining the consensus sequence across a read family using read collapsing. To achieve this, duplicate reads were grouped into ‘read families’ based on both start and end fragment positions, previously termed ‘endogenous barcodes’11,12, and molecular barcodes. Read families were collapsed, and a minimum requirement was set at 90% consensus between all family members for a base to be called. Without error suppression, the average background error rate was 2×10−4. Prior to applying error suppression, we determined the optimal minimum number of duplicates per read family (‘family size’). The proportion of read families retained and corresponding error rates for data with minimum family size requirements of 1, 2, 3 and 5 are shown in
Using a stringent level of error suppression (consensus required in 90% of family members, with a minimum family size of 5), with a median of 4.4 ng input, we obtained a median of 3.2×105 read families at each time point (IQR 8.7×104-6.2×105), each of which cover a locus that is mutated in that patient's cancer. Under the assumption that each such read family corresponds to a single molecule, we were thus able to probe for each sample hundreds of thousands of target molecules even though the starting material contained only ˜1300 copies of the genome.
When ctDNA levels are low, many patient-specific loci will have no mutant DNA fragments at that position (
The significance of the observed number of mutant reads was determined using a one-sided Fisher's exact test, given a contingency table with the number of mutant reads and total reads for both the sample of interest, and from the background error rate. Mutant reads were only considered as contributing to signal at a locus if there was at least one forward (F) and one reverse (R) read from PE150 sequencing data; this may suppresses sequencing artefacts, and also biases towards data from short cfDNA fragments (covered by reads in both direction) which are enriched in ctDNA13-15.
Based on known differences in error rate between base substitutions in hybrid-capture sequencing7, we assessed error rates in TAPAS data by mutation class using INVAR. Data were split into 12 classes (C>G, G>C, T>G, A>C, C>A, G>T, T>C, A>G, T>A, A>T, C>T, T>C), which showed differences in error rate by class both before and after error suppression (
We developed an algorithm to detect ctDNA based on splitting the read families from each sample into the 12 classes; a P-value was derived for each error class separately using Fisher's Exact Test, and P-values were combined using the Empirical Brown's method, an extension of Fisher's method that can be used to combine dependent P-values16 (Methods). To increase specificity of this approach further, we specified that mutant signal must be present in at least two mutation classes, thereby reducing the reliance of detection on an individual noisy locus or class.
All patients were sequenced with the same sequencing panel, and since 99.9% of the variants were unique to each patient, all other patients could be used to determine the false positive for ctDNA detection and thus set the P-value threshold for the cohort for each detection algorithm. This approach leverages the redundant sequencing performed that would otherwise be discarded, and utilises the multiple samples sequenced from each individual to exclude germline variants (Methods). As such, sequencing data was split into ‘patient-specific’ and ‘non-patient-specific’ based on whether the locus was mutated in the patient's tumour. Significance thresholds for detection were determined empirically, using Receiver Operating Characteristic (ROC) analysis on patient-specific (test) and non-patient-specific (controls) samples using the OptimalCutpoints package in R, which maximises classification accuracy. In accordance with the present invention, either ROC analysis can be used to identify an optimal threshold based on maximising both sensitivity and specificity, or specificity can be fixed at a certain level, e.g. 99.5%, and the sensitivity interrogated.
To assess the sensitivity of INVAR-TAPAS, a spike-in dilution experiment was generated in duplicate with 3.7 ng per library using plasma DNA from a patient for whom 2,743 mutations were covered in the TAPAS panel. Using error suppression with endogenous barcodes, we first applied INVAR without splitting reads into mutation classes, and detected a sample with an expected mutant allele fraction of 1.9×10−6 (
Next, a second spike-in dilution experiment was made in duplicate with up to 50 ng input cfDNA, and molecular barcodes were used. For this experiment, DNA was pooled from 6 patients, and serially diluted in healthy individual DNA (Methods). The patients' cfDNA pool comprised a total of 9,636 patient-specific mutations. 50 ng input DNA corresponds to the cfDNA in 3.0 mL plasma from this cohort (median cfDNA concentration 5,160 copies/mL). Using INVAR without analysis by class, we detected the expected 3 ppm mutant allele fraction spike-in sample, at an observed allele fraction of 9 ppm (
INVAR was then applied by splitting samples into 12 mutation classes, as described above. By leveraging differences in error rate between mutation classes, significant detection was achieved down to 0.3 ppm (
To test the sensitivity of this approach with a smaller panel design, a subset of between 50 and 5,000 mutations were randomly sampled in silico alongside the BRAF V600 locus, and detection of ctDNA using INVAR by mutation class was repeated iteratively (Methods). BRAF V600 was included in each sampled panel to simulate a panel design for BRAFmut patients. The sensitivity achieved for each number of mutations is shown in
The start and end positions of the reads were used to determine fragment size distributions. Error-suppressed data from all plasma samples were combined, and the distribution of fragments was calculated (
Given these findings, we aimed to enrich for mutant signal through in silico size selection. Based on the size ranges that showed enrichment of ctDNA, data were size-selected in silico for reads within the size ranges of 115-190 bp, 250-400 bp and 440-460 bp. These relatively wide ranges were chosen in order to minimise loss of rare mutant alleles, as the size distributions of mutant and wild-type fragments were mostly overlapping. Overly stringent size selection may result in drop out of rare mutant molecules, which becomes increasingly problematic as ctDNA levels approach parts per million. In principle, with more input DNA and further sequencing, narrower filters may be applied to generate a stronger enrichment factor. When applied to plasma samples and spike-in dilutions, size selection produced a median enrichment of 6.3% in ctDNA relative to wild-type while retaining 93.7% of mutant reads. The extent of enrichment post-size-selection was related to the starting mutant allele fraction of the sample, and followed an exponential relationship with decreasing mutant allele fraction (
Across the cohort, ctDNA mutant allele fractions were compared against volumetric CT imaging data, which showed a Pearson correlation of 0.67 (P=0.0002;
Across all time points, there was a Pearson correlation of 0.86 between ctDNA and serum lactate dehydrogenase, a prognostic marker used for melanoma patients (P=2.2×10−1s;
Following initiation of systemic therapy, 3 out of 10 patients' ctDNA declined to levels below 10 ppm. We found that patients whose ctDNA fell to less than 10 ppm had 24 months longer overall survival compared to patients with higher levels of residual ctDNA (median 954 vs. 229 days; log-rank test P=0.009;
Despite the limited DNA input mass used for library preparation (median 4.4 ng per library, 1320 haploid genomes), 40% of plasma samples had significantly-detected allele fractions that were below the theoretical limit of detection (with 95% sensitivity) using a perfect single-locus assay
When ctDNA levels are sufficiently high, resistance mutations may be identified de novo, and clonal evolution may be monitored through changes in the allele fractions of mutations9. An example from one patient (MR1022) is shown in
A combination of multiplexed deep sequencing of thousands of tumour-derived mutations and integration of variant reads enabled us to detect ctDNA down to 0.3 ppm. Through the characterisation of error rates and fragmentation patterns from cfDNA sequencing data, we optimised our workflow for hybrid-capture sequencing of cfDNA. In this study, we analysed a large number of mutations for each patient by using all mutations identified by exome or targeted sequencing of baseline tumour biopsies. This enabled sensitive analysis despite limited DNA input amounts into library preparation, approximately 10-fold lower than were used for other high-sensitivity amplicon and hybrid-capture methods6,8. Error suppression was used to reduce background sequencing errors, and in silico size selection was used to enhance for mutation signal. By generating a large amount of patient-specific reads overlying known tumour mutations for each patient, TAPAS compensates for small input amounts, and for data losses incurred by error suppression and size selection, while still retaining sufficient mutant reads for highly sensitive detection. As a result, both high sensitivity (below parts per million) and high specificity (>99.5%) were achieved.
INVAR-TAPAS leverages differences in error rates between mutation classes to detect rare mutant alleles while efficiently using the available data. Detection by mutation class, followed by the combination of each test statistic, allowed each class to contribute to the overall signal based on its background error rate. We used a method for combining P-values for correlated datasets16 to account for dependence between mutation classes. Here we used analysis by 12 mutation classes; a larger dataset might enable analysis based on a greater number of sequence subsets such as by tri-nucleotide context or by individual locus, which may improve resolution into error rates further still.
Using error-suppressed data, size profiles were visualised for both mutant and wild-type reads with minimal confounding error from PCR and/or sequencing. We confirm that ctDNA is enriched in short plasma cfDNA fragments, and provide evidence for enrichment of mutant DNA in di-nucleosomal DNA, which may have contributed to previous findings of longer mutant DNA in the plasma of cancer patients20,21. We applied size selection to our data, which was first demonstrated in the field of non-invasive pre-natal testing22 (where fetal DNA fragments are shorter than maternal fragments23), and is beginning to be used experimentally on cancer patient samples″. Fan et al.22 highlight the challenge of retaining mutant molecules with size selection, which we confirm is largely due to highly overlapping size distributions of mutant and wild-type fragments. In the current study, we opted for a relaxed size selection to retain a large fraction of the starting mutant molecules, and demonstrated that a relaxed cut-off can provide benefit, particularly when the mutant fraction is very low (in the range of 1 ppm mutant allele fraction and lower). With greater sequencing depth and DNA input, more stringent filtering can provide further enrichment.
INVAR-TAPAS leverages knowledge of tumour-derived mutations, which requires analysis of an initial sample with high tumour content. This method has potential utility for monitoring disease recurrence post-treatment, particularly after surgery where the tumour tissue DNA can be obtained for sequencing. We showed in one example where this method detected as little as 1.3 cm3 residual disease, with 9.1 ppm ctDNA; this mutant allele fraction observed is consistent with predicted allele fractions for given tumour volumes from a previously described model6, and indicates that INVAR-TAPAS may theoretically identify lesions at the limit of detection of CT detection. Earlier detection of relapse or disease progression with a high-sensitivity approach may facilitate earlier initiation of adjuvant therapy or change of therapy. For guiding subsequent therapy, we demonstrate that mutations can be identified de novo, though the sensitivity of this is directly proportional to the number of molecules analysed at that locus which can be limiting. Signal may be further integrated across multiple longitudinal samples to enhance identification in the context of limited input DNA. One advantage of the present approach is that low level signal in a previous sample can provide evidence to support mutation detection in a later sample. Thus, each longitudinal sample supports another.
This tailored approach may be carried out using different types of input data from the plasma, and different lists of mutations to inform analysis. Tumour-derived mutations can be identified using exome sequencing as demonstrated here, but also across smaller focused panels or larger scales such as whole genome. In this cohort of 10 melanoma patients, exome sequencing was sufficient to identify hundreds to thousands of mutations per patient. Based on known mutation rates of cancer types24, exome sequencing would also suffice for many cancer types with relatively high mutation rates, for example: lung, bladder, oesophageal, or colorectal cancers. For cancers with a mutation rate of ˜1 per megabase or less24, whole-genome sequencing of tumours for mutation profiling would be desirable: for ovarian and brain cancers, this would result in thousands of mutations identified per patient.
To generate data for INVAR, we used targeted sequencing with patient-specific panels (e.g. TAPAS), which provided deep sequencing of a large number of mutations, but required development of patient-specific sequencing panels. This is cost-effective for generating data for INVAR from longitudinal samples, as they can be analysed with the same TAPAS panel. In a different implementation, whole-exome or whole-genome sequencing without the design of patient-specific panels may produce similar data suitable for INVAR. While reducing the complexity of the workflow, with this approach much of the sequencing data would not cover tumour-mutated loci (and thus would not be informative for INVAR), resulting in fewer patient-specific read families available for INVAR unless more sequencing is performed. As sequencing costs decrease, and tumour sequencing becomes more frequent with the advent of personalised oncology, we suggest that integration of variant reads from error-suppressed sequencing of plasma cfDNA would provide a highly sensitive means of treatment monitoring, disease surveillance and detection of residual disease.
To achieve high sequencing depth at defined loci that were mutated in the patients' tumours, tailored sequencing panels were designed based on single nucleotide variants (SNVs) identified in sequencing of fresh frozen or FFPE tumour biopsies from 48 patients with Stage II-IV melanoma. Mutation calling was performed on all tumour biopsies, and variant calls were filtered to exclude common SNP sites, repeat regions, and loci with signal in the patient's matched germline DNA (Methods).
Mutation profiles were assessed in fresh frozen tumour biopsy sequencing (
We found that background error rates in plasma from hybrid-capture sequencing varies between trinucleotide context using error-suppressed data with a minimum family size threshold of 2 (
The use of trinucleotide contexts enables the determination of background error rates down to 1 in 10 million through aggregation of read families across contexts; to achieve the same accuracy of background-error estimation on a per-locus level would require a vast number of samples to be sequenced. The use of trinucleotide contexts enables maximum retention of read families following error-suppression (
When ctDNA levels are low, many patient-specific loci will have no mutant DNA fragments at that position. Therefore, to overcome sampling error, all patient-specific read-families were aggregated and analysed together using INVAR. For each sample, mutant read and total read families were aggregated by trinucleotide context, and the proportion determined:
The significance of the observed number of mutant reads was determined using a one-sided Fisher's exact test for each context, generating a vector of P-values for each sample. The length of each P-value vector varied between samples, since the number of contexts represented in each patient differed based on the mutation profile of that patient. In order to account for this, and to take into consideration that in the minimal residual disease (MRD) setting only a small number of molecules would be present, we combined the P-values from the 6 most significant trinucleotide contexts per sample. This was performed on both test samples and control samples, and controls were used to determine a P-value cut-off with 97.5% specificity.
The high mutation rate in cutaneous melanoma is almost entirely attributable to the abundance of cytidine to thymidine (C>T) transitions, which are characteristic of UV-induced mutations (Hodis et al., 2012). We confirm this mutational signature in our data (
In melanoma, CC>TT mutations provide an opportunity to achieve ultra-low error rates, since any stochastic PCR/sequencing error would have to occur twice in series. CC>TT mutations can be aggregated as their own mutation class, whereas individual indels each would have a separate error profile. Thus, CC>TT mutations may be sufficiently prevalent in the data to allow the interrogation of a sufficient number of molecules to take advantage of the low noise profile. We are currently generating a script to identify mutant reads with CC>TT in adjacent bases from data error-suppressed with a minimum family size of 2. These mutations can be treated as a separate class with their own error profile for INVAR.
In order to optimise INVAR for detection of residual disease, we generated a spike-in dilution series using a mixture of patient cfDNA and healthy individual cfDNA, and characterised the signal occurring at the lowest dilutions. For this experiment, cfDNA was pooled from 6 patients to create a theoretical patient with a total of 9,636 patient-specific mutations. This pool was then serially diluted in healthy individual DNA (Methods).
Histograms of mutant allele fractions of individual patient-specific mutations for the dilution experiment are shown in
At the lowest levels of residual disease, ctDNA would be found in individual mutant molecules at individual loci. It would be highly unlikely for many mutant molecules to be concentrated entirely at one locus without other loci being represented, and this is supported by our data (
Using this approach, the probability of mis-genotyping a SNP (expected AF=50%) by focusing on loci with 2 or fewer molecules in 50 total molecules is 1×10−12 (2 or fewer successes out of 50; p=0.5). This is reduced further by the prior SNP filtering carried out at the tumour sequencing stage, whereby loci are filtered based on common SNP databases (i.e. 1000 Genomes ALL, EUR).
In addition, we also placed a lower limit on the number of mutant reads per locus. Mutant reads were only considered as contributing to signal at a locus if there was at least one F and one R read at a locus. Given that we sequenced with PE150, requiring overlapping F and R mutant read support served a dual purpose of suppressing sequencing artifacts, and selected for mutant reads from short cfDNA fragments (supported by reads in both directions), which are slightly enriched for ctDNA (
Together, these above parameters focus the INVAR algorithm on aggregating signal from mutant molecules most likely to originate from the tumour that were randomly sampled in the context of MRD.
We assessed the representation of mutations in plasma at time points where ctDNA was high. We found a correlation between tumour exome AF and plasma AF (
Patient-specific sequencing provides an opportunity to leverage such tumour prior information. As such, INVAR signal per locus was weighted by the tumour AF prior to aggregation of signal by mutation context. This was performed by dividing both the number of mutant read families and total read families at that locus by 1-tumour allele fraction. This places greater weight on loci more likely to contain true signal in plasma.
The raw number of mutant families per locus is shown in
Next, we applied INVAR to exome sequencing data to demonstrate its generalisability to non-individualised sequencing data. Plasma exome sequencing was carried out on a subset of samples from patients with Stage IV disease.
For exome sequencing data, we did not use molecular barcodes in order to demonstrate that INVAR may be applied to existing exome data, where the use of molecular barcodes is less frequent. Given that INVAR aims to target many loci, the families of interest are spread across multiple genomic regions, and thus the likelihood of clashing endogenous barcodes is reduced. This probability is further reduced by the reduced number of families per locus obtained by exome sequencing. We pooled 3-6 exome libraries per lane of HiSeq 4000 (60-100M reads per sample).
The number of mutant reads at MRD-filtered loci pre- and post-tumour weighting are shown in
The aggregation of signal across trinucleotide contexts as opposed to calling individual loci enables INVAR to potentially be generalised to plasma sequencing data without a priori tumour knowledge. This may have applicability in patients without tumour sequencing available, though the trade-off would be expected to be a lower sensitivity and reduced ability to quantify ctDNA levels due to the abundance of loci that would never contribute any true mutant signal.
Initially, we used TAPAS data and applied error-suppression with a minimum family size of 2, as before. Next, all bases with >=50 read families in the data were identified, and the mutant signal at each was determined at each position.
In order to focus on mutant signal arising only from ctDNA, the top 100 frequently mutated genes in public exomes were excluded (Shyr et al., 2014), as was the mitochondrial chromosome and recurrently mutated families of genes identified from Shyr et al. (2014; Supplementary Methods).
INVAR was carried out across all bases with sufficient families on a spike-in dilution experiment. Following locus blacklisting (i.e. filtering out of certain loci on the basis of having a higher locus-specific error rate) and after applying an MRD-filter (for 1F+1R MRD signal only), we show preliminary evidence for the use of INVAR in an untargeted manner (
Patient cohort. Samples were collected from patients enrolled on the MelResist (REC 11/NE/0312), AVAST-M (REC 07/Q1606/15, ISRCTN81261306)30 and LUCID (REC 14/WM/1072) studies. Consent to enter the studies was taken by a research/specialist nurse or clinician who was fully trained regarding the research. MelResist is a translational study of response and resistance mechanisms to systemic therapies of melanoma, including BRAF targeted therapy and immunotherapy, in patients with stage IV melanoma. AVAST-M is a randomised control trial which assessed the efficacy of bevacizumab in patients with stage IIB-III melanoma at risk of relapse following surgery; only patients from the observation arm were selected for this analysis. LUCID is a prospective and observational study of stage I-IIIB non-small cell lung cancer patients (NSCLC) who are planning to undergo radical treatment (surgery or radiotherapy+/− chemotherapy) with curative intent. The Cambridge Cancer Trials Unit-Cancer Theme coordinated all studies, and demographics and clinical outcomes were collected prospectively.
Sample collection and processing. Fresh frozen tumour biopsies prior to treatment were collected from patients with stage IV cutaneous melanoma. Formalin-fixed paraffin-embedded (FFPE) tumour tissue was obtained for the AVAST-M and LUCID (from surgery) trials. For patients on the AVAST-M study, plasma samples were collected within 12 weeks of tumour resection, with a subsequent sample after 3 months, where available. Patients on the LUCID study had one plasma and matched buffy coat sample taken pre-surgery. Longitudinal samples were collected during treatment of patients with stage IV melanoma as part of the MelResist study. Peripheral blood samples were collected at each clinic visit in S-Monovette 9 mL EDTA tubes. For plasma collection, samples were centrifuged at 1600 g for 10 minutes within an hour of the blood draw, and then an additional centrifugation of 20,000 g for 10 minutes was carried out. All aliquots were stored at −80° C.
Tissue and plasma extraction and quantification. FFPE samples were sectioned into up to 8 μm sections, and one H&E stained slide was generated, which was outlined for tumour regions by a histopathologist. Marked tumour regions were macrodissected, and DNA extraction was performed using the QIAamp DNA FFPE Tissue Kit using the standard protocol, except with incubation at 56° C. overnight and 500 rpm agitation on a heat block. DNA was eluted twice using 20 μL ATE buffer each time with centrifugation at full speed. Following extraction, DNA repair was performed using the NEBNext® FFPE DNA Repair Mix as per the manufacturer's protocol. Fresh frozen tissue biopsies were first homogenised prior to DNA extraction, which was performed as follows: up to 30 mg of each fresh frozen tissue biopsy sample was combined with 600 μL RLT buffer, then placed in a Precellys CD14 tube (Bertin Technologies) and homogenised at 6,500 rpm for two bursts of 20 seconds separated by 5 seconds.
Subsequently, the Qiagen AllPrep extraction kit as per the manufacturer's protocol.
Genomic DNA was extracted from up to 1 mL whole blood or buffy coat using the Gentra Puregene Blood Kit (Qiagen) as per the manufacturer's protocol. Samples were eluted in two rounds of 70 μL buffer AE and incubated for 3 minutes before centrifugation. Up to 4 mL of plasma was extracted using the QIAsymphony (Qiagen) with a QIAamp protocol. DNA was eluted in 90 μl elution buffer and stored at −80° C. Plasma samples were extracted using the QIAsymphony instrument (Qiagen) using the 2-4 mL QIAamp protocol. For each QIAsymphony batch, 24 samples were extracted, which included a positive and negative control.
Following extraction of fresh frozen, FFPE and genomic DNA, eluted DNA concentration was quantified using a Qubit fluorimeter with a dsDNA broad range assay (ThermoFisher Scientific). To quantify cell-free DNA concentration of plasma DNA eluates, digital PCR was carried out using a Biomark HD (Fluidigm) with a Taq-man probe for the housekeeping gene RPP30 (Sigma Aldrich). 55 PCR cycles were used. The RPP30 assay was 65 bp in length. The estimated number of RPP30 DNA copies per μl of eluate was used to determine the cell-free DNA concentration in the original sample.
Tumour library preparation. FFPE tumour tissue DNA samples (up to 150 ng) and buffy coat DNA samples (75 ng) were sheared to a length of 150 bp, using the Covaris LE 220 (Covaris, Mass., USA). The standard Covaris protocol for a final fragment length of 150 bp and an input volume of 15 μl using the 8 microTUBE-15 AFA Beads Strip V2 was used. After the shearing, the fragmentation pattern was verified using a Bioanalyser (Agilent).
Sequencing libraries were prepared using the ThruPLEX DNA-seq kit (Rubicon). 100 ng and 50 ng sheared tumour and buffy coat DNA, respectively, were used and the protocol was carried out according to the manufacturer's instructions. The number of amplification cycles was varied during library preparation according to the manufacturer's recommendations. Library concentration was determined using qPCR with the Illumina/ROX low Library Quantification kit (Roche). Library fragment sizes were determined using a Bioanalyser (Agilent). After library preparation, exome capture was performed with The TruSeq Exome Library Kit (Illumina), using a 45 Mbp exome baitset. Three libraries were multiplexed in one capture reaction and 250 ng of each library was used as input. For compatibility with ThruPLEX libraries, the protocol was altered by adding 1 μl of i5 and i7 TruSeq HT xGen universal blocking oligos (IDT) during each hybridisation step. To compensate for the increased hybridisation volume, the volume of CT3 buffer was adjusted to 51 μl. Two rounds of hybridisations were carried out, each lasting for 24 hours. Library QC was performed using qPCR and Bioanalyser, as above. Samples were multiplexed and sequenced with a HiSeq 4000 (Illumina).
Fresh frozen tumour biopsies and matched buffy coat library preparation was performed as described by Varela et al.31 using the SureSelectXT Human All Exon 50 Mb (Agilent) bait set. Samples were multiplexed and sequenced with a HiSeq 2000 (Illumina).
Tumour mutation calling. For fresh frozen tumour biopsies, mutation calling was performed as described by Varela et al.31. For FFPE tumour biopsies, mutation calling was performed with Mutect2 with the default settings: --cosmic v77/cosmic.vcf and --dbsnp v147/dbsnp.vcf. To maximise the number of mutations retained, variants achieving Mutect2 pass (LUCID and AVAST-M samples) OR tumour LOD >5.3 were retained (AVAST-M samples). Mutation calls were filtered as follows:
In addition, for FFPE data in the melanoma cohort, the filter for C/A errors proposed by Costello et al.32 was applied to suppress C/A artefacts. As a result, we generated patient-specific mutation lists for 64 patients with stage II-IV melanoma and stage I-IIIA lung cancer. A median of 625 (IQR 411-1076) and 388 (IQR 230-600) patient-specific mutations were identified per patient with melanoma and lung cancer, respectively (
Plasma library preparation. Cell-free DNA samples were vacuum concentrated at 30° C. using a SpeedVac (ThemoFisher) prior to library preparation where required. The median input into the library was 1652 haploid genomes (IQR 900-3013). Whole genome library preparation for plasma cell-free DNA was performed using the Rubicon ThruPLEX Tag-Seq kit. The number of PCR amplification cycles during the ThruPLEX protocol was varied between 7-15 cycles, as recommended by the manufacturer. Following amplification and sample barcoding, libraries were purified using AMPure XP beads (Beckman Coulter) at a 1:1 ratio. Library concentration was determined using the Illumina/ROX low Library Quantification kit (Roche). Library fragment sizes were determined using a Bioanalyser (Agilent). For the stage IV melanoma cohort, library preparation and sequencing was run in duplicate to assess the technical reproducibility of the experimental and computational method, showing a correlation between IMAF values generated by the INVAR pipeline of 0.97 (Pearson's r, p-value <2.2×10−16). For the early-stage cohorts, input cell-free DNA material was not split and was instead prepared and sequenced as a single sample per time point.
Custom hybrid-capture panel design and plasma sequencing. Following mutation calling, custom hybrid-capture sequencing panels were designed using Agilent SureDesign software. Between 5 and 20 patients were grouped together per panel in this implementation. Baits were designed with 4-5× density and balanced boosting for melanoma patients and 1× density and balanced boosting for lung cancer patients. 95.5% of the variants had baits successfully designed; bait design was not reattempted for loci that had failed. Custom panels ranged in size between 1.26-2.14 Mb with 120 bp RNA baits. For each panel, mutation classes and tumour allele fractions are shown in
Libraries were captured either in single or 3-plex (to a total of 1000 ng capture input) using the Agilent SureSelectXT protocol, with the addition of i5 and i7 blocking oligos (IDT) as recommended by the manufacturer for compatibility with ThruPLEX libraries33. Custom Agilent SureSelectXT baits were used, with 13 cycles of post-capture amplification. Post-capture libraries were purified with AMPure XP beads at a 1:1.8 ratio, then were quantified and library fragment size was determined using a Bioanalyser (Agilent).
Exome capture sequencing of plasma. For exome sequencing of plasma, the Illumina TruSeq Exome capture protocol was followed. Libraries generated using the Rubicon ThruPLEX protocol (as above) were pooled in 3-plex, with 250 ng input for each library. Libraries underwent two rounds of hybridisation and capture in accordance with the protocol, with the addition of i5 and i7 blocking oligos (IDT) as recommended by the manufacturer for compatibility with ThruPLEX libraries. Following target enrichment, products were amplified with 8 rounds of PCR and purified using AMPure XP beads prior to QC.
Plasma sequencing data processing. Cutadapt v1.9.1 was used to remove known 5′ and 3′ adaptor sequences specified in a separate FASTA of adaptor sequences. Trimmed FASTQ files were aligned to the UCSC hg19 genome using BWA-mem v0.7.13 with a seed length of 19. Error-suppression was carried out on ThruPLEX Tag-seq library BAM files using CONNOR34. The consensus frequency threshold -f was set as 0.9 (90%), and the minimum family size threshold -s was varied between 2 and 5 for characterisation of error rates. For custom capture and exome sequencing data, a minimum family size of 2 was used. For sWGS and bloodspot analysis, a minimum family size of 1 was used.
To leverage signal across multiple time points, error-suppressed BAM files could be combined using ‘samtools view -ubS-| samtools sort -’ prior to further data processing. In the early-stage melanoma cohort (AVAST-M), where samples were available at both 3 and 6 month time points post-surgery, BAM files were merged prior to analysis.
Low-depth whole-genome sequencing of plasma. For WGS, 30 libraries were sequenced per lane of HiSeq 4000, achieving a median of 0.6× deduplicated coverage per sample. For these libraries, since the number of informative reads (IR) would limit sensitivity before background errors would become limiting, we used error-suppression with family size 1 for this particular setting. Error rates per trinucleotide were compared between WGS and custom hybrid-capture sequencing data for family size 1, showing a Pearson r of 0.91. WGS data underwent data processing (Supplementary Methods) except the minimum depth at a locus was set to 1, and patient-specific outlier-suppression (Supplementary Methods) was not used because loci with signal vs. loci without signal would only give allele fractions of 0 or 1 given a depth of 0.6×.
Cell-free DNA extraction from dried blood spots. 50 μl fresh (or thawed frozen) whole blood was collected from patients on the MelResist study on to Whatman™ FTA™ Classic Cards, and was allowed to air dry. 50 μl fresh whole blood was obtained from an ovarian cancer xenograft mouse model immediately following its sacrifice, and similarly applied Whatman™ FTA™ Classic Cards, and allowed to air dry. Blood spot card samples were stored at room temperature inside a re-sealable plastic bag. DNA was extracted from the card using the QIAamp DNA Investigator kit, using the manufacturer's recommended extraction protocol for FTA and Guthrie cards, which is conventionally used for assessment for inherited genetic conditions in neonates from gDNA. Three 3 mm punches were made from the blood spot, and carrier RNA was added to Buffer AL as per the manufacturer's recommendation. Blood spot DNA (which we hypothesised contained both cell-free DNA and gDNA) was eluted in 25 μl water, which was reapplied to the membrane and re-eluted.
Size-selection and library preparation of blood spot cell-free DNA. Blood spot DNA eluates contain a low concentration of cell-free DNA, among a large background of gDNA (
A right-sided size-selection was performed on DNA eluates prior to library preparation using AMPure XP beads (Beckman Coulter) in order to remove long gDNA fragments. For this purpose, we adapted a published protocol for right sided size selection that is conventionally used for DNA library size-selection prior to Next Generation Sequencing36. Following optimisation of bead:sample ratios for cell-free DNA fragment sizes, we used a bead:sample ratio of 1:1 to remove contaminating gDNA. The supernatant was retained as part of the right-sided selection protocol. A second size-selection step used a 3:1 bead:sample ratio to capture all remaining fragments, and the size-selected DNA was eluted in 20 μl water. Blood spot eluates were concentrated to 10 ul volume using a vacuum concentrator (SpeedVac). Next, Rubicon Tag-Seq library preparation was carried out, and libraries underwent QC using Bioanalyser (Agilent) and qPCR (as described above). Libraries were submitted for whole-genome sequencing on a HiSeq4000 (Illumina) and the INVAR analysis pipeline was used (Supplementary Methods).
Survival analysis for resected stage II-III melanoma cohort. Disease-free interval (DFI), and overall survival were calculated from the date of randomisation of the AVAST-M trial to the date of first recurrence or date of death, respectively9. Kaplan-Meier analysis was used to generate survival curves for differences between DFI and OS in patients with detected ctDNA vs. non-detected levels and compared using a Cox proportional hazards model to obtain hazard ratios and 95% CIs.
Imaging. CT imaging was acquired as part of the standard of care from each patient of the stage IV melanoma cohort and was examined retrospectively. Slice thickness was 5 mm in all cases. All lesions with a largest diameter greater than −5 mm were outlined slice by slice on the CT images by an experienced operator, under the guidance of a radiologist, using custom software written in MATLAB (Mathworks, Natick, Mass.). The outlines were subsequently imported into the LIFEx software37 in NifTI format for processing. Tumour volume was then reported by LIFEx as an output parameter from its texture based processing module.
Circulating tumour DNA (ctDNA) can be robustly detected in plasma when multiple copies are present; however, when samples have few copies of tumour DNA, analysis of individual mutation loci can result in false negatives even if assays have perfect analytical performance due to sampling noise (
Detection of ctDNA is limited by the amount of DNA, which we quantify as the number of haploid genomes analysed (hGA). In terms of sequencing data, hGA is equivalent to the average unique sequencing coverage. In methods such as shallow whole-genome sequencing (sWGS), often <1 hGA of DNA is analysed (less than 1× coverage), and although this is often generated from nanogram (ng) amounts of DNA, it could in principle be generated from picograms of DNA. Other methods generate many thousands-fold sequencing depth, which may represent duplicate reads of the same molecules if DNA input is few ng or less. Another determinant of analytical sensitivity is the number of tumour-mutated loci that are analysed2,5-7. Sensitivity for detecting ctDNA is limited by the total number of ‘informative reads’ (IR), which we define as the sum of all reads covering loci with patient-specific mutations. This is equivalent to the product of the number of mutations and the average unique depth (across the mutated loci). We therefore plot these two variables in a two-dimensional space (
To obtain information from a large number of mutations per patient, we sequenced plasma DNA using custom capture panels, whole-exome sequencing (WES) or whole-genome sequencing (WGS). In analysing sequencing data, ctDNA detection algorithms have previously relied on identification of individual mutations, which uses limited information inefficiently: any signal that does not pass a mutation calling threshold is discarded and lost. Studies have highlighted the potential advantages of aggregating signal across multiple loci to detect DNA from transplanted organs11 or diluted tumour DNA5. In cancer monitoring, multiple mutations per patient have been previously analysed3,5,6,12,13, but detection was performed for each mutation separately. To efficiently use sequencing information from plasma we developed INtegration of VAriant Reads (INVAR). INVAR uses prior information from tumour sequencing to guide analysis and aggregate signal across 102-104 mutated loci in the patient's cancer (
To identify patient-specific mutations, we performed tumour sequencing of 45 patients with stage II-IV melanoma and 19 patients with stage I-IIIA NSCLC. After identifying tumour mutations (Methods), we generated patient-specific mutation lists (
Increasing the number of informative reads addresses sampling error, by either increasing the input (hGA) or the mutations analysed. To reduce likelihood of false positive detection at high IR, the background error must be lower than the reciprocal of IR. As part of the INVAR workflow (
Previous studies have shown relationships between tumour allele fraction and plasma allele fraction13,15, as well as showing size differences between mutant and wild-type cell-free DNA fragments16-18. To effectively use sequencing information, INVAR enriches for ctDNA signal through probability weighting based on ctDNA fragment sizes and the tumour allele fraction of each mutation locus (
14.2 Analytical Performance in Positive and Negative Controls
We evaluated the analytical performance of INVAR by analysis of sequencing from a custom capture panel, in a dilution series of plasma from one melanoma patient (stage IV), for whom we had identified 5,073 mutations through exome sequencing (Supplementary Methods), diluted into plasma from two healthy control volunteers to expected IMAF as low as 3.6×10−7, and analysed in replicates. Without error-suppression, the lowest dilution detected with analytical specificity >0.85 (Methods) had an expected ctDNA concentration of 3.6×10−5. At this concentration, 2/2 replicates were detected at average IMAF of 4.7×10−5 (
The correlation between IMAF and the expected mutant fraction was 0.98 (Pearson's r, p<2.2×10−16r
We defined the analytical specificity using sequencing data from plasma DNA of patients, using non-matched mutation lists (
14.3 Application of INVAR to Detect ctDNA in Plasma of Cancer Patients
We applied INVAR to sequencing data generated using custom capture panels from 125 plasma samples from 47 stage II-IV melanoma patients, and 19 plasma samples from 19 stage I-IIIA NSCLC patients. We analysed a median of 625 mutations per patient with melanoma and 388 mutations per patient with stage I-IIIA NSCLC, resulting in up to 2.9×106 IR per sample (median 1.7×105 IR), thus analysing orders of magnitude more cell-free DNA fragments compared to methods that analyse individual or few loci (
A small number of samples had <20,000 IR, thus not achieving the high sensitivity that INVAR can in principle produce. When implementing INVAR in future practice, we suggest that such cases where ctDNA is not detected with a low IR be defined as technical failures and would be either repeated with larger DNA input/more sequencing, or by re-analysis of the tumour and normal DNA from that patient by wider-scale sequencing such as WGS (
14.4 ctDNA Monitoring to Parts Per Million and Fractions of a Cellular Genome
We detected ctDNA and quantified its levels, as indicated by IMAF values, ranging from 2.5×10−6 to 0.25 (
In patients with metastatic melanoma, IMAF showed a correlation of 0.8 with imaging (Pearson's r, P=6.7×10−10,
14.5 ctDNA Detection in Early-Stage NSCLC
We tested ctDNA detection by INVAR in plasma samples collected pre-treatment from 19 patients with newly diagnosed stage I-IIIA NSCLC (consisting of 11, 6 and 2 patients with stage I/II/IIIA, respectively). In two samples, ctDNA was not detected, but fewer than 20,000 IR were analysed (
To test INVAR in the residual disease setting, we analysed samples from 38 patients with resected stage II-III melanoma recruited in the UK AVAST-M trial (
14.7 Assessing Detection Rates with Varying IR
Using IMAF values from the clinical samples, we estimated the expected detection rates for different cohorts of patients with limited numbers of IRs, and fitted a linear model (R2=0.95) to predict the IRs that would be required to achieve different detection rates. In stage IV melanoma patients at baseline time points, ctDNA was detected in 100% of cases using 105 IR (
14.8 Sensitive Detection of ctDNA from WES and WGS
Patient-specific capture panels allow deep sequencing of patient-specific mutation lists at lower sequencing costs, but adds a time-consuming step. We hypothesised that INVAR can be leveraged to achieve increased sensitivity by aggregating informative reads also when applied to standardised workflows such as whole exome or genome sequencing. This can allow sequencing of tumour-normal material to occur in parallel to plasma sequencing, and the resulting tumour-normal data can be used for INVAR analysis on sequencing data generated from plasma cell-free DNA (
We hypothesised that ctDNA could be detected and quantified with INVAR from even smaller amounts of input data. We performed whole-genome sequencing on libraries from cell-free DNA from longitudinal plasma samples from a subset of six patients with stage IV melanoma to a mean depth of 0.6× (
These results demonstrated detection of ctDNA from untargeted sequencing data with <1 hGA and suggest that with a sufficiently large number of tumour-specific mutations, INVAR may yield high sensitivity for ctDNA detection even with minute amounts of DNA input.
14.9 Detection of ctDNA from Dried Blood Spots
We next hypothesised that by integrating mutant reads across the genome, ctDNA could be detected from limited sequencing data generated from few copies of the genome, extracted from a dried blood spot (from a single drop of blood with volume of 50 μL). Real-time PCR has previously been used to carry out fetal RHD genotyping and HIV detection using maternal dried blood spots20,21, though NGS of cell-free DNA from blood spots has not been previously described. Generating a cell-free DNA sequencing library from a blood spot is challenging due to the low number of cell-free DNA copies present, and due to the abundance of long genomic DNA (gDNA) fragments released by blood cells (
The size distribution of the DNA fragments sequenced from the blood-spot was similar to that obtained from cell-free DNA of plasma samples (
Aside from clinical utility in humans, analysis of minute amounts of blood may facilitate longitudinal ctDNA monitoring from other organisms or models, such as rodents23. Using an orthotopically xenografted ovarian tumour mouse model, 50 μL of whole blood was sampled using a dried blood spot card, and a sequencing library prepared and sequenced with sWGS (Methods). Upon alignment of sequencing reads, both human genome (tumour-derived) and mouse genome (wild-type) reads were observed, with characteristic fragmentation patterns of mutant and wild-type cell-free DNA (
We estimate the potential sensitivity of ctDNA in dried blood spots (50 μL volume) using known mutation rates for different cancer types24. If patient-specific mutation lists were generated from WGS of tumour and normal DNA from each patient (rather than WES as was used in this study), larger mutations lists for every patient would be generated. This would allow WGS data from blood spots to generate 1-2 orders of magnitude greater IR per sample, and correspondingly lower limits for ctDNA detection (
Integration of variant reads provides a method to overcome the inherent limitations of sampling noise to detect ctDNA in samples that contain much fewer than one copy of the cancer genome, by combining signal across multiple mutations that have been identified in the patient's tumour (
We showed that INVAR can be applied flexibly to NGS data generated using patient-specific capture panels (
The INVAR pipeline takes error-suppressed BAM files, a BED file of patient-specific loci, and a CSV file indicating the tumour allele fraction of each mutation and which patient it belongs to. It is optimised for a cluster running Slurm. The workflow is shown in
SAMtools mpileup 1.3.1 was used at patient-specific loci based on a BED file of mutations, with the following settings: --ff UNMAP, -q 40 (mapping quality), -Q 20 (base quality), -x, --d 10,000, then multiallelic calls were split using BCFtools 1.3.1. Next, all TSV files were annotated with 1,000 Genomes SNP data, COSMIC data, and trinucleotide context using a custom Python script. Output files were then concatenated, compressed, and read into R. First, based on prior knowledge from tumour sequencing data, all loci were annotated per patient with being either: patient-specific (present in patient's tumour) or non-patient-specific (not present in patient's tumour, or individual does not have cancer). Since each non-patient-specific sample contains the loci from multiple patients, every non-patient-specific sample may control for all other patients analysed with the same sequencing panel or method (excluding loci that are shared between individuals).
The following filters were applied to INVAR data on a per locus basis:
After data filtering, data was annotated with both locus noise filter and trinucleotide error rate. Since the locus noise filter is limited by the number of control samples and cfDNA molecules at that locus, we also assessed trinucleotide error rate. Trinucleotide error rates were determined from the region up to 10 bp either side of every patient-specific locus (excluding patient-specific locus itself), and data was pooled by trinucleotide context. After pooling data in this manner, a median of 3.0×108 informative reads (or deduplicated reads) per trinucleotide context were analysed. Trinucleotide error rate was calculated as a mismatch rate for each specific mutation context. If a trinucleotide context had zero mutant deduplicated reads, the error rate was set to the reciprocal of the number of IR/deduplicated reads in that context.
In addition, each data point was annotated with the cfDNA fragment size of that read using a custom Python script. Then, to eliminate outlier signal that was not consistent with the remainder of that patient's loci, we performed patient-specific outlier suppression (
Patient-specific sequencing data consists of informative reads at multiple known patient-specific loci, providing the opportunity to compare mutant allele fractions across loci as a means of error-suppression. The distribution of signal across loci potentially allows for the identification of noisy loci not consistent with the overall signal distribution. Each locus was tested for the probability of having observed mutant reads given the average signal across all loci (
Based on the ctDNA level of the sample, the binomial probability of observing each individual locus given the IMAF of that sample was calculated. Loci with a Bonferroni corrected P-value <0.05 (corrected for the number of loci interrogated) were excluded in that sample, thereby suppressing outliers. As a result of outlier-suppression, background noise was reduced to 33% control samples, while retaining 96.1% of signal in patient samples (
We developed a statistical method to model the number of mutant reads at multiple patient-specific loci, incorporating prior information available from patient-specific sequencing, such as the background error of the trinucleotide context, the tumour allele fraction at the locus, and fragment length. This approach aggregates signal across multiple patient-specific mutations following error-suppression. For each locus, we test the significance of the number of mutant reads given the trinucleotide error rate of that context. Trinucleotide error rates were used instead of locus-specific error rates in order to determine a more accurate estimation of background error rates to 10−7 (
Tumour allele fractions and trinucleotide error rates were considered as follows: Denote AFi as the tumour mutant allele fraction at locus i, ei as the background error in the context of locus i, and let p be an estimate of ctDNA content in that sample for the INVAR algorithm. A random read at locus i can be observed to be mutant either if it arose from a mutant molecule, or an incorrectly sequenced wild type DNA molecule. This occurs with probability qi:
q
i
=AF
i·(1−ei)·p+(1−AFi)·ei·p+ei·(1−p) (1)
Testing for the presence of ctDNA is now equivalent to testing the statistical hypothesis H0:p=0. Assuming the number of observed mutant reads is independent between loci, the following likelihood function can be produced:
where Mij is the indicator for a mutation in read j of locus i, and Ri is the number of reads in locus i. The above method allows weighting of signal by tumour allele fraction, which we confirm influences plasma mutation representation in patient samples with early stage and advanced disease (
Each sequencing read provides fragment size information (
Assuming that given the read length and mutation status are independent given the source of the read (mutant or wild-type DNA), we can factor the likelihood as follows:
where zij is the indicator that read j of locus i came from ctDNA, pk(lij)=P(lij|zij=k), and gi=AFi·(1−ei)+(1−AFi)·ei. The above method weights the signal based on both fragment length of mutant and wild-type reads, though in this implementation of INVAR, we set the weight of all wild-type size bins to be equal, thereby neglecting size information from wild-type reads.
Lastly, a score is generated for each sample through aggregation of signal across all patient-specific loci in that sample using the Generalized Likelihood Ratio test (GLRT). The GLRT directly compares the likelihood under the null hypothesis against the likelihood under the maximum likelihood estimate of p:
The higher the value of the likelihood ratio, the greater evidence for ctDNA presence in a sample. Classification of samples was performed based on comparison of likelihood ratios between patient and control samples.
Other patients were used to control for one another at non-shared loci (
In order to determine a threshold for the likelihood ratio (LR) based on controls accurately, reads from each control sample were iteratively resampled with replacement 10 times, and the GLRT script was run. To minimise the risk of any patient-specific contamination of signal at non-patient-specific control loci (through de novo mutations overlapping with patient-specific sites), only samples with patient-specific IMAF <1% were used as controls for determination of the cut-point. Based on the LR distribution in patient controls and patient samples, the cut-off for LR was determined for each cohort using the ‘OptimalCutpoints’ package in R3, maximising sensitivity and specificity using the ‘MaxSnSp’ setting. Based on the LRs per cohort, an analytical specificity was determined for each cohort (
26 healthy individuals' cfDNA from plasma were analysed using the stage IV melanoma and stage I-IIIA NSCLC custom capture panels. These samples were treated as ‘patient’ samples, and so had no influence on the filters in the pipeline, and were not used for the determination of LR thresholds. After determination of the LR thresholds (described above), the LRs from healthy individuals' samples were assessed for false positive detection of ctDNA. For each of these cohorts, the clinical specificity values in healthy individuals were determined (
Estimation of ctDNA Content Per Sample for Likelihood Ratio Determination
In this section we derive an Expectation Maximization (EM) algorithm to estimate p as part of the INVAR method. If we treat the tumour of origin zij as a latent variable, and assume that it is known, the joint likelihood of Z, M (mij is the indicator for a mutation in read j of locus i), L (lij is the length of read i of locus j), AF (AFi is the tumour allele fraction at locus i), e (ei is the background error in the context of locus i) can be written as:
Where gi=AFi·(1−ei)+(1−AFi)·ei. The log-likelihood is linear in zij, so taking the expectation of the likelihood amounts simply to replacing the zij with their expectation at stage l, zijl=E(zij|mij,lij,pl), where pl is the best estimate of p at iteration l. We can thus use EM to find a maximum likelihood estimate for p, by iteratively maximizing the likelihood with respect to p, and taking the expectation of the likelihood with respect to zij. An estimate for pl is obtained by taking the derivative with respect to pl and equating it to zero:
The above is simply the expected proportion of reads from ctDNA at stage l. Bayes' theorem can be used to compute 4:
By substituting the respective probabilities we obtain:
The algorithm proceeds by alternating the maximization of p, and the expectation of the zij.
Size-weighting with INVAR depends on first having a known distribution of sizes of mutant and wild-type reads against which to perform weighting. In order to estimate the read length distribution with the greatest accuracy, we used all wild type and mutant reads from all samples in that cohort, leaving out the one sample being tested, and we used kernel density estimation to smooth the respective probabilities.
The size distributions from each of the studied cohorts are shown in
To estimate the probability that a read is of length l, given that the cell of origin of wild type, P(L=l|z=0), we used all of the wild-type reads from each pooled dataset. For both of the data sets, we used the R function “density”, with a Gaussian kernel, to smooth the estimated probabilities, and obtained a density estimate {circumflex over (f)}(l|Z=z). Finally, to estimate P(L=l|z=z), we integrated the respective density:
P(L=l|Z=z)=∫l−0.5l+0.5{circumflex over (f)}(t|Z=z)dt
Smoothing the size distribution estimates is important in datasets where data is sparse to avoid assigning too large a weight to any given mutant fragment.
The number of informative reads (IR) for a sample is the product of the number of mutations targeted (i.e. length of the mutation list) and the number of haploid genomes analysed by sequencing (hGA, equivalent to the deduplicated coverage following read-collapsing). Thus, the limit of detection for every sample can be calculated based on 1/IR (with adjustment for sampling mutant molecules based on binomial probabilities). For non-detected samples, the 1/IR value provides an estimate for the upper limit of ctDNA in that sample; this allows quantification of samples even if no mutant molecules are present, and is utilised in
To quantify ctDNA across multiple mutation loci, we calculated an ‘integrated mutant allele fraction’, as follows:
Plasma DNA from one patient with a total of 5,073 patient-specific variants was serially diluted 10-fold each step in a pool of plasma cfDNA from 11 healthy individuals (Seralab) to give a dilution series spanning 1-100,000×. Library preparation was performed, as described in Methods, with 50 ng input per dilution. In order to interrogate a sufficiently large number of molecules in the dilution series to assess sensitivity, the lowest dilution (100,000×) was generated in triplicate. The healthy control cfDNA pools were included as control samples for the determination of locus error rate to identify and exclude potential SNP loci (
Estimated Detection Rates with Fewer Informative Reads
Based on the IMAFs of detected samples, detection rates can be estimated where fewer IR were achieved with a perfectly sensitive assay. For a given number of IR (r), the 95% limit of detection for ctDNA (p) can be determined as follows:
p=1−elog(1−0.95)/r
Thus, for each entry in a vector of IR values (102, 103 . . . 107), the detection rates for cancer were calculated per cohort, and are plotted in
For all cohorts outlined below we applied an INVAR score threshold of 0.1 for ctDNA detection. Samples below that threshold are shown as not detected (ND). Samples with an INVAR score lower than 0.1 and a total IR of less than 20,000 were classified as not evaluable (not shown in the plot).
A total of 90 samples was analysed with INVAR. 4 cases were classified as “not evaluable” (not shown) as they did not yield enough informative reads. Out of the remaining 86 samples, ctDNA was detected in 60 samples, yielding to a total detection rate of 70%. Currently we are still blinded to the stage of these patients but the cohort is heavily biased towards early stage patients where 60% are stage I patients (33% with stage IA and 27% with stage IB).
INVAR was applied to 39 plasma samples, of which seven were not evaluable (not shown). Of the remaining 32 samples, 22 reached an INVAR score greater than 0.1 (69%). INVAR was also applied to urine samples. Of the 23 analysed samples, six were not evaluable (not shown). Of the remaining 17 samples, six reached an INVAR score greater than 0.1 (35%). Samples stem from different disease subtypes.
Urine samples were collected prior to surgery (mean 8.6, range 0-35 days pre-surgery). From the same urine sample, we isolated both urine supernatant (USN) and urine cell pellet (UCP), as described below.
30-50 ml urine was collected in a 50 ml falcon tube, and 0.5M EDTA was added within an hour of collection (pH 8.0; 600 μl for 30 ml, final concentration 10 mM. For larger volumes of urine, the volume of EDTA was adjusted accordingly). After gentle inversion, the sample was spun at 2,400 g for 10 minutes. Subsequently, ˜3.6 ml aliquots of supernatant were transferred into a separate cryotube. An additional 1 ml of supernatant was transferred to a separate microfuge tube and then returned to the original falcon tube containing the urine cell pellet (UCP). The tube was agitated and the remaining liquid was transferred to a sterile 2 ml microfuge tube, which was spun at 13,300 rpm for 10 minutes and the supernatant was discarded leaving a dry UCP for storage.
Urines samples were extracted using the QIAsymphony platform (Qiagen, Germany). Up to 4 ml of urine were extracted and eluted in 60 uL.
Detection of plasma and cerebrospinal fluid (CSF) in glioblastoma patients using INVAR. In CSF samples, 2 samples were excluded from analysis due to insufficient number of informative reads (not shown). All remaining 6 samples were detected at INVAR scores greater than 0.1. In plasma 11 out of 12 samples reached this INVAR threshold (92%).
Lumbar puncture was performed immediately prior to craniotomy for tumour debulking. After sterile field preparation, the thecal sac was cannulated between the L3 and L5 intervertebral spaces using a 0.61 mm gauge lumbar puncture needle, and 10 ml of CSF was removed. After collection, CSF samples were immediately placed on ice and then rapidly transferred to a pre-chilled centrifuge for processing. Samples were centrifuged at 1500 g at 4C for 10 minutes, the supernatant was removed and further centrifuged at 20,000 g for 10 minutes, and aliquoted into 2 mL microtubes for storage at −80C (Sarstedt, Germany).
Fluids were extracted using the QIAsymphony platform (Qiagen, Germany). Up to 8 mL of CSF was used for the extraction. DNA from CSF samples was eluted in 90 uL, and further concentrated down to 30 uL using a Speed-Vac concentrator (Eppendorf). The sample then followed the normal library preparation protocol that was also used for the plasma samples.
INVAR was applied to breast cancer samples. 34 out of 35 samples reached an INVAR score of 0.1 or greater and were classified as detected (97%). Tumour mutations in this study were identified across the entire genome rather than just the exonic regions. Samples were taken longitudinally from 7 patients and fluctuations in ctDNA concentration across time courses likely represent response to treatment.
The results are shown in
14. Mouliere F, Rosenfeld N. Circulating tumor-derived DNA is shorter than somatic DNA in plasma. Proc Natl Acad Sci 2015; 112(11):201501321.
Shyr C, Tarailo-Graovac M, Gottlieb M, Lee J J, van Karnebeek C, Wasserman W W. FLAGS, frequently mutated genes in public exomes. BMC Medical Genomics. 2014; 7:64. doi:10.1186/s12920-014-0064-y.
All references cited herein are incorporated herein by reference in their entirety and for all purposes to the same extent as if each individual publication or patent or patent application was specifically and individually indicated to be incorporated by reference in its entirety.
The specific embodiments described herein are offered by way of example, not by way of limitation. Any sub-titles herein are included for convenience only, and are not to be construed as limiting the disclosure in any way.
| Number | Date | Country | Kind |
|---|---|---|---|
| 1803596.4 | Mar 2018 | GB | national |
| 1819134.6 | Nov 2018 | GB | national |
| Filing Document | Filing Date | Country | Kind |
|---|---|---|---|
| PCT/EP2019/055610 | 3/6/2019 | WO | 00 |
