RNA 25. SCI文章中估计组织浸润免疫细胞和基质细胞群的群体丰度(MCP-counter)

点击关注，桓峰基因

今天来介绍一个利用基因表达估计组织浸润免疫细胞和基质细胞群的群体丰度的软件包——MCP-counter，亲试，非常好用。

桓峰基因的教程不但教您怎么使用，还会定期分析一些相关的文章，学会教程只是基础，但是如果把分析结果整合到文章里面才是目的，觉得我们这些教程还不错，并且您按照我们的教程分析出来不错的结果发了文章记得告知我们，并在文章中感谢一下我们哦！

公司英文名称：Kyoho Gene Technology (Beijing) Co.,Ltd.

桓峰基因公众号推出转录组分析和临床预测模型教程，有需要生信的老师可以联系我们！首选看下转录分析教程整理如下：

RNA 1. 基因表达那些事–基于 GEO

RNA 2. SCI文章中基于GEO的差异表达基因之 limma

RNA 3. SCI 文章中基于T CGA 差异表达基因之 DESeq2
RNA 4. SCI 文章中基于TCGA 差异表达之 edgeR
RNA 5. SCI 文章中差异基因表达之 MA 图
RNA 6. 差异基因表达之-- 火山图 (volcano)
RNA 7. SCI 文章中的基因表达——主成分分析 (PCA)
RNA 8. SCI文章中差异基因表达–热图 (heatmap)

RNA 9. SCI 文章中基因表达之 GO 注释

RNA 10. SCI 文章中基因表达富集之–KEGG
RNA 11. SCI 文章中基因表达富集之 GSEA

RNA 12. SCI 文章中肿瘤免疫浸润计算方法之 CIBERSORT

RNA 13. SCI 文章中差异表达基因之 WGCNA

RNA 14. SCI 文章中差异表达基因之蛋白互作网络 (PPI)

RNA 15. SCI 文章中的融合基因之 FusionGDB2

RNA 16. SCI 文章中的融合基因之可视化

RNA 17. SCI 文章中的筛选 Hub 基因 (Hub genes)

RNA 18. SCI 文章中基因集变异分析 GSVA

RNA 19. SCI 文章中无监督聚类法（ConsensusClusterPlus）

RNA 20. SCI 文章中单样本免疫浸润分析（ssGSEA）
RNA 21. SCI 文章中单基因富集分析

RNA 22. SCI 文章中基于表达估计恶性肿瘤组织的基质细胞和免疫细胞（ESTIMATE）

RNA 23. SCI文章中表达基因模型的风险因子关联图（ggrisk）

RNA 24. SCI文章中基于TCGA的免疫浸润细胞分析 (TIMER)

RNA 25. SCI文章中估计组织浸润免疫细胞和基质细胞群的群体丰度(MCP-counter)

临床预测模型整理如下：

Topic 1. _临床_标志物生信分析常规思路

Topic 2. 生存分析之 Kaplan-Meier

Topic 3. SCI文章第一张表格–基线表格

Topic 4. _临床_预测模型构建 Logistic 回归

Topic 5. 样本量确定及分割

Topic 6 计数变量泊松回归

Topic 7. _临床_预测模型–Cox回归

Topic 8. _临床_预测模型-Lasso回归

Topic 9. SCI 文章第二张表—单因素回归分析表

Topic 10. 单因素 Logistic 回归分析—单因素分析表格

Topic 11. SCI中多元变量筛选—单/多因素表
Topic 12 _临床_预测模型—列线表 (Nomogram)

Topic 13. _临床_预测模型—一致性指数 (C-index)

Topic 14. _临床_预测模型之校准曲线 (Calibration curve)
Topic 15. _临床_预测模型之决策曲线 (DCA)
Topic 16. _临床_预测模型之接收者操作特征曲线 (ROC)

Topic 17. 临床预测模型之缺失值识别及可视化

Topic 18. 临床预测模型之缺失值插补方法

前言

我们引入了微环境细胞种群计数器(MCP-counter)方法，它允许鲁棒性8种免疫细胞和2种基质细胞异质群体绝对丰度的定量来自转录组数据的组织。我们提供了体外mRNA混合物和体外免疫组化数据定量地支持我们方法估计的有效性。此外，我们演示了MCP计数器克服了以前提出的计算方法的一些限制或弱点。MCP-counter是用于绘制人体健康组织和非造血人体免疫浸润的全局图肿瘤和总结了与肺部总生存率相关的微环境患者分层腺癌，结直肠癌和乳腺癌。

