原理

这是一篇发表GB上题目为《scINSIGHT for interpreting single-cell gene expression from biologically heterogeneous data》 的文章，利用 Integrate NMF 的方法解决单细胞批次整合的问题

越来越多的 scRNA-seq 数据强调需要综合分析来解释单细胞样本之间的异同。尽管已经开发了不同的批次效应去除方法，但没有一种方法适用于来自多种生物条件的异质单细胞样品。 我们提出了一种方法 scINSIGHT，用于学习在不同生物条件中常见或特定于不同生物条件的基因表达模式，即以联合建模和解释来自生物异质来源的单细胞样本中的基因表达模式。

因此，问题就转换为了矩阵分解的最优化问题

Codes

# 准备前期数据，准备表达矩阵
S1 <- matrix(runif(50000,0,2), 500,100)
S2 <- matrix(runif(60000,0,2), 500,120)
S3 <- matrix(runif(80000,0,2), 500,160)
S4 <- matrix(runif(75000,0,2), 500,150)
data = list(S1, S2, S3, S4)
# 为表达矩阵命名 
sample = c("sample1", "sample2", "sample3", "sample4")
# 设置不同的处理条件
condition = c("control", "activation", "control", "activation")
names(data) = sample
names(condition) = sample

# 将data
scINSIGHTx <- create_scINSIGHT(data, condition)

进行iNMF矩阵分解

# 读入数据
object=scINSIGHTx
K = seq(5,15,2)
K_j = 2
LDA = c(0.001, 0.01, 0.1, 1, 10)
thre.niter = 500
thre.delta = 0.01
num.cores = 1
B = 5
out.dir = NULL
method = "increase"

其中，参数解释如下：

1. K Number of common gene pathways，解释为共同基因通路的数量
2. K_j Number of dataset-specific gene pathways，解释为每个sample特有的基因通路的数量
3. LDA Regularization parameters

# 将每个sample特有的基因通路的数量存入对象中object
object@parameters[["K_j"]] = K_j

# 对LDA参数进行排序，选择最小的存入对象object
LDA = sort(LDA)
index_lda = which.min(abs(LDA-0.01))
lda0 = LDA[[index_lda]]
object@parameters[["lda"]] = lda0
  
# 对共同基因通路的数量进行排序
n_K = length(K)
K = sort(K)

# 提取object中的单细胞表达矩阵
cnt_list = object@norm.data
# 做转置，行为细胞，列为基因
cnt_list = lapply(cnt_list, function(x){t(x)})

# 设置标签
uLabels = unique(object@condition)
labels = sapply(object@condition, function(condition){which(uLabels==condition)-1})

其中 labels 表示如下：

labels

然后对表达矩阵进行iNMF分解

# 设置种子数
seeds = 1:B
  
# 读入C++脚本
Rcpp::sourceCpp('E:/iNMF_BCD_Increase.cpp')
  
