任意细胞亚群的差异分析

2021-05-18 12:38:00 浏览数 (1)

分享是一种态度

我们以 seurat 官方教程为例:

代码语言:javascript复制
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') ,我们后面的教程都是依赖于这个 文件哦!

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

代码如下:

代码语言:javascript复制

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 ) 

如下所示,

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

代码语言:javascript复制

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) 

如下所示:

任意划分亚群再找差异

代码语言:javascript复制
# 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

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

代码语言:javascript复制
> 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

然后再找差异:

代码语言:javascript复制
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 ) 

如下所示:


0 人点赞