RNAseq|免疫浸润也杀疯了，cibersoert？xCELL？ESTIMATE？你常用哪一个

肿瘤免疫浸润分析是一个文献中的网红分析内容，分析软件有很多，本次先介绍一下cibersort ，xCELL 和 ESTIMATE ，这几款软件在文章中的出镜率都很高 。

那都是做免疫浸润的，我们需要了解这么多软件吗？下面简单介绍下使用场景 ，可以根据自己的场景适配分析软件和内容。

（1）Cibersort：基于反卷积算法，可以得到22种免疫细胞的丰度结果。

（2）xCELL：基于标记基因集合的ssGSEA算法，可以得到5个大类的64种细胞类型的丰度结果，含有非免疫的。

（3）ESTIMATE：可以得到整体的基质细胞和免疫细胞评分。

零 准备数据

通过XENA（UCSC Xena）下载TCGA-SKCM的表达量矩阵 ，记得同时下载对应的probeMap文件，方便将Ensembl_ID转为常见的基因symbol。

详细可见 TCGA数据挖掘 | Xena - TCGA数据下载

数据处理

建议使用Rstudio建立project，每个项目一个project 。TCGA | 以项目方式管理代码数据  以及 数据读取存储

在project下新建一个data文件夹，存放XENA下载的数据。读取data文件夹中的数据

library(tidyverse)library(openxlsx)library(estimate)#1）读取表达量数据，，或者给绝对路径data <- read.table("../data/TCGA-SKCM.htseq_fpkm.tsv",sep = "\t" , header = T,#row.names = "Ensembl_ID",                   stringsAsFactors = FALSE ,                   check.names = FALSE)#2）读取probeMap文件，转换Ensembl_IDprobeMap <- read.table("../data/gencode.v22.annotation.gene.probeMap",sep = "\t" , header = T,                       stringsAsFactors = FALSE ,                       check.names = FALSE)# 3）转为基因symbolexp <- data %>%  inner_join(probeMap, by = c("Ensembl_ID" = "id")) %>%  select(gene , starts_with("TCGA") ) #4）相同探针/基因，取表达量最大expr <- exp %>%  mutate(rowMean =rowMeans(.[grep("TCGA", names(.))])) %>%  arrange(desc(rowMean)) %>%  distinct(gene,.keep_all = T) %>%  select(-rowMean) %>%  column_to_rownames(var = "gene")expr[1:4,1:4]
#       TCGA-EE-A2GJ-06A TCGA-EE-A2GI-06A TCGA-WE-A8ZM-06A TCGA-DA-A1IA-06A#MT-CO2         13.45460         12.86170         13.42982         12.99413#MT-CO3         12.67183         11.81769         13.10443         12.62491#MT-ND4         13.21440         12.37920         13.10761         12.44728#MT-CO1         13.14420         12.23376         12.93742         12.69075

一 CIBERSORT 方法

CIBERSORT 是基于支持向量回归（SVR：support vector regression）的原理对人类免疫细胞亚型的表达矩阵进行去卷积的工具。

提供了含有22种免疫细胞亚型的基因表达特征集文件-LM22.txt。

CIBERSORT的整体原理图如下所示

1， 输入数据要求

在原文"Profiling tumor infiltrating immune cells with CIBERSORT "中提到了对于数据的要求

（1）表达量矩阵不能有负值，不能有缺失值 ，不要log转化 ；

（2）RNA-Seq表达量的话，FPKM和TPM都可以；

2，数据转化

CIBERSORT需要2个输入文件，一个是表达矩阵文件，一个是前面提到的LM22.txt（22种免疫细胞的基因表达特征数据）。

因为下载的是log2(FPKM+1)后的FPKM，这里根据CIBERSORT的要求转为FPKM。另外就是需要将行名变成一列 。

#########转为fpkm##########fpkm <- 2^expr - 1#对应列名fpkm <- fpkm %>%   rownames_to_column("gene")#保存成输入文件write.table(fpkm,"SKCM.fpkm.txt",sep = "\t" ,quote = F, row.names=F)

3，CIBERSORT运行

CIBERSORT主要是由doPerm ，CoreAlg 和 CIBERSORT 三个函数构成，可以source的方式直接加载函数。https://rdrr.io/bioc/tidybulk/src/R/cibersort.R

## 加载函数source("cibersort.R")## 指定基因特征集sig_matrix <- "LM22.txt"## 指定表达矩阵mixture_file = 'SKCM.fpkm.txt'
#运行cibersort 此步骤比较慢cibersort_res <- CIBERSORT(sig_matrix, mixture_file, perm=10, QN=TRUE)# output CIBERSORT resultscibersort_res[1:4,1:4]
#                B cells naive B cells memory Plasma cells T cells CD8#TCGA-EE-A2GJ-06A   0.006688073     0.00000000  0.000000000   0.3108573#TCGA-EE-A2GI-06A   0.057898624     0.00000000  0.065030989   0.2741003#TCGA-WE-A8ZM-06A   0.002775345     0.00000000  0.000000000   0.0000000#TCGA-DA-A1IA-06A   0.000000000     0.04971243  0.003565012   0.0000000

