Persistent Homology Identifies Pathways Associated with Hepatocellular Carcinoma from Peripheral Blood Samples

Abstract

Topological data analysis (TDA) methods have recently emerged as powerful tools for uncovering intricate patterns and relationships in complex biological data, demonstrating their effectiveness in identifying key genes in breast, lung, and blood cancer. In this study, we applied a TDA technique, specifically persistent homology (PH), to identify key pathways for early detection of hepatocellular carcinoma (HCC). Recognizing the limitations of current strategies for this purpose, we meticulously used PH to analyze RNA sequencing (RNA-seq) data from peripheral blood of both HCC patients and normal controls. This approach enabled us to gain nuanced insights by detecting significant differences between control and disease sample classes. By leveraging topological descriptors crucial for capturing subtle changes between these classes, our study identified 23 noteworthy pathways, including the apelin signaling pathway, the IL-17 signaling pathway, and the p53 signaling pathway. Subsequently, we performed a comparative analysis with a classical enrichment-based pathway analysis method which revealed both shared and unique findings. Notably, while the IL-17 signaling pathway was identified by both methods, the HCC-related apelin signaling and p53 signaling pathways emerged exclusively through our topological approach. In summary, our study underscores the potential of PH to complement traditional pathway analysis approaches, potentially providing additional knowledge for the development of innovative early detection strategies of HCC from blood samples.

1. Introduction

Primary liver cancer is one of the most common forms of cancer worldwide, with more than 905,000 new cases diagnosed and over 830,000 cancer-related deaths reported each year [1]. Hepatocellular carcinoma (HCC) is the most prevalent type of primary liver cancer accounting for over 75% of all cases, thus representing a significant global health concern [1]. One of the major challenges in HCC diagnosis is the lack of noticeable symptoms in the early stages of the disease, often leading to patients being diagnosed only when the disease has reached an advanced stage [2].

Diagnosing HCC and predicting treatment responses can be achieved using various biomarkers [3]. The World Health Organization (WHO) defines a biomarker as “any substance, structure or process that can be measured in the body or its products and [that can] influence or predict the incidence of outcome or disease” [4]. The conventional biomarker for early HCC detection, α-fetoprotein (AFP), is no longer recommended due to its limited sensitivity and specificity [5,6]. While newer biomarkers like des-γ-carboxy-prothrombin (DCP), osteopontin, interleukin-6 (Il-6), and Golgi protein 73 (GP73) have shown promise, none have achieved sufficient sensitivity or specificity [7]. Combining multiple biomarkers is a current approach showing potential to enhance HCC detection and monitoring [3]. Nevertheless, there is still room and a need for more accurate and reliable biomarkers for effective HCC prognosis, diagnosis, and treatment. This is especially true at the level of functional pathways involving many genes, which still poses a challenge in current research.

Topological data analysis (TDA) is a modern interdisciplinary field that combines algebraic topology, computational geometry, and statistical learning to analyze large volumes of high-dimensional data and extract relevant information [8]. It offers advantages such as independence from the choice of metric, robustness against noise, and a qualitative approach to understanding the ‘shape’ of data [9,10]. It can also integrate and visualize multidimensional data, potentially providing a better understanding of relationships between clinical samples where traditional data analysis methods may fall short [11].

TDA has proven to be a versatile tool across diverse applications, including digital data classification [12], multidimensional data derived from magnetic resonance imaging (MRI) scans of osteoarthritis patients [13], and the integration of computed tomography (CT) and MRI scans with clinical data from patients with traumatic brain injury [14]. Its efficacy extends to the analysis of mammogram scans/images [15,16] and it has demonstrated utility in genomics, particularly in the analysis of DNA copy number aberrations [17]. In the realm of infectious diseases, TDA has been applied to whole-blood RNA sequencing (RNA-seq) gene expression data [18] and to measurements of inflammatory biomarkers in blood plasma [19] from coronavirus disease 2019 (COVID-19) patients. A TDA-based approach also demonstrated its capability to distinguish idiopathic pulmonary fibrosis from multiple other diseases based on peripheral blood mononuclear cell (PBMC) gene expression data [20].

In short, TDA holds immense potential as a reliable tool for analyzing multidimensional data, such as gene expression levels obtained by RNA-seq. However, a specific implementation applying this method to functional pathways in order to capture topological changes based on the set of genes involved in a functional pathway is still lacking. To address this, our research utilizes persistent homology (PH) to identify significant pathways as potential prognostic biomarkers.

Figure 1 visually outlines our research pipeline, illustrating the process of evaluating alterations in genome-wide correlation of gene expression levels between control and disease sample classes. With this approach, we first identify a set of topological descriptors that show significant changes between HCC samples and normal controls. Then, we apply this pipeline to multiple pathway-specific subsets of genes to assess changes in the correlation of gene expression specifically within these pathways, mirroring the evaluation conducted on the full (genome-wide) gene expression data. Thus, we identify those cellular pathways that reflect the global changes in HCC samples as potential biomarker pathways. In addition to the pipeline sketch in Figure 1, we have included Algorithm 1 in Section 3 to further enhance the presentation of our proposed method.

In the case of analyzed HCC expression data from peripheral blood, taken from a study by Han et al. [21], we hypothesized that PH could facilitate early HCC detection by identifying such potential biomarker pathways. Specifically, we analyzed RNA-seq gene expression data from PBMCs of both HCC patients and normal controls, with 17 samples in each class. For each class, we first mapped the gene expression profile of each sample to a multidimensional space to form a point cloud. Then, we utilized PH to explore various topological descriptors associated with a persistence diagram in order to extract relevant features from the shape of the point cloud. Finally, we compared the descriptors obtained from the two sample classes to identify relevant differences.

Our analysis unveiled several significant cellular pathways whose expression patterns differ between peripheral blood from HCC patients and that from healthy controls. We evaluated our findings by conducting a comparative assessment with the frequently used enrichment-based pathway analysis, revealing pathways that were common or unique to both methods. The literature support and discussion of our results clearly demonstrate the benefit of PH for detecting functional pathways. Our findings reveal that previously unidentified pathways can be uncovered through the application of TDA using PH, and a set of key topological descriptors that effectively capture changes in gene expression between healthy and disease samples can be defined. This approach proves to be a promising tool that can be applied to other datasets and diseases, hopefully contributing to a deeper understanding of the complexities inherent in biological data and systems.

