前面我們已經學習了單細胞轉錄組分析的：使用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之間的距離更明顯
從上面的圖可以看出不同樣本其實是有批次效應的，下一講我會介紹兩種去批次效應的方法。
因水平有限，有錯誤的地方，歡迎批評指正！