#paper# METAFlux: Characterizing metabolism from bulk and single-cell RNA-seq data# https://www./articles/s41467-023-40457-w
#github:https://github.com/KChen-lab/METAFlux?tab=readme-ov-file
devtools::install_github('KChen-lab/METAFlux')library(METAFlux)library(Seurat)

load("C:/Users/tq199/Downloads/adj_scRNA.RData")

#我們單獨(dú)提取一個(gè)組的數(shù)據(jù)進(jìn)行演示SD_scRNA <- subset(adj_scRNA, orig.ident=='SD_scRNA')mean_exp=calculate_avg_exp(myseurat = SD_scRNA, #需要分析的單細(xì)胞seurat obj,單細(xì)胞obj是normalization 的                           myident = 'celltype',#celltype分組                           n_bootstrap=3, #抽樣次數(shù)                           seed=1)#隨機(jī)數(shù)
#計(jì)算代謝反應(yīng)分?jǐn)?shù)scores<-calculate_reaction_score(data=mean_exp)#計(jì)算這個(gè)樣本中每種celltype占比，每種celltype的數(shù)量除以總的細(xì)胞數(shù)#需要注意，這里celltype fraction順序要與score里面celltype的順序一致fr =  round(table(SD_scRNA$celltype)/nrow(SD_scRNA@meta.data),2)fr = as.vector(fr)
#計(jì)算代謝通量flux=compute_sc_flux(num_cell = 10, #celltype數(shù)量,比如這里我的數(shù)據(jù)中有10種celltype                     fraction = fr, #celltype比例，和必須是1                     fluxscore=scores,                     medium = human_blood)#

#代謝model，可以查看代謝反應(yīng)ID，代謝公式，代謝通路！data("human_gem")human_gem <- human_gem
#選擇感興趣的代謝反應(yīng)，提取數(shù)據(jù)進(jìn)行可視化#比如這里看glucose的代謝反應(yīng)，我們可以通過nutrient_lookup_files文件查找glucose對應(yīng)的IDdata("nutrient_lookup_files")glucose_ID <- nutrient_lookup_files$EQUATION[grepl("glucose", nutrient_lookup_files$EQUATION)]select_ID <- subset(nutrient_lookup_files, nutrient_lookup_files$EQUATION %in% glucose_ID)#這樣我們就可以得知glucose IDglucose<-data.frame(glucose=flux[grep("HMR_9034",rownames(flux)),])
#if needed we can apply cubic root normalization to normalize the scorescbrt <- function(x) {  sign(x) * abs(x)^(1/3)}
glucose=cbrt(glucose)
#ggplot做一個(gè)箱線圖展示library(ggplot2)glucose$celltype=rownames(glucose)long_glucose=reshape2::melt(glucose,id.vars='celltype')ggplot(long_glucose,aes(y=-value,x=celltype))+  geom_boxplot()+  ggtitle("Glucose uptake level for different cell types")+  xlab("")+ylab("Glucose uptake scores")+   theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

flux1 <- rowMeans(flux) %>% as.data.frame()colnames(flux1) <- "flux"
celltype <- c("Endo","Fib","Musc","Ker","APCs","Other","Tcell","Mast","LY","Mela")#順序和mean-exp順序一樣celltype <- paste0(celltype," ")flux_celltype <- paste("celltype", seq(length(celltype)),"")for(i in 1:length(celltype)){
  rownames(flux1) <-  gsub(flux_celltype[i], celltype[i], rownames(flux1))
}
flux1$celltype <- sub(" .*$", "", rownames(flux1))flux1 <- flux1[!grepl("external_medium", flux1$celltype), ]flux1$pathway <- rep(human_gem$SUBSYSTEM, length(celltype))
library(tidyr)library(dplyr)library(pheatmap)
df_summary <- flux1 %>%  group_by(celltype, pathway) %>%  summarize(mean_value = mean(abs(flux)))
df_summary_w <- df_summary  %>%  pivot_wider(names_from = pathway, values_from = mean_value) %>%  as.data.frame()

rownames(df_summary_w) <- df_summary_w[,1]df_summary_w <- df_summary_w[,-1]

#heatmap mapal <- colorRampPalette(RColorBrewer::brewer.pal(11,"RdBu"))(256)pheatmap::pheatmap(t(df_summary_w),cluster_cols = F,color = rev(mapal),scale = "row")

#load your seurat object,多樣本數(shù)據(jù)循環(huán)分析# obj.list <- SplitObject(seurat, split.by = "orig.ident")# for (i in c(1:length(obj.list))){#   sc<-obj.list[[i]]#   mean_exp=calculate_avg_exp(myseurat = sc,myident = 'celltype',n_bootstrap=50,seed=1)#   scores<-calculate_reaction_score(data=mean_exp)#   fr =  round(table(sc$celltype)/nrow(sc@meta.data),2)#   fr = as.vector(fr)#   flux=compute_sc_flux(num_cell = length(unique(sc$celltype)),#                        fraction = fr,#                        fluxscore=scores,medium = human_blood)#   saveRDS(flux,paste0(unique(sc$orig.ident),"_flux",".rds"))# }