scMetabolism分析细胞代谢

包安装

https://github.com/wu-yc/scMetabolism

第一步 分离对比对象

1
2
3
4
HSPC <- subset(immune, group == 'HSPC')
CPRC <- subset(immune, group == 'CRPC')
countexp_HSPC <-sc.metabolism.Seurat(obj = HSPC, method = "VISION", imputation = F, ncores = 12, metabolism.type = "KEGG")
countexp_CPRC <-sc.metabolism.Seurat(obj = CPRC, method = "VISION", imputation = F, ncores = 12, metabolism.type = "KEGG")

第二步 获得代谢矩阵

1
2
etabolism.matrix_HSPC <- countexp_HSPC@assays$METABOLISM$score
etabolism.matrix_CPRC <- countexp_CPRC@assays$METABOLISM$score

第三步 代谢差异分析

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
f_DEmetabolism_01 <- function(matrixN, matrixE){
res <- data.frame(row.names = rownames(matrixN))
res[['meanN']] = 0
res[['meanE']] = 0
res[['p-value']] = 1
for(i in 1:nrow(matrixN)){
test <- t.test(unlist(matrixN[i,]),unlist(matrixE[i,]))
res[i, 'p-value'] = test$p.value
res[i, 'meanN'] = test$estimate[1]
res[i, 'meanE'] = test$estimate[2]
}
res[['def']] <- with(res, meanE - meanN)
res
}
f_metaN2cellN <- function(scRNA, metaN, typeN, matrixN=F){
res <- scRNA[[metaN]]
idx <- which(res[[1]] == typeN)
if(matrixN){
return(gsub('-','.',rownames(res)[idx]))
}else{
return(rownames(res)[idx])
}
}
f_DEmetabolism <- function(countexpN, countexpE, metaN, typeN){
matrixN <- countexpN@assays$METABOLISM$score
matrixE <- countexpE@assays$METABOLISM$score
matrixN <- matrixN[,f_metaN2cellN(countexpN, metaN, typeN, matrixN = T)]
matrixE <- matrixE[,f_metaN2cellN(countexpE, metaN, typeN, matrixN = T)]
f_DEmetabolism_01(matrixN, matrixE)
}

第四步 代谢差异展示

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
require(ggrepel)
f_ggplot2_ti <- function(p, title, subtitle=''){
if (subtitle == ''){
(p + ggtitle(title) + theme(plot.title = element_text(hjust = 0.5)))
}else{
(p + ggtitle(title, subtitle=subtitle) +
theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5)))
}
}
f_DEmetabolism_Volcano <- function(countexpN, countexpE, metaN, typeN, delta = 1, Threshold_p = 0.05, lc_rep=5){
DEmm <- f_DEmetabolism(countexpN, countexpE, metaN, typeN)
col_vector = rep(rgb(108, 200, 228, maxColorValue = 255), nrow(DEmm))
Threshold = with(DEmm, delta*sd(abs(def)))

col_vector[with(DEmm, `p-value` < Threshold_p & def > Threshold)] = rgb(226, 61, 75, maxColorValue = 255)
col_vector[with(DEmm, `p-value` < Threshold_p & def < -Threshold)] = rgb(232, 168, 71, maxColorValue = 255)

DEmm[['Name']] <- rownames(DEmm)
DEmm[['col']] <- col_vector
DEmm[['-log10(P)']] <- with(DEmm, -log10(`p-value`))

lc_tp_logFC <- DEmm$def
lc_tp_logFC[with(DEmm, `p-value` >= Threshold_p)] = 0
lc_idx <- order(lc_tp_logFC)[c(1:lc_rep, (nrow(DEmm)+1-lc_rep):nrow(DEmm))]

res <- ggplot() + geom_point(aes(def, `-log10(P)`, col=I(col)), data = DEmm)
res <- res + theme_bw() + theme(panel.grid=element_line(colour=NA))
res <- res + geom_hline(yintercept=-log10(Threshold_p), linetype="longdash")
res <- res + geom_vline(xintercept=c(Threshold, -Threshold), linetype="longdash")
res <- res + geom_text_repel(data=DEmm[lc_idx,],aes(def,`-log10(P)`,label=Name), force=T, max.overlaps=Inf)
f_ggplot2_ti(res, typeN)
}

第五步 全部展示

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
require(export)
f_ppt_output <- function(fileName, image){
for (p in image){
graph2ppt(x=p, file=fileName, width=9, aspectr=1.5, append = TRUE)
}
}
pl <- list()
for (tN in unique(unlist(countexp_HSPC[['immune_type_2']]))){
pl[[tN]] <- f_DEmetabolism_Volcano(countexp_HSPC, countexp_CPRC, 'immune_type_2', tN)
}
f_ppt_output('HSPC vs CRPC.pptx', pl)

input.pathway<- rownames(countexp_HSPC@assays$METABOLISM$score)
p1 <- DotPlot.metabolism(obj = countexp_HSPC, pathway = input.pathway, phenotype = "immune_type_2", norm = "y") + DotPlot.metabolism(obj = countexp_CPRC, pathway = input.pathway, phenotype = "immune_type_2", norm = "y")
ggsave(p1, filename = 'HSPC vs CRPC.pdf', dpi = 1200, width = 24, height = 36, device = 'pdf', limitsize = FALSE)

