How to make multiple boxplots by two different groups in one graph?

Question

The part of dataset is like this:

  Treatment   Status     gene1    gene2
1      Both Deceased  3.1934860 63.8697194   
2      Both Deceased  0.0000000 11.3436426   
3     Chemo Deceased  7.2186817 35.0621681   
4      Both Deceased  7.2186817 23.7185255   
5     Chemo Deceased  0.8049256 17.7083638   
6     Chemo Censored  0.8250437  0.8250437   
7     Chemo Censored  3.4136505 23.895533   
8     Radio Censored  0.9428735  4.7143673   
9      None Censored  3.3001750 10.7255686   

I want to make compare each gene expression in "deceased" vs "censored" for each treatment. I only could make one gene expression for now, which is like this:

ggboxplot(df, x="Treatment", y= "gene1", fill = "Status")

Is there any way I can combine two genes' boxplots in one graph? Or any other better way to show these genes expression level difference between deceased vs censored in each group?

Answer 1

using the data from jay.sf you can try a 'ggplot'. I'm using the tidyverse, but this is not required.

library(tidyverse)
dat %>% 
  as_tibble() %>%
  gather(gene, mRNA, -Treatment, -Status) %>% 
  ggplot(aes(Status, mRNA, fill =gene)) + 
   geom_boxplot() +
   facet_wrap(~Treatment, ncol = 2, scales = "free_y")

and with facet_grid you can add significance levels automatically

dat %>% 
  as_tibble() %>%
  gather(gene, mRNA, -Treatment, -Status) %>% 
  ggplot(aes(gene, mRNA, fill =gene)) + 
   geom_boxplot(show.legend = F) +
   ggbeeswarm::geom_beeswarm(show.legend = F) +
   ggsignif::geom_signif(comparisons = list(c("gene1", "gene2"))) +
   facet_grid(Status~Treatment, scales = "free_y") 

Answer 2

We may use boxplot() in base R, where we need to use reshape() first to get a long format.

boxplot(gene ~ Status + time + Treatment, 
        reshape(cbind(id=rownames(dat), dat), 4:5, sep="", direction="long"), 
        border=1:2)

However, this yields a quite crowded plot. We could do separate boxplots for e.g. each treatment group using sapply().

par(mfrow=c(2, 2))
sapply(unique(dat$Treatment), function(x) {
  boxplot(value ~ Status + gene, 
          reshape(cbind(id=rownames(dat[dat$Treatment == x, ]), dat[dat$Treatment == x, ]), 
                  4:5, sep="", direction="long", v.names="value", timevar="gene"), 
          at=c(1:2, 4:5),
          main=x,
          border=1:2)
})

Result

Data

dat <- structure(list(Treatment = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 
3L, 3L, 3L, 4L, 4L, 4L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 
4L, 4L), .Label = c("Both", "Chemo", "None", "Radio"), class = "factor"), 
    Status = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L
    ), .Label = c("Censored", "Deceased"), class = "factor"), 
    gene1 = c(2.83185327992901, 5.21658677992433, 9.36719279899948, 
    1.77809421116808, 6.39453760571561, 3.08376117126782, -1.99524072673447, 
    0.380722587753265, -0.947148460332481, 1.73014054712629, 
    0.855919162512028, 0.501667581598007, 0.0638735169737497, 
    10.1712355237258, 5.34317645471502, -7.96626158445742, -0.0781613844302278, 
    5.59930916967042, -0.725717330717595, 0.492793009977729, 
    -0.546677404630108, 0.290301979542245, 2.83540215865274, 
    -1.25738031049913), gene2 = c(6.97361394841868, -6.86012827859373, 
    -0.193731972798249, -5.64669185350061, -20.6664537342379, 
    32.5477488386544, 12.6210452154023, 6.56845245925654, 13.5491140544121, 
    -2.9113829554538, 2.90958200298303, -6.56806056188421, 50.2577234864485, 
    17.0734922804668, 49.0769939658538, -2.0186433516603, 32.3823429023035, 
    17.7654319738005, 12.2884241568455, 21.7600566866782, 19.68978862329, 
    -12.6277420840716, 27.555120882401, 17.5164450232983)), row.names = c(3L, 
23L, 13L, 44L, 34L, 50L, 90L, 67L, 62L, 100L, 95L, 96L, 132L, 
144L, 124L, 174L, 171L, 168L, 196L, 205L, 207L, 233L, 229L, 212L
), class = "data.frame")