標(biāo)題：Seurat - Guided Clustering Tutorial
網(wǎng)址：https://satijalab.org/seurat/v3.2/pbmc3k_tutorial.html

部分內(nèi)容轉(zhuǎn)載來自：caokai001 鏈接：http://www.reibang.com/p/7e4c205a576e
部分內(nèi)容轉(zhuǎn)載來自： 劉小澤 鏈接：http://www.reibang.com/p/b46b6b6d344f

一、設(shè)置Seurat對象

對于本教程，我們將分析可從10X Genomics免費(fèi)獲得的外周血單個核細(xì)胞（PBMC）數(shù)據(jù)集究西。Illumina NextSeq 500上已對2700個單細(xì)胞進(jìn)行了測序绵脯。原始數(shù)據(jù)可在此處找到窒百。

我們從讀取數(shù)據(jù)開始仿滔。該Read10X函數(shù)從10X讀取cellranger管道的輸出月趟，返回唯一的分子識別（UMI）計數(shù)矩陣。該矩陣中的值表示在每個單元格（列）中檢測到的每個特征（即基因嫌拣；行）的分子數(shù)。

接下來呆躲，我們使用計數(shù)矩陣創(chuàng)建一個Seurat對象异逐。該對象用作包含單個單元格數(shù)據(jù)集的數(shù)據(jù)（如計數(shù)矩陣）和分析（如PCA或聚類結(jié)果）的容器。有關(guān)Seurat對象結(jié)構(gòu)的技術(shù)討論插掂，請查看我們的GitHub Wiki灰瞻。例如，計數(shù)矩陣存儲在中pbmc[["RNA"]]@counts辅甥。

library(dplyr)
library(Seurat)
library(patchwork)

# Load the PBMC dataset
pbmc.data <- Read10X(data.dir = "./10x_Rawdata/pbmc3k_filtered_gene_bc_matrices/filtered_gene_bc_matrices/hg19/")
# Initialize the Seurat object with the raw (non-normalized data).
#每個基因至少在3個細(xì)胞中表達(dá)酝润，每個細(xì)胞中至少有200個基因表達(dá)。
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
> pbmc
An object of class Seurat 
13714 features across 2700 samples within 1 assay 
Active assay: RNA (13714 features, 0 variable features)

# Lets examine a few genes in the first thirty cells
>pbmc.data[c("CD3D", "TCL1A", "MS4A1"), 1:30]
## 3 x 30 sparse Matrix of class "dgCMatrix"
##                                                                    
## CD3D  4 . 10 . . 1 2 3 1 . . 2 7 1 . . 1 3 . 2  3 . . . . . 3 4 1 5
## TCL1A . .  . . . . . . 1 . . . . . . . . . . .  . 1 . . . . . . . .
## MS4A1 . 6  . . . . . . 1 1 1 . . . . . . . . . 36 1 2 . . 2 . . . .

.矩陣中的值表示0（沒有檢測到分子）璃弄。由于scRNA-seq矩陣中的大多數(shù)值為0要销，因此Seurat盡可能使用稀疏矩陣表示。這樣可以為Drop-seq / inDrop / 10x數(shù)據(jù)節(jié)省大量內(nèi)存并節(jié)省速度夏块。

>dense.size <- object.size(as.matrix(pbmc.data))
>dense.size
## 709591472 bytes
>sparse.size <- object.size(pbmc.data)
>sparse.size
## 29905192 bytes
>dense.size/sparse.size
## 23.7 bytes

二蕉陋、標(biāo)準(zhǔn)的預(yù)處理流程

以下步驟包含Seurat中scRNA-seq數(shù)據(jù)的標(biāo)準(zhǔn)預(yù)處理工作流程。這些代表基于質(zhì)量控制指標(biāo)拨扶，數(shù)據(jù)歸一化和縮放以及高度可變特征的檢測來選擇和過濾細(xì)胞。

2.1 質(zhì)控和選擇細(xì)胞進(jìn)行進(jìn)一步分析

通過Seurat茁肠，您可以輕松地探索質(zhì)量控制指標(biāo)并根據(jù)任何用戶定義的條件過濾單元患民。社區(qū)常用的一些質(zhì)量控制指標(biāo)包括

• 每個細(xì)胞中檢測到的獨(dú)特基因的數(shù)量。
  • 低質(zhì)量的細(xì)胞或空液滴通常很少有基因
  • 細(xì)胞雙峰或多重峰可能顯示異常高的基因計數(shù)
• 同樣垦梆，細(xì)胞內(nèi)檢測到的分子總數(shù)（與獨(dú)特基因高度相關(guān)）
• 映射到線粒體基因組的讀段的百分比
  • 低質(zhì)量/瀕死細(xì)胞通常表現(xiàn)出廣泛的線粒體污染
  • 我們使用PercentageFeatureSet函數(shù)計算線粒體QC指標(biāo)匹颤，該函數(shù)計算源自一組功能的計數(shù)的百分比
  • 我們使用所有MT-以線粒體基因開頭的基因

> pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
> head(pbmc@meta.data,5)
                 orig.ident nCount_RNA nFeature_RNA percent.mt
