SenSet identifies cell-type specific aging markers in the human lung

Several senescence gene sets have been published to date (GO:0090398↗, Fridman, SenMayo, and CellAge (Ashburner et al, 2000; Fridman and Tainsky, 2008; Saul et al, 2022; Avelar et al, 2020)). We first examined the extent of overlap between these gene sets. We observed that the pairwise overlap between them is relatively small compared with the size of each individual gene set, with the highest overlap of 34 genes shared between the GO and CellAge sets. The union U of all sets contains 501 unique marker genes, of which 434 were detected in the HLCA.

We sought to identify a subset of U that demonstrates greater sensitivity for senescence. The PUc estimator (Sakai and Shimizu, 2019) constructs a model of healthy cells based on data from (non-smoker) young and middle-aged individuals (groups \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathscr{Y}}}}$$\end{document}Y, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathscr{M}}}}$$\end{document}M, respectively) and applies this model to identify cells that are senescent in the aged group (group \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathscr{A}}}}$$\end{document}A, Figs. 1 and 2A). PUc accounts for potential covariate shifts in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathscr{A}}}}$$\end{document}A, which may arise due to other aging processes and hallmarks, including inflammation or epigenetic alterations (Li et al, 2023a; Wang et al, 2022a). The advantage of this approach is that it allows for a direct comparison of SnCs against non-SnCs within the same group, the oldest age group \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathscr{A}}}}$$\end{document}A. This addresses the challenge of age-related confounding factors that may arise when comparing older with younger individuals. While PUc identifies cells that generally deviate from the healthy (young) profile, we hypothesize that a significant proportion of these non-healthy cells identified by PUc are indeed senescent. We denote these cells with a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(-)$$\end{document}(−) superscript to signify that they belong to the negative (senescent) class. Similarly, a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(+)$$\end{document}(+) superscript will denote non-senescent cells for that group.

Cell types showed considerable heterogeneity in the expression of SenSet markers (Fig. 3E,F). From all SenSet genes, 87 were upregulated in alveolar macrophages⁽⁻⁾ and 84 in tracheobronchial smooth muscle cells (TSM)⁽⁻⁾, representing the highest numbers among the 22 cell types considered. The finding for TSM cells aligns with the GO analysis performed earlier. Conversely, basal⁽⁻⁾ cells and type 1 fibroblasts⁽⁻⁾ showed a downregulation of 43 and 34 genes, respectively.

Type II pneumocytes and fibroblasts are crucial structural cell types in the lung that have been implicated in senescence (Yao et al, 2021; Hayflick and Moorhead, 1961; Lin and Xu, 2020). In fibroblasts⁽⁻⁾, 19 SenSet genes were upregulated. Type II pneumocytes⁽⁻⁾ also show an upregulation of 19 SenSet genes (different set), and 1 downregulated gene, CTNNB1. We found substantial overlap in upregulated genes between basal⁽⁻⁾ cells and type II⁽⁻⁾ pneumocytes, fibroblasts⁽⁻⁾, respectively, with TNFRSF1A, CITED2, and ZFP36 in common across all three. For instance, 9 genes were upregulated in both fibroblasts⁽⁻⁾ and basal⁽⁻⁾ cells, and 9 genes were also upregulated in both type II pneumocytes⁽⁻⁾ and basal⁽⁻⁾ cells (Fig. 3G).

Basal cells represent bona fide stem cells of the lung, and stem cell exhaustion, recognized as a hallmark of aging, has been associated with senescence (Ruzankina and Brown, 2007). Among 34 downregulated SenSet genes in fibroblasts⁽⁻⁾, 19 of these were also downregulated in basal⁽⁻⁾ cells (Fig. 3H).

Several genes were found to be upregulated in one cell type and downregulated in the other. For instance, LMNA, encoding for lamin A protein, known to be downregulated in more immature/undifferentiated cells (Sehgal et al, 2013), was downregulated in basal⁽⁻⁾ cells, but upregulated in both fibroblasts⁽⁻⁾ and type II pneumocytes⁽⁻⁾.

To assess whether SenSet reliably captures senescence across diverse biological settings, we analyzed three independent single-cell and one bulk RNA-seq dataset.

