Non-invasive multi-cancer detection using DNA hypomethylation of LINE-1 retrotransposons

Targeting primate-specific LINE-1 elements reveals plasma DNA-methylation patterns genome-wide

We developed a PCR-based targeted bisulfite method coupled to deep sequencing to detect methylation patterns of L1PA elements. We used sodium bisulfite chemical conversion to achieve base-pair resolution analysis and designed a multiplexed PCR based on 8 amplicons covering L1PAs (Fig. 1A, Table S1, Fig. S1A). We detected thousands of L1PA elements scattered throughout the genome as shown by the genomic hits obtained from a healthy plasma, an ovarian tumor, and a uveal melanoma tumor sequenced at high depth (Fig. 1B, Table S2). We observed similar profiles for the three samples, as well as for healthy and cancer plasmas with standard coverage (Fig. S1B-E). This demonstrated the robustness of the approach. Overall, the estimated number of L1PA targets is about 30-40,000 elements per genome including half of the human specific copies (L1HS) and many copies of the other L1PA subfamilies (Fig. 1C, Table S2). This represents an estimate of 87-120,000 CpG sites.

• Download figure
• Open in new tab

Fig 1. Targeting primate-specific LINE-1 elements reveals plasma DNA-methylation patterns genome-wide

A. CpG density along the structure of a human specific LINE-1 (L1HS) element, which contains 95 CpG. The DIAMOND assay targets 30 CpG. Each target amplicon is highlighted by a black bar below the structure. The number of CpG sites detected per amplicon is displayed in blue. B. L1PA copy number hit by uniquely and/or randomly mapped reads, obtained from a healthy plasma versus ovarian (OVC, top track) or uveal melanoma (UVM, middle track) tumor tissue samples ‘deep sequenced’ (54M, 44M or 46M reads respectively) over the distribution of L1PA elements annotated in the genome (RepeatMasker on hg38, grey bottom track). C. Histogram summarizing the most represented sub-families of L1 targeted by the DIAMOND assay in the 3 ‘deep sequenced’ samples, in descending order (sum of copies across the 3 samples). The colors highlight the relative contribution of L1PA copies hit by reads uniquely mapped, randomly mapped or both. D. Methylation pattern observed across the 8 regions targeted along the L1 element in the healthy plasma sample ‘deep-sequenced’. Metaplot showing the average methylation levels at each CpG position. Amplicon limits are delineated with grey dotted lines. The dark line marks the end of the 5’UTR. Average levels per amplicon are indicated.

Following deep sequencing, reads are traditionally mapped back to the genome. However, the majority of sequencing reads from repetitive sequences are assigned randomly during mapping steps and are subsequently lost for classical differentially methylated region (DMR) calling ⁴⁹. We, thus, developed a new computational pipeline to accurately align repetitive sequencing data without using a reference genome (Fig. S1F). To perform this, we clustered all good quality reads based on their similarity, extracted representative sequences from the largest clusters and used them for multiple sequence alignment. We then aligned all the reads back onto this custom database. Using such reference-free method, we preserved the majority of our data and could extract the informative CpG sites agnostically. We selected sites with a CG/TG content ≥ 20% including at least 5% of CG to ensure that the position of interest carries some DNA methylation marks. This selection was done on healthy samples to avoid biases related to cancer hypomethylation. We retrieved 35 CpG positions covered by our panel including two additional CpGs with respect to the L1HS consensus annotations, located within amplicon 2 (Fig. S2A-B).

As expected, the 5’ end of the L1 copies targeted is heavily methylated ^42,50, particularly within the 2^nd amplicon. We also observed quite high levels in both the 5^th amplicon (69% in average, Fig. 1D), which covers part of the ORFI, and the last two CpGs of amplicon 8, which is located immediately upstream of the 3’UTR. Amplicon 3, which has the lowest methylation levels within the 5’ end, displayed sequencing data with atypical distributions and showed less robust performances (not shown). Hence, we further eliminated it from the rest of the study, resulting in a total of 30 CpG positions analyzed. Overall, this reference-free method retrieved methylated sites contained by the youngest LINE-1 elements present in the human genome allowing us to study their DNA-methylation levels and motifs from minute amount of DNA such as plasma cfDNA.

L1PA hypomethylation is detectable from plasma DNA in multiple forms of cancer

