量化免疫浸润时CIBERSORT的注意事项

CIBERSORT的注意事项

CIBERSORT的使用，我们举了几个例子，但是还有一些疑问，只要有疑问在，就感觉心里不踏实。
比如：
1.CIBERSORT需要什么样的数据？
是芯片数据还是测序数据，是什么格式的，需要取log么？
2.CIBERSORT输出的数据是表示什么？
在进过一系列运算后我们会得到一个矩阵

行是样本，总共25列，前面22列代表22个免疫细胞。那么表格内的数值是什么意思？
3.上面的表格中，最后还有三列

这里的p值，相关性系数怎么算出来的？RMSE是什么啊？

带着这些问题，我们就探索一下吧。

CIBERSORT究竟在干什么？

我们能够对CIBERSORT进行探索，在于他除了有网页版本之外，还提供了一个脚本。
这个脚本在登陆官网后可以申请获得。
我们先来展示一下，在R语言里面如何运行CIBERSORT
你需要三个东西：
第一，这个脚本。
第二，基因特征文本
第三，表达量矩阵。

前两个都可以在官网下载。
最后一个是我们自己的数据，行是基因，列是样本。

然后运行就可以了

## 加载函数
source("CIBERSORT/CIBERSORT.R") 
## 指定基因特征集
sig_matrix <- "CIBERSORT/LM22.txt"
## 指定表达矩阵
mixture_file = 'exprMat.txt'
## 运行
res_cibersort <- CIBERSORT(sig_matrix, mixture_file, perm=10, QN=TRUE)

其中nperm给的是置换的次数，QN如果是芯片设置为T，如果是测序就设置为F，测序数据最好是TPM
因为我对这个小函数进行了一点改造，所以会展示计算过程，就是缓解一些等待的焦虑

最后就得到了我们之前说的结果。就这么简单。
当我们运行第一句的时候

source("CIBERSORT/CIBERSORT.R") 

环境中出现了三个函数

CIBERSORT是主函数，CoreAlg用来进行SVM支持向量机分析，doPerm置换检验，产生背景r值。

我们可以打开CIBERSORT函数看看

CIBERSORT <- function(sig_matrix, mixture_file, perm=0, QN=TRUE){
  print("Modified by Guozi")
  print(paste0("Started at:",Sys.time()))
  #read in data
  X <- read.table(sig_matrix,header=T,sep="t",row.names=1,check.names=F)
  Y <- read.table(mixture_file, header=T, sep="t", check.names=F)
  Y <- Y[!duplicated(Y[,1]),]
  rownames(Y)<-Y[,1]
  Y<-Y[,-1]
  X <- data.matrix(X)
  Y <- data.matrix(Y)

#order
X <- X[order(rownames(X)),]
Y <- Y[order(rownames(Y)),]

P <- perm #number of permutations

#anti-log if max < 50 in mixture file
if(max(Y) < 50) {Y <- 2^Y}

#quantile normalization of mixture file
if(QN == TRUE){
tmpc <- colnames(Y)
tmpr <- rownames(Y)
Y <- preprocessCore::normalize.quantiles(Y)
colnames(Y) <- tmpc
rownames(Y) <- tmpr
}

#intersect genes
Xgns <- row.names(X)
Ygns <- row.names(Y)
YintX <- Ygns %in% Xgns
Y <- Y[YintX,]
XintY <- Xgns %in% row.names(Y)
X <- X[XintY,]

#standardize sig matrix
X <- (X – mean(X)) / sd(as.vector(X))

# 矩阵X是LM22矩阵
# 矩阵Y是待预测的矩阵
# 然后对
#empirical null distribution of correlation coefficients
if(P > 0) {nulldist <- sort(doPerm(P, X, Y)$dist)}

#print(nulldist)

header <- c(‘Mixture’,colnames(X),“P-value”,“Correlation”,“RMSE”)
#print(header)

output <- matrix()
itor <- 1
mixtures <- dim(Y)[2]
pval <- 9999

#iterate through mixtures
print(paste0(length(colnames(Y)),” Samples in total”))
while(itor <= mixtures){

print(paste0(“Dealing with sample:”,itor))
y <- Y[,itor]

#standardize mixture
y <- (y – mean(y)) / sd(y)

#run SVR core algorithm
result <- CoreAlg(X, y)

#get results
w <- result$w
mix_r <- result$mix_r
mix_rmse <- result$mix_rmse

#calculate p-value
if(P > 0) {pval <- 1 – (which.min(abs(nulldist – mix_r)) / length(nulldist))}

#print output
out <- c(colnames(Y)[itor],w,pval,mix_r,mix_rmse)
if(itor == 1) {output <- out}
else {output <- rbind(output, out)}

itor <- itor + 1

}

#save results
## write.table(rbind(header,output), file=“CIBERSORT-Results.txt”, sep=“t”, row.names=F, col.names=F, quote=F)

#return matrix object containing all results
obj <- rbind(header,output)
obj <- obj[,-1]
obj <- obj[-1,]
obj <- matrix(as.numeric(unlist(obj)),nrow=nrow(obj))
rownames(obj) <- colnames(Y)
colnames(obj) <- c(colnames(X),“P-value”,“Correlation”,“RMSE”)
print(paste0(“Finished at:”,Sys.time()))
obj
}

