NMF算法简介
非负矩阵分解(Non-negative Matrix Factorization, NMF)的思想可以描述为，对于任意给定的一个非负矩阵V，NMF算法能够找到一个非负矩阵W和一个非负矩阵H，使得非负矩阵V≈W*H 成立 ，从而将一个非负的矩阵分解为左右两个非负矩阵的乘积。
利用矩阵分解来解决实际问题的分析方法很多，如PCA(主成分分析)、ICA(独立成分分析)、SVD(奇异值分解)等。在这些方法中，原始的大矩阵V被近似分解为低秩的V=WH形式。这些方法的共同特点是，即使输入的初始矩阵元素是全正的，传统的秩削减算法也不能保证原始数据的非负性，因子W和H中的元素往往含有负值元素。从数学和计算的观点看，分解结果中存在负值没有问题，但负值元素在实际问题中往往是没有意义的。NMF约束了原始矩阵V和分解矩阵W、H的非负性，这就意味着只能通过特征的相加来实现原始矩阵V的还原，最终导致的结果是：非负性会引发稀疏，非负性会使计算过程进入部分分解。
非负矩阵分解在单细胞数据分析中的主要应用：
1. 非负矩阵分解在做亚群细分和提取feature的时候是一个非常有效的工具。
2. 多样本整合，LIGER中用的是非负矩阵分解
3. CellChat做细胞通讯也用到了非负矩阵
4. 空间转录组去卷积工具SPOTlight也用到了NMF
使用NMF对scRNAseq的细胞进行分群操作
1. 读取数据，重新构建Seurat对象，保存矩阵和metadata信息，去除降维信息
⚠️这里使用pbmc数据集来演示NMF做细胞分群，实际上NMF在对亚群细分和提取feature时候效果更好，对亚群定义很有帮助。
library(NMF)
library(Seurat)
pbmc <- readRDS("pbmc.rds")
pbmc <- CreateSeuratObject(pbmc@assays$RNA@counts,
meta.data = pbmc@meta.data)
2. NMF
重新进行标准化和归一化，注意设置do.center = F，这样就不会得到负值
在做nmf时，rank值可以选的比目的预期的细胞类型/细胞状态稍微大一些的值，因为分解的一些因子会富集到线粒体核糖体等噪音，而不会落到一个具体的细胞亚群上面。（这里选了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（提取分解得到的每个因子）
⚠️：被非负矩阵分解出来的一些marker基因，一般是很有生物学意义的基因（因此适合于亚群细分和提取feature，便于细胞亚群注释）
# 每个因子提取30个
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. 选择用于后续分析的因子，使用NMF运行的结果进行降维和聚类
### 选择用于后续分析的因子
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. 降维聚类结果可视化
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")
注意：
如果细胞数比较多，在做NMF时，用于做非负矩阵分解的高变基因选择6000个左右比较合适，没有必要用全部基因来做。（这里演示用了2000个）
NMF运行很慢，在做大群定义结果和Seurat相差也不大，因此做大群定义的时候没有必要使用NMF。
