劉小澤寫于19.10.26-27，更新于2020-06-23，最近會(huì)繼續(xù)整理這本資料
為何取名叫“交響樂(lè)”毯辅？因?yàn)閱渭?xì)胞分析就像一個(gè)大樂(lè)團(tuán)，需要各個(gè)流程的協(xié)同配合

分析單細(xì)胞數(shù)據(jù)煞额，常見(jiàn)的一個(gè)名稱就是SingleCellExperiment或者sce思恐，那么這次就來(lái)認(rèn)識(shí)一下這個(gè)基礎(chǔ)知識(shí)點(diǎn)

前言

我們首先要對(duì)單細(xì)胞分析的流程有一個(gè)大概的認(rèn)識(shí)：

上半場(chǎng)分析：

這個(gè)差不多屬于固定的流程了

• 首先是數(shù)據(jù)導(dǎo)入

• 緊接著進(jìn)行質(zhì)量控制和歸一化，這個(gè)步驟包含了矯正實(shí)驗(yàn)因素以及其他影響因素膊毁，目的是從raw count得到clean count

• 有了clean count胀莹，就要挑取導(dǎo)致表達(dá)量差異的基因們（用于后面的降維）
  因?yàn)榻稻S的過(guò)程其實(shí)就是去繁存簡(jiǎn)，每個(gè)基因的變化都對(duì)整體有影響婚温，對(duì)細(xì)胞來(lái)講都是一個(gè)變化維度描焰。但這么多維度我們看不了也處理不了，需要盡可能保證真實(shí)差異的前提下減少維度的數(shù)量，因此需要挑出那些更能代表整體差異的基因

  
  

• 降維

下半場(chǎng)分析：

這個(gè)就可以分出很多分支荆秦，例如

• 分群（clustering）：意在探索如何把scRNA數(shù)據(jù)集給拆分掉篱竭，變成一小塊一小塊的數(shù)據(jù)膜毁，這每一小塊基因變化都是相似的
• 差異分析（differential expression）：不同組的細(xì)胞之間表達(dá)量差異是如何產(chǎn)生的
• 數(shù)據(jù)集的整合(Integrating datasets)：scRNA數(shù)據(jù)集越來(lái)越多莱衩，數(shù)據(jù)集之間的比較也日益顯著
• 處理大型數(shù)據(jù)（large scale data）：這部分僅僅依靠?jī)?nèi)存存儲(chǔ)的數(shù)據(jù)是不夠的，還有更好的辦法
• 還有：軌跡推斷衡奥、細(xì)胞周期推斷等等特定需求

知道了大體的流程瓤介，那么就先看看必須了解的數(shù)據(jù)結(jié)構(gòu)

這是單細(xì)胞分析中的非常常用的S4對(duì)象吕喘，里面包羅萬(wàn)象，但依然有據(jù)可循刑桑。那么它是如何組織的氯质？存儲(chǔ)了什么內(nèi)容？這就是我們這次要探索的任務(wù)漾月。內(nèi)容來(lái)自：https://osca.bioconductor.org/data-infrastructure.html

首先要上一張圖

這張圖會(huì)不斷反復(fù)查看

圖中最核心的部分病梢，是藍(lán)色的data部分；另外還有綠色的基因注釋信息feature metadata梁肿、橙色的細(xì)胞注釋信息cell metadata蜓陌。除了這三大件，還會(huì)包含一些下游分析結(jié)果吩蔑，比如PCA钮热、tSNE降維結(jié)果就會(huì)保存在紫色的Dimension Reductions

這個(gè)SingleCellExperiment對(duì)象來(lái)自SingleCellExperimentR包，現(xiàn)在Bioconductor上的70多個(gè)單細(xì)胞相關(guān)的R包都使用了這個(gè)對(duì)象烛芬，可以說(shuō)是單細(xì)胞領(lǐng)域的通用貨幣隧期，包含了： gene-by-cell expression data、 per-cell metadata 赘娄、 per-gene annotation

開(kāi)始發(fā)揮想象力了哈仆潮！

