search

ixxmu / mp_duty

抓取网络文章到github issues保存
https://archives.duty-machine.now.sh/
115 stars 30 forks source link

基于singleR和人工审查的单细胞亚群命名比较 #2567

Closed ixxmu closed 2 years ago

ixxmu commented 2 years ago

https://mp.weixin.qq.com/s/347TOy8D4zZhlcddhjdpJQ

github-actions[bot] commented 2 years ago

基于singleR和人工审查的单细胞亚群命名比较 by 生信技能树

我们一直在安排学徒和实习生做单细胞数据集图表复现的腾讯会议分享，其中2021的100次单细胞分享，仅仅是30个成功录屏，已经编辑上传到b站：

  • https://www.bilibili.com/video/BV14b4y177EJ【生信技能树】单细胞文献分享计划11月
    • https://www.bilibili.com/video/BV1zY411L7Zg【生信技能树】单细胞文献分享计划1月
    • https://www.bilibili.com/video/BV1734y117f7【生信技能树】单细胞文献分享计划12月

因为大家拿到的都是我的标准代码，所以流程就是我五年前摸索成型的基于单细胞表达量矩阵的降维聚类分群啦，其中基于umap的不同亚群可视化其常见基因后，需要人工审核不同基因代表的生物学亚群然后确定名字这一点一直被大家诟病。因为确实工作量有点大，还需要背诵大量基因。比如 PBMC 的第一层次降维聚类分群后通常是可视化如下所示基因：

th=theme(axis.text.x = element_text(angle = 45, 
                                    vjust = 0.5, hjust=0.5))  
genes_to_check =list(
  Tcells=c('CD3D', 'CD3E', 'CD4','CD8A'),
  Bcells=c(  'CD19', 'CD79A', 'MS4A1'  ),
  plasma=c(    'CD38','TNFRSF17','IGHG1','IGHG4' ),
  Myeloids=c(   'CD68', 'CD163', 'CD14',   'S100A9', 'S100A8'),
  NK=c(      'KLRB1','NCR1') 
)
p_all_markers=DotPlot(sce.all, 
                      features = genes_to_check,
                      scale = T,assay='RNA' )+ th +
  theme(axis.text.x=element_text(angle=45,hjust = 1))
p_all_markers
ggsave('check_paper_markers.pdf',
       height = 8,width = 10)

这样的话， 肿瘤第一层次分群，肿瘤基质细胞里面的成纤维内皮周细胞平滑肌细胞的细分，PBMC里面的Myeloids继续细分，都需要背诵或者收集整理大量的基因。有了大量已知基因的背景铺垫才有可能更好的给我们的单细胞降维聚类分群结果进行合理的命名！

恰好最近学徒就咨询了中大五院泌尿外科赵亮宇博士后所在的中大五院/上海交通大学附属第一人民医院联合课题组在国际顶尖期刊《Nature Communications》发表了关于人阴茎海绵体单细胞的研究论文里面的关于成纤维内皮周细胞平滑肌细胞的细分问题。

《Nature Communications》

文章的第一层次降维聚类分群主要是成纤维，内皮，周细胞，平滑肌细胞的细分：

成纤维，内皮，周细胞，平滑肌细胞的细分

可以看到， 其实周细胞，平滑肌细胞在第一层次虽然umap上面可以说是泾渭分明，但是两个亚群的top基因并不是非此即彼的关系而是高低强弱的关系。这样就很难通过是背诵基因的方式来进行第一层次的分群，文章里面选取了 KCNJ8和ACTA2进行代表。

  • VWF (endothelial cells)
  • PDGFRA (fibroblast)
  • KCNJ8 (pricytes)
  • ACTA2 (smooth muscle cells)
  • S100B (Schwann cells)
  • CD163 (macro- phages)
  • CD3E (T cells).

这个时候很难判断中大五院泌尿外科赵亮宇博士后做这个分析的时候使用的是基于singleR还是人工审查的单细胞亚群命名方式。因为后续在继续对成纤维，内皮，周细胞，平滑肌细胞的细分的时候，文章里面的平滑肌细胞的细分里面又一次出现了周细胞，这个时候的RGS5基因是我们比较熟悉的了。

平滑肌细胞的细分里面又一次出现了周细胞

我看到了这个文章的表达量矩阵（GSE206528）是公开的可以获取的，所以去下载了玩了一下。一般来说单细胞转录组测序数据走完cellranger的定量流程，每个样品就会得到3个表达量矩阵文件（barcodes.tsv.gz，matrix.mtx.gz，genes.tsv.gz或者features.tsv.gz），然后就可以走seurat流程进行单细胞降维聚类分群。这样的基础分析，也可以看基础10讲：