We first tested the DIAMOND approach on methylation controls, cancer cell lines and tissue samples. The overall methylation levels demonstrated an extensive L1PA hypomethylation specifically in cancer samples, including colorectal (CRC), ovarian (OVC), breast (BRC) and uveal melanoma (UVM) cancer cell lines as well as OVC, BRC and UVM tumors compared to healthy white blood cells and healthy tissues collected adjacent to ovarian tumors (Fig. 2A, Table S3). Next, we tested a cohort of 473 plasma samples including 123 healthy plasma controls and plasma samples from patients with 6 different types of cancer, covering metastatic (M+) and localized (M0) stages (Table S4. This includes colorectal and ovarian cancers in which a substantial rate of L1 hypomethylation has previously been reported ^51,52. We detected a statistically significant L1PA hypomethylation in cfDNA of metastatic colorectal cancer (CRC M+), breast cancer (BRC M+) and uveal melanoma (UVM M+) samples as well as in locally advanced ovarian cancers (OVC M0, stages III) and localized gastric cancers (GAC M0) (Fig. 2B, Table S3). The global methylation was not significantly different in metastatic non-small cell lung cancers (LC M+) nor in localized stages of breast cancer (BRC M0). However, focusing strictly on global methylation levels provides only part of the information.

• Download figure
• Open in new tab

Fig 2. L1PA hypomethylation is detectable from plasma DNA in multiple forms of cancer

A. Global DNA methylation of fully methylated (healthy white blood cells (WBC) DNA treated with SssI) and unmethylated (whole-genome amplified healthy WBC DNA) controls, cancer cell lines or tissues. Ovarian healthy tissues were collected next to ovarian tumors. The global methylation levels for each sample correspond to the percentage of CG dinucleotides at each CpG site averaged by the number of CpG sites. Statistical differences between controls, cell lines or tissues and healthy white blood cells (WBC) were computed using Mann–Whitney U test (p_{Fully_meth.} = 9.97e-07, p_Unmeth. = 1.86e-06, p_{CRC_Cells} = 0.266, p_{OVC_Cells} = 1.20e-02, p_{BRC_Cells} = 2.47e-06, p_{UVM_Cells} = 6.77e-05, p_{Healthy_OVC_Tissues} = 0.063, p_{OVC_Tissues} = 8e-03, p_{BRC_Tissues} = 4.10e-04, p_{UVM_Tissues} = 0.88, Table S3). B. Global DNA methylation in cancer plasma including metastatic stages (M+) and non-metastatic stages (M0) as well as healthy donor plasmas. Statistical differences between each cancer subgroup and healthy samples were computed using Mann–Whitney U test (p_{CRC_M+} = 1.27e-29, p_{BRC_M+} = 3.79e-19, p_{UVM_M+} = 8.29e-06, p_{LC_M+} = 0.655, p_{OVC_M0} = 1.94e-05, p_{GAC_M0} = 4.28e-08, p_{BRC_M0} = 9.10e-01, Table S3). Black dotted lines represent the median. C. Methylation level at each targeted CpG sites (x-axis), for each healthy sample (y-axis) depicted as a heatmap. CpG numbers are indicated. The metaplot represents the average methylation levels of the population. Amplicon numbers are indicated. D. Differential methylation levels between healthy samples and patients for each type of cancer represented as metaplots. E. Proportion of methylation motifs, called haplotypes, for each amplicon (mean centered per amplicon). Only the most important features are represented (see Fig. 3F and methods). Blue arrows highlight the most abundant haplotype in each amplicon. F. Mean centered abundance of the most important haplotypes with the highest co-methylation patterns (mostly fully methylated or fully unmethylated molecules) in cancer subgroups compared to healthy donors. Statistical significances were computed using Mann–Whitney U test on raw haplotype proportions (Table S5).

We further computed the levels of methylation at each CpG target (n=30) for these plasma samples and observed specific patterns of methylation along the L1 structure, which are robustly conserved among the 123 healthy donors (Fig. 2C). When considering all cancer samples together, we observed a clear difference with the methylation of healthy samples. We detected a steady hypomethylation through all CpG targets except for the two sites within amplicon 8 (Fig. 2D). This is also true for metastatic colorectal cancers (CRC M+), breast cancers (BRC M+) and uveal melanoma (UVM M+). Clear hypomethylation is also observable for localized gastric (GAC M0) and ovarian (OVC M0) cancers, in particular at amplicon #1, #4 and #6, while the differences are less striking for localized breast cancers (BRC M0) and metastatic non-small cell lung cancers (LC M+). The distinction between most cancers and healthy samples were dependent on multiple CpG positions belonging to different amplicons along L1s, as shown by PCA analysis (Fig. S2C). The least discriminating positions were located within amplicon 8, which is consistent with the metaplots shown in Fig. 2D.

Next, we analyzed the motifs of methylation at the molecule level, which provide a more detailed signal. These haplotypes correspond to true patterns of methylation of adjacent CpGs, detected for each amplified DNA molecule. This was achieved by the incorporation of unique molecular identifiers (UMIs) into the library (Fig. S1A). Based on the combination of the 30 CpG targets divided into their 7 amplicons, we extracted a total of 372 unique features (Fig. S2D). We observed highly robust representation profiles of haplotypes among the 123 healthy samples (Fig. 2E). For most amplicons, the fully methylated molecules were the most represented, as expected for healthy controls. However, we observed a high proportion of totally unmethylated haplotypes in amplicon #6 and #7. This can be explained by the fact that older L1 copies are often truncated in 5’ and less regulated by DNA methylation, leading to the capture of molecules with lower DNA methylation in 3’. Nevertheless, several intermediate patterns were also among the most important features and were found to be differentially represented in healthy and cancer samples (Table S5). Fully methylated haplotypes were significantly under-represented in most cancer subgroups and in most amplicons (Fig. 2F). On the contrary, fully unmethylated haplotypes were over-represented in most cancer subgroups and in most amplicons. This is also well illustrated by the PCA analysis shown in Fig. S2E, underlining the contribution of highly methylated haplotypes towards the healthy group versus the lowly methylated haplotypes separating cancer samples (middle panel). This separation involves haplotypes from all amplicons (right panel).

Next, we compared the methylation profiles of tumor and plasma paired samples (OVC = 10, Fig. S3A-F; BRC = 16, Fig. S3G-K) by calculating the correlation between their methylation differences relative to the mean methylation of healthy donor plasmas. We observed a better correlation with methylation haplotype portions than with single CG methylation features (Fig. S3B-C and Fig. S3H-I, Table S6). These results demonstrate that L1 hypomethylation can robustly be observed from cancer plasma DNA at the level of single CpG sites and more importantly at the level of methylation haplotypes.

L1PA hypomethylation-based classifiers recognize samples from multiple forms of cancer

We then trained classification models using random forests, with the 30 features corresponding to the levels of methylation at each CpG target or the 372 features corresponding to the proportions of haplotypes, or both, and assessed their performances to automatically separate healthy from tumor plasmas. By testing all cancer samples without cancer-type specification, the methylation of L1PA elements showed an extremely good ability to discriminate between healthy and tumor plasmas, with an overall area under the curve (AUC) of 94-95% for the 3 types of features (Fig. 3A, 3C, S4A-B). Next, we trained distinct models to estimate the performances for each cancer type and/or dissemination stage (M0 vs M+). These models were extremely performant in metastatic colorectal and breast cancers but also stage III ovarian cancers and localized gastric cancers, with nearly perfect classifications and AUCs between 98-100% (Fig. 3B-C, S4A-B). Additionally, we observed excellent performances for metastatic lung cancers and uveal melanoma and more importantly for localized stages of breast cancer (AUC_{BRC_M0} = 92-95%). These models provide very good sensitivities at 99% specificity (Fig. 3D), in particular for CRC M+, BRC M+, OVC M0, GAC M0 and BRC M0. The latter is one of the most difficult cancer to detect non invasively, as reported in previous liquid biopsy multi-cancer tests ^11,39,40.

• Download figure
• Open in new tab

Fig 3. L1PA hypomethylation-based classifiers recognize samples from multiple forms of cancers

A-B. Receiver Operating Characteristic (ROC) curves obtained for plasma samples classification using single-CpG methylation levels (n=30) or haplotype proportions (n=372) with the ‘all cancers’ model (A) or the ‘cancer-types’ models (B). All classifications include 5000 stratified random repetitions of learning on 60% of the samples and testing on the 40% left, with undersampling for classes equilibrium (results with and without undersampling are presented in Fig. S3A-B). N_{CRC_M+} = 75, N_{BRC_M+} = 97, N_{LC_M+} = 50, N_{UVM_M+} = 55, N_{OVC_M+} = 4 (included only in ‘all cancers’ testing), N_{OVC_M0} = 18, N_{GAC_M0} = 27, N_{BRC_M0} = 23 tested versus 123 healthy donors. ROC curves shown are obtained by averaging the sensitivity and specificity of each repetition of learning. C-D. Performances for classifiers using single CpG methylation levels (grey), haplotype proportions (black) or both (white) presented as AUCs (C) or sensitivities at 99% specificity (D). Average AUCs are computed from the 5000 AUCs generated by each repetition of learning. Bars indicate 95% CI. E-F. Importance (mean decrease in impurity) of the features used by the classifiers depicted as clustered heatmaps. The features correspond to the CpG targets (E) or the haplotypes (F). Only the most important haplotypes (feature importance level >1%) are shown. G. ROC curves obtained for plasma samples classification with the 3-stage model, using haplotype features. H-I. Performances for the 3-stage classifiers using single CpG methylation levels (grey) or haplotype proportions (black) presented as AUCs (H) or sensitivities at 99% specificity (I). Early stages (I/II, N=31), locally advanced stages (III, N=30) and metastatic stages (IV, N=281) J. Cancer detection rates with the methylation-based DIAMOND assay (haplotypes and CG methylation) vs common recurrent mutations for samples assessed in previous studies (^13,14,55,56) or with NGS (Table S6).

Overall, we observed similar results using single-CpG methylation levels, haplotype proportions or both features. This can be explained by the high correlation observed between the 2 types of features (Fig. S4C). We compared these performances for each subgroup with the classification rates extracted from the model ‘all’, testing all cancer samples together. We observed similar AUCs for most cancers, but overall, they were better classified with their ‘expert’ cancer specific models (Fig. S4D-E). This shows that certain specificities of cancer type may confuse the current ‘all’ model and affect the sensitivities at 99% specificities.

Subsequently, we evaluated the importance of the features used by our classifiers (Fig. 3E-F). CpG positions displayed different patterns in the various cancer subgroups that can be informative for distinct cancer types or stages (Fig. 3E). Nonetheless, we identified features which are common to many types of cancer such as most CpGs of amplicon 1 and the first CpG of amplicon 6. Other features seemed to be characteristic of specific subgroups, such as CG7-14 which are the most important features for sorting localized stages of breast cancer (BRC M0) or CG15-18, in particular CG17, which are part of the top features for metastatic breast cancers (BRC M+). Besides their dissemination status, the different important features detected in the BRC M0 and M+ subgroups may result from the fact that they were composed of different breast cancer subtypes, with 100% of hormone dependent (HR+ HER2-) cancers in the BRC M+ subgroup and a mixture of subtypes, including also triple negative (TNBC) and HER2+ cancers, in the BRC M0 subgroup (Table S7). Indeed, we observed that slightly different hypomethylation patterns in HR+ HER2-BRC versus TNBC (Fig. S4F-H). We observed that haplotype features showed consistent patterns with important CG features (Fig. 3F). However, haplotypes provide a more detailed view of the methylation patterns with a strong importance of the most methylated or non-methylated molecules. We still observed that some methylation intermediates are important for cancer detection (ex: in amplicon #1 in CRC M+ and GAC M0, #2 in BRC M0, #4 in BRC M+ and other subgroups, #5 in LC M+ and UVM M+, #7 in OVC, #8 in BRC M0 and LC M+). BRC M+ and M0 subgroups cluster together, showing that even if they are not highly similar, they are closer to each other than to other cancer types. These breast cancer specificities are worth exploring further in the future. Overall, this suggests that L1PA methylation alterations detected from cfDNA vary in different types and stages of cancer.

To estimate the ability of DIAMOND to detect cancer at early stages of the disease, we build classifiers for 3 stage classes gathering all cancer types: early stages (I/II, N=31), locally advanced stages (III, N=30) and metastatic stages (IV, N=281). Classifications were highly performant for all 3 stage categories (AUC_Early = 95%, AUC_Adv. = 97%, AUC_Meta = 95%; Fig. 3G-H, S4I, S4D) with a mean sensitivity of 70% for early stages, (Sen_Early = 70%, Sen_Adv. = 90%, Sen_Meta = 69%; Fig. 3I, S4E).

Next, we analyzed whole genome bisulfite sequencing (WGBS) from healthy and cancer plasma from two recently published studies ^53,54 to evaluate if we can retrieve L1PA hypomethylation signal with non-targeted methods (Fig. S5A, see methods). Consistent with our findings, we observed statistically significant hypomethylation in cancer samples (Fig. S5B-H). These whole genome approaches were nevertheless not as efficient as the DIAMOND targeted approach to automatically detect cancer from L1PA hypomethylation as they cover less well the regions of interest (Fig. S5A, S5F, S5I). Indeed, they involve mapping on a reference genome leading to loss of data. Moreover, L1PA methylation measured with DIAMOND largely outperforms methods based on the detection of mutations. In comparison, the identification of the same tumor samples via the detection of frequent recurrent mutations, which is commonly used in the clinic, does not exceed 57% for ovarian cancer (Table S8), 38% for colon cancer ⁵⁵ and 52% for metastatic breast cancer ^13,14 (Fig. 3J). We particularly achieved remarkable performance on the cohort of 27 localized gastric cancers with a detection rate of 95% of true positive as compared to 12% for mutation screening ⁵⁶. This is mostly due to the fact that methylation changes occur in virtually all cancer patients, unlike recurrent mutations.

Multi-cancer classification performances are reproducible on an independent cohort

To validate the DIAMOND approach, we tested a second independent cohort consisting of 214 patients affected with the same types of cancers as in the first cohort, excluding uveal melanoma and non-metastatic gastric cancers, along with 60 healthy donors (Fig. 4A). First, we confirmed that the methylation patterns along the L1 structure were highly reproducible between healthy donors from cohorts 1 and 2, at the level of single-CpG targets (Fig. 4B) but also for haplotype proportions (Fig. 4C, Table S9). While methylation at single-CpG within cancer subgroups showed slightly more variability (Fig. S6A), global methylation levels were quite reproducible between the two cohorts, showing similar distributions and no statistical differences (Fig. 4D, Table S10), except for non-metastatic ovarian cancers. There was an important heterogeneity among the OVC M0 samples of cohort 2, which clustered into two distinct groups, while cohort 1 was more homogeneous (Fig. S6B). Notably, no correlation was found with available clinico-histopathological parameters (age, staging, CA125 level, mutational status, treatment or response to therapy). Differential haplotype proportions between healthy and cancer subgroups were also mostly conserved (Fig. S6C, Table S11). Overall, the method showed good reliability with the 7-amplicon panel used and good robustness in detecting L1 methylation levels and changes.

• Download figure
• Open in new tab

Fig 4. Multi-cancer classification performances are reproducible on an independent cohort

A. Number of patients and healthy donors (HD) in the discovery cohort (C1) and in the validation cohort (C2) for each cancer type and dissemination stage (non-metastatic: M0 vs metastatic: M+, NA: stage not available). Generated using Servier Medical Art. B. Methylation level at each targeted CpG sites (x-axis), for each healthy sample (y-axis) from C1 vs C2, depicted as a heatmap. No clustering is done on the data, which comes ordered by targeted CpG site on the x-axis (amplicon numbers are indicated). The metaplots represent the average levels for donors of C1 versus C2 at each CpG site. C. Mean centered abundance of the most important haplotypes, with the highest co-methylation patterns, in healthy donors from C1 vs C2 (Statistical differences computed using Mann–Whitney U test are available in Table S9). D. Comparison of the global levels of methylation in C1 vs C2. Methylation levels are calculated as explained previously in Fig. 2. The p-values are computed using Mann–Whitney U test (p_{CRC_M+} = 0.680, p_{OVC_M+} = 0.816, p_{BRC_M+} = 0.783, p_{LC_M+} = 0.596, p_Healthy = 0.316, p_{OVC_M0} = 4.74e-05, p_{BRC_M0} = 0.132, Table S10). Black doted lines represent the median. E. Performances for validation classifiers using haplotype features presented as ROC curves, AUCs and sensitivities at 99% specificity obtained with the ‘all’ cancers model or the ‘expert’ models for cancer subgroups. All classifications include 5000 stratified random repetitions of learning on the whole discovery cohort and testing on the whole validation cohort without undersampling. ROC curves shown are obtained by averaging the sensitivity and specificity of each repetition of learning. Average AUCs are computed from the 5000 AUCs generated by each repetition of learning. Bars indicate 95% CI. F. Performances for 3-stage ‘expert’ classifiers: early stages (I/II, N_C1=31, N_C2=38), locally advanced stages (III, N_C1=30, N_C2=54) and metastatic stages (IV, N_C1=281, N_C2=113) presented as mean ROC curves, AUCs, or sensitivities at 99% specificity. G. Performances for integrated models (‘all’ or ‘stack’) when training for the specific group tested (seen) versus when not training for this subgroup (blind). These classifications have been performed on the whole sample set including C1 and C2.

Since age-related changes in DNA methylation have been described ^57,58 and that the healthy donors included in the study are younger overall than the cancer patients (Fig. S7A), we have investigated whether there was an effect on the methylation patterns we studied. We found a significant but very small effect which appeared much smaller that the effect of disease status (Fig. S7B-C). This small effect was tending towards an increase in methylation with age (Fig. S7D), and we observed similar patterns and differences between healthy and cancer samples when adjusting for the age (Fig. S7E-G, Table S12). Furthermore, we observed similar performances in age-matched and non-age-matched cohorts extracted from C2 (Fig. S7H-J), demonstrating that age is not a confounding factor.

To validate our classifiers, we trained models on the entire first cohort and evaluated them on the second set of independent samples. The results showed excellent classification performances with an overall AUC of 88% when testing all cancers together with no annotations of their histological types and AUC between 88%-100% for the cancer subgroup ‘expert’ models with sensitivities at 99% specificity between 49%-100% (Fig. 4E) including 55% for localized breast cancers. It was, however, lower for metastatic lung cancer with a mean sensitivity of 49%. We observed that haplotype models were more robust compared to single-CpG methylation rates (Fig. S6D-G). This could be explained by the fact that haplotypes capture true methylation patterns at the molecule level, enabling to discard noise caused by experimental variability for example. Next, we applied the same validation method, training on C1 and testing on C2, for the 3-stage ‘expert’ classifiers and observed great classification performances with a mean AUC of 99% and a mean sensitivity of 79% for early stages (Fig. 4F). When comparing the performances for each subgroup with the classification rates extracted from the model ‘all’, we observed similar AUCs for most cancers, but lower rates for the BRC M0 subgroup and the early-stage group (Fig. S8A), which is mostly composed of BRC M0 (Fig. S8B).

To further investigate the generalization of our marker, we have also tested the ability of DIAMOND to detect cancers of a type or subgroup for which the model is not trained. To do so, we have combined the data from cohorts C1 and C2 and trained on the whole data set, systematically removing one cancer subgroup or one cancer type and testing specifically for this subgroup or type (see Methods, Fig. 4G, S8C). We observed performances similar to when the subgroup is included in the train set, which demonstrates the universality of L1PA hypomethylation. Classification performance is noticeably lower for the BRC M0 subgroup when the model still performs well for OVC M0 and GAC M0. This may be because the BRC M0 subgroup contains the earliest stages of both cohorts (Fig. S8B) and the model is less well trained for these stages when removing it. To improve this, we have developed another integrated model that is trained for all cancer types and stages together but also incorporates the probabilities from the ‘expert’ models, defining profiles for various stages and dissemination status. This ‘stack’ model provides a unique prediction for cancer versus healthy status for each sample at once, independently of their cancer type / dissemination status, and allows to balance for the cancer characteristics used to train the model by including them all. We have then extracted the performances per cancer subgroups out of this integrated model, which reached better performances to classify BRC M0 and early-stages samples (see Methods, Fig. 4G, S8A).

This demonstrates the robustness of cancer detection by probing L1PA hypomethylation from plasma DNA with the DIAMOND assay and its ability to detect early stages. Finally, we evaluated if the extent of L1PA hypomethylation was associated with survival. In the validation cohort, we compared patients with high methylation to those with low methylation levels (above or below the median threshold - Fig. S8D) and observed that more pronounced hypomethylation is clearly associated with shorter survival (Fig. S8E). This effect is thought to mainly reflect the tumor burden, which directly impacts the fraction of tumor DNA circulating in the blood, but it remains an interesting non-invasive tumor marker.

DIAMOND data contain signal to infer the tumor burden, which improves cancer detection

When looking at the overall methylation rates, we detected significantly more hypomethylation for more advanced stages of the disease, in particular in metastatic stages compared to localized stages (Fig. 5A, Table S13). However, there was no significant differences between metastatic tumor tissues and primary tissues (Fig. 5B, Table S14), which confirmed that L1PA methylation alteration is an early event in carcinogenesis ^59,60 and also affects early-stage cancers. The differences observed in the blood (Fig. 5A) reflect the ctDNA fraction, which is known to be directly influenced by the tumor burden and the stage of the disease ^1,9. This demonstrates the quantitative potential of DIAMOND, which could help quantify the tumor burden and monitor the disease.

• Download figure
• Open in new tab

Fig 5. DIAMOND data contain signal to infer the tumor burden, which improves cancer detection

A-B. Comparison of the average levels of methylation observed in localized vs metastatic plasma samples (A, p_Early/Adv. = 0.327, p_Adv./Meta. = 3.14e-11, p_Early/Meta. = 2.82e-14; p_{BRC_M0/M+} = 1.3e-18, p_{OVC_M0/M+} = 0.006, p_{GAC_M0/M+} = 0.005, Table S13) or in primary vs metastatic tissues (B, p_OVC = 0.257, p_UVM = 0.820, Table S14). C. L1PA unique hits obtained for 15 breast cancer cell lines compared to the distribution of Cytoscan probes distributed throughout the human genome. D. Genome-wide z-score for all cancer (N = 564) vs healthy plasma samples (N = 120, 63 of the total 183 HDs are used as references to compute the z-score and are not displayed here, p = 1.21e-20). E. Genome-wide z-score by cancer subgroups vs healthy samples. The p-values are computed using Mann–Whitney U test (p_{CRC_M+} = 2.05e-18, p_{BRC_M+} = 1.01e-18, p_{UVM_M+_}= 0.169, p_{LC_M+} = 0.769, p_{OVC_M+_}= 1.84e-11, p_{GAC_M+} = 0.003, p_{BRC_M0} = 5.12e-17, p_{OVC_M0} = 1.09e-12, p_{GAC_M0} = 8.40e-06, Table S15) F. Correlation analysis for genome-wide z-score versus global methylation (r_overall = −0.62, p = 1.25e-69). G. Performances of the 2-step model incorporating CNA with DNA methylation analysis (Classification is done as follow: Proba_Cancer ≤ Threshold C1 AND GZ-score ≤ 121: prediction = Healthy; Proba_Cancer > Threshold C1 OR GZ-score > 121: prediction = Cancer, see methods).

A recognized marker to non-invasively estimate the tumor burden and the fraction of ctDNA is the aneuploidy or copy number alterations (CNA) ^61,62, a hallmark of cancer genomes ⁶³. Given that DIAMOND hits are dispersed throughout the genome (Fig. 1B), we investigated the possibility to use our data to perform CNA analysis. The mFast-SeqS approach had previously used a PCR-based L1PA targeting as a prescreening tool to estimate the fraction of ctDNA ^64,65. This was done on native DNA whereas our data resulted from bisulfite-treated DNA.

We first tested this approach on 15 breast cancer cell lines that were also assessed by CytoScan HD microarrays for aneuploidy. DIAMOND provided an average of 78 000 uniquely mappable reads per cell line, corresponding to around 10 000 L1PA copies precisely located in the genome. These L1PA hits homogeneously overlapped with regions covered by CytoScan probes along the genome (Fig. 5C). We computed z-scores, quantifying copy number alterations, at the level of chromosome arms as previously described (^65,66, see methods) and obtained similar results to those found with CytoScan arrays (Fig. S9A). We observed low alteration scores for the normal-like breast cell line HTERT-HME1 and good correlations between the 2 methods for the majority of the cell lines (Fig. S9B).

Next, we computed genome-wide z-scores in healthy and cancer plasma samples and observed high alteration scores specifically in cancer samples (Fig. 5D). Cancer subgroups z-scores mirrored global hypomethylation profiles (Fig. 5E, 2B and 4D, Table S15), both reflecting tumor burden and ctDNA fractions available. However, global methylation rates and z-scores were only moderately anti-correlated (Fig. 5F), demonstrating that these are partially independent markers that can provide distinct signals (Fig. S9C).

To obtain a final classification labelling each sample as healthy or cancer, we used a 2-step categorization incorporating CNA analysis, which improved cancer detection. We used the probability of the cancer prediction provided by the methylation-based validation models (‘expert’ and ‘all’ models), applying a threshold identified on the discovery cohort, followed by a reclassification of samples which presented a z-score > 121, as cancer. This cut-off value was deduced from a cross-validation applied on C1 (see methods and Fig. S9C). This classifier achieved high sensitivities with specificities between 97-100% for the ‘expert’ models and 93% for the integrated ‘all’ model in 4 distinct cancer types (BRC, CRC, LC, OVC, Fig. 5G) and could be applied as is in the clinic. This was particularly promising for localized breast cancer with a sensitivity of 100% with both model types.

Materials

Preparation of DNA from cell lines and tissues and cfDNA

Isolation of DNA from cell lines and healthy white blood cells (buffy coats) was performed using the QIAamp DNA Mini Kit or QIAamp DNA Blood Mini Kit (Qiagen) according to the manufacturer’s instructions. DNA from cryopreserved and formalin-fixed paraffin embedded (FFPE) tumor tissues was extracted using a classical phenol chloroform protocol and the NucleoSpin® FFPE DNA kit (Macherey Nagel), respectively. cfDNA was extracted from 2 ml of plasma using the automated QIAsymphony Circulating DNA kit (Qiagen), the Maxwell RSC ccfDNA LV plasma kit (Promega) or manual QIAamp circulating nucleic acid kit (Qiagen), according to the manufacturer’s instructions, and eluted in 60 μl, 75 μl or 36 μl, respectively. We verified that the extraction method did not impact our results (Fig. S10A-B). Isolated DNA was quantified by Qubit® 2.0 Fluorometer using dsDNA HS Assay Kit (Thermo Fisher Scientific) according to the manufacturer’s instructions and stored at −20°C until use.

Bisulfite conversion

We used sodium bisulfite-based chemical conversion to achieve base-pair resolution analysis, which is crucial to address methylation levels at single CpG dinucleotides and the co-methylation of multiple CpG sites to determine methylation haplotypes (methylation state of successive CpG sites). Bisulfite treatment of the isolated genomic DNA (up to 200 ng) from the cancer tissues, cancer cell lines and buffy coats was performed using an EZ DNA Methylation-Gold Kit™ (Zymo Research, CA, USA), following the manufacturer’s instructions. Bisulfite treatment of cfDNA (isolated from 2 ml of plasma) was performed using the Zymo EZ DNA Methylation-Lightning Kit™ (Zymo Research, CA, USA), according to the manufacturer’s instructions. Bisulfite-treated DNA was stored at −80°C and further used to build a sequencing library. We have compared the methylation profiles obtained with bisulfite conversion or enzymatic conversion (NEBNext® Enzymatic Methyl-seq Conversion Module, ⁷³) followed by amplification with the DIAMOND targets and deep sequencing. We observed similar methylation profiles (Fig. S11). We have tested various starting quantities of DNA (100ng, 10ng, 5ng) of white blood cells extracted from healthy donors (buffy coats). We have also compared cfDNA (10ng) and DNA from MCF7 breast cancer cell lines (100ng) know to display hypomethylation at L1PA elements. We observed that the methylation profiles along the 30 CpG targets are similar when starting with 100ng (Fig. S11A) but, in our hands, the bisulfite conversion seems more robust with decreasing stating DNA quantities (Fig. S11B-C). We thus observe some differences with the cfDNA samples at 10ng (Fig. S11D). We could still detect very similar hypomethylation profiles in MCF7 breast cancer cell lines (Fig. S11E). This shows that EM-seq could be used to profile L1PA methylation changes but the same conversion should be used throughout the study.

Primer design

Eight primer pairs were designed using the LINE-1 Human Specific (L1HS) consensus sequence from Repbase (Fig. 1A). Although 5’UTR (promoter region) is CpG-rich and common target for methylation quantitation, L1PA copies are frequently 5’-truncated. Therefore, primers were also designed for ORFI and ORFII to target more L1PA elements and improve the sensitivity of our assay. All primers were designed for plus strand of bisulfite converted DNA, using the MethPrimer or PyroMark Softwares. Targeted regions contained 2-7 CpG targets and ranged from 101bp to 150bp, to better capture cfDNA fragments, which have a mean size of 167bp ⁷⁴, (Table S1). Primers were methylation-independent, encompassing 0 to 2 CpGs (none toward the 5′ end), and were degenerated to target both the methylated and unmethylated states. They contained Fluidigm universal CS (common sequence) tags at their 5′ ends. We incorporated a 16 N (random nucleotides) as unique molecular identifiers (UMI) between the target-specific sequence and the CS2 in the reverse primers for signal deconvolution to detect true low frequency alterations and for reducing errors. As LINE-1 hold thousands of copies per genome, a high number of distinct UMIs is essential for unique barcoding of each target molecule. The 16 N stretch between the target-specific sequence and the CS1 in forward primers was used to increase diversity of sequencing libraries and improve sequencing quality. All primers were obtained from Eurogentec (RP-cartridge purification method).

Preparation of targeted bisulfite sequencing libraries

To limit batch effects, sequencing and library preparation batches contained both cancer samples and healthy donors. We also specifically processed the validation cohort C2 with equilibrated number of healthy and cancer samples distributed in 5 experimental batches. Sequencing libraries were prepared using three PCR steps (Fig. S1A): 1) target-specific linear amplification for UMI assignment, 2) target-specific exponential amplification and 3) barcoding PCR for sample identification. Each library was prepared in two individual reactions (due to the overlap of amplicon 2 with other primers), including: I) Multiplex PCR amplification of 7 probes (Amplicon 1, 3, 4, 5, 6, 7, 8), and II) Single PCR amplification of amplicon 2.

