實驗記錄3：用R包Seurat進行QC、PCA分析與t-SNE聚類

版本信息：

Seurat v2.0不是3.0！現(xiàn)在Seurat更新了3.0版本，下載也是默認的3.0，這篇記錄只適用于用2.0的。

梗概

  1. 將Cellranger中的基因表達矩陣filtered_gene_bc_matrices用于分析。
  2. 進行質量控制（QC），以刪除異常細胞；
  3. 標準化與歸一化，消除技術噪音與批次效應；
  4. 主成分分析（PCA）與挑選
  5. t-SNE聚類

參考網(wǎng)站：https://satijalab.org/seurat/pbmc3k_tutorial.html
（注意！?。‖F(xiàn)在這個網(wǎng)站會自動跳轉到3.0版本）
Seurat的安裝：R中運行install.packages("Seurat")

上次結果：

經過Cellranger的數(shù)據(jù)整理之后，得到：

  • Filtered gene-barcode matrices MEX: /outs/filtered_gene_bc_matrices
    此輸出結果應為基因-細胞的表達矩陣，用Seurat包進行后續(xù)分析。

Seurat是一種R包，設計用于QC，分析和探索單細胞RNA-seq數(shù)據(jù)。 Seurat旨在使用戶能夠從單細胞轉錄組測量中識別和解釋異質性來源，并整合不同類型的單細胞數(shù)據(jù)。

運行R，并且加載這兩個包

library(Seurat)
library(dplyr)

讀取數(shù)據(jù)

spleen.data <- Read10X(data.dir = '/GRCh38/')

dim(spleen.data)
[1] 33694  1960

原始數(shù)據(jù)的基因數(shù)為33694，細胞數(shù)為1960.

比較普通與疏松矩陣的內存使用：

> dense.size <- object.size(x = as.matrix(x = spleen.data))
> dense.size
530488272 bytes

#轉化為疏松矩陣，查看大小
> sparse.size <- object.size(x = spleen.data)
> sparse.size
45955656 bytes

> dense.size/sparse.size
11.5 bytes

初始化Seurat對象：
命令CreateSeuratObject
輸入數(shù)據(jù)spleen.data
留下所有在>=3個細胞中表達的基因min.cells = 3；
留下所有檢測到>=200個基因的細胞min.genes = 200。
(為了除去一些質量差的細胞)

spleen <- CreateSeuratObject(raw.data = spleen.data, min.cells = 3, min.genes = 200, project = "10X_spleen")

spleen
An object of class seurat in project 10X_spleen 
15655 genes across 1959 samples.

剩下15655 基因和 1959 個細胞

質量控制

以下步驟包括Seurat中scRNA-seq數(shù)據(jù)的標準預處理工作流程。這些代表了Seurat對象的創(chuàng)建，基于QC指標的細胞選擇和過濾，數(shù)據(jù)標準化和縮放，以及高度可變基因的檢測。

mito.genes <- grep(pattern = "^MT-", x = rownames(x = spleen@data), value = TRUE)
percent.mito <- Matrix::colSums(spleen@raw.data[mito.genes, ])/Matrix::colSums(spleen@raw.data)
spleen <- AddMetaData(object = spleen, metadata = percent.mito, col.name = "percent.mito")
VlnPlot(object = spleen, features.plot = c("nGene", "nUMI", "percent.mito"), nCol = 3)
VlnPlot_of_spleen.png
> par(mfrow = c(1, 2))
> GenePlot(object = spleen, gene1 = "nUMI", gene2 = "percent.mito")
> GenePlot(object = spleen, gene1 = "nUMI", gene2 = "nGene")
GenePlot_of_spleen.png

過濾細胞，根據(jù)上面的兩幅圖，去除異常值，這里選擇基因數(shù)從300-5000，線粒體基因占比大于0.1的細胞。（主要看小提琴圖1和圖3）

spleen <- FilterCells(spleen, subset.names = c("nGene", "percent.mito"), low.thresholds = c(300, -Inf), high.thresholds = c(5000,0.10))

查看過濾掉剩下多少細胞：

spleen
An object of class seurat in project 10X_spleen 
 15655 genes across 1940 samples.

剩下15655個基因，1940個細胞。

數(shù)據(jù)標準化

加個log：

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

Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

spleen <- FindVariableGenes(object = spleen, mean.function = ExpMean, dispersion.function = LogVMR, x.low.cutoff = 0.0125, x.high.cutoff = 3, y.cutoff = 0.5)

Calculating gene means
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating gene variance to mean ratios
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
TEXT_SHOW_BACKTRACE environmental variable.
> length(x=spleen@var.genes)
[1] 1829
高度變異基因.png

縮放數(shù)據(jù)并刪除不需要的變體來源

您的單細胞數(shù)據(jù)集可能包含“不感興趣”的變異來源。這不僅包括技術噪音，還包括批次效應，甚至包括生物變異來源（細胞周期階段）。正如(Buettner, et al NBT，2015)中所建議的那樣，從分析中回歸這些信號可以改善下游維數(shù)減少和聚類。為了減輕這些信號的影響，Seurat構建線性模型以基于用戶定義的變量預測基因表達。這些模型的縮放得分殘差存儲在Scale.data槽中，用于降維和聚類。

我們可以消除由批次（如果適用）驅動的基因表達中的細胞 - 細胞變異，細胞比對率（由Drop-seq數(shù)據(jù)的Drop-seq工具提供），檢測到的分子數(shù)量和線粒體基因表達。對于循環(huán)細胞，我們還可以學習“細胞周期”評分（參見此處的示例）并對其進行回歸。在這個有絲分裂后血細胞的簡單例子中，我們回歸了每個細胞檢測到的分子數(shù)量以及線粒體基因含量百分比。

