Question: Error in volcano plot drawing with R.
0
gravatar for Sib
8 weeks ago by
Sib20
Sib20 wrote:

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
R • 176 views
ADD COMMENTlink modified 8 weeks ago by serovop5410 • written 8 weeks ago by Sib20

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?

ADD REPLYlink written 8 weeks ago by RamRS25k

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 ...
ADD REPLYlink written 8 weeks ago by Sib20

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

ADD REPLYlink written 8 weeks ago by cpad011212k

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

ADD REPLYlink written 8 weeks ago by Sib20

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

ADD REPLYlink written 8 weeks ago by RamRS25k

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?

ADD REPLYlink written 7 weeks ago by Sib20

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

ADD REPLYlink modified 7 weeks ago • written 7 weeks ago by cpad011212k

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

ADD REPLYlink written 7 weeks ago by cpad011212k

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

ADD REPLYlink written 8 weeks ago by cpad011212k

In fit object, there are no p-values

ADD REPLYlink written 8 weeks ago by Sib20
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: 797 users visited in the last hour