如何对应两个不同单细胞数据集的分群结果？

生信技能树

发布于 2025-02-05 05:18:26

52400

代码可运行

文章被收录于专栏：生信技能树生信技能树

运行总次数：0

代码可运行

我们生信技能树有个学徒在过来中山进行学习的时候，学到了单细胞部分，然后他做了两个同样组织样本的数据，问：我这两个不同的数据集中，怎么样比较A数据集中的比如上皮细胞亚群与B数据集中的上皮细胞亚群是不是同一种上皮细胞亚群呢？这个问题挺有意思的，而且是一些多数据集整合算法或者利用已知注释好的单细胞数据来注释新的单细胞未知细胞亚群数据的基础思想。

首先，来问问你的私人顾问人工智能大模型kimi

kimi（https://kimi.moonshot.cn/）：两个不同数据集的单细胞降维聚类分群结果如何对应？

在单细胞转录组学研究中，将两个不同数据集的降维聚类分群结果进行对应是一个常见的问题，尤其是在跨样本、跨物种或跨实验条件的研究中。以下是几种常用的方法来实现这种对应关系：

1. 数据整合（Data Integration） 数据整合是最直接的方法之一，通过将两个数据集合并到一个统一的分析框架中，消除技术变异和批次效应，从而进行统一的降维和聚类。 操作步骤：

数据预处理：对两个数据集分别进行标准化、对数转换、高变基因筛选等预处理步骤。
整合工具选择：使用如Seurat的 SCTransform、Harmony、BBKNN或Scanorama等工具进行整合。
统一降维和聚类：整合后，对合并的数据集进行降维（如PCA、t-SNE或UMAP）和聚类。
分析整合结果：通过可视化（如UMAP图）和标记基因分析，确定聚类的细胞类型。

2. 基于标记基因的对应（Marker Gene Matching） 如果不想进行数据整合，可以分别对两个数据集进行降维和聚类，然后通过标记基因来寻找对应的细胞群。 操作步骤：

分别降维和聚类：对两个数据集分别进行降维和聚类。
标记基因分析：使用FindMarkers或FindAllMarkers函数找到每个聚类的标记基因。
标记基因匹配：比较两个数据集中聚类的标记基因，找到具有相似标记基因的聚类。

3. 基于细胞类型注释的对应（Cell Type Annotation） 如果已知某些标记基因或细胞类型特征，可以直接对两个数据集的聚类结果进行细胞类型注释，然后比较注释结果。 操作步骤：

分别降维和聚类：对两个数据集分别进行降维和聚类。
细胞类型注释：使用标记基因或自动化注释工具（如SingleR、scCATCH）对聚类结果进行细胞类型注释。
比较注释结果：比较两个数据集中相同细胞类型的聚类。

4. 基于相似性度量的对应（Similarity Metrics） 如果两个数据集的细胞类型较为复杂，可以使用相似性度量（如Jaccard指数）来量化聚类之间的相似性。 操作步骤：

分别降维和聚类：对两个数据集分别进行降维和聚类。
计算相似性：计算两个数据集中聚类之间的相似性（如Jaccard指数）。
匹配聚类：根据相似性得分找到最匹配的聚类。

总结选择哪种方法取决于具体的研究需求和数据特点：

数据整合：适合需要统一分析两个数据集的情况，能够消除批次效应。
标记基因匹配：适合已知标记基因且不想进行数据整合的情况。
细胞类型注释：适合已知细胞类型特征的情况。
相似性度量：适合细胞类型复杂且需要量化比较的情况。

在实际应用中，可以结合多种方法进行综合分析，以提高对应关系的准确性。

示例数据

学员询问的两个数据集分别为：GSE167297、GSE184198

GSE167297

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE167297，这个数据样本组织类型为弥漫型胃癌，利用10X Genomics技术测序，共30,365个细胞。

GSM5101013 Pt1_normal
GSM5101014 Pt1_superficial_cancer
GSM5101015 Pt1_deep_cancer
GSM5101016 Pt2_superficial_cancer
GSM5101017 Pt2_deep_cancer
GSM5101018 Pt3_normal
GSM5101019 Pt3_superficial_cancer
GSM5101020 Pt3_deep_cancer
GSM5101021 Pt4_normal
GSM5101022 Pt4_superficial_cancer
GSM5101023 Pt4_deep_cancer
GSM5101024 Pt5_normal
GSM5101025 Pt5_superficial_cancer
GSM5101026 Pt5_deep_cancer

对此数据进行简单的读取并创建Seurat对象：

rm(list=ls())
options(stringsAsFactors = F) 
source('scRNA_scripts/lib.R')
getwd()

