Score Metabolic Pathways

Related to Figure 3.

Load Seurat objects

ren1c_obj <- readRDS("ren1c_obj.rds")
captopril_obj <- readRDS("captopril_obj.rds")
DefaultAssay(captopril_obj) <- "RNA"

ren1c_obj$cell_type <- ren1c_obj$reconciled_celltype
captopril_obj$cell_type <- captopril_obj$reconciled_celltype

ren1c_obj[["umap"]] <- ren1c_obj[["umap_harmony"]]

num_available_cores <- parallelly::availableCores()
register(MulticoreParam(workers = num_available_cores[[1]] - 4), default = TRUE)

Score metabolic pathways

calcMetabolicScores adds UCell-based pathway scores to @meta.data. UCell uses per-cell gene rankings so scores are robust to sequencing depth variation.

ren1c_obj     <- calcMetabolicScores(ren1c_obj)
captopril_obj <- calcMetabolicScores(captopril_obj)

Scores added: OXPHOS_Score, Glycolysis_Score, TCA_Score, FAO_Score, PPP_Score.

Statistical testing — Ren1^c KO model

# OXPHOS
ren1_oxphos_stats <- calcStatistics(
  ren1c_obj,
  group_var = "genotype",
  group1 = "WT",
  group2 = "KO",
  metric = "OXPHOS_Score"
)

# Glycolysis
ren1_glyc_stats <- calcStatistics(
  ren1c_obj,
  group_var = "genotype",
  group1 = "WT",
  group2 = "KO",
  metric = "Glycolysis_Score"
)

Statistical testing — captopril model

capt_oxphos_stats <- calcStatistics(
  captopril_obj,
  group_var = "condition",
  group1 = "ctrl",
  group2 = "trt",
  metric = "OXPHOS_Score"
)

capt_glyc_stats <- calcStatistics(
  captopril_obj,
  group_var = "condition",
  group1 = "ctrl",
  group2 = "trt",
  metric = "Glycolysis_Score"
)

Cohen's d effect sizes — main signatures

To compare effect magnitudes across pathways (rather than relying on p-values alone), Cohen's d is computed using the pooled-SD formula. Because the objects are enriched for renin-lineage cells, whole-object Cohen's d values are attenuated by cell-type composition. Interpret relative magnitudes across signatures, not against absolute benchmarks.

cohensD <- function(x, y) {
  nx <- length(x); ny <- length(y)
  pooled_sd <- sqrt(((nx - 1) * var(x) + (ny - 1) * var(y)) / (nx + ny - 2))
  (mean(y) - mean(x)) / pooled_sd
}

main_scores <- c("OXPHOS_Score", "Glycolysis_Score",
                 "TCA_Score", "FAO_Score",
                 "PPP_Score")

calcMainEffectSizes <- function(seurat_obj, group_var, group1, group2, model_label) {
  meta <- seurat_obj@meta.data
  g1   <- meta[meta[[group_var]] == group1, ]
  g2   <- meta[meta[[group_var]] == group2, ]
  scores_present <- main_scores[main_scores %in% colnames(meta)]
  results <- lapply(scores_present, function(s) {
    d <- cohensD(g1[[s]], g2[[s]])
    data.frame(
      model = model_label, score = s,
      group1 = group1, group2 = group2,
      mean_group1 = mean(g1[[s]]), mean_group2 = mean(g2[[s]]),
      cohens_d = d,
      magnitude = dplyr::case_when(
        abs(d) >= 0.8 ~ "large",
        abs(d) >= 0.5 ~ "medium",
        abs(d) >= 0.2 ~ "small",
        TRUE          ~ "negligible"
      )
    )
  })
  do.call(rbind, results)
}

main_effect_capt  <- calcMainEffectSizes(captopril_obj, "condition", "ctrl", "trt", "captopril")
main_effect_ren1c <- calcMainEffectSizes(ren1c_obj,     "genotype",  "WT",   "KO",  "ren1c")
main_effect_all   <- rbind(main_effect_capt, main_effect_ren1c)

fwrite(main_effect_all,
       paste0(parent_dir, "main_metabolic_signatures_cohens_d.csv"),
       row.names = FALSE, sep = ",")

Captopril model summary (ctrl → trt, positive d = elevated in treated):

Glycolysis shows the largest transcriptional shift, consistent with HIF-driven metabolic reprogramming.

K-means metabolic zone clustering

Produces manuscript Figures 3g and 3m (cells clustered by OXPHOS, glycolysis, and balance scores into k=4 zones ranked by mean glycolysis score; KO and captopril-treated cells are enriched in the highest glycolytic cluster).

calcMetabolicClustering and plotClustering are defined in metabolic_analysis_functions.R. Features used: OXPHOS_Score, Glycolysis_Score, Metabolic_Balance (OXPHOS − Glycolysis), and Metabolic_Ratio (OXPHOS / Glycolysis), scaled before clustering. Clusters are ranked 1–4 by mean glycolysis (rank 1 = highest). Enrichment of KO/treated cells toward higher-glycolysis clusters is tested with a Wilcoxon test on per-cell Glycolysis_Rank values; per-cluster Fisher's exact tests are BH-corrected.

source("metabolic_analysis_functions.R")

ren1_clustering <- calcMetabolicClustering(
    ren1c_obj,
    group_var    = "genotype",
    model_name   = "Ren1 KO",
    group_levels = c("WT", "KO"),
    method       = "kmeans",
    k            = 4
)

capt_clustering <- calcMetabolicClustering(
    captopril_obj,
    group_var    = "condition",
    model_name   = "Captopril",
    group_levels = c("ctrl", "trt"),
    method       = "kmeans",
    k            = 4
)

ren1_cluster_plots <- plotClustering(ren1c_obj,     ren1_clustering, "genotype",  "Ren1 KO")
capt_cluster_plots <- plotClustering(captopril_obj, capt_clustering, "condition", "Captopril")

for (plt_name in names(ren1_cluster_plots)) {
    ggsave(paste0(parent_dir, "ren1c_metabolic-clustering_", plt_name, "_", Sys.Date(), ".svg"),
           plot = ren1_cluster_plots[[plt_name]], units = "mm", dpi = 300)
    ggsave(paste0(parent_dir, "captopril_metabolic-clustering_", plt_name, "_", Sys.Date(), ".svg"),
           plot = capt_cluster_plots[[plt_name]], units = "mm", dpi = 300)
}