Making a volcano plot look less cluttered

Question

I have generated a volcano plot with a differential expression file.

Code for inputing file:

library(tidyr)
final_tumor <- df %>% pivot_longer(., -c(Feature.ID, Feature.Name), names_to = c("set",".value"), names_pattern = "(.+)_(.+)")

# A tibble: 80 x 6
Feature.ID Feature.Name set       Mean.Counts Log2.fold.change Adjusted.p.value
<fct>      <fct>        <chr>           <dbl>            <dbl>            <dbl>
  1 a          A            Cluster.1    0.000961            0.292           1     
2 a          A            Cluster.2    0.000902            0.793           1     
3 a          A            Cluster.3    0.00181             1.46            0.758 
4 a          A            Cluster.4    0.000642            0.269           1     
5 b          B            Cluster.1    0.000320            1.95            0.910 
6 b          B            Cluster.2    0.00180             4.77            0.154 
7 b          B            Cluster.3    0                   2.19            1     
8 b          B            Cluster.4    0                   1.66            1     
9 c          C            Cluster.1    0.00128            -2.01            0.0467
10 c          C            Cluster.2    0.00632             0.352           1     
# … with 70 more rows

Here the output of head(final_tumor)

# A tibble: 6 x 6
Feature.ID        Feature.Name set       Mean.Counts Log2.fold.change Adjusted.p.value
<fct>             <fct>        <chr>           <dbl>            <dbl>            <dbl>
  1 ENSG00000227232.5 WASH7P       Cluster.1     0                   1.50            1    
2 ENSG00000227232.5 WASH7P       Cluster.2     0                   1.73            1    
3 ENSG00000227232.5 WASH7P       Cluster.3     0                   1.77            1    
4 ENSG00000227232.5 WASH7P       Cluster.4     0.00114             4.30            0.293
5 ENSG00000227232.5 WASH7P       Cluster.5     0                   2.15            1    
6 ENSG00000227232.5 WASH7P       Cluster.6     0                   1.22            1 

And here is the output of tail(final_tumor):

# A tibble: 6 x 6
Feature.ID        Feature.Name set        Mean.Counts Log2.fold.change Adjusted.p.value
<fct>             <fct>        <chr>            <dbl>            <dbl>            <dbl>
  1 ENSG00000210196.2 MT-TP        Cluster.6       0.0699          -0.202           0.790  
2 ENSG00000210196.2 MT-TP        Cluster.7       0.0801           0.0386          1      
3 ENSG00000210196.2 MT-TP        Cluster.8       0.0711           0.0875          1      
4 ENSG00000210196.2 MT-TP        Cluster.9       0.0152          -2.31            0.00127
5 ENSG00000210196.2 MT-TP        Cluster.10      0.0147          -2.30            0.00612
6 ENSG00000210196.2 MT-TP        Cluster.11      0.122            0.762           1

Here is the code for generating volcano plot:

library(ggplot2)
library(ggrepel)
ggplot(final_tumor, aes(x = Log2.fold.change,y = -log10(Adjusted.p.value), label = Feature.Name))+ 
  geom_point()+
  geom_text_repel(data = subset(final_tumor, Adjusted.p.value < 0.05), 
                  aes(label = Feature.Name))+
  geom_text_repel(data = subset(final_tumor, Feature.Name == "Casp14"), 
                  aes(label = Feature.Name), color = "red")

Generating the ggplot with all the data was taking forever, it was suggested to select random rows of the dataframe:

rep_final = final_tumor[sample(1:nrow(final_tumor), size = 1000), ]

Then I did:

ggplot(rep_final, aes(x = Log2.fold.change,y = -log10(Adjusted.p.value), label = Feature.Name))+ 
  geom_point()+ 
  geom_text_repel(data = subset(rep_final, Adjusted.p.value < 0.05), aes(label = Feature.Name))+ 
  geom_text_repel(data = subset(final_tumor, Feature.Name == "Casp14"), aes(label = Feature.Name), color = "red")

Everything is very tight and cluttered. I want to show the top differentially expressed genes. How do I go about simplifying the plot?

Is there a way to do this on EnhancedVolcano, pulling out a certain gene?

