Seurat 指导聚类教程

参照官网教程 用了自己的一批真实的数据，总共有7038个细胞。以下是cellranger count跑出来的标准结果。

我们从读取数据开始。Read10X函数从10X读取cellranger流程的输出，返回UMI计数矩阵。矩阵中的值表示每个特征(即基因;在每个细胞(列)中检测到的。

接下来，我们使用count矩阵创建一个Seurat对象。该对象充当一个容器，其中包含单细胞数据集的数据(如计数矩阵)和分析(如PCA或聚类结果)。

library(Seurat)
library(dplyr)

读取数据：

real_10x_data = Read10X(data.dir = "\\example_G48E2L1\\filtered_feature_bc_matrix")
Read10X(data.dir = NULL, gene.column = 2, unique.features = TRUE)

data.dir 参数包含矩阵的目录。包含matrix.mtx，gene.tsv（或features.tsv）和barcodes.tsv。为了加载多个数据目录，可以给出一个向量或命名向量。如果给定了命名向量，则cell barcode 名称将以该名称为前缀。

gene.column 参数指定基因在哪一列。features.tsv或gene.tsv用于基因名称的tsv;默认是2，表示第二列是基因名，我们来看一下features.tsv，包含3列：

unique.features 参数默认为TRUE，表示 使features name unique。

如果features.csv表明数据具有多个数据类型，则返回一个包含每种类型数据的稀疏矩阵的列表。否则将返回一个包含表达式数据的稀疏矩阵。

使用原始数据(非规范化数据)初始化Seurat对象。

CreateSeuratObject(counts, project = "SeuratProject", assay = "RNA",
  min.cells = 0, min.features = 0, names.field = 1,
  names.delim = "_", meta.data = NULL)

注意:在以前的版本(<3.0)中，该函数还接受一个参数来设置“检测到的”特征(基因)的表达阈值。为了简化初始化过程/假设，删除了此功能。如果您仍然希望为特定的数据集设置这个阈值，那么只需在调用此函数之前对输入表达式矩阵进行筛选即可。

> G48E2L1 <-CreateSeuratObject(counts = real_10x_data, project = "G48E2L1",)
> G48E2L1
An object of class Seurat
33538 features across 7038 samples within 1 assay
Active assay: RNA (33538 features)

可以发现7038samples 和网页版报告一致，33538 features也和features.csv数量一致。

count matrix 数据长什么？

例如，计数矩阵存储在G48E2L1[["RNA"]]@counts中。

只看几个基因：

>real_10x_data[c("CD3D", "TCL1A","MS4A1"),1:30]
3 x 30 sparse Matrix of class "dgCMatrix"
   [[ suppressing 30 column names ‘AAACCTGAGGAATGGA’, ‘AAACCTGAGGCTCAGA’, ‘AAACCTGCAGCCTATA’ ... ]]
CD3D  6 . 7 3 3 12 5 1 6 1 . 5 2 5 7 3 4 14 4 . 1 4 2 1 7 2 3 5 13 7
TCL1A . . . . .  . . . . . . . . . . . .  . . . . . . . . . . .  . .
MS4A1 . . . . .  . . . . . . . . . . . .  . . . . . . . . . . .  . .

点. 矩阵中的值表示0(未检测到分子)。由于scRNA-seq矩阵中的大多数值都是0，所以Seurat尽可能使用稀疏矩阵表示。这为Drop-seq/inDrop/10x数据节省了大量内存和速度。

如果需要查看稀疏矩阵的空间大小（个人理解），这些可以忽略。

> dense.size = object.size(as.matrix(real_10x_data))
> dense.size
1891208224 bytes
> spare.size <- object.size(real_10x_data)
> spare.size
179737888 bytes
> dense.size/spare.size
10.5 bytes

标准的预处理流程：

以下步骤包含Seurat中scRNA-seq数据的标准预处理工作流。这些代表细胞的选择和过滤基于QC指标，数据归一化和缩放，并检测高度可变的特征。

QC和选择细胞进行进一步分析
Seurat允许您轻松地探索QC指标，并根据任何用户定义的标准过滤单元格。大家通常使用的一些QC指标包括

• 在每个细胞中检测到的独特基因的数量
  * 低质量的细胞或空液滴通常只有很少的基因
  * 细胞双重或多胞胎可能表现出异常高的基因计数
• 同样，在一个细胞内检测到的分子总数(与独特的基因密切相关)
• reads map到线粒体基因组的比例
  * 低质量/濒死细胞常表现出广泛的线粒体污染
  * 我们使用PercentageFeatureSet函数计算线粒体QC指标，该函数计算来自一组特征的计数的百分比
  * 我们使用所有以MT-开头的基因作为一组线粒体基因

"[[" 操作符可以向对象元数据添加列。新增一列储存QC统计最好不过了。
G48E2L1[["percent.mt"]] <- PercentageFeatureSet(G48E2L1, pattern = "^MT-")

QC指标存储在哪里？

> head(G48E2L1@meta.data,5)
                 orig.ident nCount_RNA nFeature_RNA percent.mt
AAACCTGAGGAATGGA    G48E2L1       8276         2922   6.911551
AAACCTGAGGCTCAGA    G48E2L1       2765         1253  10.886076
AAACCTGCAGCCTATA    G48E2L1       8529         2841   8.136945
AAACCTGGTCAGAATA    G48E2L1       9102         3003   5.207647
AAACCTGGTTAAAGAC    G48E2L1       3539         1721  13.139305

在官方的样例中，可视化QC指标，并进行cell过滤。

• 过滤unique feature counts 大于2500或少于200的细胞
• 过滤线粒体比例大于5%的细胞

先来可视化看一下：

VlnPlot()是Seurat中用于绘制单细胞数据的小提琴图函数(基因表达、指标、PC分数等)，小提琴图用于显示数据分布及其概率密度。

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

FeatureScatter通常用于可视化 feature-feature 关系，也可以用于计算对象的任何东西，i.e. 对象数据中的列，PC分数等。 个人理解：就是用点图看两个数据之间的相关性。

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

image.png

官方教程中在这里过滤掉 2500 > nFeature_RNA >200 和percent.mt < 5的数据：

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

但是我不想过滤，本文数据没有做过滤处理。哈哈哈！

数据标准化

从数据集中删除不需要的细胞后，下一步是数据标准化。默认情况下，我们使用全局缩放归一化方法“LogNormalize”，它将每个细胞的特征表达式测量值归一化为总的表达式，再乘以一个缩放因子(默认为10,000)，对结果对数化处理。标准化的数值存储在pbmc[["RNA"]]@data中。

G48E2L1 <- NormalizeData(G48E2L1, normalization.method = "LogNormalize", scale.factor = 10000)
G48E2L1[["RNA"]]@data

高度可变特征的识别(feature selection)

接下来，我们将计算数据集中显示高细胞间差异的特征子集(i。e，它们在一些细胞中高表达，在另一些细胞中低表达)。我们和其他人发现，在下游分析中关注这些基因有助于在单细胞数据集中突出生物信号。

这里详细描述了Seurat3中的过程，并通过直接建模单细胞数据中固有的均值-方差关系改进了先前的版本，并在FindVariableFeatures函数中实现。默认情况下，我们为每个数据集返回2,000个特性。这些将用于下游分析，如PCA。

Find variable features

识别“平均变异性图”上的异常点。

FindVariableFeatures(object, ...)

如何选择****selection.method:

vst:首先，用局部多项式回归(loess)拟合对数(方差)与对数(均值)的关系。然后使用观察到的平均值和期望的方差(由拟合线给出)对特征值进行标准化。然后，在裁剪到最大值之后，根据标准化的值计算特征方差(参见clip)。max参数)。

