【单细胞测序数据分析-1】认识Seurat对象数据结构/数据格式

2021-02-04 本文已影响0人森尼啊

一文了解单细胞对象数据结构/数据格式，单细胞数据操作不迷茫。本文内容包括单细胞seurat对象数据结构，内容构成，对象的调用、操作，常见函数的应用等。

1. 单细胞seurat对象数据结构

单细胞数据结构.jpg

此图引用自周运来老师，稍作注释。原文地址https://www.jianshu.com/p/13142bf51e81。里面还有SingleCellExperiment，anndata，h5 数据结构的介绍，写的很详细，推荐阅读。

2. 数据结构及内容

Assays

默认情况下，我们是对Seurat中的RNA的Assay进行操作。可以通过@active.assay查看当前默认的assay，通过DefaultAssay()更改当前的默认assay。
结构
counts 存储原始数据，是稀疏矩阵
data存储logNormalize() 规范化的data。总表达式对每个单元格的要素表达式度量进行标准化，将其乘以比例因子（默认为10,000），并对结果进行对数转换
scale.data#存储 ScaleData()缩放后的data，此步骤需要时间久。

meta.data

元数据，对每个细胞的描述。一般的meta.data包括orig.ident, nCount_RNA, nFeature_RNA, 以及计算后的percent.mt，RNA_snn_res.0.5等

image.png
可以通过pbmc$percent.mt或pbmc[['percent.mt']]来调用、操作

reduction

降维后的每个细胞的坐标信息，包括pca，tsne，umap等

3. 对象操作

① 通过结构图上的@,$符号依次取

② 两个中括号操作，pbmc[[ ]]。

教程中，pbmc[['percent.MT']]向meta.data添加 percent.MT 这一列。
pbmc[[]]，中括号取的是上面结构图中的二级数据名称

以上两种方法的区别是？

如果取的一级结构Assays的下属内容：

无差别

 > class(pbmc[['RNA']])
 [1] "Assay"
 attr(,"package")
 [1] "Seurat"
 > class(pbmc@assays$RNA)
 [1] "Assay"
 attr(,"package")
 [1] "Seurat"

如果是一级结构meta.data里的下属内容：

返回的数据类型不同

> class(pbmc[['nCount_RNA']])
[1] "data.frame"
> class(pbmc@meta.data$nCount_RNA)
[1] "numeric"

pbmc[['nCount_RNA']] 取出来,是所有细胞的nCount_RNA，包含细胞信息，数据框

image.png
而pbmc@meta.data$nCount_RNA取出来的是单独nCount_RNA一列，是向量

image.png

如果取的一级结构里reductions的下属内容：

无差别

> class(pbmc[['pca']])
 [1] "DimReduc"
 attr(,"package")
 [1] "Seurat"
> class(pbmc@reductions$pca)
 [1] "DimReduc"
 attr(,"package")
 [1] "Seurat"

4. seurat常用函数的应用

之前写了对seurat对象的主要操作列表https://www.jianshu.com/p/b9078c71f057，这里不赘述。下面主要看数据操作。

#数据准备
library(Seurat)
pbmc <- readRDS(file = "pbmc/pbmc3k_final.rds")
set.seed(42)
pbmc$replicate <- sample(c("rep1", "rep2"), size = ncol(pbmc), replace = TRUE) #随机给样本设置两个样本，分别为rep1，rep2，添加到metadata中。replace =T，有放回抽样
DimPlot(pbmc, reduction = "umap")

为什么我画出来的和教程上的不一样？

IMG_1388.JPG

切换细胞ident，简单统计各样本细胞cluster 频率信息

#切换细胞ident，从细胞名称转为样本rep1，rep2。绘图看批次性
pbmc$CellType <- Idents(pbmc)
Idents(pbmc) <- "replicate"
DimPlot(pbmc, reduction = "umap")
Idents(pbmc) <- "CellType" #转变回来，下一步使用
#查看不同样本、细胞名称的 细胞数量
table(Idents(pbmc))
table(pbmc$replicate)
prop.table(table(Idents(pbmc)))#计算每种细胞频率
table(Idents(pbmc), pbmc$replicate)#计算每个样本、每种细胞的数量
prop.table(table(Idents(pbmc), pbmc$replicate), margin = 2)#计算每个重复、每种细胞的频率,margin=2,按照列的值 频率和为1来计算。因为列是样本，每个样本的细胞频率和为1。

