Publication-Ready Single-Cell RNA-seq Visualization Pipeline: Analysis of Harmony-Integrated Seurat Object

This R Markdown document outlines a complete, step-by-step pipeline for generating all essential and advanced single-cell RNA-seq visualizations required for a high-impact publication in cancer genomics. The methodologies are rigorously derived from the latest literature (e.g., Herrera et al., Xue et al., Borcherding et al.).

CRITICAL NOTE FOR EXECUTION: This script is designed to operate on your Harmony-integrated Seurat object, which is assumed to be named All_samples_Merged. The R code below uses a placeholder PBMC dataset to ensure the script is fully executable and demonstrates the expected output. To run this pipeline on your data, you MUST replace the placeholder data loading and processing steps in the data_prep chunk with the command to load your object: ss.integrated <- readRDS("path/to/your/All_samples_Merged.rds")object is assumed to have a dimensionality reduction named harmony and a metadata column for cell type annotation named cell_type.

Setup and Data Preparation

The following chunk loads all necessary libraries, defines publication-quality parameters, and sets up the output directory.

# Load Libraries (Ensure these are installed in your R environment)
library(Seurat)
library(ggplot2)
library(patchwork)
library(dplyr)
library(schex) # For Hexbin Plots
# library(monocle3) # For Trajectory Analysis
# library(CellChat) # For Cell-Cell Communication
# library(scRepertoire) # For TCR Clonality

# --- Define Publication Parameters ---
OUTPUT_DIR <- "./publication_figures/"
PUB_WIDTH <- 10
PUB_HEIGHT <- 8
PUB_RES <- 300 # DPI (Dots Per Inch)

if (!dir.exists(OUTPUT_DIR)) {
  dir.create(OUTPUT_DIR)
}

message(paste("Output directory created at:", OUTPUT_DIR))

Data Object Initialization and Validation

This section simulates your pre-processed, Harmony-integrated Seurat object. We load a standard dataset and process it up to the clustering and UMAP stage.

# Load the PBMC 3k dataset (Placeholder for your Harmony-integrated object)
 ss.integrated <- readRDS("/home/nabbasi/cluster_home/PHD_3rd_YEAR_Analysis/0-Seurat_RDS_OBJECT_FINAL/All_samples_Merged_with_Renamed_Clusters_final-26-10-2025.rds") # <-- REPLACE THIS LINE

I. Foundational Visualizations: Establishing Cellular Identity and Heterogeneity

These plots are essential for establishing the cellular landscape and validating cell type annotation.

1. Dimensionality Reduction Visualization (UMAP)

Purpose: Global visualization of cellular heterogeneity.

umap_cluster_plot <- DimPlot(ss.integrated, reduction = "umap", group.by = "seurat_clusters",label = TRUE, pt.size = 0.5) 

# Display plot in Rmd
print(umap_cluster_plot)

# Save high-resolution version
ggsave(paste0(OUTPUT_DIR, "UMAP_Cluster.png"), plot = umap_cluster_plot, width = PUB_WIDTH, height = PUB_HEIGHT, dpi = PUB_RES)

umap_cellline_plot <- DimPlot(ss.integrated, reduction = "umap", group.by = "orig.ident",label = TRUE, pt.size = 0.5) 

# Display plot in Rmd
print(umap_cellline_plot)

# Save high-resolution version
ggsave(paste0(OUTPUT_DIR, "umap_cellline_plot.png"), plot = umap_cellline_plot, width = PUB_WIDTH, height = PUB_HEIGHT, dpi = PUB_RES)

umap_cluster_A_plot <- DimPlot(ss.integrated, reduction = "umap", group.by = "renamed_clusters",label = F, pt.size = 0.5) 

# Display plot in Rmd
print(umap_cluster_A_plot)

# Save high-resolution version
ggsave(paste0(OUTPUT_DIR, "umap_cluster_A_plot.png"), plot = umap_cluster_A_plot, width = 12, height = 6, dpi = PUB_RES)

