单细胞Day6-marker基因和注释
1.加载数据
直接加载之前的obj.Rdata
rm(list = ls())
library(Seurat)
load("obj.Rdata")
p1 = DimPlot(seu.obj,reduction = "umap",label = T)
p1
2.计算makergene
和差异基因里面的上调基因有点类似,某个基因在某一簇细胞里表达量都很高,在其他簇表达量很低,那么这个基因就是这簇细胞的象征。通常会设置 only.pos = TRUE,只关注上调。
2.1 FindAllMarkers
library(dplyr)
f = "markers.Rdata"
if(!file.exists(f)){
markers <- FindAllMarkers(seu.obj, only.pos = TRUE,min.pct = 0.25)
save(markers,file = f)
}
load(f)
FindAllMarkers:计算出所有细胞簇的marker基因,比如计算0这一簇时,就是0 vs 剩下的1-10,1-10全部算一个组。
min.pct = 0.25:一个基因至少要在25%的细胞中表达。还是以0这一簇为例,就是说要么在0这一簇里25%以上的细胞中表达,要么在剩下的1-10簇里(1-10全部算一个组)25%以上的细胞中表达。这样可以去掉一些表达比例太少的基因,比例太少的话不能代表这个簇的整体情况,即使logFC很大也没有意义。
导出结果:
Calculating cluster 0
For a (much!) faster implementation of the Wilcoxon Rank Sum Test,
(default method for FindMarkers) please install the presto package
--------------------------------------------
install.packages('devtools')
devtools::install_github('immunogenomics/presto')
--------------------------------------------
After installation of presto, Seurat will automatically use the more
efficient implementation (no further action necessary).
This message will be shown once per session
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03s
Calculating cluster 1
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02s
Calculating cluster 2
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02s
Calculating cluster 3
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02s
Calculating cluster 4
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=06s
Calculating cluster 5
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=00s
Calculating cluster 6
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=12s
Calculating cluster 7
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=04s
Calculating cluster 8
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=08s
Calculating cluster 9
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=05s
Calculating cluster 10
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=15s
Calculating cluster 11
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03s
2.2 取每个细胞簇logFC排名前n的marker基因
mks = markers %>% group_by(cluster) %>% top_n(n = 3, wt = avg_log2FC) #n=3就是每个簇取前3
g = unique(mks$gene) #去重复,有一些基因在不同簇中都被计为marker基因,排名都比较靠前。但有重复的话后面画图就会出错
3.makergene的五种可视化
3.1 热图
library(ggplot2)
DoHeatmap(seu.obj, features = g) + NoLegend()+
scale_fill_gradientn(colors = c("#2fa1dd", "white", "#f87669"))
3.2 气泡图
红色代表高表达,气泡的大小代表基因在某一簇里的表达比例
DotPlot(seu.obj, features = g,cols = "RdYlBu") +
RotatedAxis()
3.3 小提琴图
VlnPlot(seu.obj, features = g[1:6])
#features = g[1:6]代表画出整个列表的前6个基因,不能画太多
3.4 Featureplot
本质上还是umap图,只是换了一种着色方式,一张图对应一个基因,展示了该基因在所有细胞里的表达情况。
FeaturePlot(seu.obj, features = g[1:6])
#features = g[1:6]代表画出整个列表的前6个基因
3.5 峰峦图/山脊图
就是密度图,横坐标是基因表达量,最多的其实还是0,这种图用的比较少。
RidgePlot(seu.obj, features = g[1:6])
#features = g[1:6]代表画出整个列表的前6个基因
4.注释亚群
手动注释是金标准,但费时费力费眼睛,能注释到什么程度基本上取决于背景知识。自动注释一气呵成,方便,但也有明显的局限性:① 没有把组织、疾病类型等背景知识考虑在内 ②现成的参考数据只有人类和小鼠的 ③ 切换不同的参考数据就会给出不同的注释结果
4.1 手动注释
手动注释需要自行查找每种可能的细胞类型已知的marker基因,可以从文献里查,可以从数据库里查。
常用的三个数据库:
http://biocc.hrbmu.edu.cn/CellMarker
https://panglaodb.se/
https://www.gsea-msigdb.org/gsea/msigdb/human/genesets.jsp?collection=C8
存储在txt文件,每一行放一个细胞名和一个marker基因名,如果一种细胞有多个marker基因就拆成多行。
a = read.table("my_markers.txt",sep = ",")
gt = split(a[,2],a[,1]) #拆分为列表的格式
DotPlot(seu.obj, features = gt,cols = "RdYlBu") +
RotatedAxis()
新建一个anno.txt,两列数据,第一列是每簇细胞,第二列是给它注释的名字
不用按照0-11的顺序填,先看每一行又大又红的点,对应的细胞是什么。没有又大又红的只能退而求其次,看最大的点或者是最红的点。有时最大的和最红的不是同一个,就为难了,没有标准答案,我习惯是优先按照最红的,除非最红的太小了才按照最大的。。。比如1这一簇,最红的是CD8 T,最大的是NK,注释为CD8 T。
#给seurat object 重新规定Idents,RenameIdents实现修改分类
new.cluster.ids <- celltype$V2
names(new.cluster.ids) <- levels(seu.obj)
sce <- RenameIdents(seu.obj, new.cluster.ids)
save(sce,file = "sce.Rdata")
p2 <- DimPlot(sce, reduction = "umap",
label = TRUE, pt.size = 0.5) + NoLegend()
p1+p2
手动注释要检查是否合理,比如T细胞和NK细胞应该是在一起的,B细胞应该是单独的。
4.2 自动注释
最常用的自动注释R包是SingleR。它使用的参考数据是celldex里的
library(celldex)
library(SingleR)
ls("package:celldex")
每个参考数据的介绍要看它的帮助文档(例如?BlueprintEncodeData)
可以下载对应的Rdata并放在工作目录下。
f = "ref_BlueprintEncode.RData"
if(!file.exists(f)){
ref <- celldex::BlueprintEncodeData()
save(ref,file = f)
}
ref <- get(load(f))
切换时参考数据注意,有两处”BlueprintEncodeData”,都换成别的数据名字
输入数据是单细胞表达矩阵;注释完就会给出每个簇对应的细胞类型;然后给表达矩阵改分类、画图。
library(BiocParallel)
scRNA = seu.obj
test = scRNA@assays$RNA$data
pred.scRNA <- SingleR(test = test,
ref = ref,
labels = ref$label.main,
clusters = [email protected])
pred.scRNA$pruned.labels
new.cluster.ids <- pred.scRNA$pruned.labels
names(new.cluster.ids) <- levels(scRNA)
scRNA <- RenameIdents(scRNA,new.cluster.ids)
p3 <- DimPlot(scRNA, reduction = "umap",label = T,pt.size = 0.5) + NoLegend()
p2+p3
共有 0 条评论