送你个对象
今天是生信星球陪你的第476天
大神一句话,菜鸟跑半年。我不是大神,但我可以缩短你走弯路的半年~
就像歌儿唱的那样,如果你不知道该往哪儿走,就留在这学点生信好不好~
这里有豆豆和花花的学习历程,从新手到进阶,生信路上有你有我!
豆豆写于19.10.26-27
复杂的分析往往需要大量的数据,在R语言中,理解这些东西就需要~对象!
分析单细胞数据,常见的一个对象就是SingleCellExperiment
或者sce
,那么这次就来认识一下这个基础知识点
这是单细胞分析中的非常常用的S4对象,里面包罗万象,但依然有据可循。那么它是如何组织的?存储了什么内容?这就是我们这次要探索的任务。内容来自:https://osca.bioconductor.org/data-infrastructure.html
首先要上一张图
这个对象一定要不断反复查看

图中最核心的部分,是蓝色的data部分;另外还有绿色的基因注释信息feature metadata、橙色的细胞注释信息cell metadata。除了这三大件,还会包含一些下游分析结果,比如PCA、tSNE降维结果就会保存在紫色的Dimension Reductions
开始发挥想象力了哈!
如果把这个sce
对象比作一艘货船,上面就会装载许多集装箱slots(理解为对象接口)
。每个箱子slot
都是独立的,箱子里面包含的东西不一样,比如有的装食物,有的装砖头,各有各的特征。就sce
来说,有的接口必须是数值型矩阵结构,有的就需要数据框结构。
核心部分-assays
创建一个sce
对象很容易,只需要一个assays
就可以,这是一个列表,其中包含了许多表达数据,例如原始数据count或者其他标准化处理过的数据,行是基因,列是样本
可以构建一个10个基因,3个细胞的矩阵【一定要是矩阵】
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中都可以 -
counts(sce)
:它实现的东西和上面一样,只不过这个方法只适合提取counts
这个名称的矩阵
counts(sce)
# 或者assay(sce, "counts")
## 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还能扩展
既然有了核心,那么就可以根据它进行多向拓展,这也是它强大的一个原因。
使用标准函数扩展
之前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
也是能有两种提取方法:
logcounts(sce)
# assay(sce, "logcounts")
## 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
通过对比这个logcounts
和counts
数据,就能发现为什么要做normalization这一步:原始矩阵中1、2表达量差不多,但和3差别很大,很有可能是细胞3本身测序深度就比较高,因此得到的reads数也多;进行归一化以后,应该就去除了样本文库差异,结果看到1、2、3之间也变得可比了
到目前为止,我们总共得到了:
assays(sce)
## List of length 2
## names(2): counts logcounts
自定义扩展assays
这里的自定义指的是,我们不使用某个R包的某个函数,而是根据自己的想法,去根据原始矩阵得到一个新的矩阵
# 例如自己创建一个新的counts_100矩阵,然后依旧是通过这个名称进行访问
counts_100 <- assay(sce, "counts") + 100
assay(sce, "counts_100") <- counts_100
看一下结果【注意:新增用的是assay
单数,查看结果用的是assays
复数】
assays(sce)
## List of length 3
## names(3): counts logcounts counts_100
小结
再回到第一张图,看看assays
那里,是不是画了深蓝和浅蓝?这也是为了更好地表达:assays
可以包含多个矩阵;构建sce对象可以一次一次加入新的矩阵,也可以用列表的形式,一次加入多个矩阵
列的注释信息:colData
之前有了”核心“——表达矩阵信息,那么其次重要的就是添加注释信息,这部分来介绍列的注释,针对的就是实验样本、细胞。这部分信息将会保存在colData
中,它的主体是样本,于是将行名设定为样本,列名设为注释信息(如:批次、作者等等),对应上面图中的橙色部分。
接下来先设置一个细胞批次注释信息:
cell_metadata <- data.frame(batch = c(1, 1, 2))
rownames(cell_metadata) <- paste0("cell_", 1:3)
有了注释信息,怎么把它加入sce对象呢?
两种方法:一种是直接构建,一种是后续添加
-
直接构建:
sce <- SingleCellExperiment(assays = list(counts = counts_matrix),
colData = cell_metadata) -
后续添加
colData(sce) <- DataFrame(cell_metadata)
然后看看sce对象添加了什么:
可以看到colData
增加了之前设置的batch
信息
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(1): batch
## reducedDimNames(0):
## spikeNames(0):
## altExpNames(0):
加入了sce对象以后,怎么获取它呢?
colData(sce)
## DataFrame with 3 rows and 1 column
## batch
## <numeric>
## cell_1 1
## cell_2 1
## cell_3 2
这个注释信息只能自己手动添加吗?
答案是no!有一些包可以自己计算,并且帮你添加进去。例如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$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):
这样看的colnames
中就只剩两个细胞了
行的注释信息: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):
看到rownames
结果就剩2个基因
除了以上3大块,还有一些重要的”集装箱“需要介绍,这些前面用
附
标识
附1:对样本进行归一化:“sizeFactors“
这里面装了根据样本文库计算的文库大小因子,是一个数值型向量,用于后面的归一化
sce <- scran::computeSumFactors(sce)
sce <- scater::normalize(sce)
sizeFactors(sce)
## [1] 0.503 0.623 1.874
前面提到的:assays
(primary data), colData
(sample metadata), rowData
/rowRanges
(feature metadata), and sizeFactors
。其实这其中前三个都来自于SummarizedExperiment
这个对象。基于这个对象,还建立了一些新的对象接口,例如下面👇的:
附2:降维结果: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
看下全部的结果都包含什么:
reducedDims(sce)
## List of length 2
## names(2): PCA TSNE
调取方法也十分类似assay
,数据矩阵存储在assays
,而调用是assay
;这里的降维结果存储在reducedDims
,调用是reducedDim
这个降维信息只能自己手动添加吗?
答案也是no!和前面的counts_100
加到assays
的思路一样,我们也可以自己计算,而不用现成的函数,最后加到reducedDims
这个接口。
例如,进行UMAP降维,虽然可以用scater::runUMAP()
,但依然可以自己处理。
比如这里使用uwot
包,不过这个包只能计算,不能添加到sce对象,需要手动添加:
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
最后再来看一下:
reducedDims(sce)
## List of length 3
## names(3): PCA TSNE UMAP_uwot
附3:metadata 接口
虽然前面介绍了许多接口,但还是有很多DIY的,不能直接导入它们,不过我们仍然需要这些信息,于是medata接口诞生。它可以存储任意类型的数据,只要给它一个名字。
例如,我们有几个感兴趣的基因(比如是高变化基因),现在想要把它保存在sce
中,方便以后使用:
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"
附4:spike-in信息
使用isSpike
来添加:
isSpike(sce, "ERCC") <- grepl("^ERCC-", rownames(sce))
# 结果就会在sce中添加:
## spikeNames(1): ERCC
spikeNames(sce)
## [1] "ERCC"
table(isSpike(sce, "ERCC"))
# 就能看存在多少Spike-in
点击底部的“阅读原文”,获得更好的阅读体验哦😻
初学生信,很荣幸带你迈出第一步
🤓生信星球 🌎~ 一个不拽术语、通俗易懂的生信知识平台
请关注“恒诺新知”微信公众号,感谢“R语言“,”数据那些事儿“,”老俊俊的生信笔记“,”冷🈚️思“,“珞珈R”,“生信星球”的支持!