CytoTRACE是一款基于单细胞计数矩阵推测细胞间活性和细胞间相对分化状态的一款软件，这里我们结合实际代码例子，解释CytoTRACE的数学原理：
step 1:
这一步先将预处理后的单细胞矩阵进行表达基因数量的检测，而处理的方法是过滤NA值和低表达的基因，利用行总和大于0进行筛选，并计算细胞相关性矩阵
# 载入相应的数据
load('C:/Users/Downloads/CytoTRACE/data/marrow_10x_expr.rda')
# mat 为13517行3427列的矩阵
mat = marrow_10x_expr
batch = NULL
enableFast = TRUE
ncores = 1
subsamplesize = 1000
range01 <- function(x){(x-min(x))/(max(x)-min(x))}
#inputs
a1 <- mat
a2 <- batch
#Checkpoint: NAs and poor quality genes
#过滤NA值和低表达的基因，利用行总和大于0进行筛选
pqgenes <- is.na(rowSums(mat>0)) | apply(mat, 1, var) == 0
num_pqgenes <- length(which(pqgenes == TRUE))
mat <- mat[!pqgenes,]
# 定义子采样的数据集大小
size <- subsamplesize
# 定义子采样的组别
chunk <- round(ncol(mat)/size)
subsamples <- split(1:ncol(mat), sample(factor(1:ncol(mat) %% chunk)))
这里的子采样分为三组0，1，2；相当于将细胞分为0，1，2三组
step 2：
这一步做基因表达矩阵的标准化，对基因进行过滤并构建马尔可夫矩阵
batches <- parallel::mclapply(subsamples, mc.cores = min(chunk, ncores), function(subsample){
# 分批次对单细胞矩阵进行log2标准化
## 若batch = NULL 则忽略
mat <- mat[,subsample]
batch <- batch[subsample]
if(max(mat)<50){
mat <- 2^mat - 1
}
#Checkpoint: ERCC standards
#去除线粒体基因的表达
if(length(grep("ERCC-", rownames(mat)))>0){
mat <- mat[-grep("ERCC-", rownames(mat)),]
}
#Checkpoint: Sequencing depth normalization
#对每个细胞的测序深度进行标准化
mat <- t(t(mat)/apply(mat, 2, sum))*1000000
#Checkpoint: NAs and poor quality cells
# 去除NA和低质量的细胞，每个细胞列的和小于零定义为低质量细胞
pqcells <- is.na(apply(mat>0, 2, sum)) | apply(mat>0, 2, sum) <= 10
num_pqcells <- length(which(pqcells == TRUE))
mat <- mat[,!pqcells]
#Checkpoint: log2-normalize
# 将表达矩阵进行log2运算
mat <- log(mat+1,2)
mat <- data.matrix(mat)
# 矫正批次效应，若batch = NULL则忽略
#Calculate pre-batch corrected gene counts
counts <- apply(mat>0, 2, sum)
#Checkpoint: Batch correction
if(ncol(a1) == length(a2)){
#filter poor quality cells from batch vector
batch <- batch[!pqcells]
#Run Combat
# 用Combat去除批次效应
suppressMessages(mat <- sva::ComBat(mat, batch, c()))
mat <- data.matrix(mat)
#Replace negative values after batch correction with zeroes for compatibility with downstream steps
# 将单细胞矩阵中的负值替换成 0
mat[which(mat<0)] <- 0
}
#Rescale each single cell with gene counts to convert relative transcript abundances to absolute RNA content prior to cell lysis (credit: Census, Qiu et al., 2017)
# 用基因计数重新调整每个单细胞，以在细胞裂解前将相对转录本丰度转换为绝对 RNA 含量
## 这里的RNA丰度是绝对丰度，将之前log2运算后的mat进行绝对丰度的计算
census_normalize <- function(mat, counts) {
xnl <- 2^data.matrix(mat) - 1
rs <- apply(xnl, 2, sum)
rnorm <- t(t(xnl) * counts/rs)
A <- log(rnorm+1,2)
return(A)
}
mat2 <- census_normalize(mat, counts)
#Function to identify the most variable genes
# 寻找高方差的基因
mvg <- function(matn) {
A <- matn
n_expr <- rowSums(A > 0);
A_filt <- A[n_expr >= 0.05 * ncol(A),];
vars <- apply(A_filt, 1, var);
means <- apply(A_filt, 1, mean);
# 计算离散度disp
disp <- vars / means;
# 筛选出1000个离散度大的基因(从大到小排序)
# last_disp 为按离散度从大到小排序的第1000个基因的离散度
last_disp <- tail(sort(disp), 1000)[1];
# 筛选出离散度大于 last_disp 的基因
A_filt <- A_filt[disp >= last_disp,];
return(A_filt)
}
#Filter out cells not expressing any of the 1000 most variable genes
# 将离散度大的1000个基因中没有表达的基因删除掉
mat2.mvg <- mvg(mat2)
rm1 <- colSums(mat2.mvg) == 0
mat2 <- mat2[, !rm1]
counts <- counts[!rm1]
#Calculate similarity matrix
# 计算马尔可夫矩阵，即将相关性矩阵通过对角线取 0 等操作转换为马尔可夫矩阵
similarity_matrix_cleaned <- function(similarity_matrix){
D <- similarity_matrix
# 设定阈值来控制网络边的生成
cutoff <- mean(as.vector(D))
# 将对角线元素替换成 0
diag(D) <- 0;
# 将负数换为 0
D[which(D < 0)] <- 0;
# 小于cutoff的换为 0
D[which(D <= cutoff)] <- 0;
Ds <- D
D <- D / rowSums(D);
D[which(rowSums(Ds)==0),] <- 0
return(D)
}
# 利用基因的方差信息筛选出来的那1000个基因的RNA绝对定量矩阵计算马尔可夫矩阵
D <- similarity_matrix_cleaned(HiClimR::fastCor(mvg(mat2)))
# mat2 为RNA绝对定量矩阵，counts 为pre-batch corrected gene counts可以理解为每个细胞的测序深度(每个细胞中基因表达counts总和），D为马尔可夫矩阵
return(list(mat2 = mat2,counts = counts, D = D))
})
其中census_normalize函数转化被证实可以检测到更多的差异基因，而返回结果：mat2 为RNA绝对定量矩阵，counts 为pre-batch corrected gene counts可以理解为每个细胞的测序深度(每个细胞中基因表达counts总和），D为马尔可夫矩阵
D矩阵如下：
D矩阵
HiClimR::fastCor(mvg(mat2)) 计算的相关性
HiClimR::fastCor(mvg(mat2))
step 3：
利用nnls回归得到真实的GCS，，结合期望马尔可夫矩阵的信息，对真实的GCS进行矫正，使之更加合理
#Prepare for downstream steps
# 准备下面分析的基因表达矩阵
## 合并每个组别（0，1，2）对应的单细胞表达矩阵
mat2 <- do.call(cbind, lapply(batches, function(x) x$mat2))
# 合并每个组别（0，1，2）每个细胞对应的测序深度
counts <- do.call(c, lapply(batches, function(x) x$counts))
filter <- colnames(a1)[-which(colnames(a1) %in% colnames(mat2))]
#Calculate gene counts signature (GCS) or the genes most correlated with gene counts
# 对利用单细胞表达矩阵mat2的每一行(每个基因在不同细胞中的表达量向量)，与细胞测序深度向量 counts 计算相关性
ds2 <- sapply(1:nrow(mat2), function(x) ccaPP::corPearson(mat2[x,],counts))
names(ds2) <- rownames(mat2)
# 提取相关性高的前200个基因并分别在每个细胞中求集合平均数
gcs <- apply(mat2[which(rownames(mat2) %in% names(rev(sort(ds2))[1:200])),],2,mean)
# 定义细胞测序深度并命名
samplesize <- unlist(lapply(lapply(batches, function(x) x$counts), length))
gcs2 <- split(gcs, as.numeric(rep(names(samplesize), samplesize)))
# 提取每个组别（0，1，2）相关性矩阵
D2 <- lapply(batches, function(x) x$D)
## gcs2为每个组别（0，1，2）与基因计数最正相关的前 200 个基因在每个细胞中的几何平均表达量向量 (GCS)；D2 为每一个组别（0，1，2）对应的相关性矩阵
#Regress gene counts signature (GCS) onto similarity matrix
# nnls回归的目的是得到期望的GCS，A 为马尔可夫矩阵，b为GCS
## 解读 regressed 函数，将会在后面 cytotrace 函数中出现
regressed <- function(similarity_matrix_cleaned, score){
out <- nnls::nnls(similarity_matrix_cleaned,score)
# 计算期望的GCS
score_regressed <- similarity_matrix_cleaned %*% out$x
return(score_regressed)
}
#Apply diffusion to regressed GCS using similarity matrix
# 利用马尔可夫随机过程对期望的GCS向量进行迭代，优化矫正期望的GCS向量
## 解读 diffused 函数，将会在后面 cytotrace 函数中出现
diffused <- function(similarity_matrix_cleaned, score, ALPHA = 0.9){
vals <- score
# 向量化，v_prev 和v_curr 为向量
v_prev <- rep(vals);
v_curr <- rep(vals);
for(i in 1:10000) {
v_prev <- rep(v_curr);
v_curr <- ALPHA * (similarity_matrix_cleaned %*% v_curr) + (1 - ALPHA) * vals;
diff <- mean(abs(v_curr - v_prev));
# 迭代至误差最小
if(diff <= 1e-6) {
break;
}
}
return(v_curr)
}
# 解读 cytotrace 函数
cytotrace <- parallel::mclapply(1:length(D2), mc.cores = ncores, function(i) {
# 分别每个组别（0，1，2）的马尔可夫矩阵和每个组别（0，1，2）的GCS做nnls回归，完成每个组别（0，1，2）所对应相关性矩阵的矫正
gcs_regressed <- regressed(D2[[i]], gcs2[[i]])
# 利用马尔可夫矩阵对期望的GCS向量进行矫正，是期望的GCS向量达到平稳
gcs_diffused <- diffused(D2[[i]], gcs_regressed)
# 按照GCS的大小进行排序
cytotrace <- rank(gcs_diffused)
}
)
cytotrace <- cytotrace_ranked <- unlist(cytotrace)
# 标准化GCS的值到0到1之间
cytotrace <- range01(cytotrace)
所得到的cytotrace为一列长度为细胞数量（3427个）的向量
gcs2代表GCS，即每个细胞中基因表达的几何平均数
gcs2
nnls的计算模型为：
与基因计数最正相关的前 200 个基因在每个细胞中的几何平均表达被定义为基因计数特征向量 (GCS)
在我们的例子中，建立了马尔可夫矩阵与GCS表达的线性关系，相关性矩阵为 A，GCS为 b
并且score_regressed <- similarity_matrix_cleaned %*% out$x得到的是期望的GCS
利用马尔可夫随机过程和对期望的GCS向量进行矫正，使期望的GCS达到平稳
step 4：
筛选结果
#Calculate genes associated with CytoTRACE
# 计算每一行基因 (每个基因在不同细胞中的表达量向量) 与 cytotrace （姑且命名为cytotrace）计算相关性
cytogenes <- sapply(1:nrow(mat2),
function(x) ccaPP::corPearson(mat2[x,], cytotrace))
names(cytogenes) <- rownames(mat2)
#Final steps
# 对结果进行筛选
names(cytotrace) <- names(cytotrace_ranked) <- names(gcs) <- names(counts) <- colnames(mat2)
cytotrace <- cytotrace[colnames(a1)]
cytotrace_ranked <- cytotrace_ranked[colnames(a1)]
gcs <- gcs[colnames(a1)]
counts <- counts[colnames(a1)]
mat2 <- t(data.frame(t(mat2))[colnames(a1),])
names(cytotrace) <- names(cytotrace_ranked) <- names(gcs) <- names(counts) <- colnames(mat2) <- colnames(a1)
总结
该软件的核心是利用细胞间相关性构建马尔可夫矩阵，并利用nnls建立马尔可夫矩阵与GCS（与基因计数最正相关的前 200 个基因在每个细胞中的几何平均表达被定义为基因计数特征向量）的回归关系，得到期望的GCS，利用马尔可夫随机过程对期望的GCS向量进行矫正，使期望的GCS向量达到平稳，那么GCS的值被定义为每个细胞的活性
