提供了编程的基础技术教程

zl程序教程

您现在的位置是：首页 >  IT要闻

当前栏目

在你的髓系里面加上中性粒吧

2023-02-18 16:30:34 时间

最近肿瘤单细胞数据分析里面中性粒细胞亚群比较火爆，不过绝大部分10x技术产出的单细胞转录组数据里面其实很难区分出来中性粒细胞亚群，具体原因大家很容易去10x的官网看到。不过，我回顾了一下不久前的一个学徒汇报的公共数据挖掘，恰好里面的髓系免疫细胞里面就有这个中性粒细胞亚群。

通常我们拿到了肿瘤相关的单细胞转录组的表达量矩阵后的第一层次降维聚类分群通常是：

  • immune (CD45+,PTPRC),
  • epithelial/cancer (EpCAM+,EPCAM),
  • stromal (CD10+,MME,fibo or CD31+,PECAM1,endo)

参考我前面介绍过 CNS图表复现08—肿瘤单细胞数据第一次分群通用规则，这3大单细胞亚群构成了肿瘤免疫微环境的复杂。绝大部分文章都是抓住免疫细胞亚群进行细分，包括淋巴系（T,B,NK细胞）和髓系（单核，树突，巨噬，粒细胞）的两大类作为第二次细分亚群。但是也有不少文章是抓住stromal 里面的fibo 和endo进行细分，并且编造生物学故事的。

我们前面已经对免疫细胞里面的髓系和B细胞细分亚群进行了简单的介绍：

但是之前的髓系免疫细胞细分的时候其实并没有中性粒细胞亚群，它也一直不在我们分享的单细胞数据处理代码里面。是时候补充起来了，首先我对单细胞转录组表达量矩阵进行了标准的降维聚类分群，然后提取了里面的髓系免疫细胞后继续进行标准的降维聚类分群，如下所示：

标准的降维聚类分群

根据我的生物学背景，以及对它们的每个单细胞亚群的理解，我给出来了如下所示的细胞亚群名字：

celltype=data.frame(ClusterID=0:13,
                    celltype=0:13 )   

#celltype[celltype$ClusterID %in% c( 9,18 ),2]='mast'  
celltype[celltype$ClusterID %in% c(  11  ),2]='pDC'  
celltype[celltype$ClusterID %in% c( 12  ),2]='DC1'  
celltype[celltype$ClusterID %in% c( 3 ),2]='DC2'   
celltype[celltype$ClusterID %in% c( 13 ),2]='DC3'   

celltype[celltype$ClusterID %in% c( 1,2,4,5,7 ),2]='Mac'     
celltype[celltype$ClusterID %in% c( 0),2]='mono'  
celltype[celltype$ClusterID %in% c( 6 ),2]='neutrophils'

非常的漂亮，最后得到的髓系免疫细胞各个单细胞亚群及其数量如下所示：

> as.data.frame(sort(table(Idents(sce.all))))
         Var1 Freq
1         DC3   17
2         DC1   27
3         pDC   38
4 neutrophils  176
5         DC2  237
6        mono  533
7         Mac 1218

可以看到，主要是巨噬细胞和单核细胞，因为第一层次降维聚类分群里面的肥大细胞非常清晰的跟髓系其它单细胞亚群区分开来，所以第二层次并没有考虑肥大细胞。

分配好群之后，可以看到每个单细胞亚群各自的高表达量基因非常清晰：

各自的高表达量基因非常清晰

其中DC2这个亚群的异质性稍微有一点点大，这个我做了六百多个肿瘤单细胞转录组数据集发现绝大部分都是这样，就先忽略它。

值得注意的是这个时候中性粒细胞亚群已经是非常清晰可见啦，而且确实细胞质量堪忧：

细胞质量堪忧

而且我把这些髓系单细胞亚群的基因也提取出来了，做成为了代码给大家：

th=theme(axis.text.x = element_text(angle = 45, 
                                    vjust = 0.5, hjust=0.5))  
myeloids = list(
  Mac=c("C1QA","C1QB","C1QC","SELENOP","RNASE1","DAB2","LGMN","PLTP","MAF","SLCO2B1"),
  mono=c("VCAN","FCN1","CD300E","S100A12","EREG","APOBEC3A","STXBP2","ASGR1","CCR2","NRG1"),
  neutrophils =  c("FCGR3B","CXCR2","SLC25A37","G0S2","CXCR1","ADGRG3","PROK2","STEAP4","CMTM2" ),
  pDC = c("GZMB","SCT","CLIC3","LRRC26","LILRA4","PACSIN1","CLEC4C","MAP1A","PTCRA","C12orf75"),
  DC1 = c("CLEC9A","XCR1","CLNK","CADM1","ENPP1","SNX22","NCALD","DBN1","HLA-DOB","PPY"),
  DC2=c( "CD1C","FCER1A","CD1E","AL138899.1","CD2","GPAT3","CCND2","ENHO","PKIB","CD1B"),
  DC3 =  c("HMSD","ANKRD33B","LAD1","CCR7","LAMP3","CCL19","CCL22","INSM1","TNNT2","TUBB2B")
)
p <- DotPlot(sce.all, features = myeloids,
             assay='RNA' ,group.by = 'celltype' )  +th

p
ggsave(plot=p, filename="check_myeloids_marker_by_celltype.pdf")

蛮方便的，大家可以复制粘贴我的代码去自己的单细胞项目里面的髓系免疫细胞里面赶快看看哦。

值得注意的是，因为这个数据集的特殊性，我这里并没有把单核细胞区分成为经典的CD14单核和 非经典的CD16单核哦。

当然了，这个数据集还可以被制作成为了singleR的数据库：

av_myeloids <-AverageExpression(sce.all,
                       group.by =c( "celltype"),
                       assays = "RNA")[[1]]
head(av_myeloids)  
ref_sce=SingleCellExperiment::SingleCellExperiment(assays=list(counts=av_myeloids))
library(scater)
ref_sce=scater::logNormCounts(ref_sce)
library(SingleCellExperiment)
logcounts(ref_sce)[1:4,1:4]
colData(ref_sce)$Type=colnames(av_myeloids) 
ref_sce_myeloids = ref_sce 
save(av_myeloids, ref_sce_myeloids,file = 'av_myeloids_for_singleR.Rdata')

只要你拿到了我的av_myeloids_for_singleR.Rdata 文件，以后遇到了髓系免疫细胞就可以很容易去注释啦。

testdata <- GetAssayData(sce.all, slot="data")
testdata[1:4,1:4]
dim(testdata)
library(SingleR)
load(file = 'av_myeloids_for_singleR.Rdata')
pred <- SingleR(test=testdata, ref=ref_sce, 
                labels=ref_sce$Type
) 
sce.all$singleR=pred$labels
DimPlot(sce.all, reduction = "umap",repel = T,
        group.by = "singleR",label = T)
ggsave('umap-singleR.pdf',width = 12,height = 6)
gplots::balloonplot(table(Idents(sce.all),pred$labels))

不过，让我诧异的是，如果是使用SingleR针对每个细胞进行相关性注释，会有一些模糊地带，如下所示：

有一些模糊地带

这就需要对SingleR软件有所理解啦，如果是按照亚群来注释虽然会很准确，但是仁者见仁智者见智啦。

另外，如果你需要我的av_myeloids_for_singleR.Rdata 文件，可以发邮件找我申请哈，来信需要注明身份，你遇到的困难，另外我也不保证时效性哈。

相关文章