前幾天看到周運(yùn)來老師的一句話，放在這兒，與諸君共勉：

看到才能想到，想到才能做到，做到才能得到，得到才能失去，只有失去了才知道適不適合自己

好了，切入正題，今天想給大家分享的是入門級(jí)的經(jīng)典的Seurat分析流程，采用的數(shù)據(jù)為10X Genomics官方提供的外周血單個(gè)核細(xì)胞(peripheral blood mononuclear cell, PBMC)的測序數(shù)據(jù)。放一個(gè)下載鏈接：點(diǎn)擊即可下載，建議電腦操作!。這里我們下載的實(shí)際上就是可以直接利用Seurat進(jìn)行分析的數(shù)據(jù)。經(jīng)過層層解壓，你會(huì)得到這三個(gè)文件：

文件內(nèi)容.jpg

下面就正式開始利用Seurat包進(jìn)行分析。我個(gè)人將Seurat的上游分析分為幾個(gè)部分，分別是數(shù)據(jù)質(zhì)控、特征選擇、降維聚類、差異基因篩選與可視化。
Seurat包目前已經(jīng)更新到4.0的版本了，還好它已經(jīng)托管到了CRAN下面，可以直接用install.packages()進(jìn)行下載。

install.packages('Seurat')
library(Seurat)

1. 首先使用Read10X()讀入數(shù)據(jù)：

#保證那三個(gè)文件位于這個(gè)工作目錄下
pbmc <- Read10X(data.dir = 'C:/Users/DELL/Downloads/filtered_feature_bc_matrix')

當(dāng)然對(duì)于這個(gè)函數(shù)，你?Read10X一下，你可能會(huì)發(fā)現(xiàn)一些你以后可能會(huì)用到的參數(shù)，例如其中的gene.column和cell.column參數(shù)，對(duì)著你的features.tsv和barcodes.tsv文件，看看它是啥意思。

2. 然后就是創(chuàng)建SeuratObject對(duì)象：

pbmc <- CreateSeuratObject(counts = pbmc, project = "pbmc3k", min.cells = 3, min.features = 200)

實(shí)際上在這里是有一小步的質(zhì)控的，就在于min.cells和min.features兩個(gè)參數(shù)上，這兩個(gè)參數(shù)的意思是：過濾掉表達(dá)基因種類小于200個(gè)的細(xì)胞和被表達(dá)在不足3個(gè)細(xì)胞中的基因。

3. 對(duì)線粒體基因進(jìn)行統(tǒng)計(jì)

看你自己項(xiàng)目的實(shí)際需求，比如你做單細(xì)胞核測序，一般來說就應(yīng)該要做這一步，當(dāng)然你也可以用同樣的正則匹配去做核糖體基因等。

mt.genes <- grep(pattern = "^MT-", x = rownames(pbmc), value = T)

注意，我們做的人的數(shù)據(jù)，人的線粒體基因一般是以MT開頭的，我們在pbmc這個(gè)文件的行名中匹配這個(gè)模式，找打線粒體基因，value這個(gè)參數(shù)，設(shè)置為T就是返回線粒體基因的名稱，否則返回所在的行號(hào)。
對(duì)每個(gè)細(xì)胞計(jì)算其線粒體基因的含量，在這里我列出3中方法，當(dāng)然也只是在看教程而已：

pbmc[['percent.mt']] <- PercentageFeatureSet(object = pbmc, pattern = "^MT-")
pbmc <- PercentageFeatureSet(object = pbmc, pattern = "^MT-", col.name = "percent.mt")
pbmc <- PercentageFeatureSet(object = pbmc, features = mt.genes, col.name = "percent.mt")

你會(huì)發(fā)現(xiàn)，其實(shí)我們也可以自己給它一個(gè)基因的集合，讓它去進(jìn)行計(jì)算，特別是對(duì)于那些不太方便用正則表達(dá)式進(jìn)行展示的基因集合。

4. 數(shù)據(jù)質(zhì)控

這個(gè)部分需要根據(jù)你自己的數(shù)據(jù)選擇合適自己的參數(shù)，這里我和官方文檔設(shè)置的一樣。

pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nCount_RNA < 2500 & percent.mt < 5)

核心意思在于過濾后只剩下：表達(dá)基因超過200種，少于2500種和線粒體基因含量小于5%的細(xì)胞。

