2023-04-26

第1章 seurat包简介

第1.1节什么是seurat

seurat是R语言下的单细胞RNA测序(scRNA-seq)数据分析包。它能够处理和分析来自不同平台(包括10x Genomics、DroNc-seq、Smart-seq等)的scRNA-seq数据。seurat在数据整理、质量控制、可视化、粒度聚类、差异表达检测和空间定位等方面都有出色的表现。它是单细胞分析领域使用最广泛的R包,能够有效支持研究人员快速分析scRNA-seq数据,挖掘数据背后的生物学意义。

第1.2节 seurat包能做什么

seurat包能够轻松完成以下 scRNA-seq 数据分析步骤:

1)数据导入与预处理:从多个文件格式(包括mtx、csv、rds等)导入数据,对数据进行过滤、标准化、归一化和质量控制等预处理。

2)可视化:丰富的绘图功能,包括散点图、小提琴图、热图等,可以直观展示数据质量与特征。

3)降维与聚类:使用PCA、TSNE、UMAP等对数据进行降维,并采用FindClusters等函数进行聚类。

4)差异表达基因检测:使用FindMarkers和其他函数检测不同细胞类型的差异表达基因,得到标记基因。

5)空间定位与区域检测:可以结合生物学知识,在空间中定位和检测已知的细胞区域。

6)细胞类型注释与预测:建立分类模型,对新样本进行细胞类型预测和注释。

7)高级分析:包括DNARanges、TransferData、CellCycleScoring等,实现更加深入的分析。

8)与其他包无缝集成:seurat可以与scran、scater、scTSS、Signac等包无缝集成,实现更加强大和全面的分析。

除此之外,seurat还在不断更新,推出更多实用功能,以满足研究人员的分析需求。它已经成为单细胞分析领域的基石型包,值得每个人精通与运用。

第1.3节如何安装seurat包

要安装seurat包,需要首先安装R语言。然后在R控制台中输入以下命令:

> install.packages("seurat")

这将使用CRAN仓库安装seurat的最新版本。

如果需要安装开发版本,可以从GitHub安装:

> devtools::install_github("satijalab/seurat")

seurat还依赖一些其他包,安装过程中会自动安装依赖包。如果安装过程中报错,很有可能是缺少依赖包而导致的。可以在未来进行补充安装:

> install.packages(c("Matrix", "graphics", "ggplot2", "dplyr", "reshape2", "leiden", "Rcpp"))

seurat支持macOS、Windows、Linux等主流操作系统。但由于Windows环境下编译依赖包的难度较大,安装过程中更容易报错。如果在Windows下无法解决安装错误,可以考虑安装RStudio Desktop或使用Linux子系统。

安装完成后,就可以在R中加载seurat包,开始单细胞数据的探索之旅!

第1.4节 seurat包版本更新日志

seurat是活跃开发的R包,会不断推出新的版本。每个版本都会修复已知bug,优化现有功能,并添加新特性。要访问最新功能,需要及时更新seurat包。

可以在seurat GitHub仓库的Release页面查看详细的版本更新日志。以3.1.4版本为例:

Fixed a bug that prevented DimPlot from correctly ordering cells when reduction = “umap”
Fixed a bug in the default print method for SeuratAssays objects
Added getSparsity() function to calculate sparsity metrics for a Seurat assay
Added name checking in SelectFeatrues and related functions
……

通读更新日志可以了解每个版本修复的bug、修改的函数和新增的功能特性。这对学习seurat包和利用新功能非常有帮助。

当seurat推出新版本后,可以使用R命令install.packages(“seurat”)或devtools::install_github(“satijalab/seurat”)安装更新。

安装新版本将覆盖原有版本,无需额外操作。由于新版本的推出速度较快,如果长时间未更新,升级后运行的代码可能会报错,需要根据更新日志进行相应修改。

第2章初识seurat对象，读入与访问数据

第2.1节 seurat对象的基本元素

一个Seurat对象包含以下基本要素:

1. assays:存储表达矩阵和其他生物学数据。例如,默认的exprs矩阵存储基因表达数值。

2. pca:保存主成分PCA结果,用于降维和可视化。

3. pcs. lui:保存LUI scores,实现局部可集聚降维。

4. tmap: TSNE或UMAP结果,用于二维空间的可视化。

5. graph: 记录细胞之间的相似性或连接,用于聚类。

6. meta:包含样本元数据,如 pretty_names、nUMI、percent.mt等。

