单细胞转录组基础分析五：细胞再聚类

单细胞测序技术是近年最大的生命科学突破之一，相关文章频繁发表于各大顶级期刊，然而单细胞数据的分析依然是大家普遍面临的障碍。本专题将针对10X Genomics单细胞转录组数据演示各种主流分析，包括基于Seurat的基础分析、以及基于clusterProfiler、Monocle、SingleR等R包的延伸分析。不足之处请大家批评指正，欢迎添加Kinesin微信交流探讨！

单细胞数据分析中，一般需要对可以细分的细胞再聚类，比如本次分析中的T细胞群体可以细分为Navie T cells、CD8+ T cells、Treg cells、Tmemory cells等。细胞子集的提取使用subset函数，主要依据meta.data的分类项目选择，也可以自定义细胞barcode提取。

提取细胞子集

library(Seurat)
library(monocle)
library(tidyverse)
library(patchwork)
rm(list=ls())
dir.create("subcluster")
scRNA <- readRDS("scRNA.rds")
##提取细胞子集
Cells.sub <- subset(scRNA@meta.data, celltype=="T_cells")
scRNAsub <- subset(scRNA, cells=row.names(Cells.sub))

提重新降维聚类

因为再聚类的细胞之间差异比较小，所以聚类函数FindClusters()控制分辨率的参数建议调高到resolution = 0.9。

##PCA降维
scRNAsub <- FindVariableFeatures(scRNAsub, selection.method = "vst", nfeatures = 2000)
scale.genes <-  rownames(scRNAsub)
scRNAsub <- ScaleData(scRNAsub, features = scale.genes)
scRNAsub <- RunPCA(scRNAsub, features = VariableFeatures(scRNAsub))
ElbowPlot(scRNAsub, ndims=20, reduction="pca")
pc.num=1:10
##细胞聚类
scRNAsub <- FindNeighbors(scRNAsub, dims = pc.num) 
scRNAsub <- FindClusters(scRNAsub, resolution = 0.9)
table(scRNAsub@meta.data$seurat_clusters)
metadata <- scRNAsub@meta.data
cell_cluster <- data.frame(cell_ID=rownames(metadata), cluster_ID=metadata$seurat_clusters)
write.csv(cell_cluster,'subcluster/cell_cluster.csv',row.names = F)
##非线性降维
#tSNE
scRNAsub = RunTSNE(scRNAsub, dims = pc.num)
embed_tsne <- Embeddings(scRNAsub, 'tsne')
write.csv(embed_tsne,'subcluster/embed_tsne.csv')
plot1 = DimPlot(scRNAsub, reduction = "tsne") 
ggsave("subcluster/tSNE.pdf", plot = plot1, width = 8, height = 7)
ggsave("subcluster/tSNE.png", plot = plot1, width = 8, height = 7)
#UMAP
scRNAsub <- RunUMAP(scRNAsub, dims = pc.num)
embed_umap <- Embeddings(scRNAsub, 'umap')
write.csv(embed_umap,'subcluster/embed_umap.csv') 
plot2 = DimPlot(scRNAsub, reduction = "umap") 
ggsave("subcluster/UMAP.pdf", plot = plot2, width = 8, height = 7)
ggsave("subcluster/UMAP.png", plot = plot2, width = 8, height = 7)
#合并tSNE与UMAP
plotc <- plot1+plot2+ plot_layout(guides = 'collect')
ggsave("subcluster/tSNE_UMAP.pdf", plot = plotc, width = 10, height = 5)
ggsave("subcluster/tSNE_UMAP.png", plot = plotc, width = 10, height = 5)

Cluster差异分析

diff.wilcox = FindAllMarkers(scRNAsub)
all.markers = diff.wilcox %>% select(gene, everything()) %>% subset(p_val<0.05)
top10 = all.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
write.csv(all.markers, "subcluster/diff_genes_wilcox.csv", row.names = F)
write.csv(top10, "subcluster/top10_diff_genes_wilcox.csv", row.names = F)

SingleR细胞鉴定

Subcluster的细胞同样可以使用SingleR鉴定细胞类型。使用的时候注意调整参考数据库和分类标签，以便鉴定结果更有针对性。上节使用SingleR时使用的参考数据库是人类主要细胞图谱（HumanPrimaryCellAtlasData），采用分类标签是主分类标签（label.main）；这次建议使用人类免疫细胞数据（MonacoImmuneData），分类标签采用精细分类标签（label.fine）。希望详细了解SingleR的朋友可以到github看看：

https://github.com/dviraran/singler

##细胞类型鉴定
library(SingleR)
refdata <- MonacoImmuneData()
testdata <- GetAssayData(scRNAsub, slot="data")
clusters <- scRNAsub@meta.data$seurat_clusters
cellpred <- SingleR(test = testdata, ref = refdata, labels = refdata$label.fine, 
                     method = "cluster", clusters = clusters, 
                     assay.type.test = "logcounts", assay.type.ref = "logcounts")
celltype = data.frame(ClusterID=rownames(cellpred), celltype=cellpred$labels, stringsAsFactors = F)
write.csv(celltype,"subcluster/celltype_singleR.csv",row.names = F)
scRNAsub@meta.data$celltype = "NA"
for(i in 1:nrow(celltype)){
scRNAsub@meta.data[which(scRNAsub@meta.data$seurat_clusters == celltype$ClusterID[i]),'celltype'] <- celltype$celltype[i]}
p1 = DimPlot(scRNAsub, group.by="celltype", label=T, label.size=5, reduction='tsne')
p2 = DimPlot(scRNAsub, group.by="celltype", label=T, label.size=5, reduction='umap')
p3 = plotc <- p1+p2+ plot_layout(guides = 'collect')
ggsave("subcluster/tSNE_celltype.pdf", p1, width=7 ,height=6)
ggsave("subcluster/UMAP_celltype.pdf", p2, width=7 ,height=6)
ggsave("subcluster/celltype.pdf", p3, width=10 ,height=5)
ggsave("subcluster/celltype.png", p3, width=10 ,height=5)

获取帮助

本教程的目的在于把常用的单细胞分析流程串起来，适合有一定R语言基础的朋友参考。分析原理和代码我没有详细解释，官网有详细的教程和权威的解释，翻译好的中文教程也有多个版本，有兴趣的可以搜索一下。如果您学习本教程有一定困难，可以点击 “阅读原文” 找到作者联系方式，获取帮助。

往期回顾

把tcga大计划的CNS级别文章标题画一个词云

单细胞基因组学前沿会议推荐（*冷泉港亚洲学术会议）

并不是所有的批次效应都可以被矫正