本文首发于“生信补给站”公众号 https://mp.weixin.qq.com/s/Seo6MZicBgBOfhz7equ80g
肿瘤免疫浸润分析是一个文献中的网红分析内容，分析软件有很多，本次先介绍一下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_ID probeMap <- read.table("../data/gencode.v22.annotation.gene.probeMap",sep = "\t" , header = T, stringsAsFactors = FALSE , check.names = FALSE) # 3）转为基因symbol exp <- 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 results cibersort_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")