7. reducers:记录其他的选择性降维方式,如LLE、Isomap等。

8. markers:存储callingcard分析的结果。

9. spatial:用于定位样品中的空间检测工具结果。

10. clusters:存储FindClusters()等函数生成的聚类标签。

11. scale.data:标准化或归一化表达矩阵的缩放因子,用于逆转缩放。

12. data.info:记录原始数据集的信息,如矩阵行列数、 Kuramoto主要聚类数等。

13. reduction:记录降维技术的名称,如”pca”、“tsne”、“umap”等。

14. genes:存储选择进行下游分析的基因或特征的名称。

15. min.genes:记录 FindClusters()等函数聚类时确定的最低基因表达检测阈值。

16. tmp.*:存储中间结果,以”tmp.”开头。例如 FindAllMarkers()的中间输出tmp.markers。

17. features:选定特征的详细信息,如ChooseAssayFeature()输出。

18. scale.factors:存储执行ScaleData()的缩放系数。

19. embeddings:存储执行Embeddings()后的坐标。

20. cell.names:存储细胞名称,默认为空。

21. project.name:存储项目名称,默认为”seurat”。

22. misc:存储其他结果,如feature.loadings()的输出。

23. @ backingfile:存储Seurat对象构建过程中使用的R数据文件。

24. ident:存储FindClusters等函数生成的聚类标签名称,如”Cluster 1”、“Cluster 2”等。

......

除此之外,Seurat对象动态添加新功能产生的结果,这使得对象变得非常丰富和强大。熟练理解和运用各个要素,是掌握seurat的关键所在。

第2.2节从CSV、RDS等文件读入数据

seurat支持从多种格式读入单细胞数据,主要包括:

1.最常用的方式: 从包含matrix.mtx，barcodes.tsv，genes.tsv三个文件的文件夹读取数据

> data <- Read10X(data.dir = "hg19/")
> sobj <- CreateSeuratObject(data,project = "pbmc3k", min.cells = 3, min.features = 200)

> #批量读取文件夹可采用以下方式
> #1.获取当前目录下以1或2结尾的所有文件夹名
> folders = list.files('./', pattern = '[12]$')
> #2.批量读取data
> dataList = lapply(folders, function(folder) {
+   CreateSeuratObject(counts = Read10X(folder), project = folder)
+ })
> #3.合并所有data，将多个Seurat对象合并成一个大的Seurat对象
> data.big <- merge(dataList[[1]], 
+                  y = c(dataList[[2]],dataList[[3]],dataList[[4]],dataList[[5]],
+                        dataList[[6]],dataList[[7]],dataList[[8]]), 
+                  add.cell.ids = folders, 
+                  project = "DCM")
> #4.可以使用 table() 函数统计 data.big 中 orig.ident 列的不同取值的个数
> table(data.big$orig.ident)

2. CSV:逗号分隔值文件,使用read.csv()读入。

data <- read.csv("data.csv")

3. RDS:R数据序列化文件,使用readRDS()读入。

data <- readRDS("data.rds")

4. mtx:矩阵市场交换文件,使用Read10X()读入10x Genomics数据。

data <- Read10X("path/to/data/")

5. loom:Loom格式文件,使用ReadLoom()读入。

data <- ReadLoom("data.loom")

6. hdf5:HDF5格式文件,使用Read10X_h5()读入。

data <- Read10X_h5("path/to/data.h5")

7. anndata:用anndata包读入后转为Seurat,使用ImportData()。

data <- anndata::read_h5ad("data.h5ad")
sobj <- Seurat::ImportData(data)

除此之外,seurat还支持从其他数据集格式进行转化和导入。多种格式支持使得seurat可以处理各种常见的单细胞数据。

导入数据后,得到的对象(如data)中包含表达矩阵以及样本元数据。我们可以查看数据维度和头几行,判断是否导入成功。下一步就可以创建Seurat对象开始分析了。

第2.3节创建seurat对象

从数据文件导入数据后,需要将其封装在Seurat对象中,才能使用seurat提供的各种分析功能。在seurat中,我们使用Seurat对象存储和表示单细胞数据。要创建一个Seurat对象,可以使用CreateSeuratObject()函数:

> library(tidyverse)
> library(Seurat)
> library(patchwork)
> data <- Read10X(data.dir = "hg19/")  #读入数据
> sobj <- CreateSeuratObject(data)
## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')

该函数需要输入 scRNA-seq数据,支持矩阵(matrix)、CSV、RDS等多种格式。创建成功后,sobj将是一个Seurat对象,存储原始数据并具有分析功能。

