NMF算法簡介
非負(fù)矩陣分解(Non-negative Matrix Factorization, NMF)的思想可以描述為量蕊，對(duì)于任意給定的一個(gè)非負(fù)矩陣V痰洒，NMF算法能夠找到一個(gè)非負(fù)矩陣W和一個(gè)非負(fù)矩陣H哨查，使得非負(fù)矩陣V≈W*H 成立 回俐，從而將一個(gè)非負(fù)的矩陣分解為左右兩個(gè)非負(fù)矩陣的乘積吹艇。
利用矩陣分解來解決實(shí)際問題的分析方法很多麦牺，如PCA(主成分分析)钮蛛、ICA(獨(dú)立成分分析)、SVD(奇異值分解)等剖膳。在這些方法中魏颓，原始的大矩陣V被近似分解為低秩的V=WH形式。這些方法的共同特點(diǎn)是吱晒，即使輸入的初始矩陣元素是全正的甸饱，傳統(tǒng)的秩削減算法也不能保證原始數(shù)據(jù)的非負(fù)性，因子W和H中的元素往往含有負(fù)值元素仑濒。從數(shù)學(xué)和計(jì)算的觀點(diǎn)看叹话，分解結(jié)果中存在負(fù)值沒有問題，但負(fù)值元素在實(shí)際問題中往往是沒有意義的躏精。NMF約束了原始矩陣V和分解矩陣W渣刷、H的非負(fù)性，這就意味著只能通過特征的相加來實(shí)現(xiàn)原始矩陣V的還原矗烛，最終導(dǎo)致的結(jié)果是：非負(fù)性會(huì)引發(fā)稀疏辅柴，非負(fù)性會(huì)使計(jì)算過程進(jìn)入部分分解。
非負(fù)矩陣分解在單細(xì)胞數(shù)據(jù)分析中的主要應(yīng)用：
1. 非負(fù)矩陣分解在做亞群細(xì)分和提取feature的時(shí)候是一個(gè)非常有效的工具瞭吃。
2. 多樣本整合碌嘀，LIGER中用的是非負(fù)矩陣分解
3. CellChat做細(xì)胞通訊也用到了非負(fù)矩陣
4. 空間轉(zhuǎn)錄組去卷積工具SPOTlight也用到了NMF
使用NMF對(duì)scRNAseq的細(xì)胞進(jìn)行分群操作
1. 讀取數(shù)據(jù)，重新構(gòu)建Seurat對(duì)象歪架，保存矩陣和metadata信息股冗，去除降維信息
??這里使用pbmc數(shù)據(jù)集來演示NMF做細(xì)胞分群，實(shí)際上NMF在對(duì)亞群細(xì)分和提取feature時(shí)候效果更好和蚪，對(duì)亞群定義很有幫助止状。
library(NMF)
library(Seurat)
pbmc <- readRDS("pbmc.rds")
pbmc <- CreateSeuratObject(pbmc@assays$RNA@counts,
meta.data = pbmc@meta.data)
2. NMF
重新進(jìn)行標(biāo)準(zhǔn)化和歸一化烹棉，注意設(shè)置do.center = F，這樣就不會(huì)得到負(fù)值
在做nmf時(shí)怯疤，rank值可以選的比目的預(yù)期的細(xì)胞類型/細(xì)胞狀態(tài)稍微大一些的值浆洗，因?yàn)榉纸獾囊恍┮蜃訒?huì)富集到線粒體核糖體等噪音，而不會(huì)落到一個(gè)具體的細(xì)胞亞群上面集峦。（這里選了12）
pbmc <- NormalizeData(pbmc) %>% FindVariableFeatures() %>% ScaleData(do.center = F)
vm <- pbmc@assays$RNA@scale.data
saveRDS(vm, file = "pbmc_vm.rds")
res <- nmf(vm, 12, method = "snmf/r") #很慢
save(res, file = "pbmc_nmf_res.rda")
3. Extract cluster top loading features（提取分解得到的每個(gè)因子）
??：被非負(fù)矩陣分解出來的一些marker基因伏社，一般是很有生物學(xué)意義的基因（因此適合于亞群細(xì)分和提取feature，便于細(xì)胞亞群注釋）
# 每個(gè)因子提取30個(gè)
fs <- extractFeatures(res, 30L)
fs <- lapply(fs, function(x) rownames(res)[x])
fs <- do.call("rbind", fs)
rownames(fs) <- paste0("cluster", 1:12)
write.csv(t(fs), "pb,c_NMF_TopGenes.csv")
DT::datatable(t(fs))
4. 選擇用于后續(xù)分析的因子聪黎，使用NMF運(yùn)行的結(jié)果進(jìn)行降維和聚類
### 選擇用于后續(xù)分析的因子
s.f <- 1:12 # 因子 1 主要是線粒體和核糖體
## 降維
cell1 <- colnames(pbmc)
cell2 <- colnames(coef(res))
cells <- intersect(cell1, cell2)
pbmc <- pbmc[,cells]
pbmc <- RunPCA(pbmc, verbose = F)
pbmc@reductions$nmf <- pbmc@reductions$pca
pbmc@reductions$nmf@cell.embeddings <- t(coef(res)[,cells])
pbmc@reductions$nmf@feature.loadings <- basis(res)
pbmc <- RunUMAP(pbmc, reduction='nmf', dims=s.f)
## 基于NMF降維矩陣的聚類
pbmc <- FindNeighbors(pbmc, reduction='nmf', dims=s.f) %>% FindClusters()
## 基于因子最大載荷分類
pbmc$cluster <- apply(NMF::coefficients(res)[s.f,], 2, which.max)
5. 降維聚類結(jié)果可視化
p1 <- DimPlot(pbmc, label = T) + ggtitle("Clustered by Louvain")
p2 <- DimPlot(pbmc, group.by = "cluster", label = T) + ggtitle("Clustered by max loading")
pc <- p1|p2
ggsave("pbmc_NMF_Cluster.pdf", pc, width = 10, height = 5)
FeaturePlot(pbmc, features = c("CD8A","CD8B","CD3D","CD4","GZMA","NKG7"), ncol = 3)
saveRDS(pbmc, "pbmc_NMF.rds")
注意：
如果細(xì)胞數(shù)比較多，在做NMF時(shí)烘跺，用于做非負(fù)矩陣分解的高變基因選擇6000個(gè)左右比較合適湘纵，沒有必要用全部基因來做。（這里演示用了2000個(gè)）
NMF運(yùn)行很慢滤淳，在做大群定義結(jié)果和Seurat相差也不大梧喷，因此做大群定義的時(shí)候沒有必要使用NMF。