With this study, and the application of TDA for the identification of disease-associated pathways, we intend to address a major limitation of the most commonly used pathway analysis strategies: the assumption that alterations in the biological activity of a pathway depend on the overexpression or underexpression of a significant number of the genes involved in the pathway. However, for both signaling pathways and metabolic pathways, the output can depend on other factors such as coordinated feedback and timing [22,23,24].

We hypothesized that such effects could be better captured by TDA rather than an enrichment analysis approach based on differential expression of individual genes, and thus TDA might offer additional clues for pathway analysis, complementing traditional approaches rather than replacing them.

This motivated our exploration of TDA, specifically PH, to identify key pathways associated with HCC from peripheral blood samples. Our aim is to demonstrate how PH can provide valuable insights into intricate patterns and relationships in complex biological data, particularly in the context of HCC. By analyzing RNA-seq data, we seek to bridge gaps in current pathway analysis approaches.

The results of this pilot study showcase that PH-based methods can complement traditional pathway analysis approaches and provide a promising tool that could be applicable to various datasets and diseases, thus contributing to a deeper understanding of complexities inherent in biological data and systems. Our key contributions include the application of PH to analyze RNA-seq data from peripheral blood samples, the identification of topological descriptors capable of capturing changes in gene expression between healthy and disease samples, the identification of significant cellular pathways distinguishing HCC and healthy samples which include pathways known to be associated with HCC, such as the apelin signaling pathway, the IL-17 signaling pathway, and the p53 signaling pathway, as well as previously unidentified pathways. In addition, comparative analysis with the frequently used enrichment-based pathway analysis method reveals pathways that are common or unique to both methods, which underscores that they can be considered as complementing each other.

The paper is structured as follows: this section described the biomedical context and the objectives of the study. In Section 2, we present a basic introduction to topological data analysis. Section 3 details the materials and methods employed in the study. In Section 4, our results are presented, followed by a discussion in Section 5. Finally, Section 6 encompasses the conclusions and future research directions.

2. Topological Data Analysis: Background

In this section, we will provide a brief overview of homology (in the mathematical sense of topology, not biology or evolution) and persistent homology. For detailed formal definitions and further mathematical terms, please refer to Appendix A.

2.1. Simplicial Complex and Homology

Topology is a mathematical discipline that explores the structural relationships and connections between various regions within a space, similar to observing individual points or composed objects and determining how they are connected and what makes them different from one another. Having in mind a rubber sheet that can be stretched and deformed, but not torn or punctured, topology can be used for understanding which shapes can be transformed into one another through stretching without breaking [25]. TDA builds on these concepts to decipher the essential features of a dataset, even if it is composed of a limited set of data points. It thus helps to determine key characteristics of a complex shape using just a few points.

Homology is a mathematical technique that utilizes algebra to study the complex shapes of objects, starting from simple building blocks called $k$-simplices (like points, lines, triangles, and more), as illustrated in Figure 2. These blocks are assembled to create more complex shapes, like networks.

Usually, specific nomenclature is used to refer to the $k$-simplices of the lowest dimensions, $k$: “vertex” (point) for a 0-simplex, “edge” (two points joined by a line) for a 1-simplex, “triangle” (three points with joining edges and their enclosed space) for a 2-simplex, and “tetrahedron” (a solid triangular pyramid) for a 3-simplex (see Figure 2).

Simplicial complexes are systematic collections of these building blocks, as illustrated in Figure 3a–c. The dimension of a simplicial complex is the maximum dimension of any of its building blocks (simplices). For example, a network can be considered as a 1-dimensional simplicial complex with nodes (points, $k=0$) and edges (lines, $k=1$).

Various approaches are available for constructing simplicial complexes from initially unconnected data points, such as those we obtain when considering each sample as a point in a multidimensional space (see Section 3.1). We will use the Vietoris–Rips (VR) [25] complex method, which is commonly preferred in PH calculations for TDA due to its computational advantages. Essentially, the VR method assembles data points into higher-order objects if balls (circles) of a given threshold radius $r$, centered at the data points, overlap each other: i.e., their pairwise distances do not exceed $2r$, as illustrated in Figure 3d. For more details and an alternative method, see Appendix A.

Some key characteristics of simplicial complexes can be expressed by means of a set of topological invariants called Betti numbers. The $k$-th Betti number of a simplicial complex $K$ is denoted as ${\beta }_k\left(K\right)$, with $k=0$ measuring the number of 0-dimensional holes (connected components), $k=1$ measuring the number of 1-dimensional or “circular” holes (1-cycles), and $k=2$ measuring the number of 2-dimensional holes (“cavities” or “voids”). A formal definition of ${\beta }_k\left(K\right)$ is provided in Appendix A.

The “holes” (of different dimensions) measured by these Betti numbers—i.e., the connected components, the circular holes, and the cavities—are referred to as “topological features” of a simplicial complex.

2.2. Persistent Homology

For a point cloud representing a set of data points in a high-dimensional space, classical homology, as outlined above, fails to give meaningful information. This limitation arises because when constructing simplicial complexes from individual data points, a dataset’s topology can vary significantly depending on the choice of the distance threshold parameter $r,$ as can be seen in Figure 3d. Thus, determining a single representative threshold for the entire dataset can be a challenging or even an impossible task. To address this issue, the PH method offers an alternative approach.

Indeed, PH considers multiple or all values of the thresholding parameter $r$. As $r$ increases, initially disconnected simplices can merge to form new, larger simplices (like edges formed from individual data points). These new simplices are then added to the current simplicial complex, resulting in a sequence of subcomplexes known as a “filtered simplicial complex” (see Figure 3d and Appendix A for more details).

