SenPred: a single-cell RNA sequencing-based machine learning pipeline to classify deeply senescent dermal fibroblast cells for the detection of an in vivo senescent cell burden

Primary human dermal fibroblasts (HDFs) were a kind donation from anonymous healthy patients under standard ethical practice, reference LREC No. 09/HO704/69. HDFs were cultured in DMEM with 4 mM L-glutamine (41,966–029, Life Technologies) supplemented with 10% foetal bovine serum (1001/500, Biosera), in the absence of antibiotics, and maintained at 37 °C, 5% CO₂, 95% humidity. Cells were seeded at 4000 cells/cm² for passages 1–19, 7500 cells/cm² for passage 20–25, and 10,000 cells/cm² for passage 26 onwards. HDFs were serially passaged and were classified as early proliferative (EP) up to passage 13, and Deeply Senescent (DS) beyond passage 39 + 3. All cells were routinely tested for mycoplasma and were shown to be negative.

Dermal gels were generated by suspending HDFs in collagen 1 (354,236, Corning) supplemented with Matrigel (7,341,100, Corning), FCS (1001/500, Biosera), 10X DMEM (11,430–030, Invitrogen), and pH balanced with 4 M NaOH in the absence of antibiotics. HDF-gel mix was added dropwise to transwells (CC405, Corning) to give a final cell density of 10,000 HDFs per gel. Nine gels were constructed per condition. Dermal gels were maintained at 37 °C, 5% CO₂, and 95% humidity.

For 2D samples, cells were grown in 6-well plates at 7000/cm² for 72 h. Cells were collected, and resuspended in HBSS (14,025,092, Thermofisher), and immediately sorted into single cell oil droplet suspension. To extract the cells from 3D dermal gels, gels were incubated with collagenase D (Sigma) at 37 °C for 1 h. Collagenase activity was subsequently inhibited by HBSS addition. Gels were strained and centrifuged to capture cell pellets, which were resuspended in HBSS before sorting into single-cell oil droplet suspension.

Library generation and sequencing were performed using the v2 Chromium Single Cell 3′ Kit (10X genomics). Samples were sequenced using the Illumina NextSeq 500 High Output Run Sequencing platform using paired-end sequencing and aligned to the human genome (GRCh37) using Cell Ranger (10 × Genomics).

For external model testing, data was acquired from the published papers listed in Table 1.

The datasets from Chan et al. and Teo et al. were used to apply and test the 2D SenPred pipeline on in vitro scRNA-seq experiments, whereas the datasets from Tabib et al., Solé-Boldo et al., Ganier et al., and Sikkema et al. were used to test the 3D SenPred pipeline. The number of cells analysed for the in vitro studies was as follows:

The number of samples used for each in vivo study after QC and any additional sample details are listed in Tables 2, 3, 4, and 5. The ethics approval and data acquisition of human tissue were acquired by the respective authors within the in vivo references listed in Table 1.

Additional preprocessing steps were carried out on both the Tabib et al. [19] and Solé-Boldo et al. [20] datasets, to combine donors and generate a single Seurat Object for each. Cell level filtering and normalisation of raw data from Tabib et al. [19] and Solé-Boldo et al. [20] were performed using Seurat R package version 3.2.2 in R version 4.0.3 [26]. To remove low-quality cells, the data was filtered using the following metrics:

Data was normalised using ScTransform [27], and the sctransform normalized counts for each sample were then integrated [26], using the top 3000 most variable genes in the SelectIntegrationFeatures() function. Then, the FindIntegrationAnchors() function was applied with default parameters to identify the anchors across the different samples in each dataset. The anchors were then passed to the IntegrateData() function to create a single normalized and integrated Seurat object for each dataset.

All subsequent data analysis was performed using R (version 4.2.2) in Rstudio (2022.07.2 Build 576), running under macOS Monterey 12.0.1. Matrix, gene, and barcode.gz data files were loaded and converted into a Seurat object [28]. Unless otherwise stated, default parameters from Seurat V4.3.0.1 were applied. To remove low-quality cells, the data was filtered using the following metrics (Supplementary Fig. S1):

Data was normalised using NormalizeData(), which normalises individual gene expression of a cell to total gene expression of this cell and multiplies by the default scale factor (10,000). This data is then log transformed.

