血液免疫细胞比例估算准确度如何？

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

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

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

而且肿瘤微环境里面的免疫微环境越来越受大家关注，很自然的大家都喜欢拿肿瘤单细胞数据集的降维聚类分群结果去看常规的bulk肿瘤转录组测序矩阵里面的肿瘤微环境细胞亚群比例构成。前面我们虽然介绍了 使用MuSiC以及MuSiC2来根据单细胞转录组结果推断bulk转录组细胞比例，但是仅仅是使用了其示范数据，而且是糖尿病领域，对绝大部分肿瘤研究小伙伴来说有点陌生，所以我们接下来以Seurat里面的单细胞降维聚类分群结果和bulk转录组，结合起来继续说明如何 使用MuSiC以及MuSiC2来根据单细胞转录组结果推断bulk转录组细胞比例。

因为是血液免疫，所以免疫细胞亚群进行细分就是，包括淋巴系（T,B,NK细胞）和髓系（单核，树突，巨噬，粒细胞）的两大类。

血液免疫细胞常规bulk转录组数据（GEO数据集）

/geo/query/acc.cgi?acc=GSE74246

样品队列是 13 normal hematopoietic cell types and 3 acute myeloid leukemia cell types ，每个病人都有很多细胞亚群的纯化后的常规bulk转录组数据：

GSM1915582 2596-Bcell
GSM1915583 2596-CD4Tcell
GSM1915584 2596-CD8Tcell
GSM1915585 2596-CLP
GSM1915586 2596-Ery
GSM1915587 2596-NKcell

而且作者提供了 1.9M的 GSE74246_RNAseq_All_Counts.txt.gz 表达量矩阵文件下载，下载链接是：

• /geo/series/GSE74nnn/GSE74246/suppl/GSE74246_RNAseq_All_Counts.txt.gz

血液免疫细胞的单细胞转录组（seurat-data和TENxPBMCData包）

们会推荐 satijalab/seurat-data ，它内置了很多数据集，如果你还没有下面的seurat-data包和pbmc3k对象 ，就自己去下载：

install.packages('devtools')
devtools::install_github('satijalab/seurat-data')
library(SeuratData) #加载seurat数据集  
getOption('timeout')
options(timeout=10000)
InstallData("pbmc3k")  
data("pbmc3k")  
sce <- pbmc3k.final  

但是因为它这个包仍然是在GitHub，所以很多人下载失败， 另外它数据也是在外网，很多人网络有问题无法下载。但是我随时搜索到了一个包，TENxPBMCData ，蛮有意思的，安装和加载：

# BiocManager::install('TENxPBMCData') # (286 KB)
library(TENxPBMCData)
args(TENxPBMCData) 

因为这个包，TENxPBMCData 并不大，所以它其实也是需要在线下载数据的。目前它有如下所示的10个数据：

1. pbmc68k
2. frozen_pbmc_donor_a
3. frozen_pbmc_donor_b
4. frozen_pbmc_donor_c
5. pbmc33k
6. pbmc3k
7. pbmc6k
8. pbmc4k
9. pbmc8k
10. pbmc5k-CITEseq

每个数据集都有自己的代号：

function (dataset = c("pbmc4k", "pbmc68k", "frozen_pbmc_donor_a", 
    "frozen_pbmc_donor_b", "frozen_pbmc_donor_c", "pbmc33k", 
    "pbmc3k", "pbmc6k", "pbmc8k", "pbmc5k-CITEseq"), as.sparse = TRUE) 

运行MuSiC计算比例

首先读取前面的 1.9M的 GSE74246_RNAseq_All_Counts.txt.gz 表达量矩阵文件：

rm(list = ls())
library(MuSiC)
library(ggplot2)
library(xbioc)
library(dplyr)
library(gridExtra)
library(grid)
library(MuSiC2)

a=data.table::fread('GSE74246_RNAseq_All_Counts.txt.gz',data.table = F)
length(unique(a[,1]))
rownames(a)=a[,1]
a=a[,-1]
phe = as.data.frame(stringr::str_split(colnames(a),'[-_]',simplify = T))
head(phe)
table(phe$V2)


kp = phe$V2 %in% c('Bcell', 'CD4Tcell', 'CD8Tcell' , 'Mono','NKcell' )
a=a[,kp]
phe=phe[kp,]

这个时候因为数据集GSE74246涉及到的免疫细胞种类实在是有点多，我们仅仅是拿PBMC3K单细胞数据集里面的即可：

# install.packages('devtools')
# devtools::install_github('satijalab/seurat-data')
library(SeuratData) #加载seurat数据集  
library(Seurat)
getOption('timeout')
options(timeout=10000)
#InstallData("pbmc3k")  
data("pbmc3k")  
sce <- pbmc3k.final  
table(Idents(sce))

有了PBMC3K单细胞数据集，而且有了常规的bulk转录组测序矩阵，很容易 使用MuSiC以及MuSiC2来根据单细胞转录组结果推断bulk转录组细胞比例。

library(scater)
pbmc3k.sce <- SingleCellExperiment(
  assays = list(counts =  sce@assays$RNA@counts), 
  colData = data.frame(cellType=Idents(sce),sampleID=colnames(sce))
)
pbmc3k.sce
colData(pbmc3k.sce)
unique(pbmc3k.sce$cellType) 


# music2 deconvolution
set.seed(1234)
# MuSiC
bulk.mtx = as.matrix(a)
prop_music=music_prop(bulk.mtx = bulk.mtx, sc.sce = pbmc3k.sce, 
                      clusters = 'cellType', samples = 'sampleID',
                      select.ct = unique(pbmc3k.sce$cellType) [1:5],
                      verbose = F)$Est.prop.weighted
pheatmap::pheatmap(prop_music)
pheatmap::pheatmap(prop_music,display_numbers = T)

可以很清楚的看到，这个细胞比例的推断还是蛮成功的 :

因为这个常规的bulk转录组测序矩阵本来就是每个人的血液的各种免疫细胞亚群纯化后的，所以纯化后的样品去卷积推测其单细胞亚群比例就基本上也是比较纯的。