第六步 基因热图-获取代谢相关基因

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
require(KEGGREST)
require(stringr)
f_KEGG_name2id <- function(keggL, nameL){
res <- NULL
for (i in 1:length(nameL)){
fuck <- gsub(pattern = '\\(', replacement = '\\\\(', x = nameL[[i]])
fuck <- gsub(pattern = '\\)', replacement = '\\\\)', x = fuck)
idx <- str_detect(keggL, fuck)
tmp <- keggL[idx]
if (length(tmp) == 1){
tmp <- names(tmp)
tmp <- substr(x =tmp, start = 6, stop=15)
res <- c(res, tmp)
}
if (length(tmp) > 1){
print(tmp)
}
}
res
}
f_KEGG_id2symbol <- function(KEGGid){
res <- NULL
for (hsa_id in KEGGid){
hsa_info <- keggGet(hsa_id)
hsa_info <- hsa_info[[1]]$GENE
hsa_info <- hsa_info[seq(from = 2, to = length(hsa_info), by = 2)]
hsa_info <- strsplit(hsa_info,split = ';')
hsa_info <- unlist(lapply(hsa_info, function(X){X[1]}))
res <- c(res, hsa_info)
}
res <- unique(res)
res <- res[str_detect(res, pattern = '\\] \\[', negate = T)]
res
}
MMPath <- keggGet('hsa01100')
MMPath <- names(MMPath[[1]]$REL_PATHWAY)
MMgene <- f_KEGG_id2symbol(MMPath)

第六步 绘制基因热图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
f_metaG2G <- function(metaG, matrixN=F){
res <- list()
alltype <- unique(metaG[[1]])
for(type in alltype){
res[[type]] <- rownames(metaG)[metaG[[1]] == type]
if (matrixN){
res[[type]] <- gsub('-','.',res[[type]])
}
}
res
}
f_matrix_groupMean <- function(matrixA, group, matrixN = T){
res <- data.frame(row.names = rownames(matrixA))
group <- f_metaG2G(group, matrixN = matrixN)
for(name in names(group)){
res[[name]] <- rowMeans(matrixA[,group[[name]]])
}
res
}
require(reshape2)
f_matrix_heatmap <- function(dfA){
# 标准化一下
# dfA <- as.data.frame(scale(dfA))
# 转换前,先增加一列ID列,保存行名字
dfA$df_ID <- rownames(dfA)
dfm <- melt(dfA, na.rm = T, id.vars = c('df_ID'))
p <- ggplot(dfm, aes(x=variable, y=df_ID))
p <- p + geom_tile(aes(fill=value))
p <- p + theme(axis.text.x=element_text(angle=45, hjust=1, vjust=1))
p <- p + scale_fill_gradient(low = "blue", high = "pink")
p <- p + xlab("samples") + theme_bw() + theme(panel.grid.major = element_blank()) + theme(legend.key=element_blank())
p
}

mono <- subset(immune, immune_type_2 %in% c('TAM','monocytes', 'S100A8+ monocytes'))
mono <- sc.metabolism.Seurat(obj = mono, method = "VISION", imputation = F, ncores = 12, metabolism.type = "KEGG")
a <- f_matrix_groupMean(mono@assays$METABOLISM$score, mono[['patient_id']])
p <- f_matrix_heatmap(a)
ggsave(p, filename = 'mono_HSPC vs CRPC_heatmap.pdf', dpi = 1200, width = 8, height = 24, device = 'pdf', limitsize = FALSE)

第七步 绘制交叉热图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
f_metaG2G <- function(metaG, matrixN=F){
res <- list()
alltype <- unique(metaG[[1]])
for(type in alltype){
res[[type]] <- rownames(metaG)[metaG[[1]] == type]
if (matrixN){
res[[type]] <- gsub('-','.',res[[type]])
}
}
res
}
f_matrix_groupMean <- function(matrixA, group, matrixN = T){
res <- data.frame(row.names = rownames(matrixA))
group <- f_metaG2G(group, matrixN = matrixN)
for(name in names(group)){
if (length(group[[name]]) == 1){
res[[name]] <- matrixA[,group[[name]]]
}else{
res[[name]] <- rowMeans(matrixA[,group[[name]]])
}
}
res
}
require(reshape2)
f_matrix_heatmap <- function(dfA){
# 标准化一下
# dfA <- as.data.frame(scale(dfA))
# 转换前,先增加一列ID列,保存行名字
dfA$df_ID <- rownames(dfA)
# dfA$df_ID <- factor(x = dfA$df_ID, ordered = T)
dfm <- melt(dfA, na.rm = T, id.vars = c('df_ID'))
dfm$variable <- factor(x = as.character(dfm$variable), ordered = T)
p <- ggplot(dfm, aes(x=variable, y=df_ID))
p <- p + geom_tile(aes(fill=value))
p <- p + scale_fill_gradient(low = "white", high = "red")
p <- p + xlab("samples") + theme_bw() + theme(panel.grid.major = element_blank()) + theme(legend.key=element_blank())
p <- p + theme(axis.text.x=element_text(angle=45, hjust=1, vjust=1))
p
}
1
2
3
4
5
immune[['tmp2']] <- paste(as.character(immune[['immune_type_2']][[1]]), as.character(immune[['patient_id']][[1]]))

tmp2 <- f_matrix_groupMean(immune@assays$METABOLISM$score[unlist(sR),], immune[['tmp2']])
p2 <- f_matrix_heatmap(tmp2)
ggsave(p2, filename = 'HSPC vs CRPC_heatmap_2.pdf', dpi = 1200, width = 24, height = 24, device = 'pdf', limitsize = FALSE)

scMetabolism分析细胞代谢
https://b.limour.top/1539.html
Author
Limour
Posted on
February 9, 2022
Licensed under