整合scRNA数据后,那么就是聚类处理,选择合适的聚类参数是识别celltype的重点。更多知识分享请到 https://zouhua.top/。
加载R包
library(dplyr)
library(Seurat)
library(ggplot2)
# 设置参数
rm(list = ls())
options(stringsAsFactors = F)
options(future.globals.maxSize = 4000 * 1024^2)
colors <- c("#A6CEE3", "#1F78B4", "#08306B", "#B2DF8A", "#006D2C", "#8E0152",
"#DE77AE", "#CAB2D6", "#6A3D9A", "#FB9A99", "#E31A1C", "#B15928",
"#619CFF","#FF67A4","#00BCD8", "#EE2B2B", "#2D6BB4")
读取数据
data.integrated <- readRDS("../../Result/RDS/data.n8.sct.rds")
聚类前可视化
# 若原数据没有PCA结果,则运行下面部分
data.integrated <- RunPCA(object = data.integrated, verbose = FALSE)
data.integrated <- RunTSNE(data.integrated, reduction = "pca", dims = 1:30)
data.integrated <- RunUMAP(data.integrated, reduction = "pca", dims = 1:30)
cowplot::plot_grid(
DimPlot(data.integrated,
reduction = "umap",
label = TRUE,
label.size = 6),
DimPlot(data.integrated,
reduction = "umap",
label = TRUE,
label.size = 6,
split.by = "Batch")+NoLegend(),
ncol = 2)
查看PC的metagenes(热图展示)
DimHeatmap(data.integrated,
reduction = "pca",
dims = 1:9,
cells = 500,
balanced = TRUE)
print(x = data.integrated[["pca"]],
dims = 1:10,
nfeatures = 5)
ElbowPlot(object = data.integrated,
ndims = 40)+
geom_vline(xintercept = 30, linetype=2)+
geom_hline(yintercept = 2, linetype=2)
聚类,设置不同的分辨率
# Determine the K-nearest neighbor graph
data.integrated <- FindNeighbors(object = data.integrated,
dims = 1:30)
# Determine the clusters for various resolutions
data.integrated <- FindClusters(object = data.integrated,
resolution = c(0.4, 0.6, 0.8, 1.0, 1.4))
查看不同分辨率下簇的分布情况
-
簇之间明显区分;
-
簇没有跨区域分散;
Idents(object = data.integrated) <- "integrated_snn_res.0.4"
DimPlot(data.integrated,
reduction = "umap",
label = TRUE,
label.size = 6)
Idents(object = data.integrated) <- "integrated_snn_res.0.6"
DimPlot(data.integrated,
reduction = "umap",
label = TRUE,
label.size = 8)
Idents(object = data.integrated) <- "integrated_snn_res.1.4"
DimPlot(data.integrated,
reduction = "umap",
label = TRUE,
label.size = 6)
注意:如没有integration步骤,则使用RNA_snn_res.0.6。
根据聚类结果选择resolution参数
Idents(object = data.integrated) <- "integrated_snn_res.0.4"
saveRDS(data.integrated, "../../Result/RDS/data.all.cluster04.rds", compress = TRUE)
Reference
参考文章如引起任何侵权问题,可以与我联系,谢谢。
网友评论