Limma Volcano Plot Highlight Top Genes by Color
2
0
Entering edit mode
dmm204 • 0
@dmm204-15412
Last seen 5.5 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. 

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             

 
limma microarray volcanoplot color • 6.9k views
ADD COMMENT •
0
Entering edit mode

You can also try out EnhancedVolcano (Bioconductor): https://github.com/kevinblighe/EnhancedVolcano

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.

Kevin

ADD REPLY •
0
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.

ADD REPLY •
0
Entering edit mode
@james-w-macdonald-5106
Last seen 2 hours 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)
ADD COMMENT •
0
Entering edit mode
@gordon-smyth
Last seen 8 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")

 
ADD COMMENT •

Login before adding your answer.

Similar Posts
Loading Similar Posts
Traffic: 957 users visited in the last hour
Content Search
Users
Tags
Badges
Help About
FAQ
Access RSS
API
Stats

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

Powered by the version 2.3.6