非负矩阵分解（NMF）的python应用 ,基因分析，以胰腺癌上皮亚簇为实例，NMF在癌症研究中的优势

在提供的两篇研究文献中，非负矩阵分解（NMF）的应用主要与胰腺癌的分子亚型分析相关，具体计算逻辑结合文献内容及NMF基本原理总结如下：

一、文献中的NMF应用场景

在第二篇文献《Integrative analysis of spatial and single-cell transcriptome data from human pancreatic cancer》中，提到NMF用于分析胰腺癌的bulk转录组数据，以识别分子亚型。例如：

• 文中指出胰腺癌上皮亚簇 Ep_PIFO 的标记基因与先前研究中通过NMF鉴定的 Signature 10 高度相关，该Signature与纤毛功能相关。
• NMF在bulk数据中通过分解基因表达矩阵，将胰腺癌分为不同分子亚型（如经典型、基底样型等），为单细胞亚簇的功能注释提供参考。

二、NMF的核心计算原理（结合文献场景）

1. 数据输入与目标

• 输入数据：胰腺癌bulk转录组的基因表达矩阵 ( \mathbf{V} )（行=基因，列=样本，元素为非负表达值）。
• 分解目标：将 ( \mathbf{V} ) 近似分解为两个非负矩阵的乘积，即：
  [
  \mathbf{V} \approx \mathbf{W} \times \mathbf{H}
  ]
  其中，( \mathbf{W} ) 为 基矩阵（行=基因，列=亚型，反映各亚型的特征基因表达模式），( \mathbf{H} ) 为 系数矩阵（行=亚型，列=样本，反映每个样本在各亚型上的权重）。

2. 非负性约束的意义

• 基因表达值本质为非负数，NMF的非负约束确保分解结果符合生物学意义（如亚型特征基因的表达量不会为负）。
• 相比PCA等方法，NMF更适合捕捉基因表达的“组合性”特征（如不同亚型的基因共表达模块）。

3. 迭代优化过程

• 初始化：随机生成非负矩阵 ( \mathbf{W} ) 和 ( \mathbf{H} )。
• 迭代更新：通过优化算法（如乘法更新规则、梯度下降）最小化重构误差，常用损失函数为：
  [
  \text{Loss} = |\mathbf{V} - \mathbf{W}\mathbf{H}|F^2 \quad (\text{ Frobenius范数})
  ]
  每次迭代中，按以下规则更新矩阵元素（保持非负）：
  [
  W{ik} \leftarrow W_{ik} \times \frac{(\mathbf{V}\mathbf{H}^{T)_{ik}}{(\mathbf{W}\mathbf{H}\mathbf{H}}T){ik}}, \quad H{kj} \leftarrow H_{kj} \times \frac{(\mathbf{W}^{T\mathbf{V})_{kj}}{(\mathbf{W}}T\mathbf{W}\mathbf{H})_{kj}}
  ]
• 终止条件：损失函数收敛或达到预设迭代次数。

4. 亚型识别与解释

• 分解后，( \mathbf{H} ) 中权重最高的亚型即代表样本的主要分子类型（如经典型、基底样型）。
• 通过分析 ( \mathbf{W} ) 中高权重基因的功能，可注释各亚型的生物学特征（如Ep_PIFO关联的纤毛功能相关基因）。

三、文献中NMF与单细胞分析的整合

• 验证单细胞亚簇的生物学意义：通过NMF在bulk数据中定义的亚型特征，可验证单细胞亚簇是否对应已知分子亚型。例如，Ep_PIFO簇的标记基因（如与纤毛相关的基因）在NMF的Signature 10中高度富集，提示其功能特异性。
• 桥接单细胞与bulk数据：NMF分解得到的亚型权重（( \mathbf{H} )）可与单细胞亚簇的比例进行关联分析，揭示肿瘤内异质性与整体亚型的关系。

四、NMF在癌症研究中的优势

1. 鲁棒性：对高维、稀疏的基因表达数据具有较好的降维效果，且不受噪声影响。
2. 可解释性：分解结果直接对应生物学上的“亚型”或“模块”，便于功能注释。
3. 兼容性：可与单细胞分析结合，从bulk水平验证单细胞亚簇的临床意义（如预后相关性）。

总结

文献中NMF的计算核心是通过非负矩阵分解将bulk转录组数据降维为生物学可解释的分子亚型，其迭代优化过程确保分解结果符合基因表达的非负性和功能模块性。该方法在胰腺癌研究中用于定义分子亚型，并与单细胞亚簇的功能注释和预后分析相结合，为解析肿瘤异质性提供了关键工具。

在单细胞转录组和癌症研究中，非负矩阵分解（NMF）常用于分析基因表达数据以识别分子亚型或细胞状态。以下结合文献中的应用场景，提供Python实现NMF的完整流程，重点关注其在胰腺癌和食管癌研究中的应用逻辑。

一、Python实现NMF的核心库与数据准备

1. 关键库导入

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.decomposition import NMF
from sklearn.preprocessing import normalize
from sklearn.metrics import silhouette_score
import seaborn as sns
import scanpy as sc  # 单细胞数据处理

