Single-Cell RNA-Seq Analysis: Immune Compartment from CELLxGENE

Introduction

Background

This analysis examines a subset of publicly available single-cell RNA sequencing dataset from the CZ CELLxGENE Discover platform. The dataset (ID: 27c51769-800a-4a44-a885-74ce413b51f5) contains immune compartment cells, providing an opportunity to characterize immune cell heterogeneity at single-cell resolution.

Single-cell RNA sequencing enables the identification of distinct cell populations, rare cell types, and cell-type-specific gene expression patterns that are masked in bulk RNA-seq approaches.

Research Questions

The primary goals of this analysis are to:

1. Load and quality control the CELLxGENE H5AD dataset
2. Identify and characterize major immune cell types present in the dataset
3. Examine cell type composition across samples/conditions
4. Identify marker genes that define each cell population
5. Perform differential expression and pathway enrichment analyses

---

Methods

Data Acquisition and Loading

The dataset was downloaded from CELLxGENE Discover using the dataset ID. We load the H5AD file using Python’s anndata library through reticulate and convert it to a Seurat object.

# Define RDS file path
rds_file <- "immune_compartment_subset_seurat.rds"

# Check if RDS file exists
if (file.exists(rds_file)) {
  cat("Loading existing Seurat object from RDS file...\n")
  seurat_obj <- readRDS(rds_file)
  cat("Seurat object loaded successfully!\n")
  print(seurat_obj)
  original_cell_count <- ncol(seurat_obj)
  
}

## Loading existing Seurat object from RDS file...
## Seurat object loaded successfully!
## An object of class Seurat 
## 38333 features across 10000 samples within 1 assay 
## Active assay: RNA (38333 features, 0 variable features)
##  2 layers present: counts, data

# Set up Python and load anndata
use_virtualenv("r-reticulate", required = TRUE)
py_install("anndata", envname = "r-reticulate")
ad <- import("anndata", convert = FALSE)
py_install(c("scanpy", "h5py", "numpy", "pandas"), envname = "r-reticulate")

# Load the H5AD file
adata <- ad$read_h5ad("27c51769-800a-4a44-a885-74ce413b51f5.h5ad")

cat("H5AD file loaded successfully!\n\n")

# Display basic information
cat("Dataset dimensions:\n")
cat("Number of cells:", py_to_r(adata$n_obs), "\n")
cat("Number of genes:", py_to_r(adata$n_vars), "\n\n")

cat("Available cell metadata (obs) columns:\n")
obs_cols <- py_to_r(adata$obs$columns$to_list())
cat(paste(obs_cols, collapse = ", "), "\n\n")

cat("Available gene metadata (var) columns:\n")
var_cols <- py_to_r(adata$var$columns$to_list())
cat(paste(var_cols, collapse = ", "), "\n")

# SUBSET THE DATA 
# ============================================
# Set subset size (adjust as needed)
SUBSET_SIZE <- 10000  # Use 10,000 cells for faster analysis

# Get total number of cells
total_cells <- py_to_r(adata$n_obs)

if (total_cells > SUBSET_SIZE) {
  cat("Subsetting data to", SUBSET_SIZE, "cells for computational efficiency...\n")
  
  # Set seed for reproducibility
  set.seed(42)
  
  # If cell_type exists, do stratified sampling to maintain cell type diversity
  obs_cols <- py_to_r(adata$obs$columns$to_list())
  
  if ("cell_type" %in% obs_cols) {
    cat("Performing stratified sampling by cell_type...\n")
    
    # Get cell types
    cell_types <- py_to_r(adata$obs['cell_type'].to_numpy())

    
    # Calculate proportional sample size per cell type
    cell_type_counts <- table(cell_types)
    cell_type_props <- cell_type_counts / sum(cell_type_counts)
    
    # Sample proportionally from each cell type
    sampled_indices <- c()
    for (ct in names(cell_type_counts)) {
      ct_indices <- which(cell_types == ct)
      n_sample <- min(ceiling(SUBSET_SIZE * cell_type_props[ct]), length(ct_indices))
      sampled_indices <- c(sampled_indices, sample(ct_indices, n_sample))
    }
    
    # Trim to exact subset size if needed
    if (length(sampled_indices) > SUBSET_SIZE) {
      sampled_indices <- sample(sampled_indices, SUBSET_SIZE)
    }
    
    # Convert to Python indices (0-based)
    py_indices <- reticulate::np_array(sampled_indices - 1, dtype = "int32")
    
  } else {
    # Random sampling if no cell_type column
    cat("Performing random sampling...\n")
    sampled_indices <- sample(1:total_cells, SUBSET_SIZE)
    py_indices <- reticulate::np_array(sampled_indices - 1, dtype = "int32")
  }
  
  # Subset the AnnData object
  adata <- adata[py_indices,]
  
  cat("Subset complete!\n")
  cat("New dataset size:", py_to_r(adata$n_obs), "cells\n\n")
} else {
  cat("Dataset is small enough, using all cells.\n\n")
}

cat("Available cell metadata (obs) columns:\n")
obs_cols <- py_to_r(adata$obs$columns$to_list())
cat(paste(obs_cols, collapse = ", "), "\n\n")

cat("Available gene metadata (var) columns:\n")
var_cols <- py_to_r(adata$var$columns$to_list())
cat(paste(var_cols, collapse = ", "), "\n")

# Try to get raw counts first, then fall back to X
if (py_has_attr(adata, "raw") && !py_is_null_xptr(adata$raw)) {
  cat("Using raw.X for counts\n")
  counts_mat <- adata$raw$X
  gene_names <- py_to_r(adata$raw$var_names$to_list())
} else if ("counts" %in% py_to_r(adata$layers$keys())) {
  cat("Using layers['counts'] for counts\n")
  counts_mat <- adata$layers["counts"]
  gene_names <- py_to_r(adata$var_names$to_list())
} else {
  cat("Using X matrix for counts\n")
  counts_mat <- adata$X
  gene_names <- py_to_r(adata$var_names$to_list())
}

cell_names <- py_to_r(adata$obs_names$to_list())

# Convert to character vectors to ensure proper names
gene_names <- as.character(gene_names)
cell_names <- as.character(cell_names)
cat("Gene names length:", length(gene_names), "\n")
cat("Cell names length:", length(cell_names), "\n")

# Handle sparse matrix efficiently
cat("Converting sparse matrix to R sparse format...\n")

if (inherits(counts_mat, "scipy.sparse.csr.csr_matrix") || 
    inherits(counts_mat, "scipy.sparse.csc.csc_matrix")) {
  
  # Convert to CSC format if not already
  if (inherits(counts_mat, "scipy.sparse.csr.csr_matrix")) {
    counts_mat <- counts_mat$tocsc()
  }
  
  # Extract sparse matrix components
  i <- py_to_r(counts_mat$indices)
  p <- py_to_r(counts_mat$indptr)
  x <- py_to_r(counts_mat$data)
  dims <- py_to_r(counts_mat$shape)
  
 cat("Sparse matrix dimensions (cells x genes):", dims[1], "x", dims[2], "\n")
  
  # Create R sparse matrix (dgCMatrix) WITHOUT names first
  # scipy stores as cells x genes
  counts_temp <- sparseMatrix(
    i = i,
    p = p,
    x = x,
    dims = dims,
    index1 = FALSE,
    dimnames = NULL  # Don't set names during creation
  )
 # Transpose to genes x cells for Seurat
  counts <- t(counts_temp)
  
  # Now assign names AFTER transposition
  rownames(counts) <- gene_names
  colnames(counts) <- cell_names
  
  cat("Sparse matrix conversion complete!\n")
  cat("Final matrix dimensions:", nrow(counts), "genes x", ncol(counts), "cells\n")
  
  # Clean up
  rm(counts_temp)
  
} else {
  cat("Matrix is dense format\n")
  counts_dense <- py_to_r(counts_mat)
  counts <- t(counts_dense)
  rownames(counts) <- gene_names
  colnames(counts) <- cell_names
  rm(counts_dense)
}

# Extract metadata
cat("Extracting metadata...\n")
metadata <- py_to_r(adata$obs)
rownames(metadata) <- cell_names

# Create Seurat object
cat("\nCreating Seurat object...\n")
seurat_obj <- CreateSeuratObject(
  counts = counts,
  meta.data = metadata,
  project = "Immune_Compartment_Subset",
  min.cells = 0,
  min.features = 0
)

cat("\nSeurat object created successfully!\n")
print(seurat_obj)

# Save original cell count
original_cell_count <- ncol(seurat_obj)

# Clear large objects to free memory
rm(counts_mat, counts, i, p, x, adata)
gc()

# ANNOTATE GENES WITH HUGO SYMBOLS
  cat("\n=== Annotating Ensembl IDs with Gene Symbols ===\n")
  
  # Get current gene names (Ensembl IDs)
  ensembl_ids <- rownames(seurat_obj)
  cat("Total genes in dataset:", length(ensembl_ids), "\n")
  
  # Check if they look like Ensembl IDs
  cat("First 5 gene IDs:", head(ensembl_ids, 5), "\n")
  
  # Connect to Ensembl BioMart
  cat("Connecting to Ensembl BioMart (human)...\n")
  ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl")
  
  # Query BioMart for gene symbols
  cat("Querying gene symbols...\n")
  gene_annotations <- getBM(
    attributes = c('ensembl_gene_id', 'external_gene_name', 'gene_biotype'),
    filters = 'ensembl_gene_id',
    values = ensembl_ids,
    mart = ensembl
  )
  
  cat("Retrieved annotations for", nrow(gene_annotations), "genes\n")
  
  # Create a mapping dataframe
  gene_map <- data.frame(
    ensembl_id = ensembl_ids,
    stringsAsFactors = FALSE
  )
  
  # Merge with annotations
  gene_map <- left_join(gene_map, gene_annotations, 
                        by = c("ensembl_id" = "ensembl_gene_id"))
  
  # For genes without symbols, use Ensembl ID
  gene_map$gene_symbol <- ifelse(
    is.na(gene_map$external_gene_name) | gene_map$external_gene_name == "",
    gene_map$ensembl_id,
    gene_map$external_gene_name
  )
  
  # Handle duplicates by making gene symbols unique
  gene_map$gene_symbol_unique <- make.unique(gene_map$gene_symbol)
  
  cat("\nAnnotation summary:\n")
  cat("Genes with symbols:", sum(!is.na(gene_map$external_gene_name) & 
                                   gene_map$external_gene_name != ""), "\n")
  cat("Genes without symbols (kept as Ensembl ID):", 
      sum(is.na(gene_map$external_gene_name) | gene_map$external_gene_name == ""), "\n")
  
  # Add gene annotations to Seurat object metadata
  # Add feature-level metadata to the RNA assay
  # Ensure rownames of gene_map match the assay
