Seurat数据结构（一）

在做单细胞分析的过程中，我们通常要基于Seurat数据之上做多样化的分析：针对某一子集，重新画图等等。所以，我们首先需要对Seurat的数据结构有一个基本的了解。

下面这个图是对于Seurat对象的总结，可以看出是个树状的数据结果，二级结构主要包含Assays、meta.data、active.assay、active.ident、graphs、reduction、commands等。

图片来源：https://www.jianshu.com/p/13142bf51e81

下面我们一个个介绍一下这些基本数据对象。

=======Assays=======

默认情况下，Seurat对象是一个叫RNA的Assay。在我们处理数据的过程中，做整合（integration），或者做变换（SCTransform），或者做去除污染（SoupX），或者是融合velocity的数据等，都会生成新的相关的Assay，用于存放这些处理之后的矩阵。

在之后的处理中，我们可以根据情况使用指定Assay下的数据。不指定Assay使用数据的时候，Seurat调用的是Default Assay下的内容。我们可以通过对象名@active.assay查看当前Default Assay，通过DefaultAssay函数更改当前Default Assay。

Assay数据中，counts为raw原始数据，data为normalized（归一化），scale.data matrix为scaled（标准化的数据矩阵）。

调用方法：PBMC@assays$RNA@data：调用normalized后数据。

=====meta.data======

这个结构也叫元数据，包含对每个细胞的描述。一般计算的nFeature_RNA等信息就以metafeature的形式存在Seurat对象的metadata中。计算的分类信息一般以RNA_snn_res.x（x指使用的resolution）存放在metadata中。

=====active.assay和active.ident======

前者是查看当前使用的assays，后者是查看当前的使用分群方式。

=====reduction======

用于储存降维之后的每个细胞的坐标信息。

每个细胞在PC轴上的坐标：pbmc@[email protected]

每个基因对每个PC轴的贡献度（loading值）：pbmc@[email protected]

下面，我们介绍几个常用的函数。

获取全部基因ID：rownames(object)

获取整个object的细胞ID：Cells(object)或者colnames(object)

按照idents获取部分细胞ID：WhichCells(object, idents = c(1, 2))

行数和列数：ncol(object) nrow(object)

每个cell的类型和class：Idents(object)  levels(object)

按照基因表达获取部分细胞ID：

WhichCells(object, expression = gene1> 1)

WhichCells(object, expression = gene1> 1, slot = "counts")

提取降维之后的坐标信息：

Embeddings(object =object[["pca"]])

Embeddings(object =object[["umap"]])

感觉和这样子类似：pbmc@[email protected]

如果要提取包含部分细胞的对象：

按照细胞ID提取：

subset(x = object, cells = cells)

按照idents提取：

subset(x = object, idents = c(1, 2))

或者前面用过的WhichCells(object, idents = 1)

subset(x=object, idents = "root")  #对细胞簇重新命名后为root

想要排除1、2细胞类型，可以这样：

subset(pbmc, idents = c(1, 2), invert =TRUE)

按照meta.data中设置过的stim信息提取：

subset(x = object, stim == "Ctrl")

按照某一个resolution下的分群提取：

subset(x = object, RNA_snn_res.2 == 2)

当然还可以根据某个基因等的表达量来提取：

subset(x = object, gene1 > 1)

subset(x = object, gene1 > 1，slot ="counts")

其实和前面whichcell的提取思想类似，只不过一个提取的cell的列表，一个提取的是seurat的子集。

在子集操作之上，我们就可以进行其它的函数等分析了。

例如：

查看每个聚类包含多少细胞：

table(Idents(pbmc))，table(pbmc$RNA_snn_res.0.5)等

每个聚类细胞数占比：

prop.table(table(Idents(pbmc)))，prop.table(table(pbmc$RNA_snn_res.0.5))

提取表达矩阵：

as.matrix(pbmc@assays$RNA@counts)  

或者

as.matrix(GetAssayData(pbmc, slot = "counts")）

提取某种细胞的表达矩阵：（在上面的基础上和whichcell结合）

as.matrix(GetAssayData(pbmc, slot = "counts")[, WhichCells(pbmc, ident = "1")])

计算平均表达量（用AverageExpression函数）：

cluster.averages <- AverageExpression(pbmc)

修改聚类后的因子水平（这些基本是常规的R操作了）：

Idents(pbmc) <- factor(Idents(pbmc), levels= c(1,2,3,4,9,8,7,6,5,0))