PH involves tracking when (i.e., at what threshold parameter $r$) new topological features (e.g., circular holes) emerge (referred to as ‘birth’) and when they disappear (referred to as ‘death’) by recording the corresponding distance thresholds $r$. This process captures how the topological properties of the complexes evolve as $r$ increases.

Various methods in computational topology allow the summary and visualization of PH and multi-scale topological features, including persistence barcodes [26], persistence diagrams [27,28], persistence landscapes [29,30], persistence images [31], and persistence curves [32], each offering insights into data topology for different applications. In Figure 4, we provide an illustrative example for a given point cloud (Figure 4a).

Both persistence barcodes (Figure 4b) and persistence diagrams (Figure 4c) offer a clear representation of the lifespan (persistence) of each feature. Essentially, a persistence barcode portrays each feature as a bar, with the length of the bar reflecting the feature’s lifespan, while a persistence diagram visually represents each topological feature as a point on a graph, with the $x$- and $y$-axes denoting the distance threshold parameter $r$ at the feature’s birth and death, respectively. In both visualization techniques, the color of the bar (in the persistence barcode) or point (in the persistence diagram) can indicate the dimension $k$ of the corresponding topological feature.

Longer bars in the persistence barcode (or points located farther from the diagonal in the persistence diagram) are typically considered robust features of the dataset, while shorter bars in the persistence barcode (or points near the diagonal in the persistence diagram) are often interpreted as noise or less significant features [9]. Thus, PH is robust against data noise, making it a stable representation, where small perturbations in the data have only a minor effect on the barcode or persistence diagram.

Additionally, the evolution of Betti numbers as $r$ increases can be visualized as Betti curves (see Figure 4d). For example, consider the topological features shown at three thresholds: $r=0.25, r=0.375,$ and $r=0.5$. At $r=0.25, { \beta }_0=6$, because there are only six connected components; the two closest points at the bottom left corner of Figure 4a have merged into one connected component, while ${\beta }_1=0$ due to the absence of 1-dimensional circular holes, i.e., 1-cycles (see also Figure 3d, third panel). At $r=0.375, { \beta }_0=4$, because there are only four connected components; the three closest points at the top right corner of Figure 4a have merged into one connected component, while ${\beta }_1=0$ because there still is no 1-cycle (see also Figure 3d, fourth panel), whereas at $r=0.5, { \beta }_0=1,$because the threshold radius$r$is large enough to merge all points into a single connected component, while ${\beta }_1=1$ due to the presence of a 1-cycle (see also Figure 3d, last panel).

3. Materials and Methods

As the research pipeline was previously outlined in Figure 1 and a more detailed implementation provided in Algorithm 1. Detailed information about the expression and pathway datasets used, the topological descriptors explored, and the methods employed in this study are explained in this section.

Algorithm 1: Persistent Homology Analysis of Gene Expression Data

Input:   Gene expression dataset (control and disease samples). Output:   Persistence diagrams, persistence barcodes, and Betti curves. Significant topological descriptors indicating differences between sample classes. Significant KEGG pathways exhibiting topological changes similar to those observed for the genome-wide gene expression dataset. PartitionDataset(GeneExpressionDataset) i. Separate the gene expression dataset into control (healthy) and disease (HCC) sample classes. For each sample class: ConstructPointCloud(GeneExpressionClass) i. Create a multidimensional point cloud, where each point represents the gene expression levels of a specific sample from the GeneExpressionClass. ii. Measure distances between points (samples) based on correlation coefficients. ComputePersistentHomology(PointCloud) i. Identify persistent features using the point cloud. ii. Generate persistence diagrams and persistence barcodes representing birth and death of topological features. ExtractTopologicalDescriptors(PersistenceHomology) i. Measure all the topological descriptors defined in Section 3 from the persistence homology computed for the sample class. AssessDifferences(TopologicalDescriptorsControl, TopologicalDescriptorsDisease) i. Randomly permute the class labels of the samples and recompute topological descriptors (steps 1 and 2). Perform this random permutation 2000 times to obtain a distribution of random values for each topological descriptor. ii. Compare the observed/actual values of each descriptor with their corresponding values obtained from permuted/randomized class labels to determine statistical significance. IdentifySignificantDescriptors(TopologicalDescriptor) i. Identify topological descriptors showing significant differences between control and disease classes. ApplyWorkflowToSubsets(GeneExpressionDataset, PathwaySubsetGeneSets) i. For the pathway-specific analysis, the focus is on subsets of genes associated with each particular pathway. ii. Apply steps 1–4 above to each KEGG pathway gene set, limiting the gene expression dataset to the subset of genes involved in the given pathway. iii. Identify pathways where the topological descriptors, specifically those that showed significance in the genome-wide gene expression dataset at step 4, also demonstrate significance.

3.1. PBMC Gene Expression Data

Gene expression dataset. The data were obtained from a study by Han et al. [21], which employed RNA-seq to profile gene expression patterns in PBMC from 17 patients with hepatocellular carcinoma and 17 healthy controls of similar age. Of 30,474 genes with discernible gene expression levels in the PBMC samples, as selected by the original authors, we used only those 26,575 genes that have a corresponding HGNC gene symbol.

RNA-seq data were generated by Han et al. based on the following experimental technique: RNA transcripts (molecules) extracted from cells are converted to their DNA equivalent and subsequently fractured into smaller pieces which are sequenced. The obtained sequence “reads” are then aligned/mapped to the reference genome or transcriptome. Since the number of RNA transcripts of a gene present in the cells determines the number of obtained sequence reads for that gene, gene expression levels can be estimated from the number of reads which map to the corresponding genomic regions.

Data pre-processing and normalization. As suggested by a recent validation by Abrams et al. [33], we employed the Transcripts Per Million (TPM) method [34] to normalize the RNA-seq gene expression data using ‘edgeR’ [35] (R package version 3.42.4, R version 4.3.0). Afterwards, a log2 transformation was carried out on the resultant TPM values.

Data representation for genome-wide gene expression. Each sample was treated as a data point (sample vector) within a 26,575-dimensional space, defined by the expression levels of the individual genes in that sample. Thus, the set of HCC or control samples form a ’point cloud’ to which PH can be applied.

