Hello,

I need help regarding making UMAP showing clusters which will be labeled by groups control(n=3) and treatment(n=3). I have attached UMAP image 1 control vs 1 Treatment for your reference

I am working with scRNAseq data from 3 treatments and 3 controls. I am able to integrate them, create a Seurat file and obtain nice clusters. However I am struggling to group the 3 controls and the 3 treatments together to compare gene expression between patients and controls. I am unable to find out how to show cell clusters by treatments vs control. Please review my code-

 # For 3 controls vs 3 Treatments- New code analysis
library(Seurat)
library(SeuratObject)
library("Seurat")
library("sctransform")
library("dplyr")
library("RColorBrewer")
library("ggthemes")
library("ggplot2")
library("cowplot")
library("data.table")



ctrl1st2.data <-Read10X(data.dir = "D:/30-649959914/01_analysis/cellranger_count/ST2-CTRL/filtered_feature_bc_matrix")

ctrl2st2.data <-Read10X(data.dir = "D:/30-649959914/01_analysis/cellranger_count/ST-2-CTRL-1/filtered_feature_bc_matrix")


ctrl3st2.data <-Read10X(data.dir = "D:/30-649959914/01_analysis/cellranger_count/ST-2-CTRL-2/filtered_feature_bc_matrix")

control.list = list() # First create an empty list to hold the Seurat objects
control.list[[1]] = CreateSeuratObject(counts = ctrl1st2.data, 
                                       project = "ctrl1st2"
)

control.list[[2]] = CreateSeuratObject(counts = ctrl2st2.data, 
                                       project = "ctrl2st2")

control.list[[3]] = CreateSeuratObject(counts = ctrl3st2.data, 
                                       project = "ctrl3st2"
)


controldisease <- merge(x=control.list[[1]], y=c(control.list[[2]],control.list[[3]]), add.cell.ids = c("ctrl1st2","ctrl2st2","ctrl3st2"), project="CSHL")

controldisease[["percent.mt"]] <- PercentageFeatureSet(object = controldisease, pattern = "^MT-")

 p1 <- VlnPlot(object = controldisease, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

 plot(p1)

 p2 <- FeatureScatter(object = controldisease, feature1 = "nCount_RNA", feature2 = "percent.mt")

plot(p2)
controldisease <- subset(x =  controldisease, subset = nFeature_RNA > 200 & nFeature_RNA < 8000 & percent.mt < 20)
 p3 <- VlnPlot(object = controldisease, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

plot(p3)

p4 <- FeatureScatter(object = controldisease, feature1 = "nCount_RNA", feature2 = "percent.mt")

plot(p4)                               


treatment1st2.data <-Read10X(data.dir = "D:/30-649959914/01_analysis/cellranger_count/ST2-treatment-7/filtered_feature_bc_matrix")
treatment2st2.data <-Read10X(data.dir = "D:/30-649959914/01_analysis/cellranger_count/ST-2-treatment-7/filtered_feature_bc_matrix")
treatment3st2.data <-Read10X(data.dir = "D:/30-649959914/01_analysis/cellranger_count/ST--2-treatment-7/filtered_feature_bc_matrix")

treatment.list = list() # First create an empty list to hold the Seurat objects
treatment.list[[1]] = CreateSeuratObject(counts = treatment1st2.data, 
                                   project = "treatment1st2"
)

treatment.list[[2]] = CreateSeuratObject(counts = treatment2st2.data, 
                                   project = "treatment2st2")

treatment.list[[3]] = CreateSeuratObject(counts = treatment3st2.data, 
                                   project = "treatment3st2"
)


treatmentdisease <- merge(x=treatment.list[[1]], y=c(treatment.list[[2]],treatment.list[[3]]), add.cell.ids = c("treatment1st2","treatment2st2","treatment3st2"), project="CSHL")

treatmentdisease[["percent.mt"]] <- PercentageFeatureSet(object = treatmentdisease, pattern = "^MT-")

p5 <- VlnPlot(object = treatmentdisease, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

plot(p5)

p6 <- FeatureScatter(object = treatmentdisease, feature1 = "nCount_RNA", feature2 = "percent.mt")

plot(p6)
treatmentdisease <- subset(x =  treatmentdisease, subset = nFeature_RNA > 200 & nFeature_RNA < 8000 & percent.mt < 20)
p7 <- VlnPlot(object = treatmentdisease, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

plot(p7)

p8 <- FeatureScatter(object = treatmentdisease, feature1 = "nCount_RNA", feature2 = "percent.mt")

plot(p8)                               

diseasest2comp <- c(controldisease,treatmentdisease)
diseasest2comp@meta.data$condition <- sapply(diseasest2comp@meta.data$orig.ident, function(ita) ifelse(ita %in% control_samples,"Control","treatment7")
controldisease$treatmentdisease <- "CONTROL"
)treatmentdisease$treatmentdisease     <- "treatment7"

diseasest2comp <- lapply(X = diseasest2comp, FUN = function(x) {
  x <- NormalizeData(x)
  x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 3000)
})
features <- SelectIntegrationFeatures(object.list = diseasest2comp)
diseasest2comp <- lapply(X = diseasest2comp, FUN = function(x) {
  x <- ScaleData(x, features = features, verbose = FALSE)
  x <- RunPCA(x, features = features, verbose = FALSE)
})

disease.anchors <- FindIntegrationAnchors(object.list = diseasest2comp, anchor.features = features, reduction = "rpca")
disease.combined <- IntegrateData(anchorset = disease.anchors)
DefaultAssay(disease.combined) <- "integrated"
disease.combined <- ScaleData(disease.combined, verbose = FALSE)
disease.combined <- RunPCA(disease.combined, npcs = 30, verbose = FALSE)
disease.combined <- RunUMAP(disease.combined, reduction = "pca", dims = 1:30)
disease.combined <- FindNeighbors(disease.combined, reduction = "pca", dims = 1:30)
disease.combined <- FindClusters(disease.combined, resolution = 0.5)

p9 <- DimPlot(disease.combined, reduction = "umap", raster=FALSE)
p10 <- DimPlot(disease.combined, reduction = "umap",  label = TRUE,
              repel = TRUE, raster=FALSE)
p9 + p10