To subset fibroblasts from whole skin data for the Tabib [19] and Solé-Boldo [20] datasets, pre-defined metadata was used. For Ganier et al. [22], dot plots were generated to visualise the mean cluster expression of a pre-defined list of genes from the original paper, using the Seurat function DotPlot(). To identify lung fibroblasts in Sikkema et al., data was subset for pre-defined fibroblasts, from healthy donors, where the age was reported. Donors were removed if they had less than 100 fibroblasts recorded.

The top 2000 variable genes were identified using the FindVariableFeatures function [26], and these were subsequently scaled to ensure a mean expression of 0 and variance of 1 across all cells. K-nearest neighbour algorithms were used based on PCA distance, followed by the Louvain algorithm to determine clusters. The resolution of all cluster-based functions was determined systematically using the Clustree package [29]. These clusters were visualised on UMAP plots. To identify differentially expressed genes (DEGs) in the clusters, FindAllMarkers default Seurat function was implemented. Genes must be detected in > 25% of cells in at least 1 cluster and have a logFC greater than 0.25.

The ScPred package developed by the Powell group was used to build ML models to classify senescence (https://github.com/powellgenomicslab/scPred/)↗ [30]. Unless otherwise stated, default parameters of ScPred V1.9.2 were used. Data was separated into training (80%) and testing (20%) datasets. ScPred ML models were then built on training data using all principal components, using a variety of methods including Support Vector Machines with Radial Basis Function Kernel (SVM), Mixture Discriminant Analysis (MDA), K-nearest neighbours (KNN), and Generalised Linear Model (GLM). The models were evaluated based on their classification of the testing dataset through a variety of methods, including ROC, sensitivity, and specificity, and the best performing model (MDA) was taken forwards. This utilised the CRAN package mda_0.5–4. For the MDA model, a fivefold cross-validation was used. These models could then be applied to classify senescence in alternative datasets.

Monocle3 V1.3.1 was used to analyse the trajectory of the PDL50 cells [31–33]. All default parameters were used unless otherwise stated. The earliest principal node was calculated using the function get_earliest_principal_node, where the node most surrounded by EP cells was selected. The cells were then ordered based on pseudotime.

NABA_MATRISOME human gene set was obtained from Petrov et al., and the median log normalised gene count expression was calculated for the 2D, 3D, Tabib et al. [19] and Solé-Boldo et al. [20] datasets. Genes which had a median count of zero for all conditions were removed, and the final list of 147 matrisome genes is reported in Supplementary Table S2. The gene counts were converted into a heatmap using heatmap.2(), and Ward hierarchical clustering was applied.

Normalised gene counts were extracted from the Tabib and Solé-Boldo data, using a list of pre-defined senescence genes: CDKN1A, CDKN2A, IGFBP3, SERPINE1, IGF1, CXCL12, IL6, MMP3 and TP53. The Pearson correlation coefficient was calculated between all pairs of genes, and Ward hierarchical clustering was applied to generate both x- and y-axis dendrograms. The results were then visualised on a heatmap using heatmap.2().

Early Proliferative (EP – P10 − P12) and Deeply Senescent (DS – P39 + 3) FLF1 HDF cell pellets were lysed with QIAzol (217,004, Qiagen miRNeasy Mini Kit), periodically vortexed for 5 min and snap frozen at -80C. Total RNA extraction was performed using the miRNeasy Mini Kit (217,004, Qiagen), following the manufacturer’s guidelines (“Quick-Start protocol miRNeasy Mini Kit” guide, Qiagen). Following this, RNA ethanol precipitation was performed to increase sample purity, where 30 mL of RNA was mixed with 75 μl of 100% ethanol (34,852, Sigma Aldrich), 3 μl of sodium acetate 3 M at pH 5.2 (R1181, Thermofisher Scientific), and 1 μl GlycoBlueTM Coprecipitant (AM9515, Invitrogen), and incubated overnight at − 20 °C. The sample was centrifuged at 12,000 g for 10 min at 4 °C, resuspended in 500 μl of 80% EtOH, and centrifuged a second time at 12,000 g for 10 min at 4 °C. The pellet was dried for 2 h at RT and resuspended in 20 μl RNAse-free distilled water (217,004, Qiagen miRNeasy Mini Kit).

Precipitated RNA was converted to cDNA using the SuperScriptTM III Reverse Transcriptase kit, as per the manufacturer’s instructions (18,080,044, Invitrogen). cDNA concentration was measured using the Qubit 2.0 Fluorometer (Thermofisher Scientific), and absolute quantification of cDNA was carried out using SYBR Green I dye (S7563, Thermofisher Scientific) on the Applied Biosystems QuantStudio Absolute Q Digital PCR System. Digital PCR was carried out using the Absolute Q DNA digital PCR Master Mix (A52490, Thermofisher Scientific), with a QuantStudio Absolute Q MAP16 plate kit (A52865, Thermofisher Scientific), as per manufacturer’s instructions. For analysis, absolute cDNA counts were normalised to 1 ng/ml of cDNA loaded.

In an effort to develop a ML model of senescence classification, we cultured adult human dermal fibroblasts (HDFs) to replicative senescence (passage 39) and allowed them to reach a deeply senescent state over a further 3 weeks in culture (DS; P39 + 3). Single-cell RNA sequencing (scRNA-seq) was carried out on early proliferative (EP) and deeply senescent (DS) HDFs cultured in a two-dimensional monolayer (extensive model validation previously reported in Wallis et al. [34]) (Fig. 1a). Using k-nearest neighbours (KNN), seven independent clusters were identified in the EP and DS HDF dataset (Fig. 1b). Overlaying these clusters with the originating cell populations showed distinct clustering of EP and DS fibroblasts, highlighting the discrete mRNA profiles of the two populations, together with inherent heterogeneity within both the EP and DS fibroblast populations (Fig. 1c). Given that the EP and DS states were transcriptionally distinct, we proceeded to build a ML model to distinguish the two conditions.

The scRNA-seq data was split into training (80%) and testing (20%) datasets, for model development and evaluation. Four ML models were explored, namely support vector machine (SVM; Fig. 1d), mixture discriminant analysis (MDA; Fig. 1e), k-nearest neighbours (KNN) and generalised linear model (GLM) (Supplementary Table 1). The MDA confusion matrix was the most accurate and highlighted a marginally higher percentage of true positive DS cells than the SVM model (0.49% higher). Both the MDA and SVM models have an ROC score of 1, and a sensitivity of 0.999, but the MDA model had higher specificity than the SVM model (1 versus 0.99, respectively). KNN and GLM were also tested but had marginally lower evaluation metrics (KNN ROC: 0.998, Sens: 0.978, Spec: 0.992; and GLM ROC: 0.999, Sens: 0.996, Spec: 0.994, respectively). These evaluation metrics are displayed in Additional File 1: Table S1. For this reason, the MDA model was selected for subsequent work. Comparing the known classification of the testing dataset (Fig. 1f), with the MDA model predictions (Fig. 1g) shows a striking prediction accuracy of 99.79%. Although intra-experimental, this testing data is previously ‘unseen’, suggesting model overfitting is unlikely and instead, the two classification outcomes are transcriptionally distinct. This work demonstrates that it is possible to build an MDA ML model of senescence classification, which can accurately predict fibroblast deep senescence in a two-dimensional monolayer. Therefore, MDA machine learning models are used throughout this work.

To test the MDA ML model’s performance, we used an independent publicly available dataset from Chan et al. (GEO—#GSE175533↗), which included scRNA-seq data for increasing population doubling levels (PDLs) of WI-38 fibroblasts [16]. These fibroblasts ranged from PDL25 to PDL50, the latter of which Chan et al. described as being in an early senescent state. The transcriptional differences of the PDL50 senescent cells were apparent when the cells were clustered and plotted onto a UMAP (Fig. 2a), with the majority of the PDL46 and PDL50 cells clustering away from the main bulk of the population. Applying our MDA ML model to the WI-38 fibroblast dataset shows an enrichment for DS predictions in the UMAP locations of the early senescent PDL50 cells (Fig. 2b). However, only 10.08% of the PDL50 cells are predicted to be DS (Fig. 2c). To further investigate this, the PDL50 cells were isolated and clustered independently, producing eight clusters (Fig. 2d). Intriguingly, cluster five was largely enriched for DS-predicted cells, suggesting that this subpopulation of cells was transcriptionally similar to the DS cells in the original HDF model.

We hypothesised that the low percentage of DS prediction in the PDL50 cells could be due to temporal differences in the two experiments. The DS cells used to build the MDA ML model had undergone three weeks of culture after reaching replicative senescence, to establish and deepen their phenotype. Conversely, the PDL50 WI-38 fibroblast cells from Chan et al. were sequenced immediately upon reaching their Hayflick limit (Fig. 2e). To test this hypothesis, we applied Monocle3 trajectory analysis to the PDL50 cells [33]. Cluster five (the cluster enriched for cells which were predicted to be DS—Fig. 2f), fell at the end of the predicted trajectory, suggesting that the model was detecting deeply senescent but not early senescent fibroblasts (Fig. 2g) within the PDL50 population.

Cells at PDL50 likely represent a heterogeneous population, where some cells have reached the Hayflick limit earlier and deepened their senescent phenotype. To confirm this, we examined the expression of known cell cycle arrest markers CDKN1A (p21) and CDKN2A (p16/Arf). Alcorta et al. reported an initial increase of p21 during the early stages of cell cycle arrest of human diploid fibroblasts, with a gradual increase of p16 expression as the senescence phenotype deepens [35]. Intriguingly, this was also apparent in the PDL50 cells, with those cells predicted to be DS having lower expression levels of CDKN1A (Fig. 2h) and higher expression levels of CDKN2A mRNA (Fig. 2i) than the remaining PDL50 population. It must be noted that an increase of CDKN2A expression could be due to increased expression of either p16 or Arf alternative reading frames. However, using sequence-specific primers, the DS HDFs yielded a significant increase in p16 mRNA and a significant decrease in Arf compared to their proliferative counterparts with senescence (Additional File 1: Fig. S2). It is therefore likely that there is an increased expression of p16 but not Arf mRNA in the PDL50 DS predicted cells. This increase supports our hypothesis that the MDA ML model is able to detect deeply senescent but not early senescent fibroblasts.

To build a ML model which could encompass early replicative senescent fibroblasts (ES), the PDL50 cells from Chan et al. [16] were incorporated into the model as an additional ES classifier. First, clustering the ES, EP, and DS cells revealed a more successive pattern, with the ES cells bridging the gap between the EP and DS cells (Fig. 3a–b). This was confirmed using Monocle3 trajectory analysis (Fig. 3c). This is perhaps unsurprising, as we have previously identified the ES cells to contain a heterogeneous mix with a small proportion of DS cells. The MDA confusion matrix highlights more than 90.2% true positives for all intra-experimental predictions, with the lowest percentage of true positives being for the ‘early senescent’ or ES cells (Fig. 3d), which could again be explained by the heterogeneity of this population.

Applying this model, the WI-38 cells at PDL25 through to PDL50 were able to be classified as either EP, ES, or DS (Fig. 3e; original clustering in Fig. 2a). Unsurprisingly, the PDL50 cells had the highest predicted percentage of ES cells (66.6%). The percentage of predicted PDL50 cells then incrementally decreased with decreasing PDL number, with 15.7% of PDL25 cells being predicted to be ES (Fig. 3f). This suggests that even at the earlier passages, the population of fibroblasts is heterogeneous and that individual cells are on their own trajectories. Finally, for completeness, this ML predictive model for ES, EP, and DS cells was tested on the EP and DS HDF cells alone, identifying a small percentage (10.06%) of the EP cells as early senescent (Fig. 3g). However, within the DS cells, only 0.39% were predicted to be ES, indicating that the DS cells are a more homogenous population. Intriguingly, testing the ML model on two alternative senescence triggers (oncogene-induced senescence (OIS) and paracrine senescence) [10] suggests that the time spent in senescence is perhaps more important than the trigger (Additional File 1: Fig. S3a–c). The ES classifier was able to detect increasing senescence in the OIS and paracrine groups compared to the proliferating cells (27.9% and 32.9% more, respectively), but no DS cells were detected (Additional File 1: Fig. S3c).

To investigate whether the EP/ES/DS model can detect senescence in vivo, the models were tested on two independent publicly available whole skin scRNA-seq datasets from Tabib et al. [19] (data kindly provided by corresponding author) and Solé-Boldo et al. [20] (#GSE130973↗). Clustering the fibroblasts from Tabib et al. [19] revealed five clusters (Fig. 4a), which did not appear to stratify based on the age of the donors (Fig. 4b), or by their predicted senescence state (Fig. 4c). Analysis of the second whole skin scRNA-seq data from Solé-Boldo et al. [20] generated similar findings. Clustering these fibroblasts formed 11 distinct clusters (Fig. 4d), and these did appear to have some age-associated influences (Fig. 4e). Particularly, the 53-year-old donor clustered alone, which could suggest fibroblast abnormality for this individual, or perhaps issues with tissue isolation and processing. Intriguingly, this abnormality was also reflected in the percentage predictions of EP, ES, and DS cells, with the 53-year-old donor having a very high prediction of ES cells (73.9%) (Fig. 4f–g).

Combining the percentage of EP, ES and DS predictions from both the Tabib [19] and Solé-Boldo [20] datasets revealed a remarkably high prevalence of ES and DS predicted cells, regardless of donor age (Fig. 4g). The combination of DS and ES predicted cells makes up a minimum of 44.18% of the fibroblasts (24-year-old donor), and a maximum of 78.02% of the fibroblasts (53-year-old donor). Evidence within the literature suggests that senescent fibroblasts within human skin in vivo make up approximately 8–15% of the fibroblast population [36–39]. However, this literature relies on single markers which are not necessarily specific to senescence, or sufficient to detect all heterogeneous senescent populations, and making these predictions is confounded by a lack of ground truth in vivo [2]. Despite this, the large disparity between the model prediction and literature reports suggests that the model is not accurately able to identify senescent fibroblasts in vivo.

In recent years, the limitations of culturing cells in a two-dimensional monolayer have been widely reported. Most importantly, cells cultured in a monolayer do not replicate the complexities and depth of the extracellular matrix in vivo, which can have wide implications on cell growth and signalling processes [12–15, 40, 41]. Consequently, three-dimensional organotypic models have been developed in an effort to recapitulate in vivo cell behaviour more accurately. For this reason, we built a 3D fibroblast matrigel-collagen dermal gel containing either EP or DS HDFs, which underwent 10X scRNA-seq (Fig. 5a). To determine whether the transcriptome of fibroblasts grown in 3D or 2D was more similar to fibroblasts in vivo, a heatmap was constructed comparing the core matrisome (structural ECM proteins) between these groups (Additional File 1: Fig. S4) [42]. The core matrisome proteins measured are listed in Additional File 2; Supplementary Table S2. Hierarchical clustering reveals two initial branched clusters, where the fibroblasts grown in 2D cluster alone, and the fibroblasts grown in a 3D gel clustered alongside the in vivo fibroblasts from both Tabib et al. [19] and Solé-Boldo et al. [20], suggesting fibroblasts grown in a 3D gel are more similar to those in vivo based on their expression of ECM components.

Clustering of the fibroblasts grown in a 3D gel revealed six distinct groups, and a clear separation of EP and DS cells (Fig. 5b–c), highlighting that the transcriptomic differences between EP and DS cells were maintained in three-dimensional culture. Next, a MDA ML model was built to classify the cells from the dermal gels as EP or DS (called 3D SenPred), and intra-experimental testing revealed over 97.92% true positives (Fig. 5d). Promisingly, applying this ML model to the dataset from Tabib et al. [19] revealed a more physiological level of senescence prediction, with the maximum percentage of senescent fibroblasts being predicted in the 63-year-old (21.34%) (Fig. 5e). This was recapitulated when the 3D SenPred model was applied to the Solé-Boldo [20] dataset, with the highest predicted senescent cell burden being 32.79% in the 69-year-old donor (Fig. 5f).

Mine et al. reported that papillary fibroblasts change with age compared to reticular fibroblasts, including increasing secretion of keratinocyte growth factor, increased age-related contraction of collagen gels, and decreased population doublings and clonogenic capacity with age [43], suggesting that papillary fibroblasts are likely to be the fibroblast subtype to undergo replicative senescence in vivo. Therefore, we investigated the proportion of senescent cells within the four fibroblast subtypes defined by Solé-Boldo et al. [20]: secretory reticular, pro-inflammatory, secretory papillary, and mesenchymal. Promisingly, the secretory papillary fibroblasts from Solé-Boldo et al. [20] had the highest percentage of senescence prediction (40.25%), providing increased confidence that the ML model is detecting senescent cells in vivo (Additional File 1: Fig. S5).

Characterisation of senescence classically relies of the use of a selected panel of markers. Therefore, we were interested to explore how the expression of a panel of 9 commonly used senescence markers (CXCL12, CDKN1A (p21), CDKN2A (p16/Arf), IGF1, IGFBP3, IL-6, MMP3, SERPINE1, TP53) compared on a per cell basis. We started by using the scRNAseq data for the in vitro DS3D and EP3D cells because these have a state of “ground truth”. Using this data, we generated Pearson correlation coefficients for the comparison of normalised gene counts for each of the senescence markers, displaying the resulting correlation coefficients in a heatmap (Additional File: Fig. S6a–b). The highest positive correlation coefficient in the DS cells was when IGFBP3 and CXCL12 were compared (0.39), and this compared to 0.28 in the EP3D cells (Additional File: Fig. S6a–b). 1 1

In parallel, we generated Pearson correlation coefficients for the cells that the 3D SenPred model predicted as DS or EP within the Tabib [19] and Solé-Boldo [20] datasets (Fig. 5g and Additional File 1: Fig. S6C, respectively). Again, the highest correlation in the predicted DS cells was between IGFBP3 and CXCL12 (0.35), compared to 0.22 in the predicted EP cells. The Pearson correlation coefficient between CDKN1A and CDKN2A expression in the DS and EP predicted cells were 0.0046 and − 0.0022, respectively (individual correlation plots provided in Additional File 1: Fig. S6d–e), and between IGFBP3 and TP53 the correlation coefficient was 0.0225 and − 0.0126, respectively (individual correlation plots provided in Additional File 1: Fig. S6f–g) [22]. Further, 49.38% of gene pairs have a Pearson correlation coefficient of zero or less than zero (i.e. a negative correlation) in the DS predicted cells. Overall, this exercise revealed the absence of a strong correlation between pairs of commonly used senescence markers at the single cell level both in vivo and even in vitro where we have a state of ground truth, and emphasises how using preselected markers might only detect a subset of senescent cells. Taken together, this data further highlights the value of using SenPred, a holistic machine learning model built using thousands of genes, to detect the likelihood of a cell being senescent based on more nuanced changes, with the potential for improved identification of heterogeneous senescent cells whilst removing the likelihood of single-marker based false positives or false negatives.

Next, we expanded the in vivo datasets to incorporate a further 12 whole skin scRNA-seq dataset donors recently generated by Ganier et al. [22]. First, the fibroblasts were identified in line with Ganier et al. [22] (Additional File 1: Fig. S7). We then used the 3D SenPred model to predict the percentage of DS fibroblasts. Finally, we combined a total of 23 donor samples from the three independent in vivo datasets: Tabib [19], Solé-Boldo [20], and Ganier [22] to explore how the SenPred output changes relative to donor age (Fig. 5H). We observed that the percentage of DS predicted fibroblasts is significantly positively correlated with age (Pearson = 0.471, p value = 0.023), indicating that SenPred generated from 3D fibroblasts successfully predicts increasing percentages of deeply senescent fibroblast cells with age. Crucially, SenPred generated from 2D fibroblasts does not show any significant trend with age (Pearson = − 0.26, p = 0.23), further confirming the advantage of utilising data generated with cells grown in a 3D environment (Additional File 1: Fig. S8a). Unsurprisingly, the trend is also lost when the model is tested on whole skin (Additional File 1: Fig. S8b) or on fibroblasts isolated from lung tissue (Additional File 1: Fig. S8c) [24], highlighting the model’s tissue and cell type specificity in vivo.

In summary, the SenPred pipeline allows the generation of an unbiased predictive model of dermal fibroblast senescence. Using scRNA-seq data from cells with a known senescence classification, ML models can be generated and evaluated both intra- and inter-experimentally (Fig. 5i). In the future, these models could be generated in multiple tissues, trigger and cell-type-specific contexts of interest. The models can then be applied to classify senescence within unknown datasets. Importantly, using scRNA-seq data generated from cells grown in 3D allows the detection of senescent cells in vivo.

Not applicable.

Not applicable.

Authors AD, LJW and DAG are Unilever employees. The remaining authors declare that they have no competing interests.

Package dependencies

All packages and versions used in this code are available at the following URL: https://github.com/bethk-h/SenPred_HDF↗ [47].

Data Availability

The scRNAseq data generated as part of this work has been uploaded to GEO #GSE282425↗ [48].

Code Availability

All code for this project is accessible via the following URL: https://github.com/bethk-h/SenPred_HDF↗ [47].