5. 降維聚類

pbmc <- NormalizeData(object = pbmc)
#=============findvariablefeatures==============#
pbmc <- FindVariableFeatures(object = pbmc, nfeatures = 2000)
#=================scale data=================#
pbmc <- ScaleData(object = pbmc, features = rownames(pbmc))
#=================PCA analysis===================#
pbmc <- RunPCA(object = pbmc, features = VariableFeatures(pbmc), verbose = F)

至于什么是NormalizeData和ScaleData，我在單細(xì)胞分析中的NormalizeData()與ScaleData()區(qū)別在哪兒？寫過一點(diǎn)自己的感悟。這里為什么ScaleData()是針對(duì)于全部的基因，官方也給出了解釋，確實(shí)這樣會(huì)在這一步消耗更多的時(shí)間，特別是在大數(shù)據(jù)集的時(shí)候，但是如果我們在這里不對(duì)高變基因以外的基因進(jìn)行Scale，就會(huì)在繪制heatmap時(shí)出現(xiàn)不必要的麻煩。
下一個(gè)重點(diǎn)在于我們選擇多少個(gè)主成分進(jìn)行下游的分析，Seurat給了兩種圖形顯示讓我們對(duì)我們的數(shù)據(jù)有直觀的感受。

pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
JackStrawPlot(pbmc, dims = 1:15)
ElbowPlot(pbmc)

進(jìn)行UMAP和TSNE聚類分析，兩者的差異在我的Literature Review(1): A short but comprehensive comparison about tSNE, UMAP and PCA這篇帖子里有過分享，當(dāng)然在這里我都做了，結(jié)果會(huì)儲(chǔ)存在SeuratObject中的不同位置，彼此并不會(huì)相互影響，這也是SeuratObject的一個(gè)突出優(yōu)勢。

#=============UMAP & TSNE analysis================#
pbmc <- FindNeighbors(object = pbmc, reduction = "pca", dims = 1:10)
pbmc <- FindClusters(object = pbmc, resolution = 0.5)
pbmc <- RunUMAP(object = pbmc, dims = 1:10)
pbmc <- RunTSNE(object = pbmc, dims = 1:10)

注意，這里FindClusters()函數(shù)中的resolution參數(shù)非常重要，這個(gè)值越大，聚出來的cluster也就越多。

6. 類型注釋

一般來說，單細(xì)胞的上游分析，最費(fèi)勁的地方就在于調(diào)整你的cluster數(shù)，聚出合適的類，進(jìn)行比較準(zhǔn)確的注釋。在這里我就直接套用官方的聚類注釋結(jié)果，直接進(jìn)行注釋了。

cluster_cell_type <- c("Naive CD4 T", "Memory CD4 T", "CD14+ Mono", "B", "CD8 T", "FCGR3A+ Mono", 
                     "NK", "DC", "Platelet")
names(cluster_cell_type) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, cluster_cell_type)
pbmc[['cell_type']] <- pbmc@active.ident

注意，pbmc[[]]這種形式是直接在meta.data中添加元素（列），是在細(xì)胞的水平上，這在我們前面的統(tǒng)計(jì)線粒體基因處也用到了。當(dāng)然有很多種不同的將注釋結(jié)果映射到每個(gè)細(xì)胞上的方法，后面有時(shí)間再一起總結(jié)、學(xué)習(xí)！

#TSNE
DimPlot(object = pbmc, 
        reduction = 'tsne',
        cols = c('#FF0033','#CC3399','#993366','#009933', '#0066CC','#003399','#9966CC','#99CC33','#003399'),
        label = T,
        label.box = T,
        pt.size = 1)
#UMAP
DimPlot(object = pbmc, 
        reduction = 'umap',
        cols = c('#FF0033','#CC3399','#993366','#009933', '#0066CC','#003399','#9966CC','#99CC33','#003399'),
        label = T,
        label.box = T,
        pt.size = 1)

tSNE.png

UMAP.png

看來，兩種不同的聚類方式繪制出的圖還是很不一樣的哈！最后，我想說，Seurat最好的教程是它的官網(wǎng)Seurat，沒有誰比創(chuàng)造它的人更清楚它是怎么工作的了，我們不過在重復(fù)它自己給出的教程罷了。

今天也是摸魚的一天！