rownames(gene_map) <- gene_map$ensembl_id
seurat_obj[["RNA"]] <- AddMetaData(
  object = seurat_obj[["RNA"]],
  metadata = gene_map %>%
    dplyr::select(ensembl_id, gene_symbol = gene_symbol_unique, gene_biotype)
)

  # Update rownames to use gene symbols
  cat("\nUpdating gene names in Seurat object...\n")
  
  # Get all assay data
  counts_data <- GetAssayData(seurat_obj, slot = "counts")
  rownames(counts_data) <- gene_map$gene_symbol_unique
  
  # Create new assay with updated names
  seurat_obj[["RNA"]] <- CreateAssayObject(counts = counts_data)
  
  # Re-add metadata
  seurat_obj@meta.data <- metadata
  
  cat("Gene annotation complete!\n")
  cat("First 10 gene names:", head(rownames(seurat_obj), 10), "\n\n")
  
  # SAVE SEURAT OBJECT TO RDS
  # ============================================
  cat("Saving Seurat object to RDS file:", rds_file, "\n")
  saveRDS(seurat_obj, file = rds_file)
  cat("RDS file saved successfully!\n")
  cat("Next time, the analysis will load this file directly.\n\n")
  
  print(seurat_obj)

---

Results

Gene Annotation Summary

# Display gene annotation information
# Get current gene names (Ensembl IDs)
  ensembl_ids <- rownames(seurat_obj)
  cat("Total genes in dataset:", length(ensembl_ids), "\n")

## Total genes in dataset: 38333

  # Check if they look like Ensembl IDs
  cat("First 5 gene IDs:", head(ensembl_ids, 5), "\n")

## First 5 gene IDs: ENSG00000286448 LINC00115 FAM41C SAMD11 NOC2L

  # Connect to Ensembl BioMart
  cat("Connecting to Ensembl BioMart (human)...\n")

## Connecting to Ensembl BioMart (human)...

  ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl")
  
  # Query BioMart for gene symbols
  cat("Querying gene symbols...\n")

## Querying gene symbols...

  gene_annotations <- getBM(
    attributes = c('ensembl_gene_id', 'external_gene_name', 'gene_biotype'),
    filters = 'ensembl_gene_id',
    values = ensembl_ids,
    mart = ensembl
  )
  
  cat("Retrieved annotations for", nrow(gene_annotations), "genes\n")

## Retrieved annotations for 11029 genes

  # Create a mapping dataframe
  gene_map <- data.frame(
    ensembl_id = ensembl_ids,
    stringsAsFactors = FALSE
  )
  
  # Merge with annotations
  gene_map <- left_join(gene_map, gene_annotations, 
                        by = c("ensembl_id" = "ensembl_gene_id"))
  
  # For genes without symbols, use Ensembl ID
  gene_map$gene_symbol <- ifelse(
    is.na(gene_map$external_gene_name) | gene_map$external_gene_name == "",
    gene_map$ensembl_id,
    gene_map$external_gene_name
  )
  
  # Handle duplicates by making gene symbols unique
  gene_map$gene_symbol_unique <- make.unique(gene_map$gene_symbol)
  
  gene_info <- gene_map
head(gene_info)

##        ensembl_id external_gene_name gene_biotype     gene_symbol
## 1 ENSG00000286448                 NA       lncRNA ENSG00000286448
## 2       LINC00115                 NA         <NA>       LINC00115
## 3          FAM41C                 NA         <NA>          FAM41C
## 4          SAMD11                 NA         <NA>          SAMD11
## 5           NOC2L                 NA         <NA>           NOC2L
## 6          KLHL17                 NA         <NA>          KLHL17
##   gene_symbol_unique
## 1    ENSG00000286448
## 2          LINC00115
## 3             FAM41C
## 4             SAMD11
## 5              NOC2L
## 6             KLHL17

if ("ensembl_id" %in% colnames(gene_info)) {
  
  cat("Gene Annotation Summary:\n")
  cat("Total genes:", nrow(gene_info), "\n")
  cat("Genes with HUGO symbols:", sum(gene_info$ensembl_id != gene_info$gene_symbol), "\n")
  cat("Genes kept as Ensembl ID:", sum(gene_info$ensembl_id == gene_info$gene_symbol), "\n\n")
  
  # Show biotype distribution if available
  if ("gene_biotype" %in% colnames(gene_info)) {
    cat("Gene biotype distribution:\n")
    print(table(gene_info$gene_biotype))
    
    # Visualize biotype distribution
    biotype_counts <- as.data.frame(table(gene_info$gene_biotype))
    colnames(biotype_counts) <- c("Biotype", "Count")
    biotype_counts <- biotype_counts %>% arrange(desc(Count)) %>% head(10)
    
    ggplot(biotype_counts, aes(x = reorder(Biotype, Count), y = Count)) +
      geom_bar(stat = "identity", fill = "steelblue") +
      coord_flip() +
      labs(title = "Top 10 Gene Biotypes in Dataset",
           x = "Gene Biotype", y = "Number of Genes") +
      theme_minimal()
  }
  
  # Display example genes
  cat("\nExample annotated genes:\n")
  gene_info %>%
    dplyr::select(ensembl_id, gene_symbol, gene_biotype) %>%
    head(20) %>%
    kable() %>%
    kable_styling(bootstrap_options = c("striped", "hover"))
}

## Gene Annotation Summary:
## Total genes: 38333 
## Genes with HUGO symbols: 0 
## Genes kept as Ensembl ID: 38333 
## 
## Gene biotype distribution:
## 
##                           artifact                          IG_V_gene 
##                                 17                                  2 
##                    IG_V_pseudogene                             lncRNA 
##                                  3                               9613 
##                           misc_RNA               processed_pseudogene 
##                                  9                                 54 
##                     protein_coding                         pseudogene 
##                                429                                  1 
##                             scaRNA                             snoRNA 
##                                  1                                 12 
##                                TEC   transcribed_processed_pseudogene 
##                                841                                 13 
##     transcribed_unitary_pseudogene transcribed_unprocessed_pseudogene 
##                                  3                                 19 
##                 unitary_pseudogene             unprocessed_pseudogene 
##                                  1                                 11 
## 
## Example annotated genes:

ensembl_id	gene_symbol	gene_biotype
ENSG00000286448	ENSG00000286448	lncRNA
LINC00115	LINC00115	NA
FAM41C	FAM41C	NA
SAMD11	SAMD11	NA
NOC2L	NOC2L	NA
KLHL17	KLHL17	NA
PLEKHN1	PLEKHN1	NA
PERM1	PERM1	NA
HES4	HES4	NA
ISG15	ISG15	NA
AGRN	AGRN	NA
RNF223	RNF223	NA
C1orf159	C1orf159	NA
TTLL10	TTLL10	NA
TNFRSF18	TNFRSF18	NA
TNFRSF4	TNFRSF4	NA
SDF4	SDF4	NA
B3GALT6	B3GALT6	NA
UBE2J2	UBE2J2	NA
SCNN1D	SCNN1D	NA

Quality Control

# Calculate QC metrics if not already present
if (!"nCount_RNA" %in% colnames(seurat_obj@meta.data)) {
  seurat_obj[["nCount_RNA"]] <- colSums(seurat_obj@assays$RNA@counts)
  seurat_obj[["nFeature_RNA"]] <- colSums(seurat_obj@assays$RNA@counts > 0)
}

# Calculate mitochondrial percentage
seurat_obj[["percent.mt"]] <- PercentageFeatureSet(seurat_obj, pattern = "^MT-")

# Calculate ribosomal percentage
seurat_obj[["percent.rb"]] <- PercentageFeatureSet(seurat_obj, pattern = "^RP[SL]")

# Visualize QC metrics
p1 <- VlnPlot(seurat_obj, 
              features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), 
              ncol = 3, 
              pt.size = 0.1) & 
  theme(axis.title.x = element_blank(), 
        plot.title = element_text(size = 10))

p1

hist(seurat_obj@meta.data$nFeature_RNA, breaks = 100)

plot1 <- FeatureScatter(seurat_obj, feature1 = "nCount_RNA", feature2 = "percent.mt") +
  geom_hline(yintercept = 20, linetype = "dashed", color = "red", linewidth = 1) +
  ggtitle("QC: UMI vs Mitochondrial %") +
  theme_minimal()

plot2 <- FeatureScatter(seurat_obj, feature1 = "nCount_RNA", feature2 = "nFeature_RNA") +
  ggtitle("QC: UMI vs Gene Count") +
  theme_minimal()

plot1 + plot2

qc_stats <- data.frame(
  Metric = c("Total Cells (Subset)", 
             "Median Genes/Cell", 
             "Median UMIs/Cell", 
             "Median % MT",
             "Mean % MT"),
  Value = c(ncol(seurat_obj),
            round(median(seurat_obj$nFeature_RNA)),
            round(median(seurat_obj$nCount_RNA)),
            round(median(seurat_obj$percent.mt), 2),
            round(mean(seurat_obj$percent.mt), 2))
)

kable(qc_stats, 
      caption = "Quality Control Summary Statistics (Subset)",
      align = 'lr') %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), 
                full_width = FALSE,
                position = "left") %>%
  row_spec(0, bold = TRUE, color = "white", background = "#4CAF50")

Quality Control Summary Statistics (Subset)

Metric	Value
Total Cells (Subset)	10000.00
Median Genes/Cell	1004.00
Median UMIs/Cell	2413.00
Median % MT	4.15
Mean % MT	5.17

Explore Existing Metadata

cat("Available metadata columns:\n")