Seurat对象保存了样本的特征基因/代谢物表达矩阵,以及用于后续分析的诸多中间结果。这使得分析过程中的任何一步都可以回顾、审视或重新运行。事实上,Seurat对象可视为一个完整的单细胞学习环境。

该函数将导入的数据(如data)转化为Seurat对象sobj。要成功创建对象,data必须包含:

1. 数据矩阵: 表达值矩阵,行为基因/特征,列表为细胞。

2. 细胞名称: 行名为细胞的唯一ID。

3. sample.id:(可选) 记录样本唯一ID。

4. genes:(可选) 基因或特征名称。

除了必要的数据元素外,也可以提供单细胞测序平台或提取方法等元数据。CreateSeuratObject()会将所有信息整合,创建一个功能完备的Seurat对象进行下游分析。

第2.4节访问和修改seurat对象

创建Seurat对象后,我们可以使用各个要素访问原始数据和添加的元数据,

> exprs_matrix <-GetAssayData(sobj, slot = "counts") #访问表达矩阵
> rownames(exprs_matrix)      #访问基因名称
> sobj@meta.data     #访问元数据

这使分析过程中的每一步都处于可控状态,方便调试和跟踪。

我们可以使用”$“访问Seurat对象的各个要素,例如:

> sobj$pca  
> sobj$tSNE 
> sobj$meta

此外,也可以使用Aliased Shortcut方式访问,如:

> pca(sobj)  
> tsne(sobj)

要修改Seurat对象,也有多种方式:

1.直接覆盖要素:sobj$pca <- new_pca

2.使用Alter标志进行更新:sobj[[“pca”]] <- Alter(sobj[[“pca”]], new_pca)

3.使用seurat内置函数进行更新:PCA(sobj, new_pca)

4.重新运行命令来覆盖结果:RunPCA(sobj)

这使我们可以轻松对Seurat对象内的任意一步进行重新分析,而不用完全重新构建对象。

Seurat对象集成了单细胞数据分析的全部过程,熟练掌握创建、访问和修改Seurat对象的方法,是使用seurat包的基础。要深入理解Seurat对象,需要耐心学习并逐步掌握其中的各个要素与功能。

第3章数据预处理

第3.1节质量控制与可视化

seurat v4提供丰富的QC度量和可视化功能,可以深入了解样本和数据质量。可视化工具包括DimPlot，FeatureScatter，VlnPlot。主要的QC可视化PLOT有:

1. FeatureScatter:绘制特征-特征散点图,视觉检查相关性。

> sobj <- PercentageFeatureSet(sobj, pattern = "^MT-", col.name = "percent.mito")
> FeatureScatter(sobj, feature1 = "nFeature_RNA", feature2 = "nCount_RNA")

> FeatureScatter(sobj, feature1 = "nCount_RNA", feature2 = "percent.mito")

2. VlnPlot:绘制特征的小提琴图,检查分布。

> VlnPlot(sobj, c("nFeature_RNA", "nCount_RNA", "percent.mito"))

3. RidgelinePlot:绘制特征的脊线图,检查密度和分布。

> RidgePlot(sobj, c("nCount_RNA", "percent.mito"))
## Picking joint bandwidth of 139
## Picking joint bandwidth of 0.153

这些可视化工具可以让我们直观判断数据和细胞的质量,识别异常值和缺失值,进而指导过滤和预处理的方案。熟练运用QC可视化工具是seurat v4高级用户的必备技能。

第3.2节过滤低质量细胞

过滤低质量和异常细胞对单细胞数据的分析至关重要。

以下代码保留了基因数大于200的细胞，并保留了线粒体基因百分比小于5%的细胞

> sobj <- subset(sobj,  subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mito < 5)

也可以使用多个条件过滤细胞：

> min_genes <- 200
> max_genes <- 6000
> min_UMI_count <- 1000
> max_UMI_count <- 50000
> max_mito_pct <- 5
> 
> sobj <- subset(sobj, subset = nFeature_RNA >= min_genes & nFeature_RNA <= max_genes & nCount_RNA >= min_UMI_count & nCount_RNA <= max_UMI_count & percent.mito <= max_mito_pct)

适当的过滤可以最大限度地去除低质量和异常细胞,这是后续分析的重要一步。

第3.3节识别multiplet细胞