To further refine and validate SenSet, we utilized a highly complex ex vivo human 3D tissue culture model based on precision-cut lung slices (PCLS). PCLS recapitulates the complexity of the lung environment in situ, enabling the study of various lung cell types and lineages present in the parenchymal region, along with the extracellular matrix (ECM) within the lungs’ native 3D architecture at high temporal and spatial resolution. We have previously shown that this model recapitulates the onset and progression of lung injury and disease and enables the testing of candidate therapeutics in viable ex vivo human tissue affected by disease (Uhl et al, 2015; Alsafadi et al, 2017, 2020).

Given that our SenSet signature was based on non-smokers, we next aimed to investigate whether SenSet enrichment is altered upon smoking. We performed DE testing on a few cell types of interest to evaluate if senescence markers were enriched in smokers for the youngest \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathscr{Y}}}}$$\end{document}Y and oldest age groups \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathscr{A}}}}$$\end{document}A. The analysis showed differences in senescence enrichment between smokers of varying ages (Fig. 8B). Specifically, we found that most markers were upregulated in basal cells from smokers (around 80% of the genes across all lists except for SenMayo at 50%). For other cell types, including CD4/CD8-positive, alpha-beta T cells and AT2 cells, we found that aged smokers showed an upregulation of a higher number of senescence markers than young smokers, across all marker lists, with SenSet containing most such genes. Lastly, for AT2 cells, around 35% of all genes were downregulated and 40% upregulated in the young, while around 60% were upregulated in the older age group.

The (core) Human Lung Cell Atlas (HLCA) (Sikkema et al, 2023) was downloaded from https://data.humancellatlas.org/hca-bio-networks/lung/atlases/lung-v1-0↗. Counts were already normalized. The HLCA harmonizes scRNA-seq data from 14 datasets, encompassing 106 individuals aged between 10 and 76 years. We removed one individual for whom the age was not available. While five levels of annotation are available in the data, we used the finest level, assigning one of 50 cell types to over 500,000 cells for the analysis in this study. The dataset also contains individuals with a smoking history, including 19 former and 28 active smokers. Smoking status is not available for 8 individuals, and these were not included in the analysis.

The protocols for preparing human PCLS and inducing senescence in PCLS through treatment with bleomycin and doxorubicin have been published on protocols.io: https://www.protocols.io/view/human-derived-precision-cut-lung-slices-hpcls-agar-36wgq3xr5lk5/v1↗https://www.protocols.io/view/senescence-induction-by-dna-damage-in-pcls-e6nvwdwpwlmk/v1↗. Lungs from healthy donors (ages 20–78) were collected, and the lower left lobes were inflated with 2.5% agarose in DMEM/F-12 with HEPES (Thermo Fisher Scientific Cat #12400024). After 45 min on ice, the lobes were sliced, and 1-cm diameter cores were extracted. Precision-cut lung slices (PCLS) were prepared using a Compresstome® to obtain 300-μm thick tissue slices with a diameter of 1 cm. The PCLS were then placed in a 24-well plate containing 1 mL of DMEM/F-12 medium supplemented with 1% FBS, 1% Penicillin–Streptomycin (Merck MilliporeSigma, Sigma-Aldrich Cat #P00781↗), and 0.3 μg/mL Amphotericin B solution (Merck MilliporeSigma, Sigma-Aldrich Cat #A2942) and incubated for 24 h at 37 °C with 5% CO₂ (day −1).

