量化免疫浸润时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化之后,就会帮我们去logif(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 %inXgns
Y <- Y[YintX,]
XintY <- Xgns %inrow.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=1else
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语言就可以用了。
那么讲这么长有什么用呢?
是这样的:
这个函数使用起来这么简单,只要我们替换掉基因特征数据,就可以迁移到其他物种,
比如,我们想要量化小鼠的组织,那怎么办呢,只要找到小鼠免疫细胞的这个特征集就行了。
还好,这个事情,不需要我们做,有人已经做好了,下次我们就讲这个。

生物信息学

如何使用R语言进行肿瘤样本纯度计算?

2020-5-22 14:45:40

生物信息学

生信分析时常用的数据格式

2020-5-28 21:06:26

声明 本网站部分文章源于互联网,出于传递更多信息和学习之目的转载,并不保证内容正确或赞同其观点。
如转载稿涉及失效、版权等问题,请立即联系管理员;我们会予以修改、删除相关文章,请留言反馈
Notice: When your legal rights are being violated, please send an email to: [email protected]
2 条回复 A文章作者 M管理员
  1. 498601851

    您好,请问为什么会出现收捲时出错: Model is empty!
    Error: no more error handlers available (recursive errors?); invoking ‘abort’ restart 报错呀?

    • 落落

      您好,我也出现了同样问题,请问解决了么

个人中心
购物车
优惠劵
今日签到
有新私信 私信列表
搜索