Step1:Integration of normal CD4Tcells to project L1

load libraries————————————

Read Seurat object with load as you save it with save() function

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

1. Subset Normal CD4⁺ T Cells for Reference

# Subset normal CD4+ T cells from merged object
reference_cd4 <- subset(
  All_samples_Merged,
  subset = cell_line %in% c("CD4Tcells_lab", "CD4Tcells_10x")
)


DefaultAssay(reference_cd4) <- "RNA"  

2. Normalize & Integrate Reference (Accommodate Multiple Donors)

ref_list <- SplitObject(reference_cd4, split.by = "cell_line")

# 🔴 Run SCTransform on each subset to ensure consistent variable features
ref_list <- lapply(ref_list, function(x) {
  x <- SCTransform(x, verbose = FALSE)
  # 🔴 Store top 3000 variable features explicitly
  VariableFeatures(x) <- head(VariableFeatures(x), 3000)
  return(x)
})

# Run PCA on each SCT-normalized subset (required for RPCA)
ref_list <- lapply(ref_list, function(x) {
  x <- RunPCA(x, assay = "SCT", verbose = FALSE)
  return(x)
})

# Select integration features
ref_features <- SelectIntegrationFeatures(object.list = ref_list, nfeatures = 3000)

# 🔴 REMOVE TCR/TRAV/TRBV GENES FROM INTEGRATION FEATURES
ref_features <- ref_features[!grepl("^TRBV|^TRAV", ref_features)]  # 🔴 CHANGE

# Prepare SCT integration
ref_list <- PrepSCTIntegration(object.list = ref_list, anchor.features = ref_features)

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

# Find integration anchors using RPCA reduction
ref_anchors <- FindIntegrationAnchors(
  object.list = ref_list,
  anchor.features = ref_features,
  normalization.method = "SCT",
  reduction = "rpca"
)

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

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

# Integrate data
reference_integrated <- IntegrateData(anchorset = ref_anchors, normalization.method = "SCT")

[1] 1
[1] 2

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

# 🔴 Store SCT variable features explicitly inside object
VariableFeatures(reference_integrated) <- ref_features   # 🔴 IMPORTANT CHANGE

# 🔴 Set default assay to SCT for mapping L1 later
DefaultAssay(reference_integrated) <- "SCT"  # 🔴 CHANGE

3. Clustering & Dimensionality Reduction

reference_integrated <- FindClusters(reference_integrated, resolution = 0.2)

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 8610
Number of edges: 291335

Running Louvain algorithm...

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

Maximum modularity in 10 random starts: 0.9143
Number of communities: 6
Elapsed time: 0 seconds

ElbowPlot(reference_integrated, ndims = 50)

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

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

# 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")

Save the mapped query object (Sézary cell lines projected onto reference trajectory):

saveRDS(reference_integrated, file = "sezary_cell_lines_mapped_to_cd4_reference_integrated_before_Monocle3_03-09-2025.rds")

FeaturePlots

library(Seurat)
library(ggplot2)
library(patchwork)

# ---- Define marker lists for differentiation states ----
marker_list <- list(
  Tnaive = c("CCR7","SELL","LEF1","TCF7","IL7R","CD27"),
  Tcm    = c("CCR7","SELL","CD27","IL7R","BCL2","TCF7"),
  Tem    = c("CCR6","CXCR3","GZMK","PRF1","IFNG","CD45RO"),
  Temra  = c("GZMB","PRF1","KLRG1","CX3CR1","CD45RA"),
  Tex    = c("PDCD1","CTLA4","LAG3","TIGIT","TOX","ENTPD1"),
  CD4CTL = c("GZMB","PRF1","NKG7","KLRG1","CX3CR1") # CD4 cytotoxic markers
)

# ---- Keep only markers present in the dataset ----
marker_list <- lapply(marker_list, function(x) x[x %in% rownames(reference_integrated)])

# ---- Compute module scores one by one ----
for (state in names(marker_list)) {
  reference_integrated <- AddModuleScore(
    reference_integrated,
    features = list(marker_list[[state]]),
    name = state
  )
}

# Module scores are named Tnaive1, Tcm1, Tem1, Temra1, Tex1, CD4CTL1

# ---- Plot module scores individually with lightblue → red gradient and labels ----
plots <- lapply(names(marker_list), function(state) {
  FeaturePlot(
    reference_integrated,
    features = paste0(state,"1"),
    reduction = "umap",
    cols = c("lightblue","red"),
    label = TRUE
  ) + ggtitle(paste(state, "Markers")) + theme(plot.title = element_text(hjust = 0.5))
})

# ---- Display all plots in a grid ----
wrap_plots(plots, ncol = 3)

FeaturePlots

# ---- Load libraries ----
library(Seurat)
library(ggplot2)
library(patchwork)

# ---- Define markers for differentiation states ----
tnaive_markers <- c("CCR7", "SELL", "LEF1", "TCF7", "IL7R", "CD27", "CD45RA")
tcm_markers    <- c("CCR7", "SELL", "CD45RO", "IL7R", "CD27")
tem_markers    <- c("CCR6", "CXCR3", "GZMK", "PRF1", "IFNG", "CD45RO")
temra_markers  <- c("GZMB", "PRF1", "KLRG1", "CX3CR1", "CD45RA")
tex_markers    <- c("PDCD1", "CTLA4", "LAG3", "TIGIT", "TOX", "ENTPD1")
cd4ctl_markers <- c("GZMB", "PRF1", "NKG7", "KLRG1", "CX3CR1")  # CD4 cytotoxic markers