## Available metadata columns:

print(colnames(seurat_obj@meta.data))

##  [1] "tissue_ontology_term_id"                 
##  [2] "tissue_type"                             
##  [3] "assay_ontology_term_id"                  
##  [4] "disease_ontology_term_id"                
##  [5] "cell_type_ontology_term_id"              
##  [6] "self_reported_ethnicity_ontology_term_id"
##  [7] "development_stage_ontology_term_id"      
##  [8] "sex_ontology_term_id"                    
##  [9] "donor_id"                                
## [10] "suspension_type"                         
## [11] "grade"                                   
## [12] "author_cell_type"                        
## [13] "batch"                                   
## [14] "is_primary_data"                         
## [15] "cell_type"                               
## [16] "assay"                                   
## [17] "disease"                                 
## [18] "sex"                                     
## [19] "tissue"                                  
## [20] "self_reported_ethnicity"                 
## [21] "development_stage"                       
## [22] "observation_joinid"                      
## [23] "nCount_RNA"                              
## [24] "nFeature_RNA"                            
## [25] "percent.mt"                              
## [26] "percent.rb"

if ("cell_type" %in% colnames(seurat_obj@meta.data)) {
  cat("\n=== Cell Type Distribution (in subset) ===\n")
  cell_type_table <- sort(table(seurat_obj$cell_type), decreasing = TRUE)
  print(head(cell_type_table, 20))
}

## 
## === Cell Type Distribution (in subset) ===
## 
##                                                               macrophage 
##                                                                     1665 
##                          effector memory CD8-positive, alpha-beta T cell 
##                                                                     1050 
##                          effector memory CD4-positive, alpha-beta T cell 
##                                                                      894 
##                                                      natural killer cell 
##                                                                      825 
##                                                            memory B cell 
##                                                                      761 
## CD4-positive, CD25-positive, CCR4-positive, alpha-beta regulatory T cell 
##                                                                      712 
##                                                             naive T cell 
##                                                                      662 
##                                                          IgG plasma cell 
##                                                                      658 
##                                                         exhausted T cell 
##                                                                      647 
##                               CD8-positive, alpha-beta regulatory T cell 
##                                                                      503 
##                                              conventional dendritic cell 
##                                                                      450 
##                                                 T follicular helper cell 
##                                                                      431 
##                                                           cycling T cell 
##                                                                      238 
##                                                       cycling macrophage 
##                                                                      180 
##                                                                mast cell 
##                                                                      166 
##                                              plasmacytoid dendritic cell 
##                                                                       92 
##                                                   myeloid dendritic cell 
##                                                                       50 
##                                               CD4-positive helper T cell 
##                                                                       16

if ("tissue" %in% colnames(seurat_obj@meta.data)) {
  cat("\n=== Tissue Distribution ===\n")
  print(table(seurat_obj$tissue))
}

## 
## === Tissue Distribution ===
## 
## breast 
##  10000

if ("disease" %in% colnames(seurat_obj@meta.data)) {
  cat("\n=== Disease Status ===\n")
  print(table(seurat_obj$disease))
}

## 
## === Disease Status ===
## 
##                                              breast mucinous carcinoma 
##                                                                    104 
##                                              breast apocrine carcinoma 
##                                                                     82 
## invasive tubular breast carcinoma || invasive lobular breast carcinoma 
##                                                                    153 
##                                       invasive ductal breast carcinoma 
##                                                                   4429 
##                                                       breast carcinoma 
##                                                                    206 
##                                      invasive lobular breast carcinoma 
##                                                                     87 
##                                       triple-negative breast carcinoma 
##                                                                    685 
##                                           metaplastic breast carcinoma 
##                                                                     13 
##                                         HER2 positive breast carcinoma 
##                                                                     94 
##                               estrogen-receptor positive breast cancer 
##                                                                    576 
##                                                          breast cancer 
##                                                                   3571

if ("donor_id" %in% colnames(seurat_obj@meta.data)) {
  cat("\n=== Number of Donors ===\n")
  cat(length(unique(seurat_obj$donor_id)), "unique donors\n")
}

## 
## === Number of Donors ===
## 136 unique donors

Filter cells based on QC metrics

seurat_obj <- subset(
  seurat_obj,
  subset = nFeature_RNA >= 200 &
           nFeature_RNA <= 5000 &
           percent.mt <= 5 &
           nCount_RNA >= 500 &
           nCount_RNA <= 30000
)

Normalization and Feature Selection

seurat_obj <- NormalizeData(seurat_obj, 
                            normalization.method = "LogNormalize", 
                            scale.factor = 10000)

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

top10 <- head(VariableFeatures(seurat_obj), 10)

cat("Top 10 most variable genes:\n")

## Top 10 most variable genes:

cat(paste(top10, collapse = ", "), "\n\n")

## IGHGP, IGLC3, IGLC2, IGHA1, IGHM, ZNF71, TPSAB1, SCRT2, PTGDS, TPSB2

plot1 <- VariableFeaturePlot(seurat_obj) + 
  theme_minimal() +
  ggtitle("Variable Feature Selection")

plot2 <- LabelPoints(plot = plot1, 
                     points = top10, 
                     repel = TRUE)
plot2

# Scale only variable features
seurat_obj <- ScaleData(seurat_obj, features = VariableFeatures(seurat_obj))
cat("Data normalized and scaled.\n")

## Data normalized and scaled.

Dimensionality Reduction

seurat_obj <- RunPCA(seurat_obj, 
                     features = VariableFeatures(object = seurat_obj), 
                     verbose = FALSE)

print(seurat_obj[["pca"]], dims = 1:5, nfeatures = 5)

## PC_ 1 
## Positive:  FTL, SSR4, PSAP, GADD45B, CTSD 
## Negative:  NPIPB9, TKTL1, PLA2G3, SCNN1A, NRG3 
## PC_ 2 
## Positive:  FAM9C, CORO2B, APC2, PLEKHG6, ARHGAP39 
## Negative:  POF1B, CT45A1, OR10Q1, PASD1, PMP2 
## PC_ 3 
## Positive:  CST3, CD68, CD14, C1QA, FCGR2A 
## Negative:  GZMA, FKBP11, KLRB1, MZB1, PIM2 
## PC_ 4 
## Positive:  HM13-AS1, SCRT2, ZNF71, FBLN1, INPP5J 
## Negative:  ENSG00000229986, SMIM17, NOG, SLC2A10, FDXACB1 
## PC_ 5 
## Positive:  LGALS1, CSTB, NEAT1, TYMP, MIF 
## Negative:  TAMALIN, CD83, NR4A3, DUSP5, BIRC3

VizDimLoadings(seurat_obj, dims = 1:2, reduction = "pca")

DimPlot(seurat_obj, reduction = "pca") + NoLegend()

DimHeatmap(seurat_obj, dims = 1, cells = 500, balanced = TRUE)

ElbowPlot(seurat_obj, ndims = 30) +
  geom_vline(xintercept = 15, linetype = "dashed", color = "red", linewidth = 1) +
  ggtitle("Elbow Plot - Selecting Number of PCs") +
  theme_minimal()

Clustering and UMAP

n_pcs <- 15  # Reduced for subset

seurat_obj <- FindNeighbors(seurat_obj, dims = 1:n_pcs)
seurat_obj <- FindClusters(seurat_obj, resolution = 0.5, verbose = FALSE)
seurat_obj <- RunUMAP(seurat_obj, dims = 1:n_pcs, verbose = FALSE)

n_clusters <- length(unique(Idents(seurat_obj)))

DimPlot(seurat_obj, 
        reduction = "umap", 
        label = TRUE, 
        label.size = 6,
        pt.size = 0.8) +
  ggtitle(paste0("UMAP - ", n_clusters, "UMAP: Cluster Assignments ")) +
  theme_minimal()

if ("cell_type" %in% colnames(seurat_obj@meta.data)) {
  DimPlot(seurat_obj, 
          reduction = "umap", 
          group.by = "cell_type",
          label = TRUE,
          repel = TRUE,
          pt.size = 0.5) +
    ggtitle("UMAP - Cell Type Annotations") +
    theme_minimal() +
    theme(legend.position = "right")
}

Cell Type Analysis

# Cluster composition
cluster_table <- as.data.frame(table(Idents(seurat_obj)))
colnames(cluster_table) <- c("Cluster", "Count")
cluster_table$Percentage <- round(cluster_table$Count / sum(cluster_table$Count) * 100, 2)

ggplot(cluster_table, aes(x = Cluster, y = Count, fill = Cluster)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = paste0(Count, "\n(", Percentage, "%)")), 
            vjust = -0.5, size = 3) +
  labs(title = "Cells per Cluster", x = "Cluster", y = "Number of Cells") +
  theme_minimal() +
  theme(
    legend.position = "none",
    axis.text.x = element_text(angle = 45, hjust = 1)  # slant x labels
  )

# If cell_type exists, show composition
if ("cell_type" %in% colnames(seurat_obj@meta.data)) {
  Idents(seurat_obj) <- seurat_obj$cell_type
  
  cell_counts <- table(seurat_obj$cell_type)
  cell_props <- prop.table(cell_counts) * 100
  
  prop_df <- data.frame(
    Cell_Type = names(cell_props),
    Count = as.numeric(cell_counts),
    Percentage = round(as.numeric(cell_props), 2)
  ) %>%
    arrange(desc(Percentage))
  
  kable(prop_df, 
        caption = "Cell Type Distribution (Subset)",
        col.names = c("Cell Type", "Cell Count", "Percentage (%)"),
        align = 'lrr') %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed"), 
                  full_width = FALSE) %>%
    row_spec(0, bold = TRUE, color = "white", background = "#2196F3")
  
  ggplot(prop_df %>% head(15), 
         aes(x = reorder(Cell_Type, Percentage), y = Percentage, fill = Cell_Type)) +
    geom_bar(stat = "identity", color = "black", size = 0.3) +
    coord_flip() +
    theme_minimal() +
    theme(legend.position = "none") +
    labs(title = "Cell Type Distribution",
         x = "Cell Type",
         y = "Percentage (%)") +
    geom_text(aes(label = paste0(Percentage, "%")), hjust = -0.1, size = 3)
}

