Step 11: Reactive Astrocyte Trajectory and lncRNA Discovery

Modeling astrocyte reactivity in PD using Slingshot pseudotime

Setup: Environment and Data

In this step, we shift from neuronal trajectory analysis to reactive astrocyte trajectory analysis.

The previous neuron-based analysis showed that the combined PD damage score was strongly cluster-dependent and did not cleanly separate PD from control. Therefore, in this module, we focus on astrocyte populations, especially the cluster annotated as Reactive Astrocyte.

The biological goal is to model a transition from:

Astrocyte
→ reactive astrocyte
→ stressed / inflammatory astrocyte state

This trajectory should be interpreted as an astrocyte reactivity trajectory, not a neuronal degeneration trajectory.

library(Seurat)
library(tidyverse)
library(patchwork)
library(SingleCellExperiment)
library(slingshot)
library(tradeSeq)

# Load validated annotated object
import_path <- "/Users/yoshimurasouhei/Downloads/010_school/4年生/bioinfomaticsリサーチクラークシップ/PD_2026/scripts/SO_08_Validated.RDS"

SO <- readRDS(import_path)

SO

An object of class Seurat 
47773 features across 43448 samples within 1 assay 
Active assay: RNA (47773 features, 2000 variable features)
 3 layers present: counts, data, scale.data
 2 dimensional reductions calculated: pca, umap

Step 11.1: Check Current Cell Type Annotation

First, we confirm the cell type identities stored in the Seurat object.

Step 11.1: Check Current Cell Type Annotation

First, we confirm the cell type identities stored in the Seurat object.

# Show only the summary, not every single cell identity
table(Idents(SO))

            Unassigned              Astrocyte        Oligodendrocyte 
                 14949                   7927                   4209 
         Stressed Cell Microglia / Macrophage                   OPCs 
                  3352                   3248                   2703 
      Endothelial Cell      Inhibitory Neuron  Fibroblast / Pericyte 
                  1777                   1057                   1045 
    Reactive Astrocyte       Unassigned Glial                 Neuron 
                  1025                    596                    992 
     Excitatory Neuron       T Cell / NK Cell             Fibroblast 
                   356                    124                     88 

colnames(SO@meta.data)

 [1] "orig.ident"        "nCount_RNA"        "nFeature_RNA"     
 [4] "Sample"            "condition"         "Sex"              
 [7] "scDblFinder.score" "scDblFinder.class" "percent.mt"       
[10] "merge_group"       "S.Score"           "G2M.Score"        
[13] "Phase"             "old.ident"         "RNA_snn_res.0.5"  
[16] "seurat_clusters"  

# Optional: run manually only if you want to inspect every cell identity.
# This output is very long, so it is disabled during rendering.
Idents(SO)

DimPlot(
  SO,
  reduction = "umap",
  label = TRUE,
  repel = TRUE,
  pt.size = 0.4
) +
  NoLegend() +
  labs(title = "Annotated Cell Types")

Step 11.2: Extract Astrocyte and Reactive Astrocyte Populations

We include both baseline astrocytes and reactive astrocytes so that Slingshot can infer a possible reactivity trajectory.

target_astrocytes <- c(
  "Astrocyte",
  "Reactive Astrocyte"
)

SO_astro <- subset(SO, idents = target_astrocytes)

SO_astro

An object of class Seurat 
47773 features across 8952 samples within 1 assay 
Active assay: RNA (47773 features, 2000 variable features)
 3 layers present: counts, data, scale.data
 2 dimensional reductions calculated: pca, umap

p1 <- DimPlot(
  SO_astro,
  reduction = "umap",
  label = TRUE,
  repel = TRUE,
  pt.size = 0.5
) +
  labs(title = "Astrocyte and Reactive Astrocyte Subset")

p2 <- DimPlot(
  SO_astro,
  reduction = "umap",
  group.by = "condition",
  pt.size = 0.5,
  cols = c("C" = "dodgerblue", "PD" = "firebrick")
) +
  labs(title = "Condition Distribution")

p1 | p2

table(Idents(SO_astro), SO_astro$condition)

                    
                        C   PD
  Astrocyte          5273 2654
  Reactive Astrocyte   12 1013

Step 11.3: Re-cluster Astrocytes

Because astrocytes are heterogeneous, we re-run variable feature selection, PCA, UMAP, and clustering within the astrocyte subset.

SO_astro <- FindVariableFeatures(
  SO_astro,
  selection.method = "vst",
  nfeatures = 2000
)

