Andy Color distribution

This notebook describes a data analysis pipeline for Andy Raptis - a student who is working with Magda Julkowska at the remote internship, to describe the color characteristics of flowers.

Load data:

First - let’s load a set of data that Andy has collected during his internship:

list.files(pattern = ".csv")

 [1] "Misc.1.csv"                 "Misc.2InsidePetal.csv"      "Misc.2OutsidePetal.csv"     "Misc.3.csv"                
 [5] "Misc.4.csv"                 "Misc.5.csv"                 "Misc.6.csv"                 "Sunflower.3.csv"           
 [9] "Sunflower.4.csv"            "Sunflower.5.csv"            "Sunflower.6.csv"            "Sunflower.Control.csv"     
[13] "Sunflower.InsidePetal.csv"  "Sunflower.OutsidePetal.csv" "Sunflower.Petal.csv"       

Sunflower1 <- read.csv("Sunflower.Control.csv")
Sunflower2 <- read.csv("Sunflower.3.csv")
Sunflower1

Sunflower2

Let’s remove all of the empty columns:

Sunflower1 <- Sunflower1[,1:4]
Sunflower2 <- Sunflower2[,1:4]
Sunflower1

Sunflower2

OK - now let’s remove empty rows with NAs:

Sunflower1 <- na.omit(Sunflower1)
Sunflower2 <- na.omit(Sunflower2)
Sunflower1

Sunflower2

Awesome - 1st column is actually called “Pixel value” and I would like to rename it as such:

colnames(Sunflower1)[1] <- "Pix.values"
colnames(Sunflower2)[1] <- "Pix.values"
Sunflower1

Sunflower2

OK - so now what I would like to do is compare the distribution of pixel values of two sunflowers - to do that - I would need to merge the frames.

To be able to distinguish them - I have to add another collumn called “sample ID”:

Sunflower1$SampleID <- "Sunflower01"
Sunflower2$SampleID <- "Sunflower03"
Sunflower1

Sunflower2

OK - now let’s merge the two data sets together:

library(reshape2)
Sunflower <- rbind(Sunflower1, Sunflower2)
Sunflower

to make sure we have two data sets correctly merged - I will check if my unique SampleIDs in merged data frame are indeed what I expect:

unique(Sunflower$SampleID)

[1] "Sunflower01" "Sunflower03"

OK - looks good - moving on to plotting:

library(ggpubr)

Loading required package: ggplot2

library(ggsci)
library(ggplot2)

Red <- ggplot(data = Sunflower, aes(x = Pix.values, y = Red, colour = SampleID, group = SampleID)) 
Red <- Red + geom_line(alpha = 0.5) + scale_color_aaas() 
Red

Now - let’s compare another color dimension:

Green <- ggplot(data = Sunflower, aes(x = Pix.values, y = Green, colour = SampleID, group = SampleID)) 
Green <- Green + geom_line(alpha = 0.5) + scale_color_aaas() 
Green

Blue <- ggplot(data = Sunflower, aes(x = Pix.values, y = Blue, colour = SampleID, group = SampleID)) 
Blue <- Blue + geom_line(alpha = 0.5) + scale_color_aaas() 
Blue

OK - this all looks good. And we can even plot all of these into one graph, but before we need to melt the dataset:

SunM <- melt(Sunflower, id = c("SampleID", "Pix.values"))
SunM

RGB_plot <- ggplot(data = SunM, aes(x = Pix.values, y = value, colour = SampleID, group = SampleID)) 
RGB_plot <- RGB_plot + geom_line(alpha = 0.5) 
RGB_plot <- RGB_plot + scale_color_aaas() 
RGB_plot <- RGB_plot + facet_grid(~ variable)
RGB_plot

Now - we can see what channel gives us most difference between the two samples we would like to compare. However - where are the differences statistically significant?

In order to calculate this - we would have to change the values in our tables to actual data.

How the table looks like for each RGB variable - it represents the % of pixels.

What we need is to have the data of points that will have that specific value. Let’s do an example here below:

SunM

So I dont want to have any of the pixels that have value of 0 - so let’s remove them:

SunPix <- subset(SunM, SunM$value > 0)
SunPix

Then - I would like to have 134 pixels that is having value of 140 in Red, same for 145 and so on…. but how do I do it?

Maybe let’s first multiply everything by a thousand so we have only the WHOLE observations:

SunPix$obs <- SunPix$value * 1000
SunPix

OK - now we need to re-calculate all of this into individual observations. Let’s start with the first one to begin with:

temp <- SunPix[1,1:3]
temp

reps <- as.numeric(SunPix$obs[1])
temp2 <- do.call("rbind", replicate(reps, temp, simplify = FALSE))
SunCount <- temp2
SunCount

cool - now let’s loop it for all of the rows in the SunPix dataframe:

for(i in 2:nrow(SunPix)){
  temp <- SunPix[i,1:3]
  temp
  reps <- as.numeric(SunPix$obs[i])
  temp2 <- do.call("rbind", replicate(reps, temp, simplify = FALSE))
  SunCount <- rbind(SunCount, temp2)
}

head(SunCount)

dim(SunCount)

[1] 599962      3

OK - now that we have our data in correct format, let’s plot it again, but using slighltly different methods:

library(ggsci)
library(ggpubr)

color_plot <- ggviolin(SunCount, x = "SampleID", y = "Pix.values", color = "SampleID", fill = "SampleID", 
                          facet.by = "variable",
                          desc_stat = "mean_sd", add = "boxplot", add.params = list(fill = "white"),
                          xlab="", ylab= "number of pixels")
color_plot <- color_plot + stat_compare_means(aes(label = after_stat(p.signif)), method = "aov", hide.ns = F)
color_plot <- color_plot + rremove("legend") + scale_color_d3("category10") + scale_fill_d3("category10")
color_plot <- color_plot + theme(axis.text.x = element_text(angle = 90))
color_plot