參考鏈接：送你個(gè)對(duì)象

在seurat 以及Monocle 逐漸都采用 SingleCellExperiment或者sce 進(jìn)行存儲(chǔ)數(shù)據(jù)遭贸， 這是單細(xì)胞分析中的非常常用的S4對(duì)象 .

https://osca.bioconductor.org/data-infrastructure.html

整體的sce 對(duì)象內(nèi)部數(shù)據(jù)結(jié)果如圖：

image.png

這個(gè)對(duì)象主要包括rowData/rowRanges制轰， ,assays,colData 意述，reduceDims 幸冻，metaData 等常用的接口。

• rowData/rowRanges ： 對(duì)基因進(jìn)行注釋舔株，比如基因別名/ 基因坐標(biāo)等等

• 最核心是assays 部分：里面含有基因-細(xì)胞表達(dá)矩陣,可以包含多個(gè)數(shù)組首启，比如TPM,logcount,count,RPKM 等等。

• colData ： 主要對(duì)細(xì)胞進(jìn)行注釋驻龟，比如說批次信息温眉，有利于后面進(jìn)行批次校正

• reduceDims ： 可以包含多種降維后細(xì)胞坐標(biāo)，如PCA,T-sne迅脐，U-MAP.行名為細(xì)胞名

• metaData : 可以添加任何其他類型的結(jié)果芍殖。比如說添加幾個(gè)標(biāo)記基因的基因名。添加為list 結(jié)構(gòu)谴蔑，支持多個(gè)向量豌骏，比如說三個(gè)細(xì)胞系的標(biāo)記基因都存儲(chǔ)到metaData 中。

? 對(duì)于存儲(chǔ)好的SingleCellExperiment或者sce 對(duì)象隐锭，如何提取數(shù)據(jù)呢窃躲，我們直接使用上面的接口函數(shù)。

? 比如說提取基因-細(xì)胞：assay(sce,1) -- 括號(hào)中1代表第一個(gè)數(shù)組钦睡，因?yàn)閍ssay 可以保存多個(gè)數(shù)組蒂窒。類似比如reduceDim(sce,"PCA") 提取PCA降維坐標(biāo)。

1.核心部分-assays

創(chuàng)建

sce 中可以存儲(chǔ)多個(gè)array. 比如 原始數(shù)據(jù)count或者其他標(biāo)準(zhǔn)化處理過的數(shù)據(jù)荞怒，行是基因洒琢，列是樣本 。

如下：

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) 

有了這個(gè)褐桌，就可以用一個(gè)list構(gòu)建出SingleCellExperiment對(duì)象衰抑。當(dāng)然，這個(gè)list中可以包括任意個(gè)矩陣

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")：這個(gè)方法是最通用的辦法荧嵌，而且其中的counts可以換成其他的名稱呛踊，只要是出現(xiàn)在之前的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

添加多個(gè)數(shù)組（通過標(biāo)準(zhǔn)函數(shù)擴(kuò)充）

之前assays中只有原始表達(dá)矩陣啦撮，其實(shí)還能根據(jù)它擴(kuò)展到歸一化矩陣谭网，例如使用一些R包的函數(shù)對(duì)包裝的矩陣進(jìn)行操作：

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 就從一個(gè)存儲(chǔ)原始矩陣的counts 赃春，又?jǐn)U增了歸一化矩陣的logcounts 愉择。同理，這個(gè)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

添加多個(gè)數(shù)組（通過自定義擴(kuò)充）

# 例如自己創(chuàng)建一個(gè)新的counts_100矩陣，然后依舊是通過這個(gè)名稱進(jìn)行訪問
counts_100 <- assay(sce, "counts") + 100
assay(sce, "counts_100") <- counts_100  

看一下結(jié)果【注意：新增用的是assay單數(shù)薄辅，查看結(jié)果用的是assays復(fù)數(shù)】 ,于是可以sce里面看到有三個(gè)assay 數(shù)組要拂。

assays(sce)
## List of length 3
## names(3): counts logcounts counts_100

2.列的注釋信息：colData

