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

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

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

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

获取数据：（.tar.gz）

读取数据

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

代码语言：javascript代码运行次数：0运行复制

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 = "pbmck",  = , min.features = 200)
pbmc

输出结果：

代码语言：javascript代码运行次数：0运行复制

## An object of class Seurat 
## 1714 features across 2700 samples within 1 assay 
## Active assay: RA (1714 features, 0 variable features)
## 1 layer present: counts

数据质量控制和筛选

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

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

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

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

工具和方法：

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

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

代码语言：javascript代码运行次数：0运行复制

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

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

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

结果图：

代码语言：javascript代码运行次数：0运行复制

# FeatureScatter is typically used to visualize feature-feature relatihips, but can be used
# for anything calculated by the object,  columns in object metadata, PC scores etc.plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RA", feature2 = "")
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RA", feature2 = "nFeature_RA")
plot1 + plot2

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

代码语言：javascript代码运行次数：0运行复制

pbmc <- ormalizeData(pbmc,  ="Logormalize", scale.factor = 10000)
#pbmc <- ormalizeData(pbmc)

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

代码语言：javascript代码运行次数：0运行复制

pbmc <- FindVariableFeatures(pbmc,  = "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[["RA"]]$scale.data。默认情况下，仅缩放可变特征。你可以指定参数来缩放其他功能features。

代码语言：javascript代码运行次数：0运行复制

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

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

代码语言：javascript代码运行次数：0运行复制

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

结果：

代码语言：javascript代码运行次数：0运行复制

# PC_ 1 
# Positive: CST, TYROBP, LST1, AIF1, FTL 
# egative: MALAT1, LTB, IL2, IL7R, CD2 
# PC_ 2 
# Positive: CD79A, MS4A1, TCL1A, HLA-DQA1, HLA-DQB1 
# egative: KG7, PRF1, CST7, GZMB, GZMA 
# PC_  
# Positive: HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPB1 
# egative: PPBP, PF4, SDPR, SPARC, GG11 
# PC_ 4 
# Positive: HLA-DQA1, CD79B, CD79A, MS4A1, HLA-DQB1 
# egative: VIM, IL7R, S100A6, IL2, S100A8 
# PC_ 5 
# Positive: GZMB, KG7, S100A8, FGFBP2, GLY 
# egative: LTB, IL7R, CKB, VIM, MS4A7

代码语言：javascript代码运行次数：0运行复制

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

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

代码语言：javascript代码运行次数：0运行复制

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

代码语言：javascript代码运行次数：0运行复制

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

聚类

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

代码语言：javascript代码运行次数：0运行复制

pbmc <- Findeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)
# Look at cluster IDs of the first 5 cells
head(Idents(pbmc), 5)

代码语言：javascript代码运行次数：0运行复制

## AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1 
## 2  2 1 
## AAACCGTGTATGCG-1 
## 6 
## Levels: 0 1 2  4 5 6 7 8

非线性降维

非线性降维方法：

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

代码语言：javascript代码运行次数：0运行复制

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：设置基因在体中的最小表达比例。
  • ：设置最小对数差异阈值。
  • 集成presto软件包显著提升 DE 分析速度，尤其适用于大型数据集。
  • 若不使用presto，可调整以下参数以加快分析：

代码语言：javascript代码运行次数：0运行复制

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

代码语言：javascript代码运行次数：0运行复制

FeaturePlot(pbmc, features = c("MS4A1", "GLY", "CDE", "CD14", "FCER1A", "FCGRA", "LYZ", "PPBP", "CD8A"))

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

代码语言：javascript代码运行次数：0运行复制

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

总结： 在本教程中，我们以外周血单核细胞（PBMC）数据集为例，完整展示了单细胞转录组分析的基本流程。从数据预处理、特征选择到降维聚类，再到可视化与差异基因分析，通过这些分析，我们能够更深入地理解单细胞层面的生物学特性，挖掘细胞间的异质性。

本文参与 腾讯云自媒体同步曝光计划，分享自。原始发表：2025-01-09，如有侵权请联系 cloudcommunity@tencent 删除基础可视化数据工具函数