library(SCpubr)
umap_cluster_plot <- SCpubr::do_DimPlot(ss.integrated, reduction = "umap", group.by = "seurat_clusters",label = TRUE, pt.size = 0.5, legend.position = "right") 

# Display plot in Rmd
print(umap_cluster_plot)

# Save high-resolution version
ggsave(paste0(OUTPUT_DIR, "UMAP_Cluster-scpubr.png"), plot = umap_cluster_plot, width = PUB_WIDTH, height = PUB_HEIGHT, dpi = PUB_RES)

library(SCpubr)
umap_cellline_plot <- SCpubr::do_DimPlot(ss.integrated, reduction = "umap", group.by = "orig.ident",label = TRUE, pt.size = 0.5, legend.position = "right") 

# Display plot in Rmd
print(umap_cellline_plot)

# Save high-resolution version
ggsave(paste0(OUTPUT_DIR, "UMAP_cellline-scpubr.png"), plot = umap_cellline_plot, width = PUB_WIDTH, height = PUB_HEIGHT, dpi = PUB_RES)

library(SCpubr)
umap_cluster_A_plot <- do_DimPlot(ss.integrated, reduction = "umap", group.by = "renamed_clusters",label = F, pt.size = 0.5, legend.position = "right") 

# Display plot in Rmd
print(umap_cluster_A_plot)

# Save high-resolution version
ggsave(paste0(OUTPUT_DIR, "umap_cluster_A_plot-scpubr.png"), plot = umap_cluster_A_plot, width = 12, height = 6, dpi = PUB_RES)

2. Feature Plot: Gene Expression Mapping for Cell Type Confirmation

Purpose: Visualizes the expression level of a key gene across the UMAP space.

3. Marker Gene Expression: Dot Plot and Heatmap for Cluster Annotation

Purpose: Systematically compares the expression of top differentially expressed genes (DEGs) across all cell types for annotation confirmation.

library(dplyr)

# Precise blacklist for uninformative genes
blacklist_patterns <- c(
  "^TRAV", "^TRBV", "^TRGV", "^TRDV", "^TRBC", "^TRAC", "^TRDC", "^TRGC", # TCR
  "^IGH", "^IGK", "^IGL", "^IGJ",                                         # Ig genes
  "^RPL", "^RPS",                                                         # ribosomal
  "^MT-",                                                                 # mitochondria
  "^HBA", "^HBB", "^HB[ABZ]",                                             # hemoglobins
  "^NEAT1$", "^MALAT1$",                                                  # optional lncRNAs
  "^XIST$"                              )

blacklist_regex <- paste(blacklist_patterns, collapse = "|")

# Preview which markers will be removed
to_remove <- ss.integrated.markers %>%
  filter(grepl(blacklist_regex, gene, ignore.case = TRUE))
message("Rows to remove: ", nrow(to_remove))
head(to_remove$gene)

[1] "TRAV17"      "TRAV9-2"     "RPL22L1"     "NEAT1"       "TRAV38-2DV8" "TRGV2"      

# Filter markers (keep important metabolic/proliferation genes)
SS_markers_filtered <- ss.integrated.markers %>%
  filter(!grepl(blacklist_regex, gene, ignore.case = TRUE))

top5 <- SS_markers_filtered %>%
  filter(p_val_adj < 0.05) %>%  
  group_by(cluster) %>%
  slice_max(n = 5, order_by = avg_log2FC, with_ties = FALSE)

# Dot Plot
dot_plot_markers <- DotPlot(ss.integrated, features = unique(top5$gene)) +
  RotatedAxis() +
  ggtitle("Top Marker Gene Expression Across Cell Types")

# Display plot in Rmd
print(dot_plot_markers)

# Save high-resolution version
ggsave(paste0(OUTPUT_DIR, "DotPlot_Markers.png"), plot = dot_plot_markers, width = 26, height = 6, dpi = PUB_RES)

