测试文件及代码地址：点击下载 1. 读入并检查数据 library (Seurat) library (dplyr) pbmc <- readRDS ( './pbmcrenamed.rds' ) pbmc # An object of class Seurat # 22254 features across 900 samples within 2 assays # Active assay: RNA (20254 features, 0 variable features) # 2 layers present: counts, data # 1 other assay present: integrated # 3 dimensional reductions calculated: pca, umap, tsne DimPlot (pbmc) 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" unique (pbmc $ group) # [1] "C57" "AS1" "P3" DimPlot (pbmc, split.by = 'group' ) 2. 差异分析 pbmc $ celltype.group <- paste (pbmc $ celltype, pbmc $ group, sep = "_" ) pbmc $ celltype <- Idents (pbmc) Idents (pbmc) <- "celltype.group" table ( Idents (pbmc)) # # VSMC_C57 T cell_C57 Fibro_C57 EC_C57 # 63 32 48 59 # Myeloid cells_C57 B cell_C57 Mono_C57 Macro_C57 # 44 39 21 3 # Mono_AS1 T cell_AS1 Macro_AS1 Neut_AS1 # 49 45 25 28 # Fibro_AS1 Myeloid cells_AS1 B cell_AS1 EC_AS1 # 16 52 46 18 # VSMC_AS1 Macro_P3 Neut_P3 Mono_P3 # 8 72 30 # B cell_P3 T cell_P3 Fibro_P3 VSMC_P3 # 15 23 36 29 # EC_P3 Myeloid cells_P3 # 23 4 mydeg <- FindMarkers (pbmc, ident.1 = 'VSMC_AS1' , ident.2 = 'VSMC_C57' , verbose = FALSE , test.use = 'wilcox' , min.pct = 0.1 ) head (mydeg) # p_val avg_log2FC pct.1 pct.2 p_val_adj # Cd24a 6.327111e-07 4.931989 0.500 0.016 0.01281493 # Spta1 9.387127e-07 4.618852 0.375 0.000 0.01901269 # Lum 9.387127e-07 9.785188 0.375 0.000 0.01901269 # Gda 9.387127e-07 5.350293 0.375 0.000 0.01901269 # Isg20 6.651476e-06 4.669678 0.500 0.032 0.13471900 # Hbb-bt 7.937909e-06 4.454887 0.500 0.032 0.16077441 3. 解放生产力 通过循环自动计算差异基因 cellfordeg <- levels (pbmc $ celltype) cellfordeg # [1] "VSMC" "T cell" "Macro" "B cell" # [5] "Fibro" "Myeloid cells" "Mono" "Neut" # [9] "EC" for (i in 1 : length (cellfordeg)){ ident1 = paste0 (cellfordeg[i], "_P3" ) print (ident1) ident2 = paste0 (cellfordeg[i], "_AS1" ) print (ident2) CELLDEG <- FindMarkers (pbmc, ident.1 = ident1 , ident.2 = ident2, verbose = FALSE ) write.csv (CELLDEG, paste0 ( "../result_seuratV5/" ,cellfordeg[i], ".CSV" )) } # [1] "VSMC_P3" # [1] "VSMC_AS1" # [1] "T cell_P3" # [1] "T cell_AS1" # [1] "Macro_P3" # [1] "Macro_AS1" # [1] "B cell_P3" # [1] "B cell_AS1" # [1] "Fibro_P3" # [1] "Fibro_AS1" # [1] "Myeloid cells_P3" # [1] "Myeloid cells_AS1" # [1] "Mono_P3" # [1] "Mono_AS1" # [1] "Neut_P3" # [1] "Neut_AS1" # [1] "EC_P3" # [1] "EC_AS1" list.files ( "../result_seuratV5/" ) # [1] "B cell.CSV" "EC.CSV" "Fibro.CSV" # [4] "Macro.CSV" "Mono.CSV" "Myeloid cells.CSV" # [7] "Neut.CSV" "pbmc.rds" "T cell.CSV" # [10] "VSMC.CSV" 4. 可视化方法 library (dplyr) # CELLDEG为EC现在是"Neut_P3"和"Neut_AS1"的组间差异基因 CELLDEG <- read.table ( "../result_seuratV5/Neut.CSV" , header = T, sep = "," , row.names = 1 ) head (CELLDEG) # p_val avg_log2FC pct.1 pct.2 p_val_adj # Sptlc2 1.731937e-06 -1.6468877 0.222 0.750 0.03507865 # Ifnar1 2.754659e-06 -1.6050075 0.125 0.607 0.05579286 # Igkc 4.157411e-06 -0.5147022 0.042 0.429 0.08420420 # Slbp 4.318635e-05 -2.9794327 0.069 0.393 0.87469624 # Rabep2 5.839358e-05 -4.0548058 0.000 0.214 1.00000000 # Itsn1 8.990620e-05 -4.3010964 0.014 0.250 1.00000000 top10 <- CELLDEG %>% top_n ( n = 10 , wt = avg_log2FC) %>% row.names top10 # [1] "Trappc11" "Fmo2" "Nqo1" "Wdr7" "Clec4n" "Stfa3" # [7] "Cxcl1" "Csta2" "Tnfsf9" "Stfa1" pbmc <- ScaleData (pbmc, features = rownames (pbmc)) pbmc_sub <- subset (pbmc, idents = c ( "Neut_P3" , "Neut_AS1" )) DoHeatmap (pbmc_sub, features = top10, size= 3 ) VlnPlot (pbmc_sub, features = top10) FeaturePlot (pbmc_sub , features = top10) DotPlot (pbmc_sub , features = top10) 5. 提取表达量，用ggplot2 DIY一个箱线图 # mymatrix <- as.data.frame (pbmc @ assays $ RNA @ data) mymatrix2 <- t (mymatrix) %>% as.data.frame mymatrix2 $ celltype <- pbmc $ celltype mymatrix2[, ncol (mymatrix2) + 1 ] <- mymatrix2[, ncol (mymatrix2) + 1 ] <- pbmc $ group colnames (mymatrix2)[ ncol (mymatrix2)] <- "group" library (ggplot2) p1 <- ggplot2 :: ggplot (mymatrix2, aes ( x= celltype, y= Thbs1, fill= group)) + geom_boxplot ( alpha= 0.7 ) + scale_y_continuous ( name = "Expression" ) + scale_x_discrete ( name= "Celltype" ) + scale_fill_manual ( values = c ( 'deepskyblue' , 'orange' , 'pink' )) p1 # 另一种提取方法,借助log1p函数提取 Idents (pbmc) <- colnames (pbmc) mymatrix <- log1p ( AverageExpression (pbmc, verbose = FALSE ) $ RNA) mymatrix2 <- t (mymatrix) %>% as.data.frame mymatrix2 $ celltype <- pbmc $ celltype mymatrix2[, ncol (mymatrix2) + 1 ] <- pbmc $ group colnames (mymatrix2)[ ncol (mymatrix2)] <- "group" library (ggplot2) p2 <- ggplot2 :: ggplot (mymatrix2, aes ( x= celltype, y= Thbs1, fill= group)) + geom_boxplot ( alpha= 0.7 ) + scale_y_continuous ( name = "Expression" ) + scale_x_discrete ( name= "Celltype" ) + scale_fill_manual ( values = c ( 'deepskyblue' , 'orange' , 'pink' )) p2 library (patchwork) p1 | p2