Deconvolution and detection of rare DNA in plasma

This application claims the priority benefit of PCT/US2017/067871 filed on Dec. 21, 2017 which claims priority benefit to U.S. Provisional Patent Application Ser. Nos. 62/437,512 and 62/438,401 filed on Dec. 21, 2016 and Dec. 22, 2016, respectively, the entire contents of which are hereby incorporated by reference.

This invention was made with government support under Grant Number R01GM097253 awarded by the National Institutes of Health. The government has certain rights in the invention.

Some embodiments described herein relate to compositions and methods for detecting a target nucleic acid in a sample. For example, some embodiments may be used for non-invasive detection of tumors or organ damage. Other embodiments may be used, for example, for detection of fetal aneuploidy at a very early stage of pregnancy.

CpG methylation in mammalian genomes is a relatively stable epigenetic modification, which can be transmitted across cell division (Wigler et al. (1981)) through DNMT1, and dynamically established, or removed by DNMT3 A/B and TET proteins. Due to the processivity of some of these enzymes, physically adjacent CpG sites on the same DNA molecules can share similar methylation status, although discordant CpG methylation has also been observed, especially in cancer cells. The theoretical framework of linkage disequilibrium (Slatkin (2008)), which was developed to model the coordinated segregation of adjacent genetic variants on human chromosomes among human populations, can be applied to the analysis of CpG co-methylation in cell populations. A number of studies related to the concepts of methylation haplotypes, epi-alleles, or epi-haplotypes have been reported, albeit at small numbers of genomic regions or limited numbers of cell/tissue types. Recent data production efforts, especially by large consortia such as the NIH RoadMap Epigenomics project (Bernstein et al. (2010)) and the EU Blueprint Epigenome project (Jones et al. (2005)) have produced a large number of whole-genome, base-resolution bisulfite sequencing data sets for many tissue and cell types. These public data sets, in combination with additional WGBS data generated, allowed full-genome characterization of local coupled CpG methylation across the largest set of human tissue types available to date to be performed, and annotate these blocks of co-methylated CpGs as a distinct set of genomic features.

DNA methylation is cell-type specific, and the pattern can be harnessed for deconvoluting the relative cell composition of heterogeneous samples, such as different white blood cells in whole blood (Houseman et al. (2016)), fetal components in maternal cell-free DNA (Sun et al. (2015)), or circulating tumor DNA in plasma (Sun et al. (2015)). Most of these recent efforts rely on the methylation level of individual CpG sites, and are fundamentally limited by the technical noise and sensitivity in measuring single CpG methylation. Very recently, Lehmann-Werman et al. demonstrated a superior sensitivity with multi-CpG haplotypes in detecting tissue-specific signatures in circulating DNA (Lehmann-Werman et al. (2016)). The markers in that study were discovered from Infinium 450k methylation array data, which represent only a very limited fraction of the human genome.

Some embodiments described herein provide accurate tissue-of-origin mapping based on comparing the patterns and abundance of methylation haplotypes against a reference set of human reference tissues and provide accurate quantitative estimation of the cancer DNA fraction.

Some embodiments are described in the following numbered paragraphs:

1. A method for detecting the presence of one or more nucleic acids indicative of a health condition, tissue of origin, germ layer of origin or organ of origin in a mixture of nucleic acids comprising: performing methylation analysis on a sample comprising a plurality of nucleic acids; and determining whether said sample includes a plurality of methylation haplotype blocks indicative of the presence one or more nucleic acids indicative of a health condition, tissue of origin, germ layer of origin, organ of origin or any combination thereof, wherein said methylation haplotype blocks comprise a plurality of methylation sites for which the methylation status is coordinated. 2. The method of Paragraph 1, wherein said methylation analysis is performed on cell-free DNA. 3. The method of any one of Paragraphs 1 and 2, wherein said methylation analysis is performed on cell-free DNA in a blood sample. 4. The method of any one of Paragraphs 1-3 wherein said plurality of methylation haplotype blocks comprises at least 2, at least 3, at least 4, at least 5, at least 10, at least 20, at least 40, at least 50, at least 100, at least 200, at least 300, at least 400, at least 500 or more than 500 methylation haplotype blocks. 5. The method of any one of Paragraphs 1-4, wherein said health condition is a tumor. 6. The method of Paragraph 5, further comprising determining whether said sample includes a plurality of methylation haplotype blocks indicative of the presence of one or more nucleic acids indicative of a normal tissue or normal organ corresponding to the tissue or organ of origin of said tumor. 7. The method of any one of Paragraphs 1-6, wherein said health condition is fetal aneuploidy. 8. The method of any one of Paragraphs 1-7, wherein said sample is a blood sample. 9. The method of any one of Paragraphs 1-8, further comprising quantitating the level of said one or more nucleic acids indicative of a health condition, tissue of origin, organ of origin or any combination thereof in said sample. 10. The method of any one of Paragraphs 1-9, wherein said methylation analysis is performed using a technique selected from the group consisting of bisulfite methylation analysis, reduced representation bisulfite sequencing, WGBS, BSPP, micro-droplet PCR, selector probe based methods, and MeDiP. 11. The method of any one of Paragraphs 1-10, further comprising determining a methylated haplotype load or unmethylated haplotype load for each methylation haplotype block, wherein said methylation haplotype load comprises the normalized fraction of methylated haplotypes at different lengths and the unmethylated haplotype load comprises the normalized fraction of unmethylated haplotypes at different lengths. 12. The method of any one of Paragraphs 1-11, wherein said methylation haplotype blocks have an average size of 95 bp. 13. The method of any one of Paragraphs 1-12, wherein said methylation haplotype blocks have a minimum of 3 CpGs per block. 14. The method of any one of Paragraphs 1-13 further comprising quantifying the level of said plurality of methylation haplotype blocks indicative of the presence of one or more nucleic acids indicative of a health condition, tissue of origin, germ layer of origin, organ of origin or any combination thereof in said sample. 15. The method of Paragraph 6 further comprising quantifying the level of said plurality of methylation haplotype blocks indicative of the presence of a tumor in said sample and quantifying the level of said plurality of methylation haplotype blocks indicative of the presence of one or more nucleic acids indicative of a normal tissue or normal organ corresponding to the tissue or organ of origin of said tumor in said sample. 16. A method of identifying methylation haplotype blocks comprising: determining methylation haplotypes in a plurality of nucleic acid segments; combining the methylation haplotypes and calculating methylation linkage disequilibrium on the combined methylation haplotypes; and partitioning each segment into a plurality of methylation haplotype blocks, wherein said methylation haplotype blocks comprise a plurality of methylation sites for which the methylation status is coordinated. 17. The method of Paragraph 16, wherein said methylation haplotypes are determined for a portion of a genome. 18. The method of any one of Paragraphs 16 and 17, wherein said methylation haplotypes are determined across a whole genome. 19. The method of any one of Paragraphs 16-18, wherein said methylation haplotype blocks are defined as the genomic region in which the r 2 value of two adjacent CpG sites is no less than 0.5. 20. The method of any one of Paragraphs 16-19, wherein said methylation haplotype blocks are identified in nucleic acids from a tumor tissue. 21. The method of any one of Paragraphs 16-20, wherein said methylation haplotype blocks are identified in nucleic acids from a known type of tissue. 22. The method of any one of Paragraphs 16-21, wherein said methylation haplotype blocks are identified in nucleic acids from a fetus. 23. The method of any one of Paragraphs 16-22, wherein said methylation haplotype blocks are identified in nucleic acids from an embryonic stem cell. 24. The method of any one of Paragraphs 16-23, wherein said methylation haplotype blocks are identified in nucleic acids from a known germ layer. 25. The method of any one of Paragraphs 16-24, wherein said methylation haplotype blocks have an average size of 95 bp. 26. The method of any one of Paragraphs 16-25, wherein said methylation haplotype blocks have a minimum of 3 CpGs per block. 27. The method of any one of Paragraphs 16-26, wherein said methylation haplotype blocks are identified in nucleic acid regions from a Whole Genome Bisulfite Sequencing analysis. 28. The method of any one of Paragraphs 16-27, wherein said methylation haplotype blocks are identified in data sets from methylation analysis of ENCODE cell lines or tissue samples. 29. The method of any one of Paragraphs 16-28, wherein said methylation haplotype blocks are identified in data sets from methylation analysis of Infinium HumanMethylation450 BeadChip (HM450K). 30. The method of any one of Paragraphs 16-29, further comprising calculating the pairwise correlation coefficient of adjacent CpG methylation levels across different sample sets for block partitioning. 31. The method of any one of Paragraphs 16-30, further comprising determining methylation haplotype load for each methylation haplotype block, wherein said methylation haplotype load comprises the normalized fraction of methylated haplotypes at different lengths. 32. The method of any one of Paragraphs 16-30, further comprising determining unmethylated haplotype load for each methylation haplotype block, wherein said unmethylated haplotype load comprises the normalized fraction of unmethylated haplotypes at different lengths.

