# 1.1 数据结构 > 分析单细胞数据，常见的一个名称就是`SingleCellExperiment`或者`sce` \ > 那么这次就来认识一下这个基础知识点 ## 前言 我们首先要对单细胞分析的流程有一个大概的认识： ### **上半场分析：** 这个差不多属于固定的流程了 * 首先是数据导入 * 紧接着进行质量控制和归一化，这个步骤包含了矫正实验因素以及其他影响因素，目的是从raw count得到clean count * 有了clean count，就要挑取导致表达量差异的基因们（用于后面的降维） 因为降维的过程其实就是去繁存简，每个基因的变化都对整体有影响，对细胞来讲都是一个变化维度。但这么多维度我们看不了也处理不了，需要尽可能保证真实差异的前提下减少维度的数量，因此需要挑出那些更能代表整体差异的基因 * 降维 ![](/files/-MCvnuUDDdWc9eNeCarp) ### **下半场分析：** 这个就可以分出很多分支，例如 * 分群（clustering）：意在探索如何把scRNA数据集给拆分掉，变成一小块一小块的数据，这每一小块基因变化都是相似的 * 差异分析（differential expression）：不同组的细胞之间表达量差异是如何产生的 * 数据集的整合(Integrating datasets)：scRNA数据集越来越多，数据集之间的比较也日益显著 * 处理大型数据（large scale data）：这部分仅仅依靠内存存储的数据是不够的，还有更好的办法 * 还有：轨迹推断、细胞周期推断等等特定需求 > 知道了大体的流程，那么就先看看必须了解的数据结构 这是单细胞分析中的非常常用的S4对象，里面包罗万象，但依然有据可循。**那么它是如何组织的？存储了什么内容？**这就是我们这次要探索的任务。 ## 首先要上一张图 > 这张图会不断反复查看 ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-23-131252.png) 图中最核心的部分，是蓝色的data部分；另外还有绿色的基因注释信息feature metadata、橙色的细胞注释信息cell metadata。除了这三大件，还会包含一些下游分析结果，比如PCA、tSNE降维结果就会保存在紫色的Dimension Reductions 这个`SingleCellExperiment`对象来自`SingleCellExperiment`R包，现在Bioconductor上的70多个单细胞相关的R包都使用了这个对象，**可以说是单细胞领域的通用货币**，包含了： gene-by-cell expression data、 per-cell metadata 、 per-gene annotation > **开始发挥想象力了哈！** 如果把这个`SingleCellExperiment`对象比作一艘货船，上面就会装载许多集装箱`slots（理解为对象结构）` 。每个箱子`slot`都是独立的，箱子里面包含的东西不一样，比如有的装食物，有的装砖头，各有各的特征。就`sce`来说，有的接口必须是数值型矩阵结构，有的就需要数据框结构。 ## 核心部分-`assays` **创建一个`sce`对象很容易，只需要一个`assays`就可以**，这是一个列表，其中包含了许多表达数据，例如原始数据count或者其他标准化处理过的数据，行是基因，列是样本 可以构建一个10个基因，3个细胞的矩阵【一定要是矩阵】 ```r 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中可以包括任意个矩阵** ```r # BiocManager::install('SingleCellExperiment') 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`这个名称的矩阵 ```r 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包的函数对包装的矩阵进行操作： ```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`也是能有两种提取方法： ```r 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之间也变得可比了 到目前为止，我们总共得到了： ```r assays(sce) ## List of length 2 ## names(2): counts logcounts ``` ### **自定义扩展assays** 这里的自定义指的是，我们不使用某个R包的某个函数，而是根据自己的想法，去根据原始矩阵得到一个新的矩阵 ```r # 例如自己创建一个新的counts_100矩阵，然后依旧是通过这个名称进行访问 counts_100 <- assay(sce, "counts") + 100 assay(sce, "counts_100") <- counts_100 ``` 看一下结果【注意：新增用的是`assay`单数，查看结果用的是`assays`复数】 ```r assays(sce) ## List of length 3 ## names(3): counts logcounts counts_100 ``` ### **小结** **再回到第一张图，看看`assays`那里，是不是画了深蓝和浅蓝？**这也是为了更好地表达：`assays`可以包含多个矩阵；构建sce对象可以一次一次加入新的矩阵，也可以用列表的形式，一次加入多个矩阵 ## 列的注释信息：`colData` 之前有了”核心“——表达矩阵信息，那么其次重要的就是添加注释信息，这部分来介绍列的注释，针对的就是实验样本、细胞。这部分信息将会保存在`colData`中，它的主体是样本，于是将行名设定为样本，列名设为注释信息（如：批次、作者等等），对应上面图中的橙色部分。 ### **接下来先设置一个细胞批次注释信息：** ```r cell_metadata <- data.frame(batch = c(1, 1, 2)) rownames(cell_metadata) <- paste0("cell_", 1:3) ``` ### **有了注释信息，怎么把它加入sce对象呢？** 两种方法：一种是直接构建，一种是后续添加 * 直接构建： ```r sce <- SingleCellExperiment(assays = list(counts = counts_matrix), colData = cell_metadata) ``` * 后续添加 ```r colData(sce) <- DataFrame(cell_metadata) ``` 然后看看sce对象添加了什么： 可以看到`colData`增加了之前设置的`batch`信息 ```r 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对象以后，怎么获取它呢？** ```r colData(sce) ## DataFrame with 3 rows and 1 column ## batch ## ## cell_1 1 ## cell_2 1 ## cell_3 2 ``` 或者直接看结果信息： ```r sce$batch ## [1] 1 1 2 ``` ### **这个注释信息只能自己手动添加吗？** 答案是no！有一些包可以自己计算，并且帮你添加进去。例如`scater`包的`calculateQCMetrics()`就会帮你计算几十项细胞的质量信息，结果依然是使用`colData`调用注释结果信息 ```r sce <- scater::addPerCellQC(sce) colData(sce)[, 1:5] ## DataFrame with 3 rows and 5 columns ## batch sum detected percent_top_50 percent_top_100 ## ## cell_1 1 80 10 100 100 ## cell_2 1 88 10 100 100 ## cell_3 2 309 10 100 100 ``` 还能任意添加 ```r sce$more_stuff <- runif(ncol(sce)) colnames(colData(sce)) sce$more_stuff ## 0.04094262 0.41983904 0.51959352 ``` ### **既然colData可以包含这么多的注释信息，那么怎么从中选取一部分呢？** `colData`的一个常用操作就是取子集，看下面操作： 例如想从`colData`中选择批次信息，和数据框取子集是类似的 ```r sce$batch # 或者colData(sce)$batch ## [1] 1 1 2 ``` 然后如果再想取批次中属于第一个批次的信息，就可以继续按照数据框的思路取子集： ```r 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(0): ## colnames(2): cell_1 cell_2 ## colData names(9): batch sum ... total more_stuff ## reducedDimNames(0): ## altExpNames(0): ``` 这样看的`colnames`中就只剩两个细胞了 ## 行的注释信息：`rowData/rowRanges` 既然样本有注释信息，那么同样的，基因也有自己的注释，它就存放在`rowData`或者`rowRanges`中，这两个的区别就是： * `rowData`：是一个数据框的结构，它就存储了核心`assays`矩阵的基因相关信息 它返回的结果就是这样： ```r # 一开始rowData(sce)是空的，可以添加 sce <- scater::addPerFeatureQC(sce) rowData(sce) ## DataFrame with 10 rows and 2 columns ## mean detected ## ## gene_1 16.0000 100 ## gene_2 14.3333 100 ## gene_3 16.0000 100 ## gene_4 18.6667 100 ## gene_5 15.3333 100 ## gene_6 16.6667 100 ## gene_7 17.6667 100 ## gene_8 13.0000 100 ## gene_9 16.6667 100 ## gene_10 14.6667 100 ``` * `rowRanges`：也是基因相关，但是它是GRange对象，存储了基因坐标信息，例如染色体信息、起始终点坐标 它返回的结果一开始是空的： ```r rowRanges(sce) ## GRangesList object of length 10: ## $gene_1 ## GRanges object with 0 ranges and 0 metadata columns: ## seqnames ranges strand ## ## ------- ## seqinfo: no sequences ## ## $gene_2 ## GRanges object with 0 ranges and 0 metadata columns: ## seqnames ranges strand ## ## ------- ## seqinfo: no sequences ## ## $gene_3 ## GRanges object with 0 ranges and 0 metadata columns: ## seqnames ranges strand ## ## ------- ## seqinfo: no sequences ## ## ... ## <7 more elements> ``` ### **怎么按行取子集？** 同样类似于数据框，可以按位置、名称取子集： ```r 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** 这里面装了根据样本文库计算的文库大小因子，是一个数值型向量，用于后面的归一化 ```r sce <- scran::computeSumFactors(sce) sizeFactors(sce) ## [1] 0.5031447 0.5534591 1.9433962 # 或者手动添加 sizeFactors(sce) <- scater::librarySizeFactors(sce) sizeFactors(sce) ## cell_1 cell_2 cell_3 ## 0.5031447 0.5534591 1.9433962 ``` 前面提到的： `assays` (primary data), `colData` (sample metadata), `rowData`/`rowRanges` (feature metadata), and `sizeFactors` 。其实这其中前三个都来自于`SummarizedExperiment`这个对象。基于这个对象，还建立了一些新的对象接口，例如下面👇的： ### **附2：降维结果：reducedDims** 存储了原始矩阵的降维结果，可以通过PCA、tSNE、UMAP等得到，它是一个数值型矩阵的list，行名是原来矩阵的列名（就是细胞、样本），它的列就是各种维度信息 它和`assays`一样，也可以包含许多降维的结果，例如用`scater`包计算PCA： ```r sce <- scater::logNormCounts(sce) sce <- scater::runPCA(sce) # 这个算法是利用了sce对象的归一化结果logcounts(sce) reducedDim(sce, "PCA") # PC1 PC2 # cell_1 -0.6897959 -0.4080860 # cell_2 0.8493039 -0.2145128 # cell_3 -0.1595080 0.6225988 # attr(,"percentVar") # [1] 67.07301 32.92699 ``` **除了PCA，tSNE的结果也是存储在这里：** ```r sce <- scater::runTSNE(sce, perplexity = 0.1) reducedDim(sce, "TSNE") ## [,1] [,2] ## cell_1 35.64851 5694.164 ## cell_2 -4951.13619 -2815.978 ## cell_3 4915.48768 -2878.186 ``` **看下全部的结果都包含什么：** ```r 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对象，需要手动添加： ```r u <- uwot::umap(t(logcounts(sce)), n_neighbors = 2) reducedDim(sce, "UMAP_uwot") <- u ## [,1] [,2] ## cell_1 -0.3952368 -0.03182602 ## cell_2 0.2347576 0.60432243 ## cell_3 0.1604792 -0.57249641 ## attr(,"scaled:center") ## [1] 7.596012 19.251450 ``` 最后再来看一下： ```r reducedDims(sce) ## List of length 3 ## names(3): PCA TSNE UMAP_uwot ``` ### **附3：metadata 接口** 虽然前面介绍了许多接口，但还是有很多DIY的，不能直接导入它们，不过我们仍然需要这些信息，于是medata接口诞生。它可以存储任意类型的数据，只要给它一个名字。 例如，我们有几个感兴趣的基因（比如是高变化基因），现在想要把它保存在`sce`中，方便以后使用： ```r my_genes <- c("gene_1", "gene_5") metadata(sce) <- list(favorite_genes = my_genes) metadata(sce) ## $favorite_genes ## [1] "gene_1" "gene_5"s ``` 既然是一个列表，就意味着支持扩展： ```r 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' is deprecated`： ```r isSpike(sce, "ERCC") <- grepl("^ERCC-", rownames(sce)) # 结果就会在sce中添加： ## spikeNames(1): ERCC spikeNames(sce) ## [1] "ERCC" table(isSpike(sce, "ERCC")) # 就能看存在多少Spike-in ``` **现在推出了：`alternative Experiments`** ```r spike_counts <- cbind(cell_1 = rpois(5, 10), cell_2 = rpois(5, 10), cell_3 = rpois(5, 30)) rownames(spike_counts) <- paste0("spike_", 1:5) spike_se <- SummarizedExperiment(list(counts=spike_counts)) spike_se ## class: SummarizedExperiment ## dim: 5 3 ## metadata(0): ## assays(1): counts ## rownames(5): spike_1 spike_2 spike_3 spike_4 spike_5 ## rowData names(0): ## colnames(3): cell_1 cell_2 cell_3 ## colData names(0): ``` 然后通过`altExp()`加入进来，它的操作和`assays`或`reducedDims()`类似： ```r altExp(sce, "spike") <- spike_se altExps(sce) ## List of length 1 ## names(1): spike # 提取数据也是类似的 sub <- sce[,1:2] # retain only two samples. altExp(sub, "spike") ## class: SummarizedExperiment ## dim: 5 2 ## metadata(0): ## assays(1): counts ## rownames(5): spike_1 spike_2 spike_3 spike_4 spike_5 ## rowData names(0): ## colnames(2): cell_1 cell_2 ## colData names(0): ``` ### **附5：对列添加label** 使用`colLabels()` ，尤其在非监督聚类过程中对细胞添加label，进行分组 ```r colLabels(sce) <- LETTERS[1:3] colLabels(sce) ## [1] "A" "B" "C" ``` 例如`scran::findMarkers`就是通过`colLabels()`来提取细胞信息的 #### 小结 `SingleCellExperiment`对象兼容性很好，可以用于多个scRNA的R包的输入数据 中间分析的每一步都会向这个对象的`assays`, `colData`, `reducedDims` 等模块添加信息 这个对象方便了未来数据传输与协作，这本书的后续也会基于这个对象进行探讨 --- # Agent Instructions: Querying This Documentation If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question. Perform an HTTP GET request on the current page URL with the `ask` query parameter: ``` GET https://jieandze1314.osca.top/01/01-1.md?ask= ``` The question should be specific, self-contained, and written in natural language. The response will contain a direct answer to the question and relevant excerpts and sources from the documentation. Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.