GEO代码分析流程 - 5. 差异分析 - 作图 - 火山图、热图
作者:互联网
5. 差异分析 - 作图 - 火山图、热图
rm(list = ls()) load(file = "step1output.Rdata") load(file = "step4output.Rdata") #1.火山图---- library(dplyr) library(ggplot2) dat = deg p <- ggplot(data = dat, aes(x = logFC, y = -log10(P.Value))) + geom_point(alpha=0.4, size=3.5, aes(color=change)) + ylab("-log10(Pvalue)")+ scale_color_manual(values=c("blue", "grey","red"))+ geom_vline(xintercept=c(-logFC_t,logFC_t),lty=4,col="black",lwd=0.8) + geom_hline(yintercept = -log10(P.Value_t),lty=4,col="black",lwd=0.8) + theme_bw() p if(T){ #自选基因 for_label <- dat%>% filter(symbol %in% c("TRPM3","SFRP1")) } if(F){ #p值最小的10个 for_label <- dat %>% head(10) } if(F) { #p值最小的前3下调和前3上调 x1 = dat %>% filter(change == "up") %>% head(3) x2 = dat %>% filter(change == "down") %>% head(3) for_label = rbind(x1,x2) } volcano_plot <- p + geom_point(size = 3, shape = 1, data = for_label) + ggrepel::geom_label_repel( aes(label = symbol), data = for_label, color="black" ) volcano_plot ggsave(plot = volcano_plot,filename = paste0(gse,"volcano.png")) #2.差异基因热图---- load(file = 'step2output.Rdata') if(F){ #全部差异基因 cg = deg$probe_id[deg$change !="stable"] length(cg) }else{ #取前30上调和前30下调 x=deg$logFC[deg$change !="stable"] names(x)=deg$probe_id[deg$change !="stable"] cg=c(names(head(sort(x),30)),names(tail(sort(x),30))) length(cg) } n=exp[cg,] dim(n) #作热图 library(pheatmap) annotation_col=data.frame(group=group_list) rownames(annotation_col)=colnames(n) library(ggplotify) heatmap_plot <- as.ggplot(pheatmap(n,show_colnames =F, show_rownames = F, scale = "row", #cluster_cols = F, annotation_col=annotation_col)) heatmap_plot ggsave(heatmap_plot,filename = paste0(gse,"heatmap.png")) load("pca_plot.Rdata") library(patchwork) (pca_plot + volcano_plot +heatmap_plot)+ plot_annotation(tag_levels = "A")
标签:load,filter,head,作图,library,dat,file,热图,GEO 来源: https://www.cnblogs.com/xiaogaobugao/p/16675834.html