Integrated multiomics signatures to optimize the accurate diagnosis of lung cancer

Lung cancer remains a major health concern worldwide^1,2, underscoring the urgent need for effective strategies that can reduce lung cancer-related mortality. Early diagnoses are vital for improving clinical outcomes by facilitating timely interventions³; however, the acknowledged noninvasive lung cancer diagnosis method, which uses low-dose computed tomography (LDCT) imaging^4,5,6, has a high false-positive rate, limiting its clinical applicability⁷. Hence, there is still an unmet clinical need for noninvasive, accurate, and cost-effective tests to facilitate the early diagnoses of lung cancer.

The identifications of 5-methylcytosine (5mC) and 5-hydroxymethylcytosine (5hmC) modifications from circulating cell-free DNA (cfDNA) have been proposed as promising methods for the early diagnoses of lung cancer^{8,9,10,11,12,13,14}. However, most published studies evaluating the diagnostic value of 5mC and 5hmC biomarkers have limitations, including small sample size of lung cancer cases, the absence of independent external validation sets, and restricted reproducibility due to the use of various sequencing methods. Besides, the single-dimensional information would restrict the diagnostic accuracy, leading to a high lung cancer-related mortality. The potential of multiomics data in highlighting the importance of comprehensive molecular profiling of cancer has been suggested in recent studies and offers valuable insights into cancer pathogenesis and therapeutic strategies^15,16,17. Hence, constructing a multiomics biomarker in multiple large, well-characterized cohorts holds promise in improving the diagnostic accuracy and clinical applicability.

Here, the plasma cfDNA end motif profiles in 5mC- and 5hmC-enriched regions were analyzed based on a multicenter, prospective, large cohort of indeterminate pulmonary nodules (IPLs). An integrated multiomics signature was then developed by integrating end motif, CT-based radiomic and clinical factors to optimize the lung cancer diagnosis process. To our knowledge, the potential value of the multiomics model for lung cancer detection is underexplored, and by conducting this study, we aimed to provide a more accurate, robust, noninvasive, and cost-effective test for diagnosing lung cancer and promoting the precision management for patients with IPLs.

Participant characteristics

The clinical characteristics of all 2032 participants (Fig. 1) were summarized in Table 1, with detailed histopathological results listed in Supplementary Data 1. Among 1718 participants from Shanghai Pulmonary Hospital, 1019 (59.3%) were females and 699 (40.7%) were males with a mean age of 55.9 years. The majority of the IPLs manifested as subsolid nodules, including 22.4% pure-GGOs and 42.3% part-solid nodules with a mean (±SD) solid component size of 15.9 ( ± 13.5) mm. After randomization, the training, validation and internal test sets included 1030, 344, and 344 participants, respectively. The external test set included 314 participants, consisting of 197 (62.7%) females and 117 (37.3%) males with a mean age of 56.2 years. The clinical characteristics of the internal test set and the external test set were compared, and the participants on the external test set had fewer pure-solid nodules (p = 0.011) and smaller solid component sizes (p = 0.012).

A participants selecting from Shanghai Pulmonary Hospital; B participants selecting from external 4 hospitals.

Full size image

Fragmentomic model establishment

The median final numbers of mapped fragments (paired reads) derived from the 5mC- and 5hmC-sequencing data were 12.7 million and 14.1 million, respectively, and no systematic bias or obvious batch effect induced by hospital regions were detected both for the 5mC- and 5hmC- sequencing data (Supplementary Fig. 1). The 4-mer (bp) end-motifs profiles were first identified from 5mC- and 5hmC- sequencing data, and then, extended to 6 bp end motifs profiles (Fig. 2A). After feature dimensionality reduction and selection, 27 4 bp end motifs and 11 6 bp end motifs from the 5mC-sequencing data, 11 4 bp end motifs and 11 6 bp end motif from the 5hmC-sequencing data were identified. The numbers of remaining features and the feature names observed after the processing steps are detailed in Supplementary Table 1. Hierarchical clustering analyses of the end motifs selected from 5mC-sequencing data is illustrated in Supplementary Fig. 2A for the training set and Fig. 2B for the rest sets, while the hierarchical clustering analyses of the end motifs selected from 5hmC-sequencing data is illustrated in Supplementary Fig. 2B, C.

A Schematic of the process for determining the first 6-nucleotide sequence (i.e., a 6-mer end motif) on each 5’ fragment end of cfDNA relative to the hg19 reference genome; B Hierarchical clustering analyses of the selected 6 bp end motifs derived from 5mC-sequencing data; Receiver operating characteristic analyses of the epigenomic models on the validation set C internal test set D and external test set (E); F Bar charts showing the TFs identified by the 6-mer end motifs selected from 5mC-sequencing data; G Bar charts showing the regulatory target genes by these identified TFs; H The top 15 most enriched GO terms based on the target genes. The area under the receiver operating characteristic curves are compared via the DeLong’s test. All statistical tests were two-sided, with p < 0.05 indicative of a statistically significant difference. 4bp-5mC, the model established by the 4-mer end motifs selected from 5mC-sequencing data; 6bp-5mC, the model established by the 6-mer end motifs selected 5mC-sequencing data; 4bp-5hmC, the model established by the 4-mer end motifs selected from 5hmC-sequencing data; 6bp-5hmC, the model established by the 6-mer end motifs selected from 5hmC-sequencing data; 5mC, 5-methylcytosine; 5hmC, 5-hydroxymethylcytosine; ROCs, receiver operating characteristic curve; AUC area under the ROCs curve, CI confidence interval, Sens sensitivity, Spec specificity, PPV positive predictive value, NPV negative predictive value, Accur accuracy, TFs transcription factors, GO gene ontology. Source data are provided as a Source Data file.

Full size image

