J
Jay2Echo
V1
2023/02/13阅读：165主题：默认主题
[单细胞|R]单细胞自动注释与手动注释
背景
上次分享了单细胞聚类分群的部分，本次一起来学习后续的基于Seurat的自动注释与手动注释。以及通过小提琴图、气泡图等多种形式进行Marker表达量展示。
[单细胞|R]单细胞降维聚类
Just Do It!
基于SingleR进行细胞自动注释
#载入参考数据集
ref_use <- HumanPrimaryCellAtlasData()
#SingleR注释
pred <- SingleR(test=as.matrix(pbmc@assays$RNA@data), #输入表达矩阵
ref=ref_use, #参考数据
labels=ref_use$label.fine, #标签列
clusters = pbmc$seurat_clusters, #细胞聚类信息（只有method = "cluster"才使用）
method = "cluster") #注释类型，按照cluster还是按照cell进行注释
#保存SingleR注释结果,需要根据自己情况设置文件位置
write.table(pred,"hg19/Singler_out.xls",sep='\t',quote = F)
注：
-
本次使用的参考数据集是HumanPrimaryCellAtlasData()，是人类参考数据集。 -
此参考数据集有37个大类，157个小类 -
通过unique(ref_use$label.main)查看所有大类
-
-
通过unique(ref_use$label.fine)查看所有小类 -
SingleR（）参数意义： -
目前当'cluster='参数被指定时，method="cluster"参数可以缺省，即省略不写。
-
-
通过“head(pred)”查看自动注释结果。 -
通过“unique(pred$pruned.labels)”查看自动注释得到的细胞类别。 -
发现原始Cluster有0~7,8种，自动注释结果只有7种，查看输出的结果文件，发现Cluster0与Cluster2均被识别为同一种细胞类型。故后续需要选择差异基因中高表达基因作为Marker进行手动注释。 
-
基于CellMarker数据库进行手动注释
-
通过FindAllMarkers（）函数可得每个Cluster中的差异表达基因。利用差异表达基因Top10表达热图从每个Cluster中寻找高表达的基因作为Marker进行后续手动注释。
-
基于Marker进入CellMarker官网进行手动注释。 -
该网站可以实现两种注释逻辑： -
通过Marker匹配细胞类型 -
通过细胞类型匹配Marker
-
-
通过手动注释提醒我们，需要关注原始数据的物种及取样组织信息。
通过Marker匹配细胞类型
-
通过点击首页右上角“search”按钮，进入下一个功能页面。 -
依次选择物种、Marker、点击Submit即可
通过细胞类型匹配Marker
-
通过依次选择物种、组织部位、细胞类型之后可以得到以词云形式呈现的所有Markers信息。 -
同样词云中字体越大表示被报道次数越多，越可信。 
自动及手动结果汇总
#细胞注释
new.cluster.ids <- c("Memory CD4+ T","Naive CD4+ T","CD14+ Mono","B cell","CD8+ T",
"FCGR3A+ Mono","NK", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids) #修改Idents
#在metadata中，添加Celltype信息
pbmc$celltype <- Idents(pbmc)
head(pbmc@meta.data)
#降维图
DimPlot(pbmc, #seurat对象(数据)
reduction = 'umap', #降维类型，umap，tsne，pca
group.by = 'celltype', #分组名称
pt.size = 1, #图中点的大小
#split.by = 'seurat_clusters', #拆分绘图的名称
label = T) #图中是否添加label信息
#保存分析结果,需要根据自己情况设置文件位置
saveRDS(pbmc, file = "hg19/pbmc3k_final.rds")
注：
-
创建new.cluster.ids对象时细胞类型需要按照Cluster0~7顺序指定。
marker展示
每个Cluster前Top10差异基因表达量热图
#DoHeatmap热图
DoHeatmap(pbmc, #seurat对象(数据)
features = top10$gene, #绘制的基因
label = F, #图中是否添加label信息
slot = "scale.data", #绘图使用的表达矩阵
group.by = 'seurat_clusters', #分组名称
group.bar = T) #是否展示bar
单Marker--CD79A表达情况展示
#FeaturePlot图
FeaturePlot(pbmc, #seurat对象(数据)
features = "CD79A", #绘制的基因
reduction = "umap", #降维类型
pt.size = 1, #图中点的大小
cols = c("lightgrey","red"), #点的颜色设置
label = T, #图中是否添加label信息
order = T) #是否按照表达量顺序对细胞进行排序
单Marker--CD8A表达小提琴图展示
#小提琴图
VlnPlot(pbmc, #seurat对象(数据)
features = "CD8A", #绘制的基因
pt.size = 0.5, #图中点的大小
group.by = "seurat_clusters", #分组名称
cols = c("#FF6700","#9370DB","#89CFF0", #绘图颜色
"#F8D568","#00AD43","#BA160C",
"#FF91AF","#A6A6A6"))
单Marker--CD8A表达气泡图展示
# 气泡图
DotPlot(pbmc,features = "CD8A")
多Marker展示汇总
-
同时观测多个Marker在不同Cluster中表达量情况的需求也很常见，直接上代码。 -
通过比较单个与多个Marker的结果图发现，两者代码完全一致，只是作图函数中features参数由原来的单个marker替换为一个多markers的列表。
#Featureplot图2
cell_marker <- c("IL7R","CCR7","S100A4","CD14","LYZ","MS4A1","CD8A","FCGR3A","MS4A7","GNLY","NKG7","PPBP")
FeaturePlot(pbmc,features = cell_marker)
#小提琴图2
VlnPlot(pbmc, #seurat对象(数据)
features = cell_marker, #绘制的基因
pt.size = 0.5, #图中点的大小
group.by = "seurat_clusters", #分组名称
stack = T, #是否将所有基因的图片拼合一起
same.y.lims = T, #是否设置相同的y轴
flip = F) + #是否对坐标轴翻转
NoLegend() #不绘制图例
#气泡图2
DotPlot(pbmc, #seurat对象(数据)
features = cell_marker, #绘制的基因
group.by = "seurat_clusters") #分组名称
数据及All_Code
数据
-
见本文引用的第一篇推文末百度网盘
All_Code
############## SingleR注释 #########################
#载入参考数据集
ref_use <- HumanPrimaryCellAtlasData()
# 此参考数据集有37个大类，157个小类
# 通过unique(ref_use$label.main)查看所有大类
# 通过unique(ref_use$label.fine)查看所有小类
#SingleR注释
pred <- SingleR(test=as.matrix(pbmc@assays$RNA@data), #输入表达矩阵
ref=ref_use, #参考数据
labels=ref_use$label.fine, #标签列
clusters = pbmc$seurat_clusters, #细胞聚类信息（只有method = "cluster"才使用）
method = "cluster") #注释类型，按照cluster还是按照cell进行注释
head(pred)
# 查看注释后的细胞类型
# unique(pred$pruned.labels)
#保存SingleR注释结果,需要根据自己情况设置文件位置
write.table(pred,"hg19/Singler_out.xls",sep='\t',quote = F)
#细胞注释
new.cluster.ids <- c("Memory CD4+ T","Naive CD4+ T","CD14+ Mono","B cell","CD8+ T",
"FCGR3A+ Mono","NK", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids) #修改Idents
#在metadata中，添加Celltype信息
pbmc$celltype <- Idents(pbmc)
head(pbmc@meta.data)
#降维图
DimPlot(pbmc, #seurat对象(数据)
reduction = 'umap', #降维类型，umap，tsne，pca
group.by = 'celltype', #分组名称
pt.size = 1, #图中点的大小
#split.by = 'seurat_clusters', #拆分绘图的名称
label = T) #图中是否添加label信息
#保存分析结果,需要根据自己情况设置文件位置
saveRDS(pbmc, file = "hg19/pbmc3k_final.rds")
#DoHeatmap热图
DoHeatmap(pbmc, #seurat对象(数据)
features = top10$gene, #绘制的基因
label = F, #图中是否添加label信息
slot = "scale.data", #绘图使用的表达矩阵
group.by = 'seurat_clusters', #分组名称
group.bar = T) #是否展示bar
#FeaturePlot图
FeaturePlot(pbmc, #seurat对象(数据)
features = "CD79A", #绘制的基因
reduction = "umap", #降维类型
pt.size = 1, #图中点的大小
cols = c("lightgrey","red"), #点的颜色设置
label = T, #图中是否添加label信息
order = T) #是否按照表达量顺序对细胞进行排序
#小提琴图
VlnPlot(pbmc, #seurat对象(数据)
features = "CD8A", #绘制的基因
pt.size = 0.5, #图中点的大小
group.by = "seurat_clusters", #分组名称
cols = c("#FF6700","#9370DB","#89CFF0", #绘图颜色
"#F8D568","#00AD43","#BA160C",
"#FF91AF","#A6A6A6"))
# 气泡图
DotPlot(pbmc,features = "CD8A")
#Featureplot图2
cell_marker <- c("IL7R","CCR7","S100A4","CD14","LYZ","MS4A1","CD8A","FCGR3A","MS4A7","GNLY","NKG7","PPBP")
FeaturePlot(pbmc,features = cell_marker)
#小提琴图2
VlnPlot(pbmc, #seurat对象(数据)
features = cell_marker, #绘制的基因
pt.size = 0.5, #图中点的大小
group.by = "seurat_clusters", #分组名称
stack = T, #是否将所有基因的图片拼合一起
same.y.lims = T, #是否设置相同的y轴
flip = F) + #是否对坐标轴翻转
NoLegend() #不绘制图例
#气泡图2
DotPlot(pbmc, #seurat对象(数据)
features = cell_marker, #绘制的基因
group.by = "seurat_clusters") #分组名称
其他
后续会更新富集分析相关内容，敬请期待！
-
推文多平台同步发布，公众号内容食用更佳 -
更多内容，请关注微信公众号“生信矿工” -
如有意见或建议可以在评论区讨论
作者介绍
J
Jay2Echo
V1