EnhancedVolcano(final_tumor, lab = as.character(final_tumor$FeatureName), x = 'Log2.fold.change', y = 'Adjusted.p.value', xlim = c(-8,8), title = 'Tumor', pCutoff = 10e-5, FCcutoff = 1.5, pointSize = 3.0, labSize = 3.0) 

Answer 1

The element that is clustering your plot the most is the ggrepel labels, so that's what you need to work on.

To get to a readable number of labels, I would not remove random genes: you risk losing important info, be it genes that could be of interest to you or some of the most regulated. Incidentally that seems to happen on your example: you try plotting in red the label for Casp14, but it does not show on the plot, probably because the gene was removed in your sampling.

What you can do however is, as dc37 suggests, increase the thresholds for the labels (you can probably play around with the thresholds to see when it does become readable). This should also help with the plotting: you mention that generating the plot was taking a long time, but with the number of genes you are looking at (~10k-30k?), the geom_points() should not be a real problem (maybe a few seconds), however the ggrepel is the most computationnally expensive so that's what you should try to reduce.

This should get you started:

ggplot(final_tumor, aes(x = Log2.fold.change,y = -log10(Adjusted.p.value), label = Feature.Name))+ 
  geom_point()+ 
  geom_text_repel(data = subset(final_tumor, Adjusted.p.value < 0.001 & abs(Log2.fold.change) > 1 ), aes(label = Feature.Name))+ 
  geom_text_repel(data = subset(final_tumor, Feature.Name == "Casp14"), aes(label = Feature.Name), color = "red")

Notice everything plots the "final_tumor" object and not the "rep_final". Feel free to be even more stringent with the cutoffs if it's still taking too long to plot.

A bit outside of the score of your question, but I also wanted to mention that you should never be trying to label too many genes on a volcano plot: this type of graph is useful to visualize your data's p-value/FC distribution and highlight a subset of genes in this context (most regulated, genes regulated in another experiment, pathway of interest...), but if you label too many of them it will never be readable.

If you really want to be able to indentify all the genes on your graph, I would try an interactive visualization (like with plotly). I'm not going to detail the code here but you could look in that direction! (However, as with ggrepel, you might need to remove some genes so that it does not "take forever").

Have fun :)

Answer 2

As we discussed in comments in your previous post: Pulling out a certain gene in a volcano plot, in order to have less labeled genes on your volcano plot, you can specify more stringent cutoff.

For example, you can decide to display only genes with an absolute fold change superior to 2 by selecting the appropriate subset of your data:

ggplot(rep_final, aes(x = Log2.fold.change,y = -log10(Adjusted.p.value), label = Feature.Name))+ 
  geom_point()+ 
  geom_text_repel(data = subset(rep_final, abs(Log2.fold.change) > 2), aes(label = Feature.Name))+ 
  geom_text_repel(data = subset(final_tumor, Feature.Name == "Casp14"), aes(label = Feature.Name), 
                  color = "red", size = 5)

Alternatively, you can maybe decide that a p value inferior to 0.05 is not enough stringent and decide to labeled only genes with an adjusted p value inferior to 0.00001 by doing:

ggplot(rep_final, aes(x = Log2.fold.change,y = -log10(Adjusted.p.value), label = Feature.Name))+ 
  geom_point()+ 
  geom_text_repel(data = subset(rep_final, Adjusted.p.value < 0.00001), aes(label = Feature.Name))+ 
  geom_text_repel(data = subset(final_tumor, Feature.Name == "Casp14"), aes(label = Feature.Name), 
                  color = "red", size = 5)

Or, you can want to do a combination of both criteria by doing:

ggplot(rep_final, aes(x = Log2.fold.change,y = -log10(Adjusted.p.value), label = Feature.Name))+ 
  geom_point()+ 
  geom_text_repel(data = subset(rep_final, abs(Log2.fold.change) > 2 & Adjusted.p.value < 0.00001), aes(label = Feature.Name))+ 
  geom_text_repel(data = subset(final_tumor, Feature.Name == "Casp14"), aes(label = Feature.Name), 
                  color = "red", size = 5)

At the end, I think it depends of what you are trying to demonstrate by presenting this volcano plot. Do you want to highlight some particular genes ? or do you want to display genes that are the most differentially expressed ? If so, how many ?

All of this depends of the scientific question underlying your experiments and why you came to make this volcano plot.

Hope this couple of example helps you to figure it out how to customize your volcano plot.