---

Differential Expression Analysis

Find Cluster Markers

# Find markers for all clusters
cat("Finding differentially expressed genes for all clusters...\n")

## Finding differentially expressed genes for all clusters...

cat("This may take a few minutes...\n\n")

## This may take a few minutes...

all.markers <- FindAllMarkers(
  seurat_obj, 
  only.pos = TRUE, 
  min.pct = 0.25, 
  logfc.threshold = 0.25
)

cat("Found", nrow(all.markers), "significant markers across all clusters\n\n")

## Found 12688 significant markers across all clusters

# Filter for highly significant markers
significant.markers <- all.markers %>%
  filter(p_val_adj < 0.05, avg_log2FC > 0.5)

cat("Highly significant markers (p_adj < 0.05, log2FC > 0.5):", 
    nrow(significant.markers), "\n")

## Highly significant markers (p_adj < 0.05, log2FC > 0.5): 7638

# Show top markers per cluster
top_markers_per_cluster <- all.markers %>%
  group_by(cluster) %>%
  arrange(desc(avg_log2FC)) %>%
  slice_head(n = 5)

cat("\nTop 5 markers per cluster:\n")

## 
## Top 5 markers per cluster:

top_markers_per_cluster %>%
  dplyr::select(cluster, gene, avg_log2FC, pct.1, pct.2, p_val_adj) %>%
  kable(digits = 3, caption = "Top Markers by Cluster") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

Top Markers by Cluster

cluster	gene	avg_log2FC	pct.1	pct.2	p_val_adj
mast cell	CPA3	11.320	0.965	0.003	0
mast cell	TPSAB1	11.124	0.991	0.019	0
mast cell	TPSB2	11.102	0.929	0.014	0
mast cell	GATA2	10.189	0.858	0.005	0
mast cell	SLC18A2	10.083	0.726	0.003	0
macrophage	RNASE1	6.142	0.561	0.029	0
macrophage	GPNMB	5.918	0.599	0.026	0
macrophage	APOE	5.892	0.802	0.142	0
macrophage	APOC1	5.858	0.749	0.074	0
macrophage	SDC3	5.768	0.316	0.007	0
CD4-positive helper T cell	SPON1	5.142	0.250	0.006	0
CD4-positive helper T cell	TRAV8-6	5.067	0.250	0.007	0
CD4-positive helper T cell	TRAV8-2	4.887	0.250	0.015	0
CD4-positive helper T cell	TRAV9-2	4.746	0.333	0.010	0
CD4-positive helper T cell	TRBV2	4.735	0.250	0.009	0
natural killer cell	KLRF1	5.630	0.275	0.009	0
natural killer cell	TRDC	4.498	0.475	0.025	0
natural killer cell	KLRC1	4.037	0.418	0.032	0
natural killer cell	XCL1	3.831	0.545	0.053	0
natural killer cell	XCL2	3.738	0.577	0.064	0
myeloid dendritic cell	AOC1	8.664	0.311	0.003	0
myeloid dendritic cell	LAD1	8.288	0.733	0.005	0
myeloid dendritic cell	CCL19	7.892	0.778	0.023	0
myeloid dendritic cell	TREML1	7.498	0.311	0.003	0
myeloid dendritic cell	SLCO5A1	7.436	0.689	0.007	0
plasmacytoid dendritic cell	CLEC4C	10.172	0.486	0.001	0
plasmacytoid dendritic cell	PTCRA	8.663	0.689	0.006	0
plasmacytoid dendritic cell	LILRA4	8.462	0.689	0.008	0
plasmacytoid dendritic cell	PTGDS	8.319	0.527	0.014	0
plasmacytoid dendritic cell	VASH2	8.176	0.284	0.001	0
memory B cell	BANK1	7.935	0.722	0.014	0
memory B cell	FCRLA	7.171	0.278	0.005	0
memory B cell	MS4A1	7.022	0.759	0.020	0
memory B cell	CD19	6.406	0.253	0.010	0
memory B cell	BLK	6.257	0.370	0.010	0
CD8-positive, alpha-beta regulatory T cell	CMPK2	2.569	0.303	0.107	0
CD8-positive, alpha-beta regulatory T cell	IFIT3	2.357	0.420	0.165	0
CD8-positive, alpha-beta regulatory T cell	HERC5	2.311	0.420	0.170	0
CD8-positive, alpha-beta regulatory T cell	MX1	2.297	0.353	0.151	0
CD8-positive, alpha-beta regulatory T cell	RSAD2	2.215	0.399	0.141	0
naive T cell	TCF7	1.727	0.250	0.115	0
naive T cell	RPL3	1.661	0.785	0.516	0
naive T cell	IL7R	1.605	0.771	0.342	0
naive T cell	SELL	1.590	0.323	0.183	0
naive T cell	LDHB	1.501	0.611	0.443	0
effector memory CD4-positive, alpha-beta T cell	ANXA1	1.879	0.777	0.459	0
effector memory CD4-positive, alpha-beta T cell	RGCC	1.848	0.719	0.427	0
effector memory CD4-positive, alpha-beta T cell	IL7R	1.797	0.817	0.322	0
effector memory CD4-positive, alpha-beta T cell	KLRB1	1.772	0.334	0.131	0
effector memory CD4-positive, alpha-beta T cell	TC2N	1.761	0.334	0.091	0
effector memory CD8-positive, alpha-beta T cell	TRGC2	2.617	0.260	0.052	0
effector memory CD8-positive, alpha-beta T cell	GZMK	2.245	0.480	0.124	0
effector memory CD8-positive, alpha-beta T cell	TUBA4A	1.965	0.572	0.327	0
effector memory CD8-positive, alpha-beta T cell	SAMD3	1.953	0.273	0.103	0
effector memory CD8-positive, alpha-beta T cell	GZMH	1.947	0.293	0.093	0
IgG plasma cell	IGLV6-57	9.789	0.259	0.006	0
IgG plasma cell	IGLC3	9.539	0.456	0.031	0
IgG plasma cell	IGLV3-1	9.343	0.346	0.009	0
IgG plasma cell	IGLL5	9.300	0.311	0.005	0
IgG plasma cell	IGKV3-20	9.180	0.253	0.028	0
conventional dendritic cell	CD1C	5.513	0.379	0.014	0
conventional dendritic cell	CD1E	5.394	0.388	0.013	0
conventional dendritic cell	FCER1A	4.921	0.415	0.022	0
conventional dendritic cell	CFP	4.600	0.403	0.023	0
conventional dendritic cell	FCN1	4.413	0.482	0.046	0
CD4-positive, CD25-positive, CCR4-positive, alpha-beta regulatory T cell	FOXP3	6.206	0.428	0.012	0
CD4-positive, CD25-positive, CCR4-positive, alpha-beta regulatory T cell	RTKN2	5.765	0.351	0.014	0
CD4-positive, CD25-positive, CCR4-positive, alpha-beta regulatory T cell	IL2RA	4.556	0.364	0.039	0
CD4-positive, CD25-positive, CCR4-positive, alpha-beta regulatory T cell	TNFRSF4	4.321	0.579	0.111	0
CD4-positive, CD25-positive, CCR4-positive, alpha-beta regulatory T cell	LAYN	4.142	0.266	0.023	0
T follicular helper cell	CXCL13	4.610	0.820	0.075	0
T follicular helper cell	CHN1	4.341	0.322	0.024	0
T follicular helper cell	CD200	3.950	0.306	0.029	0
T follicular helper cell	NMB	3.623	0.278	0.064	0
T follicular helper cell	TSHZ2	3.171	0.325	0.044	0
exhausted T cell	ENSG00000237342	9.639	0.250	0.000	0
exhausted T cell	NPAS4	8.301	0.277	0.001	0
exhausted T cell	LINC00052	8.231	0.263	0.002	0
exhausted T cell	PLEKHG6	8.181	0.275	0.003	0
exhausted T cell	FAM9C	8.123	0.279	0.002	0
cycling T cell	RRM2	6.203	0.482	0.016	0
cycling T cell	ASPM	5.762	0.376	0.015	0
cycling T cell	MKI67	5.669	0.525	0.022	0
cycling T cell	TYMS	5.363	0.546	0.029	0
cycling T cell	AURKB	5.190	0.355	0.014	0
cycling macrophage	CD1A	5.327	0.298	0.010	0
cycling macrophage	PCLAF	4.920	0.617	0.018	0
cycling macrophage	TROAP	4.899	0.298	0.006	0
cycling macrophage	ESCO2	4.689	0.362	0.010	0
cycling macrophage	FOXM1	4.553	0.277	0.008	0

Marker Gene Heatmap

# Select top markers for heatmap
top10_markers <- all.markers %>%
  filter(p_val_adj < 0.05) %>%
  group_by(cluster) %>%
  arrange(desc(avg_log2FC)) %>%
  slice_head(n = 10) %>%
  ungroup()

if (nrow(top10_markers) > 0) {
  DoHeatmap(seurat_obj, features = top10_markers$gene) + 
    ggtitle("Top 10 Marker Genes per Cluster") +
    theme(axis.text.y = element_text(size = 8))
} else {
  cat("No significant markers found for heatmap\n")
}

Violin Plots for Key Markers

# Select top 2 markers per cluster for violin plots
top_genes_violin <- all.markers %>%
  filter(p_val_adj < 0.05) %>%
  group_by(cluster) %>%
  arrange(desc(avg_log2FC)) %>%
  slice_head(n = 2) %>%
  pull(gene) %>%
  unique() %>%
  head(12)  # Limit to 12 for visualization

if (length(top_genes_violin) > 0) {
  VlnPlot(seurat_obj, features = top_genes_violin, ncol = 4, pt.size = 0)
} else {
  cat("No markers available for violin plots\n")
}

Feature Plots for Top Markers

# Feature plots showing spatial distribution of top markers
top_genes_feature <- all.markers %>%
  filter(p_val_adj < 0.05) %>%
  group_by(cluster) %>%
  arrange(desc(avg_log2FC)) %>%
  slice_head(n = 1) %>%
  pull(gene) %>%
  head(9)  # Top marker from each cluster, max 9