不过，这个文章并不是每个单细胞样品3个文件，而是独立的一个文件，所以批量读取作者上传的（GSE206528）每个样品的gene_matrix.csv.gz文件，构建成为了基于Seurat的单细胞对象：

rm(list=ls())
options(stringsAsFactors = F) 
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(dplyr)

###### step1:导入数据 ######   
library(data.table)
dir='GSE206528_RAW/' 
samples=list.files( dir )
samples
# samples = head(samples,10)

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]
  ct[nrow(ct),1:4]
  rownames(ct)=ct[,1]
  ct=ct[,-1]
  ct[1:4,1:4]
  sce =CreateSeuratObject(counts =  ct ,
                          project =   gsub('_gene_matrix.csv.gz','',gsub('^GSM[0-9]*_','',pro) ) ,
                          min.cells = 5,
                          min.features = 300 )
  return(sce)
})
names(sceList)  


samples
# gsub('^GSM[0-9]*','',samples)

sce.all=merge(x=sceList[[1]],
              y=sceList[ -1 ],
              add.cell.ids = gsub('_gene_matrix.csv.gz','',gsub('^GSM[0-9]*_','',samples) )  )

as.data.frame(sce.all@assays$RNA@counts[1:10, 1:2])
head(sce.all@meta.data, 10)
table(sce.all$orig.ident)

然后走harmony的整合，以及单细胞降维聚类分群：

如下所示：

 

这些基因当然是不够的，所以需要加入 成纤维，内皮，周细胞，平滑肌细胞的细分相关基因列表，然后人工审查，我肉眼进行了如下所示的分类：


# 需要自行看图，定细胞亚群： 

celltype=data.frame(ClusterID=0:25 ,
                    celltype= 0:25) 
#定义细胞亚群   
celltype[celltype$ClusterID %in% c( 10,22,24),2]='Mac' 
celltype[celltype$ClusterID %in% c( 11,23 ),2]='Tcells'  

celltype[celltype$ClusterID %in% c( 0,2,5,13,15,17,20 ),2]='Endo'  

celltype[celltype$ClusterID %in% c(  1,3,4,6,8,12,15,16,18,21 ),2]='Fibro'  
celltype[celltype$ClusterID %in% c( 7,14 ),2]='Pericyte'   
celltype[celltype$ClusterID %in% c(  9,19,25 ),2]='SMC'  

人工审查然后进行命名确实在很多亚群都非常纠结，尤其是 周细胞，平滑肌细胞 ，它们在7,14以及 9,19,25 都有RGS5基因表达量，仅仅是这个基因压根就没办法把大家合理的区分 ：

都有RGS5基因表达量

如果我们使用singleR就简单很多了，当然了前提是自己制作好了  ref_sce ，这里省略（步骤略微有点复杂）：

testdata <- GetAssayData(sce.all, slot="data")
testdata[1:4,1:4]
dim(testdata)
library(SingleR)
pred <- SingleR(test=testdata, ref=ref_sce, 
                labels=ref_sce$Type
)
as.data.frame(table(pred$labels))
head(pred) 
labels=pred$labels
table(labels)  
sce.all$singleR=labels
DimPlot(sce.all, reduction = "umap",repel = T,
        group.by = "singleR",label = T)
ggsave('umap-by-singleR.pdf')



pdf('singleR-vs-RNA_snn_res.0.8.pdf',width = 12)
gplots::balloonplot(table(sce.all$RNA_snn_res.0.8,sce.all$singleR))
dev.off()

load('../3-cell/phe-by-markers.Rdata')
pdf('singleR-vs-RNA_snn_res.0.8.pdf',width = 12)
gplots::balloonplot(table(phe$celltype,sce.all$singleR))
dev.off()

可以看到如果使用singleR进行单细胞亚群注释，其实都不需要走降维聚类分群流程了，仅仅是最开始读取表达量矩阵即可。

让我们来看看两个方法的注释结果差异吧，如下所示：

两个方法的注释结果差异

在判断 内皮和成纤维的时候，两个方法基本上没有差异，因为确实它们的基因很特异性，如下所示：

LUM和DCN的成纤维特异性

但是在判断周细胞和平滑肌细胞的时候，两个方法出现了冲突，但是这个冲突是可以容忍的，是第一层次降维聚类分群目前不可能完美避开的。

基于singleR和人工审查的单细胞亚群命名都有自己的合理性，亲爱的读者朋友们，你们会如何选择呢？


写在文末

我在《生信技能树》，《生信菜鸟团》，《单细胞天地》的大量推文教程里面共享的代码都是复制粘贴即可使用的， 有任何疑问欢迎留言讨论，也可以发邮件给我，详细描述你遇到的困难的前因后果给我，我的邮箱地址是 jmzeng1314@163.com

如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：

We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.

十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你。