脚本更新---NMF识别单细胞数据中的基因模块

原创

追风少年i

发布于 2024-12-03 08:00:41

36300

代码可运行

运行总次数：0

代码可运行

作者，Evil Genius

今日更新脚本，NMF在单细胞数据中的分析

NMF在单细胞分析中的主要运用

降维和特征提取：NMF可以用于从高维的基因表达矩阵中提取低维特征，帮助我们发现细胞亚群、基因模块等隐藏的结构。

细胞亚群的识别：通过对单细胞数据进行NMF分解，可以识别细胞的潜在亚群，并将基因表达模式与这些亚群关联起来。

基因模块的发现：NMF可以用于识别在单细胞数据中表达模式相似的基因，这些基因可能在生物学上有共同的功能或调控机制。

数据的噪声过滤和降维：NMF通过分解矩阵的方式提取低秩近似，可以有效地去除数据中的噪声，保留数据的主要结构。

NMF在单细胞分析中的具体应用

1. 单细胞数据的降维与特征提取

NMF的一个主要用途是将高维的单细胞RNA-seq数据分解为低秩近似。通过将基因表达矩阵V分解为W（基矩阵）和H（系数矩阵），可以识别出隐藏在数据中的主要模式：

W代表基因的特征模式。
H代表细胞的特征表示。

这种方式不仅能够降维，还能帮助我们识别哪些基因在特定细胞类型或状态下表现突出。

library(NMF)

# 假设 scRNA-seq 数据是一个基因（行）- 细胞（列）的矩阵
set.seed(123)
V <- matrix(abs(rnorm(1000)), nrow = 100, ncol = 10)  # 假设100个基因，10个细胞

# 使用NMF对数据进行分解
nmf_results <- nmf(V, rank = 3, nrun = 10)

# 提取W和H矩阵
W <- basis(nmf_results)  # 基因模式
H <- coef(nmf_results)   # 细胞模式

# 可视化W和H
heatmap(W, main = "Gene Feature Matrix (W)")
heatmap(H, main = "Cell Feature Matrix (H)")

通过对W和H的解读，研究人员可以理解哪些基因在不同细胞亚群中具有高表达，并且可以进一步分析这些基因是否与特定的生物学功能或路径有关。

2. 细胞亚群的识别

NMF可以用于识别细胞群体。在单细胞数据中，细胞通常具有不同的状态或类型，而这些类型或状态的差异可以通过NMF分解的H矩阵（细胞模式）来揭示。

例如，通过将H矩阵（细胞特征矩阵）进行聚类，可以发现不同的细胞亚群。

# 假设H是细胞的特征表示
# 可以使用层次聚类或K-means聚类来识别细胞亚群
library(pheatmap)
library(cluster)

# 聚类细胞模式H
cell_clusters <- kmeans(t(H), centers = 3)  # 假设我们希望分成3个细胞亚群

# 绘制热图
pheatmap(t(H), cluster_rows = TRUE, cluster_cols = TRUE, annotation_col = as.data.frame(cell_clusters$cluster))

这样，H矩阵的列可以对应不同的细胞群体，而基于这些群体的表达模式，可以进一步分析细胞类型的特征。

3. 基因模块的发现

NMF也可以用来发现基因模块，即在不同细胞群体或条件下共同表达的基因集合。通过查看W矩阵的行，可以识别在多个细胞群体中共同表达的基因。

# 假设W是基因的特征矩阵
# 你可以识别哪些基因在不同的细胞群体中表现出高表达
# 提取具有高表达模式的基因
top_genes <- apply(W, 1, function(x) mean(x) > 0.5)
genes_of_interest <- rownames(W)[top_genes]

# 打印出这些基因
print(genes_of_interest)

这些基因可能在不同的细胞类型或状态下有相似的表达模式，可能代表具有相似生物学功能的基因。

4. 通过NMF去除噪声并提高数据质量

NMF通过对原始数据的低秩近似分解，可以有效地去除噪声并保留数据的主要特征。在处理单细胞RNA-seq数据时，噪声和技术变异可能会影响分析结果，而NMF能够提供更稳定和一致的特征。

例如，通过NMF分解，细胞和基因的模式会被提取并重建，从而减少不相关的变异。

# 重构数据
reconstructed_V <- W %*% H

# 可视化重构后的数据与原始数据的比较
heatmap(V, main = "Original Data")
heatmap(reconstructed_V, main = "Reconstructed Data")

这种重建能够帮助去除一些与生物学无关的技术噪声，使得下游分析（如聚类、差异表达分析）更加准确。

5. 单细胞转录组数据的多样性分析

NMF还可以用于研究细胞的多样性和变异性。在分析不同细胞群体的基因表达时，NMF能够揭示细胞群体之间的相似性和差异性，并帮助发现潜在的细胞类型或转录状态。

例如，可以使用consensus_nmf方法将多次NMF分解的结果结合起来，得到更加稳定的细胞群体和基因模块。

# 假设我们已经使用多次NMF分解
nmf_results <- nmf(V, rank = 3, nrun = 50)  # 多次分解以获得共识

# 计算共识矩阵
consensus_matrix <- consensus(nmf_results)

# 可视化共识矩阵
heatmap(consensus_matrix, main = "Consensus Matrix")

通过共识NMF方法，可以获得更加稳定和一致的分解结果，从而提高数据分析的可靠性。

NMF方法在实际项目中的运用

我们首先要理解整体的架构