SO_astro <- ScaleData(SO_astro)
SO_astro <- RunPCA(SO_astro, npcs = 50, verbose = FALSE)

ElbowPlot(SO_astro, ndims = 50) +
  labs(title = "Astrocyte Elbow Plot")

SO_astro <- FindNeighbors(SO_astro, dims = 1:20, verbose = FALSE)
SO_astro <- FindClusters(SO_astro, resolution = 0.3, verbose = FALSE)
SO_astro <- RunUMAP(SO_astro, dims = 1:20, verbose = FALSE)

p1 <- DimPlot(
  SO_astro,
  reduction = "umap",
  group.by = "seurat_clusters",
  label = TRUE,
  repel = TRUE,
  pt.size = 0.5
) +
  labs(title = "Astrocyte Subclusters")

p2 <- DimPlot(
  SO_astro,
  reduction = "umap",
  group.by = "condition",
  pt.size = 0.5,
  cols = c("C" = "dodgerblue", "PD" = "firebrick")
) +
  labs(title = "Condition Distribution in Astrocytes")

p1 | p2

Step 11.4: Validate Astrocyte and Reactive Astrocyte Markers

We check canonical astrocyte markers and reactive astrocyte markers.

Astrocyte markers:

AQP4, GFAP, ALDH1L1, SLC1A3, S100B

Reactive astrocyte / inflammatory markers:

CHI3L1, TNC, C3, SERPINA3, LCN2, HSPB1, VIM

FeaturePlot(
  SO_astro,
  features = c(
    "AQP4", "GFAP", "ALDH1L1", "SLC1A3", "S100B",
    "CHI3L1", "TNC", "C3", "SERPINA3", "LCN2", "HSPB1", "VIM"
  ),
  ncol = 4,
  pt.size = 0.5
)

Warning: The following requested variables were not found: LCN2

Step 11.5: Define Astrocyte Reactivity Modules

Instead of using a broad neuronal PD damage score, we define astrocyte-specific reactivity scores.

astro_identity_genes <- c(
  "AQP4", "GFAP", "ALDH1L1", "SLC1A3", "S100B"
)

reactive_astrocyte_genes <- c(
  "CHI3L1", "TNC", "C3", "SERPINA3", "LCN2", "HSPB1", "VIM"
)

inflammatory_genes <- c(
  "IL6", "CXCL10", "CCL2", "CXCL8", "NFKBIA", "TNFAIP3", "SOCS3"
)

stress_genes <- c(
  "HSPA1A", "HSPA1B", "HSPB1", "DDIT3", "FOS", "JUN", "ATF3"
)

ecm_remodeling_genes <- c(
  "TNC", "COL1A1", "COL1A2", "COL4A1", "COL4A2", "MMP2", "MMP9"
)

astro_identity_present <- intersect(astro_identity_genes, rownames(SO_astro))
reactive_astro_present <- intersect(reactive_astrocyte_genes, rownames(SO_astro))
inflammatory_present <- intersect(inflammatory_genes, rownames(SO_astro))
stress_present <- intersect(stress_genes, rownames(SO_astro))
ecm_present <- intersect(ecm_remodeling_genes, rownames(SO_astro))

astro_identity_present

[1] "AQP4"    "GFAP"    "ALDH1L1" "SLC1A3"  "S100B"  

reactive_astro_present

[1] "CHI3L1"   "TNC"      "C3"       "SERPINA3" "HSPB1"    "VIM"     

inflammatory_present

[1] "IL6"     "CXCL10"  "CCL2"    "CXCL8"   "NFKBIA"  "TNFAIP3" "SOCS3"  

stress_present

[1] "HSPA1A" "HSPA1B" "HSPB1"  "DDIT3"  "FOS"    "JUN"    "ATF3"  

ecm_present

[1] "TNC"    "COL1A1" "COL1A2" "COL4A1" "COL4A2" "MMP2"   "MMP9"  

SO_astro <- AddModuleScore(
  SO_astro,
  features = list(astro_identity_present),
  name = "Astro_identity_score"
)

SO_astro <- AddModuleScore(
  SO_astro,
  features = list(reactive_astro_present),
  name = "Reactive_astro_score"
)

SO_astro <- AddModuleScore(
  SO_astro,
  features = list(inflammatory_present),
  name = "Inflammatory_score"
)

SO_astro <- AddModuleScore(
  SO_astro,
  features = list(stress_present),
  name = "Stress_score"
)

