setwd("D:/geo/GSE70689/")

DEG = read.csv("CBD_LPS-LPS.csv", header = T, row.names = 1)

#全部的数据做火山图
vol = data.frame("log2FC"=as.numeric(DEG$logFC),"p"=as.numeric(DEG$adj.P.Val))
rownames(vol) = rownames(DEG)
dim(vol)
#[1] 19709    2
vol = na.omit(vol)
dim(vol)
#[1] 18709    2

#标记分组，先将所有的标记为没有差异的，然后通过循环将其标记为上调UP，和下调Down
vol$SI = "NotDEG"
for(i in 1:length(vol[,1])){
  if(vol[i,1]> 1 & vol[i,2]< 0.05){
    vol[i,3] = "UP"
  }else if(vol[i,1]< -1 & vol[i,2]< 0.05){
    vol[i,3] = "Down"
  }
}
vol$SI=as.factor(vol$SI)

print(length(rownames(vol[vol$SI=="UP",])))
# 111
print(length(rownames(vol[vol$SI=="Down",])))
# 28

library(ggplot2)
pdf("Vol.pdf", width=8, height=8)
layer1 <- ggplot(
  vol,
  aes(x=log2FC,y=-1*log10(p))
  )+
  geom_point(aes(color=SI),alpha=0.4)+ 
  xlim(-5,5) + ylim(0,8)+ #设置x和y轴的范围
  scale_color_manual(values =c("green","blue", "red"))+ 
  labs(title="mRNA",x="log2(fold change)", y="-log10(p value)")+
  geom_hline(yintercept=1.30103,linetype=3)+
  geom_vline(xintercept=c(-1,1),linetype=3)
  #theme(legend.position="none",title=element_text(size=16,color="black",hjust=0.2,lineheight=0.2),axis.title.x=element_text(size=14,face="bold",color="black",hjust=0.5))
layer1
dev.off()

#换一种风格
pdf("Vol2.pdf", width=8, height=8)
p<-ggplot( vol, 
           aes(x = log2FC, 
               y = -log10(p), 
               colour=SI)) +
  geom_point(alpha=0.4, size=3.5) +
  scale_color_manual(values=c("#546de5", "#d2dae2","#ff4757"))+
  geom_vline(xintercept=c(-1,1),lty=4,col="black",lwd=0.8) +
  geom_hline(yintercept = 1.30103,lty=4,col="black",lwd=0.8) +
  labs(x="log2(fold change)",
       y="-log10 (p-value)")+
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5), 
        legend.position="right", 
        legend.title = element_blank()
  )
p
dev.off()


#在火山图中标记按P值排序,取前10的基因
library("ggrepel")

mRNA =DEG[order(DEG$adj.P.Val),]
mRNA = mRNA[1:10,]

TopTenGenes = row.names(mRNA)
TopTenGenes
#[1] "Mt2"           "Vegfa"         "Mmp23"         "Aqp9"         
#[5] "2210008I11Rik" "2810428C21Rik" "Glrp1"         "Slc6a9"       
#[9] "Mtbp"          "Trib3" 
# 也可以将TopTenGenes改为所需要展示的基因的名称
# 如 TopTemGenes = c("noc4l","p53")


vol$label = ""
for(i in 1:nrow(vol)){
  if(row.names(vol)[i] %in% TopTenGenes){
    vol[i,'label']=row.names(vol)[i]
  }
}
pdf("Vol2_label.pdf", width=8, height=8)
p<-ggplot( vol, 
           aes(x = log2FC, 
               y = -log10(p), 
               color=SI)) +
  geom_point(alpha=0.4, size=3.5) +
  scale_color_manual(values=c("#546de5", "#d2dae2","#ff4757"))+
  geom_vline(xintercept=c(-1,1),lty=4,col="black",lwd=0.8) +
  geom_hline(yintercept = 1.30103,lty=4,col="black",lwd=0.8) +
  labs(x="log2(fold change)",
       y="-log10 (p-value)")+
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5), 
        legend.position="right", 
        legend.title = element_blank()
  )
p+geom_text_repel(data = vol, aes(x = vol$log2FC, 
                                      y = -log10(vol$p), 
                                      label = label),
                  size = 3,box.padding = unit(0.5, "lines"),
                  point.padding = unit(0.8, "lines"), 
                  segment.color = "#000000", 
                  show.legend = FALSE)
dev.off()