1. 读入数据

差异表达基因来自limma分析结果。

# read the file
data <- read.csv("diff_expr_genes.csv",row.names=1)
# sorting
data = diff_genes[order(diff_genes$adj.P.Val),]
head(data)
class(data)

colnames(data)
dim(data)

# 如果之前保存的是R对象
save(nrDEG_limma_voom, file = 'diff_expr_genes')
load('diff_expr_genes')
data <- nrDEG_limma_voom

2. ggplot作火山图

# 颜色设定
data$color <- ifelse(data$adj.P.Val<0.05 & abs(data$logFC)>= 1,
                     ifelse(data$logFC > 1,'red','blue'),'gray')

#install.packages("ggplot2")
#install.packages("ggrepel")

# 导入包
library(ggplot2)
library(ggrepel) # labels 不重叠

#tiff(filename = "volcano.tif",width = 720, height = 720)
rownames(data) <- data$symbol

# 设定要标出的基因名
data$sign <- ifelse(data$adj.P.Val < 1e-10 & abs(data$logFC) >4,
                    rownames(data),NA)

p <- ggplot(data, aes(logFC, -log10(adj.P.Val), col = color)) +
  geom_point() +
  theme_bw() +
  scale_color_manual(values = color) +
  labs(x="log2 (fold change)",y="-log10 (q-value)") +
  geom_hline(yintercept = -log10(0.05), lty=4,col="grey",lwd=0.6) +
  geom_vline(xintercept = c(-1, 1), lty=4,col="grey",lwd=0.6) +
  theme(legend.position = "none",
        panel.grid=element_blank(),
        axis.title = element_text(size = 18),
        axis.text = element_text(size = 14)) +
  geom_text_repel(aes(label = sign), size = 3)
print(p)

#dev.off()
# ggsave保存
ggsave("ggsave_volcano.tif",dpi=100,plot=p,device = 'tiff')

#device:    "eps", "ps", "tex" (pictex), "pdf", "jpeg", 
#"tiff", "png", "bmp","svg" or "wmf" (windows only)