Network-Community Analysis of Cellular Senescence

• The paper introduces an integrative network analysis using bulk RNA-seq to isolate a minimal, predictive gene set for cellular senescence.
• It employs eigenvector centrality and Louvain modularity to rank genes and reveal markers like Fabp3 and Cd59a with robust statistical validation.
• Validation via permutation tests and single-cell RNA-seq confirms the method’s cross-tissue relevance and practical implications for senescence research.

Cellular senescence is defined as a highly regulated, stable growth-arrest state, typically induced by diverse cellular insults and characterized by complex, multi-hallmark transcriptional remodeling. Identifying the minimal set of genes most predictive and mechanistically informative of senescence remains a significant challenge due to the genetic and phenotypic heterogeneity inherent to cell populations. Network-community analysis provides a framework to address this complexity through integrative feature selection and modular decomposition of gene expression data, enabling discovery of robust gene markers and their functional relationships within cellular senescence programs (Sabalic et al., 2024).

1. Experimental Design and Data Processing

The analytical workflow begins with low-input bulk RNA sequencing (RNA-seq) of sorted muscle-resident cell types—macrophages, fibro-adipogenic progenitors (FAPs), and satellite cells—under three cell states (basal, non-senescent injured, and senescent), sampled from both young and aged mice at days 3 and 7 after muscle injury. Each of the distinct 36 cell-condition combinations (3 cell types × 3 states × 2 ages × 2 time points) is profiled from three biological replicates, resulting in 108 libraries and expression quantification for 46,078 genes. After replicate averaging and stringent dropout filter application (removing genes with zero expression in at least two out of three replicates per condition), the retained dataset comprises 36 samples across 28,603 expressed genes.

State-specific gene selection is performed using the Kolmogorov–Smirnov test, Bonferroni-corrected for multiple comparisons; specifically, genes whose expression distributions in the senescent group significantly differ from the union of non-senescent and basal groups (p < 0.05) are retained, yielding 142 genes for subsequent network analysis.

2. Construction of Co-Expression Networks

The primary analytical object is a weighted, undirected network in which each node represents one of the 36 cell-condition samples. Edges are weighted by the Pearson correlation coefficient computed over the expression vectors of the 142 selected state-discriminative genes. Mathematically, for samples $i$ and $j$ with expression profiles $x_i, x_j \in \mathbb{R}^N$ ($N=142$ genes), the edge weight is

$$w_{ij} = \mathrm{corr}(x_i, x_j) = \frac{\sum_{k=1}^N (x_{i,k}-\bar{x}_i)\,(x_{j,k}-\bar{x}_j)} {\sqrt{\sum_{k}(x_{i,k}-\bar{x}_i)^2}\,\sqrt{\sum_{k}(x_{j,k}-\bar{x}_j)^2}}$$

No binarization or hard thresholding is performed; the network remains fully connected, with edge strengths distributed over $[-1, 1]$.

3. Eigenvector Centrality Feature Selection

To prioritize genes most relevant for distinguishing senescent from non-senescent and basal states, a supervised feature ranking based on eigenvector centrality is introduced. Each of the 142 genes forms a node within a fully connected gene-gene adjacency matrix $A$, constructed as the outer product of two supervised measures:

• Fisher score $f_u$ quantifies each gene’s mean expression separation between senescent and non-senescent/basal states, normalized by variance,

$$f_u = \frac{(\mu_{u,1}-\mu_{u,2})^2}{\sigma_{u,1}^2+\sigma_{u,2}^2}$$

• Mutual information $m_u$ quantifies the dependency of gene $u$’s expression on cell state labels.

The adjacency entry between genes $u$ and $v$ is $A_{uv} = f_u m_v$. The eigenproblem $A x = \lambda x$ is solved for the largest eigenvalue $\lambda$ and the associated non-negative eigenvector $x$. The centrality score $x_u$ for gene $u$ thus reflects both its individual discriminative power and its connectivity to other discriminative genes. Genes are ranked in descending order of $x_u$, producing an ordered candidate list $(G_1, G_2, \ldots, G_{142})$.

4. Community Detection and Community-State Alignment

For each prefix of the ranked gene list—the top $N$ genes for $N=1,2,\ldots,142$—a new co-expression network is generated over the 36 cell-conditions, rewiring edge weights as the Pearson correlation across the restricted gene subset. The Louvain modularity optimization algorithm is then applied to partition the network into communities, maximizing the modularity functional

