相比较SVR而言，这里有另外一种解决单细胞反卷积的方法，nnls（非负最小二乘）
文章链接：《Bulk tissue cell type deconvolution with multi-subject single-cell expression reference》

核心思想

1. X_jp 代表给定tissue的 sample j 中gene g的mRNA分子数
2. X_jpc 代表给定tissue的 sample j ，某一cell c中gene g的mRNA分子数
3. C^k_j 代表cell type k
4. 
   m^k_j 代表cell type k的细胞数

这个式子代表给定tissue的 sample j ，cell type k 中gene g 的平均相对表达量
上式可以改编为

1.

代表cell type k中细胞总mRNA分子数的平均值（除以cell type k的总细胞数记作平均）

2. 

   记作total number of cells in the bulk tissue

3. 

   记作proportion of cells from cell type k

4. 

   记作 bulk tissue 中 gene g 的相对表达量

因此基于上述关系我们容易得知：

这样的正比关系也暗示我们可以用线性模型来衡量这种关系，因此我们的目标就是估计 p^k_j 来表示细胞成分

代码分析

1. 读取数据并运算

library(MuSiC)

GSE50244.bulk.eset = readRDS('C:/Users/lenovo/Downloads/GSE50244bulkeset.rds')
GSE50244.bulk.eset

EMTAB.eset = readRDS('C:/Users/lenovo/Downloads/EMTABesethealthy.rds')
EMTAB.eset


XinT2D.eset = readRDS('C:/Users/lenovo/Downloads/XinT2Deset.rds')
XinT2D.eset


Est.prop.GSE50244 = music_prop(bulk.eset = GSE50244.bulk.eset, sc.eset = EMTAB.eset, clusters = 'cellType',
                               samples = 'sampleID', select.ct = c('alpha', 'beta', 'delta', 'gamma',
                                                                   'acinar', 'ductal'), verbose = F)

对应函数music_prop分步解析：

# 读入数据
bulk.eset = GSE50244.bulk.eset ## 读入bulk 的数据
sc.eset = EMTAB.eset  ## 读入scRNA 的数据
markers = NULL
clusters = 'cellType'
samples = 'sampleID'
select.ct =  c('alpha', 'beta', 'delta',
           'gamma','acinar', 'ductal')   ## cell type 的名称
cell_size = NULL
ct.cov = FALSE
verbose = TRUE
iter.max = 1000
nu = 0.0001
eps = 0.01
centered = FALSE
normalize = FALSE

其中bulk.eset代表的是不同的bulk seq的数据，该例子中一共89个sample：

bulk.eset

本项目的scRNA seq数据如下：
首先是sample的情况

sc.eset sampleID

1097代表的是总共的细胞数目，形如1，2，3，.....这样的属于不同的sample（全为 1 的属于同一个sample，全为 2 的属于另一个sample，以此类推）
其次是cell type的情况：

sc.eset cellTypeID

1097代表的是总共的细胞数目，形如1，2，3，.....这样的属于不同的cell type（全为 1 的属于同一个cell type，全为 2 的属于另一个cell type，以此类推）

接下来对数据进行处理：

# bulk 中选择基因表达不全为0的基因
bulk.gene = rownames(bulk.eset)[rowMeans(exprs(bulk.eset)) != 0]
bulk.eset = bulk.eset[bulk.gene, , drop = FALSE]

# 把这部分的基因定义为scRNA的markers
sc.markers = bulk.gene
# 对scRNA 数据进行标准化，并计算 Sigma 和 S (M.S)
sc.basis = music_basis(sc.eset, non.zero = TRUE, markers = sc.markers, clusters = clusters, samples = samples, select.ct = select.ct, cell_size = cell_size, ct.cov = ct.cov, verbose = verbose)

# sc.basis$Disgn.mtx 代表着单细胞不同 cell type 的基因表达矩阵，取bulk和scRNA基因的交集
cm.gene = intersect(rownames(sc.basis$Disgn.mtx), bulk.gene)
m.sc = match(cm.gene, rownames(sc.basis$Disgn.mtx))
m.bulk = match(cm.gene, bulk.gene)

