1.1 seurat包概述

seurat包是一款R语言的包,主要用于分析单细胞RNA测序数据。它由Satija实验室开发,于2015年首次发布,现已更新到3.0版本。seurat包实现了单细胞RNAseq数据分析的主要步骤,包括质控、归一化、可视化、聚类、差异表达等,并有丰富的图形化界面,十分易于使用。

seurat包的基础是AnnotatedChip类,它是一个S4类,存储了单细胞RNAseq数据及所有相关分析结果。我们可以通过AnnotatedChip对象对数据进行操作,实现各种分析及结果访问。seurat包通过定义AnnotatedChip类的子类,如AnnotatedSCE类,可以支持多种单细胞数据格式,大大增加其兼容性和应用性。

seurat包的主要功能包括:

1. 从原始gene表达矩阵构建AnnotatedChip对象,或从已有格式如SingleCellExperiment对象进行转换。

2. 用于质控的功能,如过滤低质量细胞,检测高变基因等。

3. 用于数据标准化的方法,如缩放,归一化等。

4. 用于数据降维和可视化的工具,如PCA,tSNE,UMAP等。

5. 用于数据聚类的方法,如Louvain算法,K-means等。

6. 用于差异表达分析的功能, findsMarkers()等。

7. 支持进行GO/KEGG富集分析,蛋白质相互作用网络分析等。

8. 提供丰富的绘图功能,如特征作图,等高线图,热图等。

9. 支持对象的导出与共享,便于重复分析或与他人交流。

1.2 AnnotatedChip类介绍

AnnotatedChip类是seurat包的基础,它是一个S4类,用于存储单细胞RNAseq数据及所有相关分析结果。一个AnnotatedChip对象包含多个域(slots),存储不同类型的信息。

AnnotatedChip对象包含以下主要域:

1. assays:存储数据集的表达矩阵及相关信息。

2. meta_data: 存储细胞元数据,如细胞类型,count数等。

3. active.assay:设置当前使用的表达矩阵,这里是RNA。

4. active.ident:细胞ID,对应列名。

5. graphs:网络图信息,用于存储PCA,tSNE等结果。

6. reductions:用于存储降维结果,如PCA,UMAP的模型对象。

7. images:存储绘图 Output,用于再现分析过程。

8. project.name:项目名称。

9. version:seurat对象的版本信息。

10. commands:存储构建seurat对象使用的所有命令,以支持再现分析过程。

11. tools:存储分析中使用的其他软件包版本信息。

通过这些域的存储,一个AnnotatedChip对象包含了单细胞数据集从原始表达矩阵到聚类与差异表达分析的全部信息,并支持这些信息的提取与再分析。这使得工作流程高度连贯且可重复。

我们可以通过多种方式访问这些域,如@运算符,slot()函数和中括号运算符[]。同时,seurat包定义了许多S4泛型函数用于操作AnnotatedChip对象,如 subset()进行子集选择,pull()提取计数数据,add_test_data()添加新的域等。这些函数与域的存储机制共同构成了AnnotatedChip对象的操作接口。

AnnotatedChip对象还具有Validity方法,用于检查对象各域的正确性与完整性。这可以确保结果的可靠性和稳定性。通过设置as和contains方法,AnnotatedChip对象具有较强的扩展性,我们可以构建子类来支持多种数据类型。

1.3 seurat工作流程概览

seurat包的工作流程主要基于AnnotatedChip对象的操作与分析。一个标准的seurat工作流程通常包括以下步骤:

1. 构建AnnotatedChip对象:从原始表达矩阵或其他格式构建AnnotatedChip对象。这一步将数据载入seurat包,建立分析基础。

2. 数据集质控:包括过滤低质量细胞,检测高变基因等。通过设置min.cells和min.genes参数控制数据集的深度与规模。

3. 数据标准化:使用缩放(ScaleData())或归一化(NormalizeData())方法减少样本间的系统误差。