AAACATACAACCAC-1     pbmc3k       2419          779  3.0177759
AAACATTGAGCTAC-1     pbmc3k       4903         1352  3.7935958
AAACATTGATCAGC-1     pbmc3k       3147         1129  0.8897363
AAACCGTGCTTCCG-1     pbmc3k       2639          960  1.7430845
AAACCGTGTATGCG-1     pbmc3k        980          521  1.2244898

在下面的示例中，我們可視化QC指標(biāo)托猩，并使用它們來過濾單元格印蓖。

• 過濾 unique feature counts 超過2500或少于200的細(xì)胞
• 過濾線粒體 > 5％的細(xì)胞

> VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

image.png

plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
#CombinePlots(plots = list(plot1,plot2))

image.png

過濾

pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
>pbmc
An object of class Seurat 
13714 features across 2638 samples within 1 assay 
Active assay: RNA (13714 features, 0 variable features)

之前是13714 features across 2700 samples ，過濾掉了62個細(xì)胞

三京腥、 預(yù)處理之歸一化

關(guān)于Seurat歸一化原理赦肃，可以看這一篇：https://www.biorxiv.org/content/biorxiv/early/2019/03/18/576827.full.pdf

標(biāo)準(zhǔn)化是一個比較簡單的過程，使用的是"logNormalize", 就是將每個基因的表達(dá)量對該細(xì)胞總表達(dá)量進(jìn)行平衡公浪，然后乘以一個因子(scale factor, 默認(rèn)值為10,000), 然后中進(jìn)行對數(shù)轉(zhuǎn)換他宛。

默認(rèn)是進(jìn)行一個全局的LogNormalize操作：log1p(value/colSums[cell-idx] *scale_factor) ，其中log1p指的是log(x + 1) 欠气，得到的結(jié)果存儲在：pbmc[["RNA"]]@data

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

四厅各、鑒定差異基因HVGs（highly variable features）

我們在Seurat3的流程詳細(xì)描述在這里，并通過直接建模在單細(xì)胞數(shù)據(jù)固有的均方差關(guān)系改進(jìn)了以前的版本预柒，并在實(shí)現(xiàn)FindVariableFeatures功能队塘。V3默認(rèn)選擇2000個差異基因袁梗。這些將用于下游分析，例如PCA憔古。

如果手動計算遮怜，相當(dāng)于計算變異系數(shù)，而不是標(biāo)準(zhǔn)差投放。因?yàn)闃?biāo)準(zhǔn)差對于數(shù)據(jù)呈現(xiàn)對稱分布奈泪，及單峰情況體現(xiàn)差異比較好，但是對于雙峰情況灸芳，變異系數(shù)cv 可能更好涝桅。

pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
SLM <<- FindVariableFeatures(SLM, selection.method = "mvp", mean.cutoff = c(0.1,100),dispersion.cutoff=c(-1,Inf))
##"mvp"這種寫法，最后鑒定出來了4000多個差異基因烙样，不受默認(rèn)值影響冯遂。
# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(pbmc), 10)

• V3橫坐標(biāo)范圍設(shè)定參數(shù)改成：mean.cutoff；縱坐標(biāo)改成：dispersion.cutoff
• 鑒定差異基因的算法包含三種：vst(默認(rèn))谒获、mean.var.plot蛤肌、dispersion
  • vst：首先利用loess對 log(variance) 和log(mean) 擬合一條直線，然后利用觀測均值和期望方差對基因表達(dá)量進(jìn)行標(biāo)準(zhǔn)化批狱。最后根據(jù)保留最大的標(biāo)準(zhǔn)化的表達(dá)量計算方差
  • mean.var.plot: 首先利用mean.function和 dispersion.function分別計算每個基因的平均表達(dá)量和離散情況裸准；然后根據(jù)平均表達(dá)量將基因們分散到一定數(shù)量（默認(rèn)是20個）的小區(qū)間（bin）中，并且計算每個bin中z-score
  • dispersion：挑選最高離散值的基因
• V3默認(rèn)選擇2000個差異基因赔硫，檢查方法也不同（V3用VariableFeatures(sce)檢查炒俱，V2用sce@var.genes檢查）

# plot variable features with and without labels
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
> plot1 + plot2 
Error in grid.Call(C_convert, x, as.integer(whatfrom), as.integer(whatto),  : 
  Viewport has zero dimension(s)
此外: Warning messages:
1: Transformation introduced infinite values in continuous x-axis 
2: Transformation introduced infinite values in continuous x-axis 

有問題，未實(shí)現(xiàn)圖片：

紅色的HVGs 有2000個爪膊，從上面 FindVariableFeatures 权悟，命令參數(shù)就可以知道.

image.png

五、標(biāo)準(zhǔn)化

從上圖中推盛，我們可以看到很多基因都表現(xiàn)出很大的差異性峦阁，但是這些差異性并不一定都是事實(shí)，有些可能是因?yàn)閷?shí)驗(yàn)耘成，或者數(shù)據(jù)采集的方式導(dǎo)入的榔昔。因此人們想盡量的消除這些因素帶來的差異。這些差異有可能是因?yàn)閎atch-effect瘪菌， 或者cell alignment rate件豌，或者其它。這些因素都可以寫入metadata中控嗜，然后通過ScaleData函數(shù)來消除茧彤。

