Problem in running R-code

Question

[sinha.puja]

Hi All,

I wrote a code to plot a matrix graph but after running the code, plot doesn't match with the matrix data. I don't able to figure out the problem. If anyone knows how to resolve it then it would be highly appreciated. R-code is:

# install.packages("ggcorrplot")
# install.packages("ggplot2")

library(ggplot2)
library(ggcorrplot)

data <- read.csv(file="C:/Users/sinha.puja/Desktop/ChromHMM/practice.csv", sep=",")
rnames <- data[,1]                                       
mat_data <- data.matrix(data[,2:ncol(data)])              
rownames(mat_data) <- rnames                             
mat_data <- expand.grid(mark = c("H3K27me3", "H3K9me3", "H3K36me3"), 
                                 state = c("TssP", "PcRepr", "Heter/low", "EnhP/low", "ConHeter"))
mat_data$Probability <- rnorm(nrow(mat_data))

ggplot(data = mat_data, aes(x = mark, y = state)) +
           geom_tile(aes(fill = Probability), colour="grey80") + 
           scale_fill_gradient(low="white", high="red")

The input dataset is as follows:

state (Emission order)  H3K27me3    H3K9me3 H3K36me3
TssP    0.525472473 0.013872603 0.012417431
PcRepr  0.028550739 0.004967125 0.000949458
Heter/low   0.001659524 0.002706355 0.000396299
EnhP/low    0.004981548 0.183992447 0.002332633
ConHeter    0.118632661 0.82344082  0.395105767

Please help me in plotting the matrix without changing the input values.

[Friederike]

I don't understand what you mean by "plot doesn't match matrix code".

It would probably be helpful if you:

• shared the image (you can check this post here on How to add images)
• showed the content of mat_data (e.g. head(mat_data))
• pointed out the discrepancy

[Friederike]

Hi, please make use of the formatting buttons that are available. I've edited your post to properly show the code portions; in the future you can do that yourself by using the button next to the quotation marks.

[sinha.puja]

Hi, Thank you for your effort. Below is the image published from the article, similar one I want to create with my data. If that makes sense what I want to illustrate and its clear too then please help me in plotting similar type of image.

I will highly appreciate.

[Friederike]

I do not know what the values represent that you have in your matrix, but since you're complaining that the rnorm values aren't it, what's wrong with this:

> test <- data.frame(state = c("Tss","PcRepr","Heter","Enhp","Con"),
                   Emission = c(0.52547247, 0.028550739,0.001659524,0.004981548,0.118632661),
                   H3K27me3 = c(0.013872603, 0.004967125, 0.002706355, 0.183992447, 0.82344082),
                   H3K9me3 = c(0.012417431, 0.000949458, 0.000396299, 0.002332633, 0.395105767)
)
> test.long <- reshape2::melt(test, id.vars = c("Emission","state"))
> ggplot(test.long, aes(x = variable, y = state, fill = value)) + geom_tile()

[sinha.puja]

Hi, I have modified the codes with the values as mentioned below:

test <- `data.frame`(state = c("Tss","PcRepr","Heter","Enhp","Con"),
                   H3K27me3 = c(0.52547247, 0.028550739,0.001659524,0.004981548,0.118632661),
                   H3K9me3   = c(0.013872603, 0.004967125, 0.002706355, 0.183992447, 0.82344082),
                   H3K36me3 = c(0.012417431, 0.000949458, 0.000396299, 0.002332633, 0.395105767)
)
> **test.long <- reshape2::melt(test, id.vars = c("Emission","state"))**

I do not know how to edit the step bold above and also in ggplot step, x = mark has to be used. Kindly help me in the test.long step which has marked with bold. Sorry for not presenting it properly.

[Friederike]

since you no longer have a column named "Emission", did you try what happens by just removing that entry from the reshape command? If you want to keep the column with the histone mark labels named "mark", you'll have to add variable.name = "mark" to the reshape2:melt command.

Answer 1

[zx8754]

Because we are using random data - rnorm, we are plotting output of below lines, which is overwriting your input file data object:

mat_data <- expand.grid(mark = c("H3K27me3", "H3K9me3", "H3K36me3"), 
                        state = c("TssP", "PcRepr", "Heter/low", "EnhP/low", "ConHeter"))

mat_data$Probability <- rnorm(nrow(mat_data))

[sinha.puja]

Hi,

Thank you for replying. Could you please suggest me what function I should use instead of rnorm so that the data do not get overwritten.

Thanks!!

[Friederike]

well, what do you want to plot? Emission order? H3K27me3 values?

[sinha.puja]

I want to plot the emission order values by stating chromatin states ("TssP", "PcRepr", "Heter/low", "EnhP/low", "ConHeter") on Y-axis and mark ("H3K27me3", "H3K9me3", "H3K36me3") on X-axis with the above mentioned emission values. And also the plotted graph shows gradient color bar indicating the title probability with color gradient ranging from 0.00, 0.25, 0.50, 0.75 & 1.00.

[Friederike]

Please provide an example plot that looks similar to what you want, either a sketch or something that's been published. The geom_tile approach doesn't make sense to me if you want to plot three (!) types of values: emission, histone marks, and probability?

You will probably also have to reformat your data.frame; here's my best guess at what you might be after:

> test <- data.frame(state = c("Tss","PcRepr","Heter","Enhp","Con"),
                   Emission = c(0.52547247, 0.028550739,0.001659524,0.004981548,0.118632661),
                   H3K27me3 = c(0.013872603, 0.004967125, 0.002706355, 0.183992447, 0.82344082),
                   H3K9me3 = c(0.012417431, 0.000949458, 0.000396299, 0.002332633, 0.395105767)
)
> test.long <- reshape2::melt(test, id.vars = c("Emission","state"))
> test.long$random <- rnorm(n = nrow(test.long)) ## disclaimer: no idea what this is supposed to be, just trying to get some values here for coloring
> test.long
      Emission  state variable       value      random
1  0.525472470    Tss H3K27me3 0.013872603  0.15772817
2  0.028550739 PcRepr H3K27me3 0.004967125 -0.86393629
3  0.001659524  Heter H3K27me3 0.002706355 -0.07175823
4  0.004981548   Enhp H3K27me3 0.183992447  0.15042093
5  0.118632661    Con H3K27me3 0.823440820 -0.26406319
6  0.525472470    Tss  H3K9me3 0.012417431  0.31408124
7  0.028550739 PcRepr  H3K9me3 0.000949458  1.50120646
8  0.001659524  Heter  H3K9me3 0.000396299 -2.09819795
9  0.004981548   Enhp  H3K9me3 0.002332633 -0.12286316
10 0.118632661    Con  H3K9me3 0.395105767  0.61224557

## this will ignore the actual emission values because they're the same
## for either Con, Enhp, Heter etc.
> ggplot(test.long, aes(x = variable, y = value, fill = random)) + 
   geom_bar(stat = "identity", position = position_dodge()) + 
   facet_wrap(~state, scales = "free_y") + ylab("Histone mark value") +
  theme_bw(base_size = 16)