UMI assignment for multiplex reaction was performed using Platinum™ Multiplex PCR kit Master Mix (Thermofisher, Life Technologies SAS) in a 25 μL reaction containing 1x Platinum™ Multiplex PCR Master Mix, 0.01-0.06 µM mix of reverse primers and up to 5 ng bisulfite-converted DNA at the following thermocycling conditions: 95°C for 5 min followed by 1 cycle at 95°C for 30 s, 58°C for 90 s, 72°C for 40 s. UMI assignment for single reaction was performed using Hot Star Taq Plus DNA Polymerase (Qiagen) in a 25 μL reaction containing 1x Taq PCR Buffer, 0.65 U Hot Star Taq (5U/μL), 0.2 μM dNTPs, 1.5 mM MgCl2, 0.1 μM amplicon 2 reverse primer, up to 4 ng of bisulfite-converted DNA at the following thermocycling conditions: 95°C for 10 min followed by 1 cycle at 94°C for 60 s, 58°C for 30 s, 72°C for 40 s. To ensure complete removal of the reverse primers and dNTPs, each 25 μL reaction was treated with 50U of Exonuclease I and 10U of FastAP Thermosensitive Alkaline Phosphatase (Thermo Fisher Scientific) at 37°C for 1 h and heat-inactivated at 80°C for 15 min.

