人人都需掌握的单细胞分群聚类分析&Marker基因的可视化

rm(list=ls())
options(stringsAsFactors = F)
getwd()
library(Seurat)
library(ggplot2)
#install.packages(c('clustree'))
library(clustree)
library(cowplot)
library(dplyr)
library(patchwork)

library(data.table)
dat=fread( "GSE182434_raw_count_matrix.txt.gz",data.table = F)
dim(dat)
dat[1:4,1:4]
rownames(dat)=dat[,1]
dat=dat[,-1]
dat[1:4,1:4]
annotation = fread( "GSE182434_cell_annotation.txt.gz",data.table = F)
annotation[1:4,1:4]
table(annotation$Patient,annotation$CellType)

#筛选条件①、四个样本 ：DLBCL002 、DLBCL007、 DLBCL008、DLBCL111
#筛选条件②、CD8+细胞
sample = c("DLBCL002","DLBCL007","DLBCL008","DLBCL111")
CellType = "T cells CD8"
annotation = annotation[annotation$Patient %in% sample,]
annotation = annotation[annotation$CellType %in% CellType,]
dim(annotation)  #4个样本，一共5千个CD8细胞
dat = dat[,annotation$ID]
dat[1:5,1:5]
dim(dat)
View(dat)
View(annotation)

library(Seurat)
sce <- CreateSeuratObject(counts = dat, 
                          project = "sce", #项目名称
                          min.cells = 3, 
                          min.features = 200)
sce 
View(sce$nCount_RNA)
View(sce$nFeature_RNA)
View(sce$orig.ident)

sce$orig.ident<-annotation$Patient #修改sce 矩阵样本的分组信息
View(sce$orig.ident)

sce = FindVariableFeatures(sce,selection.method = "vst", nfeatures = 2000)#它们在某些细胞中高表达，而在其他细胞中低表达，在下游分析中关注这些基因有助于突出单细胞数据集中的生物信号。
#在 Seurat 中的过程，通过直接建模单细胞数据中固有的均值-方差关系来改进以前的版本，并在FindVariableFeatures()函数中实现。默认情况下，为每个数据集返回 2,000 个特征基因用于下游分析。
鉴别前10的高变基因。
top10 <- head(VariableFeatures(sce), 10)
画出没有标签的高变基因图
plot1 <- VariableFeaturePlot(sce)
加入前十个标签
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2

sce = NormalizeData(sce, normalization.method = "LogNormalize", scale.factor = 10000)
sce = FindVariableFeatures(sce)
sce = ScaleData(sce, vars.to.regress = c("nFeature_RNA", "percent_mito"))
sce = RunPCA(sce, npcs = 20)#PCA降维，对高变基因降维
print(sce [["pca"]], dims = 1:5, nfeatures = 5)#展示一部分结果
VizDimLoadings(sce , dims = 1:2, reduction = "pca")##点图形式展示
DimPlot(sce, reduction = "pca")#投影的降维图
DimHeatmap(sce, dims = 1:15, cells = 500, balanced = TRUE)#DimHeatmap()允许轻松探索数据集中异质性的主要来源，并且在尝试决定要包括哪些 PC 以进行进一步的下游分析时非常有用。单元格和特征都根据其 PCA 分数排序。

#Seurat 提供了多种非线性降维技术，例如 tSNE 和 UMAP，以可视化和探索这些数据集。这些算法的目标是学习数据的底层流形，以便将相似的单元格放在低维空间中。上面确定的基于图形的集群中的单元格应该在这些降维图上共同定位。作为 UMAP 和 tSNE 的输入，我们建议使用相同的 PC 作为聚类分析的输入。
sce = RunTSNE(sce, npcs = 20)
sce = Runtsne(sce, dims = 1:10)
sce = RunUMAP(sce, dims = 1:10)
sce = FindNeighbors(sce, dims = 1:20, k.param = 60, prune.SNN = 1/15)  #这个dims要和PCA的npcs对应起来。
DimPlot(sce,reduction = "umap",label=T,group.by = "orig.ident")
ggsave(filename="no_harmony_DimPlot.pdf")