# Run iNMF BCD
# 以increase为例子
if(method == "increase")
{
    res_all = mclapply(seeds, function(seed){
    res = iNMF_BCD_Increase(cnt_list, labels, K_j, K = K, lda1 = lda0, lda2 = lda0, 
                            eps = 1e-5, hre.niter, thre.delta, seed, TRUE)
    gc()
    return(res)
}, mc.cores = num.cores)

# 输入函数 iNMF_BCD_Increase 的变量解释如下
## cnt_list 为不同 sample 和 condition 下的表达矩阵
## labels 为各sample对应的 condition，该例子为 0 和 1
## K_j 为每个sample特有的基因通路的数量
## K 为共同基因通路的数量
## lda1 为最优的LDA参数
## lda2 为最优的LDA参数

函数 iNMF_BCD_Increase 是利用C++写的，我们分段来解读下它们的功能

算法步骤：

1. 初始化 V ，Wι1，Wι2和Hj矩阵，随机的非负矩阵
2. 推断各个矩阵的数值，
   
3. 计算 loss function

step 1 初始化相关矩阵

  Rcpp::Environment baseEnv("package:base");
  Rcpp::Function setSeed = baseEnv["set.seed"];
  setSeed(seed);

  // Initialize output
  // 这里的 count_list 对应不同sample的单细胞表达矩阵 cnt_list ，这里有四个sample，因此 L = 4
  int L = count_list.size();
  // 这里Label用的是0，1表示，因此max(Label) = 1，所以 J = 2
  int J = max(Label)+1;
  
  // 该例子 K 为 5，7 ，9，11，13，15；n_K = 6；K_1 = 5
  int n_K = K.size();
  K.sort();
  int K_1 = K[0];
  
  // 建立四个向量X(L), W_2(L), W_1(L), H(J)
  std::vector<arma::mat> X(L), W_2(L), W_1(L), H(J);

  // 初始化X(L)，X(L)里面存储着各个sample的单细胞表达矩阵
  for(int i=0;i<L;i++)
  {
    mat temp = count_list[i];
    X[i] = temp;
  }
  
  // N 统计的是 sample 1 的基因数量，本例中 N = 500
  int N = X[0].n_cols;
  
  // 随机生成矩阵 V ，行数为 K_1 = 5，列 N = 500
  mat V = randu<mat>(K_1,N);
  
  //创建数值型向量 M_l(L)，里面存储的是 X(L) 每个sample的细胞数
  NumericVector M_l(L);
  for(int i=0; i<L; i++)
  {
    M_l[i] = X[i].n_rows;
  }
  
  // W_1 行为每个sample细胞数，列为 K_l = 2；W_2 行为每个sample细胞数，列为 K_1 = 5；两个均为随机产生
  for(int i=0; i<L; i++)
  {
    W_2[i] = randu<mat>(M_l[i], K_1);
    W_1[i] = randu<mat>(M_l[i], K_l);
  }
  
  // 随机创建 H[i] 矩阵，行为 K_l = 2，列为 N = 500，这里 J = 2，表示有两种 condition
  for(int i=0; i<J; i++)
  {
    H[i] = randu<mat>(K_l, N);
   // H[i] 的行代表每个sample的细胞数
    for(int k=0; k<K_l; k++)
    {
      double n_H = sqrt(sum(H[i].row(k)%H[i].row(k)));
      H[i].row(k) /= n_H;
      for(int m=0; m<L; m++)
      {
        if(Label[m]==i)
        {
          W_1[m].col(k) *= n_H;
        }
      }
    }
  }

接下来就是定义loss function 进行迭代优化iNMF的结果，这是迭代的全部代码，接下来一步一步看

step2 参数推断：

估算出 V 矩阵的值

// Iterations
    double Loss_1 = 0.0;
    double Loss_2 = 1000000.0;

    start = clock(); // Record the start time of iterations
    int iter; // Record the number of iterations

    for(iter=0; iter<thre_niter; iter++)
    {
      // 每 20 倍的迭代次数输出一次
      if(iter % 20 == 0)
      {
        Rcpp::Rcout<<iter<<std::endl;
      }

      // 定义残差矩阵 Res_v(L)
      std::vector<arma::mat> Res_v(L);
      for(int i=0; i<L; i++)
      {
        // X[i] = W_1[i]*H[Label[i]] + W_2[i]*V + E[i]
        // 那么 Res_v[i] = W_2[i]*V + E[i]
        Res_v[i] = X[i] - W_1[i]*H[Label[i]];
      }
  
     // K_2 为 5，7 ，9，11，13，15；每次循环代表其中一个数
      for(int k=0; k<K_2; k++)
      {
        // 初始化行向量 V_a 为零向量，长度为 N = 500
        rowvec V_a = zeros<rowvec>(N);
        double V_b = 0;
        for(int i=0;i<L;i++)
        {
          V_a += W_2[i].col(k).t()*(Res_v[i]-W_2[i]*V)/M_l[i];
          V_b += sum(W_2[i].col(k)%W_2[i].col(k))/M_l[i];
        }
        // 对 V 矩阵的每一列进行迭代
        V.row(k) = max(V.row(k) + V_a/V_b, eps*ones<rowvec>(N));
      }

其中 V_a 代表算式为

作为分子部分；
V_b 代表的算式是

作为分母部分
对于最终的迭代算式为：

V.row(k) = max(V.row(k) + V_a/V_b, eps*ones<rowvec>(N));

由于V矩阵的初始值为 0，因此 V 矩阵每一列的迭代为 0 + V_a/V_b

估算出 W_1 和 W_2 矩阵的值

     // 这里 loop = TRUE，因此只看 TRUE 部分的
      if(loop) // Whether to solve W_2 and W_1 at the same time
      {
        for(int i=0; i<L; i++)
        {
           // X[i] = W_1[i]*H[Label[i]] + W_2[i]*V + E[i]
           //那么 Res_2 = W_2[i]*V + E[i] 
          mat Res_2 = X[i] - W_1[i]*H[Label[i]];
        
          for(int k=0; k<K_2; k++)
          {
             // Res_2 - W_2[i]*V  代表残差矩阵  E[i]
             // 对 W_2 矩阵的列进行遍历，并更新 W_2[i] 的列
            W_2[i].col(k) = max(W_2[i].col(k) + (Res_2 - W_2[i]*V)*(V.row(k).t())/sum(V.row(k)%V.row(k)), eps*ones<colvec>(M_l[i]));
          }
          // Res_h 代表 W_1[i]*H[Label[i]] + E[i]
          mat Res_h = X[i] - W_2[i]*V;
          for(int k=0; k<K_l; k++)
          {
            // Res_h - (1 + lda1)*W_1[i]*H[Label[i]]  代表残差矩阵  E[i]，只不过加了惩罚系数 lda1
            // 对 W_1 矩阵的列进行遍历，并更新 W_1[i] 的列
            W_1[i].col(k) = max(W_1[i].col(k) + ((Res_h - (1 + lda1)*W_1[i]*H[Label[i]])*(H[Label[i]].row(k).t()))/sum((1+lda1)*H[Label[i]].row(k)%H[Label[i]].row(k)),eps*ones<colvec>(M_l[i]));
          }
        }
      }

对于W_2矩阵：
其中，(Res_2 - W_2[i]*V)*(V.row(k).t()) 代表算式

作为分子；
sum(V.row(k)%V.row(k)) 代表算式

作为分母
对于最终的W_2的算式：

W_2[i].col(k) = max(W_2[i].col(k) + (Res_2 - W_2[i]*V)*(V.row(k).t())/sum(V.row(k)%V.row(k)), eps*ones<colvec>(M_l[i]));

其中W_2[i].col(k)代表之前随机初始化的W_2矩阵

对于W_1矩阵：
其中，(Res_h - (1 + lda1)*W_1[i]*H[Label[i]])*(H[Label[i]].row(k).t()) 代表算式

作为分子；
sum((1+lda1)*H[Label[i]].row(k)%H[Label[i]].row(k) 代表算式

作为分母，lda1 代表 λ1
对于最终的W_1的算式：

W_1[i].col(k) = max(W_1[i].col(k) + ((Res_h - (1 + lda1)*W_1[i]*H[Label[i]])*(H[Label[i]].row(k).t()))/sum((1+lda1)*H[Label[i]].row(k)%H[Label[i]].row(k)),eps*ones<colvec>(M_l[i]));

其中W_1[i].col(k)代表之前随机初始化的W_1矩阵

估算出 H 矩阵的值

      // update H
      for(int i=0; i<J; i++)
      {
        for(int k=0; k<K_l; k++)
        {
          rowvec temp_H = zeros<rowvec>(N);
          double temp_W = 0.0;
          for(int m=0; m<L; m++)
          {
            if(Label[m]==i)
            {  
              // 计算 H 矩阵分子的前半部分
              temp_H += (W_1[m].col(k).t()*(X[m]-W_2[m]*V-(1+lda1)*W_1[m]*H[i]))/M_l[m];
              // 计算 H 矩阵分母部分
              temp_W += (1+lda1)*sum(W_1[m].col(k)%W_1[m].col(k))/M_l[m];
            }
          }
          rowvec s_H = zeros<rowvec>(N);
          for(int j=0; j<J; j++)
          {
            if(j==i)
            {
              continue;
            }
            for(int t=0; t<K_l; t++)
            {
              // 计算 H 矩阵分子的后半部分
              s_H += H[j].row(t);
            }
          }

          H[i].row(k) = max(H[i].row(k) + (temp_H - (lda2/4)*s_H)/temp_W, eps*ones<rowvec>(N));

          double n_H = sqrt(sum(H[i].row(k)%H[i].row(k)));
          H[i].row(k) /= n_H;
          for(int m=0; m<L; m++)
          {
            if(Label[m]==i)
            {
              W_1[m].col(k) *= n_H;
            }
          }
        }
      }

其中 H 矩阵的最终算式为

H[i].row(k) = max(H[i].row(k) + (temp_H - (lda2/4)*s_H)/temp_W, eps*ones<rowvec>(N));
double n_H = sqrt(sum(H[i].row(k)%H[i].row(k)));
H[i].row(k) /= n_H;

(temp_H - (lda2/4)*s_H) 代表算式

作为分子，而temp_H计算方式如下：

for(int m=0; m<L; m++)
{
   if(Label[m]==i)
   {
      // temp_H + 表示自加
      temp_H += (W_1[m].col(k).t()*(X[m]-W_2[m]*V-(1+lda1)*W_1[m]*H[i]))/M_l[m];
   }
}

算式为：

lda1 代表 λ1
而(lda2/4)*s_H)的算式为：

其中lda2代表 λ2
而s_H的计算方式如下：

for(int j=0; j<J; j++)
{
   if(j==i)
   {
       continue;
   }
   for(int t=0; t<K_l; t++)
   {
      s_H += H[j].row(t);
   }
}

算式为

temp_W代表算式

作为分子
而temp_W的计算方式如下：

for(int m=0; m<L; m++)
{
   if(Label[m]==i)
   {
      // temp_W + 表示自加
       temp_W += (1+lda1)*sum(W_1[m].col(k)%W_1[m].col(k))/M_l[m];
   }
}

其中lda1 代表 λ

step 3 计算 loss 和矫正

      // Calculate the value of Objective function
      Loss_1 = Loss_2;
      Loss_2 = 0;

      for(int i=0;i<L;i++)
      {
        // 计算迭代后的 loss，即 X[i]-W_1[i]*H[Label[i]]-W_2[i]*V = E[i]
        Loss_2 += pow(norm(X[i]-W_1[i]*H[Label[i]]-W_2[i]*V,"fro"),2)/M_l[i];
        Loss_2 += lda1*pow(norm(W_1[i]*H[Label[i]],"fro"),2)/M_l[i];
      }


      for(int i=0; i<J; i++)
      {
        for(int j=0;j<J; j++)
        {
          if(j==i)
          {
            continue;
          }

          Loss_2 += lda2/2*accu(H[i]*(H[j].t()));
        }
      }


      if((Loss_1-Loss_2)<thre_delta)
      {
        break;
      }

    }

    end = clock(); // Record the end time of iterations

    // temporary variables
    std::vector<arma::mat> oW_2(L), oW_1(L), oH(J);
    mat oV = V;


    for(int i=0; i<L; i++)
    {
      // 将 W_1 和 W_2 矩阵传给 oW_1 和 oW_2
      oW_2[i] = W_2[i];
      oW_1[i] = W_1[i];
    }


    for(int j=0; j<J; j++)
    {
      // 将 H 矩阵传给 oH
      oH[j] = H[j];
    }

    // Convert values that less than eps to zero
    oV.elem(find(oV<=eps)).zeros();

    for(int l=0;l<L;l++)
    {
      // 将矩阵 oW_1 和 oW_2 中元素小于eps改为 0
      oW_1[l].elem(find(oW_1[l]<=eps)).zeros();
      oW_2[l].elem(find(oW_2[l]<=eps)).zeros();
    }

    for(int l=0;l<J;l++)
    {
      // 将矩阵 oH 中元素小于eps改为 0
      oH[l].elem(find(oH[l]<=eps)).zeros();
    }


    // Normalize V
    // 对 V 矩阵进行标准化
    for(int i=0;i<K_2;i++)
    {
      double n_V = sqrt(sum(oV.row(i)%oV.row(i)));
      oV.row(i) /= n_V;
      for(int l=0;l<L;l++)
      {
        oW_2[l].col(i) *= n_V;
      }
    }



    // Convert the results to several Lists (Can not directly assign the array of mat into List)
    // It means that it is illegal to write down 'Named("W_1") = W_1', W_1 needs to be converted to a List W_1r
    // 将处理好的所有矩阵进行最终的赋值
    List W_1r(L),W_2r(L),H_r(J);
    for(int i=0;i<L;i++)
    {
      W_1r[i] = oW_1[i];
      W_2r[i] = oW_2[i];
    }

    for(int j=0;j<J;j++)
    {
      H_r[j] = oH[j];
    }


    output[num] = List::create(Named("W_1") = W_1r,
                               Named("W_2") = W_2r,
                               Named("H")  =  H_r,
                               Named("V")  =  oV,
                               Named("iters") = iter,
                               Named("loss")  =  Loss_2,
                               Named("times") = (double)(end-start)/CLOCKS_PER_SEC);
  }