II. Quantitative and Comparative Visualizations

These plots are used to quantify differences between groups or visualize specific gene signatures.

4. Quantitative Analysis: Cell Type Composition Across Conditions

Purpose: Visualizes the change in the proportion of cell types across different conditions or samples (simulated here by ‘orig.ident’).

# Create a dummy 'orig.ident' column to simulate different samples/conditions
ss.integrated$orig.ident <- sample(c("L1", "L2", "L3", "L4", "L5", "L6", "L7"), size = ncol(ss.integrated), replace = TRUE)

cell_counts <- as.data.frame(table(ss.integrated$Prediction, ss.integrated$orig.ident))
names(cell_counts) <- c("Prediction", "orig.ident", "Count")

composition_plot <- ggplot(cell_counts, aes(x = orig.ident, y = Count, fill = Prediction)) +
  geom_bar(stat = "identity", position = "fill") +
  scale_y_continuous(labels = scales::percent) +
  labs(y = "Proportion of Cells", x = "Sample/Condition", fill = "Prediction") +
  ggtitle("Cell Type Composition Across Simulated Samples") +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

# Display plot in Rmd
print(composition_plot)

# Save high-resolution version
ggsave(paste0(OUTPUT_DIR, "CellComposition_BarPlot.png"), plot = composition_plot, width = PUB_WIDTH, height = PUB_HEIGHT, dpi = PUB_RES)

5. High-Resolution Density Visualization (Hexbin Plot)

Purpose: High-resolution density visualization for large datasets, overcoming point overlap (Borcherding et al. methodology).

library(schex)
# Requires the schex package
ss.integrated.hex <- make_hexbin(ss.integrated, nbins = 80, dimension_reduction = "UMAP")

# Plot density of a key marker (e.g., CD14 for Monocytes)
hexbin_plot_CD74 <- hexbin_plot(ss.integrated.hex, col = "CD74") +
  ggtitle("Hexbin Plot of CD14 Expression Density") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

# Display plot in Rmd
print(hexbin_plot_CD74)

# Save high-resolution version
ggsave(paste0(OUTPUT_DIR, "HexbinPlot_CD74_Density.png"), plot = hexbin_plot_CD74, width = PUB_WIDTH, height = PUB_HEIGHT, dpi = PUB_RES)

III. Mechanistic and Trajectory Visualizations (Advanced Templates)

These advanced analyses require specialized packages and often additional data (TCR, CNV). The code below provides the structure and key functions.

6. Functional State Mapping: Pathway Activity Score (Module Score)

Purpose: Maps the functional state (e.g., proliferation, exhaustion) of cells onto the UMAP space.

# NOTE: This is a template using Seurat's AddModuleScore (similar to UCell/AUCell)
# We will use a dummy gene set for "T Cell Exhaustion" for demonstration purposes.
exhaustion_genes <- list(c("PDCD1", "CTLA4", "LAG3", "TIGIT", "HAVCR2"))

ss.integrated <- AddModuleScore(ss.integrated, features = exhaustion_genes, name = "Exhaustion_Score")

pathway_plot <- FeaturePlot(ss.integrated, features = "Exhaustion_Score1", pt.size = 0.5, reduction = "umap") +
  ggtitle("T Cell Exhaustion Pathway Activity Score") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

# Display plot in Rmd
print(pathway_plot)

# Save high-resolution version
ggsave(paste0(OUTPUT_DIR, "PathwayActivity_Exhaustion.png"), plot = pathway_plot, width = PUB_WIDTH, height = PUB_HEIGHT, dpi = PUB_RES)

7. Pseudotime Trajectory Analysis (Monocle3 Template)

Purpose: Orders cells along a developmental or disease progression path (e.g., malignant T-cell evolution).

# NOTE: Requires monocle3 installation and the object to be converted.

# 1. Convert Seurat object to Monocle3 CellDataSet
# cds <- as.cell_data_set(ss.integrated)

