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. 

Thanks! 

-Deena 

> 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: pd.hugene.2.0.st
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
Normalizing
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 <- contrasts.fit(fit, 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

locale:
[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             

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)

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")