JP2024515565A - Cell-free DNA sequencing data analysis methods to investigate nucleosome protection and chromatin accessibility - Google Patents

Abstract

In one aspect, the present disclosure provides a computer-implemented method for enhancing sequence read data from a cell-free DNA sample to predict cell type. The method includes receiving sequence read data including a plurality of fragment reads, each fragment read having a fragment length and a GC content indicating the percentage of bases that are G or C in the fragment read. A computing system determines a GC bias value for each fragment read based on the fragment length and GC content of the fragment read. The sequence read data and the GC bias value are used to generate a GC bias adjusted genome coverage distribution. Based on the genome coverage distribution, predict the cell type. The method can be utilized to assess cell subtypes and phenotypes based on cell-free DNA present in a biological sample, for example, for cancer diagnosis, monitoring, and precision therapy.

Description

CROSS-REFERENCE TO RELATED APPLICATIONS This application claims the benefit of U.S. Patent Application No. 63/172,590, filed April 8, 2021, and U.S. Patent Application No. 63/276,378, filed November 5, 2021, the entire disclosures of which are incorporated herein by reference.

GOVERNMENT LICENSE RIGHTS STATEMENT This invention was made with government support under grants CA228944, CA264383, CA237746, CA097186, CA234715, CA076930, and HL007093 awarded by the National Institutes of Health, and W81XWH-21-1-0513, W81XWH-18-1-0406, and W81XWH-17-1-0380 awarded by the United States Army Medical Research and Development Command. The Government has certain rights in this invention.

Background Metastatic cancer is a terminal cancer and often leads to cancer-related death. Treatment options at the time of diagnosis of new or recurrent metastatic cancer are often based on clinical diagnosis from the primary tumor. In addition, molecular changes in tumors, such as genetic alterations or phenotypic changes, may appear during metastatic progression or the development of resistance to treatment. For example, in breast cancer, conversion of hormone receptors is frequently observed during the development of resistance to targeted treatment. Therefore, it is important to classify tumor subtypes and identify patterns of transcriptional regulation that drive tumor phenotypic changes during treatment. This type of research has crucial implications for studying mechanisms of resistance to treatment and informing clinical treatment decisions to provide patients with life-prolonging treatments and care.

Current approaches to subtyping solid tumors involve taking tissue biopsies and applying imaging techniques such as immunohistochemistry to assess cellular phenotype. However, obtaining surgical biopsies for disease monitoring from patients with metastatic cancer is often difficult, especially for advanced cancers, due to the location and/or number of metastatic sites. Furthermore, investigating molecular changes in tumors is particularly challenging, as repeat biopsies are cumbersome and are not considered standard of care. Thus, accurate subtyping is crucial to address subtype switching during tumor recurrence and resistance to treatment, but characterizing tumor evolution and subtype plasticity during treatment is difficult, thereby illustrating a major limitation of current treatment strategies and precision medicine for patients with metastatic cancer.

To illustrate, breast cancer is one of the most common causes of cancer, accounting for 23% of cancer diagnoses and 14% of cancer-related deaths among women worldwide. Targeted therapy is guided by tumor subtype, including expression of three hormone receptors: ER, PR, and HER2. In approximately 15% of cases, breast cancer tumors undergo hormone subtype switching during tumor recurrence or as a mechanism of resistance to endocrine therapy. However, clinical determination of tumor subtype remains limited to tissue biopsies, which are not routinely taken in advanced disease or repeatedly obtained during the course of treatment.

Similarly, prostate cancer is the second leading cause of cancer mortality among men, with an estimated 33,000 deaths in the United States in 2020. Castration-resistant prostate cancer (CRPC) is a stage in which resistance to androgen deprivation therapy of the disease occurs and progresses to metastatic CRPC (mCRPC), a stage in which there is no curative treatment and which is invariably fatal. mCRPC is understood to include multiple distinct subtype lineages and molecular subtypes that are generally classified by specific genomic or epigenetic alterations. Prostate cancer can be categorized by phenotypic features, including a range of transdifferentiated pathologies, including neuroendocrine (NE) cancer, low androgen-regulated disease states (ARlowPC), and double-negative prostate cancer (DNPC, AR negative NE negative). These phenotypic subtypes are particularly important because they are increasingly frequent as they emerge in the setting of strong androgen deprivation therapy resistance and are indicative of distinctive prognoses for patients with advanced prostate cancer. In the new era of precision medicine in prostate cancer, several newer treatment options have emerged that are guided by specific molecular alterations in the tumor genome (currently PARP inhibitors, immune checkpoint blockade, and specific DNA repair defects). Accurate molecular classification defined by the relevant molecular landscape will be crucial, as it will create the scope of future "precision medicine" for prostate cancer.

Therefore, there remains a need for easy and accurate methods for identifying cell or tissue types, such as cancer phenotypes, that can be routinely performed in the clinic. The present disclosure addresses these and related needs.

SUMMARY This Summary is provided to introduce a selection of concepts in a simplified form that are further described below in the Detailed Description. This Summary is not intended to identify key features of the claimed subject matter, nor is it intended to be used as an aid in determining the scope of the claimed subject matter.

In one aspect, the present disclosure provides a computer-implemented method for augmenting sequence read data from a cell-free DNA sample to predict cell type, the method comprising:
receiving, by a computing system, sequence read data, the sequence read data comprising a plurality of fragment reads, each fragment read having a fragment length and a GC content indicating a percentage of bases in the fragment read that are G or C;
determining, by the computing system, a GC bias value for each fragment read based on the fragment length and the GC content of the fragment read;
generating, by the computing system, a genomic coverage distribution adjusted for GC bias using the sequence read data and the GC bias value; and predicting, by the computing system, the cell type based on the genomic coverage distribution.

In one embodiment, predicting the cell type based on the genome coverage distribution includes predicting a cell phenotype. In one embodiment, predicting the cell phenotype includes predicting a tissue type, a cancer type, or a cancer subtype. In one embodiment, predicting the cell phenotype includes predicting expression of one or more genes of interest. In one embodiment, determining the GC bias value based on the fragment length and the GC content of the fragment reads includes counting the number of observed reads of each combination of fragment length and GC content to determine a GC count for the sequence read data, dividing the GC count by the corresponding GC frequency in a GC frequency matrix to determine the GC bias for each fragment length, normalizing the average GC bias for each fragment length to determine an approximate GC bias value, and smoothing the approximate GC bias value to determine the GC bias value. In one embodiment, the GC occurrence matrix stores the occurrence of each GC content for each fragment length of a plurality of fragment lengths in a mappable region of a reference genome. In one embodiment, the plurality of fragment lengths includes fragment lengths ranging from a short length threshold to a long length threshold. In one embodiment, the short length threshold is in the range of 10-20 base pairs and the long length threshold is in the range of 450-550 base pairs. In one embodiment, the short length threshold is 15 base pairs and the long length threshold is 500 base pairs. In one embodiment, the method further comprises determining genomic regions of interest with respect to a cell type and filtering the genomic regions of interest to identify sites that provide information about the cell type. In one embodiment, determining the genomic regions of interest comprises determining an average mappability within a fixed size window around each genomic region of interest and discarding genomic regions of interest with an average mappability below a predetermined threshold. In one embodiment, filtering the genomic region of interest to identify sites that provide information about the cell type includes determining sites that have differential signals between the first cell type and the second cell type.

In one embodiment, the step of generating the genome coverage distribution includes determining fragment midpoints within a window around each site that provides information about a cell type, assigning a weight to each fragment read based on the inverse of the GC bias value for each fragment read, using the weighted fragment reads to determine a GC-corrected midpoint coverage profile, excluding positions that overlap with an excluded region, determining a mean profile based on determining an average of the GC-corrected midpoint coverage profiles for all sites, smoothing the mean profile to generate a smoothed mean profile, and normalizing the smoothed mean profile by dividing by the average of the surrounding coverage to determine a normalized mean profile. In one embodiment, the excluded regions include one or more regions that are within the encoding integrated GRCh38 excluded list, centromeres, gaps in the human genome assembly, modified patches, alternative haplotypes, regions of zero mappability, or have coverage at least 10 standard deviations above the mean. In one embodiment, predicting the cell type based on the genome coverage distribution comprises generating one or more features based on the genome coverage distribution, providing the one or more features as input to a classifier model, and determining the cell type based on an output of the classifier model. In one embodiment, the one or more features comprise an average of coverage within a first predefined window around each site providing information about a cell type, an average of coverage within a second predefined window around each site providing information about a cell type, the second predefined window having a different size than the first predefined window, and an amplitude of the genome coverage distribution around each site providing information about a cell type. In one embodiment, the first predefined window is larger than the second predefined window. In one embodiment, the width of the first predefined window is in the range of 1800-2200 base pairs and the width of the second predefined window is in the range of 40-80 base pairs. In one embodiment, the width of the first predefined window is 2000 base pairs and the width of the second predefined window is 60 base pairs. In one embodiment, the amplitude of the genome coverage distribution around each site that provides information about cell type is determined by trimming the genome coverage distribution to a window containing 10 peaks, performing a fast Fourier transform on the window of the genome coverage distribution, and determining the magnitude of the 10th frequency. In one embodiment, the classifier model comprises a logistic regression model, an artificial neural network, a decision tree, a support vector machine, or a Bayesian network.

In another aspect, the disclosure provides a method for determining a chromatin accessibility profile for a cell of interest from a sample comprising cell-free DNA derived from said cell of interest, the method comprising:
obtaining sequence read data from the cell-free DNA;
receiving, by a computing system, sequence read data, the sequence read data comprising a plurality of fragment reads, each fragment read having a fragment length and a GC content indicating a percentage of bases in the fragment read that are G or C;
determining, by the computing system, a GC bias value for each fragment read based on the fragment length and the GC content of the fragment read;
generating, by the computing system, a GC bias adjusted genome coverage distribution using the sequence read data and the GC bias value; and determining a chromatin accessibility profile from the genome coverage distribution.

In one embodiment, the method further comprises determining the cell phenotype of interest based on the chromatin occupancy profile. In one embodiment, determining the cell phenotype comprises determining a tissue type, a cancer type, a cancer subtype, an aggressive phenotype of a malignant tumor, and/or a drug responsiveness phenotype. In one embodiment, the method further comprises performing one or more steps of a computer-implemented method described herein.

In another aspect, the disclosure provides a method for determining the cell type of a cell of interest from a sample comprising cell-free DNA derived from said cell of interest, the method comprising:
obtaining sequence read data generated from said sample comprising cell-free DNA;
performing a computer-implemented method as described herein; and determining the cell type of the cell of interest based on a prediction provided by the computing system.

In one embodiment, the step of determining the cell type comprises determining a cell phenotype. In one embodiment, determining the cell phenotype comprises determining a tissue type, a cancer type, a cancer subtype, a high-grade phenotype of a malignant tumor, and/or a drug-responsive phenotype. In one embodiment, determining the cell phenotype comprises determining the expression of one or more genes of interest.

In another aspect, the disclosure provides a method of detecting the presence of cancer cells in a subject, comprising:
obtaining sequence read data generated from a sample comprising cell-free DNA obtained from the subject;
performing a computer-implemented method as described herein; and determining the presence of cancer cells in the subject based on the prediction provided by the computing system.

In one embodiment, the method is performed multiple times over time, with each performance further characterizing the cancer cell(s) detected in the subject and determining a cancer subtype or phenotype of the detected cancer cell(s) based on the prediction provided by the computing system. In one embodiment, the method is performed multiple times over time, and further includes detecting a change in the phenotype of the detected cancer cell(s) over time. In one embodiment, the subject receives a cancer treatment between performances of the method, and the method further includes determining the responsiveness of the cancer cell(s) to the treatment.

In another aspect, the present disclosure provides a method for determining a cancer subtype of a target cancer cell from a sample comprising cell-free DNA derived from said target cancer cell, said method comprising:
obtaining sequence read data generated from said sample comprising cell-free DNA;
22. The method according to claim 5, further comprising: performing a computer-implemented method according to any one of claims 5 to 21; and determining the cell type of a cell of origin based on a predicted cancer subtype provided by the computing system.

In one embodiment, the sample is obtained from a subject with cancer. In one embodiment, the cancer is characterized as metastatic breast cancer. In one embodiment, determining the cancer subtype comprises determining whether the cancer is ER+ or ER-. In one embodiment, determining the cancer subtype comprises determining whether the cancer is PR+ or PR-. In one embodiment, determining the cancer subtype comprises determining whether the cancer is HER2+ or HER2-. In one embodiment, determining the cancer subtype comprises:
whether the cancer is ER+ or ER-;
determining two or all of: whether the cancer is PR+ or PR-; and whether the cancer is HER2+ or HER2-.

In one embodiment, the cancer is characterized as metastatic prostate cancer. In one embodiment, determining the cancer subtype comprises determining whether the cancer is AR+ (ARPC) or AR-. In one embodiment, determining the cancer subtype comprises determining whether the cancer is ARPC or AR-low. In one embodiment, determining the cancer subtype comprises determining whether the cancer has a neuroendocrine prostate cancer (NEPC) phenotype signature. In one embodiment, determining the cancer subtype comprises determining whether the cancer is biendocrine. In one embodiment, determining the cancer subtype comprises determining whether the cancer is biendocrine.
whether the cancer is AR+ (ARPC) or AR-;
Whether the cancer is AR-low or ARPC;
whether the cancer has a neuroendocrine prostate cancer (NEPC) phenotypic signature;
Whether the cancer is AR-low or NEPC;
This includes determining two or all of whether the cancer is bisecretory, ARPC, or NEPC.

In one embodiment, the cancer is characterized as lung cancer. In one embodiment, determining the cancer subtype includes determining whether the cancer is small cell lung cancer (SCLC) or non-small cell lung cancer (NSCLC). In one embodiment, the method further includes determining whether the NSCLC is adenocarcinoma or squamous cell carcinoma. In one embodiment, the sequence read data is generated from a panel of genomic targets. In one embodiment, the panel of genomic targets includes transcription factor binding sites (TFBS) of one or more transcription factors associated with SCLC. In one embodiment, the one or more transcription factors associated with SCLC include one or more of ASLC, NEUROD1, POU2F3, REST, etc., and the method includes determining the nucleosome occupancy of the TFBS. In one embodiment, the TFBS is retained in the panel if it is proximal to a transcription start site of a gene associated with lung cancer, identified, for example, by ChIP-seq data. In one embodiment, the panel of genomic targets includes transcription start sites (TSSs) of one or more markers associated with lung cancer, and the method includes determining the nucleosome occupancy of the TSSs.

In some embodiments of any method aspect described herein relating to cancer detection or characterization, the sample is obtained from a subject. The method may further include administering an effective treatment to the subject based on the determined cancer subtype. In one embodiment, the method further includes performing the method on a plurality of samples obtained from the subject at a plurality of separate time points after an initial diagnosis of cancer. In one embodiment, the sequence read data is generated by ultra-low pass whole genome sequencing. In one embodiment, the sequence read data is generated by a chromatin accessibility assay. In one embodiment, the sequence read data is generated by ATAC-seq. In one embodiment, the sequence read data is generated by ChIP-seq. In one embodiment, the sequence read data is generated by a DNase sensitivity assay. In one embodiment, the sequence read data is generated by a cut-and-run assay. In one embodiment, the cut-and-run assay incorporates an affinity reagent that targets one or more of the post-translational modifications H3K27ac, H3K4me1, and H3K27ac.

In some embodiments of any of the method aspects described herein, the method may further include generating the sequence read data.

In some embodiments of any of the method aspects described herein, the sequence read data comprises sequence read data generated from a panel of genomic targets. In some embodiments of any of the method aspects described herein, the panel of genomic targets comprises transcription factor binding sites (TFBSs) of one or more transcription factors associated with the cancer type of interest. In some embodiments of any of the method aspects described herein, the method comprises determining the nucleosome occupancy of the TFBSs. In some embodiments of any of the method aspects described herein, the TFBSs are identified, such as by ChIP-seq data, and are retained in the panel if they are proximal to a transcription start site of a gene associated with the cancer type of interest. In some embodiments of any of the method aspects described herein, the panel of genomic targets comprises transcription start sites (TSSs) of one or more markers associated with the cancer type of interest, and the method comprises determining the nucleosome occupancy of the TSSs.

In some embodiments of any of the method aspects described herein, the sample may be blood, plasma, serum, or the like.

In another aspect, the present disclosure provides a computer-implemented method for augmenting sequence read data from a cell-free DNA sample to predict cell type, the method comprising:
receiving, by a computing system, sequence read data, the sequence read data including a plurality of fragment reads, each fragment read having a fragment length;
determining, by the computing system, a variability in fragment sizes for at least one gene associated with a cell type; and predicting, by the computing system, the cell type based on the variability in fragment sizes for the at least one gene.

In one embodiment, the step of determining the fragment size variability comprises determining a fragment size coefficient of variation. In one embodiment, the step of predicting the cell type based on the genome coverage distribution comprises predicting a cell phenotype. In one embodiment, predicting the cell phenotype comprises predicting a cancer subtype. In one embodiment, predicting the cell phenotype comprises predicting a cancer subtype of prostate cancer. In one embodiment, predicting the cancer subtype comprises distinguishing between ARPC and NEPC.

In one embodiment, predicting the cell type based on the variability of fragment sizes comprises:
generating one or more features based on the variability in fragment sizes;
providing the one or more features as an input to a classifier model; and determining the cell type based on an output of the classifier model.

In one embodiment, generating the one or more features based on the fragment size variability includes generating a log2 fold change value of the coefficient of variation of fragment size in a first cell type relative to a second cell type. In one embodiment, the log2 fold change value predicts at least one of gene expression and gene transcription activity between the first cell type and the second cell type. In one embodiment, the first cell type is an ARPC cell and the second cell type is a NEPC cell. In one embodiment, the classifier model includes a logistic regression model, an artificial neural network, a decision tree, a support vector machine, or a Bayesian network.

In another aspect, the disclosure provides a method for determining a cell type of a cell of interest from a sample comprising cell-free DNA derived from said cell of interest, comprising:
obtaining sequence read data generated from said sample comprising cell-free DNA;
The present invention provides a method comprising: performing a computer-implemented method described herein (e.g., relating to predicting the cell type based on the variability of the fragment sizes); and determining the cell type of the cell of interest based on the prediction provided by the computing system.

In one embodiment, the step of determining the cell type comprises determining a cell phenotype. In one embodiment, determining the cell phenotype comprises determining a cancer subtype. In one embodiment, determining the cancer subtype comprises distinguishing between ARPC and NEPC.

In another aspect, the disclosure provides a method of detecting the presence of cancer cells in a subject, comprising:
obtaining sequence read data generated from a sample comprising cell-free DNA obtained from the subject;
performing a computer-implemented method described herein (e.g., relating to predicting the cell type based on the variability of the fragment sizes); and determining the presence of cancer cells in the subject based on the prediction provided by the computing system.

In one embodiment, the method is performed multiple times over time, with each performance further characterizing the cancer cell(s) detected in the subject and determining a cancer subtype or phenotype of the detected cancer cell(s) based on the prediction provided by the computing system. In one embodiment, the method is performed multiple times over time, with the method further comprising detecting changes in the detected cancer cell phenotype over time. In one embodiment, the subject receives a cancer treatment between performances of the method, and the method further comprises determining the responsiveness of the cancer cells to the treatment.

In another aspect, the present disclosure provides a method for determining a cancer subtype of a target cancer cell from a sample comprising cell-free DNA derived from said target cancer cell, comprising:
obtaining sequence read data generated from said sample comprising cell-free DNA;
performing a computer-implemented method described herein (e.g., relating to predicting the cell type based on the variability of the fragment sizes); and determining the cell type of a cell of origin based on a predicted cancer subtype provided by the computing system.

In one embodiment, the sample is obtained from a subject having cancer. In one embodiment, the cancer is characterized as metastatic prostate cancer. In one embodiment, the determining the cancer subtype comprises determining whether the cancer is AR+ (ARPC) or AR-. In one embodiment, the determining the cancer subtype comprises determining whether the cancer is ARPC or AR-low prostate cancer (ARLPC). In one embodiment, the determining the cancer subtype comprises determining whether the cancer has a neuroendocrine prostate cancer (NEPC) phenotype signature. In one embodiment, the sample is obtained from a subject and the method further comprises administering an effective treatment to the subject based on the determined cancer subtype.

In one embodiment, the method further comprises performing the method on a plurality of samples obtained from the subject at a plurality of distinct time points after an initial diagnosis of cancer.

In one embodiment, the sequence read data is generated by ultra-low pass whole genome sequencing. In one embodiment, the sequence read data is generated by a chromatin accessibility assay. In one embodiment, the sequence read data is generated by ATAC-seq. In one embodiment, the sequence read data is generated by ChIP-seq. In one embodiment, the sequence read data is generated by a DNase sensitivity assay. In one embodiment, the sequence read data is generated by a cut-and-run assay. In one embodiment, the cut-and-run assay incorporates affinity reagents that target post-translational modifications to one or more of H3K27ac, H3K4me1, and H3K27ac. In one embodiment, the method further comprises the step of generating the sequence read data.

In one embodiment, the sequence read data is generated from a panel of genomic targets. In one embodiment, the panel of genomic targets includes transcription factor binding sites (TFBSs) for one or more transcription factors associated with the cancer type of interest. In one embodiment, the method includes determining the nucleosome occupancy of the TFBSs. In one embodiment, the TFBSs are identified, such as by ChIP-seq data, and are retained in the panel if they are proximal to a transcription start site of a gene associated with the cancer type of interest. In one embodiment, the panel of genomic targets includes transcription start sites (TSSs) for one or more markers associated with the cancer type of interest, and the method includes determining the nucleosome occupancy of the TSSs. In one embodiment, the sample is blood, plasma, or serum.

The foregoing aspects and many of the attendant advantages of the present invention will become more readily appreciated as the same become better understood by reference to the following detailed description taken in conjunction with the accompanying drawings.

FIG. 1 is a flow chart illustrating a non-limiting example of an embodiment of a method for predicting cancer subtype according to various aspects of the present disclosure.

FIG. 2 is a flow chart illustrating a non-limiting example of an embodiment of a procedure for determining and filtering sites that provide information about a tissue, cell type, cancer type, or cancer subtype of interest to identify sites that provide information specific to a cancer subtype, in accordance with various aspects of the present disclosure.

FIG. 3 is a flow chart illustrating a non-limiting example of an embodiment of a procedure for determining a GC occurrence matrix for a genome according to various aspects of the present disclosure.

FIG. 4 is a flow chart illustrating a non-limiting example of an embodiment of a procedure for determining a GC bias value for sequence read data using a GC occurrence matrix according to various aspects of the disclosure.

FIG. 5 is a flow chart illustrating a non-limiting example of an embodiment of a procedure for generating nucleosome profiles of sequence read data for sites that provide subtype-specific information using GC bias values, in accordance with various aspects of the present disclosure.

FIG. 6 is a block diagram illustrating aspects of an exemplary computing device suitable for use as a computing device of the present disclosure.

7A and 7B illustrate the Griffin framework for cfDNA nucleosome profiling to predict cancer subtypes and tumor phenotypes. FIG. 7A is a diagram of accessible (left panel) and inaccessible (right panel) sites, such as TFBS. Nucleosomes (gray) are arranged in an organized manner around accessible sites (boxes; left panel), but not around inaccessible sites (right panel). These nucleosomes protect DNA from degradation upon release into peripheral blood. The fragments protected from plasma are sequenced and aligned, leading to a coverage profile that reflects nucleosome protection in the cell of origin. FIG. 7B is a schematic diagram showing the Griffin workflow for cfDNA nucleosome profiling analysis. cfDNA whole genome sequencing (WGS) data with coverage ≥ 0.1x are aligned to the established hg38 genome. (1) For each sample, a fragment-based GC bias is computed for each fragment size. (2) Sites of interest are selected from any assay. Paired-end reads aligned to each site are collected, fragment midpoint coverage is counted, and corrected for GC bias to create a coverage profile. (3) Coverage profiles from all sites in a group (e.g., open chromatin for tumor subtypes) are averaged to create a composite coverage profile. The composite profile is normalized using the surrounding region (-5kb to +5kb). (4) Extract three features from the composite coverage profile: central coverage (coverage between -30 bp and +30 bp of the site; "a"), mean coverage (between -1 kb and +1 kb; "b"), and amplitude (calculated using a fast Fourier transform (FFT), "c"). 7A and 7B illustrate the Griffin framework for cfDNA nucleosome profiling to predict cancer subtypes and tumor phenotypes. FIG. 7A is a diagram of accessible (left panel) and inaccessible (right panel) sites, such as TFBS. Nucleosomes (gray) are arranged in an organized manner around accessible sites (boxes; left panel), but not around inaccessible sites (right panel). These nucleosomes protect DNA from degradation upon release into peripheral blood. The fragments protected from plasma are sequenced and aligned, leading to a coverage profile that reflects nucleosome protection in the cell of origin. FIG. 7B is a schematic diagram showing the Griffin workflow for cfDNA nucleosome profiling analysis. cfDNA whole genome sequencing (WGS) data with coverage ≥ 0.1x are aligned to the established hg38 genome. (1) For each sample, a fragment-based GC bias is computed for each fragment size. (2) Sites of interest are selected from any assay. Paired-end reads aligned to each site are collected, fragment midpoint coverage is counted, and corrected for GC bias to create a coverage profile. (3) Coverage profiles from all sites in a group (e.g., open chromatin for tumor subtypes) are averaged to create a composite coverage profile. The composite profile is normalized using the surrounding region (-5kb to +5kb). (4) Extract three features from the composite coverage profile: central coverage (coverage between -30 bp and +30 bp of the site; "a"), mean coverage (between -1 kb and +1 kb; "b"), and amplitude (calculated using a fast Fourier transform (FFT), "c"). 7A and 7B illustrate the Griffin framework for cfDNA nucleosome profiling to predict cancer subtypes and tumor phenotypes. FIG. 7A is a diagram of accessible (left panel) and inaccessible (right panel) sites, such as TFBS. Nucleosomes (gray) are arranged in an organized manner around accessible sites (boxes; left panel), but not around inaccessible sites (right panel). These nucleosomes protect DNA from degradation upon release into peripheral blood. The fragments protected from plasma are sequenced and aligned, leading to a coverage profile that reflects nucleosome protection in the cell of origin. FIG. 7B is a schematic diagram showing the Griffin workflow for cfDNA nucleosome profiling analysis. cfDNA whole genome sequencing (WGS) data with coverage ≥ 0.1x are aligned to the established hg38 genome. (1) For each sample, a fragment-based GC bias is computed for each fragment size. (2) Sites of interest are selected from any assay. Paired-end reads aligned to each site are collected, fragment midpoint coverage is counted, and corrected for GC bias to create a coverage profile. (3) Coverage profiles from all sites in a group (e.g., open chromatin for tumor subtypes) are averaged to create a composite coverage profile. The composite profile is normalized using the surrounding region (-5kb to +5kb). (4) Extract three features from the composite coverage profile: central coverage (coverage between -30 bp and +30 bp of the site; "a"), mean coverage (between -1 kb and +1 kb; "b"), and amplitude (calculated using a fast Fourier transform (FFT), "c").

8A-8G illustrate that Griffin GC bias correction improves detection of tissue-specific accessibility from cfDNA. FIG. 8A graphically illustrates a collection of 10,000 GRHL2 binding sites and the GC content in the surrounding 2 kb region. The mean GC content (line) and interquartile range (shading) are shown. FIG. 8B graphically illustrates that cfDNA GC bias is specific to each sample and each fragment length. Computed GC bias for cfDNA from a healthy donor sample (HD_46; dashed shading) and cfDNA from a metastatic breast cancer sample (MBC_315; solid shading) are shown for various fragment sizes. FIG. 8C graphically illustrates the combined coverage profiles of 10,000 GRHL2 binding sites before and after GC correction, shown for HD_46 (dashed line) and MBC_315 (solid line). Before GC correction, the "center coverage" values are higher due to the effect of GC bias, obscuring the differential signal between samples. After GC correction, the center coverage values of the MBC samples are lower, consistent with GRHL2 activity being increased in breast cancer but not in immune cells constituting healthy donor samples. FIG. 8D graphically illustrates the combined coverage profiles of 10,000 LYL1 sites before and after GC correction for two MBC samples using deep WGS (9-25×, orange), two healthy donors (17-20×, green), and 191 MBC samples using ULP-WGS (0.1-0.3×, blue). The median +/- IQR of the 191 ULP-WGS samples is shaded blue. The smaller "core coverage" in the healthy donor samples, corresponding to greater site accessibility, is predicted to be due to LYL1 being a transcription factor associated with hematopoiesis. FIG. 8E graphically illustrates the correlation of cfDNA tumor fraction with central coverage for LYL1, shown for ULP-WGS (0.1-0.3×, n=191) and WGS (9-25×, n=2) of MBC and healthy donor samples (17-20×, n=2). cfDNA contains a mixture of tumor and blood cells; therefore, central coverage values (lower representing increased accessibility) are expected to positively correlate with tumor fraction. After GC correction, the correlation (for MBC ULP-WGS samples) is much stronger based on Pearson's r correlation coefficient. Linear fit root mean square error (RMSE) is shown. FIG. 8F illustrates box plots showing the distribution of RMSE (linear fit between central coverage and tumor fraction in the MBC ULP-WGS dataset [0.1-0.3×, n=191]) across 377 TFs before and after GC correction. Boxed ranges represent median ± IQR, whiskers indicate range of non-outlier data (maximum is 1.5×IQR). Outliers are plotted in grey. p-values were calculated using Wilcoxon signed rank test (two-sided). FIG. 8G illustrates box plots showing the distribution of mean absolute deviation (of central coverage across 215 healthy donors [1-2×WGS]) across 377 TFs before and after GC correction. Box elements are the same as in (8F). p-values were calculated using Wilcoxon signed rank test (two-sided). 8A-8G illustrate that Griffin GC bias correction improves detection of tissue-specific accessibility from cfDNA. FIG. 8A graphically illustrates a collection of 10,000 GRHL2 binding sites and the GC content in the surrounding 2 kb region. The mean GC content (line) and interquartile range (shading) are shown. FIG. 8B graphically illustrates that cfDNA GC bias is specific to each sample and each fragment length. Computed GC bias for cfDNA from a healthy donor sample (HD_46; dashed shading) and cfDNA from a metastatic breast cancer sample (MBC_315; solid shading) are shown for various fragment sizes. FIG. 8C graphically illustrates the combined coverage profiles of 10,000 GRHL2 binding sites before and after GC correction, shown for HD_46 (dashed line) and MBC_315 (solid line). Before GC correction, the "center coverage" values are higher due to the effect of GC bias, obscuring the differential signal between samples. After GC correction, the center coverage values of the MBC samples are lower, consistent with GRHL2 activity being increased in breast cancer but not in immune cells constituting healthy donor samples. FIG. 8D graphically illustrates the combined coverage profiles of 10,000 LYL1 sites before and after GC correction for two MBC samples using deep WGS (9-25×, orange), two healthy donors (17-20×, green), and 191 MBC samples using ULP-WGS (0.1-0.3×, blue). The median +/- IQR of the 191 ULP-WGS samples is shaded blue. The smaller "core coverage" in the healthy donor samples, corresponding to greater site accessibility, is predicted to be due to LYL1 being a transcription factor associated with hematopoiesis. FIG. 8E graphically illustrates the correlation of cfDNA tumor fraction with central coverage for LYL1, shown for ULP-WGS (0.1-0.3×, n=191) and WGS (9-25×, n=2) of MBC and healthy donor samples (17-20×, n=2). cfDNA contains a mixture of tumor and blood cells; therefore, central coverage values (lower representing increased accessibility) are expected to positively correlate with tumor fraction. After GC correction, the correlation (for MBC ULP-WGS samples) is much stronger based on Pearson's r correlation coefficient. Linear fit root mean square error (RMSE) is shown. FIG. 8F illustrates box plots showing the distribution of RMSE (linear fit between central coverage and tumor fraction in the MBC ULP-WGS dataset [0.1-0.3×, n=191]) across 377 TFs before and after GC correction. Boxed ranges represent median ± IQR, whiskers indicate range of non-outlier data (maximum is 1.5×IQR). Outliers are plotted in grey. p-values were calculated using Wilcoxon signed rank test (two-sided). FIG. 8G illustrates box plots showing the distribution of mean absolute deviation (of central coverage across 215 healthy donors [1-2×WGS]) across 377 TFs before and after GC correction. Box elements are the same as in (8F). p-values were calculated using Wilcoxon signed rank test (two-sided). 8A-8G illustrate that Griffin GC bias correction improves detection of tissue-specific accessibility from cfDNA. FIG. 8A graphically illustrates a collection of 10,000 GRHL2 binding sites and the GC content in the surrounding 2 kb region. The mean GC content (line) and interquartile range (shading) are shown. FIG. 8B graphically illustrates that cfDNA GC bias is specific to each sample and each fragment length. Computed GC bias for cfDNA from a healthy donor sample (HD_46; dashed shading) and cfDNA from a metastatic breast cancer sample (MBC_315; solid shading) are shown for various fragment sizes. FIG. 8C graphically illustrates the combined coverage profiles of 10,000 GRHL2 binding sites before and after GC correction, shown for HD_46 (dashed line) and MBC_315 (solid line). Before GC correction, the "center coverage" values are higher due to the effect of GC bias, obscuring the differential signal between samples. After GC correction, the center coverage values of the MBC samples are lower, consistent with GRHL2 activity being increased in breast cancer but not in immune cells constituting healthy donor samples. FIG. 8D graphically illustrates the combined coverage profiles of 10,000 LYL1 sites before and after GC correction for two MBC samples using deep WGS (9-25×, orange), two healthy donors (17-20×, green), and 191 MBC samples using ULP-WGS (0.1-0.3×, blue). The median +/- IQR of the 191 ULP-WGS samples is shaded blue. The smaller "core coverage" in the healthy donor samples, corresponding to greater site accessibility, is predicted to be due to LYL1 being a transcription factor associated with hematopoiesis. FIG. 8E graphically illustrates the correlation of cfDNA tumor fraction with central coverage for LYL1, shown for ULP-WGS (0.1-0.3×, n=191) and WGS (9-25×, n=2) of MBC and healthy donor samples (17-20×, n=2). cfDNA contains a mixture of tumor and blood cells; therefore, central coverage values (lower representing increased accessibility) are expected to positively correlate with tumor fraction. After GC correction, the correlation (for MBC ULP-WGS samples) is much stronger based on Pearson's r correlation coefficient. Linear fit root mean square error (RMSE) is shown. FIG. 8F illustrates box plots showing the distribution of RMSE (linear fit between central coverage and tumor fraction in the MBC ULP-WGS dataset [0.1-0.3×, n=191]) across 377 TFs before and after GC correction. Boxed ranges represent median ± IQR, whiskers indicate range of non-outlier data (maximum is 1.5×IQR). Outliers are plotted in grey. p-values were calculated using Wilcoxon signed rank test (two-sided). FIG. 8G illustrates box plots showing the distribution of mean absolute deviation (of central coverage across 215 healthy donors [1-2×WGS]) across 377 TFs before and after GC correction. Box elements are the same as in (8F). p-values were calculated using Wilcoxon signed rank test (two-sided). 8A-8G illustrate that Griffin GC bias correction improves detection of tissue-specific accessibility from cfDNA. FIG. 8A graphically illustrates a collection of 10,000 GRHL2 binding sites and the GC content in the surrounding 2 kb region. The mean GC content (line) and interquartile range (shading) are shown. FIG. 8B graphically illustrates that cfDNA GC bias is specific to each sample and each fragment length. Computed GC bias for cfDNA from a healthy donor sample (HD_46; dashed shading) and cfDNA from a metastatic breast cancer sample (MBC_315; solid shading) are shown for various fragment sizes. FIG. 8C graphically illustrates the combined coverage profiles of 10,000 GRHL2 binding sites before and after GC correction, shown for HD_46 (dashed line) and MBC_315 (solid line). Before GC correction, the "center coverage" values are higher due to the effect of GC bias, obscuring the differential signal between samples. After GC correction, the center coverage values of the MBC samples are lower, consistent with GRHL2 activity being increased in breast cancer but not in immune cells constituting healthy donor samples. FIG. 8D graphically illustrates the combined coverage profiles of 10,000 LYL1 sites before and after GC correction for two MBC samples using deep WGS (9-25×, orange), two healthy donors (17-20×, green), and 191 MBC samples using ULP-WGS (0.1-0.3×, blue). The median +/- IQR of the 191 ULP-WGS samples is shown in blue shading. The smaller "core coverage" in the healthy donor samples, corresponding to greater site accessibility, is predicted to be due to LYL1 being a transcription factor associated with hematopoiesis. FIG. 8E graphically illustrates the correlation of cfDNA tumor fraction with central coverage for LYL1, shown for ULP-WGS (0.1-0.3×, n=191) and WGS (9-25×, n=2) of MBC and healthy donor samples (17-20×, n=2). cfDNA contains a mixture of tumor and blood cells; therefore, central coverage values (lower representing increased accessibility) are expected to positively correlate with tumor fraction. After GC correction, the correlation (for MBC ULP-WGS samples) is much stronger based on Pearson's r correlation coefficient. Linear fit root mean square error (RMSE) is shown. FIG. 8F illustrates box plots showing the distribution of RMSE (linear fit between central coverage and tumor fraction in the MBC ULP-WGS dataset [0.1-0.3×, n=191]) across 377 TFs before and after GC correction. Boxed ranges represent median ± IQR, whiskers indicate range of non-outlier data (maximum is 1.5×IQR). Outliers are plotted in grey. p-values were calculated using Wilcoxon signed rank test (two-sided). FIG. 8G illustrates box plots showing the distribution of mean absolute deviation (of central coverage across 215 healthy donors [1-2×WGS]) across 377 TFs before and after GC correction. Box elements are the same as in (8F). p-values were calculated using Wilcoxon signed rank test (two-sided). 8A-8G illustrate that Griffin GC bias correction improves detection of tissue-specific accessibility from cfDNA. FIG. 8A graphically illustrates a collection of 10,000 GRHL2 binding sites and the GC content in the surrounding 2 kb region. The mean GC content (line) and interquartile range (shading) are shown. FIG. 8B graphically illustrates that cfDNA GC bias is specific to each sample and each fragment length. Computed GC bias for cfDNA from a healthy donor sample (HD_46; dashed shading) and cfDNA from a metastatic breast cancer sample (MBC_315; solid shading) are shown for various fragment sizes. FIG. 8C graphically illustrates the combined coverage profiles of 10,000 GRHL2 binding sites before and after GC correction, shown for HD_46 (dashed line) and MBC_315 (solid line). Before GC correction, the "center coverage" values are higher due to the effect of GC bias, obscuring the differential signal between samples. After GC correction, the center coverage values of the MBC samples are lower, consistent with GRHL2 activity being increased in breast cancer but not in immune cells constituting healthy donor samples. FIG. 8D graphically illustrates the combined coverage profiles of 10,000 LYL1 sites before and after GC correction for two MBC samples using deep WGS (9-25×, orange), two healthy donors (17-20×, green), and 191 MBC samples using ULP-WGS (0.1-0.3×, blue). The median +/- IQR of the 191 ULP-WGS samples is shaded blue. The smaller "core coverage" in the healthy donor samples, corresponding to greater site accessibility, is predicted to be due to LYL1 being a transcription factor associated with hematopoiesis. FIG. 8E graphically illustrates the correlation of cfDNA tumor fraction with central coverage for LYL1, shown for ULP-WGS (0.1-0.3×, n=191) and WGS (9-25×, n=2) of MBC and healthy donor samples (17-20×, n=2). cfDNA contains a mixture of tumor and blood cells; therefore, central coverage values (lower representing increased accessibility) are expected to positively correlate with tumor fraction. After GC correction, the correlation (for MBC ULP-WGS samples) is much stronger based on Pearson's r correlation coefficient. Linear fit root mean square error (RMSE) is shown. FIG. 8F illustrates box plots showing the distribution of RMSE (linear fit between central coverage and tumor fraction in the MBC ULP-WGS dataset [0.1-0.3×, n=191]) across 377 TFs before and after GC correction. Boxed ranges represent median ± IQR, whiskers indicate range of non-outlier data (maximum is 1.5×IQR). Outliers are plotted in grey. p-values were calculated using Wilcoxon signed rank test (two-sided). FIG. 8G illustrates box plots showing the distribution of mean absolute deviation (of central coverage across 215 healthy donors [1-2×WGS]) across 377 TFs before and after GC correction. Box elements are the same as in (8F). p-values were calculated using Wilcoxon signed rank test (two-sided).

Figures 9A and 9B illustrate that Griffin enables accurate cancer detection and tissue of origin prediction. Figure 9A illustrates receiver operating characteristic (ROC) curves for logistic regression classification of cancer vs. healthy controls in three datasets: the DELFI dataset (Cristiano, S. et al. Genome-wide cell-free DNA fragmentation in patients with cancer. Nature 570, 385-389 (2019)), the LUCAS dataset, and the LUCAS validation dataset (Mathios, D. et al. Detection and characterization of lung cancer using cell-free DNA fragmentomes. Nat Commun 12, 5060 (2021)). For each dataset, the performance of both the original low-pass (1-2x) WGS and the ultra-low-pass (0.1x) WGS generated by in-silico downsampling is shown. Logistic regression was performed on the top PCA components explaining 80% of the variance in features (center coverage, mean coverage, and amplitude) extracted from nucleosome profiles around TFBS. The ROCs for each stage of cancer and health are shown. Figure 9B illustrates box plots of AUC values for 1000 bootstrap replicates. The boxed range represents the median ± IQR, and the whiskers represent the range of non-outlier data (maximum is 1.5 × IQR). Values below the box plots indicate the median and 95% confidence interval. Figures 9A and 9B illustrate that Griffin enables accurate cancer detection and tissue of origin prediction. Figure 9A illustrates receiver operating characteristic (ROC) curves for logistic regression classification of cancer vs. healthy controls in three datasets: the DELFI dataset (Cristiano, S. et al. Genome-wide cell-free DNA fragmentation in patients with cancer. Nature 570, 385-389 (2019)), the LUCAS dataset, and the LUCAS validation dataset (Mathios, D. et al. Detection and characterization of lung cancer using cell-free DNA fragmentomes. Nat Commun 12, 5060 (2021)). For each dataset, the performance of both the original low-pass (1-2x) WGS and the ultra-low-pass (0.1x) WGS generated by in-silico downsampling is shown. Logistic regression was performed on the top PCA components explaining 80% of the variance in features (center coverage, mean coverage, and amplitude) extracted from nucleosome profiles around TFBS. The ROCs for each stage of cancer and health are shown. Figure 9B illustrates box plots of AUC values for 1000 bootstrap replicates. The boxed range represents the median ± IQR, and the whiskers represent the range of non-outlier data (maximum is 1.5 × IQR). Values below the box plots indicate the median and 95% confidence interval. Figures 9A and 9B illustrate that Griffin enables accurate cancer detection and tissue of origin prediction. Figure 9A illustrates receiver operating characteristic (ROC) curves for logistic regression classification of cancer vs. healthy controls in three datasets: the DELFI dataset (Cristiano, S. et al. Genome-wide cell-free DNA fragmentation in patients with cancer. Nature 570, 385-389 (2019)), the LUCAS dataset, and the LUCAS validation dataset (Mathios, D. et al. Detection and characterization of lung cancer using cell-free DNA fragmentomes. Nat Commun 12, 5060 (2021)). For each dataset, the performance of both the original low-pass (1-2x) WGS and the ultra-low-pass (0.1x) WGS generated by in-silico downsampling is shown. Logistic regression was performed on the top PCA components explaining 80% of the variance in features (center coverage, mean coverage, and amplitude) extracted from nucleosome profiles around TFBS. The ROCs for each stage of cancer and health are shown. Figure 9B illustrates box plots of AUC values for 1000 bootstrap replicates. The boxed range represents the median ± IQR, and the whiskers represent the range of non-outlier data (maximum is 1.5 × IQR). Values below the box plots indicate the median and 95% confidence interval. Figures 9A and 9B illustrate that Griffin enables accurate cancer detection and tissue of origin prediction. Figure 9A illustrates receiver operating characteristic (ROC) curves for logistic regression classification of cancer vs. healthy controls in three datasets: the DELFI dataset (Cristiano, S. et al. Genome-wide cell-free DNA fragmentation in patients with cancer. Nature 570, 385-389 (2019)), the LUCAS dataset, and the LUCAS validation dataset (Mathios, D. et al. Detection and characterization of lung cancer using cell-free DNA fragmentomes. Nat Commun 12, 5060 (2021)). For each dataset, the performance of both the original low-pass (1-2x) WGS and the ultra-low-pass (0.1x) WGS generated by in-silico downsampling is shown. Logistic regression was performed on the top PCA components explaining 80% of the variance in features (center coverage, mean coverage, and amplitude) extracted from nucleosome profiles around TFBS. The ROCs for each stage of cancer and health are shown. Figure 9B illustrates box plots of AUC values for 1000 bootstrap replicates. The boxed range represents the median ± IQR, and the whiskers represent the range of non-outlier data (maximum is 1.5 × IQR). Values below the box plots indicate the median and 95% confidence interval.

Figures 10A-10H illustrate that Griffin enables accurate prediction of breast cancer estrogen receptor subtypes from ultra-low-pass WGS. Figure 10A: ER+ and ER-specific open chromatin sites were selected from an assay for transposase-accessible chromatin using sequencing (ATAC-seq) data from ER+ (n=44) and ER- (n=15) breast tumors in the Cancer Genome Atlas (TCGA) (Corces, MR et al. The chromatin accessibility landscape of primary human cancers. Science 362, eaav1898 (2018)). Differential sites were identified using DESeq2 software (Love, MI, et al. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15, 550 (2014)) and the q-score and log2 fold change for each site were calculated. Sites with q-scores < ^5x10-4 and _log2 fold changes >0.5 or <-0.5 were considered differential. Figure 10B illustrates the combined coverage profile (median ± IQR) for ER+ (n=18,240) and ER- (n=19,347) specific sites, shown for MBC patients (tumor fraction ≥0.1) separated by clinical ER status (ER+, n=50; ER-, n=51). Sites shared with hematopoietic cells were excluded from this figure (Satpathy, AT et al. Massively parallel single-cell chromatin landscapes of human immune cell development and intratumoral T cell exhaustion. Nature Biotechnology 37, 925-936 (2019)). Figure 10C illustrates a coMut (Crowdis, J., He, MX, Reardon, B. & Van Allen, EM CoMut: visualizing integrated molecular information with comutation plots. Bioinformatics 36, 4348-4349 (2020)) plot showing information on 101 MBC patients with tumor fraction >0.10. The top row shows the ER status used to train and evaluate the regression model. This was the metastatic ER status obtained by IHC for the majority of patients, and primary ER status was used when metastatic ER status was unavailable. ER low (1-10% ER staining) was considered ER positive. Second row, left upward triangle contains primary ER status, right downward triangle contains metastatic ER status. ER low biopsies are shown in light blue, unknown status in white. Third row: tumor fraction, fraction of cfDNA originating from tumor, calculated using ichorcNA. Fourth row, median probability of ER positivity calculated by Griffin over 1000 bootstrap iterations. FIG. 10D is a receiver operating characteristic (ROC) curve of a logistic regression model predicting ER+ and ER- subtypes. ROC curves, accuracy and AUC are shown for all patients and for patients grouped by tumor fraction (TFx) between 0.05-0.1 and 0.1 or greater. 95% CI was obtained by bootstrapping. For patients with multiple samples, the first sample with tumor fraction >0.05 was used. FIG. 10E graphically illustrates the performance of the model on samples from the three validation cohorts. For patients with multiple time points, the first sample was used. Figure 10F graphically illustrates subtype prediction in patients separated by clinical metastatic ER status and clinical primary tumor ER status. P values were calculated using Fisher's exact test (2-sided). Figure 10G illustrates the ROC curve for predicting ER loss among patients with primary ER-positive tumors. 95% CI was obtained by bootstrapping. Figure 10H illustrates the time series for two patients (MBC1413 and MBC1099) with different subtypes and multiple cfDNA samples across multiple biopsies. ER+ predicted probability (thick grey line) and tumor fraction (thin grey line) for all cfDNA samples passing thresholds of tumor fraction >0.05 and 0.1x coverage are shown. The decision boundaries for ER+ (>0.5) and ER- (<0.5) are shown as dotted lines. Time series in months from metastatic diagnosis to death for each patient are shown. For patient MBC_1413, a metastatic biopsy (pleural fluid) was obtained on the day the metastasis was diagnosed and showed ER- disease. However, approximately 7 months later, another metastatic biopsy (liver) showed weak ER+ staining (5%). A final biopsy (pleural fluid) was obtained at approximately 12 months and again showed ER- staining. Plasma for cfDNA was removed between the second and third metastatic biopsies, with one final removal after the third biopsy. For patient MBC_1099, two ER- biopsies were obtained at 0 months (bone) and 7 months (liver). cfDNA was removed after this time point, but another biopsy (liver) showed the presence of low levels of ER+ disease between the two cfDNA removals. Figures 10A-10H illustrate that Griffin enables accurate prediction of breast cancer estrogen receptor subtypes from ultra-low-pass WGS. Figure 10A: ER+ and ER-specific open chromatin sites were selected from an assay for transposase-accessible chromatin using sequencing (ATAC-seq) data from ER+ (n=44) and ER- (n=15) breast tumors in the Cancer Genome Atlas (TCGA) (Corces, MR et al. The chromatin accessibility landscape of primary human cancers. Science 362, eaav1898 (2018)). Differential sites were identified using DESeq2 software (Love, MI, et al. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15, 550 (2014)) and the q-score and log2 fold change for each site were calculated. Sites with q-scores < ^5x10-4 and _log2 fold changes >0.5 or <-0.5 were considered differential. Figure 10B illustrates the combined coverage profile (median ± IQR) for ER+ (n=18,240) and ER- (n=19,347) specific sites, shown for MBC patients (tumor fraction ≥0.1) separated by clinical ER status (ER+, n=50; ER-, n=51). Sites shared with hematopoietic cells were excluded from this figure (Satpathy, AT et al. Massively parallel single-cell chromatin landscapes of human immune cell development and intratumoral T cell exhaustion. Nature Biotechnology 37, 925-936 (2019)). Figure 10C illustrates a coMut (Crowdis, J., He, MX, Reardon, B. & Van Allen, EM CoMut: visualizing integrated molecular information with comutation plots. Bioinformatics 36, 4348-4349 (2020)) plot showing information on 101 MBC patients with tumor fraction >0.10. The top row shows the ER status used to train and evaluate the regression model. This was the metastatic ER status obtained by IHC for the majority of patients, and primary ER status was used when metastatic ER status was unavailable. ER low (1-10% ER staining) was considered ER positive. Second row, left upward triangle contains primary ER status, right downward triangle contains metastatic ER status. ER low biopsies are shown in light blue, unknown status in white. Third row: tumor fraction, fraction of cfDNA originating from tumor, calculated using ichorcNA. Fourth row, median probability of ER positivity calculated by Griffin over 1000 bootstrap iterations. FIG. 10D is a receiver operating characteristic (ROC) curve of a logistic regression model predicting ER+ and ER- subtypes. ROC curves, accuracy and AUC are shown for all patients and for patients grouped by tumor fraction (TFx) between 0.05-0.1 and 0.1 or greater. 95% CI was obtained by bootstrapping. For patients with multiple samples, the first sample with tumor fraction >0.05 was used. FIG. 10E graphically illustrates the performance of the model on samples from the three validation cohorts. For patients with multiple time points, the first sample was used. Figure 10F graphically illustrates subtype prediction in patients separated by clinical metastatic ER status and clinical primary tumor ER status. P values were calculated using Fisher's exact test (2-sided). Figure 10G illustrates the ROC curve for predicting ER loss among patients with primary ER-positive tumors. 95% CI was obtained by bootstrapping. Figure 10H illustrates the time series for two patients (MBC1413 and MBC1099) with different subtypes and multiple cfDNA samples across multiple biopsies. ER+ predicted probability (thick grey line) and tumor fraction (thin grey line) for all cfDNA samples passing thresholds of tumor fraction >0.05 and 0.1x coverage are shown. The decision boundaries for ER+ (>0.5) and ER- (<0.5) are shown as dotted lines. Time series in months from metastatic diagnosis to death for each patient are shown. For patient MBC_1413, a metastatic biopsy (pleural fluid) was obtained on the day the metastasis was diagnosed and showed ER- disease. However, approximately 7 months later, another metastatic biopsy (liver) showed weak ER+ staining (5%). A final biopsy (pleural fluid) was obtained at approximately 12 months and again showed ER- staining. Plasma for cfDNA was removed between the second and third metastatic biopsies, with one final removal after the third biopsy. For patient MBC_1099, two ER- biopsies were obtained at 0 months (bone) and 7 months (liver). cfDNA was removed after this time point, but another biopsy (liver) showed the presence of low levels of ER+ disease between the two cfDNA removals. Figures 10A-10H illustrate that Griffin enables accurate prediction of breast cancer estrogen receptor subtypes from ultra-low-pass WGS. Figure 10A: ER+ and ER-specific open chromatin sites were selected from an assay for transposase-accessible chromatin using sequencing (ATAC-seq) data from ER+ (n=44) and ER- (n=15) breast tumors in the Cancer Genome Atlas (TCGA) (Corces, MR et al. The chromatin accessibility landscape of primary human cancers. Science 362, eaav1898 (2018)). Differential sites were identified using DESeq2 software (Love, MI, et al. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15, 550 (2014)) and the q-score and log2 fold change for each site were calculated. Sites with q-scores < ^5x10-4 and _log2 fold changes >0.5 or <-0.5 were considered differential. Figure 10B illustrates the combined coverage profile (median ± IQR) for ER+ (n=18,240) and ER- (n=19,347) specific sites, shown for MBC patients (tumor fraction ≥0.1) separated by clinical ER status (ER+, n=50; ER-, n=51). Sites shared with hematopoietic cells were excluded from this figure (Satpathy, AT et al. Massively parallel single-cell chromatin landscapes of human immune cell development and intratumoral T cell exhaustion. Nature Biotechnology 37, 925-936 (2019)). Figure 10C illustrates a coMut (Crowdis, J., He, MX, Reardon, B. & Van Allen, EM CoMut: visualizing integrated molecular information with comutation plots. Bioinformatics 36, 4348-4349 (2020)) plot showing information on 101 MBC patients with tumor fraction >0.10. The top row shows the ER status used to train and evaluate the regression model. This was the metastatic ER status obtained by IHC for the majority of patients, and primary ER status was used when metastatic ER status was unavailable. ER low (1-10% ER staining) was considered ER positive. Second row, left upward triangle contains primary ER status, right downward triangle contains metastatic ER status. ER low biopsies are shown in light blue, unknown status in white. Third row: tumor fraction, fraction of cfDNA originating from tumor, calculated using ichorcNA. Fourth row, median probability of ER positivity calculated by Griffin over 1000 bootstrap iterations. FIG. 10D is a receiver operating characteristic (ROC) curve of a logistic regression model predicting ER+ and ER- subtypes. ROC curves, accuracy and AUC are shown for all patients and for patients grouped by tumor fraction (TFx) between 0.05-0.1 and 0.1 or greater. 95% CI was obtained by bootstrapping. For patients with multiple samples, the first sample with tumor fraction >0.05 was used. FIG. 10E graphically illustrates the performance of the model on samples from the three validation cohorts. For patients with multiple time points, the first sample was used. Figure 10F graphically illustrates subtype prediction in patients separated by clinical metastatic ER status and clinical primary tumor ER status. P values were calculated using Fisher's exact test (2-sided). Figure 10G illustrates the ROC curve for predicting ER loss among patients with primary ER-positive tumors. 95% CI was obtained by bootstrapping. Figure 10H illustrates the time series for two patients (MBC1413 and MBC1099) with different subtypes and multiple cfDNA samples across multiple biopsies. ER+ predicted probability (thick grey line) and tumor fraction (thin grey line) for all cfDNA samples passing thresholds of tumor fraction >0.05 and 0.1x coverage are shown. The decision boundaries for ER+ (>0.5) and ER- (<0.5) are shown as dotted lines. Time series in months from metastatic diagnosis to death for each patient are shown. For patient MBC_1413, a metastatic biopsy (pleural fluid) was obtained on the day the metastasis was diagnosed and showed ER- disease. However, approximately 7 months later, another metastatic biopsy (liver) showed weak ER+ staining (5%). A final biopsy (pleural fluid) was obtained at approximately 12 months and again showed ER- staining. Plasma for cfDNA was removed between the second and third metastatic biopsies, with one final removal after the third biopsy. For patient MBC_1099, two ER- biopsies were obtained at 0 months (bone) and 7 months (liver). cfDNA was removed after this time point, but another biopsy (liver) showed the presence of low levels of ER+ disease between the two cfDNA removals. Figures 10A-10H illustrate that Griffin enables accurate prediction of breast cancer estrogen receptor subtypes from ultra-low-pass WGS. Figure 10A: ER+ and ER-specific open chromatin sites were selected from an assay for transposase-accessible chromatin using sequencing (ATAC-seq) data from ER+ (n=44) and ER- (n=15) breast tumors in the Cancer Genome Atlas (TCGA) (Corces, MR et al. The chromatin accessibility landscape of primary human cancers. Science 362, eaav1898 (2018)). Differential sites were identified using DESeq2 software (Love, MI, et al. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15, 550 (2014)) and the q-score and log2 fold change for each site were calculated. Sites with q-scores < ^5x10-4 and _log2 fold changes >0.5 or <-0.5 were considered differential. Figure 10B illustrates the combined coverage profile (median ± IQR) for ER+ (n=18,240) and ER- (n=19,347) specific sites, shown for MBC patients (tumor fraction ≥0.1) separated by clinical ER status (ER+, n=50; ER-, n=51). Sites shared with hematopoietic cells were excluded from this figure (Satpathy, AT et al. Massively parallel single-cell chromatin landscapes of human immune cell development and intratumoral T cell exhaustion. Nature Biotechnology 37, 925-936 (2019)). Figure 10C illustrates a coMut (Crowdis, J., He, MX, Reardon, B. & Van Allen, EM CoMut: visualizing integrated molecular information with comutation plots. Bioinformatics 36, 4348-4349 (2020)) plot showing information on 101 MBC patients with tumor fraction >0.10. The top row shows the ER status used to train and evaluate the regression model. This was the metastatic ER status obtained by IHC for the majority of patients, and primary ER status was used when metastatic ER status was unavailable. ER low (1-10% ER staining) was considered ER positive. Second row, left upward triangle contains primary ER status, right downward triangle contains metastatic ER status. ER low biopsies are shown in light blue, unknown status in white. Third row: tumor fraction, fraction of cfDNA originating from tumor, calculated using ichorcNA. Fourth row, median probability of ER positivity calculated by Griffin over 1000 bootstrap iterations. FIG. 10D is a receiver operating characteristic (ROC) curve of a logistic regression model predicting ER+ and ER- subtypes. ROC curves, accuracy and AUC are shown for all patients and for patients grouped by tumor fraction (TFx) between 0.05-0.1 and 0.1 or greater. 95% CI was obtained by bootstrapping. For patients with multiple samples, the first sample with tumor fraction >0.05 was used. FIG. 10E graphically illustrates the performance of the model on samples from the three validation cohorts. For patients with multiple time points, the first sample was used. Figure 10F graphically illustrates subtype prediction in patients separated by clinical metastatic ER status and clinical primary tumor ER status. P values were calculated using Fisher's exact test (2-sided). Figure 10G illustrates the ROC curve for predicting ER loss among patients with primary ER-positive tumors. 95% CI was obtained by bootstrapping. Figure 10H illustrates the time series for two patients (MBC1413 and MBC1099) with different subtypes and multiple cfDNA samples across multiple biopsies. ER+ predicted probability (thick grey line) and tumor fraction (thin grey line) for all cfDNA samples passing thresholds of tumor fraction >0.05 and 0.1x coverage are shown. The decision boundaries for ER+ (>0.5) and ER- (<0.5) are shown as dotted lines. Time series in months from metastatic diagnosis to death for each patient are shown. For patient MBC_1413, a metastatic biopsy (pleural fluid) was obtained on the day the metastasis was diagnosed and showed ER- disease. However, approximately 7 months later, another metastatic biopsy (liver) showed weak ER+ staining (5%). A final biopsy (pleural fluid) was obtained at approximately 12 months and again showed ER- staining. Plasma for cfDNA was removed between the second and third metastatic biopsies, with one final removal after the third biopsy. For patient MBC_1099, two ER- biopsies were obtained at 0 months (bone) and 7 months (liver). cfDNA was removed after this time point, but another biopsy (liver) showed the presence of low levels of ER+ disease between the two cfDNA removals. Figures 10A-10H illustrate that Griffin enables accurate prediction of breast cancer estrogen receptor subtypes from ultra-low-pass WGS. Figure 10A: ER+ and ER-specific open chromatin sites were selected from an assay for transposase-accessible chromatin using sequencing (ATAC-seq) data from ER+ (n=44) and ER- (n=15) breast tumors in the Cancer Genome Atlas (TCGA) (Corces, MR et al. The chromatin accessibility landscape of primary human cancers. Science 362, eaav1898 (2018)). Differential sites were identified using DESeq2 software (Love, MI, et al. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15, 550 (2014)) and the q-score and log2 fold change for each site were calculated. Sites with q-scores < ^5x10-4 and _log2 fold changes >0.5 or <-0.5 were considered differential. Figure 10B illustrates the combined coverage profile (median ± IQR) for ER+ (n=18,240) and ER- (n=19,347) specific sites, shown for MBC patients (tumor fraction ≥0.1) separated by clinical ER status (ER+, n=50; ER-, n=51). Sites shared with hematopoietic cells were excluded from this figure (Satpathy, AT et al. Massively parallel single-cell chromatin landscapes of human immune cell development and intratumoral T cell exhaustion. Nature Biotechnology 37, 925-936 (2019)). Figure 10C illustrates a coMut (Crowdis, J., He, MX, Reardon, B. & Van Allen, EM CoMut: visualizing integrated molecular information with comutation plots. Bioinformatics 36, 4348-4349 (2020)) plot showing information on 101 MBC patients with tumor fraction >0.10. The top row shows the ER status used to train and evaluate the regression model. This was the metastatic ER status obtained by IHC for the majority of patients, and primary ER status was used when metastatic ER status was unavailable. ER low (1-10% ER staining) was considered ER positive. Second row, left upward triangle contains primary ER status, right downward triangle contains metastatic ER status. ER low biopsies are shown in light blue, unknown status in white. Third row: tumor fraction, fraction of cfDNA originating from tumor, calculated using ichorcNA. Fourth row, median probability of ER positivity calculated by Griffin over 1000 bootstrap iterations. FIG. 10D is a receiver operating characteristic (ROC) curve of a logistic regression model predicting ER+ and ER- subtypes. ROC curves, accuracy and AUC are shown for all patients and for patients grouped by tumor fraction (TFx) between 0.05-0.1 and 0.1 or greater. 95% CI was obtained by bootstrapping. For patients with multiple samples, the first sample with tumor fraction >0.05 was used. FIG. 10E graphically illustrates the performance of the model on samples from the three validation cohorts. For patients with multiple time points, the first sample was used. Figure 10F graphically illustrates subtype prediction in patients separated by clinical metastatic ER status and clinical primary tumor ER status. P values were calculated using Fisher's exact test (2-sided). Figure 10G illustrates the ROC curve for predicting ER loss among patients with primary ER-positive tumors. 95% CI was obtained by bootstrapping. Figure 10H illustrates the time series for two patients (MBC1413 and MBC1099) with different subtypes and multiple cfDNA samples across multiple biopsies. ER+ predicted probability (thick grey line) and tumor fraction (thin grey line) for all cfDNA samples passing thresholds of tumor fraction >0.05 and 0.1x coverage are shown. The decision boundaries for ER+ (>0.5) and ER- (<0.5) are shown as dotted lines. Time series in months from metastatic diagnosis to death for each patient are shown. For patient MBC_1413, a metastatic biopsy (pleural fluid) was obtained on the day the metastasis was diagnosed and showed ER- disease. However, approximately 7 months later, another metastatic biopsy (liver) showed weak ER+ staining (5%). A final biopsy (pleural fluid) was obtained at approximately 12 months and again showed ER- staining. Plasma for cfDNA was removed between the second and third metastatic biopsies, with one final removal after the third biopsy. For patient MBC_1099, two ER- biopsies were obtained at 0 months (bone) and 7 months (liver). cfDNA was removed after this time point, but another biopsy (liver) showed the presence of low levels of ER+ disease between the two cfDNA removals. Figures 10A-10H illustrate that Griffin enables accurate prediction of breast cancer estrogen receptor subtypes from ultra-low-pass WGS. Figure 10A: ER+ and ER-specific open chromatin sites were selected from an assay for transposase-accessible chromatin using sequencing (ATAC-seq) data from ER+ (n=44) and ER- (n=15) breast tumors in the Cancer Genome Atlas (TCGA) (Corces, MR et al. The chromatin accessibility landscape of primary human cancers. Science 362, eaav1898 (2018)). Differential sites were identified using DESeq2 software (Love, MI, et al. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15, 550 (2014)) and the q-score and log2 fold change for each site were calculated. Sites with q-scores < ^5x10-4 and _log2 fold changes >0.5 or <-0.5 were considered differential. Figure 10B illustrates the combined coverage profile (median ± IQR) for ER+ (n=18,240) and ER- (n=19,347) specific sites, shown for MBC patients (tumor fraction ≥0.1) separated by clinical ER status (ER+, n=50; ER-, n=51). Sites shared with hematopoietic cells were excluded from this figure (Satpathy, AT et al. Massively parallel single-cell chromatin landscapes of human immune cell development and intratumoral T cell exhaustion. Nature Biotechnology 37, 925-936 (2019)). Figure 10C illustrates a coMut (Crowdis, J., He, MX, Reardon, B. & Van Allen, EM CoMut: visualizing integrated molecular information with comutation plots. Bioinformatics 36, 4348-4349 (2020)) plot showing information on 101 MBC patients with tumor fraction >0.10. The top row shows the ER status used to train and evaluate the regression model. This was the metastatic ER status obtained by IHC for the majority of patients, and primary ER status was used when metastatic ER status was unavailable. ER low (1-10% ER staining) was considered ER positive. Second row, left upward triangle contains primary ER status, right downward triangle contains metastatic ER status. ER low biopsies are shown in light blue, unknown status in white. Third row: tumor fraction, fraction of cfDNA originating from tumor, calculated using ichorcNA. Fourth row, median probability of ER positivity calculated by Griffin over 1000 bootstrap iterations. FIG. 10D is a receiver operating characteristic (ROC) curve of a logistic regression model predicting ER+ and ER- subtypes. ROC curves, accuracy and AUC are shown for all patients and for patients grouped by tumor fraction (TFx) between 0.05-0.1 and 0.1 or greater. 95% CI was obtained by bootstrapping. For patients with multiple samples, the first sample with tumor fraction >0.05 was used. FIG. 10E graphically illustrates the performance of the model on samples from the three validation cohorts. For patients with multiple time points, the first sample was used. Figure 10F graphically illustrates subtype prediction in patients separated by clinical metastatic ER status and clinical primary tumor ER status. P values were calculated using Fisher's exact test (2-sided). Figure 10G illustrates the ROC curve for predicting ER loss among patients with primary ER-positive tumors. 95% CI was obtained by bootstrapping. Figure 10H illustrates the time series for two patients (MBC1413 and MBC1099) with different subtypes and multiple cfDNA samples across multiple biopsies. ER+ predicted probability (thick grey line) and tumor fraction (thin grey line) for all cfDNA samples passing thresholds of tumor fraction >0.05 and 0.1x coverage are shown. The decision boundaries for ER+ (>0.5) and ER- (<0.5) are shown as dotted lines. Time series in months from metastatic diagnosis to death for each patient are shown. For patient MBC_1413, a metastatic biopsy (pleural fluid) was obtained on the day the metastasis was diagnosed and showed ER- disease. However, approximately 7 months later, another metastatic biopsy (liver) showed weak ER+ staining (5%). A final biopsy (pleural fluid) was obtained at approximately 12 months and again showed ER- staining. Plasma for cfDNA was removed between the second and third metastatic biopsies, with one final removal after the third biopsy. For patient MBC_1099, two ER- biopsies were obtained at 0 months (bone) and 7 months (liver). cfDNA was removed after this time point, but another biopsy (liver) showed the presence of low levels of ER+ disease between the two cfDNA removals. Figures 10A-10H illustrate that Griffin enables accurate prediction of breast cancer estrogen receptor subtypes from ultra-low-pass WGS. Figure 10A: ER+ and ER-specific open chromatin sites were selected from an assay for transposase-accessible chromatin using sequencing (ATAC-seq) data from ER+ (n=44) and ER- (n=15) breast tumors in the Cancer Genome Atlas (TCGA) (Corces, MR et al. The chromatin accessibility landscape of primary human cancers. Science 362, eaav1898 (2018)). Differential sites were identified using DESeq2 software (Love, MI, et al. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15, 550 (2014)) and the q-score and log2 fold change for each site were calculated. Sites with q-scores < ^5x10-4 and _log2 fold changes >0.5 or <-0.5 were considered differential. Figure 10B illustrates the combined coverage profile (median ± IQR) for ER+ (n=18,240) and ER- (n=19,347) specific sites, shown for MBC patients (tumor fraction ≥0.1) separated by clinical ER status (ER+, n=50; ER-, n=51). Sites shared with hematopoietic cells were excluded from this figure (Satpathy, AT et al. Massively parallel single-cell chromatin landscapes of human immune cell development and intratumoral T cell exhaustion. Nature Biotechnology 37, 925-936 (2019)). Figure 10C illustrates a coMut (Crowdis, J., He, MX, Reardon, B. & Van Allen, EM CoMut: visualizing integrated molecular information with comutation plots. Bioinformatics 36, 4348-4349 (2020)) plot showing information on 101 MBC patients with tumor fraction >0.10. The top row shows the ER status used to train and evaluate the regression model. This was the metastatic ER status obtained by IHC for the majority of patients, and primary ER status was used when metastatic ER status was unavailable. ER low (1-10% ER staining) was considered ER positive. Second row, left upward triangle contains primary ER status, right downward triangle contains metastatic ER status. ER low biopsies are shown in light blue, unknown status in white. Third row: tumor fraction, fraction of cfDNA originating from tumor, calculated using ichorcNA. Fourth row, median probability of ER positivity calculated by Griffin over 1000 bootstrap iterations. FIG. 10D is a receiver operating characteristic (ROC) curve of a logistic regression model predicting ER+ and ER- subtypes. ROC curves, accuracy and AUC are shown for all patients and for patients grouped by tumor fraction (TFx) between 0.05-0.1 and 0.1 or greater. 95% CI was obtained by bootstrapping. For patients with multiple samples, the first sample with tumor fraction >0.05 was used. FIG. 10E graphically illustrates the performance of the model on samples from the three validation cohorts. For patients with multiple time points, the first sample was used. Figure 10F graphically illustrates subtype prediction in patients separated by clinical metastatic ER status and clinical primary tumor ER status. P values were calculated using Fisher's exact test (2-sided). Figure 10G illustrates the ROC curve for predicting ER loss among patients with primary ER-positive tumors. 95% CI was obtained by bootstrapping. Figure 10H illustrates the time series for two patients (MBC1413 and MBC1099) with different subtypes and multiple cfDNA samples across multiple biopsies. ER+ predicted probability (thick grey line) and tumor fraction (thin grey line) for all cfDNA samples passing thresholds of tumor fraction >0.05 and 0.1x coverage are shown. The decision boundaries for ER+ (>0.5) and ER- (<0.5) are shown as dotted lines. Time series in months from metastatic diagnosis to death for each patient are shown. For patient MBC_1413, a metastatic biopsy (pleural fluid) was obtained on the day the metastasis was diagnosed and showed ER- disease. However, approximately 7 months later, another metastatic biopsy (liver) showed weak ER+ staining (5%). A final biopsy (pleural fluid) was obtained at approximately 12 months and again showed ER- staining. Plasma for cfDNA was removed between the second and third metastatic biopsies, with one final removal after the third biopsy. For patient MBC_1099, two ER- biopsies were obtained at 0 months (bone) and 7 months (liver). cfDNA was removed after this time point, but another biopsy (liver) showed the presence of low levels of ER+ disease between the two cfDNA removals. Figures 10A-10H illustrate that Griffin enables accurate prediction of breast cancer estrogen receptor subtypes from ultra-low-pass WGS. Figure 10A: ER+ and ER-specific open chromatin sites were selected from an assay for transposase-accessible chromatin using sequencing (ATAC-seq) data from ER+ (n=44) and ER- (n=15) breast tumors in the Cancer Genome Atlas (TCGA) (Corces, MR et al. The chromatin accessibility landscape of primary human cancers. Science 362, eaav1898 (2018)). Differential sites were identified using DESeq2 software (Love, MI, et al. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15, 550 (2014)) and the q-score and log2 fold change for each site were calculated. Sites with q-scores < ^5x10-4 and _log2 fold changes >0.5 or <-0.5 were considered differential. Figure 10B illustrates the combined coverage profile (median ± IQR) for ER+ (n=18,240) and ER- (n=19,347) specific sites, shown for MBC patients (tumor fraction ≥0.1) separated by clinical ER status (ER+, n=50; ER-, n=51). Sites shared with hematopoietic cells were excluded from this figure (Satpathy, AT et al. Massively parallel single-cell chromatin landscapes of human immune cell development and intratumoral T cell exhaustion. Nature Biotechnology 37, 925-936 (2019)). Figure 10C illustrates a coMut (Crowdis, J., He, MX, Reardon, B. & Van Allen, EM CoMut: visualizing integrated molecular information with comutation plots. Bioinformatics 36, 4348-4349 (2020)) plot showing information on 101 MBC patients with tumor fraction >0.10. The top row shows the ER status used to train and evaluate the regression model. This was the metastatic ER status obtained by IHC for the majority of patients, and primary ER status was used when metastatic ER status was unavailable. ER low (1-10% ER staining) was considered ER positive. Second row, left upward triangle contains primary ER status, right downward triangle contains metastatic ER status. ER low biopsies are shown in light blue, unknown status in white. Third row: tumor fraction, fraction of cfDNA originating from tumor, calculated using ichorcNA. Fourth row, median probability of ER positivity calculated by Griffin over 1000 bootstrap iterations. FIG. 10D is a receiver operating characteristic (ROC) curve of a logistic regression model predicting ER+ and ER- subtypes. ROC curves, accuracy and AUC are shown for all patients and for patients grouped by tumor fraction (TFx) between 0.05-0.1 and 0.1 or greater. 95% CI was obtained by bootstrapping. For patients with multiple samples, the first sample with tumor fraction >0.05 was used. FIG. 10E graphically illustrates the performance of the model on samples from the three validation cohorts. For patients with multiple time points, the first sample was used. Figure 10F graphically illustrates subtype prediction in patients separated by clinical metastatic ER status and clinical primary tumor ER status. P values were calculated using Fisher's exact test (2-sided). Figure 10G illustrates the ROC curve for predicting ER loss among patients with primary ER-positive tumors. 95% CI was obtained by bootstrapping. Figure 10H illustrates the time series for two patients (MBC1413 and MBC1099) with different subtypes and multiple cfDNA samples across multiple biopsies. ER+ predicted probability (thick grey line) and tumor fraction (thin grey line) for all cfDNA samples passing thresholds of tumor fraction >0.05 and 0.1x coverage are shown. The decision boundaries for ER+ (>0.5) and ER- (<0.5) are shown as dotted lines. Time series in months from metastatic diagnosis to death for each patient are shown. For patient MBC_1413, a metastatic biopsy (pleural fluid) was obtained on the day the metastasis was diagnosed and showed ER- disease. However, approximately 7 months later, another metastatic biopsy (liver) showed weak ER+ staining (5%). A final biopsy (pleural fluid) was obtained at approximately 12 months and again showed ER- staining. Plasma for cfDNA was removed between the second and third metastatic biopsies, with one final removal after the third biopsy. For patient MBC_1099, two ER- biopsies were obtained at 0 months (bone) and 7 months (liver). cfDNA was removed after this time point, but another biopsy (liver) showed the presence of low levels of ER+ disease between the two cfDNA removals. Figures 10A-10H illustrate that Griffin enables accurate prediction of breast cancer estrogen receptor subtypes from ultra-low-pass WGS. Figure 10A: ER+ and ER-specific open chromatin sites were selected from an assay for transposase-accessible chromatin using sequencing (ATAC-seq) data from ER+ (n=44) and ER- (n=15) breast tumors in the Cancer Genome Atlas (TCGA) (Corces, MR et al. The chromatin accessibility landscape of primary human cancers. Science 362, eaav1898 (2018)). Differential sites were identified using DESeq2 software (Love, MI, et al. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15, 550 (2014)) and the q-score and log2 fold change for each site were calculated. Sites with q-scores < ^5x10-4 and _log2 fold changes >0.5 or <-0.5 were considered differential. Figure 10B illustrates the combined coverage profile (median ± IQR) for ER+ (n=18,240) and ER- (n=19,347) specific sites, shown for MBC patients (tumor fraction ≥0.1) separated by clinical ER status (ER+, n=50; ER-, n=51). Sites shared with hematopoietic cells were excluded from this figure (Satpathy, AT et al. Massively parallel single-cell chromatin landscapes of human immune cell development and intratumoral T cell exhaustion. Nature Biotechnology 37, 925-936 (2019)). Figure 10C illustrates a coMut (Crowdis, J., He, MX, Reardon, B. & Van Allen, EM CoMut: visualizing integrated molecular information with comutation plots. Bioinformatics 36, 4348-4349 (2020)) plot showing information on 101 MBC patients with tumor fraction >0.10. The top row shows the ER status used to train and evaluate the regression model. This was the metastatic ER status obtained by IHC for the majority of patients, and primary ER status was used when metastatic ER status was unavailable. ER low (1-10% ER staining) was considered ER positive. Second row, left upward triangle contains primary ER status, right downward triangle contains metastatic ER status. ER low biopsies are shown in light blue, unknown status in white. Third row: tumor fraction, fraction of cfDNA originating from tumor, calculated using ichorcNA. Fourth row, median probability of ER positivity calculated by Griffin over 1000 bootstrap iterations. FIG. 10D is a receiver operating characteristic (ROC) curve of a logistic regression model predicting ER+ and ER- subtypes. ROC curves, accuracy and AUC are shown for all patients and for patients grouped by tumor fraction (TFx) between 0.05-0.1 and 0.1 or greater. 95% CI was obtained by bootstrapping. For patients with multiple samples, the first sample with tumor fraction >0.05 was used. FIG. 10E graphically illustrates the performance of the model on samples from the three validation cohorts. For patients with multiple time points, the first sample was used. Figure 10F graphically illustrates subtype prediction in patients separated by clinical metastatic ER status and clinical primary tumor ER status. P values were calculated using Fisher's exact test (2-sided). Figure 10G illustrates the ROC curve for predicting ER loss among patients with primary ER-positive tumors. 95% CI was obtained by bootstrapping. Figure 10H illustrates the time series for two patients (MBC1413 and MBC1099) with different subtypes and multiple cfDNA samples across multiple biopsies. ER+ predicted probability (thick grey line) and tumor fraction (thin grey line) for all cfDNA samples passing thresholds of tumor fraction >0.05 and 0.1x coverage are shown. The decision boundaries for ER+ (>0.5) and ER- (<0.5) are shown as dotted lines. Time series in months from metastatic diagnosis to death for each patient are shown. For patient MBC_1413, a metastatic biopsy (pleural fluid) was obtained on the day the metastasis was diagnosed and showed ER- disease. However, approximately 7 months later, another metastatic biopsy (liver) showed weak ER+ staining (5%). A final biopsy (pleural fluid) was obtained at approximately 12 months and again showed ER- staining. Plasma for cfDNA was removed between the second and third metastatic biopsies, with one final removal after the third biopsy. For patient MBC_1099, two ER- biopsies were obtained at 0 months (bone) and 7 months (liver). cfDNA was removed after this time point, but another biopsy (liver) showed the presence of low levels of ER+ disease between the two cfDNA removals.

Figures 11A and 11B illustrate the workflow for characterizing advanced prostate cancer by corresponding tumors and liquid biopsies from PDX models. The top panel of Figure 11A illustrates blood and tissue samples obtained from 26 patient-derived xenograft (PDX) mouse models with tumors originating from metastatic castration-resistant prostate cancer (mCRPC) with AR-positive adenocarcinoma (ARPC), neuroendocrine prostate cancer (NEPC) and AR-low non-neuroendocrine prostate cancer (ARLPC) phenotypes. Cell-free DNA (cfDNA) was extracted from pooled plasma collected from 7-10 mice and whole genome sequencing (WGS) was performed. After bioinformatics subtraction of mouse reads, pure human circulating tumor DNA (ctDNA) reads remained. ATAC-Seq and CUT&RUN (targeting H3K27ac, H3K4me1, and H3K27me3) data were generated from PDX tissues. The middle panel of FIG. 11A illustrates two separate ctDNA features analyzed at transcription factor binding sites (TFBS) and open chromatin sites throughout the genome using Griffin (see Example 1 and Doebley et al. (2021). Griffin: Framework for clinical cancer subtyping from nucleosome profiling of cell-free DNA. MedRxiv 2021.08.31.21262867 and Methods). The bottom right panel of FIG. 11A shows phenotyping using a probability model that takes into account ctDNA tumor content and is informed by PDX features, applied to 159 samples in a cohort of three patients. 11B illustrates PDX phenotypes and sequencing of mouse plasma. Inclusion status based on final average depth after mouse read subtraction (less than 3x coverage was excluded unless AR coordinate amplification signal was reliably detected; lower dotted line). Phenotype status including 6 NEPC, 18 ARPC (2 excluded), and 2 ARLPC. Average depth of coverage before and after mouse subtraction (average coverage 20.5x; upper dotted line). Percentage of cfDNA samples containing human ctDNA after mouse read subtraction. Figures 11A and 11B illustrate the workflow for characterizing advanced prostate cancer by corresponding tumors and liquid biopsies from PDX models. The top panel of Figure 11A illustrates blood and tissue samples obtained from 26 patient-derived xenograft (PDX) mouse models with tumors originating from metastatic castration-resistant prostate cancer (mCRPC) with AR-positive adenocarcinoma (ARPC), neuroendocrine prostate cancer (NEPC) and AR-low non-neuroendocrine prostate cancer (ARLPC) phenotypes. Cell-free DNA (cfDNA) was extracted from pooled plasma collected from 7-10 mice and whole genome sequencing (WGS) was performed. After bioinformatics subtraction of mouse reads, pure human circulating tumor DNA (ctDNA) reads remained. ATAC-Seq and CUT&RUN (targeting H3K27ac, H3K4me1, and H3K27me3) data were generated from PDX tissues. The middle panel of FIG. 11A illustrates two separate ctDNA features analyzed at transcription factor binding sites (TFBS) and open chromatin sites throughout the genome using Griffin (see Example 1 and Doebley et al. (2021). Griffin: Framework for clinical cancer subtyping from nucleosome profiling of cell-free DNA. MedRxiv 2021.08.31.21262867 and Methods). The bottom right panel of FIG. 11A shows phenotyping using a probability model that takes into account ctDNA tumor content and is informed by PDX features, applied to 159 samples in a cohort of three patients. 11B illustrates PDX phenotypes and sequencing of mouse plasma. Inclusion status based on final average depth after mouse read subtraction (less than 3x coverage was excluded unless AR coordinate amplification signal was reliably detected; lower dotted line). Phenotype status including 6 NEPC, 18 ARPC (2 excluded), and 2 ARLPC. Average depth of coverage before and after mouse subtraction (average coverage 20.5x; upper dotted line). Percentage of cfDNA samples containing human ctDNA after mouse read subtraction. Figures 11A and 11B illustrate the workflow for characterizing advanced prostate cancer by corresponding tumors and liquid biopsies from PDX models. The top panel of Figure 11A illustrates blood and tissue samples obtained from 26 patient-derived xenograft (PDX) mouse models with tumors originating from metastatic castration-resistant prostate cancer (mCRPC) with AR-positive adenocarcinoma (ARPC), neuroendocrine prostate cancer (NEPC) and AR-low non-neuroendocrine prostate cancer (ARLPC) phenotypes. Cell-free DNA (cfDNA) was extracted from pooled plasma collected from 7-10 mice and whole genome sequencing (WGS) was performed. After bioinformatics subtraction of mouse reads, pure human circulating tumor DNA (ctDNA) reads remained. ATAC-Seq and CUT&RUN (targeting H3K27ac, H3K4me1, and H3K27me3) data were generated from PDX tissues. The middle panel of FIG. 11A illustrates two separate ctDNA features analyzed at transcription factor binding sites (TFBS) and open chromatin sites throughout the genome using Griffin (see Example 1 and Doebley et al. (2021). Griffin: Framework for clinical cancer subtyping from nucleosome profiling of cell-free DNA. MedRxiv 2021.08.31.21262867 and Methods). The bottom right panel of FIG. 11A shows phenotyping using a probability model that takes into account ctDNA tumor content and is informed by PDX features, applied to 159 samples in a cohort of three patients. 11B illustrates PDX phenotypes and sequencing of mouse plasma. Inclusion status based on final average depth after mouse read subtraction (less than 3x coverage was excluded unless AR coordinate amplification signal was reliably detected; lower dotted line). Phenotype status including 6 NEPC, 18 ARPC (2 excluded), and 2 ARLPC. Average depth of coverage before and after mouse subtraction (average coverage 20.5x; upper dotted line). Percentage of cfDNA samples containing human ctDNA after mouse read subtraction.

Figures 12A-12G illustrate that analysis of tumor histone modifications and ctDNA revealed nucleosomal patterns consistent with transcriptional regulation in CRPC phenotype-specific genes. Figure 12A illustrates H3K27ac peak signals across 10,000 AR binding sites (left) and ASCL1 binding sites (right) among ARLPC, ARPC, and NEPC PDX tumor phenotypes. Binding sites were selected from GTRD (Yevshin et al. (2019). GTRD: a database on gene transcription regulation-2019 update. Nucleic Acids Res 47, D100-D105) (Methods). Figures 12B and 12C graphically illustrate the combined coverage profile across 1000 AR binding sites (12B) and ASCL1 binding sites (12C) in ctDNA analyzed using Griffin. Coverage profile means (lines) and 95% confidence intervals (shading) using 1000 bootstraps are shown. Regions ±150 bp are indicated with vertical dotted lines and yellow shading. Figure 12D is a heatmap of _log2 fold changes in key up- and down-regulated genes between ARPC and NEPC established by RNA-Seq, grouped by the type of histone modification that dictates translation levels (left): group 1 shows genes where the predominant PTM marks are due to H3K27ac or H3K4me1 active marks at gene promoters or putative distal enhancers, and lack H3K27me3 heterochromatin marks within the transcribed region of the gene (gene body); group 2 features H3K27me3 repressive marks that extend within the transcribed region of the gene. The middle column shows the difference in peak intensity for each of the histone modifications assayed, separated by whether they appeared upstream, in the promoter or in the transcribed region of each gene. On the right, the log ₂ fold change in coefficient of variation (CV) of fragment sizes between ARPC and NEPC lines is shown for the TSS +/- 1 KB window and within the transcribed region of each gene. Figure 12E graphically illustrates a comparison of the log ₂ fold change in mean mRNA expression and mean coefficient of variation (CV) in the promoter regions of 47 phenotypic lineage marker genes (ARPC vs. NEPC). Figure 12F (top) presents an example of the predicted ctDNA coverage profile for group 1 genes with or without H3K27ac or H3K4me1 modifications that guide active and inactive transcription, respectively. Figure 12F (bottom) ±1000 bp surrounding the promoter regions of AR and ASCL1 in ARPC and NEPC. The coverage profile mean (line) and 95% confidence interval (shading) with 1000 bootstraps are shown. A decrease in coverage reflects increased nucleosome accessibility and therefore increased transcription. The dotted line and yellow shading highlight a confined window around the transcription start site (TSS) (TSS-230bp to +170bp). Figure 12G is an illustration of a predicted ctDNA coverage profile for group 2 genes with repression of transcription caused by H3K27me3 modification within the transcribed region of the gene. The neuronal gene UNC13A has increased nucleosome staging in ctDNA of ARPC samples compared to NEPC. Figures 12A-12G illustrate that analysis of tumor histone modifications and ctDNA revealed nucleosomal patterns consistent with transcriptional regulation in CRPC phenotype-specific genes. Figure 12A illustrates H3K27ac peak signals across 10,000 AR binding sites (left) and ASCL1 binding sites (right) among ARLPC, ARPC, and NEPC PDX tumor phenotypes. Binding sites were selected from GTRD (Yevshin et al. (2019). GTRD: a database on gene transcription regulation-2019 update. Nucleic Acids Res 47, D100-D105) (Methods). Figures 12B and 12C graphically illustrate the combined coverage profile across 1000 AR binding sites (12B) and ASCL1 binding sites (12C) in ctDNA analyzed using Griffin. Coverage profile means (lines) and 95% confidence intervals (shading) using 1000 bootstraps are shown. Regions ±150 bp are indicated with vertical dotted lines and yellow shading. Figure 12D is a heatmap of _log2 fold changes in key up- and down-regulated genes between ARPC and NEPC established by RNA-Seq, grouped by the type of histone modification that dictates translation levels (left): group 1 shows genes where the predominant PTM marks are due to H3K27ac or H3K4me1 active marks at gene promoters or putative distal enhancers, and lack H3K27me3 heterochromatin marks within the transcribed region of the gene (gene body); group 2 features H3K27me3 repressive marks that extend within the transcribed region of the gene. The middle column shows the difference in peak intensity for each of the histone modifications assayed, separated by whether they appeared upstream, in the promoter or in the transcribed region of each gene. On the right, the log ₂ fold change in coefficient of variation (CV) of fragment sizes between ARPC and NEPC lines is shown for the TSS +/- 1 KB window and within the transcribed region of each gene. Figure 12E graphically illustrates a comparison of the log ₂ fold change in mean mRNA expression and mean coefficient of variation (CV) in the promoter regions of 47 phenotypic lineage marker genes (ARPC vs. NEPC). Figure 12F (top) presents an example of the predicted ctDNA coverage profile for group 1 genes with or without H3K27ac or H3K4me1 modifications that guide active and inactive transcription, respectively. Figure 12F (bottom) ±1000 bp surrounding the promoter regions of AR and ASCL1 in ARPC and NEPC. The coverage profile mean (line) and 95% confidence interval (shading) with 1000 bootstraps are shown. A decrease in coverage reflects increased nucleosome accessibility and therefore increased transcription. The dotted line and yellow shading highlight a confined window around the transcription start site (TSS) (TSS-230bp to +170bp). Figure 12G is an illustration of a predicted ctDNA coverage profile for group 2 genes with repression of transcription caused by H3K27me3 modification within the transcribed region of the gene. The neuronal gene UNC13A has increased nucleosome staging in ctDNA of ARPC samples compared to NEPC. Figures 12A-12G illustrate that analysis of tumor histone modifications and ctDNA revealed nucleosomal patterns consistent with transcriptional regulation in CRPC phenotype-specific genes. Figure 12A illustrates H3K27ac peak signals across 10,000 AR binding sites (left) and ASCL1 binding sites (right) among ARLPC, ARPC, and NEPC PDX tumor phenotypes. Binding sites were selected from GTRD (Yevshin et al. (2019). GTRD: a database on gene transcription regulation-2019 update. Nucleic Acids Res 47, D100-D105) (Methods). Figures 12B and 12C graphically illustrate the combined coverage profile across 1000 AR binding sites (12B) and ASCL1 binding sites (12C) in ctDNA analyzed using Griffin. Coverage profile means (lines) and 95% confidence intervals (shading) using 1000 bootstraps are shown. Regions ±150 bp are indicated with vertical dotted lines and yellow shading. Figure 12D is a heatmap of _log2 fold changes in key up- and down-regulated genes between ARPC and NEPC established by RNA-Seq, grouped by the type of histone modification that dictates translation levels (left): group 1 shows genes where the predominant PTM marks are due to H3K27ac or H3K4me1 active marks at gene promoters or putative distal enhancers, and lack H3K27me3 heterochromatin marks within the transcribed region of the gene (gene body); group 2 features H3K27me3 repressive marks that extend within the transcribed region of the gene. The middle column shows the difference in peak intensity for each of the histone modifications assayed, separated by whether they appeared upstream, in the promoter or in the transcribed region of each gene. On the right, the log ₂ fold change in coefficient of variation (CV) of fragment sizes between ARPC and NEPC lines is shown for the TSS +/- 1 KB window and within the transcribed region of each gene. Figure 12E graphically illustrates a comparison of the log ₂ fold change in mean mRNA expression and mean coefficient of variation (CV) in the promoter regions of 47 phenotypic lineage marker genes (ARPC vs. NEPC). Figure 12F (top) presents an example of the predicted ctDNA coverage profile for group 1 genes with or without H3K27ac or H3K4me1 modifications that guide active and inactive transcription, respectively. Figure 12F (bottom) ±1000 bp surrounding the promoter regions of AR and ASCL1 in ARPC and NEPC. The coverage profile mean (line) and 95% confidence interval (shading) with 1000 bootstraps are shown. A decrease in coverage reflects increased nucleosome accessibility and therefore increased transcription. The dotted line and yellow shading highlight a confined window around the transcription start site (TSS) (TSS-230bp to +170bp). Figure 12G is an illustration of a predicted ctDNA coverage profile for group 2 genes with repression of transcription caused by H3K27me3 modification within the transcribed region of the gene. The neuronal gene UNC13A has increased nucleosome staging in ctDNA of ARPC samples compared to NEPC. Figures 12A-12G illustrate that analysis of tumor histone modifications and ctDNA revealed nucleosomal patterns consistent with transcriptional regulation in CRPC phenotype-specific genes. Figure 12A illustrates H3K27ac peak signals across 10,000 AR binding sites (left) and ASCL1 binding sites (right) among ARLPC, ARPC, and NEPC PDX tumor phenotypes. Binding sites were selected from GTRD (Yevshin et al. (2019). GTRD: a database on gene transcription regulation-2019 update. Nucleic Acids Res 47, D100-D105) (Methods). Figures 12B and 12C graphically illustrate the combined coverage profile across 1000 AR binding sites (12B) and ASCL1 binding sites (12C) in ctDNA analyzed using Griffin. Coverage profile means (lines) and 95% confidence intervals (shading) using 1000 bootstraps are shown. Regions ±150 bp are indicated with vertical dotted lines and yellow shading. Figure 12D is a heatmap of _log2 fold changes in key up- and down-regulated genes between ARPC and NEPC established by RNA-Seq, grouped by the type of histone modification that dictates translation levels (left): group 1 shows genes where the predominant PTM marks are due to H3K27ac or H3K4me1 active marks at gene promoters or putative distal enhancers, and lack H3K27me3 heterochromatin marks within the transcribed region of the gene (gene body); group 2 features H3K27me3 repressive marks that extend within the transcribed region of the gene. The middle column shows the difference in peak intensity for each of the histone modifications assayed, separated by whether they appeared upstream, in the promoter or in the transcribed region of each gene. On the right, the log ₂ fold change in coefficient of variation (CV) of fragment sizes between ARPC and NEPC lines is shown for the TSS +/- 1 KB window and within the transcribed region of each gene. Figure 12E graphically illustrates a comparison of the log ₂ fold change in mean mRNA expression and mean coefficient of variation (CV) in the promoter regions of 47 phenotypic lineage marker genes (ARPC vs. NEPC). Figure 12F (top) presents an example of the predicted ctDNA coverage profile for group 1 genes with or without H3K27ac or H3K4me1 modifications that guide active and inactive transcription, respectively. Figure 12F (bottom) ±1000 bp surrounding the promoter regions of AR and ASCL1 in ARPC and NEPC. The coverage profile mean (line) and 95% confidence interval (shading) with 1000 bootstraps are shown. A decrease in coverage reflects increased nucleosome accessibility and therefore increased transcription. The dotted line and yellow shading highlight a confined window around the transcription start site (TSS) (TSS-230bp to +170bp). Figure 12G is an illustration of a predicted ctDNA coverage profile for group 2 genes with repression of transcription caused by H3K27me3 modification within the transcribed region of the gene. The neuronal gene UNC13A has increased nucleosome staging in ctDNA of ARPC samples compared to NEPC. Figures 12A-12G illustrate that analysis of tumor histone modifications and ctDNA revealed nucleosomal patterns consistent with transcriptional regulation in CRPC phenotype-specific genes. Figure 12A illustrates H3K27ac peak signals across 10,000 AR binding sites (left) and ASCL1 binding sites (right) among ARLPC, ARPC, and NEPC PDX tumor phenotypes. Binding sites were selected from GTRD (Yevshin et al. (2019). GTRD: a database on gene transcription regulation-2019 update. Nucleic Acids Res 47, D100-D105) (Methods). Figures 12B and 12C graphically illustrate the combined coverage profile across 1000 AR binding sites (12B) and ASCL1 binding sites (12C) in ctDNA analyzed using Griffin. Coverage profile means (lines) and 95% confidence intervals (shading) using 1000 bootstraps are shown. Regions ±150 bp are indicated with vertical dotted lines and yellow shading. Figure 12D is a heatmap of _log2 fold changes in key up- and down-regulated genes between ARPC and NEPC established by RNA-Seq, grouped by the type of histone modification that dictates translation levels (left): group 1 shows genes where the predominant PTM marks are due to H3K27ac or H3K4me1 active marks at gene promoters or putative distal enhancers, and lack H3K27me3 heterochromatin marks within the transcribed region of the gene (gene body); group 2 features H3K27me3 repressive marks that extend within the transcribed region of the gene. The middle column shows the difference in peak intensity for each of the histone modifications assayed, separated by whether they appeared upstream, in the promoter or in the transcribed region of each gene. On the right, the log ₂ fold change in coefficient of variation (CV) of fragment sizes between ARPC and NEPC lines is shown for the TSS +/- 1 KB window and within the transcribed region of each gene. Figure 12E graphically illustrates a comparison of the log ₂ fold change in mean mRNA expression and mean coefficient of variation (CV) in the promoter regions of 47 phenotypic lineage marker genes (ARPC vs. NEPC). Figure 12F (top) presents an example of the predicted ctDNA coverage profile for group 1 genes with or without H3K27ac or H3K4me1 modifications that guide active and inactive transcription, respectively. Figure 12F (bottom) ±1000 bp surrounding the promoter regions of AR and ASCL1 in ARPC and NEPC. The coverage profile mean (line) and 95% confidence interval (shading) with 1000 bootstraps are shown. A decrease in coverage reflects increased nucleosome accessibility and therefore increased transcription. The dotted line and yellow shading highlight a confined window around the transcription start site (TSS) (TSS-230bp to +170bp). Figure 12G is an illustration of a predicted ctDNA coverage profile for group 2 genes with repression of transcription caused by H3K27me3 modification within the transcribed region of the gene. The neuronal gene UNC13A has increased nucleosome staging in ctDNA of ARPC samples compared to NEPC. Figures 12A-12G illustrate that analysis of tumor histone modifications and ctDNA revealed nucleosomal patterns consistent with transcriptional regulation in CRPC phenotype-specific genes. Figure 12A illustrates H3K27ac peak signals across 10,000 AR binding sites (left) and ASCL1 binding sites (right) among ARLPC, ARPC, and NEPC PDX tumor phenotypes. Binding sites were selected from GTRD (Yevshin et al. (2019). GTRD: a database on gene transcription regulation-2019 update. Nucleic Acids Res 47, D100-D105) (Methods). Figures 12B and 12C graphically illustrate the combined coverage profile across 1000 AR binding sites (12B) and ASCL1 binding sites (12C) in ctDNA analyzed using Griffin. Coverage profile means (lines) and 95% confidence intervals (shading) using 1000 bootstraps are shown. Regions ±150 bp are indicated with vertical dotted lines and yellow shading. Figure 12D is a heatmap of _log2 fold changes in key up- and down-regulated genes between ARPC and NEPC established by RNA-Seq, grouped by the type of histone modification that dictates translation levels (left): group 1 shows genes where the predominant PTM marks are due to H3K27ac or H3K4me1 active marks at gene promoters or putative distal enhancers, and lack H3K27me3 heterochromatin marks within the transcribed region of the gene (gene body); group 2 features H3K27me3 repressive marks that extend within the transcribed region of the gene. The middle column shows the difference in peak intensity for each of the histone modifications assayed, separated by whether they appeared upstream, in the promoter or in the transcribed region of each gene. On the right, the log ₂ fold change in coefficient of variation (CV) of fragment sizes between ARPC and NEPC lines is shown for the TSS +/- 1 KB window and within the transcribed region of each gene. Figure 12E graphically illustrates a comparison of the log ₂ fold change in mean mRNA expression and mean coefficient of variation (CV) in the promoter regions of 47 phenotypic lineage marker genes (ARPC vs. NEPC). Figure 12F (top) presents an example of the predicted ctDNA coverage profile for group 1 genes with or without H3K27ac or H3K4me1 modifications that guide active and inactive transcription, respectively. Figure 12F (bottom) ±1000 bp surrounding the promoter regions of AR and ASCL1 in ARPC and NEPC. The coverage profile mean (line) and 95% confidence interval (shading) with 1000 bootstraps are shown. A decrease in coverage reflects increased nucleosome accessibility and therefore increased transcription. The dotted line and yellow shading highlight a confined window around the transcription start site (TSS) (TSS-230bp to +170bp). Figure 12G is an illustration of a predicted ctDNA coverage profile for group 2 genes with repression of transcription caused by H3K27me3 modification within the transcribed region of the gene. The neuronal gene UNC13A has increased nucleosome staging in ctDNA of ARPC samples compared to NEPC. Figures 12A-12G illustrate that analysis of tumor histone modifications and ctDNA revealed nucleosomal patterns consistent with transcriptional regulation in CRPC phenotype-specific genes. Figure 12A illustrates H3K27ac peak signals across 10,000 AR binding sites (left) and ASCL1 binding sites (right) among ARLPC, ARPC, and NEPC PDX tumor phenotypes. Binding sites were selected from GTRD (Yevshin et al. (2019). GTRD: a database on gene transcription regulation-2019 update. Nucleic Acids Res 47, D100-D105) (Methods). Figures 12B and 12C graphically illustrate the combined coverage profile across 1000 AR binding sites (12B) and ASCL1 binding sites (12C) in ctDNA analyzed using Griffin. Coverage profile means (lines) and 95% confidence intervals (shading) using 1000 bootstraps are shown. Regions ±150 bp are indicated with vertical dotted lines and yellow shading. Figure 12D is a heatmap of _log2 fold changes in key up- and down-regulated genes between ARPC and NEPC established by RNA-Seq, grouped by the type of histone modification that dictates translation levels (left): group 1 shows genes where the predominant PTM marks are due to H3K27ac or H3K4me1 active marks at gene promoters or putative distal enhancers, and lack H3K27me3 heterochromatin marks within the transcribed region of the gene (gene body); group 2 features H3K27me3 repressive marks that extend within the transcribed region of the gene. The middle column shows the difference in peak intensity for each of the histone modifications assayed, separated by whether they appeared upstream, in the promoter or in the transcribed region of each gene. On the right, the log ₂ fold change in coefficient of variation (CV) of fragment sizes between ARPC and NEPC lines is shown for the TSS +/- 1 KB window and within the transcribed region of each gene. Figure 12E graphically illustrates a comparison of the log ₂ fold change in mean mRNA expression and mean coefficient of variation (CV) in the promoter regions of 47 phenotypic lineage marker genes (ARPC vs. NEPC). Figure 12F (top) presents an example of the predicted ctDNA coverage profile for group 1 genes with or without H3K27ac or H3K4me1 modifications that guide active and inactive transcription, respectively. Figure 12F (bottom) ±1000 bp surrounding the promoter regions of AR and ASCL1 in ARPC and NEPC. The coverage profile mean (line) and 95% confidence interval (shading) with 1000 bootstraps are shown. A decrease in coverage reflects increased nucleosome accessibility and therefore increased transcription. The dotted line and yellow shading highlight a confined window around the transcription start site (TSS) (TSS-230bp to +170bp). Figure 12G is an illustration of a predicted ctDNA coverage profile for group 2 genes with repression of transcription caused by H3K27me3 modification within the transcribed region of the gene. The neuronal gene UNC13A has increased nucleosome staging in ctDNA of ARPC samples compared to NEPC. Figures 12A-12G illustrate that analysis of tumor histone modifications and ctDNA revealed nucleosomal patterns consistent with transcriptional regulation in CRPC phenotype-specific genes. Figure 12A illustrates H3K27ac peak signals across 10,000 AR binding sites (left) and ASCL1 binding sites (right) among ARLPC, ARPC, and NEPC PDX tumor phenotypes. Binding sites were selected from GTRD (Yevshin et al. (2019). GTRD: a database on gene transcription regulation-2019 update. Nucleic Acids Res 47, D100-D105) (Methods). Figures 12B and 12C graphically illustrate the combined coverage profile across 1000 AR binding sites (12B) and ASCL1 binding sites (12C) in ctDNA analyzed using Griffin. Coverage profile means (lines) and 95% confidence intervals (shading) using 1000 bootstraps are shown. Regions ±150 bp are indicated with vertical dotted lines and yellow shading. Figure 12D is a heatmap of _log2 fold changes in key up- and down-regulated genes between ARPC and NEPC established by RNA-Seq, grouped by the type of histone modification that dictates translation levels (left): group 1 shows genes where the predominant PTM marks are due to H3K27ac or H3K4me1 active marks at gene promoters or putative distal enhancers, and lack H3K27me3 heterochromatin marks within the transcribed region of the gene (gene body); group 2 features H3K27me3 repressive marks that extend within the transcribed region of the gene. The middle column shows the difference in peak intensity for each of the histone modifications assayed, separated by whether they appeared upstream, in the promoter or in the transcribed region of each gene. On the right, the log ₂ fold change in coefficient of variation (CV) of fragment sizes between ARPC and NEPC lines is shown for the TSS +/- 1 KB window and within the transcribed region of each gene. Figure 12E graphically illustrates a comparison of the log ₂ fold change in mean mRNA expression and mean coefficient of variation (CV) in the promoter regions of 47 phenotypic lineage marker genes (ARPC vs. NEPC). Figure 12F (top) presents an example of the predicted ctDNA coverage profile for group 1 genes with or without H3K27ac or H3K4me1 modifications that guide active and inactive transcription, respectively. Figure 12F (bottom) ±1000 bp surrounding the promoter regions of AR and ASCL1 in ARPC and NEPC. The coverage profile mean (line) and 95% confidence interval (shading) with 1000 bootstraps are shown. A decrease in coverage reflects increased nucleosome accessibility and therefore increased transcription. The dotted line and yellow shading highlight a confined window around the transcription start site (TSS) (TSS-230bp to +170bp). Figure 12G is an illustration of a predicted ctDNA coverage profile for group 2 genes with repression of transcription caused by H3K27me3 modification within the transcribed region of the gene. The neuronal gene UNC13A has increased nucleosome staging in ctDNA of ARPC samples compared to NEPC.

FIG. 13 illustrates hierarchical clustering of normalized composite central mean coverage in TFBS by Griffin analysis of ctDNA for 107 TFs in LuCaP PDX lines with ARPC (n=16), NEPC (n=6), and ARLPC (n=2) phenotypes. This list of TFs was initially selected as having differential expression between ARPC and NEPC from LuCaP PDX RNA-Seq analysis. Heatmap colors indicate increased accessibility (low values; light) and decreased accessibility (high values; dark) in ctDNA. TFs with increased accessibility in NEPC samples ( _log2 fold change >0.05, Mann-Whitney U test p<0.05) are shown in red text; increased accessibility in ARPC (log2 fold change <-0.05, p<0.05) are shown in blue text. FIG. 13 illustrates hierarchical clustering of normalized composite central mean coverage in TFBS by Griffin analysis of ctDNA for 107 TFs in LuCaP PDX lines with ARPC (n=16), NEPC (n=6), and ARLPC (n=2) phenotypes. This list of TFs was initially selected as having differential expression between ARPC and NEPC from LuCaP PDX RNA-Seq analysis. Heatmap colors indicate increased accessibility (low values; light) and decreased accessibility (high values; dark) in ctDNA. TFs with increased accessibility in NEPC samples ( _log2 fold change >0.05, Mann-Whitney U test p<0.05) are shown in red text; increased accessibility in ARPC (log2 fold change <-0.05, p<0.05) are shown in blue text.

Figures 14A-14G illustrate a comprehensive assessment of genome-wide ctDNA signatures for CRPC phenotyping in PDX models. Figure 14A illustrates a volcano plot of log ₂ fold changes in ATAC-Seq peak intensities between five ARPC and five NEPC lines; dotted lines demarcate sites with q-values <0.05. Figures 14B and 14C graphically illustrate combined coverage profiles at ARPC PDX tumor-specific (14B) and NEPC PDX tumor-specific (14C) open chromatin sites analyzed by Griffin. Sites from (14A) were filtered for overlap with known TFBS in 338 factors from GTRD (Yevshin et al. (2019). Nucleic Acids Res 47, D100-D105). The coverage profile mean (line) and 95% confidence interval (shading) with 1000 bootstraps are shown. The region ±150 bp is indicated with a vertical dotted line and yellow shading. Figure 14D illustrates CA of ctDNA features demonstrating grouping between ARPC and NEPC phenotypes: (left panel) Combined central coverage of significant TFBS for 74 TFs with differential accessibility among 338 factors between ARPC and NEPC. (right panel) Fragment size variability (coefficient of variation) at H3K4me1 histone modification sites (n=9,750). Figure 14E graphically illustrates the performance of supervised machine learning (XGBoost) for classification of ARPC and NEPC PDX from ctDNA using different region types (all genes, TFBS, and open regions, methods). The area under the receiver operating characteristic curve (AUC) with 95% confidence intervals (100 replicates of stratified cross-validation) is shown for the performance of all feature types. FIG. 14F is an example of a composite coverage profile at ARPC-specific (left) and NEPC-specific (right) open chromatin sites identified in 14B-14C. Shown are simulated mixtures generated using ARPC (left) mixed with healthy donor (HD) and NEPC (right) mixed with HD for various tumor fractions. FIG. 14G graphically illustrates the performance of classifying mixed samples using stochastic mixture models. For each phenotype from the PDX lineage, five ctDNA mixtures were generated, each with various sequencing coverage and tumor fraction. In total, 125 mixtures were evaluated. The average AUC across the five mixtures for each configuration is shown. Figures 14A-14G illustrate a comprehensive assessment of genome-wide ctDNA signatures for CRPC phenotyping in PDX models. Figure 14A illustrates a volcano plot of log ₂ fold changes in ATAC-Seq peak intensities between five ARPC and five NEPC lines; dotted lines demarcate sites with q-values <0.05. Figures 14B and 14C graphically illustrate combined coverage profiles at ARPC PDX tumor-specific (14B) and NEPC PDX tumor-specific (14C) open chromatin sites analyzed by Griffin. Sites from (14A) were filtered for overlap with known TFBS in 338 factors from GTRD (Yevshin et al. (2019). Nucleic Acids Res 47, D100-D105). The coverage profile mean (line) and 95% confidence interval (shading) with 1000 bootstraps are shown. The region ±150 bp is indicated with a vertical dotted line and yellow shading. Figure 14D illustrates CA of ctDNA features demonstrating grouping between ARPC and NEPC phenotypes: (left panel) Combined central coverage of significant TFBS for 74 TFs with differential accessibility among 338 factors between ARPC and NEPC. (right panel) Fragment size variability (coefficient of variation) at H3K4me1 histone modification sites (n=9,750). Figure 14E graphically illustrates the performance of supervised machine learning (XGBoost) for classification of ARPC and NEPC PDX from ctDNA using different region types (all genes, TFBS, and open regions, methods). The area under the receiver operating characteristic curve (AUC) with 95% confidence intervals (100 replicates of stratified cross-validation) is shown for the performance of all feature types. FIG. 14F is an example of a composite coverage profile at ARPC-specific (left) and NEPC-specific (right) open chromatin sites identified in 14B-14C. Shown are simulated mixtures generated using ARPC (left) mixed with healthy donor (HD) and NEPC (right) mixed with HD for various tumor fractions. FIG. 14G graphically illustrates the performance of classifying mixed samples using stochastic mixture models. For each phenotype from the PDX lineage, five ctDNA mixtures were generated, each with various sequencing coverage and tumor fraction. In total, 125 mixtures were evaluated. The average AUC across the five mixtures for each configuration is shown. Figures 14A-14G illustrate a comprehensive assessment of genome-wide ctDNA signatures for CRPC phenotyping in PDX models. Figure 14A illustrates a volcano plot of log ₂ fold changes in ATAC-Seq peak intensities between five ARPC and five NEPC lines; dotted lines demarcate sites with q-values <0.05. Figures 14B and 14C graphically illustrate combined coverage profiles at ARPC PDX tumor-specific (14B) and NEPC PDX tumor-specific (14C) open chromatin sites analyzed by Griffin. Sites from (14A) were filtered for overlap with known TFBS in 338 factors from GTRD (Yevshin et al. (2019). Nucleic Acids Res 47, D100-D105). The coverage profile mean (line) and 95% confidence interval (shading) with 1000 bootstraps are shown. The region ±150 bp is indicated with a vertical dotted line and yellow shading. Figure 14D illustrates CA of ctDNA features demonstrating grouping between ARPC and NEPC phenotypes: (left panel) Combined central coverage of significant TFBS for 74 TFs with differential accessibility among 338 factors between ARPC and NEPC. (right panel) Fragment size variability (coefficient of variation) at H3K4me1 histone modification sites (n=9,750). Figure 14E graphically illustrates the performance of supervised machine learning (XGBoost) for classification of ARPC and NEPC PDX from ctDNA using different region types (all genes, TFBS, and open regions, methods). The area under the receiver operating characteristic curve (AUC) with 95% confidence intervals (100 replicates of stratified cross-validation) is shown for the performance of all feature types. FIG. 14F is an example of a composite coverage profile at ARPC-specific (left) and NEPC-specific (right) open chromatin sites identified in 14B-14C. Shown are simulated mixtures generated using ARPC (left) mixed with healthy donor (HD) and NEPC (right) mixed with HD for various tumor fractions. FIG. 14G graphically illustrates the performance of classifying mixed samples using stochastic mixture models. For each phenotype from the PDX lineage, five ctDNA mixtures were generated, each with various sequencing coverage and tumor fraction. In total, 125 mixtures were evaluated. The average AUC across the five mixtures for each configuration is shown. Figures 14A-14G illustrate a comprehensive assessment of genome-wide ctDNA signatures for CRPC phenotyping in PDX models. Figure 14A illustrates a volcano plot of log ₂ fold changes in ATAC-Seq peak intensities between five ARPC and five NEPC lines; dotted lines demarcate sites with q-values <0.05. Figures 14B and 14C graphically illustrate combined coverage profiles at ARPC PDX tumor-specific (14B) and NEPC PDX tumor-specific (14C) open chromatin sites analyzed by Griffin. Sites from (14A) were filtered for overlap with known TFBS in 338 factors from GTRD (Yevshin et al. (2019). Nucleic Acids Res 47, D100-D105). The coverage profile mean (line) and 95% confidence interval (shading) with 1000 bootstraps are shown. The region ±150 bp is indicated with a vertical dotted line and yellow shading. Figure 14D illustrates CA of ctDNA features demonstrating grouping between ARPC and NEPC phenotypes: (left panel) Combined central coverage of significant TFBS for 74 TFs with differential accessibility among 338 factors between ARPC and NEPC. (right panel) Fragment size variability (coefficient of variation) at H3K4me1 histone modification sites (n=9,750). Figure 14E graphically illustrates the performance of supervised machine learning (XGBoost) for classification of ARPC and NEPC PDX from ctDNA using different region types (all genes, TFBS, and open regions, methods). The area under the receiver operating characteristic curve (AUC) with 95% confidence intervals (100 replicates of stratified cross-validation) is shown for the performance of all feature types. FIG. 14F is an example of a composite coverage profile at ARPC-specific (left) and NEPC-specific (right) open chromatin sites identified in 14B-14C. Shown are simulated mixtures generated using ARPC (left) mixed with healthy donor (HD) and NEPC (right) mixed with HD for various tumor fractions. FIG. 14G graphically illustrates the performance of classifying mixed samples using stochastic mixture models. For each phenotype from the PDX lineage, five ctDNA mixtures were generated, each with various sequencing coverage and tumor fraction. In total, 125 mixtures were evaluated. The average AUC across the five mixtures for each configuration is shown. Figures 14A-14G illustrate a comprehensive assessment of genome-wide ctDNA signatures for CRPC phenotyping in PDX models. Figure 14A illustrates a volcano plot of log ₂ fold changes in ATAC-Seq peak intensities between five ARPC and five NEPC lines; dotted lines demarcate sites with q-values <0.05. Figures 14B and 14C graphically illustrate combined coverage profiles at ARPC PDX tumor-specific (14B) and NEPC PDX tumor-specific (14C) open chromatin sites analyzed by Griffin. Sites from (14A) were filtered for overlap with known TFBS in 338 factors from GTRD (Yevshin et al. (2019). Nucleic Acids Res 47, D100-D105). The coverage profile mean (line) and 95% confidence interval (shading) with 1000 bootstraps are shown. The region ±150 bp is indicated with a vertical dotted line and yellow shading. Figure 14D illustrates CA of ctDNA features demonstrating grouping between ARPC and NEPC phenotypes: (left panel) Combined central coverage of significant TFBS for 74 TFs with differential accessibility among 338 factors between ARPC and NEPC. (right panel) Fragment size variability (coefficient of variation) at H3K4me1 histone modification sites (n=9,750). Figure 14E graphically illustrates the performance of supervised machine learning (XGBoost) for classification of ARPC and NEPC PDX from ctDNA using different region types (all genes, TFBS, and open regions, methods). The area under the receiver operating characteristic curve (AUC) with 95% confidence intervals (100 replicates of stratified cross-validation) is shown for the performance of all feature types. FIG. 14F is an example of a composite coverage profile at ARPC-specific (left) and NEPC-specific (right) open chromatin sites identified in 14B-14C. Shown are simulated mixtures generated using ARPC (left) mixed with healthy donor (HD) and NEPC (right) mixed with HD for various tumor fractions. FIG. 14G graphically illustrates the performance of classifying mixed samples using stochastic mixture models. For each phenotype from the PDX lineage, five ctDNA mixtures were generated, each with various sequencing coverage and tumor fraction. In total, 125 mixtures were evaluated. The average AUC across the five mixtures for each configuration is shown. Figures 14A-14G illustrate a comprehensive assessment of genome-wide ctDNA signatures for CRPC phenotyping in PDX models. Figure 14A illustrates a volcano plot of log ₂ fold changes in ATAC-Seq peak intensities between five ARPC and five NEPC lines; dotted lines demarcate sites with q-values <0.05. Figures 14B and 14C graphically illustrate combined coverage profiles at ARPC PDX tumor-specific (14B) and NEPC PDX tumor-specific (14C) open chromatin sites analyzed by Griffin. Sites from (14A) were filtered for overlap with known TFBS in 338 factors from GTRD (Yevshin et al. (2019). Nucleic Acids Res 47, D100-D105). The coverage profile mean (line) and 95% confidence interval (shading) with 1000 bootstraps are shown. The region ±150 bp is indicated with a vertical dotted line and yellow shading. Figure 14D illustrates CA of ctDNA features demonstrating grouping between ARPC and NEPC phenotypes: (left panel) Combined central coverage of significant TFBS for 74 TFs with differential accessibility among 338 factors between ARPC and NEPC. (right panel) Fragment size variability (coefficient of variation) at H3K4me1 histone modification sites (n=9,750). Figure 14E graphically illustrates the performance of supervised machine learning (XGBoost) for classification of ARPC and NEPC PDX from ctDNA using different region types (all genes, TFBS, and open regions, methods). The area under the receiver operating characteristic curve (AUC) with 95% confidence intervals (100 replicates of stratified cross-validation) is shown for the performance of all feature types. FIG. 14F is an example of a composite coverage profile at ARPC-specific (left) and NEPC-specific (right) open chromatin sites identified in 14B-14C. Shown are simulated mixtures generated using ARPC (left) mixed with healthy donor (HD) and NEPC (right) mixed with HD for various tumor fractions. FIG. 14G graphically illustrates the performance of classifying mixed samples using stochastic mixture models. For each phenotype from the PDX lineage, five ctDNA mixtures were generated, each with various sequencing coverage and tumor fraction. In total, 125 mixtures were evaluated. The average AUC across the five mixtures for each configuration is shown.

Figures 15A-15C illustrate accurate classification of NEPC phenotype from plasma in three patient cohorts using a probability model informed by PDX ctDNA features. Figure 15A graphically illustrates the receiver operating characteristic (ROC) curves for 101 mCRPC patients (DFCI cohort I) with ultra low-pass WGS (ULP-WGS) data. The optimal performance of 90.4% sensitivity (for predicting NEPC) and 97.5% specificity (for predicting ARPC) corresponds to a prediction score cutoff of 0.3314, shown by the horizontal and vertical dotted lines, respectively. Figure 15B illustrates the prediction scores for 11 plasma samples from seven patients (DFCI cohort II) with both WGS and ULP-WGS data. A 0.3314 score cutoff threshold (dotted line) was used to classify NEPC and ARPC. Tumor fraction was estimated by ichorcna from WGS data. Patients were treated for adenocarcinoma (ARPC) or had high PSA levels. Figure 15C illustrates the prediction scores for 47 plasma samples with clinical phenotypes including mixed or ambiguous phenotypes (triangles), including 26 ARPC, 5 NEPC, and 16 double-negative prostate cancer (DNPC). Scores for WGS and ULP-WGS (0.1x) for the same ctDNA samples are shown. A cutoff threshold of 0.3314 (dotted line) was used to classify NEPC and ARPC. Tumor fraction was estimated by ichorcna on WGS data. Figures 15A-15C illustrate accurate classification of NEPC phenotype from plasma in three patient cohorts using a probability model informed by PDX ctDNA features. Figure 15A graphically illustrates the receiver operating characteristic (ROC) curves for 101 mCRPC patients (DFCI cohort I) with ultra low-pass WGS (ULP-WGS) data. The optimal performance of 90.4% sensitivity (for predicting NEPC) and 97.5% specificity (for predicting ARPC) corresponds to a prediction score cutoff of 0.3314, shown by the horizontal and vertical dotted lines, respectively. Figure 15B illustrates the prediction scores for 11 plasma samples from seven patients (DFCI cohort II) with both WGS and ULP-WGS data. A 0.3314 score cutoff threshold (dotted line) was used to classify NEPC and ARPC. Tumor fraction was estimated by ichorcna from WGS data. Patients were treated for adenocarcinoma (ARPC) or had high PSA levels. Figure 15C illustrates the prediction scores for 47 plasma samples with clinical phenotypes including mixed or ambiguous phenotypes (triangles), including 26 ARPC, 5 NEPC, and 16 double-negative prostate cancer (DNPC). Scores for WGS and ULP-WGS (0.1x) for the same ctDNA samples are shown. A cutoff threshold of 0.3314 (dotted line) was used to classify NEPC and ARPC. Tumor fraction was estimated by ichorcna on WGS data. Figures 15A-15C illustrate accurate classification of NEPC phenotype from plasma in three patient cohorts using a probability model informed by PDX ctDNA features. Figure 15A graphically illustrates the receiver operating characteristic (ROC) curves for 101 mCRPC patients (DFCI cohort I) with ultra low-pass WGS (ULP-WGS) data. The optimal performance of 90.4% sensitivity (for predicting NEPC) and 97.5% specificity (for predicting ARPC) corresponds to a prediction score cutoff of 0.3314, shown by the horizontal and vertical dotted lines, respectively. Figure 15B illustrates the prediction scores for 11 plasma samples from seven patients (DFCI cohort II) with both WGS and ULP-WGS data. A 0.3314 score cutoff threshold (dotted line) was used to classify NEPC and ARPC. Tumor fraction was estimated by ichorcna from WGS data. Patients were treated for adenocarcinoma (ARPC) or had high PSA levels. Figure 15C illustrates the prediction scores for 47 plasma samples with clinical phenotypes including mixed or ambiguous phenotypes (triangles), including 26 ARPC, 5 NEPC, and 16 double-negative prostate cancer (DNPC). Scores for WGS and ULP-WGS (0.1x) for the same ctDNA samples are shown. A cutoff threshold of 0.3314 (dotted line) was used to classify NEPC and ARPC. Tumor fraction was estimated by ichorcna on WGS data. Figures 15A-15C illustrate accurate classification of NEPC phenotype from plasma in three patient cohorts using a probability model informed by PDX ctDNA features. Figure 15A graphically illustrates the receiver operating characteristic (ROC) curves for 101 mCRPC patients (DFCI cohort I) with ultra low-pass WGS (ULP-WGS) data. The optimal performance of 90.4% sensitivity (for predicting NEPC) and 97.5% specificity (for predicting ARPC) corresponds to a prediction score cutoff of 0.3314, shown by the horizontal and vertical dotted lines, respectively. Figure 15B illustrates the prediction scores for 11 plasma samples from seven patients (DFCI cohort II) with both WGS and ULP-WGS data. A 0.3314 score cutoff threshold (dotted line) was used to classify NEPC and ARPC. Tumor fraction was estimated by ichorcna from WGS data. Patients were treated for adenocarcinoma (ARPC) or had high PSA levels. Figure 15C illustrates the prediction scores for 47 plasma samples with clinical phenotypes including mixed or ambiguous phenotypes (triangles), including 26 ARPC, 5 NEPC, and 16 double-negative prostate cancer (DNPC). Scores for WGS and ULP-WGS (0.1x) for the same ctDNA samples are shown. A cutoff threshold of 0.3314 (dotted line) was used to classify NEPC and ARPC. Tumor fraction was estimated by ichorcna on WGS data.

FIG. 16 is a schematic diagram of an integrated, non-invasive, targeted sequencing assay based on cfDNA for detection of genetic alterations and prediction of key tumor epigenetic signatures in SCLC.

17A and 17B illustrate detection of transcription factor (TF) expression in SCLC models using targeted sequencing of cfDNA. FIG. 17A is a schematic of the experimental workflow for proof of concept of a negative control sample ("healthy donor") and a positive control sample ("flank tumor" from a SCLC cell model). FIG. 17B graphically illustrates the aggregated coverage across TFBS in targeted sequencing data for a healthy donor (top row) and a flank tumor (bottom row). TFBS is expected to be located at position 0 on the x-axis. Data is color-coded by predicted TF expression. cfDNA from a healthy donor is expected to reflect expression of REST, but not ASCL1, NEUROD1, or POU2F3. Systematic differences in coverage distribution as a function of TF expression were evident in the SCLC model. 17A and 17B illustrate detection of transcription factor (TF) expression in SCLC models using targeted sequencing of cfDNA. FIG. 17A is a schematic of the experimental workflow for proof of concept of a negative control sample ("healthy donor") and a positive control sample ("flank tumor" from a SCLC cell model). FIG. 17B graphically illustrates the aggregated coverage across TFBS in targeted sequencing data for a healthy donor (top row) and a flank tumor (bottom row). TFBS is expected to be located at position 0 on the x-axis. Data is color-coded by predicted TF expression. cfDNA from a healthy donor is expected to reflect expression of REST, but not ASCL1, NEUROD1, or POU2F3. Systematic differences in coverage distribution as a function of TF expression were evident in the SCLC model. 17A and 17B illustrate detection of transcription factor (TF) expression in SCLC models using targeted sequencing of cfDNA. FIG. 17A is a schematic of the experimental workflow for proof of concept of a negative control sample ("healthy donor") and a positive control sample ("flank tumor" from a SCLC cell model). FIG. 17B graphically illustrates the aggregated coverage across TFBS in targeted sequencing data for a healthy donor (top row) and a flank tumor (bottom row). TFBS is expected to be located at position 0 on the x-axis. Data is color-coded by predicted TF expression. cfDNA from a healthy donor is expected to reflect expression of REST, but not ASCL1, NEUROD1, or POU2F3. Systematic differences in coverage distribution as a function of TF expression were evident in the SCLC model. 17A and 17B illustrate detection of transcription factor (TF) expression in SCLC models using targeted sequencing of cfDNA. FIG. 17A is a schematic of the experimental workflow for proof of concept of a negative control sample ("healthy donor") and a positive control sample ("flank tumor" from a SCLC cell model). FIG. 17B graphically illustrates the aggregated coverage across TFBS in targeted sequencing data for a healthy donor (top row) and a flank tumor (bottom row). TFBS is expected to be located at position 0 on the x-axis. Data is color-coded by predicted TF expression. cfDNA from a healthy donor is expected to reflect expression of REST, but not ASCL1, NEUROD1, or POU2F3. Systematic differences in coverage distribution as a function of TF expression were evident in the SCLC model. 17A and 17B illustrate detection of transcription factor (TF) expression in SCLC models using targeted sequencing of cfDNA. FIG. 17A is a schematic of the experimental workflow for proof of concept of a negative control sample ("healthy donor") and a positive control sample ("flank tumor" from a SCLC cell model). FIG. 17B graphically illustrates the aggregated coverage across TFBS in targeted sequencing data for a healthy donor (top row) and a flank tumor (bottom row). TFBS is expected to be located at position 0 on the x-axis. Data is color-coded by predicted TF expression. cfDNA from a healthy donor is expected to reflect expression of REST, but not ASCL1, NEUROD1, or POU2F3. Systematic differences in coverage distribution as a function of TF expression were evident in the SCLC model.

18A-18C illustrate the inference of transcription factor activity using TFBS coverage distributions from SCLC patient samples for which corresponding tumor gene expression data were available. FIG. 18A graphically illustrates the aggregated coverage across TFBS in targeted sequencing data for healthy donors (top row) and patients with SCLC for which corresponding tumor tissue with gene expression data was available (bottom row). Samples are color-coded by predicted TF expression. Systematic differences in coverage distribution as a function of predicted TF expression are again evident. FIG. 18B illustrates gene expression of significant genes in selected patient samples displayed as a heat map. Cells are color-coded by Z-score and inset is log2(TPM+1). FIG. 18C illustrates the peak and trough amplitudes calculated from the coverage distribution across TFBS in each patient sample displayed as a heat map. Amplitudes are indicated by color and are also indicated inset. The magnitude of trough depth corresponds to the gene expression of key TFs in these bona fide SCLC patient samples. 18A-18C illustrate the inference of transcription factor activity using TFBS coverage distributions from SCLC patient samples for which corresponding tumor gene expression data were available. FIG. 18A graphically illustrates the aggregated coverage across TFBS in targeted sequencing data for healthy donors (top row) and patients with SCLC for which corresponding tumor tissue with gene expression data was available (bottom row). Samples are color-coded by predicted TF expression. Systematic differences in coverage distribution as a function of predicted TF expression are again evident. FIG. 18B illustrates gene expression of significant genes in selected patient samples displayed as a heat map. Cells are color-coded by Z-score and inset is log2(TPM+1). FIG. 18C illustrates the peak and trough amplitudes calculated from the coverage distribution across TFBS in each patient sample displayed as a heat map. Amplitudes are indicated by color and are also indicated inset. The magnitude of trough depth corresponds to the gene expression of key TFs in these bona fide SCLC patient samples. 18A-18C illustrate the inference of transcription factor activity using TFBS coverage distributions from SCLC patient samples for which corresponding tumor gene expression data were available. FIG. 18A graphically illustrates the aggregated coverage across TFBS in targeted sequencing data for healthy donors (top row) and patients with SCLC for which corresponding tumor tissue with gene expression data was available (bottom row). Samples are color-coded by predicted TF expression. Systematic differences in coverage distribution as a function of predicted TF expression are again evident. FIG. 18B illustrates gene expression of significant genes in selected patient samples displayed as a heat map. Cells are color-coded by Z-score and inset is log2(TPM+1). FIG. 18C illustrates the peak and trough amplitudes calculated from the coverage distribution across TFBS in each patient sample displayed as a heat map. Amplitudes are indicated by color and are also indicated inset. The magnitude of trough depth corresponds to the gene expression of key TFs in these bona fide SCLC patient samples. 18A-18C illustrate the inference of transcription factor activity using TFBS coverage distributions from SCLC patient samples for which corresponding tumor gene expression data were available. FIG. 18A graphically illustrates the aggregated coverage across TFBS in targeted sequencing data for healthy donors (top row) and patients with SCLC for which corresponding tumor tissue with gene expression data was available (bottom row). Samples are color-coded by predicted TF expression. Systematic differences in coverage distribution as a function of predicted TF expression are again evident. FIG. 18B illustrates gene expression of significant genes in selected patient samples displayed as a heat map. Cells are color-coded by Z-score and inset is log2(TPM+1). FIG. 18C illustrates the peak and trough amplitudes calculated from the coverage distribution across TFBS in each patient sample displayed as a heat map. Amplitudes are indicated by color and are also indicated inset. The magnitude of trough depth corresponds to the gene expression of key TFs in these bona fide SCLC patient samples. 18A-18C illustrate the inference of transcription factor activity using TFBS coverage distributions from SCLC patient samples for which corresponding tumor gene expression data were available. FIG. 18A graphically illustrates the aggregated coverage across TFBS in targeted sequencing data for healthy donors (top row) and patients with SCLC for which corresponding tumor tissue with gene expression data was available (bottom row). Samples are color-coded by predicted TF expression. Systematic differences in coverage distribution as a function of predicted TF expression are again evident. FIG. 18B illustrates gene expression of significant genes in selected patient samples displayed as a heat map. Cells are color-coded by Z-score and inset is log2(TPM+1). FIG. 18C illustrates the peak and trough amplitudes calculated from the coverage distribution across TFBS in each patient sample displayed as a heat map. Amplitudes are indicated by color and are also indicated inset. The magnitude of trough depth corresponds to the gene expression of key TFs in these bona fide SCLC patient samples.

19 is a series of graphs illustrating the quantification of the peak and trough amplitudes of transcription factor binding sites in sample types. Distribution of peak and trough amplitudes in TFBS calculated from the aggregated coverage distribution according to the predicted ground truth of TF expression. Pdx samples labeled "non-SCLC" are NSCLC pdx models. Patient samples labeled "non-SCLC" are either samples from patients with NSCLC (n=11) or from patients without a diagnosed malignancy (n=4). Peak and trough amplitudes of ASCL1 sites are associated with both SCLC status and ASCL1 positivity, whereas peak and trough amplitudes of NEUROD1 and POU2F3 are associated with TF positivity only. 19 is a series of graphs illustrating the quantification of the peak and trough amplitudes of transcription factor binding sites in sample types. Distribution of peak and trough amplitudes in TFBS calculated from the aggregated coverage distribution according to the predicted ground truth of TF expression. Pdx samples labeled "non-SCLC" are NSCLC pdx models. Patient samples labeled "non-SCLC" are either samples from patients with NSCLC (n=11) or from patients without a diagnosed malignancy (n=4). Peak and trough amplitudes of ASCL1 sites are associated with both SCLC status and ASCL1 positivity, whereas peak and trough amplitudes of NEUROD1 and POU2F3 are associated with TF positivity only. 19 is a series of graphs illustrating the quantification of the peak and trough amplitudes of transcription factor binding sites in sample types. Distribution of peak and trough amplitudes in TFBS calculated from the aggregated coverage distribution according to the predicted ground truth of TF expression. Pdx samples labeled "non-SCLC" are NSCLC pdx models. Patient samples labeled "non-SCLC" are either samples from patients with NSCLC (n=11) or from patients without a diagnosed malignancy (n=4). Peak and trough amplitudes of ASCL1 sites are associated with both SCLC status and ASCL1 positivity, whereas peak and trough amplitudes of NEUROD1 and POU2F3 are associated with TF positivity only.

20A and 20B graphically illustrate gene expression inference using TSS coverage distributions in flank tumor positive control samples. FIG. 20A illustrates TSS coverage distributions from targeted sequencing of cfDNA grouped by quintiles of gene expression in SCLC flank tumor models (quintiles 1-5) and blood ("B", dark blue). Shown are 1,912 TSSs corresponding to 1,213 genes selected based on low expression in whole blood and correlation of TSS coverage distribution with gene expression. TSS coverage distributions vary systematically depending on the expression of the corresponding genes. FIG. 20B illustrates receiver operating characteristic curves for prediction of gene expression below threshold (shown for thresholds 0.1, 0.5, 1.0, and 2.0) as inferred from the coverage distributions of the above or corresponding TSSs. Estimates of gene expression were calculated from the TSS coverage profile as the magnitude of the difference between the average coverage depth at positions +130 and +145 relative to the TSS minus the average depth at positions -45, -30, and -15 (shown as dotted lines in 20A). The AUC of the ROC curves is shown in parentheses for each gene expression cutoff. In this preliminary analysis, restricted to particularly variable and therefore challenging genes, the TSS coverage distribution can be used to predict, with good test properties, whether a gene is expressed above or below a certain value.

21A-21C are a series of graphs illustrating the use of aggregated coverage profiles across a large, rationally selected subset of the TSS panel to predict SCLC or NSCLC status in lung cancer Pdx models and patient samples. The graphs provide examples of aggregated TSS coverage distribution across selected gene TSSs for upregulation in NSCLC (n=396) and SCLC (n=1045) for three different samples: one healthy donor (21A), one NSCLC Pdx model (21B), and one SCLC Pdx model (21C). As shown overlaid on the NSCLC PDX model, amplitude features were calculated from each coverage distribution curve as the difference between coverage at position -45 and coverage at position +120 for the TSS, facilitating comparisons within and between samples. 21A-21C are a series of graphs illustrating the use of aggregated coverage profiles across a large, rationally selected subset of the TSS panel to predict SCLC or NSCLC status in lung cancer Pdx models and patient samples. The graphs provide examples of aggregated TSS coverage distribution across selected gene TSSs for upregulation in NSCLC (n=396) and SCLC (n=1045) for three different samples: one healthy donor (21A), one NSCLC Pdx model (21B), and one SCLC Pdx model (21C). As shown overlaid on the NSCLC PDX model, amplitude features were calculated from each coverage distribution curve as the difference between coverage at position -45 and coverage at position +120 for the TSS, facilitating comparisons within and between samples. 21A-21C are a series of graphs illustrating the use of aggregated coverage profiles across a large, rationally selected subset of the TSS panel to predict SCLC or NSCLC status in lung cancer Pdx models and patient samples. The graphs provide examples of aggregated TSS coverage distribution across selected gene TSSs for upregulation in NSCLC (n=396) and SCLC (n=1045) for three different samples: one healthy donor (21A), one NSCLC Pdx model (21B), and one SCLC Pdx model (21C). As shown overlaid on the NSCLC PDX model, amplitude features were calculated from each coverage distribution curve as the difference between coverage at position -45 and coverage at position +120 for the TSS, facilitating comparisons within and between samples.

22A and 22B are a series of graphs illustrating the use of aggregated coverage profiles across a large, rationally selected subset of the TSS panel to predict SCLC or NSCLC status in lung cancer PDX models (22A) and patient samples (22B). Aggregated coverage of SCLC-specific gene TSS (y-axis, n=1045) versus NSCLC-specific gene TSS (x-axis, n=396) in plasma samples from lung cancer PDX samples (also shown as "benign" for reference for non-cancer control patients) or plasma samples from lung cancer patients. SCLC PDXs that have transdifferentiated from adenocarcinoma are identified by the thick red line. 22A and 22B are a series of graphs illustrating the use of aggregated coverage profiles across a large, rationally selected subset of the TSS panel to predict SCLC or NSCLC status in lung cancer PDX models (22A) and patient samples (22B). Aggregated coverage of SCLC-specific gene TSS (y-axis, n=1045) versus NSCLC-specific gene TSS (x-axis, n=396) in plasma samples from lung cancer PDX samples (also shown as "benign" for reference for non-cancer control patients) or plasma samples from lung cancer patients. SCLC PDXs that have transdifferentiated from adenocarcinoma are identified by the thick red line.

FIG. 23 is a flow chart illustrating a non-limiting example of an embodiment of a method of cell (e.g., cancer, e.g., prostate cancer) subtype prediction according to aspects of the present disclosure.

DETAILED DESCRIPTION The present disclosure is based on the inventors' development of an easy and sensitive method to assess chromatin architecture from cell-free DNA (cfDNA) and to provide an accurate signal to detect and distinguish cell and/or tissue phenotypes based on the determined chromatin architecture.

Cell-free DNA (cfDNA) is released from dying cells, including tumor cells, and can be isolated from peripheral blood for biological investigation. In oncology, as a subtype of cfDNA, circulating tumor DNA (ctDNA) is released from tumor cells into the blood. The presence of ctDNA represents an opportunity for a non-invasive "liquid biopsy" solution to address the tissue accessibility challenges mentioned above. Current research and clinical attempts are focused on detecting genetic mutations in selected cancer genes from ctDNA, and have demonstrated potential clinical utility. Sequencing analysis of ctDNA to detect genomic alterations has also helped classify several subsets of tumors based on genetic differences. However, investigating tumor phenotypes from ctDNA remains challenging and remains a nascent research area.

We and others have noted that cfDNA is protected from degradation by nucleosomes and other DNA-binding proteins in the bloodstream, resulting in a coverage pattern that reflects the genome organization in the cell of origin. Genomic organization includes patterns of chromatin accessibility and transcriptional regulation, which in turn drive the differential phenotype of the cell of origin. Thus, cfDNA may provide a non-invasive means to identify tumor subtypes by analysis of tumor phenotypes beyond traditional analysis of genotypes with DNA alterations.

Although there is an intriguing possibility of performing genomic, chromatin accessibility, transcriptomic, and transcriptional regulation "multi-omics" profiling from analysis of ctDNA samples, there have been no robust tools previously available for predicting these multi-omics profiles from ctDNA, particularly for non-invasive applications to identify genomic and phenotypic signature changes during disease progression and resistance to treatment in metastatic cancer. For example, previous attempts to demonstrate nucleosome occupancy at the transcription start sites (TSS) of actively transcribed genes can be used to predict the presence of transcription for individual genes (see Ulz P, et al. Inferring expressed genes by whole-genome sequencing of plasma DNA. Nat Genet. 2016; 48 (10): 1273-1278, incorporated herein by reference in its entirety). However, this approach requires large (>75%) tumor fractions or areas of somatic copy number amplification to make predictions with any reliable accuracy. In another study, the same concept was used to assess nucleosome occupancy of transcription factor binding sites (TFBS) to predict transcription factor (TF) activity (see Ulz P, et al. Inference of transcription factor binding from cell-free DNA enables tumor subtype prediction and early detection. Nat Commun. 2019; 10 (1): 4666, incorporated herein by reference in its entirety). This study demonstrated that TFs exhibited differential signals between adenocarcinoma and NEPC ctDNA samples. However, this approach did not account for or incorporate features such as local sequence bias (e.g., GC content), somatic CNAs, ctDNA fraction, tumor fraction, and cfDNA fragment size, which are enriched for shorter fragments in cancer patients compared to healthy donors. These factors were observed to significantly affect the results, thus obscuring the actual signal from the data. Therefore, it is unlikely that such existing methods can be sensitive and robust enough to work with data from ultra-low-pass whole genome sequencing data (ULP-WGS) (0.1x), which would make them a more cost-effective and therefore accessible sequencing option in clinical settings.

The inventors have addressed the shortcomings of the art to create an easy, robust, and sensitive approach to detect and distinguish cell phenotypes. As described in more detail below, the approach is based in part on a core method called "Griffin" for investigating nucleosome protection and chromatin accessibility by quantifying cfDNA fragments around accessible sites. Unlike previous methods, Griffin takes into account fragment length-based GC correction, implementing a crucial approach to remove signal-obscuring GC bias, especially prevalent in ULP-WGS applications (e.g., WGS with coverage as low as 0.1x). This novel fragment size-aware GC bias correction approach implemented in Griffin helps maximize signal-to-noise and optimize analysis of sequence data such as ULP-WGS for cfDNA. As an initial proof of concept described in Example 1, the inventors applied the Griffin approach to cfDNA samples. Griffin achieved excellent performance in detecting tumor cfDNA in early cancer patients (AUC=0.96). The approach was then applied to samples from subjects with metastatic breast cancer (MBC) with different phenotypes (i.e., ER+ and ER-), demonstrating that accurate determination of nucleosome occupancy allows for differentiation of the two phenotypes. Specifically, 254 samples from 139 patients were analyzed, and ER subtypes were predicted with high performance (AUC=0.89), leading to valuable insights into tumor heterogeneity. This demonstration that the use of cfDNA can be used to predict hormone subtypes has immediate application in clinical diagnostics to identify subtypes and potentially switch to minimally invasive and cost-effective modalities, thereby facilitating appropriate treatment. In this particular embodiment, the Griffin approach can be used to detect breast cancer during early disease, when the efficacy and impact of treatment may be greatest.

The Griffin method offers additional advantages. Griffin is flexible to analyze any region throughout the genome that may provide information on differential chromatin accessibility between cell/tissue/cancer phenotype environments. For example, important transcriptional factors that differentiate between tumor subtypes can be predicted using Griffin by analyzing the binding sites of these transcription factors. Furthermore, Griffin can be applied to various input data generated by different assay methods, including ATAC-seq, ChIP-seq, transcription factor profiling data, cut-and-run, etc., to investigate chromatin architecture and accessibility. Moreover, significantly different from existing technologies, Griffin is able to:
・Gene expression prediction (equivalent to transcriptomics)
Transcriptional regulation, e.g., transcription factor activity (gene regulation; reguromics)
Transcription factor activity during targeted therapy (gene regulation; reguromics)
・Chromatin accessibility (epigenetics)
Chromatin modifications, e.g., H3K27 acetylation (epigenetics)
- Hematopoietic and immune cell profiling (immunology)
By enabling the analysis of multiple "omics" such as those above, numerous hypotheses can be addressed.

The Griffin approach is adaptable to existing ctDNA sequencing techniques, thus enabling scalability, adaptability, and accessibility, even from ULP-WGS data that is highly susceptible to bias and signal obscuration. Key applications of the approach include tumor (subtype) classification, identification of mixed histology/phenotype, detection of potential subtype switching (transdifferentiation) during treatment in "real-time," and prediction of biomarkers that may signal treatment resistance (e.g., ARv7 splice variants).

As described in more detail below (e.g., in Example 2), we utilized circulating tumor DNA (ctDNA) to investigate tumor phenotypes by identifying nucleosome placement patterns associated with transcriptional regulation. We performed genome-wide sequencing of ctDNA in mouse plasma from 24 patient-derived xenograft models of androgen receptor-activated prostate cancer (ARPC) and neuroendocrine prostate cancer (NEPC). Nucleosome patterns associated with transcriptional activity were reflected in ctDNA at genes, promoters, histone modifications, transcription factor binding, and regions of accessible chromatin. We identified the activity of key transcription regulators from ctDNA, including AR, ASCL1, HOXB13, HNF4G, and NR3C1, that were associated with prostate cancer phenotypes. We subsequently designed a predictive model to distinguish NEPC from ARPC with 97-100% sensitivity and 85-100% specificity for 159 plasma samples across three clinical cohorts. These results highlight the utility of combining ctDNA with the Griffin workflow to investigate molecular phenotypes and advance diagnostics in precision medicine.

In accordance with the foregoing, in one aspect, the disclosure provides a computer-implemented method for augmenting sequence read data from a cell-free DNA sample to predict cell type. In this context, the phrase "cell type prediction" is used in a general sense to refer to predicting the identity of a cell of origin (i.e., a cell that contributes DNA in a cfDNA sample) or a trait thereof. For example, the trait may be a distinguishable phenotype compared to cells having the same or similar developmental lineage, including a developmental lineage with a transformation event (i.e., for cancer cells). Alternatively, the trait may be a distinguishable developmental lineage compared to a separate developmental lineage. As described in more detail below, the method encompasses predicting or identifying different cell lineages, different tissue types, different tissue subtypes, different cancer types, different cancer subtypes (i.e., subtypes of the same cancer type), and the like. All that is required is that cell types, as a broad definition, are distinguishable by unique nucleosome occupancy and/or chromatin accessibility profiles.

The method is:
receiving sequence read data by a computing system, the sequence read data including a plurality of fragment reads, each fragment read having a fragment length and a GC content indicating a percentage of bases in the fragment read that are G or C;
determining, by a computing system, a GC bias value for each fragment read based on the fragment length and GC content of the fragment read;
generating, by a computing system, a genome coverage distribution adjusted for GC bias using the sequence read data and the GC bias value; and predicting, by the computing system, a cell type based on the genome coverage distribution.

FIG. 1 is a flow chart illustrating a non-limiting example of an embodiment of a method for predicting cell type according to various aspects of the present disclosure. Method 100 includes using the GRIFFIN technique described elsewhere herein to enable extraction of meaningful features from short cancer DNA nucleic acid sequences obtained from sequencing cell-free DNA fragments in a sample. Method 100 can be used for a variety of different types of cell type prediction, including but not limited to tissue type prediction, cell type prediction, cancer type prediction, and cancer subtype prediction.

From a start block, method 100 proceeds to subroutine block 102, where genomic regions of interest are determined and filtered to identify sites that provide information about cell type. Any suitable technique for determining and filtering sites that provide information about cell type may be used, and different techniques may be used for different cancer types, different molecular subtypes of a cancer type, different tissues, different cell types, and different assay types. A non-limiting example of an embodiment of a suitable procedure for determining and filtering sites that provide information about cell type is illustrated in FIG. 2 and described in further detail below.

In subroutine block 104, a GC occurrence matrix is determined for a combination of fragment length and GC content. For a particular sequencing technology, fragments with a certain amount of G and C bases ("GC content") are over-represented in the sequence read data. This bias is not constant, as fragments of different sizes have different GC biases. Because sequence read data from cell-free DNA fragments typically contain many short fragments of different lengths, establishing a GC occurrence matrix that specifies the expected GC content ratios for a variety of different fragment lengths allows the sequence read data to be appropriately corrected for GC bias and to obtain a meaningful signal from sequence read data that would otherwise be too noisy. A non-limiting example of a technique for determining a GC occurrence matrix is illustrated in FIG. 3 and described in more detail below.

It will be appreciated that the actions described with respect to subroutine blocks 102 and 104 can be performed on the reference genomic data before obtaining the sample or analyzing the sequence data.

At block 106, sequence read data is received. In some embodiments, the sequence read data represents sequence reads generated for a sample obtained from a subject. In some embodiments, the sequence read data may be obtained from an archive or other previously obtained sample.

In subroutine block 108, the GC occurrence matrix is used to determine GC bias values for the sequence read data. Any suitable technique can be used in subroutine block 108, including but not limited to the non-limiting examples illustrated in FIG. 4 and described in more detail below.

In subroutine block 110, the GC bias values are used to generate a genome coverage distribution of sequence read data for sites that provide information about the cell type. Again, any suitable technique can be used in subroutine block 110, including but not limited to the non-limiting examples illustrated in FIG. 5 and described in more detail below.

At block 112, features are extracted from the genome coverage distribution. Any features suitable for use in a classifier model can be extracted, depending on the type of classifier model used, the assay that generated the sequence reads, and/or the cell type (e.g., cancer type, cancer subtype, tissue, or cell type) to be detected. As a non-limiting example, for estrogen receptor (ER) subtyping for breast cancer, three features can be extracted: mean coverage, center coverage, and amplitude.

The average coverage can be extracted by determining the average coverage within a window around the informative site. The window around the informative site for determining the average coverage can be of any suitable size, including but not limited to a range of 1800-2200 bp (+/- 900 bp to +/- 1100 bp). One non-limiting example of a suitable size of window for determining the average coverage is 2000 bp (+/- 1000 bp).

The central coverage can be extracted by determining the average coverage within a smaller window around the informative site. The window around the informative site for determining the central coverage can be of any suitable size, including but not limited to a range of 40-80 bp (+/- 20 bp to +/- 40 bp). One non-limiting example of a suitable size of window for determining the average coverage is 60 bp (+/- 30 bp).

The amplitude can be extracted by trimming the genome coverage distribution to a region containing a given number of peaks (e.g., a region of +/- 960 bp containing 10 peaks), performing a fast Fourier transform, and taking the magnitude of the frequency based on the given number of peaks (e.g., the 10th frequency for a region containing 10 peaks).

At block 114, the features are provided as input to a classifier model to predict the cell subtype. Any suitable classifier model may be used. In one non-limiting example embodiment, the classifier model may be a logistic regression model.

Once the cancer subtype is predicted by the classifier model, then method 100 proceeds to an end block and ends. Of course, in some embodiments, once the cancer subtype has been determined, further actions can be taken, including, but not limited to, diagnosing the appropriate cancer, identifying a change or switch in cancer subtype, recommending a new course of treatment, modifying an existing course of treatment, or any other suitable action.

One consideration and challenge with analyzing patient-derived plasma is the presence of cfDNA released by hematopoietic cells, which results in a small ctDNA fraction (i.e., tumor fraction). Furthermore, the small patient cohorts for which information on tumor phenotype is available make supervised machine learning approaches suboptimal. Therefore, an unsupervised probabilistic model was developed to estimate the proportion of cell types contributing to individual plasma samples. One advantage of this model is that the ctDNA tumor fraction in patients is explicitly modeled.

Inputs to the model include signals generated from patient-derived xenografts (PDXs). PDXs provide an ideal resource for investigating ctDNA characteristics by comparison with matched tumors, developing new analytical tools, and validating both genetic and phenotypic features. Using ctDNA fraction estimates and these input PDX signals, the model is applied to a statistical mixture model approach to estimate mixture weighting parameters that represent the proportions of cell types. As discussed below in Example 2 and illustrated in Figures 14-15, the mixture weighting parameters can be used as a predictive score to classify cell types, such as ARPC vs. NEPC. Other cell types, such as phenotypes and subtypes, can also be modeled and predicted using this framework.

Figure 2 is a flow chart illustrating a non-limiting example of an embodiment of a procedure for determining and filtering genomic regions of interest to identify sites that provide information about a cell type, according to various aspects of the present disclosure. In some embodiments, the cell types of interest for determining and filtering sites that provide information about a cell type are different cancer types, different cancer subtypes, different tissue types, or different cell types.

From a start block, the procedure 200 proceeds to block 202, where a list of potentially informative sites in the cell type of interest is selected. Sites can be selected using available data, including but not limited to public research databases and repositories, published scientific data, and sequencing data. These data may be derived from assays including but not limited to Assay for Transposase-Accessible Chromatin (ATACs-eq), Micrococcal Nuclease (MNase-seq), DNase Hypersensitive Sites, Chromatin Immunoprecipitation (ChIP-seq), and sequencing techniques for on-target cleavage and release using nucleases (cut-and-run). From these data, any suitable comparison is used to select sites that distinguish cell types (e.g., tissue types, cell types, cancer types, or cancer subtypes), including, but not limited to, statistical hypothesis testing using two-group Mann-Whitney U (also called Wilcoxon rank sum) test or Student's t-test and multigroup Kruskal-Wallis test or analysis of variance (ANOVA). Additional filtering can be performed using fold changes between groups.

In optional block 204, the average mappability score (a metric that represents the uniqueness of the genomic sequence) within a fixed-size window around each potentially informative site is determined, and in optional block 206, sites with an average mappability score below a predefined threshold are discarded. Retaining only mappable sites limits the analysis to sites that are likely to be accurately represented in the sequence read data. Mappability can be determined based on reference data, such as a mappability score track from the UCSC genome browser. In some embodiments, the actions of optional block 204 and optional block 206 may not be performed.

In block 208, remaining sites that provide information for determining cell type are identified. Any suitable technique can be used. For example, for ER subtyping of breast cancer, Cancer Genome Atlas (TCGA) ATAC seq data can be used to identify sites with differential ATAC signal between ER-positive and ER-negative TCGA samples. Any suitable technique can be used to identify these sites. In some embodiments, a Mann-Whitney U test can be used for each site, followed by a false discovery rate (FDR) correction using the Benjamini-Hochberg procedure, which can hold all sites at an adjusted p-value (i.e., q-value) of less than 0.05. In some embodiments, the ATAC seq read counts around each site can be provided as input to the DESeq2 software, which can then identify differential sites and generate adjusted fold changes and FDR-corrected p-values for each site.

In some embodiments, sites can be further refined by examining fold changes and retaining all sites with a log2 fold change of >0.5 for the subtype of interest relative to other subtypes. For ER subtyping of breast cancer, ER positive and ER negative sites can be separated into those shared with hematopoietic cells and those not shared with hematopoietic cells using separate datasets of hematopoietic ChIP-seq peaks to generate a total of four lists of subtype-specific informative sites.

Procedure 200 then proceeds to the end block and ends.

Figure 3 is a flow chart illustrating a non-limiting example of an embodiment of a procedure for determining a GC occurrence matrix for a genome according to various aspects of the present disclosure. The technique described in Figure 3 differs from previous techniques such as those described in Benjamini & Speed, 2012 and implemented in DeepTools (Ramirez, Duendar, Diehl, Gruening, & Manke, 2014) at least because the previous techniques did not correct for fragments of different lengths and have never been shown to work for cell-free DNA sequencing data. In procedure 300, separate GC bias curves are determined for each of the different fragment lengths.

First, all mappable regions of the genome are examined (block 302). Then, for each fragment length (where a for loop is defined as between a for loop start block 304 and a for loop end block 310), the number of times each GC content is observed within fragments of that fragment length within the mappable region is counted to determine the GC frequency for the genome (block 306), and the GC frequency for that fragment length is stored in a GC frequency matrix (block 308).

After the for loop, procedure 300 proceeds to the end block and ends.

In some embodiments, various fragment lengths between the short and long length thresholds are analyzed in procedure 300. In some embodiments, the short length threshold may be in the range of 10-20 bp and the long length threshold may be in the range of 450-550 bp. In one specific non-limiting embodiment, the short length threshold may be 15 bp and the long length threshold may be 500 bp. A for loop may be run for each fragment length between the short and long length thresholds.

Figure 4 is a flow chart illustrating a non-limiting example embodiment of a procedure for determining GC bias values for sequence read data using a GC frequency matrix in accordance with various aspects of the present disclosure.

In block 402, the number of observed reads and GC content of each fragment length are counted to determine a GC count for the sequence read data.

In block 404, the GC counts are divided by the values in the GC frequency matrix to determine the GC bias for each fragment length.

At block 406, the average GC bias for each fragment length is normalized to determine an approximate GC bias value. In some embodiments, the average GC bias may be normalized to 1, resulting in an approximate GC bias value for all possible combinations of fragment size and GC content.

At block 408, the approximate GC bias values are smoothed to determine a GC bias value. In some embodiments, for each fragment size, all GC bias values for similarly sized fragments (as a non-limiting example, for a 165 bp fragment, fragments sized from 155 bp to 175 bp may be considered) may be determined. The GC bias values for similarly sized fragments may be sorted by GC content, and kernel smoothing may be performed by taking the median of the nearest neighbors to determine the GC bias value.

Procedure 400 then proceeds to the end block and ends.

Figure 5 is a flow chart illustrating a non-limiting example of an embodiment of a procedure for generating genome coverage distributions of sequence read data for cell type-specific informative sites using GC bias values according to various aspects of the present disclosure.

Procedure 500 proceeds from a start block to block 502, where fragment midpoints within a window around each site that provides cell type-specific information are determined.

At block 504, a weight is assigned to each fragment based on the appropriate GC bias value for the fragment length and the GC content (i.e., the GC bias value and GC content for the fragment length determined, for example, by procedure 400 in subroutine block 108). The weight is then based on the appropriate GC bias value. In some embodiments, the weight can be the inverse of the GC bias value (1/GC bias value). For example, if in a given sample, 165 bp fragments with a GC content of 60% have a GC bias of 2.5 (a large proportion compared to 165 bp fragments with other GC contents), then these fragments are assigned a weight of 1/2.5 = 0.4.

In block 506, the weights are used to determine a GC corrected midpoint profile.

At block 508, positions that overlap with the excluded regions are excluded. The excluded regions may be determined using any suitable technique. In some embodiments, the excluded regions may be derived from one or more lists of excluded regions. The list of excluded regions may include, but is not limited to, the Encoding Integrated GRCh38 Exclusion List, centromeres, gaps in the human genome assembly, correction patches, alternative haplotypes, regions of zero mappability, and regions of unusually high coverage (e.g., 10 standard deviations above the mean).

In block 510, the GC-corrected midpoint profiles for all sites are averaged to determine the average profile.

At block 512, the average profile is smoothed to generate a smoothed average profile. Any suitable technique for smoothing may be used. For example, in some embodiments, the average profile may be smoothed using a Savitzky-Golay filter with a window length of 165 bp and a third order polynomial.

In block 514, the smoothed average profile is normalized by dividing it by the average of the surrounding coverage. In some embodiments, a range of 9,000-11,000 bp (+/- 4,500 bp to +/- 5,500 bp), e.g., 10,000 bp (+/- 5,000 bp), is taken into account for normalization. This allows samples with different depths of sequencing coverage to be compared.

The normalized average profile can be used as the resulting genome coverage distribution.

Procedure 500 then advances to and end block and ends.

6 is a block diagram illustrating aspects of an exemplary computing device suitable for use as a computing device of the present disclosure. The techniques described above, including but not limited to those described in method 100, may be implemented in whole or in part in one or more computing systems including one or more computing devices, such as computing device 600, communicatively coupled to each other.

The exemplary computing device 600 illustrates various elements common to many different types of computing devices, including, but not limited to, desktop computing devices, laptop computing devices, server computing devices, mobile computing devices, and computing devices that are part of a cloud computing system. Although FIG. 6 is described with reference to a computing device implemented as a device on a network, the following description is applicable to servers, personal computers, mobile phones, smartphones, tablet computers, embedded computing devices, and other devices that can be used to perform some of the embodiments of the present disclosure. Some embodiments of a computing device may be implemented in or include an application specific integrated circuit (ASIC), a field programmable gate array (FPGA), or other customized device. Moreover, those skilled in the art will appreciate that the computing device 600 may be any one of any number of currently available or yet to be developed devices.

In its most basic configuration, the computing device 600 includes at least one processor 602 and a system memory 610 connected by a communication bus 608. Depending on the exact configuration and type of device, the system memory 610 may be volatile or non-volatile memory, such as read-only memory ("ROM"), random access memory ("RAM"), EEPROM, flash memory or similar memory technology. Those skilled in the art will appreciate that the system memory 610 typically stores data and/or program modules that are immediately accessible to and/or currently being operated on by the processor 602. In this regard, the processor 602 may act as the computer center of the computing device 600 by supporting the execution of instructions.

As further illustrated in FIG. 6, the computing device 600 may include a network interface 606 that includes one or more components for communicating with other devices over a network. In embodiments of the present disclosure, basic services utilizing the network interface 606 may be accessed to effect communication using a common network protocol. The network interface 606 may also include a wireless network interface configured to communicate via one or more wireless communication protocols, such as, for example, Wi-Fi, 2G, 3G, LTE, WiMAX, Bluetooth, Bluetooth low energy, etc. As will be appreciated by those skilled in the art, the network interface 606 illustrated in FIG. 6 may be one or more of the wireless interfaces or physical communication interfaces described and illustrated above with respect to particular components of the computing device 600.

In the exemplary embodiment shown in FIG. 6, the computing device 600 also includes a storage medium 604. However, the service may be accessed using a computing device that does not include a means for persisting data on a local storage medium. Thus, the storage medium 604 shown in FIG. 6 is represented by a dashed line to indicate that the storage medium 604 is optional. In any event, the storage medium 604 may be volatile or non-volatile, removable or non-removable, and implemented using any technology capable of storing information, such as, but not limited to, a hard drive, solid state drive, CD ROM, DVD, or other disk storage, magnetic cassette, magnetic tape, magnetic disk storage, etc.

Suitable implementations of a computing device including a processor 602, a system memory 610, a communication bus 608, a storage medium 604, and a network interface 606 are known and commercially available. For ease of illustration and as not important to an understanding of the claimed subject matter, FIG. 6 does not show some components that are typical of many computing devices. In this regard, the computing device 600 may include input devices such as, for example, a keyboard, a keypad, a mouse, a microphone, a touch input device, a touch screen, a tablet, and the like. Such input devices may be coupled to the computing device 600 by wired or wireless connections, including RF, infrared, serial, parallel, Bluetooth, Bluetooth low energy, USB, or other suitable connection protocols using wireless or physical connections. Similarly, the computing device 600 may also include output devices such as a display, a speaker, a printer, and the like. These devices are well known in the art and therefore will not be illustrated or further described herein.

As shown, the computer-implemented method implementing the Griffin workflow is highly adaptable to different types of input data reflecting chromatin architecture (e.g., nucleosome occupancy and chromatin accessibility). The method can be applied to various environments of analysis depending on the source and characteristics of the original cell or tissue analyzed. Thus, in another aspect, the present disclosure provides a method for determining a chromatin accessibility profile for a cell of interest from a sample containing cell-free DNA derived from the cell of interest. The method applies the Griffin data optimization workflow described in more detail above to determine a chromatin accessibility profile for the cell of interest. The method is flexible and accepts input data obtained from a variety of sequencing and capture protocols. The method includes:
obtaining sequence read data from the cell-free DNA;
receiving sequence read data by a computing system, the sequence read data including a plurality of fragment reads, each fragment read having a fragment length and a GC content indicating a percentage of bases in the fragment read that are G or C;
determining, by a computing system, a GC bias value for each fragment read based on the fragment length and GC content of the fragment read;
The method includes generating, by a computing system, a genomic coverage distribution adjusted for GC bias using the sequence read data and the GC bias value; and determining a chromatin accessibility profile from the genomic coverage distribution.

The method may further include determining a cell phenotype of interest based on the chromatin occupancy profile. For example, determining the cell phenotype may include determining the tissue type of the cell of origin, determining whether the cell is transformed (e.g., cancerous or malignant), determining a cancer type or subtype, determining a malignant aggressive phenotype, and/or determining a drug responsive phenotype. The term malignant aggressive phenotype refers to the relative aggressiveness of transformed (e.g., cancer) cells in terms of their rate of reproduction, migration, drug responsiveness, etc. The phenotype may be qualitative or may be one that can be assessed by various metrics to allow for quantitative comparison. The term "drug responsive phenotype" refers to the relative responsiveness (i.e., susceptibility or resistance) of cancer cells to cancer therapy. The metrics may be quantitative or qualitative. These determinations may be made based on sequence data optimized by the Griffin workflow using various classifiers described in more detail above. Elements of the Griffin workflow and computer-implemented methods are described in more detail above and are incorporated into the present aspects without limitation. An exemplary, non-limiting implementation of the Griffin workflow and associated classifier for determining subtypes of cancer cells having distinct phenotypes is presented in the Examples.

As shown herein, the Griffin workflow augments data from various sequencing and capture platforms to provide nucleosome accessibility profiles that can provide highly accurate insights into the nature of the cells that contribute to the ctDNA present in a biological sample. These insights allow for the detection and characterization of the cells that contribute to the ctDNA, including enabling the ability to detect certain types of cells and/or distinguish between different subtypes of cells. Thus, in certain aspects, the present disclosure also provides a method for determining or identifying the cell type of a cell of interest from a sample that contains cell-free DNA derived from the cell of interest. The method of this aspect comprises:
Obtaining sequence read data generated from a sample comprising cell-free DNA;
The method includes: performing a computer-implemented method as described in more detail above (and incorporated in this aspect in all its embodiments); and determining or identifying the cell type of the cell of interest based on the prediction provided by the computing system. The determining step can be performed based on the data augmented by the Griffin workflow by any of a number of suitable classifiers. As described above, the determining step can include determining the cell phenotype, for example, the tissue type, cancer type, cancer subtype, aggressive phenotype of a malignant tumor, drug responsiveness phenotype, or expression (or expression level) of the gene of interest.

In yet another aspect, the disclosure provides a method for detecting the presence of cancer cells in a subject, the method comprising:
Obtaining sequence read data generated from a sample comprising cell-free DNA obtained from a subject;
performing a computer-implemented method as described in more detail above (and incorporated in this aspect in all its embodiments); and determining the presence of cancer cells in the subject based on the prediction provided by the computing system.

In some embodiments, the method is performed multiple times. Thus, the method can be a method of monitoring for the presence and/or identity of cancer in a subject. The cancer cell(s) detected in the subject at each performance of the method can be further characterized. For example, the method can be used to monitor the cell(s) over time to determine the cancer subtype or phenotype of the detected cancer cell(s) based on a prediction provided by a computing system. In some embodiments, the method further includes detecting changes in the phenotype of the detected cancer cell(s) over time. For example, as described in more detail below, certain cancer types may progress from one subtype to another during the course of the disease. Cancer cells may develop and essentially switch between characterized subtypes. These changes can be associated with changes in malignancy and/or responsiveness to various treatments, all of which can be detected given the demonstrated sensitivity of the Griffin workflow. Monitoring and documenting such changes over time provides information about the requirements for modifying treatment to optimize outcomes. As a non-limiting example, non-small cell lung cancer (NSCLC) can be monitored for transdifferentiation into small cell lung cancer (SCLC). Alternatively, SCLC subtypes can be monitored for transdifferentiation into distinct subtypes. In some embodiments, the method can be performed beginning before or during treatment for the cancer. Thus, the cancer can be monitored for responsiveness to treatment or for phenotypic changes during the course of treatment. These characteristics can provide information for any appropriate adjustments to the treatment regimen. In some embodiments, the method includes performing a treatment or treatment change based on the status of the monitored cancer cells as determined by the method.

In another aspect, the present disclosure provides a method for determining a cancer subtype of a target cancer cell from a sample comprising cell-free DNA derived from the target cancer cell, the method comprising:
Obtaining sequence read data generated from a sample comprising cell-free DNA;
performing a computer-implemented method as described in more detail above (and incorporated in this aspect in all its embodiments); and determining the cell type of the target cancer cells based on the predicted cancer subtype provided by the computing system.

The sample may be a biological sample from a subject, e.g., a subject having or suspected of having cancer. Exemplary biological samples are described in more detail below. In some embodiments, the method includes obtaining a biological sample from a subject and/or generating sequence read data from the sample according to standard techniques appropriate for a desired sequencing platform and/or targeted capture technology.

As described in more detail below, the Griffin platform has been used to successfully distinguish between important cancer subtypes for a variety of different unrelated cancers, demonstrating broad applicability across cancer types in general. Thus, in some embodiments, the cancer is characterized as metastatic breast cancer. In some further embodiments, the determining step includes determining the status of the breast cancer as ER+ or ER-, which refers to expression of the estrogen receptor (ER) and whether the cancer cells respond to exposure to estrogen hormone. Because ER+ breast cancer can be addressed by administering endocrine therapy, this status may be crucial to provide information regarding the appropriate course of treatment. In other embodiments, the determining step includes determining the status of the breast cancer as PR+ or PR-, which refers to expression of the progesterone receptor (PR) and whether the cancer cells respond to exposure to progesterone hormone. Similarly, because PR+ breast cancer can also be addressed by administering appropriate hormone therapy drugs such as tamoxifen and aromatase inhibitors, this status may be crucial to provide information regarding the appropriate course of treatment. In yet another embodiment, the determining step includes determining the status of the breast cancer as HER2+ or HER2-, which refers to the expression of human epidermal growth factor receptor 2 (HER2). HER2+ breast cancer cells tend to grow faster and are more likely to spread, for example, to lymph nodes, resulting in a poor prognosis. PR+ breast cancer can also be addressed by administering an appropriate Her2-targeted therapy, such as trastuzumab or pertuzumab, so this status can be crucial to inform the appropriate course of treatment. Of course, it will be understood that the present disclosure also encompasses distinguishing embodiments that determine the expression status of multiple informative markers. For example, the method can include determining whether the cancer is ER+ or ER-; whether the cancer is PR+ or PR-; and/or whether the cancer is HER2+ or HER2-, in any combination. For example, in one embodiment, the method includes determining whether the cancer is ER+ or ER-, whether the cancer is PR+ or PR-, and whether the cancer is HER2+ or HER2-. Patients with triple-negative breast cancer (i.e., ER-, PR-, HER-) can receive neoadjuvant chemotherapy, such as carboplatin and paclitaxel, in combination with immunotherapy, such as pembrolizumab and atezolizumab.

In another embodiment, the cancer is characterized as metastatic prostate cancer. In a further embodiment, determining the prostate cancer subtype deals with determining whether the cancer expresses various markers characteristic of distinct subtypes. For example, in one embodiment, the step of cancer subtype includes determining whether the prostate cancer is AR+ (ARPC) or AR-, which refers to the state of expression of the androgen receptor. Alternatively, the step of cancer subtype includes determining whether the prostate cancer is AR+ (ARPC) or AR(low). Prostate cancers that are AR+ are often treated with androgen receptor signaling inhibitors (ARSIs), which suppress androgen receptor activity in cells. In another embodiment, the step of cancer subtype includes determining whether the prostate cancer has a neuroendocrine prostate cancer (NEPC) phenotypic signature. NEPC cells lack AR activity and have a transcriptional programming regulatory profile, including epigenetic modifications, that are distinct from CRPC cells, thus resulting in a distinct phenotype that requires alternative therapeutic intervention. In another embodiment, the step of cancer subtyping includes determining whether the prostate cancer is amphicrine, which refers to having both exocrine and neuroendocrine traits in the same cell. As demonstrated in Example 2 below, the Griffin workflow can be utilized to accurately distinguish these cell types from input sequence reads generated from ctDNA. Of course, it will be understood that the present disclosure also encompasses differentiation embodiments that determine the status of multiple features to precisely determine the prostate cancer subtype in a subject. For example, determining the cancer subtype includes determining two, three, four, or all of any combination of the following: whether the cancer is AR+ (ARPC) or AR-, whether the cancer is AR-low or ARPC, whether the cancer has a neuroendocrine prostate cancer (NEPC) phenotypic signature, whether the cancer is AR-low or NEPC, whether the cancer is amphicrine, ARPC, or NEPC.

In another embodiment, the cancer is characterized as metastatic lung cancer. In a further embodiment, determining the subtype of the lung cancer includes determining whether the cancer is small cell lung cancer (SCLC) or non-small cell lung cancer (NSCLC). If the lung cancer is NSCLC, in a further embodiment, the method further includes determining whether the NSCLC is adenocarcinoma or squamous cell carcinoma.

As described above, the input sequence read data may be generated from a variety of platforms and using a variety of techniques, including whole genome analysis. In Example 3, the inventors established whole genome analysis, but this is not necessary. Instead, the inventors designed and implemented a panel of genomic targets deemed relevant for scientific research (e.g., subtyping lung cancer cells). Thus, in some embodiments, lung cancer is further subtyped using sequence read data generated from a panel of genomic targets. In some embodiments, the panel of genomic targets includes transcription factor binding sites (TFBS) of one or more transcription factors associated with the specified subtype that is the subject of the analysis, e.g., SCLC. For example, with respect to subtyping SCLC, the one or more associated transcription factors include one or more of ASLC, NEUROD1, POU2F3, REST, etc. In such embodiments, the method includes a step of determining the nucleosome occupancy of the TFBS using any suitable technique (e.g., cut and run, etc.). TFBSs can be identified by ChIP-seq data or similar techniques known in the art. If a candidate TFBS is proximal to the transcription start site (TSS) of a gene associated with lung cancer, or a subtype of lung cancer of interest in subtyping, it can be retained in the panel. In this regard, the term proximal can mean that the TFBS is in a vicinity that has a large functional impact on the initiation of transcription at the TSS. In some cases, a functional impact or association can be established if the TSS is the closest TSS to the TFBS. In other embodiments, the panel of genomic targets includes the transcription start site (TSS) of one or more markers associated with lung cancer (or a particular subtype of lung cancer of interest). In such embodiments, the method includes determining the nucleosome occupancy of the TSS by known techniques.

The biological sample described herein may be any sample obtained from a subject that may have cell-free DNA. Illustrative non-limiting examples encompassed by the present disclosure include samples that are blood, plasma, or serum, which are particularly useful for assessing cfDNA and ctDNA from a subject. In any embodiment of the foregoing aspects of detecting or assessing cancer in a subject, the method may further include obtaining a biological sample from the subject. Additionally, for a subject determined to have cancer or cancer subtype at any time, the method may further include prescribing an appropriate treatment or appropriately aggressively treating the subject according to accepted practices in the medical field for the determined cancer based on the determination of the cancer type or subtype.

In any of the aspects described herein, the described methods can be performed multiple times to provide multiple assessments. This can be useful to provide a method for monitoring the presence of a cell type or subtype or its evolution from a source. For example, the method can be performed from sequence read data obtained from a biological sample from a subject prior to and/or at or after the subject is first diagnosed with cancer.

As noted above, the Griffin workflow is flexible and is not limited to a particular set of genomic regions of interest or a particular type of sequence data for generating coverage profiles. Exemplary non-limiting techniques for generating sequence read data include whole genome sequencing (e.g., at a depth between 0.05× coverage and 100× coverage) and chromatin accessibility assays. In some embodiments, sequence read data is generated by or regions of interest are identified using techniques such as ATAC-seq, ChIP-seq, DNase sensitivity assays, etc., known in the art. In some embodiments, sequence data is generated by cut-and-run. See, e.g., WO 2019/060907, which is incorporated herein by reference in its entirety. For example, the cut-and-run assay can incorporate the use of one or more affinity reagents (e.g., antibodies or antibody fragments) targeting H3K27ac, H3K4me1, and/or H3K27ac post-translational modifications. In some embodiments, the method includes generating sequence read data deterministically, e.g., using any of the exemplary techniques described herein or other suitable techniques known in the art.

As described above, for lung cancer subtyping, sequence read data can be generated from a panel of genomic targets. It will be understood that this targeted panel approach is applicable to other cancer types beyond lung cancer subtyping. Thus, in some embodiments of any of the above methods, the sequence read data can include sequence read data generated from a panel of genomic targets. The panel of genomic targets can be designed and assembled according to the approach described in Example 3 for lung cancer (see also FIG. 16). For example, the panel can include TFBSs of one or more transcription factors associated with the cancer type of interest. Transcription factors associated with the cancer type of interest can be readily identified in the art. TFBSs for a specified transcription factor can be determined by standard assays that establish binding sites within the genome, such as, for example, ChIP-seq data. Additionally, candidate TFBSs can be further retained based on an assessment of association with or proximity to the transcription start site (TSS) of a gene with a transcription level (on, off, high, low, etc.) associated with the relevant cancer or cancer subtype. In some embodiments, the panel of genomic targets includes the transcription start sites (TSSs) of one or more markers associated with the cancer type of interest. The panel can be constructed using any combination of TFBSs and/or TSSs. Once established, directed sequencing reads are generated from the targets. In some embodiments, the nucleosome occupancy and/or TSSs of the TFBSs are determined. The sequence read data is input into the computer-implemented Griffin method described above to facilitate appropriate subtyping or other analysis.

As disclosed in Example 2, analysis of fragment size variation using the coefficient of variation to assess standard deviation yielded a metric with a surprisingly strong correlation with gene expression/activity. This analysis can be integrated into the implementation of the Griffin workflow or can be independent of it. Thus, in another aspect, the present disclosure provides a computer-implemented method for augmenting sequence read data from cell-free DNA samples to predict cell type. The method includes:
receiving, by a computing system, sequence read data, the sequence read data including a plurality of fragment reads, each fragment read having a fragment length;
determining, by a computing system, a variability in fragment sizes for at least one gene associated with a cell type; and predicting, by the computing system, the cell type based on the variability in fragment sizes for the at least one gene.

FIG. 23 is a flow chart illustrating a non-limiting example of an embodiment for enhancing sequence read data from cell-free DNA samples for improved cell type prediction according to various aspects of the present disclosure.

In block 702, sequence read data is received by a computing system, where the sequence read data includes a plurality of fragment reads, each fragment read having a fragment length.

At block 704, the computing system determines the fragment size variability for at least one gene associated with a cell type. In some embodiments, the locations of genes whose mRNA expression and transcriptional activity are known to be associated with a given cell type, such as the 47 genes illustrated in FIG. 12D that are known to be associated with prostate cancer, can be used.

In some embodiments, the coefficient of variation of fragment sizes for fragments at positions associated with one or more genes can be determined and used as a measure of fragment size variability. The coefficient of variation (CV) has been found to be particularly useful for distinguishing cell types based on fragment size variability when analyzing fragments in genes associated with cell types. In particular, CV has been found to be less sensitive to the depth of sequencing coverage than other techniques (e.g., measures of entropy, etc.).

At block 706, the computing system predicts the cell type based on the fragment size variability for at least one gene. In some embodiments, a feature can be generated based on the fragment size variability and provided as an input to a classifier model to determine whether the feature is representative of a given cell type. In one non-limiting example, a ratio of fragment size variability in a first cell type to fragment size variability in a second cell type can be used as a feature. The classifier model can be used to determine whether a calculated feature for a given sample is closer to a feature of the first cell type or a feature of the second cell type. Any suitable classifier model can be used, including but not limited to logistic regression models, artificial neural networks, decision trees, support vector machines, and Bayesian networks.

A non-limiting example of an embodiment of the use of method 700 is described in Example 2, in which analysis of fragment size variability is used to distinguish between different types of androgen receptor pathway-activated prostate cancer (ARPC) and different types of neuroendocrine prostate cancer (NEPC) prostate cancer cell types.

Additional definitions Unless otherwise specifically defined herein, all terms used herein have the same meaning as those skilled in the art of the present invention.For definitions and technical terms, practitioners are particularly directed to Sambrook J., et al. (eds)., Molecular Cloning: A Laboratory Manual, 3rd ed., Cold Spring Harbor Press, Plainsview, New York (2001); and Ausubel, FM, et al. (eds)., Current Protocols in Molecular Biology, John Wiley & Sons, New York (2010).

Although use of the term "or" in the claims is used to mean "and/or" unless expressly stated to refer to alternatives only or that the alternatives are mutually exclusive, this disclosure supports the definition referring to alternatives only and "and/or."

In accordance with long-standing patent law, the words "a" and "an," when used in conjunction with the word "comprising" in a claim or specification, refer to one or more, unless specifically stated otherwise.

Unless the context clearly requires otherwise, throughout the description and claims, words such as "comprise," "comprising," and the like, are to be construed in an inclusive sense, as opposed to an exclusive or exhaustive sense; that is, to indicate the meaning "including but not limited to." Words in the singular or plural also encompass the plural and singular, respectively. Additionally, the words "herein," "above," and "below," and words of similar import, when used in this application, shall refer to this application as a whole and not to any particular portion of this application. The word "about" refers to a number that falls within a range of minor variations above and below the stated reference number. For example, "about" may refer to a number that falls within a range of 10%, 9%, 8%, 7%, 6%, 5%, 4%, 3%, 2%, or 1% above or below the indicated reference number.

The terms "subject," "individual," and "patient" are used interchangeably herein to refer to a mammal being evaluated for and/or receiving treatment. In certain embodiments, the mammal is a human. The terms "subject," "individual," and "patient" include, but are not limited to, individuals with cancer. Although the subject can be a human, the terms also include other mammals, particularly mammals useful as experimental models for human disease, such as mice, rats, dogs, non-human primates, and the like.

The term "treating" and grammatical variations thereof may refer to any indication of successful treatment or reversal or prevention of a disease or condition (e.g., cancer, an infectious disease, or an autoimmune disease), including any objective or subjective parameter, such as relief; remission; a decrease in symptoms or making the disease state more tolerable to the patient; a slowing of the rate of degeneration or decline; or making the end point of degeneration less severe.

The treatment or amelioration of symptoms may be based on objective or subjective parameters, including the results of a physician's examination. Thus, the term "treating" includes administering a compound or agent of the present disclosure to prevent or delay, alleviate, improve clinical outcome, reduce the presence of symptoms, improve quality of life, prolong disease-free status, stabilize, extend survival, arrest or inhibit the development of associated symptoms or conditions, or any combination thereof, of a disease or condition (e.g., cancer). The term "therapeutic effect" refers to the reduction, elimination, or prevention of a disease or condition, a symptom of a disease or condition, or a side effect of a disease or condition in a subject.

As used herein, the term "nucleic acid" or "polynucleic acid" refers to a polymer of nucleotide monomeric units or "residues", typically DNA or RNA. The nucleotide monomeric subunits, or residues, of a nucleic acid each contain a nitrogenous base (i.e., nucleobase), a five-carbon sugar, and a phosphate group. The identity of each residue is generally indicated herein with reference to the identity of the nucleobase (or nitrogenous base) structure of each residue. Standard nucleobases include adenine (A), guanine (G), thymine (T), uracil (U) (in place of thymine (T) residues in RNA), and cytosine (C). However, the nucleic acids of the present disclosure may include any modified nucleobases, nucleobase analogs, and/or non-standard nucleobases known in the art.

Disclosed are materials, compositions, and components that can be used for, can be used in conjunction with, can be used in preparation for, or are products of the disclosed methods and compositions. When combinations, subsets, interactions, groups, etc. of these materials are disclosed, it is understood that each of the various individual and collective combinations is specifically contemplated, even though specific reference to any and all single combinations and permutations of these compounds may not be expressly disclosed. This concept applies to all aspects of the disclosure, including but not limited to the steps in the described methods. Thus, specific elements of any of the foregoing embodiments can be combined or substituted with elements in other embodiments. For example, where there are various additional steps that can be performed, it is understood that each of these additional steps can be performed with any particular method step or combination of method steps of the disclosed methods, and that such combinations or subsets of combinations should be considered specifically contemplated and disclosed. Additionally, it is understood that the embodiments described herein can be performed using any suitable materials, such as those described elsewhere herein or known in the art.

Publications cited herein and the subject matter for which they are cited are hereby specifically incorporated by reference in their entirety.

The following examples are provided so as to provide those of ordinary skill in the art with a complete disclosure and description of how to make and use the invention, and are not intended to limit the scope of what the inventors regard as their invention, nor are they intended to represent that the experiments below are all or the only experiments performed.

Example 1
This example describes a study that provides proof of concept that sequence analysis applied to an embodiment of the Griffin workflow disclosed herein enhances sequence signal with sufficient power and specificity to allow for breast cancer subtyping. Elements of this work are also described in Doebley, A.-L., et al. (2021). Griffin: Framework for clinical cancer subtyping from nucleosome profiling of cell-free DNA. MedRxiv 2021.08.31.21262867, which is incorporated herein by reference in its entirety.

Introduction Accurate cancer diagnosis and subtyping is crucial to guide clinical care and precision oncology. Current methods for determining tumor subtype require tissue biopsies, which are often difficult to obtain from patients with metastatic cancer. Thus, treatment options at the time of diagnosis of recurrent or metastatic cancer are often informed by the clinical diagnosis of the primary tumor. However, during metastatic progression and in the setting of treatment resistance, molecular changes in tumors can emerge. Moreover, molecular changes are difficult to investigate, as repeat biopsies are problematic and not routine for solid tumor practice.

Cell-free DNA (cfDNA) is DNA released into the circulation by cells during apoptosis and necrosis. In patients with cancer, a portion of this cfDNA is released from tumor cells and is referred to as circulating tumor DNA (ctDNA). Analysis of ctDNA can address the challenges of tissue accessibility and has demonstrated great potential for clinical utility. Many of the current research and clinical efforts have been focused on detecting genetic alterations in ctDNA. Shallow coverage sequencing of cfDNA, including ultra-low-pass whole genome sequencing (ULP-WGS, 0.1x), provides a cost-effective and scalable solution for estimating tumor fraction (the fraction of cfDNA originating from the tumor) from analysis of genomic copy number alterations. Sequencing analysis of genomic alterations from ctDNA has helped to distinguish molecular subsets of tumors. However, these genomic alterations, including somatic mutations, do not always fully explain treatment failure or identify therapeutic targets, thereby illustrating a major limitation of precision cancer medicine.

Tumor subtypes often feature distinct transcriptional regulation that can change during resistance to treatment, resulting in distinct clinical tumor phenotypes. For example, prostate and lung cancers can undergo transdifferentiation from adenocarcinoma to small cell neuroendocrine phenotype. For metastatic breast cancer (MBC), treatment is guided based on clinical subtype, often determined by the expression of estrogen receptor (ER), progesterone receptor (PR), and human epidermal growth factor receptor 2 (HER2) in the primary tumor; endocrine therapy is prescribed for patients with ER-positive (ER+) or PR-positive (PR+) cancers, while anti-HER2 drugs are prescribed for patients with HER2-positive tumors. Patients with tumors lacking expression of all three receptors have triple-negative breast cancer (TNBC) and undergo chemotherapy. However, receptor switching between primary and metastatic disease progression has been frequently observed, including approximately 20% of patient tumors switching from ER+ to ER-negative (ER-) subtype. Furthermore, similar to the presence of intratumoral genomic heterogeneity in breast cancer, a mixture of clinical subtypes may coexist across or within metastatic lesions in the same patient, thereby presenting a major clinical challenge. Accurate subtyping of the underlying clinical phenotype and identification of transcriptional patterns that emerge during treatment are therefore crucial to investigate mechanisms of resistance and inform treatment decisions.

Recent studies have shown that computational analysis of cfDNA fragmentation patterns from genome sequencing data can reveal nucleosome occupancy in the cell of origin. DNA is protected from degradation by nucleosomes when released into peripheral blood after cell death. At accessible genomic locations, e.g., actively bound transcription factor binding sites (TFBS) and open chromatin regions, nucleosomes are arranged in an organized manner that allows the access of DNA-binding proteins (Figure 7A). This nucleosome organization results in a loss of sequencing coverage, reflecting DNA degradation at unprotected binding sites, accompanied by coverage peaks at surrounding protected locations.

The application of nucleosome profiling from cfDNA has been demonstrated for cancer detection and tumor tissue of origin prediction, including analysis of short cfDNA fragments that tend to be enriched from tumor cells. Tumor subtyping from cfDNA has been explored for prostate cancer by analyzing TFBS locations, but subtyping from cfDNA appears to be unproven for other cancers. Specifically, prediction of histological subtype for breast cancer from cfDNA has not been shown. Furthermore, current cfDNA nucleosome profiling approaches are not optimized for ULP-WGS data. Although investigating the clinical phenotype of tumors from ctDNA remains challenging due to the lack of robust computational methods, there is clear potential clinical utility for guiding treatment decisions for patients with metastatic cancer.

In this study, we developed a computational framework, called Griffin, to classify tumor subtypes from nucleosome profiling of cfDNA. Griffin overcomes the challenges of current analyses for profiling nucleosome accessibility and transcriptional regulation from analysis of standard cfDNA genomic sequencing including ULP-WGS (0.1x) coverage. Griffin uses a novel GC correction procedure that is specific to DNA fragment size and therefore unique to cfDNA sequencing data. Griffin was applied to perform high performance cancer detection and tumor tissue of origin analysis. We then demonstrated the first application of breast cancer ER subtyping from cfDNA, which showed robust classification accuracy and insights into tumor heterogeneity and prognosis, all realized by analysis of ULP-WGS data. Overall, Griffin is a generalizable framework that can detect molecular changes in transcriptional regulation and chromatin accessibility from cfDNA, potentially directing personalized treatments to improve patient outcomes.

Results Griffin Framework for Nucleosome Profiling to Predict Tumor Phenotype We developed Griffin as an analytical framework that uses a GC correction procedure to accurately profile nucleosome occupancy from cfDNA. Griffin processes fragment coverage to distinguish between accessible and inaccessible features of nucleosome protection (Figure 7A). Griffin was designed to be applied to whole genome sequencing (WGS) data of cfDNA from patients with cancer and optimized to work on ULP-WGS data to quantify nucleosome protection around sites of interest (Figure 7B). Sites of interest can be selected from a variety of chromatin-based assays, e.g., assays for transposase-accessible chromatin using sequencing (ATAC-seq), and are tailored to address specific questions, including cancer detection and tumor subtyping.

The analysis workflow begins by computing genome-wide fragment-based GC bias for each sample. Then, for the region at each site of interest, fragment midpoint coverage is computed and reweighted to remove GC bias (Methods). The reason for using midpoint coverage rather than full fragment coverage is that it results in a higher amplitude nucleosome protection signal (not shown). A composite coverage profile is then computed as the average of GC-corrected coverage across a set of sites specific to the tissue type, tumor type, transcription factor (TF), or any phenotypic comparison of interest. By examining these coverage profiles around known cancer-specific and blood-specific TFs, three quantitative features were identified that distinguish between accessible and inaccessible sites: (a) coverage within a window between -30 bp and +30 bp ("center coverage"), with smaller values representing greater accessibility, (b) coverage within a window between -1000 bp and +1000 bp ("average coverage"), and (c) global nucleosome peak amplitude calculated using a fast Fourier transform ("amplitude"). These features can be used to quantify transcription factor activity or chromatin accessibility, and can be used as features to detect cancer, determine tumor subtypes, or investigate other phenotypes of interest.

Griffin reduces GC bias, thereby enabling detection of tissue-specific accessibility A novel aspect of Griffin is the implementation of fragment-based GC bias correction. In open chromatin regions, especially in TFBS, the GC content is not uniform, thereby resulting in GC-related coverage bias (Figure 8A) (Wang, J. et al. Sequence features and chromatin structure around the genomic regions bound by 119 human transcription factors. Genome Res. 22, 1798-1812 (2012)). GC bias varies between samples and among different fragment lengths within a sample (Benjamini, Y. & Speed, TP Summarizing and correcting the GC content bias in high-throughput sequencing. Nucleic Acids Research 40, e72-e72 (2012)) (Fig. 8B), which strongly affects nucleosome accessibility predictions (Fig. 8C). To correct for this GC bias, Griffin computes a global estimated average fragment coverage ("prediction") for each sample and each fragment length using a fragment length position model (Benjamini, Y. & Speed, TP Nucleic Acids Research 40, e72-e72 (2012)) (Methods, Fig. 8B). Then, when calculating the coverage profile around the site of interest, each fragment is assigned a weight based on the global predicted coverage and GC bias for its length. This correction eliminates unexpected gains (or losses) in coverage at binding sites, thereby removing technical bias and enhancing tissue-associated accessibility signals when analyzing WGS (9-25x, Figure 8C) cancer patient cfDNA and ULP-WGS (0.1-0.3x, Figure 8D).

To test the performance of nucleosome profiling after Griffin correction for GC bias, we compared the estimated TFBS accessibility to the amount of tumor-derived DNA (i.e., tumor fraction) predicted by ichorcna for ULP-WGS data from 191 MBC cfDNA samples with tumor fraction ≥ 0.1 (Adalsteinsson, V. A. et al. Scalable whole-exome sequencing of cell-free DNA reveals high concordance with metastatic tumors. Nature Communications 8, (2017)). Tumor fraction was predicted to be negatively corrected for central coverage around tumor-specific sites and positively correlated for blood-specific sites. For the blood-specific TF LYL1, we observed that central coverage in TFBS correlated positively with tumor fraction before GC correction (Pearson's r = 0.41), as predicted, but this correlation became much stronger after GC correction (Pearson's r = 0.63, Figure 8E). For the tumor-specific TF GRHL2, a negative correlation between central coverage and tumor fraction was observed as expected (Pearson's r = -0.62, not shown). Mean coverage and amplitude features also correlated with tumor fraction, but the effect of GC bias appeared to be small (not shown). For blood-specific and cancer-specific DNase I hypersensitive sites (DHS), a similar correlation was observed between nucleosome profile features and tumor fraction after GC correction (not shown).

To quantify how GC correction reduces inter-sample signal variability, we tested the central coverage in 191 MBC cfDNA ULP-WGS samples for 377 TFs in the Gene Transcription Regulation Database (GTRD) (Ulz, P. et al. Inference of transcription factor binding from cell-free DNA enables tumor subtype prediction and early detection. Nature Communications 10, 4666 (2019); Yevshin, I., et al. GTRD: A database on gene transcription regulation - 2019 update. Nucleic Acids Research 47, D100-D105 (2019)). For each factor, we compared the variability between central coverage and tumor fraction using root mean square error (RMSE) from linear regression fits before and after GC correction. For LYL1, the RMSE decreased (from 0.062 to 0.046), indicating that the data after GC correction had less inter-sample variability (Figure 8E). Similarly, for 351 (93.1%) TFs, the RMSE decreased after GC correction, indicating that the inter-sample variability decreased after accounting for the correlation between tumor fraction and center coverage (two-tailed Wilcoxon signed rank test p = 1.0 × 10 ^-58 , test statistic = 1421, Figure 8F). Furthermore, we investigated the center coverage before and after GC correction for 377 TFs for a cohort of 215 healthy donors (Cristiano, S. et al. Genome-wide cell-free DNA fragmentation in patients with cancer. Nature 570, 385-389 (2019)). As the healthy donor samples have no tumor content, the mean absolute deviation (MAD) was evaluated for each TF to compare the inter-sample variability. We found that for 365 (96.8%) TFs, the MAD was reduced after GC correction (two-tailed Wilcoxon signed-rank test p=6.28× ^10-62 , test statistic=466, Figure 8G), indicating that sample-to-sample variability was reduced for almost all TFs. Collectively, these results demonstrate that GC correction in the Griffin framework reduces the variability in chromatin accessibility signals due to sample-to-sample GC bias and allows for improved detection of tissue-specific accessibility for ULP-WGS data.

Griffin analysis in TFBS enables accurate cancer detection To determine whether Griffin can perform cancer detection, a published WGS (1-2x) dataset of cfDNA samples from healthy donors (n=215) and cancer patients (n=208) was analyzed (Cristiano, S. et al. Nature 570, 385-389 (2019)). For 377 TFs, nucleosome profiles around the TFBS were generated using nucleosome-sized (100-200 bp) fragments and three features were extracted from each profile (center coverage, mean coverage, and amplitude), for a total of 1131 features. Using logistic regression, high performance was achieved for predicting the presence of cancer, with an area under the receiver operating characteristic curve (AUC) of 0.94 (Figure 9B). The highest performance was observed for stage IV cancers (AUC=0.99), and lower performance was observed for stage I cancers (AUC=0.93). These performances may reflect the greater tumor fraction observed in late stage cancers compared to early stage cancers. Higher performance was observed for samples with tumor fraction ≧0.05 (AUC 0.99) than for samples with undetectable tumor (0 tumor fraction, AUC=0.90). Somewhat lower performance was observed for Griffin analysis around DNase I hypersensitive sites (DHS) (AUC=0.83).

To test the ability to detect cancer at ULP-WGS coverage (0.1x), Griffin was applied to the same cfDNA data downsampled to 0.1x coverage, achieving a performance of AUC 0.89 (Figure 9B). Next, since tumor-derived DNA is enriched for fragments shorter than 150 bp (Cristiano, S. et al. Nature 570, 385-389 (2019)), we tested whether we could improve the ability of this framework to detect cancer by using only short fragments. We applied Griffin to analyze fragments of only 35-150 bp in the same TFBS and observed a decrease in performance (AUC = 0.91, not shown). Finally, the results were compared to the method by Ulz et al. Inference of transcription factor binding from cell-free DNA enables tumor subtype prediction and early detection. Nature Communications 10, 4666 (2019), which analyzed all lengths of cfDNA fragments in TFBS and achieved an AUC of 0.82 for 1-2x data and an AUC of 0.55 for downsampled data (not shown). Griffin performed better using nucleosome-sized fragments or short fragments and ULP-WGS coverage, demonstrating that Griffin can accurately detect cancer using chromatin-based assays of cfDNA and various sites from ULP-WGS in a cost-effective manner.

Griffin enables accurate prediction of breast cancer subtype from ultra-low-pass WGS Breast cancer tumor classification relies primarily on accurate clinical determination of hormone receptor status by immunohistochemistry (IHC) to quantify ER expression, but no ctDNA approach exists for this application. Therefore, our goal was to determine whether Griffin could be used to predict ER subtype status from ULP-WGS (0.1×) of cfDNA from MBC patients. We analyzed 254 samples with tumor fractions greater than 0.05 from 139 patients (Adalsteinsson, VA et al. Scalable whole-exome sequencing of cell-free DNA reveals high concordance with metastatic tumors. Nature Communications 8, (2017); Stover, DG et al. Association of Cell-Free DNA Tumor Fraction and Somatic Copy Number Alterations With Survival in Metastatic Triple-Negative Breast Cancer. Journal of Clinical Oncology JCO.2017.76.003 (2018)). First, we scrutinized the Griffin profile for key factors including ESR1, FOXA1, and GATA3, which are known to be associated with ER-positive tumors in TFBS (Albergaria, A. et al. Expression of FOXA1 and GATA-3 in breast cancer: the prognostic significance in hormone receptor-negative tumours. Breast Cancer Research 11, R40 (2009)). We observed that these TFBS were more accessible in cfDNA samples from patients with ER+ metastases compared to ER-; ER+ samples had significantly reduced central coverage after accounting for tumor fraction (ANCOVA q-value for ER status <3.38× ^10-2 , not shown). To predict ER status, a logistic regression classifier was first established using features from the Griffin profiles for all 377 TFs, achieving an accuracy of 0.71 (AUC 0.79) (not shown). TFBS features computed by the method of Ulz were also used for ER subtyping, and an accuracy of 0.53 (AUC = 0.55) was observed (not shown), likely because it was not designed for ULP-WGS data.

We then used a site selection approach that was further refined by analyzing regions of differential chromatin accessibility. Using ATAC-seq data generated from 44 ER+ and 15 ER- primary breast tumors by Cancer Genome Atlas (TCGA) (Corces, M. R. et al. The chromatin accessibility landscape of primary human cancers. Science 362, eaav1898 (2018)), we identified open chromatin sites specific for each ER subtype (Methods, Fig. 10A). ER+-specific sites (n=28,170) were enriched in the TFBSs of ESR1, PGR, FOXA1, and GATA3, whereas ER-specific sites (n=41,712) were enriched in the TFBSs of STAT3 and NFKB1 (not shown). Differences in coverage profiles were observed between ER subtype-specific sites shared with accessible chromatin and those not shared with accessible chromatin in hematopoietic cells (Satpathy, A. T. et al. Massively parallel single-cell chromatin landscapes of human immune cell development and intratumoral T cell exhaustion. Nature Biotechnology 37, 925-936 (2019)), and therefore were analyzed separately (Figure 10B).

We applied Griffin to profile nucleosome accessibility in these four sets of ER subtype-specific accessible chromatin sites, thereby extracting a total of 12 features (Fig. 10B). A logistic regression classifier was established to predict ER subtypes from these chromatin accessibility features, achieving an overall accuracy of 0.81 (AUC=0.89, n=139) (Methods, Fig. 10D). Compared with samples with a small tumor fraction (accuracy 0.69, AUC=0.75, n=38, tumor fraction 0.05-0.1), samples with a large tumor fraction performed better (accuracy 0.86, AUC=0.92, n=101, tumor fraction ≥0.1) (Fig. 10D). Repeating the analysis using only short fragments (35-150 bp) did not improve performance (accuracy 0.73, AUC = 0.81), likely due to a further reduction in fragment coverage (not shown). These results illustrate the utility of using chromatin accessibility to determine cancer subtypes from ULP-WGS data and highlight the first application of ER status prediction in breast cancer from cfDNA. This finding was validated by examining cfDNA samples from two other studies and one new dataset for patients with breast cancer (see Ahuno ST, et al. Ghana Breast Health Study Team. Circulating tumor DNA is readily detectable among Ghanaian breast cancer patients supporting non-invasive cancer genomic studies in Africa. NPJ Precis Oncol. 5 (1): 83 (2021) and Zivanovic Bujak, A., et al. Circulating tumour DNA in metastatic breast cancer to guide clinical trial enrolment and precision oncology: A cohort study. PLoS medicine 17.10 (2020): e1003363) and using the trained model on the original MBC dataset, it was able to predict ER status for all samples with tumor fraction >0.05 with an accuracy of 0.92 (AUC of 0.96). Looking only at samples with tumor fraction >0.1, the accuracy was 0.96 and the AUC was 0.98. This analysis further supports the ability of Griffin to perform accurate ER status predictions on independent data sets.

Analysis of ER status from cfDNA suggests heterogeneity of tumor subtypes To further explore ER prediction, the classification results were examined to explore the proportion of tumors for samples with incorrect predictions and the pattern of ER status of primary tumors (Figure 10C). It was observed that the majority of incorrect predictions were in ER-loss (ER+ primary and ER- metastasis) samples (Figure 10F), and that the number of ER+ predictions for these patients (5 of a total of 9 ER-loss patients were predicted as ER+) was significantly different from both the number of ER+ predictions in ER- patients who had ER- primary (2-tailed Fisher's exact test p=3.7× ^10-4 ) and the number of ER+ predictions in ER+ patients who retained ER+ in both primary and metastasis (4.3× ^10-2 , 2-tailed Fisher's exact test p=0.0183.7× ^10-4 , Figure 10F). However, despite many of the ER loss patients being incorrectly predicted to remain ER+, ROC analysis of metastatic ER status prediction among patients with ER+ primaries and tumor fraction >0.1 was found to yield an AUC of 0.74, suggesting that Griffin has reasonable ability to detect ER loss among patients with ER+ primaries (FIG. 10G). However, the overall poorer performance among ER loss patients compared to patients whose ER status did not change suggests that ER+ tumor features may remain in ER loss patients or that Griffin analysis may capture a heterogeneous mixture of ER subtypes from ctDNA.

To further assess whether this observation was due to tumor heterogeneity, other metastatic biopsies obtained from patients with ER loss were examined. In many of these cases, patients were found to have additional ER+ metastatic biopsies after the initial ER- diagnosis. Two particularly interesting cases are shown in Figure 10H. Patient MBC_1413 was initially diagnosed with ER- metastatic disease on a pleural fluid biopsy, but later underwent a second metastatic biopsy on a liver metastasis that showed ER expression in 5% of cells. The first blood draw for cell-free DNA was taken immediately after this biopsy, which interestingly predicted ER+, consistent with the ER low status of the liver biopsy. Later, a third biopsy was taken from the pleural fluid, again showing ER- disease. Immediately after this biopsy, blood was drawn for cfDNA, which was also predicted to be ER-. In a second patient, MBC_1099, the first two metastatic biopsies (bone and liver) showed ER- disease. However, when cfDNA was removed several months later, the patient's subtype was predicted to be ER+ at both time points. Interestingly, when another liver biopsy was obtained between the two cfDNA draws, this biopsy showed 5% ER+ cells, potentially explaining the prediction from cfDNA of ER+. These results suggest that in some cases where Griffin fails to detect ER loss, the prediction may detect true ER subtype heterogeneity, and that Griffin can be used to monitor subtype dynamics over the course of treatment.

Discussion This study describes the development of Griffin, a new framework and analytical tool for investigating transcriptional regulation and tumor phenotype. Griffin uses a novel, cfDNA fragment length-specific normalization of GC content bias that obscures chromatin accessibility information. We demonstrate that Griffin can be used to detect cancer with high accuracy from low-pass WGS. Furthermore, we developed an approach to perform ER subtyping in breast cancer from ULP-WGS, which shows the first ER phenotype prediction from ctDNA.

Griffin is versatile and can be used for a variety of applications related to cancer. In this disclosure, we highlight use cases for cancer detection, tissue of origin, and tumor subtype. However, Griffin can also be used for any biological comparison that can delineate transcriptional regulation and chromatin accessibility differences. The application described here uses TFBS from chromatin immunoprecipitation sequencing (ChIP-seq) and accessible chromatin sites from ATAC-seq. However, Griffin differs from existing methods because it can analyze custom sites of interest specific to any biological context. These sites can be from external sources and various assays, such as ChIP-seq, DNase I hypersensitivity, ATAC-seq, or on-target cleavage and release by nucleases (cut-and-run). As additional epigenetic data is collected by the cancer research community, including from single-cell experiments (Wu, S. J. et al. Single-cell CUT&Tag analysis of chromatin modifications in differentiation and tumor progression. Nat Biotechnol 39, 819-824 (2021); Pierce, S. E., Granja, J. M. & Greenleaf, W. J. High-throughput single-cell chromatin accessibility CRISPR screens enable unbiased identification of regulatory networks in cancer. Nat Commun 12, 2969 (2021)), Griffin will be essential to advance tumor phenotyping from liquid biopsies.

Griffin is optimized for analysis of cfDNA ULP-WGS (0.1x), whereas other nucleosome profiling methods have focused on deeper coverage sequencing. Griffin utilizes site-wide analysis, as opposed to individual loci, inspired by a similar strategy used by Ulz, P. et al. Inference of transcription factor binding from cell-free DNA enables tumor subtype prediction and early detection. Nature Communications 10, 4666 (2019). Due to the novel bias correction and versatility to analyze any set of genomic regions, Griffin demonstrates better performance compared to Ulz's method for both detecting cancer from ULP-WGS data and predicting ER status. However, Griffin is not limited to low coverage data. Increasing cfDNA sequencing coverage may enable analysis of specific gene promoters and cis-regulatory elements and may provide information about gene expression (Ulz, P. et al. Inferring expressed genes by whole-genome sequencing of plasma DNA. Nature Genetics 48, 1273-1278 (2016)). Recent studies have shown promise for cfDNA methylation and cfRNA analysis for tumor phenotyping and cancer detection (Beltran, H. et al. Circulating tumor DNA profile recognizes transformation to castration- resistant neuroendocrine prostate cancer. J Clin Invest 130, 1653-1668 (2020); Wu, A. et al. Genome-wide plasma DNA methylation features of metastatic prostate cancer. J Clin Invest 130, 1991-2000 (2020); Shen, S. Y. et al. Sensitive tumour detection and classification using plasma cell-free DNA methylomes. Nature 563, 579-583 (2018); Liu, M. C. et al. Sensitive and specific multi-cancer detection and localization using methylation signatures in cell-free DNA. Annals of Oncology 31, 745-759 (2018)). (2020); Larson, M. H. et al. A comprehensive characterization of the cell-free transcriptome reveals tissue- and subtype-specific biomarkers for cancer detection. Nature Communications 12, 2357 (2021); Kang, S. et al. CancerLocator: non-invasive cancer diagnosis and tissue-of-origin prediction using methylation profiles of cell-free DNA. Genome Biology 18, 53 (2017); Chan, K. C. A. et al. Noninvasive detection of cancer-associated genome-wide hypomethylation and copy number aberrations by plasma DNA bisulfite sequencing. Proceedings of the National Academy of Sciences 110, 18761-18768 (2013)), these specimens are difficult to isolate from clinical specimens or require specialized assays. Griffin provides a cost-effective, scalable method that requires only standard low-coverage WGS of cfDNA and can be rapidly incorporated into existing platforms to predict clinical cancer phenotypes.

A limitation of binary ER classification (ER+ or ER-) is the reduced accuracy for samples with low tumor fraction (0.05-0.1); however, patients with cfDNA tumor fraction ≥10% have a poor prognosis (Stover, D. G. et al. Association of Cell-Free DNA Tumor Fraction and Somatic Copy Number Alterations With Survival in Metastatic Triple-Negative Breast Cancer. JCO 36, 543-553 (2018)) and would be more beneficial to monitor their tumor. It may be possible to improve the performance of ER subtyping for samples with low tumor fraction using additional sequencing depth or joint analysis of multiple cfDNA time points from the same patient.

Applying Griffin to predict ER status from cfDNA in MBC patients provides interesting insights into tumor heterogeneity and potential explanations for misclassified predictions. Interestingly, for patients with ER- tumors by IHC, we observed a significant enrichment of the ER+ prediction when the primary tumor was ER+. Two patients with ER loss and ER+ prediction had both subtypes in their metastases. Importantly, this subtype heterogeneity and switching is not typically captured in a single metastasis biopsy, but these results demonstrate the feasibility of using ctDNA to monitor subtype heterogeneity over time during treatment using ER probability.

Breast cancer subtyping was focused on ER prediction because ER status has important value in predicting which endocrine treatments may be beneficial (E.B.C.T.C. Relevance of breast cancer hormone receptors and other factors to the efficacy of adjuvant tamoxifen: patient-level meta-analysis of randomised trials. The Lancet 378, 771-784 (2011)). PR expression is also determined in the clinic, and ER-/PR+ tumors are considered hormone receptor positive, but these are rare and not reproducible or of low prognostic value (Hefti, M.M. et al. Estrogen receptor negative/progesterone receptor positive breast cancer is not a reproducible subtype. Breast Cancer Research 15, R68 (2013)). In the cohort, only 2 of 139 patients (1.4%) were ER-/PR+. HER2 overexpression has important implications for prognosis and treatment decisions such as trastuzumab (Slamon, D. J. et al. Human breast cancer: correlation of relapse and survival with amplification of the HER-2/neu oncogene. Science 235, 177-182 (1987)). However, an insufficient number of open chromatin sites have been identified that are specific for distinguishing HER2 status. Alternatively, ERBB2 copy number amplification can be assessed from ctDNA genomic analysis, since ERBB2 (encoding the HER2 protein) is amplified in approximately 20% of breast cancers (Curtis, C. et al. The genomic and transcriptomic architecture of 2,000 breast tumours reveals novel subgroups. Nature 486, 346-352 (2012)). Alternatively, a model for predicting PAM50 status may be useful as it may be a better prognostic indicator than ER/PR/HER2 IHC alone (Nielsen, T. O. et al. A Comparison of PAM50 Intrinsic Subtyping with Immunohistochemistry and Clinical Prognostic Factors in Tamoxifen-Treated Estrogen Receptor-Positive Breast Cancer. Clinical Cancer Research 16, 5222-5232 (2010)).

The Griffin framework is a unique advancement of our previous methods for analyzing genomic alterations and estimating tumor fraction from ULP-WGS of cfDNA (Adalsteinsson, V. A. et al. Scalable whole-exome sequencing of cell-free DNA reveals high concordance with metastatic tumors. Nature Communications 8, (2017)). Collectively, these methods form a set of tools to establish a new paradigm for investigating both tumor genotype and phenotype from ULP-WGS of cfDNA. Griffin has the potential to reveal clinically relevant tumor phenotypes, thereby supporting therapeutic resistance studies, informing treatment decisions, and accelerating applications in cancer precision medicine.

Methods Griffin: Calculating GC bias GC content affects the efficiency of amplification and sequencing, resulting in different predicted coverage (coverage bias) for fragments with different GC content and fragment length. This is called GC bias and is unique to each sample. We calculated the GC bias of each bam file using an implementation of the method developed by Benjamini and Speed in 2012 (Benjamini, Y. & Speed, TP Summarizing and correcting the GC content bias in high-throughput sequencing. Nucleic Acids Research 40, e72-e72 (2012)) previously implemented in deepTools (Ramirez, F. et al. deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Research 44, W160-W165 (2016)). However, unlike the implementation in deepTools, which assumes that all fragments have the same length, we used a "fragment length model" that calculates separate GC bias curves for each fragment length. This is useful for cfDNA, where different samples may have different fragment size distributions and different fragment lengths have biological significance. Before performing the GC bias calculation, we identified all mappable regions of the genome using the Umap multi-read mappability track on 50 bp reads downloaded from the UCSC genome browser (Karimzadeh, M., Ernst, C., Kundaje, A. & Hoffman, MM Umap and Bismap: quantifying genome and methylome mappability. Nucleic Acids Research 46, e120-e120 (2018)) (hgdownload.soe.ucsc.edu/gbdb/hg38/hoffmanMappability/k50.Umap.MultiTrackMappability.bw). Using pybedtools (Dale, RK, Pedersen, BS & Quinlan, AR Pybedtools: a flexible Python library for manipulating genomic datasets and annotations. Bioinformatics 27, 3423-3424 (2011)), we found mappable regions (defined as mappability score = 1) for hg38 downloaded from the UCSC table browser (genome.ucsc.edu/cgi-bin/hgTables) as well as additional excluded regions with known mapping problems, including the Encode Unified Exclusion List (encodeproject.org/files/ENCFF356LFX/), centromeres, correction patches, and alternative haplotypes. We then explored all remaining regions of the genome and counted the observed GC content of all possible fragments that overlap these positions for each fragment length. The frequency of occurrence of each observed GC content for each fragment length is the "genomic GC frequency". We then developed a "griffin GC bias" pipeline to compute the GC bias in a given bam file. This pipeline obtained the bam files, the bedGraph files of valid (mappable, non-excluded) regions, and the genomic GC frequencies for these regions. For each given sample, we fetched all reads that aligned to valid regions on the autosomes using pysam (github.com/pysam-developers/pysam) (Li, H. et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics 25, 2078-2079 (2009)). We counted the number of reads observed for each length and GC content, and excluded reads with poor mapping quality (<20), duplicates, unpaired reads, and reads that were not quality controlled. These read counts are the "GC counts" of the sample. We then divided the GC count for the sample by the GC frequency for the genome to obtain the GC bias for that bam file, and normalized the average GC bias for each fragment length to 1, resulting in a GC bias value for all fragment size and GC content combinations (other than those not observed in the genome). We then smoothed the GC bias curve. For each fragment size, we obtained all GC bias values for fragments of similar length (+/- 10 bp). We sorted these values by fragment GC content to create a vector of GC bias values for fragments of similar size. We then smoothed this vector by taking the median of the k nearest neighbors (k = 5% of the vector length or 50, whichever is greater), and repeated for each possible fragment length. We then normalized to an average GC bias of 1 for each possible fragment length (excluding unobserved GC content) to generate smoothed GC bias values for all possible fragment lengths and observed GC content in the genome.

の重みを割り当て、断片中点の位置を同定した。部位を１５ｂｐのビンに分割し、各ビン内の重み付けされた断片中点を合計して、ＧＣ補正された中点カバレッジプロファイルを得た（概略図について図１ｂを参照されたい）。次に、既知のマッピング問題を有する領域とオーバーラップするビン（Ｇｒｉｆｆｉｎ：ＧＣの偏りの算出において記載されている）およびマッピングできない位置を少なくとも１つ有するビンを除外した。非常に高いカバレッジ（平均を１０標準偏差上回る）を有するビンも同定し、これらのビンを取り除いた。これを部位一覧上の全ての部位について繰り返し、全ての部位の平均を取って（それらの部位内の除外されたビンは無視する）、その部位一覧についてのカバレッジプロファイルを生成した。次いで、Ｓａｖｉｔｚｋｙ－Ｇｏｌａｙフィルターをウインドウ長１６５ｂｐおよび多項式の次数３で使用してカバレッジプロファイルを平滑化した。最後に、異なる深度の試料を比較できるようにするために、カバレッジプロファイルを平均カバレッジ１に対して正規化し、中心的領域（＋／－１０００ｂｐ）をさらなる解析のために保持した。 Griffin: Nucleosome Profiling The griffin nucleosome profiling pipeline was designed to perform nucleosome profiling around sites of interest. The pipeline takes a bam file and a list of sites, as well as a variety of other parameters described below. For a given bam file and list of sites, all reads within a window (-5000 to +5000 bp) around each site were fetched using pysam (excluding those for which quality control measures were not available). Read pairs were then filtered by fragment length to select those that fell within the fragment length range (100-200 bp unless otherwise specified). For each read pair, the GC bias for the fragment was determined and the fragment length was chosen.

The sites were assigned weights of 0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 1.6, 1.7, 1.8, 1.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, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 112, 123, 131, 142, 143, 152, 153, 164, 176, 185, 197, 198, 199, 102, 113, 123, 13 Finally, to be able to compare samples at different depths, coverage profiles were normalized to a mean coverage of 1 and central regions (+/- 1000 bp) were retained for further analysis.

Griffin: Quantification of Nucleosome Profile Features To quantify the coverage profiles, three features were extracted from each coverage profile. First, a "mean coverage" value was calculated +/- 1000 bp from the site. Second, the coverage value at the site (+/- 30 bp) was calculated. And third, the amplitude of the nucleosome peaks surrounding the site was calculated using a Fast Fourier Transform for a window of +/- 960 bp from the site (as implemented in Numpy (Harris, CR, Millman, KJ, van der Walt, SJ et al. Array programming with NumPy. Nature 585, 357-362 (2020)) and taking the amplitude of the 10th frequency term. This window and frequency were chosen by the nucleosome peak spacing observed at the active site (190 bp), resulting in approximately 10 peaks within a window of +/- 960 bp.

Early-stage cancer and healthy donor cfDNA samples - DELFI dataset Whole-genome sequencing (WGS) cfDNA from patients with various types of early-stage cancer and from healthy donors was obtained from an existing dataset published in Cristiano, S. et al. Genome-wide cell-free DNA fragmentation in patients with cancer. Nature 570, 385-389 (2019). The bam file was downloaded from EGA (dataset ID: EGAD00001005339). The data consisted of 1-2x low-pass whole-genome sequencing of 100 bp paired-end Illumina sequencing reads. For our analysis, we used a subset of samples that included 1-2x WGS of cfDNA from 208 previously untreated cancer patients and 215 healthy donors. These are the samples used in Cristiano et al.'s cancer detection analysis. cfDNA tumor fraction was estimated using ichorcNA (Adalsteinsson, VA et al. Scalable whole-exome sequencing of cell-free DNA reveals high concordance with metastatic tumors. Nature Communications 8, (2017)). All 215 healthy donors in the dataset were used to create a hg38 normal panel (PoN) with bin size of 1 mb. ichorcNA was then run on all cancer and healthy samples to estimate tumor fraction. ichorcNA_fracReadsInChrYFormale was set to 0.001. Defaults were used for all other settings.

Early stage lung cancer and healthy donor cfDNA samples - LUCAS dataset Whole genome sequencing (WGS) cfDNA from a prospective study of patients with and without lung cancer was obtained from an existing dataset published by Mathios and coworkers (Mathios, D. et al. Detection and characterization of lung cancer using cell-free DNA fragmentomes. Nat Commun 12, 5060 (2021)). The bam file was downloaded from EGA (dataset ID: EGAD00001007796). The data consisted of 1-2x low-pass whole genome sequencing of 100 bp paired-end Illumina sequencing reads. For our analysis, we used a subset of samples described in the paper as the "LUCAS" cohort and a second subset of samples described as the validation cohort. The LUCAS cohort included 158 patients with no history of cancer and no future cancer diagnosis, and 129 patients diagnosed with lung cancer within days (0-44 days) of blood draw. The validation cohort included 46 patients with cancer and 385 patients without cancer. All samples were realigned to hg38 as described in Sequence Data Processing below. Tumor fraction was determined using ichorcNA as described above, with a normal panel constructed from 54 separate non-cancer samples from this same study.

Metastatic Breast Cancer (MBC) and Healthy Donor cfDNA Samples WGS of cfDNA from patients with metastatic breast cancer (MBC) and healthy donors was obtained from an existing dataset (Adalsteinsson, VA et al. Nature Communications 8, (2017)). BAM files were downloaded from dbGaP (accession code: phs001417.v1.p1). The data consisted of approximately 0.1× ultra-low-pass whole genome sequencing (ULP-WGS) of 100 bp paired-end Illumina sequencing reads. For our analysis, we used a subset of 254 samples with WGS coverage >0.1×, tumor fraction >0.05×, and known estrogen receptor (ER) status. Of these 254 samples, 133 were ER positive (from 74 unique patients) and 121 were ER negative (from 65 unique patients). Coverage and tumor fraction metrics were obtained from supplementary data in the publication (Adalsteinsson, VA et al. Nature Communications 8, (2017)). Primary and metastatic ER status was determined by immunohistochemistry and also abstracted from diagnostic records. In addition, deep (9-25x) WGS from two MBC patients (MBC_315 and MBC_288) from the same source and deep (17-20x) WGS from two healthy donors (HD45 and HD46) from the same source were used for pipeline design and validation.

To train and evaluate the ER status classifier, information on ER status from the medical records was used and each sample was labeled as ER+ or ER-. If metastatic ER status was known, samples were labeled according to that status. If metastatic ER status was not known, samples were labeled according to primary tumor ER status (20 samples from 11 patients). ER low samples (11 samples from 6 patients) were labeled as ER positive for the purposes of the binary classifier. For three patients (MBC_1405, MBC_1406, MBC_1408), information was available on multiple metastatic biopsies with different ER status. In these cases, the last biopsy taken was used for the purposes of the binary ER status classifier.

Human subjects WGS of cfDNA samples from patients with MBC were obtained from the existing study mentioned above (Adalsteinsson, VA et al. Nature Communications 8, (2017)). Additional information including primary ER status, metastatic ER status, and survival time was abstracted from diagnostic records. Use of this data was approved by the Institutional Review Board (Dana-Farber Cancer Institute IRB protocol identifiers 05-246, 09-204, 12-431 [NCT01738438; Closure effective date 6/30/2014]).

Sequence data processing All sequencing data used in this study were realigned to the hg38 version of the human genome (downloaded from hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz). Bam files were unmapped from previous alignments using Picard SamToFastq (Picard Toolkit. (Broad Institute, 2021)). The unmapped bam files were then realigned to the human reference genome using the following procedure, following GATK best practices (DePristo, MA et al. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nature Genetics 43, 491-498 (2011)). Fastq files were realigned using BWA-MEM (Li, H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. 00, 1-3 (2013)). Files were then sorted with samtools (Danecek, P. et al. Twelve years of SAMtools and BCFtools. GigaScience 10, 1-4 (2021)), duplicates were marked with Picard, and base recalibration was performed with GATK using known polymorphisms downloaded from the following locations: console.cloud.google.com/storage/browser/genomics-public-data/resources/broad/hg38/v0/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz and ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh38p7/VCF/GATK/All_20180418.vcf.gz.

Transcription factor binding site (TFBS) selection Transcription factor binding sites (TFBS) were downloaded from the GTRD database (Yevshin, I., GTRD: A database on gene transcription regulation-2019 update. Nucleic Acids Research 47, D100-D105 (2019)). This database contains a compilation of ChIP-seq data from various sources. For our analysis, we used metacluster data (version 19.10, downloaded from gtrd.biouml.org/downloads/19.10/chip-seq/Homo%20sapiens_meta_clusters.interval.gz), which contains metapeaks observed in one or more ChIP-seq experiments. The GTRD database contains some ChIP-seq experiments for targets that are not transcription factors (TFs). These were eliminated by comparing with the list of TFs with known binding sites in the CIS-BP database (Weirauch, MT et al. Determination and Inference of Eukaryotic Transcription Factor Sequence Specificity. Cell 158, 1431-1443 (2014)) (v2.00, downloaded from cisbp.ccbr.utoronto.ca/bulk.php). Site locations were identified as the average of "start" and "end". TFs with less than 10,000 sites on autosomes were eliminated. For each remaining TF, the top 10,000 sites were selected by choosing the one with the highest "peak.count" (the number of times the peak was observed across all experiments).

DNase I hypersensitive site selection DNase I hypersensitive sites for various tissue types were downloaded from zenodo.org/record/3838751/files/DHS_Index_and_Vocabulary_hg38_WM20190703.txt.gz (Meuleman, W. et al. Index and biological spectrum of human DNase I hypersensitive sites. Nature 584, 244-251 (2020)). These sites were split by tissue type, resulting in a total of 16 sites. The "summit" column was used as the site location. Sites were sorted by the number of samples in which the site was observed ("numsamples"), and the top 10,000 most frequently observed sites for each tissue type were selected.

ATAC-seq site selection for ER subtyping. Assay for transposase accessible chromatin using sequencing (ATAC-seq) site accessibility for primary breast cancer samples from Cancer Genome Atlas (TCGA) downloaded from the TCGA ATAC-seq hub (gdc.cancer.gov/about-data/publications/ATACseq-AWG) (Corces, MR et al. The chromatin accessibility landscape of primary human cancers. Science 362, eaav1898 (2018)). A file containing raw counts for all cancer type-specific peaks was downloaded ("All cancer type-specific count matrices in raw counts"), and the file containing breast cancer-specific peaks was used ("BRCA_raw_counts.txt"). The locations of these sites and patient metadata were obtained from a supplementary table in the paper (Corces, MR et al. The chromatin accessibility landscape of primary human cancers. Science 362, eaav1898 (2018)). Autosomal sites were retained for further analysis for a total of 211,938 sites. Differentially accessible sites between ER+ (n=44) and ER- (n=15) tumors were identified using DESeq2 software (Love, MI, Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15, 550 (2014)). The software was run using the default settings described in the "quick start" guide. Differential expression experiments were performed using the "DESeq" and "results" functions, followed by log fold change shrinkage using the "lfcShrink" function. Sites with q values below ^5x10-4 were selected. Additionally, the selected sites were further filtered based on the log2 fold change between ER+ and ER- tumors. Sites with a log2 fold change >0.5 were classified as ER+ specific, while sites with a log2 fold change <-0.5 were classified as ER- specific. These site lists were further divided into sites shared with hematopoietic cells and sites not shared with hematopoietic cells. Hematopoietic sites were obtained from a database of single-cell ATAC-seq data (Satpathy, AT et al. Massively parallel single-cell chromatin landscapes of human immune cell development and intratumoral T cell exhaustion. Nature Biotechnology 37, 925-936 (2019)) (GEO accession number: GSE129785, peak file available here: ftp.ncbi.nlm.nih.gov/geo/series/GSE129nnn/GSE129785/suppl/GSE129785%5FscATAC%2DHematopoiesis%2DAll%2Epeaks%2Etxt%2Egz). Hematopoietic peaks were lifted over to hg38 using the UCSC liftover command line tool and sites that changed size during the liftover (0.2% of peaks) were discarded. BRCA ATAC-seq sites overlapping with hematopoietic sites (overlapping peaks were defined as centers of sites within 500 bp of each other). This was performed using pybedtools intersect (Dale, RK, et al. Bioinformatics 27, 3423-3424 (2011); Quinlan, AR & Hall, IM BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26, 841-842 (2010)). This resulted in a total of four differential site lists: ER-positive sites not shared with hematopoietic cells (18,240 sites), ER-positive sites shared with hematopoietic cells (9,930 sites), ER-negative sites not shared with hematopoietic cells (19,347 sites), and ER-negative sites shared with hematopoietic cells (22,365 sites).

These differential ATAC-seq site lists were then overlapped with the top 10,000 sites for each of 338 transcription factors (TFs) using pybedtools intersect. Overlapping site pairs were defined as those with less than 500 bp between the centers of the sites. Each differential ATAC-seq site list was compared with each list of TFBSs, and the total number of ATAC sites in a given list that overlapped with one or more TFBSs was recorded.

Griffin's score before and after GC correction Correlation of tumor fraction in TFBS Nucleosome profiling was performed on the top 10,000 sites for each of 377 transcription factors (TFs) for 191 MBC ULP samples with tumor fraction >0.1, with and without GC correction. For each TF, the association between center coverage and tumor fraction was modeled using scipy. stats. lineregress (Virtanen, P. et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat Methods 17, 261-272 (2020)) to generate Pearson correlation (r) and best-fit lines. Root mean square error (RMSE) was calculated from the best-fit lines. This was done both before and after GC correction, as illustrated for Lyl-1 in Figure 2e. For all 377 TFs, RMSE values before and after GC correction were compared using the Wilcoxon signed rank test (two-tailed).

Mean Absolute Deviation (MAD) in TFBS
Nucleosome profiling with and without GC correction was performed on the top 10,000 sites for each of the 377 TFs for 215 healthy donors. For each TF, the MAD of the central coverage values was calculated both before and after GC correction. For all 377 TFs, the MAD values before and after GC correction were compared using the Wilcoxon signed rank test (two-tailed).

Machine learning, bootstrapping, and performance evaluation procedures To detect cancer, predict histology, or predict ER subtype, we used logistic regression and ridge regularization (i.e., L2 standard) implemented in scikit-learn (Pedregosa, F. et al. Scikit-learn: Machine Learning in Python. Journal of Machine Learning Research 12, 2825-2830 (2011)). Bootstrapping and model fitting were performed after scaling all feature values to mean 0 and standard deviation 1. We trained our models and assessed their performance using the following bootstrapping procedure. First, n samples were selected with replacement from the total set of n samples and used as the training set. The unselected samples were used as the test set. Then, using 10-fold cross-validation on the training set, the parameter "C" (the inverse of the regularization strength) was selected from the following choices: ^10-5 , ^10-4 , ^10-3 , ^10-2 , ^10-1 , ¹⁰⁰ , ¹⁰¹ , ^102. To account for class imbalance in the data used, the "class weight" parameter was set to "balanced" to adjust the weights of the samples inversely proportional to their class occurrence frequency. The final model was trained on all the training data using the selected regularization strength. Finally, this model was tested on the test set and the performance (accuracy and AUC values) as well as the probability from each sample were recorded. A new training set was then selected and the procedure was repeated for 1000 iterations (for cancer detection and tissue of origin analysis) or 1000 iterations (for breast cancer subtyping). After completing the bootstrap iterations, AUC and accuracy were calculated from each bootstrap iteration and used to generate the mean and 95% confidence interval for each of these values. To visualize the mean ROC curve, the median probability from all bootstraps whose samples were included in the test set was used. The same median probability was used for further downstream analysis, including comut plot bar plots and time series.

Features used for cancer detection classification To detect cancer, the logistic regression approach described above was applied to establish four different models using four different sets of features extracted from pan-cancer patient samples and healthy donor samples. First, nucleosome profiling was performed on these samples (selecting fragments with lengths between 100 and 200 bp) for 377 TFs selected from the GTRD database. Three features (as described above) were extracted from each coverage profile, for a total of 1,014 features. PCA was used to reduce the dimensionality of the data and features explaining 80% of the variance were selected. These PCA components were then used as inputs for the logistic regression models.

Second, nucleosome profiling was performed on these same samples and sites, but only "short" fragments (35-150 bp) were selected to be counted in the nucleosome profile.

Third, we downsampled these samples to approximately 0.1x coverage (procedure below) and performed nucleosome profiling on the same 377 TFs by selecting fragments of 100-200 bp in length.

Fourth, we used the original (non-downsampled) samples to perform nucleosome profiling on the 16 tissue-specific DHS site lists listed above. We extracted the same three features from each site profile, for a total of 48 features.

Downsampling of pan-cancer and healthy donor cfDNA sequencing data. 1-2x WGS pan-cancer patient and healthy donor bam files aligned to hg38 were downsampled using Picard DownSampleSam. Probabilities used in DownSampleSam were calculated based on 2,463,109 read-pair targets, resulting in approximately 0.11x coverage as calculated by Picard CollectWgsMetrics. Downsampled bam files were realigned to hg19 for use in the Ulz pipeline. The realignment procedure was the same as above, but using the hg19 genome (downloaded from hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/hg19.fa.gz) and known polymorphic sites of hg19 for base recalibration (downloaded from gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg37/Mills_and_1000G_gold_standard.indels.hg37.vcf.gz and ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh37p13/VCF/GATK/All_20180423.vcf.gz).

ER status classification in the MBC cohort To predict ER status, the logistic regression approach described above was applied to features extracted from MBC patient samples. Because there were multiple samples for some patients, the bootstrapping procedure was modified so that 139 patients (rather than samples) were selected with replacement from the full set of 139 patients. For each selected patient, all samples from that patient were added to the training set (if a patient was selected multiple times, all of their samples were included multiple times). This ensured that separate samples from the same patient (biological replicates) could not appear in both the training and test sets. Samples from unselected patients were used as the test set.

Using these training and test sets, three different models were established based on three different sets of features. First, nucleosome profiling using 100-200 bp fragments was applied to 377 TFs from GTRD, extracting three features per profile, for a total of 1131 features. PCA was then used to identify components explaining 80% of the variance, as described above. Second, nucleosome profiling using 100-200 bp fragments was applied to the four ER differential ATAC seq lists, extracting three features per profile, for a total of 12 features. Finally, nucleosome profiling using 35-150 bp fragments was applied to the four ER differential ATAC seq lists, extracting three features per list, for a total of 12 features.

For model evaluation, only the first time point for each patient in the test set was included when calculating accuracy and AUC for each bootstrap replicate, preventing a small number of patients with many samples from significantly affecting the scores.

Transcription factor profiling using the pipeline from Ulz et al.

We downloaded the transcription factor profiling pipeline (Ulz, P. et al. Inference of transcription factor binding from cell-free DNA enables tumor subtype prediction and early detection. Nature Communications 10, 4666 (2019)) published by Ulz and coworkers from Github (github.com/PeterUlz/TranscriptionFactorProfiling) and ran it using the following procedure described in the paper. Since the pipeline was written for the hg19 version of the genome, we used bam files aligned to hg19. The script was modified to work in python3. Reads within each bam were trimmed to 60 bp using "trim from bam single end" modified to skip unaligned reads. ichorcna was run on the original (untrimmed) bam using the default ichorcna settings for hg19, except that the bin size was modified to 50,000 bp and no regular panel was used. Transcription factor profiling analysis was then performed on the trimmed bam using the script run_tf_analyses_from_bam.py, with the "-calccov" and "-a tf_gtrd_1000sites" options, and the ichorcna corrected depth file as the "-norm-file". This performed transcription factor profiling on 1,000 sites for each of the 504 TFs. Finally, the scoring pipeline was run. For each of the 504 TFs in the accessibility output file (Accessibility1KSitesAdjusted.txt), we used the high frequency amplitude ("HighFreqRange") as a feature in a logistic regression model using the same bootstrapping scheme as above.

Data Availability Sequencing data used in this study were obtained from dbGaP (accession phs001417.v1.p1) and EGA (dataset ID EGAD00001005339).

Code Availability The Griffin software and subtype classifier tools are available at github.com/adoebley/Griffin. Code for the analyses and machine learning models is available at github.com/adoebley/Griffin_analyses.

Example 2
Example 1 above is a proof of concept that sequence analysis applying an embodiment of the Griffin workflow can enhance sequence signals from low-pass sequencing data with sufficient power and specificity to allow breast cancer subtype determination.In this example, the application of the Griffin workflow is extended to other cancer types and uses data from alternative sequence profiling platforms.Specifically, histone modification profiling is performed using cut and run on prostate cancer cells of different subtypes.As in Example 1, the Griffin workflow provides robust signals to clearly distinguish different subtypes of prostate cancer, thereby demonstrating the power and flexibility of the analysis workflow.

1. Background Metastatic castration-resistant prostate cancer (mCRPC) describes the stage at which the disease develops resistance to androgen deprivation therapy and becomes fatal. Androgen receptor signaling inhibitors (ARSIs) designed to treat CRPC suppress androgen receptor (AR) activity and improve survival, but these therapies ultimately fail. Since ARSIs were adopted as the standard of care for mCRPC, there has been a significant increase in the incidence of treatment-resistant tumors with neuroendocrine (NE) differentiation and small cell carcinoma features. These high-grade tumors may develop through a resistance mechanism of transdifferentiation from AR-positive adenocarcinoma (ARPC) to NE prostate cancer (NEPC) that lacks AR activity. Based on AR activity and expression of the NE gene, additional phenotypes may also occur, including AR-low prostate cancer (ARLPC) and double-negative prostate cancer (DNPC; AR-null/NE-null). Differentiating between prostate cancer subtypes is of clinical importance given their differential response to therapeutic agents, but can be complicated by the need to perform biopsies to diagnose the histological status of the tumor: invasive procedures are expensive, are associated with morbidity, there are subsets of tumors for which biopsy is unavailable, and bone sites pose particular challenges regarding sample quality.

Circulating tumor DNA (ctDNA), released from tumor cells into the blood as cell-free DNA (cfDNA), is a non-invasive "liquid biopsy" solution to harness tumor molecular information. Analyzing ctDNA to detect mutations and copy number alterations has been useful to classify genomic subtypes of CRPC tumors. However, defining loss of TP53 and RB1 in NEPC does not necessarily result in NE transdifferentiation. Instead, ARPC and NEPC tumors are associated with distinct transcriptional reprogramming. Analyzing methylation of cfDNA in mCRPC to profile the epigenome shows promise for distinguishing phenotypes, but requires specialized assays such as bisulfite treatment, enzymatic digestion, or immunoprecipitation.

The majority of cfDNA represents DNA protected by nucleosomes when released into the circulation from dying cells, thus resulting in DNA fragmentation that reflects enzymatic cleavage by nonrandom nucleases. A new approach to analyze cfDNA fragmentation patterns from plasma to investigate cancer can be implemented directly from standard whole genome sequencing (WGS). cfDNA fragments have a characteristic size of 167 bp consistent with protection by a single core nucleosome octamer and histone linkers, although the size distribution can vary between healthy individuals and cancer patients. Recent studies have demonstrated that nucleosome occupancy of cfDNA at transcription start sites (TSS) and transcription factor binding sites (TFBS) can be used to infer gene expression and transcription factor (TF) activity from cfDNA. However, nucleosome positioning and spacing are dynamic in cases of active and repressed gene regulation. A detailed understanding of nucleosome organization and positioning patterns associated with transcription regulation has not been fully explored for cfDNA.

A major challenge of ctDNA analysis is the low tumor content (tumor fraction) in patient plasma samples. In contrast, plasma from patient-derived xenograft (PDX) models can contain nearly pure human ctDNA after bioinformatics elimination of mouse DNA reads. This provides an ideal resource for investigating ctDNA properties, developing new analytical tools, and validating both genetic and phenotypic features by comparison with corresponding tumors. In this study, we performed deep WGS of ctDNA from mouse plasma across 24 CRPC PDX lines with diverse phenotypes. By applying the computational method described in Example 1 to nucleosome patterns, we comprehensively explored regulatory loci, TFBS, TSS, and open chromatin sites across genes to reveal transcriptional regulation associated with mCRPC phenotypes. Finally, we designed a probabilistic model to accurately classify treatment-resistant tumors into different phenotypes and validated its performance for 159 plasma samples from three mCRPC patient cohorts. Overall, these results highlight the potential utility of ctDNA to pinpoint transcriptional regulation of tumor phenotype, with potential diagnostic applications in precision cancer medicine.

Results Comprehensive resource of representative tumors and liquid biopsies from patient-derived xenograft (PDX) models of advanced prostate cancer We used 26 models from the LuCaP PDX series of advanced prostate cancer with a well-defined mCRPC phenotype (Nguyen et al. , (2017). LuCaP Prostate Cancer Patient-Derived Xenografts Reflect the Molecular Heterogeneity of Advanced Disease and Serve as Models for Evaluating Cancer Therapeutics. The Prostate 77, 654-671). The models consisted of 18 cases classified as ARPC, 2 cases classified as AR-low and NE-negative prostate cancer (ARLPC), and 6 cases classified as NEPC (Figure 11A). For each PDX line, mouse plasma from 7-10 mice was pooled, cfDNA extracted, and deep whole genome sequencing was performed (WGS; mean 38.4x coverage, range 21-85x) (Methods, Fig. 11A). Human ctDNA constituted more than 10% of samples in 25 lines (mean 52.9%, range 10.6-96%), with a significantly higher human fraction in NEPC samples (mean 85.1%, range 77.1-96%, two-tailed Mann-Whitney U test p= ^9.6x10-4 ) (Fig. 11B). Bioinformatics subtraction of mouse sequenced reads was used to obtain nearly pure human ctDNA data (Methods). After subsequent filtering by human ctDNA sequencing coverage, 24 PDX lines were retained for further analysis (16 ARPC, 6 NEPC, 2 ARLPC; mean 20.5×, range 3.8-50.6×). Targeted cleavage and release by nuclease (CUT&RUN) was performed on corresponding tumors to profile H3K27ac, H3K4me1, and H3K27me3 histone post-translational modifications (PTMs) (Meers et al. , (2019). Peak calling by Sparse Enrichment Analysis for CUT&RUN chromatin profiling. Epigenetics & Chromatin 12, 42; Skene and Henikoff (2017). An efficient targeted nuclease strategy for high-resolution mapping of DNA binding sites. ELife 6, e21856.). It was hypothesized that nucleosome organization inferred from ctDNA reflects the transcriptional activity landscape regulated by histone PTMs (Zhou et al. (2011). Charting histone modifications and the functional organization of mammalian genomes. Nat Rev Genet 12, 7-18).

To investigate transcriptional regulation in the mCRPC phenotype from ctDNA, four different features were investigated: local promoter coverage, nucleosome positioning, fragment size analysis, and combined TFBS and open chromatin site analysis using the Griffin framework (Example 1, and Doebley et al. (2021). Griffin: Framework for clinical cancer subtyping from nucleosome profiling of cell-free DNA. MedRxiv 2021.08.31.21262867) (Figure 11A, Methods). Three different local regions within the ctDNA were analyzed: all gene promoters, and within the transcribed regions of genes, and sites of histone PTMs guided by cut-and-run analysis. Next, ctDNA was analyzed for transcription factor binding sites (TFBS) and open chromatin regions. For each transcription factor (TF), ctDNA coverage in the TFBS was aggregated into a composite profile representing inferred activity (Example 1 and Doebley et al. , 2021; Ulz et al. (2019). Inference of transcription factor binding from cell-free DNA enables tumor subtype prediction and early detection. Nature Communications 10, 4666). Similarly, features in the composite profile of subtype-specific open chromatin regions were extracted to analyze the signature of chromatin accessibility in ctDNA. In summary, a multi-omics sequencing dataset was assembled from corresponding tumors and plasma for a total of 24 PDX lines, providing a unique molecular resource and platform for developing transcriptional regulatory signatures predictive of tumor phenotype from ctDNA.

Characterization of AR and ASCL1 Transcriptional Activity in PDX Phenotypes by Analysis of Tumor Histone Modifications and ctDNA Prostate cancer phenotypes in mCRPC patients have distinct transcriptional signatures, which are also observed in the LuCaP PDX line (Labrecque et al. (2021b). RNA Splicing Factors SRRM3 and SRRM4 Distinguish Molecular Phenotypes of Castration-Resistant Neuroendocrine Prostate Cancer. Cancer Research 81, 4736-4750). Transcriptional activity in the different tumor phenotypes was further characterized by investigating epigenetic regulation via histone PTMs. Extensive peak regions for H3K4me1 (median 17,643 regions, range 1,894-64,934), H3K27ac (median 7,093, range 1610-34,047), and H3K27me3 (median 8,737, range 2,024-42,495) were identified in tumors from 24 PDX lines and an additional 9 LuCaP PDX lines when only tumors were available (a total of 25 ARPC, 2 ARLPC, and 6 NEPC) (Methods). Using unsupervised clustering and principal component analysis (PCA), we identified putative actively regulated enhancer and promoter regions (H3K27ac, H3K4me1) and gene-repressive heterochromatin marks (H3K27me3) that are ARPC-, ARLPC-, and NEPC-phenotype-specific (Soares et al. (2017). Determinants of Histone H3K4 Methylation Patterns. Molecular Cell 68, 773-785.e6).

AR and ASCL1 are two important differentially expressed TFs with known regulatory roles in ARPC and NEPC phenotypes, respectively (Brady et al. (2021). Temporal evolution of cellular heterogeneity during the progression to advanced AR-negative prostate cancer. Nat Commun 12, 3372; Cejas et al. (2021). Subtype heterogeneity and epigenetic convergence in neuroendocrine prostate cancer. Nat Commun 12, 5775; Rapa et al. (2008). Human ASH1 expression in prostate cancer with neuroendocrine differentiation. Mod Pathol 21, 700-707; Wang et al. (2020). Molecular tracing of prostate cancer lethality. Oncogene 39, 7225-7238). When we closely examined AR binding sites in ARPC tumors, we observed increased signal from nucleosomes adjacent to H3K27ac PTMs compared to other phenotypes (mean area under peak profile was 18.46 vs. 15.08 in ARLPC and 10.63 in NEPC) (FIG. 12A, Methods). In ARPC, the strongest signal of H3K27ac was also observed in nucleosome depleted regions (NDRs) (coverage reduction of 1.54 vs. 0.78 in ARLPC and 0.41 in NEPC). Conversely, in NEPC tumors, stronger signals were observed in nucleosomes with H3K27ac PTMs adjacent to ASCL1 binding sites (mean area under peak profile 62.65 vs. 29.18 in ARLPC and 10.83 in ARPC) and stronger NDR signals were observed (coverage reduction of 2.26 vs. 0.19 in ARPC and 0.37 in ARLPC). A similar trend was observed for H3K4me1 PTMs in LuCaP PDX lines.

ctDNA composite coverage profiles were analyzed in TFBS to assess nucleosome accessibility, whereby lower normalized centroid (±30 bp window) average coverage across these sites suggests greater nucleosome depletion (Methods). For AR TFBS, the strongest signal of nucleosome depletion was observed in ARPC, as indicated by the lowest average centroid coverage (average 0.64, n=16), compared with moderate signals observed in ARLPC (average 0.88, n=2) and the weakest signals observed in NEPC (average 0.95, n=6) (Figure 12B). Conversely, the combined coverage profile in ASCL1 TFBS showed the strongest nucleosome depletion in NEPC samples (mean central coverage 0.69) compared to ARLPC (0.86) and ARPC (0.88) (Figure 12C). These findings were consistent with differential binding activity by AR and ASCL1 in each phenotype from tumor tissue. Furthermore, the ctDNA coverage pattern of nucleosome depletion in ctDNA resembled nucleosome-fringed NDR with H3K27ac and H3K4me1 peak profiles exemplified when only nucleosome-sized fragments (140bp-200bp) generated by cut-and-run were analyzed (Figure 12A). Collectively, these results suggest that nucleosome depletion in ctDNA AR and ASCL1 binding sites represents active TF binding and regulatory activity in a specific prostate PDX tumor phenotype.

Nucleosome patterns in gene promoters inferred from ctDNA are consistent with transcriptional activity of phenotype-specific genes Forty-seven genes, including 12 previously established ARPC lineage markers and 35 NEPC lineage markers (Beltran et al. (2016). Divergent clonal evolution of castration-resistant neuroendocrine prostate cancer. Nature Medicine 22, 298-305; Labrecque et al. (2021b). RNA Splicing Factors SRRM3 and SRRM4 Distinguish Molecular Phenotypes of Castration-Resistant Neuroendocrine Prostate Cancer. Cancer Research 81, 4736-4750), were selected and confirmed by differential expression analysis from PDX tumor RNA-Seq data (FIG. 12D, Methods). To assess the activity of these genes from ctDNA, the TSS (±1 kb window) and ctDNA fragment size within the transcribed region of the genes were analyzed. We found that differential size variability between phenotypes correlated positively with relative expression levels (Spearman's ρ=0.844, p=9.4×10 ⁻¹⁴ , FIG. 12E , Methods). We next analyzed relative ctDNA coverage in the TSS (±1 kb) and observed no association between phenotypes. However, closer examination of ctDNA coverage patterns at promoters revealed consistent nucleosome organization for transcriptional activation and repression (Jiang and Zhang(2021). On the role of transcription in positioning nucleosomes. PLOS Computational Biology 17, e1008556; Klemm et al. (2019). Chromatin accessibility and the regulatory epigenome. Nature Reviews Genetics 20, 207-220; Oruba et al. (2020). Role of cell-type specific nucleosome positioning in inducible activation of mammalian promoters. Nat Commun 11, 1075; Ramachandran et al. (2017). Transcription and Remodeling Produce Asymmetrically Unwrapped Nucleosomal Intermediates. Molecular Cell 68, 1038-1053.e4) (Figure 12D). Therefore, genes were grouped based on differential signals of H3K27me3 histone PTMs associated with transcriptional repression or heterochromatin condensation. For 25 genes (group 1) that do not have differential H3K27me3 peaks, including AR, FOXA1, KLK3 and ASCL1, for example, nucleosome depletion was observed in TSSs for AR in ARPC samples (average coverage 0.47, n=16) and ASCL1 in NEPC samples (0.30, n=6), consistent with the presence of active PTMs (Figure 12F). In contrast, increased coverage was observed in the TSSs of AR in NEPC (1.08) and ASCL1 in ARPC (0.42), thereby supporting nucleosome depletion and inactive transcription in the absence of PTMs. For the 22 genes with differential H3K27me3 peaks (group 2), including STEAP1, CHGB and SRRM4, a relatively more consistent increase in nucleosome occupancy and phasing was observed for NE-specific genes in the TSS as well as within the transcribed regions of the genes (Figure 12G). This group of neural signaling genes, such as UNC13A and INSM1, had reduced signals of nucleosome positioning, consistent with the heterogeneous ("fuzzy") nucleosome patterns described for actively transcribed genes (Jiang and Pugh(2009). Nucleosome positioning and gene regulation: advances through genomics. Nat Rev Genet 10, 161-172; Lai and Pugh(2017). Understanding nucleosome dynamics and their links to gene expression and DNA replication. Nat Rev Mol Cell Biol 18, 548-562). Interestingly, UNC13A was repressed in ARPC tumors, but despite being expressed in NEPC tumors, did not carry the H3K27ac or H3K4me1 accessible PTM marks. These results illustrate that ctDNA analysis can reveal patterns consistent with transcriptional regulation by histone modifications for key genes that define prostate cancer phenotype.

TF activity inferred from analysis of nucleosome accessibility in ctDNA TFBS identifies key regulators of tumor phenotype To characterize lineage-defining TFs with respect to prostate tumor phenotype, nucleosome accessibility in the TFBS of PDX ctDNA was considered. 107 TFs were identified based on the intersection of 338 TFs analyzed using Griffin with 404 TFs differentially expressed between ARPC and NEPC PDX tumors (Methods). Of these TFs, 38 had accessibility in ctDNA that was significantly different between ARPC and NEPC phenotypes (2-tailed Mann-Whitney U test, Benjamini-Hochberg adjusted p<0.05). By unsupervised hierarchical clustering of the combined TFBS central coverage values for the 107 TFs, distinct groups of TFs were observed in PDX ctDNA (FIG. 13). REST had the greatest difference in accessibility as evidenced by reduced coverage within the ARPC model compared to NEPC (log ₂ fold change -0.77, adjusted p= ^5.7x10-4 ). FOXA1, and GRHL2 were significantly more accessible in ARPC (and ARLPC) samples compared to NEPC (log ₂ fold change <-0.57, adjusted p< ^1.3x10-3 ). AR, HOXB13, and NKX3-1 were more accessible in ARPC compared to NEPC (log ₂ fold change <-0.37, adjusted p< ^1.3x10-3 ), but as expected, were only moderately accessible in ARLPC. Interestingly, progesterone receptor (PGR) accessibility was also higher in ARPC (log ₂ fold change -0.33, adjusted p= ^2.6x10-3 ). A group of genes regulated in ARPC that followed a similar trend was also observed, including the glucocorticoid receptor (NR3C1) and other nuclear hormone receptors (NR2F2, RARG), pioneer factors GATA2 and GATA3, and nuclear factors HNF4G and HNF1A (log ₂ fold change <-0.10, adjusted p<0.027).

For factors with increased accessibility in NEPC models compared to ARPC and ARLPC, ASCL1 had the largest TFBS coverage difference (log ₂ fold change 0.36, adjusted p= ^5.7x10-4 , Fig. 12C, Fig. 13F). Other TFs, including RUNX1, BCL11B, POU3F2, NEUROG2, and SOX2, were also more active in NEPC (log ₂ fold change >0.06, adjusted p<0.048), but the differences were modest. HEY1, IRF1, and IKZF1 had a similar trend consistent with increased accessibility in NEPC samples, but were not significantly different from ARPC (adjusted p>0.10). NKX2-1 and CEBPA were increased in accessibility in NEPC compared to ARPC (albeit not significant at adjusted p=0.47 and 0.36, respectively), but these factors were also moderately active in ARLPC. Other notable factors, such as MYC and ETS transcription family genes (ETV4, ETV5, ETS1, ETV1), were highly accessible across all phenotypes, whereas NEUROD1, RUNX3, and TP63 were inaccessible in nearly all samples. Overall, we identified accessibility of known prostate cancer regulators, including ASCL1, NR3C1, HNF4G, HNF1A, and SOX2, not previously demonstrated by ctDNA analysis for these tumor phenotypes (Arora et al. (2013). Glucocorticoid Receptor Confers Resistance to Antiandrogens by Bypassing Androgen Receptor Blockade. Cell 155, 1309-1322; Mu et al. (2017). SOX2 promotes lineage plasticity and antiandrogen resistance in TP53- and RB1-deficient prostate cancer. Science 355, 84-88; Shukla et al. (2017). Aberrant Activation of a Gastrointestinal Transcriptional Circuit in Prostate Cancer Mediates Castration Resistance. Cancer Cell 32, 792-806.e7).

Phenotype-specific open chromatin regions in PDX tumor tissues are reflected in ctDNA profiles of nucleosome accessibility Nucleosome profiling from cfDNA sequencing analysis has been shown to be consistent with global chromatin accessibility in tumor tissues (Snyder et al. (2016). Cell-free DNA Comprises an in Vivo Nucleosome Footprint that Informs Its Tissues-Of-Origin. Cell 164, 57-68; Sun et al. (2019). Orientation-aware plasma cell-free DNA fragmentation analysis in open chromatin regions informs tissue of origin. Genome Research 29, 418-427; Ulz et al. (2019). Inference of transcription factor binding from cell-free DNA enables tumor subtype prediction and early detection. Nature Communications 10, 4666), but its application for distinguishing tumor phenotypes is limited. To provide information on phenotypic differences in chromatin accessibility, we investigated the use of ATAC-Seq data from tumor tissues for 10 LuCaP PDX lines (5 ARPC and 5 NEPC) (Cejas et al. (2021). Subtype heterogeneity and epigenetic convergence in neuroendocrine prostate cancer. Nat Commun 12, 5775). We defined an initial set of 28,765 ARPC and 21,963 NEPC differential consensus open chromatin regions, which we further restricted to overlapping TFBS for 338 TFs, resulting in 15,881 ARPC and 11,694 NEPC sites (Methods, Fig. 14A). For ARPC-specific open chromatin sites, reduced overall complex site coverage (window of +/- 1 kb) and central coverage (+/- 30 bp) was observed in ctDNA for ARPC PDX lines (mean central coverage 0.75, n=16) compared to NEPC lines (mean 0.96, n=6) and cfDNA from healthy human donors (mean 0.97, n=14) (Figure 14B). Conversely, for NEPC-specific open chromatin sites, coverage in ctDNA was reduced for NEPC lines (mean 0.89) compared to ARPC lines (mean 1.01) and healthy donor cfDNA (mean 1.00) (Figure 14C). These results confirmed that tumor tissue chromatin accessibility can be confirmed with ctDNA and that ARPC and NEPC phenotypes have distinct ctDNA complex site coverage profiles.

Comprehensive evaluation of ctDNA features across genomic contexts for CRPC phenotyping To assess the utility of ctDNA nucleosome profiling to inform prostate cancer phenotyping, a set of comprehensive genome-wide ctDNA features was systematically evaluated: fragment size, local coverage profiling, and composite site coverage profiling (FIG. 11A). By principal component analysis (PCA), distinct feature signals were observed between ARPC and NEPC phenotypes for composite TFBS coverage of TFs and fragment size variability at global sites of PTMs (FIG. 14D, Methods). In addition to these features, previously reported similar methods also included the ratio of short to long fragments and the regional coverage pattern in the TSS (maximum peak between -120 bp and 195 bp) (Cristiano et al. (2019). Genome-wide cell-free DNA fragmentation in patients with cancer. Nature 570, 385-389; Ulz et al. (2016b). Inferring expressed genes by whole-genome sequencing of plasma DNA. Nature Genetics 48, 1273-1278) (Methods).

All combinations of coverage and fragment size features were quantitatively evaluated for different genomic contexts to explore their potential for classifying ARPC and NEPC phenotypes. For each feature set, 100 iterations of stratified cross-validation were performed using a supervised machine learning classifier (XGBoost) on ctDNA samples from 16 ARPC and 6 NEPC models, and the area under the receiver operating characteristic curve (AUC) was computed (Methods). First, an established set of 10 genes associated with AR activity was evaluated (Bluemn et al. (2017). Androgen Receptor Pathway-Independent Prostate Cancer Is Sustained through FGF Signaling. Cancer Cell 32, 474-489.e6; Labrecque (2021a). The heterogeneity of prostate cancers lacking AR activity will require diverse treatment approaches. Endocrine-Related Cancer 28, T51-T66). We observed that graded nucleosome distance at H3K27ac sites and central coverage at TSS had moderate predictive performance (AUC 0.88). When considering all PTM sites, promoters, genes, TFs, and open chromatin regions, the best-performing features included the average fragment size feature at H3K4me1 sites (n = 9,750, AUC 1.0) and the promoter TSS feature (n = 17,946, AUC 1.0), and both open chromatin complex site features (AUC 1.0) (Figure 14E).

PDX ctDNA analysis informs accurate classification of ARPC and NEPC phenotypes from patient plasma using a probabilistic model An important consideration and challenge in analyzing plasma from patients is the presence of cfDNA released by hematopoietic cells, which reduces the ctDNA fraction (i.e., tumor fraction). Furthermore, the small patient cohorts for which information on tumor phenotype is available make supervised machine learning approaches suboptimal. Therefore, a probabilistic model was developed to estimate the proportion of ARPC and NEPC from individual plasma samples, taking into account the tumor fraction (Methods). PDX plasma ctDNA signals were used to inform the model, focusing on phenotype-specific open chromatin complex site features (Figures 14B and 14C). The model produces a normalized prediction score that represents the estimated signature of ARPC (low value) and NEPC (high value). This method was applied to a benchmarking dataset generated by simulating various tumor fractions and sequencing coverages using five ARPC PDX ctDNA and five NEPC PDX ctDNA samples, respectively (FIG. 14F, Methods). An AUC of 1.0 was achieved at 25× coverage down to a tumor fraction of 0.01, at 1× down to a tumor fraction of 0.2, and at 0.2× coverage, an AUC of 1.0 was achieved at a tumor fraction of 0.3, suggesting a potential upper bound performance for classification of samples with small tumor fractions in plasma (FIG. 14G).

To test the classification performance of this model on patient samples, a published dataset of ultra-low-pass whole genome sequencing (ULP-WGS) of plasma cfDNA from 101 mCRPC patients (DFCI cohort I), including 80 adenocarcinoma (ARPC) and 21 NEPC samples (mean coverage 0.52×, range 0.28-0.92×) was analyzed (Berchuck et al. (2022). Detecting Neuroendocrine Prostate Cancer Through Tissue-Informed Cell-Free DNA Methylation Analysis. Clinical Cancer Research 28, 928-938). Using a model that is unsupervised and uses parameters informed only by PDX analysis, an overall AUC of 0.96 was achieved (Figure 15A). When considering samples with large (≧0.1) and small (<0.1) tumor fractions, the model had AUCs of 0.97 and 0.76, respectively. The best overall performance was identified with 97.5% specificity (ARPC) and 90.4% sensitivity (NEPC), corresponding to a prediction score of 0.3314 (Figure 15A). In another published dataset of 11 mCRPC samples from 6 patients (DFCI cohort II) with high PSA and/or treated with ARSI (Choudhury et al. (2018). Tumor fraction in cell-free DNA as a biomarker in prostate cancer. JCI Insight 3; Viswanathan et al. , 2018), the model correctly classified patients as ARPC for 11 (100%) WGS (~20x) and 8 (73%) ULP-WGS (~0.1x) samples when using the best score cutoff (Figure 15B).

Next, we analyzed 61 clinical plasma samples from 30 CRPC patients with ARPC, NEPC, and mixed phenotypes representing typical clinical histories. ULP-WGS of cfDNA was performed and 47 samples from 30 patients (26 ARPC, 5 NEPC, and 16 mixed phenotypes) were selected for deeper WGS (mean 22.13× coverage, range 15.15× to 31.79×) based on tumor fraction and AR copy number status (Methods). For the 26 samples with ARPC clinical phenotype, all were predicted to be primarily ARPC using a score cutoff of 0.3314 (Figure 15C). For the NEPC clinical phenotype, all five were predicted to be NEPC with scores above the cutoff. We also found a negative association between patient ctDNA coverage and tumor fraction at open chromatin sites for both ARPC (Spearman's ρ = -0.93) and NEPC (Spearman's ρ = -1.00) predictions, suggesting that the observed ctDNA signal may be tumor specific. The ULP-WGS data correctly predicted 22 samples (84%) with ARPC clinical phenotype and all 5 samples (100%) with NEPC clinical phenotype (Figure 15C). The remaining 16 samples had clinical history or tumor histology reflecting mixed phenotypes, e.g., mixed NEPC and AR-positive adenocarcinoma tumors (Figure 15C). Of the 12 samples that included the presence of ARPC in the mixed clinical phenotype, 10 (83%) were classified as ARPC at the optimal score cutoff. All 3 samples with NEPC but not ARPC in the clinical phenotype were classified as NEPC by the model. Overall, 100% accuracy was achieved for WGS data in samples with unambiguous clinical phenotypes (87% for ULP-WGS). However, predictions for mixed or ambiguous phenotypes were variable, highlighting the complexities associated with classification for patients with advanced prostate cancer, where tumor heterogeneity can be observed.

Discussion The study presented here is believed to be the largest sequencing study of human ctDNA from mouse plasma in PDX models to date. Mouse plasma sequencing provides a unique opportunity to comprehensively investigate epigenetic nucleosome patterns in ctDNA from well-characterized tumor models. A computational methodology was developed and applied to construct a large number of ctDNA features, each of which was linked to transcriptional regulation across CRPC tumor phenotypes in the LuCaP PDX model. Using the features learned from PDX ctDNA, a probabilistic model was developed to accurately classify ARPC and NEPC phenotypes from patient plasma in three clinical cohorts.

The use of PDX mouse plasma overcomes the challenges of low ctDNA content or incomplete tumor knowledge when investigating patient samples and may facilitate cfDNA diagnostic development, basic cancer research, and clinical translation. Furthermore, the LuCaP ctDNA sequencing data complements the characterization of CRPC tumor phenotypes from mature tissues. In addition to supporting molecular studies of CRPC, the ctDNA data and disclosed approaches detail the potential utility of the PDX model for translational research. Although these data focus on ARPC and NEPC phenotypes, this study may serve as a framework for using PDX plasma from additional CRPC phenotypes and other cancer models.

Analysis of LuCaP PDX ctDNA sequencing data confirmed the activity of key regulators, including a set of 47 established differentially expressed gene markers, between ARPC and NEPC phenotypes. While inference of gene expression from ctDNA has been shown in proof-of-concept studies (Ulz et al. (2016b). Inferring expressed genes by whole-genome sequencing of plasma DNA. Nature Genetics 48, 1273-1278; Zhu et al. (2021). Tissue-specific cell-free DNA degradation quantifies circulating tumor DNA burden. Nature Communications 12, 2229), PDX ctDNA allows detailed analysis of nucleosome organization associated with transcriptional activity of individual genes that define tumor phenotype.

In addition to the existing molecular profiling available for these models, this study now provides characterization of histone PTMs in LuCaP PDX tumors using cut-and-run. In regions harboring these PTMs on the histone tails, predicted nucleosome patterns in ctDNA were observed, consistent with active or repressed gene transcription. We believe this is the first time that ctDNA analysis has been performed on histone PTMs, providing a blueprint for developing new methods to investigate additional epigenetic changes using PDX plasma.

While ctDNA has demonstrated regulation of key factors such as AR, HOXB13, NKX-3.1, FOXA1, and REST in CRPC (Ulz et al. (2019). Inference of transcription factor binding from cell-free DNA enables tumor subtype prediction and early detection. Nature Communications 10, 4666), this study is the first to reveal differential activity of other key factors in CRPC from ctDNA analysis. These key factors included the glucocorticoid receptor (NR3C1), nuclear factors HNF4G and HNF1A, and pioneering factors GATA2 and GATA3, all of which are associated with adenocarcinoma of the prostate (ARPC) (Arora et al. (2013). Glucocorticoid Receptor Confers Resistance to Antiandrogens by Bypassing Androgen Receptor Blockade. Cell 155, 1309-1322; Chaytor et al. , 2019; Shukla et al. , 2017). ASCL1 is a pioneer TF that has a role in neuronal differentiation and has recently been described to be active during NE transdifferentiation and in NEPC (Cejas et al. , 2021; Rapa et al. , 2008). To our knowledge, this study is the first to demonstrate ASCL1 binding site accessibility in NEPC from plasma ctDNA and provide a detailed characterization of its transcriptional activity.

This study provides a comprehensive analysis of the TFBS of 338 factors in each plasma sample without the need for chromatin immunoprecipitation or other epigenetic assays. However, no significant differences in accessibility were observed for 69 of the 107 TFs in ctDNA, which may be consistent with TF activity not necessarily correlating with its own expression level (Corces et al. (2018). The chromatin accessibility landscape of primary human cancers. Science 362). On the other hand, TFBS accessibility may not necessarily indicate true TF activity, such as the binding of multiple factors to the same locus. Moreover, this analysis was based on TFBS obtained from public databases. However, a prostate phenotype-specific TF cistrome could better guide this approach.

State-of-the-art computational methods based on existing and new concepts in ctDNA data analysis were applied to extract tumor-specific features. While other approaches have considered regions such as TSS, TFBS, and DNase hypersensitive sites (Peneder et al. (2021). Multimodal analysis of cell-free DNA whole-genome sequencing for pediatric cancers with low mutational burden. Nat Commun 12, 3230; Snyder(2016). Cell-free DNA Comprises an in Vivo Nucleosome Footprint that Informs Its Tissues-Of-Origin. Cell 164, 57-68; Ulz et al. (2016b). Inferring expressed genes by whole-genome sequencing of plasma DNA. Nature Genetics 48, 1273-1278; Ulz et al. (2019). Inference of transcription factor binding from cell-free DNA enables tumor subtype prediction and early detection. Nature Communications 10, 4666), after systematic evaluation, we focused on ctDNA signatures in open chromatin sites derived from ATAC-Seq of PDX tissues (Cejas et al. (2021). Subtype heterogeneity and epigenetic convergence in neuroendocrine prostate cancer. Nat Commun 12, 5775) was found to provide the best performance for distinguishing CRPC phenotypes. An unsupervised probabilistic model is presented that estimates the proportion of ARPC and NEPC in patient plasma using a statistical framework informed by ideal parameters from LuCaP PDX ctDNA analysis. The model does not require training on patient samples, but does require estimates of tumor fraction (ichorCNA (Adalsteinsson(2017). Scalable whole-exome sequencing of cell-free DNA reveals high concordance with metastatic tumors. Nature Communications 8) and prediction score cutoffs, determined from DFCI cohort I. The framework presented here can be extended to model multiple phenotypic classes, provided that parameters that inform these additional situations can be learned. Insights from additional datasets, such as single-cell nucleosome and accessibility profiling of PDX tumors and clinical samples (Fang et al. (2021). Comprehensive analysis of single cell ATAC-seq data with SnapATAC. Nat Commun 12, 1337; Wu et al. (2021). Single-cell CUT & Tag analysis of chromatin modifications in differentiation and tumor progression. Nat Biotechnol 39, 819-824), can improve the resolution of ctDNA analysis.

By applying the predictive model to the patient dataset with the most reliable clinical phenotype, high performance was obtained despite the use of low-depth sequencing. Notably, the performance for the DFCI cohort I was also consistent with reported phenotyping results in which ctDNA methylation in the same patients was used (Berchuck et al. (2022). Detecting Neuroendocrine Prostate Cancer Through Tissue-Informed Cell-Free DNA Methylation Analysis. Clinical Cancer Research 28, 928-938). Similarly, for the UW cohort, samples with well-defined clinical phenotypes had perfect concordance from the deep WGS data. However, for samples with mixed or unclear clinical phenotypes, a subset of cases had complex clinical and histopathological features, limiting the ability to conclusively assess the model's performance. Tumor heterogeneity and the coexistence of different molecular phenotypes are common in mCRPC, where treatment-induced phenotypic plasticity can vary within and between tumors of individual patients. Future extension of the model to predict mixed phenotypes from ctDNA requires larger studies with comprehensive assessment of tumor histology.

In summary, this study illustrates for the first time that large-scale analysis of ctDNA from PDX mouse plasma can facilitate more detailed investigation of tumor regulation. These results, together with the suite of computational methods presented herein, highlight the utility of ctDNA to investigate the transcriptional regulation of tumor phenotype and its potential diagnostic applications in cancer precision medicine.

Experimental Model and Subject Details PDX Mouse Model LuCaP patient-derived xenograft tumors (established at the University of Washington) were initiated from tumor specimens resected from men with advanced prostate cancer. The establishment and characterization of the PDX model was previously described (Lam et al. (2018). Generation of Prostate Cancer Patient-Derived Xenografts to Investigate Mechanisms of Novel Treatments and Treatment Resistance. In Prostate Cancer: Methods and Protocols, Z. Culig, ed. (New York, NY: Springer), pp. 1-27). PDX were grown in vivo in male NOD scid IL2R-gamma-null (NSG) mice (cat#005557) from Jackson Labs. Tumor harvesting to establish PDX lines was approved by the University of Washington Human Subjects Division IRB (IRB#2341) after written informed consent was obtained from the patients. Up to five mice were caged in a pathogen-free facility, provided with food and water ad libitum, and maintained on a 12-h light/dark cycle. Surgery was performed under isoflurane anesthesia, and mice were given supplemental buprenorphine sustained release (SR). PDX lines were assessed by at least two expert pathologists using histopathology and histologic phenotypic subtype annotations were orthogonally validated based on transcriptome-derived signature marker expression scores to define phenotypes (Beltran et al. (2016). Divergent clonal evolution of castration-resistant neuroendocrine prostate cancer. Nature Medicine 22, 298-305; Bluemn(2017). Androgen Receptor Pathway-Independent Prostate Cancer Is Sustained through FGF Signaling. Cancer Cell 32, 474-489.e6; Nyquist et al. (2020). Combined TP53 and RB1 Loss Promotes Prostate Cancer Resistance to a Spectrum of Therapeutics and Confers Vulnerability to Replication Stress. Cell Reports 31, 107669): Adenocarcinoma AR positive (ARPC), neuroendocrine positive (NEPC), and AR-low, neuroendocrine negative (ARLPC). Animals were sacrificed and removed from the study if the tumor began to ulcerate during LuCaP growth or if the animal's health was compromised. Resected PDX tumors (300-800 mm ³ ) were divided into pieces of approximately 50 mg to approximately 100 mg and stored at -80°C. Animal studies were approved by the Fred Hutchinson Cancer Research Center (FHCRC) IACUC (protocol 1618) and performed in accordance with NIH guidelines. For this study, blood was collected by cardiac puncture from animals bearing PDX tumors (measurable size 300-800 mm ³ ).

Human Subjects UW Cohort: Blood samples were collected from men with metastatic castration-resistant prostate cancer at the University of Washington (collected between 2014 and 2021 under University of Washington Human Subjects Division IRB protocol number CC6932). Sixty-one plasma samples from 30 patients were analyzed in this study. After initial ultra-low-pass whole genome sequencing (ULP-WGS) analysis, 47 plasma samples from 30 patients were retained for further high-depth coverage whole genome sequencing (WGS) analysis. All samples were de-identified before ctDNA analysis, and a double-blind approach was used to evaluate clinical phenotype predictions. Initial patient selection was based on information on clinical disease burden and availability of clinically derived phenotype subtype annotations. Clinical information regarding these patients was protected by IRB protocol restrictions.

DFCI Cohort I: Plasma was collected from men diagnosed with mCRPC and treated at Dana-Farber Cancer Institute (DFCI), Brigham and Women's Hospital, or Weill Cornell Medicine (WCM) between April 2003 and August 2021. All patients provided written informed consent to participate in the study and for genomic analysis of biospecimens and blood. Use of samples was approved by the DFCI IRB (#01-045 and 09-171) and the WCM (1305013903) IRBs. ULP-WGS data with a mean coverage of 0.5x (range 0.3x-0.9x) for 101 patients was previously published (Berchuck et al. (2022). Detecting Neuroendocrine Prostate Cancer Through Tissue-Informed Cell-Free DNA Methylation Analysis. Clinical Cancer Research 28, 928-938). The presence of high-grade neuroendocrine carcinoma of prostatic origin was confirmed by two genitourinary pathologists according to contemporary practice, based on histological review of available material, reinterpretation of the original reports, and incorporation of available molecular results (Epstein et al. (2014). Proposed Morphologic Classification of Prostate Cancer With Neuroendocrine Differentiation. The American Journal of Surgical Pathology 38, 756-767). Patients with ARPC (clinically annotated as PRAD) had castration-resistant prostate adenocarcinoma without pathologic evidence of neuroendocrine differentiation throughout the disease course.

DFCI Cohort II: Plasma samples in this cohort were collected from men diagnosed with and treated for mCRPC at the Dana-Farber Cancer Institute (DFCI). All patients provided written informed consent for research blood collection and analysis of clinical and genetic data (DFCI protocols #01-045 and 11-104). WGS data with an average coverage of 27× (range 11× to 44×) (Viswanathan et al. (2018). Structural Alterations Driving Castration-Resistant Prostate Cancer Revealed by Linked-Read Genome Sequencing. Cell 174, 433-447.e19) and ULP-WGS data with an average coverage of 0.13× (range 0.07× to 0.18×) (Adalsteinsson et al. (2017). Scalable whole-exome sequencing of cell-free DNA reveals high concordance with metastatic tumors. Nature Communications 8; Choudhury et al. (2018). Tumor fraction in cell-free DNA as a biomarker in prostate cancer. JCI Insight 3) were downloaded from dbGAP accession phs001417. Eleven samples from six patients had WGS and ULP-WGS matching with paired-end reads required for analysis by Griffin. Prostate-specific antigen (PSA, ng/mL) levels at time of blood collection and treatments have been previously published (Choudhury et al. (2018). JCI Insight 3). Six patients had been treated with abiraterone, enzalutamide, or bicalutamide for adenocarcinoma or had detectable levels of PSA.

Healthy donor plasma cfDNA WGS data used in this study were obtained from a previously published study. Two samples (HD45 and HD46) with 13× and 15× coverage, respectively, were obtained from dbGAP under accession phs001417 (Adalsteinsson et al. (2017). Nature Communications 8; Viswanathan et al. (2018). Cell 174, 433-447.e19). Informed consent was obtained from these donors under DFCI protocol IRB (#03-022). Thirteen healthy donor plasma cfDNA WGS data with coverage ranging from 13.5× to 27.6× (12 males: NPH002, 03, 06, 07, 12, 18, 23, 26, 33, 34, 35, 36; 1 female (mixed): NPH004) were obtained from the European Phenome Archive (EGA) under accession EGAD00001005343 (Ulz et al. (2019). Nature Communications 10, 4666).

Method Details PDX Plasma Processing Blood samples were collected from NSG mice bearing subcutaneous PDX tumors at the time of sacrifice. PDX strains were maintained in the animal care facilities at the University of Washington and FHCRC. Blood was processed for subsequent DNA isolation as described for human plasma DNA processing. Blood was collected in purple-capped EDTA tubes and processed within 4 hours. All blood samples were double spun using centrifugation at room temperature at 2500g for 10 minutes followed by a 16000g spin for 10 minutes for the plasma fraction. Seven to ten mouse plasma samples were pooled for each PDX strain. Processed plasma samples were stored in clean screw-cap cryo-microcentrifuge tubes and stored at -80°C until cfDNA isolation.

Cell-Free DNA Isolation cfDNA was isolated from plasma from PDX mice using the QIAamp Circulating Nucleic Acid Kit using the recommended protocol. Pooled plasma samples from 7-10 mice for each PDX line contained approximately 2-3 mL total plasma volume for each line. The cfDNA kit method based on filter retention does not perform any fragment size class enrichment. Carrier RNA spike-in was removed from the elution buffer. Isolated cfDNA was quantified using the Qubit dsDNA HS assay (Invitrogen) and cfDNA fragment size profiles were analyzed using Tapestation HS D5000 and HS D1000 assays (Agilent).

Cell-free DNA library preparation and sequencing NGS libraries were prepared for LuCaP PDX mouse plasma samples using 50 ng of input cfDNA. Illumina NGS sequencing libraries were prepared using the KAPA hyperprep kit, employing 9 cycles of amplification and purification using laboratory standardized SPRI beads. KAPA UDI dual index library adapters were used. Library concentrations were equilibrated, pooled for multiplexing, and sequenced at Fred Hutch Genomics Shared Resources using an Illumina HiSeq 2500 (200 cycles) and at Broad Institute Genomics Platform Walkup-Seq Services using an Illumina NovaSeq platform and S4 flow cell (300 cycles). Truncated 200-cycle FASTQ files were generated (100 bp paired-end reads) to match the Illumina HiSeq 2500 data.

Clinical patient plasma samples collected at the University of Washington (UW cohort) were submitted to Broad Institute Blood Biopsy Services. Briefly, cfDNA was extracted from 2 mL of double-spun plasma and ultra-low-pass whole genome sequencing (ULP-WGS) was performed to approximately 0.2× coverage. The ichorcNA pipeline was used to estimate tumor DNA content (i.e., tumor fraction, see below). Forty-seven samples (from 30 patients) had either ≥5% tumor fraction or ≥2% tumor fraction, AR amplification was observed in ichorcna, and were subsequently sequenced to greater depth WGS coverage (approximately 20×).

Cell-free DNA sequencing analysis and mouse subtraction All cfDNA sequencing data used in this study were realigned to the hg38 human reference genome (hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz). FASTQ files were realigned using BWA mem (Li (2013). Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. ArXiv:1303.3997 [q-Bio]), and post-alignment processing was performed according to the GATK Best Practices workflow (DePristo et al. (2011). A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nat Genet 43, 491-498).

For PDX ctDNA whole genome sequence data, mouse genome subtraction was performed following a previously described protocol (Jo et al. (2019). Impact of mouse contamination in genomic profiling of patient-derived models and best practice for robust analysis. Genome Biology 20, 231). For mouse genome subtraction, reads were aligned to a concatenated reference consisting of both the human reference genome (hg38) and the mouse reference genome (mm10, GRCm38.p6, igenomes.illumina.com.s3-website-us-east-1.amazonaws.com/Mus_musculus/NCBI/GRCm38/Mus_musculus_NCBI_GRCm38.tar.gz) using BWA mem. Read pairs where both reads aligned to the human reference genome were kept and all other read pairs were removed. The remaining reads were then realigned to the human-only reference. Finally, the GATK best practices workflow was applied to each sample. After mouse subtraction, samples with a depth below 3x were removed from downstream analysis. The mouse subtraction pipeline used in this study is available atgithub.com/GavinHaLab/PDX_mouseSubtraction.

Differential mRNA Expression Analysis RNA isolation of 102 tumors from 46 LuCaP PDX samples was performed as previously described (Labrecque (2019). Molecular profiling stratifies diverse phenotypes of treatment-refractory metastatic castration-resistant prostate cancer. J Clin Invest 129, 4492-4505). RNA concentration, purity, and integrity were assessed by NanoDrop (Thermo Fisher Scientific Inc) and Agilent TapeStation, and RNA RIN ≥ 8 was retained for library preparation. RNA-Seq libraries were constructed from 1 μg of total RNA using the Illumina TruSeq Stranded mRNA LT Sample Prep Kit according to the manufacturer's protocol. Barcoded libraries were pooled and sequenced by Illumina NovaSeq 6000 or Illumina HiSeq 2500 to generate 50 bp paired-end reads. Sequencing reads were mapped to the hg38 human genome and mm10 mouse genome using STAR. v2.7.3a (Dobin et al. (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15-21). All subsequent analyses were performed in R-4.1.0. Sequences aligned to the mouse genome and therefore derived from potential contamination of mouse tissues were removed from the analysis using XenofilteR (v1.6) (Kluin (2018). XenofilteR: computational deconvolution of mouse and human reads in tumor xenograft sequence data. BMC Bioinformatics 19, 366). Gene-level abundance was quantified using the R package GenomicAlignments summarizeOverlaps function, using mode=IntersectionStrict and restricting to primary aligned reads. refSeq gene annotations were used for transcriptome analysis. Transcript abundances were input into edgeR (Robinson et al. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139-140) and filtered for minimum expression levels using the filterByExpr function with default parameters, then limma voom was used for differential expression analysis of NEPC vs. ARPC and ARLPC vs. ARPC. Results were then filtered using a previously published list of 1,635 human transcription factors (Lambert et al. (2018). The Human Transcription Factors. Cell 172, 650-665), resulting in 514 genes with FDR<0.05 and fold change>3. Of these 514, deregulation of gene expression for 404 transcription factor genes delineated ARPC from NEPC.

On-target cleavage and release using nucleases (cut and run)
Cut-and-Run is an antibody-targeted enzyme-tethering chromatin profiling assay that releases specific protein-DNA complexes into the supernatant by controlled cleavage with micrococcal nuclease for paired-end DNA sequencing analysis. Cut-and-Run assays were performed for three histone modifications, H3K27ac, H3K4me1, and H3K27me3, according to published protocols (Skene and Henikoff(2017). An efficient targeted nuclease strategy for high-resolution mapping of DNA binding sites. ELife 6, e21856). Cut-and-Run was performed on LuCaP PDX tumors using approximately 75 mg of flash-frozen tissue pieces. Briefly, frozen tissue was thoroughly minced into small pieces and converted to smaller cell clumps using collagenase and dispase. Cell clumps were permeabilized using digitonin and rocked with target antibodies in EDTA antibody buffer. A time-sensitive micrococcal nuclease enzyme treatment was performed on ice. The released DNA was precipitated with glycogen carriers, followed by NGS library preparation using the picogram input DNA library preparation protocol.

Paired-end (50 bp) sequencing was performed and reads were aligned to the hg38 human reference assembly using bowtie2 version 2.4.2 (Langmead et al. (2019). Scaling read aligners to hundreds of threads on general-purpose processors. Bioinformatics 35, 421-432). Aligned reads were processed as described in the SEACR protocol (github.com/FredHutch/SEACR#preparing-input-bedgraph-files). Peaks were called using SEACR version 1.3 (Meers et al. (2019). Peak calling by Sparse Enrichment Analysis for CUT&RUN chromatin profiling. Epigenetics & Chromatin 12, 42) using the "string" setting and with reference to paired IgG controls. Bam files were inspected, parsed, and filtered using SAMtools (Danecek et al. (2021). Twelve years of SAMtools and BCFtools. GigaScience 10, giab008). Insert size-restricted analyses were performed by retaining fragments smaller than 120 bp (subnucleosome size) or by retaining fragments between 140 and 200 bp (nucleosome size). BigWig files were prepared using bamCoverage in deepTools 3.5.0 (Ramirez et al. (2016). deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Research 44, W160-W165) using a bin size of 10 and the "extendReads" option on the bam files. Genome-wide peak heatmaps, targeted heatmaps, and respective profiles were plotted using deepTools. The mean function of wiggletools 1.2.8. and deepTools computeMatrix were used to obtain bigWig format files for each phenotype. The coordinates of the regions yielding phenotype-specific information were obtained from diffBind v3.5.0, and the top 10,000 most significant regions (all FDR<0.05) that were differentially open between ARPC and NEPC lines were used for downstream feature analysis (see Selection of transcribed and promoter regions of genes for additional subsetting criteria applied per feature). The plotHeatmap function was used for heatmaps and profiles. The "Peak Center" option was used to extract the desired heatmaps. All these steps were performed for the H3K27ac, H3K4me1 and H3K27me3 antibodies. The area under the curve (AUC) of the scaled heatmap profiles and the peak height at the profile center were estimated using deepStats v0.4 (Richard, 2020) (equivalent profiles were scaled to 10 units).

Differential histone post-translational modification (PTM) analysis Differential PTM analysis was performed using the Diffbind version 2.16.0 package in R-4.0.1 (Ross-Innes et al. (2012). Differential oestrogen receptor binding is associated with clinical outcome in breast cancer. Nature 481, 389-393) using standard parameters (bioconductor.riken.jp/packages/3.0/bioc/html/DiffBind.html). ARPC, NEPC and ARLPC samples were grouped by phenotype defined by histopathology and transcriptome signatures as described in the "PDX mouse model" section. Samples were loaded using the dba function, reads were counted using the dba.count function and dba.count was used to count reads. Contrast and minimum member 2 were used to identify contrast as phenotype. Differential peak sites were computed using the dba.analyze function with default settings. Differential peak binding of NEPC and ARLPC was computed for ARPC samples. Unique binding sites in NEPC and ARLPC were cataloged using bedtools v2.29.2 (Quinlan and Hall (2010). BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26, 841-842). Differentially bound peaks between groups were analyzed using ChIPseeker 1.28.3 (Yu et al. (2015). ChIPseeker: an R/Bioconductor package for ChIP peak annotation, comparison and visualization. Bioinformatics 31, 2382-2383) and TxDb. Hsapiens. UCSC.hg38. was annotated using knownGene 3.2.2.

ATAC-Seq Analysis ATAC-Seq sequence data for 15 tumor samples from 10 PDX lines have been previously published, with FASTQ files available upon request (Cejas et al. (2021). Subtype heterogeneity and epigenetic convergence in neuroendocrine prostate cancer. Nat Commun 12, 5775). These lines included LuCaP PDX lines with ARPC histology (23.1, 77, 78, 81, 96) and LuCaP PDX lines with NEPC histology (2 replicates each of 49, 93, 145.1, 173.1 and 1 replicate of 145.2). Paired-end reads were aligned using bowtie2 2.4.2 (Langmead et al. (2019). Scaling read aligners to hundreds of threads on general-purpose processors. Bioinformatics 35, 421-432) and aligned against the UCSC hg38 human reference assembly with very-sensitive -k 10 settings. Peaks were called using Genrich version 0.6.1 (github.com/jsh58/Genrich). Differential binding analysis was performed using the Diffbind version 3.5.0 package in R version 4.1.0. hg38-blacklist. v2 (Amemiya et al. (2019). The ENCODE Blacklist: Identification of Problematic Regions of the Genome. Sci Rep 9, 9354) (github.com/Boyle-Lab/Blacklist) was used to exclude ENCODE blacklisted regions. AR-positive NE-null PDX samples (n=5) were compared to AR-null NE-positive PDX samples (n=5) using RNA-Seq derived phenotypes. Phenotype-specific binding sites were first selected for positive fold change open chromatin enrichment and then isolated by using Intervene 0.6.5 (Khan and Mathelier (2017). Intervene: a tool for intersection and visualization of multiple gene or genomic region sets. BMC Bioinformatics 18, 287). Regions were considered to overlap if they shared at least 1 bp. Regions with an FDR-adjusted p-value <0.05 were then subset to those overlapping at least 1 bp with 338,000 established TFBSs (338 TFs x 1,000 binding sites, see Griffin analysis for site selection) using BedTools Intersect. Only regions overlapping with established TFBSs were retained.

Griffin Analysis Griffin is a method for profiling nucleosome protection and accessibility for predefined genomic loci (see Example 1 and Doebley et al. (2021). Griffin: Framework for clinical cancer subtyping from nucleosome profiling of cell-free DNA. MedRxiv 2021.08.31.21262867). Griffin filters sites by mappability, estimates and corrects for GC bias at a per-fragment level, and generates a GC-corrected coverage profile around each site. First, griffin takes a list of sites and examines the mappability within a window (+/- 5000 bp around each site). Mappability (hg38 Umap multi-read mappability for 50 bp reads) was obtained from the UCSC genome browser (Karimzadeh et al. (2018). Umap and Bismap: quantifying genome and methylome mappability. Nucleic Acids Research 46, e120) (hgdownload.soe.ucsc.edu/gbdb/hg38/hoffmanMappability/k50.Umap.MultiTrackMappability.bw). Sites with mappability below 0.95 were excluded from further analysis. Next, for each sample, GC bias was quantified using a modified version of a previously described method (Benjamini and Speed (2012). Summarizing and correcting the GC content bias in high-throughput sequencing. Nucleic Acids Research 40, e72-e72). Briefly, for each possible fragment length and GC content, the number of reads in the bam file and the number of positions in the genome with that particular length and GC content were counted. The GC bias and GC content for each fragment length was calculated by dividing the number of observed reads by the number of observed positions in the genome of that fragment length and GC content. The GC bias for all possible GC contents at a given fragment length was then normalized to a mean bias of 1. The GC bias was then smoothed by taking the median of fragments of similar length and GC content (k-nearest neighbor smoothing) to generate a smoothed GC bias value.

Nucleosome profiling was performed for each sample after GC correction. For each mappable site of interest, fragments that aligned to a region of ±5000 bp from the site were fetched from the bam file. Fragments were filtered to remove duplicates and low quality alignments (mapping quality <20) and by fragment length. Nucleosome-sized fragments (140-250 bp) were retained. Fragments were then GC-corrected by assigning each fragment a weight of 1/GC_bias for its given fragment length and GC content to identify fragment midpoints. The number of weighted fragment midpoints within 15 bp bins across the site were counted. For composite sites, all sites of a given type (e.g., all sites for a given transcription factor) were summed to generate a single coverage profile. Individual or composite coverage profiles were normalized to the average coverage of 1 within a region of ±5000 bp surrounding the site. Finally, sites were smoothed using a Savitsky-Golay filter with a window length of 165 bp and polynomial degree 3. A window of ±1000 bp around the site was kept for plotting and feature extraction (see Griffin manuscript for further details); where sites are plotted, shading illustrates the 95% confidence interval within the sample group. Features extracted from individual or combined sites included:
a. "average center coverage", the average coverage between -30 bp and 30 bp relative to the site center;
b. "Average window coverage", the average coverage between -990 bp and 990 bp relative to the site center, and c. "Maximum height", the absolute difference between the minimum coverage in a window between -120 bp and 30 bp relative to the TSS and the maximum coverage in a window between 31 bp and 195 bp relative to the TSS.

Griffin Analysis of Selective Transcription Factor Binding Sites (TFBS) Griffin analysis of transcription factor binding sites (TFBS) was performed using the same list of TFBSs utilized in Griffin (see Example 1 and Doebley et al. (2021). MedRxiv 2021.08.31.21262867). Briefly, TFBS locations were downloaded from the Gene Transcriptional Regulation Database (GTRD) (Yevshin et al. (2019). GTRD: a database on gene transcription regulation-2019 update. Nucleic Acids Res 47, D100-D105), which contains a compilation of ChIP-Seq data from multiple sources. Metacluster data or metapeaks (version 19.10) observed in one or more ChIP-Seq experiments were used. An initial list of 1,314 transcription factors (TFs) in GTRD was compared to the CIS-BP database (Weirauch et al. (2014). Determination and inference of eukaryotic transcription factor sequence specificity. Cell 158, 1431-1443) (v2.00, downloaded from cisbp.ccbr.utoronto.ca/bulk.php). TFs from GTRD that were also present in CIS-BP and had known binding motifs were retained. The selected TF-bound genomic loci were then filtered for mappability as described above (Griffin analysis) and TFs with fewer than 10,000 highly mappable sites on autosomes were removed, resulting in 338 TFs. For each TF, sites were sorted by "peak count" value from GTRD, and the 1,000 sites (genomic loci) with the highest number of peaks were used for downstream analysis with Griffin. After crossing these 338 with the 404 differentially expressed TFs identified by RNA-Seq, 107 remained for which unsupervised hierarchical clustering of centered window means was performed (see Griffin analysis). Hierarchical clustering was performed using the Ward. D2 method with Euclidean distance and complete linkage settings; groupings were determined in the dendrogram using cutree_cols=2 for columns (LuCaP CRPC phenotype) and cutree_rows=13 for rows (TFs).

Selection of transcribed and promoter regions of genes For analysis of the transcribed and promoter regions of individual genes, Ensembl BioMart v104 (hg38) (Howe et al. (2021). Ensembl 2021. Nucleic Acids Research 49, D884-D891) was used to directly search for protein-coding transcription start (TSS) and transcription end (TES) coordinates. For analysis of promoter regions, a window of ±1000 bp relative to the TSS was considered. For analysis of the transcribed regions of genes, the region between the TSS and TES was considered. In the case of genes with multiple transcripts, the analysis was limited to the longest transcript, resulting in 19,336 regions. In the downstream analysis of LuCaP PDX cfDNA, if any lineage did not meet certain criteria for a region (containing a differential open histone modification region), that feature/region combination was excluded from the analysis, thereby allowing a variable number of regions to be considered based on that feature. These criteria included a total number of fragments in the region of at least 10 for all fragment size analyses (see below), and a non-zero number of "short" and "long" fragments for short-to-long ratios; short-to-long ratios less than 0.01 or greater than 10.0 were also excluded as outliers. Any regions without coverage within the lineage were excluded from all analyses. This resulted in a list of genes with different numbers between genome contexts and feature types.

も算出した。本試験において使用した断片サイズ解析コードおよび実装はgithub.com/GavinHaLab/CRPCSubtypingPaper/tree/main/FragmentAnalysisにおいて入手することができる。 Fragment size analysis Fragments were first filtered to remove duplicates and low quality alignments (mapping quality < 20) and by fragment length (15-500 bp). For each genomic locus/window, the fragment short-length ratio (FSLR) was computed as the ratio of short fragments (15-120 bp) to long fragments (140-250 bp). For each selected window, the mean, median absolute deviation (MAD: median (|X _i - median (X)|)) of the fragment length distribution, and the mean of the fragment length distribution were calculated.

The fragment size analysis code and implementation used in this study is available at github.com/GavinHaLab/CRPCSubtypingPaper/tree/main/FragmentAnalysis.

ctDNA Tumor-Normal Mixes and Benchmarking Mixes to assess benchmarking performance were constructed using five ARPC lines (LuCaP 35, 35CR, 58, 92, 136CR) and five NEPC lines (LuCaP 49, 93, 145.2, 173.1, 208.4) mixed with a single healthy donor plasma line (NPH004, EGAD00001005343) to 1%, 5%, 10%, 20%, and 30% tumor fraction, assuming a mean coverage of approximately 25x and 100% tumor fraction in mouse-subtracted PDX sequencing data. Chromosomal DNA was extracted with SAMtools (Danecek et al. (2021). Twelve years of SAMtools and BCFtools. GigaScience 10, giab008), and BAM files were merged using SAMtools after removing duplicates with Picard (broadinstitute.github.io/picard/). The mixtures were then downsampled to the number of reads corresponding to 1× and 0.2× using SAMtools to evaluate (ultra) low-pass WGS performance. During unsupervised benchmarking of each mixture, the healthy and LuCaP strains used in the mixture were excluded from the generation of feature distributions to ensure that the model training was not done from the lineages under investigation. The mixture pipeline used in this study is available at github.com/GavinHaLab/Admixtures_snakemake.

Binary supervised classification of ARPC and NEPC Binary classification of ARPC and NEPC subtypes using individual regions and feature combinations was performed using the Python implementation of XGBoost "XGBClassifier" with default parameters. Features included histone modification regions, promoters, and within transcribed regions of genes; fragment size mean, short-to-long ratio, and coefficient of variation (see fragment size analysis) in histone modification regions and promoters; center and window coverage in promoters, composite TFBS, and composite differential open chromatin regions identified by ATAC-Seq (see Griffin analysis); and Max Wave Height in promoters (see Griffin analysis). Stratified 6-fold cross-validation was applied, with two ARPC and one NEPC sample provided for each fold. This was repeated 100 times and performance was computed using the area under the receiver operating characteristic (ROC) curve (AUC) and 95% confidence intervals for each individual feature-region combination. Code and implementation of this method can be found at github.com/GavinHaLab/CRPCSubtypingPaper/tree/main/SupervisedLearning.

ichorcna Tumor Fraction Estimation Tumor fraction from patient plasma samples was assessed using ichorcna (Adalsteinsson et al. (2017). Scalable whole-exome sequencing of cell-free DNA reveals high concordance with metastatic tumors. Nature Communications 8) with a binSize of 1,000,000 bp and the hg19 reference genome. The default tumor fraction estimates reported by ichorcna were used. See github.com/GavinHaLab/CRPCSubtypingPaper/tree/main/ichorCNA_configuration for full configuration settings.

（式中、αは各試験試料についての腫瘍割合推定値である）を取る。最終的なモデルでは、４つの特徴を使用した：Ｇｒｉｆｆｉｎを使用してＬｕＣａＰ ＰＤＸ ＡＴＡＣ－Ｓｅｑ解析によって同定された特定の表現型（ＡＲＰＣおよびＮＥＰＣ）についての複合オープンクロマチン領域（中心カバレッジおよびウインドウ平均カバレッジ）（Ｇｒｉｆｆｉｎ解析を参照されたい）。各特徴ｉについて、観察された試料が腫瘍割合補正されたガウス分布の混合物に由来する確率が見いだされ、ここで、θはＮＥＰＣ混合物の重みである：
ｐ_ｉ（ｘ｜θ）＝θｐ（ｘ｜ｋ＝ＮＥＰＣ）＋（１－θ）ｐ（ｘ｜ｋ＝ＡＲＰＣ） Phenotype Prediction Model A probabilistic model was developed to classify mCRPC phenotype (ARPC or NEPC) in individual patient plasma ctDNA samples. This is an unsupervised generative mixture model that is not trained on a patient cohort of interest. However, the model accepts as prior information the tumor fraction pre-estimated from ichorcNA for a given patient ctDNA sample, as well as pre-computed ctDNA feature values from LuCaP PDX ctDNA and healthy donor ctDNA. For each patient ctDNA sample, the heterogeneous tumor fraction is fitted to a pure PDX LuCaP model. The predicted feature values (mean μ and standard deviation σ) from each phenotype k for feature i taken from the mean of the LuCaP PDX samples (μ _i,k ) or taken from the mean of normal panel H (μ _i,H , males only, n=14; see healthy donor cohort) are shifted, assuming a Gaussian distribution, so that the shifted values μ′ _i,k , σ′ _i,k are calculated according to the formula:

where α is the tumor fraction estimate for each test sample. The final model used four features: composite open chromatin regions (center coverage and window average coverage) for specific phenotypes (ARPC and NEPC) identified by LuCaP PDX ATAC-Seq analysis using Griffin (see Griffin analysis). For each feature i, the probability that the observed sample comes from a mixture of tumor fraction corrected Gaussians is found, where θ is the weight for the NEPC mixture:
p _i (x|θ)=θp(x|k=NEPC)+(1−θ)p(x|k=ARPC)

θは、範囲［０，１］を有し、値が大きいほど、ＮＥＰＣ表現型を有する試料の比率が大きいことを示し、ＮＥＰＣ予測スコア測定基準として使用される。この方法のコードおよび実装は、github.com/GavinHaLab/CRPCSubtypingPaper/tree/main/GenerativeMixtureModelにおいて見いだすことができる。 The θ parameters are estimated by maximizing the joint log-likelihood L for a given patient sample:

θ has the range [0,1], with higher values indicating a greater proportion of samples with the NEPC phenotype, and is used as the NEPC prediction score metric. The code and implementation of this method can be found at github.com/GavinHaLab/CRPCSubtypingPaper/tree/main/GenerativeMixtureModel.

であり、特異度は、ＡＲＰＣ試料の同定についての真陰性率

である。感度および特異度を最大にする代替閾値として、０．１０７７の場合、９５％感度が実現され、特異度は９３．８％に低下し、０．３７６９では、感度が８１．０％に低下したが、特異度が９８．８％に上昇した。 Analysis and Classification of Clinical Patient Samples After establishing feature distributions using the LuCaP PDX line and normal panels as described in the Generative model, the model was applied to three clinical patient cohorts (see human subjects for cohort information). Initial scoring using the model was performed on DFCI cohort I, consisting of 101 ULP-WGS samples with paired-end reads. Tumor fraction estimates predicted by ichorcNA in the original study (Berchuck et al. (2022). Detecting Neuroendocrine Prostate Cancer Through Tissue-Informed Cell-Free DNA Methylation Analysis. Clinical Cancer Research 28, 928-938) and tumor phenotype classification were obtained from the original study. A prediction score threshold of 0.3314 for calling NEPC yielded optimal performance in terms of sensitivity (90%) and specificity (97.5%) and was therefore selected. Here, sensitivity is the true positive rate for the identification of NEPC samples.

and specificity is the true negative rate for identifying ARPC samples.

As alternative thresholds maximizing sensitivity and specificity, 0.1077 achieved 95% sensitivity with a reduced specificity of 93.8%, and 0.3769 reduced sensitivity to 81.0% but increased specificity to 98.8%.

The model was then validated on two cohorts, starting from the already published DFCI cohort II (Adalsteinsson et al. (2017). Nature Communications 8; Choudhury et al. (2018). Tumor fraction in cell-free DNA as a biomarker in prostate cancer. JCI Insight 3; Viswanathan et al. (2018). Structural Alterations Driving Castration-Resistant Prostate Cancer Revealed by Linked-Read Genome Sequencing. Cell 174, 433-447.e19). The analysis was restricted to 11 samples from 6 patients with paired-end reads with corresponding ULP-WGS and WGS data. Tumor fraction estimates from ichorcna were obtained from the original study (Adalsteinsson et al. (2017). Nature Communications 8). Based on clinical history, all samples were considered to be adenocarcinoma (ARPC) (see human subjects). A scoring threshold of 0.3314 determined from DFCI cohort I was used for phenotype classification.

For the UW cohort of 47 samples from 30 patients, the tumor fraction of the samples was estimated using ichorcNA as described above, while clinical phenotype was determined by medical history and expert chart review. Model performance was evaluated against the corresponding ULP-WGS and WGS data for unequivocal clinical phenotypes of ARPC and NEPC. Using the selected scoring threshold of 0.3314, the fraction of correctly predicted ARPC (n=26) and NEPC (n=5) was computed. Performance was not evaluated for the remaining 16 samples with mixed histology.

Quantitative and statistical analysis Quantitative and statistical methods for high-throughput sequencing data analysis are described in the methods above. When comparing nonparametric (non-normal) distributions of values of a particular parameter in a population (using box plots or in tables), a two-tailed Mann-Whitney U test (also known as the Wilcoxon rank sum test; scipy.stats.mannwhitneyu, (Virtanen et al. (2020). SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat Methods 17, 261-272) was used to test whether any two distributions being compared were significantly different, and Benjamini-Hochberg (statsmodels.stats.multitest.fdrcorrection, statsmodels.org) correction was applied for multiple testing scenarios. All box plots have medians represented by center lines, interquartile ranges (IQR) represented by boxes, and the first quartile -1.5IQR and third quartile +1.5IQR represented by whiskers. PCA was performed in Python (sklearn.decomposition.PCA; scikit-learn.org).

Example 3
In Example 1, an embodiment of the Griffin workflow was applied to enhance sequence signal to enable accurate breast cancer subtype determination from low-pass sequencing data. In Example 2, the embodiment Griffin workflow approach was applied to identify other cancers, namely prostate cancer subtypes, and successfully leveraged data from alternative sequence profiling platforms (e.g., from the cut-and-run platform for nucleosome accessibility), demonstrating the power and flexibility of the Griffin analysis workflow for different cancers and input data. This example describes the development of a targeted sequencing panel for use with the Griffin workflow to understand the transcriptional characteristics of small cell lung cancer, non-small cell lung cancer, and other cancer types from blood ctDNA.

Background This example describes an innovative analytical assay based on the analysis of cell-free DNA that demonstrates clear transformation potential for clinical lung cancer diagnosis. Cell-free DNA circulating in the blood of cancer patients has been widely used to assess gene mutations and, more recently, to infer the activation of certain transcription factors by analysis of total genomic DNA. Cancer cells generate cell-free DNA through cell death, which is overwhelmingly nucleosomal, i.e., bound to histone octamers, protecting the DNA from degradation. The positioning of histones in the genome is influenced by chromatin components, including transcription factors and RNA polymerase complexes. Nucleosomes are displaced at transcription factor binding sites (TFBSs) that bind to cognate transcription factors (TFs) or transcription start sites (TSSs) that correspond to highly expressed genes. Thus, sophisticated analysis of cell-free DNA sequencing data reveals the location of nucleosomes in the cell of origin, which in turn reveals TF occupancy and gene expression. However, analyses performed to date have mainly used deep whole-genome sequencing of ctDNA, which is not practical or cost-effective for widespread clinical testing in lung cancer patients.

Results and Discussion As an innovation, we have identified highly informative TFBSs and TSSs that can be used to distinguish NSCLC from SCLC, or subtypes of NSCLC or SCLC, and thus facilitate the generation of high-resolution maps of TFBSs of key TFs (ASCL1, NEUROD1, POU2F3, REST) in SCLC, and nucleosome occupancy in gene TSSs that are markers of key transcriptional features of lung cancer cells, using DNA sequencing based on hybridization capture of ctDNA. Alternatively, these informative sites can be investigated with low-coverage whole genome sequencing to extract similar transcriptional features. Targeted capture panels are routinely applied in clinics to call mutations from ctDNA in blood, and applying targeted sequencing to assess transcriptional activity in cancer cells is highly feasible as a clinical test. This technology is particularly relevant and feasible for SCLC, which is responsible for approximately 30,000 deaths in the United States every year. Tissue sampling for SCLC is generally performed only once during a patient's disease course, often by transbronchial fine needle aspiration, resulting in very small amounts of tissue. Surgery is very rarely performed. However, SCLC has high levels of ctDNA compared to most other cancer types, reflecting its high metastatic potential, making the assay practical and potentially particularly useful for application to SCLC patients. There is also an increasing understanding that SCLC subtypes exist based on the activation of key transcription factors and their downstream programs (e.g., ASLC1, NEUROD1, and POU2F3, etc.) rather than mutations. Despite the clinical urgency of characterizing these subtypes, there are no established techniques for determining SCLC transcription subtypes using blood samples. The disclosed targeted assay is designed to identify SCLC transcription subtypes from ctDNA, thus providing a powerful clinical application of the assay. Furthermore, the panel is designed to call exonic genetic mutations from a panel of approximately 600 genes. Thus, the assay has broad clinical utility for correlation analysis of both mutations and transcriptional activity in clinical samples.

Another important current challenge in lung cancer management is the transdifferentiation of driver mutation-positive lung cancer to SCLC, which generally occurs after long-term disease control using targeted therapy. Transdifferentiation to SCLC is treated differently than advanced disease, but without notable histological changes. However, transdifferentiation can be significantly underdiagnosed, since it can currently only be assessed by biopsy of advanced lesions, which is often not feasible or desirable. This assay can also be applied to lung adenocarcinoma patients who have developed resistance to EGFR inhibitors to determine whether the resistance is related to the activation of the SCLC transcriptional profile. Thus, the main non-invasive applications include:
- Tumor classification of SCLC by estimating the activity of transcription factors (ASCL1, NEUROD1, POU2F3, REST) - Differentiating between SCLC and NSCLC - Differentiating between the major NSCLC histological subtypes: adenocarcinoma and squamous cell carcinoma - Estimating mixed histology - Detecting potential subtype changes during treatment in "real time" (either treatment-induced transcriptional subtype changes in SCLC or NSCLC to SCLC changes occurring as a resistance mechanism to targeted therapies such as EGFR inhibition in NSCLC).

All of these applications have great potential for investigating novel therapeutics and informing clinical treatment decisions.

A schematic of the panel design showing the generation of the capture panel is presented in FIG. 16. More specifically, the approach involved rationally designing a targeted sequencing panel for integrated detection of SCLC gene mutations, transcription factor (TF) subtype identity, and expression of key gene programs. Public mutation databases and functional mutation data were interrogated for coding mutations in approximately 600 genes associated with SCLC. For TF subtype identity, TFBSs of four key SCLC-associated TFs (ASCL1, NEUROD1, POU2F3, and REST) were targeted. For expression of other key gene programs, TSSs corresponding to the majority of protein-coding genes in the genome were targeted. To select specific sites, multiple data sources were incorporated as follows: For SCLC TFBSs, ChIP-seq data was used to identify TFBSs, resulting in 4-30k sites per factor. These candidate sites were then annotated with the distance to the nearest gene TSS. Sites that were retained were those whose nearest gene TSS was a gene known to be upregulated in SCLC cells expressing the factor of interest as determined by available RNAseq data. This resulted in approximately 400-700 SCLC-focused sites per factor. In the final probe set, we targeted a 1 kb window (500 bp on each side) symmetrically encompassing these ∼2 k sites. For TSS profiling, we started with established transcript annotations and removed TSSs corresponding to non-coding transcripts, Y chromosome genes, and multi-exon genes with low annotation confidence, resulting in approximately ∼36 k theoretically targeted TSSs. In the probe set, we targeted a region 260 bp downstream and 100 bp upstream of the TSS. The use of application-specific orthogonal chromatin profiling data to select sites is a key feature of the approach. However, please note that other types of chromatin profiling data, such as ATAC-seq, cut-and-run/TAG, DNAse-seq, histone modification ChIP-seq, etc., can be readily substituted or added to obtain the same or similar results.

A data analysis pipeline was developed to quantify nucleosome-protected cfDNA fragments in both TFBS- and TSS-captured DNA. The analysis pipeline, Griffin (described in more detail above), uses fragment-length-based GC correction to remove signal-obscuring GC bias. A fragment-size-aware GC bias correction approach helps maximize signal-to-noise and optimize analysis of captured DNA.

17A and 17B illustrate detection of transcription factor (TF) expression in SCLC models using targeted sequencing of cfDNA. FIG. 17A is a schematic of the experimental workflow for proof of concept of a negative control sample ("healthy donor") and a positive control sample ("flank tumor" from a SCLC cell model). FIG. 17B graphically illustrates aggregated coverage across TFBS in targeted sequencing data for a healthy donor (top row) and a flank tumor (bottom row). It is expected that TFBS will be located at position 0 on the x-axis. Data are coded per predicted TF expression. cfDNA from a healthy donor is expected to reflect expression of REST, but not ASCL1, NEUROD1, or POU2F3. Systematic differences in coverage distribution as a function of TF expression were evident in the SCLC model.

Figures 18A-18C illustrate the inference of transcription factor activity using TFBS coverage distributions from SCLC patient samples for which corresponding tumor gene expression data were available. Figure 18A graphically illustrates aggregated coverage across TFBS in targeted sequencing data for healthy donors (top row) and patients with SCLC for which corresponding tumor tissue with gene expression data was available (bottom row). Samples are coded by predicted TF expression. Systematic differences in coverage distribution as a function of predicted TF expression are again evident. Figure 18B illustrates gene expression of significant genes in selected patient samples displayed as a heat map. Cells are coded by Z-score and inset letters are log2(TPM+1). Figure 18C illustrates peak and trough amplitudes calculated from the coverage distribution across TFBS in each patient sample displayed as a heat map. Amplitudes are indicated by color and are also indicated in inset letters. The magnitude of trough depth corresponds to the gene expression of key TFs in these bona fide SCLC patient samples.

Figure 19 is a series of graphs illustrating the quantification of transcription factor binding site peak and trough amplitudes across sample types. Distribution of peak and trough amplitudes in TFBS calculated from aggregated coverage distributions according to predicted ground truth of TF expression. Pdx samples labeled "non-SCLC" are NSCLC pdx models. Patient samples labeled "non-SCLC" are either samples from patients with NSCLC (n=11) or from patients without diagnosed malignancies (n=4). Peak and trough amplitudes of ASCL1 sites are associated with both SCLC status and ASCL1 positivity, whereas peak and trough amplitudes of NEUROD1 and POU2F3 are associated with TF positivity only.

20A and 20B graphically illustrate gene expression inference using TSS coverage distributions in flank tumor positive control samples. FIG. 20A illustrates TSS coverage distributions from targeted sequencing of cfDNA grouped by quintiles of gene expression in SCLC flank tumor models (quintiles 1-5) and blood ("B", dark blue). Shown are 1,912 TSSs corresponding to 1,213 genes selected based on low expression in whole blood and correlation of TSS coverage distribution with gene expression. TSS coverage distributions vary systematically depending on the expression of the corresponding genes. FIG. 20B illustrates receiver operating characteristic curves for prediction of gene expression below threshold (shown for thresholds 0.1, 0.5, 1.0, and 2.0) as inferred from the coverage distributions of the above or corresponding TSSs. Estimates of gene expression were calculated from the TSS coverage profile as the magnitude of the difference between the average coverage depth at positions +130 and +145 relative to the TSS minus the average depth at positions -45, -30, and -15 (shown as dotted lines in 20A). The AUC of the ROC curves is shown in parentheses for each gene expression cutoff. In this preliminary analysis, which is restricted to particularly variable and therefore challenging genes, the TSS coverage distribution can be used to predict, with good test properties, whether a gene is expressed above or below a certain value.

Figure 21 is a series of graphs illustrating the use of aggregated coverage profiles across a large rationally selected subset of the TSS panel to predict SCLC or NSCLC status in lung cancer Pdx models and patient samples. The graphs present examples of aggregated TSS coverage distribution across selected gene TSSs for upregulation in NSCLC (n=396) and SCLC (n=1045) for three different samples: one healthy donor, one NSCLC Pdx model, and one SCLC Pdx model. As shown overlaid on the NSCLC PDX model, amplitude features were calculated from each coverage distribution curve as the difference between coverage at position -45 and coverage at position +120 for the TSS, thereby facilitating comparison within and between samples.

Figure 22 is a series of graphs illustrating the use of aggregated coverage profiles across a large rationally selected subset of the TSS panel to predict SCLC versus NSCLC status in lung cancer PDX models and patient samples. Aggregated coverage of SCLC-specific gene TSS (y-axis, n=1045) versus NSCLC-specific gene TSS (x-axis, n=396) in plasma samples from lung cancer PDX samples (also shown as "benign" for reference for non-cancer control patients) or plasma samples from lung cancer patients. SCLC PDXs that have transdifferentiated from adenocarcinoma are identified by the thick red line.

Thus, these data demonstrate that rationally designed capture panels enable the recovery of sequence data using the Griffin workflow that, when optimized with appropriate classifiers, can accurately discriminate between SCLC and NSCLC cells from both PDF models and patient samples. There are clear technical advantages throughout the disclosed workflow.

Griffin uses special, unique cfDNA sequence data normalization for nucleosome profiling and chromatin accessibility analysis. This includes GC bias correction, repetitive sequence filtering, and local coverage normalization. All of these normalization techniques are not available in existing proof-of-concept methods such as those in Ulz P, et al. Inference of transcription factor binding from cell-free DNA enables tumor subtype prediction and early detection. Nat Commun. 2019; 10 (1): 4666. Additionally, multi-omics feature extraction from Griffin for use in building machine learning classifiers to predict cancer subtypes is unique to this approach. It is expected that the use of targeted sequencing panels will provide higher resolution, especially with respect to resolution for more similar cell types, while retaining practical costs and being more easily integrable with resequencing of regions of interest to detect genetic mutations (i.e., cancer gene panel sequencing). From the Griffin output, many features can be extracted from each binding site of interest, and a machine learning classifier can be used to predict lung cancer histological subtype subtypes from cfDNA Griffin-optimized data.

While exemplary embodiments have been illustrated and described, it will be appreciated that various changes can be made without departing from the spirit and scope of the invention.
The embodiments of the invention in which an exclusive property or privilege is claimed are defined as follows:

Claims (111)

1. A computer-implemented method for augmenting sequence read data from a cell-free DNA sample for predicting cell type, comprising:
receiving, by a computing system, sequence read data, the sequence read data comprising a plurality of fragment reads, each fragment read having a fragment length and a GC content indicating a percentage of bases in the fragment read that are G or C;
determining, by the computing system, a GC bias value for each fragment read based on the fragment length and the GC content of the fragment read;
generating, by the computing system, a genomic coverage distribution adjusted for GC bias using the sequence read data and the GC bias value; and predicting, by the computing system, the cell type based on the genomic coverage distribution. The computer-implemented method of claim 1, wherein predicting the cell type based on the genome coverage distribution includes predicting a cell phenotype. The computer-implemented method of claim 2, wherein predicting the cell phenotype includes predicting a tissue type, a cancer type, or a cancer subtype. The computer-implemented method of claim 2, wherein predicting the cell phenotype includes predicting the expression of one or more genes of interest. determining the GC bias value based on the fragment length and the GC content of the fragment reads,
counting the number of observed reads for each combination of fragment length and GC content to determine a GC count for the sequence read data;
dividing the GC counts by the corresponding GC frequencies in a GC frequency matrix to determine a GC bias for each fragment length;
normalizing the average GC bias for each fragment length to determine an approximate GC bias value; and
The computer-implemented method of claim 1 , comprising smoothing the approximate GC bias value to determine the GC bias value. The computer-implemented method of claim 5, wherein the GC occurrence frequency matrix stores the occurrence frequency of each GC content for each of a plurality of fragment lengths within a mappable region of a reference genome. The computer-implemented method of claim 6, wherein the plurality of fragment lengths includes fragment lengths ranging from a short length threshold to a long length threshold. The computer-implemented method of claim 7, wherein the short length threshold is in the range of 10 to 20 base pairs and the long length threshold is in the range of 450 to 550 base pairs. The computer-implemented method of claim 8, wherein the short length threshold is 15 base pairs and the long length threshold is 500 base pairs. 2. The computer-implemented method of claim 1, further comprising: determining genomic regions of interest with respect to a cell type; and filtering the genomic regions of interest to identify sites that provide information about the cell type. determining the genomic region of interest,
11. The computer-implemented method of claim 10, comprising: determining an average mappability within a fixed-size window around each genomic region of interest; and rejecting genomic regions of interest whose average mappability is below a predetermined threshold. The computer-implemented method of claim 10, wherein filtering the genomic region of interest to identify sites that provide information about cell type includes determining sites that have differential signals between a first cell type and a second cell type. The step of generating a genome coverage distribution comprises:
determining fragment midpoints within a window around each of the sites that provides information about cell type;
assigning a weight to each fragment read based on the inverse of said GC bias value for each fragment read;
determining a GC-corrected midpoint coverage profile using the weighted fragment reads;
excluding locations that overlap the exclusion region;
determining a mean profile based on determining an average of the GC corrected midpoint coverage profiles for all sites;
11. The computer-implemented method of claim 10, comprising: smoothing the average profile to generate a smoothed average profile; and normalizing the smoothed average profile by dividing by an average of surrounding coverage to determine a normalized average profile. The computer-implemented method of claim 13, wherein the excluded regions include one or more regions within the Encoding Integrated GRCh38 Exclusion List, a centromere, a gap in the human genome assembly, a correction patch, an alternative haplotype, a region of zero mappability, or having coverage at least 10 standard deviations above the mean. predicting the cell type based on the genome coverage distribution,
generating one or more features based on the genome coverage distribution;
11. The computer-implemented method of claim 10, comprising providing the one or more features as input to a classifier model; and determining the cell type based on an output of the classifier model. The computer-implemented method of claim 15, wherein the one or more features include an average of coverage within a first predefined window around each site that provides information about a cell type, an average of coverage within a second predefined window of a different size than the first predefined window around each site that provides information about a cell type, and an amplitude of the genome coverage distribution around each site that provides information about a cell type. The computer-implemented method of claim 16, wherein the first predetermined window is larger than the second predetermined window. The computer-implemented method of claim 17, wherein the width of the first predetermined window is within a range of 1800 to 2200 base pairs, and the width of the second predetermined window is within a range of 40 to 80 base pairs. The computer-implemented method of claim 18, wherein the first predetermined window has a width of 2000 base pairs and the second predetermined window has a width of 60 base pairs. The amplitude of the genome coverage distribution around each site that provides information about cell type is
trimming the genome coverage distribution to a window containing 10 peaks;
17. The computer-implemented method of claim 16, wherein the magnitude of a 10th frequency is determined by performing a Fast Fourier Transform on the window of the genome coverage distribution; and determining a magnitude of a 10th frequency. The computer-implemented method of claim 15, wherein the classifier model comprises a logistic regression model, an artificial neural network, a decision tree, a support vector machine, or a Bayesian network. 1. A method for determining a chromatin accessibility profile for a cell of interest from a sample comprising cell-free DNA derived from said cell of interest, comprising:
obtaining sequence read data from the cell-free DNA;
receiving, by a computing system, sequence read data, the sequence read data comprising a plurality of fragment reads, each fragment read having a fragment length and a GC content indicating a percentage of bases in the fragment read that are G or C;
determining, by the computing system, a GC bias value for each fragment read based on the fragment length and the GC content of the fragment read;
generating, by the computing system, a GC bias adjusted genome coverage distribution using the sequence read data and the GC bias value; and determining a chromatin accessibility profile from the genome coverage distribution. 23. The method of claim 22, further comprising determining the cellular phenotype of interest based on a chromatin occupancy profile. 24. The method of claim 23, wherein the step of determining the cell phenotype comprises determining a tissue type, a cancer type, a cancer subtype, an aggressive phenotype of a malignant tumor, and/or a drug responsive phenotype. 23. The method of claim 22, further comprising performing one or more steps as recited in one or more of claims 5 to 21. 1. A method for determining the cell type of a cell of interest from a sample comprising cell-free DNA derived from said cell of interest, comprising:
obtaining sequence read data generated from said sample comprising cell-free DNA;
22. A method comprising: performing a computer-implemented method according to any one of claims 5 to 21; and determining the cell type of the cell of interest based on a prediction provided by the computing system. 27. The method of claim 26, wherein the step of determining the cell type comprises determining a cell phenotype. 28. The method of claim 27, wherein determining the cell phenotype comprises determining a tissue type, a cancer type, a cancer subtype, an aggressive phenotype of a malignant tumor, and/or a drug responsive phenotype. 28. The method of claim 27, wherein determining the cell phenotype comprises determining the expression of one or more genes of interest. 1. A method for detecting the presence of cancer cells in a subject, comprising:
obtaining sequence read data generated from a sample comprising cell-free DNA obtained from the subject;
22. A method comprising: performing a computer-implemented method according to any one of claims 5 to 21; and determining the presence of cancer cells in the subject based on a prediction provided by the computing system. 31. The method of claim 30, which is performed multiple times over time, and at each time further characterizing the cancer cell(s) detected in the subject to determine a cancer subtype or phenotype of the detected cancer cell(s) based on the predictions provided by the computing system. 32. The method of claim 31, further comprising the step of detecting a change in the phenotype of the detected cancer cell(s) over time, the step being performed multiple times over time. 33. The method of claim 31 or 32, wherein the subject undergoes cancer treatment during administration of the method, and the method further comprises determining the responsiveness of the cancer cell(s) to the treatment. 1. A method for determining a cancer subtype of a target cancer cell from a sample containing cell-free DNA derived from the target cancer cell, comprising:
obtaining sequence read data generated from said sample comprising cell-free DNA;
22. A method comprising: performing a computer-implemented method according to any one of claims 5 to 21; and determining the cell type of a cell of origin based on a predicted cancer subtype provided by the computing system. The method of claim 34, wherein the sample is obtained from a subject having cancer. The method of any one of claims 30 to 35, wherein the cancer is characterized as metastatic breast cancer. The method of claim 36, wherein determining the cancer subtype comprises determining whether the cancer is ER+ or ER-. The method of claim 36, wherein determining the cancer subtype comprises determining whether the cancer is PR+ or PR-. The method of claim 36, wherein determining the cancer subtype comprises determining whether the cancer is HER2+ or HER2-. determining the cancer subtype,
whether the cancer is ER+ or ER-;
37. The method of claim 36, comprising determining two or all of: whether the cancer is PR+ or PR-; and whether the cancer is HER2+ or HER2-. The method of any one of claims 30 to 35, wherein the cancer is characterized as metastatic prostate cancer. The method of claim 41, wherein determining the cancer subtype comprises determining whether the cancer is AR+ (ARPC) or AR-. The method of claim 41, wherein determining the cancer subtype includes determining whether the cancer is ARPC or AR-low. 42. The method of claim 41, wherein determining the cancer subtype comprises determining whether the cancer has a neuroendocrine prostate cancer (NEPC) phenotype signature. 42. The method of claim 41, wherein determining the cancer subtype comprises determining whether the cancer is bisecretory. determining the cancer subtype,
whether the cancer is AR+ (ARPC) or AR-;
Whether the cancer is AR-low or ARPC;
whether the cancer has a neuroendocrine prostate cancer (NEPC) phenotypic signature;
Whether the cancer is AR-low or NEPC;
42. The method of claim 41, comprising determining two or all of whether the cancer is bisecretory, ARPC, or NEPC. The method of any one of claims 30 to 35, wherein the cancer is characterized as lung cancer. 48. The method of claim 47, wherein determining the cancer subtype comprises determining whether the cancer is small cell lung cancer (SCLC) or non-small cell lung cancer (NSCLC). 49. The method of claim 48, further comprising determining whether the NSCLC is an adenocarcinoma or a squamous cell carcinoma. The method of any one of claims 47 to 49, wherein the sequence read data is generated from a panel of genomic targets. 51. The method of claim 50, wherein the panel of genomic targets comprises transcription factor binding sites (TFBS) of one or more transcription factors associated with SCLC. 52. The method of claim 51, wherein the one or more transcription factors associated with SCLC include one or more of ASLC, NEUROD1, POU2F3, REST, etc., and the method includes determining the nucleosome occupancy of the TFBS. The method of claim 52, wherein the TFBS is retained in the panel if it is proximal to the transcription start site of a gene associated with lung cancer, as identified, for example, by ChIP-seq data. 51. The method of claim 50, wherein the panel of genomic targets includes transcription start sites (TSSs) of one or more markers associated with lung cancer, and the method includes determining nucleosome occupancy of the TSSs. The method of any one of claims 30 to 54, wherein the sample is obtained from a subject, and the method further comprises administering an effective treatment to the subject based on the determined cancer subtype. 56. The method of any one of claims 35 to 55, further comprising performing the method on a plurality of samples obtained from the subject at a plurality of distinct time points after an initial diagnosis of cancer. The method of any one of claims 22 to 56, wherein the sequence read data is generated by ultra-low pass whole genome sequencing. The method of any one of claims 22 to 56, wherein the sequence read data is generated by a chromatin accessibility assay. The method according to any one of claims 22 to 56, wherein the sequence read data is generated by the ATAC-seq method. The method according to any one of claims 22 to 56, wherein the sequence read data is generated by the ChIP-seq method. The method of any one of claims 22 to 56, wherein the sequence read data is generated by a DNase sensitivity assay. The method of any one of claims 22 to 56, wherein the sequence read data is generated by a cut-and-run assay. 63. The method of claim 62, wherein the cut-and-run assay incorporates an affinity reagent that targets a post-translational modification to one or more of H3K27ac, H3K4me1, and H3K27ac. The method of any one of claims 1 to 63, further comprising the step of generating the sequence read data. The method of any one of claims 1 to 64, wherein the sequence read data comprises sequence read data generated from a panel of genomic targets. 66. The method of claim 65, wherein the panel of genomic targets comprises transcription factor binding sites (TFBS) for one or more transcription factors associated with a cancer type of interest. The method of claim 66, comprising determining the nucleosome occupancy of the TFBS. The method of claim 66, wherein the TFBS is retained in the panel if it is proximal to the transcription start site of a gene associated with the cancer type of interest, as identified, for example, by ChIP-seq data. 67. The method of claim 66, wherein the panel of genomic targets comprises transcription start sites (TSSs) of one or more markers associated with the cancer type of interest, and the method comprises determining nucleosome occupancy of the TSSs. The method of any one of claims 22 to 64, wherein the sample is blood, plasma, or serum. 1. A computer-implemented method for augmenting sequence read data from a cell-free DNA sample for predicting cell type, comprising:
receiving, by a computing system, sequence read data, the sequence read data including a plurality of fragment reads, each fragment read having a fragment length;
determining, by the computing system, a variability in fragment sizes for at least one gene associated with a cell type; and predicting, by the computing system, the cell type based on the variability in fragment sizes for the at least one gene. 72. The computer-implemented method of claim 71, wherein the step of determining the variability of fragment sizes includes determining a coefficient of variation of fragment sizes. The computer-implemented method of claim 71, wherein predicting the cell type based on the genome coverage distribution includes predicting a cell phenotype. The computer-implemented method of claim 73, wherein predicting the cell phenotype includes predicting a cancer subtype. The computer-implemented method of claim 74, wherein predicting the cell phenotype includes predicting a cancer subtype of prostate cancer. The computer-implemented method of claim 75, wherein predicting the cancer subtype includes distinguishing between ARPC and NEPC. predicting the cell type based on the variability of the fragment sizes,
generating one or more features based on the variability in fragment sizes;
72. The computer-implemented method of claim 71, comprising providing the one or more features as input to a classifier model; and determining the cell type based on an output of the classifier model. 78. The computer-implemented method of Claim 77, wherein generating the one or more features based on the variability in fragment size comprises generating a log ₂ fold change value of a coefficient of variation in fragment size in a first cell type relative to a second cell type. 80. The computer-implemented method of claim 78, wherein the log ₂ fold change value predicts at least one of gene expression and gene transcription activity between the first cell type and the second cell type. The computer-implemented method of claim 78, wherein the first cell type is an ARPC cell and the second cell type is a NEPC cell. The computer-implemented method of claim 77, wherein the classifier model comprises a logistic regression model, an artificial neural network, a decision tree, a support vector machine, or a Bayesian network. 1. A method for determining the cell type of a cell of interest from a sample comprising cell-free DNA derived from said cell of interest, comprising:
obtaining sequence read data generated from said sample comprising cell-free DNA;
82. A method comprising: performing a computer-implemented method according to any one of claims 71 to 81; and determining the cell type of the cell of interest based on a prediction provided by the computing system. 83. The method of claim 82, wherein the step of determining the cell type comprises determining a cell phenotype. The method of claim 83, wherein determining the cell phenotype comprises determining a cancer subtype. 85. The method of claim 84, wherein determining the cancer subtype comprises distinguishing between ARPC and NEPC. 1. A method for detecting the presence of cancer cells in a subject, comprising:
obtaining sequence read data generated from a sample comprising cell-free DNA obtained from the subject;
82. A method comprising: performing a computer-implemented method according to any one of claims 71 to 81; and determining the presence of cancer cells in the subject based on a prediction provided by the computing system. 87. The method of claim 86, wherein the method is performed multiple times over time, and during each of the multiple times, the cancer cell(s) detected in the subject are further characterized to determine a cancer subtype or phenotype of the detected cancer cell(s) based on the predictions provided by the computing system. The method of claim 87, further comprising a step of detecting changes in the detected cancer cell phenotype over time, the step being performed multiple times over time. The method of claim 87 or 88, wherein the subject receives a cancer treatment during the administration of the method, and the method further comprises a step of determining the responsiveness of the cancer cells to the treatment. 1. A method for determining a cancer subtype of a target cancer cell from a sample containing cell-free DNA derived from the target cancer cell, comprising:
obtaining sequence read data generated from said sample comprising cell-free DNA;
82. A method comprising: performing a computer-implemented method according to any one of claims 71 to 81; and determining the cell type of a cell of origin based on a predicted cancer subtype provided by the computing system. The method of claim 90, wherein the sample is obtained from a subject having cancer. The method of any one of claims 86 to 91, wherein the cancer is characterized as metastatic prostate cancer. The method of claim 92, wherein determining the cancer subtype comprises determining whether the cancer is AR+ (ARPC) or AR-. The method of claim 92, wherein determining the cancer subtype includes determining whether the cancer is ARPC or AR-low prostate cancer (ARLPC). 93. The method of claim 92, wherein determining the cancer subtype comprises determining whether the cancer has a neuroendocrine prostate cancer (NEPC) phenotype signature. The method of any one of claims 86 to 95, wherein the sample is obtained from a subject, and the method further comprises administering an effective treatment to the subject based on the determined cancer subtype. The method of any one of claims 86 to 96, further comprising performing the method on a plurality of samples obtained from the subject at a plurality of distinct time points after an initial diagnosis of cancer. The method of any one of claims 82 to 97, wherein the sequence read data is generated by ultra-low pass whole genome sequencing. The method of any one of claims 82 to 97, wherein the sequence read data is generated by a chromatin accessibility assay. The method according to any one of claims 82 to 97, wherein the sequence read data is generated by the ATAC-seq method. The method according to any one of claims 82 to 97, wherein the sequence read data is generated by the ChIP-seq method. The method of any one of claims 82 to 97, wherein the sequence read data is generated by a DNase sensitivity assay. The method of any one of claims 82 to 97, wherein the sequence read data is generated by a cut-and-run assay. The method of claim 103, wherein the cut-and-run assay incorporates an affinity reagent that targets a post-translational modification to one or more of H3K27ac, H3K4me1, and H3K27ac. The method of any one of claims 71 to 104, further comprising generating the sequence read data. The method of any one of claims 71 to 105, wherein the sequence read data is generated from a panel of genomic targets. 107. The method of claim 106, wherein the panel of genomic targets comprises transcription factor binding sites (TFBS) for one or more transcription factors associated with a cancer type of interest. The method of claim 107, comprising determining the nucleosome occupancy of the TFBS. The method of claim 107, wherein the TFBS is retained in the panel if it is proximal to the transcription start site of a gene associated with the cancer type of interest, as identified, for example, by ChIP-seq data. 107. The method of claim 107, wherein the panel of genomic targets comprises transcription start sites (TSSs) of one or more markers associated with the cancer type of interest, and the method comprises determining nucleosome occupancy of the TSSs. The method of any one of claims 82 to 110, wherein the sample is blood, plasma, or serum.