MCP计数器方法的目的和发展。比较MCP-counter估计和cibersort估计的白细胞分数, MCP计数器的开发和验证框架。微环境细胞种群计数器（MCP-count）方法，该方法允许从转录组数据对组织中八种免疫细胞和两种基质细胞种群的绝对丰度进行量化。离体免疫组织化学数据支持该方法的有效性。因此此MCP-counter可以用于绘制人类健康组织和非造血人类肿瘤的免疫浸润的全局图。MCP-counter也提供R包。从基因表达矩阵中，它为每个样本生成 CD3+ T 细胞、CD8+ T 细胞、细胞毒性淋巴细胞、NK 细胞、B 淋巴细胞、源自单核细胞（单核细胞谱系）的细胞、髓样树突细胞、中性粒细胞以及内皮细胞和成纤维细胞。MCP-counter 是“单样本”分数，因为它们是在每个样本上独立计算的。然后，这些分数可用于直接比较队列中样本中相应细胞类型的丰度。通过使用石蜡包埋组织切片上免疫组织化学细胞定量对 MCP 计数器进行了定量验证。结果说明了它在 47 种健康组织类型和 32 种非血液系统恶性肿瘤中评估组织浸润的成功应用。

软件安装

通过加载devtools软件包进行安装，如下：

if(!require(devtools))install.packages("devtools")
if(!require(curl))install.packages("curl")library(devtools)
library(curl)
if(!require(MCPcounter))install_github("ebecht/MCPcounter",ref="master", subdir="Source")

数据读取

我们需要准备三个文件，包括表达矩阵，探针文件，以及对应基因文件，下面我们分析看看文件格式。

1. 探针文件

这个文件我们可以从官方直接获取，文件不是很大，直接下载也不会出问题。

library(curl)
library(MCPcounter)
probesets = read.table(curl("http://raw.githubusercontent.com/ebecht/MCPcounter/master/Signatures/probesets.txt"),sep = "\t", stringsAsFactors = FALSE, colClasses = "character")head(probesets)
##             V1      V2
## 1 1557733_a_at T cells
## 2   1558971_at T cells
## 3 1558972_s_at T cells
## 4   1566448_at T cells
## 5  204777_s_at T cells
## 6    206485_at T cellswrite.csv(probesets, file = "probeset.csv")

2. 基因文件

这个同样可以从官网下载，如下：

genes = read.table(curl("http://raw.githubusercontent.com/ebecht/MCPcounter/master/Signatures/genes.txt"),sep = "\t", stringsAsFactors = FALSE, header = TRUE, colClasses = "character",check.names = FALSE)head(genes)
##   HUGO symbols Cell population  ENTREZID      ENSEMBL ID
## 1         CD28         T cells       940 ENSG00000178562
## 2         CD3D         T cells       915 ENSG00000167286
## 3         CD3G         T cells       917 ENSG00000160654
## 4          CD5         T cells       921 ENSG00000110448
## 5          CD6         T cells       923 ENSG00000013725
## 6    CHRM3-AS2         T cells 100506915 ENSG00000233355
write.csv(genes, file = "genes.csv")

3. 基因表达矩阵

我们从GEO上获得表达矩阵之后，可以通过标准化后进行分析，行是样本名称，列是基因名称，如下：

GSE_expr = read.table("E:/代理/吴隆亮/data/GSE32918.rpkm.mat.xls", header = T, row.names = 1,sep = "\t")
GSE_expr = GSE_expr[, 1:100]
head(GSE_expr[, 1:5])
##       GSM814863 GSM814864 GSM814865 GSM814866 GSM814867
## 7A5   10.506953 10.359537 10.344342 10.433970  8.422291
## A1BG   8.783805  8.378983  8.703535  8.355771  8.880121
## A1CF   7.346380  7.385903  7.247202  7.341750  7.534278
## A26A1  7.844565  8.145302  8.114154  7.950246  8.421318
## A26B1  8.212963  7.554923  7.926285  7.861465  8.066015
## A26C3  8.168278  8.167457  8.279593  8.025435 11.741347
dat <- scale(GSE_expr)

