R版本的NMF分析中，我們介紹了相關(guān)概念，可查看（（視頻教程）NMF: 單細(xì)胞非負(fù)矩陣分解R語(yǔ)言版-數(shù)據(jù)降維及MetaProgram分析）。但是R語(yǔ)言版本有一個(gè)嚴(yán)重缺陷就是太耗費(fèi)時(shí)間了，而cnmf卻很快，成了替代的有力工具！

1-cnmf簡(jiǎn)介
cNMF（Consensus NMF） 是一個(gè)從單細(xì)胞 RNA 測(cè)序（scRNA-Seq）數(shù)據(jù)中推斷基因表達(dá)程序的流程。以一個(gè)計(jì)數(shù)矩陣（N 個(gè)細(xì)胞×G 個(gè)基因）作為輸入，并生成一個(gè)（K×G）的基因表達(dá)程序（GEPs）矩陣和一個(gè)（N×K）矩陣，該矩陣指明了數(shù)據(jù)中每個(gè)細(xì)胞對(duì)每個(gè)程序的使用情況。cNMF是基于python的分析工具，輸入要求是單細(xì)胞RNA測(cè)序數(shù)據(jù),可以應(yīng)用于任何類型的單細(xì)胞數(shù)據(jù),包括非腫瘤樣本。官網(wǎng)有詳細(xì)的教程，原理感興趣的可以閱讀paper。本貼僅為數(shù)據(jù)演示，無(wú)任何生物學(xué)意義。視頻解說(shuō)教程已上傳微信VIP，請(qǐng)自行下載！

github：https://hub.whtrys.space/dylkot/cNMF
原文paper：https://elifesciences.org/articles/43803

2-創(chuàng)建分析環(huán)境

#終端中運(yùn)行：創(chuàng)建分析環(huán)境（最好單獨(dú)創(chuàng)建一個(gè)環(huán)境，避免庫(kù)的沖突）
conda create -n cnmf_env python=3.10
conda activate cnmf_env
#安裝cNMF
pip install cnmf -i https://pypi.tuna.tsinghua.edu.cn/simple
#Successfully installed anndata-0.11.4 array-api-compat-1.11.2 cnmf-1.7.0 contourpy-1.3.2 cycler-0.12.1 exceptiongroup-1.2.2 
#fonttools-4.57.0 h5py-3.13.0 joblib-1.5.0 kiwisolver-1.4.8 legacy-api-wrap-1.4.1 llvmlite-0.44.0 matplotlib-3.10.1 natsort-8.4.0 
#networkx-3.4.2 numba-0.61.2 numpy-2.2.5 packaging-25.0 palettable-3.3.3 pandas-2.2.3 patsy-1.0.1 pillow-11.2.1 pynndescent-0.5.13 
#pyparsing-3.2.3 python-dateutil-2.9.0.post0 pytz-2025.2 pyyaml-6.0.2 scanpy-1.11.1 scikit-learn-1.5.2 scipy-1.15.2 seaborn-0.13.2 
#session-info2-0.1.2 six-1.17.0 statsmodels-0.14.4 threadpoolctl-3.6.0 tqdm-4.67.1 typing-extensions-4.13.2 tzdata-2025.2 umap-learn-0.5.7
#也可以運(yùn)行一下help，看看是否正常(因?yàn)槲抑爸苯影惭b在別的環(huán)境，help出錯(cuò)，說(shuō)明環(huán)境和包有沖突，所以最好單獨(dú)創(chuàng)建環(huán)境)
cnmf --help

3-Prepare for cNMF analysis
如果你的數(shù)據(jù)是scanpy分析的，直接將需要分析的對(duì)象保存為h5ad文件即可。如果你的數(shù)據(jù)是R語(yǔ)言seurat分析的，提取表達(dá)矩陣即可，保存為txt用于分析。

#在R中：
setwd('./data_analysis/NMF_singlecell/')
library(Seurat)
epi_count <-  GetAssayData(Epi, slot = 'counts') %>% as.data.frame()
epi_count <- epi_count[rowSums(epi_count) > 0,]#將表達(dá)量全部為0的基因去除
epi_count <- as.data.frame(t(epi_count)) # 轉(zhuǎn)置數(shù)據(jù)框，行為細(xì)胞，列為基因
write.table(epi_count,file = "epi_count.txt",quote = F,sep = "\t",row.names = T,col.names = T)
#scanpy
count_adat_fn = './counts.h5ad'
sc.write(count_adat_fn, Epi_dat)

4-cNMF分析步驟（python環(huán)境）
前期準(zhǔn)備工作完成之后，就可以安心進(jìn)行分析了，這里我們先演示python環(huán)境中分析，數(shù)據(jù)就直接使用adata h5ad文件了。cnmf的分析總的來(lái)說(shuō)有5大步驟，對(duì)應(yīng)五個(gè)函數(shù)：
prepare：運(yùn)行準(zhǔn)備；
factorize：分解矩陣； ·
combine：將多個(gè)k值迭代結(jié)果合并； ·
k_selection_plot：選擇最佳k值；
·consensus：聚類及獲取結(jié)果

