image
忙里偷闲，整理了一套绘图代码，代码300行，效果如上图所示：其实就是把单细胞转录组中常见的气泡图根据分组信息展开了，代码目前适用于分组后level数目为2/3/4这三种情况。
如果不分组，就是普通的气泡图，之前已经讲过几次了：
- 单细胞分析实录(9): 展示marker基因的4种图形（二）
- 单细胞转录组绘图视频教程
绘图完整流程如下
1. seurat标准流程
library(Seurat)
library(tidyverse)
library(xlsx)
library(harmony)
### 导入数据 ########################################
testanno=readRDS("testanno.rds")
testmat=readRDS("testmat.rds")
### 标准流程&质控 ###################################
testseu = CreateSeuratObject(counts = testmat)
testseu[["percent.mt"]] <- PercentageFeatureSet(testseu, pattern = "^mt-")
testseu@meta.data$sample=testseu@meta.data$orig.ident
VlnPlot(testseu,features = c("nCount_RNA", "nFeature_RNA", "percent.mt"),group.by = "sample",pt.size = 0)
testseu <- subset(testseu, subset = nCount_RNA < 40000 &
nFeature_RNA < 5000 & nFeature_RNA > 500 &
percent.mt < 10)
testseu <- NormalizeData(testseu, normalization.method = "LogNormalize", scale.factor = 10000)
testseu <- FindVariableFeatures(testseu, selection.method = "vst", nfeatures = 2000)
testseu <- ScaleData(testseu)
testseu <- RunPCA(testseu, npcs = 50, verbose = FALSE)
### harmony去批次 ###################################
testseu=testseu %>% RunHarmony("sample", plot_convergence = TRUE)
### 降维&聚类 #######################################
testseu <- testseu %>%
RunUMAP(reduction = "harmony", dims = 1:20) %>%
FindNeighbors(reduction = "harmony", dims = 1:20) %>%
FindClusters(resolution = 0.5)
### 添加注释 ########################################
testseu@meta.data$CB=rownames(testseu@meta.data)
testseu@meta.data=testseu@meta.data%>%inner_join(testanno,by="CB")
rownames(testseu@meta.data)=testseu@meta.data$CB
DimPlot(testseu,reduction = "umap",group.by = "seurat_clusters",pt.size = 1,label = T,repel = T,label.size = 6)+
DimPlot(testseu,reduction = "umap",group.by = "celltype",pt.size = 1,label = T,repel = T,label.size = 6)
rm(list = c("testanno", "testmat"))
2. 准备差异基因
### 找每一群细胞的marker基因 ###################################################
testseu@meta.data$celltype=factor(testseu@meta.data$celltype,
levels = sort(unique(testseu@meta.data$celltype)))
Idents(testseu) = "celltype"
markerdf=FindAllMarkers(testseu,logfc.threshold = 0.8,only.pos = T)
markerdf=markerdf %>% group_by(cluster) %>% top_n(100,wt = avg_log2FC)
markerdf=as.data.frame(markerdf)
markerdf=markerdf%>%arrange(cluster,desc(avg_log2FC))
write.csv(markerdf,file = "celltype_marker_log2fc0.8_top100.csv",quote = F,row.names = F)
# 之后从这个表格中挑出想画图的基因，另存为plotgene.xlsx
plotgene=read.xlsx("plotgene.xlsx",sheetIndex = 1,header = T)
### 接下来为了得到最终图的效果，人为构建分组信息 ###############################
# 实际分析中，分组信息是不能人为随意构建的，而是课题样本内在具有的信息。
# 比如，本示例数据有7个样本: Arep1 Arep2 Arep3 Brep1 Brep2 Brep3 Brep4
# 显然可以分为两组: Arep和Brep
# 这里人为设定了3种分组，分别包含2/3/4个水平，以此来测试代码的适用性，这三种
# 应该涵盖了绝大多数用户的需求。
testseu@meta.data$group1=sample(c("groupA","groupB"),dim(testseu@meta.data)[1],replace = T)
testseu@meta.data$group2=sample(c("groupA","groupB","groupC"),dim(testseu@meta.data)[1],replace = T)
testseu@meta.data$group3=sample(c("groupA","groupB","groupC","groupD"),dim(testseu@meta.data)[1],replace = T)
testseu@meta.data$group3[testseu@meta.data$group3 == "groupD" &
testseu@meta.data$celltype == "cell_B"] = "groupC"
3. 绘制分组气泡图
# 必要的输入：
# seurat对象，meta数据框中有表示细胞类型的列和表示分组的列；
# 提供细胞类型的顺序；
# 表示marker基因的数据框，至少有两列，一列为cluster，一列为gene，基因名无重复。
source("bubble_split.R")
bubble_split(seu.obj = testseu,celltype.column = "celltype",
celltype.order = levels(Idents(testseu)),
group.column = "group1",deg = plotgene,
rect_color = "red",rect_size = 1)
ggsave("bubble_分两组.pdf",width = 40,height = 7,units = "cm")
bubble_split(seu.obj = testseu,celltype.column = "celltype",
celltype.order = levels(Idents(testseu)),
group.column = "group3",deg = plotgene,
rect_color = "black",rect_size = 0.5)
ggsave("bubble_分四组.pdf",width = 40,height = 11,units = "cm")
获取代码
淘，
宝店-
铺TOP生物信息
