单细胞非负矩阵分解分析python版(cNMF)学习

2024-09-15 22:35:37 浏览数 (1)

前置一个推文，老师的推文已经详细讲解了非负矩阵分解的算法原理~ 如果对算法原理感兴趣的可以点击以下链接~

单细胞天地： https://mp.weixin.qq.com/s/-sdYyBG_zB6Lhi9vHkpKBw

笔者之前也写过一个帖子，有兴趣的朋友可以点击去看一看~ https://mp.weixin.qq.com/s/3zySnfkflHfitqh4p4chsQ

接下来笔者会根据个人理解，将非负矩阵分解(Non-negative Matrix Factorization，NMF)、一致性聚类(consensus clustering)、主成分分析(Principal Component Analysis，PCA) 做一个的类比解释。

(硬蹭一下黑神话悟空/西游记hhhh~)

小时候的我们看了很多遍西游记，但那时候懵懵懂懂，看到悟空的时候我们一般会按照主观印象对它进行评价，首先可能会说这是只猴子，接下来可能会说这只猴子长的挺有趣的，最后可能会说这只猴子还挺聪明，这种强烈按照主要、次要等有秩序的特征评价方式就类似于PCA的方法。

后来随着我们长大，不再以主观第一印象去评价事物了。这时候我们再评价悟空的时候评价方式可能会略有不同，我们可能还是会说悟空是只猴子，但我们也会提到打妖怪的时候每次都悟空去打，以及悟空的沟通能力都值得我们所有人学习等等，这些评价的话语是看了很多集的故事而总结出来的，没有十分显著的描述词汇，但诉说的语句又都是正确。这种方式就类似于一致性聚类的方法，它通过了频繁的抽样把矩阵中的信息分成多个聚类，这些聚类内部是非常稳定的，不同聚类之间互相独立，组合在一起可以完整的描述矩阵的特色，但是每一个聚类不存在十分显著的代表特征(相对于非负矩阵分解的方法) 。

再后来人到了中年，我们的知识水平又到了进一步的提高。这时候我们再评价悟空的时候又有了新的变化。我们会说悟空是兼具了英勇气概，谋略才智，强大力量，仁爱正义等特征的一个人物。这些词汇都是对悟空这个角色特点的高度概括，单独一个词汇就能代表悟空的一部分但都无法去准确完整的描述悟空。这种方式就类似于非负矩阵分解的方法，同样非负矩阵分解会把矩阵分成不同的聚类，但此时提取的聚类是按照这个聚类群的关键特征所提取，具有高度的“代表性”。这些具有代表性的特征互相之间是“平等”的，都可以代表矩阵，但是单独一个特征又无法完整的说明整个矩阵的特点，此时的聚类内部可能没有像一致性聚类分析(相对而言)那样稳定，不同聚类之间也是相互独立。

以上就是三种方式的类比解释啦~

接下来的分析流程采用了python版的非负矩阵分解-cNMF

github上的图挺形象的，从左到右把细胞信息转化成基因矩阵然后进行分解提取。

步骤流程
R语言部分

1、导入

代码语言：javascript复制
rm(list=ls())
options(stringsAsFactors = F) 
library(Seurat)
library(ggplot2)
library(qs)
library(BiocParallel)
register(MulticoreParam(workers = 4, progressbar = TRUE)) 
sc_data <- qread("./sub_data.qs")
dim(sc_data)
# [1] 28269  1164

一共有28269个基因，1164个细胞

2、提取数据

代码语言：javascript复制
subdat <- as.data.frame(sc_data@assays$RNA@counts)
subdat <- subdat[rowSums(subdat) > 0,]
subdat <- subdat[!str_detect(rownames(subdat),"^MT-"),]
subdat <- as.data.frame(t(subdat)) # 转置成行为样本，列为基因的矩阵
write.table(subdat,file = "subdat.txt",quote = F,sep = "t",
            row.names = T,col.names = T)

保存成一个txt文件，其中去除了线粒体基因，同时使用者也可以根据自己需求去除一些基因信息。

python部分

1、环境部署

代码语言：javascript复制
# 环境构建/cnmf包安装
conda create -n cNMF_env python 3.7
conda activate cNMF_env
pip install cnmf -i https://mirrors.aliyun.com/pypi/simple

# 接下来需要使用cd进入工作文件夹

2、cNMF运行

./res是分析导出的文件夹名

cNMF_res分析之后数据文件夹名称

subdat.txt是纳入的分析文件

k后面的数字是不同的聚类数，需要提前预设

--n-iter后面的数字代表每个指定的 k 值（即不同成分个数）下运行的迭代次数为 100 次。

--worker-index 0 --total-workers 1 ，这里的workers代表多少线程，index后面的数字代表数据分成几部分运行，如果选择了workers为5，那么需要重复运行factorize步骤5次(index后面的数字要从0开始到4结束)

