作者，追风少年i
我们来复习复习吧。
2025年，希望霉运散去，一切向好。
前段时间有粉丝觉得NMF之前的代码都太复杂了，看看能不能来个简单版本的。
最近的文章也都不打算收费了，大家有个3块5块的打赏一下，买瓶绿茶喝一喝。
这一篇满足一下
library(Seurat)
library(NMF)
library(pheatmap)
library(RColorBrewer)
library(dplyr)
library(Matrix)
library(ggplot2)
library(pheatmap)
library(RColorBrewer)
library(viridis)
sdata = readRDS('sc.hsa.rds')
这里的rds我已经做过了基础分析，不过为了完整，再做一遍
sdata <- NormalizeData(sdata, normalization.method = "LogNormalize", scale.factor = 10000)
sdata <- FindVariableFeatures(sdata, selection.method = "vst")
sdata <- ScaleData(object = sdata, vars.to.regress = "percent_mito")
如果是多样本进行NMF呢？V5矩阵分开了，需要采用之前的合并策略,不过当下还是推荐V4那种合并成大矩阵的策略。
data <- as.matrix(sdata@assays$RNA@scale.data)
data[data<0] <- 0
运行NMF,一般需要提前选择数量，提取每个因子前50个基因作为特征值。
set.seed(123)
rk <- 2:20
res <- nmf(data, rank = rk, nrun = 20, seed=1)
plot(res)
这里采用了10
res <- nmf(data, rank = rk, nrun = 10, seed=1)
s <- extractFeatures(res, 50L)
plot(s)
可以看到获取的是基因的坐标值。
对每个特征的因子基因进行打分，获取NMF分数矩阵，这里采用了Seurat自带的函数AddModuleScore。
data1 <- c()
for (i in 1:as.numeric(rk)) {
data1 <- rbind(data1, as.matrix(sdata@assays$RNA@scale.data)[s[[i]],])
sdata <- AddModuleScore(object = sdata, features = list(rownames(data)[s[[i]]]), name = paste("NMF", i, sep = "_"))
}
my_palette<-colorRampPalette(c("dodgerblue3","white","Tomato2"))(n=19)
data1[data1>3] <- 3
data1[data1< -3] <- -3
x <- pheatmap(data1,cluster_cols=T, cluster_rows=F, color=my_palette, scale="none",
fontsize=5,treeheight_row=0,treeheight_col=0, cellheight = 1.2,cellwidth = 0.08,
show_rownames=F,show_colnames=F, clustering_method = "ward.D",
border_color=NA)
下面是文章的例图
选择特征值变化较大的因子
##### Filter programs ######
Pg_mat <- sdata@meta.data[ ,paste('Pg',"_",seq(1:rk), rep('1',rk), sep="")]
SD <- apply(Pg_mat, 2, sd)
SD
data2 <- c()
for (i in which(SD>0.1)) {
data2 <- rbind(data2, as.matrix(sdata@assays$RNA@scale.data)[s[[i]],])
}
data2[data2>3] <- 3
data2[data2< -3] <- -3
x <- pheatmap(data2,cluster_cols=T, cluster_rows=F, color=my_palette, scale="none",
fontsize=5,treeheight_row=0,treeheight_col=0, cellheight = 1.2,cellwidth = 0.08,
show_rownames=F,show_colnames=F, clustering_method = "ward.D",
border_color=NA)
expression program clustering,计算相关性系数
####### Jaccard Index #####
jac <- function(x, y) {
inter <- intersect(x, y)
total <- union(x, y)
similarity <- length(inter)/length(total)
return(similarity)
}
Mat <- matrix(0, ncol = ncol(Pgs), nrow = ncol(Pgs))
rownames(Mat) <- colnames(Pgs)
colnames(Mat) <- colnames(Pgs)
for (i in 1:ncol(Pgs)) {
for (j in 1:ncol(Pgs)) {
ss <- jac(as.vector(Pgs[,i]), as.vector(Pgs[,j]))
Mat[i,j] <- ss*100
}
}
custom_magma <- c(colorRampPalette(c("white", rev(magma(323, begin = 0.15))[1]))(10), rev(magma(323, begin = 0.18)))
#custom_magma<-colorRampPalette(c("white","blue"))(n=299)
Mat[Mat>40] <- 40
x <- pheatmap(as.matrix(Mat), cluster_cols=T, cluster_rows=T,
clustering_distance_rows="euclidean", color=custom_magma,
fontsize=2,treeheight_row=0,treeheight_col=30,
cellheight = 1.5,cellwidth = 1.5,show_rownames=T,
show_colnames=F, border_color = "NA")
其生物学意义，要根据自己的项目而定
生活很好，有你更好
