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 = scRNA@active.ident)
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
