Metabolic Gene Sets and Functions
Related to Figure 3.
Five curated gene sets cover the major metabolic axes relevant to vascular smooth muscle cell (VSMC) and renin-lineage cell biology: OXPHOS, Glycolysis, TCA cycle, fatty acid oxidation (FAO), and the pentose phosphate pathway (PPP). MSigDB Hallmark metabolic pathways are loaded alongside for cross-reference. calcMetabolicScores and calcStatistics are defined here and called from both the main scoring analysis (5B) and the FAO decomposition (5C).
OXPHOS
The OXPHOS gene set covers all five electron transport chain complexes, mitochondrially-encoded subunits, and accessory factors involved in mitochondrial biogenesis, import, and membrane integrity:
oxphos_genes <- c(
# Complex I (NADH dehydrogenase)
"Ndufa1", "Ndufa2", "Ndufa3", "Ndufa4", "Ndufa5", "Ndufa6", "Ndufa7",
"Ndufa8", "Ndufa9", "Ndufa10", "Ndufa11", "Ndufa12", "Ndufa13",
"Ndufb1", "Ndufb2", "Ndufb3", "Ndufb4", "Ndufb5", "Ndufb6", "Ndufb7",
"Ndufb8", "Ndufb9", "Ndufb10", "Ndufb11",
"Ndufs1", "Ndufs2", "Ndufs3", "Ndufs4", "Ndufs5", "Ndufs6", "Ndufs7", "Ndufs8",
"Ndufc1", "Ndufc2",
"Ndufv1", "Ndufv2", "Ndufv3",
# Complex II (Succinate dehydrogenase)
"Sdha", "Sdhb", "Sdhc", "Sdhd",
# Complex III (Cytochrome bc1)
"Uqcrc1", "Uqcrc2", "Uqcrfs1", "Uqcrb", "Uqcrh", "Uqcrq",
"Cyc1", "Uqcr10", "Uqcr11",
# Complex IV (Cytochrome c oxidase)
"Cox4i1", "Cox4i2", "Cox5a", "Cox5b", "Cox6a1", "Cox6a2", "Cox6b1", "Cox6b2",
"Cox6c", "Cox7a1", "Cox7a2", "Cox7b", "Cox7c", "Cox8a",
"Cox10", "Cox11", "Cox15", "Cox17", "Cox7a2l", "Cox6cp3", "Cox7b2", "Cox8c",
# Complex V (ATP synthase)
"Atp5f1a", "Atp5f1b", "Atp5f1c", "Atp5f1d", "Atp5f1e",
"Atp5mc1", "Atp5mc2", "Atp5mc3",
"Atp5me", "Atp5mf", "Atp5mg", "Atp5mj", "Atp5mk", "Atp5ml", "Atp5mo", "Atp5mp",
"Atp5pb", "Atp5pd", "Atp5pf", "Atp5po",
"Atp5mc1p5",
# Mitochondrially-encoded subunits
"mt-Atp6", "mt-Atp8", "mt-Co1", "mt-Co2", "mt-Co3", "mt-Cyb",
"mt-Nd1", "mt-Nd2", "mt-Nd3", "mt-Nd4", "mt-Nd4l", "mt-Nd5", "mt-Nd6",
# Accessory and regulatory
"Phb2", "Fdx1", "Nqo2", "Gpx4", "Etfa", "Etfb", "Etfdh",
"Tomm70a", "Hccs", "Vdac2", "Vdac1", "Vdac3", "Dld", "Bckdha", "Nnt",
"Casp7", "Abcb7", "Atp6ap1", "Fxn", "Ogdh", "Cpt1a", "Acaa2", "Idh3g",
"Sucla2", "Atp6v0e", "Hsd17b10", "Grpel1", "Timm17a", "Mrps12", "Polr2f",
"Aifm1", "Tcirg1", "Pdk4", "Pdhx", "Mrpl15", "Timm8b", "Timm10", "Timm13",
"Timm9", "Slc25a4", "Slc25a5", "Slc25a3", "Ech1", "Mtx2", "Atp6v0b", "Mtrr",
"Atp6v1f", "Mrpl35", "Atp6v1g1", "Atp6v1c1", "Iscu", "Mrps15", "Mgst3",
"Timm50", "Acadsb", "Retsat", "Decr1", "Atp6v1h", "Aldh6a1", "Idh3a",
"Slc25a11", "Mpc1", "Mrps11", "Pdhb", "Oxa1l", "Afg3l2", "Bdh2", "Ndufab1",
"Pmpca", "Lrpprc", "Atp6v1d", "Opa1", "Immt", "Rhot1", "Slc25a12", "Dlst",
"Mrps30", "Cyb5a", "Suclg1", "Mrps22", "Htra2", "Slc25a20", "Hadha", "Hadhb",
"Echs1", "Mrpl11", "Mrpl34", "Acaa1a", "Idh3b", "Mtrf1", "Rhot2", "Dlat",
"Supv3l1", "Mfn2", "Tomm22", "Pdp1", "AK157302", "Gm10053", "Acadm", "Acat1",
"Aco2", "Alas1", "Prdx3", "Atp1b1", "Atp6v0c", "Cs", "Phyh", "Atp6v1e1",
"Acadvl", "Eci1", "Cyb5r3", "Fh1", "Glud1", "Got2", "Gpi1", "Hspa9",
"Idh1", "Idh2", "Ldha", "Ldhb", "Maob", "Mdh2", "Mdh1", "Oat", "Pdha1",
"Por", "Surf1", "Bax", "Atp12a", "Atp4a", "Atp4b",
"Atp6v0a1", "Atp6v0a2", "Atp6v0a4", "Atp6v0d1", "Atp6v0d2", "Atp6v0e1",
"Atp6v0e2", "Atp6v1a", "Atp6v1b1", "Atp6v1b2", "Atp6v1c2", "Atp6v1e2",
"Atp6v1g2", "Atp6v1g3", "Lhpp", "Ndufa4l2", "Ppa1", "Ppa2",
"Sdhcp5", "Uqcr10p1", "Uqcrhl"
)
Glycolysis
Core glycolytic enzymes plus HIF-responsive glucose transporters, lactate dehydrogenases, pyruvate dehydrogenase regulators, and glycogen metabolism genes:
glycolysis_genes <- c(
# Core glycolysis
"Hk1", "Hk2", "Hk3",
"Gpi1",
"Pfkl", "Pfkm", "Pfkp",
"Aldoa", "Aldob", "Aldoc",
"Tpi1",
"Gapdh", "Gapdhs",
"Pgk1", "Pgk2",
"Pgam1", "Pgam2",
"Eno1", "Eno2", "Eno3",
"Pkm", "Pklr",
"Ldha", "Ldhb", "Ldhc",
# Glycogen metabolism
"Gys1", "Gys2",
"Pygb", "Pygl",
"Agl",
# Extended glycolysis-related genes
"Vcan", "Vegfa", "Dsc2", "Lct", "Tff3", "Hmmr", "Ier3", "Gpc4", "Gpc3",
"Gclc", "Il13ra1", "G6pdx", "Chpf", "Efna3", "Dld", "Pfkfb1", "Ext2",
"Aldh7a1", "Kif2a", "Stc1", "Irs2", "Cxcr4", "Rragd", "Ankzf1", "Nsdhl",
"Capn5", "Elf3", "Gpc1", "Kif20a", "Bik", "Taldo1", "Tpst1", "Cited2",
"Casp6", "Slc37a4", "Lhx9", "Stc2", "Ppp2cb", "Cldn3", "Qsox1", "Artn",
"Bpnt1", "Cth", "Tpbg", "Rbck1", "Med24", "Hs2st1", "Psmc4", "Hax1",
"Plod2", "Gnpda1", "Homer1", "Sdc1", "Sdc3", "Sdc2", "Chst4", "Dpysl4",
"Zfp292", "Slc25a10", "Ero1a", "Slc25a13", "Gne", "Hs6st2", "Copb2",
"Nasp", "Gal3st1", "B4galt2", "Fut8", "Pmm2", "Ak3", "Aldh9a1", "Angptl4",
"Paxip1", "Ndufv3", "Chst2", "Arpp19", "Miox", "Cldn9", "Sdhc", "Pgls",
"Rpe", "Polr3k", "Rars1", "Nanp", "Pkp2", "Ppfia4", "Gmppa", "Galk2",
"Fam162a", "Chpf2", "Slc35a3", "Ecd", "B3gat3", "Abcb6", "Gale", "Ddit4",
"Depdc1a", "Cog2", "Lhpp", "B3gat1", "Chst1", "Nol3", "Cyb5a", "B4galt4",
"Cacna1h", "Isg20", "Chst12", "Sap30", "Akr1a1", "Srd5a3", "Chst5", "Egln3",
"Ndst3", "Tktl1", "Slc16a3", "Gpr87", "Glrx", "Glce", "Kdelr3", "Me2",
"B3gnt3", "B3galt6", "Cln6", "Ugp2", "Pdk3", "Alg1", "B4galt7", "Spag4",
"Xylt2", "Gmppb", "Eno1b", "Agrn", "Ak4", "Ang", "Cd44", "Cdk1", "Cenpa",
"Col5a1", "P4ha2", "Ext1", "Aurka", "Dcn", "Egfr", "Fbp2", "Fkbp4", "Gfpt1",
"B4galt1", "Galk1", "Got1", "Got2", "Hspa5", "Gusb", "Idh1", "Idua",
"Igfbp3", "Stmn1", "Mertk", "Met", "Mif", "Me1", "Mdh2", "Mdh1", "Mpi",
"Mxi1", "P4ha1", "Pam", "Pcx", "Pgm2", "Phka2", "Ppia", "Prps1", "Sod1",
"Sox9", "Tgfa", "Gfus", "Txn1", "Vldlr", "Hdlbp", "Adora2b", "Nt5e",
"Plod1", "Tgfbi", "Gck", "Gpi", "Pfkfb2", "Pfkfb3", "Pfkfb4", "Ppp2r5d",
"Arnt", "Entpd5", "Gpd1", "Hif1a", "Htr2a", "Igf1", "Ins", "Insr", "Myc",
"P2rx7", "Prkaa1", "Prkaa2", "Mlxipl", "Ldhal6b", "Fbp1", "G6pc1", "Pck1",
"Pdha1", "Pdha2", "Slc2a1", "Slc2a2", "Slc2a3", "Slc2a4", "Dlat", "Pdhx",
"Mpc1", "Slc2a5", "Pdhb", "Mpc2"
)
TCA cycle
Core Krebs cycle enzymes plus assembly factors and regulatory proteins (37 genes):
tca_cycle_genes <- c(
"Cs", "Aco1", "Aco2",
"Idh1", "Idh2", "Idh3a", "Idh3b", "Idh3g",
"Ogdh", "Dlst", "Dld",
"Suclg1", "Suclg2", "Sucla2",
"Sdha", "Sdhb", "Sdhc", "Sdhd",
"Fh1", "Mdh1", "Mdh2",
"Trap1", "Isca2", "Sdhaf2", "Sirt3", "Iscu", "Sdhaf4",
"Nfs1", "Acat1", "Sdhaf3", "Isca1", "Lyrm4", "Fxn",
"Mrps36", "Sdhaf1", "Nnt"
)
FAO
Fatty acid transport, mitochondrial β-oxidation enzymes, and peroxisomal acyl-CoA oxidases (27 genes):
fatty_acid_oxidation_genes <- c(
# Fatty acid transport
"Slc27a1", "Slc27a2", "Slc27a3", "Slc27a4", "Slc27a5",
"Cd36",
"Cpt1a", "Cpt1b", "Cpt1c",
"Cpt2",
"Slc25a20",
# Mitochondrial beta-oxidation
"Acadl", "Acadm", "Acads", "Acadvl",
"Echs1",
"Hadha", "Hadhb",
"Acaa1", "Acaa2",
"Eci1", "Eci2",
"Decr1",
# Peroxisomal
"Acox1", "Acox2", "Acox3"
)
PPP
Pentose phosphate pathway, oxidative and non-oxidative phases (16 genes):
pentose_phosphate_genes <- c(
"G6pdx", "G6pd2", "Pgls", "Pgd",
"Rpe", "Rpia", "Tkt", "Taldo1",
"Shpk", "Prps2", "Rbks", "Pgm2", "Dera", "Prps1", "Prps1l3", "Prps1l1"
)
Hallmark gene sets
Hallmark metabolic pathways are loaded from MSigDB for cross-reference with the curated scores. Mouse gene symbols are requested directly (db_species = "MM"):
hallmark_all <- msigdbr(species = "Mus musculus", db_species = "MM")
metabolism_hallmarks <- c(
"HALLMARK_OXIDATIVE_PHOSPHORYLATION",
"HALLMARK_GLYCOLYSIS",
"HALLMARK_FATTY_ACID_METABOLISM",
"HALLMARK_ADIPOGENESIS",
"HALLMARK_CHOLESTEROL_HOMEOSTASIS",
"HALLMARK_XENOBIOTIC_METABOLISM",
"HALLMARK_BILE_ACID_METABOLISM",
"HALLMARK_HEME_METABOLISM",
"HALLMARK_REACTIVE_OXYGEN_SPECIES_PATHWAY"
)
hallmark_gene_lists <- list()
for (pathway in metabolism_hallmarks) {
hallmark_gene_lists[[pathway]] <- hallmark_all %>%
dplyr::filter(gs_name == pathway) %>%
dplyr::pull(gene_symbol)
}
calcMetabolicScores
Scores all five curated pathways plus nine Hallmark sets with AddModuleScore_UCell, applies KNN smoothing over the UMAP neighborhood, and derives Metabolic_Ratio (OXPHOS/Glycolysis) and Metabolic_Balance (OXPHOS − Glycolysis).
calcMetabolicScores <- function(seurat_obj) {
seurat_obj <- AddModuleScore_UCell(
obj = seurat_obj, features = list(oxphos_genes),
name = "OXPHOS_Score", assay = "RNA", storeRanks = TRUE
)
seurat_obj <- AddModuleScore_UCell(
obj = seurat_obj, features = list(glycolysis_genes),
name = "Glycolysis_Score", assay = "RNA"
)
seurat_obj <- AddModuleScore_UCell(
obj = seurat_obj, features = list(tca_cycle_genes),
name = "TCA_Score", assay = "RNA"
)
seurat_obj <- AddModuleScore_UCell(
obj = seurat_obj, features = list(fatty_acid_oxidation_genes),
name = "FAO_Score", assay = "RNA"
)
seurat_obj <- AddModuleScore_UCell(
obj = seurat_obj, features = list(pentose_phosphate_genes),
name = "PPP_Score", assay = "RNA"
)
for (pathway_name in names(hallmark_gene_lists)) {
seurat_obj <- AddModuleScore_UCell(
obj = seurat_obj,
features = list(hallmark_gene_lists[[pathway_name]]),
name = pathway_name, assay = "RNA"
)
}
# UCell prepends "signature_1" to the name argument; rename to expected column names
colnames(seurat_obj@meta.data)[colnames(seurat_obj@meta.data) == "signature_1OXPHOS_Score"] <- "OXPHOS_Score"
colnames(seurat_obj@meta.data)[colnames(seurat_obj@meta.data) == "signature_1Glycolysis_Score"] <- "Glycolysis_Score"
colnames(seurat_obj@meta.data)[colnames(seurat_obj@meta.data) == "signature_1TCA_Score"] <- "TCA_Score"
colnames(seurat_obj@meta.data)[colnames(seurat_obj@meta.data) == "signature_1FAO_Score"] <- "FAO_Score"
colnames(seurat_obj@meta.data)[colnames(seurat_obj@meta.data) == "signature_1PPP_Score"] <- "PPP_Score"
for (pathway_name in names(hallmark_gene_lists)) {
clean_name <- gsub("signature_1", "", pathway_name)
old_col <- paste0("signature_1", clean_name)
if (old_col %in% colnames(seurat_obj@meta.data))
colnames(seurat_obj@meta.data)[colnames(seurat_obj@meta.data) == old_col] <- clean_name
}
# KNN smoothing across UMAP neighbors
score_names <- colnames(seurat_obj@meta.data)[
grep("_Score|^HALLMARK_", colnames(seurat_obj@meta.data))
]
k <- ceiling(sqrt(ncol(seurat_obj)))
seurat_obj <- SmoothKNN(seurat_obj, reduction = "umap",
signature.names = score_names,
k = k, suffix = paste0("_kNN", k), up.only = TRUE)
smoothed_cols <- colnames(seurat_obj@meta.data)[grep("_kNN", colnames(seurat_obj@meta.data))]
# Rename HALLMARK columns to title-case to avoid name collision with smoothed copies
colnames(seurat_obj@meta.data)[colnames(seurat_obj@meta.data) == "HALLMARK_OXIDATIVE_PHOSPHORYLATION"] <- tools::toTitleCase(tolower("HALLMARK_OXIDATIVE_PHOSPHORYLATION"))
colnames(seurat_obj@meta.data)[colnames(seurat_obj@meta.data) == "HALLMARK_GLYCOLYSIS"] <- tools::toTitleCase(tolower("HALLMARK_GLYCOLYSIS"))
colnames(seurat_obj@meta.data)[colnames(seurat_obj@meta.data) == "HALLMARK_FATTY_ACID_METABOLISM"] <- tools::toTitleCase(tolower("HALLMARK_FATTY_ACID_METABOLISM"))
colnames(seurat_obj@meta.data)[colnames(seurat_obj@meta.data) == "HALLMARK_ADIPOGENESIS"] <- tools::toTitleCase(tolower("HALLMARK_ADIPOGENESIS"))
colnames(seurat_obj@meta.data)[colnames(seurat_obj@meta.data) == "HALLMARK_CHOLESTEROL_HOMEOSTASIS"] <- tools::toTitleCase(tolower("HALLMARK_CHOLESTEROL_HOMEOSTASIS"))
colnames(seurat_obj@meta.data)[colnames(seurat_obj@meta.data) == "HALLMARK_XENOBIOTIC_METABOLISM"] <- tools::toTitleCase(tolower("HALLMARK_XENOBIOTIC_METABOLISM"))
colnames(seurat_obj@meta.data)[colnames(seurat_obj@meta.data) == "HALLMARK_BILE_ACID_METABOLISM"] <- tools::toTitleCase(tolower("HALLMARK_BILE_ACID_METABOLISM"))
colnames(seurat_obj@meta.data)[colnames(seurat_obj@meta.data) == "HALLMARK_HEME_METABOLISM"] <- tools::toTitleCase(tolower("HALLMARK_HEME_METABOLISM"))
colnames(seurat_obj@meta.data)[colnames(seurat_obj@meta.data) == "HALLMARK_REACTIVE_OXYGEN_SPECIES_PATHWAY"] <- tools::toTitleCase(tolower("HALLMARK_REACTIVE_OXYGEN_SPECIES_PATHWAY"))
# Rename smoothed columns: "_Score_kNNk" -> "_Score"; HALLMARK strips kNN suffix
for (col in smoothed_cols) {
new_name <- if (grepl("^HALLMARK_", col)) {
gsub(paste0("_kNN", k, "$"), "", col)
} else {
gsub(paste0("_Score_kNN", k, "$"), "_Score", col)
}
colnames(seurat_obj@meta.data)[grep(col, colnames(seurat_obj@meta.data))] <- new_name
}
seurat_obj$Metabolic_Ratio <- seurat_obj$OXPHOS_Score / seurat_obj$Glycolysis_Score
seurat_obj$Metabolic_Ratio[is.infinite(seurat_obj$Metabolic_Ratio)] <- NA
seurat_obj$Metabolic_Balance <- seurat_obj$OXPHOS_Score - seurat_obj$Glycolysis_Score
return(seurat_obj)
}
calcStatistics
Runs a Wilcoxon rank-sum test comparing a single metric column between group1 and group2 cells, both overall and stratified by cell type. Cell-type comparisons require at least five cells per group; p-values are BH-adjusted:
calcStatistics <- function(seurat_obj, group_var, group1, group2, metric = "OXPHOS_Score") {
data <- seurat_obj@meta.data
overall_data <- data %>% dplyr::filter(get(group_var) %in% c(group1, group2))
wilcox_overall <- wilcox.test(
overall_data[[metric]][overall_data[[group_var]] == group1],
overall_data[[metric]][overall_data[[group_var]] == group2],
alternative = "two.sided"
)
cell_types <- unique(data$cell_type)
celltype_results <- list()
for (ct in cell_types) {
ct_data <- data %>%
dplyr::filter(cell_type == ct, get(group_var) %in% c(group1, group2))
group1_vals <- ct_data[[metric]][ct_data[[group_var]] == group1]
group2_vals <- ct_data[[metric]][ct_data[[group_var]] == group2]
if (length(group1_vals) > 5 && length(group2_vals) > 5) {
test_result <- wilcox.test(group1_vals, group2_vals, alternative = "two.sided")
celltype_results[[ct]] <- data.frame(
cell_type = ct,
n_group1 = length(group1_vals),
n_group2 = length(group2_vals),
median_group1 = median(group1_vals, na.rm = TRUE),
median_group2 = median(group2_vals, na.rm = TRUE),
mean_group1 = mean(group1_vals, na.rm = TRUE),
mean_group2 = mean(group2_vals, na.rm = TRUE),
p_value = test_result$p.value,
statistic = test_result$statistic
)
}
}
celltype_df <- if (length(celltype_results) > 0) {
df <- dplyr::bind_rows(celltype_results)
df$p_adj <- p.adjust(df$p_value, method = "BH")
df$significant <- df$p_adj < 0.05
df
} else NULL
list(
overall_pvalue = wilcox_overall$p.value,
overall_statistic = wilcox_overall$statistic,
celltype_results = celltype_df
)
}