Limma Volcano Plot Highlight Top Genes by Color
Entering edit mode
dmm204 • 0
Last seen 4.0 years ago

Hi Everyone, 

I am trying to highlight (by change of color) the top DEG in my volcano plot from differential expression analysis using the limma package. However, the codes that I am using, change all of the points, instead of just the highlighted points, in the plot to the indicated color.  Can anymore tell me how to rework my command so only the highlighted genes are a different color? I know its probably something silly with my coding, but any help would be appreciated! I have listed the command codes and session info below. 



> library(oligo)
> library(limma)
> getwd()

[1] "/Users/dmaurer/Desktop/Butterfield Lab Data/Microarray 09-021 Patient Raw Files/CEL files - Patient #/Raw Files"
> Celfiles<-list.celfiles("/Users/dmaurer/Desktop/Butterfield Lab Data/Microarray 09-021 Patient Raw Files/CEL files - Patient #/Raw Files")
> Rawdata<-read.celfiles(Celfiles)
Loading required package:
Loading required package: RSQLite
Loading required package: DBI
Platform design info loaded.
Reading in : iDC_10.CEL
Reading in : iDC_11.CEL
Reading in : iDC_12.CEL
Reading in : iDC_13.CEL
Reading in : iDC_14.CEL
Reading in : iDC_15.CEL
Reading in : iDC_16.CEL
Reading in : iDC_17.CEL
Reading in : iDC_18.CEL
Reading in : iDC_19.CEL
Reading in : iDC_2.CEL
Reading in : iDC_20.CEL
Reading in : iDC_21.CEL
Reading in : iDC_22.CEL
Reading in : iDC_23.CEL
Reading in : iDC_24.CEL
Reading in : iDC_25.CEL
Reading in : iDC_26.CEL
Reading in : iDC_27.CEL
Reading in : iDC_28.CEL
Reading in : iDC_29.CEL
Reading in : iDC_3.CEL
Reading in : iDC_30.CEL
Reading in : iDC_31.CEL
Reading in : iDC_32.CEL
Reading in : iDC_33.CEL
Reading in : iDC_34.CEL
Reading in : iDC_35.CEL
Reading in : iDC_4.CEL
Reading in : iDC_5.CEL
Reading in : iDC_6.CEL
Reading in : iDC_7.CEL
Reading in : iDC_8.CEL
Reading in : iDC_9.CEL
Reading in : mDC_10.CEL
Reading in : mDC_11.CEL
Reading in : mDC_12.CEL
Reading in : mDC_13.CEL
Reading in : mDC_14.CEL
Reading in : mDC_15.CEL
Reading in : mDC_16.CEL
Reading in : mDC_17.CEL
Reading in : mDC_18.CEL
Reading in : mDC_19.CEL
Reading in : mDC_2.CEL
Reading in : mDC_20.CEL
Reading in : mDC_21.CEL
Reading in : mDC_22.CEL
Reading in : mDC_23.CEL
Reading in : mDC_24.CEL
Reading in : mDC_25.CEL
Reading in : mDC_26.CEL
Reading in : mDC_27.CEL
Reading in : mDC_28.CEL
Reading in : mDC_29.CEL
Reading in : mDC_3.CEL
Reading in : mDC_30.CEL
Reading in : mDC_31.CEL
Reading in : mDC_32.CEL
Reading in : mDC_33.CEL
Reading in : mDC_34.CEL
Reading in : mDC_35.CEL
Reading in : mDC_4.CEL
Reading in : mDC_5.CEL
Reading in : mDC_6.CEL
Reading in : mDC_7.CEL
Reading in : mDC_8.CEL
Reading in : mDC_9.CEL
Reading in : TMM2_10.CEL
Reading in : TMM2_11.CEL
Reading in : TMM2_12.CEL
Reading in : TMM2_13.CEL
Reading in : TMM2_14.CEL
Reading in : TMM2_15.CEL
Reading in : TMM2_16.CEL
Reading in : TMM2_17.CEL
Reading in : TMM2_18.CEL
Reading in : TMM2_19.CEL
Reading in : TMM2_2.CEL
Reading in : TMM2_20_.CEL
Reading in : TMM2_21.CEL
Reading in : TMM2_22.CEL
Reading in : TMM2_23.CEL
Reading in : TMM2_24.CEL
Reading in : TMM2_25.CEL
Reading in : TMM2_26.CEL
Reading in : TMM2_27.CEL
Reading in : TMM2_28.CEL
Reading in : TMM2_29.CEL
Reading in : TMM2_3.CEL
Reading in : TMM2_30.CEL
Reading in : TMM2_31.CEL
Reading in : TMM2_32.CEL
Reading in : TMM2_33.CEL
Reading in : TMM2_34.CEL
Reading in : TMM2_35.CEL
Reading in : TMM2_4.CEL
Reading in : TMM2_5.CEL
Reading in : TMM2_6.CEL
Reading in : TMM2_7.CEL
Reading in : TMM2_8.CEL
Reading in : TMM2_9.CEL
> eset<-rma(Rawdata)
Background correcting
Calculating Expression
> sample <-factor(rep(c("iDC","mDC","TMM2"), each =34))
> design.mat <-model.matrix (~0 +sample)
> colnames(design.mat) <- levels(sample)

