3_Human_data_analysis

# --------------------------
# Load required libraries
# --------------------------
library(Seurat)
library(SeuratObject)
library(SeuratWrappers)
library(future)
library(future.apply)
library(tidyverse)
library(RCurl)
library(dplyr)
library(paletteer)
library(patchwork)
library(readxl)
library(ggpubr)
library(ggrepel)

# --------------------------
# Load Seurat object
# --------------------------
# The RData file should contain "human_obj" (Seurat object with metadata + RNA assay)
load('human_obj.RData')

head(human_obj@meta.data, 3)   # Quick check

# --------------------------
# Extract CD276 expression by sample
# --------------------------
cd276_expr <- FetchData(human_obj, vars = c("CD276", "sample_ID"))

# Compute per-sample average CD276 expression
cd276_by_sample <- cd276_expr %>%
  group_by(sample_ID) %>%
  summarise(mean_CD276 = mean(CD276, na.rm = TRUE)) %>%
  arrange(desc(mean_CD276))

summary(cd276_by_sample$mean_CD276)   # Get min/median/max stats

# Define high vs low groups using median as cutoff
cd276_median <- median(cd276_by_sample$mean_CD276, na.rm = TRUE)
CD276_High_sample_ID <- as.character(cd276_by_sample[cd276_by_sample$mean_CD276 > cd276_median, ]$sample_ID)
CD276_Low_sample_ID  <- as.character(cd276_by_sample[cd276_by_sample$mean_CD276 <= cd276_median, ]$sample_ID)

# --------------------------
# Define signature genes
# --------------------------
# Derived from mouse markers → filtered to keep only consistent human trends
Gr_MDSC_signatures <- c("CXCR2", "ACTB", "TMSB4X", "S100A9", "IFITM2", "CD52", "CSF3R",
                        "S100A8", "TYROBP", "IL1B", "GMFG", "LCP1", "CORO1A", "JAML",
                        "SNRPG", "PFN1", "COX6A1", "RGS2", "RBX1", "COX7A2", "NEDD8",
                        "CST3", "RNF213")

M_MDSC_signature <- c("CCL3", "CXCL2", "CSTB", "CCL4", "IL1RN", "ACOD1", "SQSTM1", "F10",
                      "ZEB2", "CTSB", "RGCC", "CLEC4E", "PLEK", "IL1A", "TNF", "BRI3",
                      "SOD2", "NLRP3", "ATP6V0C", "SLC31A2", "BASP1", "EIF5", "IFRD1",
                      "FNIP2", "TNFRSF1B", "EEA1", "EHD1", "RALGDS", "SEMA4D", "CD63",
                      "MAPKAPK2", "EIF4A1", "CCRL2", "HMOX1", "TXNRD1", "GADD45B",
                      "ATP6V1C1", "NPC1", "ABCC1", "PTS", "MAFF", "RAB20", "LGALS3", "SPP1")

CAF_signatures <- c("COL4A1", "COL4A2", "COL18A1", "ACTA2", "TAGLN", "MYH11", "MYL9", 
                    "POSTN", "MMP11", "SPARC", "BGN", "NME2", "INHBA", "FN1", "IGFBP7", "CTSZ")

# --------------------------
# Subset and analyze MDSCs
# --------------------------
MDSC_obj <- subset(human_obj, subset = cellstate %in% c(
  "Polymorphonuclear myeloid-derived suppressor cells (PMN-MDSC)",
  "Mononuclear myeloid-derived suppressor cells (M-MDSC)"
))

# Add simplified labels
MDSC_obj$sample_ID <- MDSC_obj$`GSM ID`
MDSC_obj$cellstate_short <- as.character(MDSC_obj$cellstate)
MDSC_obj@meta.data[MDSC_obj$cellstate_short == "Polymorphonuclear myeloid-derived suppressor cells (PMN-MDSC)", ]$cellstate_short <- "Gr-MDSC"
MDSC_obj@meta.data[MDSC_obj$cellstate_short == "Mononuclear myeloid-derived suppressor cells (M-MDSC)", ]$cellstate_short <- "M-MDSC"

