對(duì)應(yīng)原版教程第8章http://bioconductor.org/books/release/OSCA/overview.html
細(xì)胞的測(cè)序結(jié)果一般都有1萬(wàn)多的基因。從表達(dá)水平來(lái)看，其中大部分的信息可能都是較為“普通”的。如何從中抽取出能夠最大程度細(xì)胞間異質(zhì)性的基因是本小節(jié)的主要目的薄翅。

筆記要點(diǎn)

1沙兰、背景知識(shí)
1.1 為什么要挑選特定的基因
1.2 怎么挑HVG
2、基因表達(dá)的方差分解
2.1 基于log-counts的mean-variance擬合
2.2 基于normalization counts的CV2變異系數(shù)
3翘魄、如何根據(jù)計(jì)算指標(biāo)篩選高變基因
3.1 基于Top X排名
3.2 3.2 其它方法
4鼎天、確定高變基因之后sce對(duì)象的取舍
4.1 僅保留高變基因信息（不建議）
4.2 標(biāo)記高變基因，降維設(shè)置subset.row=參數(shù)（建議）
5熟丸、補(bǔ)充：關(guān)于“技術(shù)誤差”的進(jìn)一步分解

1、背景知識(shí)

1.1 為什么要挑選特定的基因

• 單細(xì)胞數(shù)據(jù)分析的主要在于考慮細(xì)胞（群）之間的異質(zhì)性伪节，最直接的方式就是通過(guò)相同基因在不同細(xì)胞間的差異表達(dá)來(lái)描述光羞。
• 不同細(xì)胞間的基因表達(dá)波動(dòng)差異可能會(huì)由技術(shù)誤差和生物水平的異質(zhì)性兩方面導(dǎo)致，而我們的目的就是找出主要由后者原因?qū)е碌母咦兓?HVG, highly variable genes)怀大，減少其它不太相關(guān)基因的信息干擾纱兑，用于下一步的降維。

1.2 怎么挑HVG

• 根據(jù)上面描述化借，最簡(jiǎn)單的想法就是計(jì)算基因在不同細(xì)胞間表達(dá)值的方差潜慎，取方差最大的基因即可；
• 但基因表達(dá)的總方差里有多少是技術(shù)誤差蓖康，有多少是生物水平差異是需要分解的铐炫；
• 而且需要考慮到mean-variance的影響，即基因表達(dá)均值越大蒜焊，可能受到越大的技術(shù)誤差倒信。

2、基因表達(dá)的方差分解

• 通過(guò)統(tǒng)計(jì)方法泳梆，計(jì)算基因的方差在多大程度上是由于不同類(lèi)型細(xì)胞差異表達(dá)造成的鳖悠；
• 下面介紹的兩種方法有個(gè)統(tǒng)一的假設(shè)就是：大部分的基因表達(dá)波動(dòng)是由于技術(shù)誤差造成的。
• 示例數(shù)據(jù)集如下优妙，已經(jīng)質(zhì)控乘综，標(biāo)準(zhǔn)化處理。具體步驟見(jiàn)原教程

sce.pbmc
#class: SingleCellExperiment 
#dim: 33694 3985 
#assays(2): counts logcounts

sce.416b
#class: SingleCellExperiment 
#dim: 46604 185 
#assays(2): counts logcounts

2.1 基于log-counts的mean-variance擬合

• 繪制基因平均表達(dá)量（log轉(zhuǎn)換）和方差的點(diǎn)圖套硼，擬合一條曲線卡辰，近似表示該表達(dá)量的基因受技術(shù)誤差導(dǎo)致的方差。

library(scran)
dec.pbmc <- modelGeneVar(sce.pbmc)
fit.pbmc <- metadata(dec.pbmc)
plot(fit.pbmc$mean, fit.pbmc$var, xlab="Mean of log-expression",
     ylab="Variance of log-expression")
curve(fit.pbmc$trend(x), col="dodgerblue", add=TRUE, lwd=2)

• 如下圖所示邪意，高于擬合曲線的點(diǎn)到擬合曲線的豎直距離表示有意義的生物水平造成的基因表達(dá)方差看政。

  
  
• 具體計(jì)算結(jié)果如下--bio = total - tech。（至于負(fù)數(shù)的結(jié)果本身是沒(méi)有意義的抄罕，可忽略）

dec.pbmc[order(dec.pbmc$bio, decreasing=TRUE),] 

image.png