在单细胞RNA测序分析过程中，识别并处理multiplet细胞（多个细胞被错误地捕获并一起测序）是关键的一步。multiplet细胞会导致数据中出现技术偏差，可能会影响分析结果的准确性和可靠性。在本节中，我们将讨论识别multiplet细胞的重要性、常用方法和实际应用。

3.3.1 multiplet细胞的影响

multiplet细胞在单细胞RNA测序数据中引入了噪声，可能导致以下问题：

1.影响细胞聚类：multiplet细胞的表达谱是多个细胞表达谱的混合，可能导致错误的聚类结果。

2.错误的生物学解释：由于multiplet细胞的表达谱不是单个细胞的真实表达谱，它们可能导致错误的生物学解释和结论。

3.降低后续分析的准确性：multiplet细胞的存在可能影响差异表达分析、轨迹推断和其他分析步骤的准确性。

3.3.2 识别multiplet细胞的方法

多种方法可用于识别multiplet细胞，以下是一些常用方法：

1.基于表达谱的方法：通过比较不同细胞之间的基因表达模式，可以识别具有异常表达谱的细胞。这些异常表达谱可能暗示着它们是multiplet细胞。例如，使用主成分分析（PCA）或t分布邻域嵌入（t-SNE）等降维方法，可以在低维空间中观察到multiplet细胞聚集在多个细胞类别之间的边界区域。

2.使用有标签数据：在某些情况下，可以通过有标签数据来识别multiplet细胞。例如，使用单细胞哈希（Cell Hashing）技术，可以通过标记特定细胞并在测序过程中检测这些标签来识别multiplet细胞。

3.使用第三方工具：许多第三方工具可以用于识别multiplet细胞，如DoubletFinder、Scrublet和DoubletDecon等。这些工具通常通过比较细胞之间的相似性或预测潜在的multiplet细胞来工作。

3.3.3 实际应用

在实际分析中，可以使用以下步骤来识别并处理multiplet细胞：

1.数据预处理：在分析单细胞数据之前，首先进行基本的数据预处理，如过滤低质量细胞、标准化、归一化和log变换等。

2.选择合适的方法：根据数据类型和实验设计，选择最适合的方法来识别multiplet细胞。这可能包括基于表达谱的方法、使用有标签数据或利用第三方工具。

3.应用所选方法：将所选方法应用于预处理后的数据，以识别潜在的multiplet细胞。确保仔细检查识别出的细胞，并评估它们是否可能是真正的multiplet细胞。

4.移除multiplet细胞：在识别出multiplet细胞后，将它们从数据中移除。这将有助于提高后续分析的准确性和可靠性。

5.重新分析数据：在移除multiplet细胞后，重新分析数据以获得更准确的结果。这可能包括重新进行细胞聚类、差异表达分析和轨迹推断等。

总之，在单细胞RNA测序分析过程中，识别并处理multiplet细胞是关键的一步，有助于提高分析结果的准确性和可靠性。在实际应用中，研究人员需要根据自己的数据类型和实验设计，选择并应用合适的方法来识别multiplet细胞。

综合使用多个方法可以更加准确识别Multiplet,获得高质量的单细胞数据用于后续研究。Multiplet的识别与过滤是单细胞分析比较难度较大的一块,需要不断实践与总结。

第3.4节标准化、归一化和log变换

单细胞RNA测序数据往往需要进行预处理,以消除技术效应和批效应,获得更加可比较的表达值。

常用的预处理方法包括:

1. 标准化:使表达值在相同尺度上进行比较。seurat提供ScaleData()进行标准化。

> all.genes <- rownames(sobj)
> sobj <- ScaleData(sobj, features = all.genes)
## Centering and scaling data matrix

2. 归一化:使reads数量在所有细胞之间相同或处于相似范围。seurat提供NormalizeData()进行归一化。 NormalizeData() 默认方式是 LogNormalize，其他方法有CLR,RC。LogNormalize的做法是，每个细胞，每个基因的count除以该细胞的总counts，乘以scale.factor(默认是10000，就好像所有单细胞总共有10k UMI),并对获得的值进行log1p转换，进行自然对数转换。归一化的数据存储在”RNA”assay的seurat_obj[[‘RNA’]]@data中

> sobj <- NormalizeData(sobj)

这2种预处理可以单独使用,也可以结合使用。具体选择取决于数据的质量与方差。预处理后,可以使用预处理前的表达矩阵数据sobj$data与sobj$scale.data对比,选择最佳的预处理方案。

3.高变异基因的选择