# Standard Seurat preprocessing
DefaultAssay(MDSC_obj) <- "RNA"
MDSC_obj <- NormalizeData(MDSC_obj) %>% ScaleData() %>% FindVariableFeatures(nfeatures = 3000)
MDSC_obj <- RunPCA(MDSC_obj) %>% FindNeighbors(dims = 1:15) %>% RunUMAP(dims = 1:15) %>% FindClusters(resolution = 0.1)

# Visualize UMAP
DimPlot(MDSC_obj, group.by = "cellstate_short")

# Add module scores for Gr-MDSC and M-MDSC signatures
MDSC_obj <- AddModuleScore(MDSC_obj,
                           features = list(Gr_MDSC_signatures, M_MDSC_signature),
                           name = c("Gr_feature", "M_feature"),
                           ctrl = 5)

colnames(MDSC_obj@meta.data)[42:43] <- c("Gr-MDSC signature", "M-MDSC signature")

# Violin plots of signature scores split by CD276 group
pdf('MDSC_obj_signature.pdf', height = 4, width = 5)
VlnPlot2(MDSC_obj,
         features = c("Gr-MDSC signature", "M-MDSC signature"),
         group.by = "CD276_Tumor_group",
         stat.method = "wilcox.test",
         cols = c("#A666E1", "#EEAAEE"),
         pt = FALSE) +
  labs(y = "Signature module score") +
  theme(axis.title.y = element_text(size = 12))
dev.off()

# --------------------------
# Subset and analyze CAFs
# --------------------------
CAF_obj <- subset(human_obj, subset = cellstate %in% c("Fibroblast"))
DefaultAssay(CAF_obj) <- "RNA"
CAF_obj <- NormalizeData(CAF_obj) %>% ScaleData() %>% FindVariableFeatures(nfeatures = 3000)
CAF_obj <- RunPCA(CAF_obj) %>% FindNeighbors(dims = 1:20) %>% RunUMAP(dims = 1:20) %>% FindClusters(resolution = 0.1)

# CAF signature scoring
CAF_obj <- AddModuleScore(CAF_obj,
                          features = list(CAF_signatures),
                          name = "CAF_feature",
                          ctrl = 5)

colnames(CAF_obj@meta.data)[45] <- "CAF_signature"

pdf('CAF_obj_signature.pdf', height = 4, width = 2.5)
VlnPlot2(CAF_obj,
         features = "CAF_signature",
         group.by = "CD276_Tumor_group",
         stat.method = "wilcox.test",
         cols = c("#A666E1", "#EEAAEE"),
         pt = FALSE) +
  labs(y = "Signature module score") +
  theme(axis.title.y = element_text(size = 12))
dev.off()


# Differential expression in CAFs between CD276 high vs low tumors
Idents(CAF_obj) <- CAF_obj$CD276_Tumor_group
DEG_all_fibroblast <- FindMarkers(CAF_obj,
                                  ident.1 = "CD276 high",
                                  ident.2 = "CD276 low",
                                  logfc.threshold = 0.1,
                                  test.use = "MAST")

# Volcano plot for CAF DEGs, highlighting ECM remodeling genes
Plot_genes <- c("COL4A1", "COL4A2", "COL18A1", "ACTA2", "TAGLN", "MYH11", "MYL9", "MMP11", "BGN", "IGFBP7")

pdf('EnhancedVolcano_CAF.pdf', height = 8, width = 6)
EnhancedVolcano(DEG_all_fibroblast,
                lab = rownames(DEG_all_fibroblast),
                x = 'avg_log2FC',
                y = 'p_val_adj',
                title = NULL,
                pCutoff = 0.05,
                FCcutoff = 0.1,
                selectLab = Plot_genes,
                labSize = 5,
                max.overlaps = Inf,
                xlim = c(-5, 5),
                pointSize = 3.5,
                colAlpha = 1,
                drawConnectors = TRUE,
                boxedLabels = TRUE,
                col = c('grey30', 'grey30', 'grey30', 'red3')) +
  theme_classic()
