1. rm(list = ls()) 
  2. load(file = "step1output.Rdata")
  3. load(file = "step4output.Rdata")

1.火山图

library(dplyr)
library(ggplot2)
dat  = deg[!duplicated(deg$symbol),]

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

for_label <- dat%>% 
  filter(symbol %in% c("HADHA","LRRFIP1"))

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_number,"_volcano.png"))

2.差异基因热图

load(file = 'step2output.Rdata')
# 表达矩阵行名替换
exp = exp[dat$probe_id,]
rownames(exp) = dat$symbol
if(F){
  #全部差异基因
  cg = dat$symbol[dat$change !="stable"]
  length(cg)
}else{
  #取前10上调和前10下调
  x=dat$logFC[dat$change !="stable"] 
  names(x)=dat$symbol[dat$change !="stable"] 
  cg=names(c(head(sort(x),10),tail(sort(x),10)))
  length(cg)
}
n=exp[cg,]
dim(n)

#差异基因热图
library(pheatmap)
annotation_col=data.frame(group=Group)
rownames(annotation_col)=colnames(n) 
heatmap_plot <- pheatmap(n,show_colnames =F,
                         scale = "row",
                         #cluster_cols = F, 
                         annotation_col=annotation_col,
                         breaks = seq(-3,3,length.out = 100)
) 
heatmap_plot
ggsave(heatmap_plot,filename = paste0(gse_number,"_heatmap.png"))

3.感兴趣基因的箱线图

g = c(head(cg,3),tail(cg,3))
library(tidyr)
library(tibble)
library(dplyr)
dat = t(exp[g,]) %>% 
  as.data.frame() %>% 
  rownames_to_column() %>% 
  mutate(group = Group)

pdat = dat%>% 
  pivot_longer(cols = 2:(ncol(dat)-1),
               names_to = "gene",
               values_to = "count")

pdat$gene = factor(pdat$gene,levels = cg,ordered = T)
pdat$change = ifelse(pdat$gene %in% head(cg,10),"down","up")
library(ggplot2)
library(paletteer)
box_plot = ggplot(pdat,aes(gene,count))+
  geom_boxplot(aes(fill = group))+
  #scale_fill_manual(values = c("blue","red"))+
  scale_fill_paletteer_d("basetheme::minimal")+
  geom_jitter()+
  theme_bw()+
  facet_wrap(~change,scales = "free")
box_plot
ggsave(box_plot,filename = paste0(gse_number,"_boxplot.png"))

4.感兴趣基因的相关性

library(corrplot)
M = cor(t(exp[g,]))
pheatmap(M)

my_color = rev(paletteer_d("RColorBrewer::RdYlBu"))
my_color = colorRampPalette(my_color)(10)
corrplot(M, type="upper",
         method="pie",
         order="hclust", 
         col=my_color,
         tl.col="black", 
         tl.srt=45)
library(cowplot)
cor_plot <- recordPlot()

5.拼图

load("pca_plot.Rdata")
library(patchwork)
library(ggplotify)
(pca_plot + volcano_plot +as.ggplot(heatmap_plot))/box_plot

plot_grid(cor_plot,heatmap_plot$gtable)