耐心一点的可以一句一句运行探索这个函数
这个事情，我们常常做。

sig_matrix <- "CIBERSORT/LM22.txt"
mixture_file = 'exprMat.txt'
#read in data
## 读入signature数据
X <- read.table(sig_matrix,header=T,sep="t",row.names=1,check.names=F)
## 读入表达量数据
Y <- read.table(mixture_file, header=T, sep="t", check.names=F)
## 去除重复基因
Y <- Y[!duplicated(Y[,1]),]
## 第一列转行名
rownames(Y)<-Y[,1]
Y<-Y[,-1]
## 数据框转为矩阵
X <- data.matrix(X)
Y <- data.matrix(Y)

#order
## 行名排序
X <- X[order(rownames(X)),]
Y <- Y[order(rownames(Y)),]

这些操作就是分别读入表达矩阵和基因特征，最终产生基因数目一样的两个矩阵。
X现在是这样的

Y是这样的

接下来就是设置一下置换次数

## 设置置换次数
## #number of permutations
perm <- 1000
P <- perm 

下面这一步很重要，这个函数会自动检测你的数据。
如果你的数据被log话之后，他会自动帮你去掉log,
所以，CIBERSORT接受的数据是不去掉log的，如果你一不小心取了log，也不要害怕，他会帮你再去掉。
当然检测数据是否被log化，我们有更好的方法，GEO已经提供了。
来完成你的生信作业，这是最有诚意的GEO数据库教程

#anti-log if max < 50 in mixture file
## 这里说明，当他检测数据已经被log化之后，就会帮我们去log化
if(max(Y) < 50) {Y <- 2^Y}

为什么要这么做，因为他是要以基因集为依据来推断，我们看看基因特征的表达量发现。
这些数据有的很大，可以看出是没有取过log的。

接下里要对数据进行quantile 矫正
理解 Quntile Normalization

使用的是这个函数preprocessCore::normalize.quantile

QN=T
if(QN == TRUE){
  tmpc <- colnames(Y)
  tmpr <- rownames(Y)
  Y <- preprocessCore::normalize.quantiles(Y)
  colnames(Y) <- tmpc
  rownames(Y) <- tmpr
}

需要强调的是，只有芯片是需要这样做的，这跟芯片产生的过程有关。
而RNAseq数据，论文中推荐的是TPM

然后把两个数据的基因取交集，然后以交集的基因再去调整两个数据
用到的函数是%in%

Xgns <- row.names(X)
Ygns <- row.names(Y)
YintX <- Ygns %in% Xgns
Y <- Y[YintX,]
XintY <- Xgns %in% row.names(Y)
X <- X[XintY,]

现在X和Y的行都是517行

接下来他们对基因特征的数据进行了z-score转换
z-score的标准化究竟怎么弄？

X <- (X - mean(X)) / sd(as.vector(X))

好了准备工作做完了，现在正式开始

核心函数CoreAlgd

这个函数是主要计算的函数，他接受基因集以及一个样本的数据，返回一个列表

CoreAlg <- function(X, y){

#try different values of nu
svn_itor <- 3

res <- function(i){
if(i==1){nus <- 0.25}
if(i==2){nus <- 0.5}
if(i==3){nus <- 0.75}
model<-e1071::svm(X,y,type=“nu-regression”,kernel=“linear”,nu=nus,scale=F)
model
}

if(Sys.info()[‘sysname’] == ‘Windows’) out <- parallel::mclapply(1:svn_itor, res, mc.cores=1) else
out <- parallel::mclapply(1:svn_itor, res, mc.cores=svn_itor)

nusvm <- rep(0,svn_itor)
corrv <- rep(0,svn_itor)

#do cibersort
t <- 1
while(t <= svn_itor) {
weights = t(out[[t]]$coefs) %*% out[[t]]$SV
weights[which(weights<0)]<-0
w<-weights/sum(weights)
u <- sweep(X,MARGIN=2,w,‘*’)
k <- apply(u, 1, sum)
nusvm[t] <- sqrt((mean((k – y)^2)))
corrv[t] <- cor(k, y)
t <- t + 1
}

#pick best model
rmses <- nusvm
mn <- which.min(rmses)
model <- out[[mn]]

#get and normalize coefficients
q <- t(model$coefs) %*% model$SV
q[which(q<0)]<-0
w <- (q/sum(q))

mix_rmse <- rmses[mn]
mix_r <- corrv[mn]

newList <- list(“w” = w, “mix_rmse” = mix_rmse, “mix_r” = mix_r)

}

我们可以测试这个函数的结果

