Volcano plot in R

Question

This question has also been asked on biostars

How can I reproduce this volcano plot?

I'm only able to do the traditional one, I'm kind knew too these field.

Answer 1

This plot is clearly done using core R functions. There are smoother alternatives how to make a pretty volcano plot (like ggplot with example here), but if you really wish to, here is my attempt to reproduce it :

I obviously had to generate data since I do not have the expression data from the figure, but the procedure will be about the same with the real data. Here is the R script :

### define interesting genes
interesting_genes <- c('Il17s', 'Casz1', 'Rorc', 'Il22')
fc_interesting_genes <- c(-1.8, -1.3, -0.9, -0.85)
pval_interesting_genes <- log10(c(0.0005, 0.01, 1e-4, 4e-05))

### define colours (grey, blue, red)
palette <- c(rgb(0, 0, 0, max = 255, alpha = 65),
    rgb(0, 0, 255, max = 255, alpha = 125),
    rgb(255, 0, 0, max = 255, alpha = 125))

### generate some expression like data, fold change using normal distr, and tranformed uniform for pvalues
fold_changes <- c(rnorm(1000, 0, 2), fc_interesting_genes)
pvalues <- c(log10(runif(1000) / 3^(1 + abs(fold_changes[1:1000]))), pval_interesting_genes)

# find deferentially expressed genes to colour
left_biased <- fold_changes < -1
right_biased <- fold_changes > 1

point_type = 20
point_size = 1.5

### plot corpus
plot(fold_changes, pvalues,
     pch = point_type, cex = point_size, col = palette[1],
     xlab = "Fold Change (log2)", ylab = "P value",
     xlim = c(-6,6), ylim = c(-5, 0), yaxt="n")
aty <- axTicks(2)
labels <- sapply(aty,function(i)
            as.expression(bquote(10^ .(i)))
          )
labels[6] <- 1
axis(2,at=aty,labels=labels)

### plot deferentially expressed genes in colour + genes on interest with black
points(fold_changes[left_biased], pvalues[left_biased],
       pch = point_type, cex = point_size, col = palette[2])
points(fold_changes[right_biased], pvalues[right_biased],
       pch = point_type, cex = point_size, col = palette[3])
points(fc_interesting_genes, pval_interesting_genes, pch = point_type, cex = point_size)

### add text
labelx <- c(-4, -4, 4, 4)
labely <- c(-1.5, -0.5, -0.5, -1.5)
labels <- c('Il17s', 'Casz1', 'Rorc', 'Il22')
text(labelx, labely, labels)
text(-4, -4, "down\n194 genes")
text(4, -4, "up\n156 genes")

### add lines between text and data points
for(i in 1:4){
    lines(c(labelx[i], fc_interesting_genes[i]), c(labely[i], pval_interesting_genes[i]))
}

Answer 2

Without a minimal reproducible example (MWE), I can't reproduce the plot but I would suggest using existing volcanoplot functions such as the limma package on Bioconductor. Here is a typical workflow for a differential expression analysis that produces a violinplot:

#install package
source("https://bioconductor.org/biocLite.R")
biocLite("limma")

#load package
library("limma")

#define phenotype/experimental conditions for samples
group = rep(c("Treatment","Control"), 
c(number_of_treated_samples,number_of_controls))
design = model.matrix(~group);
colnames(design) = c("Control","TreatmentvsControl")

#fit empirical Bayes model
fit <- lmFit(NormalisedExpression, design)
fit <- eBayes(fit)

#summarise significant results
tt <- topTable(fit, coef="TreatmentvsControl",
                adjust="BH",n=nrow(NormalisedExpression))

#plot results
volcanoplot(fit, coef="TreatmentvsControl")

This plot can be customised in a similar manner to base R plots by passing the relevant arguments as shown in the limma documentation. This is the recommended plot format that readers in the field will be familiar with. While I do not recommend deviating from this standard, it is possible to reverse axes in R and the source code for the plot function can be retrieved from R by entering volcanoplot (without any arguments) if you wish to modify it.

Answer 3

Easy to visualize volcano plot using Python. Check this tutorial for inverted volcano plot