How to highlight certain genes in a Volcano Plot in R
2
1
Entering edit mode
5.6 years ago
ishackm ▴ 110

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 • 9.5k views
ADD COMMENT
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

Sorry i dont understand. Can you clarify please?

ADD REPLY
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

Thanks but i still get the same error.

Is there another way?

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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?

ADD REPLY
0
Entering edit mode

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

ADD REPLY
1
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY
2
Entering edit mode
5.6 years ago
Prakash ★ 2.2k

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"
)
ADD COMMENT
2
Entering edit mode
5.6 years ago
1769mkc ★ 1.2k
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)
ADD COMMENT
1
Entering edit mode

Thank you all, it finally works.

ADD REPLY

Login before adding your answer.

Traffic: 1810 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6