I am new to the use of microarray analysis; I have been looking at the use of lumi for the analysis of BeadArray chips.
I have used the tutorial by Du, Feng, Kibbe and Lin (2010) (http://www.bioconductor.org/packages//2.7/bioc/vignettes/lumi/inst/doc/lumi.pdf) to learn the process. In that vignette, the authors use the packages lumi and limma to normalize the data and then fit the linear method (creating an object 'fit') that can identify the most expressed genes. The data come from the example.lumi dataset contained in the package lumiBarnes.
I then used the tutorial by Gillespie (http://bioinformatics.knowledgeblog.org/2011/06/21/volcano-plots-of-microarray-data/) to learn how to draw a volcano plot. It essentially plots the fold change against the -log of the p-value. When I applied this approach to the data of the lumi tutorial, however, I did not obtain the classical 'V' shaped profile but rather a '√' shape. This 'shape alteration' might be inherent to the data, but I reckon it codued be do to the fitting method applied. Both the fitting method and the plotting procedure are based on the topTable function and this is derived from the model.matrix function, thus if the analysis has produced the observed skewness in the data is the model.matrix that I need to modify. In the vignette by Gillespie, the contrast.matrix and contrasts.fit(fit, contrast.matrix) are used instead of model.matrix.
My questions are then: (a) how can I customize the contrast.matrix? (b) how can I draw a volcano plot of the data? (maybe I simply need to change the fitting method to obtain the expected 'V' shape)
The code follows.
### # lumi analysis library(lumi) library(lumiBarnes) library(limma) data("example.lumi") lumi_t <-lumiT(example.lumi) lumi_n <- lumiN(lumi_t, method='rsn') dataMatrix <- exprs(lumi_n) probeList <- rownames(dataMatrix) sampleType <- c('100:0', '95:5', '100:0', '95:5') design <- model.matrix(~ factor(sampleType)) colnames(design) <- c('100:0', '95:5-100:0') fit <- lmFit(dataMatrix, design) fit <- eBayes(fit) # print most expressed genes topTable(fit, coef='95:5-100:0', adjust='fdr', number=10) # volcano plot gene_list <- topTable(fit, coef='95:5-100:0', number=1000000, sort.by="logFC") plot(gene_list$logFC, -log10(gene_list$P.Value))