mean.var.plot (mvp):首先，使用一个函数计算每个特征的平均表达式(mean.function)和离散度(diffusion .function)。接下来，根据每个bin的平均表达式将特征划分为number .bin (默认 20)，并计算每个bin内的离散度z-score。这样做的目的是识别变量特征，同时控制可变性和平均表达之间的强烈关系。

dispersion(disp):选择分散值最高的基因

G48E2L1 <- FindVariableFeatures(G48E2L1, selection.method = "vst", nfeatures = 2000)

找出10个差异最大的基因：

top10 <- head(VariableFeatures(G48E2L1),10)
# 绘制带有和不带有标签的变量特性
plot1 <- VariableFeaturePlot(G48E2L1)
plot2 <- LabelPoints(plot = plot1, points = top10)
CombinePlots(plots = list(plot1, plot2))

数据的缩放（Scaling the data）

接下来，我们应用一个线性变换(“scaling”)，这是一个标准的预处理步骤，比PCA等降维技术更重要。

ScaleData函数功能:

• 改变每个基因的表达，使细胞间的平均表达为0
• 测量每个基因的表达，使细胞间的差异为1
  * 这一步给予下游分析同等的权重，这样高表达基因就不会占主导地位
• 其结果存储在G48E2L1[["RNA"]]@scale.data中