The predictive models were then established accordingly, with hyperparameter results summarized in Supplementary Table 2. As illustrated in Fig. 2C–E, for 5mC-sequencing data, the model developed by 6 bp end motifs (6bp-5mC model) produced higher AUCs of 0.805, 0.826 on the validation set and external test set, respectively, outperforming the model developed by 4 bp end motifs (4bp-5mC model), although no significant differences (p > 0.05 for all). The AUC of the 6bp-5mC model on the internal test set was slightly lower than that of the 4bp-5mC model but still did not significantly differ (p = 0.100). Also, the 6bp-5mC model maintained higher specificity and PPV than those of the 4bp-5mC model across all datasets, and higher accuracy on the external test set. For 5hmC-sequencing data, the predictive model developed with the 6 bp end motifs (6bp-5hmC model) achieved greater AUCs than did the model developed with the 4 bp end motifs (4bp-5hmC model) on the external test set (p = 0.399; Fig. 2E). The 6bp-5hmC model also exhibited greater accuracy. Hence, we selected the 6bp-5mC model and the 6bp-5hmC model for subsequent analyses. Upon downsampling the fragments to 25%, 50%, and 75%, the AUC values gradually increased, and the 6 bp end motifs-based models maintained the best performance and robustness (Supplementary Fig. 2D, E), revealing their appropriate applications. The hyperparameter of models derived from 25%, 50%, and 75% sequencing fragments are summarized in Supplementary Table 2.

For the functional analysis of 6 bp end motifs incorporated into the 6bp-5mC model, 131 transcription factors (TFs) were identified (Fig. 2F), with the top 20 TFs illustrated in Fig. 2G. The most enriched biological processes were cancer-associated pathways (Fig. 2H), such as the epithelial cell proliferation. The similar analysis was performed for 6 bp end motifs incorporated into the 6bp-5hmC model (Supplementary Fig. 2F, G), while the most enriched pathways illustrated in Supplementary Fig. 2H. The impacts of each end motif on the 6bp-5mC model and 6bp-5hmC model predictions for malignancy risks of IPLs are visualized in Supplementary Fig. 3.

Radiomic model establishment

The construction schemes of the radiomic models are illustrated in Fig. 3A, B. The trainable parameters of the DL-radiomics model are described in Supplementary Method 5, 6, while the hyperparameter results for the C-radiomics model are listed in Supplementary Table 2. The DL-radiomics model achieved AUCs of 0.842, 0.842 and 0.855 on the validation set (Fig. 3C), internal test set (Fig. 3D) and external test set (Fig. 3E), respectively, outperformed the C-radiomics model, especially on the external test set (p = 0.012), revealing the superiority of the deep learning technology. The attribution scores of all 64 deep learning-based radiomics features obtained from the training set is illustrated in Fig. 3F, and the DL-feature 07, 52, and 25 had the largest contributions to the DL-radiomics model’s output (Fig. 3G). Nevertheless, the definite meanings of DL-features are still difficult to explain. To further enhance the interpretability of the DL-radiomics model, the biological basis of it was evaluated. A total of 116 patients from Radiogenomics dataset were identified, and classified as low-risk (DL-radiomics model score<47.8, the cutoff value calculated based on the maximum Youden index on the training set) or high-risk (DL-radiomics model score≥47.8). As illustrated in Fig. 3H, 218 differentially expressed genes, consisting of 89 upregulated genes and 129 downregulated genes, were found between the low and high-risk subgroups. The tumorigenesis-related pathways, such as the HIF-1 signaling pathways and transcriptional misregulation in cancer, were significantly upregulated (Fig. 3I), whereas the pathways related to tumor immunity, such as the B-cell receptor signaling pathway and the T-cell receptor signaling pathway, were significantly downregulated (Fig. 3J). Further, high-risk patients classified by the DL-radiomics model were tend to show fewer activated CD8 T cells, fewer myeloid-derived suppressor cells (MDSCs), and fewer T follicular helper cells in their tumor microenvironments, although the difference were not significant (Fig. 3K).

Construction procedures utilized for the DL-radiomics model A and the C-radiomics model B; Receiver operating characteristics analyses of the DL-radiomics model and C-radiomics model on the validation set C, the internal test set D, and the external test set E; F, The attribution score distribution of each deep learning-based radiomics feature on the training set; G The impact of the deep learning-based radiomics features included in the DL-radiomics models on the risk probabilities output for the training set; H Volcano plot for differentially expressed genes between the low-risk and high-risk subgroups classified by the DL-radiomics model; I GSEA conducted for upregulated genes; J GSEA conducted for downregulated genes; K ssGSEA conducted for the low-risk and high-risk subgroups classified by the DL-radiomics model, The center of the box denotes the 50th percentile, the bounds of the box contain the 25th to 75th percentiles, the whiskers mark the maximum and minimum values, and values beyond these upper and lower whiskers are considered outliers and are marked with dots. n = 116 biologically independent samples were analyzed. All statistical tests were two-sided, with p < 0.05 indicative of a statistically significant difference. DL-radiomics, the deep learning-based radiomics model score, while C-radiomics referred to the classic radiomics model. ROCs, receiver operating characteristic analyses; AUC, area under the ROCs curve; CI, confidence interval; DL, deep learning; KEGG, Kyoto Encyclopedia of Genes and Genomes; MDSC, myeloid derived suppressor cell; GSEA, gene set enrichment analyses; ssGSEA, single sample gene set enrichment analyses. Source data are provided as a Source Data file.

Full size image

Multiomics model establishment