• 引申兩個(gè)參數(shù)：
  1允蚣、擬合中的可能存在的過(guò)擬合問(wèn)題：這里主要針對(duì)基因表達(dá)量高，且方差大的少數(shù)點(diǎn)（基因）呆贿。造成曲線的膨脹,很大概率上高估了技術(shù)誤差嚷兔，擬合曲線表現(xiàn)為高蹺的尾巴森渐。可通過(guò)modelGeneVar()函數(shù)的density.weights參數(shù)設(shè)置(默認(rèn)為T(mén))冒晰。如下圖同衣，一般紅線的擬合情況是我們期望看到的(density.weights=FALSE)
  
  
  2、針對(duì)多批次的高變基因的挑選壶运，最直接的方法就是對(duì)每個(gè)批次單獨(dú)擬合好啰，再挑選。具體實(shí)現(xiàn)是在擬合函數(shù)中漓帅，設(shè)置block參數(shù)指定批次的設(shè)計(jì)
  modelGenemodelGeneVar(sce.416b, "ERCC", block=sce.416b$block)
  

2.2 基于normalization counts的CV2變異系數(shù)

• 變異系數(shù)(Coefficient of variation)也是均值水平歸一化的方差比較方法筷凤；
• 與2.1不同的是，這里使用的表達(dá)量是未經(jīng)log轉(zhuǎn)換的normalized counts棵癣；
• 但同樣是擬合不同表達(dá)量水平下辕翰，技術(shù)誤差造成的變異系數(shù)大小，然后取高于擬合曲線的點(diǎn)狈谊，認(rèn)為是能夠反應(yīng)生物水平差異的基因喜命。

dec.cv2.pbmc <- modelGeneCV2(sce.pbmc)
fit.cv2.pbmc <- metadata(dec.cv2.pbmc)
plot(fit.cv2.pbmc$mean, fit.cv2.pbmc$cv2, log="xy")
curve(fit.cv2.pbmc$trend(x), col="dodgerblue", add=TRUE, lwd=2)

• 具體計(jì)算結(jié)果如下--ratio = total / trend。注意ratio指標(biāo)是用除法河劝，而不是2.1里描述的減法壁榕。因?yàn)镃V2=(標(biāo)準(zhǔn)差/均值)2，需要排除分母均值的干擾赎瞎，才能使不同基因的CV2指標(biāo)具有可比性护桦。

dec.cv2.pbmc[order(dec.cv2.pbmc$ratio, decreasing=TRUE),]

• 如下圖結(jié)果，較2.1的方法煎娇，由于未進(jìn)行l(wèi)og壓縮二庵，更有可能挑選出均值水平較低的高變基因。

  
  

關(guān)于選擇2.1還是2.2計(jì)算基因的篩選指標(biāo)缓呛，如教程所說(shuō)催享，二者都能很好的捕獲基因的生物水平造成的差異性。但后面計(jì)算細(xì)胞間的距離是基于log-counts的方法哟绊，因此使用方法2.1可能會(huì)更具有統(tǒng)一性因妙。

3、如何根據(jù)計(jì)算指標(biāo)篩選高變基因

3.1 基于Top X排名

• 最簡(jiǎn)單的方法就是選擇上一步計(jì)算的高變基因指標(biāo)最顯著的一批基因票髓；
• 這會(huì)引出第二個(gè)問(wèn)題：選多少個(gè)高變基因合適攀涵？一般的選擇范圍可能從500到5000不等；
• 一方面洽沟，選出很多高變基因可能會(huì)包含部分不感興趣的無(wú)關(guān)基因的干擾以故；另一方面，選擇的高變基因過(guò)少裆操，又可能漏掉一些重要基因的信息怒详。
• 因此炉媒，最好先對(duì)測(cè)序的細(xì)胞群體的異質(zhì)性有一個(gè)大致評(píng)估。異質(zhì)性高的話昆烁，高變基因就定義多一點(diǎn)吊骤；相反，則定義少一點(diǎn)高變基因静尼。在此基礎(chǔ)上白粉，選一個(gè)大致的數(shù)量，跑接下來(lái)的流程鼠渺；如果發(fā)現(xiàn)不合適再回到這里修改鸭巴。而不是浪費(fèi)太多時(shí)間，糾結(jié)最佳取值這個(gè)難題系冗。

#基于log-count的計(jì)算指標(biāo)
dec.pbmc <- modelGeneVar(sce.pbmc)
hvg.pbmc.var <- getTopHVGs(dec.pbmc, n=1000)
#getTopHVGs(dec.pbmc, prop=0.1)
str(hvg.pbmc.var)
#chr [1:1000] "LYZ" "S100A9" "S100A8" "HLA-DRA" "CD74" "CST3" "TYROBP" "IGKC" "CCL5" "HLA-DRB1" ...