ScaleData 可以通過metadata信息，消除很多其他因素影響疆栏。

5.1 縮放的標(biāo)準(zhǔn)

? 應(yīng)用線性變換縮放數(shù)據(jù)曾掂，這是一個標(biāo)準(zhǔn)的預(yù)處理步驟惫谤，比PCA等降維技術(shù)更重要。使用ScaleData函數(shù):
（1）使每個基因在所有細(xì)胞間的表達(dá)量均值為0
（2）使每個基因在所有細(xì)胞間的表達(dá)量方差為1

縮放操作給予下游分析同等的權(quán)重珠洗，這樣高表達(dá)基因就不會占主導(dǎo)地位溜歪，其結(jié)果存儲在pbmc[["RNA"]]@scale.data中，這里默認(rèn)ScaleData 進(jìn)行了中心化许蓖，標(biāo)準(zhǔn)化蝴猪。

## 全局縮放
all.genes <- rownames(pbmc)
pbmc <- ScaleData(object = pbmc, features = all.genes)
> head(pbmc[["RNA"]]@scale.data[, 1:3])

##               AAACATACAACCAC AAACATTGAGCTAC AAACCGTGCTTCCG
## AL627309.1       -0.06452497    -0.06452497    -0.06452497
## AP006222.2       -0.03726409    -0.03726409    -0.03726409
## RP11-206L10.2    -0.04153370    -0.04153370    -0.04153370
## RP11-206L10.9    -0.03734170    -0.03734170    -0.03734170
## LINC00115        -0.08343360    -0.08343360    -0.08343360
## NOC2L            -0.32126699    -0.32126699    -0.32126699

縮放是Seurat工作流程中必不可少的步驟，但僅限于將用作PCA輸入的基因膊爪。因此自阱，默認(rèn)設(shè)置ScaleData是僅對先前確定的變量特征執(zhí)行縮放（默認(rèn)為2,000）。為此米酬，請忽略features上一個函數(shù)調(diào)用中的參數(shù)沛豌，即

>pbmc <- ScaleData(pbmc)   #默認(rèn)為2,000個高邊基因

這樣操作的結(jié)果中降維和聚類不會被影響，但是只對HVGs基因進(jìn)行標(biāo)準(zhǔn)化赃额，下面畫的熱圖依然會受到全部基因中某些極值的影響加派。所以為了熱圖的結(jié)果，還是對所有基因進(jìn)行標(biāo)準(zhǔn)化比較好跳芳。

ScaleData的操作過程是怎樣的?

看一下幫助文檔就能了解：
ScaleData這個函數(shù)有兩個默認(rèn)參數(shù)：do.scale = TRUE和do.center = TRUE芍锦，然后需要輸入進(jìn)行scale/center的基因名（默認(rèn)是取前面FindVariableFeatures的結(jié)果）。
scale和center這兩個默認(rèn)參數(shù)應(yīng)該不陌生了飞盆，center就是對每個基因在不同細(xì)胞的表達(dá)量都減去均值娄琉；
scale就是對每個進(jìn)行center后的值再除以標(biāo)準(zhǔn)差（就是進(jìn)行了一個z-score的操作）

再看一個來自BioStar的討論：https://www.biostars.org/p/349824/

問題1：為什么在NormalizeData()之后，還需要進(jìn)行ScaleData操作桨啃？

前面NormalizeData進(jìn)行的歸一化主要考慮了測序深度/文庫大小的影響（reads *10000 divide by total reads, and then log），但實(shí)際上歸一化的結(jié)果還是存在基因表達(dá)量分布區(qū)間不同的問題檬输；而ScaleData做的就是在歸一化的基礎(chǔ)上再添zero-centres （mean/sd =》 z-score）照瘾，將各個表達(dá)量放在了同一個范圍中，真正實(shí)現(xiàn)了”表達(dá)量的同一起跑線“丧慈。

另外析命，新版的ScaleData函數(shù)還支持設(shè)定回歸參數(shù)，如（nUMI逃默、percent.mito）鹃愤，來校正一些技術(shù)噪音

問題2： FindVariableGenes() or RunPCA() or FindCluster()這些參數(shù)是基于歸一化Normalized_Data還是標(biāo)準(zhǔn)化Scaled_Data ？

所有操作都是基于標(biāo)準(zhǔn)化的數(shù)據(jù)Scaled_Data 完域，因?yàn)檫@個數(shù)據(jù)是針對基因間比較的软吐。例如有兩組基因表達(dá)量如下:

g1 10 20 30 40 50
g2 20 40 60 80 100

雖然它們看起來是g2的表達(dá)量是g1的兩倍，但真要降維聚類時吟税，就要看scale的結(jié)果：

> scale(c(10,20,30,40,50))
           [,1]
[1,] -1.2649111
[2,] -0.6324555
[3,]  0.0000000
[4,]  0.6324555
[5,]  1.2649111
attr(,"scaled:center")
[1] 30
attr(,"scaled:scale")
[1] 15.81139
> scale(c(20,40,60,80,100))
           [,1]