SO_astro <- AddModuleScore(
  SO_astro,
  features = list(ecm_present),
  name = "ECM_score"
)

FeaturePlot(
  SO_astro,
  features = c(
    "Astro_identity_score1",
    "Reactive_astro_score1",
    "Inflammatory_score1",
    "Stress_score1",
    "ECM_score1"
  ),
  ncol = 3,
  pt.size = 0.5
)

VlnPlot(
  SO_astro,
  features = c(
    "Reactive_astro_score1",
    "Inflammatory_score1",
    "Stress_score1",
    "ECM_score1"
  ),
  group.by = "condition",
  pt.size = 0
)

Step 11.6: Check Reactivity by Astrocyte Subcluster

We check whether reactivity is associated with specific astrocyte subclusters.

VlnPlot(
  SO_astro,
  features = c(
    "Reactive_astro_score1",
    "Inflammatory_score1",
    "Stress_score1",
    "ECM_score1"
  ),
  group.by = "seurat_clusters",
  pt.size = 0
)

table(SO_astro$seurat_clusters, SO_astro$condition)

   
       C   PD
  0 2115    6
  1 1015  619
  2 1131    2
  3  558  458
  4    6  904
  5  220  675
  6    6  754
  7    0  243
  8  217    3
  9   17    3

astro_cluster_summary <- SO_astro@meta.data %>%
  mutate(cluster = as.character(seurat_clusters)) %>%
  group_by(cluster) %>%
  summarise(
    n_cells = n(),
    mean_reactive = mean(Reactive_astro_score1, na.rm = TRUE),
    mean_inflammatory = mean(Inflammatory_score1, na.rm = TRUE),
    mean_stress = mean(Stress_score1, na.rm = TRUE),
    mean_ecm = mean(ECM_score1, na.rm = TRUE),
    percent_PD = mean(condition == "PD", na.rm = TRUE),
    percent_C = mean(condition == "C", na.rm = TRUE),
    .groups = "drop"
  ) %>%
  arrange(mean_reactive)

astro_cluster_summary

# A tibble: 10 × 8
   cluster n_cells mean_reactive mean_inflammatory mean_stress mean_ecm
   <chr>     <int>         <dbl>             <dbl>       <dbl>    <dbl>
 1 0          2121       -0.435          -0.0212       -0.333   -0.178 
 2 6           760       -0.429          -0.0380       -0.286   -0.174 
 3 1          1634       -0.357          -0.0263       -0.0947  -0.160 
 4 9            20       -0.148          -0.0244        0.0847  -0.0897
 5 2          1133        0.0775          0.0114        0.187   -0.0213
 6 8           220        0.179          -0.0105        0.414    0.161 
 7 3          1016        0.225           0.0133        0.221    0.0957
 8 5           895        0.240          -0.000901      0.112    0.0953
 9 7           243        0.400           0.0107        0.946    0.0391
10 4           910        0.713           0.105         0.327    0.290 
# ℹ 2 more variables: percent_PD <dbl>, percent_C <dbl>

Step 11.7: Convert Astrocytes to SingleCellExperiment

Slingshot requires a SingleCellExperiment object.

sce_astro <- as.SingleCellExperiment(SO_astro)

reducedDims(sce_astro)$UMAP <- Embeddings(SO_astro, "umap")

sce_astro$cluster <- as.character(SO_astro$seurat_clusters)
sce_astro$condition <- SO_astro$condition
sce_astro$Reactive_astro_score <- SO_astro$Reactive_astro_score1
sce_astro$Inflammatory_score <- SO_astro$Inflammatory_score1
sce_astro$Stress_score <- SO_astro$Stress_score1
sce_astro$ECM_score <- SO_astro$ECM_score1

Step 11.8: Choose Root Cluster

For astrocyte reactivity trajectory, the root should ideally be:

control-enriched
low reactive astrocyte score
low inflammatory score
baseline astrocyte-like

astro_root_summary <- SO_astro@meta.data %>%
  mutate(cluster = as.character(seurat_clusters)) %>%
  group_by(cluster) %>%
  summarise(
    n_cells = n(),
    mean_reactive = mean(Reactive_astro_score1, na.rm = TRUE),
    mean_inflammatory = mean(Inflammatory_score1, na.rm = TRUE),
    mean_stress = mean(Stress_score1, na.rm = TRUE),
    percent_PD = mean(condition == "PD", na.rm = TRUE),
    percent_C = mean(condition == "C", na.rm = TRUE),
    .groups = "drop"
  ) %>%
  filter(n_cells >= 20) %>%
  arrange(mean_reactive, mean_inflammatory, desc(percent_C))

