How to plot UMAP for control vs treatment scRNA seq samples which will depict clearly control and treatment clusters
2
0
Entering edit mode
2.0 years ago

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

enter image description here

scRNAseq Seurat • 1.8k views
ADD COMMENT
3
Entering edit mode
2.0 years ago
EagleEye 7.5k

Use split.by from DimPlot

ADD COMMENT
0
Entering edit mode

Thank you EagleEye for your reply. I have tried with your suggestion but did not work out.

ADD REPLY
1
Entering edit mode
2.0 years ago
Soheil ▴ 110

split.by should work to separate your umap based on conditions, but if you want to separate/color each sample you can create a new column in your object using the cell names. Then you can use the split.by/group.by (or both) to split/color your umap using that column. sth like this:

disease.combined$sample <- sapply(strsplit(colnames(disease.combined), split = "_"), "[", 1)
ADD COMMENT

Login before adding your answer.

Traffic: 1948 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6