这样就得到了每个样本的22种免疫亚型的浸润结果，后续可以进行箱线图，小提琴图，柱形图等的可视化，以及各种分组（分子分型，riskscore高低，临床分期等）的比较。

4，ERROR：出现如下报错，按照提示 不存在哪个包就安装哪个包即可。

二 xCELL 方法

xCELL是基于标记基因集合的ssGSEA算法，可以得到如下图所示的5个大类的64种细胞类型的丰度结果，包括适应性和先天免疫细胞，造血祖细胞，上皮细胞和细胞外基质细胞（含有非免疫的细胞类型）。

1， 输入数据要求

根据xCELL包的github官网https://github.com/dviraran/xCell中提到，注意rowname需要是 gene ，而不是像cibersort中gene需要是首列。

2，运行xCELL

可以按照github使用 rawEnrichmentAnalysis ，transformScores 和 spillOver分步骤运算。

这里可以使用作者包装好的pipeline的方式运行（推荐）。

#devtools::install_github('dviraran/xCell')library(xCell)fpkm2 <- fpkm %>%   column_to_rownames("gene")fpkm2[1:4,1:4]
xCell  = xCellAnalysis(fpkm2,                        rnaseq = TRUE, # RNA-seq                       parallel.sz = 1, #线程数 ，window建议改为1                       parallel.type = "SOCK") dim(xCell)xCell[1:4,1:4]#           TCGA-EE-A2GJ-06A TCGA-EE-A2GI-06A TCGA-WE-A8ZM-06A TCGA-DA-A1IA-06A#aDC            2.673902e-01     2.845002e-01     1.058243e-02     8.379095e-02#Adipocytes     1.363339e-18     3.572483e-18     2.347432e-02     0.000000e+00#Astrocytes     0.000000e+00     0.000000e+00     7.729962e-02     2.331681e-18#B-cells        1.118969e-01     2.076587e-01     3.362450e-19     8.289009e-03

三步运行的方式参照https://github.com/dviraran/xCell

这样就得到了每个样本的64种细胞亚型的结果（含有非免疫），以及ImmuneScore，StromaScore 和 MicroenvironmentScore 。

后续同样可以进行箱线图，小提琴图，柱形图等的可视化，以及各种分组（分子分型，riskscore高低，临床分期等）的比较。

三 ESTIMATE 方法

ESTIMATE可以利用表达谱数据来估计肿瘤样本的基质分数(stromal score )和免疫分数(immune score)，算法使用的是ssGSEA。

1， 载入R包 数据

library(utils)rforge <- "http://r-forge.r-project.org"#install.packages("estimate", repos=rforge, dependencies=TRUE)library(estimate)

2，运行ESTIMATE

filterCommonGenes(input.f="SKCM.fpkm.txt",                    output.f="SKCM.gct",                    id="GeneSymbol")#[1] "Merged dataset includes 10221 genes (191 mismatched)." 
estimateScore(input.ds = "SKCM.gct",               output.ds="SKCM_estimate_score.gct",                platform="illumina") #转录组数据与芯片数据计算过程不同的地方是platform是illumina#[1] "1 gene set: StromalSignature  overlap= 139"#[1] "2 gene set: ImmuneSignature  overlap= 141"
estimate=read.table("SKCM_estimate_score.gct",skip = 2,header = T,stringsAsFactors = F,                  check.names = F)rownames(estimate)=estimate[,1]estimate=t(estimate[,3:ncol(estimate)])head(estimate)
#                 StromalScore ImmuneScore ESTIMATEScore#TCGA.EE.A2GJ.06A   -290.24425  1712.12627      1421.882#TCGA.EE.A2GI.06A   -429.93780  1635.89999      1205.962#TCGA.WE.A8ZM.06A   -320.94643  -804.57717     -1125.524#TCGA.DA.A1IA.06A  -1332.98050   -66.10589     -1399.086#TCGA.D3.A51H.06A     72.56161  3397.46556      3470.027#TCGA.XV.A9VZ.01A   -540.02737  -690.88707     -1230.914

注意：如果是芯片数据， plotform = 'affymetrix'；如果是二代测序数据， platform = 'illumina' 

3，保存数据

保存所有数据，以备后面使用

save(expr,fpkm,cibersort_res,xCell,estimate,file = "RNAseq.SKCM.RData")