CytoTRACE是一款基于單細胞計數(shù)矩陣推測細胞間活性和細胞間相對分化狀態(tài)的一款軟件揽碘，這里我們結合實際代碼例子龄章，解釋CytoTRACE的數(shù)學原理：
step 1:
這一步先將預處理后的單細胞矩陣進行表達基因數(shù)量的檢測，而處理的方法是過濾NA值和低表達的基因，利用行總和大于0進行篩選，并計算細胞相關性矩陣
# 載入相應的數(shù)據(jù)
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,]
# 定義子采樣的數(shù)據(jù)集大小
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){
# 分批次對單細胞矩陣進行l(wèi)og2標準化
## 若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和低質(zhì)量的細胞护蝶，每個細胞列的和小于零定義為低質(zhì)量細胞
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
# 將表達矩陣進行l(wèi)og2運算
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)
# 用基因計數(shù)重新調(diào)整每個單細胞翩迈，以在細胞裂解前將相對轉(zhuǎn)錄本豐度轉(zhuǎn)換為絕對 RNA 含量
## 這里的RNA豐度是絕對豐度持灰，將之前l(fā)og2運算后的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 等操作轉(zhuǎn)換為馬爾可夫矩陣
similarity_matrix_cleaned <- function(similarity_matrix){
D <- similarity_matrix
# 設定閾值來控制網(wǎng)絡邊的生成
cutoff <- mean(as.vector(D))
# 將對角線元素替換成 0
diag(D) <- 0;
# 將負數(shù)換為 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函數(shù)轉(zhuǎn)化被證實可以檢測到更多的差異基因喂链，而返回結果：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個基因并分別在每個細胞中求集合平均數(shù)
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）與基因計數(shù)最正相關的前 200 個基因在每個細胞中的幾何平均表達量向量 (GCS)；D2 為每一個組別（0纵揍，1顿乒，2）對應的相關性矩陣
#Regress gene counts signature (GCS) onto similarity matrix
# nnls回歸的目的是得到期望的GCS，A 為馬爾可夫矩陣泽谨，b為GCS
## 解讀 regressed 函數(shù)淆游，將會在后面 cytotrace 函數(shù)中出現(xiàn)
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向量進行迭代，優(yōu)化矯正期望的GCS向量
## 解讀 diffused 函數(shù)隔盛，將會在后面 cytotrace 函數(shù)中出現(xiàn)
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 函數(shù)
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向量達到平穩(wěn)
gcs_diffused <- diffused(D2[[i]], gcs_regressed)
# 按照GCS的大小進行排序
cytotrace <- rank(gcs_diffused)
}
)
cytotrace <- cytotrace_ranked <- unlist(cytotrace)
# 標準化GCS的值到0到1之間
cytotrace <- range01(cytotrace)
所得到的cytotrace為一列長度為細胞數(shù)量（3427個）的向量
gcs2代表GCS，即每個細胞中基因表達的幾何平均數(shù)
gcs2
nnls的計算模型為：
與基因計數(shù)最正相關的前 200 個基因在每個細胞中的幾何平均表達被定義為基因計數(shù)特征向量 (GCS)
在我們的例子中搜骡，建立了馬爾可夫矩陣與GCS表達的線性關系拂盯，相關性矩陣為 A，GCS為 b
并且score_regressed <- similarity_matrix_cleaned %*% out$x得到的是期望的GCS
利用馬爾可夫隨機過程和對期望的GCS向量進行矯正记靡，使期望的GCS達到平穩(wěn)
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（與基因計數(shù)最正相關的前 200 個基因在每個細胞中的幾何平均表達被定義為基因計數(shù)特征向量）的回歸關系团驱，得到期望的GCS，利用馬爾可夫隨機過程對期望的GCS向量進行矯正空凸，使期望的GCS向量達到平穩(wěn)嚎花，那么GCS的值被定義為每個細胞的活性
