Important questions about creating specific plots for microarray datasets with ggplot2

[svlac...@eie.gr]

Dear ggplot2 group,

Firstly, i would like to wish to everybody  Merry Christmas and a happy new year !!  My name is Efstathios-Iason Vlachavas, I'm from Greece and I'm currently working as a PhD candidate at National Hellenic Research Foundation under the team of Metabolic Engineering & Bioinformatics, in Greece (http://www.eie.gr/nhrf/institutes/ibrb/programmes/metabolicengineering-en.html). I'm using R/Bioconductor about one year, but i recently discovered and have a first look in ggplot2 capabilities and documentation !! In detail, because im my current project im dealing mainly with microarray technology, & besides my main analysis, i'm also trying to build some pipelines regarding exploratory data analysis. Thus, i want to produce more elegant and higher quality graphics and plots regarding both my project analysis and for my pipelines--so,

As i have searched in the updated ggplot2 cheatsheet, 

http://www.rstudio.com/wp-content/uploads/2015/12/ggplot2-cheatsheet-2.0.pdf,

and also in other pipelines(i.e. RNA-seq analysis from Bioconductor http://www.bioconductor.org/help/workflows/rnaseqGene/), 

i have mainly experienced problems about creating some specific plots of interest:

In detail, a first "naive" function i created, was to deal with the output of function topTable from limma, and create a volcano plot, like below:

volcano.r <- function(object, FC_cutoff_1, FC_cutoff_2, pval_cutoff) { # where object is the output of topTable

  

  results <- object

  results$threshold <- "1"

  a <- subset(results, adj.P.Val < pval_cutoff & logFC > FC_cutoff_1)

  results[rownames(a), "threshold"] <- "2"

  b <- subset(results, adj.P.Val < pval_cutoff & logFC < FC_cutoff_2 )

  results[rownames(b), "threshold"] <- "3"

  

  volcano <- ggplot(data=results, aes(x=logFC, y=-1*log10(adj.P.Val), colour=threshold,      label=SYMBOL))

  

  volcano <- volcano + geom_point() + scale_color_manual(values=c("black","red","green"),   labels=c("Not Significant", "Upregulated", "Downregulated"), name="key/Legend") 

  

  volcano2 <- volcano + theme(legend.text=element_text(size=16))

  print(volcano2)

  ggsave(filename="eleos.jpeg", plot=volcano2)

}

(im posting the code in order to illustrate my main question below)

In this specific case, here i can distinguish the aesthetics regarding the data frame, both for the x and the y axis. However, my main difficulty was to try to create some general boxplots, histograms and Q-Q plots regarding for instance a specific microarray dataset in R, as in this case i don't deal with data frames, but mainly with ExpressionSets or matrices. Thus, i could not understand how to appropriately define the aesthetics, because here all the rows are the genes, and the columns represent the different samples. So,

1) Is there a tutorial or book regarding high-throughput data analysis with ggplot2, especially with examples for microarray data ?

2) Or is it possible to provide me with a small idea or a brief explanation on this matter, or where i could look to deal with my issue, as my main goal is to learn how to create these specific plots for microarray expression sets(boxplots, histograms and quantile-quantile plots) ??

Thank you in advance for your time and patience, and i hope i don't disturb anyone with any naive questions--but any suggestions or help on this matter would be beneficial for a newby like me !!

Best Regards,

Efstathios

[Brandon Hurr]

Efstathios, 

It's been a while, but I have done some exploratory visualizations of microarray and RNAseq data using ggplot2. To answer your questions:

1) to my knowledge such a thing does not exist, but there might be some things out there in personal scripts and gists.

2) It sounds like you just need to transpose your array. If you have a matrix as you describe just wrap it in as.data.frame(t("matrix goes here")). My understanding is that you'll want some sort of normalization of your matrix before plotting. 

If you have more specific questions related to plotting, could you dput() a subset of your dataset and provide a better description of what you want?

B

--
--
You received this message because you are subscribed to the ggplot2 mailing list.
Please provide a reproducible example: https://github.com/hadley/devtools/wiki/Reproducibility
 
To post: email ggp...@googlegroups.com
To unsubscribe: email ggplot2+u...@googlegroups.com
More options: http://groups.google.com/group/ggplot2