> sobj <- FindVariableFeatures(sobj, selection.method = "vst", nfeatures = 2000)
> gc()  #释放内存
> top10 <- head(VariableFeatures(sobj), 10)
> plot1 <- VariableFeaturePlot(sobj)
> LabelPoints(plot = plot1, points = top10, repel = TRUE)
## When using repel, set xnudge and ynudge to 0 for optimal results
## Warning: Transformation introduced infinite values in continuous x-axis

##             used  (Mb) gc trigger   (Mb)  max used   (Mb)
## Ncells   3479891 185.9    6388828  341.3   6388828  341.3
## Vcells 102905979 785.2  335896301 2562.7 326150178 2488.4

4.也可以使用sctransform包进行标准化处理，单个命令将替换NormalizeData、ScaleData和FindVariableFeatures。

> 
> sobj <- SCTransform(sobj, vars.to.regress = "percent.mito", verbose = FALSE, conserve.memory=TRUE)
> gc()  #释放内存
> #> sobj <- SCTransform(sobj, vars.to.regress = "percent.mito", verbose = FALSE, conserve.memory=TRUE)
> #Will not return corrected UMI because residual type is not set to 'pearson'

这个警告信息并不会影响到SCTransform函数的执行。实际上，即使不返回校正后的UMI值，您仍然可以继续进行后续的单细胞分析流程。这个警告信息是为了让您了解到，在使用非Pearson残差方法时，SCTransform不会返回校正后的UMI值。

预处理使表达值在不同样本之间具有可比性,这是进行下游分析的前提。选择适当的预处理方法可以最大限度地消除技术差异,获得真实的生物学信号。但过度预处理也会造成信息损失,需谨慎操作。

第4章集群分析

4.1 主成分分析(PCA)

PCA是一种无监督的线性降维方法,通过正交变换将高维数据映射到低维空间,最大限度地保留数据的方差。使用RunPCA()在seurat中运行PCA:

> sobj <- RunPCA(sobj, features = VariableFeatures(object = sobj))
## PC_ 1 
## Positive:  FTL, FTH1, CST3, AIF1, FCER1G, S100A4, TYROBP, COTL1, LGALS1, S100A6 
##     OAZ1, LST1, ACTB, PSAP, SAT1, S100A11, TYMP, TMSB4X, CTSS, SERPINA1 
##     SPI1, LYZ, IFITM3, FCN1, CFD, RP11-290F20.3, TIMP1, HLA-DRB1, CFP, HLA-DPA1 
## Negative:  MALAT1, IL32, LTB, IL7R, CCL5, CTSW, GZMA, NKG7, CD247, CD2 
##     GZMK, CST7, ACAP1, GZMH, HOPX, LINC00926, BEX2, FGFBP2, TCL1A, GNLY 
##     NCR3, MYC, STK17A, AQP3, SAMD3, SLC2A3, SPON2, ZAP70, PRF1, MS4A1 
## PC_ 2 
## Positive:  B2M, NKG7, GZMA, CTSW, GNLY, PRF1, FGFBP2, CCL5, SPON2, GZMB 
##     CST7, AKR1C3, GZMH, MALAT1, CLIC3, FCGR3A, HOPX, CD247, IL32, XCL2 
##     ACTB, TMSB4X, IFITM2, CTSC, UBB, XCL1, CCL4, TTC38, TUBA1B, ARPC5L 
## Negative:  S100A8, S100A9, LYZ, CD14, LGALS2, FCN1, MS4A6A, RBP7, FTL, GSTP1 
##     S100A12, TYROBP, S100A6, TYMP, CST3, FOLR3, ID1, FCGR1A, LGALS3, IL8 
##     ASGR1, PLBD1, RETN, LGALS1, ALDH2, FCGRT, CFP, APOBEC3A, NCF2, MNDA 
## PC_ 3 
## Positive:  B2M, NKG7, GZMA, CD74, GNLY, PRF1, HLA-DPB1, HLA-DPA1, CTSW, HLA-DRA 
##     MALAT1, FGFBP2, SPON2, HLA-DQA1, HLA-DRB1, GZMB, HLA-DQB1, CST7, HLA-DQA2, HLA-DRB5 
##     CLIC3, PLAC8, AKR1C3, GZMH, HOPX, UBB, CD1C, CD247, XCL2, VIM 
## Negative:  PPBP, GNG11, SPARC, PF4, AP001189.4, ITGA2B, CD9, SDPR, CLU, GP9 
##     TREML1, NRGN, LY6G6F, CMTM5, TUBB1, RGS18, F13A1, RP11-367G6.3, GP1BA, C6orf25 
##     GPX1, CA2, SCGB1C1, RUFY1, CLDN5, SEPT5, CLEC1B, ITGB3, HIST1H2AC, HGD 
## PC_ 4 
## Positive:  CD74, HLA-DQA1, HLA-DQB1, HLA-DQA2, HLA-DRA, HLA-DPB1, HLA-DRB1, HLA-DPA1, HLA-DRB5, CD1C 
##     CD79A, MS4A1, HLA-DMB, CD79B, HLA-DMA, LINC00926, CLEC10A, TCL1A, FCER1A, PLD4 
##     ITM2C, PPP1R14A, IRF8, LTB, SERPINF1, BLNK, P2RX5, MZB1, TNFRSF13B, TUBA1B 
## Negative:  FCGR3A, NKG7, FTL, S100A4, RP11-290F20.3, FCER1G, TYROBP, LST1, GZMA, AIF1 
##     S100A8, SERPINA1, S100A6, CFD, GNLY, FGFBP2, S100A9, SPON2, PRF1, MS4A7 
##     CTSS, CTSW, S100A11, FTH1, CST7, FCN1, PSAP, GZMH, IFITM2, HMOX1 
## PC_ 5 
## Positive:  LYZ, NKG7, S100A9, S100A8, GNLY, GPX1, GZMA, CD14, SPON2, FGFBP2 
##     LGALS2, GZMB, PRF1, GSTP1, CTSW, CST7, AKR1C3, CLIC3, CCL3, GZMH 
##     XCL2, MS4A6A, CCL5, CCL4, ID1, FCN1, HOPX, XCL1, TTC38, IGFBP7 
## Negative:  FCGR3A, MALAT1, RP11-290F20.3, MS4A7, LTB, IFITM2, CKB, LST1, HMOX1, SIGLEC10 
##     AIF1, TIMP1, TMSB4X, IL7R, SERPINA1, COTL1, FCER1G, VMO1, MT-CO2, PPM1N 
##     CYTIP, CD79B, CEBPB, FAM110A, IFITM3, CTD-2006K23.1, GSTA4, LYN, B2M, PTGES3
> VizDimLoadings(sobj, dims = 1:2, reduction = "pca")

