今天我們來學習一個簡單的數(shù)學內(nèi)容，NMF

NMF算法簡介

NMF是什么？

圖片

非負矩陣分解(Non-negative Matrix Factorization, NMF)本質(zhì)上說是一種矩陣分解的方法，對于任意給定的一個非負矩陣V，NMF算法能夠?qū)ふ业揭粋€非負矩陣W和一個非負矩陣H，使得 V≈W*H 成立 ，從而將一個非負的矩陣分解為左右兩個非負矩陣的乘積。

NMF的特點？

NMF最重要的特點是非負性約束。

矩陣分解的方法有很多種，如奇異值分解 (singular value decomposition, SVD) 、獨立成分分析 (independent component analysis, ICA) 、主成分分析 (principal component analysis, PCA) 等。這些方法的共同特點是，即使初始矩陣 V 元素是非負的，分解出來的因子 W 和 H 中的元素往往含有負值元素。從計算科學的角度來看，分解出來的因子 W 和 H 中的元素含有負值元素并沒有問題， 但負值元素通常是無法解釋的。NMF約束了原始矩陣V和分解矩陣W、H的非負性，這就意味著只能通過特征的相加來實現(xiàn)原始矩陣V的還原，最終導(dǎo)致的結(jié)果是：

• 非負性會引發(fā)稀疏
• 非負性會使計算過程進入部分分解

給大家對比一下PCA與NMF分解圖像的效果：

圖片.png

PCA分解之后，每個主成分（PC）都會保留與其他PC不正交的全局特征，并且PC保留的特征是遞減的。

圖片.png

NMF分解之后，每個因子保留的都是局部特征，它們的權(quán)重是基本平等的。通過這張圖可以看出，很多因子能與面部特征一一對應(yīng)起來，例如鼻子、眼睛、嘴巴都能找到相應(yīng)的因子。

NMF在單細胞研究中的優(yōu)勢

單細胞研究避免不了要回答兩個問題：組織中有哪些細胞類型，每個細胞類型又有哪些表達模式？NMF解決這類問題具有天然的優(yōu)勢，因為它分解的因子很容易與細胞類型或表達模式對應(yīng)起來。Github上有很多基于NMF和其變種算法的單細胞分析工具，我比較喜歡的有單細胞整合分析工具liger和空間轉(zhuǎn)錄組去卷積工具SPOTlight。應(yīng)用NMF分析方法發(fā)表的高分文章也有很多，我給大家介紹一篇，更多的文章請自己搜索。

Chen, YP., Yin, JH., Li, WF. et al. Single-cell transcriptomics reveals regulators underlying immune cell diversity and immune subtypes associated with prognosis in nasopharyngeal carcinoma. Cell Res 30, 1024–1042 (2020). https://doi.org/10.1038/s41422-020-0374-x

基礎(chǔ)NMF包的安裝與用法簡介

安裝NMF基礎(chǔ)包

BiocManager::install('Biobase')

nmf函數(shù)簡介

NMF包通過nmf()函數(shù)實現(xiàn)矩陣分解，它的用法及重要參數(shù)如下：

nmf(x, rank, method, seed, nrun, ...)

rank值一般是要通過測試評估后確定的，但是分析單細胞數(shù)據(jù)這是一個很難完成的工作，5000個細胞的測試時間可能超過10個小時。替代辦法是使用經(jīng)驗或先驗知識指定，可以嘗試略多于細胞類型或細胞狀態(tài)（細胞亞群再聚類時）的一個數(shù)值，例如我在本帖的PBMC數(shù)據(jù)分解中就指定為rank=10。因為NMF一般是從隨機數(shù)開始，通過迭代算法收斂誤差的方法求出最優(yōu)W和H矩陣，所以seed不同最后的結(jié)果也不同。為了減少seed的影響求得最優(yōu)解，常規(guī)的辦法是通過nrun參數(shù)設(shè)置運行100-200次矩陣分解選取最優(yōu)值，也可以使用特殊的算法選擇一個最佳的seed（設(shè)置seed='nndsvd'或seed='ica'），這樣運行一次也能得到最優(yōu)解。下面我們測試一下不同方法的運行時間：

library(Biobase)

10X PBMC數(shù)據(jù)實測

測試數(shù)據(jù)

測試數(shù)據(jù)來源于10x genomics官網(wǎng)的示例數(shù)據(jù)集。

圖片

數(shù)據(jù)下載鏈接：https://cf.10xgenomics.com/samples/cell-exp/3.0.2/5k_pbmc_v3_nextgem/5k_pbmc_v3_nextgem_filtered_feature_bc_matrix.h5

