In this study, a total of 1,168 participants were enrolled from Xijing Hospital (March 2019 to December 2019) and Tangdu Hospital (March 2019 to February 2020) of Fourth Military Medical University (FMMU) and Changhai Hospital of Second Military Medical University (March 2019 to September 2019), including patients with HBV-related HCC (n=541), non-HBV-related HCC (n=30), patients with CHB (n=125), patients with LC (n=176) and healthy controls (n=296). The demographic and clinicopathological characteristics were summarized in Supplementary Table 1. Patients with HCC, CHB or LC were diagnosed based on the American Association for the Study of Liver Diseases (AASLD) guidelines [22]. HC was collected from individuals who underwent routine physical examinations. The inclusion criteria and exclusion criteria of patients were described in the Supplementary information. Tumor stage was determined according to the Barcelona Clinic Liver Cancer (BCLC) staging system. The HCC patients from BCLC stage 0-C were selected and enrolled in this study. All HCC patients included in this study were subjected to surgical removal for confirmation of pathological diagnosis and then followed up every 6 months. The vascular invasion, lymph node metastasis, and bile duct invasion were confirmed by histopathology. Finally, follow-up data was successfully collected from 388 HCC patients, with the final follow-up time in December 2023. The median follow-up duration was 29.9 months, ranging from 3 to 57 months, and 263 patients were dead of HCC during the follow-up. Overall survival (OS) was considered as the time span from diagnosis to death or last follow-up. Recurrence-free survival (RFS) was considered as the time from surgery to disease recurrence, death, or the last follow-up. Diseasefree survival (DFS) was defined as the time from surgery to the occurrence of disease recurrence, the development of a second primary cancer, death or the last follow-up. This study was approved by the three Hospital Review Boards and all participants provided written informed consent.

Sample collection was performed following the acquisition of diagnostic information. Peripheral blood samples (5 mL for each) from all participants were collected 1 day before treatment and plasma was separated within 2 hours using a standard two-step centrifugation protocol and then stored at –80°C.

The QIAamp Circulating Nucleic Acid kit (Qiagen, USA) was utilized for DNA extraction from the plasma according to the manufacturer’s instructions [19]. The DNA was quantified using Qubit 4.0 (Thermo Fisher, USA). Subsequently, the whole-genome sequencing (WGS) libraries were constructed using plasma cfDNA (10–20 ng for each sample) with the NEB ultra v2 kit (NEB, US) as previously described [23]. The number of PCR cycles during library construction was 12. Then, the WGS libraries were hybridized with homemade biotinylated mtDNA capture probes as previously described [17]. Finally, all captured mtDNA libraries were sequenced by the Illumina Hiseq XTen platform using paired-end runs with 2×150 cycles (PE 150). To minimize bias resulting from systematic errors in the experiment, we ensured that samples of different disease types are included in the same batch of experiments, and standardized experimental procedures were followed. And, all experiments were carried out by two individuals in a centralized laboratory, utilizing the same instruments across all samples.

After sequencing, Illumina adapters and bases with quality scores below 30 were trimmed in raw sequencing data using Fastp software (v.0.20.0) [24]. Furthermore, to minimize contaminations from nuclear sequences of mitochondrial origin (NUMTs), qualified reads were aligned to both the revised Cambridge reference sequence (rCRS) and the human genome reference (hg19) using BWA software (v.0.7.1521). Subsequently, reads that only mapped to rCRS were retained for further analyses. Picard Tools software (v.1.119) was applied for sorting the reads and the duplicated reads were removed. The Genome Analysis Toolkit 4 software (v.3.2-2) was employed for local realignment. The final mtDNA coverage depths ranged from 2,469 X to 3,383 X for capture-based mtDNA sequencing data (Supplementary Table 2).

Based on our established method, the average sequencing depth of mtDNA and the average sequencing depth of the nuclear reference gene were used to calculate the relative mtDNA CN. And the following formula was used to calculate the mtDNA copy: 2×average depth of mtDNA/average depth of nuclear reference gene [21].