The multivariable logistic regression analyses revealed that age, radiological solid component size, the DL-radiomics model score, the 6bp-5mC model score, and the 6bp-5hmC model score were significant predictors for malignancy risk of IPLs (p < 0.001 for all; Table 2). All single-omics model score (DL-radiomics, 6bp-5mC, and 6bp-5hmC) showed weak correlations with clinical factors, with the absolute value of Pearson correlation coefficients≤0.10 for all (Fig. 4A–C), suggesting that both radiomic and fragmentomic features provide additional information to clinical data for lung cancer detection. A clinical model and 4 multiomics models were then established, with the parameter information detailed in Supplementary Table 2. The clinic-RadmC model, which integrated the clinical factors, DL-radiomics model with the 6bp-5mC model, yielded significantly improved AUCs of 0.883 on the validation set (Fig. 4D), 0.889 on the internal test set (Fig. 4E), and 0.923 on the external test set (Fig. 4F), outperforming the models who only integrated clinical factors with the DL-radiomics model, or with the 6bp-5mC model (p < 0.050 for all; Supplementary Table 3). Interestingly, we found that after adding the diagnostic information from the 5hmC-sequencing data, the clinic-Rad(h)mC model did not significantly outperform the clinic-RadmC model (p > 0.050 for all; Supplementary Table 3). Based on the cutoff value of 78.0, 77.9, 80.5, 82.7, and 82.7 for the clinical, clinic-Radiomics, clinic-mC, clinic-RadmC and clinic-Rad(h)mC model, comparable performance metrics were observed between the clinic-RadmC model and the clinic-Rad(h)mC model (Fig. 4D–F). A calibration curve analysis revealed that the predicted probabilities of the clinic-RadmC model and the actual malignancy risk probabilities had good agreement (Hosmer-Lemeshow test, p > 0.050 for all datasets; Supplementary Fig. 4A–C). DCA suggested that the clinic-RadmC model achieved better net clinical benefits than did the clinic-Radiomics, clinic-mC and clinical model, but similar clinical benefits with the clinic-Rad(h)mC model (Supplementary Fig. 4D–F).

Pearson correlation coefficients analyses for the features included in the multiomics model on validation set A, internal test set B, and external test set C. Receiver operating characteristic analyses and performance metrics for the models on validation set D, internal test set E, and external test set F. All statistical tests were two-sided, with p < 0.05 indicative of a statistically significant difference, and *** denotes p < 0.001; ** denotes p < 0.01; * denotes p < 0.05. Size, the radiological solid component size of pulmonary nodules; Radiomics, the deep learning-based radiomics model score; 6bp-5mC, the model score established by the 6-mer end motifs selected from 5mC-equencing data; 6bp-5hmC, the model score established by the 6-mer end motifs selected from 5hmC-sequencing data; clinical, the model established by the age and radiological solid component size of pulmonary nodule; clinic-Radiomics, the model established by combining clinical variables with the DL-radiomics model sore; clinic-mC, the model established by combining clinical variables with the 6bp-5mC model score; clinic-RadmC, the model established by combining clinical variables, DL-radiomics model score with the 6bp-5mC model score; clinic-Rad(h)mC, the model established by combining clinical variables, the DL-radiomics model score, the 6bp-5mC model score with the 6bp-5hmC model score. ROCs, receiver operating characteristics analyses; AUC area under the ROCs curve, Sens sensitivity, Spec specificity, PPV positive predictive value, NPV negative predictive value, Accur accuracy. Source data are provided as a Source Data file.

Full size image

In further comparisons with single-omics models, the clinic-RadmC outperformed the 6 bp end motifs-based model did (p < 0.010 for all datasets; Supplementary Fig. 5A–C). It kept the better performance than the DL-radiomics model, especially on the external test set (p < 0.001; Supplementary Fig. 5C). DCA suggested that the net clinical benefits derived by the clinic-RadmC model was better than all single-omics models (Supplementary Fig. 5D–F), as well as the IDI and NRI analyses (Supplementary Table 4). The impact of each feature on the malignancy probability output by the clinic-RadmC model is illustrated in Fig. 5A–C, with the DL-radiomics model score harboring the largest contribution, followed by the 6bp-5mC score. Additionally, 2 example participants with similar radiological characteristics, the SHAP maps of which are illustrated in Fig. 5D, E, could be classified into different subgroups with distinct malignancy risk levels based on the probability output by the clinic-RadmC model, enabling to select high-risk participants for early intervention.

SHAP analyses results obtained for ranking the impact of the continuous features included in the clinic-RadmC models on the risk probabilities output for the validation set A, internal test set B, and external test set C; 2 example participants with similar clinic-radiological characteristics in the low-risk D and high-risk subgroups E, respectively. Each participant was represented by a single dot on each feature flow. The horizontal position of the dot was determined by the SHAP value of that feature, and dots were accumulated along each feature row to show density values. DL-radiomics, deep learning-based radiomic model; 6bp-5mC, the model established by the 6-mer end motifs selected from 5mC-sequencing data; solid size, the radiological solid component size of pulmonary nodules; SHAP, Shapley additive explanations. Source data are provided as a Source Data file.

Full size image

Performance evaluation of the models in various clinical scenarios

The diagnostic performance of the clinic-RadmC model was further evaluated after adjusting for clinical and radiological characteristics. The clinic-RadmC model yielded significantly improved AUCs over clinical model for participants aged≤50 years on test sets (p < 0.050 for all; Fig. 6A & Supplementary Table 5). For participants aged>50 years who were identified as high-risk population for lung cancer¹⁸, the clinic-RadmC model outperformed the other models (p < 0.050 for all; p = 0.070 for vs. clinic-Radiomics on the internal test set) but performed comparably with the clinic-Rad(h)mC model (p > 0.050 for all; Fig. 6B & Supplementary Table 5).

Predictive performance achieved by the models in subgroups stratified by age A, B radiological image C–E, and nodule size F–H. Pure-GGO comprised nodules with only GGO, and part-solid nodules consisted of GGOs and solid components, whereas pure-solid nodules had only solid components without GGOs. Subcentimeter pulmonary nodules were defined as the nodules with solid component size≤10 mm, and large nodules were defined as those with 15 mm≤solid component size≤30 mm, whereas pulmonary massed were defined as those with solid component size>30 mm. 6bp-5mC, the model established by the 6-mer end motifs selected from 5mC-sequencing data; DL-radiomics, the deep learning-based radiomic model score; clinic-mC, the model established by combining clinical variables with the 6bp-5mC model score; clinic-Radiomics, the model established by combining clinical variables with the DL-radiomics model sore; clinic-Rad(h)mC, the model established by combining clinical variables, the DL-radiomics model score, the 6bp-5mC model score with the 6bp-5hmC model score; clinic-RadmC, the model established by combining clinical variables, the DL-radiomics model score with the 6bp-5mC model score. AUCs areas under the receiver operating characteristics curves; GGO ground-glass opacity. Source data are provided as a Source Data file.