FIGS. 1 ( a )- 1 ( e ) . Identification and characterization of human methylation haplotype blocks (MHBs). FIG. 1 ( a ) Schematic overview of data collection, generation, and analysis. FIG. 1 ( b ) An example of MHB at the promoter of the gene APC. FIG. 1 ( c ) Smooth scatterplots of methylation linkage disequilibrium decay of adjacent CpG sites over larger distances. 500,000 adjacent CpG loci in MHB regions were randomly sampled and their corresponding r 2 values within different sets of cell types were plotted to demonstrate the differential decay characteristics. The dotted lines at high linkage (r 2 =0.9) is where stem and progenitor cells (10 WGBS samples), normal adult tissue cells (49 WGBS samples), and primary tumor (6 WGBS samples) cells were found to have 94.8%, 91.2% and 87.8% of CpGs respectively. FIG. 1 ( d ) Co-localization analysis of MHBs with known genomic features. Genome distribution (left) and CpG-island relationships (right). FIG. 1 ( e ) . Enrichment of MHBs in known genomic features. Bootstrap random sampling regions with same size for 10,000 times to estimate empirical statistical significance and enrichment factor (fold-change).

FIG. 2 . Comparison of methylation haplotype load with four metrics used in the literature. Five patterns of methylation haplotype combinations are used to illustrate the difference between methylation frequency, methylation entropy, epi-polymorphism and methylation haplotype load. Methylation haplotype load can discriminate all the five patterns while other metrics cannot.

FIGS. 3 ( a )- 3 ( b ) . Heatmaps of tissue specific MHBs. FIG. 3 ( a ) Heatmap of MHL values for tissue specific MHBs selected by MHL. FIG. 3 ( b ) Heatmap of uMHL values for tissue specific MHBs selected by uMHL.

FIG. 4 . Comparison of signal to noise ratio between average methylation frequency (AMF) and methylation haplotype load (MHL) metrics at tissue specific differentially methylated loci (Lokk et al. (2014)). For most loci, the MHL metric had a much higher signal to noise ratio than the AMF metric even though these loci were selected using AMF for another data set.

FIG. 5 . Heatmap of the MHL in plasma samples in test sets and primary cancer tissues for colon cancer markers (left) and lung cancer markers (right). Cancer markers were identified by comparison of the respective cancer tissues against background (normal plasma set aside for features selection).

FIG. 6 . Boxplots of average MHL in plasma samples in test sets and cancer tissues for colon cancer markers (left) and lung cancer markers (right). The cancer markers were significantly (Two Sample one-sided T-Test) different between the test sets of normal plasma (NP) and cancer plasma (CCP=colon cancer plasma, LCP=lung cancer plasma). The marker regions were highly methylated in tumor tissues (CCT=colon cancer tissue, LCT=lung cancer tissue).

FIG. 7 . Simulated standard curves for estimation of tumor fraction. The mean for each cancer fraction from 20 simulations plotted with standard deviations as error bars. A fitted linear model on the standard curve for colon cancer had an adjusted r 2 of 0.9621 and the fitted linear model for lung cancer had an adjusted r 2 of 0.9573.

FIG. 8 . The tumor fraction between test sets of cancer and normal plasma samples. Cancer plasma (n=30 colon cancer, n=29 lung cancer) had significantly elevated tumor fractions (Two Sample t-test with unequal variance) when compared with normal plasma (n=23) using either colon cancer markers (left) or lung cancer markers (right).

FIGS. 9 ( a )- 9 ( c ) . Results of MARS-based binary classification for test plasma data sets. The set of tissue-specific markers derived from an independent set of normal tissue data was able to segregate each class of plasma samples from the other classes in a binary fashion. FIG. 9 ( a ) Normal vs other. FIG. 9 ( b ) Colon vs other. FIG. 9 ( c ) Lung vs other.

FIG. 10 . Confusion matrix of MARS-based feature selection and PLSDA prediction of tissue-of-origin in the test plasma data set. The set of lung and colon tissue specific markers identified from an independent data set was able to correctly classify 74% of the test plasma samples.

FIGS. 11 ( a )- 11 ( c ) . Characteristics of MHB in human genome. FIG. 11 ( a ) Distribution of MHB sizes. FIG. 11 ( b ) Distribution of CpG density (CpGs/bp) in MHB regions. FIG. 11 ( c ) Co-localization of MHB with known genomics features breaking down based on CpG density. MHBs were placed in quartiles where the CpGs/bp of MHBs within quartiles were as follows (0, 0.46), (0.046, 0.096), (0.096, 0.155), and (0.155, 0.6). The 1 st quartile (MHBs with lowest CpG density) was mostly CGI shelves or shores, and was enriched for LAD, LOCK, and enhancers.

FIGS. 12 ( a )- 12 ( b ) . Validation of MHB with Illumina 450k methylation array and RRBS data. FIG. 12 ( a ) Pearson correlation coefficient (r 2 ) versus absolute LD r 2 . FIG. 12 ( b ) The Pearson correlation coefficient r 2 in RRBS and HM450K were significantly higher in overlapped MHBs with WGBS compared with the MHBs without overlapping with WGBS MHBs. IN: denotes RRBS or HM450K regions within MHB. OUT: denotes RRBS or HM450K regions beyond MHB regions.

FIG. 13 . The distribution of incidence of cancer-associated high-methylated haplotypes (caHMH) in colon cancer (CRC) and lung cancer (LC) plasma samples. Y-axis denotes the frequency of caHMH and x-axis denotes the incidence (number of samples) of the caHMH in cancer plasmas. A majority of caHMH are patient specific.

Unless defined otherwise, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art pertinent to the methods and systems described. As used herein, the following terms and phrases have the meanings ascribed to them unless specified otherwise. All patents, applications, published applications and other publications referred to herein are incorporated by reference in their entirety.

“A,” “an,” and “the” include plural referents, unless the context clearly indicates otherwise. For example, “a nucleic acid” as used herein is understood to represent one or more nucleic acids. As such, the terms “a” (or “an”), “one or more,” and “at least one” can be used interchangeably herein.

“About” is used herein to mean approximately, roughly, around, or in the region of. When the term “about” is used in conjunction with a numerical range, it modifies that range by extending the boundaries above and below the values set forth. In general, the term “about” is used herein to modify a numerical value above and below the stated value by a deviation of ±10% and preferably ±5%.

“Comprise,” “comprises,” “comprising,” “include,” “includes,” and “including” are interchangeable and not intended to be limiting. Furthermore, where the description of one or more embodiments uses the term “comprising,” those skilled in the art would understand that, in some specific instances, the embodiment or embodiments can be alternatively described using the language “consisting essentially of” and/or “consisting of.”

The following abbreviations are used throughout the application: Methylation Haplotype Block (MHB); Methylation Haplotype Load (MHL); Unmethylated Haplotype Load (uMHL); Group Specific Index (GSI); Circulating cell-free DNA (cfDNA); Reduced Representation Bisulfite Sequencing (RRBS); single-cell Reduced Representation Bisulfite Sequencing (scRRBS); Whole Genome Bisulfite Sequencing (WGBS); The Cancer Genome Atlas project (TCGA); The Encyclopedia of DNA Elements (ENCODE); Gene Expression Omnibus (GEO); Lung Cancer (LC); Colorectal or colon cancer (CRC); cancer associated High Methylation Haplotype (caHMH); tissue-specific Methylation Haplotype Block regions (tsMHB): Colorectal or colon cancer tissue (CCT); colorectal or colon cancer plasma (CCP); lung cancer tissue (LCT); lung cancer plasma (LCP); normal plasma (NP).

“Amplification” refers to any known procedure for obtaining multiple copies of a target nucleic acid or its complement, or fragments thereof. The multiple copies may be referred to as amplicons or amplification products. Amplification, in the context of fragments, refers to production of an amplified nucleic acid that contains less than the complete target nucleic acid or its complement, e.g., produced by using an amplification oligonucleotide that hybridizes to, and initiates polymerization from, an internal position of the target nucleic acid. Known amplification methods include, for example, replicase-mediated amplification, polymerase chain reaction (PCR), reverse transcription polymerase chain reaction (RT-PCR), ligase chain reaction (LCR), strand-displacement amplification (SDA), and transcription-mediated or transcription-associated amplification.

“Complementary” means that a contiguous nucleic acid base sequence is capable of hybridizing to another base sequence by standard base pairing (hydrogen bonding) between a series of complementary bases. Complementary sequences may be completely complementary (i.e. no mismatches in the nucleic acid duplex) at each position in an oligomer sequence relative to its target sequence by using standard base pairing (e.g., G:C, A:T or A:U pairing) or sequences may contain one or more positions that are not complementary by base pairing (e.g., there exists at least one mismatch or unmatched base in the nucleic acid duplex), but such sequences are sufficiently complementary because the entire oligomer sequence is capable of specifically hybridizing with its target sequence in appropriate hybridization conditions (i.e. partially complementary). Contiguous bases in an oligomer are typically at least 80%, preferably at least 90%, and more preferably completely complementary to the intended target sequence.