#導(dǎo)入庫(kù)
from cnmf import cNMF
cnmf_obj = cNMF(output_dir="./cnmf_test/", name="Epi_cnmf")#初始化cnm obj 。
#output_dir：結(jié)果保存文件夾,默認(rèn)是當(dāng)前文件夾，這里我們是當(dāng)前文件夾下的cnmf_test文件夾
#name:此次分析的名稱，所有文件前綴，避免多個(gè)分析混淆！
#step1：
cnmf_obj.prepare(counts_fn='./counts.h5ad', components=[2,3,4,5,6,7,8,9,10,11,12], n_iter=20, seed=14, num_highvar_genes=2000)
#counts_fn: 矩陣文件，可以是txt制表符分割文件也可以是h5ad文件
#components: k值，一個(gè)數(shù)組，這個(gè)一般設(shè)置多個(gè)，類似于R中方nmf，多少自行設(shè)置
#step2：
cnmf_obj.factorize(worker_i=0, total_workers=1)#total_workers 并行線程數(shù)，cnm快得很，所以1個(gè)work通常就夠了
#step3：
#將每個(gè)k值得迭代結(jié)果合并
cnmf_obj.combine()
#step4：
#計(jì)算每個(gè)選擇K的穩(wěn)定性和誤差，選擇合適的k值。
cnmf_obj.k_selection_plot()

image.png

#step5：
#（正常情況下，step5運(yùn)行兩次）
#針對(duì)給定的 K 值計(jì)算共識(shí)解。首先可以在不進(jìn)行任何異常值過(guò)濾的情況下運(yùn)行它，以查看結(jié)果如何。
cnmf_obj.consensus(k=4, density_threshold=2.00, show_clustering=True, close_clustergram_fig=False)

image.png

#根據(jù)上面的圖，一致性不錯(cuò)，但還是有一些異常值，根據(jù)直方圖過(guò)濾，這里設(shè)置density_threshold為0.1，可以過(guò)濾最大異常值
cnmf_obj.consensus(k=4, density_threshold=0.10, show_clustering=True, close_clustergram_fig=False)

image.png

5-輸出結(jié)果
輸出的結(jié)果文件都在我們第一步prepare的output_dir文件夾
最終結(jié)果文件如下:(k=4,density_threshold=0.10)
文件名也是有特點(diǎn)的，cnmf name+。。。+k+dt。。。
Z-score unit gene expression program matrix : Epi_cnmf.gene_spectra_score.k_4.dt_0_1.txt
TPM unit gene expression program matrix： Epi_cnmf.gene_spectra_tpm.k_4.dt_0_1.txt
Usage matrix： Epi_cnmf.usages.k_4.dt_0_1.consensus.txt
另外還有兩個(gè)圖，一個(gè)是選擇最佳k值折線圖，一個(gè)是聚類驗(yàn)證圖。
我們可以看看結(jié)果是什么形式的：

gene_spectra_score = pd.read_csv('./cnmf_test/Epi_cnmf/Epi_cnmf.gene_spectra_score.k_4.dt_0_1.txt',sep='\t',index_col=0)
gene_spectra_score
圖片
gene_spectra_tpm = pd.read_csv('./cnmf_test/Epi_cnmf/Epi_cnmf.gene_spectra_tpm.k_4.dt_0_1.txt',sep='\t',index_col=0)
gene_spectra_tpm
#每個(gè)GEP中基因權(quán)重
圖片
Epi_cnmf_usages = pd.read_csv('./cnmf_test/Epi_cnmf/Epi_cnmf.usages.k_4.dt_0_1.consensus.txt',sep='\t',index_col=0)
#每個(gè)細(xì)胞中，每個(gè)GEP的強(qiáng)弱
Epi_cnmf_usages.head()
圖片
可以在UMAP中可視化GEP，如果你的數(shù)據(jù)是scanpy，那么就可以直接將Epi_cnmf_usages與adata合并可視化。
Epi_dat.obs = pd.merge(left=Epi_dat.obs, right=Epi_cnmf_usages, how='left', left_index=True, right_index=True)
sc.pl.umap(Epi_dat, color=Epi_cnmf_usages.columns,ncols=3, vmin=0, vmax=1)

image.png

cnmf的運(yùn)行得到了和我們之前講過(guò)的的R nmf相對(duì)應(yīng)的結(jié)果（也就是分解得到的兩個(gè)矩陣）。針對(duì)這個(gè)結(jié)果，可能想做一些具體的應(yīng)用，比如說(shuō)降維，比如說(shuō)metaprograme等等?？傊甤nmf還是很快的，終于不用怕NMF分析費(fèi)時(shí)間的窘境了。

覺(jué)得我們分享有些用的，點(diǎn)個(gè)贊再走唄！