Full size image

For subsolid nodules, the clinic-RadmC model did not perform better in pure-GGOs (p > 0.050 for all; p = 0.005 for vs. the DL-radiomics model on the internal test set; Fig. 6C and Supplementary Table 6), but performed better than the other models on the external test set for part-solid nodules (p < 0.050 for all; Fig. 6D & Supplementary Table 6). However, the clinic-RadmC model exhibited excellent performance for pure-solid nodules (Fig. 6E), outperforming the single-omics models (p < 0.050 for all; Supplementary Table 6). Our findings revealed the increasing applicability of the clinic-RadmC model along with the increase in solid components size of IPLs.

Diagnosing subcentimeter pulmonary nodules remains challenging^19,20 in clinical practice, and compared with the clinical model, the clinic-RadmC model yielded higher AUCs of 0.776, and 0.937 for predicting the malignancy risks of subcentimeter IPLs on test sets (p < 0.001; Fig. 6F & Supplementary Table 7). In addition, the current guidelines still rely on the clinical risk models, and large IPLs (size≥15 mm) satisfy the 10% threshold for investigation, whereas the British Thoracic Society Guidelines suggest a broad range of options. The clinic-RadmC model achieved improved AUCs of 0.988, and 0.936 on the internal and the external test set, when compared with the 6bp-5mC model (p = 0.002 for all), the clinic-mC model (p = 0.003 for all), and the clinical model did (p < 0.001 for all; Fig. 6G & Supplementary Table 7). Similar trends concerning the superiority of the clinic-RadmC model were observed for pulmonary massed with size>30 mm (Fig. 6H & Supplementary Table 7).

Support for guiding individual management process

The impacts of the models on the malignancy risk reclassification for IPLs were evaluated on the combined test sets (n = 658 consisting of 137 benign and 521 malignant samples). Based on the cutoff value of 82.7 on the training set, the clinic-RadmC model correctly classified 112 samples as low-risk and 456 samples as high-risk, leading to a 78.2% reduction in unnecessary invasive procedures (112/137; Fig. 7A). When compared with the clinical model (Fig. 7B), the clinic-RadmC model could correctly reclassified 48 benign samples as low-risk, and 202 malignant samples as high-risk, reducing the unnecessary invasive procedures for benign diseases by 35.0% (48/137), and avoiding the delayed treatment for cancer by 38.8% (202/521). When compared with single-omics model consisted of the DL-radiomics and 6bp-5mC model, the clinic-RadmC model could provide a reduction of unnecessary invasive procedures for benign cases by 13.1% ~ 16.1% (18/137; 22/137), and provide avoidance of potential delayed treatment for lung cancer by 3.1% ~ 20.3% (16/521; 106/521; Fig. 7C, D). The superior reclassification performances of the clinic-RadmC model to the clinic-Radiomics or clinic-mC model were also observed, and it could reduce the unnecessary invasive procedures for benign IPLs by 10.9% ~ 21.9%, and prevent delayed treatments for lung cancer by 3.8% ~ 13.8% (Fig. 7E, F).

A Confusion matrices illustrating the predicted outcomes generated by the multiomics model in comparison with the actual outcomes, as well as between the multiomics model and the clinical model B the DL-radiomics model C the 6bp-5mC model D the clinic-Radiomics model E and the clinic-mC model F with emphasis placed on the patients ruled in and ruled out. The dotted lines demarcate the corresponding cutoff values of the different models. The number labeled with ^* refer to cancer cases misclassified as low-risk samples by the x-axis model but correctly reclassified as high-risk samples by the multiomics model on y-axis, whereas the number labeled with ^# refer to benign cases misclassified as high-risk samples by the x-axis model but correctly reclassified as low-risk samples by the multiomics model on the y-axis. Source data are provided as a Source Data file.

Full size image

In this multi-institutional study involving a large population, we established the clinic-RadmC model by integrating clinical, fragmentomic with radiomic features and demonstrated that this multiomics model provided additional diagnostic information for lung cancer. We further revealed its superior performance to single-omics models and clinical model in predicting the malignancy risks of IPLs, with an execution time of less than one second, offering a fast, effective and noninvasive method for lung cancer detection and facilitating individual management.

Currently, chest LDCT is still a low-financial-cost and highly-sensitive method for detecting lung cancer, but 96% of IPLs are diagnosed as benign after follow-up scans; as many as 20% of IPLs are diagnosed as benign after surgical resection, and 38% after biopsy^21,22. In previously published studies, radiomic analyses enable to noninvasively and longitudinally capture the radiographic phenotype and characterize the underlying pathophysiology of a tumor, providing an improved disease diagnosis accuracy^23,24. A dizzying array of possibilities, such as circulating tumor DNA, DNA methylation, have also been investigated as potential diagnostic biomarkers for IPLs^25,26, allowing for the optimization of early cancer diagnoses²⁶. The fragmentation patter of cfDNA represents another class of liquid biopsy-based biomarkers in oncology, and hold promise for cancer detection^{11,27,28,29,30}. Among the emergent class of cfDNA fragment features, end motif-based approach act as an effective method for identifying various cancers^9,31,32,33, outperforming other cfDNA fragment metrics³⁴. End motif also has the advantage that only a small number of DNA molecules required to achieve its plateau of performance comparing to other tumor markers such as copy number aberration and tumor-derived single-nucleotide variants^27,31. Recent study has indicated that cfDNA motif-based method held promise in stage I lung adenocarcinoma detection⁹. These preliminary findings suggested the end motifs-based approach outperformed other cfDNA fragment metrics in identifying cancers, shedding light on the diagnosis of lung cancer at early stage in our study. In this context, it is conceivable that molecular biomarkers combined with radiomic might improve the cancer detection process through the noninvasive characterization of cancer biology. In 2021, Kammer et al.³⁵ developed a biomarker model by combining clinical variables, blood-based biomarkers with radiomic features, which first demonstrated the notable strength of combined biomarkers in hastening the early diagnosis of lung cancer. In 2023, PulmonSeek Plus, established by combining clinical data, CT imaging data, and cfDNA methylation biomarkers, achieved excellent accuracy¹⁷, suggesting the superiority of the integrative model. Nevertheless, the diagnostic value of plasma cfDNA end motif profiles in 5mC- or 5hmC-enriched regions specifically for lung cancer is still underexplored, and studies evaluating the multiomics signatures for diagnosing lung cancer are mostly limited by relatively small sample sizes.

