Step2:L2 projection of ReferenceCD4Tcells

1. load libraries

2. Read object with all samples

All_samples_Merged <- readRDS("../../0-Seurat_RDS_OBJECT_FINAL/All_samples_Merged_with_STCAT_and_renamed_FINAL.rds")

3. Read object with all Reference_CD4Tcells

 reference_integrated <- readRDS("Step1_sezary_cell_lines_mapped_to_cd4_reference_integrated_before_Query_Projection_03-09-2025.rds")

Trajectory and Pseudotime with Monocle3

reference_integrated$pseudotime[!is.finite(reference_integrated$pseudotime)] <- NA


max_finite <- max(reference_integrated$pseudotime, na.rm = TRUE)
reference_integrated$pseudotime[!is.finite(reference_integrated$pseudotime)] <- max_finite

summary(reference_integrated$pseudotime)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.00000 0.01969 0.19290 0.93865 0.90736 5.63665 

table(is.finite(reference_integrated$pseudotime))

TRUE 
8610 

FeaturePlot(
  reference_integrated,
  features = "pseudotime",
  reduction = "umap",
  cols = c("lightblue", "red"),
  label = TRUE
) + ggtitle("UMAP of Integrated CD4⁺ T Cells by Pseudotime") +
  theme(plot.title = element_text(hjust = 0.5))

NA
NA

Trajectory and Pseudotime with Monocle3

# Visualize UMAP colored by original donor (cell_line)
DimPlot(reference_integrated, group.by = "Prediction", reduction = "umap") +
  ggtitle("UMAP of Integrated CD4⁺ T Cells")

# Visualize UMAP colored by original donor (cell_line)
DimPlot(reference_integrated, group.by = "predicted.celltype.l2", reduction = "umap") +
  ggtitle("UMAP of Integrated CD4⁺ T Cells")

3. Subset L2 cells

# Subset L2 cells
L2 <- subset(All_samples_Merged, subset = cell_line == "L2")


# Use SCT assay for both
DefaultAssay(reference_integrated) <- "SCT"
DefaultAssay(L2) <- "SCT"


# Make sure both have variable features set
if (length(VariableFeatures(reference_integrated)) == 0) {
  reference_integrated <- FindVariableFeatures(reference_integrated, assay = "SCT", selection.method = "vst", nfeatures = 3000)
}
if (length(VariableFeatures(L2)) == 0) {
  L2 <- FindVariableFeatures(L2, assay = "SCT", selection.method = "vst", nfeatures = 3000)
}

# Find anchors
anchors <- FindTransferAnchors(
  reference = reference_integrated,
  query = L2,
  normalization.method = "SCT",
  reference.reduction = "pca",
  dims = 1:18
)

[1] "Given reference assay has multiple sct models, selecting model with most cells for finding transfer anchors"

reference_integrated <- RunUMAP(
  reference_integrated,
  reduction = "pca",
  dims = 1:18,
  assay = "SCT",
  return.model = TRUE,   # For Seurat v5, save.model = TRUE deprecated; use return.model
)

Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

# Map query onto reference and transfer pseudotime
mapped_L2 <- MapQuery(
  anchorset = anchors,
  query = L2,
  reference = reference_integrated,
  refdata = list(
    pseudotime = reference_integrated$pseudotime,
    seurat_clusters = reference_integrated$seurat_clusters,
    Prediction = reference_integrated$Prediction,
    predicted.celltype.l2 = reference_integrated$predicted.celltype.l2
  ),
  reference.reduction = "pca",
  reduction.model = "umap"
)

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

  |                                                  | 0 % ~calculating  
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=01s  

Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

# Convert numeric pseudotime vector into a 1-row matrix with rowname
pseudo_mat <- matrix(reference_integrated$pseudotime, nrow = 1)
colnames(pseudo_mat) <- colnames(reference_integrated)  # cells
rownames(pseudo_mat) <- "pseudotime"                     # feature name

# Transfer numeric pseudotime
pseudotime_transfer <- TransferData(
  anchorset = anchors,
  refdata = pseudo_mat,
  weight.reduction = L2[["pca"]],
  dims = 1:18
)

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

# Extract the data from the assay
pseudotime_vector <- GetAssayData(pseudotime_transfer, slot = "data")["pseudotime", ]