4. 降维与可视化:使用PCA、tSNE、UMAP等方法对数据进行降维,并通过特征图、散点图等可视化结果。这可以直观判断样本间的异质性与相似性。

5. 数据聚类:通过FindClusters()方法和相应的聚类算法,如Louvain算法,对数据进行聚类并获得聚类标签。

6. 差异表达分析:通过FindMarkers()方法检测不同聚类或细胞类型之间的差异表达基因。

7. GO/KEGG富集分析:通过EnrichR包等对差异表达基因进行GO或KEGG富集分析,发现重要的生物学路径或进程。

8. 蛋白质相互作用网络分析:通过STRINGdb包构建差异表达基因对应的蛋白质相互作用网络,发现重要的功能模块。

9. 结果导出与共享:通过ExportToCDS(), ExportToLoom()或SaveRDS()方法导出AnnotatedChip对象,以支持结果分享、重复分析或Meta分析。

这一工作流程覆盖了单细胞RNAseq分析的主要过程。研究人员可以根据自己的数据集与研究问题选择部分步骤,或根据需要加入其他扩展步骤。seurat包提供的丰富方法与完备功能,可以支持大部分标准分析并具备定制化的灵活性。

2.1 slots域详解

AnnotatedChip对象包含以下主要域:

1. assays:存储数据集的表达矩阵及相关信息。这里只有一个RNA域,对应RNA表达数据。它包含:

- counts:原始表达计数矩阵,行为基因,列为细胞。

- data: 与counts相同,用于兼容以前的seurat版本。 数据进行了VST变换后,存储在data域,而非scale_data域。

- scale.data: 标准化后的表达矩阵。

- key:表达矩阵的标识符,此处为”rna_“。

- assay.orig:原始表达矩阵,此处为空。

- var_features:变异基因信息。如果进行筛选，例如选定了2000个变异基因,存储在var_features域

- meta.features:包含每个特征的统计量,如平均值、方差、标准化方差等,用于过滤低变异特征。

- misc:杂项,此处为空。

2. meta_data: 存储细胞元数据,如细胞类型,count数等。meta_data域包含更丰富的信息,如percent.mito,聚类标签等

- orig.ident:样品的原始ID,此处全部为”SeuratProject”。

- nCount_RNA:每个细胞的UMI计数,在RNA表达矩阵中。

- nFeature_RNA:每个细胞检测到的特征数,在RNA表达矩阵中。

- percent.mito:每个细胞线粒体特征的百分比,用于过滤低质量细胞。

- RNA_snn_res.0.5:构建RNA_snn网络图时resolution为0.5时的聚类结果。

- seurat_clusters:Seurat的FindClusters命令产生的聚类结果,resolution同为0.5。

这个域汇总了样品和表达数据的基本信息,以及两种不同参数下的聚类结果。

3. active.assay:设置当前使用的表达矩阵,这里是RNA。

4. active.ident:细胞ID,对应列名。active.ident设置为聚类标签,而非原始ID

5. graphs:网络图信息,用于存储PCA,tSNE等结果。

- RNA_nn:最近邻网络图。

- RNA_snn:共享最近邻网络图。

它们包含节点(细胞)与边(相似性)的信息。

6. reductions:用于存储降维结果,如PCA,UMAP的模型对象。

- cell.embeddings: 细胞在低维空间的坐标。

- feature.loadings:每个维度对原始特征的贡献度。只有PCA包含。

- assay.used:使用的表达矩阵,这里是RNA。

- global:是否是全局分析,tSNE和UMAP为TRUE。

- stdev:标准差。只有PCA包含。

- key:降维 techniuqes的名称,如PC_、tSNE_ 和UMAP_。

- jackstraw:置换检验的结果,评估每个维度的显著性。只有PCA包含。

- misc:其他信息,如PCA的总方差。

7. images:存储绘图 Output,用于再现分析过程。

8. project.name:项目名称。

9.misc:杂项,此处为空。

