本文涉及：亞群合并，亞群提取卿嘲，亞群細(xì)分，亞群抽樣
一般來講夫壁，不同亞群的細(xì)胞具有不同的功能拾枣，因此亞群分析對于研究十分重要，會(huì)影響后續(xù)分析結(jié)果盒让。
進(jìn)行細(xì)胞亞群的分群的依據(jù)：
- 細(xì)胞異質(zhì)性（每個(gè)細(xì)胞都有獨(dú)特的表達(dá)模式和功能梅肤，都有自己特有的基因）；
- 細(xì)胞共性（同一類型的細(xì)胞都有近似的表達(dá)模式）邑茄；
- 生物學(xué)基礎(chǔ)知識(shí)（基于已有的知識(shí)姨蝴，對細(xì)胞進(jìn)行分類鑒定）
單細(xì)胞聚類很少能一次性得到符合預(yù)期的結(jié)果，需要對結(jié)果不斷調(diào)整肺缕，如細(xì)胞亞群數(shù)目調(diào)整左医。如前文中所言，Seurat工具可通過調(diào)整 FindClusters 函數(shù)中的resolution參數(shù)進(jìn)行調(diào)整細(xì)胞亞群數(shù)目同木，一般在0.1-1之間炒辉，值越大，亞群數(shù)目越多泉手，但是亞群數(shù)目過多黔寇，后續(xù)分析越耗力耗神。另一種方法是對感興趣的細(xì)胞亞群進(jìn)行亞群細(xì)分或者對相近的細(xì)胞亞群進(jìn)行合并斩萌，進(jìn)而來增加或減少亞群數(shù)目缝裤。
然而在單細(xì)胞數(shù)據(jù)分析中，一般初步的亞群分類結(jié)果可能將同類細(xì)胞分成若亞類颊郎”锓桑基于分子標(biāo)記，我們可以將其歸為同類細(xì)胞的亞類重新合并姆吭。
下面以 seurat 官方教程為例：
## 魔幻清空
rm(list = ls())
## 加載R包
library(Seurat)
library(ggplot2)
library(patchwork)
library(dplyr)
load(file = 'basic.sce.pbmc.Rdata')
levels(Idents(pbmc)) #查看細(xì)胞亞群
# [1] "Naive CD4 T" "CD14+ Mono" "Memory CD4 T" "B" "CD8 T" "FCGR3A+ Mono"
# [7] "NK" "DC" "Platelet"
## UMAP可視化
DimPlot(pbmc, reduction = 'umap',
label = TRUE, pt.size = 0.5) + NoLegend()
1. 亞群合并
假如我們的生物學(xué)背景知識(shí)不夠榛做，不需要把T細(xì)胞分成 "Naive CD4 T" , "Memory CD4 T" , "CD8 T", "NK" 這些亞群，因此可以將這些亞群合并為T細(xì)胞這個(gè)大的亞群：
## 'PTPRC', 'CD3D', 'CD3E','CD4', 'CD8A','FOXP3','KLRD1', ## Tcells
## 'CD19', 'CD79A', 'MS4A1' , # B cells
## 'IGHG1', 'MZB1', 'SDC1', # plasma
## 'CD68', 'CD163', 'CD14', 'CD86','CCL22', 'S100A4', ## DC(belong to monocyte)
## 'CD68', 'CD163', 'MRC1', 'MSR1', 'S100A8', ## Macrophage (belong to monocyte)
## 'PRF1', 'NKG7', 'KLRD1', ## NKcells
## 'PPBP', ##Platelet
sce=pbmc ##為防止影響pbmc本來的數(shù)據(jù)内狸，接下來將以sce變量進(jìn)行
genes_to_check = c('PTPRC', 'CD3D', 'CD3E', 'PRF1' , 'NKG7',
'CD19', 'CD79A', 'MS4A1' ,
'CD68', 'CD163', 'CD14',
'FCER1A', 'PPBP' )
DotPlot(sce, group.by = 'seurat_clusters',
features = unique(genes_to_check)) + RotatedAxis()
接下來我們對亞群進(jìn)行合并检眯，目的只能區(qū)分 B、DC昆淡、Mono锰瘸、Platelet、T 這5個(gè)細(xì)胞亞群昂灵。
方法一：使用 RenameIdents 函數(shù)
levels(Idents(sce)) #查看細(xì)胞亞群避凝，與上述結(jié)果一致
***根據(jù)levels(Idents(sce)) 順序重新賦予對應(yīng)的 B舞萄、DC、Mono管削、Platelet倒脓、T
這5個(gè)細(xì)胞亞群名稱，順序不能亂
new.cluster.ids <- c("T", "Mono", "T",
"B", "T", "Mono",
"T", "DC", "Platelet")
names(new.cluster.ids) <- levels(sce)
sce <- RenameIdents(sce, new.cluster.ids)
levels(sce) #查看是否已改名
# [1] "T" "Mono" "B" "DC" "Platelet"
DimPlot(sce, reduction = 'umap',
label = TRUE, pt.size = 0.5) + NoLegend()
方法二：使用unname函數(shù)配合向量：
cluster2celltype <- c("0"="T",
"1"="Mono",
"2"="T",
"3"= "B",
"4"= "T",
"5"= "Mono",
"6"= "T",
"7"= "DC",
"8"= "Platelet")
sce[['cell_type']] = unname(cluster2celltype[sce@meta.data$seurat_clusters])
DimPlot(sce, reduction = 'umap', group.by = 'cell_type',
label = TRUE, pt.size = 0.5) + NoLegend()
方法三：使用數(shù)據(jù)框
(n=length(unique(sce@meta.data$seurat_clusters))) #獲取亞群個(gè)數(shù)
celltype=data.frame(ClusterID=0:(n-1),
celltype='unkown') #構(gòu)建數(shù)據(jù)框
## 判斷亞群ID屬于那類細(xì)胞
celltype[celltype$ClusterID %in% c(0,2,4,6),2]='T'
celltype[celltype$ClusterID %in% c(3),2]='B'
celltype[celltype$ClusterID %in% c(1,5),2]='Mono'
celltype[celltype$ClusterID %in% c(7),2]='DC'
celltype[celltype$ClusterID %in% c(8),2]='Platelet'
## 重新賦值
sce@meta.data$celltype = "NA"
for(i in 1:nrow(celltype)){
sce@meta.data[which(sce@meta.data$seurat_clusters == celltype$ClusterID[i]),'celltype'] <- celltype$celltype[i]}
table(sce@meta.data$celltype)
# B DC Mono Platelet T
# 342 32 642 14 1608
DimPlot(sce, reduction = 'umap', group.by = 'celltype',
label = TRUE, pt.size = 0.5) + NoLegend()
## 查看多種方法結(jié)果
table(sce@meta.data$cell_type,
sce@meta.data$celltype)
# B DC Mono Platelet T
# B 342 0 0 0 0
# DC 0 32 0 0 0
# Mono 0 0 642 0 0
# Platelet 0 0 0 14 0
# T 0 0 0 0 1608
## 可視化一下
p1=DimPlot(sce, reduction = 'umap', group.by = 'celltype',
label = TRUE, pt.size = 0.5) + NoLegend()
p2=DotPlot(sce, group.by = 'celltype',
features = unique(genes_to_check)) + RotatedAxis()
p1+p2
最后總結(jié)一下含思，合并亞群其實(shí)就是亞群重命名崎弃，實(shí)現(xiàn)的基本原理就是將分析過程聚類亞群0-8重新命名，核心技術(shù)是R語言中匹配替換功能茸俭。
當(dāng)然這僅僅是示例吊履，實(shí)際分析過程中可能不會(huì)這么簡單安皱，會(huì)略微復(fù)雜调鬓，但是萬變不離其宗，掌握好R語言基本功是硬道理酌伊。
2. 亞群提取
在單細(xì)胞分析過程中腾窝，我們通常會(huì)對感興趣的細(xì)胞亞群進(jìn)行亞群細(xì)分，這樣可以把一個(gè)亞群或者多個(gè)亞群提取出來居砖，然后再進(jìn)行亞群細(xì)分虹脯。
亞群細(xì)分有兩種方法：
第一種，調(diào)整FindClusters函數(shù)中的resolution參數(shù)使亞群數(shù)目增多奏候；
第二種循集，將此亞群提取出來，再重新降維聚類分群蔗草。
前面咒彤，我們將幾個(gè)亞群合并為T細(xì)胞這個(gè)大亞群。接下來我們看看如何對這部分細(xì)胞亞群進(jìn)行亞群細(xì)分咒精。
## 重新讀取數(shù)據(jù)
rm(list = ls())
library(Seurat)
library(ggplot2)
library(patchwork)
library(dplyr)
load(file = 'basic.sce.pbmc.Rdata')
sce=pbmc
首先定位到感興趣的亞群
levels(sce)
# [1] "Naive CD4 T" "CD14+ Mono" "Memory CD4 T" "B" "CD8 T" "FCGR3A+ Mono"
# [7] "NK" "DC" "Platelet"
genes_to_check = c('PTPRC', 'CD3D', 'CD3E',
'CD4','IL7R','NKG7','CD8A')
p1=DimPlot(sce, reduction = 'umap', group.by = 'seurat_clusters',
label = TRUE, pt.size = 0.5) + NoLegend()
p2=DotPlot(sce, group.by = 'seurat_clusters',
features = unique(genes_to_check)) + RotatedAxis()
p1+p2
假設(shè)這個(gè)時(shí)候镶柱，我們想提取CD4的T細(xì)胞，那么根據(jù)上文聚類0/2/4/6均為T細(xì)胞模叙，其中0和2表達(dá)CD4相對4/6較高歇拆，但是其實(shí)示例里面的CD4的T細(xì)胞并不怎么高表達(dá)CD4，在此不深究范咨，繼續(xù)向下走故觅。
提取指定單細(xì)胞亞群
cd4_sce1 = sce[,sce@meta.data$seurat_clusters %in% c(0,2)]
cd4_sce2 = sce[, Idents(sce) %in% c( "Naive CD4 T" , "Memory CD4 T" )]
cd4_sce3 = subset(sce,seurat_clusters %in% c(0,2))
## 較簡單，核心原理就是R里取子集的3種策略：邏輯值渠啊，坐標(biāo)逻卖，名字
重新降維聚類分群
sce=cd4_sce1
sce <- NormalizeData(sce, normalization.method = "LogNormalize", scale.factor = 1e4)
sce <- FindVariableFeatures(sce, selection.method = 'vst', nfeatures = 2000)
sce <- ScaleData(sce, vars.to.regress = "percent.mt")
sce <- RunPCA(sce, features = VariableFeatures(object = sce))
sce <- FindNeighbors(sce, dims = 1:10)
sce <- FindClusters(sce, resolution = 1 )
head(Idents(sce), 5)
table(sce$seurat_clusters)
sce <- RunUMAP(sce, dims = 1:10)
DimPlot(sce, reduction = 'umap')
genes_to_check = c('PTPRC', 'CD3D', 'CD3E', 'FOXP3',
'CD4','IL7R','NKG7','CD8A')
DotPlot(sce, group.by = 'seurat_clusters',
features = unique(genes_to_check)) + RotatedAxis()
# 亞群水平
p1=DimPlot(sce, reduction = 'umap', group.by = 'seurat_clusters',
label = TRUE, pt.size = 0.5) + NoLegend()
p2=DotPlot(sce, group.by = 'seurat_clusters',
features = unique(genes_to_check)) + RotatedAxis()
p1+p2
可以看到，這次重新降維聚類分群已經(jīng)是非常的勉強(qiáng)昭抒， 各細(xì)胞亞群之間的界限并不是很清晰评也。僅僅提取出來CD4的T細(xì)胞進(jìn)行細(xì)分亞群炼杖，而結(jié)果又多出來了一個(gè)CD8 T細(xì)胞亞群，令人費(fèi)解盗迟。嘗試將resolution改成0.5再試一次坤邪，結(jié)果如下：
這次分群比上次要好點(diǎn)，但還是難以接受罚缕，掌握技能而已艇纺，不糾結(jié)，后面再深入探查邮弹。
單細(xì)胞分多少群合適
我的細(xì)胞到底分多少個(gè)群是合適的黔衡？每個(gè)細(xì)胞都是獨(dú)一無二的，也就是最小可以以單細(xì)胞為單位腌乡，但是高通量的單細(xì)胞測序技術(shù)是數(shù)以萬計(jì)的細(xì)胞盟劫。那么，我的細(xì)胞到底分多少個(gè)群是合適的与纽？
那我們來試試clustree侣签，首先依舊是讀取我們的數(shù)據(jù)
## 重新讀取數(shù)據(jù)
rm(list = ls())
library(Seurat)
library(ggplot2)
library(patchwork)
library(dplyr)
load(file = 'basic.sce.pbmc.Rdata')
sce=pbmc
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5) #直接使用了官方resolution = 0.5
## 實(shí)際上resolution 是可以多次調(diào)試的,執(zhí)行不同resolution 下的分群
library(clustree)
sce <- FindClusters(
object = sce,
resolution = c(seq(.1,1.6,.2))
)
clustree(sce@meta.data, prefix = "RNA_snn_res.")
如圖所示，可以看到不同的resolution 急迂，分群的變化情況影所。紅框圈出來的對應(yīng)我們使用的 resolution = 0.5 ，聚類9個(gè)亞群僚碎。根據(jù)動(dòng)態(tài)分群的樹猴娩，可見3對應(yīng)的B細(xì)胞亞群，無論怎么樣調(diào)整參數(shù)勺阐，都很難細(xì)分了卷中，同樣還有7對應(yīng)DC亞群，和8對應(yīng)的Platelet亞群也是很難再細(xì)分啦皆看。但是T細(xì)胞和monocyte還有進(jìn)一步細(xì)分的可能性仓坞！
p1=DimPlot(sce, reduction = 'umap', group.by = 'RNA_snn_res.0.5',
label = TRUE, pt.size = 0.5) + NoLegend()
p2=DimPlot(sce, reduction = 'umap',group.by = 'RNA_snn_res.1.5',
label = TRUE, pt.size = 0.5) + NoLegend()
p1+p2
這個(gè)時(shí)候就需要根據(jù)研究目的判斷一下了，即使有分群的可能性腰吟，但不代表一定要進(jìn)行細(xì)分亞群无埃，一味的追求細(xì)分亞群是沒有意義的！
比如前面的CD4的T細(xì)胞亞群細(xì)分：
load(file = 'sce.cd4.subset.Rdata')
#先執(zhí)行不同resolution 下的分群
library(Seurat)
library(clustree)
sce <- FindClusters(
object = sce,
resolution = c(seq(.1,1.6,.2))
)
clustree(sce@meta.data, prefix = "RNA_snn_res.")
這時(shí)候分群混亂毛雇，相互交錯(cuò)嫉称，個(gè)人覺得沒啥意義。分群是為了下一步分析灵疮，切莫作繭自縛织阅。
3. 單細(xì)胞亞群抽樣
單細(xì)胞的gsva, 細(xì)胞通訊，轉(zhuǎn)錄因子震捣，擬時(shí)序, inferCNV這些分析荔棉，特別的消耗計(jì)算資源闹炉，因?yàn)閷?shí)際項(xiàng)目中每個(gè)細(xì)胞亞群都是過萬的細(xì)胞，希望可以將這些單細(xì)胞亞群進(jìn)行抽樣润樱，使得其細(xì)胞數(shù)量一致渣触。
rm(list = ls())
library(Seurat)
library(ggplot2)
library(patchwork)
library(dplyr)
## 讀取數(shù)據(jù)
load(file = 'basic.sce.pbmc.Rdata')
DimPlot(pbmc, reduction = 'umap',
label = TRUE, pt.size = 0.5) + NoLegend()
sce=pbmc
features= c('IL7R', 'CCR7','CD14', 'LYZ', 'IL7R', 'S100A4',"MS4A1", "CD8A",'FOXP3',
'FCGR3A', 'MS4A7', 'GNLY', 'NKG7',
'FCER1A', 'CST3','PPBP')
DoHeatmap(subset(sce ),
features = features,
size = 3
)
如圖，細(xì)胞數(shù)量較少的亞群在熱圖中占很少壹若。
方法一：subset函數(shù)
table(Idents(sce))
# Naive CD4 T CD14+ Mono Memory CD4 T B CD8 T FCGR3A+ Mono NK
# 709 480 429 342 316 162 154
# DC Platelet
# 32 14
## 可見最少的一個(gè)亞群Platelet細(xì)胞只有14個(gè)細(xì)胞嗅钻，因此我們對每個(gè)亞群抽取15個(gè)細(xì)胞繪制熱圖
DoHeatmap(subset(sce, downsample = 15),
features = features, size = 3)
方法二：自定義函數(shù)
# 每個(gè)細(xì)胞亞群抽10
allCells=names(Idents(sce))
allType = levels(Idents(sce))
choose_Cells = unlist(lapply(allType, function(x){
cgCells = allCells[Idents(sce)== x ]
cg=sample(cgCells,10)
cg
}))
cg_sce = sce[, allCells %in% choose_Cells]
cg_sce
as.data.frame(table(Idents(cg_sce)))
# Var1 Freq
# 1 Naive CD4 T 10
# 2 CD14+ Mono 10
# 3 Memory CD4 T 10
# 4 B 10
# 5 CD8 T 10
# 6 FCGR3A+ Mono 10
# 7 NK 10
# 8 DC 10
# 9 Platelet 10
轉(zhuǎn)自：https://mp.weixin.qq.com/s/mF5rziWHlCIDCIBJloSTqA
可配合此視頻學(xué)習(xí)：https://www.bilibili.com/video/BV19Q4y1R7cu?p=1