At 24 h (day 0), 72 h (day 2), and 120 h (day 4), the medium was replaced with fresh medium containing treatments diluted in DMEM/F-12 with 0.1% FBS, 1% Penicillin–Streptomycin, and 0.3 μg/mL Amphotericin B. The PCLS were treated under the following conditions: in PBS as the control, with 15 μg/mL bleomycin (Fresenius Kabi Cat #10361), with DMSO diluted at 1:100,000 in medium (Merck Millipore Sigme, Sigma-Aldrich, Cat #D2438), or with 0.1 μM doxorubicin hydrochloride (Merck MilliporeSigma, Sigma-Aldrich, Cat #D1515-10MG), originally dissolved in DMSO at 10 mM. After 168 h (day 6), the supernatants were collected and frozen at −80 °C for future multiplex immunoassay analysis by Luminex^TM. The PCLS were snap-frozen in liquid nitrogen for future protein extraction and snRNA-seq using the 10x Genomics^TM platform. In addition, PCLS were fixed for 30 min in 4% formaldehyde diluted in PBS (Life Technologies Cat #28908), washed in PBS, and embedded in paraffin. Separate PCLS were fixed for 30 min in the fixative solution from the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document}β-galactosidase staining kit (Cell Signaling Technology Cat #9860) and then washed in PBS.

Peritumor control tissue from three non-chronic lung diseases (N-CLD) patients and one COPD patient were obtained from the CPC-M bioArchive at the Comprehensive Pneumology Center (CPC Munich, Germany). Patients were two males and two females and had a mean age of 72.5 years old. Human lung tissue was filled with 3% of low gelling temperature agarose in DMEM/F-12 (Thermo Scientific, USA) with phenol red supplemented with 0.1% FCS, 1% P/S, and 1% amphotericin B and kept at 4 °C for at least 1 h. In total, 500 μm PCLS were generated using either a vibratome HyraxV50 (Zeiss, Germany) or 7000smz-2 Vibratome (Campden Instruments, England). The day after slicing, fresh medium was added, and PCLS were exposed to ionizing radiation using the RS225 X-ray cabinet (Xstrahl, Camberley, UK). Dose was calculated according to exposure time (30 Gray (Gy) = 12 min 24 sec) at 195 kV and 15 mA. Then, PCLS were kept in culture for up to 5 days at 37 °C, 5% CO₂, and medium was changed every 2–3 days.

All studies involving human participants were conducted in accordance with the ethical standards of the institutional research committee and with the 1964 Helsinki Declaration and its later amendments. Informed consent was obtained from all participants prior to inclusion in the study. The study also adhered to the principles outlined in the U.S. Department of Health and Human Services Belmont Report.

The study was approved by the local ethics committee of the Ludwig-Maximilians University of Munich, Germany (Ethic vote 19-630). Written informed consent was obtained for all study participants.

Paraffin-embedded PCLS were sliced at a thickness of 4 μm. One slide was subjected to H&E staining, and another slide was used for p21 and Ki67 immunostaining. The slides were rehydrated through a series of baths in xylene, followed by 100%, 95%, 85% ethanol, and finally water. The slides were then incubated at 105 °C for 20 min in 1× DAKO high pH buffer (Agilent Technologies, Dako Target Retrieval Solution pH 9 10×, Cat #S236784-2). Then the slides were washed in buffer A from Duolink® In Situ Wash Buffers (MilliporeSigma, Sigma-Aldrich Cat #DUO82049-20L), followed by incubation in 300 mM glycine for 30 min, and then in PBS containing 0.1% Tween and 0.5% Triton for 15 min.

The slides were then incubated with a blocking solution consisting of 2% BSA, 0.1% Triton, and 0.1% Tween at 37 °C for 45 min, followed by an overnight incubation at 4 °C with the primary antibody diluted in blocking solution (anti-p21 antibody [EPR362] Abcam, Cat #ab109520, 1:500 and anti-Ki67 antibody (8D5), Cell Signaling, Cat #9449, 1:600). After washing in buffer A, the slides were incubated for 1 h at room temperature with a 1:1000 dilution of the secondary antibody, anti-rabbit IgG 647; Biotium, Cat #20047 or anti-mouse IgG 750; Biotum, Cat #20463. Following further washes in buffer B from Duolink® In Situ Wash Buffers, the nuclei were stained with DAPI. The slides were then washed in buffer B and mounted with Fluoroshield Mounting Medium (Abcam, Cat #ab104135). Images were captured using an IX83 Olympus microscope, acquiring the entire PCLS area at 20x magnification. Images were analyzed using Olympus OlyVIA 3.1 (Olympus Soft Imaging Solutions GmbH).

Quantification of p21 and Ki67-positive cells was performed using Fiji (version 2.9.0). Automatic quantification of positive cell number identified by a nuclear p21 or Ki67 signal overlapping with the DAPI signal and normalized to total cell count by image based on the DAPI signal.

No blinding was performed during the experiments or outcome assessment.

Four PCLS were sonicated three times at 35% amplitude for 10 sec in 200 μL of buffer (TPER buffer, Thermo Scientific, Cat #78510, supplemented with Halt^TM Protease and Phosphatase Inhibitor Cocktail, EDTA-free (100X), Thermo Scientific, Cat #1861281, to a final concentration of 1×) and kept on ice. The samples were then homogenized for 15 sec using a Thermo Fisher Scientific homogenizer and centrifuged at 300 × g for 5 min at 4 °C. The supernatants were transferred to new 1.5-mL tubes and centrifuged at 10,000 rpm for 10 min at 4 °C. The supernatants were gently collected, and protein concentrations were quantified in triplicate using the Pierce detergent Compatible Bradford Assay kit (Thermo Scientific, Cat #23246), with 150 μL of reagent and 5 μL of sample per well. A standard curve was generated using the Prediluted Protein Assay Standards BSA Set (Thermo Scientific, Cat #23208). Absorbance was measured at 595 nm using SpectraMax ABS Plus Molecular Devices with SoftMaxPro 7.2.

Protein extracts were denatured in 1× Protein Loading Buffer (Li-COR, Cat #928-40004) containing 0.1 M dithiothreitol (DTT, Sigma-Aldrich, Cat #43816) for 10 min at 95 °C. A total of 10 μg of denatured protein was loaded onto Criterion TGX long shelf-life Precast Gels (4-15%, Bio-Rad, Cat #5671083) using 1× Tris/Glycine/SDS buffer (Bio-Rad, Cat #1610772). The proteins were then transferred onto an Odyssey Nitrocellulose Membrane (Li-COR, Cat #926-31092).

The membrane was then blocked for 1 h at room temperature in a 1:1 mixture of 1× TBS and Intercept Blocking Buffer (Li-COR, Cat #927-70001). The blocked membrane was incubated overnight at 4 °C with p21 antibody (Abcam, Cat #ab109520) diluted 1:1000 or GAPDH antibody (Abcam, Cat #ab9485) diluted 1:2500 in a 1:1 mixture of 1× TBST and Intercept Blocking Buffer.

The membrane was then washed three times for 10 min each with TBST and incubated with secondary antibodies: anti-rabbit red (IRDye® 680RD Donkey anti-Rabbit IgG, Li-COR, Cat #926-68073), diluted 1:20,000 in a 1:1 mixture of 1X TBST and Intercept Blocking Buffer for 1 h at room temperature.

The images were captured using Odyssey CLx Imager LICORbio and analyzed with Image Studio version 6.0 software.

PCLS supernatants were harvested after 6 days in culture and stored at −80 °C. The multiplex immunoassay was performed using the Luminex^TM platform, following the manufacturer’s instructions (Bio-Techne).

PCLS were snap-frozen in cryotubes using liquid nitrogen after 6 days of culture and stored in liquid nitrogen. Nuclei isolation from frozen PCLS follows the protocol at 10.17504/protocols.io.36wgqndb5gk5/v1. Four PCLS per experimental condition were then used for snRNA-seq, following the manufacturer’s instructions (10x Genomics^TM). Some minor individual differences in total-counts-per-cell were found across conditions (Fig. 6A).

Given the lack of labeled data, supervised algorithms such as SenCID (Tao et al, 2024) cannot be applied to classify SnCs in the HLCA dataset. This limitation necessitates the use of learning paradigms that do not rely on fully labeled training sets. One such approach is positive-unlabeled (PU) learning.

PU learning is a variant of binary classification where the goal is to distinguish between positive and negative samples, with the restriction that only positive samples are seen during the training phase (Bekker and Davis, 2020). More concretely, the real data distribution is a mixture given by\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p(x)=\alpha {p}_{+}(x)+(1-\alpha ){p}_{-}(x)$$\end{document}p(x)=αp+(x)+(1−α)p−(x)where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document}α is the mixture proportion, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${p}_{+},{p}_{-}$$\end{document}p+,p− are the probability density functions of the positive and negative samples, respectively. In traditional binary classification, we are given training data \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\mathscr{D}}}}}_{{{{\rm{tr}}}}}\sim p(x)$$\end{document}Dtr~p(x) and test data \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\mathscr{D}}}}}_{{{{\rm{te}}}}}\sim p(x)$$\end{document}Dte~p(x). However, in PU learning, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\mathscr{D}}}}}_{{{{\rm{tr}}}}}\sim {p}_{+}(x)$$\end{document}Dtr~p+(x) which makes it more challenging. Many methods have been proposed, which typically assume some form of smoothness or separability of the classes or treat the negative samples as noise (Scott, 2015; Garg et al, 2021; Kiryo et al, 2017).