if (length(top_genes_feature) > 0) {
  FeaturePlot(seurat_obj, features = top_genes_feature, ncol = 3)
} else {
  cat("No markers available for feature plots\n")
}

Dot Plot for Marker Genes

# Dot plot showing marker expression across clusters
top_markers_dotplot <- all.markers %>%
  filter(p_val_adj < 0.05) %>%
  group_by(cluster) %>%
  arrange(desc(avg_log2FC)) %>%
  slice_head(n = 3) %>%
  pull(gene)

if (length(top_markers_dotplot) > 0) {
  DotPlot(seurat_obj, features = unique(top_markers_dotplot)) + 
    coord_flip() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
    ggtitle("Marker Gene Expression Across Clusters")
} else {
  cat("No markers available for dot plot\n")
}

Known Immune Markers

immune_markers <- c(
  "CD3D", "CD3E", "CD4", "IL7R", "CD8A", "CD8B",
  "NKG7", "GNLY", "MS4A1", "CD79A", "CD14", "LYZ", "MZB1", "JCHAIN", "IGKC", "FCER1A", "CST3", "CLEC9A", "CD19", "NCAM1", "CCR7", "CD3G", "FCGR3A")

immune_markers_present <- immune_markers[immune_markers %in% rownames(seurat_obj)]

if (length(immune_markers_present) > 0) {
  FeaturePlot(seurat_obj, 
              features = immune_markers_present[1:min(9, length(immune_markers_present))], 
              ncol = 3,
              pt.size = 0.5) &
    theme_minimal()
}

Pairwise Cluster Comparisons

# Example: Compare two specific clusters if they exist
clusters <- levels(Idents(seurat_obj))

if (length(clusters) >= 2) {
  cat("Comparing cluster", clusters[1], "vs cluster", clusters[2], "\n")
  
  cluster_comparison <- FindMarkers(
    seurat_obj, 
    ident.1 = clusters[1], 
    ident.2 = clusters[2],
    min.pct = 0.25
  )
  
  # Add gene names as column
  cluster_comparison$gene <- rownames(cluster_comparison)
  
  # Show top upregulated in cluster 1
  cat("\nTop genes upregulated in cluster", clusters[1], ":\n")
  cluster_comparison %>%
    filter(avg_log2FC > 0) %>%
    arrange(p_val_adj) %>%
    head(10) %>%
    dplyr::select(gene, avg_log2FC, pct.1, pct.2, p_val_adj) %>%
    kable(digits = 3) %>%
    kable_styling(bootstrap_options = c("striped", "hover"))
  
  # Show top upregulated in cluster 2
  cat("\nTop genes upregulated in cluster", clusters[2], ":\n")
  cluster_comparison %>%
    filter(avg_log2FC < 0) %>%
    arrange(p_val_adj) %>%
    head(10) %>%
    dplyr::select(gene, avg_log2FC, pct.1, pct.2, p_val_adj) %>%
    kable(digits = 3) %>%
    kable_styling(bootstrap_options = c("striped", "hover"))
}

## Comparing cluster mast cell vs cluster macrophage 
## 
## Top genes upregulated in cluster mast cell :
## 
## Top genes upregulated in cluster macrophage :

	gene	avg_log2FC	pct.1	pct.2	p_val_adj
CD74	CD74	-3.948	0.531	0.975	0
HLA-DRA	HLA-DRA	-4.458	0.442	0.959	0
FTL	FTL	-3.083	0.982	0.987	0
AIF1	AIF1	-5.411	0.097	0.874	0
HLA-DRB1	HLA-DRB1	-3.746	0.451	0.917	0
TYROBP	TYROBP	-2.231	0.717	0.939	0
PSAP	PSAP	-2.843	0.584	0.916	0
HLA-DPA1	HLA-DPA1	-3.681	0.389	0.891	0
HLA-DPB1	HLA-DPB1	-3.482	0.363	0.898	0
ANXA2	ANXA2	-4.135	0.115	0.811	0

Marker Gene Summary Table

# Create comprehensive marker table
marker_summary <- all.markers %>%
  group_by(cluster) %>%
  summarise(
    n_markers = n(),
    n_significant = sum(p_val_adj < 0.05),
    top_gene = gene[which.max(avg_log2FC)],
    max_log2FC = max(avg_log2FC)
  )

marker_summary %>%
  kable(caption = "Marker Gene Summary by Cluster", digits = 2) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

Marker Gene Summary by Cluster

cluster	n_markers	n_significant	top_gene	max_log2FC
mast cell	444	259	CPA3	11.32
macrophage	1227	1224	RNASE1	6.14
CD4-positive helper T cell	801	130	SPON1	5.14
natural killer cell	370	225	KLRF1	5.63
myeloid dendritic cell	1083	569	AOC1	8.66
plasmacytoid dendritic cell	988	539	CLEC4C	10.17
memory B cell	356	255	BANK1	7.94
CD8-positive, alpha-beta regulatory T cell	200	109	CMPK2	2.57
naive T cell	180	123	TCF7	1.73
effector memory CD4-positive, alpha-beta T cell	412	286	ANXA1	1.88
effector memory CD8-positive, alpha-beta T cell	274	154	TRGC2	2.62
IgG plasma cell	154	153	IGLV6-57	9.79
conventional dendritic cell	1345	1345	CD1C	5.51
CD4-positive, CD25-positive, CCR4-positive, alpha-beta regulatory T cell	308	173	FOXP3	6.21
T follicular helper cell	369	164	CXCL13	4.61
exhausted T cell	448	323	ENSG00000237342	9.64
cycling T cell	1311	770	RRM2	6.20
cycling macrophage	2418	2289	CD1A	5.33

# Export top markers to CSV
top_markers_export <- all.markers %>%
  filter(p_val_adj < 0.05) %>%
  arrange(cluster, desc(avg_log2FC))

write.csv(top_markers_export, "cluster_markers.csv", row.names = FALSE)
cat("\nAll significant markers exported to 'cluster_markers.csv'\n")

## 
## All significant markers exported to 'cluster_markers.csv'

---

Clinical Metadata Analysis

Disease Subtype Comparisons

# Visualize cell distribution by disease subtype
if ("disease" %in% colnames(seurat_obj@meta.data)) {
  p1 <- DimPlot(seurat_obj, reduction = "umap", group.by = "disease", pt.size = 0.3) +
    ggtitle("UMAP by Disease Subtype") +
    theme(legend.position = "right", legend.text = element_text(size = 8))
  
  p2 <- DimPlot(seurat_obj, reduction = "umap", group.by = "cell_type", pt.size = 0.3) +
    ggtitle("UMAP by Cell Type") +
    theme(legend.position = "right", legend.text = element_text(size = 6))
  
  p1 | p2
}

# Cell type composition across disease subtypes
if ("disease" %in% colnames(seurat_obj@meta.data) && 
    "cell_type" %in% colnames(seurat_obj@meta.data)) {
  
  disease_celltype <- seurat_obj@meta.data %>%
    group_by(disease, cell_type) %>%
    summarise(count = n(), .groups = "drop") %>%
    group_by(disease) %>%
    mutate(percentage = count / sum(count) * 100)
  
  # Stacked bar chart
  ggplot(disease_celltype, aes(x = disease, y = percentage, fill = cell_type)) +
    geom_bar(stat = "identity") +
    coord_flip() +
    labs(title = "Cell Type Composition by Disease Subtype",
         x = "Disease", y = "Percentage", fill = "Cell Type") +
    theme_minimal() +
    theme(legend.text = element_text(size = 7))
  
  # Focus on immune cell differences
  # Calculate cell type proportions per disease
  celltype_props <- seurat_obj@meta.data %>%
    group_by(disease, cell_type) %>%
    summarise(count = n(), .groups = "drop") %>%
    group_by(disease) %>%
    mutate(total = sum(count),
           proportion = count / total) %>%
    dplyr::select(disease, cell_type, proportion) %>%
    pivot_wider(names_from = disease, values_from = proportion, values_fill = 0)
  
  celltype_props %>%
    head(20) %>%
    kable(digits = 3, caption = "Cell Type Proportions by Disease") %>%
    kable_styling(bootstrap_options = c("striped", "hover"), font_size = 10) %>%
    scroll_box(width = "100%")
}

Cell Type Proportions by Disease

cell_type	breast mucinous carcinoma	breast apocrine carcinoma	invasive tubular breast carcinoma || invasive lobular breast carcinoma	invasive ductal breast carcinoma	breast carcinoma	invasive lobular breast carcinoma	triple-negative breast carcinoma	metaplastic breast carcinoma	HER2 positive breast carcinoma	estrogen-receptor positive breast cancer	breast cancer
mast cell	0.06	0.000	0.141	0.019	0.000	0.036	0.000	0.000	0.000	0.024	0.014
macrophage	0.29	0.000	0.121	0.139	0.468	0.255	0.127	0.667	0.000	0.137	0.202
natural killer cell	0.10	0.154	0.128	0.082	0.137	0.145	0.038	0.000	0.176	0.097	0.035
myeloid dendritic cell	0.03	0.000	0.000	0.010	0.000	0.000	0.006	0.000	0.000	0.016	0.005
memory B cell	0.05	0.038	0.007	0.070	0.000	0.036	0.058	0.000	0.029	0.040	0.044
effector memory CD4-positive, alpha-beta T cell	0.04	0.423	0.148	0.089	0.048	0.018	0.101	0.000	0.059	0.056	0.065
effector memory CD8-positive, alpha-beta T cell	0.02	0.115	0.268	0.115	0.024	0.127	0.038	0.167	0.206	0.097	0.045
IgG plasma cell	0.06	0.000	0.007	0.063	0.000	0.091	0.081	0.000	0.118	0.073	0.156
conventional dendritic cell	0.20	0.077	0.060	0.073	0.073	0.127	0.017	0.000	0.059	0.016	0.027
T follicular helper cell	0.06	0.000	0.000	0.050	0.000	0.000	0.069	0.000	0.000	0.056	0.039
exhausted T cell	0.06	0.000	0.007	0.051	0.065	0.018	0.176	0.167	0.088	0.081	0.132
cycling T cell	0.01	0.000	0.000	0.021	0.040	0.000	0.029	0.000	0.000	0.032	0.037
cycling macrophage	0.02	0.000	0.007	0.016	0.024	0.000	0.012	0.000	0.000	0.008	0.021
naive T cell	0.00	0.115	0.040	0.061	0.000	0.055	0.046	0.000	0.206	0.121	0.028
CD4-positive, CD25-positive, CCR4-positive, alpha-beta regulatory T cell	0.00	0.077	0.040	0.088	0.048	0.018	0.092	0.000	0.029	0.081	0.088
CD8-positive, alpha-beta regulatory T cell	0.00	0.000	0.027	0.040	0.073	0.036	0.090	0.000	0.000	0.048	0.039
CD4-positive helper T cell	0.00	0.000	0.000	0.002	0.000	0.000	0.012	0.000	0.000	0.000	0.000
plasmacytoid dendritic cell	0.00	0.000	0.000	0.010	0.000	0.036	0.009	0.000	0.029	0.016	0.023