之前有了”核心“——表達(dá)矩陣信息，那么其次重要的就是添加注釋信息站楚，這部分來介紹列的注釋脱惰，針對(duì)的就是實(shí)驗(yàn)樣本、細(xì)胞窿春。這部分信息將會(huì)保存在colData中拉一，它的主體是樣本，于是將行名設(shè)定為樣本旧乞，列名設(shè)為注釋信息（如：批次蔚润、作者等等），對(duì)應(yīng)上面圖中的橙色部分尺栖。

手動(dòng)創(chuàng)建colData

cell_metadata <- data.frame(batch = c(1, 1, 2))
rownames(cell_metadata) <- paste0("cell_", 1:3)

## 直接創(chuàng)建
sce <- SingleCellExperiment(assays = list(counts = counts_matrix),
                         colData = cell_metadata)
## 后續(xù)添加
colData(sce) <- DataFrame(cell_metadata)                   

通過函數(shù)直接創(chuàng)建

比如 scater包的calculateQCMetrics()就會(huì)幫你計(jì)算幾十項(xiàng)細(xì)胞的質(zhì)量信息嫡纠，結(jié)果依然是使用colData調(diào)用注釋結(jié)果信息

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 (某列/某行)

操作和數(shù)據(jù)框一樣。

### 提取某列
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

對(duì)基因進(jìn)行注釋延赌，比如說基因的別名及其基因區(qū)間. rowData / rowRanges

• rowData：是一個(gè)數(shù)據(jù)框的結(jié)構(gòu)除盏，它就存儲(chǔ)了核心assays矩陣的基因相關(guān)信息

它返回的結(jié)果就是這樣：

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：也是基因相關(guān)，但是它是GRange對(duì)象挫以，存儲(chǔ)了基因坐標(biāo)信息者蠕，例如染色體信息、起始終點(diǎn)坐標(biāo)

  它返回的結(jié)果就是這樣：

  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
  

查詢部分基因注釋信息

同樣類似于數(shù)據(jù)框掐松，可以按位置踱侣、名稱取子集：

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. 降維結(jié)果：reducedDims

存儲(chǔ)了原始矩陣的降維結(jié)果，可以通過PCA大磺、tSNE抡句、UMAP等得到，它是一個(gè)數(shù)值型矩陣的list杠愧，行名是原來矩陣的列名（就是細(xì)胞待榔、樣本），它的列就是各種維度信息

它和assays一樣殴蹄，也可以包含許多降維的結(jié)果，例如用scater包計(jì)算PCA：

sce <- scater::runPCA(sce)
# 這個(gè)算法是利用了sce對(duì)象的歸一化結(jié)果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的結(jié)果也是存儲(chǔ)在這里：

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 查看所有的存儲(chǔ)的降維類型

reducedDims(sce)
## List of length 3
## names(3): PCA TSNE UMAP_uwot

5.metadata 接口

我們之前幾個(gè)類型，存儲(chǔ)的行和列都是和基因數(shù)目及其細(xì)胞數(shù)目有關(guān)系绑嘹，比如PCA 結(jié)果,行名都是細(xì)胞稽荧。假設(shè)我們需要存儲(chǔ)一些HVGs基因怎么辦呢，所有需要一個(gè)通用接口工腋。

創(chuàng)建接口：

my_genes <- c("gene_1", "gene_5")
metadata(sce) <- list(favorite_genes = my_genes)

metadata(sce)
## $favorite_genes
## [1] "gene_1" "gene_5"s

擴(kuò)充接口：

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 信息

創(chuàng)建：

使用isSpike來添加：

isSpike(sce, "ERCC") <- grepl("^ERCC-", rownames(sce))

# 結(jié)果就會(huì)在sce中添加：
## spikeNames(1): ERCC

spikeNames(sce)
## [1] "ERCC"

table(isSpike(sce, "ERCC"))
# 就能看存在多少Spike-in

附：對(duì)樣本進(jìn)行歸一化：sizeFactors 接口

這里面裝了根據(jù)樣本文庫計(jì)算的文庫大小因子姨丈，是一個(gè)數(shù)值型向量畅卓，用于后面的歸一化

sce <- scran::computeSumFactors(sce)
sce <- scater::normalize(sce)
sizeFactors(sce)
## [1] 0.503 0.623 1.874