参数说明

关于参数起始很简单，就是读入三个文件以及一个筛选类型featuresType，这个参数有三个选择，根据自己的实际实况进行选择即可，细节如下：

expression： matrix or data.frame with features in rows and samples in columns
featuresType： type of identifiers for expression features. Defaults to “affy133P2_probesets” for Affymetrix Human Genome 133 Plus 2.0 probesets. Other options are “HUGO_symbols” (Official gene symbols), “ENTREZ_ID” (Entrez Gene ID)or “ENSEMBL_ID” (ENSEMBL Gene ID)
probesets： Data frame of probesets transcriptomic markers and corresponding cell populations. Fetched from github by a call to read.table by default, but can also be a data.frame
genes：Data frame of genes transcriptomic markers (HUGO symbols or ENTREZ_ID) and corresponding cell populations. Fetched from github by a call to read.table by default, but can also be a data.frame

例子实操

三个文件都准备好后，我们就可以运行结果了得到结果了！

results = MCPcounter.estimate(dat, featuresType = "HUGO_symbols", probesets = probesets,genes = genes)
results[, 1:5]
##                           GSM814863    GSM814864   GSM814865   GSM814866
## T cells                  0.32714921  0.437527214  0.31284051  0.45884846
## CD8 T cells             -0.74049749 -0.765199781 -0.78409775 -0.83310670
## Cytotoxic lymphocytes    0.04693324  0.050832675  0.06533680  0.02455273
## B lineage                0.38729205  0.367077204  0.39510092  0.41235912
## NK cells                -0.63445267 -0.608239131 -0.63877276 -0.59221126
## Monocytic lineage        1.22024519  1.215150052  1.18869549  1.19793489
## Myeloid dendritic cells -0.38268746 -0.449054615 -0.39392798 -0.41873187
## Neutrophils             -0.52048129 -0.444310031 -0.48737525 -0.47471500
## Endothelial cells        0.01874103 -0.003792511  0.02700539  0.03138468
## Fibroblasts              1.18510310  1.185551407  1.20075863  1.19891693
##                           GSM814867
## T cells                  0.30414823
## CD8 T cells             -0.65260954
## Cytotoxic lymphocytes   -0.05112254
## B lineage                1.05286167
## NK cells                -0.69801742
## Monocytic lineage        0.84725859
## Myeloid dendritic cells  0.15522807
## Neutrophils             -0.57250680
## Endothelial cells       -0.23137362
## Fibroblasts              0.49508066
write.csv(results, "MCPcounter.csv")
save(results, file = "MCPcounter.rdata")

绘制热图

heatmap(as.matrix(results), col = colorRampPalette(c("blue", "white", "red"))(100))

免疫细胞差异分析

首先将宽数据转为长数据，再进行绘制箱线图，具体可以参考桓峰基因公众号绘图教程 FigDraw 6. SCI 文章绘图之箱线图 (Boxplot)

library(reshape)
results = as.data.frame(results)
results$Cell = rownames(results)
group = read.table("group.xls", header = T, sep = "\t")
group = group[group$Sample %in% colnames(GSE_expr), ]
cell_long <- merge(group, melt(results, id = "Cell"), by.x = "Sample", by.y = "variable")
library(ggplot2)
library(tibble)
library(ggpubr)
ggplot(cell_long, aes(Cell, value, fill = group)) + geom_boxplot(outlier.shape = 21,color = "black") + theme_bw() + labs(x = "Cell Type", y = "Estimated Proportion") +theme(legend.position = "top") + theme(axis.text.x = element_text(angle = 80,vjust = 0.5)) + stat_compare_means(aes(group = group, label = ..p.signif..),method = "kruskal.test")

结果解读

MCP-counter和cibersort最大的区别就是cibersort是计算淋巴细胞的比例，而MCP-counter是一个决定计数方法。

跟着桓峰基因教程学习，铸造成功的您，记得关注我们哦，如果觉得自己搞生信有难度，快来联系桓峰基因吧！！！

References:

Becht, E., Giraldo, N.A., Lacroix, L. et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol 17, 218 (2016).