2. 数据准备（以胰腺癌bulk转录组为例）

假设输入数据为基因表达矩阵（行=基因，列=样本），非负且已标准化：

# 加载bulk转录组数据（示例数据，实际需替换为真实数据）
# 数据格式：行为基因，列为样本，值为log2(TPM+1)
expression_matrix = pd.read_csv("pancreatic_cancer_expression.csv", index_col=0)# 确保数据非负（NMF的基本要求）
if (expression_matrix < 0).any().any():expression_matrix = expression_matrix - expression_matrix.min().min()print("数据已调整为非负矩阵")# 示例：若数据为单细胞矩阵，需先提取非免疫细胞或上皮细胞群体
# adata = sc.read_h5ad("single_cell_data.h5ad")
# non_immune_cells = adata[adata.obs['cell_type'] != 'immune_cell']
# expression_matrix = non_immune_cells.X.toarray()

二、NMF模型构建与参数优化

1. 模型训练与分解

def run_nmf(matrix, n_components, random_state=42):"""执行NMF分解并返回基矩阵和系数矩阵"""nmf = NMF(n_components=n_components,init='nndsvd',       # 推荐初始化方法，提高收敛速度solver='mu',          # 乘法更新规则，适用于大规模数据random_state=random_state,max_iter=1000,       # 最大迭代次数tol=1e-4             # 收敛阈值)W = nmf.fit_transform(matrix)  # 系数矩阵 (样本×亚型)H = nmf.components_         # 基矩阵 (亚型×基因)return W, H, nmf# 示例：分解胰腺癌表达矩阵为k=2种亚型（经典型+基底样型）
W, H, nmf = run_nmf(expression_matrix, n_components=2)# 计算重构误差（评估分解质量）
reconstruction_error = nmf.reconstruction_err_
print(f"NMF重构误差: {reconstruction_error}")

2. 确定最佳亚型数量（k值）

# 方法1：基于重构误差随k的变化
def find_optimal_k(matrix, max_k=10):errors = []for k in range(1, max_k+1):_, _, nmf = run_nmf(matrix, k)errors.append(nmf.reconstruction_err_)plt.figure(figsize=(8, 6))plt.plot(range(1, max_k+1), errors, 'o-')plt.xlabel('亚型数量 (k)')plt.ylabel('重构误差')plt.title('NMF重构误差随k的变化')plt.grid(True)plt.show()# 肘部法则：选择误差下降趋缓的k值elbow_k = np.argmin(np.diff(errors) > 0) + 1return elbow_koptimal_k = find_optimal_k(expression_matrix)
print(f"最佳亚型数量: {optimal_k}")# 方法2：结合生物学先验（如文献中已知胰腺癌分为2-3种亚型）
optimal_k = 2  # 假设根据文献预设k=2

三、NMF结果分析与生物学解释

1. 解析基矩阵（H）：识别亚型特征基因

def identify_subtype_markers(H, gene_names, n_top_genes=20):"""从基矩阵中提取各亚型的 top 标记基因"""subtype_markers = {}for i in range(H.shape[0]):# 按基因权重排序gene_indices = np.argsort(H[i, :])[::-1]top_genes = gene_names[gene_indices][:n_top_genes]subtype_markers[f"Subtype_{i+1}"] = top_genesreturn subtype_markers# 示例：假设gene_names为表达矩阵的行索引
gene_names = expression_matrix.index.values
subtype_markers = identify_subtype_markers(H, gene_names)# 可视化亚型特征基因热图
plt.figure(figsize=(12, 8))
top_genes = np.unique([gene for markers in subtype_markers.values() for gene in markers])
heatmap_data = expression_matrix.loc[top_genes].Tsns.heatmap(heatmap_data,cmap="virdis",annot=False,fmt=".1f",cbar_kws={"label": "标准化表达量"}
)
plt.title("胰腺癌NMF亚型特征基因热图")
plt.show()

2. 解析系数矩阵（W）：样本亚型分类

# 将样本分配到主要亚型
sample_subtypes = np.argmax(W, axis=1)
subtype_proportions = pd.Series(sample_subtypes).value_counts(normalize=True)
print("各亚型样本比例:", subtype_proportions)# 结合临床数据分析（如文献中的肿瘤分期、生存数据）
# 假设clinical_data为包含样本临床信息的DataFrame
# clinical_data['NMF_Subtype'] = sample_subtypes
# 后续可进行生存分析、亚型与临床特征的相关性分析

四、NMF与单细胞数据的整合分析

1. 单细胞亚簇与NMF亚型的关联

def correlate_sc_subclusters_with_nmf(sc_expression,      # 单细胞基因表达矩阵 (细胞×基因)nmf_h,              # NMF基矩阵 (亚型×基因)sc_subclusters      # 单细胞亚簇标签
):"""计算单细胞亚簇与NMF亚型的相关性"""# 计算单细胞亚簇的平均基因表达cluster_avg_expr = {}for cluster in np.unique(sc_subclusters):cluster_cells = sc_expression[sc_subclusters == cluster]cluster_avg_expr[cluster] = np.mean(cluster_cells, axis=0)# 转换为DataFrame便于计算cluster_avg_df = pd.DataFrame(cluster_avg_expr).Tcluster_avg_df.columns = sc_expression.columns# 计算亚簇与NMF亚型的相关性correlation_matrix = pd.DataFrame()for i in range(nmf_h.shape[0]):subtype_genes = np.argsort(nmf_h[i, :])[::-1][:50]  # 取亚型前50个特征基因subtype_expr = cluster_avg_df[subtype_genes].mean(axis=1)correlation_matrix[f"Subtype_{i+1}"] = subtype_expr# 计算相关性系数correlation_coeffs = pd.DataFrame(np.corrcoef(correlation_matrix.T),index=correlation_matrix.index,columns=correlation_matrix.index)return correlation_coeffs# 示例：假设sc_data为单细胞表达矩阵，sc_clusters为细胞亚簇标签
# correlation = correlate_sc_subclusters_with_nmf(sc_data, H, sc_clusters)
# sns.heatmap(correlation, annot=True, cmap="coolwarm")
# plt.title("单细胞亚簇与NMF亚型的相关性")

2. 验证单细胞亚簇与NMF亚型的一致性

以文献中胰腺癌的Ep_PIFO亚簇为例，验证其与NMF的Signature 10相关性：

# 假设ep_pifo_genes为Ep_PIFO亚簇的标记基因
ep_pifo_genes = ["PIFO", "KRT6B", "TRIM54", ...]  # 替换为实际标记基因# 计算Ep_PIFO基因在NMF亚型中的富集程度
pifo_gene_indices = [np.where(gene_names == gene)[0][0] for gene in ep_pifo_genes if gene in gene_names]
pifo_enrichment = np.mean(H[:, pifo_gene_indices], axis=1)# 可视化Ep_PIFO与各亚型的关联
plt.bar(range(1, len(pifo_enrichment)+1), pifo_enrichment)
plt.xlabel("NMF亚型")
plt.ylabel("Ep_PIFO基因平均权重")
plt.title("Ep_PIFO亚簇标记基因在NMF亚型中的富集")
plt.show()

五、进阶应用：结合空间转录组数据

1. NMF亚型的空间分布分析

def analyze_nmf_spatial_distribution(spatial_expression,  # 空间转录组表达矩阵 (spot×基因)nmf_h,               # NMF基矩阵spot_coordinates     # spot的空间坐标 (x, y)
):"""分析NMF亚型在空间上的分布"""# 计算每个spot的亚型权重spot_subtype_weights = NMF(n_components=nmf_h.shape[0]).fit_transform(spatial_expression)# 转换为DataFrame并添加坐标spot_subtype_df = pd.DataFrame(spot_subtype_weights, columns=[f"Subtype_{i+1}" for i in range(nmf_h.shape[0])])spot_subtype_df["x"] = spot_coordinates[:, 0]spot_subtype_df["y"] = spot_coordinates[:, 1]# 可视化亚型空间分布（以Subtype_1为例）plt.figure(figsize=(10, 8))plt.scatter(spot_subtype_df["x"], spot_subtype_df["y"], c=spot_subtype_df["Subtype_1"], cmap="viridis",s=50,alpha=0.8)plt.colorbar(label="Subtype_1权重")plt.title("NMF亚型在空间转录组中的分布")plt.show()

六、代码应用场景与文献关联

1. 胰腺癌分子亚型分析
   在第二篇文献中，NMF用于将胰腺癌分为经典型和基底样型，对应单细胞亚簇Ep_TRIM54（经典型）和Ep_KRT6A（基底样型）。通过上述代码可复现：

   • 分解bulk转录组数据，识别亚型特征基因（如经典型的GATA6，基底样型的KRT6B）。
   • 验证单细胞亚簇的标记基因是否在NMF亚型中富集（如Ep_PIFO与纤毛相关Signature 10的关联）。

2. 食管癌免疫微环境分析
   第一篇文献中虽未明确提及NMF，但可类比应用：

   • 对食管癌bulk转录组进行NMF分解，识别与免疫抑制相关的亚型（如富含Treg、TEx细胞的亚型）。
   • 结合单细胞数据，分析NMF亚型与免疫细胞亚群（如tDC、Treg）的空间共定位关系。

七、注意事项

1. 数据预处理：

   • 确保输入数据非负，若为原始计数数据，建议使用log(CPM+1)或标准化处理。
   • 对于高维数据，可先通过PCA降维（保留80%方差）再进行NMF，提高计算效率。

2. 生物学验证：

   • NMF亚型需结合标记基因的功能注释（如GO/KEGG富集分析），避免纯数学分解导致的生物学误解。
   • 参考文献中的已知亚型（如胰腺癌的经典-基底样二分法），调整k值以匹配生物学先验。

3. 单细胞-空间整合：

   • 若需复现文献中Ep_VGLL1亚簇的过渡特性，可通过NMF亚型权重与单细胞亚簇比例的相关性分析，验证其“中间状态”。

通过上述代码框架，可实现从NMF模型构建到生物学解释的全流程分析，紧密结合两篇文献中胰腺癌和食管癌的研究场景，为肿瘤微环境的分子亚型解析提供量化工具。