Error in volcano plot drawing with R.

Question

[Sib]

I want to draw a volcano plot. I don't know where I make a mistake that I get an error. These are the codes:

library(Biobase)
library(GEOquery)
library(ggplot2)
library(pheatmap)
library(limma)
library(plyr)
gset<-getGEO("GSE116959",GSEMatrix = TRUE,AnnotGPL = TRUE,destdir ="data/")
gset<-gset[[1]]
gr<-c(rep("Tumor",3),("Normal"),rep("Tumor",7),("Normal"),rep("Tumor",21),rep("Normal",4)
      ,rep("Tumor",5),("Normal"),rep("Tumor",17),rep("Normal",2),("Tumor"),("Normal"),rep("Tumor",2),("Normal"),("Tumor"))
gr<-as.factor(gr)
gset$description<-gr
design<-model.matrix(~description+0,gset)
colnames(design)<-levels(gr)
fit<-lmFit(gset,design)
contmatrix<-makeContrasts(Tumor-Normal,levels = design)
fit2<-contrasts.fit(fit,contmatrix)
fit2<-eBayes(fit2,0.01)
tT <- topTable(fit2, adjust="fdr", sort.by="B", number=Inf)
volcanoplot(fit, coef = 1, style = "p-value", highlight = 0,  hl.col="blue",
            xlab = "Log2 Fold Change", ylab = NULL, pch=16, cex=0.35)
