肿瘤微环境（TME）大家都比较熟悉了，肿瘤周围的组织、免疫细胞、血管、细胞外基质等元素共同形成了“肿瘤微环境”。我们讲了很多工具可以推断bulk转录组的微环境组成，目录是：
- estimate的两个打分值本质上就是两个基因集的ssGSEA分析
- 针对TCGA数据库全部的癌症的表达量矩阵批量运行estimate
- 不同癌症内部按照estimate的两个打分值高低分组看蛋白编码基因表达量差异
- 使用CIBERSORT算法推断全部tcga样品的免疫细胞比例
这些工具都是依据肿瘤病人的转录组测序表达量矩阵进行的分析，也有几百篇类似的数据挖掘文章了，它们总是喜欢落脚到estimate或者CIBERSORT结果的预后意义。
虽然这些工具肿瘤微环境（TME）的思路是ok的，绝大部分的肿瘤相关的单细胞转录组研究我介绍过 CNS图表复现08—肿瘤单细胞数据第一次分群通用规则，这个第一次分群规则是 ：
- immune (CD45 ,PTPRC),
- epithelial/cancer (EpCAM ,EPCAM),
- stromal (CD10 ,MME,fibo or CD31 ,PECAM1,endo)
这3大单细胞亚群构成了肿瘤免疫微环境的复杂。绝大部分文章都是抓住免疫细胞亚群进行细分，包括淋巴系（T,B,NK细胞）和髓系（单核，树突，巨噬，粒细胞）的两大类作为第二次细分亚群。但是也有不少文章是抓住stromal 里面的fibo 和endo进行细分，并且编造生物学故事的。
而TCGA等公共数据库数据库的转录组测序数据毕竟是bulk转录组测序，病人的肿瘤样品里面虽然是混合了各种各样的肿瘤微环境里面的基质细胞和免疫细胞，但是在数据层面被混杂成为了一个样品，并不是单细胞测序，所以并没有细胞比例信息。这个时候如果想推断细胞比例就需要借助于算法啦，有一个2019的综述文章：《Comprehensive evaluation of transcriptome-based cell-type quantification methods for immuno-oncology》比较了常见的免疫细胞比例推断工具的表现。
但是这些工具都是在单细胞还没有流行开来的时候开发工具，虽然说它们确实是可以跑通流程，也可以进行细胞亚群比例推断。而现在各个疾病研究领域的单细胞转录组公开数据多如牛毛，我们自己对单细胞转录组数据的降维聚类分群和命名后的信息，如果可以用来推断bulk转录组细胞比例会更加精准。下面我们就介绍一下使用MuSiC以及MuSiC2来根据单细胞转录组结果推断bulk转录组细胞比例。
首先是安装MuSiC以及MuSiC2
比较麻烦的是MuSiC以及MuSiC2目前都是在github的包，所以安装起来有点麻烦，代码如下所示：
代码语言：javascript复制# BiocManager::install('TOAST')
# devtools::install_github( repo = "crhisto/xbioc")
# install devtools if necessary
if (!"devtools" %in% rownames(installed.packages())) {
install.packages('devtools')
}
# install MuSiC package
if (!"MuSiC" %in% rownames(installed.packages())) {
devtools::install_github('xuranw/MuSiC')
}
# install the MuSiC2 package
if (!"MuSiC2" %in% rownames(installed.packages())) {
devtools::install_github('Jiaxin-Fan/MuSiC2')
}
MuSiC以及MuSiC2目前都是在github的包，大家完全是可以下载其github的源代码压缩包，里面也有示例数据。分别有 bulk-eset.rds 和 EMTABesethealthy.rds
认识示例数据（bulk转录组矩阵和单细胞矩阵）
示例数据文件是 bulk-eset.rds 和 EMTABesethealthy.rds ，可以在MuSiC以及MuSiC2的github源代码找到。然后分别加载后熟悉一下：
代码语言：javascript复制## Bulk RNA-seq data
#benchmark.eset = readRDS("https://xuranw.github.io/MuSiC/data/bulk-eset.rds")
benchmark.eset = readRDS("bulk-eset.rds")
benchmark.eset
table(benchmark.eset$group)
bulk.control.mtx = exprs(benchmark.eset)[, benchmark.eset$group == 'healthy']
bulk.case.mtx = exprs(benchmark.eset)[, benchmark.eset$group == 't2d']
bulk.case.mtx[1:4,1:4]
可以看到的是bulk转录组矩阵存储成为了 ExpressionSet 对象，需要自己去熟悉它。如下所示：
代码语言：javascript复制> benchmark.eset
ExpressionSet (storageMode: lockedEnvironment)
assayData: 10000 features, 200 samples
element names: exprs
protocolData: none
phenoData
sampleNames: 1 10 ... 200 (200 total)
varLabels: sampleID group
varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation:
> table(benchmark.eset$group)
healthy t2d
100 100
> bulk.case.mtx[1:4,1:4]
101 102 103 104
SFT2D1 16791 8362 17888 12894
TRAPPC13 8583 7823 6674 6626
C9orf89 6745 4659 8482 4790
RPS6 218932 142579 190503 160069
这样的转录组测序矩阵大家很容易自己拿到，看起来就是普普通通的counts矩阵，虽然bulk转录组矩阵存储成为了 ExpressionSet 对象，但是后续在使用MuSiC以及MuSiC2需要的都是从 ExpressionSet 对象里面拿到的普普通通的counts矩阵即可哦。
然后熟悉一下单细胞表达量矩阵：
代码语言：javascript复制## scRNA-seq data
#seger.sce = readRDS("https://xuranw.github.io/MuSiC/data/EMATBsce_healthy.rds")
seger.sce = readRDS("EMTABesethealthy.rds")
seger.sce
# sc.sce SingleCellExperiment for single cell data
tail(sort(table(seger.sce$cellType)))
# 需要熟练掌握它们的对象，：一些单细胞转录组R包的对象
# https://mp.weixin.qq.com/s/lpoHhZqi-_ASUaIfpnX96w
seger.sce <- SingleCellExperiment(
assays = list(counts = exprs(seger.sce) ),
colData = seger.sce@phenoData@data
)
colData(seger.sce)
值得玩味的是作者给出来的示例数据里面的单细胞表达量矩阵居然也是一个 ExpressionSet 对象，其实对它的MuSiC以及MuSiC2来说是不合法的，我给它修改成为了SingleCellExperiment对象， 就合情合理啦。
代码语言：javascript复制> seger.sce
ExpressionSet (storageMode: lockedEnvironment)
assayData: 25453 features, 1097 samples
element names: exprs
protocolData: none
phenoData
sampleNames: AZ_A10 AZ_A11 ... HP1509101_P9 (1097 total)
varLabels: sampleID SubjectName cellTypeID cellType
varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation:
> tail(sort(table(seger.sce$cellType)))
delta gamma acinar ductal beta alpha
59 75 112 135 171 443
> seger.sce
class: SingleCellExperiment
dim: 25453 1097
metadata(0):
assays(1): counts
rownames(25453): SGIP1 AZIN2 ... KIR2DL2 KIR2DS3
rowData names(0):
colnames(1097): AZ_A10 AZ_A11 ... HP1509101_P8 HP1509101_P9
colData names(4): sampleID SubjectName cellTypeID cellType
reducedDimNames(0):
mainExpName: NULL
altExpNames(0):
可以看到其实单细胞转录组数据分析，就是对各种各样的R对象的理解而已。
运行MuSiC并且理解结果
其实就是一个平平无奇的 music_prop 函数而已：
代码语言：javascript复制# MuSiC
bulk.mtx = exprs(benchmark.eset)
prop_music=music_prop(bulk.mtx = bulk.mtx, sc.sce = seger.sce,
clusters = 'cellType', samples = 'sampleID',
select.ct = c('acinar','alpha', 'beta', 'delta', 'ductal','gamma'),
verbose = F)$Est.prop.weighted
prop_all = cbind('proportion'=c(prop_music),
'celltype'=rep(colnames(prop_music),
each=nrow(prop_music)),
'sampleID'=rep(rownames(prop_music),times=ncol(prop_music)),
'Method'='MuSiC')
prop_all=as.data.frame(prop_all)
prop_all$proportion=as.numeric(as.character(prop_all$proportion))
就得到了需要分解的单细胞亚群（ delta gamma acinar ductal beta alpha ）在每个样品的比例情况；
运行MuSiC2并且理解结果
也是一个平平无奇的函数 music2_prop_t_statistics ：
代码语言：javascript复制
# music2 deconvolution
set.seed(1234)
library(scater)
est = music2_prop_t_statistics(bulk.control.mtx = bulk.control.mtx,
bulk.case.mtx = bulk.case.mtx,
sc.sce = seger.sce,
clusters = 'cellType',
samples = 'sampleID',
select.ct = c('acinar','alpha','beta','delta','ductal','gamma'),
n_resample=20, sample_prop=0.5,cutoff_c=0.05,cutoff_r=0.01)
est.prop = est$Est.prop
# plot estimated cell type proportions
prop_all = cbind('proportion'=c(est.prop), 'sampleID'=rep(rownames(est.prop),times=ncol(est.prop)), 'celltype'=rep(colnames(est.prop), each=nrow(est.prop)))
prop_all = as.data.frame(prop_all)
prop_all$proportion = as.numeric(as.character(prop_all$proportion))
prop_all$group = ifelse(prop_all$sampleID %in% seq(from=1, to=100, by=1), 'Healthy', 'T2D')
就得到了需要分解的单细胞亚群（ delta gamma acinar ductal beta alpha ）在每个样品的比例情况，可以看到它的区别居然是把需要分解的bulk转录组矩阵根据其表型分开了一下。
如果你对两个函数（music_prop和music2_prop_t_statistics）有疑惑，可以去细看官方文档：https://xuranw.github.io/MuSiC/articles/pages/MuSiC2.html
用法还是蛮简单的，自己准备好单细胞转录组矩阵以及bulk转录组矩阵即可，运行的速度也是很快；
有了需要分解的单细胞亚群（ delta gamma acinar ductal beta alpha ）在每个样品的比例情况，就可以很轻松的可视化，或者分组比较细胞比例差异：
细胞比例差异
可以看到在糖尿病组里面的 acinar单细胞亚群应该是统计学显著的比例上升了。
