#-----------单细胞亚群合并与提取----------
#>>进行细胞亚群的分群的依据：#>细胞异质性（每个细胞都有独特的表达模式和功能，都有自己特有的基因）；#>细胞共性（同一类型的细胞都有近似的表达模式）；#>生物学基础知识（基于已有的知识，对细胞进行分类鉴定）
#https://cloud.tencent.com/developer/article/1842672
## 魔幻清空rm(list = ls())## 加载R包library(Seurat)library(ggplot2)library(patchwork)library(dplyr)
getwd()setwd("C:\\Users\\Desktop\\Seurat_example")
## 导入上节保存的数据load(file = \\\\'basic.sce.pbmc.Rdata\\\\')levels(Idents(pbmc)) #查看细胞亚群 # [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()
#亚群合并#假如我们的生物学背景知识不够，不需要把T细胞分成 "Naive CD4 T" , "Memory CD4 T" , "CD8 T", "NK" 这些亚群，#因此可以将这些亚群合并为T细胞这个大的亚群：
## \\\\'PTPRC\\\\', \\\\'CD3D\\\\', \\\\'CD3E\\\\',\\\\'CD4\\\\', \\\\'CD8A\\\\',\\\\'FOXP3\\\\',\\\\'KLRD1\\\\', ## T cells## \\\\'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\\\\', ## NK cells## \\\\'PPBP\\\\', ##Platelet
sce=pbmc ##为防止影响pbmc本来的数据，接下来将以sce变量进行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()
sce$CellType DotPlot(sce, group.by = \\\\'CellType\\\\', features = unique(genes_to_check)) + RotatedAxis()
##>接下来我们对亚群进行合并，目的只能区分 B、DC、Mono、Platelet、T 这5个细胞亚群。#方法一：使用 RenameIdents 函数levels(Idents(sce)) #查看细胞亚群，与上述结果一致## 根据levels(Idents(sce)) 顺序重新赋予对应的 B、DC、Mono、Platelet、T 这5个细胞亚群名称，顺序不能乱new.cluster.ids "B", "T", "Mono", "T", "DC", "Platelet")names(new.cluster.ids) sce levels(sce) #查看是否已改名# [1] "T" "Mono" "B" "DC" "Platelet"
DimPlot(sce, reduction = \\\\'umap\\\\', label = TRUE, pt.size = 0.5) + NoLegend()
#方法二：使用unname函数配合向量：
cluster2celltype "1"="Mono", "2"="T", "3"= "B", "4"= "T", "5"= "Mono", "6"= "T", "7"= "DC", "8"= "Platelet")sce[[\\\\'cell_type\\\\']] = unname(cluster2celltype[[email protected]$seurat_clusters])DimPlot(sce, reduction = \\\\'umap\\\\', group.by = \\\\'cell_type\\\\', label = TRUE, pt.size = 0.5) + NoLegend()
#>注意使用 @ and $符号：[email protected]$seurat_clusters
sce$seurat_clusters[email protected]
#方法三：使用数据框(n=length(unique([email protected]$seurat_clusters))) #获取亚群个数celltype=data.frame(ClusterID=0:(n-1), celltype=\\\\'unkown\\\\') #构建数据框class(celltype)
## 判断亚群ID属于那类细胞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\\\\'
#>>>----------%in%----------------#>%in% is value matching and "returns a vector of the positions of (first) matches #>of its first argument in its second"#1:10 %in% 3:7#[1] FALSE FALSE TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE#is same as #sapply(1:10, function(a) any(a == 3:7))#[1] FALSE FALSE TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE
## 重新赋值[email protected]$celltype = "NA"for(i in 1:nrow(celltype)){ [email protected][which([email protected]$seurat_clusters == celltype$ClusterID[i]),\\\\'celltype\\\\'] table([email protected]$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()
## 查看多种方法结果table([email protected]$cell_type, [email protected]$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+p2p1+p2+p1+p2#p1+p2+p1+p2+p1+p2+p1+p2
#最后总结一下，合并亚群其实就是亚群重命名，#实现的基本原理就是将分析过程聚类亚群0-8重新命名，核心技术是R语言中匹配替换功能。
###----------------------------#####--------亚群提取------------####
#在单细胞分析过程中，我们通常会对感兴趣的细胞亚群进行亚群细分，#这样可以把一个亚群或者多个亚群提取出来，然后再进行亚群细分。#亚群细分有两种方法：#@>第一种，调整FindClusters函数中的resolution参数使亚群数目增多；#@>第二种，将此亚群提取出来，再重新降维聚类分群。
#上文我们将几个亚群合并为T细胞这个大亚群。#接下来我们看看如何对这部分细胞亚群进行亚群细分。
## 重新读取数据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
#提取指定单细胞亚群cd4_sce1 = sce[,[email protected]$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种策略：逻辑值，坐标，名字
#重新降维聚类分群sce=cd4_sce1sce sce sce sce sce sce head(Idents(sce), 5)table(sce$seurat_clusters) sce 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
#可以看到，这次重新降维聚类分群已经是非常的勉强， 各细胞亚群之间的界限并不是很清晰。#仅仅提取出来CD4的T细胞进行细分亚群，而结果又多出来了一个CD8 T细胞亚群，令人费解。#尝试将resolution改成0.5再试一次，结果如下：#这次分群比上次要好点，但还是难以接受，掌握技能而已，不纠结，后面再深入探查。
##----------------######-----########--------单细胞分多少群合适---###
#过滤双细胞的方法有很多种，一种比较直接粗暴的方法就是把细胞基因表达数目超过一定阈值的细胞去掉（比如PBMC细胞，阈值为2500），不过不同的样品阈值不同；另外一种方式就是通过算法来去掉双细胞，现在去掉双细胞的工具有很多种，比如DoubletFinder(https://www.cell.com/cell-systems/fulltext/S2405-4712(19)30073-0)、scrublet(https://github.com/AllonKleinLab/scrublet)、DoubletDecon(https://github.com/EDePasquale/DoubletDecon)、DoubletDetection(https://github.com/JonathanShor/DoubletDetection.git) 等等，python和R语言的工具都有，可以根据自己需要进行工具选择。
#那我们来试试clustree，首先依旧是读取我们的数据## 重新读取数据rm(list = ls())library(Seurat)library(ggplot2)library(patchwork)library(dplyr)load(file = \\\\'basic.sce.pbmc.Rdata\\\\')sce=pbmcpbmc pbmc
## 实际上resolution 是可以多次调试的,执行不同resolution 下的分群library(clustree)sce object = sce, resolution = c(seq(.1,1.6,.2)))clustree([email protected], prefix = "RNA_snn_res.") #这个很重要#如图所示，可以看到不同的resolution ，分群的变化情况。#红框圈出来的对应我们使用的 resolution = 0.5 ，聚类9个亚群。#根据动态分群的树，可见3对应的B细胞亚群，无论怎么样调整参数，都很难细分了，
