回溯亚群细分效果到初始分群

分享是一种态度

• 1 加载R包
• 2 读取Seurat object
• 3 提取子集并做亚群分析
• 4 修改亚群命名，便于区分初始分群序号
• 5 回溯亚群细分效果到初始分群

1 加载R包

代码语言：javascript复制

library(Seurat)
library(ggplot2)

2 读取Seurat object

代码语言：javascript复制

# Remove variables
rm(list = ls())

# Set workdir
workdir <- 'F:/R_Language/R_Practice/scRNA_Seq_column'
setwd(workdir)

# Load seurat object variable
data.filt <- readRDS(file = "data/Raw_data/Cluster_all_again.Rdata")
print(table(Idents(data.filt)))

代码语言：javascript复制

## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16 
## 406 399 352 340 332 283 247 208 184 181 159  81  56  47  41  37

代码语言：javascript复制

DT::datatable(data.filt@meta.data)

3 提取子集并做亚群分析

代码语言：javascript复制

sub_cells <- subset(data.filt, cluster == 1)
DT::datatable(sub_cells@meta.data)

代码语言：javascript复制

# Determine the K-nearest neighbor graph
sub_cells <- FindNeighbors(sub_cells, reduction="pca", dims = 1:20, verbose = F)
# Determine the clusters for various resolutions
sub_cells <- FindClusters(sub_cells, resolution = 0.5, verbose = F)

# t-SNE and Clustering
sub_cells <- RunTSNE(sub_cells, reduction = "pca", dims = 1:20, verbose = F) # check_duplicates = FALSE, Remove duplicates before running TSNE
sub_cells <- RunUMAP(sub_cells, reduction = "pca", dims = 1:20, verbose = F)

4 修改亚群命名，便于区分初始分群序号

代码语言：javascript复制

## Rename sub cells
metadata <- sub_cells@meta.data
metadata$cluster <- Idents(sub_cells)
metadata$sub_cluster <- paste0('M', metadata$cluster)
sub_cells@meta.data <- metadata
table(sub_cells$sub_cluster)

代码语言：javascript复制

## 
##  M0  M1  M2 
## 170 151  85

5 回溯亚群细分效果到初始分群

代码语言：javascript复制

## Mapping to old clusters
metadata <- data.filt@meta.data
metadata$cell_type <- NA
metadata$cell_type <- sub_cells@meta.data[match(rownames(data.filt@meta.data), rownames(sub_cells@meta.data)), 'sub_cluster']
print(table(metadata$cell_type))

代码语言：javascript复制

## 
##  M0  M1  M2 
## 170 151  85

代码语言：javascript复制

DT::datatable(metadata)

代码语言：javascript复制

celltype_names <- NULL
for(i in 1:dim(metadata)[1]){
  # print(i)
  sub_data <- metadata[i,]
  if(is.na(sub_data$cell_type)){
    # print('Change value')
    sub_data$cell_type <- sub_data$cluster
    celltype_names <- c(celltype_names, sub_data$cell_type)
  }else{
    celltype_names <- c(celltype_names, sub_data$cell_type)
  }
}

print(table(celltype_names))

代码语言：javascript复制

## celltype_names
##  10  11  12  13  14  15  16   2   3   4   5   6   7   8   9  M0  M1  M2 
## 181 159  81  56  47  41  37 399 352 340 332 283 247 208 184 170 151  85

代码语言：javascript复制

metadata$cell_type <- celltype_names
metadata$cell_type <- factor(metadata$cell_type)
metadata$cell_type <- factor(metadata$cell_type, levels = c(names(table(sub_cells$sub_cluster)), 2:nlevels(Idents(data.filt))))
print(table(metadata$cell_type))

代码语言：javascript复制

## 
##  M0  M1  M2   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16 
## 170 151  85 399 352 340 332 283 247 208 184 181 159  81  56  47  41  37

代码语言：javascript复制

data.filt@meta.data <- metadata 

DimPlot(data.filt, reduction='tsne', label=T, label.size=5, pt.size=1, group.by='cluster')   ggtitle('original celltype')

代码语言：javascript复制

DimPlot(data.filt, reduction='tsne', label=T, label.size=5, pt.size=1, group.by='cell_type')   ggtitle('new celltype')

往期回顾

对细胞簇重命名

OSCA单细胞数据分析笔记10—Marker gene detection

NC单细胞文章复现(七)：Gene expression signatures(2)

多组学分析肺结核队列的记忆T细胞状态

---

---

---

如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程

• 生信爆款入门-2021第4期
• 数据挖掘（GEO,TCGA,单细胞）2021第4期
• 明码标价之共享96线程384G内存服务器

detection expression geo marker object

0 人点赞