I followed the code at https://rpubs.com/mjulkowska/Cowpea06SideView2022 for this 7-side view image analysis.
getwd()
## [1] "C:/Users/Julkowska Lab/Desktop/R codes by Maryam/20220120_Whole_Pheno_3tomatos_continouts_salt_soil_Images"
my_files <- list.files(pattern = "traits.csv")
length(my_files)
## [1] 4
data <- read.csv(my_files[1])
my_data <- data[,c(31, 18:25,27:30)]
my_data
#Let’s split the ROI into specific information - including the timestamp!
for(i in 1:nrow(my_data)){
my_data$timestamp[i] <- strsplit(my_data$roi[i], "_")[[1]][1]
}
my_data
#install.packages("doBy")
library(doBy)
sum_data <- summaryBy(area + convex_hull_area + solidity + perimeter + width + height ~ timestamp, data = my_data, FUN = function(x) c(mean = mean(x), sum=sum(x), max=max(x)))
sum_data
let’s keep only the information we find valuable:
colnames(sum_data)
## [1] "timestamp" "area.mean" "area.sum"
## [4] "area.max" "convex_hull_area.mean" "convex_hull_area.sum"
## [7] "convex_hull_area.max" "solidity.mean" "solidity.sum"
## [10] "solidity.max" "perimeter.mean" "perimeter.sum"
## [13] "perimeter.max" "width.mean" "width.sum"
## [16] "width.max" "height.mean" "height.sum"
## [19] "height.max"
sum_clean <- sum_data[,c(1, 3, 5, 8, 11, 16, 19)]
sum_clean
Great - before moving forward - let’s add also info on DAY:
for(i in 1:nrow(sum_clean)){
sum_clean$day[i] <- strsplit(sum_clean$timestamp[i], "-")[[1]][1]
sum_clean$time[i] <- strsplit(sum_clean$timestamp[i], "-")[[1]][2]
}
sum_clean <- sum_clean[,c(1,8:9,2:7)]
sum_clean
#I have pot# 1 to 58 with #17 and 18 have never been imaged.
pots <- 1:58
nots <- c(17,18)
pots <- subset(pots, !(pots %in% nots))
pots
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 19 20 21 22 23 24 25 26 27
## [26] 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
## [51] 53 54 55 56 57 58
sum_clean$pot.no <- pots
sum_clean
CP06 <- sum_clean
OK - now let’s do next day but much quicker:
data <- read.csv(my_files[2])
my_data <- data[,c(31, 18:25,27:30)]
for(i in 1:nrow(my_data)){
my_data$timestamp[i] <- strsplit(my_data$roi[i], "_")[[1]][1]
}
sum_data <- summaryBy(area + convex_hull_area + solidity + perimeter + width + height ~ timestamp, data = my_data, FUN = function(x) c(mean = mean(x), sum=sum(x), max=max(x)))
sum_clean <- sum_data[,c(1, 3, 5, 8, 11, 16, 19)]
for(i in 1:nrow(sum_clean)){
sum_clean$day[i] <- strsplit(sum_clean$timestamp[i], "-")[[1]][1]
sum_clean$time[i] <- strsplit(sum_clean$timestamp[i], "-")[[1]][2]
}
sum_clean <- sum_clean[,c(1,8:9,2:7)]
sum_clean
pots <- 1:58
nots <- c(17,18)
pots <- subset(pots, !(pots %in% nots))
pots
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 19 20 21 22 23 24 25 26 27
## [26] 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
## [51] 53 54 55 56 57 58
sum_clean$pot.no <- pots
sum_clean
CP06 <- rbind(CP06, sum_clean)
data <- read.csv(my_files[3])
my_data <- data[,c(31, 18:25,27:30)]
for(i in 1:nrow(my_data)){
my_data$timestamp[i] <- strsplit(my_data$roi[i], "_")[[1]][1]
}
sum_data <- summaryBy(area + convex_hull_area + solidity + perimeter + width + height ~ timestamp, data = my_data, FUN = function(x) c(mean = mean(x), sum=sum(x), max=max(x)))
sum_clean <- sum_data[,c(1, 3, 5, 8, 11, 16, 19)]
for(i in 1:nrow(sum_clean)){
sum_clean$day[i] <- strsplit(sum_clean$timestamp[i], "-")[[1]][1]
sum_clean$time[i] <- strsplit(sum_clean$timestamp[i], "-")[[1]][2]
}
sum_clean <- sum_clean[,c(1,8:9,2:7)]
sum_clean
pots <- 1:58
nots <- c(17,18)
pots <- subset(pots, !(pots %in% nots))
pots
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 19 20 21 22 23 24 25 26 27
## [26] 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
## [51] 53 54 55 56 57 58
sum_clean$pot.no <- pots
sum_clean
CP06 <- rbind(CP06, sum_clean)
data <- read.csv(my_files[4])
my_data <- data[,c(31, 18:25,27:30)]
for(i in 1:nrow(my_data)){
my_data$timestamp[i] <- strsplit(my_data$roi[i], "_")[[1]][1]
}
sum_data <- summaryBy(area + convex_hull_area + solidity + perimeter + width + height ~ timestamp, data = my_data, FUN = function(x) c(mean = mean(x), sum=sum(x), max=max(x)))
sum_clean <- sum_data[,c(1, 3, 5, 8, 11, 16, 19)]
for(i in 1:nrow(sum_clean)){
sum_clean$day[i] <- strsplit(sum_clean$timestamp[i], "-")[[1]][1]
sum_clean$time[i] <- strsplit(sum_clean$timestamp[i], "-")[[1]][2]
}
sum_clean <- sum_clean[,c(1,8:9,2:7)]
sum_clean
pots <- 1:58
nots <- c(17,18)
pots <- subset(pots, !(pots %in% nots))
pots
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 19 20 21 22 23 24 25 26 27
## [26] 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
## [51] 53 54 55 56 57 58
sum_clean$pot.no <- pots
sum_clean
CP06 <- rbind(CP06, sum_clean)
####Decoding the pot information using the FW data We will load the data containing fresh and dry weight, as well as all of the decoding information:
decode <- read.csv("Total.FW.3accessions.continous.salt.csv")
decode
FW_data <- decode
decode <- decode[,1:4]
decode <- subset(decode, !(decode$pot.no %in% nots))
CP06 <- CP06[,c(1:4,10,4:9)]
CP06
CP06_decoded <- merge(CP06, decode, all = TRUE)
unique(CP06_decoded$day)
## [1] "2022.02.11" "2022.02.16" "2022.02.14" "2022.02.18"
CP06_decoded
Visualize the data over time: For this - we need to rename the date into a numeric day of salt exposure:
CP06_decoded$day <- gsub("2022.02.11", 28, CP06_decoded$day)
CP06_decoded$day <- gsub("2022.02.14", 31, CP06_decoded$day)
CP06_decoded$day <- gsub("2022.02.16", 33, CP06_decoded$day)
CP06_decoded$day <- gsub("2022.02.18", 35, CP06_decoded$day)
CP06_decoded$day <- as.numeric(as.character(CP06_decoded$day))
CP06_decoded
Now - let’s visualize the area.sum over time for each pot, and divide colours by treatments:
unique(CP06_decoded$Condition)
## [1] "control" "salt"
CP06_decoded$Condition <- factor(CP06_decoded$Condition, levels = c("control", "salt"))
library(ggplot2)
library(ggpubr)
library("ggsci")
CP06_decoded$day <- as.factor(CP06_decoded$day)
Area_lgraph_CP06 <- ggplot(data=CP06_decoded, aes(x= day, y=area.sum, group = pot.no, color = Condition))
Area_lgraph_CP06 <- Area_lgraph_CP06 + geom_line(alpha = 0.1)
Area_lgraph_CP06 <- Area_lgraph_CP06 + stat_summary(fun.data = mean_se, geom="ribbon", linetype=0, aes(group= Condition), alpha=0.3)
Area_lgraph_CP06 <- Area_lgraph_CP06 + stat_summary(fun=mean, aes(group= Condition), size=0.7, geom="line", linetype = "dashed")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
Area_lgraph_CP06 <- Area_lgraph_CP06 + stat_compare_means(aes(group = Condition), label = "p.signif", method = "t.test", hide.ns = T)
Area_lgraph_CP06 <- Area_lgraph_CP06 + ylab("Shoot Size (7 x SV pixels)") + xlab("Days After Stress") + scale_color_jco()
Area_lgraph_CP06
Let’s save the data too into a clean and separate csv file:
write.csv(CP06_decoded, "Maryam_3tomatoes_continuous_salt_Clean_data.csv", row.names = FALSE)
Visualize the correlation between FW and pixels at final day:
last_day <- subset(CP06_decoded, CP06_decoded$day == 35)
last_day
unique(last_day$Condition)
## [1] control salt
## Levels: control salt
unique(FW_data$Condition)
## [1] "control" "salt"
FW_data$Condition <- gsub(" ", "", FW_data$Condition)
last_day_FW <- merge(last_day, FW_data, by=c("pot.no", "Accession", "Condition"), all=TRUE)
last_day_FW
last_day_C <- subset(last_day_FW, last_day_FW$Condition == "control")
last_day_S <- subset(last_day_FW, last_day_FW$Condition == "salt")
FW_Area_C <- ggscatter(last_day_C, x = "area.sum", y = "Total.FW.g.x",rug = TRUE) + stat_cor()
FW_Area_C
FW_Area_S <- ggscatter(last_day_S, x = "area.sum", y = "Total.FW.g.x",rug = TRUE) + stat_cor()
FW_Area_S
I would like to continue the code to see growth for varoius accessions in my experiment…
Area_lgraph_CP06 <- ggplot(data=CP06_decoded, aes(x= day, y=area.sum, group = pot.no, color = Condition))
Area_lgraph_CP06 <- Area_lgraph_CP06 + geom_line(alpha = 0.1) + facet_wrap(~ Accession)
Area_lgraph_CP06 <- Area_lgraph_CP06 + stat_summary(fun.data = mean_se, geom="ribbon", linetype=0, aes(group= Condition), alpha=0.3)
Area_lgraph_CP06 <- Area_lgraph_CP06 + stat_summary(fun=mean, aes(group= Condition), size=0.7, geom="line", linetype = "dashed")
Area_lgraph_CP06 <- Area_lgraph_CP06 + stat_compare_means(aes(group = Condition), label = "p.signif", method = "t.test", hide.ns = T)
Area_lgraph_CP06 <- Area_lgraph_CP06 + ylab("Shoot Size (7 x SV pixels)") + xlab("Days After Stress") + scale_color_aaas()
Area_lgraph_CP06
####to change the color to visualize better the difference between Accessions