PU learning is a natural fit for our task. Because SnCs cannot yet be labeled with certainty from single-cell transcriptomes, we treat the population of SnCs as negative samples within an unlabeled set. The objective is therefore to train a PU classifier that separates SnCs (negative) from healthy cells (positive). Subsequently, differential expression (DE) analysis between these two groups can be used to identify marker genes.

To build the labeled training set and the unlabeled test set, we rely on the age of the donors, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\mathscr{D}}}}}_{{{\rm{tr}}}}$$\end{document}Dtr models young individuals, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\mathscr{D}}}}}_{{{\rm{te}}}}$$\end{document}Dte models older ones. To accommodate the PU learning setup, we require the following assumptions to hold: (a) the young group contains few to no SnCs, (b) the negative (non-healthy) class in the senior group consists mostly of SnCs, and (c) healthy cells from older donors follow the same distribution as those from young donors. The last assumption is often violated, as older individuals suffer from other aging hallmarks not related to senescence, such as inflammation, epigenetic alterations, or mitochondrial dysfunction (Li et al, 2023a; Wang et al, 2022a; Amorim et al, 2022). Stated differently:\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${p}_{+}^{{{{\rm{young}}}}}\left(x\right)\,\ne\, {p}_{+}^{{{{\rm{senior}}}}}\left(x\right).$$\end{document}p+youngx≠p+seniorx.