step 4 依据stability进行筛选

# C++ 进行iNMF分解
if(method == "increase")
  {
    res_all = mclapply(seeds, function(seed){
      res = iNMF_BCD_Increase(cnt_list, labels, K_j, K = K, lda1 = lda0, lda2 = lda0, eps = 1e-5,
                              thre.niter, thre.delta, seed, TRUE)
      gc()
      return(res)
    }, mc.cores = num.cores)
    res_parallel = list()

    # 该例子 K 为 5，7 ，9，11，13，15；n_K = 6
    for(i in 1:n_K){
      res_name = paste0("K_",as.character(K[[i]]))
     ## 获得每一个 seed 的结果
      res_parallel[[res_name]] = lapply(seeds, function(seed){res_all[[seed]][[i]]})
      names(res_parallel[[res_name]]) = sapply(seeds, function(seed){paste0("seed_",seed)})
      if(!is.null(out.dir)){
        saveRDS(res_parallel[[res_name]], file = paste0(out.dir, "res-k", K[[i]], "-lda", lda0, ".rds"))
      }

   # 获得较为稳定的结果
      object@parameters[["stability"]][i] = get_stability_StrictAndWeight_k(res_parallel[[res_name]], nk = 20)
    }
    names(object@parameters[["stability"]]) = sapply(1:n_K,function(i){paste0("K_",as.character(K[[i]]))})
  }

由于在默认参数下，共同基因通路的数量 K 为为 5，7 ，9，11，13，15，因此要筛选出最佳的K值。所以上述代码的目的是根据每一个 K值对应各个 seed 的之间聚类的情况进行筛选，对与每一个 K值，每一个seed的结果应该趋于一致