前幾天看到周運(yùn)來(lái)老師的一句話搏明，放在這兒，與諸君共勉：
看到才能想到闪檬，想到才能做到星著，做到才能得到，得到才能失去粗悯，只有失去了才知道適不適合自己
好了虚循，切入正題，今天想給大家分享的是入門(mén)級(jí)的經(jīng)典的Seurat分析流程样傍，采用的數(shù)據(jù)為10X Genomics官方提供的外周血單個(gè)核細(xì)胞(peripheral blood mononuclear cell, PBMC)的測(cè)序數(shù)據(jù)横缔。放一個(gè)下載鏈接：點(diǎn)擊即可下載，建議電腦操作!铭乾。這里我們下載的實(shí)際上就是可以直接利用Seurat進(jìn)行分析的數(shù)據(jù)剪廉。經(jīng)過(guò)層層解壓，你會(huì)得到這三個(gè)文件：
當(dāng)然炕檩，這三個(gè)文件是什么斗蒋，可以看看我前面的文章我眼中的barcodes.tsv.gz/features.tsv.gz/matrix.mtx.gz捌斧。
下面就正式開(kāi)始利用Seurat包進(jìn)行分析。我個(gè)人將Seurat的上游分析分為幾個(gè)部分泉沾，分別是數(shù)據(jù)質(zhì)控捞蚂、特征選擇、降維聚類(lèi)跷究、差異基因篩選與可視化姓迅。
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ù)的意思是：過(guò)濾掉表達(dá)基因種類(lèi)小于200個(gè)的細(xì)胞和被表達(dá)在不足3個(gè)細(xì)胞中的基因觉增。
3. 對(duì)線粒體基因進(jìn)行統(tǒng)計(jì)
看你自己項(xiàng)目的實(shí)際需求，比如你做單細(xì)胞核測(cè)序翻斟，一般來(lái)說(shuō)就應(yīng)該要做這一步逾礁，當(dāng)然你也可以用同樣的正則匹配去做核糖體基因等。
mt.genes <- grep(pattern = "^MT-", x = rownames(pbmc), value = T)
注意访惜，我們做的人的數(shù)據(jù)敞斋，人的線粒體基因一般是以MT開(kāi)頭的，我們?cè)趐bmc這個(gè)文件的行名中匹配這個(gè)模式疾牲，找打線粒體基因植捎，value這個(gè)參數(shù)，設(shè)置為T(mén)就是返回線粒體基因的名稱阳柔，否則返回所在的行號(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)
核心意思在于過(guò)濾后只剩下：表達(dá)基因超過(guò)200種低滩，少于2500種和線粒體基因含量小于5%的細(xì)胞召夹。
5. 降維聚類(lèi)
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ū)別在哪兒恕沫？寫(xiě)過(guò)一點(diǎn)自己的感悟监憎。這里為什么ScaleData()是針對(duì)于全部的基因，官方也給出了解釋?zhuān)_實(shí)這樣會(huì)在這一步消耗更多的時(shí)間婶溯，特別是在大數(shù)據(jù)集的時(shí)候鲸阔，但是如果我們?cè)谶@里不對(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聚類(lèi)分析叙身，兩者的差異在我的Literature Review(1): A short but comprehensive comparison about tSNE, UMAP and PCA這篇帖子里有過(guò)分享死讹，當(dāng)然在這里我都做了，結(jié)果會(huì)儲(chǔ)存在SeuratObject中的不同位置曲梗，彼此并不會(huì)相互影響，這也是SeuratObject的一個(gè)突出優(yōu)勢(shì)妓忍。
#=============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è)值越大世剖，聚出來(lái)的cluster也就越多定罢。
6. 類(lèi)型注釋
一般來(lái)說(shuō)，單細(xì)胞的上游分析旁瘫，最費(fèi)勁的地方就在于調(diào)整你的cluster數(shù)祖凫，聚出合適的類(lèi)，進(jìn)行比較準(zhǔn)確的注釋酬凳。在這里我就直接套用官方的聚類(lèi)注釋結(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)
看來(lái)攘蔽，兩種不同的聚類(lèi)方式繪制出的圖還是很不一樣的哈！最后呐粘，我想說(shuō)满俗，Seurat最好的教程是它的官網(wǎng)Seurat转捕，沒(méi)有誰(shuí)比創(chuàng)造它的人更清楚它是怎么工作的了谋国，我們不過(guò)在重復(fù)它自己給出的教程罷了顿苇。
今天也是摸魚(yú)的一天！
