往期回顾：
在前面的课程中我们已经举行过“单样本数据分析”、“多样本数据整合”、“细胞类型解释”等内容的学习，信任各人现在已经能够对单细胞测序数据分析流程及Seurat对象的根本布局拥有了肯定的相识。这一讲紧张领导各人举行组间差异的盘算及可视化方法的学习，这部分内容能够资助科研工作者直接证明该数据集的前期试验计划，从前期枯燥的数据预处置惩罚走向文章中的Figure！
视频教程：
保姆级教程 《手把手教你做单细胞测序数据分析》（六）——组间差异分析及可视化
（B站同步播出，先看一遍视频再跟着代码一起利用，发起每个视频至少看三遍）
代码：
测试数据与第四讲多样本整合雷同：
读入并查抄数据
library(Seurat)## Attaching SeuratObjectlibrary(dplyr)## ## 载入程辑包：'dplyr'## The following objects are masked from 'package:stats':## ## filter, lag## The following objects are masked from 'package:base':## ## intersect, setdiff, setequal, unionpbmc <- 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)## 1 other assay present: integrated## 3 dimensional reductions calculated: pca, umap, tsneDimPlot(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" " 3"DimPlot(pbmc,split.by = 'group')
差异分析
pbmc$celltype.group <- paste(pbmc$celltype, pbmc$group, sep = "_")pbmc$celltype <- Idents(pbmc)Idents(pbmc) <- "celltype.group"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 1.4046048 0.500 0.016 0.01281493## Spta1 9.387127e-07 0.3391453 0.375 0.000 0.01901269## Lum 9.387127e-07 3.8953383 0.375 0.000 0.01901269## Gda 9.387127e-07 0.6064680 0.375 0.000 0.01901269## Isg20 6.651476e-06 1.4016408 0.500 0.032 0.13471900## Hbb-bt 7.937909e-06 4.3779094 0.500 0.032 0.16077441解放生产力 通过循环自动盘算差异基因
cellfordeg<-levels(pbmc$celltype)for(i in 1:length(cellfordeg)){ CELLDEG <- FindMarkers(pbmc, ident.1 = paste0(cellfordeg,"_P3"), ident.2 = paste0(cellfordeg,"_AS1"), verbose = FALSE) write.csv(CELLDEG,paste0(cellfordeg,".CSV"))}list.files()## [1] "B cell.CSV" "EC.CSV" ## [3] "Fibro.CSV" "Macro.CSV" ## [5] "Mono.CSV" "Myeloid cells.CSV" ## [7] "Neut.CSV" "pbmcrenamed.rds" ## [9] "T cell.CSV" "VSMC.CSV" ## [11] "组间差异分析及可视化.html" "组间差异分析及可视化.Rmd" ## [13] "组间差异分析及可视化_files" "组间差异及可视化.r"差异分析解果解读：
可视化方法
library(dplyr)top10 <- CELLDEG %>% top_n(n = 10, wt = avg_log2FC) %>% row.names()top10## [1] "Thbs1" "Acta2" "Myl9" "Tagln" "Ccn2" "lvap" ## [7] "Igfbp7" "Ifi27l2a" "Dcn" "Gdf15"pbmc <- ScaleData(pbmc, features = rownames(pbmc))## Centering and scaling data matrixDoHeatmap(pbmc,features = top10,size=3)
Idents(pbmc) <- "celltype"VlnPlot(pbmc,features = top10,split.by = 'group',idents = 'EC')## The default behaviour of split.by has changed.## Separate violin plots are now plotted side-by-side.## To restore the old behaviour of a single split violin,## set split.plot = TRUE.## ## This message will be shown once per session.FeaturePlot(pbmc,features = top10,split.by = 'group')#DotPlot(pbmc,features = top10,split.by ='group')#默认只有两种颜色DotPlot(pbmc,features = top10,split.by ='group',cols = c('blue','yellow','pink'))
提取表达量，用ggplot2 DIY一个箱线图
####提取表达量#######mymatrix <- as.data.frame(pbmc@assays$RNA@data)mymatrix2<-t(mymatrix)%>%as.data.frame()mymatrix2[,1]<-pbmc$celltypecolnames(mymatrix2)[1] <- "celltype"mymatrix2[,ncol(mymatrix2)+1]<-pbmc$groupcolnames(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#########另一种提取方法########Idents(pbmc) <- colnames(pbmc)mymatrix <- log1p(AverageExpression(pbmc, verbose = FALSE)$RNA)mymatrix2<-t(mymatrix)%>%as.data.frame()mymatrix2[,1]<-pbmc$celltypecolnames(mymatrix2)[1] <- "celltype"mymatrix2[,ncol(mymatrix2)+1]<-pbmc$groupcolnames(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'))
###比较一下两种方法，发现并没有差异library(patchwork)p1|p2往期回顾
单细胞数据分析系列教程：
B站视频，先看一遍视频再去看推送利用，发起至少看三遍：https://www.bilibili.com/video/BV1S44y1b76Z/
单细胞测序根本数据分析保姆级教程，代码部分整理在往期推送之中：
手把手教你做单细胞测序数据分析（一）——绪论
手把手教你做单细胞测序数据分析（二）——各类输入文件读取
手把手教你做单细胞测序数据分析（三）——单样天职析
手把手教你做单细胞测序数据分析（四）——多样本整合
手把手教你做单细胞测序数据分析（五）——细胞类型解释
欢迎关注同名公众号~ |
