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

<pre><code>#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 &lt;- lmFit(NormalisedExpression, design)
fit &lt;- eBayes(fit)

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

#plot results
volcanoplot(fit, coef="TreatmentvsControl")
</code></pre>

<p>This plot can be customised in a similar manner to base R plots by passing the relevant arguments as shown in the <a href="https://www.bioconductor.org/packages/devel/bioc/manuals/limma/man/limma.pdf" rel="nofollow noreferrer">limma documentation</a>. 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 <code>volcanoplot</code> (without any arguments) if you wish to modify it.</p>

CC BY-SA 4.0
2018-05-17T08:58:40.483Z
2018-05-22T02:47:26.903Z
2018-05-22T02:47:26.903Z
Tom Kelly ケリー・トム
<p>Easy to visualize volcano plot using Python. Check <a href="https://reneshbedre.github.io/blog/volcano.html" rel="nofollow noreferrer">this tutorial</a> for inverted volcano plot</p>

CC BY-SA 4.0
2018-12-06T19:25:05.160Z
2018-12-07T15:13:09.043Z
2018-12-07T15:13:09.043Z
Ren
<p>This plot is clearly done using core R functions. There are smoother alternatives how to make a pretty volcano plot (like <code>ggplot</code> with example <a href="https://www.biostars.org/p/309962/" rel="noreferrer">here</a>), but if you really wish to, here is my attempt to reproduce it :</p>

<p><a href="https://i.sstatic.net/Key5z.png" rel="noreferrer"><img src="https://i.sstatic.net/Key5z.png" alt="volcano plot"></a></p>

<p>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 :</p>

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

### define colours (grey, blue, red)
palette &lt;- 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 &lt;- c(rnorm(1000, 0, 2), fc_interesting_genes)
pvalues &lt;- c(log10(runif(1000) / 3^(1 + abs(fold_changes[1:1000]))), pval_interesting_genes)

# find deferentially expressed genes to colour
left_biased &lt;- fold_changes &lt; -1
right_biased &lt;- fold_changes &gt; 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 &lt;- axTicks(2)
labels &lt;- sapply(aty,function(i)
            as.expression(bquote(10^ .(i)))
          )
labels[6] &lt;- 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 &lt;- c(-4, -4, 4, 4)
labely &lt;- c(-1.5, -0.5, -0.5, -1.5)
labels &lt;- 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]))
}
</code></pre>

CC BY-SA 3.0
2018-04-18T19:48:47.460Z
2018-04-18T19:48:47.460Z
Kamil S Jaron
<p><em>This question has also been asked on <a href="https://www.biostars.org/p/309962/" rel="nofollow noreferrer">biostars</a></em></p>

<p>How can I reproduce this volcano plot?</p>

<p><a href="https://i.sstatic.net/KD9Eh.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/KD9Eh.png" alt="Volcano plot"></a></p>

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

<rna-seq><r><visualization>
Volcano plot in R
# Volcano plot in R

(Asked by: username_0 on 2018-04-18)

*This question has also been asked on [biostars](https://www.biostars.org/p/309962/)*

How can I reproduce this volcano plot?

[![Volcano plot](https://i.sstatic.net/KD9Eh.png)](https://i.sstatic.net/KD9Eh.png)

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

**Tags:** rna-seq, r, visualization

---

## 3 answers

### Answer by username_3 on 2018-04-18. Score: 13 | Accepted answer

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 by username_1 on 2018-05-17. Score: 3

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 by username_2 on 2018-12-06. Score: 0

Easy to visualize volcano plot using Python. Check [this tutorial](https://reneshbedre.github.io/blog/volcano.html) for inverted volcano plot


---

