手把手教你做单细胞测序数据分析（六）——组间差异分析及可视化

往期回顾：

在前面的课程中我们已经进行过“单样本数据分析”、“多样本数据整合”、“细胞类型注释”等内容的学习，相信大家现在已经能够对单细胞测序数据分析流程及Seurat对象的基本结构拥有了一定的了解。这一讲主要带领大家进行组间差异的计算及可视化方法的学习，这部分内容能够帮助科研工作者直接证明该数据集的前期试验设计，从前期枯燥的数据预处理走向文章中的Figure！

视频教程：

保姆级教程 《手把手教你做单细胞测序数据分析》（六）——组间差异分析及可视化

（B站同步播出，先看一遍视频再跟着代码一起操作，建议每个视频至少看三遍）

代码：
测试数据与第四讲多样本整合相同：

读入并检查数据

library(Seurat)
## Attaching SeuratObject
library(dplyr)
## 
## 载入程辑包：'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
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)
##  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')

差异分析

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[i],"_P3"), ident.2 = paste0(cellfordeg[i],"_AS1"), verbose = FALSE)
  write.csv(CELLDEG,paste0(cellfordeg[i],".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"     "Plvap"   
##  [7] "Igfbp7"   "Ifi27l2a" "Dcn"      "Gdf15"
pbmc <- ScaleData(pbmc, features =  rownames(pbmc))
## Centering and scaling data matrix
DoHeatmap(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$celltype
colnames(mymatrix2)[1] <- "celltype"

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

#########另一种提取方法########
Idents(pbmc) <- colnames(pbmc)
mymatrix <- log1p(AverageExpression(pbmc, verbose = FALSE)$RNA)
mymatrix2<-t(mymatrix)%>%as.data.frame()
mymatrix2[,1]<-pbmc$celltype
colnames(mymatrix2)[1] <- "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'))

###比较一下两种方法，发现并没有差异
library(patchwork)
p1|p2

本系列其他课程

手把手教你做单细胞测序数据分析（一）——绪论

手把手教你做单细胞测序数据分析（二）——各类输入文件读取

手把手教你做单细胞测序数据分析（三）——单样本分析

手把手教你做单细胞测序数据分析（四）——多样本整合

手把手教你做单细胞测序数据分析（五）——细胞类型注释

手把手教你做单细胞测序数据分析（六）——组间差异分析及可视化

手把手教你做单细胞测序数据分析（七）——基因集富集分析

欢迎关注同名公众号~