“Configured to” denotes an actual arrangement of a nucleic acid sequence configuration of a referenced oligonucleotide. For example, a primer that is configured to generate a specified amplicon from a target nucleic acid has a nucleic acid sequence that hybridizes to the target nucleic acid or a region thereof and can be used in an amplification reaction to generate the amplicon. Also as an example, an oligonucleotide that is configured to specifically hybridize to a target nucleic acid or a region thereof has a nucleic acid sequence that specifically hybridizes to the referenced sequence under stringent hybridization conditions.

“Configured to specifically hybridize to” means that an oligonucleotide is designed to have a nucleic acid sequence that can hybridize with a target nucleic acid or region thereof. The oligonucleotide is designed to function as a component of an assay for amplification and detection of a target nucleic acid (or a region thereof) in a sample, and therefore is designed to hybridize with a target nucleic acid (or a region thereof) in the presence of other nucleic acids that may be found in testing samples.

“Fragment” refers to a piece of contiguous nucleic acid that contains fewer nucleotides than the complete nucleic acid.

“Hybridization” or “annealing” refer to the base-pairing interaction of one nucleic acid with another nucleic acid (typically an antiparallel nucleic acid) that results in formation of a duplex or other higher-ordered structure (i.e. a hybridization complex). The primary interaction between the antiparallel nucleic acid molecules is typically base specific, e.g., A/T and G/C. It is not a requirement that two nucleic acids have 100% complementarity over their full length to achieve hybridization. Nucleic acids hybridize due to a variety of well characterized physio-chemical forces, such as hydrogen bonding, solvent exclusion, base stacking and the like. An extensive guide to the hybridization of nucleic acids is found in Tijssen (1993) Laboratory Techniques in Biochemistry and Molecular Biology—Hybridization with Nucleic Acid Probes part I chapter 2, “Overview of principles of hybridization and the strategy of nucleic acid probe assays,” (Elsevier, New York), as well as in Ausubel (Ed.) Current Protocols in Molecular Biology, Volumes I, II, and III, 1997, which is incorporated by reference.

“Nucleic acid” or “nucleic acid molecule” refers to a multimeric compound comprising two or more covalently bonded nucleosides or nucleoside analogs having nitrogenous heterocyclic bases, or base analogs, where the nucleosides are linked together by phosphodiester bonds or other linkages to form a polynucleotide. Nucleic acids include RNA, DNA, or chimeric DNA-RNA polymers or oligonucleotides, and analogs thereof. A nucleic acid backbone can be made up of a variety of linkages, including one or more of sugar-phosphodiester linkages, peptide-nucleic acid bonds, phosphorothioate linkages, methylphosphonate linkages, or combinations thereof. Sugar moieties of the nucleic acid can be ribose, deoxyribose, or similar compounds having known substitutions (e.g. 2′-methoxy substitutions and 2′-halide substitutions). Nitrogenous bases can be conventional bases (A, G, C, T, U) or analogs thereof (e.g., inosine, 5-methylisocytosine, isoguanine). A nucleic acid can comprise only conventional sugars, bases, and linkages as found in RNA and DNA, or can include conventional components and substitutions (e.g., conventional bases linked by a 2′-methoxy backbone, or a nucleic acid including a mixture of conventional bases and one or more base analogs). Nucleic acids can include “locked nucleic acids” (LNA), in which one or more nucleotide monomers have a bicyclic furanose unit locked in an RNA mimicking sugar conformation, which enhances hybridization affinity toward complementary sequences in single-stranded RNA (ssRNA), single-stranded DNA (ssDNA), or double-stranded DNA (dsDNA). Nucleic acids can include modified bases to alter the function or behavior of the nucleic acid (e.g., addition of a 3′-terminal dideoxynucleotide to block additional nucleotides from being added to the nucleic acid). Synthetic methods for making nucleic acids in vitro are well known in the art although nucleic acids can be purified from natural sources using routine techniques. Nucleic acids can be single-stranded or double-stranded.

“Primer” refers to an enzymatically extendable oligonucleotide, generally with a defined sequence that is designed to hybridize in an antiparallel manner with a complementary, primer-specific portion of a target nucleic acid. A primer can initiate the polymerization of nucleotides in a template-dependent manner to yield a nucleic acid that is complementary to the target nucleic acid when placed under suitable nucleic acid synthesis conditions (e.g. a primer annealed to a target can be extended in the presence of nucleotides and a DNA/RNA polymerase at a suitable temperature and pH). Suitable reaction conditions and reagents are known to those of ordinary skill in the art. A primer is typically single stranded for maximum efficiency in amplification, but may alternatively be double stranded. If double stranded, the primer is generally first treated to separate its strands before being used to prepare extension products. The primer generally is sufficiently long to prime the synthesis of extension products in the presence of the inducing agent (e.g. polymerase). Specific length and sequence will be dependent on the complexity of the required DNA or RNA targets, as well as on the conditions of primer use such as temperature and ionic strength. Preferably, the primer is about 5-100 nucleotides. Thus, a primer can be, e.g., 5, 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, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95 or 100 nucleotides in length. A primer does not need to have 100% complementarity with its template for primer elongation to occur; primers with less than 100% complementarity can be sufficient for hybridization and polymerase elongation to occur. A primer can be labeled if desired. The label used on a primer can be any suitable label, and can be detected by, for example, spectroscopic, photochemical, biochemical, immunochemical, chemical, or other detection means. A labeled primer therefore refers to an oligomer that hybridizes specifically to a target sequence in a nucleic acid, or in an amplified nucleic acid, under conditions that promote hybridization to allow selective detection of the target sequence.

“Sample preparation” refers to any steps or methods that prepare a sample for subsequent sequencing, amplification, and/or detection of target nucleic acids present in the sample. Sample preparation may include any known method of concentrating components, such as nucleic acids, from a larger sample volume. Sample preparation may include physical disruption and/or chemical lysis of cellular components to release intracellular components into a substantially aqueous or organic phase and removal of debris, such as by using filtration, centrifugation or adsorption. Sample preparation may include use of a nucleic acid oligonucleotide that selectively or non-specifically captures a target nucleic acid and separates it from other sample components.

“Sequencing” refers to any known procedure, method, or technology for determining the precise order of the nucleosides or nucleoside analogs of a target nucleic acid molecule, or its complement, or fragments thereof. Sequencing, in the context of fragments, refers to determining the precise order of nucleosides or nucleotides within a nucleic acid molecule that contains less bases than the complete target nucleic acid molecule e.g., determined by sequencing amplicons produced by using an amplification oligonucleotide that hybridizes to, and initiates polymerization from, an internal position of the target nucleic acid. Known sequencing methods include, for example, whole-genome sequencing as well as targeted sequencing wherein only subset of genes or regions of the genome are isolated and sequenced.

Some embodiments described herein allow a quantitative deconvolution of biological samples (for example, human plasma) that contain mixed DNA molecules, based on comparing the patterns and abundance of methylation haplotypes against a reference set of human tissues. Detection and quantification of low-abundant cancer DNA in plasma, and simultaneous mapping to the tissue of origin has been successfully demonstrated. With the high sensitivity and the capability of detecting rare species in heterogeneous samples. Some embodiments described herein are, for example, suitable for non-invasive detection of tumor or other organ damage (diabetes, stroke etc.) in a subject's plasma, urine, stool, or cerebrospinal fluid, and detection of fetal aneuploidy in maternal blood at a very early stage of pregnancy.

Some embodiments described herein provide methods of creating a database of patient specific methylation haplotypes and unmethylated haplotypes. Some embodiments described herein provide methods of preparing data of comparative methylation haplotypes and unmethylated haplotypes. Some embodiments described herein include further steps for methods of treatment of diagnosed diseases or conditions associated with methylation haplotypes and unmethylated haplotypes.

Some embodiments described herein utilize the co-methylation status of multiple adjacent CpG sites in single DNA molecules rather than the methylation status of individual loci and or their local average. Some embodiments described herein suppress the stochastic noise from the low coverage methylation data and therefore have high sensitivity and specificity in both detecting circulating tumor DNA and mapping the tissue-of-origin of the tumor. Additionally, some embodiments described herein provide a set of MHL based biomarkers and a statistical approach to the detection.

For the quantitative estimation of the level of cancer DNA fractions, a strategy based on a tissue-specific MHL sampling technique which is distinctly different from existing approaches, such as quadratic programming, was developed. Due to the low level of cancer and fetal DNA fragments in the blood (plasma), some embodiments are applicable to even low coverage data based on methylation sequencing, such as RRBS and WGBS, whereas other methods would fail due to the difficulties in dealing with missing values.

For non-invasive cancer diagnosis, current mature commercial assays in the market use bronchial fluid (EPI PROLUNG BL Reflex assay), or can only be used for specific cancer type detection (EPI PROCOLON blood assay for colon cancer) without origin-tissue prediction. Some embodiments described herein provide accurate tissue-of-origin mapping based on comparing the patterns and abundance of methylation haplotypes against a reference set of human reference tissues and provide accurate quantitative estimation of the cancer DNA fraction.