保存文件名：pbmc.h5

基于PCA分解的降維聚類

library(Seurat)

圖片

圖片

基于NMF分解的降維聚類

## 高變基因表達矩陣的分解
# pbmc大體可分成T，B，NK，CD14+Mono，CD16+Mono，DC，Platelet等類型，考慮冗余后設(shè)置rank=10
vm <- pbmc@assays$RNA@scale.data
res <- nmf(vm, 10, method = "snmf/r", seed = 'nndsvd') 
runtime(res)
#    用戶     系統(tǒng)     流逝 
#1063.147   78.019 1139.831 

## 分解結(jié)果返回suerat對象
pbmc@reductions$nmf <- pbmc@reductions$pca
pbmc@reductions$nmf@cell.embeddings <- t(coef(res))    
pbmc@reductions$nmf@feature.loadings <- basis(res)  

## 使用nmf的分解結(jié)果降維聚類
set.seed(219)
pbmc.nmf <- RunUMAP(pbmc, reduction = 'nmf', dims = 1:10) %>% 
  FindNeighbors(reduction = 'nmf', dims = 1:10) %>% FindClusters()

## 結(jié)果可視化  
p <- DimPlot(pbmc.nmf, label = T) + ggsci::scale_color_igv()
ggsave("pbmc_nmf.png", p, width = 9, height = 6)
p <- FeaturePlot(pbmc.nmf, features = c('CD3D', 'CD3E', 'MS4A1', 'CD79A', 'GNLY', 'NKG7', 'CD14', 
                                        'FCGR3A', 'PPBP', 'FCER1A', 'CD4', 'CD8A'), ncol = 4)
ggsave("pbmc_nmf_markers.png", p, width = 12, height = 8)

圖片.png

圖片

NMF因子可解釋性探索

對比PCA分析的結(jié)果，NMF雖然毫不遜色，但是它的運行時間更長，我們?yōu)槭裁匆肗MF呢？一個很重要的原因是NMF的因子可解釋性更強，每個因子貢獻度最大的基因基本代表了某種或某個狀態(tài)細胞的表達模式，相比差異分析得到marker基因更有代表性。

NMF因子與細胞類型的關(guān)系

## 人工定義細胞類型
pbmc.nmf$celltype <- pbmc.nmf$seurat_clusters
pbmc.nmf$celltype <- recode(pbmc.nmf$celltype,
                            '1' = "B cells", 
                            '5' = "NKs", 
                            '10' = "CD8+ T", 
                            '6' = "CD8+ T", 
                            '4' = "CD4+ T",
                            '9' = "CD4+ T", 
                            '0' = "CD4+ T", 
                            '3' = "CD4+ T", 
                            '7' = "CD4+ T", 
                            '2' = "CD14+ Mono",
                            '8' = "CD14+ Mono", 
                            '11' = "CD16+ Mono", 
                            '12' = "DCs", 
                            '13' = "Platelet", 
                            '14' = "Unknown")
p <- DimPlot(pbmc.nmf, group.by = 'celltype', label = T, label.size = 3) + ggsci::scale_color_npg(alpha = 0.6)
ggsave("pbmc_nmf_celltype.png", p, width = 9, height = 6)

## 查看細胞因子上的荷載
tmp <- data.frame(t(coef(res)), check.names = F)
colnames(tmp) <- paste0("factor", 1:10)
pbmc.nmf <- AddMetaData(pbmc.nmf, metadata = tmp)
p <- FeaturePlot(pbmc.nmf, features = paste0("factor", 1:10), ncol = 4)
ggsave("pbmc_nmf_factors.png", p, width = 12, height = 8)

## 查看細胞主成分上的荷載
p <- FeaturePlot(pbmc.nmf, features = paste0("PC_", 1:12), ncol = 4)
ggsave("pbmc_nmf_PCs.png", p, width = 12, height = 8)

人工鑒定的細胞類型

圖片

細胞在因子上的值

圖片

細胞在PC軸上的值

圖片.png

對比上下兩張圖，很容易發(fā)現(xiàn)NMF的因子比PCA的PC軸解釋性更強。

提取celltype的signatures

## 提取每個因子貢獻度最大的20個基因
f <- extractFeatures(res, 20L)
f <- lapply(f, function(x) rownames(res)[x])
f <- do.call("rbind", f)
DT::datatable(t(f))

圖片

NMF的10個因子中，很容易發(fā)現(xiàn)我們對細胞進行分類的marker基因，是不是很神奇？

生活很好，等你超越