Hence, our study was conducted to establish a multiomics model by integrating clinical, radiomic and end motifs features from 5mC-sequencing data, aiming to provide a more accurate, stable, cost-effective and noninvasive tool for optimizing lung cancer diagnosis. Several strengths of our study might be emphasized. First, 2032 participants from 5 different clinical centers were prospectively enrolled, representing the largest sample size reported in this field to date. Second, a deep learning technique was used to extract radiomic features and establish the radiomic model. The construction of a classic radiomic model still requires the ROIs in CT images to be manually segmented in a layer-by-layer manner. Compared with the classic radiomic, deep learning techniques reduce the burden of data processing by extracting radiomic features directly from the labeled ROIs and can also achieve better and more stable performance. In addition, the existing DL-based radiomics model reported in previously published studies yielded a pooled AUC of 0.88 ~ 0.94, with a sensitivity of 63.0% ~ 84.0% and a specificity of 83% ~ 90%^{36,37,38,39,40,41}. The DL-radiomics models in our study achieved AUCs around 0.85 and maintained comparable performance with that of previously DL-based radiomics models, indicating the appropriate application of the DL-radiomics model for further analyses in our study. Third, our study added the information to support the diagnostic value of the cfDNA 6 bp end motifs profiles in 5mC- and 5hmC-enriched regions for lung cancer. Mechanistically, the cfDNA end motifs profiles might be globally affected by the changes in chromatin accessibility and nuclease activity⁴², and chromatin accessibility is important for the fragmentation of cfDNA, which contributes to their preferred ends⁴³. Accordingly, we assume that the abnormal epigenetic changes in cancer genomes might influence the chromatin accessibility and DNA endonucleases activity, leading to the changes of cfDNA end motifs^9,27, and the 5mC and 5hmC modifications would exert different effects on cfDNA end motif profile, since they are 2 distinct epigenetic marks in human genomes. We investigated the performance of 4 bp and 6 bp end motifs derived from 5mC-sequencing data, and found that both of them yielded good performance (>0.80) in distinguishing malignant and benign nodules. Compared with similar assays previously reported^{9,10,11,30,32,33,34}, both 4bp- and 6bp-5mC model kept the comparable and robust performance for lung cancer detection. Besides, the incorporation of the multidimensional epigenomics features into multiomics models might comprehensively reveal more epigenetic characteristics, such as the previously reported efficacy of multidimensional “omics” features^11,44. 5hmC biomarkers have been used as notable tools for lung cancer detection in our previous studies^13,14. However, in this study, after methylation information has been incorporated, adding a diagnostic score form 5hmC-sequencing data could not provide additional improvement for predicting malignancy risk of IPLs. DCA, NRI and IDI analyses further validated the similar benefits derived from the clinic-RadmC model and the clinic-Rad(h)mC model. These finding suggested that it is not the case that the more -omics information, the more accurate the lung cancer diagnosis, and the combination of information derived from clinical, radiomic, and cfDNA 5mC-sequencing data could provide sufficient performance, serving as a more cost-effective method for diagnosing lung cancer in clinical practice. Fourth, machine learning methods, such as, SVM, deep learning, and logistic analyses, were used to construct models in our study, and the interpretability analyses of the predictive models were conducted in our study. Translating AI research to the optimization of the clinical workflow is the ultimate goal of this state-of-the-art technology, which is also the benchmark indicator for measuring the clinical value of an AI model. We conducted a functional analysis for the end-motifs, exploring the importance levels of DL-based radiomic features via DeepLIFT, and ranked the importance of each feature in the clinic-RadmC model via SHAP. We further explored the biological basis of the DL-radiomics model via the published Radiogenomics dataset. These analyses might help to enhance the credibility and applicability of our predictive models in clinical practice.

We further compared the predictive performance of the clinic-RadmC model with that of a similar model, namely PulmoSeek Plus. First, despite the similar AUC, the samples size used to establish the clinic-RadmC model included many more than 543 samples used to establish the PulmoSeek Plus model, ensuring the robustness of the clinic-RadmC model. Second, the cfDNA methylation biomarkers in the PulmoSeek Plus model were identified via high throughput DNA bisulfite sequencing technology^45,46, which might be complex and costly, and might limit the generalizability of this approach. Besides, the CT image features included in the PulmoSeek Plus model were evaluated manually, which might impose a greater workload on radiologists. In contrast, we applied the low-depth cfDNA sequencing technology to identify fragment patterns on 5mC- and 5hmC-enriched regions, and radiomic features were identified directly via an “end-to-end” deep learning algorithm, which improved the clinical cost-effectiveness and convenience of the clinic-RadmC model. Third, an explainability analysis of the clinic-RadmC model were performed to elucidate the mechanism underlying the prediction process from the perspective of genomic biology, and could provide insight into its diagnostic mechanism and enable broad adoption of the model. Similar explorations were rarely conducted on validated models in previously public studies. Accordingly, the clinic-RadmC model could serve as a more efficient, credible, stable, and cost-effective signature in terms of lung cancer diagnoses.