---
You received this message because you are subscribed to the Google Groups "ggplot2" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ggplot2+u...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[svlac...@eie.gr]

Dear Brandon, thank you for your answer, and please excuse me for the late response, but the time zone here in Greece is different. So, i will create a "fake" expression set, like an affymetrix dataset i have in R, in order to post more appropriately my questions :

> set.seed(81)

> mat <- matrix(data=sample(1:100), nrow=12136, ncol=60)

pdat <- data.frame(Condition=rep(c("Normal", "Diseased"),times=30))

new.set <- new("ExpressionSet", exprs=mat)

 pData(new.set) <-pdat

> head(pData(new.set))

  Condition

1    Normal

2  Diseased

3    Normal

4  Diseased

5    Normal

6  Diseased

> dim(new.set)

Features  Samples 

   12136       60

Thus, lets suppose that i want to produce some exploratory plots for this "aritficial" microarray dataset--in detail, histogram and boxplot plots. Now, i understand very naively, that in order to use the ggplot() function, i have to use somehow:

for istance, ggplot(data, aes(x=, y=) + geom_histogram()----where data is a data frame, and x represents the log2 intensities(assume we have these above) and y the sample names. But the main problem, is when i try to create a data frame to include these two columns. the samples are 60, whereas the features-genes are 12.136 !! So, my main problem is how to move or which transformations do i have to perform in order to apply the same logic, to create the histogram above or a boxplot  ??

Best,

Efstathios

thank you for your answers 

[svlac...@eie.gr]

Just a small correction----in the possible histogram with ggplot2, maybe the column names aren't necessary, as possibly there are too many, so perhaps the code above for ggplot2 is more appropriate for the boxplot implementation !!

[Brandon Hurr]

What packages are you using? 

set.seed(81)

mat <- matrix(data=sample(1:100), nrow=12136, ncol=60)

#Warning message:

#In matrix(data = sample(1:100), nrow = 12136, ncol = 60) :

  data length [100] is not a sub-multiple or multiple of the number of rows [12136]

new.set <- new("ExpressionSet", exprs=mat)

#Error in getClass(Class, where = topenv(parent.frame())) : 

  “ExpressionSet” is not a defined class

pData(new.set) <-pdat

#Error: object 'pdat' not found

[svlac...@eie.gr]

the R package affy

[Ista Zahn]

I'm not sure I've understood this correctly, but probably you want
something like

library(tidyr)
x <- as.data.frame(new.set)
x <- as.data.frame(new.set)
x$sample <- factor(1:nrow(x))
x <- gather(x, feature, intensity, -Condition, -sample)
ggplot(x, aes(x = sample, y = intensity, fill = Condition)) +
geom_boxplot() + coord_flip()

Best,
Ista

[svlac...@eie.gr]

Dear Ista,

thank you for your answer and code provided !! i have never used tidyr in the past, so i will search it---just to pinpoint some important parts----

inside the gather function, with the - symbol you remove these variables ?

Also the same logic applies to the creation of a histogram plot ??

Best,

Efstathios

[Ista Zahn]

Hi Efstathios,

You don't have to use the tidyr package, but you do need to reshape
the data before plotting. You could do that with stats::reshape, with
reshape2::melt, data.table::melt, or tidyr::gather or any other means
you prefer. For details on the arguments to tidyr::gather see ?gather.

Histograms will need somewhat different treatment, perhaps

ggplot(x, aes(x = intensity, fill = Condition)) + geom_histogram() +
facet_wrap(~sample)

Best,
Ista

[svlac...@eie.gr]

Dear Ista,

thank you again for your answer---i will search your instructions, as i understand that anyway it is necessary to perform the above kind of transformation---anyway,

i used your code, but i removed from the boxplot the function coord_flip in order for the boxplot to look right :

[Brandon Hurr]

Yes, Ista's approach is better than what I was coming up with for summarizing across all samples. 

My approach was to extract the phenotype and matrix from the expressionSet object, transpose it, filter it for each gene and then plot it.  

https://gist.github.com/bhive01/aa15ac57546bc64b2f08

Then I ran into a plotting issue. 

Worked with mtcars. I don't get it. 

Anyway, this is helpful for checking the overall signal differences between samples, but I think what you're after is a comparison of differences between Normal and Diseased for each gene, but this would be impossible to read and take ages to plot. You need to filter your dataset for things that are different by some statistical test (maybe an output of the affy analysis) and then plot those that have interesting differences. I never had enough samples to really differentiate treatments, but if you have 30 and 30 you should be able to easily test for differences. 

[adam.l...@pnc.com]

Is tidyr not available on CRAN? I do not find it.




Adam

The contents of this email are the property of PNC. If it was not addressed to you, you have no legal right to read it. If you think you received it in error, please notify the sender. Do not forward or copy without permission of the sender. This message may be considered a commercial electronic message under Canadian law or this message may contain an advertisement of a product or service and thus may constitute a commercial electronic mail message under US law. You may unsubscribe at any time from receiving commercial electronic messages from PNC at http://pages.e.pnc.com/globalunsub/
PNC, 249 Fifth Avenue, Pittsburgh, PA 15222; pnc.com

[Ista Zahn]

https://cran.rstudio.com/web/packages/tidyr/index.html

[svlac...@eie.gr]

Dear Brandon,

 thank you again for your detailed reply and your consideration on this matter !! basically, my main goal of the 2 mentioned plots, as with others--is basically to have a brief view after normalization and non-specific intensity filtering, about any "unusual" abnormalities, in any sample, or to detect any "weird" distributions, than to inspect for differences between studied groups---although your idea is nice---but i believe that in my case, maybe the group variable could be used as a color indicator for instance !!

Best, 

Efstathios

[svlac...@eie.gr]

Dear Ista,

one last question---

regarding the histogram creation, can i somehow include a different color for each different sample ? or here it can't be visualized ?

Moreover, can also include a "trendline" to the histogram ? maybe with + geom_density() ? or i have to change also more arguments ?

Best,

Efstathios

[Ista Zahn]

On Tue, Jan 5, 2016 at 2:50 PM, <svlac...@eie.gr> wrote:
> Dear Ista,
>

> one last question---
>
> regarding the histogram creation, can i somehow include a different color
> for each different sample ? or here it can't be visualized ?

There is no technical reason you cannot do it, but you won't be able
to distinguish 60 different colors so there isn't much point.

> Moreover, can also include a "trendline" to the histogram ? maybe with +
> geom_density() ? or i have to change also more arguments ?

Try it and see.

Best,
Ista

[svlac...@eie.gr]

Dear Ista,

please excuse me to return with one more question, but i stuck into an issue:

in detail, when i used your above code with one of my actual affymetrix microarray datasets---the main difference is that in my pheno Data i have more variables than my above artificial example, such as

> head(pData(eset.4))

              Disease      Location Meta_factor       Study

St_1_WL57.CEL  Normal sigmoid_colon           0 hgu133plus2

St_2_WL57.CEL  Cancer sigmoid_colon           0 hgu133plus2

St_N_EC59.CEL  Normal sigmoid_colon           0 hgu133plus2

St_T_EC59.CEL  Cancer sigmoid_colon           0 hgu133plus2

St_N_EJ58.CEL  Normal         cecum           0 hgu133plus2

St_T_EJ58.CEL  Cancer         cecum           0 hgu133plus2

So the variable "Disease" is similar to the condition above--but when i used the code as above:

x <- as.data.frame(eset.4) # eset.4 my actual expression set

x$sample <- factor(1:nrow(x))

Then first i got this warning:

> x <- gather(x, feature, intensity, -Disease, -sample)

Warning message:

attributes are not identical across measure variables; they will be dropped 

and then when i typed:

> ggplot(x, aes(x=sample, y=intensity, fill=Disease)) + geom_boxplot()

i got the following very weird image !!

Any ideas on this matter ??

Thank you,

Efstathios

[Ista Zahn]

Presumably you need to exclude Location, Meta_factor, and Study from the gather call.

Best,
Ista

[svlac...@eie.gr]

Dear Ista,

thank you for your indication. I played a bit with gather function, and when i removed all the mentioned labels it worked. Maybe from the beggining, to avoid this i could create first a matrix of only the expression values, and then coerce to a data.frame in order to avoid further arguments. But one last thing that struggles me a bit : is it possible to visualize also sample names, instead of just the number of the samples ? I could not find any way to solve this !!

Best,

Efstathios