Therefore, we face a covariate shift (Dharani. et al, 2019). In covariate shift, the dependence of the response variable \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y$$\end{document}y (i.e., senescence status) on gene expression \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x$$\end{document}x is the same for both training and test sets; however, the input distributions may not be:\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{array}{c}{p}_{{{{\rm{tr}}}}}(y| x)={p}_{{{{\rm{te}}}}}(y| x)\\ {p}_{{{{\rm{tr}}}}}(x)\,\ne\, {p}_{{{{\rm{te}}}}}(x).\end{array}$$\end{document}ptr(y∣x)=pte(y∣x)ptr(x)≠pte(x).

To address this, the PU learner needs to accommodate a covariate shift. Here, we rely on the formulation of Sakai & Shimizu which propose an unbiased risk estimator for covariate shift adaptation on PU learning, termed PUc (Sakai and Shimizu, 2019).

PUc assumes we are given three sets of samples: labeled training data \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{{x}_{i}^{{{\rm{Ptr}}}}\}_{i}\sim {p}_{{{\rm{tr}}}}(x{| }y=1)$$\end{document}{xiPtr}i~ptr(x∣y=1), unlabeled training data \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{{x}_{j}^{{{\rm{Utr}}}}\}_{j}\sim {p}_{{{\rm{tr}}}}(x)$$\end{document}{xjUtr}j~ptr(x), and unlabeled test data \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{{x}_{k}^{{{\rm{Ute}}}}\}_{k}\sim {p}_{{{\rm{te}}}}(x)$$\end{document}{xkUte}k~pte(x). When covariate shift occurs, the PU risk on the test distribution \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${p}_{{{\rm{te}}}}$$\end{document}pte differs from the PU risk on the train distribution \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${p}_{{{\rm{tr}}}}$$\end{document}ptr. PUc addresses this by importance-weighting, where the ratio between test and train densities \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$w(x):={p}_{{{{\rm{te}}}}}(x)/{p}_{{{{\rm{tr}}}}}(x)$$\end{document}w(x):=pte(x)/ptr(x) is used to weight each sample during the computation of the risk (Sugiyama et al, 2007, 2012). In this case, the PUc risk becomes\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${R}_{{{{\rm{PUc}}}}}\left(g\right):=\alpha {{\mathbb{E}}}_{x\sim {p}_{{{{\rm{tr}}}}}\left(x|y=1\right)}\left[\widetilde{{{{\mathscr{\ell }}}}}\left(g\left(x\right)\right)w\left(x\right)\right]+{{\mathbb{E}}}_{x\sim {p}_{{{{\rm{tr}}}}}\left(x\right)}\left[{{{\mathscr{\ell }}}}\left(-g\left(x\right)\right)w\left(x\right)\right],$$\end{document}RPUcg:=αEx~ptrx∣y=1\ell~gxwx+Ex~ptrx\ell−gxwx,where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g$$\end{document}g is a classifier and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathscr{\ell }}}}$$\end{document}\ell is a loss function with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widetilde{{{{\mathscr{\ell }}}}}(x):={{{\mathscr{\ell }}}}(x)-{{{\mathscr{\ell }}}}(-x)$$\end{document}\ell~(x):=\ell(x)−\ell(−x). The PUc risk on training data can be shown to be an unbiased estimator of the PU risk on test data. The mixture proportion \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document}α is estimated from prior knowledge. Similar to the original work, we employ a linear-in-parameter classifier with a Gaussian kernel basis function. Concretely, we represent the decision function\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g\left(x\right)=\mathop{\sum }_{{{{\mathscr{\ell }}}}=1}^{b}{\alpha }_{{{{\mathscr{\ell }}}}}{\phi }_{{{{\mathscr{\ell }}}}}\left(x\right)\,{{{\rm{with}}}}\,{\phi }_{{{{\mathscr{\ell }}}}}(x)=\exp (-{{||x}-{c}_{{{{\mathscr{\ell }}}}}{||}}^{2}/(2{\sigma }^{2})),$$\end{document}gx=∑\ell=1bα\ellϕ\ellxwithϕ\ell(x)=exp(−∣∣x−c\ell∣∣2/(2σ2)),where the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${c}_{{{{\mathscr{\ell }}}}}$$\end{document}c\ell are up to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b$$\end{document}b centers randomly sampled from the unlabeled test data and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma$$\end{document}σ is selected by cross-validation. For more details on PUc, please refer to Sakai and Shimizu (2019).