#基于log-count的計(jì)算指標(biāo)
dec.pbmc <- modelGeneVar(sce.pbmc, )
hvg.pbmc.var <- getTopHVGs(dec.pbmc, n=1000)
str(hvg.pbmc.var)
#chr [1:1000] "PPBP" "PRTFDC1" "HIST1H2AC" "FAM81B" "PF4" "GNG11" "KIAA1211" "HIST2H2BE" ...

3.2 其它方法

根據(jù)教程奕扣，再簡(jiǎn)單列舉幾種其它方法

3.2.1 基于FDR顯著性

• 根據(jù)顯著性的P值篩選高變基因薪鹦，具體是設(shè)置getTopHVGs()函數(shù)的fdr.threshold參數(shù)掌敬。

getTopHVGs(dec.pbmc, fdr.threshold=0.05)

• 對(duì)于log-count指標(biāo)，零假設(shè)就是bio<=0
• 對(duì)于CV2指標(biāo)池磁，零假設(shè)就是ratio<=1
• 根據(jù)設(shè)定的閾值奔害，可以控制假陽(yáng)性高變基因的數(shù)目。

3.2.2 取全部合格的高變基因

getTopHVGs(dec.pbmc, var.threshold=0)

• 對(duì)于log-count指標(biāo)地熄，取所有bio > 0的基因华临；
• 對(duì)于CV2指標(biāo)，取所有ratio > 1的基因端考。
• 不做推薦雅潭，因?yàn)榭赡芎蟹浅６嗟脑胍艋?/li>

之后教程還舉出一種方法是預(yù)先定義一組基因集作為作為高變基因。雖然與高變基因最初定義不符却特，但在某些情況下還是挺有意思的扶供，就不展開(kāi)了(其實(shí)沒(méi)怎么看明白)，可參看原教程裂明。

4椿浓、確定高變基因之后sce對(duì)象的取舍

確定高變基因后，后續(xù)降維操作主要根據(jù)能夠表示細(xì)胞間生物水平差異的高變基因展開(kāi)

4.1 僅保留高變基因信息（不建議）

dec.pbmc <- modelGeneVar(sce.pbmc)
chosen <- getTopHVGs(dec.pbmc, prop=0.1)
str(chosen)
sce.pbmc.hvg <- sce.pbmc[chosen,]
dim(sce.pbmc.hvg)
#[1] 1274 3985

4.2 標(biāo)記高變基因闽晦，降維設(shè)置subset.row=參數(shù)（建議）

# Performing PCA only on the chosen HVGs.
library(scater)
sce.pbmc <- runPCA(sce.pbmc, subset_row=chosen)
reducedDimNames(sce.pbmc)

• 小技巧：可以通過(guò)rowSubset()函數(shù)扳碍，將高變基因信息儲(chǔ)存到sce對(duì)象里(可以添加多次)

rowSubset(sce.pbmc,"hvg")=chosen
colnames(rowData(sce.pbmc))
rowData(sce.pbmc)[,2][rowData(sce.pbmc)[,3]]

5、補(bǔ)充：關(guān)于“技術(shù)誤差”的進(jìn)一步分解

• 如筆記一開(kāi)始所說(shuō)可將基因表達(dá)值的總方差分解為計(jì)數(shù)誤差和生物水平的誤差仙蛉。
  -是實(shí)際上笋敞，前者還包括由于轉(zhuǎn)錄擾動(dòng)造成的總體細(xì)胞表達(dá)水平的變化；可稱之為“uninteresting” biological variation荠瘪；例如管家基因house keeping gene液样，一般可忽略振亮，統(tǒng)稱為技術(shù)誤差。
• 但在某些情況下鞭莽，例如某些特定類(lèi)型細(xì)胞的特異基因受到調(diào)控作用顯著上調(diào)坊秸，即出現(xiàn)了一批（均值大、方差也大）的高變基因(hvg)澎怒，會(huì)造成擬合曲線的膨脹褒搔。即提高了該均值水平下的擬合技術(shù)誤差的標(biāo)準(zhǔn)，從而可能漏掉潛在的高變基因喷面。
• 針對(duì)此種情況星瘾，可使用外參轉(zhuǎn)錄本擬合單純的技術(shù)誤差。前提假設(shè)就是外參轉(zhuǎn)錄本的含量不受細(xì)胞內(nèi)的生物過(guò)程調(diào)控惧辈。
• 如果數(shù)據(jù)集中沒(méi)有外參轉(zhuǎn)錄本信息琳状，那么可使用泊松分布近似擬合技術(shù)誤差曲線。
• 相關(guān)函數(shù)如下盒齿，具體使用可參考原教程念逞。

modelGeneVarWithSpikes(sce.416b, "ERCC")
modelGeneVarByPoisson(sce.pbmc)