單細(xì)胞學(xué)習(xí)入門教程，里面的東西可以反復(fù)琢磨，標(biāo)準(zhǔn)分析基本上就這么多內(nèi)容。教程中提到了大量的資源介紹，參數(shù)含義解讀，甚至高級部分的簡潔版算法介紹，結(jié)果解讀，可視化方法，以及數(shù)據(jù)變換后使用什么函數(shù)來提取對應(yīng)的變換后的數(shù)據(jù)。
1、數(shù)據(jù)介紹
此次分析所有數(shù)據(jù)來源于Peripheral Blood Mononuclear Cells (PBMC)，使用10X Genomics， Illumina NextSeq 500進(jìn)行測序，共獲得2700個單細(xì)胞?？梢栽谶@里進(jìn)行下載：
[https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz]
數(shù)據(jù)下載下來后進(jìn)行解壓，主要包含以下三個文件：
2、安裝Seurat包
安裝Seurat包
install.packages('Seurat')
library('Seurat')
packageVersion("Seurat")
相關(guān)包下載
install.packages(c('dplyr','patchwork'))
3、加載包和讀取數(shù)據(jù)
library(dplyr)
library(Seurat)
library(patchwork)
# Load the PBMC dataset, 讀取數(shù)據(jù)主要使用Read10X這個函數(shù)
pbmc.data <- Read10X(data.dir = "data/filtered_gene_bc_matrices/hg19")
# Initialize the Seurat object with the raw (non-normalized data).
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)
count matrix長什么樣？
# 首先看看三個基因的count值
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，即no molecules detected。當(dāng)然，這個地方還有另外一種含義就是這個基因是真的沒有表達(dá)。
由于單細(xì)胞測序數(shù)據(jù)中大多數(shù)的值都為0，因此，seurat使用一個稀疏矩陣來保存測序得到的count matrix，這樣有利于數(shù)據(jù)存儲空間的節(jié)省。
我們來看看使用稀疏矩陣和使用0來存儲兩種方式的大小對比
dense.size <- object.size(as.matrix(pbmc.data))
dense.size
sparse.size <- object.size(pbmc.data)
sparse.size
dense.size/sparse.size
23.7 bytes
dense為轉(zhuǎn)換為0后的matrix存儲大小，709591472 bytes
as.matrix(pbmc.data)[c("CD3D", "TCL1A", "MS4A1"), 1:30]
AAACATACAACCAC AAACATTGAGCTAC AAACATTGATCAGC AAACCGTGCTTCCG
CD3D 4 0 10 0
TCL1A 0 0 0 0
MS4A1 0 6 0 0
AAACCGTGTATGCG AAACGCACTGGTAC AAACGCTGACCAGT AAACGCTGGTTCTT
CD3D 0 1 2 3
TCL1A 0 0 0 0
MS4A1 0 0 0 0
AAACGCTGTAGCCA AAACGCTGTTTCTG AAACTTGAAAAACG AAACTTGATCCAGA
CD3D 1 0 0 2
TCL1A 1 0 0 0
MS4A1 1 1 1 0
AAAGAGACGAGATA AAAGAGACGCGAGA AAAGAGACGGACTT AAAGAGACGGCATT
CD3D 7 1 0 0
TCL1A 0 0 0 0
MS4A1 0 0 0 0
4、單細(xì)胞數(shù)據(jù)分析預(yù)處理
預(yù)處理主要包括基于QC指標(biāo)的細(xì)胞和基因過濾，數(shù)據(jù)標(biāo)準(zhǔn)化和歸一化，高變基因選擇。
4.1首先是QC來篩選高質(zhì)量的細(xì)胞
一般篩選條件有三個：
1.每個細(xì)胞中檢測到的唯一基因數(shù)
低質(zhì)量的細(xì)胞或者空的droplet液滴通常含有很少的基因
Cell doublets雙胞體或多胞體含有很高的異常的gene counts
2.線粒體基因含量比例
低質(zhì)量或者死亡細(xì)胞含有很高的線粒體基因
MT-開頭的基因認(rèn)為是線粒體基因
# The [[ operator can add columns to object metadata. This is a great place to stash QC stats
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
# 查看QC指標(biāo)
# Show QC metrics for the first 5 cells
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
過濾前三個指標(biāo)可視化
# Visualize QC metrics as a violin plot
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
每個特征之間的相關(guān)性
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)
4.2數(shù)據(jù)標(biāo)準(zhǔn)化
過濾了不想要的cells之后，下一步就是標(biāo)準(zhǔn)化。標(biāo)準(zhǔn)化方法為：LogNormalize。標(biāo)準(zhǔn)化的數(shù)據(jù)存在pbmc[["RNA"]]@data中。
歸一化和標(biāo)準(zhǔn)化的區(qū)別
首先去除樣本/細(xì)胞效應(yīng)：因為不同樣本或者細(xì)胞的測序數(shù)據(jù)量不一樣，所以同樣的一個基因在不同細(xì)胞，哪怕看到的表達(dá)量是一樣的，但是背后細(xì)胞整體測序數(shù)據(jù)量的差異其實反而說明了這個基因在不同細(xì)胞表達(dá)量其實是有差異的。
然后去除基因效應(yīng)，這個主要是在繪制熱圖的時候會需要使用，因為個別基因表達(dá)量超級高，在熱圖里面一枝獨秀，實際上我們并不會關(guān)心不同基因的表達(dá)量高低，我們僅僅是想看指定基因在不同細(xì)胞的高低而已，這樣的話，就把該基因的表達(dá)量在不同細(xì)胞的數(shù)值，進(jìn)行z-score這樣的轉(zhuǎn)換。這樣處理后的表達(dá)矩陣，就可以進(jìn)行后續(xù)的降維聚類分群。
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
Performing log-normalization
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
4.3識別高變基因
高表達(dá)基因：在一些細(xì)胞中高表達(dá)，在另外的細(xì)胞中低表達(dá)。
使用均值與方差之間的關(guān)系，來挑選高變基因，便于進(jìn)行特征選擇篩選生物學(xué)相關(guān)基因，默認(rèn)返回前2000個高變基因進(jìn)入下游分析，如PCA。
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
# Identify the 10 most highly variable genes
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
4.4識別高變基因
在數(shù)據(jù)降維之前，我們使用[ScaleData()進(jìn)行數(shù)據(jù)歸一化。歸一化的功能：使得每一個基因在所有cell中的表達(dá)均值為0，方差為1。
# 默認(rèn)只對2千個高變基因進(jìn)行數(shù)據(jù)歸一化
pbmc <- ScaleData(pbmc)
# 也可以使用全部的基因進(jìn)行歸一化，但是耗時會比較久
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
使用2000個高變基因和所有的基因進(jìn)行歸一化對后續(xù)的PCA和聚類分析并不會有太大影響，有影響的地方是后面繪制熱圖的時候 DoHeatmap()，因為繪制熱圖需要輸入歸一化后的基因表達(dá)值。
4.5識別高變基因
接下來，我們對歸一化后的數(shù)據(jù)進(jìn)行PCA分析。默認(rèn)的，只是用前面決定的高變基因進(jìn)行PCA分析，也可以使用features參數(shù)設(shè)置用戶自己選擇的基因進(jìn)行PCA分析。
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
# Examine and visualize PCA results a few different ways
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
DimPlot(pbmc, reduction = "pca")
此外，[DimHeatmap()]允許我們探索數(shù)據(jù)中的最初的異質(zhì)性，然后決定下游使用多少個PC進(jìn)行分析。
DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)
4.6 確定數(shù)據(jù)的維度
這里我們需要選擇出主成分的數(shù)目，用于后續(xù)細(xì)胞分類。這里定義的“維度”并不代表細(xì)胞類型的數(shù)目，而是對細(xì)胞分類時需要用到的一個參數(shù)。
Seurat對細(xì)胞進(jìn)行聚類主要基于他們的PCA打分，每一個PC代表一個綜合特征，它綜合了數(shù)據(jù)中相關(guān)基因表達(dá)的一些信息。前幾的主成分代表了一個數(shù)據(jù)集的穩(wěn)定的綜合的信息。那么我們應(yīng)該選擇多少個PC進(jìn)行分析呢？
我們使用了一個樣本隨機(jī)擾動的方法：the JackStraw procedure
我們隨機(jī)選擇一個數(shù)據(jù)的subset（默認(rèn)為數(shù)據(jù)得1%）進(jìn)行PCA分析，返回PCA打分，然后重復(fù)這個過程，形成一個打分的零分布。
# NOTE: This process can take a long time for big datasets, comment out for expediency.
# More approximate techniques such as those implemented in ElbowPlot() can be used to
# reduce computation time
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
# 對PC打分進(jìn)行可視化
JackStrawPlot(pbmc, dims = 1:15)
還可以根據(jù)PC的貢獻(xiàn)度進(jìn)行排序繪制一個散點圖,碎石圖
ElbowPlot(pbmc)
根據(jù)上圖，我們可以看到在第10個PC處有一個拐點。
4.7 聚類
我們使用k最近鄰法來對細(xì)胞進(jìn)行聚類（K-nearest neighbor (KNN) graph），具有相似基因表達(dá)模式的細(xì)胞之間繪制邊，然后嘗試將這張圖劃分為高度相互關(guān)聯(lián)的“準(zhǔn)派系”或“群體”。
PhenoGrap中，我們基于PCA空間的歐氏距離，構(gòu)建一個KNN圖，然后根據(jù)兩個細(xì)胞的局部鄰居的共享重疊部分(Jaccard相似性)來優(yōu)化它們之間的邊權(quán)值。這一步使用FindNeighbors()函數(shù)，使用前面決定的10個PCs作為輸入。
為了對細(xì)胞聚類，我們接下來使用模塊化優(yōu)化技術(shù)如the Louvain algorithm (default) or SLM [SLM, Blondel et al., Journal of Statistical Mechanics]，迭代地將細(xì)胞分組。這一步主要使用FindClusters()，其中有個函數(shù)resolution會隨著值變大，分類數(shù)變多。對于單細(xì)胞數(shù)據(jù)3k個細(xì)胞，resolution值為0.4-1.2就可以返回一個比較好的聚類結(jié)果。對于較大的數(shù)據(jù)集，最佳分辨率通常會增加。
pbmc <- FindNeighbors(pbmc, dims = 1:10)
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: 96033
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8720
Number of communities: 9
Elapsed time: 0 seconds
這里聚成了9個類。
# Look at cluster IDs of the first 5 cells
head(Idents(pbmc), 5)
4.8 非線性降維(UMAP/tSNE)
Seurat提供幾個非線性降維方法，如tSNE和UMAP，來可視化和探索這些數(shù)據(jù)集。這些算法的目標(biāo)是學(xué)習(xí)數(shù)據(jù)的底層流形，以便將相似的細(xì)胞放在低維空間中。
# If you haven't installed UMAP, you can do so via reticulate::py_install(packages =
# 'umap-learn')
pbmc <- RunUMAP(pbmc, dims = 1:10)
# note that you can set `label = TRUE` or use the LabelClusters function to help label
# individual clusters
DimPlot(pbmc, reduction = "umap", label=T)
保存結(jié)果下次使用
dir.create("output")
saveRDS(pbmc, file = "output/pbmc_tutorial.rds")
4.9 差異表達(dá)分析
這一步主要使用FindMarkers和[FindAllMarkers]
# find all markers of cluster 1
cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)
head(cluster1.markers, n = 5)
# find all markers distinguishing cluster 5 from clusters 0 and 3
cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
head(cluster5.markers, n = 5)
# find markers for every cluster compared to all remaining cells, report only the positive ones
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)
cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)
可視化差異表達(dá)基因：VlnPlot、RidgePlot()、CellScatter()，
VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
image-20210404151249117.png
# you can plot raw counts as well
VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
image-20210404151855546.png
這里作者展示的marker基因都是根據(jù)他的背景知識來的，他事先就知道這些個基因是哪些細(xì)胞類型里面特異表達(dá)的基因，可以用來鑒定細(xì)胞類型。
如果是我們，可能就比較困難，需要看大量的文獻(xiàn)，或者借助已知marker的數(shù)據(jù)庫。
繪制9個cluster的差異top20的基因的熱圖
top20 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 20, wt = avg_log2FC)
DoHeatmap(pbmc, features = top20$gene) + NoLegend()
3.10 細(xì)胞類型注釋
這個地方也是整個單細(xì)胞數(shù)據(jù)分析里面比較難的部分，有很多種方法可以進(jìn)行細(xì)胞類型注釋。這里作者采用了已知的canonical markers進(jìn)行注釋。 NovoCM輔助細(xì)胞定義+SingleR細(xì)胞定義。
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()
保存最終的結(jié)果
saveRDS(pbmc, file = "output/pbmc3k_final.rds")
總結(jié)
綜合以上內(nèi)容，去掉可視化部分，標(biāo)準(zhǔn)分析步驟有
pbmc.counts <- Read10X(data.dir = "pbmc3k/filtered_gene_bc_matrices/hg19/")
pbmc <- CreateSeuratObject(counts = pbmc.counts)
pbmc <- NormalizeData(object = pbmc)
pbmc <- FindVariableFeatures(object = pbmc)
pbmc <- ScaleData(object = pbmc)
pbmc <- RunPCA(object = pbmc)
pbmc <- FindNeighbors(object = pbmc)
pbmc <- FindClusters(object = pbmc)
pbmc <- RunTSNE(object = pbmc)
pbmc <- RunUMAP(object = pbmc)
DimPlot(object = pbmc, reduction = "tsne")
