单细胞入门(2)-经典案例分析
一、文件预处理
1.相关文件下载——单样本
(1)R包:
Seurat、dplyr、libr.。
Seurat介绍:
Analysis, visualization, and integration of Visium HD spatial datasets with Seurat • Seurathttps://satijalab.org/seurat/articles/pbmc3k_tutorial.html
(2)数据来源:
10X Genomics公司提供的外周血单核细胞(Peripheral Blood Mononuclear Cells,PBMC)数据集,有2700个单细胞在illuminate NextSeq500上测序,https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gzhttps://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz
(3)数据库
CellMarker:CellMarker2.0
PanglaoDB:
Single Cell Portal - Broad Institute:
Single Cell Expression Atlas - EMBL-EBI:Home < Single Cell Expression Atlas < EMBL-EBI
Mouse Cell Atlas:MCA | Mouse Cell Atlas
CancerSEA:CancerSEA - Database Commons
Mouse Cell Atlas:MCA | Mouse Cell Atlas
scRNASeqDB:scRNASeqDB
注意:由于软件更新等原因,部分名称可能有一小部分的不同,大家根据提示选用正确的数据。加油偶。
2.读取数据
(1)读取数据(设置seurat对象)
该函数从 10X 读取 cellranger 管道的输出,返回唯一的分子识别 (unique molecular identified,UMI) 计数矩阵。此矩阵中的值表示在每个单元格(列)中检测到的每个特征(即 gene;row)的分子数。
使用 count 矩阵创建一个对象。该对象用作容器,其中包含单个单元数据集的数据(如计数矩阵)和分析(如 PCA 或聚类结果)。
getwd() #获取工作路径
#BiocManager::install("Seurat")
library(Seurat)
library(dplyr)
library(patchwork) #加载包rm(list = ls())#删除环境变量
#导入数据,读取数据,是一些tsv、mtx文件
pbmc.data <- Read10X(data.dir = "./pbmc3k_filtered_gene_bc_matrices/filtered_gene_bc_matrices/hg19")
#创建seurat对象
pbmc <- CreateSeuratObject(counts = pbmc.data,project = "pbmc3k",min.cells = 3,min.features = 200)
project ="pbmc3k":pbmc3k相当于对样本进行标记了;
min.cells = 3:一个基因至少在3个细胞中表达才能被保留;min.features = 200:这个细胞表达多少个基因。根据样本特征设置,过滤低质量的细胞。
#读取稀疏矩阵,一个
matrix_data <- read.table("xxx",seq = "\t",head = T.row.names = 1)
dim(matrix_data)
seurat_obj <- CreateSeuratObject(counts = matrix_data)###读取RDS5文件
pbmc <- readRDS5("xxx")
saveRDS5 <- (pbmc, "pbmc.rds") #保存
(2)标准预处理工作流程
以下步骤包括 Seurat 中 scRNA-seq 数据的标准预处理工作流程。这些代表了基于 QC 指标的细胞选择和过滤、数据标准化和缩放以及高度可变特征的检测
QC 和选择用于进一步分析的细胞
Seurat 允许您轻松探索 QC 指标并根据任何用户定义的标准筛选细胞。常用的一些 QC 指标包括
- 在每个细胞中检测到的唯一基因的数量。低质量的细胞或空液滴通常具有很少的基因;细胞双联体或多联体可能表现出异常高的基因计数
- 在细胞内检测到的分子总数(与独特基因密切相关)
- 映射到线粒体基因组的 reads 百分比
- 低质量/垂死细胞通常表现出广泛的线粒体污染
- 我们使用该函数计算线粒体 QC 指标,该函数计算源自一组特征的计数百分比
PercentageFeatureSet()
- 我们使用以 开头的所有基因集作为一组线粒体基因
MT-
- 独特基因的数量和总分子的数量在
CreateSeuratObject()
。可以在对象元数据中找到
# 获取 RNA 测序数据的计数矩阵(Seurat >= 5.0.0)
counts_matrix <- GetAssayData(pbmc, assay = "RNA", layer = "counts")
#质控及数据可视化,添加线粒体基因的百分比
pbmc[['percemt.mt']] <- PercentageFeatureSet(pbmc,pattern = 'MT-')
head(pbmc@meta.data)
# orig.ident nCount_RNA nFeature_RNA percemt.mt
# AAACATACAACCAC-1 pbmc3k 2419 779 3.0177759
# AAACATTGAGCTAC-1 pbmc3k 4903 1352 3.7935958
# AAACATTGATCAGC-1 pbmc3k 3147 1129 0.8897363
# AAACCGTGCTTCCG-1 pbmc3k 2639 960 1.7430845
# AAACCGTGTATGCG-1 pbmc3k 980 521 1.2244898
# AAACGCACTGGTAC-1 pbmc3k 2163 781 1.6643551
(3)可视化
QC 指标,并使用这些指标来筛选单元格。
- 我们过滤具有超过 2,500 个或少于 200 个唯一特征计数的单元格
- 我们过滤线粒体计数为 >5% 的细胞
#数据可视化,可以用小提琴图表示
VlnPlot(pbmc,features = c('nCount_RNA','nFeature_RNA','percemt.mt'),ncol = 3) #小提琴图# 或者FeatureScatter 通常用于可视化两个特征之间的关系
plot1 <- FeatureScatter(pbmc,feature1 = 'nCount_RNA',feature2 = 'nFeature_RNA')
plot2 <- FeatureScatter(pbmc,feature1 = 'nCount_RNA',feature2 = 'percemt.mt')
plot1 + plot2#质控,去除低质量的数据
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percemt.mt < 5)
ncol = 3:每行放的图数量。
(4)数据规范化
从数据集中删除不需要的单元格后,下一步是规范化数据。默认情况下,我们采用全局缩放归一化方法“LogNormalize”,该方法通过总表达式对每个单元格的特征表达式测量值进行归一化,将其乘以比例因子(默认为 10,000),然后对结果进行对数转换。在 Seurat v5 中,标准化值存储在 pbmc[["RNA"]]$data
中。
#数据归一化
pbmc <- NormalizeData(pbmc,normalization.method = 'LogNormalize',scale.factor = 10000)
Normalize相当于给每个细胞标准化。
(5)寻找高变基因
我们计算在数据集中表现出高细胞间差异的特征子集(即,它们在某些细胞中高表达,而在其他细胞中低表达)。在下游分析中关注这些基因有助于突出单细胞数据集中的生物信号。
在 Seurat 中通过直接建模单细胞数据中固有的均值-方差关系来改进以前的版本,并在函数中实现。默认情况下,我们为每个数据集返回 2000 。这些将用于下游分析,如 PCA。FindVariableFeatures()
#寻找特征/高变基因
pbmc <- FindVariableFeatures(pbmc,selection.method = 'vst',nfeatures = 2000)# 指定使用的特征选择方法为 vst(方差稳定变换)。vst 是一种常用的方法,用于识别在不同样本或条件下具有较高方差的特征。选取2000个基因
# VariableFeatures(pbmc) -- 可变特征列表,高表达的10个基因
top10 <- head(VariableFeatures(pbmc),10)
# [1] "PPBP" "S100A9" "IGLL5" "LYZ" "GNLY" "FTL" "PF4" "FTH1" "FCER1A" "GNG11"
plot3 <- VariableFeaturePlot(pbmc)
plot4 <- LabelPoints(plot = plot1,points = top10,repel = TRUE,xnudge = 0,ynudge = 0)
plot3 + plot4
(6)缩放数据
接下来,我们应用线性变换(“缩放”),这是 PCA 等降维技术之前的标准预处理步骤。功能:ScaleData()
- 改变每个基因的表达,使细胞间的平均表达量为 0
- 缩放每个基因的表达,使细胞间的方差为 1。此步骤在下游分析中给予同等的权重,因此高表达基因不会占主导地位
- 其结果存储在
pbmc[["RNA"]]$scale.data
- 默认情况下,仅缩放可变特征。
- 可以指定参数来缩放其他功能
features
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
Scale相当于每一种基因在所有细胞中的标准化,为了防止后续的细胞分群当中由于某一个基因表达值过于离群,从而造成细胞分群的差异/误差。
二、细胞亚群分析
在获取高质量的 RNA-seq 数据矩阵后,我们需要对单细胞数据进行降维分析和细胞聚类。目前,降维方法有很多,例如:主成分分析(PCA, Principal Component Analysis)、t-SNE(t-Distributed Stochastic Neighbor Embedding)、UMAP(Uniform Manifold Approximation and Projection)、FIt-SNE(Fast Interpolation-based t-SNE)、其他新兴降维方法。
这些方法的核心思想是将 高维数据映射到低维(通常是二维)空间,使得相似的细胞能够聚集在一起,便于进行可视化和进一步分析。
区别:
1.PCA 标准化、线性降维
- 主成分分析(PCA, Principal Component Analysis)
- PC 是主成分(Principal Component)的缩写
#数据缩放 – 标准化
# 缩放数据
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc,features = all.genes)# 线性降维,初步降维
pbmc <-RunPCA(pbmc,features = VariableFeatures(object = pbmc))
print(pbmc[["pca"]],dims = 1:5,nfeatures = 5)
## 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
(2)可视化
# 可视化。Seurat提供可视化细胞和定义PCA,包括功能的几种有用的方法ViziDimReduction(),DimPlot()和DimHeatmap()
VizDimLoadings(pbmc,dims = 1:2,reduction = "pca")
DimPlot(pbmc,reduction = "pca")
DimHeatmap(pbmc,dims = 1:9,cells = 500,balanced = TRUE)
特别是,可以轻松探索数据集中异质性的主要来源,并且在尝试决定要包含哪些 PC 以进行进一步的下游分析时非常有用。像元和特征都根据其 PCA 分数进行排序。设置为数字会在频谱的两端绘制“极端”单元格,从而显著加快大型数据集的绘制速度。虽然显然是一种监督分析,但我们发现这是探索相关特征集的宝贵工具。DimHeatmap()
cells
DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
15个pc的表现:前边几个可以,可以分离的比较清楚的,随着pc数量的增加,后续展示的不是那么的泾渭分明,这些是要舍弃的,但是数据可以适当多选一些,不要选的太少
(3)确定数据集的“维数”
为了克服 scRNA-seq 数据任何单个特征中的广泛技术噪音,Seurat 根据细胞的 PCA 分数对细胞进行聚类,每台 PC 基本上代表一个“元特征”,该元特征结合了相关特征集中的信息。因此,顶部主成分表示数据集的稳健压缩。但是,我们应该选择包含多少个组件呢?10?20?100?
高阶的pca分析,选出适合进行下一步分析的数据:
JackStraw
函数用于评估主成分分析(PCA)中特征的显著性。通过重抽样,可以生成多个随机数据集,并计算每个特征在这些随机数据集中的统计量(如 p 值),并将结果储存在pbmc对象中。重抽样分析的结果可以用于以下方面:①特征选择:根据重抽样分析得到的 p 值,可以确定哪些特征在 PCA中具有显著的贡献。②模型评估:通过比较原始数据和重抽样数据的结果,可以评估PCA模型的稳定性和可靠性。
ScoreJackStraw 函数用于计算每个主成分的 JackStraw 得分。JackStraw 得分是一种用于评估主成分中特征的显著性的指标。得分结果可以用于以下方面:①特征选择:根据得分的大小,可以确定哪些主成分中的特征具有较高的显著性,从而可以选择这些特征进行进一步的分析或建模。②主成分解释:得分可以帮助解释每个主成分的生物学意义,通过查看得分较高的特征,可以了解主成分所代表的生物学过程或细胞状态。
#确定数据维度,高阶的pca分析,选出适合进行下一步分析的数据
pbmc<-JackStraw(pbmc,num.replicate = 100)
pbmc<-ScoreJackStraw(pbmc,dims = 1:20)# 可视化处理
JackStrawPlot(pbmc,dims = 1:20)
ElbowPlot(pbmc) #散点图
通过可视化的方法告诉我们应该选取哪些pc进行下游的分析,在虚线之下的维度可以不要,认为它变化不够显著,不够为我们后续的分群提供参考。
与上图一样,前面几个PC,他的下降趋势很快,到后边比较平缓,大概10个pc之后后边的变化就不太有了,后边可以决定在find neighbor的时候,选取前10个pc即可。
识别数据集的真实维度——对用户来说可能具有挑战性/不确定性。因此,我们建议考虑这三种方法。第一个是更受监督的,探索 PC 以确定异质性的相关来源,例如可以与 GSEA 结合使用。第二个实现了基于随机空模型的统计测试,但对于大型数据集很耗时,并且可能无法返回清晰的 PC 截止值。第三种是常用的启发式方法,可以即时计算。在这个例子中,所有三种方法都产生了相似的结果,但我们可能有理由选择 PC 7-12 之间的任何东西作为截止点。
我们在这里选择了 10 个,但鼓励用户考虑以下几点:
树突状细胞和 NK 爱好者可能会认识到与 PC 12 和 13 密切相关的基因定义了罕见的免疫亚群(即 MZB1 是浆细胞样 DC 的标志物)。然而,这些群体非常罕见,在没有先验知识的情况下,很难将这种规模的数据集与背景噪声区分开来。
鼓励用户使用不同数量的 PC(10、15 甚至 50 台!)重复下游分析。正如您将观察到的,结果通常不会有太大差异。
建议在选择此参数时将错误偏高。例如,仅使用 5 台 PC 执行下游分析会对结果产生显着的不利影响。
2.聚类分析
降维后,我们可以利用不同的聚类算法(如 Louvain 聚类、K-means 聚类)对细胞进行分类。例如,在 肿瘤组织 中,我们可以通过降维分析将免疫细胞、上皮细胞、基质细胞等不同类型的细胞区分开来。
Seurat 应用了基于图形的聚类方法,以 (Macosko et al) 的初始策略为基础。重要的是,驱动聚类分析的距离指标(基于以前确定的 PC)保持不变。然而,我们将蜂窝距离矩阵划分为聚类的方法已经有了显着改进。我们的方法深受最近手稿的启发,这些手稿将基于图的聚类方法应用于 scRNA-seq 数据 [SNN-Cliq, Xu 和 Su, Bioinformatics, 2015] 和 CyTOF 数据 [PhenoGraph, Levine et al., Cell, 2015]。简而言之,这些方法将单元格嵌入到图形结构中——例如 K 最近邻 (KNN) 图,在具有相似特征表达模式的单元格之间绘制边缘,然后尝试将该图划分为高度互连的“准小圈子”或“社区”。
与PhenoGraph一样,我们首先根据PCA空间中的欧几里得距离构建一个KNN图,并根据它们局部邻域中的共享重叠(Jaccard相似性)来细化任意两个单元之间的边权重。此步骤使用函数执行,并将先前定义的数据集维度(前 10 台 PC)作为输入。FindNeighbors()
为了对单元进行聚类,我们接下来应用模块化优化技术,如 Louvain algorithm(默认)或 SLM [SLM, Blondel et al., Journal of Statistical Mechanics],以迭代方式将单元分组在一起,以优化标准模块化函数。该函数实现此过程,并包含一个 resolution 参数,该参数设置下游聚类的“粒度”,值增加会导致集群数量增加。我们发现将此参数设置在 0.4-1.2 之间通常对于大约 3K 细胞的单细胞数据集会返回良好的结果。对于较大的数据集,最佳分辨率通常会提高。可以使用该函数找到集群FindClusters()
Idents()
# 聚类细胞
# 建立KNN图,并基于其局部领域中的共享重叠细化任意两个单元之间的边权重
pbmc<-FindNeighbors(pbmc,dims = 1:10)#对细胞进行聚类,应用模块化优化技术,Louvain算法(默认)或SLM
# 以迭代方式将细胞分组在一起,目标是优化标准模块化函数
pbmc<-FindClusters(pbmc,resolution = 0.5)
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
resolution越大,获得的细胞类群越多,根据数据反复修改。
运行非线性降维 (UMAP/tSNE)
Seurat 提供了几种非线性降维技术,例如 tSNE (t-Distributed Stochastic Neighbor Embedding)和 UMAP(Uniform Manifold Approximation and Projection),用于可视化和探索这些数据集。这些算法的目标是学习数据集中的底层结构,以便将相似的单元放在低维空间中。因此,在上面确定的基于图形的聚类中分组在一起的单元格应该在这些降维图上共定位。
虽然我们和其他人经常发现 tSNE 和 UMAP 等 2D 可视化技术是探索数据集的宝贵工具,但所有可视化技术都有局限性,无法完全代表底层数据的复杂性。特别是,这些方法旨在保持数据集中的局部距离(即确保具有非常相似基因表达谱的细胞共定位),但通常不会保留更多的全局关系。我们鼓励用户利用 UMAP 等技术进行可视化,但要避免仅根据可视化技术得出生物学结论。
# 运行非线性降维(UMAP/tSNE)
pbmc <- RunUMAP(pbmc,dims = 1:10)
DimPlot(pbmc,reduction = "umap")
pbmc <- RunTSNE(pbmc,dims = 1:10)
DimPlot(pbmc,reduction = 'tsne')saveRDS(pbmc, file = "../output/pbmc_tutorial.rds")
umap适合大型的数据,计算比较快;
tSNE通常情况下适合小型的数据比较好,计算比较慢,分离效果好;
3.发现差异表达的特征(聚类生物标志物)
Seurat 可以帮助您找到通过差异表达 (DE) 定义簇的标记物。默认情况下,与所有其他像元相比,它识别单个聚类(在 中指定)的正标记和负标记。 为所有集群自动执行此过程,但您也可以测试集群组彼此之间或所有单元之间的对比。ident.1
FindAllMarkers()
在 Seurat v5 中,我们使用 presto 软件包(如此处所述,可在此处安装)来显著提高 DE 分析的速度,尤其是对于大型数据集。对于不使用 presto 的用户,您可以查看此函数 () 的文档以探索 and 参数,可以增加这些参数以提高 DE 测试的速度。?FindMarkers
min.pct
logfc.threshold
# 寻找差异表达的特征(簇生物标志物)
# findmarkers为所有集群自动执行此过程,也可以测试集群组之间的对比,或针对所有单元格进行测试
# 默认情况下,ident.1与所有其他细胞相比,他识别单个簇的阳性和阴性标记。
# min.pct参数要求在两组细胞中的任何一组中以最小百分比检测到一个特征
# 而 thresh.test 参数要求一个特征在两组之间差异表达(平均)一定量
# 寻找cluster2的所有markers
cluster2.markers <- FindMarkers(pbmc,ident.1 = 2,min.pct = 0.25)
head(cluster2.markers,n=5)
# 寻找cluster5中与cluster0和cluster3n不同的所有markers
cluster5.markers <- FindMarkers(pbmc,ident.1 = 5,ident.2=c(0,3),min.pct = 0.25)
head(cluster5.markers,n=5)# 找出每个细胞簇的标记物,与所有剩余的细胞进行比较,只报告阳性细胞
pbmc.markers <- FindAllMarkers(pbmc,only.pos = TRUE,min.pct = 0.25,logfc.threshold =0.25 )
pbmc.markers %>%group_by(cluster) %>%slice_max(n=2,order_by = avg_log2FC) #选出高表达的一些基因,管道符%>%,
Findmarkers可以指定一类类群和二类类群,细胞来找到高表达的marker基因,ident.1/ident.2细胞编号在上一步的图中,=5代表可以找出class5中的一些mark基因。
FindAllMarkers:通常用它寻找所有的marker基因,帮助我们寻找出每一个细胞类型中对应的高表达基因
%>%
是 magrittr 包引入的管道操作符,后被 dplyr 等包广泛采用,其核心功能是:
将左侧对象作为右侧函数的第一个参数传递,实现链式编程。
Seurat 有几种差异表达测试,可以使用 test.use 参数进行设置,
cluster0.markers <- FindMarkers(pbmc, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)
我们提供了几种用于可视化标记表达的工具。 (显示跨集群的表达概率分布)和 (在 tSNE 或 PCA 图上可视化特征表达)是我们最常用的可视化。我们还建议探索 、 和 作为查看数据集的其他方法。VlnPlot()
FeaturePlot()
RidgePlot()
CellScatter()
DotPlot()
# 还有用于可视化标记表达的工具,Vlnplot(显示跨集群的表达概率分布)和FeaturePlot()(在tSNE或PCA图上可视化特征表达)是最常用的可视化,建议探索RidgePlot(),CellScatter(),和DotPlot()作为查看数据集的其他方法
VlnPlot(pbmc,features = c("MS4A1","CD79A"))
# 可以绘制行数据
VlnPlot(pbmc,features = c("NKG7","PF4"),slot = "counts",log = TRUE)FeaturePlot(pbmc,features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP","CD8A"))
#DoHeatmap()为给定的细胞和特征生成一个表达热图。在这种情况下,我们为每个集群绘制前 20 个标记(或所有标记,如果小于 20)。
pbmc.markers%>%group_by(cluster)%>%top_n(n=10,wt=avg_log2FC)->top10
DoHeatmap(pbmc,features = top10$gene)+NoLegend()top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10,wt = avg_log2FC)
DoHeatmap(pbmc, features = top10$gene)+ NoLegend()
VlnPlot(pbmc, features = top10$gene[1:20],pt.size=0)
当确定高变基因在哪些细胞中高表达之后,这些细胞类群对应的生物学名称就明了了
4.为集群分配单元类型标识
- 通过文献检索,查看细胞类群的来源(类群名称)
- 通过cellmarker网站
在这个数据集的情况下,我们可以使用规范标记轻松地将无偏聚类与已知细胞类型相匹配:
# 将细胞类型标识分配给集群
new.cluster.ids<-c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T", "FCGR3A+ Mono","NK", "DC", "Platelet")
names(new.cluster.ids)<-levels(pbmc)
pbmc<-RenameIdents(pbmc,new.cluster.ids)
DimPlot(pbmc,reduction = "umap",label = TRUE,pt.size = 0.5)+NoLegend()saveRDS(pbmc,file = "./pbmc3k_final.rds")
- RDS形式存储suerat对象
- 之后可以直接读取(readRDS)
参考资料:
七龙珠 |召唤一份单细胞数据库汇总
可能是最全的单细胞数据库汇总!张泽民团队开发的排名35!
第三讲 单样本分析_哔哩哔哩_bilibili
3. Raw data processing — Single-cell best practices
https://github.com/satijalab/seurat/
Analysis, visualization, and integration of Visium HD spatial datasets with Seurat • Seurathttps://github.com/satijalab/seurat/blob/e44cb2ce21aff3cbd94c2ddf24f9a1db6745e121/vignettes/pbmc3k_tutorial.Rmd