美文网首页基因组数据绘图
复现Immunity文章图表:分组雷达图展示富集结果

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

作者: KS科研分享与服务 | 来源:发表于2023-06-24 16:26 被阅读0次

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

    image.png

    (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)
    
    
    image.png

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

    相关文章

      网友评论

        本文标题:复现Immunity文章图表:分组雷达图展示富集结果

        本文链接:https://www.haomeiwen.com/subject/viybydtx.html