--components 5 , 这里面的数字是根据组分数量与稳定性和误差的关系图得到的

--local-density-threshold 2 这里的2是不设定阈值，这要根据聚类图右上角的密度图而定。

代码语言：javascript复制
# python代码运行
# prepare步骤 
cnmf prepare --output-dir ./res 
--name cNMF_res 
-c ./subdat.txt 
-k 2 3 4 5 6 7 8 9 10 11 12 13 
--n-iter 100 --seed 123

# factorize步骤
cnmf factorize --output-dir ./res 
--name cNMF_res 
--worker-index 0 --total-workers 1

# combine步骤
cnmf combine --output-dir ./res --name cNMF_res

# k_selection_图步骤
cnmf k_selection_plot --output-dir ./res 
--name cNMF_res

# 一致性聚类
cnmf consensus --output-dir ./res 
--name cNMF_res 
--components 4 
--local-density-threshold 2 
--show-clustering

这个图展示了组分数量与稳定性和误差的关系, 一般选择下降最快的曲线的前一个点，这里是选择(该图产生于k_selection_图步骤)。

聚类图，可以看出中间聚类和边上的颜色色差很大，说明聚类效果很好。右上角的密度图也没有异常值。

github示例数据，此时右上角的密度图在0值之后出现了一些小的柱状图，这时候就需要设定阈值进行过滤。

比如按照这个图应该是把--local-density-threshold 设置在了0.15左右。

再次进入R语言部分

本次暂时只需要用到cNMF_res.gene_spectra_score.k_4.dt_2_0.txt 文件,重点关注含有score字样的文件。

与其他老师使用分解样本之后再聚类的方式不同，笔者并不将样本分解而是直接对细胞群进行NMF处理。

如果按照其他老师的方法还需要关注含有consensus字样的文件。

1、将文件数据进行预处理

代码语言：javascript复制
exprSet <- read.delim("cNMF_res.gene_spectra_score.k_4.dt_2_0.txt",header = T,row.names = 1)
rownames(exprSet) <- c("C1","C2","C3","C4")

# 提取每一行的前50个基因
top_genes <- apply(exprSet, 1, function(x) {
  names(sort(x, decreasing = TRUE)[1:50])
})

expr <- t(exprSet)
# 转换为矩阵用于绘图
top_genes_matrix <- sapply(top_genes, function(genes) {
  expr[genes, ]
})
total_genes <- t(top_genes_matrix)
# scale data
total_genes <- scale(total_genes)

2、Complexheatmap绘图可视化

代码语言：javascript复制
#定义样本分组信息 
clusters <- colnames(total_genes)  

# 定义列注释
col_ha = HeatmapAnnotation(
  Clusters = clusters,
  simple_anno_size = unit(2, 'mm'),
  col = list(Clusters = c('C1' = '#1F77B4',  # 蓝色
                       'C2' = '#FF7F0E',  # 橙色
                       'C3' = '#2CA02C',  # 绿色
                       'C4' = '#D62728')  # 红色
  ),
  show_annotation_name = FALSE
)

# 热图绘制
library(ComplexHeatmap)# 设置热图参数
p <- Heatmap(
  as.matrix(total_genes),
  cluster_rows = FALSE,  # 行不聚类
  cluster_columns = FALSE,  # 列不聚类
  show_row_names = FALSE,  # 不显示行名
  show_column_names = TRUE,  # 显示列名
  top_annotation = col_ha,  # 添加列注释
  heatmap_legend_param = list(title = "Score"),  # 热图的图例标题
  col = colorRampPalette(c("#313695", "#4575b4", "white", "firebrick3", "#a50026"))(100)
)
p

# 添加注释信息
genes <- c("S100A10",           
           "CD1B",
           "IGFBP7",
           "WBP5",
           "RAMP1"
           )
genes <- as.data.frame(genes)

# 设置输出设备，指定文件名和尺寸
pdf("heatmap_with_annotations.pdf", width = 6, height = 8)
# 添加基因信息
draw(p   rowAnnotation(link = anno_mark(at = which(rownames(total_genes) %in% genes$genes),                                          labels = genes$genes, 
                                   labels_gp = gpar(fontsize = 10)))
)
# 关闭设备
dev.off()
参考资料：

1、cNMF：

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

2、Top生物信息：

https://mp.weixin.qq.com/s/GRim4kuvZmCxrxnosLreag https://mp.weixin.qq.com/s/GjtHyWP11ttPd7PJuKsiPQ

3、生信技能树：

https://mp.weixin.qq.com/s/RQCRUolM9HpZSodam6B8dQ

致谢：感谢曾老师以及生信技能树团队全体成员。

注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台(欢迎交流)。更多内容可关注公众号：生信方舟

- END -

数据分析单细胞 非负矩阵分解

0 人点赞