1. rm(list = ls()) 
  2. load(file = "step1output.Rdata")
  3. load(file = "step4output.Rdata")
  4. #1.火山图----
  5. library(dplyr)
  6. library(ggplot2)
  7. dat  = deg[!duplicated(deg$symbol),]
  9. p <- ggplot(data = dat, 
  10.             aes(x = logFC, 
  11.                 y = -log10(P.Value))) +
  12.   geom_point(alpha=0.4, size=3.5, 
  13.              aes(color=change)) +
  14.   ylab("-log10(Pvalue)")+
  15.   scale_color_manual(values=c("blue", "grey","red"))+
  16.   geom_vline(xintercept=c(-logFC_t,logFC_t),  ##竖线
  17.              lty=4,col="black",lwd=0.8) +
  18.   geom_hline(yintercept = -log10(P.Value_t),  ##横线
  19.              lty=4,col="black",lwd=0.8) +
  20.   theme_bw()
  21. p
  23. for_label <- dat%>% 
  24.   filter(symbol %in% c("HADHA","LRRFIP1"))
  26. volcano_plot <- p +
  27.   geom_point(size = 3, shape = 1, data = for_label) +
  28.   ggrepel::geom_label_repel(          ##加label的图层，加标签的图层
  29.     aes(label = symbol),
  30.     data = for_label,
  31.     color="black"
  32.   )
  33. volcano_plot
  34. ggsave(plot = volcano_plot,filename = paste0(gse_number,"_volcano.png"))

2.差异基因热图——

  1. load(file = 'step2output.Rdata')

行名替换

  1. exp = exp[dat$probe_id,] 
  2. rownames(exp) = dat$symbol
  3. if(F){
  4.   #全部差异基因
  5.   cg = dat$symbol[dat$change !="stable"]
  6.   length(cg)
  7. }else{
  8.   #取前10上调和前10下调
  9.   x=dat$logFC[dat$change !="stable"] 
  10.   names(x)=dat$symbol[dat$change !="stable"] 
  11.   cg=names(c(head(sort(x),10),tail(sort(x),10)))
  12.   length(cg)
  13. }
  14. n=exp[cg,]
  15. dim(n)

差异基因热图

  1. library(pheatmap)
  2. annotation_col=data.frame(group=Group)
  3. rownames(annotation_col)=colnames(n) 
  4. heatmap_plot <- pheatmap(n,show_colnames =F,
  5.                          scale = "row",
  6.                          #cluster_cols = F, 
  7.                          annotation_col=annotation_col,
  8.                          breaks = seq(-3,3,length.out = 100)
  9. ) 
  10. heatmap_plot
  11. ggsave(heatmap_plot,filename = paste0(gse_number,"_heatmap.png"))

感兴趣基因的箱线图

  1. g = c(head(cg,3),tail(cg,3))
  2. library(tidyr)
  3. library(tibble)
  4. library(dplyr)
  5. dat = t(exp[g,]) %>% 
  6.   as.data.frame() %>% 
  7.   rownames_to_column() %>% 
  8.   mutate(group = Group)
  10. pdat = dat%>% 
  11.   pivot_longer(cols = 2:(ncol(dat)-1),
  12.                names_to = "gene",
  13.                values_to = "count")
  15. pdat$gene = factor(pdat$gene,levels = cg,ordered = T)
  16. pdat$change = ifelse(pdat$gene %in% head(cg,10),"down","up")
  17. library(ggplot2)
  18. library(paletteer)
  19. box_plot = ggplot(pdat,aes(gene,count))+
  20.   geom_boxplot(aes(fill = group))+
  21.   #scale_fill_manual(values = c("blue","red"))+
  22.   scale_fill_paletteer_d("basetheme::minimal")+
  23.   geom_jitter()+
  24.   theme_bw()+
  25.   facet_wrap(~change,scales = "free")
  26. box_plot
  27. ggsave(box_plot,filename = paste0(gse_number,"_boxplot.png"))
  29. # 4.感兴趣基因的相关性----
  30. library(corrplot)
  31. M = cor(t(exp[g,]))   ##默认求列的相关性
  32. pheatmap(M)
  34. my_color = rev(paletteer_d("RColorBrewer::RdYlBu"))
  35. my_color = colorRampPalette(my_color)(10)
  36. corrplot(M, type="upper",
  37.          method="pie",
  38.          order="hclust", 
  39.          col=my_color,
  40.          tl.col="black", 
  41.          tl.srt=45)
  42. library(cowplot)
  43. cor_plot <- recordPlot()

感兴趣基因的相关性

  1. # 4.感兴趣基因的相关性----
  2. library(corrplot)
  3. M = cor(t(exp[g,]))   ##默认求列的相关性
  4. pheatmap(M)
  6. my_color = rev(paletteer_d("RColorBrewer::RdYlBu"))
  7. my_color = colorRampPalette(my_color)(10)
  8. corrplot(M, type="upper",
  9.          method="pie",
  10.          order="hclust", 
  11.          col=my_color,
  12.          tl.col="black", 
  13.          tl.srt=45)
  14. library(cowplot)
  15. cor_plot <- recordPlot()

拼图

  1. load("pca_plot.Rdata")
  2. library(patchwork)
  3. library(ggplotify)
  4. (pca_plot + volcano_plot +as.ggplot(heatmap_plot))/box_plot
  6. plot_grid(cor_plot,heatmap_plot$gtable)