原理

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

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

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

Codes

# 準(zhǔn)備前期數(shù)據(jù)，準(zhǔn)備表達(dá)矩陣
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)
# 為表達(dá)矩陣命名 
sample = c("sample1", "sample2", "sample3", "sample4")
# 設(shè)置不同的處理?xiàng)l件
condition = c("control", "activation", "control", "activation")
names(data) = sample
names(condition) = sample

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

進(jìn)行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，解釋為每個(gè)sample特有的基因通路的數(shù)量
3. LDA Regularization parameters

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

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

# 提取object中的單細(xì)胞表達(dá)矩陣
cnt_list = object@norm.data
# 做轉(zhuǎn)置，行為細(xì)胞，列為基因
cnt_list = lapply(cnt_list, function(x){t(x)})

# 設(shè)置標(biāo)簽
uLabels = unique(object@condition)
labels = sapply(object@condition, function(condition){which(uLabels==condition)-1})

其中 labels 表示如下：

labels

然后對(duì)表達(dá)矩陣進(jìn)行iNMF分解

# 設(shè)置種子數(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 下的表達(dá)矩陣
## labels 為各sample對(duì)應(yīng)的 condition，該例子為 0 和 1
## K_j 為每個(gè)sample特有的基因通路的數(shù)量
## K 為共同基因通路的數(shù)量
## lda1 為最優(yōu)的LDA參數(shù)
## lda2 為最優(yōu)的LDA參數(shù)

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

算法步驟：

1. 初始化 V ，Wι1，Wι2和Hj矩陣，隨機(jī)的非負(fù)矩陣
2. 推斷各個(gè)矩陣的數(shù)值，
   
3. 計(jì)算 loss function

step 1 初始化相關(guān)矩陣

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

  // Initialize output
  // 這里的 count_list 對(duì)應(yīng)不同sample的單細(xì)胞表達(dá)矩陣 cnt_list ，這里有四個(gè)sample，因此 L = 4
  int L = count_list.size();
  // 這里L(fēng)abel用的是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];
  
  // 建立四個(gè)向量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)里面存儲(chǔ)著各個(gè)sample的單細(xì)胞表達(dá)矩陣
  for(int i=0;i<L;i++)
  {
    mat temp = count_list[i];
    X[i] = temp;
  }
  
  // N 統(tǒng)計(jì)的是 sample 1 的基因數(shù)量，本例中 N = 500
  int N = X[0].n_cols;
  
  // 隨機(jī)生成矩陣 V ，行數(shù)為 K_1 = 5，列 N = 500
  mat V = randu<mat>(K_1,N);
  
  //創(chuàng)建數(shù)值型向量 M_l(L)，里面存儲(chǔ)的是 X(L) 每個(gè)sample的細(xì)胞數(shù)
  NumericVector M_l(L);
  for(int i=0; i<L; i++)
  {
    M_l[i] = X[i].n_rows;
  }
  
  // W_1 行為每個(gè)sample細(xì)胞數(shù)，列為 K_l = 2；W_2 行為每個(gè)sample細(xì)胞數(shù)，列為 K_1 = 5；兩個(gè)均為隨機(jī)產(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);
  }
  
  // 隨機(jī)創(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] 的行代表每個(gè)sample的細(xì)胞數(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;
        }
      }
    }
  }

接下來(lái)就是定義loss function 進(jìn)行迭代優(yōu)化iNMF的結(jié)果，這是迭代的全部代碼，接下來(lái)一步一步看

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)代表其中一個(gè)數(shù)
      for(int k=0; k<K_2; k++)
      {
        // 初始化行向量 V_a 為零向量，長(zhǎng)度為 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];
        }
        // 對(duì) V 矩陣的每一列進(jìn)行迭代
        V.row(k) = max(V.row(k) + V_a/V_b, eps*ones<rowvec>(N));
      }

其中 V_a 代表算式為

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

作為分母部分
對(duì)于最終的迭代算式為：

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]
             // 對(duì) W_2 矩陣的列進(jìn)行遍歷，并更新 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]，只不過(guò)加了懲罰系數(shù) lda1
            // 對(duì) W_1 矩陣的列進(jìn)行遍歷，并更新 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]));
          }
        }
      }

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

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

作為分母
對(duì)于最終的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)代表之前隨機(jī)初始化的W_2矩陣

對(duì)于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
對(duì)于最終的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)代表之前隨機(jī)初始化的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)
            {  
              // 計(jì)算 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];
              // 計(jì)算 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++)
            {
              // 計(jì)算 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計(jì)算方式如下：

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的計(jì)算方式如下：

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的計(jì)算方式如下：

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 計(jì)算 loss 和矯正

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

      for(int i=0;i<L;i++)
      {
        // 計(jì)算迭代后的 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
    // 對(duì) V 矩陣進(jìn)行標(biāo)準(zhǔn)化
    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
    // 將處理好的所有矩陣進(jìn)行最終的賦值
    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進(jìn)行篩選

# C++ 進(jìn)行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]]))
     ## 獲得每一個(gè) 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]]))})
  }

由于在默認(rèn)參數(shù)下，共同基因通路的數(shù)量 K 為為 5，7 ，9，11，13，15，因此要篩選出最佳的K值。所以上述代碼的目的是根據(jù)每一個(gè) K值對(duì)應(yīng)各個(gè) seed 的之間聚類(lèi)的情況進(jìn)行篩選，對(duì)與每一個(gè) K值，每一個(gè)seed的結(jié)果應(yīng)該趨于一致