[1,] -1.2649111
[2,] -0.6324555
[3,]  0.0000000
[4,]  0.6324555
[5,]  1.2649111
attr(,"scaled:center")
[1] 60
attr(,"scaled:scale")
[1] 31.62278

二者標(biāo)準(zhǔn)化后的分布形態(tài)是一樣的（只是中位數(shù)不同）凹耙，而我們后來聚類也是看數(shù)據(jù)的分布姿现，分布相似聚在一起。所以歸一化差兩倍的兩個基因肖抱，根據(jù)標(biāo)準(zhǔn)化的結(jié)果依然聚在了一起

問題3： ScaleData()在scRNA分析中一定要用嗎备典？

實(shí)際上標(biāo)準(zhǔn)化這個過程不是單細(xì)胞分析固有的，在很多機(jī)器學(xué)習(xí)意述、降維算法中比較不同feature的距離時經(jīng)常使用提佣。當(dāng)不進(jìn)行標(biāo)準(zhǔn)化時，有時feature之間巨大的差異（就像??問題2中差兩倍的基因表達(dá)量荤崇，實(shí)際上可能差的倍數(shù)會更多）會影響分析結(jié)果拌屏，讓我們誤以為它們的數(shù)據(jù)分布相差很遠(yuǎn)。

很多時候天试，我們會混淆normalization和scale：那么看一句英文解釋：

Normalization "normalizes" within the cell for the difference in sequenicng depth / mRNA throughput. 主要著眼于樣本的文庫大小差異
Scaling "normalizes" across the sample for differences in range of variation of expression of genes . 主要著眼于基因的表達(dá)分布差異

另外槐壳，看scale()函數(shù)的幫助文檔就會發(fā)現(xiàn)，它實(shí)際上是對列進(jìn)行操作：

image

那么具體操作中是要先對表達(dá)矩陣轉(zhuǎn)置喜每，然后scale务唐，最后轉(zhuǎn)回來

t(scale(t(dat[genes,])))

問題4：RunCCA()函數(shù)需要?dú)w一化數(shù)據(jù)還是標(biāo)準(zhǔn)化數(shù)據(jù)？

通過看這個函數(shù)的幫助文檔：

RunCCA(object, object2, group1, group2, group.by, num.cc = 20, genes.use,
  scale.data = TRUE, rescale.groups = FALSE, ...)

顯示scale.data = TRUE带兜，那么就需要標(biāo)準(zhǔn)化后的數(shù)據(jù)

刪除不需要的數(shù)據(jù)

使用ScaleData函數(shù)從單細(xì)胞數(shù)據(jù)集中刪除不需要的變量枫笛。例如，我們可以“regress out”與細(xì)胞周期階段或線粒體污染相關(guān)的異質(zhì)性（去除不需要的變量刚照，即校正協(xié)變量刑巧，校正線粒體基因比例的影響）。如下：

# 刪除不需要的數(shù)據(jù)无畔，校正線粒體基因引起的影響
# 這個步驟非常耗時
pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")

我們對矩陣進(jìn)行了細(xì)胞水平質(zhì)控啊楚，但是如何對基因進(jìn)行質(zhì)控，當(dāng)然我們創(chuàng)建seurat 對象時候浑彰，已經(jīng)設(shè)定一個cutoff. 基因必須早多少個細(xì)胞中表達(dá). 但是如果直接用這些基因進(jìn)行聚類和降維（T-sne/UMAP 可視化 可能引起偏差恭理，所以需要先進(jìn)行一下PCA 初步選取主要的主成分，對應(yīng)著權(quán)重不同的基因郭变，再用這些基因進(jìn)行聚類颜价，降維可視化效果更好）

六、PCA

接下來诉濒，我們對縮放后的數(shù)據(jù)執(zhí)行PCA周伦。默認(rèn)情況下，僅將先前確定的變量特征用作輸入未荒，但是features如果您希望選擇其他子集专挪，則可以使用arguments進(jìn)行定義。默認(rèn)使用標(biāo)準(zhǔn)化的數(shù)據(jù)

pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))

# Examine and visualize PCA results a few different ways
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
## PC_ 1 
## Positive:  CST3, TYROBP, LST1, AIF1, FTL 
## Negative:  MALAT1, LTB, IL32, IL7R, CD2 
## PC_ 2 
## Positive:  CD79A, MS4A1, TCL1A, HLA-DQA1, HLA-DQB1 
## Negative:  NKG7, PRF1, CST7, GZMB, GZMA 
## PC_ 3 
## Positive:  HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPB1 
## Negative:  PPBP, PF4, SDPR, SPARC, GNG11 
## PC_ 4 
## Positive:  HLA-DQA1, CD79B, CD79A, MS4A1, HLA-DQB1 
## Negative:  VIM, IL7R, S100A6, IL32, S100A8 
## PC_ 5 
## Positive:  GZMB, NKG7, S100A8, FGFBP2, GNLY 
## Negative:  LTB, IL7R, CKB, VIM, MS4A7

Seurat提供了多種手段來查看PCA分析的結(jié)果，其中包括PrintPCA, VizPCA, PCAPlot, 以及PCHeatmap

【一維可視化】 對每一個維度狈蚤，權(quán)重較大的基因進(jìn)行可視化

VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")

image.png