Two categories of ccf-mtDNA fragmentomic features including fragment-related and 5’ end-related features were extracted from the mtDNA capture based sequencing (mtDNA Cap-seq) data and detailed methods were described in the Supplementary information.

A total of 210 patients with HCC and 240 non-HCC individuals (120 CHB/LC and 120 HC) from Xijing Hospital were randomly divided into a training cohort and an internal validation cohort at a ratio of 7:3. A total of 210 patients with HCC and 240 non-HCC individuals (120 CHB/LC and 120 HC) from Tangdu Hospital and 121 patients with HCC and 117 non-HCC individuals (61 CHB/LC, and 56 HC) from Changhai Hospital were used as external validation cohorts 1 and 2, respectively.

The training cohort was used to construct the HCC detection model utilizing fragment-related, 5’ end-related and all fragmentomic features, respectively, using the random forest (RF) algorithm from R package ‘randomForest’ (v.4.6-14) with the number of trees set to 500 and six features randomly selected for each tree. The ten-fold cross validation was used to evaluate the model performance in the training cohort. Consequently, the HD score for each patient was determined from the HCC detection model. A detection model integrating the HCC detection (HD) score and AFP was formulated through logistic regression in the R package ‘stats’ (v.4.1.2). Internal and external validation cohorts were used for evaluating the performance of HCC detection model based on all fragmentomic features. The performance of our model was also evaluated in a non-HBV-related HCC cohort.

The LASSO-Cox regression analysis, using CoxPHSurvivalAnalysis model in ‘scikit-survival’ package (v.0.22.1), was utilized to develop the HCC prognosis prediction model based on all ccf-mtDNA fragmentomic features in the training HCC cohort, which consisted of 197 patients with HCC from Xijing Hospital. A five-fold cross validation was conducted to identify the minimum lambda value, resulting in 20 features with non-zero coefficients. Subsequently, the HPP score (range: 0 to 1) was calculated as the sum of values weighted by the coefficients from the Cox hazard model. Based on the median HCC prognosis prediction (HPP) score of HCC patients in the training cohort, all HCC patients were divided into HPP score ≤0.65 and HPP score >0.65 groups. Next, a combined HCC prognosis prediction model integrating the HPP score and clinical characteristics (BCLC stage, Child-Pugh class and AFP level) was formulated through multivariate Cox regression. The predictive ability of the HCC prognosis prediction model was evaluated in the validation HCC cohort, which consisted of 191 patients with HCC from Tangdu Hospital.

Totally, based on our sequencing and analysis methods, it will take 5–7 working days to accomplish one processing pipeline (including cfDNA extraction, library construction, next generation sequencing and data analysis) and thus obtain the fragmentomic data of ccf-mtDNA, with cost of 30–40 US dollars for each sample.

Statistical analyses were performed using GraphPad Prism v.9.5.1 (GraphPad Software). Detailed methods were described in the Supplementary information. All P values reported were two-tailed and a significance level of 0.05 was used.

A total of 1168 participants including 541 patients with HBV-related HCC, 30 patients with non-HBV-related HCC, 125 patients with CHB, 176 patients with HBV-related LC and 296 HC were recruited from three hospitals into the training and validation cohorts (Fig. 1). As shown in Supplementary Figure 1, our data indicated that patients with HCC were significantly older than HC and patients with CHB/LC (both P<0.05). Further analysis revealed a higher male percentage in the HCC group compared to the CHB/LC and HC groups. In addition, no significant difference in the demographic and clinicopathological characteristics was observed among the Xijing, Tangdu and Changhai cohorts (all P>0.05, Supplementary Table 1).