10. version:seurat对象的版本信息。

11. commands:存储构建seurat对象使用的所有命令,以支持再现分析过程。

- ScaleData:标准化表达矩阵。

- NormalizeData:对表达矩阵进行归一化。

- FindVariableFeatures:找到高变异特征。

- RunPCA:进行PCA分析。

- RunTSNE:进行tSNE分析。

- RunUMAP:进行UMAP分析。

- JackStraw:进行置换检验。

- ScoreJackStraw:给置换检验结果打分。

- FindNeighbors:构建邻居网络图。

- FindClusters:进行聚类分析。

每个命令都包含name、time.stamp、assay.used、call.string和params等域,记录命令名称、运行时间、使用的表达矩阵、命令语句和参数等信息。

11. tools:存储分析中使用的其他软件包版本信息。

2.2 data域:原始计数矩阵

data域包含原始的表达计数矩阵,它使用dgCMatrix类表示,包含以下slot:

1. i:非零元索引,表示计数矩阵中非零元的行和列索引。

2. p: 指向i中每列第一个非零元的索引,用于快速提取某一列的数据。

3. x: 非零元对应的计数值。

4. Dim: 矩阵的行数(特征数)和列数(样品数)。

5. Dimnames: 行名和列名,对应特征名和样品名。

这个域存储了数字基因表达计数矩阵的原始数据。这是进行下游分析的基础,数据的质量直接决定结果的准确性。

为了获得高质量的数据,需要在分析前进行严格的质控。主要步骤如下:

1. 去除低质量样品:过滤UMI计数太低或检测特征太少的样品。

2. 去除线粒体和ри波体RNA:这些RNA的表达丰度通常很高,会影响后续分析。

3. 选择高质量特征:过滤低变异或高缺失率特征。

4. 归一化:消除样品间的系统差异,常用方法有LogNormalize、CLR等。

5. 特征选择:选择变异较大、对样品分类有较高区分度的特征,以减小维度。

6. 标准化:将特征表达值缩放到相近范围,便于不同特征间的比较。

经过上述处理,我们可以获得一个高质量的表达矩阵,作为Seurat下游分析的输入。这需要对AnnotatedChip对象中各域的信息与处理结果进行综合判断,选择最优的参数与流程。

2.3 meta_data域:样本元数据

meta_data域包含样品的元数据,以data.frame的形式存储,主要包含以下信息:

1. orig.ident:样品的原始ID,用于在不同数据集间建立对应关系。

2. nCount_RNA:RNA表达矩阵中每个样品的UMI计数,用于过滤低质量样品。

3. nFeature_RNA:RNA表达矩阵中每个样品检测到的特征数,也用于过滤低质量样品。

4. percent.mito:每个样品线粒体特征的百分比,用于过滤线粒体高表达样品。

5. 测序相关信息:比如测序深度、reads质量等,用于判断数据质量。

6. 实验相关信息:比如时间点、处理组别、来源等,用于后续分析与解读。

7. 聚类/降维结果:包含不同参数下产生的结果,用于比较选择最佳参数。

这个域为下游分析提供了重要的上下文信息。通过其中的测序信息和统计量,我们可以进行严格的质控,选择高质量样品。这直接决定后续结果的准确性。

利用实验设计信息,我们可以找到不同处理组别之间的差异表达特征,进行差异表达分析与解读。这需要将meta_data与表达数据进行综合分析,找到生物学意义上有价值的信息。

聚类或降维结果的比较可以帮助我们确定最佳的算法及参数。但这需要根据实验设计和生物知识进行判断,简单地选择可解释性高的结果是不够的。

另外,通过orig.ident,我们可以将meta_data与其他数据集的元数据相结合,进行更加全面深入的生物信息挖掘与解析。

2.4 scale_data域:规范化表达矩阵

scale_data域包含标准化后的表达矩阵,它使用矩阵的形式存储,包含以下内容:

