单细胞S4 对象学习
参考链接：送你个对象
在seurat 以及Monocle 逐渐都采用 SingleCellExperiment或者sce 进行存储数据， 这是单细胞分析中的非常常用的S4对象 .
https://osca.bioconductor.org/data-infrastructure.html
整体的sce 对象内部数据结果如图：
image.png
这个对象主要包括rowData/rowRanges， ,assays,colData ，reduceDims ，metaData 等常用的接口。
-
rowData/rowRanges ： 对基因进行注释，比如基因别名/ 基因坐标等等
-
最核心是assays 部分：里面含有基因-细胞表达矩阵,可以包含多个数组，比如TPM,logcount,count,RPKM 等等。
-
colData ： 主要对细胞进行注释，比如说批次信息，有利于后面进行批次校正
-
reduceDims ： 可以包含多种降维后细胞坐标，如PCA,T-sne，U-MAP.行名为细胞名
-
metaData : 可以添加任何其他类型的结果。比如说添加几个标记基因的基因名。添加为list 结构，支持多个向量，比如说三个细胞系的标记基因都存储到metaData 中。
对于存储好的SingleCellExperiment或者sce 对象，如何提取数据呢，我们直接使用上面的接口函数。
比如说提取基因-细胞：assay(sce,1) -- 括号中1代表第一个数组，因为assay 可以保存多个数组。类似比如reduceDim(sce,"PCA") 提取PCA降维坐标。
1.核心部分-assays
创建
sce 中可以存储多个array. 比如 原始数据count或者其他标准化处理过的数据，行是基因，列是样本 。
如下：
counts_matrix <- data.frame(cell_1 = rpois(10, 10),
cell_2 = rpois(10, 10),
cell_3 = rpois(10, 30))
rownames(counts_matrix) <- paste0("gene_", 1:10)
counts_matrix <- as.matrix(counts_matrix)
有了这个，就可以用一个list构建出SingleCellExperiment对象。当然，这个list中可以包括任意个矩阵
sce <- SingleCellExperiment(assays = list(counts = counts_matrix))
>sce
## class: SingleCellExperiment
## dim: 10 3
## metadata(0):
## assays(1): counts
## rownames(10): gene_1 gene_2 ... gene_9 gene_10
## rowData names(0):
## colnames(3): cell_1 cell_2 cell_3
## colData names(0):
## reducedDimNames(0):
## spikeNames(0):
## altExpNames(0):
查看：
assay(sce, "counts")：这个方法是最通用的办法，而且其中的counts可以换成其他的名称，只要是出现在之前的list中都可以 。注意是assay不是assays.
assay(sce, "counts")
# 或者counts(sce)
## cell_1 cell_2 cell_3
## gene_1 7 9 35
## gene_2 7 6 38
## gene_3 10 14 32
## gene_4 7 9 32
## gene_5 19 19 48
## gene_6 8 7 26
## gene_7 10 10 28
## gene_8 4 10 26
## gene_9 10 9 37
## gene_10 6 16 26
添加多个数组（通过标准函数扩充）
之前assays中只有原始表达矩阵，其实还能根据它扩展到归一化矩阵，例如使用一些R包的函数对包装的矩阵进行操作：
sce <- scran::computeSumFactors(sce)
sce <- scater::normalize(sce)
> sce
## class: SingleCellExperiment
## dim: 10 3
## metadata(1): log.exprs.offset
## assays(2): counts logcounts
## rownames(10): gene_1 gene_2 ... gene_9 gene_10
## rowData names(0):
## colnames(3): cell_1 cell_2 cell_3
## colData names(0):
## reducedDimNames(0):
## spikeNames(0):
## altExpNames(0):
这样，assays 就从一个存储原始矩阵的counts ，又扩增了归一化矩阵的logcounts 。同理，这个logcounts也是能有两种提取方法：
assay(sce, "logcounts")
# logcounts(sce)
## cell_1 cell_2 cell_3
## gene_1 3.90 3.95 4.30
## gene_2 3.90 3.41 4.41
## gene_3 4.38 4.55 4.18
## gene_4 3.90 3.95 4.18
## gene_5 5.28 4.98 4.73
## gene_6 4.08 3.61 3.89
## gene_7 4.38 4.09 3.99
## gene_8 3.16 4.09 3.89
## gene_9 4.38 3.95 4.37
## gene_10 3.69 4.74 3.89
添加多个数组（通过自定义扩充）
# 例如自己创建一个新的counts_100矩阵，然后依旧是通过这个名称进行访问
counts_100 <- assay(sce, "counts") + 100
assay(sce, "counts_100") <- counts_100
看一下结果【注意：新增用的是assay单数，查看结果用的是assays复数】 ,于是可以sce里面看到有三个assay 数组。
assays(sce)
## List of length 3
## names(3): counts logcounts counts_100
2.列的注释信息：colData
之前有了”核心“——表达矩阵信息，那么其次重要的就是添加注释信息，这部分来介绍列的注释，针对的就是实验样本、细胞。这部分信息将会保存在colData中，它的主体是样本，于是将行名设定为样本，列名设为注释信息（如：批次、作者等等），对应上面图中的橙色部分。
手动创建colData
cell_metadata <- data.frame(batch = c(1, 1, 2))
rownames(cell_metadata) <- paste0("cell_", 1:3)
## 直接创建
sce <- SingleCellExperiment(assays = list(counts = counts_matrix),
colData = cell_metadata)
## 后续添加
colData(sce) <- DataFrame(cell_metadata)
通过函数直接创建
比如 scater包的calculateQCMetrics()就会帮你计算几十项细胞的质量信息，结果依然是使用colData调用注释结果信息
sce <- scater::calculateQCMetrics(sce)
colData(sce)[, 1:5]
## DataFrame with 3 rows and 5 columns
## batch is_cell_control total_features_by_counts
## <numeric> <logical> <integer>
## cell_1 1 FALSE 10
## cell_2 1 FALSE 10
## cell_3 2 FALSE 10
## log10_total_features_by_counts total_counts
## <numeric> <integer>
## cell_1 1.04139268515822 88
## cell_2 1.04139268515822 109
## cell_3 1.04139268515822 328
查询所有的colData
colData(sce)
## DataFrame with 3 rows and 1 column
## batch
## <numeric>
## cell_1 1
## cell_2 1
## cell_3 2
查询部分的colData (某列/某行)
操作和数据框一样。
### 提取某列
sce$batch
# 或者colData(sce)$batch
## [1] 1 1 2
### 提取某行
sce[, sce$batch == 1]
## class: SingleCellExperiment
## dim: 10 2
## metadata(0):
## assays(1): counts
## rownames(10): gene_1 gene_2 ... gene_9 gene_10
## rowData names(7): is_feature_control mean_counts ... total_counts
## log10_total_counts
## colnames(2): cell_1 cell_2
## colData names(10): batch is_cell_control ...
## pct_counts_in_top_200_features pct_counts_in_top_500_features
## reducedDimNames(0):
## spikeNames(0):
## altExpNames(0):
3.行注释信息：rowData/rowRanges
对基因进行注释，比如说基因的别名及其基因区间. rowData / rowRanges
-
rowData：是一个数据框的结构，它就存储了核心assays矩阵的基因相关信息
它返回的结果就是这样：
rowData(sce)[, 1:3]
## DataFrame with 10 rows and 3 columns
## is_feature_control mean_counts log10_mean_counts
## <logical> <numeric> <numeric>
## gene_1 FALSE 17 1.25527250510331
## gene_2 FALSE 17 1.25527250510331
## gene_3 FALSE 18.6666666666667 1.29373075692248
## gene_4 FALSE 16 1.23044892137827
## gene_5 FALSE 28.6666666666667 1.47226875192525
## gene_6 FALSE 13.6666666666667 1.16633142176653
## gene_7 FALSE 16 1.23044892137827
## gene_8 FALSE 13.3333333333333 1.15634720085992
## gene_9 FALSE 18.6666666666667 1.29373075692248
## gene_10 FALSE 16 1.23044892137827
-
rowRanges：也是基因相关，但是它是GRange对象，存储了基因坐标信息，例如染色体信息、起始终点坐标它返回的结果就是这样：
rowRanges(sce) ## GRangesList object of length 10: ## $gene_1 ## GRanges object with 0 ranges and 0 metadata columns: ## seqnames ranges strand ## <Rle> <IRanges> <Rle> ## ------- ## seqinfo: no sequences
查询部分基因注释信息
同样类似于数据框，可以按位置、名称取子集：
sce[c("gene_1", "gene_4"), ]
# 或者 sce[c(1, 4), ]
## class: SingleCellExperiment
## dim: 2 3
## metadata(0):
## assays(1): counts
## rownames(2): gene_1 gene_4
## rowData names(7): is_feature_control mean_counts ... total_counts
## log10_total_counts
## colnames(3): cell_1 cell_2 cell_3
## colData names(10): batch is_cell_control ...
## pct_counts_in_top_200_features pct_counts_in_top_500_features
## reducedDimNames(0):
## spikeNames(0):
## altExpNames(0):
4. 降维结果：reducedDims
存储了原始矩阵的降维结果，可以通过PCA、tSNE、UMAP等得到，它是一个数值型矩阵的list，行名是原来矩阵的列名（就是细胞、样本），它的列就是各种维度信息
它和assays一样，也可以包含许多降维的结果，例如用scater包计算PCA：
sce <- scater::runPCA(sce)
# 这个算法是利用了sce对象的归一化结果logcounts(sce)
reducedDim(sce, "PCA")
## PC1 PC2
## cell_1 -0.639 -0.553
## cell_2 0.982 -0.123
## cell_3 -0.343 0.677
## attr(,"percentVar")
## [1] 65.7 34.3
除了PCA，tSNE的结果也是存储在这里：
sce <- scater::runTSNE(sce, perplexity = 0.1)
reducedDim(sce, "TSNE")
## [,1] [,2]
## cell_1 -5664 -542
## cell_2 3306 -4642
## cell_3 2359 5184
查看全部的维度类型: 注意是reduceDims 不是reduceDim. 类似assays与assay 差别。
reducedDims(sce)
## List of length 2
## names(2): PCA TSNE
除此之外：还可以添加UMAP降维信息
u <- uwot::umap(t(logcounts(sce)), n_neighbors = 2)
reducedDim(sce, "UMAP_uwot") <- u
reducedDim(sce, "UMAP_uwot")
# [,1] [,2]
## cell_1 0.453 -0.464
## cell_2 -0.115 0.633
## cell_3 -0.339 -0.169
## attr(,"scaled:center")
## [1] 8.69 -2.22
用reduceDims 查看所有的存储的降维类型
reducedDims(sce)
## List of length 3
## names(3): PCA TSNE UMAP_uwot
5.metadata 接口
我们之前几个类型，存储的行和列都是和基因数目及其细胞数目有关系，比如PCA 结果,行名都是细胞。假设我们需要存储一些HVGs基因怎么办呢，所有需要一个通用接口。
创建接口：
my_genes <- c("gene_1", "gene_5")
metadata(sce) <- list(favorite_genes = my_genes)
metadata(sce)
## $favorite_genes
## [1] "gene_1" "gene_5"s
扩充接口：
your_genes <- c("gene_4", "gene_8")
metadata(sce)$your_genes <- your_genes
metadata(sce)
## $favorite_genes
## [1] "gene_1" "gene_5"
##
## $your_genes
## [1] "gene_4" "gene_8"
6.spike-in 信息
创建：
使用isSpike来添加：
isSpike(sce, "ERCC") <- grepl("^ERCC-", rownames(sce))
# 结果就会在sce中添加：
## spikeNames(1): ERCC
spikeNames(sce)
## [1] "ERCC"
table(isSpike(sce, "ERCC"))
# 就能看存在多少Spike-in
附：对样本进行归一化：sizeFactors 接口
这里面装了根据样本文库计算的文库大小因子，是一个数值型向量，用于后面的归一化
sce <- scran::computeSumFactors(sce)
sce <- scater::normalize(sce)
sizeFactors(sce)
## [1] 0.503 0.623 1.874