# 筛选scRNA矩阵的交集基因
D1 = sc.basis$Disgn.mtx[m.sc, ]; 
M.S = colMeans(sc.basis$S, na.rm = T);

# 定义 bulk 表达矩阵的相对表达量
Yjg = relative.ab(exprs(bulk.eset)[m.bulk, ]); N.bulk = ncol(bulk.eset);

Sigma = sc.basis$Sigma[m.sc, ]

# 对scRNA的cell type表达情况进行筛选
valid.ct = (colSums(is.na(Sigma)) == 0)&(colSums(is.na(D1)) == 0)&(!is.na(M.S))
D1 = D1[, valid.ct]
M.S = M.S[valid.ct]
Sigma = Sigma[, valid.ct]
  
# 初始化
Est.prop.allgene = NULL
Est.prop.weighted = NULL
Weight.gene = NULL
r.squared.full = NULL
Var.prop = NULL

# 对bulk 表达矩阵的每个 sample 进行遍历
for(i in 1:N.bulk){
  if(sum(Yjg[, i] == 0) > 0){
    # 筛选非零表达的基因
    D1.temp = D1[Yjg[, i]!=0, ];
    Yjg.temp = Yjg[Yjg[, i]!=0, i];
    Sigma.temp = Sigma[Yjg[,i]!=0, ];
  }else{
    D1.temp = D1;
    Yjg.temp = Yjg[, i];
    Sigma.temp = Sigma;
  }

  # 对每个 sample 的bulk表达矩阵 （向量）与 scRNA 表达矩阵做 nnls
  lm.D1.weighted = music.iter(Yjg.temp, D1.temp, M.S, Sigma.temp, iter.max = iter.max,
                              nu = nu, eps = eps, centered = centered, normalize = normalize)
  Est.prop.allgene = rbind(Est.prop.allgene, lm.D1.weighted$p.nnls)
  Est.prop.weighted = rbind(Est.prop.weighted, lm.D1.weighted$p.weight)
  weight.gene.temp = rep(NA, nrow(Yjg)); weight.gene.temp[Yjg[,i]!=0] = lm.D1.weighted$weight.gene;
  Weight.gene = cbind(Weight.gene, weight.gene.temp)
  r.squared.full = c(r.squared.full, lm.D1.weighted$R.squared)
  Var.prop = rbind(Var.prop, lm.D1.weighted$var.p)
}
colnames(Est.prop.weighted) = colnames(D1)
rownames(Est.prop.weighted) = colnames(Yjg)
colnames(Est.prop.allgene) = colnames(D1)
rownames(Est.prop.allgene) = colnames(Yjg)
names(r.squared.full) = colnames(Yjg)
colnames(Weight.gene) = colnames(Yjg)
rownames(Weight.gene) = cm.gene
colnames(Var.prop) = colnames(D1)
rownames(Var.prop) = colnames(Yjg)

## 其中 Est.prop.weighted 表示细胞组分

在music_basis内置函数中，有两个重要的参数，一个是 Sigma ，另一个是 M.S

1. Sigma

Sigma <- sapply(unique(clusters), function(ct){
   #  apply 计算某个 cell type 在不同 sample 中每个基因平均相对表达量的方差
   apply(
        ##  sapply 计算的是某个 cell type 在不同 sample 中每个基因的平均相对表达量
        sapply(unique(samples), function(sid){
            y = exprs(x)[,clusters %in% ct & samples %in% sid, drop = FALSE]
            rowSums(y)/sum(y)
   }), 1, var, na.rm = TRUE)
 })

# 以其中一个细胞类型 'delta' 为例
test = sapply(unique(samples), function(sid){
   y = exprs(x)[,clusters %in% 'delta' & samples %in% sid, drop = FALSE]
   rowSums(y)/sum(y)
})

test(cell type 'delta')