如果把這個(gè)SingleCellExperiment對(duì)象比作一艘貨船，上面就會(huì)裝載許多集裝箱slots（理解為對(duì)象結(jié)構(gòu)） 遣臼。每個(gè)箱子slot都是獨(dú)立的性置，箱子里面包含的東西不一樣，比如有的裝食物揍堰，有的裝磚頭鹏浅，各有各的特征。就sce來(lái)說(shuō)屏歹，有的接口必須是數(shù)值型矩陣結(jié)構(gòu)隐砸，有的就需要數(shù)據(jù)框結(jié)構(gòu)。

核心部分-assays

創(chuàng)建一個(gè)sce對(duì)象很容易蝙眶，只需要一個(gè)assays就可以季希，這是一個(gè)列表，其中包含了許多表達(dá)數(shù)據(jù)，例如原始數(shù)據(jù)count或者其他標(biāo)準(zhǔn)化處理過(guò)的數(shù)據(jù)胖眷，行是基因武通，列是樣本

可以構(gòu)建一個(gè)10個(gè)基因，3個(gè)細(xì)胞的矩陣【一定要是矩陣】

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è)矩陣

# 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):

為了得到這個(gè)對(duì)象中的矩陣境析，可以用兩種方式獲得：

• assay(sce, "counts")：這個(gè)方法是最通用的辦法囚枪，而且其中的counts可以換成其他的名稱，只要是出現(xiàn)在之前的list中都可以
• counts(sce)：它實(shí)現(xiàn)的東西和上面一樣劳淆，只不過(guò)這個(gè)方法只適合提取counts這個(gè)名稱的矩陣

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還能擴(kuò)展

既然有了核心链沼，那么就可以根據(jù)它進(jìn)行多向拓展，這也是它強(qiáng)大的一個(gè)原因沛鸵。

使用標(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也是能有兩種提取方法：

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

通過(guò)對(duì)比這個(gè)logcounts和counts數(shù)據(jù)乱豆，就能發(fā)現(xiàn)為什么要做normalization這一步：原始矩陣中1、2表達(dá)量差不多吊趾，但和3差別很大宛裕，很有可能是細(xì)胞3本身測(cè)序深度就比較高，因此得到的reads數(shù)也多论泛；進(jìn)行歸一化以后揩尸，應(yīng)該就去除了樣本文庫(kù)差異屁奏，結(jié)果看到1岩榆、2、3之間也變得可比了

到目前為止了袁，我們總共得到了：

assays(sce)

## List of length 2
## names(2): counts logcounts

自定義擴(kuò)展assays

這里的自定義指的是朗恳，我們不使用某個(gè)R包的某個(gè)函數(shù)湿颅，而是根據(jù)自己的想法载绿，去根據(jù)原始矩陣得到一個(gè)新的矩陣

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

看一下結(jié)果【注意：新增用的是assay單數(shù)油航，查看結(jié)果用的是assays復(fù)數(shù)】

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

小結(jié)

再回到第一張圖崭庸，看看assays那里，是不是畫了深藍(lán)和淺藍(lán)？這也是為了更好地表達(dá)：assays可以包含多個(gè)矩陣怕享；構(gòu)建sce對(duì)象可以一次一次加入新的矩陣执赡，也可以用列表的形式，一次加入多個(gè)矩陣

列的注釋信息：colData

之前有了”核心“——表達(dá)矩陣信息函筋，那么其次重要的就是添加注釋信息沙合，這部分來(lái)介紹列的注釋，針對(duì)的就是實(shí)驗(yàn)樣本跌帐、細(xì)胞首懈。這部分信息將會(huì)保存在colData中，它的主體是樣本谨敛，于是將行名設(shè)定為樣本究履，列名設(shè)為注釋信息（如：批次、作者等等）脸狸，對(duì)應(yīng)上面圖中的橙色部分最仑。

接下來(lái)先設(shè)置一個(gè)細(xì)胞批次注釋信息：

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

有了注釋信息，怎么把它加入sce對(duì)象呢炊甲？

兩種方法：一種是直接構(gòu)建泥彤，一種是后續(xù)添加

• 直接構(gòu)建：

  sce <- SingleCellExperiment(assays = list(counts = counts_matrix),
                             colData = cell_metadata)
  

• 后續(xù)添加

  colData(sce) <- DataFrame(cell_metadata)
  

然后看看sce對(duì)象添加了什么：

可以看到colData增加了之前設(shè)置的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對(duì)象以后，怎么獲取它呢蜜葱？

colData(sce)
## DataFrame with 3 rows and 1 column
##            batch
##        <numeric>
## cell_1         1
## cell_2         1
## cell_3         2

或者直接看結(jié)果信息：

sce$batch
## [1] 1 1 2

這個(gè)注釋信息只能自己手動(dòng)添加嗎全景？

答案是no！有一些包可以自己計(jì)算牵囤，并且?guī)湍闾砑舆M(jìn)去爸黄。例如scater包的calculateQCMetrics()就會(huì)幫你計(jì)算幾十項(xiàng)細(xì)胞的質(zhì)量信息，結(jié)果依然是使用colData調(diào)用注釋結(jié)果信息

