前言
最近的文章：《Gene signature extraction and cell identity recognition at the single-cell level with Cell-ID》，作者在这篇文章中提出了一种新的降维方式和寻找每个cell中特征基因的算法。作者仍未传统的PCA，ICA，t-SNE和UMAP降维后，仅仅只是在cell cluster（细胞亚群水平）中寻找基因的特征基因
而作者利用MCA降维以后，可以针对每一个单独的细胞寻找它们的特征基因
大致步骤
软件Cell ID工作原理大致分为三步：
step 1：
将单细胞的计数矩阵进行标准化和对数转换，使矩阵内的值介于0-1之间，此时记为M矩阵
step 2：
利用MCA对M矩阵进行降维，降至低维空间，此时称作MCA空间
step 3：
在MCA空间计算每个cell中，基因与cell之间的距离，并展开排序，排名越靠前，越代表这些基因是该cell的特征基因
作者通过MCA降维后，在低维空间里面分别获得了cell和基因在低维空间的坐标矩阵，那么基因与cell的距离，就用低维空间每个基因和每个cell坐标的欧氏距离表示
数学原理
MCA降维的具体原理：
假设标准化后的单细胞矩阵为 MN,P，N代表N cells；P代表 P genes；而 mn,p 代表cell n，gene p的表达水平，Mp 第p个基因的列向量
定义：
x + = f + (m) : R->R[0,1]
x - = f - (m) : R->R[0,1]；其中 f - (m) = 1 - f + (m)
将其标准化到0—1之间，这种编码模式来自于文献：Fuzzy coding
经过Fuzzy coding后的矩阵为：
因此标准化后新构成的 MN,K 矩阵如下：
K=2P
然后计算频率矩阵：
其中， NP代表总的样本数，X代表上述X矩阵，FN,K代表X矩阵第N行，第K列的数值
接下来定义两个对角矩阵：
定义相关频率矩阵 SNK：
其中：
Dr = [drn,n]
Dc = [dck,k]
F为上述FN,K矩阵
然后对相关频率矩阵 SNK 进行SVD分解：
SVD分解的意义是寻找cells和基因在低维空间的坐标
[1]. 那么对于cell在低维空间的坐标定义为：其中Φ包含了每个cell在低维空间的坐标
[2]. 那么对于基因在低维空间的坐标定义为：其中G包含了每个基因在低维空间的坐标
每一个细胞的特征基因定义为：
其中，
代码分析
1.MCA降维
读取数据以及对单细胞表达矩阵进行标准化
library(CelliD)
BaronMatrix <- readRDS(url("https://storage.googleapis.com/cellid-cbl/BaronMatrix.rds"))
BaronMetaData <- readRDS(url("https://storage.googleapis.com/cellid-cbl/BaronMetaData.rds"))
SegerMatrix <- readRDS(url("https://storage.googleapis.com/cellid-cbl/SegerstolpeMatrix.rds"))
SegerMetaData <- readRDS(url("https://storage.googleapis.com/cellid-cbl/SegerstolpeMetaData2.rds"))
data("HgProteinCodingGenes")
BaronMatrixProt <- BaronMatrix[rownames(BaronMatrix) %in% HgProteinCodingGenes,]
SegerMatrixProt <- SegerMatrix[rownames(SegerMatrix) %in% HgProteinCodingGenes,]
# Create Seurat object and remove remove low detection genes
Baron <- CreateSeuratObject(counts = BaronMatrixProt, project = "Baron", min.cells = 5, meta.data = BaronMetaData)
Seger <- CreateSeuratObject(counts = SegerMatrixProt, project = "Segerstolpe", min.cells = 5, meta.data = SegerMetaData)
# Library-size normalization, log-transformation, and centering and scaling of gene expression values
Baron <- NormalizeData(Baron)
Baron <- ScaleData(Baron, features = rownames(Baron))
MCA降维函数RunMCA()源代码分析
# 引入cpp代码
library(Rcpp)
Rcpp::sourceCpp('mca.cpp')
library(stringr)
# 从seurat对象中获得单细胞表达矩阵
data_matrix <- as.matrix(GetAssayData(Baron, "data"))
X = data_matrix
nmcs = 50
features = NULL
reduction.name = "MCA"
# preprocessing matrix
# 预处理表达矩阵
X <- as.matrix(X)
X <- X[rowVars(X) != 0, ]
X <- X[str_length(rownames(X)) > 0, ]
X <- X[!duplicated(rownames(X)), ]
cellsN <- colnames(X)
featuresN <- rownames(X)
# 运行cpp的 MCAStep1 函数
MCAPrepRes <- MCAStep1(X)
# 对由 MCAStep1 函数计算出的相关频率矩阵 S[NK] 进行SVD分解
SVD <- irlba::irlba(A = MCAPrepRes$Z, nv = nmcs + 1, nu = 1)[seq(3)]
# 利用 MCAStep2 函数计算低维空间每个cell和每个基因的坐标
## MCAPrepRes$Z 代表相关频率矩阵 S[NK]，
MCA <- MCAStep2(Z = MCAPrepRes$Z, V = SVD$v[, -1], Dc = MCAPrepRes$Dc)
component <- paste0(reduction.name, "_", seq(ncol(MCA$cellsCoordinates)))
colnames(MCA$cellsCoordinates) <- component
colnames(MCA$featuresCoordinates) <- component
rownames(MCA$cellsCoordinates) <- cellsN
rownames(MCA$featuresCoordinates) <- featuresN
MCA$stdev <- SVD$d[-1]
class(MCA) <- "MCA"
其中cpp脚本如下：
#include <RcppArmadillo.h>
using namespace Rcpp;
using namespace RcppArmadillo;
//[[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
List MCAStep1(NumericMatrix X) {
# 读入preprocessing后的单细胞表达矩阵
arma::mat AM = arma::mat(X.begin(), X.rows(), X.cols(), true);
# 计算 min(Mp)，详细见MCA降维的数学原理
arma::colvec rmin = arma::min(AM,1);
# 计算 max(Mp)，详细见MCA降维的数学原理
arma::colvec rmax = arma::max(AM,1);
# 计算 max(Mp) - min(Mp)，详细见MCA降维的数学原理
arma::colvec range = (rmax -rmin);
# 计算 m[np] - min(Mp)，详细见MCA降维的数学原理
AM.each_col() -= rmin;
# 分子分母相除，详细见MCA降维的数学原理
AM.each_col() /= range;
# 计算 x[-] = 1 - x[+]，并将两个矩阵合并
arma::mat FM = join_cols(AM, 1 - AM);
# 清除变量AM
AM.clear();
# 计算 NF 总的样本数
long total = arma::accu(FM);
# 计算 D[r] 矩阵
arma::rowvec colsum = arma::sum(FM,0);
# 计算 D[c] 矩阵
arma::colvec rowsum = arma::sum(FM,1);
# 计算 D[r]^(-1/2) 矩阵
FM.each_row() /= sqrt(colsum);
# 计算 D[c]^(-1/2) 矩阵
FM.each_col() /= sqrt(rowsum);
arma::colvec Dc = 1/(sqrt(rowsum/total));
# 此时的 “FM” 即为相关频率矩阵 S[NK] ，命名为 Z
return List::create(Named("Z") = wrap(FM),
Named("Dc") = wrap(Dc));
}
//[[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
List MCAStep2(NumericMatrix Z, NumericMatrix V, NumericVector Dc) {
# 读入相应数据
arma::mat AV = arma::mat(V.begin(), V.rows(), V.cols(), true);
arma::mat AZ = arma::mat(Z.begin(), Z.rows(), Z.cols(), true);
arma::colvec ADc = arma::colvec(Dc);
arma::mat FeaturesCoordinates = AZ * AV;
int AZcol = AZ.n_cols;
AZ.clear();
FeaturesCoordinates.each_col() %= ADc;
ADc.clear();
# 计算低维空间cell和基因的坐标矩阵
return List::create(Named("cellsCoordinates") = wrap(std::sqrt(AZcol) * AV),
Named("featuresCoordinates") = wrap(FeaturesCoordinates.head_rows(FeaturesCoordinates.n_rows/2)));
}
2.计算cell与基因距离：
# 进行MCA降维
Baron <- RunMCA(Baron)
# 定义维度
dims = 1:50
# 获得cell在低维空间的坐标矩阵
cells <- rownames(Embeddings(Baron, "mca"))
cellsEmb <- Embeddings(Baron, "mca")[cells, dims]
# 获得gene在低维空间的坐标矩阵
features <- rownames(Loadings(Baron, "mca"))
genesEmb <- Loadings(Baron, "mca")[features, dims]
CellGeneDistance <- pairDist(cellsEmb, genesEmb)
### 备注 pairDist函数
pairDist <- function(x, y) {
z <- fastPDist(y, x)
rownames(z) <- rownames(y)
colnames(z) <- rownames(x)
return(z)
}
其中fastPDist是cpp函数：
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
using namespace arma;
//[[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
NumericMatrix fastPDist(NumericMatrix Ar, NumericMatrix Br) {
int m = Ar.nrow(),
n = Br.nrow(),
k = Ar.ncol();
# 分别读入cell和基因在低维空间的坐标矩阵
arma::mat A = arma::mat(Ar.begin(), m, k, false);
arma::mat B = arma::mat(Br.begin(), n, k, false);
# 计算欧式距离
arma::colvec An = sum(square(A),1);
arma::colvec Bn = sum(square(B),1);
arma::mat C = -2 * (A * B.t());
C.each_col() += An;
C.each_row() += Bn.t();
return wrap(sqrt(C));
}
然后根据每个cell与每个基因在低维空间的坐标来计算欧式距离并进行排序，那么对于每个cell来说，排名靠前的基因为该细胞的特征基因