> DimPlot(sobj, reduction = "pca")

> DimHeatmap(sobj, dims = 1, cells = 500, balanced = TRUE)

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

> print(sobj[["pca"]], dims = 1:5, nfeatures = 5)
## PC_ 1 
## Positive:  FTL, FTH1, CST3, AIF1, FCER1G 
## Negative:  MALAT1, IL32, LTB, IL7R, CCL5 
## PC_ 2 
## Positive:  B2M, NKG7, GZMA, CTSW, GNLY 
## Negative:  S100A8, S100A9, LYZ, CD14, LGALS2 
## PC_ 3 
## Positive:  B2M, NKG7, GZMA, CD74, GNLY 
## Negative:  PPBP, GNG11, SPARC, PF4, AP001189.4 
## PC_ 4 
## Positive:  CD74, HLA-DQA1, HLA-DQB1, HLA-DQA2, HLA-DRA 
## Negative:  FCGR3A, NKG7, FTL, S100A4, RP11-290F20.3 
## PC_ 5 
## Positive:  LYZ, NKG7, S100A9, S100A8, GNLY 
## Negative:  FCGR3A, MALAT1, RP11-290F20.3, MS4A7, LTB

DimPlot()绘制PCA图可视化降维结果:,reduction =“pca”指定绘制PCA图。可以查看样本在低维空间的分布,识别聚集模式与分隔。

PCA检验可通过ElbowPlot()选择主成分数;通过JackStrawPlot()检验PC的统计显著性。这些工具可以更科学合理地选择主成分,使其最大限度保留数据信息。

DimHeatmap():观察同一聚类内细胞间的距离heatmap。热点代表聚类核心,冷点可能被错误聚类。

4.2 tSNE t-stochastic neighbor embedding

tSNE是一种非线性降维方法,通过随机近邻嵌入保留局部数据结构,将高维数据映射到二维或三维。使用RunTSNE()在seurat中运行tSNE:

> sobj <- RunTSNE(sobj, dims = 1:30)   #使用前30个维度
> DimPlot(sobj, reduction = "tsne")

DimPlot()绘制tSNE图,reduction = “tsne”指定绘制tSNE图。tSNE图可识别局部聚类模式,但全局结构可能失真。