列 1，2，3...代表不同的sample，行代表不同的基因，这张表说明 细胞类型 'delta' 在不同的sample中的基因平均相对表达量；相对基因表达量的计算为 cell type A 中 属于 sample 1 的所有细胞，所组成的表达矩阵中行元素求和（行为基因）再除以全部元素的求和，rowSums(y)/sum(y)

Sigma[1:20.1:5]

Sigma 代表每个 cell type 在所有 sample 中相对基因表达量均方差

另外

y = exprs(x)[,clusters %in% 'delta' & samples %in% '1', drop = FALSE]

y(sample 1 cell type delta)

y 表示 sample 1 中属于 cell type delta 的细胞表达矩阵，行表示基因，列表示单个的细胞

2. M.S

S <- sapply(unique(clusters), function(ct){
    # 计算 某个 cell type 在所有 sample 中表达谱特征值的均值
    my.rowMeans(
        ## sapply 计算的是某个 cell type 在不同 sample 中表达谱的特征值
       sapply(unique(samples), function(sid){
            y = exprs(x)[, clusters %in% ct & samples %in% sid, drop = FALSE]
            sum(y)/ncol(y)
 }), na.rm = TRUE)
})

# 以其中一个细胞类型 "delta" 为例
test = my.rowMeans(sapply(unique(samples), function(sid){
  y = exprs(x)[, clusters %in% "delta" & samples %in% sid, drop = FALSE]
  sum(y)/ncol(y)
}), na.rm = TRUE)

test(cell type 'delta')

列 1，2，3... 代表不同的sample，里面的值代表个 cell type 在不同 sample 中表达谱的特征值，这个特征值利用的是 cell type A 中 属于 sample 1 的所有细胞，所组成的表达矩阵中全部元素求和再除以细胞个数，sum(y)/ncol(y)
因此，

S

S代表每个 cell type 在每个 sample 中的表达谱特征值，所以 M.S（sc.basis$S 即为上述的S） 代表的每个 cell type 在所有 sample 中表达谱特征值的均值

M.S = colMeans(sc.basis$S, na.rm = T)  # sc.basis$S 即为上述的 S

M.S

另外

y = exprs(x)[,clusters %in% 'delta' & samples %in% '1', drop = FALSE]

y(sample 1 cell type delta)

y 表示 sample 1 中属于 cell type delta 的细胞表达矩阵，行表示基因，列表示单个的细胞

分析music.iter和music.basic（music.iter的内置函数）的源码：

# 对于函数 music.iter
## 读入数据
Y = Yjg.temp   # 某个 bulk 的表达矩阵
D = D1.temp    # scRNA 的表达矩阵
S = M.S
Sigma = Sigma.temp
iter.max = 1000
nu = 0.0001
eps = 0.01
centered = FALSE
normalize = FALSE


k = ncol(D); # number of cell types
# 计算 bulk 和 scRNA 基因的交集，得到相应的表达矩阵
common.gene = intersect(names(Y), rownames(D))
Y = Y[match(common.gene, names(Y))];
D = D[match(common.gene, rownames(D)), ]
Sigma = Sigma[match(common.gene, rownames(Sigma)), ]
X = D
Y = Y*100
## nnls计算
lm.D = music.basic(Y, X, S, Sigma, iter.max = iter.max, nu = nu, eps = eps)

#<-----------------分割线----------------->
# 对于函数music.basic
k = ncol(X)
## nnls 计算某个 bulk (Y 矩阵) 的表达矩阵和 scRNA (X 矩阵) 的表达矩阵
lm.D = nnls(X, Y)
r = resid(lm.D)

## 定义基因的权重，lm.D$x 代表的是 nnls 回归系数
weight.gene = 1/(nu + r^2 + colSums((lm.D$x*S)^2*t(Sigma)))

## 响应变量为 bulk 表达矩阵乘开根号后的weight.gene
Y.weight = Y*sqrt(weight.gene)
## 决策变量为 scRNA 表达矩阵对应元素乘根号后的weight.gene
D.weight = sweep(X, 1, sqrt(weight.gene), '*')