We also evaluated the clinical applicability of our established clinic-RadmC model by investigating the reclassification performance of it, since an ideal diagnostic tool for the clinical setting would be expected to facilitate individual management of IPLs. Our findings revealed that the clinic-RadmC model could not only achieve high sensitivity of 0.895 and high specificity of 0.806 in diagnosing lung cancer on the external test set, but also reduce the unnecessary invasive procedures by 10.9% ~ 35% for benign IPLs, and avoid the delayed treatment by 3.1% ~ 38.8% for lung cancer. Hence, this multiomics assay allow to use easy-accessible information from screening pipeline, including clinical, radiological and blood-derived features, to diagnose lung cancer more accurately and promote precise treatment for IPLs. Meanwhile, we also found that the clinic-RadmC model might also show potential value in confirming the suspected lung cancer cases. It showed superior performance in diagnosing lung cancer among patients aged>50 years, who were identified as high-risk population for lung cancer, as well as among those with radiological pure-solid IPLs, those with large IPLs (size≥15 mm), and pulmonary mass (size≥30 mm), indicating the great potential of it in clinical application.

Nevertheless, several limitations of our study must be acknowledged. First, the role of all involved centers as tertiary care centers might have generated some selection bias within our sample population, since less suspicious SSNs that were likely to be stable over time were rarely referred to us. Future studies investigating the performance of the clinic-RadmC model for less suspicious pulmonary nodules are warranted. Second, the present study was only conducted in clinical centers from China. Owing to differences among disease spectra, further studies recruiting diverse non-Asian populations are needed to demonstrate the generalizability. Third, we only investigated the diagnostic value of plasma cfDNA end motif profile in our study, and future studies are of still necessity to investigate the potential value of different fragment patterns or biomarkers derived from 5mC- or 5hmC-enriched regions for cancer detection. Fourth, although the clinic-RadmC model enables the integration of clinically convenient plasma sources with LDCT data to diagnose lung cancer, its clinical benefit and cost-effectiveness in practical applications still require further study if more people are to be reached.

In conclusion, the clinic-RadmC model constructed in this study attained superior predictive performance to that of single-omics models, and could serve as a more effective, robust and noninvasive tool for optimizing lung cancer diagnoses, providing a promising avenue for improving lung cancer screening programs and reducing the global burden of this prevalent cancer.

Ethical statements

This study was conducted according to the Declaration of Helsinki, and registered in the Chinese Clinical Trial Registry (ChiCTR2100041971). The ethical committees of all involved centers, including Shanghai Pulmonary Hospital, Ningbo Hwa Mei Hospital, The First Affiliated Hospital of Nanchang University, the Affiliated Hospital of Zunyi Medical College, and The First Hospital of Lanzhou University, approved this study (L21-022, L21-022-1, L21-022-2). Written consents were obtained from participants prior to the enrollment, and they were not compensated for participation.

Participants enrollment

As shown in Fig. 1A, participants who had IPLs were prospectively recruited at the Department of Thoracic Surgery of Shanghai Pulmonary Hospital between March 2021 and December 2021 based on the following inclusion criteria: (1) aged ≥ 18 years; (2) did not have a history of malignant tumors within 5 years prior to their diagnoses; (3) did not have definite signs of tumor metastasis; (4) did not receive any anticancer treatments before their diagnoses; and (5) signed an informed consent form. Participants who had multiple resected IPLs were excluded, and those without definite pathological outcomes were also excluded to ensure the robustness of our study. Participants whose samples had inadequate DNA after extraction, whose sequencing quality control tests failed (spike-in DNA > 10 reads with capture efficiency <0.5, and sequencing reads <one million) or whose chest CT images could not be extracted were excluded. In total, 1718 participants were enrolled and randomly classified into training set, validation set, and internal test set at a ratio of 6:2:2. Further, participants from 4 clinical centers (Ningbo Hwa Mei Hospital, The First Affiliated Hospital of Nanchang University, the Affiliated Hospital of Zunyi Medical College, and The First Hospital of Lanzhou University) were recruited from September 2021 to December 2022 (Fig. 1B), and 314 participants were ultimately enrolled, serving as an independent external test set. The disease outcomes of all participants were surgery-, or biopsy-proven, and the baseline characteristics were collected via clinical records.

Plasma sample collection and processing

