细胞的基因表达谱决定了细胞的功能和细胞的状态，在单细胞研究中，我们一般通过做差异基因分析来寻找我们感兴趣的疾病特异性表达的基因，从而研究疾病的发生发展。但差异基因分析存在统计假设检验的大前提，所以可能会在某些情况下筛选不到我们感兴趣的疾病特异性表达程序（比如样本数较多、样本分组较多、异质性较大，可能很难通过差异基因的方法观观察到多个样本、分组之间共有的基因表达程序）。

cNMF（Consensus Non-negative Matrix factorization ）是一个用于从单细胞 RNA-Seq (scRNA-Seq) 数据推断基因表达程序的分析管道。它采用count矩阵(N 细胞 * G 基因) 作为输入，并生成一个 (K x G)基因表达程序 (GEPs)矩阵和一个(N x K) 在指定数据中每个细胞的每个程序的使用情况矩阵。

已经有非常多的高影响因子文章通过cNMF来研究特定状态下（如疾病、发育阶段等）的基因表达程序，比如2022年7月发表在《Nature Genetics》题为Single-nucleus and spatial transcriptome profiling of pancreatic cancer identifies multicellular dynamics associated with neoadjuvant treatment的文章中，研究者使用共识非负矩阵分解（cNMF）从头学习了跨不同肿瘤的恶性细胞和CAF的反复表达程序。通过稳定性和错误选择cNMF程序的数量，重点关注来自多个患者的细胞之间共享的那些，并用其前200个进行注释加权基因。最终确定了4个反映细胞谱系（经典细胞、鳞状细胞、基底细胞样、间充质、腺泡样、神经内分泌样和神经样祖细胞（NRP））或细胞状态（cycling-S、cycling-G2/M、MYC、干扰素、肿瘤坏死因子）的恶性细胞程序/核因子κB（TNF-NFκB），核糖体和粘合剂）和四个CAF程序（肌成纤维细胞祖细胞，嗜神经细胞，免疫调节剂和粘合剂）。

安装

cNMF需要 scikit-learn的版本大于1.0，并且依赖scanpy。安装完这些依赖包后可使用pip安装cNMF。

pip install cnmf

使用方法

cNMF的运行分为5个步骤：

步骤 1 ：- 标准化输入矩阵并准备运行参数

cnmf prepare --output-dir ./example_data --name example_cNMF \
             -c ./example_data/counts_prefiltered.txt \
             -k 5 6 7 8 9 10 11 12 13 \
             --n-iter 100 --seed 14 --numgenes 2000

路径结构：

• --output-dir - 所有结果将被放置到的输出目录。默认: .，即当前目录。
• --name - 将创建子目录 output_dir/name ，并且所有输出文件都将以 name 作为前缀。默认: cNMF

输入数据：

• -c - 细胞 x 基因count文件的路径。这应该是以制表符分隔的文本文件或以 h5ad 格式保存的 Scanpy 对象
• --tpm [可选] - 预先计算TPM标准化数据。如果未提供，将自动计算TPM标准化数据。如果需要特定的标准化，这会很有帮助。这些可以以与count文件相同的格式加载。默认: None
• --genes-file [可选] - 用于分解步骤的过度分散基因列表。如果没有提供，将自动计算过度分散的基因，并且可以通过下面的 --numgenes 参数设置要使用的基因数量。默认: None

参数：

• -k - 将测试 cNMF 的 K 值列表，以空格分割
• --n-iter - 为每个 K 运行的 NMF 迭代次数。默认: 100
• --total-workers - 指定可以并行使用的工作线程数量，默认: 1
• --seed - 将用于为每个 NMF 冲服生成单独种子的主种子，默认: None
• --numgenes - 将用于运行分解的最高方差基因的数量。去除低方差基因有助于放大信号，是正确推断数据中程序的重要因素。然而，不用担心，最后程序谱会重新拟合以包含所有基因的估计值，甚至包括那些未包含在高方差集中的基因。默认: 2000
• --beta-loss - NMF 的损失函数, 可选 frobenius, kullback-leibler, itakura-saito中的一个，默认: frobenius
• --densify -- 输入为密集矩阵的标识符，即输入数据并不稀疏，不建议用于大多数单细胞 RNA-Seq 数据，特别是10x平台的数据。默认: False

这个命令生成一个经过过滤和归一化的矩阵，用于运行矩阵分解。它首先将数据划分为一组过度分散的基因，这些基因可以作为输入文件提供或在此处计算。虽然将为输入count文件中的所有基因计算最终程序谱，但如果仅在一组高方差基因上运行，则分解速度要快得多，并且可以找到更好的模式。还可以提供每个细胞归一化的输入文件，以便可以根据该归一化来计算最终的基因表达程序。

此外，此命令还将要运行的特定分解作业分配给不同的workers（通过--total-workers参数指定 ）。

在上面的示例中，假设不使用并行化（--total-workers 1），因此所有作业都分配给单个workers。

注意：输入矩阵不应包含任何总计数为 0 的细胞或基因。此外，如果任何细胞最终的过度分散基因计数为 0，则可能会导致错误。需要在运行 cNMF 之前过滤掉计数较低的细胞和基因。