4.3 UMAP Uniform Manifold Approximation and Projection

UMAP是一种非线性降维方法,结合全局和局部维度缩减技术,生成二维或三维布局。使用RunUMAP()在seurat中运行UMAP:

> sobj <- RunUMAP(sobj, dims = 1:30) #使用1-30维数据 
## 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
## 20:37:46 UMAP embedding parameters a = 0.9922 b = 1.112
## 20:37:46 Read 2638 rows and found 30 numeric columns
## 20:37:46 Using Annoy for neighbor search, n_neighbors = 30
## 20:37:46 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 20:37:46 Writing NN index file to temp file C:\Users\Administrator\AppData\Local\Temp\RtmpeoxNlr\file55c16151050
## 20:37:46 Searching Annoy index using 1 thread, search_k = 3000
## 20:37:46 Annoy recall = 100%
## 20:37:47 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 20:37:47 Initializing from normalized Laplacian + noise (using irlba)
## 20:37:47 Commencing optimization for 500 epochs, with 117352 positive edges
## 20:38:04 Optimization finished
> DimPlot(sobj, reduction = "umap")

DimPlot()绘制UMAP图,reduction = “umap”指定绘制UMAP图。UMAP可同时保留全局和局部结构,生成更连续更可信的低维布局,较tSNE有改进。

PCA,tSNE和UMAP是seurat中常用的降维和可视化方法。通过这些图可以识别样本间的聚集模式和群落结构,为后续的细胞注释和聚类分析提供依据。

理解这3种方法的原理,熟练掌握seurat中运行和绘图的工具,是单细胞数据分析的基本技能。

4.4 测试集群效果

选择合适的聚类分辨率和降维空间是确保聚类效果的关键。我们可以通过以下方法评估不同参数下的聚类效果.

Elbow plot:绘制聚类数与SSE的散点图,选择肘点处的聚类数。但仅当resolution在合适范围时有效。

JackStrawPlot()功能提供了一个可视化工具，用于将每个 PC 的 p 值分布与统一分布（虚线）进行比较。“重要”PC 将显示在虚线上方。另一种启发式方法生成”肘部图

> # JackStrawPlot()功能提供了一个可视化工具，用于将每个 PC 的 p 值分布与统一分布（虚线）进行比较。"重要"PC 将显示在虚线上方。
> sobj <- JackStraw(sobj, num.replicate = 100)
> sobj <- ScoreJackStraw(sobj, dims = 1:20)
> JackStrawPlot(sobj, dims = 1:15)
## Warning: Removed 21941 rows containing missing values (`geom_point()`).

另一种启发式方法生成”肘部图”：[ElbowPlot()]根据每个（函数）解释的方差百分比对原则组件进行排名。在此示例中，我们可以观察到 PC9-10 周围的”肘部”，表明大多数真实信号在前 10 个 PC 中被捕获。

> 
> ElbowPlot(sobj)

4.5 FindClusters函数进行集群

在降维空间识别的样本聚集结构为细胞注释和聚类提供参考。FindClusters()函数用于在seurat对象上执行细胞聚类。

> sobj <- FindNeighbors(sobj, dims = 1:10)
## Computing nearest neighbor graph
## Computing SNN
> sobj <- FindClusters(sobj, resolution = 0.5)
> head(Idents(sobj), 5)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 2638
## Number of edges: 90496
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8845
## Number of communities: 12
## Elapsed time: 0 seconds
## AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1 
##                4                3                7                5 
## AAACCGTGTATGCG-1 
##                4 
## Levels: 0 1 2 3 4 5 6 7 8 9 10 11

reduction指定在哪个降维空间进行聚类,resolution调整聚类分辨率,范围0-1,值越大聚类数越多。FindClusters()将在聚类标签”seurat_clusters”中返回聚类结果。使用DimPlot观察聚类结构:

> VizDimLoadings(sobj, dims = 1:2, reduction = "pca")

> DimPlot(sobj, reduction = "umap")

> DimPlot(sobj, reduction = "pca")

> DimPlot(sobj, reduction = "tsne")

> DimPlot(sobj, split.by = "seurat_clusters")

第5章鉴定不同ially expressed features

5.1 FindMarkers函数

FindMarkers()函数用于在不同细胞类型或群落之间鉴定差异表达的特征(基因或蛋白)。该函数基于差异性检验(如Wilcoxon秩和检验),计算每个特征在相比群之间的差异性评分。然后选择评分最高的特征作为群间的差异标记物。

