11  VEXAS Syndrome

12 Clonal Dominance in VEXAS Syndrome

In this manuscript(Molteni et al. 2025)the researchers examined the mechanisms driving clonal dominance in VEXAS syndrome through:

  • Immunophenotypic characterization of hematopoiesis in nine male patients
  • Single-cell transcriptomics analysis
  • Development of humanized models using base editing to insert the causative UBA1 mutation in healthy hematopoietic stem and progenitor cells (HSPCs)

Key findings: * VEXAS patients show pervasive inflammation across all hematopoietic lineages * Patient HSPCs are skewed toward myelopoiesis and acquire senescence-like programs * Humanized models successfully recapitulated disease hallmarks including proteostatic defects and inflammatory signatures * Competitive transplantations revealed that UBA1-mutant cells are more resilient to inflammatory environments while wild-type cells become exhausted * This dynamic leads to clonal dominance, impaired hematopoiesis, and eventually bone marrow failure * The mechanism involves “poisoning” of wild-type HSPCs by inflammatory factors produced by UBA1-mutant progeny

The study provides novel insight into the pathogenesis of VEXAS syndrome and establishes models for preclinical testing of potential therapeutic strategies targeting either inflammatory signaling or the resilience mechanisms of mutant clones.

```{r} # VEXAS Syndrome Data Analysis # Based on: “Mechanisms of hematopoietic clonal dominance in VEXAS syndrome” (Nature Medicine, 2025)

13 Load required libraries

library(ggplot2) library(dplyr) library(tidyr) library(viridis) library(pheatmap) library(pheatmap) library(gridExtra) # Note: ComplexHeatmap was removed as it requires installation

14 Set a consistent theme for plots

theme_set(theme_minimal() + theme(text = element_text(size = 12), plot.title = element_text(size = 14, face = “bold”), legend.position = “bottom”))

15 —- 1. VEXAS Patient Cohort Characteristics —-

16 Create patient data based on Extended Data Table 1

patient_data <- data.frame( PatientID = paste0(“PT”, 1:9), Age = c(71, 71, 75, 66, 72, 74, 70, 71, 84), Mutation = c(“Met41Val”, “Met41Thr”, “Met41Leu”, “Met41Thr”, “Met41Val”, “Splice”, “Met41Leu”, “Met41Thr”, “Met41Val”), MDS_at_diagnosis = c(TRUE, TRUE, FALSE, TRUE, FALSE, TRUE, FALSE, TRUE, FALSE), Steroid_treatment = c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE), Vacuolization = c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE) )

17 Mutation distribution

mutation_plot <- ggplot(patient_data, aes(x = Mutation, fill = Mutation)) + geom_bar() + scale_fill_viridis_d() + labs(title = “Distribution of UBA1 Mutations in VEXAS Patients”, y = “Count”, x = “Mutation”) + theme(legend.position = “none”)

18 Clinical feature distribution

clinical_features <- patient_data %>% select(MDS_at_diagnosis, Steroid_treatment, Vacuolization) %>% pivot_longer(cols = everything(), names_to = “Feature”, values_to = “Present”) %>% group_by(Feature, Present) %>% summarise(Count = n(), .groups = “drop”) %>% filter(Present == TRUE)

clinical_features\(Feature <- factor(clinical_features\)Feature, levels = c(“MDS_at_diagnosis”, “Steroid_treatment”, “Vacuolization”), labels = c(“MDS at diagnosis”, “Steroid treatment”, “Vacuolization”))

clinical_plot <- ggplot(clinical_features, aes(x = Feature, y = Count, fill = Feature)) + geom_col() + scale_fill_viridis_d() + labs(title = “Clinical Features in VEXAS Patients”, y = “Number of Patients”, x = ““) + theme(legend.position =”none”) + geom_text(aes(label = Count), vjust = -0.5)

19 —- 2. Simulate Hematopoietic Cell Population Data —-

20 Create simulated data for hematopoietic cell populations based on Fig. 1

set.seed(123)

21 Cell types and their expected distributions in controls vs VEXAS

cell_types <- c(“HSC”, “MPP”, “MLP”, “GMP”, “CMP”, “MEP”, “PreBNK”, “ETPs”) control_mean <- c(15, 25, 20, 5, 8, 12, 10, 5) vexas_mean <- c(8, 15, 5, 15, 18, 8, 3, 28) control_sd <- rep(3, length(cell_types)) vexas_sd <- rep(4, length(cell_types))

22 Generate data for 7 controls and 19 VEXAS patients (as mentioned in the paper)

control_data <- lapply(1:7, function(i) { data.frame( Sample = paste0(“Control_”, i), CellType = cell_types, Count = pmax(0, rnorm(length(cell_types), mean = control_mean, sd = control_sd)), Group = “Control” ) }) %>% bind_rows()

vexas_data <- lapply(1:19, function(i) { data.frame( Sample = paste0(“VEXAS_”, i), CellType = cell_types, Count = pmax(0, rnorm(length(cell_types), mean = vexas_mean, sd = vexas_sd)), Group = “VEXAS” ) }) %>% bind_rows()

cell_population_data <- bind_rows(control_data, vexas_data)

23 Calculate summary statistics

cell_summary <- cell_population_data %>% group_by(Group, CellType) %>% summarise( Mean = mean(Count), SD = sd(Count), .groups = “drop” )

24 Plot cell population differences

cell_plot <- ggplot(cell_summary, aes(x = CellType, y = Mean, fill = Group)) + geom_bar(stat = “identity”, position = position_dodge()) + geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), position = position_dodge(width = 0.9), width = 0.25) + scale_fill_manual(values = c(“Control” = “#3498db”, “VEXAS” = “#e74c3c”)) + labs(title = “HSPC Population Distribution in VEXAS vs. Controls”, x = “Cell Type”, y = “Cell Count”) + theme(axis.text.x = element_text(angle = 45, hjust = 1))

25 Calculate myeloid:lymphoid ratio (based on Fig. 1e)

myeloid_cells <- c(“GMP”, “CMP”, “ETPs”) lymphoid_cells <- c(“MLP”, “PreBNK”)

ratio_data <- cell_population_data %>% mutate(LineageType = case_when( CellType %in% myeloid_cells ~ “Myeloid”, CellType %in% lymphoid_cells ~ “Lymphoid”, TRUE ~ “Other” )) %>% filter(LineageType %in% c(“Myeloid”, “Lymphoid”)) %>% group_by(Sample, Group, LineageType) %>% summarise(TotalCount = sum(Count), .groups = “drop”) %>% pivot_wider(names_from = LineageType, values_from = TotalCount) %>% mutate(Ratio = Myeloid / Lymphoid)

ratio_plot <- ggplot(ratio_data, aes(x = Group, y = Ratio, fill = Group)) + geom_boxplot(outlier.shape = NA) + geom_jitter(width = 0.2, alpha = 0.5) + scale_fill_manual(values = c(“Control” = “#3498db”, “VEXAS” = “#e74c3c”)) + labs(title = “Myeloid:Lymphoid HSPC Ratio”, y = “Ratio”, x = ““) + theme(legend.position =”none”)

26 —- 3. Simulate UBA1 Mutation Distribution Across Cell Types —-

27 Create simulated data for UBA1 mutation distribution based on Fig. 1h

cell_lineages <- c(“CD34+”, “CD13+”, “CD3+”, “CD19+”, “CD56+”) mutation_mean <- c(85, 80, 15, 10, 45) mutation_sd <- c(10, 15, 8, 5, 20)

28 Generate data for patients

mutation_data <- data.frame( CellType = rep(cell_lineages, each = 7), Patient = rep(paste0(“PT”, 1:7), length(cell_lineages)), PercentMutation = pmin(100, pmax(0, rnorm(length(cell_lineages) * 7, rep(mutation_mean, each = 7), rep(mutation_sd, each = 7)))) )

29 Plot UBA1 mutation distribution

mutation_plot_cells <- ggplot(mutation_data, aes(x = CellType, y = PercentMutation, color = Patient)) + geom_boxplot(outlier.shape = NA, fill = “lightgray”, alpha = 0.3) + geom_jitter(width = 0.2) + scale_color_viridis_d() + labs(title = “UBA1 Mutation Distribution Across Cell Lineages”, x = “Cell Type”, y = “Percentage of UBA1 Mutant Cells”) + theme(axis.text.x = element_text(angle = 45, hjust = 1))

30 —- 4. Simulate Gene Expression Data Based on Fig. 4 —-

31 Create simulated gene expression data for hallmark pathways in UBA1mut vs UBA1wt

hallmark_pathways <- c( “IFNα response”, “IFNγ response”, “Inflammatory response”, “TNF signaling via NF-κB”, “Complement”, “Apoptosis”, “Heme metabolism”, “p53 pathway”, “Oxidative phosphorylation”, “MYC targets v1”, “E2F targets”, “G2M checkpoint” )

32 Simulated NES scores for different cell types

cell_clusters <- c(“HSC/MPP”, “GMP”, “Monocytes”, “B cells”, “Erythroid”)

33 Create a matrix of NES scores

set.seed(123) nes_matrix <- matrix( c( # HSC/MPP rnorm(6, 3, 0.5), # Inflammatory pathways positive rnorm(3, 1, 0.5), # Metabolism pathways slightly positive rnorm(3, -2, 0.5), # Cell cycle pathways negative

# GMP
rnorm(6, 4, 0.5),  # Inflammatory pathways strongly positive
rnorm(3, 2, 0.5),  # Metabolism pathways positive
rnorm(3, -3, 0.5), # Cell cycle pathways strongly negative

# Monocytes
rnorm(6, 5, 0.5),  # Inflammatory pathways very strongly positive
rnorm(3, 1, 0.5),  # Metabolism pathways slightly positive
rnorm(3, -2, 0.5), # Cell cycle pathways negative

# B cells
rnorm(6, 2, 0.5),  # Inflammatory pathways moderately positive
rnorm(3, 0.5, 0.5), # Metabolism pathways neutral/slightly positive
rnorm(3, -1, 0.5),  # Cell cycle pathways slightly negative

# Erythroid
rnorm(6, 2.5, 0.5), # Inflammatory pathways positive
rnorm(3, 3, 0.5),   # Metabolism pathways strongly positive (esp. heme)
rnorm(3, -1.5, 0.5) # Cell cycle pathways negative

), nrow = length(hallmark_pathways), ncol = length(cell_clusters) )

rownames(nes_matrix) <- hallmark_pathways colnames(nes_matrix) <- cell_clusters

34 Convert to a data frame for plotting

nes_data <- as.data.frame(nes_matrix) %>% mutate(Pathway = hallmark_pathways) %>% pivot_longer(cols = -Pathway, names_to = “CellType”, values_to = “NES”)

35 Assign pathway categories for coloring

nes_data <- nes_data %>% mutate(Category = case_when( Pathway %in% c(“IFNα response”, “IFNγ response”, “Inflammatory response”, “TNF signaling via NF-κB”, “Complement”, “Apoptosis”) ~ “Inflammatory”, Pathway %in% c(“Heme metabolism”, “p53 pathway”, “Oxidative phosphorylation”) ~ “Metabolic”, Pathway %in% c(“MYC targets v1”, “E2F targets”, “G2M checkpoint”) ~ “Cell Cycle”, TRUE ~ “Other” ))

36 —- 5. Create a summary of the competitive transplantation experiment (Fig. 6) —-

37 Generate data for UBA1mut percentages at different input ratios

transplant_conditions <- c(“0% UBA1mut”, “25% UBA1mut”, “40% UBA1mut”, “100% UBA1mut”) cell_types_transplant <- c(“HSPCs”, “Myeloid”, “Lymphoid”)

38 Expected UBA1mut percentages in output for each condition and cell type

expected_percentages <- matrix( c( 0, 0, 0, # 0% input 25, 45, 5, # 25% input
80, 85, 15, # 40% input 100, 100, 100 # 100% input ), nrow = length(transplant_conditions), ncol = length(cell_types_transplant), byrow = TRUE )

39 Create data frame with simulated measurements (5 mice per group)

transplant_data <- data.frame() for (i in 1:length(transplant_conditions)) { for (j in 1:length(cell_types_transplant)) { for (mouse in 1:5) { transplant_data <- rbind(transplant_data, data.frame( Condition = transplant_conditions[i], InputPercentage = c(0, 25, 40, 100)[i], CellType = cell_types_transplant[j], Mouse = paste0(“Mouse_”, mouse), OutputPercentage = pmin(100, pmax(0, rnorm(1, expected_percentages[i, j], 8))) )) } } }

40 Calculate means and SD for plotting

transplant_summary <- transplant_data %>% group_by(Condition, InputPercentage, CellType) %>% summarise( Mean = mean(OutputPercentage), SD = sd(OutputPercentage), .groups = “drop” )

41 Plot the competitive transplantation results

transplant_plot <- ggplot(transplant_summary, aes(x = factor(InputPercentage), y = Mean, fill = CellType)) + geom_bar(stat = “identity”, position = position_dodge()) + geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), position = position_dodge(width = 0.9), width = 0.25) + scale_fill_viridis_d() + labs(title = “UBA1mut Percentage in Output Cells After Competitive Transplantation”, x = “Input Percentage of UBA1mut HSPCs”, y = “Output Percentage of UBA1mut Cells”) + geom_hline(yintercept = c(0, 25, 40, 100), linetype = “dashed”, color = “darkgray”) + theme(legend.position = “bottom”)

42 —- 6. Create a heatmap visualization of inflammatory cytokines (Fig. 6p) —-

43 Simulated cytokine concentrations from Fig. 6p

cytokines <- c(“IL1RA”, “IL1β”, “IL18”, “CXCL10”, “IL8”, “IFNγ”, “CCL4”, “TNFα”, “IFNα”) conditions <- c(“UBA1mut”, “UBA1wt”)

44 Log-scaled concentrations

set.seed(123) cytokine_values <- matrix( c( rlnorm(3, 8, 0.5), # IL1RA (high in mutant) rlnorm(3, 6, 0.5), # IL1β (high in mutant) rlnorm(3, 5.5, 0.5), # IL18 (high in mutant) rlnorm(3, 5, 0.5), # CXCL10 (high in mutant) rlnorm(3, 4, 0.5), # IL8 (moderate in mutant) rlnorm(3, 3, 0.5), # IFNγ (moderate in mutant) rlnorm(3, 3.5, 0.5), # CCL4 (moderate in mutant) rlnorm(3, 2, 0.5), # TNFα (low in mutant) rlnorm(3, 1, 0.5), # IFNα (low in mutant)

rlnorm(4, 5, 0.5),    # IL1RA (lower in WT)
rlnorm(4, 2, 0.5),    # IL1β (lower in WT)
rlnorm(4, 2.5, 0.5),  # IL18 (lower in WT)
rlnorm(4, 1.5, 0.5),  # CXCL10 (lower in WT)
rlnorm(4, 3, 0.5),    # IL8 (similar in WT)
rlnorm(4, 2, 0.5),    # IFNγ (lower in WT)
rlnorm(4, 2.5, 0.5),  # CCL4 (lower in WT)
rlnorm(4, 1.5, 0.5),  # TNFα (lower in WT)
rlnorm(4, 1, 0.5)     # IFNα (similar in WT)

), ncol = 9 )

cytokine_data <- data.frame( Condition = rep(conditions, c(3, 4)), Mouse = c(paste0(“Mut_”, 1:3), paste0(“WT_”, 1:4)) )

for (i in 1:length(cytokines)) { cytokine_data[[cytokines[i]]] <- cytokine_values[, i] }

45 Convert to long format for plotting

cytokine_long <- cytokine_data %>% pivot_longer(cols = -c(Condition, Mouse), names_to = “Cytokine”, values_to = “Concentration”) %>% group_by(Condition, Cytokine) %>% summarise( Mean = mean(Concentration), SD = sd(Concentration), .groups = “drop” )

46 Calculate fold changes

cytokine_wide <- cytokine_long %>% pivot_wider(id_cols = Cytokine, names_from = Condition, values_from = Mean) %>% mutate(FoldChange = UBA1mut / UBA1wt, LogFC = log2(FoldChange))

47 Cytokine bar plot

cytokine_plot <- ggplot(cytokine_long, aes(x = Cytokine, y = Mean, fill = Condition)) + geom_bar(stat = “identity”, position = position_dodge()) + geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), position = position_dodge(width = 0.9), width = 0.25) + scale_fill_manual(values = c(“UBA1mut” = “#e74c3c”, “UBA1wt” = “#3498db”)) + labs(title = “Inflammatory Cytokine Levels in VEXAS Model”, x = ““, y =”Concentration (pg/ml)“) + theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position =”bottom”) + scale_y_log10()

48 —- 7. Visualize the model of VEXAS clonal dominance —-

49 Create schematic data to represent the model in Extended Data Fig. 7

set.seed(123)

50 Time points for simulation

time_points <- 0:100

51 Initial cell counts

initial_wt <- 100 initial_mut <- 10

52 Simulation parameters

wt_growth_rate <- 0.05 mut_growth_rate <- 0.04 inflammation_effect_on_wt <- 0.08 inflammation_effect_on_mut <- 0.02 inflammation_rate <- 0.03

53 Simulate cell populations over time

simulation_data <- data.frame( Time = time_points, WildType = numeric(length(time_points)), Mutant = numeric(length(time_points)), Inflammation = numeric(length(time_points)) )

simulation_data\(WildType[1] <- initial_wt simulation_data\)Mutant[1] <- initial_mut simulation_data$Inflammation[1] <- 10

for (i in 2:length(time_points)) { # Update inflammation level - increases with mutant cell number simulation_data\(Inflammation[i] <- simulation_data\)Inflammation[i-1] + inflammation_rate * simulation_data$Mutant[i-1] / 100

# Wild type cells grow but are more sensitive to inflammation simulation_data\(WildType[i] <- simulation_data\)WildType[i-1] * (1 + wt_growth_rate - inflammation_effect_on_wt * simulation_data$Inflammation[i-1] / 100)

# Mutant cells grow slower but are less sensitive to inflammation simulation_data\(Mutant[i] <- simulation_data\)Mutant[i-1] * (1 + mut_growth_rate - inflammation_effect_on_mut * simulation_data$Inflammation[i-1] / 100)

# Safety checks to prevent negative values simulation_data\(WildType[i] <- max(simulation_data\)WildType[i], 0) simulation_data\(Mutant[i] <- max(simulation_data\)Mutant[i], 0) }

54 Convert to long format for plotting

simulation_long <- simulation_data %>% select(Time, WildType, Mutant) %>% pivot_longer(cols = c(WildType, Mutant), names_to = “CellType”, values_to = “Count”)

55 Calculate percentages

simulation_percentages <- simulation_data %>% mutate( Total = WildType + Mutant, WildTypePercent = 100 * WildType / Total, MutantPercent = 100 * Mutant / Total ) %>% select(Time, WildTypePercent, MutantPercent, Inflammation) %>% pivot_longer(cols = c(WildTypePercent, MutantPercent), names_to = “CellType”, values_to = “Percentage”)

56 Plot the simulation

model_counts_plot <- ggplot(simulation_long, aes(x = Time, y = Count, color = CellType)) + geom_line(size = 1.2) + scale_color_manual(values = c(“WildType” = “#3498db”, “Mutant” = “#e74c3c”), labels = c(“Wild-type cells”, “UBA1mut cells”)) + labs(title = “VEXAS Clonal Dominance Model: Absolute Cell Counts”, x = “Time”, y = “Cell Count”) + theme(legend.title = element_blank())

model_percent_plot <- ggplot(simulation_percentages %>% filter(CellType %in% c(“MutantPercent”)), aes(x = Time, y = Percentage)) + geom_line(size = 1.2, color = “#e74c3c”) + geom_area(alpha = 0.3, fill = “#e74c3c”) + labs(title = “VEXAS Clonal Dominance: UBA1mut Percentage”, x = “Time”, y = “Percentage of UBA1mut cells”) + scale_y_continuous(limits = c(0, 100))

inflammation_plot <- ggplot(simulation_data, aes(x = Time, y = Inflammation)) + geom_line(size = 1.2, color = “#9b59b6”) + geom_area(alpha = 0.3, fill = “#9b59b6”) + labs(title = “Inflammatory Milieu in VEXAS Model”, x = “Time”, y = “Inflammation Level”) + scale_y_continuous(limits = c(0, max(simulation_data$Inflammation)))

57 —- 8. Create an integrated figure showing the key findings —-

59 VEXAS Treatment Options

59.1 AML-directed therapies in VEXAS patients.