Triple-Negative vs ER-Positive Comparison

# Compare TNBC to ER-positive breast cancer
if ("disease" %in% colnames(seurat_obj@meta.data)) {
  
  # Subset to TNBC and ER+ samples
  tnbc_cells <- seurat_obj@meta.data$disease == "triple-negative breast carcinoma"
  erpos_cells <- seurat_obj@meta.data$disease == "estrogen-receptor positive breast cancer"
  
  if (sum(tnbc_cells) > 0 && sum(erpos_cells) > 0) {
    cat("Comparing TNBC (n=", sum(tnbc_cells), ") vs ER+ (n=", sum(erpos_cells), ")\n\n")
    
    # Create a binary classification
    seurat_obj$disease_comparison <- ifelse(
      tnbc_cells, "TNBC",
      ifelse(erpos_cells, "ER+", "Other")
    )
    
    # Visualize
    p <- DimPlot(seurat_obj, reduction = "umap", group.by = "disease_comparison", 
            cells.highlight = list(
              "TNBC" = colnames(seurat_obj)[tnbc_cells],
              "ER+" = colnames(seurat_obj)[erpos_cells]
            ),
            cols.highlight = c("red", "blue")) +
      ggtitle("TNBC vs ER-Positive Breast Cancer")
    print(p)
    
    # Cell type differences
    comparison_celltype <- seurat_obj@meta.data %>%
      filter(disease_comparison %in% c("TNBC", "ER+")) %>%
      group_by(disease_comparison, cell_type) %>%
      summarise(count = n(), .groups = "drop") %>%
      group_by(disease_comparison) %>%
      mutate(percentage = count / sum(count) * 100) %>%
      dplyr::select(disease_comparison, cell_type, percentage) %>%
      pivot_wider(names_from = disease_comparison, values_from = percentage, values_fill = 0) %>%
      mutate(difference = `TNBC` - `ER+`) %>%
      arrange(desc(abs(difference)))
    
    cat("\nCell type enrichment differences (TNBC vs ER+):\n")
    comparison_celltype %>%
      head(15) %>%
      kable(digits = 2, caption = "Top Cell Type Differences") %>%
      kable_styling(bootstrap_options = c("striped", "hover"))
    
    # Plot differences
    top_diff_celltypes <- comparison_celltype %>%
      head(10) %>%
      pivot_longer(cols = c("TNBC", "ER+"), names_to = "Disease", values_to = "Percentage")
    
    ggplot(top_diff_celltypes, aes(x = reorder(cell_type, difference), y = Percentage, fill = Disease)) +
      geom_bar(stat = "identity", position = "dodge") +
      coord_flip() +
      labs(title = "Top Cell Types with Differential Abundance (TNBC vs ER+)",
           x = "Cell Type", y = "Percentage") +
      theme_minimal()
  }
}

## Comparing TNBC (n= 346 ) vs ER+ (n= 124 )

## 
## Cell type enrichment differences (TNBC vs ER+):

Donor-Level Analysis

# Analyze cell composition per donor
if ("donor_id" %in% colnames(seurat_obj@meta.data)) {
  
  # Cells per donor
  donor_counts <- as.data.frame(table(seurat_obj$donor_id))
  colnames(donor_counts) <- c("Donor", "Cells")
  donor_counts <- donor_counts %>% arrange(desc(Cells))
  
  cat("Cells per donor (top 20):\n")
  donor_counts %>%
    head(20) %>%
    kable() %>%
    kable_styling(bootstrap_options = c("striped", "hover"))
  
  # Distribution
  ggplot(donor_counts, aes(x = Cells)) +
    geom_histogram(bins = 30, fill = "steelblue", color = "white") +
    labs(title = "Distribution of Cell Counts per Donor",
         x = "Number of Cells", y = "Number of Donors") +
    theme_minimal()
}

## Cells per donor (top 20):

Immune Checkpoint Gene Expression

# Key immune checkpoint genes
checkpoint_genes <- c("PDCD1", "CD274", "PDCD1LG2", "CTLA4", "LAG3", 
                      "HAVCR2", "TIGIT", "BTLA", "CD27", "ICOS")

# Check which are present
present_checkpoints <- checkpoint_genes[checkpoint_genes %in% rownames(seurat_obj)]

if (length(present_checkpoints) > 0) {
  cat("Found", length(present_checkpoints), "checkpoint genes:", 
      paste(present_checkpoints, collapse = ", "), "\n\n")
  
  # Feature plots
  if (length(present_checkpoints) <= 9) {
    FeaturePlot(seurat_obj, features = present_checkpoints, ncol = 3)
  } else {
    FeaturePlot(seurat_obj, features = head(present_checkpoints, 9), ncol = 3)
  }
  
  # Dot plot by cell type
  if ("cell_type" %in% colnames(seurat_obj@meta.data)) {
    DotPlot(seurat_obj, features = present_checkpoints, group.by = "cell_type") +
      coord_flip() +
      theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8)) +
      ggtitle("Immune Checkpoint Expression by Cell Type")
  }
  
  # Compare checkpoint expression across disease subtypes
  if ("disease" %in% colnames(seurat_obj@meta.data) && length(present_checkpoints) >= 1) {
    # Calculate average expression per disease
    checkpoint_expr <- FetchData(seurat_obj, vars = c(present_checkpoints, "disease"))
    
    checkpoint_avg <- checkpoint_expr %>%
      group_by(disease) %>%
      summarise(across(all_of(present_checkpoints), mean, na.rm = TRUE)) %>%
      pivot_longer(cols = -disease, names_to = "gene", values_to = "avg_expression")
    
    # Heatmap of checkpoint expression by disease
    checkpoint_matrix <- checkpoint_avg %>%
      pivot_wider(names_from = gene, values_from = avg_expression) %>%
      column_to_rownames("disease") %>%
      as.matrix()
    
    pheatmap(t(checkpoint_matrix),
             scale = "row",
             main = "Immune Checkpoint Expression by Disease Subtype",
             fontsize_row = 10,
             fontsize_col = 8)
  }
}

## Found 10 checkpoint genes: PDCD1, CD274, PDCD1LG2, CTLA4, LAG3, HAVCR2, TIGIT, BTLA, CD27, ICOS

Cytotoxic and Regulatory T Cell Analysis

# Focus on regulatory vs cytotoxic T cells
if ("cell_type" %in% colnames(seurat_obj@meta.data)) {
  
  # Identify T cell subsets
  tcell_types <- grep("T cell|regulatory", seurat_obj$cell_type, value = TRUE, ignore.case = TRUE)
  
  if (length(tcell_types) > 0) {
    cat("T cell types in dataset:\n")
    cat(paste(unique(tcell_types), collapse = "\n"), "\n\n")
    
    # Subset to T cells
    tcell_subset <- subset(seurat_obj, subset = cell_type %in% tcell_types)
    
    if (ncol(tcell_subset) > 0) {
      # UMAP of just T cells
      DimPlot(tcell_subset, reduction = "umap", group.by = "cell_type", label = TRUE) +
        ggtitle("T Cell Subsets") +
        theme(legend.position = "right")
      
      # Key T cell markers
      tcell_markers <- c("CD3D", "CD4", "CD8A", "FOXP3", "IL2RA", "GZMB", "PRF1", "IFNG")
      present_tcell_markers <- tcell_markers[tcell_markers %in% rownames(seurat_obj)]
      
      if (length(present_tcell_markers) > 0) {
        cat("Found T cell markers:", paste(present_tcell_markers, collapse = ", "), "\n")
        
        # Dot plot
        DotPlot(tcell_subset, features = present_tcell_markers, group.by = "cell_type") +
          coord_flip() +
          theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
          ggtitle("T Cell Marker Expression")
      }
      
      # Compare T cell proportions across disease subtypes
      if ("disease" %in% colnames(tcell_subset@meta.data)) {
        tcell_disease <- tcell_subset@meta.data %>%
          group_by(disease, cell_type) %>%
          summarise(count = n(), .groups = "drop") %>%
          group_by(disease) %>%
          mutate(proportion = count / sum(count))
        
        ggplot(tcell_disease, aes(x = disease, y = proportion, fill = cell_type)) +
          geom_bar(stat = "identity") +
          coord_flip() +
          labs(title = "T Cell Subset Proportions by Disease",
               x = "Disease", y = "Proportion") +
          theme_minimal() +
          theme(legend.text = element_text(size = 8))
      }
    }
  }
}

## T cell types in dataset:
## cycling T cell
## exhausted T cell
## effector memory CD8-positive, alpha-beta T cell
## CD4-positive, CD25-positive, CCR4-positive, alpha-beta regulatory T cell
## mast cell
## effector memory CD4-positive, alpha-beta T cell
## naive T cell
## CD8-positive, alpha-beta regulatory T cell
## CD4-positive helper T cell 
## 
## Found T cell markers: CD3D, CD4, CD8A, FOXP3, IL2RA, GZMB, PRF1, IFNG

Macrophage Polarization Analysis

