单细胞数据中：NMF在Nature杂志文献中的新用法（张泽民院士团队）

NMF已被广泛应用于单细胞和空间表达数据的分析中，用于提取基因表达程序。但是在这篇来自张泽民院士团队的文献《Cross-tissue multicellular coordination and its rewiring in cancer》（2025年5月28号发表在Nature杂志）中，里面使用NMF做了一点不一样的应用，来看看！

在本研究中，作者开发了一个工具CoVarNet，其将NMF应用于频率矩阵以解析细胞共现程序，具体采用 nsNMF 方法，秩的范围为2至20，该功能由R语言的NMF包实现。为确保结果的稳健性，这篇研究进行了30次运行以获得共识输出，该输出包含k个因子及其在每个样本中的活性。具体而言，在泛组织分析中，每个因子中排名前十的细胞亚群被用作单个细胞模块网络中的共现节点候选。

特异性相关的亚群对：CoVarNet 利用皮尔逊相关系数来评估任意两个细胞亚群是否共现。

CoVarNet工作流程示意图：CoVarNet通过分析不同样本间细胞亚群频率的协变性来识别共现的细胞模块网络。

工具网址如下：

https://github.com/QiangShiPKU/CoVarNet?tab=readme-ov-file

文献应用

单细胞分辨率下的人类泛组织转录组图谱构建：

从35个人类组织的706份健康样本中，共获得2,293,951个高质量细胞用于下游分析（图1a，扩展数据图1及补充表1）。利用经典标志物，注释了八个主要的细胞类型（图1b），该图谱揭示了细胞组成中显著的跨组织变异性。

CoVarNet 通过利用细胞丰度的协变性，识别出12个具有不同细胞组成、组织偏好和空间组织的细胞模块。

软件安装

安装：

## 使用西湖大学的 Bioconductor镜像
options(BioC_mirror="https://mirrors.westlake.edu.cn/bioconductor")
options("repos"=c(CRAN="https://mirrors.westlake.edu.cn/CRAN/"))
devtools::install_github(repo = "https://github.com/QiangShiPKU/CoVarNet")
library(CoVarNet)

我这里使用的本地安装：

# bash命令本地安装
R CMD INSTALL -l /usr/local/software/miniconda3/envs/R4.4/lib/R/library CoVarNet-main

显示安装成功：

* installing *source* package ‘CoVarNet’ ...
** using staged installation
** R
** data
*** moving datasets to lazyload DB
** byte-compile and prepare package for lazy loading
✔ Setting active project to
"CoVarNet-main".
** help
*** installing help indices
** building package indices
** testing if installed package can be loaded from temporary location
** testing if installed package can be loaded from final location
** testing if installed package keeps a record of temporary installation path
* DONE (CoVarNet)

学习一个教程：https://htmlpreview.github.io/?https://github.com/QiangShiPKU/CoVarNet/blob/main/vignette/tutorial_discovery.html

step1：输入数据

输入数据为单细胞注释后的数据，该文件至少应包含三列："sampleID"（样本ID）、"majorCluster"（每个细胞所属的大类细胞类型）和"subCluster"（细胞所属的亚群）。

示例数据为上面文献的泛组织单细胞图谱，该数据可在data/example_discovery.rds中获取。网址如下：

https://github.com/QiangShiPKU/CoVarNet/tree/main/data

rm(list=ls())
library(CoVarNet)

meta <- readRDS("CoVarNet-main/data/example_discovery.rds")
dim(meta)
meta[1:5,]
table(meta$majorCluster)
table(meta$subCluster)
colnames(meta)

table(meta$cellSort)
meta <- meta[meta$cellSort %in% c("Mix", "Total"), ]  
rt_sp <- names(table(meta$sampleID))[table(meta$sampleID) >= 100]
meta <- meta[meta$sampleID %in% rt_sp, ]

上面包含了我们需要的三列：sampleID，majorCluster，subCluster，有7个大类，76个亚群，原始670个样本。

需要注意的点：

• 仅纳入未经分选的样本。这样就可以排除因实验性细胞类型富集操作而导致某些细胞类型缺失的样本。
• 过滤掉细胞数量少于 100 的样本。因为细胞数量不足可能导致细胞频率计算出现错误。

step2：计算细胞亚群频率矩阵

对于符合条件的样本，计算相应大类细胞类型中细胞亚群的频率。为了减轻不同大类细胞类型中细胞亚群数量差异所带来的影响，采用归一化方法对频率矩阵进行校正。

freq_calculate：会自动识别上面的三列 sampleID，majorCluster，subCluster

mat_fq_raw <- freq_calculate(meta)
mat_fq_norm <- freq_normalize(mat_fq_raw,normalize="minmax")
dim(mat_fq_norm)
mat_fq_norm[1:5,1:5]
colSums(mat_fq_raw[,-1])[1:5] 

得到的矩阵：行为76个亚群，列为670个样本，里面的值为每个样本中不同细胞亚群的相对比例。

step3：对频率矩阵进行 NMF 分解

非负矩阵分解（NMF）通过识别一组元亚群，并根据其权重对细胞亚群进行排序。为确定 NMF 分析的最优分解秩（参数 K），这里使用共表型相关系数（cophenetic correlation coefficient ）作为评估指标。将秩为 k 时的共表型相关系数记为 ρₖ，并针对当前情境建立了一套自动化确定最优秩的流程，其判定标准如下：

