Seurat
Seurat整合流程與原理
使用CCA分析將兩個(gè)數(shù)據(jù)集降維到同一個(gè)低維空間，因?yàn)镃CA降維之后的空間距離不是相似性而是相關(guān)性娃循，所以相同類(lèi)型與狀態(tài)的細(xì)胞可以克服技術(shù)偏倚重疊在一起炕檩。
-
CCA降維之后細(xì)胞在低維空間有了可以度量的“距離”，MNN(mutual nearest neighbor)算法以此找到兩個(gè)數(shù)據(jù)集之間互相“距離”最近的細(xì)胞捌斧，Seurat將這些相互最近鄰細(xì)胞稱(chēng)為“錨點(diǎn)細(xì)胞”笛质。我們用兩個(gè)數(shù)據(jù)集A和B來(lái)說(shuō)明錨點(diǎn)，假設(shè)：
- A樣本中的細(xì)胞A3與B樣本中距離最近的細(xì)胞有3個(gè)（B1,B2,B3）
- B樣本中的細(xì)胞B1與A樣本中距離最近的細(xì)胞有4個(gè)（A1,A2,A3,A4）
- B樣本中的細(xì)胞B2與A樣本中距離最近的細(xì)胞有2個(gè)（A5,A6）
- B樣本中的細(xì)胞B3與A樣本中距離最近的細(xì)胞有3個(gè)（A1,A2,A7）
從圖中可以看出捞蚂，A3與B1是“相互“最近鄰的細(xì)胞妇押，而A3與B2、B3不是相互最近鄰細(xì)胞姓迅，A3+B1就是A敲霍、B兩個(gè)數(shù)據(jù)集中的錨點(diǎn)之一。實(shí)際數(shù)據(jù)中丁存，兩個(gè)數(shù)據(jù)集之間的錨點(diǎn)可能有幾百上千個(gè)肩杈。
-
理想情況下，相同類(lèi)型和狀態(tài)的細(xì)胞才能構(gòu)成配對(duì)錨點(diǎn)細(xì)胞解寝，但是異常的情況也會(huì)出現(xiàn)扩然，如C圖中query數(shù)據(jù)集中黑色的細(xì)胞團(tuán)。它在reference數(shù)據(jù)集沒(méi)有相同類(lèi)型的細(xì)胞编丘，但是它也找到了錨點(diǎn)配對(duì)細(xì)胞（紅色連線）与学。Seurat會(huì)通過(guò)兩步過(guò)濾這些不正確的錨點(diǎn)：
- 在CCA低維空間找到的錨點(diǎn)彤悔，返回到基因表達(dá)數(shù)據(jù)構(gòu)建的高維空間中驗(yàn)證，如果它們的轉(zhuǎn)錄特征相似性高則保留索守，否則過(guò)濾此錨點(diǎn)晕窑。
- 檢查錨點(diǎn)細(xì)胞所在數(shù)據(jù)集最鄰近的30個(gè)細(xì)胞，查看它們重疊的錨點(diǎn)配對(duì)細(xì)胞的數(shù)量卵佛，重疊越多分值越高杨赤，代表錨點(diǎn)可靠性更高。
經(jīng)過(guò)層層過(guò)濾剩下的錨點(diǎn)細(xì)胞對(duì)截汪，可以認(rèn)為它們是相同類(lèi)型和狀態(tài)的細(xì)胞疾牲，它們之間的基因表達(dá)差異是技術(shù)偏倚引起的。Seurat計(jì)算它們的差異向量衙解，然后用此向量校正這個(gè)錨點(diǎn)錨定的細(xì)胞子集的基因表達(dá)值阳柔。校正后的基因表達(dá)值即消除了技術(shù)偏倚，實(shí)現(xiàn)了兩個(gè)單細(xì)胞數(shù)據(jù)集的整合蚓峦。
Seurat工作流程圖
Seurat示例代碼
pbmc.list <- SplitObject(pbmc, split.by="batch")
for(i in 1:length(pbmc.list)) {
pbmc.list[[i]] <- NormalizeData(pbmc.list[[i]],verbose=FALSE)
pbmc.list[[i]] <- FindVariableFeatures(pbmc.list[[i]],selection.method="vst",verbose=FALSE)
}
AllBatch.anchors <- FindIntegrationAnchors(object.list = pbmc.list, dims = 1:20,k.filter=80)
pbmc <- IntegrateData(anchorset = AllBatch.anchors, dims = 1:20)
DefaultAssay(pbmc) <- "integrated"
pbmc <- ScaleData(pbmc, verbose = T, do.center=F,split.by = 'batch')
pbmc <- RunPCA(pbmc, npcs = 50, verbose = T)
pbmc <- RunUMAP(pbmc, reduction = "pca", dims = 1:20,
reduction.name = 'seurat_umap', reduction.key = 'seurat_umap')
DimPlot(pbmc, reduction='seurat_umap', group.by='batch', pt.size=0.3,label=F)
參考
Tutorial：https://satijalab.org/seurat/articles/integration_mapping.html
Paper：https://www.cell.com/cell/fulltext/S0092-8674(19)30559-8
LIGER
LIGER原理
- LIGER是從多個(gè)批次的數(shù)據(jù)中挑出一組共同的潛在因子舌剂，這些因子代表了每個(gè)數(shù)據(jù)集中相同的生物“信號(hào)”，但同時(shí)也保留了數(shù)據(jù)集的差異暑椰。為的是用這些因子聯(lián)合識(shí)別細(xì)胞類(lèi)型霍转，以及識(shí)別和保留細(xì)胞類(lèi)型特異性的差異。
- LIGER首先采用綜合非負(fù)矩陣分解（iNMF）來(lái)學(xué)習(xí)低維空間一汽，首先對(duì)每個(gè)批次的單細(xì)胞基因表達(dá)矩陣經(jīng)過(guò)iNMF分解避消，可以得到三個(gè)矩陣：因子(行) x 細(xì)胞 (列)的H矩陣，數(shù)據(jù)集特有基因(行) x 因子 (列)的V矩陣召夹，以及數(shù)據(jù)集共有基因(行) x 因子 (列)的W矩陣岩喷。
iNMF原理
- 我們知道NMF的優(yōu)化目標(biāo)是在迭代中優(yōu)化以下誤差
NMF優(yōu)化目標(biāo)
- 而當(dāng)我們有多個(gè)數(shù)據(jù)集并想將其整合起來(lái)的時(shí)候，就需要使用iNMF戳鹅。設(shè)有k個(gè)待整合的數(shù)據(jù)集均驶，iNMF對(duì)NMF的形式進(jìn)行如下修正：
iNMF的優(yōu)化目標(biāo)
我們可以假定待整合的不同數(shù)據(jù)集之間必然存在著共有的一些特征，也必然存在著一些獨(dú)特的特征枫虏。因此我們可以把原始W矩陣寫(xiě)成W矩陣與Vk矩陣的疊加妇穴，新的W矩陣則相當(dāng)于不同數(shù)據(jù)集之間共有的特征，而Vk矩陣則是第k個(gè)數(shù)據(jù)集獨(dú)特的特征隶债。另外引入λ正則化項(xiàng)腾它，λ的大小反映了數(shù)據(jù)集之間的相似性。當(dāng)λ等于0時(shí)死讹，VH的二范數(shù)對(duì)目標(biāo)函數(shù)沒(méi)有影響瞒滴，因此在優(yōu)化的過(guò)程中我們得到的W矩陣（也就是共有的特征）會(huì)趨向于0，而數(shù)據(jù)集獨(dú)特的特征則會(huì)使得Error項(xiàng)盡可能的小，因此數(shù)據(jù)集之間的差異（也就是數(shù)據(jù)集獨(dú)特的特征）就會(huì)盡可能的大妓忍。同理虏两，當(dāng)λ趨于正無(wú)窮大時(shí)，減小VH的二范數(shù)對(duì)優(yōu)化目標(biāo)函數(shù)的權(quán)重極大世剖，在這種情況下定罢，W矩陣，也就是數(shù)據(jù)集共有的特征旁瘫，就會(huì)盡可能的大祖凫。
因此在進(jìn)行非負(fù)矩陣分解的過(guò)程中，我們需要選擇合適的參數(shù)K和λ酬凳，注意惠况，這里的K指的是W矩陣的維度。在最理想的狀態(tài)下宁仔，數(shù)據(jù)集中的每一個(gè)真實(shí)存在的cluster將占據(jù)W矩陣的一個(gè)緯度稠屠，此時(shí)數(shù)據(jù)之間分散得越開(kāi)，H矩陣的KL散度就越大翎苫。因此參數(shù)K的選擇可以基于H矩陣的KL散度進(jìn)行優(yōu)化完箩，當(dāng)KL散度剛剛飽和時(shí)的K應(yīng)當(dāng)是最優(yōu)的K。同理拉队，一般而言，我們希望在進(jìn)行數(shù)據(jù)整合之后阻逮，不同數(shù)據(jù)集之間能夠盡可能均勻的map在一起粱快，這種均勻的程度可以用Alignment score進(jìn)行衡量，隨著λ的增大叔扼，我們假定的數(shù)據(jù)集之間的相似性越高事哭，通常Alignment score就會(huì)越高。因此參數(shù)λ的選擇可以基于聚類(lèi)的Alignment score進(jìn)行優(yōu)化瓜富，當(dāng)Alignment score剛剛飽和時(shí)的λ應(yīng)當(dāng)是最優(yōu)的λ鳍咱。但必須要指出的是，這種參數(shù)的優(yōu)化僅僅基于通常無(wú)先驗(yàn)知識(shí)的情況与柑，我們更建議以此為參考谤辜，根據(jù)我們對(duì)數(shù)據(jù)的理解和先驗(yàn)知識(shí)選擇合適的參數(shù)。
- iNMF分解矩陣后的“因子”价捧，相當(dāng)于PCA降維之后的"PC"軸丑念；但是iNMF的因子比 PCA的PC軸可解釋性更強(qiáng)，在單細(xì)胞數(shù)據(jù)分析中结蟋，每個(gè)因子常常對(duì)應(yīng)一種細(xì)胞類(lèi)型脯倚。V矩陣和W矩陣中l(wèi)oading值靠前的基因，往往具有顯著的生物學(xué)意義。
- Liger的聚類(lèi)算法也比較有意思推正，它并不基于單個(gè)細(xì)胞的坐標(biāo)（label）進(jìn)行聚類(lèi)恍涂，而是用單個(gè)細(xì)胞周?chē)鷎個(gè)鄰居細(xì)胞的坐標(biāo)（label）為該細(xì)胞賦予一個(gè)新的坐標(biāo)，再根據(jù)曼哈頓距離或其他的范數(shù)距離進(jìn)行聚類(lèi)植榕。Liger作者的觀點(diǎn)是單個(gè)細(xì)胞被錯(cuò)誤分配坐標(biāo)（label）的概率要遠(yuǎn)大于其周?chē)鷎個(gè)鄰居細(xì)胞坐標(biāo)都被系統(tǒng)性地錯(cuò)誤分配的概率要高得多再沧，因此使用周?chē)鷎個(gè)鄰居細(xì)胞的坐標(biāo)來(lái)替代該細(xì)胞本身的坐標(biāo)會(huì)使聚類(lèi)的結(jié)果更加魯棒。
LIGER工作流程圖
LIGER示例代碼
library(liger)
library(SeuratWrappers)
pbmc <- NormalizeData(pbmc)
pbmc <- FindVariableFeatures(pbmc)
pbmc <- ScaleData(pbmc, split.by = "batch", do.center = FALSE)
pbmc <- RunOptimizeALS(pbmc, k = 20, lambda = 5, split.by = "batch")
pbmc <- RunQuantileNorm(pbmc, split.by = "batch")
pbmc <- RunUMAP(pbmc, dims = 1:ncol(pbmc[["iNMF"]]), reduction = "iNMF",
reduction.name = 'liger_umap', reduction.key = 'liger_umap')
DimPlot(pbmc, group.by = "batch", label=F, reduction = 'liger_umap')
參考
Tutorial：https://htmlpreview.github.io/?https://github.com/satijalab/seurat.wrappers/blob/master/docs/liger.html
Paper：https://www.cell.com/cell/fulltext/S0092-8674%2819%2930504-5
https://www.plob.org/article/24538.html
Harmony
Harmony原理
- Harmony需要輸入低維空間的坐標(biāo)值（embedding）内贮，一般使用PCA的降維結(jié)果产园。Harmony導(dǎo)入PCA的降維數(shù)據(jù)后，會(huì)采用soft k-means clustering算法將細(xì)胞聚類(lèi)夜郁。常用的聚類(lèi)算法僅考慮細(xì)胞在低維空間的距離什燕，但是soft clustering算法會(huì)考慮我們提供的校正因素。
- 例如竞端，細(xì)胞c2距離cluster1有點(diǎn)遠(yuǎn)屎即，本來(lái)不能算作cluster1的一份子；但是c2和cluster1的細(xì)胞來(lái)自不同的數(shù)據(jù)集事富，因?yàn)槲覀兤谕煌臄?shù)據(jù)集融合技俐，所以破例讓它加入cluster1了。
- 聚類(lèi)之后统台，先計(jì)算每個(gè)cluster內(nèi)各個(gè)數(shù)據(jù)集的細(xì)胞的中心點(diǎn)雕擂，然后根據(jù)這些中心點(diǎn)計(jì)算各個(gè)cluster的中心點(diǎn)。最后通過(guò)算法讓cluster內(nèi)的細(xì)胞向中心聚集贱勃，實(shí)在收斂不了的離群細(xì)胞就過(guò)濾掉井赌。調(diào)整之后的數(shù)據(jù)重復(fù)：聚類(lèi)—計(jì)算cluster中心點(diǎn)—收斂細(xì)胞—聚類(lèi)的過(guò)程，不斷迭代直至聚類(lèi)效果趨于穩(wěn)定贵扰。
Harmony工作流程圖
Harmony示例代碼
library(harmony)
pbmc <- NormalizeData(pbmc)
pbmc <- FindVariableFeatures(pbmc)
pbmc <- ScaleData(pbmc, split.by = "batch", do.center = FALSE)
pbmc <- RunPCA(pbmc, npcs = 50, verbose = FALSE)
pbmc <- RunHarmony(pbmc, group.by.vars = "batch")
pbmc <- RunUMAP(pbmc, reduction = "harmony", dims = 1:20,
reduction.name = 'harmony_umap', reduction.key = 'harmony_umap')
DimPlot(pbmc, reduction='harmony_umap', label=F, group.by = 'batch')
參考
Tutorial：http://htmlpreview.github.io/?https://github.com/immunogenomics/harmony/blob/master/docs/SeuratV3.html
Paper：https://www.nature.com/articles/s41592-019-0619-0
BEER
BEER的主要步驟包含兩步
- 首先對(duì)每個(gè)批次的單細(xì)胞矩陣進(jìn)行降維仇穗，對(duì)相似的10個(gè)細(xì)胞進(jìn)行分組，產(chǎn)生一個(gè)代表性的表達(dá)譜戚绕，然后找出兩個(gè)批次中所有相似的細(xì)胞組（Mutual Nearest, MN）纹坐。
- 將兩個(gè)批次整合，然后進(jìn)行PCA舞丛。對(duì)于一對(duì)MN耘子，如果沒(méi)有批次效應(yīng)，那么他們?cè)谀硞€(gè)PC上應(yīng)該具有相似的值球切；BEER刪除那些在一對(duì)MN具有低相關(guān)性的PC拴还，以達(dá)到去除批次效應(yīng)的目的。
BEER工作流程圖
BEER示例代碼
source("../BEER-master/BEER.R")
DATA <- pbmc@assays$RNA@counts
DATA.BATCH <- pbmc$batch
RM1=read.table('../data/Apoptosis.txt',sep='\t',header=FALSE)[,1]
RM2=read.table('../data/CellCycle.txt',sep='\t',header=FALSE)[,1]
RM3=read.table('../data/Ribosome.txt',sep='\t',header=FALSE)[,1]
RMG=c(as.character(RM1),as.character(RM2),as.character(RM3))
mybeer <- BEER(DATA, DATA.BATCH, GNUM=30, PCNUM=50, ROUND=1, GN=2000, SEED=1, COMBAT=TRUE, RMG=RMG)
mybeer <- ReBEER(mybeer, GNUM=30, PCNUM=50, ROUND=1, SEED=1, RMG=RMG)
PCUSE <- mybeer$select
PCUSE <- PCUSE[which(PCUSE <= 20)]
# PCUSE=.selectUSE(mybeer, CUTR=0.8, CUTL=0.8, RR=0.5, RL=0.5)
mybeer$seurat@meta.data <- pbmc@meta.data
pbmc0 <- mybeer$seurat
pbmc0 <- RunUMAP(object = pbmc0, reduction='pca',dims = PCUSE, check_duplicates=FALSE,
reduction.name = "beer_umap",reduction.key = "beer_umap",
min.dist = 0.3, n.neighbors = 200L)
DimPlot(pbmc0, reduction='beer_umap', group.by='batch', pt.size=0.1) + NoLegend()
參考
Tutorial：https://github.com/jumphone/BEER
Paper：https://www.nature.com/articles/s41421-019-0114-x
參考
- https://mp.weixin.qq.com/s/XUm-1Xi4StA6N7sz9CA1IA
- https://mp.weixin.qq.com/s/txnz-wfBk5vLpBiuGBnNTQ
- https://mp.weixin.qq.com/s?__biz=MzIyMzMwNDQ2MA==&mid=2247484164&idx=1&sn=9da8c20591b95cc6663cea2c5cb195bd&chksm=e8210d67df568471612b6426a599e618d0158e734a0f8c3fafa4b1eaa2b7adaec373a30a2ee6&scene=178&cur_album_id=1593332494622359552#rd