（1）选择要分析的基因，一般我们选择高变基因（或者高度变异的前多少个基因，一般5000为主）
（2）矩阵的选择，一般要选择经过scale后的矩阵，并且将负值设置为零来调整居中的表达式矩阵。
（3）conducted NMF runs with various K values (K = 5,6,7,8,9,10), resulting in 45 NMF programs for each sample.
（4）识别共识转录程序
（5）Spearman correlation was performed to determine the interrelationships between different programs.
（6）基于阈值0.8对相关矩阵进行二值化，并对二值化后的矩阵进行层次聚类分析。将不同K运行中相似的程序进行合并，计算合并后程序中细胞和基因的平均归一化系数。根据每个合并程序的系数得分，选择前150个基因作为特征基因。为了识别不同样本中反复出现的NMF程序（“meta程序”），计算Jaccard相似性矩阵，该矩阵表示不同程序的特征基因的重叠。

我们在实际中运用一下，大家要分析好单细胞数据，注释好，获取rds文件。

library(NMF)
library(Seurat)
library(plyr)
library(dplyr)
library(future)
library(gtools)
library(ggplot2)
library(cowplot)
library(data.table)
library(tidyverse)
library(future.apply)
library(RColorBrewer)

获取scale的矩阵

# NMF ======================================================================

# generate input data --------------------

nmf_input <- lapply(allexpmc, function(i){
  
  data <- NormalizeData(i)
  data <- FindVariableFeatures(data, nfeatures = 7000)
  data <- ScaleData(data, do.center = F)
  
  as.matrix(GetAssayData(data,
            slot = "scale.data", assay = "RNA"))
})

识别program，合并并选择特征基因

# 安装和加载NMF包
install.packages("NMF")
library(NMF)

# 生成示例数据
set.seed(123)
V <- data

# 定义K值列表
K_values <- 5:10

# 存储所有K值的W和H矩阵
W_list <- list()
H_list <- list()

# 执行不同K值的NMF分解并计算W和H矩阵
for (K in K_values) {
  cat("Running NMF for K =", K, "\n")
  
  # 执行NMF分解，nrun=10进行10次分解以增加稳定性
  nmf_result <- nmf(V, rank = K, nrun = 10)
  
  # 提取W和H矩阵（基因模式和细胞模式）
  W <- basis(nmf_result)  # 基因模式
  H <- coef(nmf_result)   # 细胞模式
  
  # 将W和H矩阵保存到列表中
  W_list[[as.character(K)]] <- W
  H_list[[as.character(K)]] <- H
}

# 合并所有K值的W和H矩阵，计算平均归一化系数
W_avg <- matrix(0, nrow = nrow(V), ncol = length(K_values))  # 初始化一个矩阵来存储W的平均系数
H_avg <- matrix(0, nrow = ncol(V), ncol = length(K_values))  # 初始化一个矩阵来存储H的平均系数

# 计算每个基因和细胞的平均系数
for (i in 1:length(K_values)) {
  W_avg[, i] <- rowMeans(W_list[[i]], na.rm = TRUE)  # 基因的平均系数
  H_avg[, i] <- rowMeans(H_list[[i]], na.rm = TRUE)  # 细胞的平均系数
}

# 对W和H矩阵进行归一化
W_avg_norm <- scale(W_avg, center = FALSE, scale = apply(W_avg, 2, max))  # 按列最大值归一化
H_avg_norm <- scale(H_avg, center = FALSE, scale = apply(H_avg, 2, max))  # 按列最大值归一化

# 计算每个基因的得分，并选择前150个基因
gene_scores <- rowSums(W_avg_norm)
top_genes_idx <- order(gene_scores, decreasing = TRUE)[1:150]  # 选择得分最高的150个基因

# 查看选择的前150个基因的名称
gene_names <- paste("Gene", 1:nrow(V))  # 假设基因名称为Gene1, Gene2, ..., Gene100
top_genes <- gene_names[top_genes_idx]
print(top_genes)

拿到consensus_matrix，分析signature relationship

# 基因特征矩阵（W_avg_norm）的热图
pheatmap(W_avg_norm, 
         main = "Gene Feature Matrix (W) - Average Normalized Coefficients",
         scale = "none",  # 不对数据进行缩放，因为我们已经进行了归一化
         cluster_rows = TRUE,  # 对基因进行聚类
         cluster_cols = TRUE,  # 对K值进行聚类
         show_rownames = TRUE, 
         show_colnames = TRUE,
         color = colorRampPalette(c("blue", "white", "red"))(50))  # 设置颜色梯度

# 细胞特征矩阵（H_avg_norm）的热图
pheatmap(H_avg_norm, 
         main = "Cell Feature Matrix (H) - Average Normalized Coefficients",
         scale = "none",  # 不对数据进行缩放
         cluster_rows = TRUE,  # 对细胞进行聚类
         cluster_cols = TRUE,  # 对K值进行聚类
         show_rownames = TRUE,
         show_colnames = TRUE,
         color = colorRampPalette(c("blue", "white", "red"))(50))  # 设置颜色梯度

####选择前150个基因并进行热图展示
# 选择前150个基因的表达矩阵（W_avg_norm）进行热图展示
W_top150 <- W_avg_norm[top_genes_idx, ]  # 选择前150个基因对应的行

# 绘制选出的前150个基因的热图
pheatmap(W_top150, 
         main = "Top 150 Genes - Average Normalized Coefficients",
         scale = "none",  # 不进行额外的缩放
         cluster_rows = TRUE,  # 对基因进行聚类
         cluster_cols = TRUE,  # 对K值进行聚类
         show_rownames = TRUE,
         show_colnames = TRUE,
         color = colorRampPalette(c("blue", "white", "red"))(50))  # 设置颜色梯度