# ---- Keep only markers present in the dataset ----
tnaive_markers <- tnaive_markers[tnaive_markers %in% rownames(reference_integrated)]
tcm_markers    <- tcm_markers[tcm_markers %in% rownames(reference_integrated)]
tem_markers    <- tem_markers[tem_markers %in% rownames(reference_integrated)]
temra_markers  <- temra_markers[temra_markers %in% rownames(reference_integrated)]
tex_markers    <- tex_markers[tex_markers %in% rownames(reference_integrated)]
cd4ctl_markers <- cd4ctl_markers[cd4ctl_markers %in% rownames(reference_integrated)]

# ---- Function to generate patchwork feature plots ----
plot_markers_grid <- function(marker_vector, state_name) {
  plots <- lapply(marker_vector, function(gene){
    FeaturePlot(reference_integrated, features = gene, reduction = "umap",
                cols = c("lightblue", "red"), label = TRUE) +
      theme(plot.title = element_text(hjust = 0.5))
  })
  wrap_plots(plots) + plot_annotation(title = paste(state_name, "Marker Expression"))
}

# ---- Generate grids for each state ----
tnaive_plot <- plot_markers_grid(tnaive_markers, "Tnaive")
tcm_plot    <- plot_markers_grid(tcm_markers, "Tcm")
tem_plot    <- plot_markers_grid(tem_markers, "Tem")
temra_plot  <- plot_markers_grid(temra_markers, "Temra")
tex_plot    <- plot_markers_grid(tex_markers, "Tex")
cd4ctl_plot <- plot_markers_grid(cd4ctl_markers, "CD4 CTL")

# ---- Display plots ----
tnaive_plot

tcm_plot

tem_plot

temra_plot

tex_plot

cd4ctl_plot

# ---- Compute module scores for each state ----
reference_integrated <- AddModuleScore(reference_integrated, features = list(tnaive_markers), name = "Tnaive_Score")
reference_integrated <- AddModuleScore(reference_integrated, features = list(tcm_markers), name = "Tcm_Score")
reference_integrated <- AddModuleScore(reference_integrated, features = list(tem_markers), name = "Tem_Score")
reference_integrated <- AddModuleScore(reference_integrated, features = list(temra_markers), name = "Temra_Score")
reference_integrated <- AddModuleScore(reference_integrated, features = list(tex_markers), name = "Tex_Score")
reference_integrated <- AddModuleScore(reference_integrated, features = list(cd4ctl_markers), name = "CD4CTL_Score")

# ---- FeaturePlot for module scores ----
score_plots <- list(
  FeaturePlot(reference_integrated, features = "Tnaive_Score1", reduction = "umap", cols = c("lightblue","red"), label = TRUE) + ggtitle("Tnaive Module Score") + theme(plot.title = element_text(hjust=0.5)),
  FeaturePlot(reference_integrated, features = "Tcm_Score1", reduction = "umap", cols = c("lightblue","red"), label = TRUE) + ggtitle("Tcm Module Score") + theme(plot.title = element_text(hjust=0.5)),
  FeaturePlot(reference_integrated, features = "Tem_Score1", reduction = "umap", cols = c("lightblue","red"), label = TRUE) + ggtitle("Tem Module Score") + theme(plot.title = element_text(hjust=0.5)),
  FeaturePlot(reference_integrated, features = "Temra_Score1", reduction = "umap", cols = c("lightblue","red"), label = TRUE) + ggtitle("Temra Module Score") + theme(plot.title = element_text(hjust=0.5)),
  FeaturePlot(reference_integrated, features = "Tex_Score1", reduction = "umap", cols = c("lightblue","red"), label = TRUE) + ggtitle("Tex Module Score") + theme(plot.title = element_text(hjust=0.5)),
  FeaturePlot(reference_integrated, features = "CD4CTL_Score1", reduction = "umap", cols = c("lightblue","red"), label = TRUE) + ggtitle("CD4 CTL Module Score") + theme(plot.title = element_text(hjust=0.5))
)

wrap_plots(score_plots, ncol = 2)

NA
NA

Save the mapped query object (Sézary cell lines projected onto reference trajectory):

saveRDS(reference_integrated, file = "sezary_cell_lines_mapped_to_cd4_reference_integrated_before_Monocle3_03-09-2025.rds")

4. Trajectory and Pseudotime with Monocle3

library(monocle3)
library(SeuratWrappers)
library(Matrix)

cds <- as.cell_data_set(reference_integrated)
cds <- cluster_cells(cds, reduction_method = "UMAP")
cds <- learn_graph(cds, use_partition = TRUE)

  |                                                                                                            
  |                                                                                                      |   0%
  |                                                                                                            
  |======================================================================================================| 100%

naive_markers <- c("CCR7", "SELL", "LEF1", "TCF7", "CD45RA", "PTPRC")
naive_markers <- naive_markers[naive_markers %in% rownames(cds)]

# Extract log-normalized expression or fallback to counts log-transformed
if("logcounts" %in% assayNames(cds)) {
  expr_mat <- assay(cds, "logcounts")
} else {
  expr_mat <- log1p(assay(cds, "counts"))
}

naive_score <- Matrix::colMeans(expr_mat[naive_markers, , drop = FALSE])
threshold <- quantile(naive_score, 0.95)
root_cells <- names(naive_score[naive_score > threshold])

cds <- order_cells(cds, root_cells = root_cells)
reference_integrated$pseudotime <- pseudotime(cds)

plot_cells(cds, color_cells_by = "pseudotime", show_trajectory_graph = TRUE)

# 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")

Save the mapped query object (Sézary cell lines projected onto reference trajectory):

saveRDS(reference_integrated, file = "sezary_cell_lines_mapped_to_cd4_reference_integrated_before_Query_Projection_03-09-2025.rds")