定量免疫浸润在单细胞研究中的应用

最近听说定量免疫浸润很火，于是我也报名参加了果子老师的课程，跑了几个R包，你比如说xCell啊，GSVA啊，MCPcounter啊，ConsensusTME，ImmuneSubtypeClassifier，当然还有大名鼎鼎的CIBERSORT。安装R包跑实例文件这件事是我经常做的，但是很少做的这么系统。相比于跑R包，对我的更大难点在于理解定量免疫浸润这件事。

免疫浸润工具箱

也有国产在线的：http://bioinfo.life.hust.edu.cn/ImmuCellAI#!/

定量免疫浸润在说一件什么事呢？

从平均人说起

平均人（average man）以人体测量结果平均值为基础建立的人体模型。其身体结构的各个尺寸均与特定总体的平均值相对应，可反映特定群体身体尺寸的平均状况。

我们知道，人有性别/年龄/身高/胖瘦/嗓门大小。。。。等等属性，每个人不同，一般来看东北人和广东人还是有区别的，虽然我们不需要一个一个地比较。

平均人就好比传统的bulk RNA测序，拿到一块组织去测序，得到一个样本各个基因的表达情况。所以我们拿到bulk RNA的结果，其准确度相当于我们对歪果仁的模糊印象。

定量免疫浸润说的是我们凭着这个模糊印象，可以大致推断出歪果仁某几种特殊类型的比例（或在某方面的分值）。这个比例或者分值，在一定程度上也可以反应该地方人口的异质性：你看，本来只是一个平均值，但是里面却蕴含着可爱的异质性。

定量免疫浸润（Immune infiltration）

首先，“浸润”是啥意思？

严阵《牡丹园记》：“这蒙蒙的绿意，这团团的红雾，真像刚滴到宣纸上的水彩一样，慢慢地浸润开来。”

这里面的“浸润”很好的诠释了其“逐渐渗透，引申为积久而发生作用”的意义。所以“浸润”有逐渐渗透的意义。

我们提取临床癌症组织去测序(bulk的)，并不是纯肿瘤的结果（肿瘤有纯的吗？来一打），在测序结果中我们会发现，一些属于血管的基因，属于免疫系统的基因也有表达。这时候我们就知道，测的是个混合体，有肿瘤细胞，有免疫细胞，有血管，有细胞外基质。

定量免疫浸润通过简单的基因表达矩阵数据，将免疫细胞类型含量计算出来。随着算法的发展，到后来也不仅是免疫细胞了，可以是很多细胞类型。

好吧，我已经被浸润了

基本的套路是计算每个样本在不同细胞类型中的分值，主要的算法分为两类，包括:

基于Marker gene的算法, 如GSE(V)A
基于 deconvolution的算法，如CIBERSORT

这俩算法前者大家尽管也不懂，但是没有那么吓人，第二个反卷积(deconvolution)真的是听着就头大。那么，到底什么是反卷积？

在理解这个反卷积之前我们来看看在数学上卷积是什么：

卷积就像加减乘除一样是一种运算，只是低等数学用不到，就没教，它是指：把二元函数 U(x,y) = f(x)g(y) 卷成一元函数 V(t) ，俗称降维打击。

在我们的bulk RNA的例子中，就是把平均值看作单个细胞表达值与某个函数（和检测方法有关）的卷积。获得样本bulk RNA 表达量的过程，就是一个卷积的过程，而反卷积，是卷积运算的逆运算。

鉴于我们bulk RNA 的数据已经做的很全了，TCGA数据库里面有大量的已经卷积好的表达量，我们可以用反卷积的方法，看看样本中到底免疫细胞有着怎样的比重（分布）。而这，也这不也是在说明样本的异质性吗？

比如一个定量免疫浸润的套路是这样的：

公共数据库下载数据
基于表达谱做定量免疫浸润（如CIBERSORT）
得到样本细胞的异质性
统计（富集分析）
验证

在众多工具中CIBERSORT的应用是最为广泛的，所以这篇文献是有必要打印出来慢慢评鉴的。在这篇文献中，我们也看到了一种把bulk RNA和scRNA结合起来的桥梁：

Determining Cell Type Abundance and Expression From Bulk Tissues With Digital Cytometry