astro_root_summary

# A tibble: 10 × 7
   cluster n_cells mean_reactive mean_inflammatory mean_stress percent_PD
   <chr>     <int>         <dbl>             <dbl>       <dbl>      <dbl>
 1 0          2121       -0.435          -0.0212       -0.333     0.00283
 2 6           760       -0.429          -0.0380       -0.286     0.992  
 3 1          1634       -0.357          -0.0263       -0.0947    0.379  
 4 9            20       -0.148          -0.0244        0.0847    0.15   
 5 2          1133        0.0775          0.0114        0.187     0.00177
 6 8           220        0.179          -0.0105        0.414     0.0136 
 7 3          1016        0.225           0.0133        0.221     0.451  
 8 5           895        0.240          -0.000901      0.112     0.754  
 9 7           243        0.400           0.0107        0.946     1      
10 4           910        0.713           0.105         0.327     0.993  
# ℹ 1 more variable: percent_C <dbl>

start_cluster <- astro_root_summary$cluster[1]

start_cluster

[1] "0"

Step 11.9: Run Slingshot on Astrocytes

sce_astro <- slingshot(
  sce_astro,
  clusterLabels = "cluster",
  reducedDim = "UMAP",
  start.clus = start_cluster
)

Step 11.10: Visualize Astrocyte Trajectory

umap_astro <- reducedDims(sce_astro)$UMAP

plot(
  umap_astro,
  col = as.factor(sce_astro$condition),
  pch = 16,
  asp = 1,
  xlab = "UMAP_1",
  ylab = "UMAP_2",
  main = "Slingshot Trajectory on Astrocyte UMAP"
)

lines(SlingshotDataSet(sce_astro), lwd = 2, col = "black")

legend(
  "topright",
  legend = levels(as.factor(sce_astro$condition)),
  col = seq_along(levels(as.factor(sce_astro$condition))),
  pch = 16
)

Step 11.11: Extract Astrocyte Pseudotime and Evaluate Lineages

Slingshot may return multiple lineages.
Therefore, we do not automatically choose the lineage with the largest number of cells.
Instead, we evaluate which lineage best represents astrocyte reactivity.

A biologically meaningful astrocyte reactivity lineage should show:

positive correlation with Reactive_astro_score
positive correlation with Stress_score
positive correlation with ECM_score
PD cells shifted toward later pseudotime

astro_pseudotime_mat <- slingPseudotime(sce_astro)

head(astro_pseudotime_mat)

                              Lineage1 Lineage2  Lineage3 Lineage4 Lineage5
GSM6459689_AAACGAAAGTGCCGAA-1       NA       NA 25.583213       NA       NA
GSM6459689_AACAAAGAGGTTGCCC-1       NA       NA 25.355640       NA       NA
GSM6459689_AACAAAGTCGAGATGG-1       NA       NA 25.340390       NA       NA
GSM6459689_AACAAGATCGTAACTG-1 8.433124 8.411299  8.361028 8.513503 8.837918
GSM6459689_AACAGGGAGAACTTCC-1       NA       NA 24.819177       NA       NA
GSM6459689_AACCTGACACTTCATT-1       NA       NA 24.601876       NA       NA

# Number of cells assigned to each lineage
lineage_cell_counts <- colSums(!is.na(astro_pseudotime_mat))

lineage_cell_counts

Lineage1 Lineage2 Lineage3 Lineage4 Lineage5 
    4798     5198     4722     4634     3691 