# Add to metadata
mapped_L2$pseudotime <- pseudotime_vector

#  Visualize pseudotime
FeaturePlot(mapped_L2,
            features = "pseudotime",
            reduction = "ref.umap",
            cols = c("lightblue","red"),
            label = TRUE) +
  ggtitle("Pseudotime Mapping of L2 onto Reference") +
  theme(plot.title = element_text(hjust = 0.5))

#  Visualize transferred cell type annotations
DimPlot(mapped_L2, group.by = "Prediction", reduction = "ref.umap") +
  ggtitle("Mapped L2: Prediction Annotations") +
  theme(plot.title = element_text(hjust = 0.5))

DimPlot(mapped_L2, group.by = "predicted.celltype.l2", reduction = "ref.umap") +
  ggtitle("Mapped L2: predicted.celltype.l2 Annotations") +
  theme(plot.title = element_text(hjust = 0.5))

Mapping of L2 on Reference

library(ggplot2)

# Prepare reference data
ref_df <- data.frame(Embeddings(reference_integrated, "umap"))
ref_df$dataset <- "Reference"

# Prepare L2 query data
query_df <- data.frame(Embeddings(mapped_L2, "ref.umap"))
query_df$dataset <- "L2"
query_df$pseudotime <- mapped_L2$pseudotime  # if pseudotime exists

# Plot reference in grey, L2 cells colored by pseudotime
ggplot() +
  geom_point(data = ref_df, aes(x = umap_1, y = umap_2),
             color = "grey80", size = 0.5) +
  geom_point(data = query_df, aes(x = refUMAP_1, y = refUMAP_2, color = pseudotime),
             size = 1) +
  scale_color_gradient(low = "lightblue", high = "red") +
  theme_classic() +
  ggtitle("L2 cells projected onto reference UMAP") +
  theme(plot.title = element_text(hjust = 0.5))

NA
NA

Mapping of L2 on Reference (Azimuth)

library(ggplot2)

# Reference coordinates (grey background)
ref_df <- data.frame(Embeddings(reference_integrated, "umap"))
ref_df$dataset <- "Reference"

# Query (L2) coordinates
query_df <- data.frame(Embeddings(mapped_L2, "ref.umap"))
query_df$dataset <- "L2"

# Make sure annotations exist in mapped_L2
query_df$celltype <- mapped_L2$predicted.celltype.l2  # or "Prediction"

# Plot reference in grey and L2 cells colored by annotation
ggplot() +
  geom_point(data = ref_df, aes(x = umap_1, y = umap_2),
             color = "grey80", size = 0.5) +
  geom_point(data = query_df, aes(x = refUMAP_1, y = refUMAP_2, color = celltype),
             size = 1) +
  scale_color_brewer(palette = "Set1") +
  theme_classic() +
  ggtitle("L2 cells projected onto reference UMAP with annotations") +
  theme(plot.title = element_text(hjust = 0.5))

Mapping of L2 on Reference (STCAT)

library(ggplot2)

# Reference coordinates (grey background)
ref_df <- data.frame(Embeddings(reference_integrated, "umap"))
ref_df$dataset <- "Reference"

# Query (L2) coordinates
query_df <- data.frame(Embeddings(mapped_L2, "ref.umap"))
query_df$dataset <- "L2"

# Make sure annotations exist in mapped_L2
query_df$celltype <- mapped_L2$Prediction  # or "Prediction"

# Plot reference in grey and L2 cells colored by annotation
ggplot() +
  geom_point(data = ref_df, aes(x = umap_1, y = umap_2),
             color = "grey80", size = 0.5) +
  geom_point(data = query_df, aes(x = refUMAP_1, y = refUMAP_2, color = celltype),
             size = 1) +
  scale_color_brewer(palette = "Set1") +
  theme_classic() +
  ggtitle("L2 cells projected onto reference UMAP with annotations") +
  theme(plot.title = element_text(hjust = 0.5))

4. Subset L2 Projection

# 1. Define output directory
out_dir <- "results/L2_projection"
if(!dir.exists(out_dir)) dir.create(out_dir, recursive = TRUE)

save(mapped_L2, file = file.path(out_dir, "L2_mapped_on_reference_with_pseudotime-3-8-2025.Robj"))