all.genes <- rownames(G48E2L1)
G48E2L1 <- ScaleData(G48E2L1, features = all.genes)

线性降维

接下来，我们对缩放的数据执行PCA。默认情况下，只使用前面确定的变量特性作为输入，但是如果您希望选择不同的子集，可以使用features参数来定义。

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

Seurat提供了几种有用的方法来可视化细胞和定义PCA的特性，包括VizDimReduction、DimPlot和DimHeatmap

检查和可视化PCA结果的几种不同的方法

> print(G48E2L1[["pca"]], dims=1:5, nfeatures = 5)
PC_ 1
Positive:  MALAT1, NEAT1, MT-ND6, XIST, AC004687.1
Negative:  RAN, H2AFZ, PRDX1, SLC25A5, PFN1
PC_ 2
Positive:  EEF1A1, BTF3, C1QBP, UNG, RPL5
Negative:  MKI67, CENPF, TOP2A, ASPM, UBE2C
PC_ 3
Positive:  MIR155HG, PTPRK, LIF, XIST, ANK3
Negative:  S100A4, LTB, S100A11, TMSB4X, FTL
PC_ 4
Positive:  CLDND1, BTG1, ARL6IP1, SELL, STAT1
Negative:  MIF, RPS18, RPS2, GZMB, SRM
PC_ 5
Positive:  GZMB, CCL5, CD8A, CD8B, CTSW
Negative:  HIST1H4C, NME4, CORO1B, STAT1, CDK1

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

DimPlot(G48E2L1, reduction = "pca")

特别是，DimHeatmap可以方便地探索数据集中主要的异构来源，并且在决定哪些PCs可以用于进一步的下游分析时非常有用。细胞和特征都是根据它们的PCA分数排序的。将cells设置为一个数字，可以绘制光谱两端的“极端”细胞，这极大地加快了绘制大型数据集的速度。虽然这显然是一个监督分析，但我们发现这是一个有价值的工具，用于探索相关的特征集。

DimHeatmap(G48E2L1, dims = 1, cells = 500, balanced = TRUE)

DimHeatmap(G48E2L1, dims = 1:15, cells = 500, balanced = TRUE)

image.png

确定数据的“维数”

为了克服scRNA-seq数据中单个特征中大量的技术噪声，Seurat根据他们的PCA评分将细胞分组，每个PC实质上代表一个“元特征”，它将跨相关特征集的信息组合在一起。因此，最主要的组件代表了数据集的健壮压缩。但是，我们应该选择包含多少个主成分? 10 个？ 20个? 100个？

在Macosko et al文章中，我们实现了一个重采样测试的灵感来自JackStraw程序。我们随机排列数据的一个子集(默认为1%)并重新运行PCA，构造一个特征得分的“null distribution”，然后重复这个过程。我们认为最“significant” 的PC是那些具有丰富的低p值特征的。

# Note：对于大型数据集，此过程可能需要很长时间，为了方便起见，请注释掉。
# 可以使用ElbowPlot()中实现的更近似的技术来减少计算时间
# 计算时间
G48E2L1 <- JackStraw(G48E2L1, num.replicate = 100)
G48E2L1 <- ScoreJackStraw(G48E2L1, dims = 1:20)

JackStrawPlot函数提供了一个可视化工具，用于用均匀分布(虚线)比较每个PC的p-values分布。“显著的”PCs将显示出一个低p值(虚线以上的实线)的强富集特性。在这种情况下，在最初的10-12个PCs之后，重要性似乎急剧下降。

JackStrawPlot(G48E2L1, dims = 1:15)

另一种启发式方法生成“Elbow plot”:根据各成分解释的方差百分比对主要成分进行排序(ElbowPlot 函数)。在这个例子中，我们可以观察到PC9-10周围的一个拐点（“elbow”），这表明大部分真实信号是在前10个pc中捕获的。

ElbowPlot(G48E2L1)

