Entering edit mode
10 months ago
Researcher
▴
30
Hi, Thanks to the community for solving my previous queries regarding sc-RNA seq analysis. I have performed differential expression analysis for T cells, of two different conditions "car" and "acar". I will provide my R code. Kindly let me know, if i have done my workflow correctly.
pbmc.data <- Read10X(data.dir = "/home/dell/filtered_feature_bc_matrix1")
pbmccar <- CreateSeuratObject(counts = pbmc.data, project = "car", min.cells = 3, min.features=200)
pbmc.data <- Read10X(data.dir = "/home/dell/filtered_feature_bc_matrix")
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "acar", min.cells = 3, min.features=200)
pbmc.combined <- merge(pbmccar, y = pbmc, add.cell.ids = c("car", "acar"), project = "car-t")
head(colnames(pbmc.combined))
table(pbmc.combined$orig.ident)
pbmc.combined <- NormalizeData(pbmc.combined)
pbmc.combined <- FindVariableFeatures(pbmc.combined)
pbmc.combined <- ScaleData(pbmc.combined)
pbmc.combined <- RunPCA(pbmc.combined)
pbmc.combined <- FindNeighbors(pbmc.combined, dims = 1:30, reduction = "pca")
pbmc.combined <- FindClusters(pbmc.combined, resolution = 2, cluster.name = "uclusters")
pbmc.combined <- RunUMAP(pbmc.combined, dims = 1:30, reduction = "pca", reduction.name = "umap.u")
DimPlot(pbmc.combined, reduction = "umap.u", group.by = c("orig.ident", "seurat_clusters"))
pbmc.combined <- IntegrateLayers(object = pbmc.combined, method = CCAIntegration, orig.reduction = "pca", new.reduction = "integrated.cca",
verbose = FALSE)
# re-join layers after integration
pbmc.combined[["RNA"]] <- JoinLayers(pbmc.combined[["RNA"]])
pbmc.combined <- FindNeighbors(pbmc.combined, reduction = "integrated.cca", dims = 1:30)
pbmc.combined <- FindClusters(pbmc.combined, resolution = 1)
pbmc.combined <- RunUMAP(pbmc.combined, dims = 1:30, reduction = "integrated.cca")
# Visualization
DimPlot(pbmc.combined, reduction = "umap", group.by = c("orig.ident", "singleR.labels"))
#annotation using singleR
library(SingleR)
library(celldex)
ref <- celldex::HumanPrimaryCellAtlasData()
ref
View(as.data.frame(colData(ref)))
pbmc_counts <- GetAssayData(pbmc.combined, slot = 'counts')
pred <- SingleR(test = pbmc_counts,
ref = ref,
labels = ref$label.main)
pbmc.combined$singleR.labels <- pred$labels[match(rownames(pbmc.combined@meta.data), rownames(pred))]
DimPlot(pbmc.combined, reduction = "umap", split.by = "orig.ident")
Idents(pbmc.combined) <- "singleR.labels"
aggregatepbmc <- AggregateExpression(pbmc.combined, group.by = c("singleR.labels", "orig.ident"), return.seurat = TRUE)
pbmc.combined$cell_type <- paste(pbmc.combined$singleR.labels, pbmc.combined$orig.ident, sep = "-")
Idents(pbmc.combined) <-"cell_type"
response <- FindMarkers(pbmc.combined, ident.1 = "T_cells-car", ident.2 = "T_cells-acar", verbose= FALSE)
genest <- response[response$p_val < 0.05,]
so the genest file contains differentially expressed genes between car and acar conditions of T-cells, right?