单细胞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