# Ensure 'cell_type' exists
if ("cell_type" %in% colnames(seurat_obj@meta.data)) {

  # Identify macrophage populations
  mac_cells <- grep("macrophage", seurat_obj$cell_type, value = TRUE, ignore.case = TRUE)

  if (length(mac_cells) > 0) {
    cat("Macrophage populations:", paste(unique(mac_cells), collapse = ", "), "\n\n")

    # Subset macrophages
    mac_subset <- subset(seurat_obj, subset = cell_type %in% mac_cells)

    if (ncol(mac_subset) > 0) {
      cat("Total macrophages:", ncol(mac_subset), "\n\n")

      # Define polarization markers
      m1_markers <- c("CD86", "IL1B", "TNF", "NOS2", "CXCL10", "IL6")
      m2_markers <- c("CD163", "CD206", "MRC1", "ARG1", "IL10", "TGFB1", "CCL18")

      # Filter markers present in data
      present_m1 <- m1_markers[m1_markers %in% rownames(mac_subset)]
      present_m2 <- m2_markers[m2_markers %in% rownames(mac_subset)]

      cat("M1-like markers found:", paste(present_m1, collapse = ", "), "\n")
      cat("M2-like markers found:", paste(present_m2, collapse = ", "), "\n\n")

      all_mac_markers <- c(present_m1, present_m2)

      # Feature plots (paginated if >9 markers)
      if (length(all_mac_markers) > 0) {
        for (i in seq(1, length(all_mac_markers), by = 9)) {
          FeaturePlot(mac_subset, features = all_mac_markers[i:min(i+8, length(all_mac_markers))], ncol = 3)
        }
      }

      # DotPlot by disease
      if ("disease" %in% colnames(mac_subset@meta.data)) {
        print(
          DotPlot(mac_subset, features = all_mac_markers, group.by = "disease") +
            coord_flip() +
            theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8)) +
            ggtitle("Macrophage Polarization Markers by Disease")
        )
      }

      # Add M1/M2 module scores
      if (length(present_m1) > 0) {
        mac_subset <- AddModuleScore(mac_subset, features = list(present_m1), name = "M1_score")
      }
      if (length(present_m2) > 0) {
        mac_subset <- AddModuleScore(mac_subset, features = list(present_m2), name = "M2_score")
      }

      # FeaturePlot for M1/M2 scores
      if (all(c("M1_score1", "M2_score1") %in% colnames(mac_subset@meta.data))) {
        print(FeaturePlot(mac_subset, features = c("M1_score1", "M2_score1"), ncol = 2))
        mac_subset$M1_M2_ratio <- mac_subset$M1_score1 - mac_subset$M2_score1
      }

      # Extract macrophage scores
      mac_scores <- FetchData(mac_subset, vars = c("M1_score1", "M2_score1", "M1_M2_ratio", "disease", "cell_type"))

      # Scatter plot: M1 vs M2 by disease
      if ("disease" %in% colnames(mac_scores)) {
        p_scatter <- ggplot(mac_scores, aes(x = M1_score1, y = M2_score1, color = disease)) +
          geom_point(alpha = 0.5) +
          geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "black") +
          annotate("text", x = Inf, y = -Inf, label = "M2-skewed", hjust = 1.1, vjust = -0.5, size = 4, color = "blue") +
          annotate("text", x = -Inf, y = Inf, label = "M1-skewed", hjust = -0.1, vjust = 1.5, size = 4, color = "red") +
          labs(title = "M1 vs M2 Polarization by Disease", x = "M1 Score", y = "M2 Score") +
          theme_minimal() +
          theme(legend.position = "right")
        print(p_scatter)

        # Violin plot of M1/M2 ratio by disease
        p_violin <- ggplot(mac_scores, aes(x = disease, y = M1_M2_ratio, fill = disease)) +
          geom_violin(alpha = 0.7) +
          geom_boxplot(width = 0.1, outlier.shape = NA, fill = "white") +
          geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
          coord_flip() +
          labs(title = "M1/M2 Balance Across Disease Subtypes",
               subtitle = "Positive = M1-skewed, Negative = M2-skewed", x = "Disease", y = "M1-M2 Score Difference") +
          theme_minimal() +
          theme(legend.position = "none")
        print(p_violin)

        # Density plot
        p_density <- ggplot(mac_scores, aes(x = M1_M2_ratio, fill = disease)) +
          geom_density(alpha = 0.4) +
          geom_vline(xintercept = 0, linetype = "dashed", color = "black", size = 1) +
          labs(title = "Distribution of M1/M2 Balance Across Disease Subtypes",
               x = "M1-M2 Score Difference", y = "Density") +
          theme_minimal() +
          theme(legend.position = "right")
        print(p_density)
      }

      # Summary table
      mac_summary <- mac_scores %>%
        group_by(disease) %>%
        summarise(
          n_macrophages = n(),
          mean_M1_score = mean(M1_score1, na.rm = TRUE),
          mean_M2_score = mean(M2_score1, na.rm = TRUE),
          mean_ratio = mean(M1_M2_ratio, na.rm = TRUE),
          .groups = "drop"
        ) %>%
        mutate(
          polarization = case_when(
            mean_ratio > 0.1 ~ "M1-skewed",
            mean_ratio < -0.1 ~ "M2-skewed",
            TRUE ~ "Mixed"
          )
        ) %>%
        arrange(desc(mean_ratio))

      print(
        mac_summary %>%
          kable(digits = 3, caption = "Macrophage Polarization by Disease Subtype") %>%
          kable_styling(bootstrap_options = c("striped", "hover"))
      )
    }
  }
}

## Macrophage populations: macrophage, cycling macrophage 
## 
## Total macrophages: 1019 
## 
## M1-like markers found: CD86, IL1B, TNF, NOS2, CXCL10, IL6 
## M2-like markers found: CD163, MRC1, ARG1, IL10, TGFB1, CCL18

## <table class="table table-striped table-hover" style="margin-left: auto; margin-right: auto;">
## <caption>Macrophage Polarization by Disease Subtype</caption>
##  <thead>
##   <tr>
##    <th style="text-align:left;"> disease </th>
##    <th style="text-align:right;"> n_macrophages </th>
##    <th style="text-align:right;"> mean_M1_score </th>
##    <th style="text-align:right;"> mean_M2_score </th>
##    <th style="text-align:right;"> mean_ratio </th>
##    <th style="text-align:left;"> polarization </th>
##   </tr>
##  </thead>
## <tbody>
##   <tr>
##    <td style="text-align:left;"> invasive lobular breast carcinoma </td>
##    <td style="text-align:right;"> 14 </td>
##    <td style="text-align:right;"> 0.276 </td>
##    <td style="text-align:right;"> -0.187 </td>
##    <td style="text-align:right;"> 0.462 </td>
##    <td style="text-align:left;"> M1-skewed </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> invasive ductal breast carcinoma </td>
##    <td style="text-align:right;"> 502 </td>
##    <td style="text-align:right;"> 0.082 </td>
##    <td style="text-align:right;"> -0.055 </td>
##    <td style="text-align:right;"> 0.137 </td>
##    <td style="text-align:left;"> M1-skewed </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> triple-negative breast carcinoma </td>
##    <td style="text-align:right;"> 48 </td>
##    <td style="text-align:right;"> -0.039 </td>
##    <td style="text-align:right;"> -0.068 </td>
##    <td style="text-align:right;"> 0.030 </td>
##    <td style="text-align:left;"> Mixed </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> metaplastic breast carcinoma </td>
##    <td style="text-align:right;"> 4 </td>
##    <td style="text-align:right;"> -0.126 </td>
##    <td style="text-align:right;"> -0.108 </td>
##    <td style="text-align:right;"> -0.018 </td>
##    <td style="text-align:left;"> Mixed </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> breast cancer </td>
##    <td style="text-align:right;"> 322 </td>
##    <td style="text-align:right;"> -0.158 </td>
##    <td style="text-align:right;"> -0.105 </td>
##    <td style="text-align:right;"> -0.053 </td>
##    <td style="text-align:left;"> Mixed </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> estrogen-receptor positive breast cancer </td>
##    <td style="text-align:right;"> 18 </td>
##    <td style="text-align:right;"> -0.303 </td>
##    <td style="text-align:right;"> -0.168 </td>
##    <td style="text-align:right;"> -0.135 </td>
##    <td style="text-align:left;"> M2-skewed </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> invasive tubular breast carcinoma &amp;#124;&amp;#124; invasive lobular breast carcinoma </td>
##    <td style="text-align:right;"> 19 </td>
##    <td style="text-align:right;"> -0.239 </td>
##    <td style="text-align:right;"> -0.037 </td>
##    <td style="text-align:right;"> -0.202 </td>
##    <td style="text-align:left;"> M2-skewed </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> breast mucinous carcinoma </td>
##    <td style="text-align:right;"> 31 </td>
##    <td style="text-align:right;"> 0.000 </td>
##    <td style="text-align:right;"> 0.215 </td>
##    <td style="text-align:right;"> -0.215 </td>
##    <td style="text-align:left;"> M2-skewed </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> breast carcinoma </td>
##    <td style="text-align:right;"> 61 </td>
##    <td style="text-align:right;"> -0.107 </td>
##    <td style="text-align:right;"> 0.285 </td>
##    <td style="text-align:right;"> -0.392 </td>
##    <td style="text-align:left;"> M2-skewed </td>
##   </tr>
## </tbody>
## </table>

Tumor Microenvironment Signatures

# Define gene signatures for TME analysis
tme_signatures <- list(
  Cytotoxicity = c("GZMA", "GZMB", "PRF1", "GNLY", "NKG7"),
  Exhaustion = c("PDCD1", "HAVCR2", "LAG3", "TIGIT", "TOX"),
  Activation = c("CD69", "CD25", "IL2RA", "ICOS", "CD40LG"),
  Proliferation = c("MKI67", "TOP2A", "PCNA", "CDK1"),
  Immunosuppression = c("TGFB1", "IL10", "ARG1", "IDO1", "CD274")
)

# Calculate module scores for each signature
for (sig_name in names(tme_signatures)) {
  genes <- tme_signatures[[sig_name]]
  present_genes <- genes[genes %in% rownames(seurat_obj)]
  
  if (length(present_genes) > 0) {
    seurat_obj <- AddModuleScore(
      seurat_obj,
      features = list(present_genes),
      name = paste0(sig_name, "_score")
    )
  }
}

