7a-volcano-plot.R

getwd()

## [1] "D:/Bioinformatics learning/nanostring/volcano plot"

setwd("D:/Bioinformatics learning/nanostring/volcano plot")
# read data
data = read.csv("D:/Bioinformatics learning/nanostring/volcano plot/7a unstim vs stim.csv", header = T)
head(data)

##      Gene Unstim   Stim Log2FC     pvalue
## 1     A2M   7.71  13.27  -0.78 0.44014713
## 2   ABCB1 237.03  67.69   1.81 0.00014744
## 3    ABL1 207.97 122.60   0.76 0.01713104
## 4     ADA 375.75 309.05   0.28 0.21649210
## 5 ADORA2A 122.55 135.95  -0.15 0.06876018
## 6   AICDA   4.04  89.89  -4.48 0.00148318

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

## ✔ ggplot2 3.3.6     ✔ purrr   0.3.4
## ✔ tibble  3.1.7     ✔ dplyr   1.0.9
## ✔ tidyr   1.2.0     ✔ stringr 1.4.0
## ✔ readr   2.1.2     ✔ forcats 0.5.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

library(ggplot2)
library(ggrepel)
#create a column Expression in data 
data = data %>% 
  mutate(Expression = 
           case_when(Log2FC >=1 & pvalue < 0.05 ~ "Up regulated",
                     Log2FC <= -1 & pvalue < 0.05 ~ "Down regulate",
                     TRUE ~ "Unchanged"))
head(data)

##      Gene Unstim   Stim Log2FC     pvalue    Expression
## 1     A2M   7.71  13.27  -0.78 0.44014713     Unchanged
## 2   ABCB1 237.03  67.69   1.81 0.00014744  Up regulated
## 3    ABL1 207.97 122.60   0.76 0.01713104     Unchanged
## 4     ADA 375.75 309.05   0.28 0.21649210     Unchanged
## 5 ADORA2A 122.55 135.95  -0.15 0.06876018     Unchanged
## 6   AICDA   4.04  89.89  -4.48 0.00148318 Down regulate

# Plot volcano plot
volcano1 = ggplot(data, aes(Log2FC,-log(pvalue, 10))) + geom_point(aes(color = Expression), size = 1/5)+
  xlab(expression("log"[2]*"FC")) + ylab(expression(-"log"[10]*"plavue")) +
  scale_color_manual(values = c("dodgerblue3", "gray50", "firebrick3")) + 
  guides(color = guide_legend(override.aes = list(size = 1.5))) + theme_classic()
volcano1

#adding label to top 10 genes

up = data %>% 
    filter(Expression == "Up regulated") %>%
    arrange(pvalue, desc(abs(Log2FC))) %>%
    head(10)

down = data %>%
    filter(Expression == "Down regulate") %>%
    arrange(pvalue, desc(abs(Log2FC))) %>%
    head(10)
top_genes_pvalue = bind_rows(up, down)
                                    
volcano = volcano1 + geom_text_repel(data = top_genes_pvalue, 
                                      mapping = aes(Log2FC, -log(pvalue,10), label = Gene), max.overlaps  =20,
                                      size = 1)
volcano

ggsave("volcano.tiff", dpi =1200)

## Saving 7 x 5 in image