單細胞分析實錄(5): Seurat標準流程

TOP生物資訊發表於2021-01-05

前面我們已經學習了單細胞轉錄組分析的：使用Cell Ranger得到表達矩陣和doublet檢測，今天我們開始Seurat標準流程的學習。這一部分的內容，網上有很多帖子，基本上都是把Seurat官網PBMC的例子重複一遍，這回我換一個資料集，細胞型別更多，同時也會加入一些實際分析中很有用的技巧。

1. 匯入資料，建立Seurat物件

library(Seurat)
library(tidyverse)
testdf=read.table("test_20210105.txt",header = T,row.names = 1)
test.seu=CreateSeuratObject(counts = testdf)

看一下長什麼樣子

> test.seu
An object of class Seurat 
33538 features across 6746 samples within 1 assay 
Active assay: RNA (33538 features)
#1 assay表示有一套資料，假如用Seurat裡面的函式去過批次效應，
#這裡會有2個assay，另外一個是去批次（整合）之後的；
#cell hashing的原始表達矩陣（行為gene+tag）生成的Seurat物件也會有2個assay，
#另外一個assay為HTO的矩陣；
#其他情況以此類推

測試資料有33538個基因，6746個細胞。除此之外，還要關注一下另外兩層資訊：test.seu@meta.data這個資料框用來儲存後設資料，每一個細胞都有多個屬性；test.seu[["RNA"]]@counts這個稀疏矩陣用來儲存原始UMI表達矩陣。

> head(test.seu@meta.data)
                   orig.ident nCount_RNA nFeature_RNA
A_AAACCCAAGGGTCACA          A       3714         1151
A_AAACCCAAGTATAACG          A       1855          816
A_AAACCCAGTCTCTCAC          A       1530          823
A_AAACCCAGTGAGTCAG          A      11145         1087
A_AAACCCAGTGGCACTC          A       2289          834
A_AAACGAAAGCCAGAGT          A       3714          990
#我這裡CB的前面人為加上了樣本來源，用下劃線連線，orig.ident是自動識別得到的

> test.seu[["RNA"]]@counts[1:4,1:4]
4 x 4 sparse Matrix of class "dgCMatrix"
            A_AAACCCAAGGGTCACA A_AAACCCAAGTATAACG A_AAACCCAGTCTCTCAC A_AAACCCAGTGAGTCAG
MIR1302-2HG                  .                  .                  .                  .
FAM138A                      .                  .                  .                  .
OR4F5                        .                  .                  .                  .
AL627309.1                   .                  .                  .                  .

2. 簡單過濾

接下來，我們根據每個細胞內部線粒體基因表達佔比、檢測到的基因數、檢測的UMI總數這三個方面來對細胞進行簡單的過濾。
先計算細胞內線粒體基因表達佔比，類似的核糖體基因(大多為RP開頭)也能這樣計算，還要注意不要將線粒體基因的MT-寫成了MT，不然就把別的基因也算進去了：

test.seu[["percent.mt"]] <- PercentageFeatureSet(test.seu, pattern = "^MT-")  #正規表示式，表示以MT-開頭；test.seu[["percent.mt"]]這種寫法會在meta.data矩陣加上一列

這裡我已經根據預先設定好的閾值過濾了，程式碼如下

test.seu <- subset(test.seu, subset = nCount_RNA > 1000 & 
                     nFeature_RNA < 5000 & 
                     percent.mt < 30 & 
                     nFeature_RNA > 600)

過濾之後的數值分佈如下，用到VlnPlot()函式，該繪圖函式裡面的feature引數可以是meta.data矩陣的某一列，也可以是某一個基因，很多文章都用這種圖展示marker gene

VlnPlot(test.seu,features = c("nCount_RNA", "nFeature_RNA", "percent.mt"),pt.size = 0)

nFeature_RNA/nCount_RNA不能太小（空液滴），不能太大（doublet、測序技術限制）， 而且閾值設定要綜合多個樣本來看，像下面這樣

一般在CD45陰性的細胞中percent.mt的閾值大一些，50%也看過幾次了

3. 標準化，消除文庫大小的影響

如何標準化：LogNormalize: Feature counts for each cell are divided by the total counts for that cell and multiplied by the scale.factor. This is then natural-log transformed using log1p.（先相除，再求對數）

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

標準化之後的矩陣儲存在test.seu[["RNA"]]@data

> test.seu[["RNA"]]@data[1:4,1:4]
4 x 4 sparse Matrix of class "dgCMatrix"
            A_AAACCCAAGGGTCACA A_AAACCCAAGTATAACG A_AAACCCAGTCTCTCAC A_AAACCCAGTGAGTCAG
MIR1302-2HG                  .                  .                  .                  .
FAM138A                      .                  .                  .                  .
OR4F5                        .                  .                  .                  .
AL627309.1                   .                  .                  .                  .

