Seurat 3.0 实例教程

2019-11-19 本文已影响0人 11的雾

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)

参数	描述	默认值
counts	未标准化的数据，如原始计数或TPMs
project	设置Seurat对象的项目名称	"SeuratProject"
assay	与初始输入数据对应的分析名称。	"RNA"
min.cells	包含至少在这么多细胞中检测到的features。将子集计数矩阵以及。要重新引入被排除的特性，请创建一个具有较低截止值的新对象。	0
min.features	包含至少检测到这些features的细胞。	0
names.field	对于每个cell的初始标识类，请从cell的名称中选择此字段。例如，如果您的cell在输入矩阵中被命名为BARCODE_CLUSTER_CELLTYPE，则设置名称。字段设置为3以将初始标识设置为CELLTYPE。	1
names.delim	对于每个cell的初始标识类，请从cell的列名中选择此分隔符。例如，如果您的cell命名为bar - cluster - celltype，则将此设置为“-”，以便将cell名称分离到其组成部分中，以选择相关字段。	"_"
meta.data	要添加到Seurat对象的其他细胞水平（cell-level）数据。应该是一个数据框，其中行是cell name，列是附加的元数据字段。	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")

Seurat 3.0 实例教程

标准的预处理流程：

数据标准化

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

数据的缩放（Scaling the data）

线性降维

确定数据的“维数”

细胞聚类（Cluster the cells）

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

找差异表达features（cluster biomarkers）

为clusters分配细胞类型标识

猜你喜欢

热点阅读