# 首先读入原始文件 library (Seurat) pbmc <- readRDS ( './pbmcrenamed.rds' ) # 查看metadata里存储了哪些信息 head (pbmc @ meta.data) # orig.ident nCount_RNA nFeature_RNA percent.mt group # AAAGGGCTCGCAGAGA-1_1 C57 15155 3673 4.988453 C57 # AACACACAGCATCCCG-1_1 C57 6300 1995 7.174603 C57 # AACAGGGTCAACACCA-1_1 C57 21299 4453 2.427344 C57 # AACCCAAAGGCCTTGC-1_1 C57 1160 688 3.706897 C57 # AACCCAAAGTGGACGT-1_1 C57 1036 387 1.833977 C57 # AACCCAAGTCCGAAGA-1_1 C57 3480 1401 2.126437 C57 # integrated_snn_res.0.025 seurat_clusters celltype.group # AAAGGGCTCGCAGAGA-1_1 0 0 VSMC_C57 # AACACACAGCATCCCG-1_1 1 1 T cell_C57 # AACAGGGTCAACACCA-1_1 4 4 Fibro_C57 # AACCCAAAGGCCTTGC-1_1 8 8 EC_C57 # AACCCAAAGTGGACGT-1_1 5 5 Myeloid cells_C57 # AACCCAAGTCCGAAGA-1_1 1 1 T cell_C57 # celltype # AAAGGGCTCGCAGAGA-1_1 VSMC # AACACACAGCATCCCG-1_1 T cell # AACAGGGTCAACACCA-1_1 Fibro # AACCCAAAGGCCTTGC-1_1 EC # AACCCAAAGTGGACGT-1_1 Myeloid cells # AACCCAAGTCCGAAGA-1_1 T cell names (pbmc @ meta.data) # [1] "orig.ident" "nCount_RNA" # [3] "nFeature_RNA" "percent.mt" # [5] "group" "integrated_snn_res.0.025" # [7] "seurat_clusters" "celltype.group" # [9] "celltype" # 查看初始分类变量 head (pbmc @ active.ident) # AAAGGGCTCGCAGAGA-1_1 AACACACAGCATCCCG-1_1 AACAGGGTCAACACCA-1_1 # VSMC T cell Fibro # AACCCAAAGGCCTTGC-1_1 AACCCAAAGTGGACGT-1_1 AACCCAAGTCCGAAGA-1_1 # EC Myeloid cells T cell # Levels: VSMC T cell Macro B cell Fibro Myeloid cells Mono Neut EC table (pbmc @ active.ident) # # VSMC T cell Macro B cell Fibro # 100 # Myeloid cells Mono Neut EC # 100 table ( Idents (pbmc)) # # VSMC T cell Macro B cell Fibro # 100 # Myeloid cells Mono Neut EC # 100 # 修改初始分类变量为group Idents (pbmc) <- "group" table (pbmc @ active.ident) # # C57 AS1 P3 # 309 287 304 table ( Idents (pbmc)) # # C57 AS1 P3 # 309 287 304 unique (pbmc @ active.ident) # [1] C57 AS1 P3 # Levels: C57 AS1 P3 unique (pbmc $ group) # [1] "C57" "AS1" "P3" # 筛选部分分组 unique (pbmc $ group)[ 1 : 2 ] # [1] "C57" "AS1" pbmc_sub <- subset (pbmc, idents = unique (pbmc $ group)[ 1 : 2 ]) # 可以看看subset前后的变化，metadata里的变量都可以这样取子集 table ( Idents (pbmc_sub)) # # C57 AS1 # 309 287 # 以group列作为分类标签进行数据展示 # 计算线粒体比例，小鼠线粒体基因开头为mt- pbmc_sub[[ "percent.mt" ]] <- PercentageFeatureSet (pbmc_sub, pattern = "^mt-" ) # VlnPlot展示细胞质量 VlnPlot (pbmc_sub, features = c ( "nFeature_RNA" , "nCount_RNA" , "percent.mt" ), ncol = 3 ) # 过滤低质量的细胞 pbmc_sub <- subset (pbmc_sub, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 & percent.mt < 10 ) # 过滤后再次展示细胞质量 VlnPlot (pbmc_sub, features = c ( "nFeature_RNA" , "nCount_RNA" , "percent.mt" ), ncol = 3 ) # 数据标准化 suppressWarnings ( suppressMessages ({ pbmc_sub <- NormalizeData (pbmc_sub, normalization.method = "LogNormalize" , scale.factor = 10000 ) # 计算高变基因 pbmc_sub <- FindVariableFeatures (pbmc_sub, selection.method = "vst" , nfeatures = 2000 ) pbmc_sub <- ScaleData (pbmc_sub, features = rownames (pbmc_sub)) # 剔除不想要的变量（如线粒体的比例） pbmc_sub <- ScaleData (pbmc_sub, vars.to.regress = "percent.mt" ) # PCA降维分析 pbmc_sub <- RunPCA (pbmc_sub, features = VariableFeatures ( object = pbmc_sub)) # 计算细胞临近关系 pbmc_sub <- FindNeighbors (pbmc_sub, dims = 1 : 20 ) pbmc_sub <- RunUMAP (pbmc_sub, dims = 1 : 20 ) })) DimPlot (pbmc_sub, label = T) # 查看[email protected]包含的信息 names (pbmc_sub @ meta.data) # [1] "orig.ident" "nCount_RNA" # [3] "nFeature_RNA" "percent.mt" # [5] "group" "integrated_snn_res.0.025" # [7] "seurat_clusters" "celltype.group" # [9] "celltype" # 查看pbmc_sub每个细胞数目 table (pbmc_sub $ celltype) # # VSMC T cell Macro B cell Fibro # 60 76 22 85 56 # Myeloid cells Mono Neut EC # 95 64 28 63 # 查看pbmc_sub初始分类变量 table ( Idents (pbmc_sub)) # # C57 AS1 # 286 263 # 修改pbmc_sub的分类变量为celltype Idents (pbmc_sub) <- "celltype" table ( Idents (pbmc_sub)) # # VSMC T cell Macro B cell Fibro # 60 76 22 85 56 # Myeloid cells Mono Neut EC # 95 64 28 63 DimPlot (pbmc_sub, label = T) # 修改pbmc_sub的分类变量为celltype后，调整分类变量里面的元素名称 library (tidyverse) pbmc_sub $ new <- pbmc_sub $ group unique (pbmc_sub $ new) # [1] "C57" "AS1" pbmc_sub $ new <- recode (pbmc_sub $ new, "C57" = "sample1" , "AS1" = "sample2" ) Idents (pbmc_sub) <- "new" table ( Idents (pbmc_sub)) # # sample1 sample2 # 286 263 DimPlot (pbmc_sub) # 根据分类变量group拆分seurat对象 pbmc_sub.list <- SplitObject (pbmc_sub, split.by = 'group' ) C57 <- pbmc_sub.list[[ 1 ]] AS1 <- pbmc_sub.list[[ 2 ]] # 查看拆分后的对象类别 class (C57) # [1] "Seurat" # attr(,"package") # [1] "SeuratObject" class (AS1) # [1] "Seurat" # attr(,"package") # [1] "SeuratObject" # 拆分后的对象都添加一个新的类别 C57 $ group <- "group1" AS1 $ group <- "group2" # 重新合并两个seurat对象 pbmc_sub_new <- merge ( x= C57, y= AS1) table ( Idents (pbmc_sub_new)) # # sample1 sample2 # 286 263 unique (pbmc_sub_new $ group) # [1] "group1" "group2" # 调整分类变量为group Idents (pbmc_sub_new) <- "group" table ( Idents (pbmc_sub_new)) # # group1 group2 # 286 263 # 小提琴图展示基因表达量，同时分类变量是我们刚刚设置的group VlnPlot (pbmc_sub_new, features = "Mndal" , pt.size = 0 )