## 提取一个样本
y <- Y[,1]
y <- (y - mean(y)) / sd(y)
## 计算
result <- CoreAlg(X, y)
str(result)

返回的结果记录了三个结果，我们可以通过再次探索CoreAlg来确定那是三个返回值
w，mix_rmse,mix_r

这一段，我写不下去了，会录制成视频更新在答疑群中。
简单来说

w实际上就是就是个人除以总数，体现在最后的22个免疫细胞的值是个率，加起来等于1
RMSE就是均方根误差(Root Mean Squared Error)，越小效果越好，解释看这个帖子
花了100个小时学习线性回归，写了个万字长文作总结。
而相关性是这样算的。
首先这个样本是有每个基因的表达量的。
然后算出来每个免疫细胞的占比后，用这个占比去把基因矩阵中的每个值都按照这个比值处理了

u <- sweep(X,MARGIN=2,w,'*')

接下来，算出预测的单个样本中每个基因的表达量

k <- apply(u, 1, sum)

这个表达量k和本身的表达量y求相关性得到相关性系数。

置换函数doPerm

那么现在只剩下p值了，p值是怎么算的呢？
使用置换检验。
就是用蒙特卡罗方法，不断地从原来的数据中跟单个样本等长的数据，模拟新样本，进行CoreAlg计算。
每次抽取出相关性系数r，当次数足够多的时候，就形成背景相关性系数分布。
我做了1000次后是这样的。

这个计算由doPerm函数完成，这个方法未来在GSEA课程的时候还会碰到。

doPerm <- function(perm, X, Y){
  itor <- 1
  Ylist <- as.list(data.matrix(Y))
  dist <- matrix()
  print(paste0("do permutations:",perm," in total"))
  while(itor <= perm){
    print(paste0("nperm:",itor))
    #random mixture
    yr <- as.numeric(Ylist[sample(length(Ylist),dim(X)[1])])
    #standardize mixture
    yr <- (yr - mean(yr)) / sd(yr)
    #run CIBERSORT core algorithm
    result <- CoreAlg(X, yr)
    mix_r <- result$mix_r
    #store correlation
    if(itor == 1) {dist <- mix_r}
    else {dist <- rbind(dist, mix_r)}
    itor <- itor + 1
  }
  newList <- list("dist" = dist)
}

这个函数做的工作刚才已经讲了，很简单。
就是通过sample命令不断产生模拟样本，进行CoreAlg计算抽取r值。

yr <- as.numeric(Ylist[sample(length(Ylist),dim(X)[1])])

有了背景r分布，每次也能产生各个细胞的比率，有均方根误差RMSE，有相关性系数，还有p值，
这样最终结果的一行就有了。
那么只要批量对每一个样本计算即可得到所有结果
接下来的CIBERSORT主函数就用这句命令，计算的p值

if(P > 0) {pval <- 1 - (which.min(abs(nulldist - mix_r)) / length(nulldist))}

我强行解释一波，每次真正的样本计算出来的r值，找到其在背景分布的位置，然后看右边红色的概率有多大。如果小于0.05，说明产生这种极端值并不是偶然的。
那么该样本的反卷积解析就是成功的，可信的。

这样批量每个样本后，就返回了我们想要的结果，以obj代替。

如果样本多，我们还能可以先把p值小的样本删掉。
既然每个细胞的值是个比率，加起来等于1，那么我们就可以画图来看了。
画图用到的技能
我喜欢的gather快要被淘汰了，迎面走来了更好用的宽长转换工具

library(dplyr)
library(tidyr)
library(tibble)
dd1 <- obj %>% 
  as.data.frame() %>% 
  rownames_to_column("sample") %>% 
  pivot_longer(cols=2:23,
               names_to= "celltype",
               values_to = "Proportion")

现在的数据是这样的

有了数据画图就很容易

library(ggplot2)
ggplot(dd1,aes(sample,Proportion,fill = celltype)) + 
  geom_bar(position = "stack",stat = "identity")+
  theme_bw()+
  guides(fill=guide_legend(ncol=1))

看纵坐标加起来就是1。
但是ssGSEA的数据，最后是富集分数，画这种图是不合适的。这个会在之后的帖子讲解。

现在CIBERSORT的输入输出都知道了。用起来会很有底气。

还能怎么改进呢？
置换那里速度太慢了，设置成1000次，就跟多了1000个样本一样，我们可以设置并行化加速。

要补充一点的是，如果研究的TCGA数据，其实大部分软件都把结果计算好了。比如CIBERSORT计算的 TCGA结果在这里
https://gdc.cancer.gov/about-data/publications/panimmune

直接去下载读入R语言就可以用了。
那么讲这么长有什么用呢？
是这样的:
这个函数使用起来这么简单，只要我们替换掉基因特征数据，就可以迁移到其他物种，
比如，我们想要量化小鼠的组织，那怎么办呢，只要找到小鼠免疫细胞的这个特征集就行了。
还好，这个事情，不需要我们做，有人已经做好了，下次我们就讲这个。