During the gene set generation step, we kept only individuals without a smoking history from the HLCA (~300,000 cells) to minimize potential confounding effects on the results. As described in the previous section, the PUc estimator requires three sets of samples: positive training samples, and unlabeled training and test samples. We estimated positive samples from individuals under the age of 30, assuming that the prevalence of SnCs in this group is minimal. The unlabeled training samples were estimated from individuals aged 30–50. Finally, the test samples were obtained from older individuals aged 50 and above, which is when covariate shift occurs. We call these three age groups Y, M, and A, respectively (Fig. 2A,E). Cell types with fewer than 50 cells in any age group were excluded from the analysis.

In PUc learning, inclusion of a middle-aged group containing both senescent and healthy cells is essential, as covariate shift between young and old samples makes it difficult to distinguish senescent from non-senescent cells in aged tissue using only these two extremes. The middle-aged group provides an intermediate anchor between training and test distributions, improving calibration and preventing conflation of gradual aging effects with true senescence.

To prepare the data for PUc learning, we first applied principal component analysis (PCA) independently for each cell type. We restricted the gene set used for PCA to the unionof all four existing senescence gene sets. By incorporating all known senescence-associated genes, we aim to achieve a “weak” separation of healthy cells and SnCs, which can be leveraged by the PUc learner. The top 10 components were used as training data for the PUc classifier. For all experiments, we set the mixture proportionto 0.9, based on the prior assumption that approximately 10% of the cells are senescent. However, the estimator was robust to this value and returned percentages in the range 0–40%. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\bf{U}}}}$$\end{document} U \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} α

Given the inherent randomness of the algorithms used, some variability across runs is expected. Nonetheless, we observed that ~80-90% of SenSet genes were consistently selected across runs with different random seeds.

Differential expression analysis is performed exclusively on the oldest age group, comparing healthy cells with SnCs. This approach is in contrast with methods that compare old and young individuals, where other aging signatures could introduce confounding factors. By directly comparing these two cell populations within the oldest age group, the analysis is specifically focused on senescence. A two-sided Wilcoxon rank-sum test was used to determine differentially expressed (DE) genes (FDR = 0.05). We tested only genes that belong to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\bf{U}}}}$$\end{document}U. FDR-adjusted P values were obtained using the Benjamini and Hochberg procedure (Benjamini and Hochberg, 1995). Cell types with fewer than 20 SnCs were not considered due to the limited sample size. We selected DE genes that were enriched in at least six cell types (either up- or downregulated), resulting in a set of 106 genes that constitute SenSet. A detailed table of all cell-specific SenSet genes and their test statistics is provided in the Appendix.

