"bioconductor-request at r-project.org" <bioconductor-request at="" r-project.org=""> wrote: Send Bioconductor mailing list submissions to bioconductor at r-project.org To subscribe or unsubscribe via the World Wide Web, visit https://stat.ethz.ch/mailman/listinfo/bioconductor or, via email, send a message with subject or body 'help' to bioconductor-request at r-project.org You can reach the person managing the list at bioconductor-owner at r-project.org When replying, please edit your Subject line so it is more specific than "Re: Contents of Bioconductor digest..." Today's Topics: 1. Heatmaps for Agilent Single Color 4x44k using limma and vsn. (Cornish, Joseph (NIH/NIAID) [F]) 2. Re: Heatmaps for Agilent Single Color 4x44k using limma and vsn. (James W. MacDonald) 3. Re: Heatmaps for Agilent Single Color 4x44k using limma and vsn. (Cornish, Joseph (NIH/NIAID) [F]) 4. Re: Heatmaps for Agilent Single Color 4x44k using limma and vsn. (Sean Davis) 5. Re: Heatmaps for Agilent Single Color 4x44k using limma and vsn. (James W. MacDonald) 6. Re: Heatmaps for Agilent Single Color 4x44k using limma and vsn. (Cornish, Joseph (NIH/NIAID) [F]) 7. analysis of HuGene2.0-st array affymetrix -aroma package (Andreia Fonseca) 8. Re: Heatmaps for Agilent Single Color 4x44k using limma and vsn. (Sean Davis) 9. Interpreting DESeq2 results (Michael Muratet) 10. identifying drosophila miRNA targets (Fiona Ingleby) 11. Re: identifying drosophila miRNA targets (James W. MacDonald) 12. Re: Limit on number of sequence files for forging a BSgenome (Blanchette, Marco) 13. Re: [R] Error in setMethod("combine"... was - Error when installing globaltest package (Shields, Rusty (IMS)) 14. Re: [R] Error in setMethod("combine"... was - Error when installing globaltest package (Martin Morgan) 15. Re: [R] Error in setMethod("combine"... was - Error when installing globaltest package (Shields, Rusty (IMS)) 16. Re: Interpreting DESeq2 results (Michael Love) 17. Re: Interpreting DESeq2 results (somnath bandyopadhyay) ---------------------------------------------------------------------- Message: 1 Date: Thu, 28 Mar 2013 13:27:23 +0000 From: "Cornish, Joseph (NIH/NIAID) [F]" <joseph.cornish@nih.gov> To: "bioconductor at r-project.org" <bioconductor at="" r-project.org=""> Subject: [BioC] Heatmaps for Agilent Single Color 4x44k using limma and vsn. Message-ID: <e68ea43fbf80174c85a92bc24efaa3b84b7bc3 at="" mlbxv01.nih.gov=""> Content-Type: text/plain I am attempting to write a script that will allow me to do bulk benchmarking of various methods for normalization and background correction. I am concerned that I am not processing the data correctly and this may be leading to the abnormal heatmaps that I am seeing. For a test case, I am using an experiment where I have two different conditions with two different states with two replicates per state. The heatmaps concern me because of the uniformity between them, the values are all completely green or completely red, there are differences between genes within a state according to the heatmap. In the examples I have seen people have suggested that I use the coefficients from the eBayes fitting of the data, however in my example the eBayes LMarray object does not have two columns in the results, only on (not sure if this is a sign of errors on my part). I have instead used the lmFit coefficients since there are column for each matrix. To further confuse someone who is new to R, the simplified example and full code generate different heatmaps for the same data and same methods (normexp for background, cyclic loess for norm). The differences being in the simplified version the gene IDs show up on the right side but in the full version, only row number appears. The color for some genes changes between states as well. My questions are: Why does the bayes model only have one column? Am I approaching the processing correctly? Why are there differences between the resulting heatmaps of each version? I know that the idea is to keep the amount of posted code to a minimum but I don't see any other way here. Here is the simplified example I have made: library(limma) library(vsn) library(gplots) set.seed(0xabcd) #from vsn manual #constants for reading agi outs from our imager COL__ <- list(G="Raw intensity (tr.mean) {A}",Gb="Background (tr.mean) {A}") ANO__ <- c("Position", "Name", "ID", "Description", "GeneName") #parse data targets <- readTargets("test.csv", sep = ",") gal <- readGAL("026806_D_20110524.gal") layout <- getLayout(gal) agi <- read.maimages(targets, columns = COL__, annotation = ANO__, green.only = TRUE) design <- model.matrix(~0+factor(c(1,1,2,2))) colnames(design) <- c("s1", "s2") #process bg <- backgroundCorrect(agi, method = "normexp") norm <- normalizeBetweenArrays(bg, method = "cyclicloess") avg <- avereps(norm, ID=norm$genes$ProbeName) #fit contmat <- makeContrasts(s2-s1, levels = design) fit_lm <- lmFit(avg$E, design) fit_be <- eBayes(contrasts.fit(fit_lm, contmat)) #results difexp <- topTable(fit_be, coef = 1, adjust = "BH") res <- decideTests(fit_be) comp <- vennDiagram(res) #plot pdf("example-heatmap.pdf") genes <- as.numeric(rownames(topTable(fit_be, n = 100))) heatmap(fit_lm$coefficients[genes,], col=redgreen(75), key = TRUE) dev.off() The full version: library(vsn) library(limma) library(gplots) set.seed(0xabcd) #from the vsn manual #constants A_SUB__ <- 1 A_NXP__ <- 2 A_MIN__ <- 3 A_VSN__ <- 1 A_QNT__ <- 2 A_LWS__ <- 3 M_SUB__ <- "subtract" M_NXP__ <- "normexp" M_MIN__ <- "minimum" M_VSN__ <- "vsn" M_LWS__ <- "cyclicloess" M_QNT__ <- "quantile" COL__ <- list(G="Raw intensity (tr.mean) {A}",Gb="Background (tr.mean) {A}") ANO__ <- c("Position", "Name", "ID", "Description", "GeneName") N_BGM__ <- 3 #number of background correction methods N_NRM__ <- 3 #number of normalization methods A_BGM__ <- c(A_SUB__, A_NXP__, A_MIN__) A_NRM__ <- c(A_VSN__, A_QNT__, A_LWS__) #utility methods #applies subtraction background correction with limma bg_subtract <- function(x){ return(list(norm = x$norm, bg = M_SUB__, agi = backgroundCorrect(x$agi, method = M_SUB__))) } #applies normexp backgroud correction with limma bg_normexp <- function(x){ return(list(norm = x$norm, bg = M_NXP__, agi = backgroundCorrect(x$agi, method = M_NXP__))) } #applies minimum background correction with limma bg_minimum <- function(x){ return(list(norm = x$norm, bg = M_MIN__, agi = backgroundCorrect(x$agi, method = M_MIN__))) } #applies VSN normalization with VSN package nb_vsn <- function(x){ return(list(norm = M_VSN__, bg = x$bg, agi = normalizeVSN(x$agi))) } #applies VSN normalization with VSN package nb_vsn <- function(x){ return(list(norm = M_VSN__, bg = x$bg, agi = normalizeVSN(x$agi))) } #applies cyclic lowess normalization with limma nb_lowess <- function(x){ return(list(norm = M_LWS__, bg = x$bg, agi = normalizeBetweenArrays(x$agi, method = M_LWS__))) } #applies quantile normalization with limma nb_quantile <- function(x){ return(list(norm = M_QNT__, bg = x$bg, agi = normalizeBetweenArrays(x$agi, method = M_QNT__))) } #generates a plot of MA values with limma pl_ma <- function(x){ name <- paste(c(x[[1]]$bg, x[[1]]$norm, "-ma.pdf"), collapse="") pdf(name) plotMA(x[[1]]$agi$E) dev.off() } #generates a plot of mean versus stdev with VS pl_meansd <- function(x){ name <- paste(c(x[[1]]$bg, x[[1]]$norm, "-meansd.pdf"), collapse="") pdf(name) meanSdPlot(x[[1]]$agi$E) dev.off() } #generate volcanoplot of log-fold versus log-odds pl_volcano <- function(x){ name <- paste(c(x[[1]]$bg, x[[1]]$norm, "-volcano.pdf"), collapse="") pdf(name) volcanoplot(x[[1]]$bayes) dev.off() } #generate plot of expression data over time pl_lines <- function(x){ name <- paste(c(x[[1]]$bg, x[[1]]$norm, "-lines.pdf"), collapse="") pdf(name) plotlines(x[[1]]$agi$E) dev.off() } #generate heatmat of final results pl_heatmap <- function(x, ntop){ name <- paste(c(x[[1]]$bg, x[[1]]$norm, "-heatmap.pdf"), collapse="") pdf(name) genes <- as.numeric(rownames(topTable(x[[1]]$bayes, n = ntop))) heatmap(x[[1]]$lm$coefficients[genes,], col = redgreen(75), key = TRUE) dev.off() } #averages replicates of probes m_repavg <- function(x){ return(list(agi = x[[1]]$agi, avg = avereps(x[[1]]$agi, ID=x[[1]]$agi$genes$ProbeName), bg = x[[1]]$bg, norm = x[[1]]$norm)) } #generate contrasts m_fit <- function(x, design){ contmat <- makeContrasts(s2-s1, levels = design) fit1 <- lmFit(x[[1]]$agi$E, design) fit2 <- eBayes(contrasts.fit(fit1, contmat)) diff <- topTable(fit2, coef = 1, adjust = "BH") res <- decideTests(fit2) comp <- vennDiagram(res) return(list(agi = x[[1]]$agi, bg = x[[1]]$bg, norm = x[[1]]$norm, ave = x[[1]]$ave, cont = contmat, lm = fit1, bayes = fit2, diff = diff, res = res, comp = comp)) } ###################################################################### ######### # Function to parse array data into limma objects # Args: # tfile: target file (see limma documentation) # tsep: limma defaults TSV, have to set manually for CSV and so on # galfile: GAL file for printer/layout # Returns data object with targets, genes and printer layout ###################################################################### ######### parse <- function(tfile, tsep, galfile){ targets <- readTargets(tfile, sep=tsep) gal <- readGAL(galfile) layout <- getLayout(gal) agi <- read.maimages(targets, columns = COL__, annotation = ANO__, green.only = TRUE) bg <- "none" norm <- "none" desg <- model.matrix(~0+factor(c(1,1,2,2))) colnames(desg) <- c("s1","s2") d <- list(agi = agi, gal = gal, layout = layout, targets = targets, bg = bg, desg = desg, norm = norm) processed <- matrix(NA, nrow = N_BGM__, ncol = N_NRM__) processed <- apply(processed, 1:2, function(x){return(d)}) return(processed) } #parse command args args <- commandArgs(trailingOnly = TRUE) #load files into matrix a <- parse(args[1], args[2], args[3]) #keepy a copy of a raw dataset raw <- a[1,1] print(raw[[1]]$desg) #perform backgroung corrections a[,A_SUB__] <- lapply(a[,A_SUB__], bg_subtract) a[,A_NXP__] <- lapply(a[,A_NXP__], bg_normexp) a[,A_MIN__] <- lapply(a[,A_MIN__], bg_minimum) #perform normalization a[A_VSN__,] <- lapply(a[A_VSN__,], nb_vsn) a[A_LWS__,] <- lapply(a[A_LWS__,], nb_lowess) a[A_QNT__,] <- lapply(a[A_QNT__,], nb_quantile) #fit model and differentiate expression a <- apply(a, 1:2, m_repavg) a <- apply(a, 1:2, m_fit, design = raw[[1]]$desg) #generate plots apply(a, 1:2, pl_ma) apply(a, 1:2, pl_meansd) apply(a, 1:2, pl_volcano) apply(a, 1:2, pl_lines) apply(a, 1:2, pl_heatmap, args[4]) [[alternative HTML version deleted]] ------------------------------ Message: 2 Date: Thu, 28 Mar 2013 10:07:32 -0400 From: "James W. MacDonald" <jmacdon@uw.edu> To: "Cornish, Joseph (NIH/NIAID) [F]" <joseph.cornish at="" nih.gov=""> Cc: "bioconductor at r-project.org" <bioconductor at="" r-project.org=""> Subject: Re: [BioC] Heatmaps for Agilent Single Color 4x44k using limma and vsn. Message-ID: <51544EA4.5080805 at uw.edu> Content-Type: text/plain; charset=ISO-8859-1; format=flowed Hi Joseph, On 3/28/2013 9:27 AM, Cornish, Joseph (NIH/NIAID) [F] wrote: > I am attempting to write a script that will allow me to do bulk benchmarking of various methods for normalization and background correction. I am concerned that I am not processing the data correctly and this may be leading to the abnormal heatmaps that I am seeing. For a test case, I am using an experiment where I have two different conditions with two different states with two replicates per state. > > The heatmaps concern me because of the uniformity between them, the values are all completely green or completely red, there are differences between genes within a state according to the heatmap. > > In the examples I have seen people have suggested that I use the coefficients from the eBayes fitting of the data, however in my example the eBayes LMarray object does not have two columns in the results, only on (not sure if this is a sign of errors on my part). I have instead used the lmFit coefficients since there are column for each matrix. > > To further confuse someone who is new to R, the simplified example and full code generate different heatmaps for the same data and same methods (normexp for background, cyclic loess for norm). The differences being in the simplified version the gene IDs show up on the right side but in the full version, only row number appears. The color for some genes changes between states as well. > > My questions are: > Why does the bayes model only have one column? Because you are computing one contrast (which is all you can do). One contrast means one test which means one column of results. In other words, the coefficient you are estimating in the contrasts.fit() step is the difference in means between the two groups, which is one value per gene. > Am I approaching the processing correctly? I tried to read through that business below, but got some weird affliction due to an apparent overdose of underscores and had to stop. > Why are there differences between the resulting heatmaps of each version? I don't know, as I worried my brain might explode if I read any more. But you should note that it is fairly unusual to plot coefficients in a heatmap, as what you are then presenting are log ratios, and most people can't understand that. Instead you might consider getting the list of genes, and then creating a heatmap using the normalized expression values for each sample instead. You then would want to scale in some way. We often sweep the mean of the controls out of the matrix (by row), so the colors can then be interpreted as differences between the treated samples and controls. Best, Jim > > I know that the idea is to keep the amount of posted code to a minimum but I don't see any other way here. > > Here is the simplified example I have made: > library(limma) > library(vsn) > library(gplots) > > set.seed(0xabcd) #from vsn manual > > #constants for reading agi outs from our imager > COL__<- list(G="Raw intensity (tr.mean) {A}",Gb="Background (tr.mean) {A}") > ANO__<- c("Position", "Name", "ID", "Description", "GeneName") > > #parse data > targets<- readTargets("test.csv", sep = ",") > gal<- readGAL("026806_D_20110524.gal") > layout<- getLayout(gal) > agi<- read.maimages(targets, columns = COL__, > annotation = ANO__, > green.only = TRUE) > design<- model.matrix(~0+factor(c(1,1,2,2))) > colnames(design)<- c("s1", "s2") > > #process > bg<- backgroundCorrect(agi, method = "normexp") > norm<- normalizeBetweenArrays(bg, method = "cyclicloess") > avg<- avereps(norm, ID=norm$genes$ProbeName) > > #fit > contmat<- makeContrasts(s2-s1, levels = design) > fit_lm<- lmFit(avg$E, design) > fit_be<- eBayes(contrasts.fit(fit_lm, contmat)) > > #results > difexp<- topTable(fit_be, coef = 1, adjust = "BH") > res<- decideTests(fit_be) > comp<- vennDiagram(res) > > #plot > pdf("example-heatmap.pdf") > genes<- as.numeric(rownames(topTable(fit_be, n = 100))) > heatmap(fit_lm$coefficients[genes,], col=redgreen(75), key = TRUE) > dev.off() > > > The full version: > library(vsn) > library(limma) > library(gplots) > set.seed(0xabcd) #from the vsn manual > > #constants > A_SUB__<- 1 > A_NXP__<- 2 > A_MIN__<- 3 > A_VSN__<- 1 > A_QNT__<- 2 > A_LWS__<- 3 > M_SUB__<- "subtract" > M_NXP__<- "normexp" > M_MIN__<- "minimum" > M_VSN__<- "vsn" > M_LWS__<- "cyclicloess" > M_QNT__<- "quantile" > COL__<- list(G="Raw intensity (tr.mean) {A}",Gb="Background (tr.mean) {A}") > ANO__<- c("Position", "Name", "ID", "Description", "GeneName") > N_BGM__<- 3 #number of background correction methods > N_NRM__<- 3 #number of normalization methods > A_BGM__<- c(A_SUB__, A_NXP__, A_MIN__) > A_NRM__<- c(A_VSN__, A_QNT__, A_LWS__) > > > #utility methods > #applies subtraction background correction with limma > bg_subtract<- function(x){ > return(list(norm = x$norm, > bg = M_SUB__, > agi = backgroundCorrect(x$agi, method = M_SUB__))) > } > > #applies normexp backgroud correction with limma > bg_normexp<- function(x){ > return(list(norm = x$norm, > bg = M_NXP__, > agi = backgroundCorrect(x$agi, method = M_NXP__))) > } > > #applies minimum background correction with limma > bg_minimum<- function(x){ > return(list(norm = x$norm, > bg = M_MIN__, > agi = backgroundCorrect(x$agi, method = M_MIN__))) > } > > #applies VSN normalization with VSN package > nb_vsn<- function(x){ > return(list(norm = M_VSN__, > bg = x$bg, > agi = normalizeVSN(x$agi))) > } > > #applies VSN normalization with VSN package > nb_vsn<- function(x){ > return(list(norm = M_VSN__, > bg = x$bg, > agi = normalizeVSN(x$agi))) > } > > #applies cyclic lowess normalization with limma > nb_lowess<- function(x){ > return(list(norm = M_LWS__, > bg = x$bg, > agi = normalizeBetweenArrays(x$agi, method = M_LWS__))) > } > > #applies quantile normalization with limma > nb_quantile<- function(x){ > return(list(norm = M_QNT__, > bg = x$bg, > agi = normalizeBetweenArrays(x$agi, method = M_QNT__))) > } > > #generates a plot of MA values with limma > pl_ma<- function(x){ > name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-ma.pdf"), collapse="") > pdf(name) > plotMA(x[[1]]$agi$E) > dev.off() > } > > #generates a plot of mean versus stdev with VS > pl_meansd<- function(x){ > name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-meansd.pdf"), collapse="") > pdf(name) > meanSdPlot(x[[1]]$agi$E) > dev.off() > } > > #generate volcanoplot of log-fold versus log-odds > pl_volcano<- function(x){ > name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-volcano.pdf"), collapse="") > pdf(name) > volcanoplot(x[[1]]$bayes) > dev.off() > } > > #generate plot of expression data over time > pl_lines<- function(x){ > name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-lines.pdf"), collapse="") > pdf(name) > plotlines(x[[1]]$agi$E) > dev.off() > } > > #generate heatmat of final results > pl_heatmap<- function(x, ntop){ > name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-heatmap.pdf"), collapse="") > pdf(name) > genes<- as.numeric(rownames(topTable(x[[1]]$bayes, n = ntop))) > heatmap(x[[1]]$lm$coefficients[genes,], > col = redgreen(75), > key = TRUE) > dev.off() > } > > #averages replicates of probes > m_repavg<- function(x){ > return(list(agi = x[[1]]$agi, > avg = avereps(x[[1]]$agi, ID=x[[1]]$agi$genes$ProbeName), > bg = x[[1]]$bg, > norm = x[[1]]$norm)) > } > > #generate contrasts > m_fit<- function(x, design){ > contmat<- makeContrasts(s2-s1, > levels = design) > fit1<- lmFit(x[[1]]$agi$E, design) > fit2<- eBayes(contrasts.fit(fit1, contmat)) > diff<- topTable(fit2, coef = 1, adjust = "BH") > res<- decideTests(fit2) > comp<- vennDiagram(res) > return(list(agi = x[[1]]$agi, > bg = x[[1]]$bg, > norm = x[[1]]$norm, > ave = x[[1]]$ave, > cont = contmat, > lm = fit1, > bayes = fit2, > diff = diff, > res = res, > comp = comp)) > } > > #################################################################### ########### > # Function to parse array data into limma objects > # Args: > # tfile: target file (see limma documentation) > # tsep: limma defaults TSV, have to set manually for CSV and so on > # galfile: GAL file for printer/layout > # Returns data object with targets, genes and printer layout > #################################################################### ########### > parse<- function(tfile, tsep, galfile){ > targets<- readTargets(tfile, sep=tsep) > gal<- readGAL(galfile) > layout<- getLayout(gal) > agi<- read.maimages(targets, columns = COL__, > annotation = ANO__, > green.only = TRUE) > bg<- "none" > norm<- "none" > desg<- model.matrix(~0+factor(c(1,1,2,2))) > colnames(desg)<- c("s1","s2") > d<- list(agi = agi, > gal = gal, > layout = layout, > targets = targets, > bg = bg, > desg = desg, > norm = norm) > > processed<- matrix(NA, nrow = N_BGM__, ncol = N_NRM__) > processed<- apply(processed, 1:2, function(x){return(d)}) > return(processed) > } > > #parse command args > args<- commandArgs(trailingOnly = TRUE) > > #load files into matrix > a<- parse(args[1], args[2], args[3]) > > #keepy a copy of a raw dataset > raw<- a[1,1] > print(raw[[1]]$desg) > > #perform backgroung corrections > a[,A_SUB__]<- lapply(a[,A_SUB__], bg_subtract) > a[,A_NXP__]<- lapply(a[,A_NXP__], bg_normexp) > a[,A_MIN__]<- lapply(a[,A_MIN__], bg_minimum) > > #perform normalization > a[A_VSN__,]<- lapply(a[A_VSN__,], nb_vsn) > a[A_LWS__,]<- lapply(a[A_LWS__,], nb_lowess) > a[A_QNT__,]<- lapply(a[A_QNT__,], nb_quantile) > > #fit model and differentiate expression > a<- apply(a, 1:2, m_repavg) > a<- apply(a, 1:2, m_fit, design = raw[[1]]$desg) > > #generate plots > apply(a, 1:2, pl_ma) > apply(a, 1:2, pl_meansd) > apply(a, 1:2, pl_volcano) > apply(a, 1:2, pl_lines) > apply(a, 1:2, pl_heatmap, args[4]) > > > [[alternative HTML version deleted]] > > _______________________________________________ > Bioconductor mailing list > Bioconductor at r-project.org > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor -- James W. MacDonald, M.S. Biostatistician University of Washington Environmental and Occupational Health Sciences 4225 Roosevelt Way NE, # 100 Seattle WA 98105-6099 ------------------------------ Message: 3 Date: Thu, 28 Mar 2013 14:22:32 +0000 From: "Cornish, Joseph (NIH/NIAID) [F]" <joseph.cornish@nih.gov> To: "James W. MacDonald" <jmacdon at="" uw.edu=""> Cc: "bioconductor at r-project.org" <bioconductor at="" r-project.org=""> Subject: Re: [BioC] Heatmaps for Agilent Single Color 4x44k using limma and vsn. Message-ID: <e68ea43fbf80174c85a92bc24efaa3b84b7bdf at="" mlbxv01.nih.gov=""> Content-Type: text/plain; charset="us-ascii" Hello James, Are you referring to function names or constants? The constants are an old habit I picked up long ago in my CS classes, there are a ton but it makes them stand out for me in vim more readily. I've come to avoid camel case since switching to linux, it makes searching for things a nightmare, it can also help if I group functions with a prefix (eg plot functions start with "pl_"). If you want I can find+replace all of the constants with their values... The differences came from using the unaveraged replicates in the larger script versus the average in the example. I'm not sure what you mean by "mean of the controls", would this simply be filtering out the control spots? I'll generate the heatmaps using the normalized expression values. I suppose in that case it would be worthwhile to filter ahead of time by p-value and fold change (the coef in topTable). _________thanks_____________ _______Joe______ -----Original Message----- From: James W. MacDonald [mailto:jmacdon@uw.edu] Sent: Thursday, March 28, 2013 10:08 AM To: Cornish, Joseph (NIH/NIAID) [F] Cc: bioconductor at r-project.org Subject: Re: [BioC] Heatmaps for Agilent Single Color 4x44k using limma and vsn. Hi Joseph, On 3/28/2013 9:27 AM, Cornish, Joseph (NIH/NIAID) [F] wrote: > I am attempting to write a script that will allow me to do bulk benchmarking of various methods for normalization and background correction. I am concerned that I am not processing the data correctly and this may be leading to the abnormal heatmaps that I am seeing. For a test case, I am using an experiment where I have two different conditions with two different states with two replicates per state. > > The heatmaps concern me because of the uniformity between them, the values are all completely green or completely red, there are differences between genes within a state according to the heatmap. > > In the examples I have seen people have suggested that I use the coefficients from the eBayes fitting of the data, however in my example the eBayes LMarray object does not have two columns in the results, only on (not sure if this is a sign of errors on my part). I have instead used the lmFit coefficients since there are column for each matrix. > > To further confuse someone who is new to R, the simplified example and full code generate different heatmaps for the same data and same methods (normexp for background, cyclic loess for norm). The differences being in the simplified version the gene IDs show up on the right side but in the full version, only row number appears. The color for some genes changes between states as well. > > My questions are: > Why does the bayes model only have one column? Because you are computing one contrast (which is all you can do). One contrast means one test which means one column of results. In other words, the coefficient you are estimating in the contrasts.fit() step is the difference in means between the two groups, which is one value per gene. > Am I approaching the processing correctly? I tried to read through that business below, but got some weird affliction due to an apparent overdose of underscores and had to stop. > Why are there differences between the resulting heatmaps of each version? I don't know, as I worried my brain might explode if I read any more. But you should note that it is fairly unusual to plot coefficients in a heatmap, as what you are then presenting are log ratios, and most people can't understand that. Instead you might consider getting the list of genes, and then creating a heatmap using the normalized expression values for each sample instead. You then would want to scale in some way. We often sweep the mean of the controls out of the matrix (by row), so the colors can then be interpreted as differences between the treated samples and controls. Best, Jim > > I know that the idea is to keep the amount of posted code to a minimum but I don't see any other way here. > > Here is the simplified example I have made: > library(limma) > library(vsn) > library(gplots) > > set.seed(0xabcd) #from vsn manual > > #constants for reading agi outs from our imager > COL__<- list(G="Raw intensity (tr.mean) {A}",Gb="Background (tr.mean) > {A}") > ANO__<- c("Position", "Name", "ID", "Description", "GeneName") > > #parse data > targets<- readTargets("test.csv", sep = ",") > gal<- readGAL("026806_D_20110524.gal") > layout<- getLayout(gal) > agi<- read.maimages(targets, columns = COL__, > annotation = ANO__, > green.only = TRUE) > design<- model.matrix(~0+factor(c(1,1,2,2))) > colnames(design)<- c("s1", "s2") > > #process > bg<- backgroundCorrect(agi, method = "normexp") > norm<- normalizeBetweenArrays(bg, method = "cyclicloess") > avg<- avereps(norm, ID=norm$genes$ProbeName) > > #fit > contmat<- makeContrasts(s2-s1, levels = design) > fit_lm<- lmFit(avg$E, design) > fit_be<- eBayes(contrasts.fit(fit_lm, contmat)) > > #results > difexp<- topTable(fit_be, coef = 1, adjust = "BH") > res<- decideTests(fit_be) > comp<- vennDiagram(res) > > #plot > pdf("example-heatmap.pdf") > genes<- as.numeric(rownames(topTable(fit_be, n = 100))) > heatmap(fit_lm$coefficients[genes,], col=redgreen(75), key = TRUE) > dev.off() > > > The full version: > library(vsn) > library(limma) > library(gplots) > set.seed(0xabcd) #from the vsn manual > > #constants > A_SUB__<- 1 > A_NXP__<- 2 > A_MIN__<- 3 > A_VSN__<- 1 > A_QNT__<- 2 > A_LWS__<- 3 > M_SUB__<- "subtract" > M_NXP__<- "normexp" > M_MIN__<- "minimum" > M_VSN__<- "vsn" > M_LWS__<- "cyclicloess" > M_QNT__<- "quantile" > COL__<- list(G="Raw intensity (tr.mean) {A}",Gb="Background (tr.mean) > {A}") > ANO__<- c("Position", "Name", "ID", "Description", "GeneName") > N_BGM__<- 3 #number of background correction methods > N_NRM__<- 3 #number of normalization methods > A_BGM__<- c(A_SUB__, A_NXP__, A_MIN__) > A_NRM__<- c(A_VSN__, A_QNT__, A_LWS__) > > > #utility methods > #applies subtraction background correction with limma > bg_subtract<- function(x){ > return(list(norm = x$norm, > bg = M_SUB__, > agi = backgroundCorrect(x$agi, method = M_SUB__))) } > > #applies normexp backgroud correction with limma > bg_normexp<- function(x){ > return(list(norm = x$norm, > bg = M_NXP__, > agi = backgroundCorrect(x$agi, method = M_NXP__))) } > > #applies minimum background correction with limma > bg_minimum<- function(x){ > return(list(norm = x$norm, > bg = M_MIN__, > agi = backgroundCorrect(x$agi, method = M_MIN__))) } > > #applies VSN normalization with VSN package > nb_vsn<- function(x){ > return(list(norm = M_VSN__, > bg = x$bg, > agi = normalizeVSN(x$agi))) } > > #applies VSN normalization with VSN package > nb_vsn<- function(x){ > return(list(norm = M_VSN__, > bg = x$bg, > agi = normalizeVSN(x$agi))) } > > #applies cyclic lowess normalization with limma > nb_lowess<- function(x){ > return(list(norm = M_LWS__, > bg = x$bg, > agi = normalizeBetweenArrays(x$agi, method = > M_LWS__))) } > > #applies quantile normalization with limma > nb_quantile<- function(x){ > return(list(norm = M_QNT__, > bg = x$bg, > agi = normalizeBetweenArrays(x$agi, method = > M_QNT__))) } > > #generates a plot of MA values with limma > pl_ma<- function(x){ > name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-ma.pdf"), collapse="") > pdf(name) > plotMA(x[[1]]$agi$E) > dev.off() > } > > #generates a plot of mean versus stdev with VS > pl_meansd<- function(x){ > name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-meansd.pdf"), collapse="") > pdf(name) > meanSdPlot(x[[1]]$agi$E) > dev.off() > } > > #generate volcanoplot of log-fold versus log-odds > pl_volcano<- function(x){ > name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-volcano.pdf"), collapse="") > pdf(name) > volcanoplot(x[[1]]$bayes) > dev.off() > } > > #generate plot of expression data over time > pl_lines<- function(x){ > name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-lines.pdf"), collapse="") > pdf(name) > plotlines(x[[1]]$agi$E) > dev.off() > } > > #generate heatmat of final results > pl_heatmap<- function(x, ntop){ > name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-heatmap.pdf"), collapse="") > pdf(name) > genes<- as.numeric(rownames(topTable(x[[1]]$bayes, n = ntop))) > heatmap(x[[1]]$lm$coefficients[genes,], > col = redgreen(75), > key = TRUE) > dev.off() > } > > #averages replicates of probes > m_repavg<- function(x){ > return(list(agi = x[[1]]$agi, > avg = avereps(x[[1]]$agi, ID=x[[1]]$agi$genes$ProbeName), > bg = x[[1]]$bg, > norm = x[[1]]$norm)) > } > > #generate contrasts > m_fit<- function(x, design){ > contmat<- makeContrasts(s2-s1, > levels = design) > fit1<- lmFit(x[[1]]$agi$E, design) > fit2<- eBayes(contrasts.fit(fit1, contmat)) > diff<- topTable(fit2, coef = 1, adjust = "BH") > res<- decideTests(fit2) > comp<- vennDiagram(res) > return(list(agi = x[[1]]$agi, > bg = x[[1]]$bg, > norm = x[[1]]$norm, > ave = x[[1]]$ave, > cont = contmat, > lm = fit1, > bayes = fit2, > diff = diff, > res = res, > comp = comp)) > } > > ###################################################################### > ######### # Function to parse array data into limma objects # Args: > # tfile: target file (see limma documentation) > # tsep: limma defaults TSV, have to set manually for CSV and so on > # galfile: GAL file for printer/layout > # Returns data object with targets, genes and printer layout > ###################################################################### > ######### > parse<- function(tfile, tsep, galfile){ > targets<- readTargets(tfile, sep=tsep) > gal<- readGAL(galfile) > layout<- getLayout(gal) > agi<- read.maimages(targets, columns = COL__, > annotation = ANO__, > green.only = TRUE) > bg<- "none" > norm<- "none" > desg<- model.matrix(~0+factor(c(1,1,2,2))) > colnames(desg)<- c("s1","s2") > d<- list(agi = agi, > gal = gal, > layout = layout, > targets = targets, > bg = bg, > desg = desg, > norm = norm) > > processed<- matrix(NA, nrow = N_BGM__, ncol = N_NRM__) > processed<- apply(processed, 1:2, function(x){return(d)}) > return(processed) > } > > #parse command args > args<- commandArgs(trailingOnly = TRUE) > > #load files into matrix > a<- parse(args[1], args[2], args[3]) > > #keepy a copy of a raw dataset > raw<- a[1,1] > print(raw[[1]]$desg) > > #perform backgroung corrections > a[,A_SUB__]<- lapply(a[,A_SUB__], bg_subtract) > a[,A_NXP__]<- lapply(a[,A_NXP__], bg_normexp) > a[,A_MIN__]<- lapply(a[,A_MIN__], bg_minimum) > > #perform normalization > a[A_VSN__,]<- lapply(a[A_VSN__,], nb_vsn) > a[A_LWS__,]<- lapply(a[A_LWS__,], nb_lowess) > a[A_QNT__,]<- lapply(a[A_QNT__,], nb_quantile) > > #fit model and differentiate expression > a<- apply(a, 1:2, m_repavg) > a<- apply(a, 1:2, m_fit, design = raw[[1]]$desg) > > #generate plots > apply(a, 1:2, pl_ma) > apply(a, 1:2, pl_meansd) > apply(a, 1:2, pl_volcano) > apply(a, 1:2, pl_lines) > apply(a, 1:2, pl_heatmap, args[4]) > > > [[alternative HTML version deleted]] > > _______________________________________________ > Bioconductor mailing list > Bioconductor at r-project.org > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: > http://news.gmane.org/gmane.science.biology.informatics.conductor -- James W. MacDonald, M.S. Biostatistician University of Washington Environmental and Occupational Health Sciences 4225 Roosevelt Way NE, # 100 Seattle WA 98105-6099 ------------------------------ Message: 4 Date: Thu, 28 Mar 2013 10:33:51 -0400 From: Sean Davis <sdavis2@mail.nih.gov> To: "Cornish, Joseph (NIH/NIAID) [F]" <joseph.cornish at="" nih.gov=""> Cc: "bioconductor at r-project.org" <bioconductor at="" r-project.org="">, "James W. MacDonald" <jmacdon at="" uw.edu=""> Subject: Re: [BioC] Heatmaps for Agilent Single Color 4x44k using limma and vsn. Message-ID: <caneavbkyv+mb61jvne0aklicxovxyvelusbhgowwr61_q9s-dw at="" mail.gmail.com=""> Content-Type: text/plain On Thu, Mar 28, 2013 at 10:22 AM, Cornish, Joseph (NIH/NIAID) [F] < joseph.cornish at nih.gov> wrote: > Hello James, > > Are you referring to function names or constants? The constants are an old > habit I picked up long ago in my CS classes, there are a ton but it makes > them stand out for me in vim more readily. I've come to avoid camel case > since switching to linux, it makes searching for things a nightmare, it can > also help if I group functions with a prefix (eg plot functions start with > "pl_"). If you want I can find+replace all of the constants with their > values... > > Hi, Joe. I think I could help to translate part of what James said. If you could boil your example down to a few lines of reproducible code that show the problem, that would be simpler to interpret and easier to comment on directly. Sometimes in the process of producing such a simplified example, I have found my own misunderstandings. It also sometimes helps to get someone local to give a second pair of eyes. > The differences came from using the unaveraged replicates in the larger > script versus the average in the example. > > I'm not sure what you mean by "mean of the controls", would this simply be > filtering out the control spots? > > No. I think he means for you to subtract the mean of the control samples (or one group of samples) for each gene from all the normalized expression values. > I'll generate the heatmaps using the normalized expression values. I > suppose in that case it would be worthwhile to filter ahead of time by > p-value and fold change (the coef in topTable). > > Yes. Sean > > _________thanks_____________ > _______Joe______ > > -----Original Message----- > From: James W. MacDonald [mailto:jmacdon at uw.edu] > Sent: Thursday, March 28, 2013 10:08 AM > To: Cornish, Joseph (NIH/NIAID) [F] > Cc: bioconductor at r-project.org > Subject: Re: [BioC] Heatmaps for Agilent Single Color 4x44k using limma > and vsn. > > Hi Joseph, > > > On 3/28/2013 9:27 AM, Cornish, Joseph (NIH/NIAID) [F] wrote: > > I am attempting to write a script that will allow me to do bulk > benchmarking of various methods for normalization and background > correction. I am concerned that I am not processing the data correctly and > this may be leading to the abnormal heatmaps that I am seeing. For a test > case, I am using an experiment where I have two different conditions with > two different states with two replicates per state. > > > > The heatmaps concern me because of the uniformity between them, the > values are all completely green or completely red, there are differences > between genes within a state according to the heatmap. > > > > In the examples I have seen people have suggested that I use the > coefficients from the eBayes fitting of the data, however in my example the > eBayes LMarray object does not have two columns in the results, only on > (not sure if this is a sign of errors on my part). I have instead used the > lmFit coefficients since there are column for each matrix. > > > > To further confuse someone who is new to R, the simplified example and > full code generate different heatmaps for the same data and same methods > (normexp for background, cyclic loess for norm). The differences being in > the simplified version the gene IDs show up on the right side but in the > full version, only row number appears. The color for some genes changes > between states as well. > > > > My questions are: > > Why does the bayes model only have one column? > > Because you are computing one contrast (which is all you can do). One > contrast means one test which means one column of results. In other words, > the coefficient you are estimating in the contrasts.fit() step is the > difference in means between the two groups, which is one value per gene. > > > Am I approaching the processing correctly? > > I tried to read through that business below, but got some weird affliction > due to an apparent overdose of underscores and had to stop. > > > Why are there differences between the resulting heatmaps of each version? > > I don't know, as I worried my brain might explode if I read any more. > But you should note that it is fairly unusual to plot coefficients in a > heatmap, as what you are then presenting are log ratios, and most people > can't understand that. > > Instead you might consider getting the list of genes, and then creating a > heatmap using the normalized expression values for each sample instead. You > then would want to scale in some way. We often sweep the mean of the > controls out of the matrix (by row), so the colors can then be interpreted > as differences between the treated samples and controls. > > Best, > > Jim > > > > > > I know that the idea is to keep the amount of posted code to a minimum > but I don't see any other way here. > > > > Here is the simplified example I have made: > > library(limma) > > library(vsn) > > library(gplots) > > > > set.seed(0xabcd) #from vsn manual > > > > #constants for reading agi outs from our imager > > COL__<- list(G="Raw intensity (tr.mean) {A}",Gb="Background (tr.mean) > > {A}") > > ANO__<- c("Position", "Name", "ID", "Description", "GeneName") > > > > #parse data > > targets<- readTargets("test.csv", sep = ",") > > gal<- readGAL("026806_D_20110524.gal") > > layout<- getLayout(gal) > > agi<- read.maimages(targets, columns = COL__, > > annotation = ANO__, > > green.only = TRUE) > > design<- model.matrix(~0+factor(c(1,1,2,2))) > > colnames(design)<- c("s1", "s2") > > > > #process > > bg<- backgroundCorrect(agi, method = "normexp") > > norm<- normalizeBetweenArrays(bg, method = "cyclicloess") > > avg<- avereps(norm, ID=norm$genes$ProbeName) > > > > #fit > > contmat<- makeContrasts(s2-s1, levels = design) > > fit_lm<- lmFit(avg$E, design) > > fit_be<- eBayes(contrasts.fit(fit_lm, contmat)) > > > > #results > > difexp<- topTable(fit_be, coef = 1, adjust = "BH") > > res<- decideTests(fit_be) > > comp<- vennDiagram(res) > > > > #plot > > pdf("example-heatmap.pdf") > > genes<- as.numeric(rownames(topTable(fit_be, n = 100))) > > heatmap(fit_lm$coefficients[genes,], col=redgreen(75), key = TRUE) > > dev.off() > > > > > > The full version: > > library(vsn) > > library(limma) > > library(gplots) > > set.seed(0xabcd) #from the vsn manual > > > > #constants > > A_SUB__<- 1 > > A_NXP__<- 2 > > A_MIN__<- 3 > > A_VSN__<- 1 > > A_QNT__<- 2 > > A_LWS__<- 3 > > M_SUB__<- "subtract" > > M_NXP__<- "normexp" > > M_MIN__<- "minimum" > > M_VSN__<- "vsn" > > M_LWS__<- "cyclicloess" > > M_QNT__<- "quantile" > > COL__<- list(G="Raw intensity (tr.mean) {A}",Gb="Background (tr.mean) > > {A}") > > ANO__<- c("Position", "Name", "ID", "Description", "GeneName") > > N_BGM__<- 3 #number of background correction methods > > N_NRM__<- 3 #number of normalization methods > > A_BGM__<- c(A_SUB__, A_NXP__, A_MIN__) > > A_NRM__<- c(A_VSN__, A_QNT__, A_LWS__) > > > > > > #utility methods > > #applies subtraction background correction with limma > > bg_subtract<- function(x){ > > return(list(norm = x$norm, > > bg = M_SUB__, > > agi = backgroundCorrect(x$agi, method = M_SUB__))) } > > > > #applies normexp backgroud correction with limma > > bg_normexp<- function(x){ > > return(list(norm = x$norm, > > bg = M_NXP__, > > agi = backgroundCorrect(x$agi, method = M_NXP__))) } > > > > #applies minimum background correction with limma > > bg_minimum<- function(x){ > > return(list(norm = x$norm, > > bg = M_MIN__, > > agi = backgroundCorrect(x$agi, method = M_MIN__))) } > > > > #applies VSN normalization with VSN package > > nb_vsn<- function(x){ > > return(list(norm = M_VSN__, > > bg = x$bg, > > agi = normalizeVSN(x$agi))) } > > > > #applies VSN normalization with VSN package > > nb_vsn<- function(x){ > > return(list(norm = M_VSN__, > > bg = x$bg, > > agi = normalizeVSN(x$agi))) } > > > > #applies cyclic lowess normalization with limma > > nb_lowess<- function(x){ > > return(list(norm = M_LWS__, > > bg = x$bg, > > agi = normalizeBetweenArrays(x$agi, method = > > M_LWS__))) } > > > > #applies quantile normalization with limma > > nb_quantile<- function(x){ > > return(list(norm = M_QNT__, > > bg = x$bg, > > agi = normalizeBetweenArrays(x$agi, method = > > M_QNT__))) } > > > > #generates a plot of MA values with limma > > pl_ma<- function(x){ > > name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-ma.pdf"), collapse="") > > pdf(name) > > plotMA(x[[1]]$agi$E) > > dev.off() > > } > > > > #generates a plot of mean versus stdev with VS > > pl_meansd<- function(x){ > > name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-meansd.pdf"), collapse="") > > pdf(name) > > meanSdPlot(x[[1]]$agi$E) > > dev.off() > > } > > > > #generate volcanoplot of log-fold versus log-odds > > pl_volcano<- function(x){ > > name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-volcano.pdf"), collapse="") > > pdf(name) > > volcanoplot(x[[1]]$bayes) > > dev.off() > > } > > > > #generate plot of expression data over time > > pl_lines<- function(x){ > > name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-lines.pdf"), collapse="") > > pdf(name) > > plotlines(x[[1]]$agi$E) > > dev.off() > > } > > > > #generate heatmat of final results > > pl_heatmap<- function(x, ntop){ > > name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-heatmap.pdf"), collapse="") > > pdf(name) > > genes<- as.numeric(rownames(topTable(x[[1]]$bayes, n = ntop))) > > heatmap(x[[1]]$lm$coefficients[genes,], > > col = redgreen(75), > > key = TRUE) > > dev.off() > > } > > > > #averages replicates of probes > > m_repavg<- function(x){ > > return(list(agi = x[[1]]$agi, > > avg = avereps(x[[1]]$agi, > ID=x[[1]]$agi$genes$ProbeName), > > bg = x[[1]]$bg, > > norm = x[[1]]$norm)) > > } > > > > #generate contrasts > > m_fit<- function(x, design){ > > contmat<- makeContrasts(s2-s1, > > levels = design) > > fit1<- lmFit(x[[1]]$agi$E, design) > > fit2<- eBayes(contrasts.fit(fit1, contmat)) > > diff<- topTable(fit2, coef = 1, adjust = "BH") > > res<- decideTests(fit2) > > comp<- vennDiagram(res) > > return(list(agi = x[[1]]$agi, > > bg = x[[1]]$bg, > > norm = x[[1]]$norm, > > ave = x[[1]]$ave, > > cont = contmat, > > lm = fit1, > > bayes = fit2, > > diff = diff, > > res = res, > > comp = comp)) > > } > > > > ###################################################################### > > ######### # Function to parse array data into limma objects # Args: > > # tfile: target file (see limma documentation) > > # tsep: limma defaults TSV, have to set manually for CSV and so on > > # galfile: GAL file for printer/layout > > # Returns data object with targets, genes and printer layout > > ###################################################################### > > ######### > > parse<- function(tfile, tsep, galfile){ > > targets<- readTargets(tfile, sep=tsep) > > gal<- readGAL(galfile) > > layout<- getLayout(gal) > > agi<- read.maimages(targets, columns = COL__, > > annotation = ANO__, > > green.only = TRUE) > > bg<- "none" > > norm<- "none" > > desg<- model.matrix(~0+factor(c(1,1,2,2))) > > colnames(desg)<- c("s1","s2") > > d<- list(agi = agi, > > gal = gal, > > layout = layout, > > targets = targets, > > bg = bg, > > desg = desg, > > norm = norm) > > > > processed<- matrix(NA, nrow = N_BGM__, ncol = N_NRM__) > > processed<- apply(processed, 1:2, function(x){return(d)}) > > return(processed) > > } > > > > #parse command args > > args<- commandArgs(trailingOnly = TRUE) > > > > #load files into matrix > > a<- parse(args[1], args[2], args[3]) > > > > #keepy a copy of a raw dataset > > raw<- a[1,1] > > print(raw[[1]]$desg) > > > > #perform backgroung corrections > > a[,A_SUB__]<- lapply(a[,A_SUB__], bg_subtract) > > a[,A_NXP__]<- lapply(a[,A_NXP__], bg_normexp) > > a[,A_MIN__]<- lapply(a[,A_MIN__], bg_minimum) > > > > #perform normalization > > a[A_VSN__,]<- lapply(a[A_VSN__,], nb_vsn) > > a[A_LWS__,]<- lapply(a[A_LWS__,], nb_lowess) > > a[A_QNT__,]<- lapply(a[A_QNT__,], nb_quantile) > > > > #fit model and differentiate expression > > a<- apply(a, 1:2, m_repavg) > > a<- apply(a, 1:2, m_fit, design = raw[[1]]$desg) > > > > #generate plots > > apply(a, 1:2, pl_ma) > > apply(a, 1:2, pl_meansd) > > apply(a, 1:2, pl_volcano) > > apply(a, 1:2, pl_lines) > > apply(a, 1:2, pl_heatmap, args[4]) > > > > > > [[alternative HTML version deleted]] > > > > _______________________________________________ > > Bioconductor mailing list > > Bioconductor at r-project.org > > https://stat.ethz.ch/mailman/listinfo/bioconductor > > Search the archives: > > http://news.gmane.org/gmane.science.biology.informatics.conductor > > -- > James W. MacDonald, M.S. > Biostatistician > University of Washington > Environmental and Occupational Health Sciences > 4225 Roosevelt Way NE, # 100 > Seattle WA 98105-6099 > > _______________________________________________ > Bioconductor mailing list > Bioconductor at r-project.org > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: > http://news.gmane.org/gmane.science.biology.informatics.conductor > [[alternative HTML version deleted]] ------------------------------ Message: 5 Date: Thu, 28 Mar 2013 10:45:21 -0400 From: "James W. MacDonald" <jmacdon@uw.edu> To: "Cornish, Joseph (NIH/NIAID) [F]" <joseph.cornish at="" nih.gov=""> Cc: "bioconductor at r-project.org" <bioconductor at="" r-project.org=""> Subject: Re: [BioC] Heatmaps for Agilent Single Color 4x44k using limma and vsn. Message-ID: <51545781.7000100 at uw.edu> Content-Type: text/plain; charset=ISO-8859-1; format=flowed Hi Joe, I don't mean to be critical. If you like to use gobtons of underscores, then rock on. I would just find it very distracting. Anyway, by 'mean of the controls', I mean the average value of the control samples. But this could be generalized in any way. For example, your contrast was s2 - s1, so you could do mat <- avg$E[<index of="" interesting="" genes="">,] mn <- fit_lm[<index of="" interesting="" genes="">,"s1"] swept <- sweep(mat, 1, mn, "-") and then do the heatmap. In this case the colors would be interpreted as the log fold difference between any sample and the mean of the s1 samples. But the general take home message here is that you can't just use the raw expression values, as you will then have a solid color to your heatmap. Other possibilities are to use the scale argument to heatmap.2() to scale by row (convert to z-scores), in which case the interpretation of the colors is less obvious. Best, Jim On 3/28/2013 10:22 AM, Cornish, Joseph (NIH/NIAID) [F] wrote: > Hello James, > > Are you referring to function names or constants? The constants are an old habit I picked up long ago in my CS classes, there are a ton but it makes them stand out for me in vim more readily. I've come to avoid camel case since switching to linux, it makes searching for things a nightmare, it can also help if I group functions with a prefix (eg plot functions start with "pl_"). If you want I can find+replace all of the constants with their values... > > The differences came from using the unaveraged replicates in the larger script versus the average in the example. > > I'm not sure what you mean by "mean of the controls", would this simply be filtering out the control spots? > > I'll generate the heatmaps using the normalized expression values. I suppose in that case it would be worthwhile to filter ahead of time by p-value and fold change (the coef in topTable). > > > _________thanks_____________ > _______Joe______ > > -----Original Message----- > From: James W. MacDonald [mailto:jmacdon at uw.edu] > Sent: Thursday, March 28, 2013 10:08 AM > To: Cornish, Joseph (NIH/NIAID) [F] > Cc: bioconductor at r-project.org > Subject: Re: [BioC] Heatmaps for Agilent Single Color 4x44k using limma and vsn. > > Hi Joseph, > > > On 3/28/2013 9:27 AM, Cornish, Joseph (NIH/NIAID) [F] wrote: >> I am attempting to write a script that will allow me to do bulk benchmarking of various methods for normalization and background correction. I am concerned that I am not processing the data correctly and this may be leading to the abnormal heatmaps that I am seeing. For a test case, I am using an experiment where I have two different conditions with two different states with two replicates per state. >> >> The heatmaps concern me because of the uniformity between them, the values are all completely green or completely red, there are differences between genes within a state according to the heatmap. >> >> In the examples I have seen people have suggested that I use the coefficients from the eBayes fitting of the data, however in my example the eBayes LMarray object does not have two columns in the results, only on (not sure if this is a sign of errors on my part). I have instead used the lmFit coefficients since there are column for each matrix. >> >> To further confuse someone who is new to R, the simplified example and full code generate different heatmaps for the same data and same methods (normexp for background, cyclic loess for norm). The differences being in the simplified version the gene IDs show up on the right side but in the full version, only row number appears. The color for some genes changes between states as well. >> >> My questions are: >> Why does the bayes model only have one column? > Because you are computing one contrast (which is all you can do). One contrast means one test which means one column of results. In other words, the coefficient you are estimating in the contrasts.fit() step is the difference in means between the two groups, which is one value per gene. > >> Am I approaching the processing correctly? > I tried to read through that business below, but got some weird affliction due to an apparent overdose of underscores and had to stop. > >> Why are there differences between the resulting heatmaps of each version? > I don't know, as I worried my brain might explode if I read any more. > But you should note that it is fairly unusual to plot coefficients in a heatmap, as what you are then presenting are log ratios, and most people can't understand that. > > Instead you might consider getting the list of genes, and then creating a heatmap using the normalized expression values for each sample instead. You then would want to scale in some way. We often sweep the mean of the controls out of the matrix (by row), so the colors can then be interpreted as differences between the treated samples and controls. > > Best, > > Jim > > >> I know that the idea is to keep the amount of posted code to a minimum but I don't see any other way here. >> >> Here is the simplified example I have made: >> library(limma) >> library(vsn) >> library(gplots) >> >> set.seed(0xabcd) #from vsn manual >> >> #constants for reading agi outs from our imager >> COL__<- list(G="Raw intensity (tr.mean) {A}",Gb="Background (tr.mean) >> {A}") >> ANO__<- c("Position", "Name", "ID", "Description", "GeneName") >> >> #parse data >> targets<- readTargets("test.csv", sep = ",") >> gal<- readGAL("026806_D_20110524.gal") >> layout<- getLayout(gal) >> agi<- read.maimages(targets, columns = COL__, >> annotation = ANO__, >> green.only = TRUE) >> design<- model.matrix(~0+factor(c(1,1,2,2))) >> colnames(design)<- c("s1", "s2") >> >> #process >> bg<- backgroundCorrect(agi, method = "normexp") >> norm<- normalizeBetweenArrays(bg, method = "cyclicloess") >> avg<- avereps(norm, ID=norm$genes$ProbeName) >> >> #fit >> contmat<- makeContrasts(s2-s1, levels = design) >> fit_lm<- lmFit(avg$E, design) >> fit_be<- eBayes(contrasts.fit(fit_lm, contmat)) >> >> #results >> difexp<- topTable(fit_be, coef = 1, adjust = "BH") >> res<- decideTests(fit_be) >> comp<- vennDiagram(res) >> >> #plot >> pdf("example-heatmap.pdf") >> genes<- as.numeric(rownames(topTable(fit_be, n = 100))) >> heatmap(fit_lm$coefficients[genes,], col=redgreen(75), key = TRUE) >> dev.off() >> >> >> The full version: >> library(vsn) >> library(limma) >> library(gplots) >> set.seed(0xabcd) #from the vsn manual >> >> #constants >> A_SUB__<- 1 >> A_NXP__<- 2 >> A_MIN__<- 3 >> A_VSN__<- 1 >> A_QNT__<- 2 >> A_LWS__<- 3 >> M_SUB__<- "subtract" >> M_NXP__<- "normexp" >> M_MIN__<- "minimum" >> M_VSN__<- "vsn" >> M_LWS__<- "cyclicloess" >> M_QNT__<- "quantile" >> COL__<- list(G="Raw intensity (tr.mean) {A}",Gb="Background (tr.mean) >> {A}") >> ANO__<- c("Position", "Name", "ID", "Description", "GeneName") >> N_BGM__<- 3 #number of background correction methods >> N_NRM__<- 3 #number of normalization methods >> A_BGM__<- c(A_SUB__, A_NXP__, A_MIN__) >> A_NRM__<- c(A_VSN__, A_QNT__, A_LWS__) >> >> >> #utility methods >> #applies subtraction background correction with limma >> bg_subtract<- function(x){ >> return(list(norm = x$norm, >> bg = M_SUB__, >> agi = backgroundCorrect(x$agi, method = M_SUB__))) } >> >> #applies normexp backgroud correction with limma >> bg_normexp<- function(x){ >> return(list(norm = x$norm, >> bg = M_NXP__, >> agi = backgroundCorrect(x$agi, method = M_NXP__))) } >> >> #applies minimum background correction with limma >> bg_minimum<- function(x){ >> return(list(norm = x$norm, >> bg = M_MIN__, >> agi = backgroundCorrect(x$agi, method = M_MIN__))) } >> >> #applies VSN normalization with VSN package >> nb_vsn<- function(x){ >> return(list(norm = M_VSN__, >> bg = x$bg, >> agi = normalizeVSN(x$agi))) } >> >> #applies VSN normalization with VSN package >> nb_vsn<- function(x){ >> return(list(norm = M_VSN__, >> bg = x$bg, >> agi = normalizeVSN(x$agi))) } >> >> #applies cyclic lowess normalization with limma >> nb_lowess<- function(x){ >> return(list(norm = M_LWS__, >> bg = x$bg, >> agi = normalizeBetweenArrays(x$agi, method = >> M_LWS__))) } >> >> #applies quantile normalization with limma >> nb_quantile<- function(x){ >> return(list(norm = M_QNT__, >> bg = x$bg, >> agi = normalizeBetweenArrays(x$agi, method = >> M_QNT__))) } >> >> #generates a plot of MA values with limma >> pl_ma<- function(x){ >> name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-ma.pdf"), collapse="") >> pdf(name) >> plotMA(x[[1]]$agi$E) >> dev.off() >> } >> >> #generates a plot of mean versus stdev with VS >> pl_meansd<- function(x){ >> name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-meansd.pdf"), collapse="") >> pdf(name) >> meanSdPlot(x[[1]]$agi$E) >> dev.off() >> } >> >> #generate volcanoplot of log-fold versus log-odds >> pl_volcano<- function(x){ >> name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-volcano.pdf"), collapse="") >> pdf(name) >> volcanoplot(x[[1]]$bayes) >> dev.off() >> } >> >> #generate plot of expression data over time >> pl_lines<- function(x){ >> name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-lines.pdf"), collapse="") >> pdf(name) >> plotlines(x[[1]]$agi$E) >> dev.off() >> } >> >> #generate heatmat of final results >> pl_heatmap<- function(x, ntop){ >> name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-heatmap.pdf"), collapse="") >> pdf(name) >> genes<- as.numeric(rownames(topTable(x[[1]]$bayes, n = ntop))) >> heatmap(x[[1]]$lm$coefficients[genes,], >> col = redgreen(75), >> key = TRUE) >> dev.off() >> } >> >> #averages replicates of probes >> m_repavg<- function(x){ >> return(list(agi = x[[1]]$agi, >> avg = avereps(x[[1]]$agi, ID=x[[1]]$agi$genes$ProbeName), >> bg = x[[1]]$bg, >> norm = x[[1]]$norm)) >> } >> >> #generate contrasts >> m_fit<- function(x, design){ >> contmat<- makeContrasts(s2-s1, >> levels = design) >> fit1<- lmFit(x[[1]]$agi$E, design) >> fit2<- eBayes(contrasts.fit(fit1, contmat)) >> diff<- topTable(fit2, coef = 1, adjust = "BH") >> res<- decideTests(fit2) >> comp<- vennDiagram(res) >> return(list(agi = x[[1]]$agi, >> bg = x[[1]]$bg, >> norm = x[[1]]$norm, >> ave = x[[1]]$ave, >> cont = contmat, >> lm = fit1, >> bayes = fit2, >> diff = diff, >> res = res, >> comp = comp)) >> } >> >> ###################################################################### >> ######### # Function to parse array data into limma objects # Args: >> # tfile: target file (see limma documentation) >> # tsep: limma defaults TSV, have to set manually for CSV and so on >> # galfile: GAL file for printer/layout >> # Returns data object with targets, genes and printer layout >> ###################################################################### >> ######### >> parse<- function(tfile, tsep, galfile){ >> targets<- readTargets(tfile, sep=tsep) >> gal<- readGAL(galfile) >> layout<- getLayout(gal) >> agi<- read.maimages(targets, columns = COL__, >> annotation = ANO__, >> green.only = TRUE) >> bg<- "none" >> norm<- "none" >> desg<- model.matrix(~0+factor(c(1,1,2,2))) >> colnames(desg)<- c("s1","s2") >> d<- list(agi = agi, >> gal = gal, >> layout = layout, >> targets = targets, >> bg = bg, >> desg = desg, >> norm = norm) >> >> processed<- matrix(NA, nrow = N_BGM__, ncol = N_NRM__) >> processed<- apply(processed, 1:2, function(x){return(d)}) >> return(processed) >> } >> >> #parse command args >> args<- commandArgs(trailingOnly = TRUE) >> >> #load files into matrix >> a<- parse(args[1], args[2], args[3]) >> >> #keepy a copy of a raw dataset >> raw<- a[1,1] >> print(raw[[1]]$desg) >> >> #perform backgroung corrections >> a[,A_SUB__]<- lapply(a[,A_SUB__], bg_subtract) >> a[,A_NXP__]<- lapply(a[,A_NXP__], bg_normexp) >> a[,A_MIN__]<- lapply(a[,A_MIN__], bg_minimum) >> >> #perform normalization >> a[A_VSN__,]<- lapply(a[A_VSN__,], nb_vsn) >> a[A_LWS__,]<- lapply(a[A_LWS__,], nb_lowess) >> a[A_QNT__,]<- lapply(a[A_QNT__,], nb_quantile) >> >> #fit model and differentiate expression >> a<- apply(a, 1:2, m_repavg) >> a<- apply(a, 1:2, m_fit, design = raw[[1]]$desg) >> >> #generate plots >> apply(a, 1:2, pl_ma) >> apply(a, 1:2, pl_meansd) >> apply(a, 1:2, pl_volcano) >> apply(a, 1:2, pl_lines) >> apply(a, 1:2, pl_heatmap, args[4]) >> >> >> [[alternative HTML version deleted]] >> >> _______________________________________________ >> Bioconductor mailing list >> Bioconductor at r-project.org >> https://stat.ethz.ch/mailman/listinfo/bioconductor >> Search the archives: >> http://news.gmane.org/gmane.science.biology.informatics.conductor > -- > James W. MacDonald, M.S. > Biostatistician > University of Washington > Environmental and Occupational Health Sciences > 4225 Roosevelt Way NE, # 100 > Seattle WA 98105-6099 > -- James W. MacDonald, M.S. Biostatistician University of Washington Environmental and Occupational Health Sciences 4225 Roosevelt Way NE, # 100 Seattle WA 98105-6099 ------------------------------ Message: 6 Date: Thu, 28 Mar 2013 15:00:49 +0000 From: "Cornish, Joseph (NIH/NIAID) [F]" <joseph.cornish@nih.gov> To: "Davis, Sean (NIH/NCI) [E]" <sdavis2 at="" mail.nih.gov=""> Cc: "bioconductor at r-project.org" <bioconductor at="" r-project.org="">, "James W. MacDonald" <jmacdon at="" uw.edu=""> Subject: Re: [BioC] Heatmaps for Agilent Single Color 4x44k using limma and vsn. Message-ID: <e68ea43fbf80174c85a92bc24efaa3b84b7c02 at="" mlbxv01.nih.gov=""> Content-Type: text/plain Hello Davis, I posted a simple example, which after fixing a problem in the longer example no produces the same results. I only posted the longer one because I wasn???t able to reproduce the results of the longer version initially. I???ve been sort of banging my head on this for a while so I missed the difference between ???agi??? and ???avg???. I???m thankful anyone bothered to look after posting the code, so thanks. If I sounded offended by the underscores, I was only admitting that I know they???re terrible. c: Thanks for clarifying the subtraction of the mean aspect. That has been one aspect that has been troubling when it comes to learning this analysis, there seem to be plenty of ways to do things but, without experience, no way to tell what is better. I will give the mean subtraction and the z-score scaling a try. Am I reading too much into the usefulness of these heatmaps? In most papers I???ve seen they seem to be there because the look nice or the attached clustering might be of some use. Thanks again! From: Davis, Sean (NCI) On Behalf Of Davis, Sean (NIH/NCI) [E] Sent: Thursday, March 28, 2013 10:34 AM To: Cornish, Joseph (NIH/NIAID) [F] Cc: James W. MacDonald; bioconductor at r-project.org Subject: Re: [BioC] Heatmaps for Agilent Single Color 4x44k using limma and vsn. On Thu, Mar 28, 2013 at 10:22 AM, Cornish, Joseph (NIH/NIAID) [F] <joseph.cornish at="" nih.gov<mailto:joseph.cornish="" at="" nih.gov="">> wrote: Hello James, Are you referring to function names or constants? The constants are an old habit I picked up long ago in my CS classes, there are a ton but it makes them stand out for me in vim more readily. I've come to avoid camel case since switching to linux, it makes searching for things a nightmare, it can also help if I group functions with a prefix (eg plot functions start with "pl_"). If you want I can find+replace all of the constants with their values... Hi, Joe. I think I could help to translate part of what James said. If you could boil your example down to a few lines of reproducible code that show the problem, that would be simpler to interpret and easier to comment on directly. Sometimes in the process of producing such a simplified example, I have found my own misunderstandings. It also sometimes helps to get someone local to give a second pair of eyes. The differences came from using the unaveraged replicates in the larger script versus the average in the example. I'm not sure what you mean by "mean of the controls", would this simply be filtering out the control spots? No. I think he means for you to subtract the mean of the control samples (or one group of samples) for each gene from all the normalized expression values. I'll generate the heatmaps using the normalized expression values. I suppose in that case it would be worthwhile to filter ahead of time by p-value and fold change (the coef in topTable). Yes. Sean _________thanks_____________ _______Joe______ -----Original Message----- From: James W. MacDonald [mailto:jmacdon@uw.edu<mailto:jmacdon@uw.edu>] Sent: Thursday, March 28, 2013 10:08 AM To: Cornish, Joseph (NIH/NIAID) [F] Cc: bioconductor at r-project.org<mailto:bioconductor at="" r-project.org=""> Subject: Re: [BioC] Heatmaps for Agilent Single Color 4x44k using limma and vsn. Hi Joseph, On 3/28/2013 9:27 AM, Cornish, Joseph (NIH/NIAID) [F] wrote: > I am attempting to write a script that will allow me to do bulk benchmarking of various methods for normalization and background correction. I am concerned that I am not processing the data correctly and this may be leading to the abnormal heatmaps that I am seeing. For a test case, I am using an experiment where I have two different conditions with two different states with two replicates per state. > > The heatmaps concern me because of the uniformity between them, the values are all completely green or completely red, there are differences between genes within a state according to the heatmap. > > In the examples I have seen people have suggested that I use the coefficients from the eBayes fitting of the data, however in my example the eBayes LMarray object does not have two columns in the results, only on (not sure if this is a sign of errors on my part). I have instead used the lmFit coefficients since there are column for each matrix. > > To further confuse someone who is new to R, the simplified example and full code generate different heatmaps for the same data and same methods (normexp for background, cyclic loess for norm). The differences being in the simplified version the gene IDs show up on the right side but in the full version, only row number appears. The color for some genes changes between states as well. > > My questions are: > Why does the bayes model only have one column? Because you are computing one contrast (which is all you can do). One contrast means one test which means one column of results. In other words, the coefficient you are estimating in the contrasts.fit() step is the difference in means between the two groups, which is one value per gene. > Am I approaching the processing correctly? I tried to read through that business below, but got some weird affliction due to an apparent overdose of underscores and had to stop. > Why are there differences between the resulting heatmaps of each version? I don't know, as I worried my brain might explode if I read any more. But you should note that it is fairly unusual to plot coefficients in a heatmap, as what you are then presenting are log ratios, and most people can't understand that. Instead you might consider getting the list of genes, and then creating a heatmap using the normalized expression values for each sample instead. You then would want to scale in some way. We often sweep the mean of the controls out of the matrix (by row), so the colors can then be interpreted as differences between the treated samples and controls. Best, Jim > > I know that the idea is to keep the amount of posted code to a minimum but I don't see any other way here. > > Here is the simplified example I have made: > library(limma) > library(vsn) > library(gplots) > > set.seed(0xabcd) #from vsn manual > > #constants for reading agi outs from our imager > COL__<- list(G="Raw intensity (tr.mean) {A}",Gb="Background (tr.mean) > {A}") > ANO__<- c("Position", "Name", "ID", "Description", "GeneName") > > #parse data > targets<- readTargets("test.csv", sep = ",") > gal<- readGAL("026806_D_20110524.gal") > layout<- getLayout(gal) > agi<- read.maimages(targets, columns = COL__, > annotation = ANO__, > green.only = TRUE) > design<- model.matrix(~0+factor(c(1,1,2,2))) > colnames(design)<- c("s1", "s2") > > #process > bg<- backgroundCorrect(agi, method = "normexp") > norm<- normalizeBetweenArrays(bg, method = "cyclicloess") > avg<- avereps(norm, ID=norm$genes$ProbeName) > > #fit > contmat<- makeContrasts(s2-s1, levels = design) > fit_lm<- lmFit(avg$E, design) > fit_be<- eBayes(contrasts.fit(fit_lm, contmat)) > > #results > difexp<- topTable(fit_be, coef = 1, adjust = "BH") > res<- decideTests(fit_be) > comp<- vennDiagram(res) > > #plot > pdf("example-heatmap.pdf") > genes<- as.numeric(rownames(topTable(fit_be, n = 100))) > heatmap(fit_lm$coefficients[genes,], col=redgreen(75), key = TRUE) > dev.off() > > > The full version: > library(vsn) > library(limma) > library(gplots) > set.seed(0xabcd) #from the vsn manual > > #constants > A_SUB__<- 1 > A_NXP__<- 2 > A_MIN__<- 3 > A_VSN__<- 1 > A_QNT__<- 2 > A_LWS__<- 3 > M_SUB__<- "subtract" > M_NXP__<- "normexp" > M_MIN__<- "minimum" > M_VSN__<- "vsn" > M_LWS__<- "cyclicloess" > M_QNT__<- "quantile" > COL__<- list(G="Raw intensity (tr.mean) {A}",Gb="Background (tr.mean) > {A}") > ANO__<- c("Position", "Name", "ID", "Description", "GeneName") > N_BGM__<- 3 #number of background correction methods > N_NRM__<- 3 #number of normalization methods > A_BGM__<- c(A_SUB__, A_NXP__, A_MIN__) > A_NRM__<- c(A_VSN__, A_QNT__, A_LWS__) > > > #utility methods > #applies subtraction background correction with limma > bg_subtract<- function(x){ > return(list(norm = x$norm, > bg = M_SUB__, > agi = backgroundCorrect(x$agi, method = M_SUB__))) } > > #applies normexp backgroud correction with limma > bg_normexp<- function(x){ > return(list(norm = x$norm, > bg = M_NXP__, > agi = backgroundCorrect(x$agi, method = M_NXP__))) } > > #applies minimum background correction with limma > bg_minimum<- function(x){ > return(list(norm = x$norm, > bg = M_MIN__, > agi = backgroundCorrect(x$agi, method = M_MIN__))) } > > #applies VSN normalization with VSN package > nb_vsn<- function(x){ > return(list(norm = M_VSN__, > bg = x$bg, > agi = normalizeVSN(x$agi))) } > > #applies VSN normalization with VSN package > nb_vsn<- function(x){ > return(list(norm = M_VSN__, > bg = x$bg, > agi = normalizeVSN(x$agi))) } > > #applies cyclic lowess normalization with limma > nb_lowess<- function(x){ > return(list(norm = M_LWS__, > bg = x$bg, > agi = normalizeBetweenArrays(x$agi, method = > M_LWS__))) } > > #applies quantile normalization with limma > nb_quantile<- function(x){ > return(list(norm = M_QNT__, > bg = x$bg, > agi = normalizeBetweenArrays(x$agi, method = > M_QNT__))) } > > #generates a plot of MA values with limma > pl_ma<- function(x){ > name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-ma.pdf"), collapse="") > pdf(name) > plotMA(x[[1]]$agi$E) > dev.off() > } > > #generates a plot of mean versus stdev with VS > pl_meansd<- function(x){ > name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-meansd.pdf"), collapse="") > pdf(name) > meanSdPlot(x[[1]]$agi$E) > dev.off() > } > > #generate volcanoplot of log-fold versus log-odds > pl_volcano<- function(x){ > name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-volcano.pdf"), collapse="") > pdf(name) > volcanoplot(x[[1]]$bayes) > dev.off() > } > > #generate plot of expression data over time > pl_lines<- function(x){ > name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-lines.pdf"), collapse="") > pdf(name) > plotlines(x[[1]]$agi$E) > dev.off() > } > > #generate heatmat of final results > pl_heatmap<- function(x, ntop){ > name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-heatmap.pdf"), collapse="") > pdf(name) > genes<- as.numeric(rownames(topTable(x[[1]]$bayes, n = ntop))) > heatmap(x[[1]]$lm$coefficients[genes,], > col = redgreen(75), > key = TRUE) > dev.off() > } > > #averages replicates of probes > m_repavg<- function(x){ > return(list(agi = x[[1]]$agi, > avg = avereps(x[[1]]$agi, ID=x[[1]]$agi$genes$ProbeName), > bg = x[[1]]$bg, > norm = x[[1]]$norm)) > } > > #generate contrasts > m_fit<- function(x, design){ > contmat<- makeContrasts(s2-s1, > levels = design) > fit1<- lmFit(x[[1]]$agi$E, design) > fit2<- eBayes(contrasts.fit(fit1, contmat)) > diff<- topTable(fit2, coef = 1, adjust = "BH") > res<- decideTests(fit2) > comp<- vennDiagram(res) > return(list(agi = x[[1]]$agi, > bg = x[[1]]$bg, > norm = x[[1]]$norm, > ave = x[[1]]$ave, > cont = contmat, > lm = fit1, > bayes = fit2, > diff = diff, > res = res, > comp = comp)) > } > > ###################################################################### > ######### # Function to parse array data into limma objects # Args: > # tfile: target file (see limma documentation) > # tsep: limma defaults TSV, have to set manually for CSV and so on > # galfile: GAL file for printer/layout > # Returns data object with targets, genes and printer layout > ###################################################################### > ######### > parse<- function(tfile, tsep, galfile){ > targets<- readTargets(tfile, sep=tsep) > gal<- readGAL(galfile) > layout<- getLayout(gal) > agi<- read.maimages(targets, columns = COL__, > annotation = ANO__, > green.only = TRUE) > bg<- "none" > norm<- "none" > desg<- model.matrix(~0+factor(c(1,1,2,2))) > colnames(desg)<- c("s1","s2") > d<- list(agi = agi, > gal = gal, > layout = layout, > targets = targets, > bg = bg, > desg = desg, > norm = norm) > > processed<- matrix(NA, nrow = N_BGM__, ncol = N_NRM__) > processed<- apply(processed, 1:2, function(x){return(d)}) > return(processed) > } > > #parse command args > args<- commandArgs(trailingOnly = TRUE) > > #load files into matrix > a<- parse(args[1], args[2], args[3]) > > #keepy a copy of a raw dataset > raw<- a[1,1] > print(raw[[1]]$desg) > > #perform backgroung corrections > a[,A_SUB__]<- lapply(a[,A_SUB__], bg_subtract) > a[,A_NXP__]<- lapply(a[,A_NXP__], bg_normexp) > a[,A_MIN__]<- lapply(a[,A_MIN__], bg_minimum) > > #perform normalization > a[A_VSN__,]<- lapply(a[A_VSN__,], nb_vsn) > a[A_LWS__,]<- lapply(a[A_LWS__,], nb_lowess) > a[A_QNT__,]<- lapply(a[A_QNT__,], nb_quantile) > > #fit model and differentiate expression > a<- apply(a, 1:2, m_repavg) > a<- apply(a, 1:2, m_fit, design = raw[[1]]$desg) > > #generate plots > apply(a, 1:2, pl_ma) > apply(a, 1:2, pl_meansd) > apply(a, 1:2, pl_volcano) > apply(a, 1:2, pl_lines) > apply(a, 1:2, pl_heatmap, args[4]) > > > [[alternative HTML version deleted]] > > _______________________________________________ > Bioconductor mailing list > Bioconductor at r-project.org<mailto:bioconductor at="" r-project.org=""> > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: > http://news.gmane.org/gmane.science.biology.informatics.conductor -- James W. MacDonald, M.S. Biostatistician University of Washington Environmental and Occupational Health Sciences 4225 Roosevelt Way NE, # 100 Seattle WA 98105-6099 _______________________________________________ Bioconductor mailing list Bioconductor at r-project.org<mailto:bioconductor at="" r-project.org=""> https://stat.ethz.ch/mailman/listinfo/bioconductor Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor [[alternative HTML version deleted]] ------------------------------ Message: 7 Date: Thu, 28 Mar 2013 15:18:59 +0000 From: Andreia Fonseca <andreia.fonseca@gmail.com> To: bioconductor <bioconductor at="" stat.math.ethz.ch=""> Subject: [BioC] analysis of HuGene2.0-st array affymetrix -aroma package Message-ID: <cag43sel74j9eypo8izygsy7=c0jye5bjw5egcohkug9mmyp_fq at="" mail.gmail.com=""> Content-Type: text/plain Dear all, I am analysing HuGene2.0-st array affymetrix using aroma and limma package. I have 2 treatments and 3 arrays per condition. I have getting 11 genes that are differentially expressed, but these loose significance after applying fdr correction. I want to apply a filter out genes with similar distributions between treatments before applying fdr correction to see I am loosing significance due to overcorrection. Can someone advise me on a function or package to do this? I have explored genefilter and I could not find a function for this. Thanks for the help. Kind regards, Andreia ---------------------------------------------------------------------- ------------------------- Andreia J. Amaral, PhD BioFIG - Center for Biodiversity, Functional and Integrative Genomics Instituto de Medicina Molecular University of Lisbon Tel: +352 217500000 (ext. office: 28253) email:andreiaamaral at fm.ul.pt ; andreiaamaral at fc.ul.pt [[alternative HTML version deleted]] ------------------------------ Message: 8 Date: Thu, 28 Mar 2013 11:47:12 -0400 From: Sean Davis <sdavis2@mail.nih.gov> To: "Cornish, Joseph (NIH/NIAID) [F]" <joseph.cornish at="" nih.gov=""> Cc: "bioconductor at r-project.org" <bioconductor at="" r-project.org="">, "James W. MacDonald" <jmacdon at="" uw.edu=""> Subject: Re: [BioC] Heatmaps for Agilent Single Color 4x44k using limma and vsn. Message-ID: <caneavb=xn9c2bjqk6gf5dptghbltgrcov+oh0t05qsmd1s92dg at="" mail.gmail.com=""> Content-Type: text/plain; charset=UTF-8 On Thu, Mar 28, 2013 at 11:00 AM, Cornish, Joseph (NIH/NIAID) [F] <joseph.cornish at="" nih.gov=""> wrote: > Hello Davis, > > I posted a simple example, which after fixing a problem in the longer example no produces the same results. I only posted the longer one because I wasn?t able to reproduce the results of the longer version initially. I?ve been sort of banging my head on this for a while so I missed the difference between ?agi? and ?avg?. I?m thankful anyone bothered to look after posting the code, so thanks. If I sounded offended by the underscores, I was only admitting that I know they?re terrible. c: > > Thanks for clarifying the subtraction of the mean aspect. That has been one aspect that has been troubling when it comes to learning this analysis, there seem to be plenty of ways to do things but, without experience, no way to tell what is better. > There is not a "right" way to do this since the only metric I have seen for the evaluation of the color scale in a heatmap is subjective and qualitative. > I will give the mean subtraction and the z-score scaling a try. Am I reading too much into the usefulness of these heatmaps? In most papers I?ve seen they seem to be there because the look nice or the attached clustering might be of some use. > Heatmaps are very good visualization tools. I would not try to attach more utility to them than that. Sean > Thanks again! > > From: Davis, Sean (NCI) On Behalf Of Davis, Sean (NIH/NCI) [E] > Sent: Thursday, March 28, 2013 10:34 AM > To: Cornish, Joseph (NIH/NIAID) [F] > Cc: James W. MacDonald; bioconductor at r-project.org > Subject: Re: [BioC] Heatmaps for Agilent Single Color 4x44k using limma and vsn. > > > > On Thu, Mar 28, 2013 at 10:22 AM, Cornish, Joseph (NIH/NIAID) [F] <joseph.cornish at="" nih.gov<mailto:joseph.cornish="" at="" nih.gov="">> wrote: > Hello James, > > Are you referring to function names or constants? The constants are an old habit I picked up long ago in my CS classes, there are a ton but it makes them stand out for me in vim more readily. I've come to avoid camel case since switching to linux, it makes searching for things a nightmare, it can also help if I group functions with a prefix (eg plot functions start with "pl_"). If you want I can find+replace all of the constants with their values... > > Hi, Joe. > > I think I could help to translate part of what James said. If you could boil your example down to a few lines of reproducible code that show the problem, that would be simpler to interpret and easier to comment on directly. Sometimes in the process of producing such a simplified example, I have found my own misunderstandings. It also sometimes helps to get someone local to give a second pair of eyes. > > The differences came from using the unaveraged replicates in the larger script versus the average in the example. > > I'm not sure what you mean by "mean of the controls", would this simply be filtering out the control spots? > > No. I think he means for you to subtract the mean of the control samples (or one group of samples) for each gene from all the normalized expression values. > > I'll generate the heatmaps using the normalized expression values. I suppose in that case it would be worthwhile to filter ahead of time by p-value and fold change (the coef in topTable). > > Yes. > > Sean > > > _________thanks_____________ > _______Joe______ > > -----Original Message----- > From: James W. MacDonald [mailto:jmacdon at uw.edu<mailto:jmacdon at="" uw.edu="">] > Sent: Thursday, March 28, 2013 10:08 AM > To: Cornish, Joseph (NIH/NIAID) [F] > Cc: bioconductor at r-project.org<mailto:bioconductor at="" r-project.org=""> > Subject: Re: [BioC] Heatmaps for Agilent Single Color 4x44k using limma and vsn. > > Hi Joseph, > > > On 3/28/2013 9:27 AM, Cornish, Joseph (NIH/NIAID) [F] wrote: >> I am attempting to write a script that will allow me to do bulk benchmarking of various methods for normalization and background correction. I am concerned that I am not processing the data correctly and this may be leading to the abnormal heatmaps that I am seeing. For a test case, I am using an experiment where I have two different conditions with two different states with two replicates per state. >> >> The heatmaps concern me because of the uniformity between them, the values are all completely green or completely red, there are differences between genes within a state according to the heatmap. >> >> In the examples I have seen people have suggested that I use the coefficients from the eBayes fitting of the data, however in my example the eBayes LMarray object does not have two columns in the results, only on (not sure if this is a sign of errors on my part). I have instead used the lmFit coefficients since there are column for each matrix. >> >> To further confuse someone who is new to R, the simplified example and full code generate different heatmaps for the same data and same methods (normexp for background, cyclic loess for norm). The differences being in the simplified version the gene IDs show up on the right side but in the full version, only row number appears. The color for some genes changes between states as well. >> >> My questions are: >> Why does the bayes model only have one column? > > Because you are computing one contrast (which is all you can do). One contrast means one test which means one column of results. In other words, the coefficient you are estimating in the contrasts.fit() step is the difference in means between the two groups, which is one value per gene. > >> Am I approaching the processing correctly? > > I tried to read through that business below, but got some weird affliction due to an apparent overdose of underscores and had to stop. > >> Why are there differences between the resulting heatmaps of each version? > > I don't know, as I worried my brain might explode if I read any more. > But you should note that it is fairly unusual to plot coefficients in a heatmap, as what you are then presenting are log ratios, and most people can't understand that. > > Instead you might consider getting the list of genes, and then creating a heatmap using the normalized expression values for each sample instead. You then would want to scale in some way. We often sweep the mean of the controls out of the matrix (by row), so the colors can then be interpreted as differences between the treated samples and controls. > > Best, > > Jim > > >> >> I know that the idea is to keep the amount of posted code to a minimum but I don't see any other way here. >> >> Here is the simplified example I have made: >> library(limma) >> library(vsn) >> library(gplots) >> >> set.seed(0xabcd) #from vsn manual >> >> #constants for reading agi outs from our imager >> COL__<- list(G="Raw intensity (tr.mean) {A}",Gb="Background (tr.mean) >> {A}") >> ANO__<- c("Position", "Name", "ID", "Description", "GeneName") >> >> #parse data >> targets<- readTargets("test.csv", sep = ",") >> gal<- readGAL("026806_D_20110524.gal") >> layout<- getLayout(gal) >> agi<- read.maimages(targets, columns = COL__, >> annotation = ANO__, >> green.only = TRUE) >> design<- model.matrix(~0+factor(c(1,1,2,2))) >> colnames(design)<- c("s1", "s2") >> >> #process >> bg<- backgroundCorrect(agi, method = "normexp") >> norm<- normalizeBetweenArrays(bg, method = "cyclicloess") >> avg<- avereps(norm, ID=norm$genes$ProbeName) >> >> #fit >> contmat<- makeContrasts(s2-s1, levels = design) >> fit_lm<- lmFit(avg$E, design) >> fit_be<- eBayes(contrasts.fit(fit_lm, contmat)) >> >> #results >> difexp<- topTable(fit_be, coef = 1, adjust = "BH") >> res<- decideTests(fit_be) >> comp<- vennDiagram(res) >> >> #plot >> pdf("example-heatmap.pdf") >> genes<- as.numeric(rownames(topTable(fit_be, n = 100))) >> heatmap(fit_lm$coefficients[genes,], col=redgreen(75), key = TRUE) >> dev.off() >> >> >> The full version: >> library(vsn) >> library(limma) >> library(gplots) >> set.seed(0xabcd) #from the vsn manual >> >> #constants >> A_SUB__<- 1 >> A_NXP__<- 2 >> A_MIN__<- 3 >> A_VSN__<- 1 >> A_QNT__<- 2 >> A_LWS__<- 3 >> M_SUB__<- "subtract" >> M_NXP__<- "normexp" >> M_MIN__<- "minimum" >> M_VSN__<- "vsn" >> M_LWS__<- "cyclicloess" >> M_QNT__<- "quantile" >> COL__<- list(G="Raw intensity (tr.mean) {A}",Gb="Background (tr.mean) >> {A}") >> ANO__<- c("Position", "Name", "ID", "Description", "GeneName") >> N_BGM__<- 3 #number of background correction methods >> N_NRM__<- 3 #number of normalization methods >> A_BGM__<- c(A_SUB__, A_NXP__, A_MIN__) >> A_NRM__<- c(A_VSN__, A_QNT__, A_LWS__) >> >> >> #utility methods >> #applies subtraction background correction with limma >> bg_subtract<- function(x){ >> return(list(norm = x$norm, >> bg = M_SUB__, >> agi = backgroundCorrect(x$agi, method = M_SUB__))) } >> >> #applies normexp backgroud correction with limma >> bg_normexp<- function(x){ >> return(list(norm = x$norm, >> bg = M_NXP__, >> agi = backgroundCorrect(x$agi, method = M_NXP__))) } >> >> #applies minimum background correction with limma >> bg_minimum<- function(x){ >> return(list(norm = x$norm, >> bg = M_MIN__, >> agi = backgroundCorrect(x$agi, method = M_MIN__))) } >> >> #applies VSN normalization with VSN package >> nb_vsn<- function(x){ >> return(list(norm = M_VSN__, >> bg = x$bg, >> agi = normalizeVSN(x$agi))) } >> >> #applies VSN normalization with VSN package >> nb_vsn<- function(x){ >> return(list(norm = M_VSN__, >> bg = x$bg, >> agi = normalizeVSN(x$agi))) } >> >> #applies cyclic lowess normalization with limma >> nb_lowess<- function(x){ >> return(list(norm = M_LWS__, >> bg = x$bg, >> agi = normalizeBetweenArrays(x$agi, method = >> M_LWS__))) } >> >> #applies quantile normalization with limma >> nb_quantile<- function(x){ >> return(list(norm = M_QNT__, >> bg = x$bg, >> agi = normalizeBetweenArrays(x$agi, method = >> M_QNT__))) } >> >> #generates a plot of MA values with limma >> pl_ma<- function(x){ >> name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-ma.pdf"), collapse="") >> pdf(name) >> plotMA(x[[1]]$agi$E) >> dev.off() >> } >> >> #generates a plot of mean versus stdev with VS >> pl_meansd<- function(x){ >> name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-meansd.pdf"), collapse="") >> pdf(name) >> meanSdPlot(x[[1]]$agi$E) >> dev.off() >> } >> >> #generate volcanoplot of log-fold versus log-odds >> pl_volcano<- function(x){ >> name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-volcano.pdf"), collapse="") >> pdf(name) >> volcanoplot(x[[1]]$bayes) >> dev.off() >> } >> >> #generate plot of expression data over time >> pl_lines<- function(x){ >> name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-lines.pdf"), collapse="") >> pdf(name) >> plotlines(x[[1]]$agi$E) >> dev.off() >> } >> >> #generate heatmat of final results >> pl_heatmap<- function(x, ntop){ >> name<- paste(c(x[[1]]$bg, x[[1]]$norm, "-heatmap.pdf"), collapse="") >> pdf(name) >> genes<- as.numeric(rownames(topTable(x[[1]]$bayes, n = ntop))) >> heatmap(x[[1]]$lm$coefficients[genes,], >> col = redgreen(75), >> key = TRUE) >> dev.off() >> } >> >> #averages replicates of probes >> m_repavg<- function(x){ >> return(list(agi = x[[1]]$agi, >> avg = avereps(x[[1]]$agi, ID=x[[1]]$agi$genes$ProbeName), >> bg = x[[1]]$bg, >> norm = x[[1]]$norm)) >> } >> >> #generate contrasts >> m_fit<- function(x, design){ >> contmat<- makeContrasts(s2-s1, >> levels = design) >> fit1<- lmFit(x[[1]]$agi$E, design) >> fit2<- eBayes(contrasts.fit(fit1, contmat)) >> diff<- topTable(fit2, coef = 1, adjust = "BH") >> res<- decideTests(fit2) >> comp<- vennDiagram(res) >> return(list(agi = x[[1]]$agi, >> bg = x[[1]]$bg, >> norm = x[[1]]$norm, >> ave = x[[1]]$ave, >> cont = contmat, >> lm = fit1, >> bayes = fit2, >> diff = diff, >> res = res, >> comp = comp)) >> } >> >> ###################################################################### >> ######### # Function to parse array data into limma objects # Args: >> # tfile: target file (see limma documentation) >> # tsep: limma defaults TSV, have to set manually for CSV and so on >> # galfile: GAL file for printer/layout >> # Returns data object with targets, genes and printer layout >> ###################################################################### >> ######### >> parse<- function(tfile, tsep, galfile){ >> targets<- readTargets(tfile, sep=tsep) >> gal<- readGAL(galfile) >> layout<- getLayout(gal) >> agi<- read.maimages(targets, columns = COL__, >> annotation = ANO__, >> green.only = TRUE) >> bg<- "none" >> norm<- "none" >> desg<- model.matrix(~0+factor(c(1,1,2,2))) >> colnames(desg)<- c("s1","s2") >> d<- list(agi = agi, >> gal = gal, >> layout = layout, >> targets = targets, >> bg = bg, >> desg = desg, >> norm = norm) >> >> processed<- matrix(NA, nrow = N_BGM__, ncol = N_NRM__) >> processed<- apply(processed, 1:2, function(x){return(d)}) >> return(processed) >> } >> >> #parse command args >> args<- commandArgs(trailingOnly = TRUE) >> >> #load files into matrix >> a<- parse(args[1], args[2], args[3]) >> >> #keepy a copy of a raw dataset >> raw<- a[1,1] >> print(raw[[1]]$desg) >> >> #perform backgroung corrections >> a[,A_SUB__]<- lapply(a[,A_SUB__], bg_subtract) >> a[,A_NXP__]<- lapply(a[,A_NXP__], bg_normexp) >> a[,A_MIN__]<- lapply(a[,A_MIN__], bg_minimum) >> >> #perform normalization >> a[A_VSN__,]<- lapply(a[A_VSN__,], nb_vsn) >> a[A_LWS__,]<- lapply(a[A_LWS__,], nb_lowess) >> a[A_QNT__,]<- lapply(a[A_QNT__,], nb_quantile) >> >> #fit model and differentiate expression >> a<- apply(a, 1:2, m_repavg) >> a<- apply(a, 1:2, m_fit, design = raw[[1]]$desg) >> >> #generate plots >> apply(a, 1:2, pl_ma) >> apply(a, 1:2, pl_meansd) >> apply(a, 1:2, pl_volcano) >> apply(a, 1:2, pl_lines) >> apply(a, 1:2, pl_heatmap, args[4]) >> >> >> [[alternative HTML version deleted]] >> >> _______________________________________________ >> Bioconductor mailing list >> Bioconductor at r-project.org<mailto:bioconductor at="" r-project.org=""> >> https://stat.ethz.ch/mailman/listinfo/bioconductor >> Search the archives: >> http://news.gmane.org/gmane.science.biology.informatics.conductor > > -- > James W. MacDonald, M.S. > Biostatistician > University of Washington > Environmental and Occupational Health Sciences > 4225 Roosevelt Way NE, # 100 > Seattle WA 98105-6099 > > _______________________________________________ > Bioconductor mailing list > Bioconductor at r-project.org<mailto:bioconductor at="" r-project.org=""> > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor > > > [[alternative HTML version deleted]] > > > _______________________________________________ > Bioconductor mailing list > Bioconductor at r-project.org > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor ------------------------------ Message: 9 Date: Thu, 28 Mar 2013 11:00:17 -0500 From: Michael Muratet <mmuratet@hudsonalpha.org> To: bioconductor at r-project.org Subject: [BioC] Interpreting DESeq2 results Message-ID: <14DFA994-3719-4FDF-BD71-4290D4149457 at hudsonalpha.org> Content-Type: text/plain; charset=us-ascii Greetings I have an experiment: > design(dse) ~ factor1 + factor2 + factor3 where factor1 has two levels, factor2 has three levels and factor3 has three levels. I extract a gene of interest from the results for each term (I've changed the indices to reflect the condition): > lapply(resultsNames(dse),function(u) results(dse,u)["gene_A",]) [["Intercept"]] baseMean log2FoldChange pvalue FDR gene_A 1596.548 10.77485 3.309439e-216 7.025442e-216 [["factor1_level2"]] baseMean log2FoldChange pvalue FDR gene_A 1596.548 0.3386776 0.1307309 0.3587438 [["factor2_level2"]] baseMean log2FoldChange pvalue FDR gene_A 1596.548 -0.6882543 0.0613569 0.1007896 [["factor2_level3"]] baseMean log2FoldChange pvalue FDR gene_A 1596.548 0.2393368 0.513216 0.6589575 [["factor3_level2"]] baseMean log2FoldChange pvalue FDR gene_A 1596.548 0.1584153 0.6423634 0.8503163 [["factor3_level3]] baseMean log2FoldChange pvalue FDR gene_A 1596.548 -1.627898 1.823141e-06 0.001409384 I want to be sure I understand the output format. Is it true that the coefficients (the vector beta) from the fit are the baseMean value scaled by the log2FoldChange? Is the true intercept value 1596.548*2^10.77485=2797274.13? mcols() tells me that the baseMean term is calculated over "all rows". The baseMean is different for different genes although it is the same for each gene across all the conditions, I'm not seeing how the rows are selected. Thanks Mike Michael Muratet, Ph.D. Senior Scientist HudsonAlpha Institute for Biotechnology mmuratet at hudsonalpha.org (256) 327-0473 (p) (256) 327-0966 (f) Room 4005 601 Genome Way Huntsville, Alabama 35806 ------------------------------ Message: 10 Date: Thu, 28 Mar 2013 16:30:12 +0000 From: Fiona Ingleby <fiona.ingleby@gmail.com> To: "bioconductor at r-project.org" <bioconductor at="" r-project.org=""> Subject: [BioC] identifying drosophila miRNA targets Message-ID: <04CBBC25-9E0F-4E2B-8008-D3EB093E2A10 at gmail.com> Content-Type: text/plain Hi everyone, I am working with mRNA data from Affy 'drosophila2' arrays and miRNA data from Affy 'mirna3' arrays. I have identified a list of differentially expressed mRNAs and miRNAs. I'm having a bit of trouble with some downstream analyses and I'm hoping someone might be able to offer some help. I would like to use my list of differentially expressed miRNAs to access online databases (e.g. miRBase, microRNA.org?) and extract the names of all the potential target mRNAs. Then I'd like to use this list of mRNAs to look through my mRNA expression data. I'm aware of packages like 'RmiR' and 'microRNA' which have built-in functions for finding miRNA targets, but as far as I can tell, 'RmiR' uses miRNA databases for humans only and 'microRNA' works with human and mouse data only. So is there a package I am unaware of (or another application of 'RmiR'/'microRNA' that I am unaware of) for looking at drosophila data? So far I have also considered the 'biomaRt' package to see if the database query function on there can help me, but I haven't had much luck. For instance, if I try an example list of miRNAs: mirna<-c("dme-miR-1002","dme-miR-312","dme-miR-973") library(biomaRt) ensembl<-useMart("ensembl",dataset="dmelanogaster_gene_ensembl") getBM(attributes="mirbase_accession",filters="mirbase_id",values=mirna ,mart=ensembl) then 'logical(0)' is returned, as if there are no records for those miRNAs - but by searching the database manually I know the records are there. Alternatively I can try: miRNA <- getBM(c("mirbase_accession","mirbase_id", "ensembl_gene_id", "start_position", "chromosome_name"), filters = c("with_mirbase"), values = list(T), mart = ensembl) which returns a table of various bits of information on miRNAs, but I cannot adapt this command to just look at my list of miRNAs of interest (ie. the 'mirna' vector above). I've included the sessionInfo() output for these at the bottom of the email, but I suspect my problem is more to do with the fact I'm not going about this the right way (as opposed to a problem with package versions and coding etc.). I'm not even sure that using 'biomaRt' will give me the information I eventually want (the target mRNAs of these miRNAs), I was just trying it out, to see what it was capable of in terms of querying these databases. So I apologise for the vagueness. Since I haven't managed to get very far by myself then it's difficult to be more specific, but I'd really appreciate it if anyone could offer some advice, even just to point me in the direction of a useful package which might have gone unnoticed by me. Many thanks, Fiona Dr Fiona C Ingleby Postdoctoral Research Fellow University of Sussex Email: F.Ingleby at sussex.ac.uk Website: fionaingleby.weebly.com > sessionInfo() R version 2.15.2 (2012-10-26) Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit) locale: [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] biomaRt_2.14.0 affy_1.36.1 Biobase_2.18.0 BiocGenerics_0.4.0 loaded via a namespace (and not attached): [1] affyio_1.26.0 BiocInstaller_1.8.3 grid_2.15.2 lattice_0.20-14 Matrix_1.0-11 MCMCglmm_2.17 [7] preprocessCore_1.20.0 RCurl_1.95-4.1 tools_2.15.2 XML_3.95-0.2 zlibbioc_1.4.0 [[alternative HTML version deleted]] ------------------------------ Message: 11 Date: Thu, 28 Mar 2013 12:43:14 -0400 From: "James W. MacDonald" <jmacdon@uw.edu> To: Fiona Ingleby <fiona.ingleby at="" gmail.com=""> Cc: "bioconductor at r-project.org" <bioconductor at="" r-project.org=""> Subject: Re: [BioC] identifying drosophila miRNA targets Message-ID: <51547322.6020803 at uw.edu> Content-Type: text/plain; charset=windows-1252; format=flowed Hi Fiona, I have a function called mirna2mrna (yeah, I know, lame function name...) in my affycoretools package that does this, based on the sanger microcosm targets data that you can download here: http://www.ebi.ac.uk/enright-srv/microcosm/cgi- bin/targets/v5/download.pl there is also a function makeHmap() that will create a heatmap with the miRNA/mRNA pairs, where the color of the cells is based on the correlation between the two RNA species (with the intent to show negative correlations, indicating that the miRNA is hypothetically causing premature degradation of the mRNA). I think the help pages for these two functions are reasonable, but please let me know if you have any questions. Best, Jim On 3/28/2013 12:30 PM, Fiona Ingleby wrote: > Hi everyone, > > I am working with mRNA data from Affy 'drosophila2' arrays and miRNA data from Affy 'mirna3' arrays. I have identified a list of differentially expressed mRNAs and miRNAs. I'm having a bit of trouble with some downstream analyses and I'm hoping someone might be able to offer some help. > > I would like to use my list of differentially expressed miRNAs to access online databases (e.g. miRBase, microRNA.org?) and extract the names of all the potential target mRNAs. Then I'd like to use this list of mRNAs to look through my mRNA expression data. I'm aware of packages like 'RmiR' and 'microRNA' which have built-in functions for finding miRNA targets, but as far as I can tell, 'RmiR' uses miRNA databases for humans only and 'microRNA' works with human and mouse data only. So is there a package I am unaware of (or another application of 'RmiR'/'microRNA' that I am unaware of) for looking at drosophila data? > > So far I have also considered the 'biomaRt' package to see if the database query function on there can help me, but I haven't had much luck. For instance, if I try an example list of miRNAs: > > mirna<-c("dme-miR-1002","dme-miR-312","dme-miR-973") > library(biomaRt) > ensembl<-useMart("ensembl",dataset="dmelanogaster_gene_ensembl") > getBM(attributes="mirbase_accession",filters="mirbase_id",values=mir na,mart=ensembl) > > then 'logical(0)' is returned, as if there are no records for those miRNAs - but by searching the database manually I know the records are there. > > Alternatively I can try: > > miRNA<- getBM(c("mirbase_accession","mirbase_id", "ensembl_gene_id", "start_position", "chromosome_name"), filters = c("with_mirbase"), values = list(T), mart = ensembl) > > which returns a table of various bits of information on miRNAs, but I cannot adapt this command to just look at my list of miRNAs of interest (ie. the 'mirna' vector above). I've included the sessionInfo() output for these at the bottom of the email, but I suspect my problem is more to do with the fact I'm not going about this the right way (as opposed to a problem with package versions and coding etc.). I'm not even sure that using 'biomaRt' will give me the information I eventually want (the target mRNAs of these miRNAs), I was just trying it out, to see what it was capable of in terms of querying these databases. So I apologise for the vagueness. Since I haven't managed to get very far by myself then it's difficult to be more specific, but I'd really appreciate it if anyone could offer some advice, even just to point me in the direction of a useful package which might have gone unnoticed by me. > > Many thanks, > > Fiona > > Dr Fiona C Ingleby > Postdoctoral Research Fellow > University of Sussex > Email: F.Ingleby at sussex.ac.uk > Website: fionaingleby.weebly.com > > >> sessionInfo() > R version 2.15.2 (2012-10-26) > Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit) > > locale: > [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8 > > attached base packages: > [1] stats graphics grDevices utils datasets methods base > > other attached packages: > [1] biomaRt_2.14.0 affy_1.36.1 Biobase_2.18.0 BiocGenerics_0.4.0 > > loaded via a namespace (and not attached): > [1] affyio_1.26.0 BiocInstaller_1.8.3 grid_2.15.2 lattice_0.20-14 Matrix_1.0-11 MCMCglmm_2.17 > [7] preprocessCore_1.20.0 RCurl_1.95-4.1 tools_2.15.2 XML_3.95-0.2 zlibbioc_1.4.0 > [[alternative HTML version deleted]] > > > > _______________________________________________ > Bioconductor mailing list > Bioconductor at r-project.org > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor -- James W. MacDonald, M.S. Biostatistician University of Washington Environmental and Occupational Health Sciences 4225 Roosevelt Way NE, # 100 Seattle WA 98105-6099 ------------------------------ Message: 12 Date: Thu, 28 Mar 2013 16:58:21 +0000 From: "Blanchette, Marco" <mab@stowers.org> To: Kasper Daniel Hansen <kasperdanielhansen at="" gmail.com=""> Cc: "bioconductor at r-project.org" <bioconductor at="" r-project.org=""> Subject: Re: [BioC] Limit on number of sequence files for forging a BSgenome Message-ID: <1CE7C98DF956594C98E0956F4747178136D97D at MBSRV02.sgc.loc> Content-Type: text/plain; charset="us-ascii" Kasper, I see your line of thought, is there a particular fasta file causing forgeBSgenomeDataPkg() to break? The answer is no. Once I reach a certain number of fasta files, adding one more contig breaks the function. For instance, taking the first 454 contigs of C. brenneri breaks while removing the last or the first fasta file from the list (keeping only 453) compile without a problem (neither the last or the first fasta files are responsible for breaking the function, the number of file is the trigger) What's even more puzzling is that the number that breaks is not a fixed number. Selecting a random selection of contigs or changing genome will change the number that triggers the function to break... However it's always around 440 files, which might be due to the size of the fasta files being all of very similar sizes. Any clues? -- Marco Blanchette, Ph.D. Stowers Institute for Medical Research 1000 East 50th Street Kansas City MO 64110 www.stowers.org Tel: 816-926-4071 Cell: 816-726-8419 Fax: 816-926-2018 On 3/27/13 8:22 PM, "Kasper Daniel Hansen" <kasperdanielhansen at="" gmail.com=""> wrote: >Marco, > >You are probably right in diagnosing the problem, but sometimes I >think I have seen FASTA files with the entire sequence on a single >line, instead of (say) 80 nucleotides and then a newline. I could >believe that a really long contig on a single line without a newline, >could cause an error like this. You could quickly check if there is a >suspicious file by > wc -l * >and look for files with #lines like 2-3. Somehow 460 seems a weird >number to fail at. > >This may not be your problem, and I am sure Herve will respond in due >time. > >Best, >Kasper > >On Wed, Mar 27, 2013 at 4:28 PM, Blanchette, Marco <mab at="" stowers.org=""> >wrote: >> Hi, >> >> Is there a maximum number of sequence files (chromosomes or contigs in >>my case) that can be fed to the forgeBSgenomeDataPkg() function? I am >>trying to build a BSgenome for C. brenneri and C. japonica available >>from EnsemblGenomes. These genomes are made from thousands of contigs >>with genes annotated to them. Currently, I get the following error when >>running "Error: Line longer than buffer size" when running on the full >>set of contigs. However, it works fine on a seed file containing a >>subset of the contigs (I can forge a genome with 450 contigs but not >>with 460!) >> >> Any suggestions will be appreciated (I can provide a toy example but I >>am not sure what would be the merit of it at this point) >> >> Thanks >> >> -- Marco Blanchette, Ph.D. >> Stowers Institute for Medical Research >> 1000 East 50th Street >> Kansas City MO 64110 >> www.stowers.org >> >> Tel: 816-926-4071 >> Cell: 816-726-8419 >> Fax: 816-926-2018 >> >> [[alternative HTML version deleted]] >> >> _______________________________________________ >> Bioconductor mailing list >> Bioconductor at r-project.org >> https://stat.ethz.ch/mailman/listinfo/bioconductor >> Search the archives: >>http://news.gmane.org/gmane.science.biology.informatics.conductor ------------------------------ Message: 13 Date: Thu, 28 Mar 2013 19:07:19 +0000 From: "Shields, Rusty (IMS)" <shieldsr@imsweb.com> To: "bioconductor at r-project.org" <bioconductor at="" r-project.org=""> Subject: Re: [BioC] [R] Error in setMethod("combine"... was - Error when installing globaltest package Message-ID: <e2e839b5bfa711438b481ad9906b0c340a6a70a8 at="" varuna.omni.imsweb.com=""> Content-Type: text/plain; charset="iso-8859-1" Thanks Martin, No apologies necessary, just happy to have some help. Here's the info on Biobase: > packageDescription("Biobase") Package: Biobase Title: Biobase: Base functions for Bioconductor Version: 2.18.0 Author: R. Gentleman, V. Carey, M. Morgan, S. Falcon Description: Functions that are needed by many other packages or which replace R functions. Suggests: tools, tkWidgets, ALL Depends: R (>= 2.10), BiocGenerics (>= 0.3.2), utils Imports: methods, BiocGenerics Maintainer: Bioconductor Package Maintainer <maintainer at="" bioconductor.org=""> License: Artistic-2.0 Collate: tools.R strings.R environment.R vignettes.R packages.R AllGenerics.R VersionsClass.R VersionedClasses.R methods-VersionsNull.R methods-VersionedClass.R DataClasses.R methods-aggregator.R methods-container.R methods-MIAxE.R methods-MIAME.R methods-AssayData.R methods-AnnotatedDataFrame.R methods-eSet.R methods-ExpressionSet.R methods-MultiSet.R methods-SnpSet.R methods-NChannelSet.R anyMissing.R rowOp-methods.R updateObjectTo.R methods-ScalarObject.R zzz.R LazyLoad: yes biocViews: Infrastructure, Bioinformatics Packaged: 2012-10-02 02:41:50 UTC; biocbuild Built: R 2.14.0; x86_64-unknown-linux-gnu; 2013-03-22 15:01:20 UTC; unix -- File: /usr/local/R-2.14.0/lib64/R/library/Biobase/Meta/package.rds -----Original Message----- From: Martin Morgan [mailto:mtmorgan@fhcrc.org] Sent: Thursday, March 28, 2013 2:52 PM To: Shields, Rusty (IMS) Cc: r-help at r-project.org Subject: Re: [R] Error in setMethod("combine"... was - Error when installing globaltest package On 3/28/2013 11:05 AM, Shields, Rusty (IMS) wrote: > Hi All, > > I posted this on the bioconductor list and didn't get a response there, so I'm hoping someone here can help. > > I don't know a heck of a lot about R, so I apologize if this seems like a trivial issue. This error comes up when trying to install the bioconductor globaltest package. > Sorry that you didn't get a response on the Bioc mailing list; I'd actually suggest returning to the original thread there and I'll see that it gets answered. My guess is that you have a version of Biobase that is too new compared to the version (2.14.0) expected in the version of Bioconductor you are using. What does packageDescription("Biobase") say? Martin > Any clues? > > Thanks! > Rusty > > -----Original Message----- > From: bioconductor-bounces at r-project.org [mailto:bioconductor- bounces at r-project.org] On Behalf Of Shields, Rusty (IMS) > Sent: Tuesday, March 26, 2013 12:34 PM > To: bioconductor at r-project.org > Subject: [BioC] Error when installing globaltest package > > Hi all, > > I've run into a problem when attempting to install the globaltest package from Bioconductor. I'm using R 2.14.0 on 64bit SLES 11. Let me know what other information you might need about my system to troubleshoot this. > > Using the method described for this installation on the Bioconductor website: > > source("http://bioconductor.org/biocLite.R") > biocLite("globaltest") > > I get and the following result, which I can't find an reference to on the list archives: > >> biocLite("globaltest") > BioC_mirror: 'http://www.bioconductor.org' > Using R version 2.14, BiocInstaller version 1.2.1. > Installing package(s) 'globaltest' > trying URL 'http://www.bioconductor.org/packages/2.9/bioc/src/contri b/globaltest_5.8.1.tar.gz' > Content type 'application/x-gzip' length 956376 bytes (933 Kb) > opened URL > ================================================== > downloaded 933 Kb > > * installing *source* package ?globaltest? ... > ** R > ** data > ** inst > ** preparing package for lazy loading > Warning in .simpleDuplicateClass(def, prev) : > A specification for class ?data.frameOrNULL? > Creating a generic function for ?sort? from package ?base? in package ?globaltest? > Creating a generic function for ?model.matrix? from package ?stats? in package ?globaltest? > Creating a generic function for ?coefficients? from package ?stats? in package ?globaltest? > Creating a generic function for ?fitted.values? from package ?stats? in package ?globaltest? > Creating a generic function for ?residuals? from package ?stats? in package ?globaltest? > Error in setMethod("combine", signature(x = "gt.result", y = "gt.result"), : > no existing definition for function ?combine? > Error : unable to load R code in package ?globaltest? > ERROR: lazy loading failed for package ?globaltest? > * removing ?/usr/local/R-2.14.0/lib64/R/library/globaltest? > > The downloaded packages are in > ?/tmp/RtmpdCwjNB/downloaded_packages? > Updating HTML index of packages in '.Library' > Making packages.html ... done > Old packages: 'caret' > Update all/some/none? [a/s/n]: > > ________________________________ > > Information in this e-mail may be confidential. It is intended only for the addressee(s) identified above. If you are not the addressee(s), or an employee or agent of the addressee(s), please note that any dissemination, distribution, or copying of this communication is strictly prohibited. If you have received this e-mail in error, please notify the sender of the error. > > [[alternative HTML version deleted]] > > > ________________________________ > > Information in this e-mail may be confidential. It is intended only for the addressee(s) identified above. If you are not the addressee(s), or an employee or agent of the addressee(s), please note that any dissemination, distribution, or copying of this communication is strictly prohibited. If you have received this e-mail in error, please notify the sender of the error. > > ______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting- guide.html > and provide commented, minimal, self-contained, reproducible code. > -- Dr. Martin Morgan, PhD Fred Hutchinson Cancer Research Center 1100 Fairview Ave. N. PO Box 19024 Seattle, WA 98109 ________________________________ Information in this e-mail may be confidential. It is intended only for the addressee(s) identified above. If you are not the addressee(s), or an employee or agent of the addressee(s), please note that any dissemination, distribution, or copying of this communication is strictly prohibited. If you have received this e-mail in error, please notify the sender of the error. ------------------------------ Message: 14 Date: Thu, 28 Mar 2013 12:26:11 -0700 From: Martin Morgan <mtmorgan@fhcrc.org> To: "Shields, Rusty (IMS)" <shieldsr at="" imsweb.com=""> Cc: "bioconductor at r-project.org" <bioconductor at="" r-project.org=""> Subject: Re: [BioC] [R] Error in setMethod("combine"... was - Error when installing globaltest package Message-ID: <51549953.6080003 at fhcrc.org> Content-Type: text/plain; charset=ISO-8859-1; format=flowed On 03/28/2013 12:07 PM, Shields, Rusty (IMS) wrote: > Thanks Martin, > > No apologies necessary, just happy to have some help. > > Here's the info on Biobase: > >> packageDescription("Biobase") > Package: Biobase > Title: Biobase: Base functions for Bioconductor > Version: 2.18.0 yes, this is the version of Biobase from the current version of R / Bioconductor, whereas you're using a relatively old R / Biocondcutor (R 2.14, Bioc 2.9, from the 'Previous versions' box at http://bioconductor.org/help/) where Biobase is at version 2.14. This doesn't bode well, as other packages are also likely at incorrect versions. The easiest solution is to upgrade to a recent R, and start again with biocLite('globaltest'). You could also try biocLite("Biobase") to get the correct version, then biocLite("globaltest"), but if other packages are also installed incorrectly... Something like source("http://bioconductor.org/biocLite.R") avail = available.packages(contriburl=contrib.url(biocinstallRepos())) inst = installed.packages(priority="NA") idx = rownames(avail) %in% rownames(inst) vers = avail[idx, "Version"] vers[vers < inst[names(vers), "Version"]] will give you the too-new packages that need to be re-installed. It's hard to know how your system got to its current state; one candidate is that .libPaths() includes a path shared by several versions of R. but using biocLite() to install packages should help avoiding that in the future. Martin > Author: R. Gentleman, V. Carey, M. Morgan, S. Falcon > Description: Functions that are needed by many other packages or which > replace R functions. > Suggests: tools, tkWidgets, ALL > Depends: R (>= 2.10), BiocGenerics (>= 0.3.2), utils > Imports: methods, BiocGenerics > Maintainer: Bioconductor Package Maintainer > <maintainer at="" bioconductor.org=""> > License: Artistic-2.0 > Collate: tools.R strings.R environment.R vignettes.R packages.R > AllGenerics.R VersionsClass.R VersionedClasses.R > methods-VersionsNull.R methods-VersionedClass.R DataClasses.R > methods-aggregator.R methods-container.R methods-MIAxE.R > methods-MIAME.R methods-AssayData.R > methods-AnnotatedDataFrame.R methods-eSet.R > methods-ExpressionSet.R methods-MultiSet.R methods-SnpSet.R > methods-NChannelSet.R anyMissing.R rowOp-methods.R > updateObjectTo.R methods-ScalarObject.R zzz.R > LazyLoad: yes > biocViews: Infrastructure, Bioinformatics > Packaged: 2012-10-02 02:41:50 UTC; biocbuild > Built: R 2.14.0; x86_64-unknown-linux-gnu; 2013-03-22 15:01:20 UTC; > unix > > -- File: /usr/local/R-2.14.0/lib64/R/library/Biobase/Meta/package.rds > > -----Original Message----- > From: Martin Morgan [mailto:mtmorgan at fhcrc.org] > Sent: Thursday, March 28, 2013 2:52 PM > To: Shields, Rusty (IMS) > Cc: r-help at r-project.org > Subject: Re: [R] Error in setMethod("combine"... was - Error when installing globaltest package > > On 3/28/2013 11:05 AM, Shields, Rusty (IMS) wrote: >> Hi All, >> >> I posted this on the bioconductor list and didn't get a response there, so I'm hoping someone here can help. >> >> I don't know a heck of a lot about R, so I apologize if this seems like a trivial issue. This error comes up when trying to install the bioconductor globaltest package. >> > > Sorry that you didn't get a response on the Bioc mailing list; I'd actually > suggest returning to the original thread there and I'll see that it gets answered. > > My guess is that you have a version of Biobase that is too new compared to the > version (2.14.0) expected in the version of Bioconductor you are using. What > does packageDescription("Biobase") say? > > Martin > > > > >> Any clues? >> >> Thanks! >> Rusty >> >> -----Original Message----- >> From: bioconductor-bounces at r-project.org [mailto:bioconductor- bounces at r-project.org] On Behalf Of Shields, Rusty (IMS) >> Sent: Tuesday, March 26, 2013 12:34 PM >> To: bioconductor at r-project.org >> Subject: [BioC] Error when installing globaltest package >> >> Hi all, >> >> I've run into a problem when attempting to install the globaltest package from Bioconductor. I'm using R 2.14.0 on 64bit SLES 11. Let me know what other information you might need about my system to troubleshoot this. >> >> Using the method described for this installation on the Bioconductor website: >> >> source("http://bioconductor.org/biocLite.R") >> biocLite("globaltest") >> >> I get and the following result, which I can't find an reference to on the list archives: >> >>> biocLite("globaltest") >> BioC_mirror: 'http://www.bioconductor.org' >> Using R version 2.14, BiocInstaller version 1.2.1. >> Installing package(s) 'globaltest' >> trying URL 'http://www.bioconductor.org/packages/2.9/bioc/src/contr ib/globaltest_5.8.1.tar.gz' >> Content type 'application/x-gzip' length 956376 bytes (933 Kb) >> opened URL >> ================================================== >> downloaded 933 Kb >> >> * installing *source* package ?globaltest? ... >> ** R >> ** data >> ** inst >> ** preparing package for lazy loading >> Warning in .simpleDuplicateClass(def, prev) : >> A specification for class ?data.frameOrNULL? >> Creating a generic function for ?sort? from package ?base? in package ?globaltest? >> Creating a generic function for ?model.matrix? from package ?stats? in package ?globaltest? >> Creating a generic function for ?coefficients? from package ?stats? in package ?globaltest? >> Creating a generic function for ?fitted.values? from package ?stats? in package ?globaltest? >> Creating a generic function for ?residuals? from package ?stats? in package ?globaltest? >> Error in setMethod("combine", signature(x = "gt.result", y = "gt.result"), : >> no existing definition for function ?combine? >> Error : unable to load R code in package ?globaltest? >> ERROR: lazy loading failed for package ?globaltest? >> * removing ?/usr/local/R-2.14.0/lib64/R/library/globaltest? >> >> The downloaded packages are in >> ?/tmp/RtmpdCwjNB/downloaded_packages? >> Updating HTML index of packages in '.Library' >> Making packages.html ... done >> Old packages: 'caret' >> Update all/some/none? [a/s/n]: >> >> ________________________________ >> >> Information in this e-mail may be confidential. It is intended only for the addressee(s) identified above. If you are not the addressee(s), or an employee or agent of the addressee(s), please note that any dissemination, distribution, or copying of this communication is strictly prohibited. If you have received this e-mail in error, please notify the sender of the error. >> >> [[alternative HTML version deleted]] >> >> >> ________________________________ >> >> Information in this e-mail may be confidential. It is intended only for the addressee(s) identified above. If you are not the addressee(s), or an employee or agent of the addressee(s), please note that any dissemination, distribution, or copying of this communication is strictly prohibited. If you have received this e-mail in error, please notify the sender of the error. >> >> ______________________________________________ >> R-help at r-project.org mailing list >> https://stat.ethz.ch/mailman/listinfo/r-help >> PLEASE do read the posting guide http://www.R-project.org/posting- guide.html >> and provide commented, minimal, self-contained, reproducible code. >> > > > -- > Dr. Martin Morgan, PhD > Fred Hutchinson Cancer Research Center > 1100 Fairview Ave. N. > PO Box 19024 Seattle, WA 98109 > > ________________________________ > > Information in this e-mail may be confidential. It is intended only for the addressee(s) identified above. If you are not the addressee(s), or an employee or agent of the addressee(s), please note that any dissemination, distribution, or copying of this communication is strictly prohibited. If you have received this e-mail in error, please notify the sender of the error. > > _______________________________________________ > Bioconductor mailing list > Bioconductor at r-project.org > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor > -- Computational Biology / Fred Hutchinson Cancer Research Center 1100 Fairview Ave. N. PO Box 19024 Seattle, WA 98109 Location: Arnold Building M1 B861 Phone: (206) 667-2793 ------------------------------ Message: 15 Date: Thu, 28 Mar 2013 20:00:24 +0000 From: "Shields, Rusty (IMS)" <shieldsr@imsweb.com> To: Martin Morgan <mtmorgan at="" fhcrc.org=""> Cc: "bioconductor at r-project.org" <bioconductor at="" r-project.org=""> Subject: Re: [BioC] [R] Error in setMethod("combine"... was - Error when installing globaltest package Message-ID: <e2e839b5bfa711438b481ad9906b0c340a6a78ad at="" varuna.omni.imsweb.com=""> Content-Type: text/plain; charset="iso-8859-1" Correcting the version of Biobase seems to have fixed it, and after that it doesn't appear that anything else is of the wrong version - based on the code that you provided below. globaltest now installs successfully through BiocLite. As for how the wrong version of Biobase was installed, I have to take the blame for that. Rather than following the instructions on the bioconductor website to install using BiocLite, I made the mistake of downloading the package source manually and installing via 'R CMD INSTALL ...'. I believe I used the biocLite method to install all of the other packages. It has been a very educational few days. Thanks again. -----Original Message----- From: Martin Morgan [mailto:mtmorgan@fhcrc.org] Sent: Thursday, March 28, 2013 3:26 PM To: Shields, Rusty (IMS) Cc: bioconductor at r-project.org Subject: Re: [BioC] [R] Error in setMethod("combine"... was - Error when installing globaltest package On 03/28/2013 12:07 PM, Shields, Rusty (IMS) wrote: > Thanks Martin, > > No apologies necessary, just happy to have some help. > > Here's the info on Biobase: > >> packageDescription("Biobase") > Package: Biobase > Title: Biobase: Base functions for Bioconductor > Version: 2.18.0 yes, this is the version of Biobase from the current version of R / Bioconductor, whereas you're using a relatively old R / Biocondcutor (R 2.14, Bioc 2.9, from the 'Previous versions' box at http://bioconductor.org/help/) where Biobase is at version 2.14. This doesn't bode well, as other packages are also likely at incorrect versions. The easiest solution is to upgrade to a recent R, and start again with biocLite('globaltest'). You could also try biocLite("Biobase") to get the correct version, then biocLite("globaltest"), but if other packages are also installed incorrectly... Something like source("http://bioconductor.org/biocLite.R") avail = available.packages(contriburl=contrib.url(biocinstallRepos())) inst = installed.packages(priority="NA") idx = rownames(avail) %in% rownames(inst) vers = avail[idx, "Version"] vers[vers < inst[names(vers), "Version"]] will give you the too-new packages that need to be re-installed. It's hard to know how your system got to its current state; one candidate is that .libPaths() includes a path shared by several versions of R. but using biocLite() to install packages should help avoiding that in the future. Martin > Author: R. Gentleman, V. Carey, M. Morgan, S. Falcon > Description: Functions that are needed by many other packages or which > replace R functions. > Suggests: tools, tkWidgets, ALL > Depends: R (>= 2.10), BiocGenerics (>= 0.3.2), utils > Imports: methods, BiocGenerics > Maintainer: Bioconductor Package Maintainer > <maintainer at="" bioconductor.org=""> > License: Artistic-2.0 > Collate: tools.R strings.R environment.R vignettes.R packages.R > AllGenerics.R VersionsClass.R VersionedClasses.R > methods-VersionsNull.R methods-VersionedClass.R DataClasses.R > methods-aggregator.R methods-container.R methods-MIAxE.R > methods-MIAME.R methods-AssayData.R > methods-AnnotatedDataFrame.R methods-eSet.R > methods-ExpressionSet.R methods-MultiSet.R methods-SnpSet.R > methods-NChannelSet.R anyMissing.R rowOp-methods.R > updateObjectTo.R methods-ScalarObject.R zzz.R > LazyLoad: yes > biocViews: Infrastructure, Bioinformatics > Packaged: 2012-10-02 02:41:50 UTC; biocbuild > Built: R 2.14.0; x86_64-unknown-linux-gnu; 2013-03-22 15:01:20 UTC; > unix > > -- File: /usr/local/R-2.14.0/lib64/R/library/Biobase/Meta/package.rds > > -----Original Message----- > From: Martin Morgan [mailto:mtmorgan at fhcrc.org] > Sent: Thursday, March 28, 2013 2:52 PM > To: Shields, Rusty (IMS) > Cc: r-help at r-project.org > Subject: Re: [R] Error in setMethod("combine"... was - Error when installing globaltest package > > On 3/28/2013 11:05 AM, Shields, Rusty (IMS) wrote: >> Hi All, >> >> I posted this on the bioconductor list and didn't get a response there, so I'm hoping someone here can help. >> >> I don't know a heck of a lot about R, so I apologize if this seems like a trivial issue. This error comes up when trying to install the bioconductor globaltest package. >> > > Sorry that you didn't get a response on the Bioc mailing list; I'd actually > suggest returning to the original thread there and I'll see that it gets answered. > > My guess is that you have a version of Biobase that is too new compared to the > version (2.14.0) expected in the version of Bioconductor you are using. What > does packageDescription("Biobase") say? > > Martin > > > > >> Any clues? >> >> Thanks! >> Rusty >> >> -----Original Message----- >> From: bioconductor-bounces at r-project.org [mailto:bioconductor- bounces at r-project.org] On Behalf Of Shields, Rusty (IMS) >> Sent: Tuesday, March 26, 2013 12:34 PM >> To: bioconductor at r-project.org >> Subject: [BioC] Error when installing globaltest package >> >> Hi all, >> >> I've run into a problem when attempting to install the globaltest package from Bioconductor. I'm using R 2.14.0 on 64bit SLES 11. Let me know what other information you might need about my system to troubleshoot this. >> >> Using the method described for this installation on the Bioconductor website: >> >> source("http://bioconductor.org/biocLite.R") >> biocLite("globaltest") >> >> I get and the following result, which I can't find an reference to on the list archives: >> >>> biocLite("globaltest") >> BioC_mirror: 'http://www.bioconductor.org' >> Using R version 2.14, BiocInstaller version 1.2.1. >> Installing package(s) 'globaltest' >> trying URL 'http://www.bioconductor.org/packages/2.9/bioc/src/contr ib/globaltest_5.8.1.tar.gz' >> Content type 'application/x-gzip' length 956376 bytes (933 Kb) >> opened URL >> ================================================== >> downloaded 933 Kb >> >> * installing *source* package ?globaltest? ... >> ** R >> ** data >> ** inst >> ** preparing package for lazy loading >> Warning in .simpleDuplicateClass(def, prev) : >> A specification for class ?data.frameOrNULL? >> Creating a generic function for ?sort? from package ?base? in package ?globaltest? >> Creating a generic function for ?model.matrix? from package ?stats? in package ?globaltest? >> Creating a generic function for ?coefficients? from package ?stats? in package ?globaltest? >> Creating a generic function for ?fitted.values? from package ?stats? in package ?globaltest? >> Creating a generic function for ?residuals? from package ?stats? in package ?globaltest? >> Error in setMethod("combine", signature(x = "gt.result", y = "gt.result"), : >> no existing definition for function ?combine? >> Error : unable to load R code in package ?globaltest? >> ERROR: lazy loading failed for package ?globaltest? >> * removing ?/usr/local/R-2.14.0/lib64/R/library/globaltest? >> >> The downloaded packages are in >> ?/tmp/RtmpdCwjNB/downloaded_packages? >> Updating HTML index of packages in '.Library' >> Making packages.html ... done >> Old packages: 'caret' >> Update all/some/none? [a/s/n]: >> >> ________________________________ >> >> Information in this e-mail may be confidential. It is intended only for the addressee(s) identified above. If you are not the addressee(s), or an employee or agent of the addressee(s), please note that any dissemination, distribution, or copying of this communication is strictly prohibited. If you have received this e-mail in error, please notify the sender of the error. >> >> [[alternative HTML version deleted]] >> >> >> ________________________________ >> >> Information in this e-mail may be confidential. It is intended only for the addressee(s) identified above. If you are not the addressee(s), or an employee or agent of the addressee(s), please note that any dissemination, distribution, or copying of this communication is strictly prohibited. If you have received this e-mail in error, please notify the sender of the error. >> >> ______________________________________________ >> R-help at r-project.org mailing list >> https://stat.ethz.ch/mailman/listinfo/r-help >> PLEASE do read the posting guide http://www.R-project.org/posting- guide.html >> and provide commented, minimal, self-contained, reproducible code. >> > > > -- > Dr. Martin Morgan, PhD > Fred Hutchinson Cancer Research Center > 1100 Fairview Ave. N. > PO Box 19024 Seattle, WA 98109 > > ________________________________ > > Information in this e-mail may be confidential. It is intended only for the addressee(s) identified above. If you are not the addressee(s), or an employee or agent of the addressee(s), please note that any dissemination, distribution, or copying of this communication is strictly prohibited. If you have received this e-mail in error, please notify the sender of the error. > > _______________________________________________ > Bioconductor mailing list > Bioconductor at r-project.org > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor > -- Computational Biology / Fred Hutchinson Cancer Research Center 1100 Fairview Ave. N. PO Box 19024 Seattle, WA 98109 Location: Arnold Building M1 B861 Phone: (206) 667-2793 ________________________________ Information in this e-mail may be confidential. It is intended only for the addressee(s) identified above. If you are not the addressee(s), or an employee or agent of the addressee(s), please note that any dissemination, distribution, or copying of this communication is strictly prohibited. If you have received this e-mail in error, please notify the sender of the error. ------------------------------ Message: 16 Date: Thu, 28 Mar 2013 22:19:57 +0100 From: Michael Love <michaelisaiahlove@gmail.com> To: Michael Muratet <mmuratet at="" hudsonalpha.org=""> Cc: bioconductor at r-project.org Subject: Re: [BioC] Interpreting DESeq2 results Message-ID: <cadqzidv-qp+=nm_goxj2aso+2maevq_i_y0mfgu1vv3o8fzy1g at="" mail.gmail.com=""> Content-Type: text/plain Hi Michael, The baseMean column is not on the log scale; it is the mean of normalized counts for a gene. The intercept from the GLM is labelled intercept in mcols(dse). Mike On Mar 28, 2013 5:00 PM, "Michael Muratet" <mmuratet at="" hudsonalpha.org=""> wrote: > Greetings > > I have an experiment: > > > design(dse) > ~ factor1 + factor2 + factor3 > > where factor1 has two levels, factor2 has three levels and factor3 has > three levels. I extract a gene of interest from the results for each term > (I've changed the indices to reflect the condition): > > > lapply(resultsNames(dse),function(u) results(dse,u)["gene_A",]) > [["Intercept"]] > baseMean log2FoldChange pvalue FDR > gene_A 1596.548 10.77485 3.309439e-216 7.025442e-216 > [["factor1_level2"]] > baseMean log2FoldChange pvalue FDR > gene_A 1596.548 0.3386776 0.1307309 0.3587438 > [["factor2_level2"]] > baseMean log2FoldChange pvalue FDR > gene_A 1596.548 -0.6882543 0.0613569 0.1007896 > [["factor2_level3"]] > baseMean log2FoldChange pvalue FDR > gene_A 1596.548 0.2393368 0.513216 0.6589575 > [["factor3_level2"]] > baseMean log2FoldChange pvalue FDR > gene_A 1596.548 0.1584153 0.6423634 0.8503163 > [["factor3_level3]] > baseMean log2FoldChange pvalue FDR > gene_A 1596.548 -1.627898 1.823141e-06 0.001409384 > > I want to be sure I understand the output format. Is it true that the > coefficients (the vector beta) from the fit are the baseMean value scaled > by the log2FoldChange? Is the true intercept value > 1596.548*2^10.77485=2797274.13? > > mcols() tells me that the baseMean term is calculated over "all rows". The > baseMean is different for different genes although it is the same for each > gene across all the conditions, I'm not seeing how the rows are selected. > > Thanks > > Mike > > Michael Muratet, Ph.D. > Senior Scientist > HudsonAlpha Institute for Biotechnology > mmuratet at hudsonalpha.org > (256) 327-0473 (p) > (256) 327-0966 (f) > > Room 4005 > 601 Genome Way > Huntsville, Alabama 35806 > > _______________________________________________ > Bioconductor mailing list > Bioconductor at r-project.org > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: > http://news.gmane.org/gmane.science.biology.informatics.conductor > [[alternative HTML version deleted]] ------------------------------ Message: 17 Date: Thu, 28 Mar 2013 17:50:04 -0400 From: somnath bandyopadhyay <genome1976@hotmail.com> Cc: "bioconductor at r-project.org" <bioconductor at="" r-project.org=""> Subject: Re: [BioC] Interpreting DESeq2 results Message-ID: <bay175-w24d794c2988d81c509ef9fcdd20 at="" phx.gbl=""> Content-Type: text/plain Will the LIMMA package work on sparse matrices? Thanks,Som. [[alternative HTML version deleted]] ------------------------------ _______________________________________________ Bioconductor mailing list Bioconductor at r-project.org https://stat.ethz.ch/mailman/listinfo/bioconductor End of Bioconductor Digest, Vol 121, Issue 29