CIBERSORT肿瘤免疫微环境分析，一文就搞定
免疫治疗是如何工作的呢？其主要通过激活人体本身的免疫系统，依靠自身的免疫机能杀灭癌细胞。目前已在多种肿瘤如黑色素瘤，非小细胞肺癌、肾癌和前列腺癌等实体瘤的治疗中展示出了强大的抗肿瘤活性，多个肿瘤免疫治疗药物已经获得美国FDA批准临床应用。
随着组学技术的进步，我们也逐渐了解到肿瘤微环境的免疫环境的复杂性和多样性以及它对免疫治疗的重要影响。通过进一步分析和了解肿瘤免疫微环境将有助于免疫治疗反应性的改善。因此很多人就开始研究肿瘤的免疫微环境，试图解析肿瘤组织中的免疫细胞构成。
什么是免疫浸润
Background
所谓肿瘤免疫微环境的概念，就是肿瘤的内部和周围往往聚集着大量的免疫细胞。这些免疫细胞之间，以及肿瘤细胞和免疫细胞之间存在着千丝万缕，免疫细胞多种多样，所以所谓免疫微环境，或者说免疫浸润的分析，本质上，就是搞清楚肿瘤组织当中免疫细胞的构成比例。
如何分析免疫浸润
Method
讲完了背景，我们就来看看，到底是如何判断免疫细胞的构成的呢？目前看来，计算免疫微环境的方法，常用的有两类：
（1）一类就是“实测法”，高精度的单细胞测序，Single cell RNA-seq。为了搞清楚构成，索性一个一个细胞去分开测序，把细胞表达都测了，通过一些标记基因，来判断到底这里面有多少类型的免疫细胞，以及每种细胞的构成比例。
（2）还有一种方法，那就是“推测法”，通过bulk RNA-seq进行推测。我们常规的普通RNA测序和芯片结果就是把所有细胞统一裂解抽提RNA，那么检测到的基因表达中，就是各种免疫细胞和肿瘤细胞混杂在一起的表达，通过一些算法，我们可以从这个混杂的表达谱中推断免疫细胞的构成比，常用的软件有CIBETSORT、TIMER、EPIC等等。
这两种方法各有有缺，单细胞最大的缺点就是贵！贵在一个样本要数万元的经费；其次是难！难在有些组织不好处理，单细胞比较难以制备。不过单细胞测序最大的优点就是信息量大，结果更接近事实情况，讲的现实一点，比较好发高分文章。Bulk RNA-seq缺点就是不准确，优点就是便宜，公共数据多，比如TCGA中大量的表达谱数据，任人取用。单细胞的分析，我们以后再说，今天就来讲讲比较亲民的bulk sequencing如何进行免疫浸润的分析。
分析软件
Software
刚才我们也提到了，分析免疫浸润的软件有很多，十几种之多，今天我们讲讲里面比较常用的一个：CIBERSORT。
为什么，这么多软件我们单独挑了CIBERSORT来讲解呢？原因就在于这是一款非常经典的软件，CIBERSORT是2015年发表于Nature Methods上的文章（IF：28.47）。
CIBERSORT 是基于线性支持向量回归（linear support vector regression）的原理对人类免疫细胞亚型的表达矩阵进行去卷积的一个R/网页版工具。多用于芯片表达矩阵，对未知混合物和含有相近的细胞类型的表达矩阵的去卷积分析优于其他方法 (LLSR,LLSR,PERT,RLR,MMAD,DSA) 。该方法是基于已知参考集，提供了22种免疫细胞亚型的基因表达特征集：LM22.
为了方便用户的使用，CIBERSORT被开发成一个网页版的工具，其网址链接：http://cibersort.stanford.edu/.不过，由于CIBERSORT初期版本不支持单细胞测序数据的分析，所以近期又升级到了CIBERSORTX，升级版的网址链接为：https://cibersortx.stanford.edu/
下面我们就来讲讲CIBERSORT的使用方法：
（1）网页版
用户只需要注册一个账号，就可获得500M存储数据和结果的空间。操作时，只需上传标准的表达矩阵文件即可分析免疫浸润；如果想要分析包含其他细胞类型的浸润比例，则需按照官网提示的格式上传相应的文件。步骤如下：
- 在网页注册账号
- 准备表达数据
- 需要分析的表达矩阵
- 参考数据集
- 上传数据，设置参数，运行CIBERSORT
（2）R语言版
看上去比网页版要求要高一些，最起码你要知道R语言怎么用。不过R语言的CIBERSORT比网页版用上去方便很多，最起码没有数据量的限制，适合大数据运算，也直接运行，非常方便的。基本步骤如下：
- 准备数据
- 需要分析的表达矩阵
- 参考数据集
- 一键运行
可以说，CIBERSORT作者把代码做的非常方便易用，只需一行命令即可运行，下面就进入到我们的实战演练环节，给大家演示以下如何使用CIBERSORT进行免疫浸润的分析。
示例操练
Implementation and Results
我们选取的测试数据呢，是TCGA数据库当中LUAD肺腺癌的数据，我们的课题是，比较不同给肿瘤分期G1、G2、G3、G4期的肺腺癌样本中，免疫微环境的差异，探寻宿主免疫在肿瘤进展过程中的作用。先来看一下我们的测试数据：
当然这个截图只是我们数据中很小的一部分，我们的数据有300多行，18000多列，分别对应300多个患者，18000多个基因。（注：测试数据可以联系客服领取）
下面就是进行分析了，在进行分析的时候我们需要准备三个文件：
- LM22.txt（可以从CIBERSORT网站下载，这个就是22种免疫细胞的参考marker基因表达）
- CIBERSORT.R（CIBERSORT源代码，从官网下载）
- .txt（基因表达谱文件，就是上面说的肺腺癌表达谱）
然后，运行CIBERSORT非常简单，三个文件放到一个文件夹中，在RStudio中设置工作路径，运行如下代码即可：
source("CIBERSORT.R")# Define LM22 fileLM22.file <- "LM22.txt"exp.file <- "datExp.processed.txt"TME.results = CIBERSORT(LM22.file, exp.file, perm = 1000, QN = TRUE)# output CIBERSORT resultswrite.table(TME.results, "TME.results.output.txt", sep = "\t", row.names = T, col.names = T, quote = F)
这个代码看上去是不是很简单呢？没错，就是这么方便。那么，输出的结果长什么样子呢？我们来看一下：
每一行一个样本，每一列一种细胞，总共有22种细胞，这里的数值代表的是免疫细胞所占的比例，比如CD8+ T Cell在第一个肿瘤样本中是0.2282，那就代表着在该样本中CD8+ T Cell占总的免疫细胞的22.82%，所以，CIBERSORT输出的是一个比例，22种免疫细胞的比例加起来就等于1。好了，到这里为止，我们获取了300多个肺腺癌样本中的免疫浸润结果，分析到此就可以结束了。
高级可视化
Advanced plots
作为对粉丝异常负责的科研猫团队，仅仅做到上面这一步肯定是不够的，我们还给大家附送了高级可视化代码，用来比较不同分组之间的免疫细胞比例差异，包括柱状图、箱线图等。近期在某友商的网站看到类似的分析代码卖到1000元高价，而我们则选择继续默默地免费分享。
这里的高级绘图，大家可以参见我们的科研绘图系列教程（更加完整）。这部分绘图的内容，不属于免疫浸润的分析范围，属于延伸出来内容，主要是为了方便大家把分析结果呈现到文章中去，所以对做图不感兴趣的同学，当然也可以跳过不看。
首先，我们用箱线图表示22种免疫细胞在所有样本种的比例，在这里看到最高的应该是Naive B cells。
# boxplotggboxplot( plot.info, x = "CellType", y = "Composition", color = "black", fill = "CellType", xlab = "", ylab = "Cell composition", main = "TME Cell composition") + theme_base() + theme(axis.text.x = element_text( angle = 90, hjust = 1, vjust = 1 ))
然后，为了比较不同分期之间的免疫细胞的差异，我们使用批量统计比较的箱线图。当然，在这里，不用我说大家都知道，为了比较不同分期之间的差异，我们还需要读入一个样本分期的表型文件。这里我就不做演示了，只把做图的代码放出来。
# box plot group by GRADEggboxplot( plot.info, x = "CellType", y = "Composition", color = "black", fill = "GRADE", xlab = "", ylab = "Cell composition", main = "TME Cell composition group by GRADE") + stat_compare_means( label = "p.signif", method = "t.test", ref.group = ".all.", hide.ns = T ) + theme_base() + theme(axis.text.x = element_text( angle = 90, hjust = 1, vjust = 1 ))
同样的，为了更加直观地感受不同细胞在各个样本中的比例差异，在这里非常适合做一个堆积的比例柱状图，为啥？因为各种细胞比例加起来等于1呀，这个我们刚才讲过了。
# Barplot of cell componment of each sample using hclustsample.index <- hclust(dist(TME.data[, TME.cells]), method = "ward.D")$ordersample.order <- TME.data$PATIENT_ID[sample.index]ggbarplot( plot.info, x = "PATIENT_ID", y = "Composition", size = 0, fill = "CellType", color = "CellType", order = sample.order) + theme_base() + theme( axis.text.x = element_text( angle = 90, hjust = 1, vjust = 1, size = 1 ), legend.position = "bottom" )
好了，以上是我们给大家演示的三种图形，当然还有很多其他的图形可以用于免疫浸润结果的可视化，这些大家在我们的科研绘图板块中学习了这么久，应该可以自行发挥啦。