Single-cell Workshop 2021 - 02 - 聚类和细胞类型鉴定

note
scRNA-seq
R
Author

Shixiang Wang

Published

July 26, 2023

学习资料：GitHub 地址

数据下载

我们将使用包含大约10K个单细胞的PBMC数据集。这个数据集是由10X基因组公司公开提供的，可以从 https://www.dropbox.com/s/wn4mgwkkzqw2pox/SC3_v3_NextGem_DI_PBMC_10K_filtered_feature_bc_matrix.h5?dl=0 下载。

准备

这个文章的内容灵感来自多个Seurat 文档。它假定你已经熟悉了最初的质控步骤。我们将从由cellranger生成的数据 x 细胞过滤矩阵开始，这是大多数分析的常见起点。

我们将专注于单细胞RNA测序分析的两个具体关键任务：聚类和对已识别聚类的注释。

安装包

install.packages("Seurat")
install.packages("hdf5r")
install.packages("clustree")
install.packages("BiocManager")
BiocManager::install("SingleR")
BiocManager::install("celldex")

加载：


Attaching package: 'dplyr'
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union
Attaching SeuratObject
Loading required package: ggraph
Loading required package: ggplot2
Loading required package: SummarizedExperiment
Loading required package: MatrixGenerics
Loading required package: matrixStats

Attaching package: 'matrixStats'
The following object is masked from 'package:dplyr':

    count

Attaching package: 'MatrixGenerics'
The following objects are masked from 'package:matrixStats':

    colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
    colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
    colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
    colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
    colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
    colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
    colWeightedMeans, colWeightedMedians, colWeightedSds,
    colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
    rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
    rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
    rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
    rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
    rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
    rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
    rowWeightedSds, rowWeightedVars
Loading required package: GenomicRanges
Loading required package: stats4
Loading required package: BiocGenerics

Attaching package: 'BiocGenerics'
The following objects are masked from 'package:dplyr':

    combine, intersect, setdiff, union
The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs
The following objects are masked from 'package:base':

    anyDuplicated, aperm, append, as.data.frame, basename, cbind,
    colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
    get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
    match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
    Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
    table, tapply, union, unique, unsplit, which.max, which.min
Loading required package: S4Vectors

Attaching package: 'S4Vectors'
The following objects are masked from 'package:dplyr':

    first, rename
The following objects are masked from 'package:base':

    expand.grid, I, unname
Loading required package: IRanges

Attaching package: 'IRanges'
The following objects are masked from 'package:dplyr':

    collapse, desc, slice
Loading required package: GenomeInfoDb
Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.

Attaching package: 'Biobase'
The following object is masked from 'package:MatrixGenerics':

    rowMedians
The following objects are masked from 'package:matrixStats':

    anyMissing, rowMedians

Attaching package: 'SummarizedExperiment'
The following object is masked from 'package:SeuratObject':

    Assays
The following object is masked from 'package:Seurat':

    Assays

Attaching package: 'celldex'
The following objects are masked from 'package:SingleR':

    BlueprintEncodeData, DatabaseImmuneCellExpressionData,
    HumanPrimaryCellAtlasData, ImmGenData, MonacoImmuneData,
    MouseRNAseqData, NovershternHematopoieticData

读入数据集并创建一个 Seurat 对象

首先，我们要读取10X测序数据并将其转换为 Seurat 对象。Seurat 对象作为一个容器，最初只包含 UMI 计数矩阵。但我们会向其添加更多的分析内容（例如PCA、聚类结果）。

我们首先通过使用 Read10X_h5 函数读取计数矩阵，该函数从10X CellRanger hdf5文件中读取计数矩阵。层次数据格式（Hierarchical Data Format，HDF5或H5）提供了数据的更压缩的表示形式。Seurat 软件包包含多个读取函数，具体取决于文件格式。

在读取数据时，我们会应用基本的质量控制，将低质量的细胞丢弃。

pbmc.data <- Read10X_h5(filename="/Users/wsx/Library/CloudStorage/OneDrive-shanghaitech.edu.cn/Public/data/SC3_v3_NextGem_DI_PBMC_10K_filtered_feature_bc_matrix.h5")
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc10k", min.cells = 3, min.features = 200)
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 & percent.mt < 15)