library(harmony)
sce <- RunHarmony(sce, "orig.ident")
names(sce@reductions)
sce = RunUMAP(sce,  dims = 1:15, reduction = "harmony")
sce = RunTSNE(sce, npcs = 20,reduction = "harmony")
sce = FindNeighbors(sce, reduction = "harmony", dims = 1:15)
DimPlot(sce,reduction = "umap",label=T,group.by = "orig.ident")
ggsave(filename="harmony_DimPlot.pdf")

#分辨率分类
#文中是分了6群
#设置不同的分辨率，观察分群效果(选择哪一个？)
for (res in c(0.01, 0.05, 0.1,0.15, 0.2, 0.3, 0.4,0.45,0.5,0.8,0.9,1)) {
  sce=FindClusters(sce, resolution = res, algorithm = 1)
}#细胞分类
colnames(sce@meta.data)
apply(sce@meta.data[,grep("RNA_snn_res",colnames(sce@meta.data))],2,table)
#0.8

p1_dim = plot_grid(ncol = 4, DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.0.01") + 
                     ggtitle("louvain_0.01"), DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.0.05") + 
                     ggtitle("louvain_0.05"), DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.0.1") + 
                     ggtitle("louvain_0.1"), DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.0.2") + 
                     ggtitle("louvain_0.2"))
p1_dim

p1_dim = plot_grid(ncol = 5, DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.0.3") + 
                     ggtitle("louvain_0.3"),DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.0.5") + 
                     ggtitle("louvain_0.5"),DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.0.8") + 
                     ggtitle("louvain_0.8"), DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.0.9") + 
                     ggtitle("louvain_0.9"),DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.1") + 
                     ggtitle("louvain_1"))
p1_dim

p2_tree = clustree(sce@meta.data, prefix = "RNA_snn_res.")
p2_tree
ggsave(plot = p2_tree, filename="Tree_diff_resolution.pdf",width = 10,height = 11)

#最后分群选择的resolution=0.4，细胞聚类为7群。
sce1 <- FindClusters(sce, resolution = 0.4) 
p1 <- DimPlot(sce1, reduction = "umap", label = TRUE, repel = TRUE)
p1 

p2 <- DimPlot(sce1, reduction = "umap", group.by = "orig.ident")
p2

#也可以单独观察每个样本中细胞的聚类情况
p3=DimPlot(sce1, reduction = "umap", split.by = "orig.ident")
p3

#利用 FindMarkers 命令，可以找到找到各个细胞类型中与其他类别的差异表达基因，作为该细胞类型的生物学标记基因。其中ident.1参数设置待分析的细胞类别，min.pct表示该基因表达数目占该类细胞总数的比例
#install.packages("metap")
library(metap)
#1_find all markers of cluster 1
cluster1.markers <- FindMarkers(sce, ident.1 = 6, min.pct = 0.25)
head(cluster1.markers, n = 5)
#利用 DoHeatmap 命令可以可视化marker基因的表达
pbmc.markers <- FindAllMarkers(sce, only.pos = TRUE, min.pct = 0.25)
#2_以cluster 6为例寻找conserved marker
DefaultAssay(sce1) <- "RNA"
nk.markers <- FindConservedMarkers(sce1, ident.1 = 6, grouping.var = "orig.ident", verbose = FALSE)
nk.markers #对每个cluster执行上述操作，就能找出所有细胞类的conserved marker
View(nk.markers)
head(nk.markers, n = 9)
#对marker在不同细胞中的丰度进行可视化
p4=FeaturePlot(sce1, features = c("RPL41", "RPL30", "RPS15A", "RPS27", "RPS12", "RPS28", "RPLP1",
"RPLP10", "RPS26"), min.cutoff = "q9")
ggsave(plot = p4, filename="FeaturePlot_conserved_marker.pdf",width = 10,height = 11)

#Seurat提供了许多方法使我们能够方便的探索感兴趣的基因在各个细胞类型中的表达情况
VlnPlot(sce1, features = c("RPL41", "RPL30")) #这里，我们随便选了两个基因