【二維可視化】 將細(xì)胞降維到二維空間

DimPlot(object = pbmc, reduction = "pca")

image

這個圖中每個點(diǎn)就是一個樣本困肩，它根據(jù)PCA的結(jié)果坐標(biāo)進(jìn)行畫圖，這個坐標(biāo)保存在了cell.embeddings中：

> head(pbmc[['pca']]@cell.embeddings)[1:2,1:2]
                     PC_1       PC_2
AAACATACAACCAC -4.7296855 -0.5184265
AAACATTGAGCTAC -0.5174029  4.5918957

其實(shí)還支持定義分組：根據(jù)pbmc@meta.data中的ident分組：

> table(pbmc@meta.data$orig.ident)

pbmc3k 
  2638 
# 當(dāng)然這里只有一個分組

Tips: 下圖以PC_1 來看脆侮，是解釋細(xì)胞分布最大的基因組合方向锌畸，所有的基因是維度，有些基因?qū)@個方向共享很大靖避，也就是說更容易分開細(xì)胞潭枣，這些基因可能更加重要。下面是熱圖1展示這些基因?qū)τ趨^(qū)分細(xì)胞的效果幻捏，可以看到將很多細(xì)胞區(qū)分成紫色和黃色.效果不錯.

可以看到5,6維度貌似區(qū)分不開了盆犁。

探索異質(zhì)性來源—DimHeatmap

每個細(xì)胞和基因都根據(jù)PCA結(jié)果得分進(jìn)行了排序，默認(rèn)畫前30個基因（nfeatures設(shè)置）篡九，1個主成分(dims設(shè)置)谐岁，細(xì)胞數(shù)沒有默認(rèn)值（cells設(shè)置）

> DimHeatmap(pbmc, dims = 1:12, cells = 500, balanced = TRUE)
# 其中balanced表示正負(fù)得分的基因各占一半

image

降維的真實(shí)目的是盡可能減少具有相關(guān)性的變量數(shù)目，例如原來有700個樣本榛臼，也就是700個維度/變量伊佃，但實(shí)際上根據(jù)樣本中的基因表達(dá)量來歸類，它們或多或少都會有一些關(guān)聯(lián)沛善，這樣有關(guān)聯(lián)的一些樣本就會聚成一小撮航揉，表示它們的”特征距離“相近。最后直接分析這些”小撮“金刁，就用最少的變量代表了實(shí)際最多的特征

既然每個主成分都表示具有相關(guān)性的一小撮變量（稱之為”metafeature“帅涂，元特征集），那么問題來了：

降維后怎么選擇合適數(shù)量的主成分來代表整個數(shù)據(jù)集尤蛮？

七媳友、確定PCA的維數(shù)

Suerat使用了jackStraw方法[4]來估計應(yīng)該使用多少個principal components的基因來進(jìn)行下游分析。

# 注意: 這可能是一個非常費(fèi)時的過程产捞。
pbmc <- JackStraw(pbmc, num.replicate = 100)  # 重復(fù)一百次
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)     # 選擇20個PC

JackStrawPlot以及PCElbowPlot函數(shù)提供給我們兩種可視化手段來查看jackStraw的結(jié)果醇锚。使用JackStrawPlot可以查看每個PC的p value，這樣方便我們通過p值進(jìn)行篩選轧葛。PCElbowPlot可以方便的查看從哪個PC開始搂抒，差分度就進(jìn)入的平臺期艇搀。

• 1.JackStrawPlot

JackStrawPlot(pbmc, dims = 1:15)

image

• 2.PCElbowPlot : 拐點(diǎn)圖 # 最大可設(shè)置到50

  ElbowPlot(pbmc,ndims = 50, reduction = "pca")   # 最大可設(shè)置到50
  
  

它是根據(jù)每個主成分對總體變異水平的貢獻(xiàn)百分比排序得到的圖尿扯，我們主要關(guān)注”肘部“的PC，它是一個轉(zhuǎn)折點(diǎn)（也即是這里的PC9-10）焰雕，說明取前10個主成分可以比較好地代表總體變化

image

但官方依然建議我們衷笋，下游分析時多用幾個成分試試（10個、15個甚至50個）矩屁，如果起初選擇的合適辟宗，結(jié)果不會有太大變化爵赵；另外推薦開始不確定時要多選一些主成分，也不要直接就定下5個成分進(jìn)行后續(xù)分析泊脐，那樣很有可能會出錯

先擇PC空幻，對于得到正確的結(jié)果具體關(guān)鍵性的作用。除了上面的方法以外容客，還可以使用其它的一些辦法秕铛，比如說使用GSEA等富集分析的辦法來選擇PC，或者說使用一個隨機(jī)模型來查看具體哪些PC會比較有效果缩挑。但是后者這一方法在實(shí)現(xiàn)起來需要消耗過多的資源但两。

在本例中，因?yàn)槭荢eurat挑選的例子供置，所以通過上面的JackStraw方法谨湘，只要把cut.off值設(shè)置在7-10之間，就可以得到差不多的結(jié)果芥丧。