lineage_eval <- lapply(seq_len(ncol(astro_pseudotime_mat)), function(i) {
  
  pt_i <- astro_pseudotime_mat[, i]
  meta_i <- SO_astro@meta.data
  meta_i$pseudotime <- pt_i
  
  meta_i <- meta_i %>%
    filter(!is.na(pseudotime))
  
  tibble(
    lineage = paste0("Lineage", i),
    lineage_number = i,
    n_cells = nrow(meta_i),
    
    cor_reactive = cor(
      meta_i$pseudotime,
      meta_i$Reactive_astro_score1,
      method = "spearman",
      use = "complete.obs"
    ),
    
    cor_inflammatory = cor(
      meta_i$pseudotime,
      meta_i$Inflammatory_score1,
      method = "spearman",
      use = "complete.obs"
    ),
    
    cor_stress = cor(
      meta_i$pseudotime,
      meta_i$Stress_score1,
      method = "spearman",
      use = "complete.obs"
    ),
    
    cor_ecm = cor(
      meta_i$pseudotime,
      meta_i$ECM_score1,
      method = "spearman",
      use = "complete.obs"
    ),
    
    mean_pt_C = mean(meta_i$pseudotime[meta_i$condition == "C"], na.rm = TRUE),
    mean_pt_PD = mean(meta_i$pseudotime[meta_i$condition == "PD"], na.rm = TRUE),
    delta_PD_minus_C = mean_pt_PD - mean_pt_C
  )
}) %>%
  bind_rows() %>%
  arrange(
    desc(cor_reactive),
    desc(cor_stress),
    desc(cor_ecm),
    desc(delta_PD_minus_C),
    desc(n_cells)
  )

lineage_eval

# A tibble: 5 × 10
  lineage  lineage_number n_cells cor_reactive cor_inflammatory cor_stress
  <chr>             <int>   <int>        <dbl>            <dbl>      <dbl>
1 Lineage2              2    5198        0.783           0.304       0.641
2 Lineage3              3    4722        0.773           0.274       0.612
3 Lineage1              1    4798        0.761           0.265       0.575
4 Lineage4              4    4634        0.721           0.332       0.627
5 Lineage5              5    3691        0.136          -0.0228      0.269
# ℹ 4 more variables: cor_ecm <dbl>, mean_pt_C <dbl>, mean_pt_PD <dbl>,
#   delta_PD_minus_C <dbl>

Step 11.12: Select the Best Astrocyte Reactivity Lineage

We select the lineage that best matches the biological direction of astrocyte reactivity.

best_lineage <- lineage_eval %>%
  as.data.frame() %>%
  dplyr::filter(
    cor_reactive > 0,
    cor_stress > 0,
    cor_ecm > 0,
    delta_PD_minus_C > 0
  ) %>%
  dplyr::arrange(
    dplyr::desc(cor_reactive),
    dplyr::desc(cor_stress),
    dplyr::desc(cor_ecm),
    dplyr::desc(delta_PD_minus_C),
    dplyr::desc(n_cells)
  ) %>%
  dplyr::slice(1)

best_lineage

   lineage lineage_number n_cells cor_reactive cor_inflammatory cor_stress
1 Lineage2              2    5198    0.7831958        0.3042907  0.6408418
    cor_ecm mean_pt_C mean_pt_PD delta_PD_minus_C
1 0.7349624  11.43553   15.73092         4.295393

best_lineage_number <- best_lineage$lineage_number[1]

best_lineage_number

[1] 2

SO_astro$astro_pseudotime <- astro_pseudotime_mat[, best_lineage_number]

summary(SO_astro$astro_pseudotime)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.     NAs 
  0.000   3.634  10.014  12.236  21.145  29.915    3754 

If no lineage satisfies all conditions, manually inspect lineage_eval and choose the most biologically reasonable lineage.

# Example manual fallback:
# best_lineage_number <- 1
# SO_astro$astro_pseudotime <- astro_pseudotime_mat[, best_lineage_number]

Step 11.13: Visualize Selected Astrocyte Pseudotime

FeaturePlot(
  SO_astro,
  features = "astro_pseudotime",
  reduction = "umap",
  pt.size = 0.5
) +
  labs(
    title = paste0(
      "Astrocyte Slingshot Pseudotime: Lineage ",
      best_lineage_number
    )
  )

Step 11.14: Validate Astrocyte Trajectory Direction

A meaningful astrocyte reactivity trajectory should show increasing reactive, stress, and ECM remodeling signatures along pseudotime.

astro_pt_df <- SO_astro@meta.data %>%
  rownames_to_column("cell_id") %>%
  filter(!is.na(astro_pseudotime))

p1 <- ggplot(
  astro_pt_df,
  aes(x = astro_pseudotime, y = Reactive_astro_score1, color = condition)
) +
  geom_point(alpha = 0.25, size = 0.5) +
  geom_smooth(method = "loess", se = TRUE) +
  theme_classic() +
  labs(
    title = "Reactive Astrocyte Score Along Pseudotime",
    x = "Astrocyte pseudotime",
    y = "Reactive astrocyte score"
  )