Some embodiments described herein relate to a set of markers for cancer detection and tissue-of-origin mapping based on an exhaustive search across the whole genome for methylation haplotype blocks (MHB) that have a higher level of methylation haplotype load (MHL) or higher level of unmethylated haplotype load (uMHL). Some embodiments described herein relate to a statistical approach that takes the information of these markers for the tissue-of-origin mapping and cancer detection. Compared with existing prediction models (such as tree-based: random forest), some embodiments described herein are suitable for sparse methylation sequencing data, which is typical in the clinical setting due to the low amounts of cell-free DNA that can be extracted from patients. Some embodiments described herein are also applicable to detection of fetal aneuploidy in maternal blood. Currently there is a commercial TISSUE OF ORIGIN test from Cancer Genetics Inc. which is an invasive assay that uses tissue biopsies from patients, as compared to non-invasive liquid biopsies that allow far less discomfort for the patients.

Human peripheral blood contains low levels of DNA molecules from other tissues or cell types, such as circulating cancer stem cells, cell-free DNA from apoptotic cancer cells in cancer patients, or fetal DNA in pregnant women. DNA methylation signals along with the DNA molecules are released into the blood simultaneously and these methylation signals are tissue-specific and can be applied to identify the tissue source of the DNA fragment. Methylation haplotypes provide high sensitivity detection of DNA molecules from DNA samples (WO2015/116837 and PCT/US2015/013562 incorporated herein by reference).

In some embodiments, DNA molecules are extracted from plasma samples, from for example cancer patients or healthy individuals, and the DNA methylation status of cell-free DNA molecules are assayed by bisulfite methylation sequencing (reduced representation bisulfite sequencing (RRBS), Meissner et al. (2005)). Note that alternative technologies, such as BSPP (Diep et al. (2012)), micro-droplet PCR (Komori et al. (2011)), Selector probes (Johansson et al. (2011)), or MeDiP (Papageorgiou et al. (2011)) can also potentially be used with some differences in the requirement of input materials and/or cost. Other methods for determining the methylation status of cell-free DNA may also be used.

In some embodiments, regardless of the specific sample preparation methods or sequencing platforms used, bisulfite sequencing reads (single-ends or paired-ends) are used as the input for the proposed analytical framework for detection of cancer and their tissue-of-origin. Methylation haplotypes and their abundance may be derived from the raw sequencing reads. Each haplotype represents the combination of binary methylation status (methylated or unmethylated) at multiple CpG sites of one sequencing read. A computational pipeline (implemented with Perl) was developed to deal with methylation haplotype and further derived metrics, such as methylation haplotype (successive methylation allele combination), methylation haplotype block (linkage disequilibrium regions), highly methylated haplotype (HMH), and methylation/unmethylated haplotype load (MHL/uMHL) (weighted methylation status by haplotype length).

To identify the most informative marker for detecting cancer or tissue mapping, the term “methylation haplotype block” (MHB) was defined to describe genomic regions in which the methylation status of CpG pairs within are in linkage disequilibrium and could be taken as co-methylation regions. A set of Whole Genome Bisulfite Sequencing (WGBS) data was compiled, such as from stem cell, cancer cell, and normal adult tissues, in order to exhaustively search the entire human genome. A total of 147,888 MHBs in the human genome were identified.

