主讲老师第四十二课：CIBERSORT 算法
接下来，我们来看看具体每种细胞的丰度或者富集水平如何计算。首先，来看下大家 在 文 章 中 时 常见 到 的 CIBERSORT 分 析 。 关 于 大 名 鼎 鼎 的 CIBERSORT 算 法（CIBERSORT.zip），这估计是目前引用次数最多的免疫细胞浸润估计工具。该算法最早发表在 2015 年 Nature Medicine 杂志上，题目是Robust enumeration of cell subsets from tissue expression profiles。大家记得在使用这些算法时候需要引用这些相应的来源文章简单介绍一下 CIBERSORT 算法的原理。CIBERSORT 算法是基于线性支持向量回归的原理对表达矩阵进行去卷机从而获得不同免疫细胞亚型的工具。对于 CIBERSORT 的算法，作者提供了在线分析和 R 语言
分析两种方法。
关于在线分析，直接进入 CIBERSORT 的官网（https://cibersort.stanford.edu/）。先用教育邮箱注册，在 menu 下拉菜单选择 upload files，上传需要分析的数据，在 CIBERSORT 里称做 mixture file。
但是，受限于教育邮箱的获取和网站上传数据存储空间大小，对于有些小伙伴而言显得并不是十分的友善，每个用户一个是存在 500Mb 的存储空间。因此，我们主要来看一下如何使用代码进行分析下面，我来给大家讲解一下 CIBERSORT 算法的使用方法及其可视化作图
1. 准备工作
在 R 语言中运行 CIBERSORT 共需要三个文件：
（1）官方提供的 22 种细胞基因集“LM22.txt”；（2）自己的表达矩阵；（3）Cibersort 代码。
很奇怪吧，为什么是需要提供代码文件？实际上，关于 Cibersort，作者已经把其运行代码进行了整个包装提供给大家，我们直接把其准备好，引用即可。
1.1 相关依赖性 R 包安装
尽管在 CIBERSORT 算法中，并没有专门的 R 包存在，但是在其注释的 R 函数中需要借助 e1071，parallel 以及 preprocessCore 的 R 包。在正式使用前，先把这几个相关的 R 包先安装好
1.2 函数的下载及准备
其次，我们需要准备一份名为“source.R”的文件，将该文件放置在分析目录下和表达矩阵 放 在 一 起 ，对 于 该 文 件 ， 我 们 可 以 直 接 从 相 应 的 网 站 上 进 行 下 载（https://content.cruk.cam.ac.uk/fmlab/sivakumar2016/Cibersort.R）。
大家把这个文件直接放在工作目录下即可，不需要打开，我们直接引用。
2. CIBERSORT 计算
接着，是基因集文件。关于官方提供的 22 种细胞基因集“LM22.txt”可以直接从文章的SupplementaryTable 1 中进行下载并整理。在 LM22 文件中，共包含了 547 个基因，这 547个基因可以用于区分人类22 种不同的免疫细胞浸润水平，包括 7 种不同的 T 细胞亚型，B细胞，浆细胞等等。关于这几个文件我已经整理好放在压缩包里，不需要再自己进行整理，直接使用即可后期运行唯一需要准备的文件是基因表达矩阵，这三个文件，大家都放到工作目录下，接着，准备相应的文件后，直接运行下面的代码即可其中 perm 表示置换次数，QN 代表结果分位数的归一化。最终，在同一文件夹下可以得到运算结果（"CIBERSORT-Results.txt"），整个的运行结果可能会稍微比较久，因此，我将运行的结果以.RData 的
形式保存下来，后续分析直接读取即可
这里讲一点，关于 log 问题，在 cibersort 中，其自动会检测是否存在 log 转换。因此，基本不需要担心这类问题，但还是建议养成习惯使用 log 转换的数据。注意 Cibersort 结果的默认文件名为 CIBERSORT-Results.txt，在同一文件夹下进行第二次运算会覆盖第一次得到的文件，建议在每一次运算之后对文件重命名。因为这个过程存在置换 100 次的情况，每次运行的结果都会有一定的差异，但基本不影响，特别是对参数 perm 的次数设置高以后，结果会趋于稳定，这个和我们前面介绍的 GSEA 分析是类似的。如果大家有时间的话，也可以把其设置为 1000，此时的结果基本已经很稳定了（这一步运算，次数越多，样本越多，会越慢）。
3.CIBERSORT 分析结果的可视化展示
接下来，我们以 CIBERSORT 分析结果为例，讲解一下免疫细胞分析结果的可视化。在 分析前，首先通过 load() 函数将前面得到的 CIBERSORT 结果读取进来：
在读取进来后，我们看到其一共有 25 列，但是，最后三列并不是相应的免疫细胞类型，而是 P 值，相关性系数等信息。因此，在分析前需要将这三列删除，取前 22 列作为细胞丰度。
同时，我们需要进一步去除了丰度全为 0 的免疫细胞亚型：
最终，一共得到 21 种免疫细胞结果。这里说明一下这个 P 值，因为它是和最标准样品进行比对得到的结果，因此计算得到的 P 值小于 0.05 的会很少，一般来说，可以不管，直接所有样品入组分析
3.1 barplot 图
接下来，进入可视化的环节。首先，自然是展示一下整个项目中患者整体的免疫细胞浸润情况。因为CIBERSORT 计算得到的是相对丰度，所以每个样品中所有细胞相加之和为 1通过 rainbow()函数获取彩虹颜色，结合 ncol(ciber.res)获得免疫细胞数量，最终创建彩虹色板，其中参数 0.7 表示颜色的透明度。接着，我们使用基础包来进行绘制柱状分布图。
设置好一块画板。接着因为柱状图，使用 barplot()函数来进行绘制：
这样，一个基本的形状就完成了。接着，我们需要把左边的 y 轴补上去，使用 axis()函数来添加轴：
这样，图的整体部分已经完成。我们需要在边上添加图注，表明每种颜色所代表的细胞类型。
最终得到整个结果。后期，大家在使用时，准备好了 cibersort 的运行结果后，直接套用这部分代码即可，除了 legend 部分需要根据自己的结果进行微调外，其他基本不需要修改。
最后，介绍一个画板相关的函数 dev.off()，他的作用是关闭当前的画板。如果画板不关闭，后期的做图会继续添加，或者出现报错。
这个图代表的是每一个列代表一个病人，所有细胞相加都等于 1。它展示的每个患者中 这 21 种免疫细胞的占比情况，经过 cibersort 计算得到的结果，展示的是各个免疫细胞在一个样本中的占有率，因此这部分结果不存在后期的均一化等等都操作。本来 x 轴显示患者 id的，但是 500 来个密密麻麻的 id，已经没有展示的必要了。一般来说，这个图展示在正文里面，然后把数据作为补充材料上传，这样就可以展示所有患者的整体免疫浸润水平。
3.2 相关性热图
相关性热图在文章里面也进行了展示，首先，绘图之前，我们需要计算一下各个免疫细胞之间的相关系数直接使用 cor()函数即可快速计算相关系数，接着，通过 round()函数保留两位小数：
下面，介绍两种不同相关性热图的绘制方法
绘制方法一：使用 pheatmap 包
直接通过调用 pheatmap::pheatmap()函数，即可完成相关性矩阵的热图
绘制方法二：使用 corrplot 包
关于这两个 R 包，都在昨天进行了介绍，大家可以再复习一下相关参数的使用。
3.3 PCA 图
随后，我们来进一步看一下示例文章中的主变量分子 BTK 的高和低表达是否会对免疫细胞浸润的水平产生差异，这时，就需要使用主成分分析看一下两组整体的免疫水平。
加载完所需要的 R 包后，就是输入数据的准备。对于输入数据，主要需要两个层面，第一个是免疫细胞表达数据 exp，第二个则是分组信息 group_list。
对于表达数据，把其转换为数据框格式即可。接着是分组数据的准备过程，我们首先需要读取 TCGA-LUAD 项目的 mRNA 表达数据，提取其中目标基因 BTK 的表达情况提取表达，并进行对数转换。
把对应病人的 id 号赋值给它。这里插播一下，向量是可以有 name 的，直接 name(向 量)就可以赋值。
但是，哪怕有了 name，它还是一维数据，和数据库、矩阵不一样。这时，我们可以使用 identical()函数来判断一下，两组数据病人的 id 号和顺序是否相同无论是名字还是顺序，两者是完全相同的，这样就可以放心的进行后面的分组。基于BTK 基因的表达情
况，根据中位值将其分成高（High）和低（Low）两组：
并且将分组变量 group_list 定义为因子格式。
通过 table()函数，快速查看两个分组中患者的数量随后，将准备好的两组数据直接输入 PCA 分析中，通过内部一系列的降维分析，最终呈现分析结果。
可以看到的是，两组之间实际上免疫水平有差异，但是差异很小，所以这也可能是为什么作者没有展示这部分结果的原因。
3.4 免疫细胞相关性分析
除了基因的表达与整体免疫特征之间的相关性分析外，文章中进一步分别计算了 BTK基因与不同免疫浸润细胞之间的相关系数以及相关性图。在此，我们使用 ggpubr 包进行绘制。
同时，我们设置了相关性分析中检验 p 值的阈值，以 P <0.05 作为显著性分析的分界线：
接下来，就是最主要的相关性分析环节了。首先，设置一个空的变量来储存分析结果：
接着对循环进行设置，对于 i 逐个取值每个免疫细胞。
然后，提取该免疫细胞的浸润水平和基因的表达值，形成一个新的数据框接着，计算 BTK 基因的表达和免疫细胞之间的相关性然后把每次的运行结果都添加到 outTab 中进行储存，便于最后整体输出随后，我们再进行一个判断，将 p 值与前面设置的阈值 0.05 相比较，如果小于 0.05，那么接下来绘制相关性图，使用 ggpubr 包的 ggscatter()函数绘制散点图。最后，将得到的图形，使用 ggsave()函数来快速保存同时，对应着把相关性的分析结果，以.txt 的格式也一起保存下来，这样，整个的数据和图就都得到了。
3.5 小提琴图
最后，我们来看下如何绘制免疫细胞分组的小提琴图
3.5.1 加载 R 包
在此，我们使用 vioplot 包进行小提琴图的绘制。
3.5.2 加载表达谱
由于涉及分组比较，首先需要将表达矩阵读取进来，并设置相应的目标基因
3.5.3 取出目标基因表达
接着，和前面的步骤一样，我们需要提取目标基因的表达信息，并对其进行命名接着，以中位数进行分组，对高于中位数和低于中位数的样本名称进行提取。
同时，使用 length()函数分别计算两组的数量。
结果显示，其中高表达组 250 个，低表达组 251 个
3.5.4 加载 CIBERSOciber 免疫丰度数据
随后，读取 cibersort 分析结果，并参照前面的方法去除丰度全为 0 的细胞，同时，根据高低样本对表达谱重新排序。
3.5.5 绘制小提琴图
接下来，我们开始小提琴图的绘制过程，首先，把整体的轮廓画出来。
接下来，使用 legend() 函数绘制在右上角绘制 high 和 low 分组的标签接着，设置循环函数，对每一个丰度不为 0 的细胞绘制小提琴图首先，提取相应的免疫细胞在两组中的表达情况。
随后，分别绘制相应的小提琴图，其中低表达样本为蓝色，高表达样本为红色最后，计算两组之间的显著性差异，并将 P 值添加到相应的位置上。
这样，关于 CIBERSORT 的常见可视化方法就介绍到这里了