1. 行为特征,列为样品。每个值表示该特征在对应样品中的标准化后的表达量。

2. Dimnames存储行名和列名,对应特征名和样品名。

标准化是为了将不同特征的表达值缩放到相近的范围,使得这些特征在后续的降维与聚类分析中具有相近的权重。主要方法有:

1. Z-score标准化:用特征在所有样品中的均值和标准差进行标准化。公式为:

Z = (X - μ) / σ

其中X为原表达值,μ为均值,σ为标准差。

2. Robust Z-score标准化:用中位数和中位数绝对偏差进行标准化,更 robust 。

3. Max-Min标准化:用特征的最大最小值将其缩放到[0,1]范围。公式为:

X’ = (X - min) / (max - min)

4. Percentile标准化:用指定的上下四分位数将表达值缩放。

标准化可以弥补不同特征原始表达范围差异过大的问题,但也会导致某些特征的变化信息 loss 。所以,在进行标准化前,我们要根据特征的生物学属性选择最适合的方法与参数。

scale_data域提供了标准化后高质量的表达矩阵,它是Seurat中大多数分析命令的输入,直接决定结果的准确性。所以,理解不同标准化方法的理论基础与效果,并根据数据集选择最优方法,是成为AnnotatedChip高级用户的关键。

2.5 tsne_mappings和umap_mappings域

这两个域分别包含t-SNE和UMAP两种降维技术产生的结果,主要包含:

1. tsne_mapping:t-SNE结果,包含以下内容:

- x:样品在t-SNE空间的x坐标。

- y:样品在t-SNE空间的y坐标。

- Dim:矩阵的维度,此处为2维。

- Dimname:样品名。

2. umap_mapping:UMAP结果,内容与tsne_mapping相似。

这两个降维技术可以将高维数据集压缩到二维或三维空间,使得我们可以直观观察和理解复杂数据之间的局部结构、聚类现象等。但是,不同技术及参数会产生不同的降维结果,我们需要根据实验设计和生物知识进行选择。

主要步骤如下:

1. 选择高质量表达矩阵作为输入:对表达数据进行严格质控和标准化。

2. 选择降维技术和参数:t-SNE更适合展示局部结构,UMAP更易于全局理解;降维结果受perplexity等参数影响较大。

3. 与聚类结果联合解释:不同的降维结果与聚类结果可以相互验证。但是简单匹配是不够的,需要考虑到生物学意义。

4. 与标记信息结合:将聚类或降维结果与实验设计、样品来源等信息结合,可以找到更有价值的生物学结论。

5. 结果重复性与稳定性分析:多次运行同模型判断结果的一致性,避免产生偶然结果。

所以,理解降维结果和生物信息之间的相互作用机制,选择适当的技术和参数,进行多角度分析与验证,是熟练理解tsne_mappings和umap_mappings域并成为AnnotatedChip高级用户的关键。

2.6 访问AnnotatedChip对象域

AnnotatedChip对象包含多个域,存储着丰富的生物信息。访问和提取这些信息是进行下游分析的基础,主要有以下方法:

1. @:直接使用@可以第一层slot,如:

> library(Seurat)
## Attaching SeuratObject
> sobj <- pbmc_small
> counts <- sobj@assays
> counts <- sobj@meta.data

提取meta.data域,获得样品元数据。

2. $:直接使用$可以访问域,如:

> counts <- sobj$RNA

3.结合$和@访问第三层slot：

> counts <- sobj$RNA@counts

提取counts域,获得原始表达矩阵。

4. slotNames():列出所有域的名称。

5. 使用[[”]] 可以访问域,并允许域名包含特殊字符(如.),和$作用类似。

> counts <- sobj[["nCount_RNA"]]
> counts <- sobj[["RNA"]]

6. VariableFeatures():提取高变异特征。

7.使用validObject()函数检查类定义中指定的所有槽是否存在，并且槽中的对象是否来自所需的类或该类的子类