Approximately 10-ml peripheral venous blood samples were collected from each participant within 2 weeks prior to the initial diagnosis via BD Vacutainer EDTA tubes (Becton, Dickinson and Company, Cat^# 367525), and processed for plasma isolation following the procedures: The blood was rested for 30 min at 4 °C and then centrifuged at 1000×g for 10 min at 4 °C; then, the supernatant was collected and centrifuged at 2000×g for 5 min at 4 °C. The plasma was frozen at -80 °C until DNA isolation. For cfDNA extraction, a VAHTS Serum/Plasma Circulating DNA Kit (Vazyme, Cat^# N902-02) was used, and the cfDNA was quantified with a Qubit fluorometer (Life Technologies). cfDNA (10-30 ng) was applied for library construction: adapter ligation was facilitated via the VAHTS Universal DNA Library Prep Kit for Illumina (Vazyme, Cat^# ND607) according to the protocol provided by the manufacturer, with spike-in controls added (0.01 pg of each amplicon per 10 ng of cfDNA). The spike-in control, including 3 distinct lambda DNA amplification products (~180 bp; one without modification and the other two with 5mC and 5hmC modifications), was prepared as follow⁴⁷: lambda DNA was PCR amplified by Taq DNA Polymerase (NEB) and purified by AMPure XP beads (Beckman Counlter) in nonoverlapping ~180 bp amplicons, with a cocktail of dATP/dGTP/dTTP and one of the following: dCTP, dmCTP or 10% dhmCTP (Zymo)/90% dCTP. The employed primers sequences were as follows: dCTP FW-5′-CGTTTCCGTTCTTCTTCGTC-3′, RV-5′-TACTCGCACCGAAAATGTCA-3′; dmCTP FW-5′-GTGGCGGGTTATGATGAACT-3′, RV-5′-CATAAAATGCGGGGATTCAC-3′; 10% dhmCTP/90% dCTP FW-5′-TGAAAACGAAAGGGGATACG-3′, RV-5′-GTCCAGCTGGGAGTCGATAC-3′.

Methylation, and hydroxymethylation profiling

For cell-free 5mC sequencing, ligated cfDNA (5–10 ng) (VAHTS Universal DNA Library Prep Kit for Illumina V3, Vazyme ND607, and QIAseq UDI Y-Adapter Kit, Qiagen) was used by the cell-free methylated DNA immunoprecipitation and high-throughput sequencing (cfMeDIP-seq)^48,49. The spike-in controls were used to roughly assess the 5mC enrichment ratio instead of the methylated DNA and unmethylated DNA spike-in controls used in the cfMeDIP-seq. We added a Y-shaped long linker sequence containing index information to the cfDNA to protect the integrity of cfDNA methylome profiling process⁵⁰, as described in Supplementary method 1. For cell-free 5hmC sequencing, a low-input whole-genome, cell-free 5hmC sequencing method based on selective chemical labeling (hMe-Seal) was developed⁵¹. The hMe-Seal method is a robust approach that uses β-glucosyltransferase (βGT) to selectively label 5hmC with biotin via an azide-modified glucose for the pull-down of 5hmC-containing DNA fragments for sequencing⁵¹. The standard hMe-Seal procedure requires micrograms of DNA. In our modified approach, cfDNA was first ligated with sequencing adapters, and 5hmC was selectively labeled with a biotin group. After the cfDNA containing 5hmC WAS captured via streptavidin beads, the final library was directly completed by PCR from the beads instead of eluting the captured DNA to minimize the sample loss induced during the purification steps. With our modified approach, the cell-free 5hmC could be readily sequenced from 1-10 ng of cfDNA. The detailed information of this method⁴⁷ has been further described in Supplementary method 2.

cfDNA sequencing data processing and end-motif profiles identification

With respect to the processing strategy applied to the cfDNA sequencing data, the 5mC- and 5hmC-captured libraries were sequenced via 150 bp paired-end runs (300 cycles) on the illumine NovaSeq6000 platforms⁵⁰. Prior to alignment, low-quality sequences were filtered using the Fastp software (version 0.21.0). BWA (version 0.7.17-r1188) was used to process sequencing reads, which were aligned to hg19 and spike-in DNA with default parameters. SAMtools (version 1.7) was used to filter the generated SAM files with parameter settings of ‘-f 3 -F 3852 -q 30’ to include high quality, primary alignment, properly paired reads, followed by converting to BAM format. Picard (version 2.25.5) was employed to sort and index filtered SAM files and to ensure the removal of duplicate reads before subsequent analyses. Bedtools (version 2.26.0) bamtobed was utilized to convert BAM alignment to BED with ‘-bedpe’ parameter. Only paired-end reads with both ends aligned to the same autosome with the correct orientation, spanning an insert size from 20 bp to 1000 bp, were used for downstream analyses. Three types of spike-in DNA sequences were included into reference; then, capture efficiency, as a quality control measurement for the 5mC and 5hmC, was calculated as counts of reads aligned to a type-specific spike-in DNA divided by counts of reads aligned to total spike-in DNA. The 4-mer (bp) end motifs profiles of the 5mC- and 5hmC-sequencing data were extracted after generating the fragment file. Then, the 4 bp end-motifs was extended to 6 bp end-motifs, and the first 6-nucleotide sequence on each 5’ end of plasma DNA fragment was identified (Fig. 2A), since the 6 bp end-motif might exhibit better diagnostic performance for cancer detection^9,29. The frequency of end motifs (4 bp: 4^4 = 256; 6 bp: 4^6 = 4096) out of the total number of motifs was calculated.

Chest CT image acquisition and nodule annotation

The thin-section chest CT images that were closest to the initial diagnosis date were extracted from the picture archiving and communication (PACS) system and imported into the 3D-slicer software (version 4.10.1), with the CT scanning parameters detailed in Supplementary method 3. Each three-dimensional region of interest (ROI) was independently annotated with a bounding box by two junior thoracic surgeons (M.M. Zhao and J.J. Deng) and a junior radiologist (T.T. Wang) under the lung (level, -450 HU; width, 1500 HU) and mediastinal (level, 40 HU; width, 400 HU) window settings. The imaging subtypes, emphysema and spiculation were recorded, and IPLs were further allocated as pure ground-glass opacity (GGO), part-solid nodules or pure-solid nodules. The pure-GGOs comprised IPLs only with GGOs, and the part-solid nodules consisted of GGOs and solid components, whereas the pure-solid nodules had only solid components. The radiological solid-component size was measured, as well. The disagreement was solved via a group discussion with a senior radiologist (X.W. Sun).

Fragmentomic analyses

The principal component analysis (PCA) was performed to evaluated the potential batch effects in 5mC- and 5hmC-sequencing data between different hospital regions. The dimensions reduction of end motif profiles were conducted on the training set as follow: (1) a Wilcoxon rank-sum test was used to select differently expressed features between benign and malignant samples, and only the feature with p < 0.01 were selected; (2) a univariate logistic analysis was performed and only the features with p < 0.01 were selected for the next step; (3) the minimum redundancy and maximum relevancy method was applied to rank each selected feature depending on its relevance to malignancy, and the top 100 most significant features were selected; (4) the 100 identified features were subjected to the least absolute shrinkage and selection operator (LASSO) model. The area under the receiver operating characteristic curve (AUC) was used as the criterion, and the optimal feature set with the maximum AUC was selected as input features. Hierarchical clustering analysis was applied to evaluate the classification performance of these identified end motifs. Then, a support vector machine (SVM) was used to construct the predictive models using the end motifs selected from 5mC-sequencing data (5mC model), and 5hmC-sequencing data (5hmC model) after multiplying the coefficient identified by LASSO, whereas hyperparameter optimization was performed via a grid-search and manual fine-tuning with 5-fold cross-validation. The models based on 4 bp and 6 bp end motif features were simultaneously established, and those with the best performance were reserved for subsequent analyses. We further evaluated the diagnostic performance of the end motif profiles obtained after gradually downsampling the 5mC- and 5hmC-sequencing fragments to 25%, 50%, and 75%. In addition, the chromatin accessibility contributes to cfDNA preferred ends⁴³, and thus, the changes in chromatin accessibility and DNA nuclease activity or the combined effect of both could have impacts on the cfDNA end motifs. Increased DNA methylation is often associated with reduced chromatin accessibility, influencing the TF binding process. Thus, we hypothesize that if the sequence pattern of the end motif resembles that of the TF binding motif, the end motif may serve as a potential binding sequence for TF. Also, epigenetic changes in the vicinity of end motif sequences may affect the TF binding⁵², leading to the transcription disturbance of TF regulatory genes. Accordingly, the functional exploration of end motifs was performed via 3 steps: (1) matching with the known motifs contained in the Homo sapiens Comprehensive Model Collection database to identify TFs (https://meme-suite.org/meme/tools/tomtom); (2) searching for genes regulated by the identified TFs in the TRRUST V2 database (https://www.grnpedia.org/trrust/); and (3) applying gene ontology (GO) for functional annotation and pathway enrichment analyses.

Radiomic analyses

A feature pyramid-DenseNet backbone network was applied to multiscale features from CT images, and the deep learning-based radiomic (DL-radiomics) model was then obtained via a fully connected network after fusing the CT features. The details of the image preprocessing and model construction are described in Supplementary method 4-6. The classic radiomic (C-radiomics) model was also established to evaluate the potential superiority of the deep learning technique, with the details of the procedures of constructing the C-radiomics model presented in Supplementary method 7.

To enhance the interpretability of the DL-radiomics model, the attribution scores of radiomics features were calculated by the deep learning importance features (DeepLIFT) method, which was a backward propagation method for gaining interpretability of deep neural networks^53,54, to explain the contribution of input feature to the model’s output. We further explored the biological basis of the DL-radiomics model by applying the Radiogenomics dataset obtained from the Gene Expression Omnibus database. Gene set enrichment analysis (GSEA) and single-sample gene set enrichment analysis (ssGSEA) were performed to reveal the heterogeneity of the genetic pathways and the infiltration of immune cells in the tumor microenvironment among patients with different DL-radiomics scores. The processes of the biological exploration are detailed in Supplementary method 8.

Multiomics analyses

A univariate logistic regression was used to evaluate the associations between the clinical factors, 5mC model score, 5hmC model score and DL-radiomics model score with the malignancy risk of IPLs on the training set, and variables with p < 0.05 were fed into the multivariate logistic regression analyses via backward stepwise selection with the Akaike information criterion. Ultimately, all variables with p < 0.05 were selected as candidates, and the correlations among them were evaluated via the Pearson correlation coefficients analyses. Then, a series of 5 predictive models involving different feature sets were derived through (1) clinical factors only (clinical model), and adding (2) the DL-radiomics model score (clinic-Radiomics model), (3) the 5mC model score (clinic-mC model), (4) the DL-radiomics model score and 5mC model score (clinic-RadmC model), and (5) the DL-radiomics model score, 5mC model score and 5hmC model score (clinic-Rad(h)mC model) to evaluate the superiority of multiomics assays to models including peer assays. Multivariable logistic regression was used to establish the multiomics models, with a fivefold cross-validation scheme used to optimize the model hyperparameters.

Support for guiding the individual management process

The cutoff values of the models were identified based on the maximum Youden index on the training set, and the participants were classified into low-risk (the model score<cutoff values) and high-risk (the model score≥cutoff values) subgroups. Then, a hypothetical clinical management scheme was determined based on model classification results by assuming that any patients “ruled out” (identified as low-risk) would not undergo additional testing beyond follow-up chest CT scans, whereas any patients “ruled in” (identified as high-risk) would proceed directly to receive invasive procedures. The stratification efficiency was analyzed with respect to the reduction in the rate of invasive procedures unnecessary for benign diseases and potential delayed treatment for lung cancers in comparison with actual care.

Statistics and reproducibility

No statistical method was used to predetermine the sample size of this study since it was a prospective, observational study. Besides, participants were selected according to the inclusion and exclusion criterions described in Fig. 1. This study was designed as a prospective, observational, and case-control study, and all participants were classified based on their histopathological results. Thus, the experiments were not randomized, and the investigators were not blinded to group allocation during outcome assessment. However, the sample operation team was blinded to the group status during the whole sequencing process, and all CT reviewers were blinded to the clinical information during the annotation process. No sex or gender-related analysis was performed, as well.

Continuous and categorical variables of enrolled participants were reported as means with standard deviations (SDs) and frequencies with percentages, respectively. The Pearson χ² test or Fisher’s exact test, and Student’s t test or Wilcoxon signed-rank test were used to compare the participants characteristics, when applicable. Univariate and multivariate logistic regression analyses were performed to assess the associations between different participants characteristics and the odds of malignant diseases. The AUC was used to evaluate the diagnostic performance of the models, and compared by the DeLong’s test, with confidence intervals (CIs) generated based on 10,000 bootstrap replicates. Performance metrics, including sensitivity, specificity, positive predictive value (PPV), negative predictive value (NPV), and accuracy, were calculated based the cutoff values derived from the training set. A calibration curve analysis and a decision curve analysis (DCA) were applied to determine the clinical usefulness, with net reclassification improvement (NRI) and integrated discrimination improvement (IDI) analyses applied to evaluate the clinical benefits of different approaches. The predictive performance of the models was further assessed in subgroups stratified by age (≤50 years, >50 years), radiological image type (pure-GGO, part-solid, pure-solid), and radiological solid-component size (subcentimeter nodules with size≤ 10 mm, large nodules with 15 mm≤size≤30 mm, and mass with size>30 mm). To facilitate the clinical use of the models, the Shapley additive explanations (SHAP) method was used to examine the importance levels of each predictor and provide consistent attribution values for features incorporated into the models⁵⁵. The statistical analyses were two-sided, and performed in SPSS (version 23.0, IBM SPSS Statistics), R software (version 4.1.3, http://www.Rproject.org) or Python (version 3.10, https://www.python.org/), with p < 0.05 indicating statistical significance.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.