作图

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

2022-07-19  本文已影响0人  Biomamba生信基地

往期回顾:

在前面的课程中我们已经进行过“单样本数据分析”、“多样本数据整合”、“细胞类型注释”等内容的学习,相信大家现在已经能够对单细胞测序数据分析流程及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)
640.png
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')
640.png

差异分析

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"

差异分析解果解读:

640.png

可视化方法

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)
640.png
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.
640.png
FeaturePlot(pbmc,features = top10,split.by = 'group')
640.png
#DotPlot(pbmc,features = top10,split.by ='group')#默认只有两种颜色
DotPlot(pbmc,features = top10,split.by ='group',cols = c('blue','yellow','pink'))
640.png

提取表达量,用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
640.png
#########另一种提取方法########
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'))
640.png
###比较一下两种方法,发现并没有差异
library(patchwork)
p1|p2
640.png

往期回顾

单细胞数据分析系列教程:

B站视频,先看一遍视频再去看推送操作,建议至少看三遍:https://www.bilibili.com/video/BV1S44y1b76Z/

单细胞测序基础数据分析保姆级教程,代码部分整理在往期推送之中:

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

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

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

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

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

欢迎关注同名公众号~

上一篇下一篇

猜你喜欢

热点阅读