FindMarkers()函数的主要参数有:

1.ident.1和ident.2:需要进行差异分析的两个细胞类型或群落的ID。

2.min.pct:要考虑为差异特征的最小百分比。默认为25%,表示在ident.1或ident.2内至少25%的细胞需要表达该特征。

3.thresh.use:选择差异性检验方法,可选”wilcox”,“t”,“bimod”。默认为”wilcox”,采用Wilcoxon秩和检验。

4. min.diff:要考虑为差异特征的最小logfold change差异度。默认为0。

5. test.use:差异性评分方法,可选”roc”,“roc01”,“wilcox”,“t”。默认为”roc”,采用ROC曲线下的面积为评分。

主要输出结果有:

1.features:差异特征的ID名单。

2. scores:每个特征的差异性评分。

3. logfc:每个特征在两组之间的log2 fold change值。

4. p.value:每个特征对应的p值。

5. percent.1和percent.2:每个特征在ident.1和ident.2中的表达百分比。

FindMarkers()函数可以高效地在不同细胞类型或群落间鉴定差异表达特征,为理解不同细胞状态或机制下调控网络的变化提供了重要工具。然而,结果的准确性与稳定性依赖于数据质量,群落划分的准确性,以及参数的合适选择。

> library(tidyverse)
> cluster2.markers <- FindMarkers(sobj, ident.1 = 2, min.pct = 0.25)
> cluster5.markers <- FindMarkers(sobj, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
> sobj.markers <- FindAllMarkers(sobj, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 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

> sobj.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)
> FeaturePlot(sobj, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", 
+                                "CD8A"))

> top10 <- sobj.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)
> DoHeatmap(sobj, features = top10$gene) + NoLegend()

## # A tibble: 24 × 7
## # Groups:   cluster [12]
##        p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene    
##        <dbl>      <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>   
##  1 1.85e- 86      1.31  0.746 0.301 6.05e- 82 0       IL7R    
##  2 6.28e- 41      1.10  0.35  0.113 2.06e- 36 0       AQP3    
##  3 2.93e-277      4.31  0.998 0.24  9.58e-273 1       S100A9  
##  4 5.95e-227      4.01  1     0.532 1.95e-222 1       LYZ     
##  5 1.41e- 37      1.05  0.422 0.158 4.60e- 33 2       CCR7    
##  6 2.54e- 13      0.847 0.388 0.238 8.33e-  9 2       LEPROTL1
##  7 0              4.29  0.933 0.042 0         3       CD79A   
##  8 9.48e-271      3.59  0.622 0.022 3.10e-266 3       TCL1A   
##  9 3.75e-200      2.84  0.936 0.208 1.23e-195 4       NKG7    
## 10 9.21e-190      3.12  0.915 0.232 3.02e-185 4       CCL5    
## # ℹ 14 more rows

5.2 绘制 volcano plot进行差异基因鉴定

Volcano plot是一种直观显示差异表达分析结果的图形方法。它同时考虑了每个特征的fold change值(效应量)和统计显著性(p值),从而鉴定同时具有较大效应量和高统计显著性的差异表达特征。

要绘制volcano plot,我们需要设置如下参数:

1. features: FindMarkers()输出的特征ID列表。

2. logfc: 每个特征的log2 fold change值。

3. p_val: 每个特征对应的p值。

4. p_cutoff: 考虑为差异特征的p值切offs,默认为0.05。

5. fc_cutoff: 考虑为差异特征的fold change阈值,默认为0。

主要代码如下:

> p_cutoff=0.05
> fc_cutoff=0
> ggplot(sobj.markers, aes(x=avg_log2FC, y=-log10(p_val))) + 
+   geom_point() +
+   xlab("log2 Fold Change") + ylab("-log10 p-value") +
+   ggtitle("Volcano Plot") +  
+   theme_bw()+
+   geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "blue", linewidth = 0.5) +
+   geom_hline(yintercept = 1.3, linetype = "dashed", color = "blue", linewidth = 0.5)

这个代码将在一张散点图中同时显示每个特征的log2 fold change值(x轴)和-log10转换后的p值(y轴)。同时添加水平和垂直参考线表示选择阈值。并用颜色高亮同时满足fold change和p值阈值要求的差异特征。

Volcano plot可以将FindMarkers()结果直观显示,辅助我们选择具有较高效应量和显著统计差异的特征。选择适当的阈值cutoffs,我们可以在该图中鉴定关键的差异表达特征。