rm(list = ls())
###############################input data 
dir_path <- "C:\\Users\\liyix\\OneDrive\\Desktop\\data\\"
dir_path_name <- list.files(pattern = ".*csv",dir_path,full.names = T, recursive = F)
#dir_path_name
volcanic <- read.csv(grep("data.csv",dir_path_name,value = T),header = T,stringsAsFactors = F)
#################################plot
library(ggplot2)
#summary(volcanic)
#head(volcanic)
#tail(volcanic)
#dim(volcanic) #[1] 20449     3
volcanic$X <- NULL
#View(volcanic)
volcanic_keep <- volcanic[,c("logFC", "ttest_pvalue")]
#dim(volcanic_keep) #[1] 20449     2
##Highlight genes that have an absolute fold changeand a p-value 
df <- volcanic_keep
colnames(df) <- c("logFC", "P.Value")
#dim(df)  #[1] 20449     2
###define Green----downregulation 
df.G <- subset(df, logFC < 0 & P.Value < 0.01) #define Green
df.G <- cbind(df.G, rep("Downregulation", nrow(df.G)))
colnames(df.G)[3] <- "Color"
#head(df.G)
#df.G$logFC <- -log2(abs(df.G$logFC))
#View(df.G)
###define Black
df.B <- subset(df, P.Value >= 0.01) #define Black
df.B <- cbind(df.B, rep("Not change", nrow(df.B)))
colnames(df.B)[3] <- "Color"
#less than zero
#head(df.B)
df.B_down <- df.B[which(df.B$logFC < 0),]
#head(df.B_down) 
#df.B_down$logFC <- -log2(abs(df.B_down$logFC))
#View(df.B_down)
#more than zero
df.B_up <- df.B[which(df.B$logFC > 0),]
#head(df.B_up) 
#df.B_up$logFC <- log2(df.B_up$logFC)
#View(df.B_up)
##define Red
df.R <- subset(df, logFC > 0 & P.Value < 0.01) #define Red
df.R <- cbind(df.R, rep("Upregulation", nrow(df.R)))
colnames(df.R)[3] <- "Color"
#head(df.R)
#dim(df.R) #[1] 29  3
#View(df.R)
#df.R$logFC <- log2(abs(df.R$logFC))
##rbind  all 
df.t <- rbind(df.G, df.B_down, df.B_up, df.R)
#View(df.t)
#str(df.t)
df.t$Color <- as.factor(df.t$Color)
dim(df.t); head(df.t);table(df.t$Color)
## [1] 20449     3
##            logFC     P.Value          Color
## 545   -0.3910112 0.009925008 Downregulation
## 1649  -0.1237625 0.006906282 Downregulation
## 7263  -0.1512422 0.006468623 Downregulation
## 10572 -0.1964464 0.002922108 Downregulation
## 11909 -0.1972898 0.006069207 Downregulation
## 12154 -0.3572661 0.005005542 Downregulation
## 
## Downregulation     Not change   Upregulation 
##             11          20409             29
##Construct the plot object
ggplot(data = df.t, aes(x = logFC, y = -log10(P.Value), color= Color )) +
  geom_point(alpha = 0.5, size = 1.75) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"),
        legend.key = element_rect(colour = "transparent", fill = "white"),
        legend.title=element_blank(),
        legend.text=element_text(size=15),
        legend.position="top",
        axis.text.x = element_text(colour="black",size=15,face="plain"),
        axis.text.y = element_text(colour="black",size=15,face="plain"),
        axis.title.x = element_text(colour="black",size=15,face="plain"),
        axis.title.y = element_text(colour="black",size=15,face="plain"),
        panel.border = element_rect(colour = "black", fill=NA, size=0.5)) +
  scale_color_manual(values = c("green", "gray", "red")) +
  labs(x=expression("logFC"), y=expression( -log[10](p.value))) +
  scale_x_continuous(breaks= seq(from = -1 , to = 1 ,by = 0.2)) +
  scale_y_continuous(breaks= seq(from = 0 , to = 4 ,by = 1)) +
geom_hline(yintercept= -log10(0.05), linetype="dashed", 
           color = "blue", size= 0.5) +
  geom_vline(xintercept = -0.5, linetype="dashed", 
             color = "black", size= 0.5) +
  geom_vline(xintercept = 0.5, linetype="dashed", 
             color = "black", size= 0.5)

###############
################output
ggsave(paste0(Sys.Date(),"-volcano_plot.tif"), plot = last_plot(), device = "tiff", path = dir_path,
       scale = 1, width = 13, height = 13, units ="cm",dpi = 300, limitsize = TRUE)