sce <- scater::addPerCellQC(sce)
colData(sce)[, 1:5]
## DataFrame with 3 rows and 5 columns
##            batch       sum  detected percent_top_50 percent_top_100
##        <numeric> <integer> <integer>      <numeric>       <numeric>
## cell_1         1        80        10            100             100
## cell_2         1        88        10            100             100
## cell_3         2       309        10            100             100

還能任意添加

sce$more_stuff <- runif(ncol(sce))
colnames(colData(sce))
sce$more_stuff
## 0.04094262 0.41983904 0.51959352

既然colData可以包含這么多的注釋信息揭鳞，那么怎么從中選取一部分呢炕贵？

colData的一個(gè)常用操作就是取子集，看下面操作：

例如想從colData中選擇批次信息野崇，和數(shù)據(jù)框取子集是類似的

sce$batch
# 或者colData(sce)$batch

## [1] 1 1 2

然后如果再想取批次中屬于第一個(gè)批次的信息称开，就可以繼續(xù)按照數(shù)據(jù)框的思路取子集：

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中就只剩兩個(gè)細(xì)胞了

行的注釋信息：rowData/rowRanges

既然樣本有注釋信息，那么同樣的乓梨，基因也有自己的注釋鳖轰，它就存放在rowData或者rowRanges中，這兩個(gè)的區(qū)別就是：

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

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

  # 一開(kāi)始rowData(sce)是空的蕴侣，可以添加
  sce <- scater::addPerFeatureQC(sce)
  rowData(sce)
  
  ## DataFrame with 10 rows and 2 columns
  ##              mean  detected
  ##         <numeric> <numeric>
  ## 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：也是基因相關(guān)，但是它是GRange對(duì)象臭觉，存儲(chǔ)了基因坐標(biāo)信息昆雀，例如染色體信息辱志、起始終點(diǎn)坐標(biāo)

  它返回的結(jié)果一開(kāi)始是空的：

  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
  ## 
  ## $gene_2
  ## GRanges object with 0 ranges and 0 metadata columns:
  ##    seqnames    ranges strand
  ##       <Rle> <IRanges>  <Rle>
  ##   -------
  ##   seqinfo: no sequences
  ## 
  ## $gene_3
  ## GRanges object with 0 ranges and 0 metadata columns:
  ##    seqnames    ranges strand
  ##       <Rle> <IRanges>  <Rle>
  ##   -------
  ##   seqinfo: no sequences
  ## 
  ## ...
  ## <7 more elements>
  

怎么按行取子集？

同樣類似于數(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):

看到rownames結(jié)果就剩2個(gè)基因

除了以上3大塊，還有一些重要的”集裝箱“需要介紹挽封，這些前面用附標(biāo)識(shí)

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

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

sce <- scran::computeSumFactors(sce)
sizeFactors(sce)
## [1] 0.5031447 0.5534591 1.9433962