A metric was defined, called methylation haplotype load (MHL) for each MHB (as described in WO2015/116837), which is the normalized fraction of methylated haplotypes at different lengths to indicate the methylation status and complexity of the methylation population. MHL is superior to other metrics (such as methylation level, methylation entropy (Xie, H. et al, (2011)), and epi-polymorphism (Landan et al, (2012)) for distinguishing different methylation patterns. In addition, MHL is bounded between 0 and 1, which allows for direct comparison of different regions across many data sets without normalization.

A metric was defined, called unmethylated haplotype load (uMHL) for each MHB, which is the normalized fraction of unmethylated haplotypes at different lengths to indicate the methylation status and complexity for the unmethylated population.

A MHL matrix for a set of human tissues (61 published samples/data sets) was built encompassing all the 147,888 identified MHBs in the entire human genome. From this matrix, tissue-specific MHBs were identified. After including additional training data subsets comprising plasma from cancer patients and healthy controls, sets of 154 MHBs for binary prediction of tissue-of-origin between colon, lung, and normal (no solid tissue) and 295 MHBs for cancer types classification (colon versus lung) non-invasively in the blood were identified (Table 5(a)-(b)). The invention demonstrated binary classification AUCs of 0.856, 0.725, and 0.751 in independent test sets from normal, colon cancer, and lung cancer plasma respectively ( FIG. 9 ( a )- 9 ( c ) ). Distinguishing between colon cancer and lung cancer plasma achieved classification accuracies of 70% and 78% for test sets comprising colon cancer and lung cancer plasma respectively using cancer type classification ( FIG. 10 ).

The disclosed method was also applicable to tracking metastasized secondary tumors. Significantly higher original tissue markers (46 and 79 for two metastasized cancer samples) in the metastasized cancers (primary colon to liver and primary breast to lung) were identified with MHL measurement and a maximum marker count approach to infer its primary cancer source.

Adjacent CpG sites in mammalian genomes tend to be co-methylated due to the processivity of enzymes responsible for adding or removing the methyl group. Yet discordant methylation patterns have also been observed, and found to be related to stochastic or uncoordinated molecular processes. The invention focused on a systematic search and investigation of regions in the human genome that exhibit highly coordinated methylation.

By examining the co-methylation patterns of multiple adjacent CpG sites, termed methylation haplotypes, in single bisulfite sequencing reads, a greedy-searching strategy was applied to define blocks of tightly coupled CpG sites, called methylation haplotype blocks (MHBs), based on 53 sets of whole genome bisulfite sequencing (WGBS) data, including 43 published sets from human adult tissues, ESC, and in vitro differentiated cell lines, as well as 10 sets from human adult tissues generated in this invention. The MHBs were further validated with 101 sets of RRBS ENCODE data, and 1,274 sets of Illumina450k methylation array data from TCGA tumor and normal samples. Globally, MHBs are enriched in, but only partially overlap with, several well-known genomic features, including CpG islands, promoters, enhancers and VMRs.

To perform quantitative analysis of the MHBs, a metric called methylation haplotype load (MHL) was computed, which covers both average methylation level and methylation complexity, and therefore is more informative than average methylation level or Shannon entropy. Using a feature selection strategy, a set of tissue-specific MHBs was identified that cluster by developmental germ-layers. Interestingly, examination of these MHBs revealed two distinct mechanisms for fate commitment during development: epigenetic silencing of pluripotent genes, such as NANOG, for mesoderm induction; and epigenetic induction (or de-suppression) of lineage-specific factors for ectoderm commitment.

Furthermore, to examine MHBs in cancers and explore its clinical utility, 162 sets of RRBS data from primary tumor tissues and matched plasma from patients with lung cancer, colon cancer, or pancreatic cancer were generated, as well as plasma controls from healthy individuals. Compared with normal tissues and stem cells, primary tumor tissues exhibit a distinct methylation pattern within MHBs, related to locally disordered methylation recently discovered in chronic lymphocytic leukemia (CLL). Importantly, a subset of blocks (Table 4) was derived that can estimate the tumor content from circulating DNA in the plasma ( FIG. 8 ). Finally, prediction of tissue-of-origin was performed on plasma using tissue specific methylation haplotype blocks (Table 5).

A manuscript entitled “Identification of methylation haplotype blocks aids in deconvolution of heterogeneous tissue samples and tumor tissue-of-origin mapping from plasma DNA” by Guo et al. (Nat Genet 2017 April; 49(4):635-642), the disclosure of which is incorporated herein by reference in its entirety, including all figures, tables, supplementary figures and supplementary tables therein or related thereto.

An exhaustive search of tissue-specific methylation haplotype blocks across the entire human genome was performed, and block-level metrics including methylation haplotype load (MHL) and unmethylated haplotype load (uMHL) were proposed, for a systematic discovery of informative markers. Applying the disclosed method's analytical framework and identified markers, accurate determination of tissue origin as well as estimation of tumor load in clinical plasma samples from patients of lung cancer (LC) and colorectal cancer (CRC) ( FIG. 1 ( a ) ) was demonstrated.

Identification and Characterization of Methylation Haplotype Blocks.

To investigate the co-methylation status of adjacent CpG sites along single DNA molecules, the concept of genetic linkage disequilibrium was extended (Slatkin (2008); Shoemaker et al. (2010)) and the r 2 metric determined to quantify the degree of coupled CpG methylation among different DNA molecules of the same samples. CpG methylation status of multiple CpG sites in single- or paired-end Illumina sequencing reads were extracted to form methylation haplotypes, and pairwise “linkage disequilibrium” of CpG methylation r 2 was calculated from the abundance of different methylation haplotypes (see Methods). The full human genome was partitioned into blocks of tightly coupled CpG methylation sites, called methylation haplotype blocks (MHBs) ( FIG. 1 ( b ) ), using a r 2 cutoff of 0.5. Using slightly different cutoff values, such as 0.3 or 0.7, for partitioning of the human genome into genetic haplotype blocks, resulted in only minor quantitative differences in the block size and number without affecting the global pattern (data not shown).

To characterize the global pattern and distribution of MHBs, starting with 51 sets of published Whole Genome Bisulfite Sequencing (WGBS) data from human primary tissues (Schultz et al. (2015); Heyn et al. (2012)), as well as the H1 human embryonic stem cells, in vitro derived progenitors (Xie, W. et al. (2013)), and human cancer cell lines (Blattler et al. (2014); Heyn et al. (2016)). An additional WGBS data set from 10 adult tissues from a single human donor was also generated and included. Across this set of 61 samples (>2000× combined genome coverage) a total of ˜55 billion methylation haplotype informative reads covering 58.2% of autosomal CpGs were identified. The uncovered CpG sites were either in regions with low mappability, or CpG sparse regions where there are too few sites within Illumina read pairs for deriving informative haplotypes. 147,888 MHBs with an average size of 95 bp and minimum 3 CpGs per block were identified, representing ˜0.5% of the human genome that tends to be tightly co-regulated on the epigenetic status at the level of single DNA molecules ( FIG. 11 ( a )- 11 ( b ) ). The majority of CpG sites within the same MHBs are near perfectly coupled (r 2˜ 1.0) regardless of the sample type. Methylation LD was found to extend further along the DNA in stem cells and progenitors, compared with normal adult tissue, both in the fraction of tightly coupled CpG pairs (94.8% versus 91.2%, P-value<2.6×10 −16 ), and the over-representation of partially coupled CpG pairs that are over 100 bp apart while the linkage was slightly decayed in primary cancer data sets (87.8%, mixture of CRC and LC), which was validated by another independent WGBS data set from kidney cancer (Chen et al. (2016)) ( FIG. 1 ( c ) ). Gene Ontology analysis show cancer loss of linkage regions was significantly associated with number of cancer related pathways and functions. This is consistent with previous observations on a smaller BSPP data set comprising 2,020 CpG islands (Shoemaker et al. (2010)) for culture cell lines and another previous report (Shao et al. (2014)). Interestingly, in tumor samples, a reduction of perfectly coupled CpG pairs was observed, which could be related to the pattern of discordant methylation recently reported in variable methylation regions (VMR) (Landau et al. (2014); Hansen et al. (2011)).

While WGBS data allowed MHBs across the entire genome to unbiasedly be identified, the 61 sets of data did not represent the full diversity of human cell/tissue types. To validate the presence of MHBs in a wider range of human tissues and culture cells, 101 published reduced representation bisulfite sequencing (RRBS) data sets from ENCODE cell lines and tissue samples were examined, as well as 637 sets of Infinium HumanMethylation450 BeadChip (HM450K) data including 11 normal human tissues from the TCGA project. The ENCODE RRBS data sets were generated with short (36 bp) Illumina sequencing reads, greatly limiting the length of methylation haplotypes that can be identified. Similarly, Illumina methylation arrays only report average CpG methylation of all DNA molecules in a sample, preventing a methylation linkage disequilibrium analysis. Therefore, the invention calculated the pairwise correlation coefficient of adjacent CpG methylation levels across different sample sets for block partitioning. Note that the presence of such correlated methylation blocks is a necessary but not sufficient condition for MHBs ( FIG. 12 ( a ) ). Nonetheless, the absence of correlated methylation blocks in these data sets would invalidate the pattern of MHBs. 23,517 and 2,212 correlated methylation blocks from ENCODE RRBS and TCGA HM450K array data, respectively, were identified, among which 8,920 and 1,258 have significant overlaps with WGBS-defined MHBs. Additionally, significantly higher correlation among the CpGs within the MHB regions compared to CpG loci outside MHBs in the HM450K and RRBS data sets was observed, further supporting the block-like organization of local CpG co-methylation across a wide variety of cells and tissues ( FIG. 12 ( b ) ). Taken together, the MHBs identified represent a distinct class of genomic features where local CpG methylation is established or removed in a highly coordinated manner at the level of single DNA molecules, presumably due to the processive activities of the related enzymes coupled with the local density of CpG dinucleotides.

Co-Localization of Methylation Haplotype Blocks with Known Regulatory Elements.

The MHBs established by 61 sets of WGBS data appear to represent a distinct type of genomic feature that partially overlaps with multiple well-documented genomic elements ( FIG. 1 ( d ) , FIG. 11 ). Among all the methylation blocks, 60,828 (41.1%) were located in intergenic regions and 87,060 (58.9%) regions in transcribed regions. These MHBs were significantly (p-value<10 −6 ) enriched in enhancers (enrichment factor=7.6), super-enhancers (enrichment factor=2.3), promoter regions (enrichment factor=14.5), CpG islands (enrichment factor=70.4), and imprinted genes (enrichment factor=54.6). In addition, modest depletion in LAD was observed (Guelen et al. (2008)) and LOCK regions (Wen et al. (2009)) (46% and 37% of the expected values), and modest enrichment in TAD (Dixon et al. (2012)). Importantly, a very strong (26-fold) enrichment in variable methylation regions (VMR) was observed (Hansen et al. (2011)) ( FIG. 1 ( e ) ), suggesting that increased epigenetic variability in a cell population or tissue can be coordinated locally among hundreds of thousands of genomic regions (Pujadas et al. (2012)). A subset of MHBs that do not overlap with CpG islands was examined, and a consistent enrichment pattern ( FIG. 1 ( e ) ) was observed, suggesting that local CpG density alone does not account for the enrichment.

Previous studies on mouse and human samples (Irizarry et al. (2009); Ziller et al. (2013)) demonstrated that dynamically methylated regions were associated with regulatory regions such as enhancer-like regions marked by H3K27ac and transcription factor binding sites. In human, 21.8% of autosomal CpGs were found to be differentially methylated across 30 human cell and tissue types (Hansen et al. (2011)). These CpGs were enriched at low to intermediate CpG density promoters. Using publicly available histone mapping data for human adult tissues, co-localization of methylation haplotype blocks with marks for active promoters (H3K4me3 with H3K27ac) was found, but not for active enhancers (Leung et al. (2015)) (no peak for H3K4me1). Meanwhile, enhancers were found to tend to overlap with CpG sparse MHBs, whereas the overlap with super-enhancers were independent of CpG density ( FIG. 11 ( c ) ). Therefore, MHBs likely capture the local coherent epigenetic signatures that are directly or indirectly coupled with transcriptional regulation.

Block-Level Analysis of Human Normal Tissues and Stem Cell Lines with Methylation Haplotype Load.

To enable quantitative analysis of the methylation patterns within individual MHBs across many samples, a single metric to define the methylated pattern of multiple CpG sites within each block is needed. Ideally this metric is not only a function of average methylation level for all the CpG sites in the block, but can also capture the pattern of co-methylation on single DNA molecules. For this purpose, methylation haplotype load (MHL) and unmethylated haplotype load (uMHL) were defined, where the first is a weighted mean of the fraction of fully methylated haplotypes and substrings at different lengths (i.e. all possible substrings) and the latter is a weighted mean of the fraction of fully unmethylated haplotypes and substrings at different lengths. Compared with other metrics used in the literature (methylation level, methylation entropy, epi-polymorphism, and haplotypes counts), MHL is capable of distinguishing blocks that have the same average methylation but various degrees of coordinated methylation ( FIG. 2 ). In addition, MHL and uMHL are bounded between 0 and 1, which allows for direct comparison of different regions across many data sets without normalization.

The invention addressed whether treating MHBs as individual genomic elements and performing quantitative analysis based on MHL and uMHL would provide an advantage over previous approaches using the weighted average methylation in genomic windows. To this end, the invention sought to identify tissue specific MHBs from a collection of human solid tissues WGBS datasets based on the MHL and uMHL. A group specific index for each MHB (see Methods) was computed using either MHL or uMHL for every tissue type and then each MHB was assigned to one or more tissues based on their respective GSI. If no other tissue type have a GSI that is at least 80% of the maximum GSI, then only the tissue type with the maximum GSI is reported. From the top 500 MHBs assigned to each tissue, the top tissue specific MHL or uMHL regions were selected using average MHL or uMHL values in blood cells of less than 0.05, and average non-blood cell values greater than 0.4. Using a maximum GSI of greater than 0.6 as a threshold resulted in the identification of 1,290 and 15,377 tissue specific MHL and uMHL regions respectively. To identify the most informative MHL or uMHL marker regions for cancer detection or tissue mapping the top 10% of the identified MHL or uMHL regions determined by GSI value are reported in Tables 1(a) and 1(b). The tissue specific regions using MHL and uMHL were visualized using heatmaps ( FIG. 3 ( a )- 3 ( b ) ). To demonstrate that MHL was a better metric than AMF, a set of MHBs overlapping with published tissue specific methylated regions (Lokk et al. (2014)) were identified. Using this set of MHBs, the performance between MHL and average methylation fraction in the MHL regions (AMF) were compared. Both MHL and AMF were able to identify tissue specific MHBs, but MHL has better signal (average value from within the tissue specific group) to noise (average value from other groups) than average methylation ( FIG. 4 ).

The human adult tissues used in this invention have various degrees of similarity amongst each other. The invention hypothesized that this is primarily defined by their developmental lineage, and that the related MHBs might reveal epigenetic insights related to germ layer speciation. All the data sets based on the three germ layers were grouped, and searched for MHBs that have differential MHL. In total 114 ectoderm-specific MHBs (99 hyper- and 15 hypo-methylated), 75 endoderm specific MHBs (58 hyper- and 17 hypo-methylated) and 31 mesoderm specific MHBs (9 hyper- and 22 hypo-methylated) were identified. The invention speculated that some of these MHBs might capture binding events of transcription factors (TF) specific to developmental germ-layers. Compared with ENCODE TFBS data (The ENCODE Project Consortium (2012)), distinctive patterns of TFs binding to layer specific MHBs were observed. For layer specific MHBs with hypo-methylation MHL, which tends to represent activation signals, 53 TF binding events in mesoderm specific MHBs, 71 in endoderm specific MHB, and 2 in ectoderm specific MHBs were identified. Gene ontology analysis showed TFs binding to mesoderm exhibit negative regulator activity, while TFs binding to endoderm exhibited positive regulator activity. For layer specific MHBs with hyper-methylation MHL, which tend to represent repressive signals, 38 TF binding events in mesoderm specific MHBs, 102 in endoderm specific MHB, and 145 in ectoderm specific MHBs were identified. Interestingly, ectoderm and endoderm shared few bounded TFs, while mesoderm tissues share multiple groups of TFs with ectoderm and endoderm. Two endoderm specific hyper-MHL regions were identified, which are related to ESRRA and NANOG. This is consistent with a previous finding that mouse ES cells differentiated spontaneously into visceral/parietal endoderm upon NANOG knock-out (Mitsui et al. (2003)). Gene ontology analysis showed that hypo-MHL regions shared by mesoderm and endoderm might have regulatory functions in the fate commitment towards multiple tissues, whereas ectoderm specific hyper-MHL regions might induce the ectoderm development by suppressing the path towards the immune lineage. These observations are indicative of two distinctive “push” and “pull” mechanisms in the transition of cell states that have been harnessed for the induction of pluripotency by over-expressing lineage specifiers (Shu et al. (2013)).

Methylation-Haplotype Based Analysis of Circulating Cell-Free DNA in Cancer Patients and Healthy Donors.

A unique aspect of methylation haplotype analysis is that the pattern of co-methylation, especially within MHBs, is robust for capturing low-frequency alleles among a heterogeneous population of molecules or cells, in the presence of biological noise or technical variability (i.e. incomplete bisulfite conversion or sequencing errors). To explore the clinical potential, the invention focused on the methylation haplotype analysis of cell-free DNA from healthy donors and cancer patients, of which various low fractions of DNA molecules were released from tumor cells and potentially carry epigenetic signatures different from blood. 4-122 ng (average 20 ng) of cell-free DNA from an average of 866 μL human plasma from 75 normal individuals and 59 cancer patients were isolated, except for four with unusually high yield due to cell lysis. Due to the limited DNA availability, scRRBS was performed (Guo et al. (2013)) on 1 to 10 ng of cfDNA from 134 plasma samples and obtained an average of 13 million paired-end 150 bp reads per sample. On average, 57.7% WGBS-defined MHBs were covered in the RRBS data set from the clinical samples.

The invention sought to detect the presence of tumor specific signatures in the plasma samples, using the methylation haplotypes identified in the reference tumor tissues and in normal samples used as the negative controls. For five lung cancer plasma samples and five colorectal cancer plasma samples, matched primary tumor tissues were also obtained, and generated RRBS data (30 million reads per sample) from 100 ng of tumor genomic DNA. The invention focused on MHBs with low MHL (i.e. genomic regions that have low or no methylation) in the blood, and determined whether cancer-associated highly methylated haplotypes (caHMH) can be detected. The invention provides in some embodiments that such haplotypes are present only in the tumor tissues and the matched plasma from the same patient, but not in whole blood or any other non-cancer samples. These highly confident tumor signatures in circulating DNA were considered. caHMH in all cancer patient plasma samples (Average=36, interquartile range (IQR)=17) were detected. These caHMHs were associated with 320 genes, some of which are known to be aberrantly methylated in human cancers such as WDR37, VAX1, SMPD1 (Table 2) The 49 additional cancer plasma samples with no matched tumor samples were analyzed, using 65 normal plasma samples as background and negative control. On average 60 (IQR=31) caHMH were identified for each cancer plasma sample. Interestingly, a significant fraction (35%) of caHMH identified with matched tumor-plasma pairs were also detected in the expanded set of cancer patient plasma samples. A majority of caHMHs were found to be individual specific while few caHMHs were present in at least 53% (16/30) and 62% (18/29) cancer plasma samples for CRC and LC ( FIG. 13 ). Improving the sampling depth, by either using more input cfDNA or reducing sample loss during analysis, will likely increase the number of caHMHs commonly observed in multiple patients.

The tumor load in cancer plasma samples was quantified, using non-negative decomposition with quadratic programming, on the RRBS data from primary cancer biopsies (LC and CRC) and the WGBS data from 10 normal tissues. (Table 3) The invention estimated that a predominant fraction, 72.0% (IQR=40%) in cancer and normal plasma samples were contributed by white blood cells, which is consistent with the levels recently reported based on shallow whole genome bisulfite sequencing (69.4%) (Sun et al. (2015)). Primary tumor and normal tissue-of-origin contributed at the similar level of 2.3% (IQR=3.7%) and 3.0% (IQR=4.4%). In contrast, the similar analysis applied to normal plasma only found residual tumor contributions (0.17%, IQR=2.9% for CRC and 1.0%, IQR=3.1% LC) to normal plasma, which were significantly lower (P=3.4×10 −5 and 5.2×10 −10 for CRC and LC, respectively) than for cancer plasma. 76.7% plasma samples from CRC patients and 89.6% from LC patients were found to have a detectible contribution from tumor tissues while only 13% and 26% normal plasmas have certain (low) tumor contribution. Therefore, circulating cell-free DNA contains a relatively stable fraction of molecules released from various normal tissues, whereas tumor cells released DNA molecules that can be more abundant in tumor cells than in normal tissues (Tables 3(a)-3(d)). The fractions of white blood cells observed are lower than what was reported previously (Sun et al. (2015)), most likely due to the inclusion of 10 normal tissue types in the deconvolution analysis.

The invention sought to use the information from normal human tissues, primary tumor biopsies, and cancer cell lines to improve the detection of cfDNA. The invention selected a subset of MHBs that show high MHL in primary cancer biopsies and low MHL in normal control plasma. A subset of MHBs that have high MHL in cancer tissues and low MHLs in normal plasma was identified for each cancer type (Table 4(a)-(b)). Cancer plasma showed significantly higher MHL in these regions than independent normal plasma (P=1.2×10 −13 and 2.2×10 −16 for CRC and LC, respectively) ( FIG. 5 and FIG. 6 ). By computationally mixing the sequencing reads from cancer tissues and whole blood samples (WBC), synthetic admixtures at various levels of tumor fraction were created to calibrate the relationship between tumor load and the MHL values in these regions ( FIG. 7 ). A fitted linear model on the standard curve for colon cancer markers had an adjusted R-squared value of 0.9621 and for lung cancer markers the adjusted R-squared value was 0.9573. Note that these MHBs were selected without using any information from the cancer plasma samples, and hence they should be applicable to other cancer plasma samples. Tumor load estimation was performed on test sets of cancer and normal plasma samples. Cancer plasma was found to have significantly higher tumor load than normal plasma (Two Sample t-test with unequal variance, P=0.005074 for colon cancer plasma versus normal and P=0.0002783 for lung cancer plasma versus normal plasma) ( FIG. 8 ).

Recent studies (Sun et al. (2015); Lehmann-Werman et al. (2016); Snyder et al. (2016)) have demonstrated that epigenetic information imbedded in cfDNA has the potential for predicting tumors tissue-of-origin. Here the invention addresses whether a MHL-based framework and a set of targets derived from whole genome data would allow for the prediction of tissue-of-origin with quantifiable sensitivity and specificity, which is crucial for future clinical applications. Training WGBS tissue data from the following tissues: colon, lung, neural, heart, liver, lung, pancreas, and stomach were compiled and a set of 15,000 MHBs selected using their Group Specific Index (GSI) was created. Note that the MHBs with a high GSI score tended to be methylated in fewer tissue types. 20 colon cancer plasma, 20 lung cancer plasma, and 30 normal plasma samples were randomly sampled to create training data sets from the RRBS plasma data sets. The remaining (10 colon, 10 lung, and 39 normal) were held out as test data sets. An ensemble MARS (Multivariate adaptive regression splines) model was generated using the training data sets (Friedman, (1991)). The invention also employed K-means clustering of the MHBs using WGBS data to generate 50 clusters and each MARS model in the ensemble only saw a feature set consisting of 3 features sampled from a cluster (3 features×150 clusters=450 features for each model). To minimize overfitting of training data sets each model selected at most 2 features to use during classification. The first classifier was an ensemble model with 100 MARS models with 154 MHB features (Table 5(a)). In a manner analogous to a random forest classifier, the scores from the resulting ensemble are averaged to compute a prediction score. By utilizing these prediction scores as a binary classifier, the invention was able to obtain a colon vs non-colon AUC of 0.725, a lung vs non-lung AUC of 0.751, and a normal vs abnormal AUC of 0.856 ( FIG. 9 ( a )- 9 ( c ) ). Note that these AUC values were obtained using an initial marker set identified from an independent training set of normal WGBS tissue samples only; no cancer tissue samples were used in the initial feature selection, suggesting that plasma samples can be segregated using only tissue specific markers.

To obtain a prediction of tissue-of-origin, the invention focused on colon and lung cancer plasma. First, the top MHBs features that were assigned to either colon, lung, or both tissues ranked by GSI were selected. An ensemble of 500 MARS models identifying 295 unique tissue specific features was made (Table 5(b)). PLSDA (partial least squares discriminant analysis) implemented in the caret R package was used to perform the classification. The classification using this model was assessed using independent test data sets and the accuracy for classifying colon cancer plasma to colon tissue was found to be 70% (7/10) and the accuracy for classifying lung cancer plasma to lung tissue was 78% (7/9) ( FIG. 10 ).

In this invention a well-established concept in population genetics, linkage disequilibrium (LD), was extended to the analysis of co-methylated CpG patterns. While the mathematical representations are identical, there are two key differences. First, traditional linkage disequilibrium (LD) was defined for human individuals in a population, whereas in this invention the analysis was performed on the diploid genome of individual cells in a heterogeneous cell population. Second, linkage disequilibrium in human populations depend on the mutation rate, frequency of meiotic recombination, effective population size and demographic history. The LD level typically decays over the range of hundreds of kilobases to megabases. In contrast, CpG co-methylation depends on DNA methytransferases and demethylases, which tend to have lower processivity, and, in the case of hemi-methyltransferases, much lower fidelity compared with DNA polymerases (Williams et al. (2011)). Therefore, methylation LD decays over much shorter distances ranging from tens to hundreds of bases, with the exception of imprinting regions. Even if longer-read sequencing methods were used no radical change of the block-like pattern presented in this work is expected, which is supported by a recent study (Saito et al. (2015)). Nonetheless, these short and punctuated blocks capture discrete entities of epigenetic regulation in individual cells widespread in the human genome. Such a phenomenon can be harnessed to improve the robustness and sensitivity of DNA methylation analysis, such as the deconvolution of data from heterogeneous samples including circulating cell-free DNA.

Epigenetic abnormalities tend to be more widespread across the genome (compared with somatic mutations), and hence enabling the integration of the sparse coverage across many loci to achieve very accurate prediction by direct counting of methylated haplotypes within the appropriate tissue-specific features.

Methods

Normal and Cancer Samples.

Ten human primary tissues were purchased from BioChain. Cancer tissue and plasma samples were collected from UCSD Moores Cancer Center and normal plasma samples were obtained from UCSD Shiley Eye center under IRB protocols approved by UCSD Human Research Protections Program (HRPP).

Generation of DNA Libraries for Sequencing.

Extracted genomic DNA were prepared for bisulfite sequencing using published protocols. For whole genome bisulfite (WGBS) and reduced representation bisulfite sequencing (RRBS), the DNA fragments were adapted to barcoded methylated adaptors (Illumina). For WGBS, the adapted DNA were converted using the EZ DNA Methylation Lightning kit (Zymo Research) and amplified for 10 cycles using iQ SYBR Green Supermix (BioRad). For RRBS, the adapted DNA were converted using the MethylCode™ Bisulfite Conversion kit (Thermo Fisher Scientific) and amplified using the PfuTurboCx polymerase (Agilent) for 12-14 cycles. Libraries were pooled and size selected using 6% TBE polyacrylamide gels. Libraries were sequenced using the Illumina HiSeq platform for paired-end 100-111 cycles, the Illumina MiSeq platform for paired-end 75 cycles, and the GAIIx (WGBS only) for single-end 36 cycles.

Methylation Haplotype Blocks (MHB).

The human genome was separated into non-overlapping “sequenceable and mappable” segments using a set of generated WGBS data from 10 tissues from a 25-year adult male individual. Mapped reads from WGBS data sets were converted into methylation haplotypes in each segment. Methylation linkage disequilibrium was calculated on the combined methylation haplotypes. Each segment was partitioned into methylation haplotype blocks (MHBs). MHBs were defined as the genomic region in which the r 2 value of two adjacent CpG sites is no less than 0.5. MHB regions inferred by WGBS data sets were also validated by bulk data of methylation level. Takai and Jones's sliding-window algorithm (Takai et al. (2002)) was applied for methylation high linkage regions in the HM450K (TCGA) and the RRBS (ENCODE) data set. Finally, simulation analysis to investigate the relationship between LD and correlation of average 5 mC of two CpG loci were conducted based on random sampling of different methylation haplotypes with 1000 individuals, and each individual sampling 10 methylation haplotypes.

Methylation Haplotype Load (MHL).

A methylated haplotype load (MHL) for each candidate region was defined, which is the normalized fraction of methylated haplotypes (MH) at different lengths:

MHL = ∑ i = 1 l ⁢ w i × P ⁡ ( MH i ) ∑ i = 1 l ⁢ w i w i = i

Unmethylated Haplotype Load (uMHL).

An unmethylated haplotype load (uMHL) for each candidate region was defined, which is the normalized fraction of unmethylated haplotypes (UMH) at different lengths:

Where 1 is the length of haplotypes, P(MH i ) or P(UMH i ) is the fraction of fully successive methylated or unmethylated haplotype with i loci respectively. For a haplotype of length L, all the sub-strings with length from 1 to L were considered in this calculation. w i is the weight for i-locus haplotype. w i =i or w i =i 2 was typically used to favor the contribution of longer haplotypes. In the present invention, w i =i was applied. Quantile normalization, standardization (scale) as well as the batch effect elimination (Johnson et al. (2007)) were applied.

Developmental Germ Layers and Tissue Specific MHB Regions.

In order to investigate the layer and tissue specific MHB regions, group specific index (see below) was applied. An empirical threshold of 0.6 was selected to filter out layer and tissue specific MHB regions. Layer specific MHB regions were selected to show the ability to distinguish between the different development layers. Tissue specific MHB regions were further used for tissue mapping and cancer diagnosis.

n indicates the number of the groups. MHL(j) denotes the average of MHL of j th group. MHL max denotes the average of MHL of highest methylated group.

Deconvolution Analysis.

The deconvolution references were constructed from normal human solid tissues, WBC, colorectal cancer tissues (CCT), and lung cancer tissues (LCT). Tissue specific MHB regions for normal human tissues were selected for brain, colon, esophagus, heart, intestine, kidney, liver, lung, and stomach using candidate features for deconvolution based on non-negative decomposition with quadratic programming (Sun et al. (2015); Houseman et al. (2012); Gong et al. (2013)). Raw MHL signals were logit transformed before deconvolution analysis on plasma samples. Samples with less than 30% whole blood content from deconvolution analysis were considered to have failed due to poor library complexity.

Tumor Load Estimation Using Tumor Specific MHBs.

Pruning and K-nearest neighbors (KNN) imputation was performed on the MHL matrix with only the RRBS plasma and RRBS tumor tissue samples removing samples with low coverage and imputed missing values. Thus 30 colon cancer plasma, 29 lung cancer plasma, 69 normal plasma, 4 colon cancer tissues, and 5 lung cancer tissues remained. The 69 healthy normal plasma samples were split into “training” and “test” sets; with 46 samples set aside for feature selection and training while the remaining 23 samples were used as a completely independent data set to test the quantitation. Tumor specific methylation haplotype blocks were identified by 2-tailed t-test with a False Discovery Rate (FDR) of 0.001 and a minimum difference cutoff of 0.3. Two sets of markers were identified separately for colon cancer and lung cancer. To calibrate the relationship between tumor load and MHL values, 20 sets of simulations in which mixed sequencing reads from cancer tissue samples and normal plasma samples at a 1:5, 1:10, 1:20, 1:100, and 0:1 ratio were generated (totaling 100 simulated data sets for each cancer tissue). The average MHL value for these regions was computed for each region and a linear regression model was generated using these values and the known cancer tissue proportions for each cancer tissue. The model was applied to the average MHL value in these regions for test sets which included 30 colon cancer plasma, 29 lung cancer plasma, and 23 normal plasma samples. The estimated tumor loads for normal versus colon cancer and for normal versus lung cancer were compared using the colon cancer and lung cancer markers respectively.

Cancer Plasma Classification.

Training WGBS data were collected to generate an MHL matrix and the top 15,000 MHBs by GSI were selected. These MHBs were clustered using K-means into 50 marker groups with the WGBS data matrix. Each cluster should have similar methylation patterns across each tissue type. In order to avoid overfitting, 20 colon cancer plasma, 20 lung cancer plasma, and 30 normal plasma samples were randomly selected to create a training data set from the RRBS plasma samples. The remaining samples (10 colon, 10 lung, and 39 normal) were held out as a test data set. For building the classifiers, an ensemble MARS (Multivariate Adaptive Regression Splines) model (Friedman, (1991)) implemented in the Earth R package (https://cran.r-project.org/web/packages/earth/index.html) was trained on the training data set using features from the top 15,000 MHBs ranked by the Group Specific Index (GSI). The ensemble MARS model was used to perform classification on the test plasma data set using the binary classifier. For cancer type classification, the features selected using MARS were used in a PLSDA (partial least squares discriminant analysis) to distinguish colon versus lung.

Data Availability

WGBS and RRBS data are available at the Gene Expression Omnibus (GEO) under accession GSE79279.

Tables:

Table 1. Top 10% tissue specific MHL and uMHL markers identified by GSI.

Table 2. Complete list of high methylated haplotype shared between matched primary tumor tissues and plasma for colon cancer (CRC) and lung cancer (LC) patients.

Table 3. Deconvolution of colon cancer (CRC), lung cancer (LC) and normal plasma into ten tissues using non-negative decomposition with quadratic programming Samples with white blood cells composition greater than 30% were used to report averages.

Table 4. Differentially methylated MHB regions between cancer tissues and normal plasma.

Table 5. The sets of cancer specific and tissue specific markers derived from MARS based features selection on training data sets.

The disclosures of all references listed herein, including the following references, are incorporated by reference herein in their entireties:

1. EP Patent No. 2,258,871

2. EP Patent No. 2,643,479

3. U.S. Pat. No. 8,026,067

4. U.S. Pat. No. 8,110,361

5. U.S. Pat. No. 8,288,100

6. U.S. Pat. No. 8,293,470

7. U.S. Pat. No. 8,586,313

8. U.S. Pat. No. 9,096,906

9. WO 2015/116837

10. Accession GSE79279, Gene Expression Omnibus (GEO)

11. The ENCODE Project Consortium. An integrated encyclopedia of DNA elements in the human genome. Nature 489:57-74 (2012).

12. Bernstein et al. The NIH Roadmap Epigenomics Mapping Consortium. Nat Biotechnol 28:1045-1048 (2010).

13. Blattler et al. Global loss of DNA methylation uncovers intronic enhancers in genes showing expression changes. Genome Biol 15:469 (2014).

14. Chen et al. Loss of 5-hydroxymethylcytosine is linked to gene body hypermethylation in kidney cancer. Cell Res 26:103-118 (2016).

15. Diep et al. Library-free methylation sequencing with bisulfite padlock probes. Nat Methods. 9(3):270-272 (2012).

16. Dixon et al. Topological domains in mammalian genomes identified by analysis of chromatin interactions. Nature 485:376-380 (2012).

17. Friedman Multivariate Adaptive Regression Splines. Ann. Of Stat. 19(1): 1-67 (1991).

18. Gong et al. DeconRNASeq: a statistical framework for deconvolution of heterogeneous tissue samples based on mRNA-Seq data. Bioinformatics 29:1083-1085 (2013).

19. Guelen et al. Domain organization of human chromosomes revealed by mapping of nuclear lamina interactions. Nature 453:948-951 (2008).

20. Guo et al. Single-cell methylome landscapes of mouse embryonic stem cells and early embryos analyzed using reduced representation bisulfite sequencing. Genome Res 23:2126-2135 (2013).

21. Hansen et al. Increased methylation variation in epigenetic domains across cancer types. Nat Genet 43:768-775 (2011).

22. Heyn et al. Distinct DNA methylomes of newborns and centenarians. Proc Natl Acad Sci USA 109:10522-10527 (2012).

23. Heyn et al. Epigenomic analysis detects aberrant super-enhancer DNA methylation in human cancer. Genome Biol 17:11 (2016).

24. Houseman et al. DNA methylation arrays as surrogate measures of cell mixture distribution. BMC Bioinformatics 13:86 (2012).

25. Houseman et al. Reference-free deconvolution of DNA methylation data and mediation by cell composition effects. BMC Bioinformatics 17:259 (2016).

26. Irizarry et al. The human colon cancer methylome shows similar hypo- and hypermethylation at conserved tissue-specific CpG island shores. Nat Genet 41:178-186 (2009).

27. Johansson et al. Targeted resequencing of candidate genes using selector probes. Nucleic Acids Res. 39(2):e8 (2011).

28. Johnson et al. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics 8:118-127 (2007).

29. Jones, P. A. & Martienssen, R. A blueprint for a Human Epigenome Project: the AACR Human Epigenome Workshop. Cancer Res 65:11241-11246 (2005).

30. Komori et al. Application of microdroplet PCR for large-scale targeted bisulfite sequencing. Genome Res. 21(10):1738-1745 (2011).

31. Landan et al. Epigenetic polymorphism and the stochastic formation of differentially methylated regions in normal and cancerous tissues. Nat Genet. 44(11):1207-1214 (2012).

32. Landau et al. Locally disordered methylation forms the basis of intratumor methylome variation in chronic lymphocytic leukemia. Cancer Cell 26:813-825 (2014).

33. Lehmann-Werman et al. Identification of tissue-specific cell death using methylation patterns of circulating DNA. Proc Natl Acad Sci USA 113, E1826-1834 (2016).

34. Leung et al. Integrative analysis of haplotype-resolved epigenomes across human tissues. Nature 518:350-354 (2015).

35. Lokk et al. DNA methylome profiling of human tissues identifies global and tissue-specific methylation patterns. Genome Biology 15(4):r54 (2014).

36. Meissner et al. Reduced representation bisulfite sequencing for comparative high-resolution DNA methylation analysis. Nucleic Acids Res. 33(18):5868-5877 (2005).

37. Mitsui et al. The homeoprotein Nanog is required for maintenance of pluripotency in mouse epiblast and ES cells. Cell 113:631-642 (2003).

38. Papageorgiou et al. Fetal-specific DNA methylation ratio permits noninvasive prenatal diagnosis of trisomy 21. Nat Med. 17(4):510-513 (2011).

39. Pujadas et al. Regulated noise in the epigenetic landscape of development and disease. Cell 148:1123-1131 (2012).

40. Saito et al. Linkage disequilibrium analysis of allelic heterogeneity in DNA methylation. Epigenetics 10:1093-1098 (2015).

41. Schultz et al. Human body epigenome maps reveal noncanonical DNA methylation variation. Nature 523:212-216 (2015).

42. Shao et al. Deciphering the heterogeneity in DNA methylation patterns during stem cell differentiation and reprogramming BMC Genomics 15:978 (2014).

43. Shoemaker et al. Allele-specific methylation is prevalent and is contributed by CpG-SNPs in the human genome. Genome Res 20:883-889 (2010).

44. Shu et al. Induction of pluripotency in mouse somatic cells with lineage specifiers. Cell 153:963-975 (2013).

45. Slatkin Linkage disequilibrium—understanding the evolutionary past and mapping the medical future. Nat Rev Genet 9:477-485 (2008).

46. Snyder et al. Cell-free DNA Comprises an In Vivo Nucleosome Footprint that Informs Its Tissues-Of-Origin. Cell 164:57-68 (2016).

47. Sun et al. Plasma DNA tissue mapping by genome-wide methylation sequencing for noninvasive prenatal, cancer, and transplantation assessments. Proc Natl Acad Sci USA 112, E5503-5512 (2015).

48. Takai et al. Comprehensive analysis of CpG islands in human chromosomes 21 and 22. Proc Natl Acad Sci USA 99:3740-3745 (2002).

49. Wen, B., Wu, H., Shinkai, Y., Irizarry, R. A. & Feinberg, A. P. Large histone H3 lysine 9 dimethylated chromatin blocks distinguish differentiated from embryonic stem cells. Nat Genet 41:246-250 (2009).

50. Wigler et al. The somatic replication of DNA methylation. Cell 24:33-40 (1981).

51. Williams et al. TET1 and hydroxymethylcytosine in transcription and DNA methylation fidelity. Nature 473:343-348 (2011).

52. Xie, H. et al. Genome-wide quantitative assessment of variation in DNA methylation patterns. Nucleic Acids Res. 39(10):4099-4108 (2011).

53. Xie, W. et al. Epigenomic analysis of multilineage differentiation of human embryonic stem cells. Cell 153:1134-1148 (2013).

54. Ziller et al. Charting a dynamic DNA methylation landscape of the human genome. Nature 500:477-481 (2013).