接下来，我们计算线粒体RNA在总RNA含量中的百分比贡献。高线粒体含量可能表明正在经历凋亡的低质量细胞。我们根据线粒体含量和每个细胞特征数的分布再次进行质量控制。特征数过低可能表明空的液滴中存在环境RNA污染。特征数过高可能是由于多个细胞被困在同一个液滴中引起的。

数据处理

经过质量控制措施以选择继续分析的细胞后，接下来的步骤包括对数据进行归一化、识别高度可变特征以及进行尺度缩放。同时，还会进行主成分分析，因为许多后续分析步骤将在较低维度空间中进行计算。

pbmc <- NormalizeData(pbmc)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
Centering and scaling data matrix
pbmc <- RunPCA(pbmc, verbose = FALSE)

前面所有步骤的基本原理都与质量控制相关，在前面的教程中有更详细的介绍。

选择主成分的数量

首要的重要决策是在后续分析中保留多少个主成分（PCs）。保留的主成分越多，信号量也会增加，但同时也会增加噪音，并且计算需求也会增加。在做出这个选择之前，请检查与每个主成分相关的基因。存在与特定细胞类型相关的基因将表明该主成分具有信息量。而存在与所有不相关基因名称相关的情况则表明相反。在外周血单个核细胞（PBMC）的背景下，我们预期主要的细胞类型将包括T细胞、B细胞、NK细胞、单核细胞等等。

检查前 10 个 PC，你会看到很多熟悉的基因名称。

print(pbmc[["pca"]], dims = 1:10, nfeatures = 5)
PC_ 1 
Positive:  LTB, IL32, TRAC, CD3D, TRBC2 
Negative:  FCN1, FGL2, CST3, IFI30, TYMP 
PC_ 2 
Positive:  MS4A1, CD79A, BANK1, IGHM, NIBAN3 
Negative:  IL32, GZMM, CD3D, CD7, CD247 
PC_ 3 
Positive:  GZMB, CLIC3, NKG7, GNLY, KLRD1 
Negative:  CCR7, LEF1, TRABD2A, TCF7, LTB 
PC_ 4 
Positive:  CD79B, MS4A1, GNLY, CD79A, LINC00926 
Negative:  LILRA4, CLEC4C, SERPINF1, TPM2, SCT 
PC_ 5 
Positive:  CDKN1C, HES4, CTSL, TCF7L2, BATF3 
Negative:  S100A12, ITGAM, VCAN, CES1, MGST1 
PC_ 6 
Positive:  CDKN1C, NRGN, PADI4, CKB, FCGR3A 
Negative:  CRIP1, GBP5, ISG15, NIBAN1, MAF 
PC_ 7 
Positive:  CLU, ITGB1, NRGN, LIMS1, GP1BB 
Negative:  MALAT1, CCR7, TXK, LINC02446, KLRK1 
PC_ 8 
Positive:  CDKN1C, VIM, S100A4, CKB, ITGB1 
Negative:  MAP3K7CL, SERPING1, GMPR, LGALS2, NEXN 
PC_ 9 
Positive:  HERC5, IFIT1, RSAD2, IFIT3, MX1 
Negative:  CLEC10A, FCER1A, ENHO, CD1C, CACNA2D3 
PC_ 10 
Positive:  TNFRSF13B, IGHG1, IGHA1, SSPN, IGFBP7 
Negative:  TCL1A, CD8A, TRGC2, GZMK, CD8B

另一种可视化信息的方法是：

VizDimLoadings(pbmc, dims = 1:2, reduction = "pca",balanced=TRUE)

这是另一种提供图示表示的方法。细胞和特征根据主成分分析得分进行排序。设置一个细胞数量有助于计算效率，因为它会忽略那些信息较少的极端细胞。强制使其平衡可在正相关和负相关之间获得相等的表示。

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

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

Jackstraw 图用于估计每个主成分所捕获的结构的显著性。它会随机对数据的子集进行排列，以建立空值分布，并根据这个空值分布估计p值。