Tips: 確定進(jìn)行聚類和降維的基因和細(xì)胞后紧阔，下一步就是聚類了，在這里需要搞情況聚類和降維作用是不一樣的娄柳。特別是T-sne是降維可視化作用寓辱，而不是用它來聚類.為什么看到T-sne 圖是不同顏色的色塊呢？常見的聚類包括kmeans/hclust .

八赤拒、細(xì)胞聚類

當(dāng)決定了使用哪些PC中的基因?qū)?xì)胞進(jìn)行分類之后秫筏，就可以使用FindClusters來對細(xì)胞聚類了。分群后的結(jié)果會保存在object@ident slot中挎挖。

Seurat 3版本提供了graph-based 聚類算法这敬，參考兩篇文章：[SNN-Cliq, Xu and Su, Bioinformatics, 2015]和 [PhenoGraph, Levine et al., Cell, 2015]. 簡言之，這些graph的方法都是將細(xì)胞嵌入一個算法圖層中蕉朵，例如：K-nearest neighbor (KNN) graph 中崔涂，每個細(xì)胞之間的連線就表示相似的基因表達(dá)模式，然后按照這個相似性將圖分隔成一個個的高度相關(guān)的小群體（名叫‘quasi-cliques’ or ‘communities’）始衅；

（Step-1）在PhenoGraph中冷蚂，首先根據(jù)PCA的歐氏距離結(jié)果，構(gòu)建了KNN圖汛闸，找到了各個近鄰組成的小社區(qū)蝙茶；然后根據(jù)近鄰中兩兩住戶（細(xì)胞）之間的交往次數(shù)（shared overlap）計算每條線的權(quán)重（術(shù)語叫Jaccard similarity） 。這些計算都是用包裝好的函數(shù)FindNeighbors() 得到的诸老，它的輸入就是前面降維最終確定的主成分隆夯。

（Step-2）得到權(quán)重后，為了對細(xì)胞進(jìn)行聚類，使用了計算模塊的算法（默認(rèn)使用Louvain蹄衷，另外還有包括SLM的其他三種）忧额，使用FindClusters() 進(jìn)行聚類。其中包含了一個resolution的選項(xiàng)愧口，它會設(shè)置一個”間隔“值睦番，該值越大，間隔越大耍属，得到的cluster數(shù)量越多抡砂。一般來說，這個值在細(xì)胞數(shù)量為3000左右時設(shè)為0.4-1.2 會有比較好的結(jié)果

pbmc <- FindNeighbors(pbmc, dims = 1:10)  ## dims 作為輸入的維度
pbmc <- FindClusters(pbmc, resolution = 0.6)  ## 分辨率參數(shù)恬涧，和cluster 個數(shù)有關(guān)系注益，如何你想獲得粗糙的分類，需要調(diào)成大一些.

# 查看一下前幾個的cluster,結(jié)果聚成3類.可能發(fā)現(xiàn)新的細(xì)胞系.
head(Idents(pbmc), 5)

## AAACATACAACCAC AAACATTGAGCTAC AAACCGTGCTTCCG AAACCGTGTATGCG AAACGCTGGTTCTT 
##              4              3              0              6              4 
## Levels: 0 1 2 3 4 5 6 7 8

九溯捆、降維后丑搔，用聚類的類別可視化【UMAP/tSNE】

使用UMAP或者tSNE的好處在于將多維的PC圖降維至2維圖像中來，這樣方便分集細(xì)胞的可視化提揍。 UMAP皺縮成團(tuán)分不開啤月，而tSNE圖則分散得多。

pbmc <- RunUMAP(pbmc, dims = 1:10)
P1 <- DimPlot(pbmc, reduction = "umap", label = TRUE)
pbmc <- RunTSNE(pbmc, dims = 1:10)
P2 <- DimPlot(pbmc, reduction = "tsne",label = TRUE)
P1 + P2

image

您可以在此時保存該對象劳跃，以便可以輕松地將其裝入谎仲，而不必重新運(yùn)行上面執(zhí)行的計算密集型步驟，也可以輕松地與協(xié)作者共享該對象刨仑。

saveRDS(pbmc, file = "../output/pbmc_tutorial.rds")

得到聚類分組結(jié)果郑诺，通常需要繼續(xù)分析，在每組里面那些標(biāo)記基因杉武。由于10x genomic 對每一個細(xì)胞都進(jìn)行標(biāo)簽化辙诞，不知道每一個標(biāo)簽對應(yīng)那種細(xì)胞，可以通過每組的標(biāo)記基因?qū)γ拷M進(jìn)行細(xì)胞類型判別（根據(jù)文獻(xiàn)說明的已知marker）,進(jìn)而發(fā)現(xiàn)新的標(biāo)記基因.

十轻抱、尋找分群marker基因

生成了tSNE圖像后飞涂，我們還需要確認(rèn)這一堆堆的細(xì)胞都分別是些什么細(xì)胞，這時祈搜，可以使用特定細(xì)胞的marker來對不同的細(xì)胞群進(jìn)行標(biāo)記较店。但是想要對它們進(jìn)行標(biāo)記，首先需要得到每個集群差異表達(dá)的基因容燕。