dev.off()

# --------------------------
# Average expression analysis
# --------------------------
genes <- c("CD276", Gr_MDSC_signatures, M_MDSC_signature, CAF_signatures)
genes_present <- genes[genes %in% rownames(human_obj)]   # keep only available genes

expr_df <- FetchData(human_obj, vars = c(genes_present, "sample_ID"))

# Compute average per sample
avg_expr <- expr_df %>%
  group_by(sample_ID) %>%
  summarise(across(all_of(genes_present), mean, na.rm = TRUE))

# Add sample information
sample_infor <- distinct(human_obj@meta.data[, c("sample_ID", "Group")], sample_ID, .keep_all = TRUE)
avg_expr <- left_join(data.frame(avg_expr), sample_infor, by = "sample_ID")

# Assign tumor groups (high vs low CD276)
avg_expr$CD276_Tumor_group <- "CD276 high"
avg_expr[avg_expr$sample_ID %in% CD276_Low_sample_ID, ]$CD276_Tumor_group <- "CD276 low"
avg_expr$CD276_Tumor_group <- factor(avg_expr$CD276_Tumor_group, levels = c("CD276 high", "CD276 low"))

# Boxplot of CD276 expression
pdf('CD276_ggboxplot.pdf', height = 4, width = 3)
ggboxplot(avg_expr, x = "CD276_Tumor_group", y = "CD276",
          color = "CD276_Tumor_group",
          palette = c("#A666E1", "#EEAAEE"),
          add = "jitter") +
  stat_compare_means(label = "p.format") +
  labs(x = NULL, y = "Expression of CD276") +
  theme(legend.position = "none")
dev.off()

# --------------------------
# Correlation analysis: CAF vs Tumor CD276
# --------------------------
# Compute tumor-level mean CD276
DefaultAssay(human_obj) <- "RNA"
cd276_vals <- GetAssayData(human_obj, assay = "RNA", layer = "data")["CD276", ]
human_obj <- AddMetaData(human_obj, metadata = cd276_vals, col.name = "CD276_expr")

tumor_cd276_by_sample <- human_obj@meta.data %>%
  group_by(sample_ID) %>%
  summarise(mean_CD276 = mean(CD276_expr, na.rm = TRUE), .groups = "drop")

# Average CAF gene expression per sample
CAF_expr_df <- FetchData(CAF_obj, vars = c(CAF_signatures, "sample_ID"))
CAF_avg_expr <- CAF_expr_df %>% group_by(sample_ID) %>%
  summarise(across(all_of(CAF_signatures), mean, na.rm = TRUE))

# Merge CAF and tumor CD276 info
merged_celltype_expr_df <- left_join(tumor_cd276_by_sample, CAF_avg_expr, by = "sample_ID")
merged_celltype_expr_df <- left_join(merged_celltype_expr_df, sample_infor, by = "sample_ID")

merged_celltype_expr_df$CD276_Tumor_group <- "CD276 high"
merged_celltype_expr_df[merged_celltype_expr_df$sample_ID %in% CD276_Low_sample_ID, ]$CD276_Tumor_group <- "CD276 low"
merged_celltype_expr_df$CD276_Tumor_group <- factor(merged_celltype_expr_df$CD276_Tumor_group,
                                                    levels = c("CD276 high", "CD276 low"))

# Example: scatterplot correlation between CD276 and COL4A1 in CAFs
g <- "COL4A1"
pdf("ggscatter_COL4A1.pdf", height = 5, width = 5)
ggscatter(merged_celltype_expr_df,
          x = "mean_CD276",
          y = g,
          add = "reg.line",
          conf.int = TRUE,
          cor.coef = TRUE,
          cor.method = "spearman",
          xlab = "Expression of CD276",
          ylab = paste0("Expression of ", g),
          add.params = list(color = "blue", fill = "lightgray")) +
  theme_classic()
dev.off()