CIBERSORT函数在果子老师的课上讲的很详细，这里仅做简单的复现，看看CIBERSORT究竟做了什么。更多介绍参见官网：https://cibersort.stanford.edu/

source("CIBERSORT/CIBERSORT.R") 
.libPaths("C:\\Users\\86158\\Documents\\R\\win-library\\3.6")

读入自带的 signature matrix（细胞类型打分矩阵），并看看他的格式如何：

sig_matrix <- "CIBERSORT/LM22.txt"
sig_matrix
sdf<- readr::read_tsv(sig_matrix)

 sdf[1:4,1:4]
# A tibble: 4 x 4
  `Gene symbol` `B cells naive` `B cells memory` `Plasma cells`
  <chr>                   <dbl>            <dbl>          <dbl>
1 ABCB4                   556.              10.7           7.23
2 ABCB9                    15.6             22.1         653.  
3 ACAP1                   215.             322.           38.6 
4 ACHE                     15.1             16.6          22.1 

colnames(sdf)
 [1] "Gene symbol"                  "B cells naive"                "B cells memory"               "Plasma cells"                
 [5] "T cells CD8"                  "T cells CD4 naive"            "T cells CD4 memory resting"   "T cells CD4 memory activated"
 [9] "T cells follicular helper"    "T cells regulatory (Tregs)"   "T cells gamma delta"          "NK cells resting"            
[13] "NK cells activated"           "Monocytes"                    "Macrophages M0"               "Macrophages M1"              
[17] "Macrophages M2"               "Dendritic cells resting"      "Dendritic cells activated"    "Mast cells resting"          
[21] "Mast cells activated"         "Eosinophils"                  "Neutrophils"

也就是有22种细胞类型。其实我们是可以根据自己的样本类型来利用scrna的数据来制作这个矩阵的。制作这个有什么用？下面再说。

读入表达谱数据。在这里了停一下，三秒禅，应该是什么格式的数据，count？TPM？

CIBERSORT接受的数据是不需要log的，如果你一不小心取了log，也不要害怕，他会帮你再去掉。

 mixture_file = 'exprMat.txt'
exp<- read.table(mixture_file,header = T)

exp[1:4,1:4]
  rownames.exprMat. LAU125 LAU355 LAU1255
1              A1BG   0.82   0.58    0.81
2              A1CF   0.00   0.01    0.00
3               A2M 247.15  24.88 2307.94
4           A2M-AS1   1.38   0.20    2.60

我们来执行核心函数CIBERSORT：

res_cibersort <- CIBERSORT(sig_matrix, mixture_file, perm=10, QN=TRUE)

        B cells naive B cells memory Plasma cells T cells CD8 T cells CD4 naive T cells CD4 memory resting T cells CD4 memory activated
LAU125     0.01146448     0.33118622  0.000000000  0.03587147                 0                  0.1465500                   0.00000000
LAU355     0.18558675     0.36063872  0.005305268  0.06538556                 0                  0.1969691                   0.07829503
LAU1255    0.04230761     0.07319686  0.014003042  0.37666902                 0                  0.0000000                   0.14757024
LAU1314    0.24241384     0.34369856  0.015805166  0.06588054                 0                  0.1580854                   0.05774793
        T cells follicular helper T cells regulatory (Tregs) T cells gamma delta NK cells resting NK cells activated  Monocytes
LAU125                 0.09940070                          0          0.00000000       0.04715444         0.00000000 0.09273837
LAU355                 0.02004924                          0          0.02277827       0.00000000         0.00000000 0.01053398
LAU1255                0.11192989                          0          0.04067900       0.00000000         0.04482746 0.02217635
LAU1314                0.04853605                          0          0.03254491       0.00000000         0.01240244 0.00000000
        Macrophages M0 Macrophages M1 Macrophages M2 Dendritic cells resting Dendritic cells activated Mast cells resting
LAU125      0.03315481     0.00000000     0.12737478              0.00000000               0.004796157        0.000000000
LAU355      0.00000000     0.00000000     0.04398640              0.00000000               0.010471658        0.000000000
LAU1255     0.02427582     0.02090543     0.06168147              0.01180571               0.000000000        0.002012619
LAU1314     0.00000000     0.00000000     0.01644681              0.00000000               0.006438391        0.000000000
        Mast cells activated Eosinophils Neutrophils P-value Correlation      RMSE
