劉小澤寫于2020.6.27
為何取名叫“交響樂(lè)”？因?yàn)閱渭?xì)胞分析就像一個(gè)大樂(lè)團(tuán)，需要各個(gè)流程的協(xié)同配合
單細(xì)胞交響樂(lè)1-常用的數(shù)據(jù)結(jié)構(gòu)SingleCellExperiment
單細(xì)胞交響樂(lè)2-scRNAseq從實(shí)驗(yàn)到下游簡(jiǎn)介
單細(xì)胞交響樂(lè)3-scRNA的細(xì)胞質(zhì)控
1 前言
scRNA數(shù)據(jù)中文庫(kù)的差異還是很常見(jiàn)的改橘，來(lái)源于：不同細(xì)胞的起始底物濃度不同，導(dǎo)致cDNA捕獲或PCR擴(kuò)增效率差異自脯。歸一化的目的就是去除細(xì)胞間與真實(shí)表達(dá)量無(wú)關(guān)的技術(shù)因素贞谓，方便后續(xù)比較。
這里需要說(shuō)明：歸一化與批次處理還是不同的叛溢。歸一化不管實(shí)驗(yàn)的批次因素塑悼，只考慮細(xì)胞中存在的技術(shù)誤差（比如測(cè)序深度），而批次處理既要考慮實(shí)驗(yàn)批次楷掉，又要考慮技術(shù)誤差（比如不同實(shí)驗(yàn)時(shí)間厢蒜、不同細(xì)胞系、不同文庫(kù)制備方法烹植、不同測(cè)序方法斑鸦、不同測(cè)序深度）。技術(shù)誤差對(duì)不同基因的影響是相似的草雕，比如有的技術(shù)誤差會(huì)影響基因長(zhǎng)度巷屿、GC含量，那么基本上所有基因都會(huì)受影響墩虹。但不同的批次產(chǎn)生的影響是不同的嘱巾，比如不同時(shí)間做的實(shí)驗(yàn)效果就是不一樣，那么基因的表達(dá)量可能也會(huì)受到這個(gè)影響诫钓。盡管有的R包可以同時(shí)處理技術(shù)誤差和批次效應(yīng)（例如zinbwave）旬昭，但大部分的分析還是各顧各的。
因此菌湃，這里只需要記得：歸一化问拘，歸的是技術(shù)誤差，而不是批次效應(yīng)即可
原文中提到：We will mostly focus our attention on scaling normalization, which is the simplest and most commonly used class of normalization strategies.
個(gè)人比較喜歡叫normalization為歸一化惧所，scale為標(biāo)準(zhǔn)化【只是一種稱呼習(xí)慣而已】骤坐。
當(dāng)然也有朋友喜歡叫normalization為標(biāo)準(zhǔn)化。其實(shí)纯路，scale與normalization的定義也不用掰扯太清楚或油，誰(shuí)包含誰(shuí)不用區(qū)分太明白。只需要知道分別做了什么事就可以了驰唬，一般來(lái)說(shuō)顶岸，應(yīng)該先進(jìn)行normalization腔彰，再進(jìn)行scale，而scale的結(jié)果會(huì)用來(lái)后續(xù)的降維和聚類
- scale：This involves dividing all counts for each cell by a cell-specific scaling factor, often called a “size factor” (Anders and Huber 2010)
The size factor for each cell represents the estimate of the relative bias in that cell, so division of its counts by its size factor should remove that bias.- Normalization "normalizes" within the cell for the difference in sequenicng depth / mRNA throughput. 主要著眼于樣本的文庫(kù)大小差異
Scaling "normalizes" across the sample for differences in range of variation of expression of genes . 主要著眼于基因的表達(dá)分布差異- 在Seurat中辖佣，一般得到原始表達(dá)矩陣并過(guò)濾后霹抛，會(huì)進(jìn)行
NormalizeData()，當(dāng)然也可以使用自己的歸一化結(jié)果（例如TPM結(jié)果卷谈，只是需要再log一下）杯拐；在挑選HVGs之后，降維處理之前世蔗，還需要用到ScaleData()進(jìn)行標(biāo)準(zhǔn)化端逼。對(duì)每個(gè)進(jìn)行center后的值再除以標(biāo)準(zhǔn)差（就是進(jìn)行了一個(gè)z-score的操作）【詳細(xì)探索可以看我之前寫的：https://mp.weixin.qq.com/s/6ioR3JE0wKg6M-YAsLBcTA】
需要說(shuō)明的是：本文的介紹基本都是基于Scater和Scran包，沒(méi)有Seurat的影子污淋，所以也看不到ScaleData的操作顶滩。大部分都是對(duì)文庫(kù)差異進(jìn)行的處理，集中體現(xiàn)在計(jì)算size factor寸爆，取log進(jìn)行數(shù)據(jù)縮放礁鲁，因此主要就是歸一化處理，而沒(méi)有看到z-score處理赁豆。這個(gè)也不用奇怪仅醇，每個(gè)包都有每個(gè)包的處理方法，也沒(méi)有一個(gè)定論說(shuō)就非要進(jìn)行Seurat的scale處理
數(shù)據(jù)準(zhǔn)備
數(shù)據(jù)依然給大家準(zhǔn)備好了魔种，sce.zeisel鏈接：https://share.weiyun.com/mNwJS8U9 密碼：g8379h
library(scRNAseq)
# 數(shù)據(jù)大概18M
sce.zeisel <- ZeiselBrainData()
rownames(sce.zeisel)
# 會(huì)看到有一些奇怪的基因名
# 將每個(gè)細(xì)胞的所有基因表達(dá)量加起來(lái)
library(scater)
sce.zeisel <- aggregateAcrossFeatures(sce.zeisel,
id=sub("_loc[0-9]+$", "", rownames(sce.zeisel)))
# 基因注釋
library(org.Mm.eg.db)
rowData(sce.zeisel)$Ensembl <- mapIds(org.Mm.eg.db,
keys=rownames(sce.zeisel), keytype="SYMBOL", column="ENSEMBL")
# 質(zhì)控（先perCellQCMetrics析二，后quickPerCellQC）
stats <- perCellQCMetrics(sce.zeisel, subsets=list(
Mt=rowData(sce.zeisel)$featureType=="mito"))
qc <- quickPerCellQC(stats, percent_subsets=c("altexps_ERCC_percent",
"subsets_Mt_percent"))
sce.zeisel <- sce.zeisel[,!qc$discard]
sce.zeisel
## class: SingleCellExperiment
## dim: 19839 2816
## metadata(0):
## assays(1): counts
## rownames(19839): 0610005C13Rik 0610007N19Rik ... Zzef1 Zzz3
## rowData names(2): featureType Ensembl
## colnames(2816): 1772071015_C02 1772071017_G12 ... 1772063068_D01
## 1772066098_A12
## colData names(10): tissue group # ... level1class level2class
## reducedDimNames(0):
## altExpNames(2): ERCC repeat
2 各種Size Factors的計(jì)算
2.1 文庫(kù)相關(guān)的size factor | library size normalization
這個(gè)的計(jì)算也是最簡(jiǎn)單最常用的
文庫(kù)大小指的就是：每個(gè)細(xì)胞中所有的基因表達(dá)量之和
進(jìn)行歸一化時(shí)，需要根據(jù)每個(gè)細(xì)胞的文庫(kù)大小各自計(jì)算一個(gè)”library size factor“务嫡，而這些factor的均值為1
lib.sf.zeisel <- librarySizeFactors(sce.zeisel)
# 說(shuō)明是對(duì)列操作
> length(lib.sf.zeisel)
[1] 2731
> ncol(sce.zeisel)
[1] 2731
# size factor 均值為1
> summary(lib.sf.zeisel)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.1725 0.5778 0.8701 1.0000 1.2716 4.0096
對(duì)size factor做個(gè)直方圖
hist(log10(lib.sf.zeisel), xlab="Log10[Size factor]", col='grey80')
文庫(kù)相關(guān)的size factor 基于的假設(shè)是：任意兩對(duì)細(xì)胞間的差異表達(dá)基因都是平衡的甲抖。也就是說(shuō)漆改，一組基因的上調(diào)表達(dá)程度心铃，勢(shì)必會(huì)被另一組基因的下調(diào)表達(dá)程度相抵消。但是在scRNA中差異基因表達(dá)并不是平衡的（也稱為：存在composition biases）挫剑，因此library size normalization可能不會(huì)得到準(zhǔn)確的歸一化結(jié)果去扣，僅僅是能用，但不是最精確樊破。
不過(guò)話又說(shuō)回來(lái)愉棱，精確的歸一化結(jié)果也不是scRNA數(shù)據(jù)探索的主要任務(wù)。還是要記得哲戚，歸一化的目的是更好地為下游展示數(shù)據(jù)做鋪墊奔滑。因此使用library size normalization對(duì)細(xì)胞分群和鑒定每群的marker基因也是足夠的。
2.2 去卷積方法計(jì)算size factor| Normalization by deconvolution
前面提到scRNA的數(shù)據(jù)存在composition biases顺少，也就是樣本間基因差異表達(dá)的不均衡性朋其。
先來(lái)理解composition biases王浴，現(xiàn)在想象一個(gè)例子：兩個(gè)細(xì)胞A、B梅猿，一個(gè)基因X在A細(xì)胞相比于B細(xì)胞表達(dá)量更高氓辣。這種高表達(dá)意味著：
- 當(dāng)細(xì)胞文庫(kù)大小一定（例如實(shí)驗(yàn)采用library quantification/文庫(kù)定量的方法），更多的測(cè)序資源偏向X基因袱蚓，勢(shì)必會(huì)降低其他本來(lái)不是差異表達(dá)基因的覆蓋度【簡(jiǎn)言之钞啸，一個(gè)盤子就這么多水果，一個(gè)異類胃口大開喇潘，吃了很多体斩，其他人本來(lái)也能吃飽，但現(xiàn)在也要挨餓】
- 當(dāng)不限制細(xì)胞文庫(kù)大小時(shí)颖低，X基因的reads或UMIs增加硕勿，也會(huì)增加總體的文庫(kù)大小，導(dǎo)致計(jì)算的library size factor增加枫甲，原本非差異基因的表達(dá)量也會(huì)由于size factor的增加而相對(duì)之前降低【簡(jiǎn)言之源武，一個(gè)異類拉高了整體平均成績(jī)，導(dǎo)致其他本來(lái)成績(jī)平平的學(xué)生看上去成績(jī)發(fā)生了下降】
這兩種情況都會(huì)導(dǎo)致一個(gè)結(jié)果：細(xì)胞A的其他非差異基因都會(huì)由于X基因的”過(guò)度表現(xiàn)“而不經(jīng)意地”下調(diào)“
那么如何去除這種composition biases呢想幻？之前在bulk RNA-Seq中研究非常透徹了粱栖。例如DESeq2使用的estimateSizeFactorsFromMatrix()，edgeR使用的calcNormFactors() 都在歸一化過(guò)程中考慮到了這點(diǎn)脏毯。它的假設(shè)是：大部分基因都不是差異基因闹究。
但是單細(xì)胞數(shù)據(jù)不能使用bulk 轉(zhuǎn)錄組的方法，因?yàn)榇嬖诖罅康牡捅磉_(dá)和0表達(dá)量食店。根據(jù)Lun, Bach, and Marioni 2016的研究渣淤，單細(xì)胞可以這么操作：
- 先混合：Pool counts from many cells to increase the size of the counts
- 后去卷積：Pool-based size factors are then “deconvolved” into cell-based factors for normalization of each cell’s expression profile.
library(scran)
set.seed(100)
clust.zeisel <- quickCluster(sce.zeisel)
table(clust.zeisel)
# clust.zeisel
# 1 2 3 4 5 6 7 8 9 10 11 12
# 215 147 237 565 162 350 276 291 175 109 103 101
# 然后去卷積
deconv.sf.zeisel <- calculateSumFactors(sce.zeisel, cluster=clust.zeisel)
summary(deconv.sf.zeisel)
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 0.1292 0.4954 0.8337 1.0000 1.3145 4.4807
quickCluster使用了 irlba 的PCA方法來(lái)進(jìn)行分群操作（注意這個(gè)分群和后面真正的分群操作不是一回事，這個(gè)頂多算是pre-clustering操作吉嫩，目的就是給去卷積鋪路）价认，并且它每次運(yùn)行結(jié)果都是隨機(jī)的，因此需要設(shè)置隨機(jī)種子來(lái)滿足重復(fù)性自娩。它得到的每群細(xì)胞都分別進(jìn)行歸一化用踩，得到size factors。
# 看一下這個(gè)數(shù)據(jù)包含了這么多的細(xì)胞類型
> table(sce.zeisel$level1class)
astrocytes_ependymal endothelial-mural
160 132
interneurons microglia
290 67
oligodendrocytes pyramidal CA1
749 938
pyramidal SS
395
# 將去卷積的size factor與之前常規(guī)方法得到的size factor進(jìn)行對(duì)比
plot(lib.sf.zeisel, deconv.sf.zeisel, xlab="Library size factor",
ylab="Deconvolution size factor", log='xy', pch=16,
col=as.integer(factor(sce.zeisel$level1class)))
abline(a=0, b=1, col="red") #截距為0忙迁，斜率為1
中間的紅線上的點(diǎn)表示去卷積和常規(guī)方法得到的size factor相等脐彩；不同顏色表示不同的細(xì)胞類型℃⑷樱可以看到去卷積的方法的確能暴露出細(xì)胞類型對(duì)文庫(kù)大小的影響惠奸，也證明了composition biases的確存在，并且在不同的細(xì)胞類型之間還很顯著恰梢。利用去卷積的size factor會(huì)得到更準(zhǔn)確的結(jié)果佛南，因?yàn)樗紤]的問(wèn)題更全面证九。
上面也說(shuō)過(guò)，使用library size normalization對(duì)細(xì)胞分群和鑒定每群的marker基因是足夠的共虑。這里使用去卷積得到更準(zhǔn)確的結(jié)果愧怜，雖然對(duì)分群改善不大，但對(duì)于每個(gè)基因的表達(dá)量數(shù)據(jù)估計(jì)與解釋還是很重要的妈拌，因?yàn)樗紤]了composition biases的影響拥坛。
2.3 使用spike-in計(jì)算size factor
它基于的假設(shè)是：外源RNA（spike-in）添加到每個(gè)細(xì)胞時(shí)的量都是一樣的。如果出現(xiàn)同一spike-in轉(zhuǎn)錄本表達(dá)量在各個(gè)細(xì)胞之間的差異尘分，問(wèn)題只能是細(xì)胞相關(guān)猜惋，例如不同細(xì)胞捕獲效率、不同細(xì)胞測(cè)序深度等等培愁。
因此著摔，計(jì)算spike-in size factor的目的也正是去除這方面的偏差。相比之前的兩種方法定续，spike-in的方法不再基于生物學(xué)背景假設(shè)（比如是否存在很多差異基因）谍咆，而是假設(shè)：
- 添加到每個(gè)細(xì)胞的spike-in的量都是一定的
- 對(duì)偏差的反應(yīng)，和內(nèi)源基因是一樣的
下面快速探索這個(gè)過(guò)程：
# 準(zhǔn)備數(shù)據(jù)
library(scRNAseq)
sce.richard <- RichardTCellData()
sce.richard <- sce.richard[,sce.richard$`single cell quality`=="OK"]
sce.richard #看到有ERCC存在
## class: SingleCellExperiment
## dim: 46603 528
## metadata(0):
## assays(1): counts
## rownames(46603): ENSMUSG00000102693 ENSMUSG00000064842 ...
## ENSMUSG00000096730 ENSMUSG00000095742
## rowData names(0):
## colnames(528): SLX-12611.N701_S502. SLX-12611.N702_S502. ...
## SLX-12612.i712_i522. SLX-12612.i714_i522.
## colData names(13): age individual ... stimulus time
## reducedDimNames(0):
## altExpNames(1): ERCC
使用computeSpikeFactors() 計(jì)算所有細(xì)胞的spike-in size factor
sce.richard <- computeSpikeFactors(sce.richard, "ERCC")
summary(sizeFactors(sce.richard))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1247 0.4282 0.6274 1.0000 1.0699 23.3161
再做一個(gè)spike-in與去卷積結(jié)果的對(duì)比圖：
to.plot <- data.frame(
DeconvFactor=calculateSumFactors(sce.richard),
SpikeFactor=sizeFactors(sce.richard),
Stimulus=sce.richard$stimulus,
Time=sce.richard$time
)
ggplot(to.plot, aes(x=DeconvFactor, y=SpikeFactor, color=Time)) +
geom_point() + facet_wrap(~Stimulus) + scale_x_log10() +
scale_y_log10() + geom_abline(intercept=0, slope=1, color="red")
發(fā)現(xiàn)每個(gè)處理的spike-in size factor與去卷積的結(jié)果呈正相關(guān)私股，表示它們都捕獲到了類似的技術(shù)偏差（測(cè)序深度摹察、捕獲效率）
同時(shí)注意到隨著時(shí)間的增加（如左上圖），spike-in factors在下降倡鲸，去卷積factors在增加供嚎，可以解釋為：隨著時(shí)間的推移，生物合成活性在增加峭状，總體RNA含量在增加克滴，而spike-in的含量是一定的，因此每個(gè)文庫(kù)中各個(gè)spike-in的比例減少（導(dǎo)致spike-in size factor結(jié)果減少）优床，而總體文庫(kù)增大（導(dǎo)致文庫(kù)size factor增加）
小結(jié)
對(duì)于歸一化方法的選擇劝赔，取決于生物學(xué)假設(shè)。大多數(shù)情況下羔巢，總體RNA含量的變化我們并不感興趣望忆，可以直接用library size或deconvolution factors當(dāng)做系統(tǒng)性偏差去除罩阵。但如果總體RNA含量變化與某個(gè)感興趣的生物學(xué)過(guò)程相關(guān)竿秆，例如細(xì)胞周期活性或T細(xì)胞激活，spike-in的歸一化方法就可以把這種變化保留下來(lái)稿壁，而去除其他系統(tǒng)性偏差幽钢。
另外，不管我們關(guān)不關(guān)心總體RNA含量傅是，對(duì)于spike-in轉(zhuǎn)錄本的歸一化都有特殊對(duì)待匪燕，不能使用針對(duì)內(nèi)源基因的方法蕾羊，比如可以利用modelGeneVarWithSpikes() 得到單獨(dú)的spike-in相關(guān)的size factor
3 得到Size Factors后進(jìn)行數(shù)據(jù)轉(zhuǎn)換
3.1 logNormCounts
log-transformation這種是最常用、最簡(jiǎn)單帽驯、最易理解的方法
一旦得到了size factors龟再，就可以利用scater包的logNormCounts()函數(shù)計(jì)算每個(gè)細(xì)胞歸一化后的結(jié)果。計(jì)算也很簡(jiǎn)單尼变，就是用某個(gè)細(xì)胞中每個(gè)基因或spike-in轉(zhuǎn)錄本的表達(dá)量除以這個(gè)細(xì)胞計(jì)算的size factor利凑，最后還進(jìn)行一個(gè)log轉(zhuǎn)換，得到一個(gè)新的矩陣：logcounts 【不過(guò)這個(gè)名字并不是很貼切嫌术，只是因?yàn)槠磳懞?jiǎn)單哀澈，真實(shí)應(yīng)該是：log-transformed normalized expression values。而不是單純字面意思取個(gè)log】
set.seed(100)
# 這里一看有pre-clustering就知道是利用的去卷積的size factor
clust.zeisel <- quickCluster(sce.zeisel)
sce.zeisel <- computeSumFactors(sce.zeisel, cluster=clust.zeisel, min.mean=0.1)
sce.zeisel <- logNormCounts(sce.zeisel)
# 最后的結(jié)果多了一項(xiàng)
assayNames(sce.zeisel)
## [1] "counts" "logcounts"
至于為什么取log度气？ 就是為了不受真實(shí)值的影響割按，而關(guān)注變化倍數(shù)。這里舉個(gè)例子磷籍，X基因在細(xì)胞A的表達(dá)量是50适荣，在B細(xì)胞的表達(dá)量是10；而Y基因在細(xì)胞A的表達(dá)量是1000院领，在B細(xì)胞的表達(dá)量是1100束凑。哪個(gè)基因更能吸引你的注意呢？我們是不是會(huì)關(guān)注變化的倍數(shù)栅盲？
另外汪诉，既然要取log，就應(yīng)該要注意存在很多0表達(dá)量的數(shù)值谈秫，因此log前還需要給表達(dá)量加上1扒寄，log1p = log(x + 1)，這個(gè)1就稱作：pseudo-count拟烫。
注意该编，每個(gè)包的log-transforming方法存在差異：
scran和scater使用的一樣，是logNormCounts硕淑，用表達(dá)量除以這個(gè)細(xì)胞計(jì)算的size factor课竣；
但Seurat使用的LogNormalize計(jì)算方法是：normalizes the feature expression measurements for each cell by the total expression, multiplies this by a scale factor (10,000 by default), and log-transforms the result，翻譯成公式就是：log1p(value/colSums[cell-idx] *scale_factor)
3.2 其他方法
其實(shí)默認(rèn)大家都在用log-transformation方法置媳，只是有時(shí)還有更特殊的需求
logNormCounts還有一個(gè)參數(shù)：downsample=T 于樟，當(dāng)使用常規(guī)流程做完后發(fā)現(xiàn)PCA分群的結(jié)果可能受到size factor的影響（即使之前通過(guò)歸一化試圖去除文庫(kù)大小的影響）。既然文庫(kù)大小影響還在拇囊，那么就試圖減弱這種影響迂曲，于是想到減少樣本數(shù)量，用更少的樣本再重新歸一化【但這種情況非常特殊】
圖中看到選取更少的樣本寥袭，結(jié)果變得清楚許多路捧，size factor的影響比之前小了很多关霸。于是可以推斷，之前PCA的結(jié)果很可能是還存在技術(shù)偏差導(dǎo)致的文庫(kù)差異
另外還有更復(fù)雜的方法可以看：sctransform 杰扫、DESeq2队寇，使用了variance stabilizing transformations方法
歡迎關(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