ρₖ₋₂ < ρₖ₋₁ < ρₖ ρₖ > ρₖ₊₁

注意：这一步骤运行很耗时间

其他参数（遵循 NMF R 包 v0.30.1 的规定）：

• method：指定 NMF 算法。
• seed：指定起始点或 seeding 方法。
• nrun：指定运行的次数。
• .options：设置运行时选项。

res <- nmf(mat_fq_norm, rank = 2:20, method = "nsNMF", seed = rep(123456, 6), .options = "vp")
save(res, file = "res.RData")
plot(res)

根据上述结果，此处设定分解秩rank K=12。

K=12
NMF_K12 <- nmf(mat_fq_norm, K, method = "nsNMF", seed = rep(77, 6), nrun = 30)

step4：CMs的特征可视化

非负矩阵分解将频率矩阵分解为：特征矩阵 Feature matrix（细胞亚群 × 细胞模块）：一系列根据权重对细胞亚群进行排序的细胞模块。系数矩阵 Coefficient matrix（细胞模块 × 样本）：每个样本中细胞模块的活性。

# NMF结果
NMF_K12@consensus[1:12,1:12]
dim(NMF_K12@consensus) # 行列都为样本
# 特征矩阵 Feature matrix（细胞亚群 × 细胞模块）：一系列根据权重对细胞亚群进行排序的细胞模块。
basis <- basis(NMF_K12) 
dim(basis )
head(basis )
# 系数矩阵 Coefficient matrix（细胞模块 × 样本）：每个样本中细胞模块的活性。
coef <- coef(NMF_K12) 
dim(coef)
coef[,1:10]

# 为确保与文章中的命名保持一致，我们在此对细胞模块进行了重命名。
colnames(basis(NMF_K12))=c("CM01","CM03","CM02","CM08","CM11","CM06","CM07","CM04","CM12","CM10","CM05","CM09")
rownames(coef(NMF_K12))=c("CM01","CM03","CM02","CM08","CM11","CM06","CM07","CM04","CM12","CM10","CM05","CM09")

热图：可视化CMs的潜在组成

展示了每个细胞模块内所有细胞亚群的权重，该权重基于 NMF 分解得到的特征矩阵获得。

## 可视化
pdf(file = "plot_NMF.heatmap.pdf",width = 18,height = 26)
gr.weight_all(NMF_K12)
dev.off()

这个图分辨率感觉被锁死了，比较模糊，看看后面有什么好办法没有。

还可以展示细胞模块内按权重排序的top细胞亚群。这些top亚群构成了细胞模块的潜在组成。

gr.weight_top(NMF_K12,num=15)

比如CM01模块主要由 MO9 Mph FOLR2组成。

可视化细胞模块在不同样本中的分布

系数矩阵被解释为细胞模块在样本中的活性。根据每个样本的主导细胞模块，所有样本被划分到互斥的细胞模块类型中。

h <- coef(NMF_K12)
max_cm <- apply(h, 2, function(x) rownames(h)[which.max(x)])
max_cm <- gsub("CM", "CMT",max_cm)
max_cm[1:10]
gr.distribution(NMF_K12,meta=meta,group="tissue") 

热图展示了细胞模块在不同组织中的活性图谱。每个样本被分配一个主导细胞模块，并且所有样本按其所属的细胞模块类型和组织类型进行排序。

step5：特异性相关的亚群对

CoVarNet 利用相关系数来评估任意两个细胞亚群是否特异性共现。在此，基于频率矩阵 mat_fq_raw 进行成对相关性检验，得到一个相关系数矩阵，记作 cor_pair。随后，我们计算亚群对的相关性显著性 pval_fdr 及其特异性 spe。对于 cor_pair 中每个满足 i<j 的元素 rᵢⱼ，其背景集定义为 {rᵢₖ, rₖⱼ | k ≠ i, j}。特异性定义为背景集中不超过 rᵢⱼ 的元素所占的比例。

此处计算所有细胞亚群对之间的相关性，并将在下一节中筛选出特异性相关的亚群对。

cor_pair <- pair_correlation(mat_fq_raw,method="pearson")
saveRDS(cor_pair,"./cor_pair.rds")
head(cor_pair)

通过特异性相关边连接共现节点以构建细胞模块网络

在每个细胞模块网络中，权重最高的细胞亚群被视为潜在节点，特异性相关的亚群对作为边进行连接，并移除孤立节点。

具体细节如下：

首先，需要确定每个细胞模块中细胞亚群数量的上限（top_n）。

其次，根据以下标准筛选特异性相关的亚群对作为边：

• 相关系数 > 参数 corr（默认值：0.2）
• 校正后p值 pval_fdr < 参数 fdr（默认值：0.05）
• 特异性 spe > 1 - ((top_n-1) * 2 - 1) / ((n_all-1) * 2 - 1)。其中 n_all 是所有细胞亚群的总数。

最后，利用 R 语言的 igraph 包为每个细胞模块构建网络。

top_n = 10
network <- cm_network(NMF_K12,cor_pair,top_n,corr=0.2,fdr=0.05)
each <- network$each
gr.igraph_each(each,Layout=layout_in_circle) 

节点为所属的大类细胞类型，边的颜色为则按其特异性成比例缩放。

step6：保存结果

最后输出每个CMs的组成并保存结果：

head(network$filter)
saveRDS(network,"./network.rds")

今天分享到这里~

友情转发：