【原】GeneNMF: 基于NMF（非負(fù)矩陣分解）的基因程序識(shí)別及可視化

TS的美夢(mèng) 2025-01-10

展開(kāi)全文

Non-negative matrix factorization（NMF）非負(fù)矩陣分解是一種popular的矩陣分解方法，它允許從一組非負(fù)數(shù)據(jù)向量中提取稀疏和有意義的特征。它非常適合分解scRNA-seq數(shù)據(jù)，有效地將大型復(fù)雜矩陣簡(jiǎn)化為幾個(gè)可解釋的基因程序。GeneNMF是一個(gè)實(shí)現(xiàn)單細(xì)胞矩陣分解和基因程序發(fā)現(xiàn)方法的軟件包。它可以直接應(yīng)用于Seurat對(duì)象，以降低數(shù)據(jù)的維數(shù)，并在多個(gè)樣本中檢測(cè)穩(wěn)定的基因程序。GeneNMF運(yùn)行分析很簡(jiǎn)單，該方法受先前將非負(fù)矩陣分解（NMF）應(yīng)用于單細(xì)胞 RNA 測(cè)序（scRNA-seq）數(shù)據(jù)工作的啟發(fā)。相比于NMF，速度也是很快的。

GeneNMF github：https://github.com/carmonalab/GeneNMF?tab=readme-ov-file

官網(wǎng)提供了詳細(xì)的教程，可以閱讀學(xué)習(xí)！本文代碼及數(shù)據(jù)已發(fā)布微信VIP群！

安裝包：

setwd('/home/tq_ziv/data_analysis/NMF_singlecell/')#===============================================================================#安裝加載R包# install.packages("GeneNMF")library(remotes)remotes::install_github("carmonalab/GeneNMF")library(GeneNMF)#R>=4.3library(Seurat)library(ggplot2)library(UCell)library(patchwork)library(Matrix)library(RcppML)library(viridis)

GeneNMF用于降維，這里我們直接演示看看大群的效果，我們的數(shù)據(jù)是已經(jīng)降維聚類注釋好的數(shù)據(jù)，只不過(guò)用的是TSNE，我們這里用UMAP看看。其實(shí)和nmf一樣，大群降維，nmf并沒(méi)有什么明顯的優(yōu)勢(shì)，seuat默認(rèn)的降維方式分群已經(jīng)就可以了！

seu <- FindVariableFeatures(uterus, nfeatures = 3000)seu <- runNMF(seu, k = 10, assay="RNA")# seu@reductions[["NMF"]]# A dimensional reduction object with key NMF_ # Number of dimensions: 10 # Projected dimensional reduction calculated:  FALSE # Jackstraw run: FALSE # Computed using assay: RNA 
#我們這個(gè)數(shù)據(jù)是已經(jīng)跑完降維并做好細(xì)胞注釋的數(shù)據(jù)，只是跑的是TSNE#我們可以這里跑以下UMAP，分別用pca和nmf降維，看看效果！seu <- RunUMAP(seu, reduction = "NMF", dims=1:10,                reduction.name = "NMF_UMAP", reduction.key = "nmfUMAP_")


seu <- RunUMAP(seu, reduction = "pca", dims=1:10,                reduction.name = "PCA_UMAP", reduction.key = "pcaUMAP_")


library(ggplot2)DimPlot(seu, reduction = "PCA_UMAP", group.by = "celltype", label=T) +   theme(aspect.ratio = 1,        axis.text = element_blank(),        axis.title = element_blank(),        axis.ticks = element_blank()) + ggtitle("PCA UMAP") + NoLegend()


DimPlot(seu, reduction = "NMF_UMAP", group.by = "celltype", label=T) +   theme(aspect.ratio = 1,        axis.text = element_blank(),        axis.title = element_blank(),        axis.ticks = element_blank()) + ggtitle("NMF UMAP") + NoLegend()

GeneNMF除了在降維上替代PCA之外，還可以識(shí)別跨多個(gè)樣本的一致NMF程序。這里我們?cè)趤喨荷线M(jìn)行演示，對(duì)上皮細(xì)胞進(jìn)行亞群聚類：

#提取亞群，重新分群Epi <- uterus[,Idents(uterus) %in% c("Unciliated epithelial cells","Ciliated epithelial cells")]DefaultAssay(Epi) <- 'RNA'Epi <- NormalizeData(Epi, normalization.method = "LogNormalize", scale.factor = 1e4) Epi <- FindVariableFeatures(Epi, selection.method = 'vst', nfeatures = 3000)Epi <- ScaleData(Epi, vars.to.regress = c("S.Score", "G2M.Score","percent.mt", "nCount_RNA"), verbose = T)
Epi <- RunPCA(Epi, npcs = 40, verbose = FALSE)Epi <- RunHarmony(object = Epi,group.by.vars = c('orig.ident'))Epi <- RunUMAP(Epi,reduction = "harmony",dims = 1:10)Epi <- FindNeighbors(Epi,reduction = 'harmony', dims = 1:10)Epi <- FindClusters(Epi,resolution = 0.5, verbose = FALSE)
DimPlot(Epi, label = T, pt.size = 1)DimPlot(Epi, label = T, pt.size = 1, split.by  = 'orig.ident')

