####gsva--分樣本進行--可
library(Seurat)
library(ggplot2)
library(dplyr)
library(pheatmap)
library(patchwork)
library(msigdbr)
library(GSVA)
setwd("/data/wanglei_lab/zhangyu/2.zy_fdy/2.TBILJN_vs_MC/tbi-20221017/4.recluster/1.Mac/Mac_GSVA/Mac_gsva_group")
Mac <- readRDS("/data/wanglei_lab/zhangyu/2.zy_fdy/2.TBILJN_vs_MC/tbi-20221017/4.recluster/1.Mac/newid_Mac/newid_Mac.rds")
head(Mac)
## 表達(dá)矩陣
expr=as.matrix(Mac@assays$RNA@data)
###預(yù)定基因集,后續(xù)gsva要求是list，所以將他轉(zhuǎn)化為list
mdb_c2 <- msigdbr(species = "Mus musculus", category = "C2")  #選取物種小鼠.GO在C5馏臭，KEGG在C2
mdb_c2 = mdb_c2 [grep("^KEGG",mdb_c2 $gs_name),] ##R語言中的grep函數(shù)可以在給定的字符串向量中搜索某個子字符串
mdb_c2$gs_name <- gsub('KEGG_','',mdb_c2$gs_name)   #去除前綴KEGG_
mdb_c2$gs_name <- tolower(mdb_c2$gs_name)  #將大寫換為小寫
mdb_c2$gs_name <- gsub('_',' ',mdb_c2$gs_name)  #將_轉(zhuǎn)化為空格
msigdbr_list = split(x = mdb_c2$gene_symbol, f = mdb_c2$gs_name) ##后續(xù)gsva要求是list蹦狂，所以將他轉(zhuǎn)化為list
head(msigdbr_list)
kegg <- gsva(expr, gset.idx.list = msigdbr_list, kcdf="Gaussian",method = "zscore",parallel.sz=1)
write.table(kegg, 'Mac_gsva_kegg.xls', row.names=T, col.names=NA, sep="\t")  ##通路和細(xì)胞的表達(dá)矩陣

library(limma)   
## limma gsva通路活性評估
de_gsva <- function(exprSet,meta,compare = NULL){
  allDiff = list()
  design <- model.matrix(~0+factor(meta))
  colnames(design)=levels(factor(meta))
  rownames(design)=colnames(exprSet)  
  fit <- lmFit(exprSet,design)
  if(length(unique(meta))==2){
    if(is.null(compare)){
      stop("there are 2 Groups,Please set  compare value")
    }
    contrast.matrix<-makeContrasts(contrasts = compare,levels = design)
    fit2 <- contrasts.fit(fit, contrast.matrix) 
    fit2 <- eBayes(fit2)
    tempOutput = topTable(fit2,adjust='fdr', coef=1, number=Inf)
    allDiff[[compare]] = na.omit(tempOutput)
    
  }else if(length(unique(meta))>2){
    for(g in colnames(design)){
      fm = ""
      for(gother in colnames(design)[which(!colnames(design) %in% g)]){
        fm = paste0(fm,"+",gother)
      } 
      fm = paste0(g,"VsOthers = ",g,"-(",substring(fm,2),")/",ncol(design)-1)
      contrast.matrix <- makeContrasts(contrasts = fm,levels=design)
      fit2 <- contrasts.fit(fit, contrast.matrix) 
      fit2 <- eBayes(fit2)
      allDiff[[g]]=topTable(fit2,adjust='fdr',coef=1,number=Inf)
    }
  }else{
    stop("error only have one group")
  }
  return(allDiff)
}
meta <- Mac@meta.data[,c("group")]
Diff =de_gsva(exprSet = kegg,meta = meta,compare = "TBI-MC")  ##"exp/ctrl"
idiff <-Diff[["TBI-MC"]]
df <- data.frame(ID = rownames(idiff), score = idiff$t )
Padj_threshold=0.01
df$label =sapply(1:nrow(idiff),function(x){
  if(idiff[x,"logFC"]>0 & idiff[x,"adj.P.Val"]<Padj_threshold){return("up")}
  else if(idiff[x,"logFC"]<0 & idiff[x,"adj.P.Val"]<Padj_threshold){return("down")}
  else{return("noSig")}
})

# 按照score排序
df$hjust = ifelse(df$score>0,1,0)
df$nudge_y = ifelse(df$score>0,-0.1,0.1)
sortdf <- df[order(df$score),]
sortdf$ID <- factor(sortdf$ID, levels = sortdf$ID)
limt = max(abs(df$score))
write.csv(df,"gsva_df.csv",row.names = F)
write.csv(sortdf,"gsva_score.csv",row.names = F)  ##給老師發(fā)這個
###篩選之后再讀入
df<-read.csv('gsva_df1.csv',header = T)
sortdf<-read.csv('gsva_score1.csv',header = T)
##再排序一下作圖才能有順序
df$hjust = ifelse(df$score>0,1,0)
df$nudge_y = ifelse(df$score>0,-0.1,0.1)
sortdf <- df[order(df$score),]
sortdf$ID <- factor(sortdf$ID, levels = sortdf$ID)
pdf("down_up_GSVA.pdf",width = 30, height = 30)
ggplot(sortdf, aes(ID, score,fill=label)) + 
  geom_bar(stat = 'identity',alpha = 0.7) + 
  scale_fill_manual(breaks=c("down","noSig","up"),
                    values = c("#008020","grey","#08519C"))+
  geom_text(data = df, aes(label = ID, y = nudge_y),
            nudge_x =0,nudge_y =0,hjust =df$hjust,
            size = 10)+
  labs(x = (""),
       y="t value of GSVA score\n")+
  scale_y_continuous(limits=c(-limt,limt))+
  coord_flip() + 
  theme_bw() + #去除背景色
  theme(panel.grid =element_blank())+
  theme(panel.border = element_rect(size = 0.6)
        #panel.border = element_blank()
  )+
  theme(plot.title = element_text(hjust = 0.5,size = 18),
        axis.text.y = element_blank(),
        axis.title = element_text(hjust = 0.5,size = 18),
        axis.line = element_blank(),
        axis.ticks.y = element_blank(),
        legend.position = limt
  )
dev.off()