spleen <-ScaleData(spleen, vars.to.regress = c("nUMI","percent.mito"))

Regressing out: nUMI, percent.mito
  |=========================================================================================| 100%
Time Elapsed:  18.0711550712585 secs
Scaling data matrix
  |=========================================================================================| 100%

PCA分析

主成分分析是什么？

主成分分析，是考察多個變量間相關性一種多元統(tǒng)計方法，研究如何通過少數(shù)幾個主成分來揭示多個變量間的內部結構，即從原始變量中導出少數(shù)幾個主成分，使它們盡可能多地保留原始變量的信息，且彼此間互不相關.通常數(shù)學上的處理就是將原來P個指標作線性組合，作為新的綜合指標。

將數(shù)據(jù)集降維，利用低階的變量去反應整體的結果。

spleen <- RunPCA(spleen, pc.genes = spleen@var.genes, do.print = TRUE, pcs.print = 1:5, genes.print = 5)

[1] "PC1"
[1] "CD69"  "CD79A" "TRAC"  "CD3D"  "MS4A1"
[1] ""
[1] "FCN1"          "LYZ"           "SERPINA1"      "CSTA"          "RP11-1143G9.4"
[1] ""
[1] ""
[1] "PC2"
[1] "CD79A"    "MS4A1"    "VPREB3"   "CD79B"    "HLA-DQB1"
[1] ""
[1] "NKG7" "CST7" "GZMA" "CD7"  "CCL5"
[1] ""
[1] ""
[1] "PC3"
[1] "TRDC"  "KLRF1" "MS4A1" "CD79B" "IRF8" 
[1] ""
[1] "IL7R" "TRAC" "CD3D" "CD2"  "CD3G"
[1] ""
[1] ""
[1] "PC4"
[1] "GIMAP7" "GZMB"   "FGFBP2" "SPON2"  "PRF1"  
[1] ""
[1] "BAG3"    "HSPD1"   "FKBP4"   "DNAJA1"  "ZFAND2A"
[1] ""
[1] ""
[1] "PC5"
[1] "UBE2C" "TYMS"  "MKI67" "TOP2A" "AURKB"
[1] ""
[1] "FCGR3A" "FGFBP2" "SPON2"  "GNLY"   "GZMB"  
[1] ""
[1] ""

PCElbowPlot(spleen)
碎石圖.jpeg

選擇了前10個PC成分

spleen <- FindClusters(spleen, reduction.type = "pca", dims.use = 1:10, resolution = 0.6, print.output = 0, save.SNN = TRUE)
PrintFindClustersParams(spleen)

Parameters used in latest FindClusters calculation run on: 2018-10-01 21:59:55
=============================================================================
Resolution: 0.6
-----------------------------------------------------------------------------
Modularity Function    Algorithm         n.start         n.iter
     1                   1                 100             10
-----------------------------------------------------------------------------
Reduction used          k.param          prune.SNN
     pca                 30                0.0667
-----------------------------------------------------------------------------
Dims used in calculation
=============================================================================
1 2 3 4 5 6 7 8 9 10

細胞聚類

spleen <- RunTSNE(spleen, dims.use = 1:10, do.fast= TRUE)
TSNEPlot(spleen)
TSNE.jpeg
> saveRDS(spleen, file = "/spleen_1.rds")

將R變量保存，利于后續(xù)的分析。

一些補充：
過濾低質量細胞：
在 scRNA-seq 分析中，有些細胞質量比較低，比如細胞處于凋亡狀態(tài),細胞中 RNA 發(fā)生降解等,這些細胞的存在會影響分析，因此我們第一步需要對細胞進行過濾。主要可分為三類:

①利用細胞檢測到的基因數(shù)或者是 reads 比對率來判斷技術噪音。
但不管是基因檢測數(shù)目還是比對率都跟實驗方法有很大相關性。 如果比對率太低,表明 RNA 可能發(fā)生了降解,或者文庫有污染或者細胞裂解不完全。

②如果實驗中加入了 spike-ins（本實驗沒有），可以通過計算比對到內源性 RNA 和外源性 RNA(spike-ins)的 reads 比例來過濾低質量細胞。
比值偏低表明細胞中的 RNA 數(shù)量較低，細胞可丟棄。但是也需要注意其實當細胞狀態(tài)不一樣，比如處于不同細胞周期時，細胞的 RNA 數(shù)量是具有很大差異的。不過我們依然認為在一大群細胞中，spike-ins比例特別高的細胞在很大概率上應該被排除在外。軟件 SinQC (Single-cell RNA-seq Quality Control)可以根據(jù)比對率和檢測到的基因數(shù)來過濾細胞。

③根據(jù)整體的基因表達譜來定義技術噪音。
比如對細胞進行聚類分析，PCA 分析等，將 outlier 細胞刪除，或者細胞表達中位值低于某一設定閾值時將該細胞過濾掉。當然這種方法也存在誤刪具有真正生物學差異的細胞,因此在刪除細胞時需要小心，可與上述另外兩種方法連用。

如果你的數(shù)據(jù)量過大，使用Seurat時內存不足，請看
實驗記錄11：海量scRNA-seq數(shù)據(jù)的質量控制、PCA、聚類

最后編輯于 ：
?著作權歸作者所有,轉載或內容合作請聯(lián)系作者
【社區(qū)內容提示】社區(qū)部分內容疑似由AI輔助生成，瀏覽時請結合常識與多方信息審慎甄別。
平臺聲明：文章內容（如有圖片或視頻亦包括在內）由作者上傳并發(fā)布，文章內容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務。

相關閱讀更多精彩內容

友情鏈接更多精彩內容