Data representation for pathway-specific gene expression. For gene sets from individual cellular pathways, we essentially used the data representation in a point cloud but defined each data point (sample vector) only by the expression levels of the genes involved in each pathway, thus reducing the 26,575-dimensional space to an $N$-dimensional space, where $N$ is the number of genes involved in the pathway.

Distance metric. To measure the distances between samples, i.e., between the data points in the point cloud, different distance metrics can be used, such as the Euclidean distance or a distance based on a correlation coefficient. In biomedical research, the similarity (or dissimilarity) of samples is usually measured in terms of the positive (or negative) correlation of their gene expression profiles. Therefore, we defined the distance between two samples as $1-\rho$, where $\rho$ is the Pearson correlation coefficient between the samples. Since $\rho$ can range from −1 (perfect anti-correlation) to 1 (perfect correlation), the distance can theoretically range from 0 to 2.

3.2. KEGG Pathways

We obtained a list of cellular pathways and their associated genes from the Kyoto Encyclopedia of Genes and Genomes (KEGG) [36] database using the R (version 4.3.0) packages, ‘AnnotationDbi’ (version 1.62.2), and ‘org.Hs.eg.db’ (version 3.17.0). For our study, we focused exclusively on the 251 metabolic pathways and signaling pathways.

3.3. PH Implementation