## 再次利用nnls计算
lm.D.weight = nnls(D.weight, Y.weight)
## 细胞组分比例为 p.weight
p.weight = lm.D.weight$x/sum(lm.D.weight$x)
p.weight.iter = p.weight
r = resid(lm.D.weight)

## 循环的目的是不停的迭代更新 weight.gene 使得 胞组分比例 p.weight 收敛
for(iter in 1:iter.max){
  weight.gene = 1/(nu + r^2 + colSums( (lm.D.weight$x*S)^2*t(Sigma)  ))
  Y.weight = Y*sqrt(weight.gene)
  D.weight = X * as.matrix(sqrt(weight.gene))[,rep(1,k)]
  ## 不断更迭scRNA的weight和bulk的weight，使得最终细胞组分 p.weight 达到稳定
  lm.D.weight = nnls(D.weight, Y.weight )
  p.weight.new = lm.D.weight$x/sum(lm.D.weight$x)
  r.new = resid(lm.D.weight)
  if(sum(abs(p.weight.new - p.weight)) < eps){
    p.weight = p.weight.new;
    r = r.new
    R.squared = 1 - var(Y - X%*%as.matrix(lm.D.weight$x))/var(Y)
    fitted = X%*%as.matrix(lm.D.weight$x)
    var.p = diag(solve(t(D.weight)%*%D.weight)) * mean(r^2)/sum(lm.D.weight$x)^2
    return(list(p.nnls = (lm.D$x)/sum(lm.D$x), q.nnls = lm.D$x, fit.nnls = fitted(lm.D), resid.nnls = resid(lm.D),
                p.weight = p.weight, q.weight = lm.D.weight$x, fit.weight = fitted, resid.weight =  Y - X%*%as.matrix(lm.D.weight$x),
                weight.gene = weight.gene, converge = paste0('Converge at ', iter),
                rsd = r, R.squared = R.squared, var.p = var.p));
  }
  p.weight = p.weight.new;
  r = r.new;
}

其中，lm.D$x 代表 nnls 模型的回归系数

由此可见，作者并不是简单的直接运用scRNA表达矩阵和bulk矩阵，利用nnls直接计算细胞组分 p.weight，而是先定义一个 weight.gene来衡量基因的贡献程度，然后再分别得到矫正后的scRNA表达矩阵 D.weight 和校正后的 bulk 表达矩阵 Y.weight （进行表达矩阵数值的矫正），利用矫正后的矩阵进行nnls的计算（bulk 矩阵的 Y.weight 为Y.weight = Y*sqrt(weight.gene)；scRNA 矩阵 D.weight 为D.weight = sweep(X, 1, sqrt(weight.gene), '*')）。通过不断迭代使得 细胞组分 p.weight 收敛

并且思考weight.gene = 1/(nu + r^2 + colSums((lm.D$x*S)^2*t(Sigma)))，其中Sigma和S=M.S分别具备不同的统计学意义。Sigma 越大代表相同 cell type 中该基因在不同sample中的波动越大；S=M.S越大代表该 cell type中所有基因从整体上偏高表达，并且S=M.S代表每个 cell type 在所有 sample 中表达谱特征值的均值，该值作用是衡量某cell type的整体表达水平，而回归系数lm.D$x代表各个cell type在bulk中的成分比例；lm.D$x*S代表了每个细胞类型在bulk中整体基因表达的贡献（结果为一列向量，向量长度与cell type数相同），他们的和就代表了bulk的整体基因表达水平；所以colSums((lm.D$x*S)^2*t(Sigma))即考虑了不同sample之间基因表达的波动，也考虑了各个cell type中基因的整体表达水平；r²越高代表nnls模型拟合越差，反之越差好；因此weight.gene = 1/(nu + r^2 + colSums((lm.D$x*S)^2*t(Sigma)))越大代表该基因对nnls回归的贡献越大，从而不断迭代使得模型收敛