Target-specific exponential amplification for multiple reaction was performed using Platinum™ Multiplex PCR kit Master Mix in a 50 μL reaction containing 1x Platinum™ Multiplex PCR Master Mix, 0.01-0.06 µM mix of forward primers, 0.2 μM CS2 reverse primer and 20 μL of purified PCR product at the following thermocycling conditions: 95°C for 5 min followed by 28 cycles at 95°C for 30 s, 58°C for 90 s, 72°C for 30 s followed by a 10 min incubation at 72°C. Target-specific exponential amplification for single reaction was performed using Hot Star Taq Plus DNA Polymerase in a 25 μL reaction containing 1x Taq PCR Buffer, 0.65 U Hot Star Taq (5U/μL), 0.2 μM dNTPs, 1.5 mM MgCl2, 0.2 μM amplicon 2 forward primer, 0.2 μM CS2 reverse primer and 8 ul of purified PCR product at the following thermocycling conditions: 95°C for 10 min, 25 cycles at 94°C for 60 s, 58°C for 30 s, 72°C for 30 s and 10 min at 72°C.

PCR products of multiplex and single reaction were pooled together after quantification by qPCR and purified using Agencourt AMPure XP (Beckman Coulter) at 1.2x ratio according to the manufacturer’s protocol. Purified DNA was eluted in 30 ul of water. Barcoding PCR was performed using universal fluidigm primers. 25 μL of purified pooled PCR product, 1x Phusion HF Buffer, 1 U Phusion Hot Start II DNA Polymerase (Thermo Fisher Scientific), 0.2 μM fluidigm primer, and 0.2 mM dNTPs were mixed in the final volume of 50 μL and amplified with the following conditions: 98 °C for 2 min, followed by 20-25 cycles of 98 °C for 10 s, 62°C for 30 s, and 72 °C for 30 s followed by a 5 min incubation at 72°C. The amplified product was purified by two consecutive AMPure XP steps using 1) a low concentration of AMPure XP beads (0.6x – 0.7x ratio) where the beads containing the larger fragments are discarded and supernatant collected (reverse purification) and 2) higher beads concentration (1.1x – 1.2x ratio) where the beads containing fragments of interest were collected and purified according to the manufacturer’s protocol. Size-selected libraries were eluted in 15 μL of low-EDTA TE buffer. The libraries were quantified with Qubit HS DNA kit (Thermo Fisher Scientific), qualified with nano-electrophoresis (TapeStation, Agilent), and pooled equimolarly for sequencing. Sequencing was performed on Illumina HiSeq rapid run mode or NovaSeq (PE 30bp, 170bp).