LAU125            0.06608013 0.004228426 0.000000000     0.5  0.01785367 1.0921730
LAU355            0.00000000 0.000000000 0.000000000     0.0  0.66799046 0.7449124
LAU1255           0.00000000 0.004148568 0.001810919     0.0  0.20761899 1.0207522
LAU1314           0.00000000 0.000000000 0.000000000     0.0  0.64988147 0.7595212

可以看出22种细胞类型在每个样本中的分布及其显著性指标：P-value ，Correlation RMSE。

我们简单探索一下这个返回的结果：

library(pheatmap)
pheatmap(res_cibersort[,1:22])

 apply(res_cibersort[,1:22] ,1,sum)
 LAU125  LAU355 LAU1255 LAU1314 
      1       1       1       1

可见它是按满分是1 来对样本打分的，也就是每个样本只有这22种细胞类型。既然如此，我们就可以：

library(ggplot2)
library(tidyverse)
allcolour=c("#DC143C","#0000FF","#20B2AA","#FFA500","#9370DB","#98FB98","#F08080","#1E90FF","#7CFC00","#FFFF00","#808000","#FF00FF","#FA8072","#7B68EE","#9400D3","#800080","#A0522D","#D2B48C","#D2691E","#87CEEB","#40E0D0","#5F9EA0","#FF1493","#0000CD","#008B8B","#FFE4B5","#8A2BE2","#228B22","#E9967A","#4682B4","#32CD32","#F0E68C","#FFFFE0","#EE82EE","#FF6347","#6A5ACD","#9932CC","#8B008B","#8B4513","#DEB887")
res_cibersort[,1:22] %>% reshape2::melt() %>%
  ggplot(aes(x=Var1,y=value,fill=Var2))+
  geom_bar(position = 'stack',stat="identity")+
  scale_fill_manual(values =allcolour )  + theme_bw()

其他工具的思路：

我们试一下这个模型在单细胞中的应用情况。在这之前我们需要有一个已经注释好的数据，以便我们看看细胞类型能不能对应上。于是，我们请出pbmc3k.final数据集，这个是可以安装使用的。

library(Seurat)
library(SeuratData)
head(pbmc3k.final@meta.data)
              orig.ident nCount_RNA nFeature_RNA seurat_annotations percent.mt RNA_snn_res.0.5 seurat_clusters
AAACATACAACCAC     pbmc3k       2419          779       Memory CD4 T  3.0177759               1               1
AAACATTGAGCTAC     pbmc3k       4903         1352                  B  3.7935958               3               3
AAACATTGATCAGC     pbmc3k       3147         1129       Memory CD4 T  0.8897363               1               1
AAACCGTGCTTCCG     pbmc3k       2639          960         CD14+ Mono  1.7430845               2               2
AAACCGTGTATGCG     pbmc3k        980          521                 NK  1.2244898               6               6
AAACGCACTGGTAC     pbmc3k       2163          781       Memory CD4 T  1.6643551               1               1

我们看到是有seurat_annotations 的细胞类型注释的。

为了模拟bulk RNA的数据我们对亚群取平均值：

avexpr <- AverageExpression(pbmc3k.final)

因为CIBERSORT直接读电脑文件，而我们又不想先把数据输出，于是改一下函数的读数据方式，再source。

source("CIBERSORT/CIBERSORT.R") 
cibe<-CIBERSORT(sig_matrix, cbind(rownames(avexpr$RNA), avexpr$RNA), perm=10, QN=TRUE)
cibe[1:4,1:4]
             B cells naive B cells memory Plasma cells T cells CD8
Naive CD4 T   0.0004745028    0.006527557   0.00000000  0.40124293
Memory CD4 T  0.0000000000    0.017245116   0.00000000  0.14626803
CD14+ Mono    0.0057763350    0.007812051   0.00000000  0.01540697
B             0.3429285094    0.450309369   0.03076117  0.00895134

同上，我们来用热图看看对应关系怎么样：