在比較時梁呈，Seurat使用一比多的辦法。如果需要對cluster 1的細(xì)胞(ident.1)進(jìn)行差異分析時缰趋，它比較的對象將是除了cluster 1以外的其它所有的細(xì)胞捧杉。也可以使用ident.2參數(shù)來傳遞需要比較的對象。 在FindMarkers函數(shù)中秘血，min.pct是指的marker基因所占細(xì)胞數(shù)的最低比例味抖，這里使用1/4。 使用max.cells.per.ident可以讓函數(shù)對細(xì)胞進(jìn)行抽樣灰粮，以加速計算仔涩。 FindAllMarkers可以一次性找到所有的高表達(dá)的marker。
FindAllMarkers()函數(shù)粘舟，它默認(rèn)會對單獨(dú)的一個細(xì)胞群與其他幾群細(xì)胞進(jìn)行比較找到正熔脂、負(fù)表達(dá)biomarker（這里的正負(fù)有點(diǎn)上調(diào)、下調(diào)基因的意思柑肴；正marker表示在我這個cluster中表達(dá)量高霞揉，而在其他的cluster中低）。

# find all markers of cluster 1
cluster1.markers <- FindMarkers(object = pbmc, ident.1 = 1, min.pct = 0.25)
head(cluster1.markers, n = 5)
##               p_val avg_logFC pct.1 pct.2     p_val_adj
## RPS12 4.371171e-118 0.5264969     1 0.990 5.994624e-114
## RPS6  1.681008e-114 0.5234111     1 0.994 2.305334e-110
## RPS27 1.719993e-107 0.5242477     1 0.991 2.358799e-103
## RPS14 1.546910e-105 0.4500125     1 0.993 2.121432e-101
## RPS25  5.848778e-98 0.5427189     1 0.973  8.021015e-94

上面這個代碼中min.pct的意思是：一個基因在任何兩群細(xì)胞中的占比最低不能低于多少晰骑，當(dāng)然這個可以設(shè)為0适秩，但會大大加大計算量.

還可以計算cluster 5與（cluster0+cluster3）的差異基因

# find all markers distinguishing cluster 5 from clusters 0 and 3
cluster5.markers <- FindMarkers(object = pbmc, ident.1 = 5, 
                                ident.2 = c(0,3), min.pct = 0.25)
head(cluster5.markers, n = 5)
##                p_val avg_logFC pct.1 pct.2     p_val_adj
## FCGR3A 3.849393e-139  2.723391 0.981 0.098 5.279058e-135
## RHOC    1.149445e-91  1.665090 0.865 0.130  1.576350e-87
## CDKN1C  6.799344e-75  1.061123 0.519 0.022  9.324621e-71
## IFITM2  1.347125e-73  1.560680 1.000 0.644  1.847448e-69
## HES4    3.616785e-69  1.107357 0.596 0.052  4.960059e-65

計算所有類和其他類的標(biāo)記基因,（對所有cluster都比較一下，only.pos = TRUE為只挑出來正表達(dá)marker）：

pbmc.markers <- FindAllMarkers(object = pbmc, only.pos = TRUE,  
                               min.pct = 0.25, logfc.threshold = 0.25)

#挑每個cluster的前2個

pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)

## # A tibble: 18 x 7
## # Groups:   cluster [9]
##        p_val avg_logFC pct.1 pct.2 p_val_adj cluster gene    
##        <dbl>     <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>   
##  1 0\.            3.72  0.969 0.136 0\.        0       S100A8  
##  2 6.30e-298     3.77  0.996 0.236 8.63e-294 0       S100A9  
##  3 7.44e- 94     0.807 0.938 0.587 1.02e- 89 1       LDHB    
##  4 3.14e- 61     0.848 0.426 0.105 4.31e- 57 1       CCR7    
##  5 1.22e- 80     0.975 0.981 0.635 1.67e- 76 2       LTB     
##  6 5.02e- 69     0.956 0.783 0.315 6.88e- 65 2       IL7R    
##  7 0\.            2.95  0.931 0.045 0\.        3       CD79A   
##  8 5.62e-214     2.46  0.617 0.024 7.71e-210 3       TCL1A   
##  9 3.86e-164     2.39  0.973 0.22  5.30e-160 4       CCL5    
## 10 2.97e-142     2.12  0.595 0.052 4.07e-138 4       GZMK    
## 11 1.66e-171     2.34  0.981 0.135 2.27e-167 5       FCGR3A  
## 12 1.68e-106     1.90  1     0.368 2.30e-102 5       LST1    
## 13 2.24e-211     3.48  0.972 0.063 3.07e-207 6       GZMB    
## 14 1.84e-127     3.52  0.943 0.138 2.52e-123 6       GNLY    
## 15 1.11e-180     2.68  0.806 0.013 1.52e-176 7       FCER1A  
## 16 5.01e- 20     1.93  1     0.548 6.87e- 16 7       HLA-DPB1
## 17 4.08e-171     5.02  1     0.012 5.60e-167 8       PF4     
## 18 6.09e- 99     5.96  1     0.026 8.35e- 95 8       PPBP

>top2 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)
> top2$gene
 [1] "LDHB"     "CCR7"     "LTB"      "AQP3"     "S100A9"   "S100A8"   "CD79A"    "TCL1A"   
 [9] "CCL5"     "GZMK"     "FCGR3A"   "LST1"     "GZMB"     "GNLY"     "TBCB"     "NDUFA2"  
