This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How to highlight certain genes in a Volcano Plot in R

Hi all,

I would like to know please how to highlight certain genes in a list with a specific colour in a Volcano Plot in R

I have the following code:

genes22= read.csv("../Python/matrix_gene_final v2.csv")
gene_list <- as.vector(genes22)

#All Genes
# MOCS vs MOSE
MOCSvsMOSE = read.csv("../Data/Comparisons/G1 - G2 MOCS vs MOSE.csv") # read csv
res = MOCSvsMOSE

# Make a basic volcano plot
with(res, plot(logFC, -log10(P.Value),pch=20, main="MOCS vs MOS", xlim=c(-8.8,8)))

# if adj.P value less than 0.05 = green
with(subset(res, adj.P.Val<.05 ), points(logFC, -log10(P.Value), pch=20, col="green"))

# if logFC value more than 1 = blue
with(subset(res, abs(logFC)>1), points(logFC, -log10(P.Value), pch=20, col="blue"))

# show the genes in the gene list as pink in the volcano plot    
with(res [rownames(res) %in% gene_list]), points(logFC, -log10(P.Value), pch=20, col="pink"))

Everything works except the final line of code:

 # show the genes in the gene list as pink in the volcano plot    
   with(res [rownames(res) %in% gene_list]), points(logFC, -log10(P.Value), pch=20, col="pink"))

I get the following error:

Error: unexpected ',' in "with(res [rownames(res) %in% gene_list]),"

The volcano plot only shows the blue and green colours.

How can I fix this please?

Many Thanks,

Ishack

r volcano plot deg
• 11,516 views • link
updated by 1769mkc 1769mkc • written by ishackm ishackm

This is purely a programming question. Hint: df [1:10,1:10] is not the same as df[1:10,1:10].

• 0 views • link
written by Ram Ram

Sorry i dont understand. Can you clarify please?

• 0 views • link
written by ishackm ishackm

I wrote two lines/statements of code as a hint so you could find the difference between the two statements and apply your learning to your code.

• 0 views • link
written by Ram Ram

Thanks but i still get the same error.

Is there another way?

• 0 views • link
written by ishackm ishackm

Yes, check if your parentheses are matched. Again, this is a simple syntax error that has nothing to do with bioinformatics.

• 0 views • link
written by Ram Ram

as Ram has been pointing out: the problem is in your last line, specifically the way you subset the data.frame -- why don't you use the same subset routine as for the other lines?

• 0 views • link
written by Friederike Friederike

Hi I have added this code and there is no error but I still don't get the orange colour 

gene_list <- scan("../Python/matrix_gene_final.csv", what="", sep="\n")
gene_list <- as.vector(gene_list) 

with(subset(res, rownames(res) %in% gene_list), points(logFC, -log10(P.Value), pch=20, col="orange"))

How can I fix this, please?

Many Thanks,

Ishack

• 0 views • link
written by ishackm ishackm

well, what do you get? Without an error or a plot it's difficult to diagnose

• 0 views • link
written by Friederike Friederike

I recommend you to use an online volcano-plot generator like this.

• 0 views • link
written by Buffo Buffo

2 answers

I wrote this code for similar purpose, see if it helps for better visualization.

library(ggplot2)
library(ggrepel)
library(dplyr)
gene_list <- as.vector(genes22)

data= subset(res, rownames(res) %in% gene_list)
data= data[which(abs(data$logFC) >=1),]
data= data[which(abs(data$P.Value) <=0.05),]
res <- res %>%
    mutate(reg = case_when(
    res$logFC >= 1 & res$P.Value <= 0.05 ~ "UP",
    res$logFC <= -1 & res$P.Value <= 0.05 ~ "DOWN",
    abs(res$logFC) <= 1 & res$P.Value >= 0.05 ~ "no_change",
    abs(res$logFC) <= 1 & res$P.Value <= 0.05 ~ "no_change",
    abs(res$logFC) > 1 & res$P.Value >0.05 ~ "no_change"
    )) %>%
  mutate(reg = factor(reg, levels = c("UP", "no_change","DOWN")))

plt <-  ggplot(res,aes(x=logFC,y=-log10(P.Value),label = rownames(res))+
          geom_point(aes(color=reg))+
          scale_color_manual(name = "Differential \n regulation",
                     values = c("DOWN" = "#058983",
                                "no_change" = "grey",
                                "UP" = "#DEB132"),
                     labels = c("UP", "No Change", "DOWN"))+
          theme(axis.text.x=element_text(size=16,color="black")
                  , axis.title.x=element_text(size=16)
                  , axis.text.y=element_text(size=16,color="black")
                  , axis.title.y=element_text(size=16)
                  , legend.justification=c(0,1)
                  , legend.box.background = element_rect(colour = "black")
                  , plot.title = element_text(hjust = 0.5))+
          geom_text_repel(
                    data          = data,
                    size          = 5,
                    nudge_y      = 30,
                    direction    = "x",
                    angle        = 0,
                    vjust        = 0,
                    segment.size  = 0.5,
                    segment.color = "black"
                    direction     = "x"
)
• 0 views • link
written by Prakash Prakash 
library(EnhancedVolcano)
res1 <- read.table("your_file_with_gene_fc_pvalue",header = TRUE,row.names = 1,sep = "\t")

EnhancedVolcano(res1,lab = rownames(res1),x = "log2FoldChange",y = "padj",
                #selectLab = c("RPL11","RPL10","RPL18A","RPL28"), [just add your list of genes you want to highlight]
                xlim = c(-6, 6),
                xlab = bquote(~Log[2]~ "fold change"),
              ylab = bquote(~-Log[10]~adjusted~italic(P)),
                transcriptPointSize = 10,
                transcriptLabSize = 10,
              border = "full",
              #legendPosition = "bottom",
              borderWidth = 1.5,
              legend=c('NS','Log2 FC','Adjusted p-value',
                       'Adjusted p-value & Log2 FC'),
              legendPosition = 'bottom',
              legendLabSize = 20,
              legendIconSize = 20,
              borderColour = "blue",
              drawConnectors = FALSE,
              #widthConnectors = 0.01,
              colConnectors = 'grey30',
              gridlines.major = FALSE,
              gridlines.minor = FALSE)
• 0 views • link
written by 1769mkc 1769mkc

Thank you all, it finally works.

• 0 views • link
written by ishackm ishackm

Log in to answer this question.

More posts like this