# 或者手動(dòng)添加
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 辅愿。其實(shí)這其中前三個(gè)都來(lái)自于SummarizedExperiment這個(gè)對(duì)象和悦。基于這個(gè)對(duì)象渠缕，還建立了一些新的對(duì)象接口鸽素，例如下面??的：

附2：降維結(jié)果：reducedDims

存儲(chǔ)了原始矩陣的降維結(jié)果，可以通過(guò)PCA亦鳞、tSNE馍忽、UMAP等得到，它是一個(gè)數(shù)值型矩陣的list燕差，行名是原來(lái)矩陣的列名（就是細(xì)胞遭笋、樣本），它的列就是各種維度信息

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

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

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

看下全部的結(jié)果都包含什么：

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

調(diào)取方法也十分類似：assay测暗，數(shù)據(jù)矩陣存儲(chǔ)在assays央串，而調(diào)用是assay；這里的降維結(jié)果存儲(chǔ)在reducedDims碗啄，調(diào)用是reducedDim

這個(gè)降維信息只能自己手動(dòng)添加嗎质和？

答案也是no！和前面的counts_100加到assays的思路一樣稚字，我們也可以自己計(jì)算饲宿，而不用現(xiàn)成的函數(shù)，最后加到reducedDims這個(gè)接口胆描。

例如瘫想，進(jìn)行UMAP降維，雖然可以用scater::runUMAP()昌讲，但依然可以自己處理国夜。

比如這里使用uwot包，不過(guò)這個(gè)包只能計(jì)算剧蚣，不能添加到sce對(duì)象支竹，需要手動(dòng)添加：

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

最后再來(lái)看一下：

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

附3：metadata 接口

雖然前面介紹了許多接口，但還是有很多DIY的鸠按，不能直接導(dǎo)入它們礼搁，不過(guò)我們?nèi)匀恍枰@些信息，于是medata接口誕生目尖。它可以存儲(chǔ)任意類型的數(shù)據(jù)馒吴，只要給它一個(gè)名字。

例如瑟曲，我們有幾個(gè)感興趣的基因（比如是高變化基因）饮戳，現(xiàn)在想要把它保存在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

既然是一個(gè)列表洞拨，就意味著支持?jǐn)U展：

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

還是可以繼續(xù)使用isSpike來(lái)添加扯罐，雖然會(huì)顯示'isSpike' is deprecated：

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

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

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

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

現(xiàn)在推出了：alternative Experiments

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):

然后通過(guò)altExp()加入進(jìn)來(lái)，它的操作和assays或reducedDims()類似：

altExp(sce, "spike") <- spike_se
altExps(sce)
## List of length 1
## names(1): spike

# 提取數(shù)據(jù)也是類似的
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：對(duì)列添加label

使用colLabels() 烦衣，尤其在非監(jiān)督聚類過(guò)程中對(duì)細(xì)胞添加label歹河，進(jìn)行分組

colLabels(sce) <- LETTERS[1:3]
colLabels(sce)
## [1] "A" "B" "C"

例如scran::findMarkers就是通過(guò)colLabels()來(lái)提取細(xì)胞信息的

小結(jié)

SingleCellExperiment對(duì)象兼容性很好，可以用于多個(gè)scRNA的R包的輸入數(shù)據(jù)

中間分析的每一步都會(huì)向這個(gè)對(duì)象的assays, colData, reducedDims 等模塊添加信息

這個(gè)對(duì)象方便了未來(lái)數(shù)據(jù)傳輸與協(xié)作花吟，這本書(shū)的后續(xù)也會(huì)基于這個(gè)對(duì)象進(jìn)行探討

歡迎關(guān)注我們的公眾號(hào)～_～　　
我們是兩個(gè)農(nóng)轉(zhuǎn)生信的小碩秸歧，打造生信星球，想讓它成為一個(gè)不拽術(shù)語(yǔ)衅澈、通俗易懂的生信知識(shí)平臺(tái)键菱。需要幫助或提出意見(jiàn)請(qǐng)后臺(tái)留言或發(fā)送郵件到jieandze1314@gmail.com

Welcome to our bioinfoplanet!