A comprehensive analysis of ccf-mtDNA fragmentomic features was conducted based on capture-based NGS data with an average sequencing depth of 2,882 X (Supplementary Table 2) among 296 HC, 125 CHB, 176 LC and 541 HCC patients from both the training and all validation cohorts. Our results indicated no difference of ccf-mtDNA fragmentomic features between the patients with CHB and LC and thus these two groups were combined for subsequent analyses (Supplementary Fig. 2). Among four groups consisting of HC, patients with CHB/LC, HCC patients at the early (EHCC) and late BCLC stage (LHCC), we observed a gradual decrease of fragment size from HC to LHCC (Fig. 2A, B). Furthermore, to investigate the fragmentation profile of ccf-mtDNA at single-base resolution, the fragment size distribution (FSD) score was determined as the overall coverage ratio of long to short ccf-mtDNA fragments covering the given site, with the median length as a cutoff. Our data revealed that both EHCC and LHCC groups showed highly variable fragmentation profile when compared with other two groups (Fig. 2C). A gradual decrease of correlation to median HC profile was observed from HC to LHCC group (Fig. 2D). Subsequently, the peak number of the fragmentation profile was calculated for all groups and new peaks were defined based on the comparison of peak location with the median HC profile. Our data showed a gradual increase in total peak number and new peak number from HC to LHCC group (Fig. 2E, F). In addition, as shown in Figure 2G, H, our results revealed a significant difference in the 5’ end base preference among four groups. Specifically, we observed an overrepresentation of C-end and an underrepresentation of G-end in the LHCC group compared to other three groups. In comparison, the T-end base preference exhibited an increased trend from CHB/LC to HC and EHCC to LHCC, whereas A-end base preference exhibited a decreased trend from CHB/LC to LHCC. Furthermore, we calculated the 5’ end motif diversity score (MDS) based on the proportion of 4-mer end base motif and observed an increased trend from CHB/LC and EHCC to LHCC to HC (Fig. 2I). When ccf-mtDNA fragmentomic features among the four groups were stratified by cohort, very similar results were obtained (Supplementary Fig. 3). Also, some of the 5’ end features such as 5’ C- and G- end base preference and 5’ MDS did not differ between the CHB/LC and EHCC groups. We also observed no obvious association between age or gender and the ccf-mtDNA fragmentomic features in all groups (Supplementary Figs. 4, 5). Altogether, these findings indicate aberrant ccf-mtDNA fragmentomic features in patients with HCC, strongly supporting the clinical application as promising detection biomarkers.

Based on fragment-related and 5’ end-related ccf-mtDNA fragmentomic features, the widely utilized RF algorithm was used to construct the HD model in the training dataset from Xijing cohort (Supplementary Fig. 6). A HD score was generated by ten-fold cross-validation for each individual to assess the performance of the HD model in distinguishing HCC from non-HCC individuals (including at-risk subjects with CHB/LC and healthy controls) and at-risk subjects with CHB/LC. As shown in Supplementary Figure 7, the HD model based on all fragmentomic features exhibited a superior area under the curve (AUC) in distinguishing HCC from non-HCC (0.9965, 95% CI 0.9937–0.9994) and CHB/LC (0.9943, 95% CI 0.9894–0.9992) when compared with models based on fragment-related or 5’ end-related features.

Furthermore, we tested the all fragmentomic features-based HD model in the internal validation dataset (n=135) from Xijing cohort. Our data revealed that both EHCC and LHCC groups exhibited significantly higher HD score than HC and CHB/LC groups (Fig. 3A). As shown in Figure 3B and Supplementary Table 3, the HD score indicated the excellent capacity for distinguishing patients with early HCC from non-HCC subjects (AUC=0.9910, 95% CI 0.9814–1.000); sensitivity=92.06%; specificity=94.44%; score cut-off=0.4968). Considering that CHB and LC in China are at a significantly higher risk of developing HCC and often complicate the diagnosis of HCC, we evaluated the performance of our HD model in distinguishing EHCC from at-risk individuals with CHB/LC. Our results showed a high AUC of 0.9903 (95% CI 0.9778–1.000), with sensitivity of 92.06% and specificity of 94.44%. As expected, our HD model exhibited very good performance for distinguishing EHCC from healthy controls (AUC=0.9916). Moreover, late HCC (BCLC stage B/C) patients were perfectly distinguished by our HD model from non-HCC, CHB/LC and healthy controls with AUC of 0.9934, 0.9921, and 0.9947, respectively.