Preprocessing of the reads

For each sample, FASTQ files containing raw sequences, composed by the following parts: CS1, forward UMI, forward primer, insert, reverse primer, reverse UMI, and CS2 (Fig. S1A) were first filtered for reads quality (average >Q20 per read) and then demultiplexed (i.e., cut using atropos v1.1.31) using forward and reverse primer sequences. FASTA files were created per primer-set, containing inserts and reverse UMIs for deduplication, as they are unique for each input DNA molecule. Inserts and reverse UMI were then filtered on expected sizes (with a tolerance of ± 5 bases for the inserts). Filtered inserts and UMIs sequences were concatenated and deduplicated using vsearch v2.15.2. Reverse UMIs were then trimmed and resulting inserts from all samples were aggregated into a single FASTA file per primer-set.

Clustering, extraction of representative sequences and global alignment

Using vsearch (with the following parameters: --cluster_fast <inputFasta> --notrunclabels -- fasta_width 0 --iddef 4 --id 0 --qmask none --clusterout_sort --consout <referenceFasta>), a clustering based on sequence identity was applied to each FASTA file, or a subset of 20 million reads randomly chosen if a given file comprised more. The 10 largest clusters’ representative sequences were isolated in separate files. Using mafft v7.508 (with the following parameters: -- globalpair --maxiterate 1000), the 10 representative sequences were aligned pairwise resulting in a reference database for each primer-set. Lastly, using mothur v1.48.0 (with the following parameters: #align.seqs(candidate=<inputFasta>, template=<referenceFasta>, align=needleman, match=1, mismatch=-1, gapopen=-1, gapextend=0)) on each primer-set FASTA file, all sequences from all samples were aligned to the corresponding reference.

CG calling, methylation levels and haplotypes extraction

To call CpG dinucleotides of interest, a sliding window of 2 bp was used on all aligned sequences to determine the distribution of dinucleotides along each amplicon target. A first threshold of ≥ 20% of CG/TG dinucleotides was used to select potential CpG sites. A second threshold was applied to eliminate dinucleotide with ≥95% TG and select position with at least 5% methylation rate. From the aligned sequences, the patterns of methylation were extracted and compiled into either average levels of methylation at each previously identified CpG sites, or proportions of methylation haplotypes for each sample.

Machine learning-based classification models

The resulting data (represented as average levels of methylation per CpG site or proportions of methylation haplotypes or both) were used to do supervised learning of statistical models using the random forest classifier algorithm ⁷⁵ from Python package scikit-learn ⁷⁶, with the following hyperparameters: n_estimators=300, criterion=’gini’, max_depth=None, min_samples_split=2, min_samples_leaf=1, min_weight_fraction_leaf=0.0, max_features=’sqrt’, max_leaf_nodes=None, min_impurity_decrease=0.0, bootstrap=True, oob_score=False, warm_start=False, class_weight=None, ccp_alpha=0.0, max_samples=None.

The rational for choosing random forest over other learning methods was driven by three main factors: 1) it is less prone to overfitting ⁷⁵; 2) it shows excellent performance even when the quantitative relationship between features and observations is biased in favor of the former, such as when using methylation haplotypes data representation ⁷⁷; 3) random forests also inherently return measures of variable importance ⁷⁵, such as mean decrease in impurity, which greatly facilitate the interpretability of model decisions. The features used to train the models were the average levels of methylation per CG site (n=30), the proportions of methylation haplotypes (i.e., the combinatorial of all the possible methylation status of CG sites within a given amplicon, n=372) or both. No additional transformation nor feature selection was performed on the data.