[17] "FCER1A"   "HLA-DPB1" "PF4"      "PPBP"

還可以使用不同的算法來尋找差異表達(dá)的基因硕舆。在test.use選項(xiàng)中設(shè)置（詳見：DE vignette ）例如：

cluster1.markers <- FindMarkers(object = pbmc, ident.1 = 0, 
                                thresh.use = 0.25, test.use = "roc", 
                                only.pos = TRUE)

想了解更多關(guān)于查找差異表達(dá)基因的信息秽荞，可以查看Seurat的vignette。

【一維可視化】標(biāo)記基因在不同cluster 表達(dá)量

拿到差異表達(dá)基因后抚官，可以使用VlnPlot來查看結(jié)果扬跋。

VlnPlot(object = pbmc, features = c("MS4A1", "CD79A"))

image

# you can plot raw UMI counts as well
VlnPlot(pbmc, features = c("MS4A1", "CD79A"), slot = "counts", log = TRUE)

image.png

## RidgePlot
## RidgePlot
RidgePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))

image.png

## DotPlot
DotPlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))

image.png

【二維可視化】標(biāo)記基因在不同細(xì)胞表達(dá)量

FeaturePlot(object = pbmc, 
            features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", 
                         "FCGR3A", "LYZ", "PPBP", "CD8A"), 
            cols = c("blue", "red"))

image

提取所有cluster 的top10 標(biāo)記基因可視化

top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
# 設(shè)置slim.col.label=TRUE將避免打印每個細(xì)胞的名稱，而只打印cluster的ID.
DoHeatmap(object = pbmc, features = top10$gene) + NoLegend()

image

目前我們通過聚類知道了凌节，將細(xì)胞分成8個類別钦听，但是并不知道具體是那種細(xì)胞。所有需要進(jìn)一步通過已知的特定細(xì)胞標(biāo)記基因?qū)luster 打上細(xì)胞類型標(biāo)簽.

十一倍奢、對分群進(jìn)行細(xì)胞類型標(biāo)記

對于示例數(shù)據(jù)彪见，我們已經(jīng)有一些比較常見的細(xì)胞marker了，很方便的就可以對分集的細(xì)胞進(jìn)行標(biāo)記了娱挨。

將原來的cluster 0-9 余指，修改成對應(yīng)的細(xì)胞名稱

new.cluster.ids <- c("Naive CD4 T", "Memory CD4 T", 
                     "CD14+ Mono", "B", "CD8 T", 
                     "FCGR3A+ Mono", 
                     "NK", "DC", "Platelet","NULL")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)

P1 <- DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
P2 <- DimPlot(pbmc, reduction = "tsne", label = TRUE, pt.size = 0.5) + NoLegend()
P1 + P2

image.png

saveRDS(pbmc, file = "./output/pbmc3k_final.rds")

當(dāng)我們修改了聚類的參數(shù)，結(jié)果會不同跷坝。比如分辨率參數(shù).

# 生成一個快照, 其實(shí)就是把pbmc@ident中的數(shù)據(jù)拷貝到pbmc@meta.data中酵镜。
# 當(dāng)我們想把pbmc@meta.data中的值賦值給pbmc@ident的話，可以使用SetAllIdent函數(shù)柴钻。
# 比如pbmc <- SetAllIdent(object=pbmc, id = 'orig.ident')
pbmc <- StashIdent(object = pbmc, save.name = "ClusterNames_0.6")

# 現(xiàn)在我們使用不同的resolution來重新分群淮韭。當(dāng)我們提高resolution的時候，
# 我們可以得到更多的群贴届。
# 之前我們在FindClusters中設(shè)置了save.snn=TRUE靠粪，所以可以不用再計算SNN了蜡吧。
# 下面就直接提高一下區(qū)分度，設(shè)置resolution = 0.8
pbmc <- FindClusters(pbmc, resolution = 0.8)

番外

分辨率左圖是0.8占键，右圖0.6

# 并排畫兩個tSNE,右邊的是用ClusterNames_0.6 (resolution=0.6) 來分組顏色的昔善。
# 我們可以看到當(dāng)resolution提高之后，CD4 T細(xì)胞被分成了兩個亞群畔乙。
plot1 <- DimPlot(object = pbmc, reduction = "umap", label = TRUE) + NoLegend()
plot2 <- DimPlot(object = pbmc, reduction = "umap", label = TRUE,
                 group.by = "ClusterNames_0.6") +NoLegend()
plot1 + plot2

image

類別不同君仆，我們又可以按照這一次聚類結(jié)果，尋找marker基因

# 我們再一次尋找不同分集中的marders牲距。
tcell.markers <- FindMarkers(object = pbmc, ident.1 = 0, ident.2 = 1)

# Most of the markers tend to be expressed in C1 (i.e. S100A4). 
# However, we can see that CCR7 is upregulated in 
# C0, strongly indicating that we can differentiate memory from naive CD4 cells.
# cols.use demarcates the color palette from low to high expression
FeaturePlot(object = pbmc, features = c("S100A4", "CCR7"), 
            cols = c("blue", "red"))

image