GSEA (Subramanian et al, 2005) was performed using the GSEApy package (Fang et al, 2023). To compute the Wasserstein distances in Fig. 8A, we sampled at random 5000 cells for cell types with too many cells to speed up computation. Scanpy was used in some parts of the analysis, including clustering and dimensionality reduction (Wolf et al, 2018).

For the cell-specific analysis, we first integrated the data using scVI (Lopez et al, 2018), focusing on the top 2000 variable genes. We used 2 hidden layers with 1000 nodes each. The dimensionality of the latent space was set to 30. Since we used raw counts for scVI, the gene likelihood was modeled as a zero-inflated negative binomial distribution. Nearest neighbors were computed in the scVI latent space, and clusters were identified via Leiden clustering (Traag et al, 2019). Clusters were manually annotated based on canonical markers of lung cell types (Travaglini et al, 2020). Clusters where no marker was significantly expressed were excluded from the cell-specific analysis. Wilcoxon rank-sum tests were used to assess whether SenSet genes were enriched in treatment samples compared to controls. For this analysis, bleomycin, doxorubicin, and irradiation cells were combined into one group.

For each gene set, we obtain a list of P values for each gene based on the DE test. In order to perform a meta-analysis, we combine these P values for each set using Pearson’s method, which emphasizes larger P values. I.e., given a set of P values \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{{p}_{i}\}_{i=1}^{n}$$\end{document}{pi}i=1n, Pearson’s method computes the statistic\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P:=-2\mathop{\sum }_{i=1}^{n}\log (1-{p}_{i}).$$\end{document}P:=−2∑i=1nlog(1−pi).

Under the null hypothesis \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${H}_{0\!\!}:{p}_{i}\sim U[0,1],i=[n]$$\end{document}H0:pi~U[0,1],i=[n], the test statistic \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P$$\end{document}P follows a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\chi }^{2}$$\end{document}χ2 distribution with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2n$$\end{document}2n degrees of freedom (Pearson, 1933). The combined P value provides an overall assessment of the enrichment of a gene set within a given cell type.

The scRNA-seq dataset of mesenchymal stem cells (MSCs) undergoing replicative senescence was downloaded from GEO (GSE200157↗). We used the two available replicas for this data (rep #2 and rep #3) at time points 1 (day 49) and time point 2 (day 86). For each sample, we filter out genes expressed in less than 50 cells and filter out cells with fewer than 100 expressed genes using Scanpy (Wolf et al, 2018). For each replica, only genes overlapping between the two time points are kept.

To rank genes for GSEA prerank, we compute the log1p fold change of average counts. I.e., if we denote the raw count vectors for genebyand, we compute the score foras \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g$$\end{document} g \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\bf{x}}}}}^{(T1)}\in {{\mathbb{R}}}^{{c}_{1}}$$\end{document} x ( ) T 1 R c 1 ∈ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{\bf{x}}}}}^{(T2)}\in {{\mathbb{R}}}^{{c}_{2}}$$\end{document} x ( ) T 2 R c 2 ∈ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g$$\end{document} g \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${g}_{{{{\rm{score}}}}}:=\frac{1}{{c}_{2}}\mathop{\sum }_{i=1}^{{c}_{2}}\log \left({{{{\bf{x}}}}}_{i}^{\left(T2\right)}+1\right)-\frac{1}{{c}_{1}}\mathop{\sum }_{i=1}^{{c}_{1}}\log \left({{{{\bf{x}}}}}_{i}^{\left(T1\right)}+1\right).$$\end{document} g score : = − . 1 c 2 1 c 1 ∑ i = 1 c 2 ∑ i = 1 c 1 log log x i T 2 + 1 x i T 1 + 1

GSEApy was then applied to all five gene sets to obtain an FDR q value.

The scRNA-seq IPF atlas was downloaded from GEO (GSE136831↗). Similar to the MSC data, we filter out genes expressed in fewer than 50 cells and cells with fewer than 100 expressed genes. This reduces the data to 147,169 IPF cells and 96,303 control cells, while leaving 33,435 genes. Also, like before, we compute gene scores for preranked GSEA based on the difference of average log1p counts. We perform preranked GSEA individually for each cell type, dropping types with fewer than 30 cells in either the healthy or fibrotic portion.

For the liver cirrhosis dataset, we also follow identical preprocessing steps to the IPF dataset. When performing preranked GSEA, we gather all cirrhotic samples and compare them against all control samples jointly.