肿瘤免疫浸润分析是一个文献中的网红分析内容,分析软件有很多,本次先介绍一下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的整体原理图如下所示
在原文"Profiling tumor infiltrating immune cells with CIBERSORT "中提到了对于数据的要求
(1)表达量矩阵不能有负值,不能有缺失值 ,不要log转化 ;
(2)RNA-Seq表达量的话,FPKM和TPM都可以;
因为下载的是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)
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种细胞类型的丰度结果,包括适应性和先天免疫细胞,造血祖细胞,上皮细胞和细胞外基质细胞(含有非免疫的细胞类型)。
根据xCELL包的github官网https://github.com/dviraran/xCell中提到,注意rowname需要是 gene ,而不是像cibersort中gene需要是首列。
可以按照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。
library(utils)
rforge <- "http://r-forge.r-project.org"
#install.packages("estimate", repos=rforge, dependencies=TRUE)
library(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")
以上就完成了几种常见的免疫浸润的分析方法