任意细胞亚群的差异分析

生信技能树jimmy

发布于 2021-05-18 12:38:00

1.1K0

发布于 2021-05-18 12:38:00

文章被收录于专栏：单细胞天地单细胞天地

分享是一种态度

我们以 seurat 官方教程为例：

rm(list = ls())
library(Seurat)
# devtools::install_github('satijalab/seurat-data')
library(SeuratData)
library(ggplot2)
library(patchwork)
library(dplyr)
load(file = 'basic.sce.pbmc.Rdata')

DimPlot(pbmc, reduction = 'umap', 
        label = TRUE, pt.size = 0.5) + NoLegend()

sce=pbmc

如果你不知道basic.sce.pbmc.Rdata 这个文件如何得到的，麻烦自己去跑一下 可视化单细胞亚群的标记基因的5个方法，自己 save(pbmc,file = 'basic.sce.pbmc.Rdata') ，我们后面的教程都是依赖于这个文件哦！

对各个细胞亚群找高表达量的标记基因

代码如下：


if (file.exists('sce.markers.all_10_celltype.Rdata')) {
  load('sce.markers.all_10_celltype.Rdata')
}else {
  sce.markers <- FindAllMarkers(object = sce, only.pos = TRUE,
                                min.pct = 0.25, 
                                thresh.use = 0.25)
  save(sce.markers,file = 'sce.markers.all_10_celltype.Rdata')
  
}

# 并且可视化它


head(sce.markers)
table(sce.markers$cluster)
# 首先挑选基因
kp=grepl('Mono',sce.markers$cluster)
table(kp)
cg_sce.markers = sce.markers [ kp ,]

# 然后挑选细胞
kp=grepl('Mono', Idents(sce ) )
table(kp)
sce=sce[,kp]
sce
table( Idents(sce ))

cg_sce.markers=cg_sce.markers[cg_sce.markers$avg_logFC>2,]
dim(cg_sce.markers)
DoHeatmap(subset(sce, downsample = 15),
          unique(cg_sce.markers$gene),
          slot = 'counts',
          size=3 )

如下所示，

对指定的两个细胞亚群找差异


levels(Idents(sce))
markers_df <- FindMarkers(object = sce, 
                          ident.1 = 'FCGR3A+ Mono',
                          ident.2 = 'CD14+ Mono',
                          #logfc.threshold = 0,
                          min.pct = 0.25)
head(markers_df)
cg_markers_df=markers_df[abs(markers_df$avg_logFC) >1,]
dim(cg_markers_df)
DoHeatmap(subset(sce, downsample = 15),
          slot = 'counts',
          unique(rownames(cg_markers_df)),size=3)

如下所示：

任意划分亚群再找差异

# drop-out
highCells= colnames(subset(x = sce, subset = FCGR3A > 1,
                           slot = 'counts')) 
highORlow=ifelse(colnames(sce) %in% highCells,'high','low')
table(highORlow)
table(Idents(sce))
sce@meta.data$highORlow=highORlow

可以看到两个方法的划分亚群对比情况：

> table(highORlow)
highORlow
high  low 
 160  482 
> table(Idents(sce))

  CD14+ Mono FCGR3A+ Mono 
         480          162 
         
> table(Idents(sce),highORlow)
              highORlow
               high low
  CD14+ Mono     15 465
  FCGR3A+ Mono  145  17

然后再找差异：

markers <- FindMarkers(sce, ident.1 = "high", 
                       group.by = 'highORlow' )
head(x = markers)
cg_markers=markers[abs(markers$avg_logFC) >1,]
dim(cg_markers)
DoHeatmap(subset(sce, downsample = 15),
          rownames(cg_markers) ,
          size=3 )

如下所示：