> validObject(sobj)
## [1] TRUE

3.1 从原始矩阵构建AnnotatedChip对象

1. 读取原始计数矩阵:使用Seurat提供的Read10x()和Read10x_h5()将10x基因组数据读取为dgCMatrix对象,其中i、j和x分别存储非零元的行索引、列索引和表达值。

2. 创建AnnotatedChip对象:使用CreateSeuratObject()将dgCMatrix对象转换为AnnotatedChip对象,主要参数为:

- raw.data:dgCMatrix对象,原始计数矩阵。

- min.cells:过滤掉在少于min.cells个细胞中检测到的基因,防止低质量基因对后续分析的影响。

- project:项目名称,会添加到对象的名称中。

- assay:指定raw.data将添加到哪个assay中,默认为”rna”。

3. 质控和过滤:根据质量控制标准过滤低质量的细胞和基因。主要依据为:

- 对细胞:检测基因的数量、UMI的数量、线粒体基因比例等。

- 对基因:缺失率、平均表达、变异性等。

4. 数据归一化:使用NormalizeData()将过滤后的表达矩阵归一化,消除技术相关的变异。常用方法有LogNormalize、CLR等。

5. 选择高变异基因:使用FindVariableFeatures()选出高变异高质量的基因,用于下游分析。主要依据为平均表达、变异系数等。

6. 标准化:使用ScaleData()将表达矩阵进行标准化,shrinkage不同基因的表达范围,使它们在降维和聚类中具有相近的权重。常用方法有Z-score、Robust Z-score等。

7. 其他分析:构建AnnotatedChip对象后,可以进行差异表达分析、降维聚类、差异路径分析等。

3.2 从SingleCellExperiment对象转换

SingleCellExperiment是Bioconductor包中用于存储单细胞数据的基本对象,它包含原始计数矩阵、基因数据、样本数据等。我们可以使用以下步骤将其转换为AnnotatedChip对象:

1. 安装BiocManager包,用于安装Bioconductor包。使用install.packages(“BiocManager”)进行安装。

2. 安装SingleCellExperiment包。使用BiocManager::install(“SingleCellExperiment”)进行安装。

3. 读取单细胞数据为SingleCellExperiment对象。使用read10x()等函数读取10x数据。

4. 创建AnnotatedChip对象。使用CreateSeuratObject()函数,主要参数为:

- object:输入的SingleCellExperiment对象。

- assay:指定导入哪个assay中的数据,默认为”counts”。

- min.cells:过滤在少于min.cells个细胞中检测到的基因。

- project:AnnotatedChip对象的名称。

5. 其他步骤与从dgCMatrix构建对象相同,包括质控、过滤、归一化、选择变量基因、标准化等。

所以,理解SingleCellExperiment与AnnotatedChip这两个对象,并掌握对象间的转换方法,可以让我们利用这两个包的优势进行单细胞数据分析。这需要我们对两者的原理、数据结构与分析方法有比较全面的理解,方能在不同工具间自由转换,选择最优的分析流程。

3.3 子类:AnnotatedSCE类

AnnotatedSCE是AnnotatedChip的子类,它继承了AnnotatedChip的所有属性与方法,并且可以直接从SingleCellExperiment对象构建,无需进行对象转换。主要步骤如下:

1. 安装SingleCellExperiment包,使用BiocManager::install(“SingleCellExperiment”)。

2. 读取单细胞数据为SingleCellExperiment对象,使用read10x()等函数。

3. 创建AnnotatedSCE对象。使用CreateSeuratObject()函数,主要参数为:

- object:输入的SingleCellExperiment对象。

- assay:指定导入哪个assay中的数据,默认为”counts”。

- min.cells:过滤在少于min.cells个细胞中检测到的基因。

- project:AnnotatedSCE对象的名称。

4. 其他步骤与AnnotatedChip对象相同,包括质控、过滤、归一化、选择变量基因、标准化等。