library(pheatmap)
pheatmap(cibe[,1:22])

虽然我们用的并不是肿瘤样本，也虽然热图的大部分都是天蓝色，但是橘黄色的色块告诉我们：对应关系还算靠谱。

cibe[,1:22] %>% reshape2::melt() %>%
  ggplot(aes(x=Var1,y=value,fill=Var2))+
  geom_bar(position = 'stack',stat="identity")+
  scale_fill_manual(values =allcolour )  + theme_bw() +theme(axis.text.x = element_text(angle = 90, hjust=1))

这个柱形图显然会得到令人困惑的结果，比如NK细胞里面为什么还会有其他21种细胞呢？这当然和pbmc3k.final做细胞注释的方法有关，也和CIBERSORT这个算法有关，他总会给一个表达谱再这22种细胞类型中找比例关系，除非自己准备一个 signature matrix。但是这些结果并不是一无用处，特别是在我们对细胞类型一无所知的时候，还是有一些参考价值的。

xCell，GSVA，MCPcounter，ConsensusTME，ImmuneSubtypeClassifier等等，我就不再一一演示了，基本上是简单的基因表达矩阵数据，将免疫细胞类型含量计算出来。不同的软件默认给定的细胞类型不同，算法有不同，但是结果都是得到每个样本在不同细胞类型中的分值。

在单细胞中的应用

基于以上的讨论，我们知道了所谓的定量免疫浸润，实在无法获得单细胞水平异质性的历史条件下，一种通过bulk RNA数据来推断(肿瘤)样本免疫细胞异质性的手段。

我们知道了做定量免疫浸润的关键是有表达谱和一个参数数据集，表达谱在我们单细胞这里是不缺的，关键在于参数数据集如何获得。定量免疫浸润的这套方法，在单细胞中至少有以下应用方向：

验证

单细胞数据科学是一个超指数发展的行业，数据是爆发式的增长的，但是数据之间并不是孤立的，之前有着丰富的数据积累，其中TCGA，GEO是比较详细的。越来越多的文章，开始利用公用数据库来挖掘有用的信息来丰富自己的研究，甚至专门有个数据挖掘的方向。在这样的背景下，免疫浸润也是揭示异质性的，那么自然地，我们会联想到：拿我们的单细胞数据通过免疫浸润与公用数据联系起来，相互验证。

细胞类型推断

既然定量免疫浸润可以得到每个样本在不同细胞类型中的分值，我们可不可以根据这个分值来推断手里表达谱的细胞类型呢？在学理上应该是可以的，反卷积的算法也许有些不当，但是如果把求平均值的过程看作卷积，也不是不可以吧。基于GSVA 的方法本质就是自定义基因集做富集分析，所以完全是可以的，关键在与marker gene 的选择。

dtangle:基于反卷积的表达谱分解确定细胞组分
 使用xcell根据表达谱推断样本组成细胞的类型
 3到11分文章解读(肿瘤免疫浸润挖掘方向)
必读|TCGA数据挖掘-肺癌肿瘤免疫浸润分析
 一文献一技术路线：再来“免疫浸润”
肿瘤免疫浸润分析工具汇总-数据挖掘新高度
 TIMER:肿瘤浸润免疫细胞分析的综合网站
 如何用转录组数据定量肿瘤浸润免疫细胞
 视频课程：TCGA数据免疫浸润的量化方法
 癌细胞浸润是什么意思？
Targeting adenosine receptor 2B in triple negative breast cancer
Identification of key biomarkers and immune infiltration in the synovial tissue of osteoarthritis by bioinformatics analysis
Profiles of Immune Infiltration and Prognostic Immunoscore in Lung Adenocarcinoma
Elisa Martini， et al.Single-Cell Sequencing of Mouse Heart Immune Infiltrate in Pressure Overload–Driven Heart Failure Reveals Extent of Immune Activation.Circulation. 2019；140：2089–2107
Single-cell immune landscape of human atherosclerotic plaques
如何通俗易懂地解释卷积？
Understanding-Convolutions/
量化免疫浸润时CIBERSORT的注意事项。
evaluation-of-methods-to-assign-cell-type-labels-to-cell-clusters-from-single-cell-rna-sequencing-data

宝藏图