简介

CIBERSORT这款软件利用反卷积的方法，利用单细胞RNA-seq的数据，提取特征后，反推Bulk-seq各类细胞成分所占比例。这款软件目前推出了CIBERSORTx，但是账号貌似有点难申请
而这次我们着重介绍下CIBERSORT的R版本
而CIBERSORT计算细胞组分的方法使用去卷积原理实现的，而且通常是线性去卷积

原理图

R

核心代码由三部分组成：CoreAlg，doPerm，CIBERSORT，其中该软件利用的是svm反卷积进行推算：

#读取文件
X <- read.table(sig_matrix,header=T,sep="\t",row.names=1,check.names=F)
Y <- read.table(mixture_file, header=T, sep="\t", row.names=1,check.names=F)
  
X <- data.matrix(X)
Y <- data.matrix(Y)
  
#order
X <- X[order(rownames(X)),]
Y <- Y[order(rownames(Y)),]

其中mixture_file是RNA-seq的数据；sig_matrix为单细胞RNA-seq数据

mixture_file

sig_matrix

紧接着对两个文件的基因取交集并且分别标准化：

#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))

#standardize mixture
y <- Y[,itor]
y <- (y - mean(y)) / sd(y)
#这里的itor遍历的是mixture_file的列元素，即y是mixture_file的列元素

注意：这里的 itor 遍历的是mixture_file的列元素，即y是mixture_file的列元素，表示的是Bulk-seq的每一个sample，而x代表单细胞数据

接下来就是整个代码的核心段了，涉及到CoreAlg：

#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<-svm(X,y,type="nu-regression",kernel="linear",nu=nus,scale=F)
    model
}
  
#输出结果
out <- mclapply(1:svn_itor, res, mc.cores=1)
##svn_itor遍历不同的nu

关于参数type:

1.If you have a classification problem, i.e., discrete label to predict, you can use C-classification and nu-classification.
2.If you have a regression problem, i.e., continuous number to predict, you can use eps-regression and nu-regression.
3.If you only have one class of the data, i.e., normal behavior, and want to detect outliers. one-classification.

这一部分res()函数通过调用svm()函数，并调试不同的nu值建立不同的模型（这里 nu 值为1-3），而对于整个svm的过程可以理解为回归问题（SVR，即y为响应变量；X为决策变量，每种细胞类型相当于每一维的数据）

这里简单介绍下SVR（支持向量机回归）：

来自微信

接下来需要将out的结果（里面囊括3个svm模型）循环赋值出来，并且定义相关的权重：

#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
 } 

##svn_itor遍历不同的nu

这里的 t 代表不同参数（nu，一共3个值）下的svm模型
这里的权重（weight）定义为：svm的coefs与SV的乘积（svm模型$一下可以找到这个两个系数），这个乘积表示为线性回归平面的法向量（权重）weight，并且权重weight小于0的都定义为0；w代表的是权重矩阵的每一个元素除以该矩阵元素的总和 ；这里的u为X矩阵（单细胞表达谱）的每一个元素减去w；定义k为u矩阵每一行求和；定义
nusvm为 sqrt((mean((k - y)^2)))；定义corrv为 corrv[t] <- cor(k, y)，即k和y的相关系数

注明：
svm的coefs指的是svm模型的系数；而SV指的是支持向量，支持向量是数据集的点，它们靠近分隔类别的平面

coefs

SV

那么怎么理解这个贡献呢？

举个例子：

1. 其中，在bulk的表达矩阵中，gene1 的表达量 E1 = β1×e1 + β2×e4，其中 e1 和 e2 分别是 cell cluster A 和 cell cluster B gene1 的表达量，同理可知其他基因的，也就是 bulk 中 gene1 的表达量由 cell cluster A 和 cell cluster B gene1 的表达量组成，而系数 β1 和 β2 理解为细胞比例
2. 从图像来看，假设 gene 多的话，单独分离开 cell cluster A 和 cell cluster B来看：
   图中的红点代表不同基因的表达量，由于响应变量代表基因在bulk中的表达量，因此纵坐标在 cell cluster A 和 cell cluster B 是一样的。显然在 cell cluster A 中的基因表达量要偏低一些，而在 cell cluster B 中表达量偏高一些；因此基因表达量偏高的细胞，其比例会偏少，因此基因表达量偏低的细胞，其比例会偏多（斜率 β1 > β2； α > θ） 注明 E1 = β1×e1 + β2×e4，bulk的表达量由对应 cell cluster A 和 cell cluster B 的基因表达量组成
3. 这种比例应该理解为SVR多元回归中拟合出来的情况，只不过借助线性方程E_bulk = β1×e_cell_A + β2×e_cell_B的加和状态，将参数β1和β2比作为细胞比例（细胞的量）乘相应的cell cluster对应基因的表达量，作为最终bulk的基因表达量

然后从 out 中选择最优模型：

#pick best model
rmses <- nusvm
##令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]  
#提取最优模型的rmses(nusvm)和corrv

那么经过CoreAlg计算出来的结果为：

result = CoreAlg(X, y)

那么相应的w，mix_r和mix_rmse为：

w <- result$w
mix_r <- result$mix_r
mix_rmse <- result$mix_rmse

其中w为：

w

然而，这只是mixture_file(Bulk-seq)的一列数据，接下来我们计算每一列的rmses(nusvm)和corrv

while(itor <= mixtures){
   
    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 
    ##这里的itor遍历的是mixture_file的列元素，即y是mixture_file的列元素
  }

最后输出的结果为bulk-seq每一个sample的w矩阵：

即每个组分所占比例

小tip

这个官方代码中quantile的标准化方式可能会报错，不妨在CIBERSORT函数下修改如下：

quantile_normalisation <- function(df){
    df_rank <- apply(df,2,rank,ties.method="min")
    df_sorted <- data.frame(apply(df, 2, sort))
    df_mean <- apply(df_sorted, 1, mean)
    
    index_to_mean <- function(my_index, my_mean){
      return(my_mean[my_index])
    }
    df_final <- apply(df_rank, 2, index_to_mean, my_mean=df_mean)
    rownames(df_final) <- rownames(df)
    return(df_final)
  }

##改动如下：
##原：
if(QN == TRUE){
    tmpc <- colnames(Y)
    tmpr <- rownames(Y)
    Y <- normalize.quantiles(Y)
    colnames(Y) <- tmpc
    rownames(Y) <- tmpr
  }

##改：
 if(QN == TRUE){
    tmpc <- colnames(Y)
    tmpr <- rownames(Y)
    Y <- quantile_normalisation(Y)
    colnames(Y) <- tmpc
    rownames(Y) <- tmpr
  }

另外，有关p值计算，该软件采用的是置换检验的方式，详见微信文章：《量化免疫浸润时CIBERSORT的注意事项》