单细胞实战之单细胞非负矩阵分解(cNMF)——入门到进阶(高级篇4）

rm(list=ls())
options(stringsAsFactors = F) 
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(data.table)
library(dplyr)
library(qs)
library(BiocParallel)
library(stringr)
register(MulticoreParam(workers = 4, progressbar = TRUE)) 

#这里加载的是seurat对象，替换自己的数据即可
sc_data <- qread("./9-CD4+T/CD4+t_final.qs")
table(Idents(sc_data))

#检查一下自己导入进来的数据
Idents(sc_data) <- "celltype"
DimPlot(sc_data,reduction = 'umap',
        label = TRUE,pt.size = 0.5) +NoLegend()

dir.create("15-NMF")
setwd("15-NMF")

subdat <- GetAssayData(object = sc_data, slot = "counts")
#subdat <- ScaleData(subdat,do.center = F)  # 一定要有do.center = F，不然存在负值
subdat <- subdat[rowSums(subdat) > 0,]
subdat <- subdat[!rowSums(is.na(subdat)), ]
subdat <- subdat[, !colSums(is.na(subdat))]

# 可以按照自己的要求过滤一些基因，也可以不过滤
#subdat <- subdat[!str_detect(rownames(subdat),"^MT-"),]
#subdat <- subdat[!str_detect(rownames(subdat),"^Rp[sl]"),]
subdat <- as.data.frame(t(subdat)) # 转置成行为样本，列为基因的矩阵
write.table(subdat,file = "subdat.txt",quote = F,sep = "\t",
            row.names = T,col.names = T)

# 环境构建/cnmf包安装
conda create -n cNMF_env python 3.7
# 如果是M1/M2芯片的Mac
CONDA_SUBDIR=osx-64 conda create -n cNMF_env python=3.7
# 激活+安装
conda activate cNMF_env
pip install cnmf -i https://mirrors.aliyun.com/pypi/simple

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

# python代码运行
# prepare步骤 
cnmf prepare --output-dir ./res \ # 输出的文件夹为res
--name cNMF_res \ # 命名为cNMF_res
-c ./subdat.txt \ # 指定输入文件
-k 2 3 4 5 6 7 8 9 10 11 12 13 \ # 指定要尝试的因子数
--n-iter 100 --seed 123 # 每个k值运行100次，并设定种子数

# factorize步骤
cnmf factorize --output-dir ./res \ # 指定输出目录，与prepare阶段一致，确保读取之前准备好的数据
--name cNMF_res \ # 与prepare保持一致，告诉程序使用哪一组预处理结果
--worker-index 0 --total-workers 1 # 当前工作进程的编号，编号从0开始。同时单核运行

# combine步骤
cnmf combine --output-dir ./res --name cNMF_res # 执行合并步骤进行结果合并

# k_selection_图
cnmf k_selection_plot --output-dir ./res \ # 生成K值评估图，指定分析相同路径
--name cNMF_res # 

# 一致性聚类
cnmf consensus --output-dir ./res \ # 执行共识分析，读取结果
--name cNMF_res \ # 与前面保持一致
--components 4 \ # 最终选定的因子数 k=4
--local-density-threshold 2 \ # 调整聚类时的密度阈值
--show-clustering # 输出可视化结果

# 主要的文件是两个：
# Z分数单位基因表达程序矩阵和TPM单元基因表达程序矩阵
exprSet_score <- read.delim("./res/cNMF_res/cNMF_res.gene_spectra_score.k_4.dt_2_0.txt",header = T,row.names = 1)
# exprSet_tpm <- read.delim("./res/cNMF_res/cNMF_res.gene_spectra_tpm.k_4.dt_2_0.txt",header = T,row.names = 1)
rownames(exprSet_score) <- c("C1","C2","C3","C4")

exprSet <- exprSet_score
# 提取每一行的前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)

library(ComplexHeatmap)

#定义样本分组信息 
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
)

# 热图绘制
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",           
           "CD7",
           "ICOS",
           "SAT1",
           "FOXP3"
           )
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()

library(ClusterGVis)
library(org.Hs.eg.db)
# gather使得数据宽变长
# cluster：是新生成的列名，将存放转换后数据框的列名（作为变量的名称）
# gene：是另一个新生成的列名，将存放转换后数据框的值（即原来列中的数据）
# 1:ncol(top_genes)：指定要“收集”的列范围
cluster_gene <-gather(as.data.frame(top_genes), cluster, gene, 1:ncol(top_genes))

Gene_ID <- bitr(cluster_gene$gene, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")#构建文件并分析
data  <- merge(Gene_ID,cluster_gene,by.x='SYMBOL',by.y='gene')
data_GO <- compareCluster(
  ENTREZID~cluster, 
  data=data, 
  fun="enrichGO", 
  OrgDb="org.Hs.eg.db",
  ont = "BP",
  pAdjustMethod = "BH",
  pvalueCutoff = 0.05,
  qvalueCutoff = 0.05
)

#除去冗杂terms
cluster_GO <- simplify(data_GO, cutoff=0.7, by="p.adjust", select_fun=min)
dotplot(cluster_GO, showCategory=5,font.size = 8)
cluster_GO_res <- data_GO@compareClusterResult