# Visualize signatures
score_names <- paste0(names(tme_signatures), "_score1")
present_scores <- score_names[score_names %in% colnames(seurat_obj@meta.data)]

if (length(present_scores) > 0) {
  FeaturePlot(seurat_obj, features = present_scores, ncol = 3)
  
  # Compare across disease subtypes
  if ("disease" %in% colnames(seurat_obj@meta.data)) {
    score_data <- FetchData(seurat_obj, vars = c(present_scores, "disease"))
    
    score_avg <- score_data %>%
      group_by(disease) %>%
      summarise(across(all_of(present_scores), mean, na.rm = TRUE)) %>%
      pivot_longer(cols = -disease, names_to = "signature", values_to = "score")
    
    ggplot(score_avg, aes(x = disease, y = score, fill = signature)) +
      geom_bar(stat = "identity", position = "dodge") +
      coord_flip() +
      labs(title = "TME Signatures by Disease Subtype",
           x = "Disease", y = "Average Score") +
      theme_minimal() +
      theme(legend.position = "right")
  }
}

Statistical Testing: Disease-Associated Markers

# Differential expression between disease subtypes
if ("disease" %in% colnames(seurat_obj@meta.data)) {
  
  # Set disease as identity
  Idents(seurat_obj) <- "disease"
  
  # Get disease categories with sufficient cells
  disease_counts <- table(seurat_obj$disease)
  diseases_to_compare <- names(disease_counts[disease_counts >= 50])
  
  if (length(diseases_to_compare) >= 2) {
    cat("Comparing diseases with >=50 cells:\n")
    cat(paste(diseases_to_compare, collapse = "\n"), "\n\n")
    
    # Example: TNBC vs invasive ductal
    if ("triple-negative breast carcinoma" %in% diseases_to_compare &&
        "invasive ductal breast carcinoma" %in% diseases_to_compare) {
      
      cat("Finding markers: TNBC vs Invasive Ductal Carcinoma\n")
      
      tnbc_markers <- FindMarkers(
        seurat_obj,
        ident.1 = "triple-negative breast carcinoma",
        ident.2 = "invasive ductal breast carcinoma",
        min.pct = 0.1,
        logfc.threshold = 0.25
      )
      
      tnbc_markers$gene <- rownames(tnbc_markers)
      
      cat("\nTop genes upregulated in TNBC:\n")
      tnbc_markers %>%
        filter(avg_log2FC > 0, p_val_adj < 0.05) %>%
        arrange(desc(avg_log2FC)) %>%
        head(15) %>%
        dplyr::select(gene, avg_log2FC, pct.1, pct.2, p_val_adj) %>%
        kable(digits = 3) %>%
        kable_styling(bootstrap_options = c("striped", "hover"))
      
      cat("\nTop genes downregulated in TNBC:\n")
      tnbc_markers %>%
        filter(avg_log2FC < 0, p_val_adj < 0.05) %>%
        arrange(avg_log2FC) %>%
        head(15) %>%
        dplyr::select(gene, avg_log2FC, pct.1, pct.2, p_val_adj) %>%
        kable(digits = 3) %>%
        kable_styling(bootstrap_options = c("striped", "hover"))
      
      # Volcano plot
      tnbc_markers <- tnbc_markers %>%
        mutate(
          significance = case_when(
            p_val_adj < 0.05 & avg_log2FC > 0.5 ~ "Up in TNBC",
            p_val_adj < 0.05 & avg_log2FC < -0.5 ~ "Down in TNBC",
            TRUE ~ "Not significant"
          )
        )
      
      ggplot(tnbc_markers, aes(x = avg_log2FC, y = -log10(p_val_adj), color = significance)) +
        geom_point(alpha = 0.6) +
        scale_color_manual(values = c("Up in TNBC" = "red", 
                                      "Down in TNBC" = "blue", 
                                      "Not significant" = "grey")) +
        geom_vline(xintercept = c(-0.5, 0.5), linetype = "dashed") +
        geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
        labs(title = "Differential Expression: TNBC vs Invasive Ductal Carcinoma",
             x = "Log2 Fold Change", y = "-Log10(Adjusted P-value)") +
        theme_minimal()
    }
  }
}

## Comparing diseases with >=50 cells:
## breast mucinous carcinoma
## invasive tubular breast carcinoma || invasive lobular breast carcinoma
## invasive ductal breast carcinoma
## breast carcinoma
## invasive lobular breast carcinoma
## triple-negative breast carcinoma
## estrogen-receptor positive breast cancer
## breast cancer 
## 
## Finding markers: TNBC vs Invasive Ductal Carcinoma
## 
## Top genes upregulated in TNBC:
## 
## Top genes downregulated in TNBC:

---

Save Results

saveRDS(seurat_obj, file = "immune_compartment_subset_seurat_final.rds")
cat("Saved Seurat object to: immune_compartment_subset_seurat_final.rds\n")

## Saved Seurat object to: immune_compartment_subset_seurat_final.rds

---

Session Information

sessionInfo()

## R version 4.4.1 (2024-06-14 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 22631)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] ggrepel_0.9.6      biomaRt_2.58.2     sctransform_0.4.1  Matrix_1.7-0      
##  [5] reticulate_1.40.0  viridis_0.6.5      viridisLite_0.4.2  pheatmap_1.0.12   
##  [9] scales_1.4.0       RColorBrewer_1.1-3 DT_0.33            kableExtra_1.4.0  
## [13] patchwork_1.3.1    lubridate_1.9.3    forcats_1.0.0      stringr_1.5.1     
## [17] dplyr_1.1.4        purrr_1.0.2        readr_2.1.5        tidyr_1.3.1       
## [21] tibble_3.2.1       ggplot2_4.0.0      tidyverse_2.0.0    Seurat_5.1.0      
## [25] SeuratObject_5.0.2 sp_2.1-4          
## 
## loaded via a namespace (and not attached):
##   [1] RcppAnnoy_0.0.22        splines_4.4.1           later_1.3.2            
##   [4] filelock_1.0.3          polyclip_1.10-7         XML_3.99-0.17          
##   [7] fastDummies_1.7.4       lifecycle_1.0.4         globals_0.16.3         
##  [10] lattice_0.22-6          MASS_7.3-61             magrittr_2.0.3         
##  [13] limma_3.60.6            plotly_4.10.4           sass_0.4.9             
##  [16] rmarkdown_2.28          jquerylib_0.1.4         yaml_2.3.10            
##  [19] httpuv_1.6.15           spam_2.11-0             spatstat.sparse_3.1-0  
##  [22] cowplot_1.1.3           pbapply_1.7-2           DBI_1.2.3              
##  [25] abind_1.4-8             zlibbioc_1.50.0         Rtsne_0.17             
##  [28] presto_1.0.0            BiocGenerics_0.50.0     rappdirs_0.3.3         
##  [31] GenomeInfoDbData_1.2.12 IRanges_2.38.1          S4Vectors_0.42.1       
##  [34] irlba_2.3.5.1           listenv_0.9.1           spatstat.utils_3.1-2   
##  [37] goftest_1.2-3           RSpectra_0.16-2         spatstat.random_3.3-2  
##  [40] fitdistrplus_1.2-2      parallelly_1.41.0       svglite_2.1.3          
##  [43] leiden_0.4.3.1          codetools_0.2-20        xml2_1.3.6             
##  [46] tidyselect_1.2.1        UCSC.utils_1.0.0        farver_2.1.2           
##  [49] matrixStats_1.4.1       stats4_4.4.1            BiocFileCache_2.12.0   
##  [52] spatstat.explore_3.3-4  jsonlite_1.8.9          progressr_0.15.1       
##  [55] ggridges_0.5.6          survival_3.7-0          systemfonts_1.3.1      
##  [58] tools_4.4.1             progress_1.2.3          ica_1.0-3              
##  [61] Rcpp_1.1.0              glue_1.7.0              gridExtra_2.3          
##  [64] xfun_0.46               GenomeInfoDb_1.40.1     withr_3.0.2            
##  [67] fastmap_1.2.0           digest_0.6.36           timechange_0.3.0       
##  [70] R6_2.5.1                mime_0.12               colorspace_2.1-1       
##  [73] scattermore_1.2         tensor_1.5              dichromat_2.0-0.1      
##  [76] spatstat.data_3.1-4     RSQLite_2.3.7           generics_0.1.3         
##  [79] data.table_1.16.0       prettyunits_1.2.0       httr_1.4.7             
##  [82] htmlwidgets_1.6.4       uwot_0.2.2              pkgconfig_2.0.3        
##  [85] gtable_0.3.6            blob_1.2.4              lmtest_0.9-40          
##  [88] S7_0.2.0                XVector_0.44.0          htmltools_0.5.8.1      
##  [91] dotCall64_1.2           Biobase_2.64.0          png_0.1-8              
##  [94] spatstat.univar_3.1-1   knitr_1.48              rstudioapi_0.16.0      
##  [97] tzdb_0.4.0              reshape2_1.4.4          curl_5.2.3             
## [100] nlme_3.1-166            cachem_1.1.0            zoo_1.8-12             
## [103] KernSmooth_2.23-24      parallel_4.4.1          miniUI_0.1.1.1         
## [106] AnnotationDbi_1.66.0    pillar_1.10.1           grid_4.4.1             
## [109] vctrs_0.6.5             RANN_2.6.2              promises_1.3.0         
## [112] dbplyr_2.5.0            xtable_1.8-4            cluster_2.1.6          
## [115] evaluate_1.0.0          cli_3.6.3               compiler_4.4.1         
## [118] rlang_1.1.4             crayon_1.5.3            future.apply_1.11.3    
## [121] labeling_0.4.3          plyr_1.8.9              stringi_1.8.4          
## [124] deldir_2.0-4            Biostrings_2.72.1       lazyeval_0.2.2         
## [127] spatstat.geom_3.3-4     RcppHNSW_0.6.0          hms_1.1.3              
## [130] bit64_4.5.2             future_1.34.0           statmod_1.5.0          
## [133] KEGGREST_1.44.1         shiny_1.9.1             highr_0.11             
## [136] ROCR_1.0-11             igraph_2.0.3            memoise_2.0.1          
## [139] bslib_0.8.0             bit_4.5.0