复现Immunity文章图表：分组雷达图展示富集结果

前不久我们写过一期雷达图的做法（雷达图展示GSEA分析结果），用的fmsb包。后来有小伙伴反应说有比这个更好的雷达图，是一篇Immunity上的文章，也是利用雷达图展示通路富集的结果，其实也是一样，只不过是分组了，我们还是用fmsb包进行复现。

（reference：Cross-tissue single-cell landscape of human monocytes and macrophages in health and disease）这里我们也复现做一个分组的雷达图。首先我们分析了一组数据的差异基因，利用上下调差异基因做了GO富集，展示相关的结果。差异分析和GO富集都是很简单的东西了，就不再赘述了。

#首先我们做一下一组数据的差异分析，分为上下条两组基因进行富集分析
setwd('D:/KS项目/公众号文章/雷达图/分组雷达图')
library(tidyverse)
library(DESeq2)
library(combinat)
library(clusterProfiler)
library(ggplot2)
Two_group <- read.csv("Two_group.csv", header = T,row.names = 1)
meta1 <- data.frame(Cancer=c("Cancer1" ,"Cancer2" ,"Cancer3"),
                    Health=c("Health1", "Health2", "Health3"))

DEGs <- DEseq_multi_group(exprSet = Two_group,meta = meta1,repNum1 = 3,
                          repNum2 = 3,test = 'Cancer', control = 'Health')
head(DEGs)#运行成功
DEGs <- na.omit(DEGs)


#我们这里做一下GO分析
df_sig <- WT_DEGs[[1]][which(WT_DEGs[[1]]$pvalue<=0.05 & abs(WT_DEGs[[1]]$log2FoldChange)>=1),]
df_sig$gene <- rownames(df_sig)
df_sig$group <- ''#新加一列
df_sig$group <- ifelse(df_sig$log2FoldChange >0,'up','down')#上下调标记

#富集
group <- data.frame(gene=df_sig$gene,group=df_sig$group)
Gene_ID <- bitr(df_sig$gene, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Mm.eg.db")

#构建文件并分析
data  <- merge(Gene_ID,group,by.x='SYMBOL',by.y='gene')
data_GO <- compareCluster(ENTREZID~group, data=data, fun="enrichGO", 
                          OrgDb="org.Mm.eg.db",ont = "BP",pAdjustMethod = "BH",
                          pvalueCutoff = 0.05,qvalueCutoff = 0.1)

data_GO_sim <- clusterProfiler::simplify(data_GO,cutoff=0.7, by="p.adjust", select_fun=min)


GO_result <- data_GO_sim@compareClusterResult
write.csv(GO_result, file = 'GO_result.csv')

我们整理一下，挑选需要的通路进行展示。复现基本上90%是完成了，但是有一点没有，就是点用渐变来实现，不知道这个包有没有办法。

df <- read.csv('GO_result.csv', header = T,row.names = 1)
df1 <- df[, c("group","Description","p.adjust")]

library(tidyr)
df1<-spread(df1, Description, p.adjust)

df1[is.na(df1)] <-  1
rownames(df1) <- df1[,1]
df1 <- df1[,-1]
df1 <- -log(df1)

df1 <- t(df1)
df1 <- as.data.frame(df1)
df1 <- arrange(df1,df1$down,df1$up)

df1 <- t(df1)
df1 <- as.data.frame(df1)

my.data <- matrix( c(rep(max(df1),ncol(df1)), 
                     rep(0, ncol(df1)), 
                     rep(-log(0.05), ncol(df1))), nrow = 3, ncol = ncol(df1), byrow=TRUE)


colnames(my.data) <- colnames(df1)
rownames(my.data) <- c("max", "min", "p")

my.data <- rbind(my.data, df1)
my.data <- my.data[c(1,2,4,5,3),]


library(fmsb)
radarchart(my.data, 
           pty = c(16,16,32),
           axistype = 1,
           pcol = c("#64299C", "#0439FD","black"), 
           pfcol = c(scales::alpha(c("#64299C", "#0439FD","black"), c(0.5, 0.5, 0.5))),
           plwd = c(3,3,3),
           plty = 1,
           cglcol = "grey60", 
           cglty = 1, 
           cglwd = 1,
           axislabcol = "grey60",
           vlcex = 0.8, 
           vlabels = colnames(colnames(my.data)),
           caxislabels = c(0, 10, 20, 30, 40),
           calcex=0.8
           )


legend(x = "bottomright", legend = c("down","up"), horiz = F,
  bty = "n", pch = 15 , col = c("#64299C", "#0439FD"),
  text.col = "black", cex = 1, pt.cex = 1.5)


legend(x = "center", legend = c("p<0.05"), horiz = TRUE,
       text.col = "white", cex = 1,bg=NULL,box.lty=0)

当然了也有其他的包做雷达图，例如ggradar或者ggplot2也能实现，感兴趣的可以自行探索，我们就不再深究了，这个展示应用还是可以的。觉得分享有用的点个赞再走呗！