选择特定细胞，取子集

#取子集的几种方法
subset(pbmc, subset = MS4A1 > 1)
subset(pbmc, subset = replicate == "rep2")
subset(pbmc, idents = c("NK", "B"))
subset(pbmc, idents = c("NK", "B"), invert = TRUE) # 提取除了NK,B细胞之外的细胞
WhichCells(pbmc, idents = "NK")# 返回是NK细胞的细胞名

提取表达矩阵GetAssayData()

GetAssayData(object, slot, assay) # slot = counts, data, scale.data
GetAssayData(object = pbmc_small[["RNA"]], slot = "data")[1:5,1:5]#出来的是稀疏矩阵，所以用as.matrix()直接转换

##①从Assay中提取
d <- as.matrix(GetAssayData(object = pbmc_small[["RNA"]], slot = "data")[1:5,1:5])
#前五个基因和前五个细胞

##②从Seurat对象中提取
f <- as.matrix(GetAssayData(object = pbmc_small, slot = "counts")[1:5,1:5])

#提取NK细胞的表达矩阵，用于其他分析
nk.raw.data <- as.matrix(GetAssayData(pbmc, slot = "counts")[, WhichCells(pbmc, ident = "NK")]) #直接提取NK细胞的原始表达矩阵

计算cluster内平均基因表达量AverageExpression()

cluster.averages <- AverageExpression(pbmc) #返回一个list元素的list。这个list元素是data.frame的list，cluster.averages$RNA$`Memory CD4 T`
head(cluster.averages[["RNA"]][, 1:5])

#将计算结果返回seurat用于下步骤计算【将空格替换成_,不然ggplot2会fail】
orig.levels <- levels(pbmc)#clusters
Idents(pbmc) <- gsub(pattern = " ", replacement = "_", x = Idents(pbmc))
orig.levels <- gsub(pattern = " ", replacement = "_", x = orig.levels)
levels(pbmc) <- orig.levels
##levels 和Ident的cluster名称里的空格都要改
cluster.averages <- AverageExpression(pbmc, return.seurat = TRUE)#计算结果存在 cluster.averages[['RNA']]@data 中, 是matrix
cluster.averages
CellScatter(cluster.averages, cell1 = "NK", cell2 = "CD8_T")
#添加样本作为细胞身份
cluster.averages <- AverageExpression(pbmc, return.seurat = TRUE, add.ident = "replicate")
CellScatter(cluster.averages, cell1 = "CD8_T_rep1", cell2 = "CD8_T_rep2")
DoHeatmap(cluster.averages, features = unlist(TopFeatures(pbmc[["pca"]], balanced = TRUE)), size = 3, draw.lines = FALSE)

对于其中TopFeatures() 不了解，查了下

#TopFeatures中balanced的参数，看了帮助文档，没明白。试一下。 
#TopFeatures适用对象就是降维后的对象，单细胞测序中pbmc[['pca']]可用,试了下pbmc[["umap"]]，报错。溢出边界。
> unlist(TopFeatures(pbmc[["pca"]], balanced = TRUE))
positive1  positive2  positive3  positive4  positive5  positive6  positive7 
"CST3"   "TYROBP"     "LST1"     "AIF1"      "FTL"     "FTH1"      "LYZ" 
positive8  positive9 positive10  negative1  negative2  negative3  negative4 
"FCN1"   "S100A9"     "TYMP"   "MALAT1"      "LTB"     "IL32"     "IL7R" 
negative5  negative6  negative7  negative8  negative9 negative10 
"CD2"      "B2M"    "ACAP1"     "CD27"   "STK17A"     "CTSW" 
> unlist(TopFeatures(pbmc[["pca"]], balanced = FALSE))
[1] "MALAT1"   "SPI1"     "IFITM3"   "SERPINA1" "LGALS2"   "CTSS"    
[7] "S100A8"   "LGALS1"   "CFD"      "FCER1G"   "TYMP"     "S100A9"  
[13] "FCN1"     "LYZ"      "FTH1"     "FTL"      "AIF1"     "LST1"    
[19] "TYROBP"   "CST3"

参考教程
https://satijalab.org/seurat/archive/v3.1/interaction_vignette.html
https://www.jianshu.com/p/1a609ee4caef