去除細胞效應和基因效應

2020 年 3 月 27 日
筆記

其實質量控制三部曲，還有一個很關鍵的點沒有講解，就是多個樣本整合，並且區分批次效應和生物學差異。但是這個點很大程度是依賴於經驗，就是說，要想搞清楚，需要寫很多自定義的代碼去一點一滴的探索，而不僅僅是流程。所以我們就只能介紹到這裡，假設大家都拿到了乾淨的表達矩陣，而且可以很肯定的說這個表達矩陣做下游分析是ok的。

那麼我們就來看看，有了乾淨的表達矩陣後下游分析的第一個分析要點就是：歸一化和標準化。

歸一化和標準化的區別

實際上口語裏面通常是沒辦法很便捷的區分這兩個概念，我查閱了大家的資料，發現基本上都是混淆在z-score轉換和0-1轉換這兩個上面。所以我這裡把歸一化和標準化替換成為去除樣本/細胞效應或者去除基因效應：

首先去除樣本/細胞效應：因為不同樣本或者細胞的測序數據量不一樣，所以同樣的一個基因在不同細胞，哪怕你看到的表達量是一樣的，但是背後細胞整體測序數據量的差異其實反而說明了這個基因在不同細胞表達量其實是有差異的。在seurat3裏面的代碼是：

sce <- NormalizeData(sce, normalization.method =  "LogNormalize", scale.factor = 10000)

默認文庫大小是1萬，默認會取log值。
然後去除基因效應，這個主要是在繪製熱圖的時候會需要使用，因為個別基因表達量超級高，在熱圖裏面一枝獨秀，實際上我們並不會關心不同基因的表達量高低，我們僅僅是想看指定基因在不同細胞的高低而已，這樣的話，就把該基因的表達量在不同細胞的數值，進行z-score這樣的轉換。在seurat3裏面的代碼是：

sce <- ScaleData(sce)

這樣處理後的表達矩陣，就可以進行後續的降維聚類分群啦，我們下期再講。

取log讓表達量離散程度降低

這個時候眼尖的朋友其實看到了，在使用NormalizeData函數的時候，有一個 normalization.method = "LogNormalize" 參數被設置了，這個是為什麼呢？

其實很簡單，原始的raw counts矩陣是一個離散型的變量，離散程度很高。有的基因表達丰度比較高，counts數為10000，有些低表達的基因counts可能10，甚至在有些樣本中為0。因為是單細胞轉錄組，drop-out現象很嚴重，其實大量的基因在很多細胞的表達量都是0，如果是10x數據，甚至會出現一個表達矩陣裏面97%的數值都是0 的現象。

如果對表達量取一下log10，發現10000變成了4，10變成了1，這樣之前離散程度很大的數據就被集中了。

有時當表達量為0時，取log會出現錯誤，可以log(counts+1)來取log值。當x=1時，所有的log系列函數值都為0。這樣原本表達量為0的值，取log後仍為0。

這也就是UCSC的XENA下載到的表達矩陣的形式。

使用GetAssayData函數查看不同形式的表達矩陣

在seurat3裏面，很方便進行各種形式的歸一化或者標準化，同樣的，也很容易查看處理前後的表達矩陣，使用GetAssayData函數即可。

This function can be used to pull information from any of the slots in the Assay class. For example, pull one of the data matrices("counts", "data", or "scale.data").

如下代碼：

# 最原始數據  GetAssayData(sce,'counts')[1:6,1:6]  # 去除了細胞測序數據量後  1/(colSums(as.data.frame(GetAssayData(sce,'counts')))[2]/10000)  GetAssayData(sce,'data')[1:6,1:6]  # z-score後  GetAssayData(sce,'scale.data')[1:6,1:6]

如下結果：

> # 最原始數據  > GetAssayData(sce,'counts')[1:6,1:6]  6 x 6 sparse Matrix of class "dgCMatrix"               6A-11 6A-13 6A-14 6A-15 6A-16 6A-17  DDX11L1          .     .     .     .     .     .  WASH7P           .     1     .     .     .     .  MIR6859-1        .     .     .     .     .     .  RP11-34P13.3     .     .     .     .     .     .  OR4G11P          .     .     .     .     .     .  OR4F5            .     .     .     .     .     .  > # 去除了細胞測序數據量後  > 1/(colSums(as.data.frame(GetAssayData(sce,'counts')))[2]/10000)        6A-13  0.009725971  > GetAssayData(sce,'data')[1:6,1:6]  6 x 6 sparse Matrix of class "dgCMatrix"               6A-11       6A-13 6A-14 6A-15 6A-16 6A-17  DDX11L1          . .               .     .     .     .  WASH7P           . 0.009678978     .     .     .     .  MIR6859-1        . .               .     .     .     .  RP11-34P13.3     . .               .     .     .     .  OR4G11P          . .               .     .     .     .  OR4F5            . .               .     .     .     .  > # z-score後  > GetAssayData(sce,'scale.data')[1:6,1:6]                    6A-11      6A-13      6A-14      6A-15      6A-16      6A-17  RP5-857K21.9 -0.8302953 -0.8302953 -0.7327066 -0.7224892 -0.7683106 -0.7692732  RP5-857K21.8 -0.5066862 -0.5688591 -0.5349974 -0.4881690 -0.3650311 -0.4524682  PRKCZ        -0.8350477 -0.5996852 -0.8350477 -0.8350477 -0.5988844 -0.8350477  TNFRSF14     -0.2550752 -0.2550752 -0.2550752 -0.2550752 -0.2550752 -0.2550752  TP73         -0.2761133 -0.2761133 -0.2761133 -0.2761133 -0.2761133 -0.2337521  HES2         -0.4569104 -0.4569104 -0.4569104 -0.4569104 -0.4569104

其實樣本/細胞效應不僅僅是文庫大小

每個細胞測序數據量的不一致是很容易理解的，但其實細胞之間還有很多其它效應，比如線粒體基因含量，ERCC含量等等，那些處理起來，其實就是深入了解我們講解seurat裏面的NormalizeData和ScaleData函數。