步骤 2：分解矩阵

接下来的 cNMF 将对上一个命令中指定的所有任务按workers索引重复运行。任务已分配给索引从 0 开始的workers，如下所示，对于对应于第一个工作程序的索引 0：

cnmf factorize --output-dir ./example_data --name example_cNMF --worker-index 0 

如果指定的workers总数超过 1 个，则需要使用单独的命令为这些workers运行命令，例如：

cnmf factorize --output-dir ./example_data --name example_cNMF \
               --worker-index 0 --total-workers 3
               
cnmf factorize --output-dir ./example_data --name example_cNMF \
               --worker-index 1 --total-workers 3
               
cnmf factorize --output-dir ./example_data --name example_cNMF \
               --worker-index 2 --total-workers 3

建议将这些命令提交给不同的处理器或机器，以便它们全部并行运行。

建议将这些命令提交给不同的处理器或机器，以便它们全部并行运行。

提示：如果一台机器上有多个核心，则默认情况下 scikit-learn 中 NMF 的实现将使用超过 1 个核心。实际应用中，在使用 GNU 并行时可能使用 2 个工作线程可以获得最佳性能。

步骤 3 ：将每个 K 的各个谱结果文件合并为一个合并文件

由于已为每个 K 的每个迭代创建了一个单独的文件，因此需要将每个 K 的迭代文件合并起来，如下所示：

cnmf combine --output-dir ./example_data --name example_cNMF

之后，可以选择删除各个谱文件，如下所示：

rm ./example_data/example_cNMF/cnmf_tmp/example_cNMF.spectra.k_*.iter_*.df.npz

步骤 4 ：通过考虑稳定性和误差之间的权衡来选择最佳K

这将迭代所有已运行的 K 值并计算稳定性和误差，然后，它输出一个 PNG 图像文件，将这种关系绘制到 output_dir/name 目录中 示例命令：

cnmf k_selection_plot --output-dir ./example_data --name example_cNMF

这会将 K 选择图输出到 example_data/example_cNMF/example_cNMF.k_selection.png。选择 K 没有普遍明确的标准，但建议使用相当稳定的最大值和/或局部稳定性最大值。

步骤 5：在所需的 K 值下获得程序及其使用的共识估计

最后一步是在首先选择性地滤除异常值后对谱进行聚类，此步骤最终输出 4 个文件：TPM单元的GEP 估计，Z-scores单元的GEP 估计（反映程序的较高使用率是否会减少或增加基因表达），非标准化 GEP 使用估计， Clustergram 诊断图（显示每次迭代结果之间的一致性程度和每个谱与其K个最近邻之间的距离直方图）

建议使用诊断图来确定过滤异常值的阈值。默认情况下，cNMF 将用于此过滤的邻居数量设置为已完成迭代次数的 30%，但这可以从命令行修改。

在实践中，建议运行此命令两次，一次使用 --local- Density-threshold 2.00 来查看平均距离的分布情况，然后第二次使用 --local- Density-threshold 设置为较小的值根据该直方图确定过滤异常值。

命令示例:

cnmf consensus --output-dir ./example_data --name example_cNMF \
               --components 10 --local-density-threshold 0.2 \
               --show-clustering

• --components -用于计算共识集群的 K 值。必须包含在准备步骤提供的选项中
• --local-density-threshold - 到 K 个最近邻居的平均距离的阈值，2.0或以上意味着不会过滤掉任何内容。默认值：0.5
• --local-neighborhood-size - 被视为局部密度过滤的最近邻居的重复百分比。例如，如果您运行 100 次重复，并将其设置为 0.3，则将使用 30 个最近邻居进行异常值检测。默认值：0.3
• --show-clustering - 控制是否输出 clustergram 图像。默认值: False

输出文件

输出目录中应该有以下结果文件：

• Z-score单元基因表达程序矩阵：example_data/example_cNMF/example_cNMF.gene_spectra_score.k_10.dt_0_01.txt
• TPM 单元基因表达程序矩阵 ：example_data/example_cNMF/example_cNMF.gene_spectra_tpm.k_10.dt_0_01.txt
• 使用矩阵：example_data/example_cNMF/example_cNMF.usages.k_10.dt_0_01.consensus.txt
• K 选择图：- example_data/example_cNMF/example_cNMF.k_selection.png
• 诊断图：- example_data/example_cNMF/example_cNMF.clustering.k_10.dt_0_01.pdf

NMF输出的基因功能集合结果可以用于各种功能分析。例如Jaccard相似性分析，比较基因功能集合与细胞类型的相似性，应用于癌症数据中可以帮助判断癌细胞的细胞来源以及研究癌细胞的异质性。另外，将功能基因集合进行富集分析（如GO和GEVA）可以更好地解读基因功能集合所对应的功能，对后期寻找疾病相关的标志基因提供更多的可能性。

参考：https://github.com/dylkot/cNMF/blob/master/Tutorials/analyze_simulated_example_data.ipynb

https://github.com/dylkot/cNMF/tree/master?tab=readme-ov-file