# 2. Learn the graph structure
# cds <- learn_graph(cds, use_partition = TRUE)

# 3. Order cells by pseudotime (requires defining a root node)
# cds <- order_cells(cds, root_cells = colnames(ss.integrated)[ss.integrated$cell_type == "CD4 T cells"][1])

# 4. Plot the pseudotime trajectory
# trajectory_plot <- plot_cells(cds, color_cells_by = "pseudotime", label_groups_by_cluster = FALSE, cell_size = 0.5) +
#   ggtitle("Pseudotime Trajectory of T-cell Differentiation")

# ggsave(paste0(OUTPUT_DIR, "Monocle3_Pseudotime_Plot.png"), plot = trajectory_plot, width = PUB_WIDTH, height = PUB_HEIGHT, dpi = PUB_RES)

8. Intercellular Communication Network (CellChat Template)

Purpose: Summarizes the overall and specific ligand-receptor interactions between cell types in the TME.

# NOTE: Requires CellChat installation and a fully annotated object.

# 1. Prepare CellChat object
# cellchat <- createCellChat(object = ss.integrated, group.by = "cell_type", assay = "RNA")
# CellChatDB <- CellChatDB.human
# cellchat@DB <- CellChatDB

# 2. Compute and aggregate the network
# cellchat <- subsetData(cellchat)
# cellchat <- computeCommunProb(cellchat)
# cellchat <- aggregateNet(cellchat)

# 3. Visualize the overall network (Circle Plot)
# net_circle_plot <- netVisual_circle(cellchat@net$count,
#                                     vertex.weight = rowSums(cellchat@net$count),
#                                     weight.scale = TRUE,
#                                     label.edge = FALSE,
#                                     title.name = "Number of Interactions")

# ggsave(paste0(OUTPUT_DIR, "CellChat_Overall_Interactions.png"), plot = net_circle_plot, width = PUB_WIDTH, height = PUB_HEIGHT, dpi = PUB_RES)

# 4. Targeted Bubble Plot (e.g., for a key Sézary axis like CCL5-ACKR1)
# target_axis_plot <- netVisual_bubble(cellchat,
#                                      sources.use = c("CD4 T cells"),
#                                      targets.use = c("CD14 Monocytes"),
#                                      pairLR.use = data.frame(ligand = "CCL5", receptor = "ACKR1"))
# ggsave(paste0(OUTPUT_DIR, "CellChat_Targeted_Axis.png"), plot = target_axis_plot, width = PUB_WIDTH, height = PUB_HEIGHT, dpi = PUB_RES)

9. Clonal and Genetic Heterogeneity Visualizations (CNV/TCR Templates)

These are crucial for cancer-specific analysis (Herrera et al. and Su et al. methodologies).

# NOTE: These require specialized packages (copykat/InferCNV, scRepertoire)
# and often raw count data or external TCR data files.

# CNV Heatmap (for Malignant Cell Identification)
# CNV_heatmap_plot <- plot(copykat.obj$cnv.scores)
# ggsave(paste0(OUTPUT_DIR, "CNV_Inference_Heatmap.png"), plot = CNV_heatmap_plot, width = 12, height = 18, dpi = PUB_RES)

# TCR Clonal Heatmap (for Tracking Clonal Expansion)
# clonal_heatmap_plot <- clonal_heatmap(combined_tcr, cloneCall = "CTstrict", group.by = "cell_type")
# ggsave(paste0(OUTPUT_DIR, "TCR_Clonal_Heatmap.png"), plot = clonal_heatmap_plot, width = 10, height = 8, dpi = PUB_RES)

Conclusion and Next Steps

This R Markdown report demonstrates the core visualization pipeline for single-cell RNA-seq data in cancer genomics. By replacing the placeholder data loading with your Harmony-integrated object and installing the necessary advanced packages, you will be able to reproduce these publication-quality figures for your Sézary syndrome study. ```