4. 找Variable基因

因為單細胞表達矩陣很稀疏（很多0），選high variable基因的目的可以找到包含資訊最多的基因（很多基因的表達差不多都是0），同時極大提升軟體執行速度

test.seu <- FindVariableFeatures(test.seu, selection.method = "vst", nfeatures = 2000)

這些基因儲存在VariableFeatures(test.seu)，有時候可能需要人為指定high variable基因，可以這樣：

VariableFeatures(test.seu)="specific genes"

5. 歸一化表達矩陣

（基於前面得到的data矩陣）
這一步之後，所有基因的表達值的分佈就差不多了，不然表達值不在一個數量級，對後續降維聚類影響挺大。新的矩陣儲存在test.seu[["RNA"]]@scale.data裡面。

test.seu <- ScaleData(test.seu, features = rownames(test.seu))

預設只對上一步選出來的基因scale，這裡調整為所有基因，是為了方便以後畫熱圖（畫熱圖一般會用scale之後的z-score）

6. 降維聚類

（基於前面得到的high variable基因的scale矩陣）

test.seu <- RunPCA(test.seu, npcs = 50, verbose = FALSE)
test.seu <- FindNeighbors(test.seu, dims = 1:30)
test.seu <- FindClusters(test.seu, resolution = 0.5)

test.seu <- RunUMAP(test.seu, dims = 1:30)
test.seu <- RunTSNE(test.seu, dims = 1:30)

Run開頭的函式降維，Find開頭的函式聚類，一般就這幾步，相對固定。PCA將原來2000維的資料降到50維，dims參數列示使用多少個主成分（一般20左右就可以了，多幾個少幾個對結果影響不大），resolution參數列達聚類的解析度，這個值大於0，一般都是在0-1範圍裡面調整，越大得到的cluster越多，這個值可以反覆調整，並不會改變降維的結果（也就是tsne、umap圖的二維座標）。

這一步之後的資料是這樣的

> test.seu
An object of class Seurat 
33538 features across 6746 samples within 1 assay 
Active assay: RNA (33538 features)
 3 dimensional reductions calculated: pca, umap, tsne
 # 幾種降維方式都會呈現出來

聚類之後test.seu@meta.data多了兩列，RNA_snn_res.0.5記錄了你用的解析度，最終的聚類結果儲存在seurat_clusters中

> head(test.seu@meta.data)
                  orig.ident nCount_RNA nFeature_RNA percent.mt RNA_snn_res.0.5
A_AAACCCAAGGGTCACA          A       3714         1151   9.585353               8
A_AAACCCAAGTATAACG          A       1855          816  12.776280               0
A_AAACCCAGTCTCTCAC          A       1530          823  14.248366              12
A_AAACCCAGTGAGTCAG          A      11145         1087   2.853297               4
A_AAACCCAGTGGCACTC          A       2289          834  15.640017               3
A_AAACGAAAGCCAGAGT          A       3714          990   5.654281               0
                  seurat_clusters
A_AAACCCAAGGGTCACA               8
A_AAACCCAAGTATAACG               0
A_AAACCCAGTCTCTCAC              12
A_AAACCCAGTGAGTCAG               4
A_AAACCCAGTGGCACTC               3
A_AAACGAAAGCCAGAGT               0

7. tsne/umap展示結果

library(cowplot)
test.seu$patient=str_replace(test.seu$orig.ident,"_.*$","")
p1 <- DimPlot(test.seu, reduction = "tsne", group.by = "patient", pt.size=0.5)
p2 <- DimPlot(test.seu, reduction = "tsne", group.by = "ident",   pt.size=0.5, label = TRUE,repel = TRUE) #後面兩個引數用來新增文字標籤
p3 <- DimPlot(test.seu, reduction = "umap", group.by = "patient", pt.size=0.5)
p4 <- DimPlot(test.seu, reduction = "umap", group.by = "ident",   pt.size=0.5, label = TRUE,repel = TRUE)

fig_tsne <- plot_grid(p1, p2, labels = c('patient','ident'),align = "v",ncol = 2)
ggsave(filename = "tsne.pdf", plot = fig_tsne, device = 'pdf', width = 30, height = 15, units = 'cm')
fig_umap <- plot_grid(p3, p4, labels = c('patient','ident'),align = "v",ncol = 2)
ggsave(filename = "umap.pdf", plot = fig_umap, device = 'pdf', width = 30, height = 15, units = 'cm')

ident表示每個細胞的標籤，聚類之後就是聚類的結果，在一些特定場景可以更換。

在umap圖中，cluster之間的距離更明顯

從上面的圖可以看出不同樣本其實是有批次效應的，下一講我會介紹兩種去批次效應的方法。

因水平有限，有錯誤的地方，歡迎批評指正！

相關文章