原理

這是一篇發(fā)表GB上題目為《scINSIGHT for interpreting single-cell gene expression from biologically heterogeneous data》 的文章，利用 Integrate NMF 的方法解決單細胞批次整合的問題

越來越多的 scRNA-seq 數(shù)據(jù)強調(diào)需要綜合分析來解釋單細胞樣本之間的異同吨拍。盡管已經(jīng)開發(fā)了不同的批次效應去除方法，但沒有一種方法適用于來自多種生物條件的異質(zhì)單細胞樣品。 我們提出了一種方法 scINSIGHT山叮，用于學習在不同生物條件中常見或特定于不同生物條件的基因表達模式，即以聯(lián)合建模和解釋來自生物異質(zhì)來源的單細胞樣本中的基因表達模式添履。

因此，問題就轉(zhuǎn)換為了矩陣分解的最優(yōu)化問題

Codes

# 準備前期數(shù)據(jù)本讥，準備表達矩陣
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矩陣分解

# 讀入數(shù)據(jù)
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"

其中珊泳，參數(shù)解釋如下：

1. K Number of common gene pathways，解釋為共同基因通路的數(shù)量
2. K_j Number of dataset-specific gene pathways拷沸，解釋為每個sample特有的基因通路的數(shù)量
3. LDA Regularization parameters

# 將每個sample特有的基因通路的數(shù)量存入對象中object
object@parameters[["K_j"]] = K_j

# 對LDA參數(shù)進行排序色查，選擇最小的存入對象object
LDA = sort(LDA)
index_lda = which.min(abs(LDA-0.01))
lda0 = LDA[[index_lda]]
object@parameters[["lda"]] = lda0
  
# 對共同基因通路的數(shù)量進行排序
n_K = length(K)
K = sort(K)

# 提取object中的單細胞表達矩陣
cnt_list = object@norm.data
# 做轉(zhuǎn)置，行為細胞撞芍，列為基因
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分解

# 設置種子數(shù)
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)

# 輸入函數(shù) iNMF_BCD_Increase 的變量解釋如下
## cnt_list 為不同 sample 和 condition 下的表達矩陣
## labels 為各sample對應的 condition秧了，該例子為 0 和 1
## K_j 為每個sample特有的基因通路的數(shù)量
## K 為共同基因通路的數(shù)量
## lda1 為最優(yōu)的LDA參數(shù)
## lda2 為最優(yōu)的LDA參數(shù)

函數(shù) iNMF_BCD_Increase 是利用C++寫的，我們分段來解讀下它們的功能

算法步驟：

1. 初始化 V 序无，Wι1验毡，Wι2和Hj矩陣衡创，隨機的非負矩陣
2. 推斷各個矩陣的數(shù)值，
   
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 統(tǒng)計的是 sample 1 的基因數(shù)量，本例中 N = 500
  int N = X[0].n_cols;
  
  // 隨機生成矩陣 V 挨稿，行數(shù)為 K_1 = 5仇轻，列 N = 500
  mat V = randu<mat>(K_1,N);
  
  //創(chuàng)建數(shù)值型向量 M_l(L)，里面存儲的是 X(L) 每個sample的細胞數(shù)
  NumericVector M_l(L);
  for(int i=0; i<L; i++)
  {
    M_l[i] = X[i].n_rows;
  }
  
  // W_1 行為每個sample細胞數(shù)奶甘，列為 K_l = 2篷店；W_2 行為每個sample細胞數(shù)，列為 K_1 = 5臭家；兩個均為隨機產(chǎn)生
  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);
  }
  
  // 隨機創(chuàng)建 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的細胞數(shù)
    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 進行迭代優(yōu)化iNMF的結(jié)果胰苏，這是迭代的全部代碼，接下來一步一步看

step2 參數(shù)推斷：

估算出 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 倍的迭代次數(shù)輸出一次
      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；每次循環(huán)代表其中一個數(shù)
      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]，只不過加了懲罰系數(shù) 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 依據(jù)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 的結(jié)果
      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"))
      }

   # 獲得較為穩(wěn)定的結(jié)果
      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]]))})
  }

由于在默認參數(shù)下胰丁，共同基因通路的數(shù)量 K 為為 5，7 喂分，9锦庸，11，13蒲祈，15甘萧，因此要篩選出最佳的K值。所以上述代碼的目的是根據(jù)每一個 K值對應各個 seed 的之間聚類的情況進行篩選梆掸，對與每一個 K值幔嗦，每一個seed的結(jié)果應該趨于一致