⚠️ 这可能需要一些时间来执行。如果时间不够，请跳过这部分，然后只专注于接下来的拐点图。

pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
JackStrawPlot(pbmc, dims = 1:20)
Warning: Removed 28000 rows containing missing values (`geom_point()`).

拐点图是经典的计算机科学方法，用于检查通过逐个添加主成分来表示数据中累积变异性。

ElbowPlot(pbmc,ndims=50)

根据这些图，我们选择使用20个主成分。但也可以为15到30之间的任何数量提出理由。

我们使用UMAP来可视化数据集。我们提前计算UMAP，以便在需要时立即使用。

pbmc <- RunUMAP(pbmc, dims = 1:20, verbose = FALSE)
Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session

聚类

Seurat默认的聚类算法采用基于图的聚类方法。它受到以前的发表文章的启发，这些文章在开发和改进这种方法方面做出了贡献。特别是，《Phenograph》论文强烈推荐给那些希望更好地理解这种方法的人。简要来说，聚类算法首先在PCA空间中计算每个细胞的K个最近邻。然后，根据共享邻居计算细胞之间的Jaccard相似性。使用Louvain算法对信息进行聚合，以便将细胞迭代地分组在一起。

pbmc <- FindNeighbors(pbmc, dims = 1:20)
Computing nearest neighbor graph
Computing SNN

第二个重要的决策是选择聚类的分辨率。较小的数值会生成混合的聚类，而较大的数值会产生过多的聚类，这些聚类可能不太具有意义。

📝 对于你自己的研究，花一些时间来优化分辨率的选择是很有帮助的。尝试使用多个候选的分辨率选项来进行所有后续分析可能会有帮助。

在这里，我们尝试三个不同的分辨率值，并使用clustree来研究其影响。在计算完聚类后，我们还保存了pbmc对象。

pbmc <- FindClusters(object = pbmc,  resolution = c(0.5, 1, 1.5),  dims.use = 1:10,  save.SNN = TRUE)
Warning: The following arguments are not used: dims.use, save.SNN
Suggested parameter: dims instead of dims.use
Warning: The following arguments are not used: dims.use, save.SNN
Suggested parameter: dims instead of dims.use
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 9729
Number of edges: 356102

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9145
Number of communities: 16
Elapsed time: 1 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 9729
Number of edges: 356102

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.8655
Number of communities: 20
Elapsed time: 1 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 9729
Number of edges: 356102

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.8279
Number of communities: 22
Elapsed time: 1 seconds
saveRDS(pbmc, file = "/Users/wsx/Library/CloudStorage/OneDrive-shanghaitech.edu.cn/Public/data/pbmc_tutorial.rds")
clustree(pbmc)

这三行表示每个分辨率值下的细胞分配情况。每行中的节点大小表示该聚类中的细胞数量。行之间的箭头显示随着分辨率的增加，聚类的分配情况如何变化。稳定的聚类可能会更改名称，但是细胞在不同的分辨率下仍然会聚集在一起。一些聚类可能会分裂成两个（或更多）子聚类。这暗示着增加分辨率。但是，如果你看到聚类之间来回跳动的情况很多，那就表示稳定性较差。在你自己的数据中，最好进行过度聚类，然后检查基因，然后手动选择要合并的聚类。

我们选择分辨率=0.5进行进一步检查：

Idents(pbmc) <- pbmc$RNA_snn_res.0.5
DimPlot(pbmc, reduction = "umap", label=TRUE)

聚类富集标记

接下来，对于每个聚类，我们想要检查在该聚类中相比其他细胞高度表达的特征（也称为标记或基因）。我们使用’roc’方法来估计每个标记的分类能力（0表示随机，1表示完美）。我们可以逐个聚类地进行这个分析，也可以同时对所有聚类进行分析。

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

⚠️ 对所有聚类执行此操作非常耗时。如果你想休息一下。在这样做之前启动它。计算时间可能长达30分钟。如果你在教程时间方面运行晚了，现在跳过这个和下面的热图。

