p_load(METAFlux, tidyverse)
### bulk Rnaseq
dir.create('./metaflux')
tpm <- read.csv('tpm.csv', row.names = 1)
data('human_gem', 'cell_medium', 'human_blood')
scores <- calculate_reaction_score(tpm)
flux <- compute_flux(mras=scores, medium = cell_medium)

### cubic root normalization to flux scores, optional
cbrt <- function(x) {sign(x) * abs(x)^(1/3)}
flux <- cbrt(flux)

### pathway level activity
pathway <- unique(unlist(human_gem$SUBSYSTEM))
pathway_score <- list()
for (i in pathway){
  path <- i
  activity_score <- c()
  for (d in 1:ncol(flux)){
    activity_score[d] <- mean(abs(flux[which(unlist(human_gem$SUBSYSTEM)==i),d]))
  } 
  pathway_score[[i]] <- activity_score
}
all_pathway_score <- as.data.frame(do.call(rbind,pathway_score))
colnames(all_pathway_score) <- colnames(tpm)
pdf('./metaflux/pathway_flux.pdf', height = 18, width = 10)
mapal <- colorRampPalette(RColorBrewer::brewer.pal(11,"RdBu"))(256)
pheatmap::pheatmap(all_pathway_score, cluster_cols = F, color = rev(mapal), scale = "row", 
                   cutree_rows = 2, gaps_col = 3, angle_col = 45)
dev.off()
write.csv(all_pathway_score, './metaflux/all_pathway_score.csv')

### differential analysis via limma
p_load(limma)
list <- gsub('_[0-9]+$','',colnames(all_pathway_score)) %>%
  factor(., levels = unique(.), ordered = F)
list <- model.matrix(~factor(list) + 0)
colnames(list) <- c('ctrl', 'treat')
df.fit <- lmFit(all_pathway_score, list)
df.matrix <- makeContrasts(treat - ctrl, levels = list)
fit <- contrasts.fit(df.fit, df.matrix)
fit <- eBayes(fit)
res <- topTable(fit, n = Inf, adjust = 'fdr')
write.csv(res, './metaflux/pathway_score_diff.csv')