单细胞Day6-marker基因和注释

1.加载数据

直接加载之前的obj.Rdata

rm(list = ls())
library(Seurat)
load("obj.Rdata")
p1 = DimPlot(seu.obj,reduction = "umap",label = T)
p1
plot_zoom4.png

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"))
heatmap.png
3.2 气泡图

红色代表高表达,气泡的大小代表基因在某一簇里的表达比例

DotPlot(seu.obj, features = g,cols = "RdYlBu") +
  RotatedAxis()
气泡图.png
3.3 小提琴图
VlnPlot(seu.obj, features = g[1:6])   
#features = g[1:6]代表画出整个列表的前6个基因,不能画太多
violin.png
3.4 Featureplot

本质上还是umap图,只是换了一种着色方式,一张图对应一个基因,展示了该基因在所有细胞里的表达情况。

FeaturePlot(seu.obj, features = g[1:6]) 
#features = g[1:6]代表画出整个列表的前6个基因
featureplot.png

3.5 峰峦图/山脊图

就是密度图,横坐标是基因表达量,最多的其实还是0,这种图用的比较少。

RidgePlot(seu.obj, features = g[1:6])  
#features = g[1:6]代表画出整个列表的前6个基因
山脊图.png

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基因就拆成多行。

image.png
a = read.table("my_markers.txt",sep = ",")
gt = split(a[,2],a[,1])   #拆分为列表的格式
DotPlot(seu.obj, features = gt,cols = "RdYlBu") +
  RotatedAxis()
注释气泡.png

新建一个anno.txt,两列数据,第一列是每簇细胞,第二列是给它注释的名字

image.png

不用按照0-11的顺序填,先看每一行又大又红的点,对应的细胞是什么。没有又大又红的只能退而求其次,看最大的点或者是最红的点。有时最大的和最红的不是同一个,就为难了,没有标准答案,我习惯是优先按照最红的,除非最红的太小了才按照最大的。。。比如1这一簇,最红的是CD8 T,最大的是NK,注释为CD8 T。

image.png
#给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
image.png

手动注释要检查是否合理,比如T细胞和NK细胞应该是在一起的,B细胞应该是单独的。

4.2 自动注释

最常用的自动注释R包是SingleR。它使用的参考数据是celldex里的

library(celldex)
library(SingleR)
ls("package:celldex")
image.png

每个参考数据的介绍要看它的帮助文档(例如?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
image.png

版权声明:
作者:Alex
链接:https://www.techfm.club/p/133971.html
来源:TechFM
文章版权归作者所有,未经允许请勿转载。

THE END
分享
二维码
< <上一篇
下一篇>>