Question: How to label custom multiple gene lists object into volcano plot?
1
gravatar for choijamtsm
7 months ago by
choijamtsm50
choijamtsm50 wrote:

Hello everyone!

I am new here and I tried to make a basic volcano plot with R.

#import data
res <- read.table("results.txt", header=TRUE)
genes <- read.table("genes.txt", header=TRUE)

# Make a basic volcano plot 
with(res, plot(log2foldchange, -log10(pvalue), pch=20, cex.lab=1.6, main="Volcano Plot", xlim=c(-8,8), col="grey"))

# Add colored points:
with(subset(res, pvalue<.05 & abs(log2foldchange)>1.5), points(log2foldchange, -log10(pvalue), pch=20, col="cyan"))

# Label points with the textxy function from the calibrate plot
library(calibrate)

#Label only 1 gene
with(subset(res, Gene == "WNT7A"), textxy(log2FoldChange, -log10(pvalue), labs=Gene, cex=.8))

#How to insert "genes" into volcano plot (overlap res)?

My question is how to add "genes" which contains multiple gene lists into res? (which should label only genes list)

Thank you so much!

rna-seq R gene • 336 views
ADD COMMENTlink modified 7 months ago by RamRS30k • written 7 months ago by choijamtsm50
1

Provide example input: res and genes objects.

Note, there is a dedicated package for volcanoes - EnhancedVolcano

ADD REPLYlink written 7 months ago by zx87549.6k

Thanks for your reply

  • Here is the input:

res <- read.table("results.txt", header=TRUE)

gene      log2foldchange    pvalue      padj
Cx3cr1      -8.039238787    1.29E-11     2.78E-116
Trem1       -5.258501816    3.01E-44       1.77E-42
Serpina3f  3.202817674  2.12E-09       2.40E-08
Fcrl1        -4.56714509    2.99E-16       5.91866E-15
AU021092 1.886580391    8.00E-03       0.031416887 and so on..

genes <- read.table("genes.txt", header=TRUE)

gene
Fcrl1
Trem1
Serpina3f and so on..

I am struggling with inject genes$gene into res volcano plot

Thank you

ADD REPLYlink modified 7 months ago • written 7 months ago by choijamtsm50

Run this:

if (!requireNamespace('BiocManager', quietly = TRUE))
    install.packages('BiocManager')

BiocManager::install('EnhancedVolcano')
require('EnhancedVolcano')

EnhancedVolcano(res,
    lab = res$gene,
    x = 'log2foldchange',
    y = 'pvalue')

More information in the vignette that I wrote: https://bioconductor.org/packages/release/bioc/vignettes/EnhancedVolcano/inst/doc/EnhancedVolcano.html

ADD REPLYlink modified 7 months ago • written 7 months ago by Kevin Blighe65k

Thank you so much for your reply. I downloaded 'EnhancedVolcano' package then run above command.

EnhancedVolcano(res,
lab = res$gene,
x = 'log2foldchange',
y = 'pvalue')

But All of the res$gene labeled on volcano plot.

Volcano plot Enhanced Volcano res$gene

Actually my goal is to label custom genes. So I've tried to run following command:

EnhancedVolcano(res,
            lab = res$gene,
            x = 'log2foldchange',
            y = 'pvalue',
            selectLab = gene_list)

But it gives me the following plot! which are only row number's shown? how can I see the gene name?

Trying to labeling custom genes

Thank you

ADD REPLYlink modified 7 months ago • written 7 months ago by choijamtsm50

But All of the res$gene labeled on volcano plot.

Volcano plot Enhanced Volcano res$gene

You simply need to alter the cut-offs for p-value and log [base 2] fold-change. Please read the manual or take a look at the vignette:

---------------------------------------

But it gives me the following plot! which are only row number's shown? how can I see the gene name?

Trying to labeling custom genes

Please show the outputs of the following:

class(res)
head(res)
head(res$gene)
class(gene_list)
head(gene_list)
ADD REPLYlink written 7 months ago by Kevin Blighe65k

Thank you,

1. I altered the cut-offs for p-value and log2 fold change and graph looks nice!

altered p-value and log2 fold change

2. Here are the outputs:

class(res)

[1] "data.frame"

head(res)

    gene log2foldchange    pvalue          padj
1    Cx3cr1      -8.039239 1.29e-118 2.780000e-116
2     Trem1      -5.258502  3.01e-44  1.770000e-42
3 Serpina3f       3.202818  2.12e-09  2.400000e-08
4     Fcrl1      -4.567145  2.99e-16  5.918660e-15
5  AU021092       1.886580  8.00e-03  3.141689e-02
6     Gbp11       2.090275  3.65e-14  6.250000e-13

head(res$gene)

[1] Cx3cr1    Trem1     Serpina3f Fcrl1     AU021092  Gbp11    
1594 Levels: 10-Mar 1190002N15Rik 1500017E21Rik 1700003F12Rik 1700009N14Rik 2610528A11Rik 2900026A02Rik 3425401B19Rik 3632451O06Rik ... Zmynd15

class(gene_list)

[1] "character"

head(gene_list)

[1] "Fam122A" "Ccnd2"   "Usp2"    "Arl2"    "Rab9B"   "Dcaf7"

My goal is to only label gene_list in res. I also tried following code but failed!

EnhancedVolcano(res,
            lab = res$gene,
            x = 'log2foldchange',
            y = 'pvalue',
            pCutoff = 10e-16,
            FCcutoff = 1.5,
            selectLab = res$gene[which(res$gene %in% gene_list)])

Thank you~

ADD REPLYlink modified 7 months ago • written 7 months ago by choijamtsm50
1

Oh, you just need to do this before running anything else:

res$gene <- as.character(res$gene)
ADD REPLYlink written 7 months ago by Kevin Blighe65k
1

It works now!! thank you so much

ADD REPLYlink written 7 months ago by choijamtsm50

Please use the formatting bar (especially the code option) to present your post better. You can use backticks for inline code (`text` becomes text), or select a chunk of text and use the highlighted button to format it as a code block. I've done it for you this time.
code_formatting

ADD REPLYlink written 7 months ago by RamRS30k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1622 users visited in the last hour