跟着Seurat 官网学单细胞转录组分析

天意生信云

发布于 2025-01-22 08:05:44

5620

在大家进行了一段时间的R语言与Linux学习后，我们开启单细胞测序数据的学习。接下来的教程中，我们将以Seurat 分析框架 为基础，从数据预处理、聚类分析到可视化的完整流程，深入讲解如何从原始数据中提取有意义的生物学信息。

我们将分析10X Genomics免费提供的外周血单核细胞（PBMC）数据集，这一数据集包含了 2,700 个单细胞，使用 Illumina NextSeq 500 进行测序，数据质量可靠，非常适合初学者学习和练习。

获取数据：（https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz）

读取数据

读取数据函数从10X 读取 cellranger 管道的输出，返回唯一的分子识别（UMI）计数矩阵。此矩阵中的值表示在每个单元格（列）中检测到的每个特征（即 gene;row）的分子数。请注意，最新版本的 cellranger 现在也使用 h5 文件格式输出，可以使用 Seurat 中的函数读取该格式。


library(dplyr)
library(Seurat)
library(patchwork)
# Load the PBMC dataset
pbmc.data <- Read10X(data.dir = "/brahms/mollag/practice/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)
## 1 layer present: counts

数据质量控制和筛选

细胞质量控制 (QC) 指标：

1、唯一基因数量：检测到的基因数过低可能表示低质量细胞或空液滴；过高可能表示细胞双联体/多联体。

2、分子总数：与唯一基因数量相关，用于评估细胞质量。

3、线粒体基因比例：线粒体 reads 比例高可能表示低质量或垂死细胞。

工具和方法：

1、使用 PercentageFeatureSet() 函数计算线粒体基因的 reads 百分比。

2、将以“MT-”开头的基因视为线粒体基因集。


# The [[ operator can add columns to object metadata. This is a great place to stash QC stats
pbmc["percent.mt"] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
# Visualize QC metrics as a violin plot
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

我们过滤具有超过 2,500 个或少于 200 个唯一特征计数的单元格

我们过滤线粒体计数为 >5% 的细胞

结果图：


# 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.plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2

从数据集中删除不需要的单元格后，下一步是规范化数据。默认情况下，我们采用全局缩放归一化方法“LogNormalize”，该方法通过总表达式对每个单元格的特征表达式测量值进行归一化，将其乘以比例因子（默认为 10,000），然后对结果进行对数转换。


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

接下来，我们计算在数据集中表现出高细胞间差异的特征子集（即，它们在某些细胞中高表达，而在其他细胞中低表达）。

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)
plot2

线性降维

拿到特征子集之后，需要进行降维之前的标准预处理步骤，函数：ScaleData()。

改变每个基因的表达，使细胞间的平均表达量为 0；缩放每个基因的表达，使细胞间的方差为 1此步骤在下游分析中给予同等的权重，因此高表达基因不会占主导地位，其结果存储在pbmc[["RNA"]]$scale.data。默认情况下，仅缩放可变特征。你可以指定参数来缩放其他功能features。

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

接下来，我们对缩放数据执行 PCA。默认情况下，只有先前确定的变量特征用作输入，但如果你想选择不同的子集，可以使用 argument 进行定义。featuresScaleData对于第一个主成分，Seurat 输出具有最多正负载和负负载的基因列表，代表数据集中单个细胞之间表现出相关性（或反相关性）的基因模块。

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")

结果：

# 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

DimPlot(pbmc,reduction = "pca") + NoLegend()

可以轻松探索数据集中异质性的主要来源，并且在尝试决定要包含哪些 PC 以进行进一步的下游分析时非常有用。像元和特征都根据其 PCA 分数进行排序。设置为数字会在频谱的两端绘制“极端”单元格，从而显著加快大型数据集的绘制速度。显然是一种监督分析，但我们发现这是探索相关特征集的宝贵工具。DimHeatmap()cells.

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

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

聚类

Seurat 应用了基于图形的聚类方法，以（Macosko et al）的初始策略为基础。驱动聚类分析的距离指标（基于以前确定的 PC）保持不变。


pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)
# Look at cluster IDs of the first 5 cells
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

非线性降维

非线性降维方法：

Seurat 提供tSNE和UMAP等方法，用于将高维数据降到低维空间，便于可视化和数据探索。
降维的目标是将相似的细胞在低维空间中定位在一起，与基于图的聚类结果一致。

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")

聚类生物标志物

差异表达（DE）分析：
- Seurat 可用于识别定义聚类的正负标记基因。
- 默认情况下，比较单个聚类与所有其他聚类的差异，也可以比较聚类组之间或与所有单元的对比。
相关功能：
- 使用FindAllMarkers()自动执行所有聚类的标记基因发现。
- 使用FindMarkers()测试特定的聚类组之间的差异。
性能优化（Seurat v5）：
- min.pct：设置基因在群体中的最小表达比例。
- logfc.threshold：设置最小对数差异阈值。
- 集成presto软件包显著提升 DE 分析速度，尤其适用于大型数据集。
- 若不使用presto，可调整以下参数以加快分析：


cluster0.markers <- FindMarkers(pbmc, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)
VlnPlot(pbmc, features = c("MS4A1", "CD79A"))

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

DoHeatmap()为给定的单元格和特征生成表达式热图。在本例中，我们将为每个聚类绘制前 20 个标记（如果少于 20 个，则绘制所有标记）。


pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE)
pbmc.markers %>%
  group_by(cluster) %>%
  dplyr::filter(avg_log2FC > 1)
pbmc.markers %>%
  group_by(cluster) %>%
  dplyr::filter(avg_log2FC > 1) %>%
  slice_head(n = 10) %>%
  ungroup() -> top10
DoHeatmap(pbmc, features = top10$gene) + NoLegend()