Damage trajectory analysis using current neuron Seurat object
Setup: Environment and Data
In this step, we start from the current neuron Seurat object generated in Step 09.
The goal is to infer a neuronal damage trajectory rather than a developmental trajectory.
We use the neuron-subclustered Seurat object as the starting point.
library(Seurat)library(tidyverse)library(patchwork)library(SingleCellExperiment)library(slingshot)library(tradeSeq)# Load current neuron Seurat objectimport_path <-"/Users/yoshimurasouhei/Downloads/010_school/4年生/bioinfomaticsリサーチクラークシップ/PD_2026/scripts/SO_09_NeuronSub.RDS"SO_neuron <-readRDS(import_path)SO_neuron
An object of class Seurat
47773 features across 1413 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 10.1: Check Metadata and Reductions
Before trajectory analysis, we confirm that the object contains UMAP, PCA, clusters, and condition metadata.
# Keep only genes present in the Seurat objectpd_damage_genes_present <-intersect(pd_damage_genes, rownames(SO_neuron))excitatory_genes_present <-intersect(excitatory_identity_genes, rownames(SO_neuron))inhibitory_genes_present <-intersect(inhibitory_identity_genes, rownames(SO_neuron))pd_damage_genes_present
For the first pseudotime test, we use all neuron subclusters.
Later, this should be repeated separately for excitatory and inhibitory neurons.
SO_traj <- SO_neuronSO_traj
An object of class Seurat
47773 features across 1413 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 10.6: Convert Seurat Object to SingleCellExperiment
Slingshot works with SingleCellExperiment.
We use the existing UMAP coordinates and Seurat clusters.
sce <-as.SingleCellExperiment(SO_traj)# Add UMAP coordinates from Seurat to SingleCellExperimentreducedDims(sce)$UMAP <-Embeddings(SO_traj, "umap")# Store cluster labelssce$cluster <-as.character(SO_traj$seurat_clusters)# Store conditionsce$condition <- SO_traj$condition# Store PD damage scoresce$PD_damage_score <- SO_traj$PD_damage_score1
Step 10.7: Run Slingshot Pseudotime
We infer lineage structure using UMAP and Seurat clusters.
Pseudotime is useful only if it reflects neuronal damage.
Therefore, PD damage score should generally increase along pseudotime.
pt_df <- SO_traj@meta.data %>%rownames_to_column("cell_id") %>%filter(!is.na(pseudotime_slingshot))ggplot(pt_df, aes(x = pseudotime_slingshot, y = PD_damage_score1, color = condition)) +geom_point(alpha =0.25, size =0.5) +geom_smooth(method ="loess", se =TRUE) +theme_classic() +labs(title ="PD Damage Score Along Pseudotime",x ="Slingshot pseudotime",y ="PD damage score" )
`geom_smooth()` using formula = 'y ~ x'
ggplot(pt_df, aes(x = pseudotime_slingshot, fill = condition)) +geom_density(alpha =0.4) +theme_classic() +labs(title ="Condition Distribution Along Pseudotime",x ="Slingshot pseudotime",y ="Density" )
Step 10.11: Fit Gene Expression Along Pseudotime
We use tradeSeq to identify genes whose expression changes along pseudotime.
For memory safety, we first test this on variable genes.
counts <-GetAssayData(SO_traj, assay ="RNA", layer ="counts")# Start with only 500 variable genes for speedgenes_use <-head(VariableFeatures(SO_traj), 500)genes_use <-intersect(genes_use, rownames(counts))counts_use <- counts[genes_use, ]# Extract pseudotime and cell weightspt <-slingPseudotime(sce)cw <-slingCurveWeights(sce)# Keep only cells with non-NA pseudotime in at least one lineagecells_keep <-rowSums(!is.na(pt)) >0# Subset counts, pseudotime, and cell weightscounts_use <- counts_use[, cells_keep]pt_use <- pt[cells_keep, , drop =FALSE]cw_use <- cw[cells_keep, , drop =FALSE]# Replace remaining NA pseudotime values with 0# and set corresponding weights to 0pt_use[is.na(pt_use)] <-0cw_use[is.na(cw_use)] <-0dim(counts_use)