To verify the applicability of our model in clinical practice, the performance of fixed HD model was then evaluated in two external validation cohorts. Very similar in internal validation cohort, a notable increase in HD score was observed from HC and CHB/LC to LHCC in both external validation cohorts (Fig. 3D–I). Furthermore, the results demonstrated remarkable discriminatory capacity of our HD model in distinguishing patients with EHCC from non-HCC in both external val idat ion cohor t 1 and 2 (AUC=0.9861 and 0.9826), as well as distinguishing patients with EHCC from CHB/LC (AUC=0.9840 and 0.9818) and healthy controls (AUC=0.9882 and 0.9835). Very similar performance was obtained in distinguishing patients with LHCC from non-HCC, CHB/LC and HC in both external validation cohorts 1 and 2. In addition, by applying HD model in a set of non-HBV-related HCC (Supplementary Table 4), including 56 HCV-related HCC and 64 HCC patients without hepatitis, we observed a pretty good performance in distinguishing them from CHB/LC and HC (all AUC>0.9385), although the power was to some extent decreased when compared with those in HBV-related HCC due to the absence of non-HBV-related HCC in the training cohort (Supplementary Fig. 8). These results clearly indicate the high performance and generalization of ccf-mtDNA fragmentomics-based HD model in HCC early detection.

To evaluate the real-world influence of ccf-mtDNA fragmentomics-based strategy in the context of HCC detection, we compared the performance of our HD model with the current screening measurement of AFP level. As shown in Figure 4A, B, in the most challenging clinical scenario of distinguishing EHCC from CHB/LC, our HD model exhibited superior performance to AFP-based method in both external validation cohort 1 (AUC=0.9840 vs. 0.8238) and 2 (AUC=0.9818 vs. 0.8407). Accordingly, HD model had the sensitivity of 89.93% and 89.61% at a score cut-off of 0.4968 in both cohorts, while AFP (cut-off=20 ng/mL) only achieved the sensitivity of 56.83% and 62.34%, which is consistent with previous reports [12]. Very similar results were observed in distinguishing LHCC from CHB/LC (Fig. 4C, D). Additionally, the combination of HD score and AFP only provided a very slight improvement in both EHCC and LHCC detection over the HD model alone in both external validation cohorts. Importantly, our HD model still demonstrated remarkable capacity in distinguishing 79 and 39 HCC patients with AFP level below 20 ng/mL (AFP-negative HCC) who would have been misclassified based on AFP alone from CHB/LC, with AUC of 0.9863 and 0.9859 in external validation cohort 1 and 2, respectively (Fig. 4E, F and Supplementary Table 5). As expected, consistent performance of HD model was achieved in distinguishing AFP-positive HCC from CHB/LC in both cohorts (Fig. 4G, H).

Furthermore, considering the previous reports that mtDNA CN in plasma cfDNA may be used as biomarker for HCC detection, we evaluated the performance of mtDNA CN in differentiating HCC from CHB/LC and HC. Similar to previous reports, both EHCC and LHCC groups exhibited a significant increase in mtDNA CN compared to HC and CHB/LC groups (Supplementary Fig. 9). When mtDNA CN was used as a classifier, an AUC of 0.5809 and 0.6278 was achieved in distinguishing EHCC patients from HC and patients with CHB/LC, respectively, while a slightly increased AUC of 0.6642 and 0.7061 was achieved in distinguishing LHCC patients from HC and patients with CHB/LC. Furthermore, no significant difference in distinguishing HCC from non-HCC was found between the combination of the HD score & mtDNA CN and the HD score alone (Fig. 4I–L). Taken together, our results highlight the utility of ccf-mtDNA fragmentomic features as much more effective alternative biomarkers for HCC detection.

