Akhilesh Experiment 01
This notebook is describing a quick pipeline for analyzing just the PS2 data collected by Akhilesh Sharma, a PhD student in Vipan Kumar’s lab of tolerant and susceptible plants exposed to a herbicide treatment.
The data was collected in BTI’s PhenoSight facility in January/February of 2025 by Akhilesh, analyzed using PhenoVation software, and implemented here in standard data processing pipeline by Magda Julkowska.
Data import and wrangling:
list.files(pattern = "Akhilesh")## [1] "20250205_Akhilesh_MMJ.html" "20250205_Akhilesh_MMJ.nb.html"
## [3] "20250205_Akhilesh_MMJ.Rmd" "20250206_Akhilesh_data.csv"
## [5] "20250213_Akhilesh_MMJ.nb.html" "20250213_Akhilesh_MMJ.Rmd"
## [7] "Akhilesh_14dR_Jan25_Analysis.TXT" "Akhilesh_14dS_Jan15_Analysis.TXT"
## [9] "Akhilesh_1dR_Jan9_Analysis.TXT" "Akhilesh_1dS_Analysis_MMJ.txt"
## [11] "Akhilesh_21dR_Jan30_Analysis.TXT" "Akhilesh_21dS_Jan22_Analysis.TXT"
## [13] "Akhilesh_7dR_Jan15_Analysis.TXT" "Akhilesh_7dS_Jan7_Analysis.TXT"
## [15] "Akhilesh_NoSprayR_Jan6_Analysis.TXT" "Akhilesh_NoSprayS_Analysis_MMJ.txt"loading the datasets:
D0R <- read.table("Akhilesh_NoSprayR_Jan6_Analysis.TXT", header = T, sep = '\t')## Warning in read.table("Akhilesh_NoSprayR_Jan6_Analysis.TXT", header = T, :
## incomplete final line found by readTableHeader on
## 'Akhilesh_NoSprayR_Jan6_Analysis.TXT'D1R <- read.table("Akhilesh_1dR_Jan9_Analysis.TXT", header = T, sep = '\t')## Warning in read.table("Akhilesh_1dR_Jan9_Analysis.TXT", header = T, sep =
## "\t"): incomplete final line found by readTableHeader on
## 'Akhilesh_1dR_Jan9_Analysis.TXT'D7R <- read.table("Akhilesh_7dR_Jan15_Analysis.TXT", header = T, sep = '\t')## Warning in read.table("Akhilesh_7dR_Jan15_Analysis.TXT", header = T, sep =
## "\t"): incomplete final line found by readTableHeader on
## 'Akhilesh_7dR_Jan15_Analysis.TXT'D14R <- read.table("Akhilesh_14dR_Jan25_Analysis.TXT", header = T, sep = '\t')## Warning in read.table("Akhilesh_14dR_Jan25_Analysis.TXT", header = T, sep =
## "\t"): incomplete final line found by readTableHeader on
## 'Akhilesh_14dR_Jan25_Analysis.TXT'D21R <- read.table("Akhilesh_21dR_Jan30_Analysis.TXT", header = T, sep = '\t')## Warning in read.table("Akhilesh_21dR_Jan30_Analysis.TXT", header = T, sep =
## "\t"): incomplete final line found by readTableHeader on
## 'Akhilesh_21dR_Jan30_Analysis.TXT'D0S <- read.table("Akhilesh_NoSprayS_Analysis_MMJ.txt", header = T, sep = '\t')## Warning in read.table("Akhilesh_NoSprayS_Analysis_MMJ.txt", header = T, :
## incomplete final line found by readTableHeader on
## 'Akhilesh_NoSprayS_Analysis_MMJ.txt'D1S <- read.table("Akhilesh_1dS_Analysis_MMJ.txt", header = T, sep = '\t')## Warning in read.table("Akhilesh_1dS_Analysis_MMJ.txt", header = T, sep = "\t"):
## incomplete final line found by readTableHeader on
## 'Akhilesh_1dS_Analysis_MMJ.txt'D7S <- read.table("Akhilesh_7dS_Jan7_Analysis.TXT", header = T, sep = '\t')## Warning in read.table("Akhilesh_7dS_Jan7_Analysis.TXT", header = T, sep =
## "\t"): incomplete final line found by readTableHeader on
## 'Akhilesh_7dS_Jan7_Analysis.TXT'D14S <- read.table("Akhilesh_14dS_Jan15_Analysis.TXT", header = T, sep = '\t')## Warning in read.table("Akhilesh_14dS_Jan15_Analysis.TXT", header = T, sep =
## "\t"): incomplete final line found by readTableHeader on
## 'Akhilesh_14dS_Jan15_Analysis.TXT'D21S <- read.table("Akhilesh_21dS_Jan22_Analysis.TXT", header = T, sep = '\t')## Warning in read.table("Akhilesh_21dS_Jan22_Analysis.TXT", header = T, sep =
## "\t"): incomplete final line found by readTableHeader on
## 'Akhilesh_21dS_Jan22_Analysis.TXT'D0RBefore fusing the files - let’s add identifiers to each of them:
D0R$ID <- "D0R"
D1R$ID <- "D1R"
D7R$ID <- "D7R"
D14R$ID <- "D14R"
D21R$ID <- "D21R"
D0S$ID <- "D0S"
D1S$ID <- "D1S"
D7S$ID <- "D7S"
D14S$ID <- "D14S"
D21S$ID <- "D21S"Let’s check if they all have same dimensions:
dim(D0R)## [1] 247 106dim(D1R)## [1] 297 106dim(D7R)## [1] 389 106dim(D14R)## [1] 349 106dim(D21R)## [1] 259 80dim(D0S)## [1] 287 62dim(D1S)## [1] 309 62dim(D7S)## [1] 397 106dim(D14S)## [1] 329 106dim(D21S)## [1] 415 106Looks good
Now - lets keep only the collumns that are meaningful for each dataset:
D0R <- D0R[,c(1:4,6,9,11:15,35,39,41,43,45,47,49,51,53,55,57,59, 106)]
colnames(D0R)## [1] "File" "Date" "Time" "Obj.No" "Obj.Size" "Fv.Fm"
## [7] "Day" "treatment" "plant" "rate" "Fq.Fm" "NPQ"
## [13] "qP" "qN" "qL" "qE" "qI" "D.no"
## [19] "D.npq" "npq.t." "ChlIdx" "AriIdx" "NDVI" "ID"D1R <- D1R[,c(1:4,6,9,11:15,35,39,41,43,45,47,49,51,53,55,57,59, 106)]
colnames(D1R)## [1] "File" "Date" "Time" "Obj.No" "Obj.Size" "Fv.Fm"
## [7] "Day" "treatment" "plant" "rate" "Fq.Fm" "NPQ"
## [13] "qP" "qN" "qL" "qE" "qI" "D.no"
## [19] "D.npq" "npq.t." "ChlIdx" "AriIdx" "NDVI" "ID"D7R <- D7R[,c(1:4,6,9,11:15,35,39,41,43,45,47,49,51,53,55,57,59, 106)]
colnames(D7R)## [1] "File" "Date" "Time" "Obj.No" "Obj.Size" "Fv.Fm"
## [7] "Day" "treatment" "plant" "rate" "Fq.Fm" "NPQ"
## [13] "qP" "qN" "qL" "qE" "qI" "D.no"
## [19] "D.npq" "npq.t." "ChlIdx" "AriIdx" "NDVI" "ID"D14R <- D14R[,c(1:4,6,9,11:15,35,39,41,43,45,47,49,51,53,55,57,59, 106)]
colnames(D14R)## [1] "File" "Date" "Time" "Obj.No" "Obj.Size" "Fv.Fm"
## [7] "Day" "treatment" "plant" "rate" "Fq.Fm" "NPQ"
## [13] "qP" "qN" "qL" "qE" "qI" "D.no"
## [19] "D.npq" "npq.t." "ChlIdx" "AriIdx" "NDVI" "ID"D21R <- D21R[,c(1:4, 6, 9, 11:15, 17, 21, 23, 25, 27, 29, 31, 33, 35, 55, 57, 59, 80)]
colnames(D21R)## [1] "File" "Date" "Time" "Obj.No" "Obj.Size" "Fv.Fm"
## [7] "Day" "treatment" "plant" "rate" "Fq.Fm" "NPQ"
## [13] "qP" "qN" "qL" "qE" "qI" "D.no"
## [19] "D.npq" "npq.t." "ChlIdx" "AriIdx" "NDVI" "ID"D0S <- D0S[,c(1:4, 6, 9, 11:15, 17, 21, 23, 25, 27, 29, 31, 33, 35, 37, 39, 41, 62)]
colnames(D0S)## [1] "File" "Date" "Time" "Obj.No" "Obj.Size" "Fv.Fm"
## [7] "Day" "treatment" "plant" "rate" "Fq.Fm" "NPQ"
## [13] "qP" "qN" "qL" "qE" "qI" "D.no"
## [19] "D.npq" "npq.t." "ChlIdx" "AriIdx" "NDVI" "ID"D1S <- D1S[,c(1:4, 6, 9, 11:15, 17, 21, 23, 25, 27, 29, 31, 33, 35, 37, 39, 41, 62)]
colnames(D1S)[7] <- "Day"
D1S$Day <- 1
D1SD7S <- D7S[,c(1:4,6,9,11:15,35,39,41,43,45,47,49,51,53,55,57,59, 106)]
colnames(D7S)## [1] "File" "Date" "Time" "Obj.No" "Obj.Size" "Fv.Fm"
## [7] "Day" "treatment" "plant" "rate" "Fq.Fm" "NPQ"
## [13] "qP" "qN" "qL" "qE" "qI" "D.no"
## [19] "D.npq" "npq.t." "ChlIdx" "AriIdx" "NDVI" "ID"D14S <- D14S[,c(1:4,6,9,11:15,35,39,41,43,45,47,49,51,53,55,57,59, 106)]
colnames(D14S)## [1] "File" "Date" "Time" "Obj.No" "Obj.Size" "Fv.Fm"
## [7] "Day" "treatment" "plant" "rate" "Fq.Fm" "NPQ"
## [13] "qP" "qN" "qL" "qE" "qI" "D.no"
## [19] "D.npq" "npq.t." "ChlIdx" "AriIdx" "NDVI" "ID"D21S <- D21S[,c(1:4,6,9,11:15,35,39,41,43,45,47,49,51,53,55,57,59, 106)]
colnames(D21S)## [1] "File" "Date" "Time" "Obj.No" "Obj.Size" "Fv.Fm"
## [7] "Day" "treatment" "plant" "rate" "Fq.Fm" "NPQ"
## [13] "qP" "qN" "qL" "qE" "qI" "D.no"
## [19] "D.npq" "npq.t." "ChlIdx" "AriIdx" "NDVI" "ID"Now - let’s fuse them:
library(reshape2)
all_data <- rbind(D0R, D1R)
all_data <- rbind(all_data, D7R)
all_data <- rbind(all_data, D14R)
all_data <- rbind(all_data, D21R)
all_data <- rbind(all_data, D0S)
all_data <- rbind(all_data, D1S)
all_data <- rbind(all_data, D7S)
all_data <- rbind(all_data, D14S)
all_data <- rbind(all_data, D21S)
all_dataComment from Akhilesh: please not include picture 6 or its data from susceptible_7day_jan7
So let’s exclude this subset of data:
unique(all_data$File)## [1] "HDR_E0001P0001N0001_GCB22060646_20250106.INF"
## [2] "HDR_E0001P0002N0001_GCB22060646_20250106.INF"
## [3] "HDR_E0001P0003N0001_GCB22060646_20250106.INF"
## [4] "HDR_E0001P0004N0001_GCB22060646_20250106.INF"
## [5] "HDR_E0001P0005N0001_GCB22060646_20250106.INF"
## [6] "HDR_E0001P0006N0001_GCB22060646_20250106.INF"
## [7] "HDR_E0001P0001N0001_GCB22060646_20250109.INF"
## [8] "HDR_E0001P0002N0001_GCB22060646_20250109.INF"
## [9] "HDR_E0001P0003N0001_GCB22060646_20250109.INF"
## [10] "HDR_E0001P0004N0001_GCB22060646_20250109.INF"
## [11] "HDR_E0001P0005N0001_GCB22060646_20250109.INF"
## [12] "HDR_E0001P0006N0001_GCB22060646_20250109.INF"
## [13] "HDR_E0001P0001N0001_GCB22060646_20250115.INF"
## [14] "HDR_E0001P0002N0001_GCB22060646_20250115.INF"
## [15] "HDR_E0001P0003N0001_GCB22060646_20250115.INF"
## [16] "HDR_E0001P0004N0001_GCB22060646_20250115.INF"
## [17] "HDR_E0001P0005N0001_GCB22060646_20250115.INF"
## [18] "HDR_E0001P0006N0001_GCB22060646_20250115.INF"
## [19] "HDR_E0001P0001N0001_GCB22060646_20250125.INF"
## [20] "HDR_E0001P0002N0001_GCB22060646_20250125.INF"
## [21] "HDR_E0001P0003N0001_GCB22060646_20250125.INF"
## [22] "HDR_E0001P0004N0001_GCB22060646_20250125.INF"
## [23] "HDR_E0001P0005N0001_GCB22060646_20250125.INF"
## [24] "HDR_E0001P0006N0001_GCB22060646_20250125.INF"
## [25] "HDR_E0001P0007N0001_GCB22060646_20250125.INF"
## [26] "HDR_E0001P0001N0001_GCB22060646_20250130.INF"
## [27] "HDR_E0001P0002N0001_GCB22060646_20250130.INF"
## [28] "HDR_E0001P0003N0001_GCB22060646_20250130.INF"
## [29] "HDR_E0001P0004N0001_GCB22060646_20250130.INF"
## [30] "HDR_E0001P0005N0001_GCB22060646_20250130.INF"
## [31] "HDR_E0001P0005N0002_GCB22060646_20250130.INF"
## [32] "HDR_E0001P0006N0001_GCB22060646_20250130.INF"
## [33] "HDR_E0001P0006N0002_GCB22060646_20250130.INF"
## [34] "HDR_10190_39933_58.INF"
## [35] "HDR_10190_39934_58.INF"
## [36] "HDR_10190_39935_58.INF"
## [37] "HDR_10190_39936_58.INF"
## [38] "HDR_10190_39937_58.INF"
## [39] "HDR_10190_39938_58.INF"
## [40] "HDR_10190_39939_58.INF"
## [41] "HDR_10190_39940_58.INF"
## [42] "HDR_10190_39941_58.INF"
## [43] "HDR_10190_39942_58.INF"
## [44] "HDR_10190_39943_58.INF"
## [45] "HDR_10190_39944_58.INF"
## [46] "HDR_10190_39945_58.INF"
## [47] "HDR_10190_39946_58.INF"
## [48] "HDR_10190_39947_58.INF"
## [49] "HDR_10190_39948_58.INF"
## [50] "HDR_10190_39949_58.INF"
## [51] "HDR_10190_39950_58.INF"
## [52] "HDR_10190_39951_58.INF"
## [53] "HDR_10190_39952_58.INF"
## [54] "HDR_10190_39953_58.INF"
## [55] "HDR_10190_39954_58.INF"
## [56] "HDR_10190_39955_58.INF"
## [57] "HDR_10190_39956_58.INF"
## [58] "HDR_10190_39957_58.INF"
## [59] "HDR_10190_39958_58.INF"
## [60] "HDR_10190_39959_58.INF"
## [61] "HDR_10190_39960_58.INF"
## [62] "HDR_10190_39961_58.INF"
## [63] "HDR_10190_39962_58.INF"
## [64] "HDR_10190_39963_58.INF"
## [65] "HDR_10190_39964_58.INF"
## [66] "HDR_10190_39965_58.INF"
## [67] "HDR_10190_39966_58.INF"
## [68] "HDR_10190_39967_58.INF"
## [69] "HDR_10190_39968_58.INF"
## [70] "HDR_10190_39969_58.INF"
## [71] "HDR_10190_39970_58.INF"
## [72] "HDR_10190_39971_58.INF"
## [73] "HDR_10190_39972_58.INF"
## [74] "HDR_10190_39973_58.INF"
## [75] "HDR_10190_39974_58.INF"
## [76] "HDR_10190_39975_58.INF"
## [77] "HDR_10190_39976_58.INF"
## [78] "HDR_10190_39977_58.INF"
## [79] "HDR_10190_39978_58.INF"
## [80] "HDR_10190_39979_58.INF"
## [81] "HDR_10190_39980_58.INF"
## [82] "HDR_10190_39981_58.INF"
## [83] "HDR_10190_39982_58.INF"
## [84] "HDR_10190_39983_58.INF"
## [85] "HDR_10190_39984_58.INF"
## [86] "HDR_10190_39985_58.INF"
## [87] "HDR_10190_39986_58.INF"
## [88] "HDR_10190_39987_58.INF"
## [89] "HDR_10190_39988_58.INF"
## [90] "HDR_10190_39989_58.INF"
## [91] "HDR_10190_39990_58.INF"
## [92] "HDR_10190_39991_58.INF"
## [93] "HDR_10190_39992_58.INF"
## [94] "HDR_10190_39993_58.INF"
## [95] "HDR_10190_39994_58.INF"
## [96] "HDR_10190_39995_58.INF"
## [97] "HDR_10190_39996_58.INF"
## [98] "HDR_10190_39933_65.INF"
## [99] "HDR_10190_39934_65.INF"
## [100] "HDR_10190_39935_65.INF"
## [101] "HDR_10190_39936_65.INF"
## [102] "HDR_10190_39937_65.INF"
## [103] "HDR_10190_39938_65.INF"
## [104] "HDR_10190_39939_65.INF"
## [105] "HDR_10190_39940_65.INF"
## [106] "HDR_10190_39941_65.INF"
## [107] "HDR_10190_39942_65.INF"
## [108] "HDR_10190_39943_65.INF"
## [109] "HDR_10190_39944_65.INF"
## [110] "HDR_10190_39945_65.INF"
## [111] "HDR_10190_39946_65.INF"
## [112] "HDR_10190_39947_65.INF"
## [113] "HDR_10190_39948_65.INF"
## [114] "HDR_10190_39949_65.INF"
## [115] "HDR_10190_39950_65.INF"
## [116] "HDR_10190_39951_65.INF"
## [117] "HDR_10190_39952_65.INF"
## [118] "HDR_10190_39953_65.INF"
## [119] "HDR_10190_39954_65.INF"
## [120] "HDR_10190_39955_65.INF"
## [121] "HDR_10190_39956_65.INF"
## [122] "HDR_10190_39957_65.INF"
## [123] "HDR_10190_39958_65.INF"
## [124] "HDR_10190_39959_65.INF"
## [125] "HDR_10190_39960_65.INF"
## [126] "HDR_10190_39961_65.INF"
## [127] "HDR_10190_39962_65.INF"
## [128] "HDR_10190_39963_65.INF"
## [129] "HDR_10190_39964_65.INF"
## [130] "HDR_10190_39965_65.INF"
## [131] "HDR_10190_39966_65.INF"
## [132] "HDR_10190_39967_65.INF"
## [133] "HDR_10190_39968_65.INF"
## [134] "HDR_10190_39969_65.INF"
## [135] "HDR_10190_39970_65.INF"
## [136] "HDR_10190_39971_65.INF"
## [137] "HDR_10190_39972_65.INF"
## [138] "HDR_10190_39973_65.INF"
## [139] "HDR_10190_39974_65.INF"
## [140] "HDR_10190_39975_65.INF"
## [141] "HDR_10190_39976_65.INF"
## [142] "HDR_10190_39977_65.INF"
## [143] "HDR_10190_39978_65.INF"
## [144] "HDR_10190_39979_65.INF"
## [145] "HDR_10190_39980_65.INF"
## [146] "HDR_10190_39981_65.INF"
## [147] "HDR_10190_39982_65.INF"
## [148] "HDR_10190_39983_65.INF"
## [149] "HDR_10190_39984_65.INF"
## [150] "HDR_10190_39985_65.INF"
## [151] "HDR_10190_39986_65.INF"
## [152] "HDR_10190_39987_65.INF"
## [153] "HDR_10190_39988_65.INF"
## [154] "HDR_10190_39989_65.INF"
## [155] "HDR_10190_39990_65.INF"
## [156] "HDR_10190_39991_65.INF"
## [157] "HDR_10190_39992_65.INF"
## [158] "HDR_10190_39993_65.INF"
## [159] "HDR_10190_39994_65.INF"
## [160] "HDR_10190_39995_65.INF"
## [161] "HDR_10190_39996_65.INF"
## [162] "HDR_E0001P0001N0001_GCB22060646_20250107.INF"
## [163] "HDR_E0001P0002N0001_GCB22060646_20250107.INF"
## [164] "HDR_E0001P0003N0001_GCB22060646_20250107.INF"
## [165] "HDR_E0001P0004N0001_GCB22060646_20250107.INF"
## [166] "HDR_E0001P0005N0001_GCB22060646_20250107.INF"
## [167] "HDR_E0001P0006N0001_GCB22060646_20250107.INF"
## [168] "HDR_E0001P0007N0001_GCB22060646_20250115.INF"
## [169] "HDR_E0001P0008N0001_GCB22060646_20250115.INF"
## [170] "HDR_E0001P0009N0001_GCB22060646_20250115.INF"
## [171] "HDR_E0001P0010N0001_GCB22060646_20250115.INF"
## [172] "HDR_E0001P0011N0001_GCB22060646_20250115.INF"
## [173] "HDR_E0001P0001N0001_GCB22060646_20250122.INF"
## [174] "HDR_E0001P0002N0001_GCB22060646_20250122.INF"
## [175] "HDR_E0001P0003N0001_GCB22060646_20250122.INF"
## [176] "HDR_E0001P0004N0001_GCB22060646_20250122.INF"
## [177] "HDR_E0001P0005N0001_GCB22060646_20250122.INF"
## [178] "HDR_E0001P0006N0001_GCB22060646_20250122.INF"all_data2 <- subset(all_data, all_data$File != "HDR_E0001P0006N0001_GCB22060646_20250107.INF")
all_data2OK - cool - now - lets make some calculations
But first - let’s remove the NA’s
all_data3 <- na.omit(all_data2)
all_data3Before going forward - let’s fix this odd issue with some rows having the right FvFm and others right all of the other values.
FvFm_data <- subset(all_data3, all_data3$Fv.Fm > 0)
notFvFm_data <- subset(all_data3, all_data3$Fv.Fm < 0)
FvFm_datanotFvFm_dataFvFm_data <- FvFm_data[,c(1:10,24)]
notFvFm_data <- notFvFm_data[,c(1:5,7:24)]
FvFm_datanotFvFm_dataall_correct <- merge(FvFm_data, notFvFm_data, by = c("File", "Date", "Time", "Obj.No", "Obj.Size", "Day", "treatment", "plant", "rate", "ID"))
all_correctThen - we need to sum all of the object size belonging to one plant - and average all of the other values
#install.packages('doBy')
#install.packages('summaryBy')
library(doBy)
sum_data <- summaryBy(data = all_correct, Obj.Size + Fv.Fm + Fq.Fm + NPQ + qP + qN + qL + qE + qI + D.no + D.npq + npq.t. + ChlIdx + AriIdx + NDVI ~ File + Date + Day + treatment + plant + rate + ID, FUN=function(x) c(sum=sum(x), median=median(x)))
sum_datanow - let’s get rid of things we dont need anymore:
data1 <- sum_data[,c(1:8,11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33, 35, 37)]
data1data1$Plant.ID <- paste(data1$treatment, data1$rate, data1$plant, sep=".")
data1Ploting data
OK - let’s try to plot it:
#install.packages('ggpubr')
#install.packages('ggsci')
library(ggplot2)
library(ggpubr)
library(ggsci)
min(data1$Obj.Size.sum)## [1] 18data1 <- subset(data1, data1$Obj.Size.sum > 100)
data1 <- subset(data1, data1$Obj.Size.sum < 431000)
want <- c("0x", "0.25x", "0.5x", "1x", "2x", "4x", "8x")
data1 <- subset(data1, data1$rate %in% want)unique(data1$rate)## [1] "0x" "0.25x" "1x" "0.5x" "2x" "8x" "4x"data1$rate <- factor(data1$rate, levels = c("0x", "0.25x", "0.5x", "1x", "2x", "4x", "8x"))
ObjSize_lgraph <- ggplot(data = data1, aes(x = Day, y = Obj.Size.sum, group = Plant.ID, color = treatment))
ObjSize_lgraph <- ObjSize_lgraph + geom_line(alpha = 0.3) + facet_wrap(~ rate, ncol = 7)
ObjSize_lgraph <- ObjSize_lgraph + stat_summary(fun.data = mean_se, geom = "ribbon", linetype = 0, aes(group = treatment), alpha = 0.3)
ObjSize_lgraph <- ObjSize_lgraph + stat_summary(fun = mean, aes(group = treatment), size = 0.7, geom = "line", linetype = "dashed") ## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.ObjSize_lgraph <- ObjSize_lgraph + stat_compare_means(aes(group = treatment), label = "p.signif", method = "t.test", hide.ns = F)
ObjSize_lgraph <- ObjSize_lgraph + labs(x = "Time After Treatment (days)", y = "Object Size (a.u.)") + scale_color_d3("category10") + theme(legend.position = "top")
ObjSize_lgraph## Warning: Computation failed in `stat_compare_means()`
## Caused by error in `purrr::map()`:
## ℹ In index: 1.
## Caused by error in `.check_npc_coord()`:
## ! '*.npc coord for x axis should be either a numeric value in [0-1] or a character strings including one of right, left, center, centre, middle## Warning: Unknown or uninitialised column: `p`.## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! argument "x" is missing, with no default## Warning: Computation failed in `stat_compare_means()`
## Computation failed in `stat_compare_means()`
## Computation failed in `stat_compare_means()`
## Caused by error in `purrr::map()`:
## ℹ In index: 1.
## Caused by error in `.check_npc_coord()`:
## ! '*.npc coord for x axis should be either a numeric value in [0-1] or a character strings including one of right, left, center, centre, middle## Warning: Unknown or uninitialised column: `p`.## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! argument "x" is missing, with no default## Warning: Unknown or uninitialised column: `p`.## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! argument "x" is missing, with no default
Let’s save it as PDF:
pdf("Object.size.pdf")
plot(ObjSize_lgraph)## Warning: Computation failed in `stat_compare_means()`
## Caused by error in `purrr::map()`:
## ℹ In index: 1.
## Caused by error in `.check_npc_coord()`:
## ! '*.npc coord for x axis should be either a numeric value in [0-1] or a character strings including one of right, left, center, centre, middle## Warning: Unknown or uninitialised column: `p`.## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! argument "x" is missing, with no default## Warning: Computation failed in `stat_compare_means()`
## Computation failed in `stat_compare_means()`
## Computation failed in `stat_compare_means()`
## Caused by error in `purrr::map()`:
## ℹ In index: 1.
## Caused by error in `.check_npc_coord()`:
## ! '*.npc coord for x axis should be either a numeric value in [0-1] or a character strings including one of right, left, center, centre, middle## Warning: Unknown or uninitialised column: `p`.## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! argument "x" is missing, with no default## Warning: Unknown or uninitialised column: `p`.## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! argument "x" is missing, with no defaultdev.off()## quartz_off_screen
## 2Ploting per treatment
OK - so Akhilesh wants me to plot all of the resistant lines together and susceptible lines together. Let’s do this.
Res_vs_Susc_lgraph <- ggplot(data = data1, aes(x = Day, y = Obj.Size.sum, group = Plant.ID, color = rate))
#Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + geom_line(alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + facet_wrap(~ treatment, ncol = 7)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun.data = mean_se, geom = "ribbon", linetype = 0, aes(group = rate), alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun = mean, aes(group = rate), size = 0.7, geom = "line", linetype = "dashed")
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_compare_means(aes(group = rate), label = "p.signif", method = "t.test", hide.ns = F)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + labs(x = "Time After Treatment (days)", y = "Object Size (a.u.)") + scale_color_d3("category10") + theme(legend.position = "top")
Res_vs_Susc_lgraph## Warning: Computation failed in `stat_compare_means()`
## Computation failed in `stat_compare_means()`
## Caused by error in `purrr::map()`:
## ℹ In index: 1.
## Caused by error in `.check_npc_coord()`:
## ! '*.npc coord for x axis should be either a numeric value in [0-1] or a character strings including one of right, left, center, centre, middle
So this graph is a bit confusing - since you cannot really see what treatments show significant differences between eachother. But hey - let’s save it!
Let’s save this as PDF file:
pdf("Res_vs_Susc_ObjectSize.pdf")
plot(Res_vs_Susc_lgraph)## Warning: Computation failed in `stat_compare_means()`
## Computation failed in `stat_compare_means()`
## Caused by error in `purrr::map()`:
## ℹ In index: 1.
## Caused by error in `.check_npc_coord()`:
## ! '*.npc coord for x axis should be either a numeric value in [0-1] or a character strings including one of right, left, center, centre, middledev.off()## quartz_off_screen
## 2There is a better way of plotting it, but for that - we need to have the 0x treatment replicated for each comparison:
sub0x <- subset(data1,data1$rate == "0x")
sub0.25x <- subset(data1,data1$rate == "0.25x")
sub0.5x <- subset(data1,data1$rate == "0.5x")
sub1x <- subset(data1,data1$rate == "1x")
sub2x <- subset(data1,data1$rate == "2x")
sub4x <- subset(data1,data1$rate == "4x")
sub8x <- subset(data1,data1$rate == "8x")
temp_sub <- rbind(sub0x, sub0.25x)
temp_sub$comparison <- "0.25x"
temp_sub2 <- rbind(sub0x, sub0.5x)
temp_sub2$comparison <- "0.5x"
data2 <- rbind(temp_sub, temp_sub2)
temp_sub <- rbind(sub0x, sub1x)
temp_sub$comparison <- "1x"
temp_sub2 <- rbind(sub0x, sub2x)
temp_sub2$comparison <- "2x"
data2 <- rbind(data2, temp_sub, temp_sub2)
temp_sub <- rbind(sub0x, sub4x)
temp_sub$comparison <- "4x"
temp_sub2 <- rbind(sub0x, sub8x)
temp_sub2$comparison <- "8x"
data2 <- rbind(data2, temp_sub, temp_sub2)
data2write.csv(data2, "Data2.csv", row.names = F)OK - now we can plot the graph where we compare everything to the respective growth at 0x:
Res_vs_Susc_lgraph <- ggplot(data = data2, aes(x = Day, y = Obj.Size.sum, group = Plant.ID, color = rate))
# Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + geom_line(alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + facet_wrap(treatment ~ comparison, ncol = 6)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun.data = mean_se, geom = "ribbon", linetype = 0, aes(group = rate), alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun = mean, aes(group = rate), size = 0.7, geom = "line", linetype = "dashed")
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + labs(x = "Time After Treatment (days)", y = "Object Size (a.u.)") + scale_color_d3("category10") + theme(legend.position = "top")
Res_vs_Susc_lgraph
Let’s save this one:
pdf("Res_vs_Susc_ObjectSize_better.pdf")
plot(Res_vs_Susc_lgraph)
dev.off()## quartz_off_screen
## 2OK lets make more graphs for all of the other traits:
data2Res_vs_Susc_lgraph <- ggplot(data = data2, aes(x = Day, y = Fv.Fm.median, group = Plant.ID, color = rate))
# Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + geom_line(alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + facet_wrap(treatment ~ comparison, ncol = 6)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun.data = mean_se, geom = "ribbon", linetype = 0, aes(group = rate), alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun = mean, aes(group = rate), size = 0.7, geom = "line", linetype = "dashed")
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + labs(x = "Time After Treatment (days)", y = "Fv/Fm (a.u.)") + scale_color_d3("category10") + theme(legend.position = "top")
Res_vs_Susc_lgraph
pdf("Res_vs_Susc_FvFm_better.pdf")
plot(Res_vs_Susc_lgraph)
dev.off()## quartz_off_screen
## 2Res_vs_Susc_lgraph <- ggplot(data = data2, aes(x = Day, y = Fq.Fm.median, group = Plant.ID, color = rate))
#Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + geom_line(alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + facet_wrap(treatment ~ comparison, ncol = 6)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun.data = mean_se, geom = "ribbon", linetype = 0, aes(group = rate), alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun = mean, aes(group = rate), size = 0.7, geom = "line", linetype = "dashed")
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + labs(x = "Time After Treatment (days)", y = "Fv'/Fm' (a.u.)") + scale_color_d3("category10") + theme(legend.position = "top")
Res_vs_Susc_lgraph
pdf("Res_vs_Susc_Fv.Fm._better.pdf")
plot(Res_vs_Susc_lgraph)
dev.off()## quartz_off_screen
## 2Res_vs_Susc_lgraph <- ggplot(data = data2, aes(x = Day, y = NPQ.median, group = Plant.ID, color = rate))
#Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + geom_line(alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + facet_wrap(treatment ~ comparison, ncol = 6)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun.data = mean_se, geom = "ribbon", linetype = 0, aes(group = rate), alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun = mean, aes(group = rate), size = 0.7, geom = "line", linetype = "dashed")
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + labs(x = "Time After Treatment (days)", y = "NPQ (a.u.)") + scale_color_d3("category10") + theme(legend.position = "top")
Res_vs_Susc_lgraph
pdf("Res_vs_Susc_NPQ_better.pdf")
plot(Res_vs_Susc_lgraph)
dev.off()## quartz_off_screen
## 2Res_vs_Susc_lgraph <- ggplot(data = data2, aes(x = Day, y = qP.median, group = Plant.ID, color = rate))
# Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + geom_line(alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + facet_wrap(treatment ~ comparison, ncol = 6)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun.data = mean_se, geom = "ribbon", linetype = 0, aes(group = rate), alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun = mean, aes(group = rate), size = 0.7, geom = "line", linetype = "dashed")
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + labs(x = "Time After Treatment (days)", y = "qP (a.u.)") + scale_color_d3("category10") + theme(legend.position = "top")
Res_vs_Susc_lgraph
pdf("Res_vs_Susc_qP_better.pdf")
plot(Res_vs_Susc_lgraph)
dev.off()## quartz_off_screen
## 2Res_vs_Susc_lgraph <- ggplot(data = data2, aes(x = Day, y = qN.median, group = Plant.ID, color = rate))
#Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + geom_line(alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + facet_wrap(treatment ~ comparison, ncol = 6)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun.data = mean_se, geom = "ribbon", linetype = 0, aes(group = rate), alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun = mean, aes(group = rate), size = 0.7, geom = "line", linetype = "dashed")
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + labs(x = "Time After Treatment (days)", y = "qN (a.u.)") + scale_color_d3("category10") + theme(legend.position = "top")
Res_vs_Susc_lgraph
pdf("Res_vs_Susc_qN._better.pdf")
plot(Res_vs_Susc_lgraph)
dev.off()## quartz_off_screen
## 2Res_vs_Susc_lgraph <- ggplot(data = data2, aes(x = Day, y = qL.median, group = Plant.ID, color = rate))
#Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + geom_line(alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + facet_wrap(treatment ~ comparison, ncol = 6)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun.data = mean_se, geom = "ribbon", linetype = 0, aes(group = rate), alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun = mean, aes(group = rate), size = 0.7, geom = "line", linetype = "dashed")
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + labs(x = "Time After Treatment (days)", y = "qL (a.u.)") + scale_color_d3("category10") + theme(legend.position = "top")
Res_vs_Susc_lgraph
pdf("Res_vs_Susc_qL._better.pdf")
plot(Res_vs_Susc_lgraph)
dev.off()## quartz_off_screen
## 2Res_vs_Susc_lgraph <- ggplot(data = data2, aes(x = Day, y = qI.median, group = Plant.ID, color = rate))
# Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + geom_line(alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + facet_wrap(treatment ~ comparison, ncol = 6)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun.data = mean_se, geom = "ribbon", linetype = 0, aes(group = rate), alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun = mean, aes(group = rate), size = 0.7, geom = "line", linetype = "dashed")
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + labs(x = "Time After Treatment (days)", y = "qI (a.u.)") + scale_color_d3("category10") + theme(legend.position = "top")
Res_vs_Susc_lgraph
pdf("Res_vs_Susc_qI._better.pdf")
plot(Res_vs_Susc_lgraph)
dev.off()## quartz_off_screen
## 2Res_vs_Susc_lgraph <- ggplot(data = data2, aes(x = Day, y = D.no.median, group = Plant.ID, color = rate))
# Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + geom_line(alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + facet_wrap(treatment ~ comparison, ncol = 6)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun.data = mean_se, geom = "ribbon", linetype = 0, aes(group = rate), alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun = mean, aes(group = rate), size = 0.7, geom = "line", linetype = "dashed")
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + labs(x = "Time After Treatment (days)", y = "Phi NO (a.u.)") + scale_color_d3("category10") + theme(legend.position = "top")
Res_vs_Susc_lgraph
pdf("Res_vs_Susc_PhiNO._better.pdf")
plot(Res_vs_Susc_lgraph)
dev.off()## quartz_off_screen
## 2Res_vs_Susc_lgraph <- ggplot(data = data2, aes(x = Day, y = D.npq.median, group = Plant.ID, color = rate))
#Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + geom_line(alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + facet_wrap(treatment ~ comparison, ncol = 6)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun.data = mean_se, geom = "ribbon", linetype = 0, aes(group = rate), alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun = mean, aes(group = rate), size = 0.7, geom = "line", linetype = "dashed")
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + labs(x = "Time After Treatment (days)", y = "phi NPQ (a.u.)") + scale_color_d3("category10") + theme(legend.position = "top")
Res_vs_Susc_lgraph
pdf("Res_vs_Susc_phiNPQ._better.pdf")
plot(Res_vs_Susc_lgraph)
dev.off()## quartz_off_screen
## 2Res_vs_Susc_lgraph <- ggplot(data = data2, aes(x = Day, y = npq.t..median, group = Plant.ID, color = rate))
#Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + geom_line(alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + facet_wrap(treatment ~ comparison, ncol = 6)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun.data = mean_se, geom = "ribbon", linetype = 0, aes(group = rate), alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun = mean, aes(group = rate), size = 0.7, geom = "line", linetype = "dashed")
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + labs(x = "Time After Treatment (days)", y = "NPQ(t) (a.u.)") + scale_color_d3("category10") + theme(legend.position = "top")
Res_vs_Susc_lgraph
pdf("Res_vs_Susc_NPQt_better.pdf")
plot(Res_vs_Susc_lgraph)
dev.off()## quartz_off_screen
## 2Res_vs_Susc_lgraph <- ggplot(data = data2, aes(x = Day, y = ChlIdx.median, group = Plant.ID, color = rate))
#Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + geom_line(alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + facet_wrap(treatment ~ comparison, ncol = 6)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun.data = mean_se, geom = "ribbon", linetype = 0, aes(group = rate), alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun = mean, aes(group = rate), size = 0.7, geom = "line", linetype = "dashed")
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + labs(x = "Time After Treatment (days)", y = "ChlIdx (a.u.)") + scale_color_d3("category10") + theme(legend.position = "top")
Res_vs_Susc_lgraph
pdf("Res_vs_Susc_ChlIdx._better.pdf")
plot(Res_vs_Susc_lgraph)
dev.off()## quartz_off_screen
## 2Res_vs_Susc_lgraph <- ggplot(data = data2, aes(x = Day, y = AriIdx.median, group = Plant.ID, color = rate))
#Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + geom_line(alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + facet_wrap(treatment ~ comparison, ncol = 6)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun.data = mean_se, geom = "ribbon", linetype = 0, aes(group = rate), alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun = mean, aes(group = rate), size = 0.7, geom = "line", linetype = "dashed")
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + labs(x = "Time After Treatment (days)", y = "AriIdx (a.u.)") + scale_color_d3("category10") + theme(legend.position = "top")
Res_vs_Susc_lgraph
pdf("Res_vs_Susc_AriIdx_better.pdf")
plot(Res_vs_Susc_lgraph)
dev.off()## quartz_off_screen
## 2Res_vs_Susc_lgraph <- ggplot(data = data2, aes(x = Day, y = NDVI.median, group = Plant.ID, color = rate))
#Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + geom_line(alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + facet_wrap(treatment ~ comparison, ncol = 6)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun.data = mean_se, geom = "ribbon", linetype = 0, aes(group = rate), alpha = 0.3)
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + stat_summary(fun = mean, aes(group = rate), size = 0.7, geom = "line", linetype = "dashed")
Res_vs_Susc_lgraph <- Res_vs_Susc_lgraph + labs(x = "Time After Treatment (days)", y = "NDVI (a.u.)") + scale_color_d3("category10") + theme(legend.position = "top")
Res_vs_Susc_lgraph
pdf("Res_vs_Susc_NDVI_better.pdf")
plot(Res_vs_Susc_lgraph)
dev.off()## quartz_off_screen
## 2Growth rate calculations:
Although this is not the best idea to calculate the growth rate just based on one top view of the plants - this is the data we got - so let’s do this
first - let’s isolate one plant from the dataset:
temp <- subset(data1, data1$Plant.ID == unique(data1$Plant.ID)[1])
# plot the data
plot(temp$Obj.Size.sum ~ temp$Day) +
# Linear regression line is added to the scatterplot using the lm() function
abline(lm(temp$Obj.Size.sum ~ temp$Day))
## integer(0)Let’s have a look at the linear model components:
# Fit a linear regression model to predict 'length' based on 'hours' in the 'temp' data frame
model <- lm(temp$Obj.Size.sum ~ temp$Day)
model##
## Call:
## lm(formula = temp$Obj.Size.sum ~ temp$Day)
##
## Coefficients:
## (Intercept) temp$Day
## 15577 5221Extract the growth rate (GR) coefficient from the linear regression model
GR <- model$coefficients[2]
GR## temp$Day
## 5221.136Calculate the coefficient of determination (R-squared) for the model
R2 <- summary(model)$r.squared
R2## [1] 0.4156404Define vector ‘names’ containing column names for the ‘growth_data’ data frame
tempnames <- c(text = "Treatment", "Plant", "Rate", "Plant.ID", "GR", "R2")
# Create an empty data frame 'growth_data' to store growth-related data
growth_data <- data.frame()
# Loop through each element 'k' in the 'names' vector
for (k in names) {
# Create a new column in 'growth_data' with the name specified by 'k' and set as character type
growth_data[[k]] <- as.character()
}
# Assign the 'treatment' value from the 'temp' data frame to the first row, first column of 'growth_data'
growth_data[1, 1] <- temp$treatment[1]
# Assign the 'plant' value from the 'temp' data frame to the first row, second column of 'growth_data'
growth_data[1, 2] <- as.character(temp$plant[1])
# Assign the 'rate' value from the 'temp' data frame to the first row, third column of 'growth_data'
growth_data[1, 3] <- as.character(temp$rate[1])
# Assign the 'Plant.ID' value from the 'temp' data frame to the first row, fourth column of 'growth_data'
growth_data[1, 4] <- as.character(temp$Plant.ID[1])
# Assign the growth rate (GR) to the first row, fifth column of 'growth_data'
growth_data[1, 5] <- GR
# Assign the coefficient of determination (R-squared, R2) to the first row, sixth column of 'growth_data'
growth_data[1, 6] <- R2
growth_dataLooks great! Now - let’s do this for all of the plants within this experiment:
# Create vector 'all_plants' containing unique 'ID' values from 'final_spline'
all_plants <- unique(data1$Plant.ID)
# Loop through each unique 'ID' in 'all_plants'
for (i in 1:length(all_plants)) {
# Subset 'final_spline_deco' to select data for the current 'ID'
temp <- subset(data1, data1$Plant.ID == all_plants[i])
# Fit a linear regression model to predict 'Area.SUM' based on 'day'
model <- lm(temp$Obj.Size.sum ~ temp$Day)
# Extract the growth rate (GR) coefficient from the linear regression model
GR <- model$coefficients[2]
# Calculate the coefficient of determination (R-squared, R2) for the model
R2 <- summary(model)$r.squared
# Populate 'growth_data' with relevant information for the current 'ID'
growth_data[i, 1] <- temp$treatment[1]
growth_data[i, 2] <- as.character(temp$plant[1])
growth_data[i, 3] <- as.character(temp$rate[1])
growth_data[i, 4] <- as.character(temp$Plant.ID[1])
growth_data[i, 5] <- GR
growth_data[i, 6] <- R2
}
growth_datagrowth_data$GR <- as.numeric(as.character(growth_data$GR))
growth_data$R2 <- as.numeric(as.character(growth_data$R2))
growth_data <- subset(growth_data, growth_data$GR > -2000)
growth_dataGreat - now we have all of this data - let’s plot it:
growth_data$Treatment <- factor(growth_data$Treatment, levels =c("resistant", "susceptible"))
growth_data$Rate <- factor(growth_data$Rate, levels =c("0x", "0.25x", "0.5x", "1x", "2x", "4x", "8x"))
GR_overall <- ggerrorplot(growth_data, x = "Rate", y = "GR", color = "Rate", fill = "Rate", facet.by = "Treatment",
desc_stat = "mean_sd", add = "jitter",
xlab="", ylab= "Growth Rate (pixel / day)", add.params = list(color = "darkgray")) + scale_color_d3("category10")
GR_overall <- GR_overall + stat_compare_means(method = "aov", label.y = 15000)
GR_overall <- GR_overall + rremove("legend") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
GR_overall
OK - this looks great - lets save this:
pdf("Growth rates_btwRates.pdf")
plot(GR_overall)
dev.off()## quartz_off_screen
## 2Save the dataset too:
write.csv(growth_data, "growth_rates.csv", row.names = F)Vipan wants us to calculate regression between dose and GR in response to increasing herbicide concentration. Let’s plot it first:
growth_data$rate2 <- growth_data$Rate
growth_datagrowth_data$rate2 <- gsub("x", "", growth_data$rate2)
growth_data$rate2 <- as.numeric(as.character(growth_data$rate2))
resistant <- subset(growth_data, growth_data$Treatment == "resistant")
susceptible <- subset(growth_data, growth_data$Treatment == "susceptible")#pdf("Resistant_plot_DoseResponse.pdf")
dose <- seq(0, 8, by = 0.25)
plot.spl <- with(resistant, smooth.spline(rate2, GR, df = 5))
# Create a scatterplot of 'Area.sum' against 'days' using data from 'temp'
plot(GR ~ rate2, data = resistant,main = "Resistant Population Dose-Response",
xlab = "Rate", ylab = "Growth Response")
# Add a blue line to the plot based on the cubic smoothing spline ('plot.spl')
lines(plot.spl, col = "blue")
# Add a red line to the plot by predicting values using the cubic smoothing spline ('plot.spl')
lines(predict(plot.spl, dose), col = "red")
#dev.off()#pdf("Susceptible_plot_DoseResponse.pdf")
plot.spl <- with(susceptible, smooth.spline(rate2, GR, df = 5))
# Create a scatterplot of 'Area.sum' against 'days' using data from 'temp'
plot(GR ~ rate2, data = susceptible,main = "Susceptible Population Dose-Response",
xlab = "Rate", ylab = "Growth Response")
# Add a blue line to the plot based on the cubic smoothing spline ('plot.spl')
lines(plot.spl, col = "blue")
# Add a red line to the plot by predicting values using the cubic smoothing spline ('plot.spl')
lines(predict(plot.spl, dose), col = "red")
#dev.off()I would rather extract the 50% of growth inhibition based on the growth rate - rather than object size at 21 days (which is not reliable due to plant growth). I would extract this data directly from the spline:
# first - let's change the dose to something ridiculously high resolution:
dose <- seq(0, 8, by = 0.01)
plot.spl <- with(resistant, smooth.spline(rate2, GR, df = 5))
resistant_predict <- predict(plot.spl, dose)
resistant_predict2 <- as.data.frame(resistant_predict)
resistant_predict2# our max rate is 1954.9510849 - let's see whats the half rate:
# 1954.9510849/2 = 977.4755
resistant_rate <- subset(resistant_predict2, resistant_predict2$y > 970)
resistant_rate <- subset(resistant_rate, resistant_rate$y < 990)
resistant_rate# There we go - the D50 for resistant is 0.55x rate# first - let's change the dose to something ridiculously high resolution:
dose <- seq(0, 8, by = 0.01)
plot.spl <- with(susceptible, smooth.spline(rate2, GR, df = 5))
susceptible_predict <- predict(plot.spl, dose)
susceptible_predict2 <- as.data.frame(susceptible_predict)
susceptible_predict2# our max rate is 3781.65440 - let's see whats the half rate:
# 3781.65440/2 = 1890.827
susceptible_rate <- subset(susceptible_predict2, susceptible_predict2$y > 1800)
susceptible_rate <- subset(susceptible_rate, susceptible_rate$y < 2000)
susceptible_rate# There we go - the D50 for susceptible is 0.09x rateOK - now let’s calculate the Herbicide Tolerance Index:
First - let’s calculate average GR for Control
library(doBy)
sum_growth <- summaryBy(GR ~ Treatment + Rate, data = growth_data)
sum_growthisolate only control samples:
C_sub <- subset(sum_growth, sum_growth$Rate == "0x")
C_sub <- C_sub[,c(1,3)]
colnames(C_sub)[2] <- "GR.0x"
C_subNOw - let’s fuse them with the GR data and merge per Treatment:
growth_data2 <- merge(growth_data, C_sub, by = "Treatment")
growth_data2Awesome - now we can calculate RELATIVE growth rates for each rate of herbicide application - we call it Herbicite Tolerance Index (HTI)
growth_data2$HTI <- growth_data2$GR / growth_data2$GR.0x
growth_data2Cool - let’s plot this now:
HTI <- ggerrorplot(growth_data2, x = "Rate", y = "HTI", color = "Rate", fill = "Rate", facet.by = "Treatment",
desc_stat = "mean_sd", add = "jitter",
xlab="", ylab= "Growth Rate (fraction of 0x)", add.params = list(color = "darkgray")) + scale_color_d3("category10")
HTI <- HTI + stat_compare_means(method = "aov", label.y = 2)
HTI <- HTI + rremove("legend") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
HTI
library(ggbeeswarm)
ObjSize_lgraph2 <- ggplot(data = data1, aes(x = treatment, y = Obj.Size.sum, color = treatment))
ObjSize_lgraph2 <- ObjSize_lgraph2 + geom_beeswarm() + facet_wrap(Day ~ rate, ncol = 6)
ObjSize_lgraph2 <- ObjSize_lgraph2 + stat_summary(fun.data = mean_se, geom = "ribbon", linetype = 0, aes(group = treatment), alpha = 0.3)
ObjSize_lgraph2 <- ObjSize_lgraph2 + stat_summary(fun = mean, aes(group = treatment), size = 0.7, geom = "line", linetype = "dashed")
#ObjSize_lgraph2 <- ObjSize_lgraph2 + stat_compare_means(aes(group = treatment), ref.group = "susceptible", label = "p.signif", method = "t.test", hide.ns = F)
ObjSize_lgraph2 <- ObjSize_lgraph2 + labs(x = "", y = "Object Size (a.u.)") + scale_color_d3("category10") + theme(legend.position = "top")
pdf("Obj.Size.Akhil.comparison.pdf", width = 10, height = 15)
plot(ObjSize_lgraph2)## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?dev.off()## quartz_off_screen
## 2write.csv(data1, "20250206_Akhilesh_data.csv", row.names = F)