> fit <- lmFit(eset, design.mat)

> contrast.mat <- makeContrasts(Diff.1 = mDC-iDC,Diff.2 = TMM2-iDC, Diff.3=TMM2-mDC, levels=design.mat)

>fit2 <-, contrast.mat)
> fit2 <- eBayes(fit2)


>volcanoplot(fit2,coef="Diff.1",highlight=10,names=fit$genes$NAME, col="red",main="mDC vs iDC") 

* This command results in all the points being red, not just the highlighted ones. 


Session Info -

R version 3.4.4 (2018-03-15)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.6

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib

[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] pd.hugene.2.0.st_3.14.1 DBI_1.0.0               RSQLite_2.1.1           limma_3.34.9            oligo_1.42.0            Biostrings_2.46.0       XVector_0.18.0         
 [8] IRanges_2.12.0          S4Vectors_0.16.0        Biobase_2.38.0          oligoClasses_1.40.0     BiocGenerics_0.24.0    

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.19               compiler_3.4.4             BiocInstaller_1.28.0       GenomeInfoDb_1.14.0        bitops_1.0-6               iterators_1.0.10           tools_3.4.4               
 [8] zlibbioc_1.24.0            digest_0.6.17              bit_1.1-14                 memoise_1.1.0              preprocessCore_1.40.0      lattice_0.20-35            ff_2.2-14                 
[15] pkgconfig_2.0.2            Matrix_1.2-14              foreach_1.4.4              DelayedArray_0.4.1         GenomeInfoDbData_1.0.0     affxparser_1.50.0          bit64_0.9-7               
[22] grid_3.4.4                 blob_1.1.1                 codetools_0.2-15           matrixStats_0.54.0         GenomicRanges_1.30.3       splines_3.4.4              SummarizedExperiment_1.8.1
[29] RCurl_1.95-4.11            affyio_1.48.0             


limma microarray volcanoplot color • 4.5k views
Entering edit mode

You can also try out EnhancedVolcano (Bioconductor):

To shade any gene or group of genes by any colour that you want, take a look at the end of the vignette. You can also label any specific gene or genes.


Entering edit mode

There's no need to print out the whole design matrix in this case, so I took it out to make your question a bit shorter and easier to read.

Entering edit mode
Last seen 2 days ago
United States

You can just add your own, using the points function.

volcanoplot(fit2,coef="Diff.1", names=fit$genes$NAME, main="mDC vs iDC") 
lp <- -log10(fit2$p.value[,"Diff.1"])
ord <- order(lp, decreasing = TRUE)[1:10]
points(fit2$coef[ord,"Diff.1"], lp[ord], pch = 16, cex = 0.45, col = 2)
Entering edit mode
Last seen 2 hours ago
WEHI, Melbourne, Australia

No, you haven't made a mistake. volcanoplot() doesn't give you a choice as to the color of the gene names on the plot. I will put add an argument to give you a choice in the next Bioconductor release.

Continuing on James MacDonald's code example, you can easily add the gene names yourself, exactly as is done by volcanoplot(), but now in red:

volcanoplot(fit2, coef="Diff.1", main="mDC vs iDC") 
lp <- -log10(fit2$p.value[,"Diff.1"])
ord <- order(lp, decreasing = TRUE)[1:10]
x <- fit2$coef[ord,"Diff.1"]
y <- lp[ord]
text(x, y, fit2$genes$NAME[ord], cex=0.8, col="red")



Login before adding your answer.

Traffic: 442 users visited in the last hour
Help About
Access RSS

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

Powered by the version 2.3.6