对用户来说，确定数据集的真实维数是一项挑战/不确定的工作。因此，我们建议考虑这三种方法。第一个是更有监督的，探索PCs以确定相关的异质性来源，并可与GSEA联合使用。第二个实现了一个基于随机空模型的统计测试，但是对于大型数据集来说非常耗时，并且可能不会返回一个明确的PC截止时间。第三种是一种常用的启发式算法，可以立即计算出来。在这个例子中，所有这三种方法都产生了相似的结果，但是我们可能有理由选择PC 7-12之间的任何一个作为截止时间。

我们在这里选择了10个，但鼓励用户考虑以下几点:

• 树突状细胞和NK细胞爱好者可能认识到，与PCs 12和PCs 13密切相关的基因定义了罕见的免疫亚群(即MZB1是浆细胞样dc的标记)。但是，这些组非常罕见，在没有预先知识的情况下，很难从这种大小的数据集的背景噪声中区分出来。

• 我们鼓励用户使用不同数量的pc(10台、15台甚至50台!)重复下游分析。正如你将观察到的，结果往往没有显著的不同。

• 我们建议用户在选择这个参数时，不要选得太高。例如，仅使用5个PCs执行下游分析会显著影响结果。

细胞聚类（Cluster the cells）

Seurat v3应用了一种基于图的集群方法，建立在(Macosko等人)的初始策略之上。重要的是，驱动聚类分析的距离度量(基于先前确定的PCs)保持不变。然而，我们将细胞距离矩阵划分成集群的方法已经得到了极大的改进。我们的方法受到最近手稿的很大启发，这些手稿将基于图的聚类方法应用于scRNA-seq数据 [SNN-Cliq, Xu and Su, Bioinformatics, 2015]和CyTOF数据 [PhenoGraph, Levine et al., Cell, 2015]。简单地说，这些方法将单元格嵌入到一个图结构中——例如k -最近邻(KNN)图，在具有相似特征表达模式的单元格之间绘制边缘，然后尝试将这个图划分为高度互连的准团或社区。

和表现型一样，我们首先在PCA空间中构造一个基于欧氏距离的KNN图，然后根据任意两个细胞在局部区域的共享重叠(Jaccard相似性)来细化它们之间的边权值。此步骤使用FindNeighbors函数执行，并将之前定义的数据集维度(前10个pc)作为输入。

为了对单元进行聚类，我们接下来应用模块化优化技术，如Louvain算法(default)或SLM [SLM, Blondel et al., Journal of Statistical Mechanics]，以迭代方式将单元分组在一起，目标是优化标准模块化函数。FindClusters函数实现这个过程，并包含一个分辨率参数，该参数设置下游集群的粒度，增加的值将导致更多的集群。我们发现，将该参数设置在0.4-1.2之间，对于3K左右的单细胞数据集通常会得到良好的结果。对于较大的数据集，最佳分辨率通常会增加。可以使用Idents函数找到集群。

>G48E2L1 <- FindNeighbors(G48E2L1, dims = 1:10)
>G48E2L1 <- FindClusters(G48E2L1, resolution = 0.5)
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 7038
Number of edges: 226547

Running Louvain algorithm...
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8435
Number of communities: 9
Elapsed time: 0 seconds

查看前5个细胞的cluster id

> head(Idents(G48E2L1),5)
AAACCTGAGGAATGGA AAACCTGAGGCTCAGA AAACCTGCAGCCTATA AAACCTGGTCAGAATA AAACCTGGTTAAAGAC
               6                1                0                0                4
Levels: 0 1 2 3 4 5 6 7 8

进行非线性降维（UMAP/tSNE）**

Seurat提供了几种非线性的降维技术，如tSNE和UMAP，以可视化和探索这些数据集。这些算法的目标是学习数据的底层流形，以便在低维空间中将相似的单元放在一起。上面所确定的基于图的集群中的单元应该在这些降维图上共同定位。作为UMAP和tSNE的输入，我们建议使用相同的PCs作为聚类分析的输入。

G48E2L1 <- RunUMAP(G48E2L1, dims = 1:10)
# 注意，您可以设置'label = TRUE'，或者使用LabelClusters函数来帮助标记各个clusters.

DimPlot(G48E2L1, reduction = "umap")

此时可以保存对象，这样就可以轻松地将其加载回来，而不必重新运行上面执行的计算密集型步骤，或者轻松地与协作者共享。

saveRDS(G48E2L1, file = "../output/G48E2L1_tutorial.rds")

