单细胞数据挖掘(5,6)-聚类、筛选marker基因及拼图(生信

本笔记来源于B站@生信技能树-jimmy；学习视频链接: 「生信技能树」单细胞数据挖掘

5. 聚类，筛选marker基因，可视化

# 5.1 聚类
# 基于上一步骤的结果
pc.num=1:20
# 基于PCA数据
scRNA <- FindNeighbors(scRNA, dims = pc.num) 
# dims参数，需要指定哪些pc轴用于分析；这里利用上面的分析，选择20
scRNA <- FindClusters(scRNA, resolution = 0.5) # resolution用于分群大小设定，在0-1之间
table(scRNA@meta.data$seurat_clusters)  #查看每个cluster有多少基因

scRNA = RunTSNE(scRNA, dims = pc.num)
DimPlot(scRNA, reduction = "tsne",label=T)
?RunTSNE
p3_1 <- DimPlot(scRNA, reduction = "tsne",label=T) +
  labs(tag = "E")
p3_1

SingleCell5_1.jpeg

# 5.2 marker gene
# 进行差异分析，一般使用标准化数据！！
scRNA <- NormalizeData(scRNA, normalization.method = "LogNormalize")
# 结果储存在"data" 这个slot 里
GetAssayData(scRNA,slot="data",assay="RNA")[1:8,1:4]

diff.wilcox = FindAllMarkers(scRNA)  ## 默认使用wilcox方法挑选差异基因; 耗时大概4-5min
# if test.use is "negbinom", "poisson", or "DESeq2", slot will be set to "counts"!!!
save(diff.wilcox, file = "diff.wilcox.Rdata") 
head(diff.wilcox)
dim(diff.wilcox) # 27002     7

SingleCell5_2.JPG

library(tidyverse)
all.markers = diff.wilcox %>% select(gene, everything()) %>%
  subset(p_val<0.05 & abs(diff.wilcox$avg_logFC) > 0.5)
# select函数将gene这列排到了前面
# An adjusted P value < 0.05and | log 2 [fold change (FC)] | > 0.5 
# were considered the 2 cutoff criteria for identifying marker genes.

dim(all.markers)
summary(all.markers)
save(all.markers,file = "../../markergene.Rdata") # 用于后续的轨迹分析

top10 = all.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
top10
top10 = CaseMatch(search = as.vector(top10$gene), match = rownames(scRNA)) 
top10
length(top10)  #130
length(unique(sort(top10))) # 125

SingleCell5_3.JPG

p3_2 <- DoHeatmap(scRNA, features = top10, group.by = "seurat_clusters")
p3_2
p3_1 | p3_2  # Fig1下图

Fig1.F,G.jpg

6. 拼图，比较

p <- (p1_1 | p1_2 | p1_3 ) /
  ((p2_1| p2_2 | p2_3) /
     (p3_1 | p3_2))
ggsave("../../Fig1.pdf", plot = p, width = 15, height = 18) 

rm(list=ls())
save(scRNA,file = "scRNA.Rdata")

Fig1_all.jpg