$$Q = \frac{1}{2m} \sum_{i,j} \Big[A_{ij} - \frac{k_i k_j}{2m}\Big] \delta(c_i, c_j)$$

where $m$ is the total edge weight, $k_i$ is the degree of node $i$, $c_i$ the community assignment, and $\delta(c_i, c_j)$ the Kronecker delta.

To quantify the alignment between detected communities and actual biological cell states (basal, non-senescent, senescent), a tailored alignment score $h$ based on the reduction in Shannon entropy is used:

$$h = \begin{cases} \frac{H_{\max} - H_{\rm tot}}{H_{\max}}, & \text{if } n=3 \ 0, & \text{otherwise} \end{cases}$$

Here, $H_{\rm tot} = H_{\rm state} + H_{\rm community}$ sums the entropies of the contingency tables.

It is observed that maximal community-state separation $(h=0.737)$ occurs using $N=10$ genes. A permutation test with 1000 randomizations yields $p=0.0019$ for $h \geq 0.737$ under the null, confirming the specificity of the result.

Rank	Gene	Essential?
1	Fabp3	Yes
2	Gm38050	No
3	Lncpint	Yes
4	Pcnp	Yes
5	Luc7l3	Yes
6	Micalcl	No
7	Ash1l	No
8	Al427809	Yes
9	Cdhr4	No
10	Cd59a	Yes

Iterative leave-one-out analysis identifies six genes (Fabp3, Lncpint, Pcnp, Luc7l3, Al427809, Cd59a) as essential; their exclusion collapses the alignment score to zero, confirming their necessity for tripartite separation.

5. Functional Interpretation and Upstream Regulatory Analysis

Among these six essential genes, different roles and levels of validation are observed. Lncpint is a p53-induced long noncoding RNA functionally linked to DNA damage response, cell cycle arrest, and senescence. Fabp3 encodes a fatty-acid binding protein implicated in lipid metabolism and morphological remodeling associated with the senescence-associated secretory phenotype (SASP). The remaining four (Pcnp, Luc7l3, Al427809, Cd59a) are either newly identified or have less established roles in senescence.

Upstream transcription factor (TF) analysis (via GRNdb and KEGG pathways) shows that all four coding genes (Fabp3, Pcnp, Luc7l3, Cd59a) are regulated by TFs directly implicated in canonical senescence pathways, including cell cycle control (p53, E2F), inflammatory/SASP signaling (NF-κB), apoptosis resistance, DNA repair, and chromatin remodeling. This mapping supports the interpretation that the identified modules bridge classical senescence hallmarks.

6. Validation in Single-Cell RNA-Seq Data

Cross-tissue and cross-platform robustness of the network-identified markers is assessed using a published single-cell RNA-seq dataset from p16^INK4a reporter mice (Omori et al., 2020), profiling fibrotic liver and kidney. Ground truth senescent identity is inferred from tdTomato fluorescence (threshold 1.5) splitting cells into p16^+ (senescent) and p16^– populations.

Dimensionality reduction via UMAP reveals cell clustering largely consistent with p16 status. Overlay of Cd59a expression demonstrates spatially precise co-localization with p16^+ macrophages/Kupffer cells and, to a lesser degree, with endothelial and plasma cells. Quantitatively, Cd59a provides superior correspondence to the ground-truth p16^+ annotations compared to canonical markers, establishing cross-tissue generalizability and practical utility of the gene set derived from network-community analysis.

7. Biological and Methodological Implications

This integrative framework refines the transcriptional definition of cellular senescence to a small, robust, and transferable gene set mapping onto multiple hallmark pathways, including cell cycle arrest, lipid metabolism, SASP, and DNA damage response. Despite the utility of cell-type stratified bulk RNA-seq, the approach is limited in resolving intra-population heterogeneity, and certain separations are confounded by biological overlap (for example, non-senescent injured versus basal cells). The method’s reliance on Pearson correlation may also limit detection of nonlinear relationships, and results may be sensitive to selected statistical thresholds.

A plausible implication is that the general ECFS and network-community pipeline can be systematically extended to other complex phenotypes given bulk (or spatial) RNA-seq data with cell-type resolution and ground-truth groupings, such as developmental trajectories, cell fate decisions, or responses to perturbations (Sabalic et al., 2024).

---

Reference: "Network-community analysis of cellular senescence" (Sabalic et al., 2024).