Expert and All cancer models

Model classifications were run 5000 times in order to estimate variance and confidence intervals. For the discovery step, in each run, as many samples from each class were randomly drawn to construct a balanced subset of the data ⁷⁸. The samples from these draws were stratified by class and split into 60% for training, 40% for evaluation. For the validation step, we trained the model on the entire cohort 1 and evaluated it on cohort 2. The true and false positive rates for all possible classification threshold were evaluated at each run, with interpolation to generate an average ROC curve with 95% confidence interval for the 5000 runs. 95% confidence interval has been calculated with the following formula: Μ±z*s/√n with M the average of the variable, z the confidence level (z∼1.96 for 95%CI), s the standard deviation, n the number of samples in the variable.

Blind models

We trained a random forest on haplotypes features, removing one cancer type or subgroup from the training set. The specific cancer type or subgroup is then assessed in the test set. We pooled together the discovery and validation cohorts, training on 2/3 of all the samples-- excluding the cancer type or subgroup to test for-- and testing on the remaining one-third of the samples. The only exception was metastatic gastric cancer: as they are made up of only 3 samples, they were systematically moved to the test set, consequently blinding the model towards GAC M+ samples. We also trained a stacked version with this setup (see below).

Stacked machine learning model

We developed a model referred to in the paper as “Stack”. This model uses one Random Forest model for each combination of cancer type and metastatic status, known as the “expert sub-model.” Each expert sub-model was trained on one-third of the healthy plasma samples and one-third of the samples matching the cancer subgroup of interest (cancer type and dissemination status). These expert sub-models were then combined into a Random Forest Stack model, which uses both the haplotype features and the probabilities output by each expert sub-model. The final Random Forest Stack model was trained on an additional one-third of the healthy and cancer plasma samples and tested on the remaining samples (which represent one third of the healthy samples and one third of each subgroup).