As tumor biomarkers may be affected by clinical characteristics, we comprehensively investigated the association of the HD score with demographic variables such as age and gender or measures of liver dysfunction in HC, CHB/LC and HCC groups gathered from two validation cohorts. Our data revealed no obvious association between the HD score and age or gender in all three groups, as well as between the HD score and measures of liver dysfunction in both CHB/LC and HCC groups (all P>0.05, Supplementary Fig. 10). Moreover, we further examined the relationship between HD score and variables of HCC progression. Notably, as shown in Figure 5A, HD score significantly increased as BCLC stage advanced in HCC patients (P for trend <0.01). In particular, HD score just showed a very slight decrease of detection power for HCC at stage 0 from CHB/LC (AUC=0.9343, sensitivity=77.78%) when compared with those at other three stages (AUC=0.9852, 0.9869 and 1.000, sensitivity=90.40%, 90.67%, and 100.0%, respectively, Fig. 5B, C). In addition, a significantly higher HD score was also observed in HCC patients with Child-Pugh class B or increasing tumor size (>2 cm) or more tumor sites (>1) when compared with corresponding control groups (all P<0.05, Fig. 5A). And patients with HCC who had the features of vascular invasion, lymph node metastasis, or bile duct invasion exhibited significantly higher HD scores compared to those without these features (all P<0.001, Supplementary Fig. 11A). In these subgroups, HD score exhibited almost equivalent detection performance (Fig. 5B, C and Supplementary Fig. 11B, C). These observations suggest that ccf-mtDNA fragmentomics-based HD model may not be remarkably affected by clinical characteristics, indicating better clinical applicability.

To evaluate the applicability of ccf-mtDNA fragmentomic features in the prognosis prediction of HCC patients undergoing surgical treatment, we assigned 197 HCC patients with complete survival information from Xijing Hospital into a training HCC cohort and 191 HCC patients from Tangdu Hospital into a validation HCC cohort. LASSO-Cox regression was employed to diminish dimensionality and construct an HCC prognosis prediction model which included 20 features (Supplementary Fig. 12). The patients with HCC were first stratified into high-risk and low-risk groups by the median HPP score of 0.65 generated from the HCC prognostic prediction model constructed based on training HCC cohort. Subsequently, Kaplan–Meier curves revealed a significantly shorter median OS in the high-risk group compared to the low-risk group (log-rank test, P<0.0001; HR 0.30, 95% CI 0.21–0.44, Fig. 6A). The AUC of HPP score for predicting 1-year, 2-year, and 3-year survival were 0.8838, 0.8558, and 0.8347, respectively, in the training HCC cohort (Fig. 6B).

As shown in Supplementary Table 6, the univariable and multivariable Cox regression analyses indicated that the HPP score was significantly associated with OS and may serve as an independent prognostic factor. Furthermore, the prognosis prediction performance of the HPP score was evaluated in validation HCC cohort. As shown in Figure 6C, the median OS time in the high-risk group was significantly shorter than that in the low-risk group (log-rank test, P<0.0001; HR 0.37, 95% CI 0.26–0.53). And the AUC of HPP score for predicting 1-year, 2-year and 3-year survival were 0.8333, 0.8145 and 0.7958, respectively (Fig. 6D). Stratified analysis indicated that significant association between HPP score and overall survival existed in all subgroups (all P<0.001, Fig. 6E), further supporting the independent role of ccf-mtDNA fragmentomic features in HCC prognosis prediction. The predictive performance of the HPP score for RFS and DFS was further evaluated. Kaplan–Meier curves demonstrated a significantly shorter median RFS and DFS in the high-risk group compared to the low-risk group in both the training (HR 0.32 and 0.31, respectively) and validation (HR 0.39 and 0.37, respectively) HCC cohorts (log-rank test, P<0.0001, Supplementary Fig. 13). Notably, the combination of HPP score and clinical characteristics significantly improved the ability to predict HCC prognosis in both the training (HR 0.12, Supplementary Fig. 14A, B) and validation HCC cohorts (HR: 0.14, Supplementary Fig. 14C, D). Altogether, these results demonstrate that ccf-mtDNA fragmentomic features can contribute to death risk stratification and prognosis prediction in patients with HCC.