找差异表达features（cluster biomarkers）

Seurat可以帮助您找到通过差异表达式定义集群的标记。默认情况下，它识别单个簇的阳性和阴性标记(在ident.1中指定)，与所有其他细胞相比较。Findallmarkers为所有集群自动化这个过程，但是您也可以测试集群组之间的相互关系，或者测试所有细胞。

min.pct参数要求至少在两组细胞中的任何一组中检测一个特性，以及thresh.test参数要求一个特性在两组之间有一定的差异(平均)。您可以将这两个值都设置为0，但是时间上有很大的增加——因为这将测试大量不太可能具有高度歧视性的特性。作为加速这些计算的另一个选项，max.cells.per.ident可以设置。这将对每个标识类进行采样，使其不具有比设置的细胞更多的细胞。虽然通常会有功率的损失，速度的增长可能是显著的，最高度差异表达的特征可能仍然会上升到顶部。

# find all markers of cluster 1
cluster1.markers <- FindMarkers(G48E2L1, ident.1 = 1, min.pct = 0.25)
head(cluster1.markers, n=5)
               p_val  avg_logFC pct.1 pct.2     p_val_adj
MT-CO1  0.000000e+00  0.8509450 0.995 0.993  0.000000e+00
ACTB    0.000000e+00 -0.8800403 1.000 0.999  0.000000e+00
EEF1A1  0.000000e+00 -0.9152317 0.974 0.992  0.000000e+00
B2M     0.000000e+00 -1.2646953 0.900 0.991  0.000000e+00
RPL28  1.044112e-307  0.6582244 0.993 0.962 3.501744e-303

找出区分cluster 5与cluster 0和cluster 3的所有标记

cluster5.markers <- FindMarkers(G48E2L1, 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
B2M    6.780927e-144 -0.6074824 0.998 0.998 2.274187e-139
ACTB   3.018783e-139 -0.4931734 1.000 1.000 1.012439e-134
PRDX1  1.243987e-134 -0.6805939 0.797 0.994 4.172084e-130
EEF1A1 9.412366e-126 -0.4098785 1.000 1.000 3.156719e-121
MT-CO1 2.831898e-125  0.5641633 0.992 0.990 9.497620e-121

找出每个cluster的标记与所有剩余的细胞相比较，只报告阳性细胞

G48E2L1.markers <- FindAllMarkers(G48E2L1, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
G48E2L1.markers %>% group_by(cluster) %>% top_n(n=2, wt = avg_logFC)

Seurat有几个关于微分表达式的测试，可以通过该测试设置。使用参数(详情请参阅我们的DE vignette)。例如，ROC测试返回任何单个标记(从0 - random到1 - perfect)的分类能力。

cluster1.markers <- FindMarkers(G48E2L1, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)

我们包括一些可视化标记表达的工具。VlnPlot(显示跨集群的表达式概率分布)和FeaturePlot(在tSNE或PCA图上可视化特性表达式)是我们最常用的可视化方法。我们还建议使用RidgePlot、CellScatter和DotPlot作为查看数据集的额外方法。

VlnPlot(G48E2L1, features = c("MS4A1", "CD79A"))

# you can plot raw counts as well
VlnPlot(G48E2L1, features = c("NKG7", "PF4"),slot = "counts", log = TRUE)

FeaturePlot(G48E2L1, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))

DoHeatmap为给定的细胞和特征生成一个表达式heatmap。在本例中，我们绘制每个集群的前20个标记(如果小于20，则绘制所有标记)。

top10 <- G48E2L1.markers %>% group_by(cluster) %>% top_n(n=10, wt = avg_logFC)
DoHeatmap(G48E2L1, features = top10$gene) + NoLegend()

image.png

为clusters分配细胞类型标识

幸运的是，在这个数据集的情况下，我们可以使用规范的标记，以方便地匹配无偏聚类到已知的细胞类型:

image.png

new.cluster.ids <- c("Naive CD4 T", "Memory CD4 T", "CD14+ Mono", "B", "CD8 T", "FCGR3A+ Mono", "NK", "DC","Platelet")
names(new.cluster.ids) <- levels(G48E2L1)
G48E2L1 <- RenameIdents(G48E2L1, new.cluster.ids)
DimPlot(G48E2L1, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()

image.png

saveRDS(G48E2L1, file = "../output/G48E2L1_final.rds")