Error in volcanoplot(fit, coef = 1, style = "p-value", highlight = 0,  : 
  No p-values found in linear model fit object

[Ram]

What have you tried? Have you looked into the fit object and the lmFit method and figured out why the object does not have p values? Is it supposed to have p values?

[Sib]

I am new to this analysis and I don't know whether It should have p-value or not But when I call fit, the results are these:

> fit
An object of class "MArrayLM"
$coefficients
                 Normal    Tumor
A_19_P00315452 7.134198 6.256601
A_19_P00315459 6.925040 7.046574
A_19_P00315482 6.000868 6.067685
A_19_P00315492 5.661245 5.743756
A_19_P00315493 7.579796 7.911303
50594 more rows ...

$rank
[1] 2

$assign
[1] 1 1

$qr
$qr
               Normal       Tumor
GSM3265243 -3.3166248  0.00000000
GSM3265244  0.0000000 -7.54983444
GSM3265245  0.0000000  0.13245324
GSM3265246  0.3015113 -0.03993615
GSM3265247  0.0000000  0.13245324
63 more rows ...

$qraux
[1] 1.000000 1.132453

$pivot
[1] 1 2

$tol
[1] 1e-07

$rank
[1] 2


$df.residual
[1] 66 66 66 66 66
50594 more elements ...

$sigma
A_19_P00315452 A_19_P00315459 A_19_P00315482 A_19_P00315492 A_19_P00315493 
     0.6351571      0.5536146      0.3861076      0.2579868      0.6370429 
50594 more elements ...

$cov.coefficients
           Normal      Tumor
Normal 0.09090909 0.00000000
Tumor  0.00000000 0.01754386

$stdev.unscaled
                  Normal     Tumor
A_19_P00315452 0.3015113 0.1324532
A_19_P00315459 0.3015113 0.1324532
A_19_P00315482 0.3015113 0.1324532
A_19_P00315492 0.3015113 0.1324532
A_19_P00315493 0.3015113 0.1324532
50594 more rows ...

$pivot
[1] 1 2

$genes
                           ID                              SPOT_ID CONTROL_TYPE    REFSEQ    GB_ACC LOCUSLINK_ID  GENE_SYMBOL                    GENE_NAME
A_19_P00315452 A_19_P00315452                       A_19_P00315452        FALSE XR_110148 XR_110148           NA LOC100130938 uncharacterized LOC100130938
A_19_P00315459 A_19_P00315459                       A_19_P00315459        FALSE              U66048           NA                                          
A_19_P00315482 A_19_P00315482                       A_19_P00315482        FALSE NR_034143 NR_034143           NA    LOC729177    uncharacterized LOC729177
A_19_P00315492 A_19_P00315492 tc|THC2614189|linc|TCONS_l2_00021703        FALSE                               NA       Q73P46                             
A_19_P00315493 A_19_P00315493                       A_19_P00315493        FALSE NR_027046 NR_027046           NA    LOC145474    uncharacterized LOC145474
               UNIGENE_ID ENSEMBL_ID                                ACCESSION_STRING     CHROMOSOMAL_LOCATION    CYTOBAND
A_19_P00315452                                           ref|XR_110148|ref|XR_111793  chr18:59273360-59273301 hs|18q21.33
A_19_P00315459                       gb|U66048|gb|U66046|tc|THC2721474|tc|THC2602433 chrX:149131107-149131166     hs|Xq28
A_19_P00315482  Hs.730673                    ref|NR_034143|gb|AK126168|tc|THC2510121   chr6:22136556-22136497   hs|6p22.3
A_19_P00315492                                  tc|THC2614189|linc|TCONS_l2_00021703 chr4:129376376-129376435   hs|4q28.2
A_19_P00315493  Hs.658241                    ref|NR_027046|gb|AK023718|tc|THC2606363  chr14:71956255-71956314  hs|14q24.2
                                                                                                                DESCRIPTION GO_ID
A_19_P00315452                        PREDICTED: Homo sapiens hypothetical LOC100130938 (LOC100130938), miscRNA [XR_110148]      
A_19_P00315459                                                                                                                   
A_19_P00315482                               Homo sapiens uncharacterized LOC729177 (LOC729177), non-coding RNA [NR_034143]      
A_19_P00315492 Q73P46_TREDE (Q73P46) Branched-chain amino acid ABC transporter, permease protein, partial (5%) [THC2614189]      
A_19_P00315493                               Homo sapiens uncharacterized LOC145474 (LOC145474), non-coding RNA [NR_027046]      
                                                                   SEQUENCE

A_19_P00315493 ATATATTCTTCAGTGTGATTTAGGTGAGTCAGGAATTTTTCCTCTGGAGACTGCCTAAAT
50594 more rows ...

$Amean
A_19_P00315452 A_19_P00315459 A_19_P00315482 A_19_P00315492 A_19_P00315493 
      6.398565       7.026914       6.056876       5.730408       7.857677 
50594 more elements ...

$method
[1] "ls"

$design
           Normal Tumor
GSM3265243      0     1
GSM3265244      0     1
GSM3265245      0     1
GSM3265246      1     0
GSM3265247      0     1
63 more rows ...

[cpad0112]

object tT will have logFC, P.value and adj.p.values. Try to plot logfc against p.value/adj.p.value. @ maryammomeni23

[Sib]

thanks for answering. I will be appreciated more if you write the codes that I should use for doing this.

[Ram]

Please explore that on your own. cpad has already given you the exact location of the p values that you need.

[Sib]

thanks. I plotted a volcano plot by ggplot2 . the codes are as below:

ggplot(tT,aes(x=logFC,y=-log10(adj.P.Val)))+geom_point()

Now, I want to change the colors of genes with |logFC|>2 and adj-P-Val<0.05. And I want to show their gene-symbols. What codes I should add to ggplots codes?

[cpad0112]

you need to create a subset of data for plotting points with different colors and labels. Refer to example code in ggplot2 section in https://digibio.blogspot.com/2018/05/volcano-plot-with-ggplot2-and-basic.html

[cpad0112]

or you can use enhanced volcano plot package by Kevin @ Sib

[cpad0112]

Did you look at the p values in fit object for coefficient 1?

[Sib]

In fit object, there are no p-values

Answer 1

[Gordon Smyth]

The error is because you ran volcanoplot on fit instead of on fit2. The eBayes function addes p-values to the object and fit2 is where you stored the eBayes output.