Mutation screening for ovarian cancer samples

Ovarian tumor genotyping was performed using the TIGER panel previously developed by Institut Curie⁷⁹, which targets 78 genes or using ‘custom NGS’ with amplicons targeting Tp53 and TSC2, which are the 2 most frequently mutated genes in this group of patients, with TruSeq library constructions for low input material (dual strand technology). After mutations identification, ovarian cancer plasmas were genotyped using custom NGS or ddPCR as previously done¹⁴. Sequencing was performed on a MiSeq V3-150 (25M) with Paired End 75bd protocol. Results are presented in Table S8.

Whole genome bisulfite sequencing analysis

To see if we could retrieve cancer-associated L1PA hypomethylation in others plasma studies, we analyzed data sets from 2 recent studies profiling cell-free DNA methylation with Whole Genome Bisulfite Sequencing (WGBS) in healthy and cancer patients (Liu et al. 2024⁵³; Gao et al. 2022⁵⁴). Liu et al. analyzed methylation at 75 617 CpGs in the whole genome (estimated with Bis-SNP after mapping with Bismark) of 17 healthy donors and 31 cancer patients. We extracted the methylation levels at CpG residing within L1PA families hit by DIAMOND (mostly L1HS-L1PA10, see Fig. 1C). To do so, we identified CpG dinucleotides covered by L1HS-L1PA10 elements, and their position within L1 based on the L1HS consensus sequence, and we selected the DIAMOND copies using Bedtools⁸⁰. Gao et al., provided FASTQ files data for 123 breast cancer and 40 healthy patients. In order to avoid breast-cancer subtype effect and age effect, we subsampled the data to generate an age matched cohort with 16 healthy and 15 HR+HER2-M+ patients. We first merged the paired-end reads using Fastp⁸¹ and mapped them on the hg38 reference genome using bismark. At this step, we followed two different approaches to obtain the methylation level at the 30 L1PA CpGs that we studied with DIAMOND. 1) Using bismark extractor, we retrieved the methylation levels of all CpGs covered and intersected it with the dinucleotides covered by DIAMOND using Bedtools. Results produced with this approach are subsequently called “Gao et al (2022) Bismark”. 2) Using mothur⁸², we aligned the reads that mapped on L1HS-L1PA10 elements on an aggregate of native and converted (CG converted to TG) reference sequences of L1HS. All read that mapped with a score lower than 51% of similarity (score computed by mothur) were rejected (this threshold was established from the percentile at 99% of similarity score distribution from 1000 random sequences 167 bp long).