AnnotatedSCE在AnnotatedChip的基础上,增加了以下优点:

1. 简化工作流程:无需进行对象转换,可以直接从SingleCellExperiment构建,简化了分析流程。

2. 利用SingleCellExperiment的优势:可以直接利用SingleCellExperiment进行质量评估、批效应校正等,获得更高质量的分析对象。

3. 结果保存:分析结果可以直接以SingleCellExperiment对象的形式保存与导出。

4. 与其他Bioconductor包兼容:可以与其他Bioconductor包无缝结合,扩展分析方法。

所以,AnnotatedSCE整合了AnnotatedChip和SingleCellExperiment两个包的优势,让我们可以直接构建和保存分析对象,扩展分析方法,简化工作流程。这需要我们对这两个对象及相关软件包有比较全面深入的理解,熟练掌握AnnotatedSCE的构建与应用方法,进行生物学假设的检验与验证。

4.1 提取数据

1.直接提取数据后转换成数据框

> expression_matrix <- sobj[["RNA"]]@counts
> expression_df <- as.data.frame(expression_matrix)
> expression_matrix <- sobj@meta.data

2.使用FetchData()函数提取meta.data和reductions数据框的数据,slot只能是 “data”, “scale.data”, “counts”其中一个

> FetchData(object = sobj, vars = 'PC_1',slot = "counts")
> FetchData(object = sobj, vars = 'nCount_RNA')
> FetchData(object = sobj, vars = 'groups')
> FetchData(object = sobj, vars = 'PPBP')
> FetchData(object = sobj, vars = c("tSNE_1","tSNE_2"),slot = "counts")
> FetchData(object = sobj, vars = c("tSNE_1","tSNE_2"),slot = "data")
> FetchData(object = sobj, vars = c("tSNE_1","tSNE_2"),slot = "scale.data")

3.使用subset():按meta.data中的因素对对象进行子集选择,过滤后assays和meta.data两个slot的矩阵都会发生变化

> sobj <- PercentageFeatureSet(sobj, pattern = "^MT-", col.name = "percent.mito")
> sobj <- subset(sobj, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mito < 5)

要熟练掌握这个子集选择方法,我们必须清晰理解AnnotatedChip对象的数据结构,尤其是metaData域所包含的信息，我们必须明确每个metaData列所包含的信息类型。

常见的数据类型包括:

- 字符型:记录样品批次、处理组别等分类信息。

- 数值型:记录样品的定量指标,如年龄、体积、浓度等。

- 逻辑型:记录样品的特征,如是否属于某组等。

4.2 在AnnotatedChip对象中添加域

AnnotatedChip对象中的信息以域(slot)的形式组织和存储。Seurat是一个R包，用于处理和分析单细胞RNA测序数据。要给一个Seurat对象添加一个名为my_data的新slot，可以使用@操作符。但首先，我们需要确保这个slot在Seurat类中是有效的。幸运的是，Seurat对象使用一个名为misc的slot来存储各种额外的用户数据。我们可以将您的my_data存储在miscslot中。

> # 为了演示方便，这里我们创建一个简单的数据框作为my_data
> my_data <- data.frame(
+   gene = c("gene1", "gene2", "gene3"),
+   value = c(1.2, 3.4, 5.6))
> # 将my_data添加到Seurat对象的misc slot中
> sobj@misc$my_data <- my_data
> # 查看misc slot中的my_data
> sobj@misc$my_data
##    gene value
## 1 gene1   1.2
## 2 gene2   3.4
## 3 gene3   5.6

注意，在添加新域前,我们必须考虑:

1. 新域的名称:应具有明确的生物学意义,反映其包含的信息类型。

2. 新域的数据类型:数值型、字符型、逻辑型等,需要与其信息内容相匹配。

3. 新域的作用:记录什么样的信息,在后续分析中扮演什么角色。

4. 新域与已有域的关系:如何与对象中其他域的数据相互作用或关联。