In this study, we employed PH as the chosen topological method. We categorized the gene expression dataset into two classes consisting of control samples only and disease samples only, and applied PH to both classes individually. We constructed the VR complexes from the corresponding sample point clouds, using the distance metric defined above. The Python version of ‘Gudhi’ [37] (version 3.7.1) (https://gudhi.inria.fr/ (accessed on 20 December 2023)) was used for the computation of PH and for generating our visualizations.

3.4. Topological Descriptors

Let $k$ represent the dimension of a $k$-simplex in an $n$-dimensional simplicial complex. We defined several key descriptors for the topological features in our analysis. Specifically, $m_k$ represents the count of $k$-dimensional topological features, $b_{i,k}$ denotes the “birth time” of the $i$-th $k$-dimensional feature (threshold radius $r$ at which the feature appears), $d_{i,k}$ denotes the “death time” of the $i$-th $k$-dimensional feature ($r$ at which the feature ceases to exist), and $p_{lk}$ corresponds to the total (summed) persistence (lifetime) of all features of dimension $k$. The topological descriptors we consider for our analysis are explicitly defined as shown in

Table 1. Topological descriptors.

Descriptor (Short Form)	Formula
Classical Euler characteristic (Classical EC)	${\sum }_{k=0}^n{(-1)}^k{m}_k$
Persistence-wise Euler characteristic (Persistence EC)	${\sum }_{k=0}^n{(-1)}^k{p}_{lk}$
Sum of persistence of $k$-dimensional features (Sum P-$k$)	$p_{lk}={\sum }_{i=1}^{m_k}(d_{i,k}-b_{i,k})$
Average persistence of $k$-dimensional features (Average P-$k$)	$\frac{1}{m_k}{\sum }_{i=1}^{m_k}(d_{i,k}-b_{i,k})$
Maximum persistence of $k$-dimensional features (Max P-$k$)	$\underset{i}{\max }\{d_{i,k}-b_{i,k}\}$
Range of persistence of $k$-dimensional features (Range P-$k$)	$\underset{i}{\max }\{d_{i,k}-b_{i,k}\}-\underset{i}{\min }\{d_{i,k}-b_{i,k}\}$
Sum of birth times of $k$-dimensional features (Sum BT-$k$)	${\sum }_{k=1}^{m_k}{b}_{i,k}$
Average birth times of $k$-dimensional features (Average BT-$k$)	$\frac{1}{m_k}{\sum }_{k=1}^{m_k}{b}_{i,k}$
Sum of death times of $k$-dimensional features (Sum DT-$k$)	${\sum }_{k=1}^{m_k}{d}_{i,k}$
Average death times of $k$-dimensional features (Average DT-$k$)	$\frac{1}{m_k}{\sum }_{k=1}^{m_k}{d}_{i,k}$

. It is worth noting that most of these descriptors apply to topological features ($k$-cycles) of specific dimension $k$ within a persistence diagram. Only two descriptors, namely the classical and persistence-wise Euler characteristics, take into account the contribution of topological features ($k$-cycles) across all possible dimensions (e.g., connected components, 1-cycles, cavities, etc.).

It is important to highlight that a value of zero for a topological descriptor of dimension $k$ typically indicates the absence of features in that particular dimension. Exceptions to this are Classical EC and Persistence EC, where a value of zero does not necessarily indicate the absence of features due to the involvement of multiple dimensions in their computation.

3.5. Differential Expression Analysis

For the classical approach to pathway analysis, differential expression analysis was performed using the ‘PyDESeq2’ [38] package (version 0.4.4) for Python (version 3.10.7). Only genes with an adjusted $p$-value $<0.05$(adjusted for multiple testing using the Benjamini–Hochberg method) and a log fold change $\ge 1$ were considered as significantly differentially expressed.

3.6. Enrichment Analysis of Pathways

The classical enrichment-based pathway analysis was carried out using the ‘Scipy’ [39] package (version 1.10.1) for Python (version 3.10.7). Essentially, it performs a hypergeometric enrichment test (one-sided Fisher’s exact test) to identify pathways containing a significantly high number of differentially expressed genes. For each pathway, it calculates a $p$-value, which in essence is the likelihood of observing at least $m$ differentially expressed genes in a pathway purely by chance, considering the total number of genes in the dataset, the total number of genes in the pathway, and the total number of differentially expressed genes. $p$-values were corrected for multiple testing (251 pathways) using the Benjamini–Hochberg method.

3.7. Significance of Topological Descriptors

To assess whether topological descriptors of the PH-based approach differ significantly between sample classes, we performed a two-tailed statistical permutation test using Python (version 3.10.7) and the Benjamini–Hochberg correction of $p$-values for multiple testing using the ‘Statsmodels’ [40] package (version 0.13.5). Essentially, 2000 random permutations of the sample labels were generated, and there were correspondingly 2000 random assignments of samples to the disease and control classes.

For each of the permutations, we applied PH and computed the values of the topological descriptors individually for the two sample classes. For each topological descriptor TD, we computed the difference between disease and control samples as: TD(disease)—TD(control). We repeated this entire process for all 2000 permutations to build a distribution of differences under the null hypothesis that there is no association (or biological distinction) between the peripheral blood samples of the two classes.

To quantify significance, for each descriptor, we calculated the empirical or permutation $p$-value, which measures the proportion of random permutations that yield a more extreme value (in either direction) than the observed (actual) value. Subsequently, to account for multiple testing, we applied the Benjamini–Hochberg procedure, which resulted in adjusted $p$-values. These adjusted $p$-values help control the false discovery rates (FDR), indicating the likelihood of false discoveries (incorrectly identifying differences) and potential for false positives (attaining statistically significant results by chance) among the descriptors. Throughout our analysis, we set the FDR threshold to be less than 0.05.

4. Results

4.1. Genome-Wide Persistent Homology Analysis

Initially, we leveraged the full (genome-wide) HCC RNA-seq gene expression dataset from the study by Han et al. [21] (see Section 3 for details) to identify key topological descriptors that effectively capture changes in genome-wide gene expression between 17 HCC samples and 17 healthy controls.

The process, outlined in Figure 1 and Algorithm 1, begins by partitioning the genome-wide gene expression data into datasets for control and disease samples. PH is then individually applied to each class.

For each sample class, the genome-wide gene expression levels of each individual sample are represented as a data point in the input point cloud, used for the construction of the VR complexes and the computations of their persistent homologies. Distances between data points are measured based on the Pearson correlation coefficient between the corresponding samples, as detailed in Section 3.1. Subsequently, topological descriptors for both sample classes were derived from the persistence diagrams obtained through PH calculations (see Section 3.4 for details about each descriptor). To visualize the topological (homological) evolution resulting from the computation of PH, we created and plotted persistence barcodes, persistence diagrams, and Betti curves for each class (see Figure 5).

We noticed that the most persistent features were the 0-dimensional features (i.e., connected components). This observation remains valid irrespective of whether the data were from the control class or the disease class. Additionally, 1-dimensional topological features (i.e., 1-cycles/circular holes) and 2-dimensional topological features (i.e., 2-cycles/voids) were also present although not as persistent as the 0-dimensional topological features. No 3-dimensional topological features (3-cycles) were detected, which allowed us to omit the calculation of topological features of higher dimensions (3 and above), resulting in a substantial reduction in computation time without compromising the quality of the output represented by the computed topological descriptors. In addition, we observed that since most of the 1-dimensional and 2-dimensional topological cycles were not very persistent (had short lifespans), the value of Persistence EC was very close to that of Sum P-0 (see

Table 2. Significant topological descriptors and differences between control and disease classes for the genome-wide gene expression dataset.

Descriptors	Dataset		Difference (Disease—Control)	$\mathbf{adj}. \mathit{p}$-Value
	Disease Class	Control Class
Persistence EC	0.6252	0.3509	0.2743	<0.00050
Sum P-0	0.6305	0.3573	0.2732	<0.00050
Average P-0	0.0371	0.0210	0.0161	<0.00050
Range P-2	0.0000	0.0003	−0.0003	0.03667
Sum BT-2	0.0000	0.0569	−0.0569	0.03667
Sum DT-2	0.0000	0.0577	−0.0577	0.03667

and Supplementary Table S1).

We then assessed through a permutation test whether the computed values of a topological descriptor significantly differed between disease and control samples. To achieve this, we first computed the values of all the considered topological descriptors for both classes and then determined the difference between the two classes for each descriptor (see Section 3.7). These differences served as our measure of the global changes in gene expression observable in HCC peripheral blood samples with respect to controls.

4.2. Choice of Relevant Topological Descriptors

To identify the subset of topological descriptors that best reflect genome-wide gene expression changes in HCC samples, we performed a two-tailed permutation test on the sample labels (random assignment of samples to disease and control classes) to evaluate the statistical significance of the observed differences in topological descriptors between peripheral blood from HCC patients and healthy controls.

Several of the topological descriptors showed highly significant differences (FDR below 0.05) between control and disease samples (see Table 2 and Figure A3, as well as Supplementary Table S1 and Figure S1 for additional details). Specifically, we have identified the persistence-wise Euler characteristic (Persistence EC), the sum of persistence of 0-dimensional features (Sum P-0), the average persistence of 0-dimensional features (Average P-0), the range of persistence of 2-dimensional features (Range P-2), the sum of birth times of 2-dimensional features (Sum BT-2), and the sum of death times of 2-dimensional features (Sum DT-2) as the descriptors meeting our statistical significance criteria. Hence, we selected them as representative of the differences in peripheral blood between HCC patients and healthy controls.

Furthermore, we made noteworthy observations regarding these descriptors. Specifically, we found that Persistence EC, Sum P-0, and Average P-0 exhibited higher values in the disease class compared to the control class, while Range P-2, Sum BT-2, and Sum DT-2 had higher values in the control class because topological features of the dimensions $k=2$ are not present in the PH of the disease class. The observation of higher values for topological descriptors in dimension $k=0$ (namely, Sum P-0 and Average P-0) indicates a greater sample heterogeneity in the disease class compared to the control class, similar to what has been observed for autism spectrum disorder in a previous study [41].

In summary, we identified a subset of six topological descriptors (see Table 2) that best reflect genome-wide gene expression changes in peripheral blood from HCC patients.

4.3. Classical Differential Gene Expression Analysis and Enrichment-Based Pathway Analysis

Before using our new PH-based approach for pathway analysis, we applied the most frequently used classical pathway analysis method, so that we would be able to compare our results to those which can be obtained using a standard approach.

First, we used the genome-wide expression data and identified 1426 differentially expressed genes (adjusted $p$-value (FDR) of less than 0.05 and log-fold change of at least 1). Among these significantly dysregulated genes, 1242 genes were up-regulated, while 184 genes were down-regulated in peripheral blood of HCC patients (see Supplementary File S2 for details).

Afterwards, we assessed the significance of 251 KEGG metabolic pathways and signaling pathways by performing a hypergeometric enrichment test (one-sided Fisher’s exact test; see Section 3) to identify pathways containing a significant number of differentially expressed genes. Only sixteen pathways exhibited a significant enrichment in differentially expressed genes with an adjusted $p$-value (FDR) of less than 0.05 (see Supplementary Table S2 for details).

4.4. Persistent Homology Analysis for Pathway-Specific Gene Sets

Subsequent to the selection of six topological descriptors indicative of genome-wide gene expression patterns in peripheral blood from HCC patients, we investigated whether there are specific cellular pathways which reflect these global changes, hypothesizing that such pathways might potentially aid in an early detection of HCC from peripheral blood.

To this end, we applied our PH pipeline to various subsets of genes involved in specific cellular pathways, taking only the corresponding gene expression data instead of the full (genome-wide) dataset. That is, using functional gene sets for 251 metabolic pathways and signaling pathways from the KEGG database [36], we applied essentially the same methodology that we had used for the genome-wide gene expression dataset (see above) to focus only on the expression levels of the sets of genes involved in individual KEGG pathways to define the sample data points in the point cloud (see Section 3.1).

Our primary objective here was to assess the significance of pathways based on whether the previously identified topological descriptors showed significant differences between HCC samples and control samples—limited to the gene expression of the corresponding pathway gene sets—because we reasoned that those pathways which yield significant changes for the same topological descriptors would best reflect the genome-wide differences between HCC samples and those from healthy controls.

First considering individual topological descriptors, we observed that Persistence EC and Average P-0 detected the same set of pathways (138 pathways), while Sum BT-2 and Sum DT-2 produced identical results (36 pathways). In addition, Sum P-0 identified 137 pathways and Range P-2 identified 181 pathways, respectively (see Supplementary File S1 for a detailed list).

Visualizations in Figure 6 illustrate the number of pathways for which the individual topological descriptors show significant differences between HCC samples and controls, as well as the overlap of these sets of pathways with the 16 pathways identified as enriched in differently expressed genes.

Since our purpose is to identify pathways that best reflect the genome-wide changes in gene expression, instead of considering pathways with respect to individual topological descriptors, we define a pathway as significant when it exhibits significant changes across all six topological descriptors, mirroring the observations made for the genome-wide gene expression dataset. This criterion was met by twenty-three pathways, as illustrated in Figure 7a and detailed in

Table 3. Significant KEGG pathways and the differences and adjusted $p$-values of individual topological descriptors between control and disease sample classes. A difference of zero with a significant $p$-value indicates that the corresponding topological features were absent in PH of both HCC and control sample classes, but not when sample classes were randomized.

Pathway (KEGG ID)	Descriptor	Difference (Disease—Control)	$\mathbf{adj}. \mathit{p}$-Value
ABC transporters (02010)	Persistence EC	1.0659	<0.0005
	Sum P-0	1.0650	<0.0005
	Average P-0	0.0626	<0.0005
	Range P-2	0.0000	<0.0005
	Sum BT-2	0.0000	0.0240
	Sum DT-2	0.0000	0.0240
Apelin signaling pathway (04371)	Persistence EC	0.3427	0.0220
	Sum P-0	0.3429	0.0220
	Average P-0	0.0202	0.0220
	Range P-2	0.0000	0.0220
	Sum BT-2	0.0383	0.0318
	Sum DT-2	0.0385	0.0318
Ascorbate and aldarate metabolism (00053)	Persistence EC	0.4311	0.0330
	Sum P-0	0.4293	0.0330
	Average P-0	0.0253	0.0330
	Range P-2	0.0000	<0.0005
	Sum BT-2	0.0000	0.0420
	Sum DT-2	0.0000	0.0420
Base excision repair (03410)	Persistence EC	0.6788	<0.0005
	Sum P-0	0.6727	<0.0005
	Average P-0	0.0396	<0.0005
	Range P-2	−0.0007	0.0220
	Sum BT-2	−0.0957	0.0283
	Sum DT-2	−0.0972	0.0283
beta-Alanine metabolism (00410)	Persistence EC	0.3086	0.0440
	Sum P-0	0.3157	0.0440
	Average P-0	0.0186	0.0440
	Range P-2	0.0000	<0.0005
	Sum BT-2	−0.0357	0.0440
	Sum DT-2	−0.0366	0.0440
Citrate cycle (TCA cycle) (00020)	Persistence EC	1.1207	0.0110
	Sum P-0	1.1070	0.0110
	Average P-0	0.0652	0.0110
	Range P-2	−0.0020	0.0110
	Sum BT-2	−0.1709	0.0477
	Sum DT-2	−0.1759	0.0477
Collecting duct acid secretion (04966)	Persistence EC	1.1848	<0.0005
	Sum P-0	1.1795	<0.0005
	Average P-0	0.0693	<0.0005
	Range P-2	0.0000	<0.0005
	Sum BT-2	0.0000	0.0040
	Sum DT-2	0.0000	0.0040
Drug metabolism—cytochrome P450 (00982)	Persistence EC	0.4315	0.0320
	Sum P-0	0.4568	0.0320
	Average P-0	0.0268	0.0320
	Range P-2	0.0000	<0.0005
	Sum BT-2	0.0000	0.0320
	Sum DT-2	0.0000	0.0320
Glycine, serine, and threonine metabolism (00260)	Persistence EC	1.6059	<0.0005
	Sum P-0	1.5992	<0.0005
	Average P-0	0.0940	<0.0005
	Range P-2	0.0000	<0.0005
	Sum BT-2	0.0000	0.0260
	Sum DT-2	0.0000	0.0260
Histidine metabolism (00340)	Persistence EC	0.8008	0.0140
	Sum P-0	0.8045	0.0140
	Average P-0	0.0473	0.0140
	Range P-2	0.0000	<0.0005
	Sum BT-2	0.0000	0.0140
	Sum DT-2	0.0000	0.0140
IL-17 signaling pathway (04657)	Persistence EC	0.3977	0.0360
	Sum P-0	0.3962	0.0360
	Average P-0	0.0233	0.0360
	Range P-2	0.0000	<0.0005
	Sum BT-2	0.0000	0.0360
	Sum DT-2	0.0000	0.0360
p53 signaling pathway (04115)	Persistence EC	0.4171	0.0165
	Sum P-0	0.4117	0.0165
	Average P-0	0.0242	0.0165
	Range P-2	0.0000	<0.0005
	Sum BT-2	0.0000	0.0385
	Sum DT-2	0.0000	0.0385
Pantothenate and CoA biosynthesis (00770)	Persistence EC	0.3973	0.0400
	Sum P-0	0.3995	0.0400
	Average P-0	0.0235	0.0400
	Range P-2	0.0000	<0.0005
	Sum BT-2	0.0000	0.0400
	Sum DT-2	0.0000	0.0400
Phosphonate and phosphinate metabolism (00440)	Persistence EC	3.2065	0.0200
	Sum P-0	3.2170	0.0154
	Average P-0	0.1892	0.0154
	Range P-2	0.0000	<0.0005
	Sum BT-2	0.0000	0.0055
	Sum DT-2	0.0000	0.0055
Porphyrin metabolism (00860)	Persistence EC	2.3922	<0.0005
	Sum P-0	2.4252	<0.0005
	Average P-0	0.1426	<0.0005
	Range P-2	0.0000	<0.0005
	Sum BT-2	0.0000	0.0100
	Sum DT-2	0.0000	0.0100
Primary bile acid biosynthesis (00120)	Persistence EC	0.3730	0.0200
	Sum P-0	0.3737	0.0200
	Average P-0	0.0220	0.0200
	Range P-2	0.0000	<0.0005
	Sum BT-2	0.0000	0.0200
	Sum DT-2	0.0000	0.0200
Protein processing in endoplasmic reticulum (04141)	Persistence EC	0.5312	<0.0005
	Sum P-0	0.5237	<0.0005
	Average P-0	0.0308	<0.0005
	Range P-2	−0.0003	0.0176
	Sum BT-2	−0.0638	0.0220
	Sum DT-2	−0.0643	0.0220
Riboflavin metabolism (00740)	Persistence EC	1.5118	<0.0005
	Sum P-0	1.5413	<0.0005
	Average P-0	0.0906	<0.0005
	Range P-2	0.0000	<0.0005
	Sum BT-2	0.0000	0.0495
	Sum DT-2	0.0000	0.0495
RNA polymerase (03020)	Persistence EC	0.3368	0.0377
	Sum P-0	0.3691	0.0330
	Average P-0	0.0218	0.0330
	Range P-2	0.0000	0.0330
	Sum BT-2	0.1141	0.0460
	Sum DT-2	0.1171	0.0460
Sulfur metabolism (00920)	Persistence EC	6.0523	<0.0005
	Sum P-0	6.0452	<0.0005
	Average P-0	0.3556	<0.0005
	Range P-2	0.0000	<0.0005
	Sum BT-2	0.0000	0.0220
	Sum DT-2	0.0000	0.0220
Synaptic vesicle cycle (04721)	Persistence EC	0.3297	<0.0005
	Sum P-0	0.3305	<0.0005
	Average P-0	0.0195	<0.0005
	Range P-2	−0.0006	0.0385
	Sum BT-2	−0.0523	0.0403
	Sum DT-2	−0.0534	0.0403
Tryptophan metabolism (00380)	Persistence EC	0.2584	0.0424
	Sum P-0	0.2530	0.0424
	Average P-0	0.0148	0.0424
	Range P-2	0.0000	<0.0005
	Sum BT-2	0.0000	0.0424
	Sum DT-2	0.0000	0.0424
Virion—Herpesvirus (03266)	Persistence EC	1.1037	0.0220
	Sum P-0	1.1240	0.0220
	Average P-0	0.0661	0.0220
	Range P-2	0.0000	<0.0005
	Sum BT-2	0.0000	0.0220
	Sum DT-2	0.0000	0.0220

.

It is important to note that for 2-dimensional features (cavities or voids), the associated topological descriptors can also be significant if their difference between HCC samples and controls is zero, indicating that these features are absent in the PH of both types of samples. When randomizing sample classes, instead, these features can appear and thus lead to a non-zero difference, such that a zero difference obtained with the correct sample classes may indeed be statistically significant.

The pathways identified as relevant by our PH method include genetic information processing pathways (RNA polymerase, protein processing in endoplasmic reticulum, and base excision repair), environmental information processing pathways (ABC transporters and apelin signaling pathway), a cellular process pathway (p53 signaling pathway), organismal system pathways (IL-17 signaling pathway and collecting duct acid secretion), and several metabolic pathways (riboflavin metabolism, citrate/TCA cycle, sulfur metabolism, ascorbate and aldarate metabolism, drug metabolism—cytochrome P450, glycine, serine and threonine metabolism, primary bile acid biosynthesis, phosphonate and phosphinate metabolism, histidine metabolism, tryptophan metabolism, beta-alanine metabolism, pantothenate and CoA biosynthesis, and porphyrin metabolism); refer to Table 3 for additional details.

4.5. Comparison of Pathways Identified by the PH Method and the Enrichment Analysis

We compared our novel approach with the commonly used enrichment-based pathway analysis method, which identifies significant pathways through the hypergeometric enrichment test of differentially expressed genes in a pathway. Notably, our results revealed mostly unique significant pathways obtained by these approaches.

Among the 23 significant pathways detected by our PH method, only the IL-17 signaling pathway was also identified by the classical enrichment-based method (see Figure 8). Remarkably, 22 pathways were unique to our PH method (see Figure 8). Although these pathways did not show a significant enrichment in differentially expressed genes, they provide valuable insights into the molecular processes that may help to detect HCC from peripheral blood samples, deserving further investigation, experimental validation, and functional studies. Similarly, the pathways uniquely identified by the hypergeometric enrichment test (see Supplementary File S1 and Table S2 for details) but not detected by our PH method should not be disregarded.

Interestingly, several pathways identified by our PH method have been previously associated with HCC, or their association has been suggested by other computational approaches. For completeness, we provide a list: genetic information processing pathways (RNA polymerase [42], protein processing in endoplasmic reticulum [43], and base excision repair [44,45]), environmental information processing pathways (ABC transporters [46,47,48] and apelin signaling pathway [49,50]), a cellular processes pathway (p53 signaling pathway [51,52,53,54]), organismal systems pathways (IL-17 signaling pathway [55,56] and synaptic vesicle cycle [57]), and several metabolism pathways (riboflavin metabolism [58], citrate/TCA cycle [59,60,61], ascorbate and aldarate metabolism [62,63], drug metabolism—cytochrome P450 [48,64,65,66,67], glycine, serine, and threonine metabolism [48,64,68], primary bile acid biosynthesis [64], histidine metabolism [64], beta-alanine metabolism [64], tryptophan metabolism [54,64,69,70], pantothenate and CoA biosynthesis [58,71,72,73], and porphyrin metabolism [58,74,75]). This suggests that our method may indeed be capable of detecting subtle, HCC-related changes from peripheral blood.

Fifteen pathways were enriched in differentially expressed genes but not detected as significant by our PH method (see Figure 8). These include complement and coagulation cascades, the PPAR signaling pathway, hematopoietic cell lineage, the cholesterol metabolism, neutrophil extracellular trap formation, osteoclast differentiation, ECM–receptor interaction, platelet activation, cytokine–cytokine receptor interaction, neuroactive ligand–receptor interaction, the phagosome, fat digestion and absorption, the glycerolipid metabolism, the JAK-STAT signaling pathway, and the calcium signaling pathway (see Supplementary Table S2 for details).

5. Discussion

Detecting the presence of a severe disease such as HCC from the peripheral blood sample of an individual necessitates the identification of key biomarkers, be they specific metabolites, proteins, genes, or pathways. Unlike previous studies that primarily focused on applications of PH for tasks such as classification, prediction, or clustering [76,77,78], our study demonstrates the effectiveness of utilizing PH to explore the global characteristics of RNA-seq expression data from peripheral blood of both HCC patients and healthy individuals to identify significant pathways.

We first identified a set of topological descriptors that showed significant differences with respect to genome-wide expression data of HCC samples and normal controls. Persistence EC assesses the topological complexity of a pathway by examining persisting features across multiple dimensions during various stages of filtration, i.e., at various threshold radii $r$. Sum P-0, which is equivalent to the sum of death times of 0-dimensional features (connected components), provides insights into their persistence, reflecting the heterogeneity of gene expression between samples. Longer persistence indicates larger distances in the data point cloud and hence lower correlation between samples. Additionally, Average P-0 captures the mean persistence of 0-dimensional features. On the other hand, Range P-2, Sum BT-2, and Sum DT-2 elucidate the evolution of 2-dimensional features, emphasizing the birth, persistence, and demise of voids or cavities within the topology of the samples’ point cloud.

We then applied the PH method to pathway-specific gene expression data rather than genome-wide data, reasoning that pathways that reflect global changes may be particularly interesting as potential biomarker pathways. Among the twenty-three identified significant pathways, the IL-17 signaling pathway, also detected by the classical enrichment analysis method, is particularly noteworthy.

Figure 9 illustrates how pathways that were either identified as significant by our PH method, or were not evidenced as such, do or do not reflect patterns observed from genome-wide gene expression data.

HCC-related samples generally showed more heterogeneity than control samples. This is suggested by 0-Betti numbers ${\beta }_0$, i.e., the numbers of connected components (but disconnected between each other), which decline much slower with increasing threshold radius $r$, indicating an overall lower pairwise correlation between samples.

The pathways identified as significant by our method displayed a pattern resembling more closely the genome-wide gene expression dataset (Figure 9a). Specifically, for significant pathways like the citrate cycle/TCA cycle (Figure 9b) and the p53 signaling pathway (Figure 9c), the 0-Betti numbers in the control class (blue curves) exhibited a sudden and faster drop than in the disease class (red curves), where the decline was more gradual. Conversely, pathways not identified by our method, such as the fructose and mannose metabolism (Figure 9d) and the phospholipase D signaling pathway (Figure 9e), showed a simultaneous but gradual decline phase for both sample classes. This observation suggests that topological features, such as the number of connected components, evolve differently as the radius increases.

For many of the 23 pathways highlighted by our method, we found indications for an experimentally validated or at least predicted association with HCC in the literature. The identified pathways, particularly those exclusive to our method, despite not being enriched in differentially expressed genes, offer a valuable resource for further exploration and experimental validation as biomarkers for peripheral blood. The importance of the p53 signaling pathway (see Figure 9c), for example, is supported both methodologically and biologically [51,52,53,54].

6. Conclusions and Future Research

Our study represents a pioneering effort in applying persistent homology (PH) analysis to identify tumor-associated pathways, particularly from PBMC samples of HCC patients. Analyzing RNA-seq data from peripheral blood, we successfully identified 23 significant pathways, including the apelin signaling pathway, the IL-17 signaling pathway, and the p53 signaling pathway. Our findings suggest that PH-based methods are complementary to classical enrichment-based pathway analysis methods and could thus be used in conjunction with current approaches to identify relevant pathways for experimental validation. We anticipate that applying our method directly to samples from affected disease tissues, rather than peripheral blood, could further advance our understanding of complex diseases.

While TDA methods, including PH, have shown promise in various applications, their use in gene expression analyses is still in its early stages and, to our best knowledge, ours is the first application to the pathway analysis problem. Looking ahead, future research could refine and expand the utility of TDA methods, particularly PH, by exploring their integration with other approaches for the comprehensive analysis of complex biomedical data. This approach holds the potential to advance our understanding of complex diseases and ultimately improve diagnostic strategies. The gained knowledge, e.g., about potential biomarker pathways, may ultimately help to develop early detection strategies that can be applied to individual samples or patients.