认识 DNA 测序的朋友应该是知道对肿瘤病人进行 WES 或者 WGS 等基因组测序后，可以得到成百上千的 somatic 突变，而这些突变都是 ATCG 碱基直接 变化，所以它们变化的组合情况就 6 种，而加上上下文碱基也就 96 种，这 96 种碱基变化的比例的特殊组成，就是 mutation signatures，而且 COSMIC 数据库上面有着 30 种已知的 signatures，我们可以把这些 signatures 当做是一个有意义的生物学功能，这样一旦我们拿到自己的突变数据， 就可以通过非负矩阵分解的方法把自己的突变数据分解为这 30 个 signatures 的组合，那么假如我们拿到的是表达矩阵呢？
由于表达矩阵通常是 bulk 转录组测序，也就是说本来就是肿瘤细胞以及其肿瘤微环境的各种其他细胞组合而成，同理我们应该是可以根据表达量推断出来他们的细胞组分，当然，这个就需要算法上面的突破啦，下面我们就先容一些相关方向的进展。
免疫微环境背景知识
低通量实验可以用免疫组化,免疫荧光和流式等方法来获得组织免疫细胞组成，而 bulk tumor 转录组数据同样可以通过算法轻松得到肿瘤免疫细胞浸润水平。两篇综述系统性的收拾整顿和比较了主流算法，
Hackl H, Charoentong P, Finotello F, et al. Computational genomics tools for dissecting tumour-immune cell interactions. Nat Rev Genet 2016;17:441-58.Finotello, Francesca, and Zlatko Trajanoski. Quantifying Tumor-Infiltrating Immune Cells from Transcriptomics Data. Cancer Immunology, Immunotherapy: CII 67, no. 7 (2018): 1031–40.
提到的算法主要是 2 大类，包括：
基于 GSEA 的半定量方法Deconvolution algorithms（分为 partial deconvolution 和 complete deconvolution）
基于 GSEA 的半定量方法
基于 GSEA 半定量的方法最常用的就是 ssGSEA(单样本 GSEA)，ssGSEA 按照单个样本表达量的绝对值排序，计算特定 gene set 的累积经验概率分布，得到 ES 值代表这些基因在单个样本中协调上下调的程度，以此来代表免疫浸润细胞丰度。xCell，MCP-counter 都是基于 ssGSEA 所开发出来的方法。
partial deconvolution
Deconvolution algorithms 可以理解为一个基因在样本中的表达量是该基因在样本中不同细胞亚群表达水平和细胞分数权重的线性组合。近几年开发的CIBERSORT 算法恰是去卷积方法的应用。
TIMER 量化 6 种免疫细胞，但是与 CIBERSORT 不同（CIBERSORT 解析结果：22 种免疫细胞相加的总占比为 100%），TIMER 没有把预测值标准化为 1，故不可以把结果解释为细胞分数或是在不同的数据集中比较。还有 EPIC，quanTIseq 等等方法。
complete deconvolution
相比较与 partial deconvolution，complete deconvolution 不仅可以估计相对细胞分数同时还能 disentangle 表达谱。非负矩阵分解（NMF）就是一个很好的理解例子，但是 NMF 是完全无监视的方法，考虑到量化免疫细胞的要求，Gaujoux 融合细胞 marker genes 的先验知识开发出半监视 NMF，有爱好了解的可以了解下 CellMix 这个 R 包。
量化工作中的 challenges
Bulk tumor 中，除了各类方法想要量化的细胞，还有一些未知的细胞。Deconvolution algorithms 其中一个挑战就是对这些 unknown tumor content 的稳健性。EPIC 和 quanTIseq 答应终极细胞分数总和低于 1，来估计未知细胞分数。
此外，多重共线性也是 deconvolution algorithm 要考虑的问题，如基因表达在 related cells 中的高度相关。CIBERSORT，TIMER，quanTIseq 采用了一系列的方法来减少共线性的影响。
今天，我们先一起来拆解一下 CIBERSORT 流程，配套 R 代码来自于 CIBERSORT R script v1.03 (last updated 07-10-2015)，最新版是 CIBERSORTX，我们下一次再解析。
相比较与基因芯片，RNA-seq 数据具有“digital” nature 和 wider dynamic range 使得量化工作更具有挑战性，目前只有 EPIC 和 quanTIseq 是针对 RNAseq 数据的。mRNA content，输入数据的统计学分布以及是否应该针对不同癌种，肿瘤类型（实体瘤和血液肿瘤）开发个性化的量化方法都是值得思索和改进的问题。
第一步，认识 LM22 和表达矩阵
代码如下，无需运行，我已经把两个矩阵留存为 input.Rdata 文件啦
rm(list = ls())options(stringsAsFactors = F)# 首先读取两个文件sig_matrix <-"LM22-ref.txt" # cibersoft 内置数据库挖掘mixture_file <- "mRNA2.txt" # 约 80M，TCGA 数据库# 两个表达矩阵需要取交集#read in dataX <- read.table(sig_matrix,header=T,sep=" ",row.names=1,check.names=F)Y <- read.table(mixture_file, header=T, sep=" ", check.names=F)Y <- Y[!duplicated(Y[,1]),]rownames(Y)<-Y[,1]Y<-Y[,-1]X <- data.matrix(X)Y <- data.matrix(Y)Y[1:4,1:4]X[1:4,1:4]dim(X)dim(Y)X <- X[order(rownames(X)),]Y <- Y[order(rownames(Y)),]#anti-log if max < 50 in mixture fileif(max(Y) < 50) {Y <- 2^Y}QN = F#quantile normalization of mixture fileif(QN == TRUE){ tmpc <- colnames(Y) tmpr <- rownames(Y) Y <- normalize.quantiles(Y) colnames(Y) <- tmpc rownames(Y) <- tmpr}#intersect genesXgns <- row.names(X)Ygns <- row.names(Y)YintX <- Ygns %in% XgnsY <- Y[YintX,]XintY <- Xgns %in% row.names(Y)X <- X[XintY,]dim(X)dim(Y)#standardize sig matrixX <- (X - mean(X)) / sd(as.vector(X))Y[1:4,1:4]X[1:4,1:4]boxplot(X[,1:4])save(X,Y,file = "input.Rdata")
从代码里面，可以看到，LM22 矩阵进行了 zscore 转换，待反卷积的表达矩阵固然并没有被 zscore，不外在后面的步骤仍旧是会 zscore 的。
第二步，搞清楚 CoreAlg 函数
这里会使用被 zscore 后的 LM22 矩阵，来使用 SVM 算法来猜测一个随机的 Y 变量。
rm(list = ls())options(stringsAsFactors = F)library(preprocessCore)library(parallel)library(e1071)source("CIBERSORT.R") load(file = "input.Rdata")Y[1:4,1:4]X[1:4,1:4]dim(X)dim(Y)# 下面的演示是为了搞清楚 CoreAlg 函数# 并不需要留存任何信息# 从表达矩阵 Y 里面，随机挑选 LM22 矩阵基因数目的表达量值Ylist <- as.list(data.matrix(Y))yr <- as.numeric(Ylist[ sample(length(Ylist),dim(X)[1]) ])# yr 这个时候是一个假设的样本#standardize mixture,就是 scale 函数yr <- (yr - mean(yr)) / sd(yr)boxplot(yr)# 每次随机挑选的 yr，都是需要走后面的流程# 一切都是默认值的支持向量机# 这里的 X 是 LM22 矩阵，不同的免疫细胞比例组合成为不同的 yr# 这里的 yr 是随机的，反推免疫细胞比例out=svm(X,yr)outout$SV# SVM 需要自行学习哦# 需要修改的参数包括：type="nu-regression",kernel="linear",nu=nus,scale=Fsvn_itor <- 3y=yrres <- 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}#Execute In a parallel way the SVMif(Sys.info()["sysname"] == "Windows") { out <- mclapply(1:svn_itor, res, mc.cores=1) }else { out <- mclapply(1:svn_itor, res, mc.cores=svn_itor)}# 运行了 Support Vector Machines，函数是 svm {e1071}out#Initiate two variables with 0nusvm <- rep(0,svn_itor)corrv <- rep(0,svn_itor)t <- 1while(t <= svn_itor) { # 得到两个向量之间矩阵乘法的权重，此时应该只得到一个数字。 # 这样做是乘以系数 # 支持向量是数据集的点，它们靠近分隔种别的平面 # 现在的问题是，我没有任何种别（离散变量，例如“运动”、“片子”），但我有一个连续变量 mySupportVectors <- out[[t]]$SV # 系数定义 myCoefficients <- out[[t]]$coefs weights = t(myCoefficients) %*% mySupportVectors # 设置权重和相关性 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 modelrmses <- nusvmcorrvmn <- which.min(rmses)mn model <- out[[mn]]# 从 nus 为 0.25,0.5,0.75 的 3 个模型里面挑选一个即可#get and normalize coefficientsq <- t(model$coefs) %*% model$SV q[which(q<0)]<-0# w 就是计算后的 22 种免疫细胞的比例w <- (q/sum(q))mix_rmse <- rmses[mn]mix_r <- corrv[mn]# 会返回这个随机的 y 的免疫细胞组成情况，就是权重 wnewList <- list("w" = w, "mix_rmse" = mix_rmse, "mix_r" = mix_r)newList# 根据对应的权重与参考集相乘u <- sweep(X,MARGIN=2,w,"*") k <- apply(u, 1, sum)plot(y,k)sqrt((mean((k - y)^2))) cor(k, y)# 通常这个猜测结果都惨不忍睹
假如你硬是要去了解 SMV 原理，可能得去购买专业书籍了：
SVM 原理
第 3 步，把 CoreAlg 函数运行 1000 次
每次随机挑选基因表达量，天生一个模拟的样本，然后使用 zscore 后的 LM22 矩阵经由 SVM 算法来猜测：
rm(list = ls())options(stringsAsFactors = F)load(file = "input.Rdata")Y[1:4,1:4]X[1:4,1:4]dim(X)dim(Y)library(preprocessCore)library(parallel)library(e1071)source("cibersort_ann.R") itor <- 1Ylist <- as.list(data.matrix(Y))dist <- matrix()# 就是把 CoreAlg 函数运行 1000 次perm=1000while(itor <= perm){ print(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)nulldist=sort(newList$dist) # 这个 nulldist 主要是用来计算 P 值if(F){ P=perm #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)save(nulldist,file = "nulldist_perm_1000.Rdata")load(file = "nulldist_perm_1000.Rdata")print(nulldist)fivenum(print(nulldist))
每次猜测的 Y 和随机天生的 y 差异都非常可怕，简直就是劫难现场！留存每次随机过程得到猜测变量和随机变量的相关性系数，固然是都是 0 四周，由于猜测的很糟糕，基本上都不相关。
> fivenum(nulldist)[1] -0.079400079 -0.009724867 0.019402797 0.056885990 0.604810742
这个 nulldist 是 1000 次的相关性系数值，只有一个用处，就是计算后面真正的样本猜测的 P 值而已。
第四步，对表达矩阵的每个样本进行 SVM 猜测
理解了前面的代码，下面代码就非常容易理解，只不过是之前的猜测变量 y 是随机的，这次是真实样本的基因表达量 y 值：
rm(list = ls())options(stringsAsFactors = F)load(file = "input.Rdata")Y[1:4,1:4]X[1:4,1:4]dim(X)dim(Y)library(preprocessCore)library(parallel)library(e1071)source("cibersort_ann.R") header <- c("Mixture",colnames(X),"P-value","Correlation","RMSE")print(header)load(file = "nulldist_perm_1000.Rdata")output <- matrix()itor <- 1mix <- dim(Y)[2]pval <- 9999P=1000# 表达矩阵的每个样本，都需要计算一下 LM22 的比例#iterate through mixwhile(itor <= mix){ ################################## ## Analyze the first mixed sample ################################## y <- Y[,itor] #标准化样本数据集 y <- (y - mean(y)) / sd(y) #执行 SVR 核心算法 result <- CoreAlg(X, y) #获得结果 w <- result$w mix_r <- result$mix_r mix_rmse <- result$mix_rmse #计算 p-value if(P > 0) {pval <- 1 - (which.min(abs(nulldist - mix_r)) / length(nulldist))} #输出 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}head(output)#save resultswrite.table(rbind(header,output), file="CIBERSORT-Results.txt", sep=" ", row.names=F, col.names=F, quote=F)#return matrix object containing all resultsobj <- 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")objsave(obj,file = "output_obj.Rdata")
第五步，可视化 cibersoft 去卷积结果
结果很轻易理解，就是每个样本的各个免疫细胞比例而已，再加上 "P.value","Correlation","RMSE" 这 3 列：