###### step1: 导入数据 ######   
# 参考：https://mp.weixin.qq.com/s/tw7lygmGDAbpzMTx57VvFw
dir='GSE167297_RAW/' 
samples=list.files( dir ,pattern = 'gz')
samples 
library(data.table)
sceList = lapply(samples,function(pro){ 
  # pro=samples[1] 
  print(pro)
  ct=fread(file.path( dir ,pro),data.table = F)
  ct[1:4,1:4]
  rownames(ct)=ct[,1]
  ct=ct[,-1]
  sce=CreateSeuratObject(counts =  ct ,
                         project = gsub('_CountMatrix.txt.gz','', gsub('GSM[0-9]*_','',pro) ) ,
                         min.cells = 5,
                         min.features = 300)
  return(sce)
})
names(sceList)  
library(stringr)
do.call(rbind,lapply(sceList, dim)) 
samples= gsub('_CountMatrix.txt.gz','', gsub('GSM[0-9]*_','',samples) ) 
samples
sce.all=merge(x=sceList[[1]], y=sceList[ -1 ],   add.cell.ids = samples  ) 
names(sce.all@assays$RNA@layers)
sce.all[["RNA"]]$counts 

# Alternate accessor function with the same result
LayerData(sce.all, assay = "RNA", layer = "counts")
sce.all <- JoinLayers(sce.all)
dim(sce.all[["RNA"]]$counts )
 
as.data.frame(sce.all@assays$RNA$counts[1:10, 1:2])
head(sce.all@meta.data, 10)
table(sce.all$orig.ident)

然后经过简单的质控，降为聚类分群，得到以下的简单注释：

GSE184198

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE184198，该数据组织样本取自一名未经治疗、原发性、非转移性胃癌（GC）患者，该患者接受了胃癌切除手术。相邻的正常胃组织样本采集自癌组织边缘5厘米以外的区域。

GSM5580154 GC tissue
GSM5580155 NT tissue

数据处理：

rm(list=ls())
options(stringsAsFactors = F) 
source('scRNA_scripts/lib.R')
getwd()

###### step1: 导入数据 ######   
# 参考：https://mp.weixin.qq.com/s/tw7lygmGDAbpzMTx57VvFw
dir='GSE184198_RAW/'
samples=list.files( dir  )
samples
sceList = lapply(samples,function(pro){ 
  # pro=samples[5] 
  print(pro)  
  tmp = Read10X(file.path(dir,pro ) ) 
  if(length(tmp)==2){
    ct = tmp[[1]] 
  }else{ct = tmp}
  print(dim(ct))
  sce =CreateSeuratObject(counts =  ct ,
                          project =  gsub('_processed','', gsub('GSM[0-9]*_','',pro) )  ,
                          min.cells = 5,
                          min.features = 300 )
  return(sce)
}) 
do.call(rbind,lapply(sceList, dim)) 
samples=gsub('_processed','', gsub('GSM[0-9]*_','',samples) ) 
samples
sce.all=merge(x=sceList[[1]], y=sceList[ -1 ], add.cell.ids = samples  ) 
names(sce.all@assays$RNA@layers)
sce.all[["RNA"]]$counts 

# Alternate accessor function with the same result
LayerData(sce.all, assay = "RNA", layer = "counts")
sce.all <- JoinLayers(sce.all)
dim(sce.all[["RNA"]]$counts )
 
as.data.frame(sce.all@assays$RNA$counts[1:10, 1:2])
head(sce.all@meta.data, 10)
table(sce.all$orig.ident)

同样的简单注释之后结果如下：

将两个数据集的注释结果对应起来：基于亚群表达相似性

1、读取前面的处理结果

# GSE167297
sce.all = readRDS('./2021-GSE167297-深层浸润和浅层-胃癌/2-harmony/sce.all_int.rds')
load('./2021-GSE167297-深层浸润和浅层-胃癌/phe.Rdata') 
identical(rownames(phe) , colnames(sce.all))
sce.all@meta.data = phe
sel.clust = "celltype"
sce.all <- SetIdent(sce.all, value = sel.clust)
table(sce.all@active.ident) 
DimPlot(sce.all,label = T) 
GSE167297 = sce.all

# GSE184198
sce.all = readRDS('./2022-GSE184198-胃癌-肿瘤干细胞-AQP5/2-harmony/sce.all_int.rds')
load('./2022-GSE184198-胃癌-肿瘤干细胞-AQP5/phe.Rdata') 
identical(rownames(phe) , colnames(sce.all))
sce.all@meta.data = phe
sel.clust = "celltype"
sce.all <- SetIdent(sce.all, value = sel.clust)
table(sce.all@active.ident) 
DimPlot(sce.all, label = T) 
GSE184198 = sce.all

2、提取亚群水平的伪bulk count表达值

# 提取亚群水平的伪bulk count表达值
table(GSE184198$celltype)
GSE184198 <- AggregateExpression(GSE184198, group.by = c("celltype" ), assays = "RNA") 
GSE184198 <- as.data.frame(GSE184198[[1]])
head(GSE184198)

table(GSE167297$celltype)
GSE167297 <-AggregateExpression(GSE167297 , group.by = c("celltype" ), assays = "RNA") 
GSE167297 <- as.data.frame(GSE167297[[1]])
head(GSE167297)

3、根据表达值计算亚群间的相似性

使用pearson计算相关性并用热图展示：

ids <- intersect(rownames(GSE167297),rownames(GSE184198))
ct=cbind(GSE167297[ids,], GSE184198[ids,] )
colnames(ct)=paste(c(rep('GSE167297',ncol(GSE167297)),rep('GSE184198',ncol(GSE184198))),colnames(ct))
dim(ct)
head(ct)

pheatmap::pheatmap(cor(log2(edgeR::cpm(ct)+1)))