一、設置Seurat對象
使用的示例數(shù)據(jù)集來自10X Genome 測序的 PBMC。
下載鏈接：https://s3-us-west-2.amazonaws.com/10x.files/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz
三個文件是“barcodes.tsv","features.tsv","matrix.mtx";
barcodes.tsv就是cell id
features.tsv就是gene id
matrix.mtx就是計數(shù)counts矩陣（UMI數(shù)量）
Read10X（）函數(shù)從10X讀取cellranger管道的輸出，返回唯一的分子識別（UMI）計數(shù)矩陣。該矩陣中的值表示在每個細胞（列）中檢測到的每個特征（即基因；行）的分子數(shù)。
【每個基因的UMI種類代表基因表達絕對量】
接下來，使用計數(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數(shù)據(jù)組
pbmc.data <- Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/")
# Initialize the Seurat object with the raw (non-normalized data).
#使用原始(非標準化數(shù)據(jù))初始化Seurat對象
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 在1次分析中對2,700個樣本進行了13714項特征（基因）分析
## Active assay: RNA (13714 features, 0 variable features)
計數(shù)矩陣中的數(shù)據(jù)長什么樣子？看一下某3行（3個基因）的前30行（前30個細胞）
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盡可能使用稀疏矩陣表示。這大大節(jié)省了Drop-Seq/inDrop/10x數(shù)據(jù)的內(nèi)存和速度。
二、標準預處理流程
1、基于質(zhì)量控制指標進行細胞過濾
2、數(shù)據(jù)標準化和縮放
3、檢測高度可變特征
1、QC和選擇細胞進行下一步分析
使用Seurat，可以瀏覽QC指標，并根據(jù)定義的標準過濾細胞。常用的一些質(zhì)量控制指標包括
1)每個細胞中檢測到的基因數(shù)nfeature（太少和太多都不要）
*低質(zhì)量細胞或空滴通常只有很少的基因
*細胞二倍體或多重體可能表現(xiàn)出異常高的基因計數(shù)
2）每個細胞內(nèi)檢測到的分子總數(shù)nCount(UMI總數(shù))：標準和上面類似
3）線粒體基因組的read比例
*低質(zhì)量/瀕死的細胞通常表現(xiàn)出廣泛的線粒體污染，出現(xiàn)大比例線粒體基因read
*將以MT-開頭的所有基因定義為線粒體基因
*PercentageFeatureSet()函數(shù)計算線粒體QC指標（計算線粒體基因的比例）
# The [[ operator can add columns to object metadata. This is a great place to stash QC stats
#[[]]運算符可以向?qū)ο笤獢?shù)據(jù)添加列?？梢杂脕泶娣臦C統(tǒng)計數(shù)據(jù)
#計算線粒體read百分比
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
QC指標存儲在Seurat的什么地方？
在CreateSeuratObject()過程中會自動計算獨特基因的數(shù)量和分子總數(shù)?？梢栽趯ο髆eta.data中找到它們
# Show QC metrics for the first 5 cells 顯示前5個細胞的QC指標
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
nCount_RNA：每個細胞的UMI總數(shù)
【關(guān)于UMI的解釋，詳見單細胞測序?qū)W習（二） - —小鄭同學— - 博客園 (cnblogs.com)】
*1個基因有幾個UMI就有幾個表達量（因為PCR擴增有偏好，用UMI來表示絕對RNA表達量）
nFeature_RNA：每個細胞所檢測到的基因數(shù)目
一個基因可以對應多個UMI，所以nFeature<nCount
接下來，將QC指標可視化，并使用這些指標來過濾細胞。
*過濾nfeature計數(shù)（基因）大于2500或小于200的細胞
*過濾線粒體計數(shù)大于5%的細胞
# Visualize QC metrics as a violin plot
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
ncol代表用幾列來顯示圖像
# FeatureScatter is typically used to visualize feature-feature relationships, but can be used
# for anything calculated by the object, i.e. columns in object metadata, PC scores etc.
# FeatureScatter通常用于可視化特征-特征關(guān)系，但可用于對象計算的任何內(nèi)容，即對象元數(shù)據(jù)中的列、PC分數(shù)等。
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
過濾細胞
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
2、標準化數(shù)據(jù)： NormalizeData()
【關(guān)于什么是標準化，詳見數(shù)據(jù)預處理之Normalization - 知乎 (zhihu.com)
對數(shù)據(jù)集進行一些變化，把數(shù)據(jù)拉到一個特定范圍里，變得更有統(tǒng)計意義，有的基因表達太高，會影響整體的分析】
默認情況下，使用全局縮放歸一化方法“LogNormalize”，該方法用總表達式標準化每個細胞的特征表達式測量值，再乘以比例因子(默認情況下為10,000)，然后對結(jié)果進行對數(shù)轉(zhuǎn)換。標準化值存儲在PBMC[[\“RNA\”]]@data中。【
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
標準化的方式和比例因子的默認值可以不寫，直接寫這個代碼也是一樣的
pbmc <- NormalizeData(pbmc)
3、識別差異基因：FindVariableFeature()
接下來，計算在數(shù)據(jù)集中表現(xiàn)出細胞間高差異的特征（基因）子集(即，它們在某些細胞中高表達，而在其他細胞中低表達)。在下游分析中關(guān)注這些基因有助于突出單細胞數(shù)據(jù)集中的生物信號。
在Seurat中的過程在這里進行了詳細描述，通過直接對單個細胞數(shù)據(jù)中固有的均值-方差關(guān)系進行建模，改進了以前的版本，并在FindVariableFeature()函數(shù)中實現(xiàn)。默認情況下，為每個數(shù)據(jù)集返回2,000個要素。這些將用于下游分析，如PCA。
#識別變異最大的前2000個基因
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
# Identify the 10 most highly variable genes 找出10個變異最大的基因
top10 <- head(VariableFeatures(pbmc), 10)
# plot variable features with and without labels 繪制帶標簽和不帶標簽的可變基因
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2
points=要標出名字的基因集
repel=true可以解決標簽重疊的問題
4、縮放數(shù)據(jù)：ScaleData()
接下來，應用線性變換(“縮放”)，可以理解為PCA降維之前的標準預處理步驟。ScaleData()函數(shù)的作用是：
*改變每個基因的表達，使每個細胞基因的平均表達值為0
*縮放每個基因的表達，使細胞間的方差為1：
**這一步在下游分析中給予同等的權(quán)重，這樣高表達的基因就不會占據(jù)主導地位?！痉讲罘磻牟煌虮磉_的離散程度】
*其結(jié)果存儲在PBMC[[\“RNA\”]]@scale.data中
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
5、線性降維：RunPCA()
接下來，對縮放后的數(shù)據(jù)執(zhí)行PCA。默認情況下，只有先前確定的可變要素用作輸入，但如果想選擇不同的子集，則可以使用Feature參數(shù)進行定義。
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
3種方法來可視化PCA的細胞和基因，包括VizDimReduce()、DimPlot()和DimHeatmap()
# Examine and visualize PCA results a few different ways通過幾種不同的方式檢查和可視化PCA結(jié)果
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
VizDimReduce()
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
DimPlot()
DimPlot(pbmc, reduction = "pca")
DimHeatmap()
DimHeatmap()可以查看數(shù)據(jù)中異質(zhì)性的主要來源，并可以確定哪些PCA維度可以用于下一步的下游分析。細胞和特征根據(jù)PCA的分數(shù)進行排序。
DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)
6、確定數(shù)據(jù)集的維度(PC=？)
為了克服在單細胞數(shù)據(jù)中在單個特征中的技術(shù)噪音，Seurat 聚類細胞是基于PCA分數(shù)的。每個PC代表著一個‘元特征’（帶有跨相關(guān)特征集的信息）。因此，最主要的主成分代表了壓縮的數(shù)據(jù)集。問題是要選多少PC呢？
方法一：JackStrawPlot
作者受JackStraw procedure 啟發(fā)。隨機置換數(shù)據(jù)的一部分子集（默認1%）再運行PCA，構(gòu)建了一個’null distribution’的特征分數(shù)，重復這一步。最終會識別出低P-value特征的顯著PCs
# 注意：對于大數(shù)據(jù)集，此過程可能需要很長時間，為方便起見，請將其注釋掉。
# 可以使用更多的#近似技術(shù)(如在ElBowPlot()中實現(xiàn)的技術(shù))來減少#計算時間。
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
JackStrawPlot()函數(shù)提供了一個可視化工具，用于比較具有均勻分布(虛線)的每臺PC的p值分布。“重要的”PC將顯示大量低p值的功能(虛線上方的實心曲線)。在這種情況下，在最初的10-12臺PC之后，重要性似乎急劇下降。
繪圖看看
JackStrawPlot(pbmc, dims = 1:15)
方法二：ElbowPlot
“ElbowPlot”：基于每個分量所解釋的方差百分比對主要成分進行排名。 在此示例中，我們可以在PC9-10周圍觀察到“elbow ”，這表明大多數(shù)真實信號是在前10臺PC中捕獲的。
ElbowPlot(pbmc)
為了識別出數(shù)據(jù)的真實維度，有三種方法：
1、用更加受監(jiān)督的方法來確定PCs的異質(zhì)性，比如可以結(jié)合GSEA來分析
2、基于隨機空模型的統(tǒng)計測試，但是對于大型數(shù)據(jù)集來說很費時間，并且可能不會返回一個清晰的 pc 連接。
3、常用的啟發(fā)式算法，可以立即計算出來。
在這個例子中三種方法均產(chǎn)生了相似的結(jié)果，以PC 7-12作為閾值。
這個例子中，作者選擇10，但是實際過程中還要考慮：
1、樹突狀細胞和NK細胞可能在PCs12和13中識別，這可能定義了罕見的免疫亞群（比如，MZB1是漿細胞樣DC的標記）。但是除非有一定的知識量，否則很難從背景噪音中發(fā)現(xiàn)。
2、用戶可以選擇不同的PCs再進行下游分析，比如選10，15，50等。結(jié)果常常差不多
3、建議在選擇該參數(shù)時候，盡量偏高一點。如果僅僅使用5PCs會對下游分析產(chǎn)生不利影響
7、聚類：FindClusters()
pbmc <- FindNeighbors(pbmc, dims = 1:10) #將之前定義的數(shù)據(jù)集維度（前10個PC作為輸入）
pbmc <- FindClusters(pbmc, resolution = 0.5)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 2638
## Number of edges: 95965
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8723
## Number of communities: 9
## Elapsed time: 0 seconds
FindClusters()函數(shù)包含一個分辨率參數(shù)，該參數(shù)設置下游集群的“粒度”，值越大，集群數(shù)量越多。將此參數(shù)設置在0.4-1.2之間通常會為大約3K個細胞的數(shù)據(jù)集返回良好的結(jié)果。對于較大的數(shù)據(jù)集，最佳分辨率通常會增加。
可以使用idents()函數(shù)找到聚類亞群ID。
#查看前5個細胞的聚類亞群ID
head(Idents(pbmc), 5)
## AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1
## 2 3 2 1
## AAACCGTGTATGCG-1
## 6
## Levels: 0 1 2 3 4 5 6 7 8
8、非線性降維(UMAP/tSNE)
Seurat提供了幾種非線性降維技術(shù)，如tSNE和UMAP，以可視化和探索這些數(shù)據(jù)集。這些算法的目標是了解數(shù)據(jù)的底層流形，以便將相似的細胞放在低維空間中。
上述基于圖形的聚類中的細胞應該在這些降維圖上共同定位。
作為UMAP和tSNE的輸入，建議使用相同的pc作為聚類分析的輸入。
umap
# 安裝umap包：reticulate::py_install(packages ='umap-learn')
pbmc <- RunUMAP(pbmc, dims = 1:10)
#可以設置`LABEL=TRUE`或使用LabelClusters函數(shù)來幫助標記單個群集
DimPlot(pbmc, reduction = "umap")
保存
saveRDS(pbmc, file = "../output/pbmc_tutorial.rds")
tsne
pbmc <- RunTSNE(pbmc, dims = 1:10)
DimPlot(pbmc, reduction = "tsne")
# 顯示聚類標簽
DimPlot(pbmc, reduction = "tsne", label = TRUE)
9、尋找亞群間差異表達基因：FindMarkers()
FindMarkers（對象，ident.1=指定亞群，Min.pct=，Thresh.test=）
Min.pct要求在兩組亞群的任一組中以最小百分比檢測到特征
【兩個亞群中，至少有一個亞群有至少Min.pct的細胞有這個表達特征】
Thresh.test要求在兩組細胞之間以一定數(shù)量(平均)差異表達特征。
【兩個亞群中，這個表達特征的平均差異要達到Thresh.test】
可以將這兩個值都設置為0，但會大大增加時間-因為這將測試大量不太可能有很大差別的特性。作為加速這些計算的另一個選項，可以設置max.cells.per.ident。這將對每個身份類進行下采樣，使其細胞數(shù)不超過設置的值。雖然通常會失去動力，但速度的提高可能會很顯著，最具差異性的功能可能仍然會上升到頂端。
發(fā)現(xiàn)亞群2的所有marker
cluster2.markers <- FindMarkers(pbmc, ident.1 = 2, min.pct = 0.25)
head(cluster2.markers, n = 5)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## IL32 2.593535e-91 1.2154360 0.949 0.466 3.556774e-87
## LTB 7.994465e-87 1.2828597 0.981 0.644 1.096361e-82
## CD3D 3.922451e-70 0.9359210 0.922 0.433 5.379250e-66
## IL7R 1.130870e-66 1.1776027 0.748 0.327 1.550876e-62
## LDHB 4.082189e-65 0.8837324 0.953 0.614 5.598314e-61
發(fā)現(xiàn)區(qū)分亞群5和亞群0、3的marker
cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
head(cluster5.markers, n = 5)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## FCGR3A 2.150929e-209 4.267579 0.975 0.039 2.949784e-205
## IFITM3 6.103366e-199 3.877105 0.975 0.048 8.370156e-195
## CFD 8.891428e-198 3.411039 0.938 0.037 1.219370e-193
## CD68 2.374425e-194 3.014535 0.926 0.035 3.256286e-190
## RP11-290F20.3 9.308287e-191 2.722684 0.840 0.016 1.276538e-186
發(fā)現(xiàn)每個亞群和其他所有細胞相比的陽性marker，每個亞群顯示兩個
# ones
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
pbmc.markers %>%
group_by(cluster) %>%
slice_max(n = 2, order_by = avg_log2FC)
## # A tibble: 18 × 7
## # Groups: cluster [9]
## p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
## <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <chr>
## 1 1.17e- 83 1.33 0.435 0.108 1.60e- 79 0 CCR7
## 2 1.74e-109 1.07 0.897 0.593 2.39e-105 0 LDHB
## 3 0 5.57 0.996 0.215 0 1 S100A9
## 4 0 5.48 0.975 0.121 0 1 S100A8
## 5 7.99e- 87 1.28 0.981 0.644 1.10e- 82 2 LTB
## 6 2.61e- 59 1.24 0.424 0.111 3.58e- 55 2 AQP3
## 7 0 4.31 0.936 0.041 0 3 CD79A
## 8 9.48e-271 3.59 0.622 0.022 1.30e-266 3 TCL1A
## 9 4.93e-169 3.01 0.595 0.056 6.76e-165 4 GZMK
## 10 1.17e-178 2.97 0.957 0.241 1.60e-174 4 CCL5
## 11 3.51e-184 3.31 0.975 0.134 4.82e-180 5 FCGR3A
## 12 2.03e-125 3.09 1 0.315 2.78e-121 5 LST1
## 13 6.82e-175 4.92 0.958 0.135 9.36e-171 6 GNLY
## 14 1.05e-265 4.89 0.986 0.071 1.44e-261 6 GZMB
## 15 1.48e-220 3.87 0.812 0.011 2.03e-216 7 FCER1A
## 16 1.67e- 21 2.87 1 0.513 2.28e- 17 7 HLA-DPB1
## 17 3.68e-110 8.58 1 0.024 5.05e-106 8 PPBP
## 18 7.73e-200 7.24 1 0.01 1.06e-195 8 PF4
Seurat有幾個差異表達測試，可以使用test.use參數(shù)進行設置(有關(guān)詳細信息，請參閱我們的DE Vignette)。例如，ROC測試返回任何單個標記的“分類能力”(范圍從0-隨機到1-完美)。
cluster0.markers <- FindMarkers(pbmc, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)
可視化marker表達的工具
VlnPlot()顯示亞群間的某基因表達分布
FeaturePlot()在tSNE或PCA圖上可視化某基因表達分布
這兩個是最常用的可視化工具。還建議研究RidgePlot()、CellScatter()和DotPlot()作為查看數(shù)據(jù)集的其他方法。
VlnPlot()顯示亞群間的某基因表達分布
VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
可視化原始的count
VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)
FeaturePlot()在tSNE或PCA圖上可視化某基因表達分布
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP",
"CD8A"))
DoHeatmap()為給定的細胞和基因生成表達式熱圖。在本例中，繪制每個群集的前10個標記(如果少于10個，則繪制所有標記)。
pbmc.markers %>%
group_by(cluster) %>%
top_n(n = 10, wt = avg_log2FC) -> top10
DoHeatmap(pbmc, features = top10$gene) + NoLegend()
10、注釋
使用規(guī)范標記將無偏聚類與已知細胞類型相匹配：
new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T", "FCGR3A+ Mono",
"NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
pt.size（字號）
保存
saveRDS(pbmc, file = "../output/pbmc3k_final.rds")