cluster.all.markers0.5 <- FindAllMarkers(pbmc, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE, min.pct = 0.25)
Calculating cluster 0
Calculating cluster 1
Calculating cluster 2
Calculating cluster 3
Calculating cluster 4
Calculating cluster 5
Calculating cluster 6
Calculating cluster 7
Calculating cluster 8
Calculating cluster 9
Calculating cluster 10
Calculating cluster 11
Calculating cluster 12
Calculating cluster 13
Calculating cluster 14
Calculating cluster 15

热图为共享高度表达标记的集群提供了很好的图形表示。注意到umap中这些集群之间的关系了吗?

cluster.all.markers0.5  %>%
    group_by(cluster) %>%
    top_n(n = 5, wt = avg_log2FC) -> top5

DoHeatmap(pbmc, features = top5$gene) + NoLegend()

接下来，我们选择一些标记，这些标记往往适用于PBMC数据集。其中许多来自cluster.all.markers0.5。

markers.to.plot <- c("CD3D", "HSPH1", "SELL", "CD14", "LYZ", "GIMAP5", "CACYBP", "GNLY", "NKG7", "CCL5", "CD8A", "MS4A1", "CD79A", "FCGR3A", "MS4A7", "S100A9", "HLA-DQA1","GPR183", "PPBP", "GNG11", "TSPAN13", "IL3RA", "FCER1A", "CST3", "S100A12")

DotPlot(pbmc, features = markers.to.plot, cols = c("blue", "red"), dot.scale = 8) +RotatedAxis()

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

VlnPlot(pbmc, features = c("NKG7", "GNLY"))

VlnPlot(pbmc, features = c("FCGR3A", "MS4A7"))

VlnPlot(pbmc, features = c("PPBP"))

VlnPlot(pbmc, features = c("FCER1A", "CST3"))   

VlnPlot(pbmc, features = c("CD8A", "CD8B", "CD3D"))

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

基于参考的聚类注释

基于参考数据的聚类注释不在Seurat软件包中进行。为此，我们将使用SingleR，并首先将Seurat数据的信息转换为SingleR的格式。

sce <- GetAssayData(object = pbmc, assay = "RNA", slot = "data")

我们将使用一篇经典论文中的数据，该论文使用分类细胞对不同的免疫细胞类型进行了比较。基于参考数据的注释允许与高度精选的数据集进行直接比较。然而，缺点是推断的结果取决于参考数据集中定义的内容。

refMonaco <- MonacoImmuneData()
snapshotDate(): 2022-10-31
see ?celldex and browseVignettes('celldex') for documentation
loading from cache
see ?celldex and browseVignettes('celldex') for documentation
loading from cache

数据格式在两个层次上显示信息:主要细胞类型和更精细的分辨率信息：

prediction_Monaco_main <- SingleR(test=sce, ref=refMonaco, clusters=Idents(pbmc), labels=refMonaco$label.main)
prediction_Monaco_fine <- SingleR(test=sce, ref=refMonaco, clusters=Idents(pbmc), labels=refMonaco$label.fine)

predicted_Monaco <- data.frame(cluster=sort(unique(Idents(pbmc))), Monaco_main= prediction_Monaco_main$labels, Monaco_fine= prediction_Monaco_fine$labels)
predicted_Monaco
   cluster     Monaco_main                  Monaco_fine
1        0       Monocytes          Classical monocytes
2        1    CD4+ T cells            Naive CD4 T cells
3        2       Monocytes          Classical monocytes
4        3    CD8+ T cells            Naive CD8 T cells
5        4    CD4+ T cells               Th1/Th17 cells
6        5         T cells           Non-Vd2 gd T cells
7        6         B cells  Non-switched memory B cells
8        7         B cells                Naive B cells
9        8        NK cells         Natural killer cells
10       9       Monocytes       Intermediate monocytes
11      10         T cells                   MAIT cells
12      11 Dendritic cells      Myeloid dendritic cells
13      12 Dendritic cells Plasmacytoid dendritic cells
14      13       Monocytes          Classical monocytes
15      14     Progenitors             Progenitor cells
16      15         B cells                 Plasmablasts