k是NMF程序關(guān)鍵的參數(shù)，為了獲取更加穩(wěn)定的programs，使用multimf（）函數(shù)自動(dòng)對(duì)樣本列表和多個(gè)k值執(zhí)行NMF。

seu.list <- SplitObject(Epi, split.by = "orig.ident")geneNMF.programs <- multiNMF(seu.list, assay="RNA",                              k=4:10,                              min.exp = 0.05,                             nfeatures=2000)
#將識(shí)別的基因程序合并geneNMF.metaprograms <- getMetaPrograms(geneNMF.programs,                                        metric = "cosine",                                        weight.explained = 0.5,                                        nMP=9,                                        max.genes=200)

#最后得到的結(jié)果geneNMF.metaprograms是一個(gè)list#check一些參數(shù)geneNMF.metaprograms$metaprograms.metrics# sampleCoverage silhouette meanSimilarity numberGenes numberPrograms# MP1      1.0000000 0.59783160          0.704          67             21# MP2      1.0000000 0.21113850          0.333          98             24# MP3      1.0000000 0.15745334          0.341          46             18# MP4      1.0000000 0.04416608          0.220          73             28# MP5      0.6666667 0.88466129          0.914          32             12# MP6      0.6666667 0.66182734          0.724          15             12# MP7      0.6666667 0.41217887          0.480          35             13# MP8      0.3333333 0.48668043          0.594          34              8# MP9      0.3333333 0.33003672          0.437         154             11

plot MP heatmap。

ht <- GeneNMF::plotMetaPrograms(geneNMF.metaprograms,similarity.cutoff = c(0.1,1))

這個(gè)默認(rèn)作圖也挺好的，但是我們?cè)谒幕A(chǔ)上修改修飾以下，讓它更接近高分文章中的狀態(tài)，為了表示對(duì)R包作者的尊重，函數(shù)還是用這個(gè)名稱！

pdf('./aaa.pdf', height = 8,width = 9)plotMetaPrograms(geneNMF.metaprograms,similarity.cutoff = c(0.1,1),                 palette = c('white','#bfd3e6','#9ebcda','#8c96c6',                             '#8c6bb1','#88419d','#810f7c','#4d004b'))dev.off()

獲取每個(gè)MP的基因：做富集分析，便于對(duì)基因程序的解釋！

MP_genes <- unlist(geneNMF.metaprograms$metaprograms.genes)MP <- rep(names(geneNMF.metaprograms$metaprograms.genes), lengths(geneNMF.metaprograms$metaprograms.genes))MP_genes <- data.frame(gene = MP_genes, MPs = MP)
#of cource, geneNMF提供了分析函數(shù)，可以直接用他的library(msigdbr)library(fgsea)top_p <- lapply(geneNMF.metaprograms$metaprograms.genes, function(program) {  runGSEA(program, universe=rownames(Epi), category = "C5", subcategory = "GO:BP")})
save(top_p, file = 'top_p.RData')

得到了MPs的基因特征，可以做基因集分?jǐn)?shù)：

mp.genes <- geneNMF.metaprograms$metaprograms.genesEpi <- AddModuleScore_UCell(Epi, features = mp.genes, ncores=4, name = "")VlnPlot(Epi, features=names(mp.genes), group.by = "orig.ident",pt.size = 0, ncol=3)VlnPlot(Epi, features=names(mp.genes), group.by = "seurat_clusters",pt.size = 0, ncol=3)#從樣本角度看，當(dāng)然我這個(gè)數(shù)據(jù)演示樣本沒(méi)有意義，有些MP在多個(gè)樣本中活躍，也有在某些樣本中突出活躍的MP#從亞群cluster的角度看，MP與cluster的活躍性顯著相關(guān)，例如MP1、MP5與cluster2明顯相關(guān)

# FeaturePlot(Epi, features = "MP1")# FeaturePlot(Epi, features = names(mp.genes),ncol=3) &#   scale_color_viridis(option="B") &#   theme(aspect.ratio = 1, axis.text=element_blank(), axis.ticks=element_blank())

integrated to seurat，使用基因程序得分表示單個(gè)細(xì)胞的特征！

#單個(gè)細(xì)胞現(xiàn)在可以用它們的基因程序得分來(lái)表示。matrix <- Epi@meta.data[,names(mp.genes)]
#dimred <- scale(matrix)dimred <- as.matrix(matrix)
colnames(dimred) <- paste0("MP_",seq(1, ncol(dimred)))#New dim reductionEpi@reductions[["MPsignatures"]] <- new("DimReduc",                                        cell.embeddings = dimred,                                        assay.used = "RNA",                                        key = "MP_",                                        global = FALSE)Epi <- RunUMAP(Epi, reduction="MPsignatures", dims=1:length(Epi@reductions[["MPsignatures"]]),               metric = "euclidean", reduction.name = "umap_MP")
library(viridis)FeaturePlot(Epi, features = names(mp.genes), reduction = "umap_MP", ncol=4) &  scale_color_viridis(option="B") &  theme(aspect.ratio = 1, axis.text=element_blank(), axis.ticks=element_blank())

DimPlot(Epi, reduction = "umap_MP", group.by = "seurat_clusters", label=T) + theme(aspect.ratio = 1,                                                                               axis.text = element_blank(),                                                                               axis.title = element_blank(),                                                                               axis.ticks = element_blank()) + ggtitle("Original cell types") + NoLegend()