Results generated by this approach are subsequently called “Gao et al (2022) Mothur”. Finally, for “Gao et al (2022) Bismark” and “Gao et al (2022) Mothur “, we refined the results to only take in account the copy of L1HS-L1PA10 targeted by DIAMOND. Further refinements were done to only consider L1HS-L1PA3 copies with the second approach and is referred as “L1HS-L1PA3”.

Survival analysis

Survival analysis has been performed with survival and survminer R packages - https://cran.r-project.org/web/packages/survival/citation.html, doi:10.32614/CRAN.package.survminer.

Copy number alterations analysis

Cytoscan HD microarrays: 250 ng of gDNA from 15 breast cell lines (1 normal-like: HTERT-HME1 and 14 cancer cell lines: MDA-MB231, MDA-MB453, HCC1569, BT20, HCC1954, HCC38, MDA-MB361, ZR 75.1, MDA-MB157, MCF7, SKBR3, HCC202, HCC70, BT474) were characterized using Affimetrix/Thermo Cytoscan HD microarrays at the Genomics facility of Institut Curie to profile aneuploidy. To compare with the z-score by chromosome arm, we calculated the mean of Weighted Log2 combining probes by chromosome arms.

DIAMOND CNA: 1) Z-score calculation: prepossessed reads were uniquely mapped on hg38 genome using Bismark (version 0.23.1). As in Belic et al. 2015 ⁶⁶, only the reads with an alignment score > 15 were kept. Resulting reads from all amplicons (excluding #2 and #3) were merged and normalized number of reads per chromosome arm (excluding sexual chromosomes X and Y) per sample were calculated with R. Next, the amplifications/deletions score was computed using the following formula: with i = a given chromosome arm, n = a given sample and controls = a set of reference samples (white blood cells from 10 healthy reference samples for the cell lines, 63 healthy plasma from C1 as a reference for cancer and healthy plasma samples). Genome-wide z-scores were computed by summing the squared z-scores of all chromosome arms. 2) Z-score threshold identification: to identify altered versus normal z-scores, we performed 5-fold cross validation of simple cutoff classification model on the discovery cohort (N_Healthy = 60, N_Cancer = 350) using the genome-wide z-score and calculated the threshold that maximize the sensitivity at 100% specificity.

2-step classification for sample labelling

First, we selected the threshold for the probability of the cancer prediction (Proba_Cancer) on the discovery cohort maximizing the sensitivity for a 99% specificity, per ‘cancer-type’ model. We applied this threshold on the Proba_Cancer computed with the validation models and reclassified samples which presented a z-score > 121, as cancer (Proba_Cancer ≤ Threshold C1 AND GZ-score ≤ 121: prediction = Healthy; Proba_Cancer > Threshold C1 OR GZ-score > 121: prediction = Cancer) see Tables S16-22.