p2 <- ggplot(
  astro_pt_df,
  aes(x = astro_pseudotime, y = Stress_score1, color = condition)
) +
  geom_point(alpha = 0.25, size = 0.5) +
  geom_smooth(method = "loess", se = TRUE) +
  theme_classic() +
  labs(
    title = "Stress Score Along Pseudotime",
    x = "Astrocyte pseudotime",
    y = "Stress score"
  )

p3 <- ggplot(
  astro_pt_df,
  aes(x = astro_pseudotime, y = ECM_score1, color = condition)
) +
  geom_point(alpha = 0.25, size = 0.5) +
  geom_smooth(method = "loess", se = TRUE) +
  theme_classic() +
  labs(
    title = "ECM Remodeling Score Along Pseudotime",
    x = "Astrocyte pseudotime",
    y = "ECM score"
  )

p4 <- ggplot(
  astro_pt_df,
  aes(x = astro_pseudotime, y = Inflammatory_score1, color = condition)
) +
  geom_point(alpha = 0.25, size = 0.5) +
  geom_smooth(method = "loess", se = TRUE) +
  theme_classic() +
  labs(
    title = "Inflammatory Score Along Pseudotime",
    x = "Astrocyte pseudotime",
    y = "Inflammatory score"
  )

(p1 | p2) / (p3 | p4)

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

Step 11.15: Condition Distribution Along Astrocyte Pseudotime

ggplot(
  astro_pt_df,
  aes(x = astro_pseudotime, fill = condition)
) +
  geom_density(alpha = 0.4) +
  theme_classic() +
  labs(
    title = "Condition Distribution Along Astrocyte Pseudotime",
    x = "Astrocyte pseudotime",
    y = "Density"
  )

Step 11.16: Cluster Distribution Along Astrocyte Pseudotime

This plot helps determine whether pseudotime reflects a transition from baseline astrocyte clusters to reactive astrocyte clusters.

p_cluster_pt <- ggplot(
  astro_pt_df,
  aes(x = as.factor(seurat_clusters), y = astro_pseudotime, fill = as.factor(seurat_clusters))
) +
  geom_boxplot(outlier.size = 0.3) +
  theme_classic() +
  labs(
    title = "Astrocyte Pseudotime by Subcluster",
    x = "Astrocyte subcluster",
    y = "Astrocyte pseudotime"
  ) +
  NoLegend()
p_cluster_pt

Step 11.17: Save Current Astrocyte Trajectory Progress

# Create results directory if it does not exist
results_dir <- "/Users/yoshimurasouhei/Downloads/010_school/4年生/bioinfomaticsリサーチクラークシップ/PD_lncRNA_Trajectory_2026/results"
scripts_dir <- "/Users/yoshimurasouhei/Downloads/010_school/4年生/bioinfomaticsリサーチクラークシップ/PD_lncRNA_Trajectory_2026/scripts"

dir.create(results_dir, showWarnings = FALSE, recursive = TRUE)
dir.create(scripts_dir, showWarnings = FALSE, recursive = TRUE)

# 1. Save Seurat object with astrocyte pseudotime
saveRDS(
  SO_astro,
  file = file.path(scripts_dir, "SO_11_ReactiveAstrocytePseudotime_Lineage2.RDS")
)

# 2. Save lineage evaluation table
write.csv(
  lineage_eval,
  file = file.path(results_dir, "Step11_lineage_evaluation.csv"),
  row.names = FALSE
)

# 3. Save astrocyte pseudotime metadata
write.csv(
  astro_pt_df,
  file = file.path(results_dir, "Step11_astrocyte_pseudotime_metadata.csv"),
  row.names = FALSE
)

# 4. Save pseudotime matrix
write.csv(
  as.data.frame(astro_pseudotime_mat),
  file = file.path(results_dir, "Step11_slingshot_pseudotime_matrix.csv"),
  row.names = TRUE
)

# 5. Save cluster summary if it exists
if (exists("astro_cluster_summary")) {
  write.csv(
    astro_cluster_summary,
    file = file.path(results_dir, "Step11_astrocyte_cluster_summary.csv"),
    row.names = FALSE
  )
}

print("Step 11 progress saved successfully.")

[1] "Step 11 progress saved successfully."

ggsave(
  filename = file.path(results_dir, "Step11_astrocyte_pseudotime_by_subcluster.png"),
  plot = p_cluster_pt,
  width = 10,
  height = 5,
  dpi = 300
)

ggsave(
  filename = file.path(results_dir, "Step11_astrocyte_pseudotime_by_subcluster.pdf"),
  plot = p_cluster_pt,
  width = 10,
  height = 5
)