GH Arabidopsis growth
This notebook describes the analysis of Arabidopsis growth in different soil types by the BTI’s Greenhouse summer intern (Ella Goggan) under supervision of Julie Bell. Data analyzed by Magda :)
list.files(pattern = "csv")## [1] "coding.csv" "GH_All-single-value-traits.csv"
## [3] "PA-single-value-traits.csv" "PB-single-value-traits.csv"
## [5] "QA-single-value-traits.csv" "QB-single-value-traits.csv"data <- read.csv("GH_All-single-value-traits.csv")
dataunique(data$day)## [1] 2 27 29 26 1 3 30 31 28 25data$doe <- paste(data$mnth, data$day, sep="_")
unique(data$doe)## [1] "8_2" "7_27" "7_29" "7_26" "8_1" "8_3" "7_30" "7_31" "7_28" "7_25"data$doe <- gsub("7_25", "0", data$doe)
data$doe <- gsub("7_26", "1", data$doe)
data$doe <- gsub("7_27", "2", data$doe)
data$doe <- gsub("7_28", "3", data$doe)
data$doe <- gsub("7_29", "4", data$doe)
data$doe <- gsub("7_30", "5", data$doe)
data$doe <- gsub("7_31", "6", data$doe)
data$doe <- gsub("8_1", "7", data$doe)
data$doe <- gsub("8_2", "8", data$doe)
data$doe <- gsub("8_3", "9", data$doe)
data$doe <- as.numeric(as.character(data$doe))
data$Plant <- paste(data$PiID, data$CameraID, data$Plant.ID, sep="_")
library(ggplot2)
library(ggsci)
Arabidopsis_graph <- ggplot(data=data, aes(x= doe, y=area, group = Plant))
Arabidopsis_graph <- Arabidopsis_graph + geom_line(alpha = 0.1)
#EVT_lgraph <- EVT_lgraph + facet_wrap(~ Treatment2)
#EVT_lgraph <- EVT_lgraph + stat_summary(fun.data = mean_se, geom="ribbon", linetype=0, aes(group= Treatment1), alpha=0.3)
#EVT_lgraph <- EVT_lgraph + stat_summary(fun=mean, aes(group= Treatment1), size=0.7, geom="line", linetype = "dashed")
#EVT_lgraph <- EVT_lgraph + stat_compare_means(aes(group = Treatment1), label = "p.signif", method = "t.test", hide.ns = T)
Arabidopsis_graph <- Arabidopsis_graph + ylab("Rosette size (pixels)") + xlab("Days After Treatment") #+ color_palette("aaas") + theme(legend.position = c(0.2, 0.8))
Arabidopsis_graph
coding <- read.csv("coding.csv")
coding$Plant <- paste(coding$PiID, coding$CameraID, coding$Plant.ID, sep="_")
library(reshape2)
coding <- coding[,c(5,4)]
data_decoded <- merge(data, coding, by="Plant")
data_decodeddata_decoded$soil.type <- factor(data_decoded$soil.type, levels = c("Cornell.Mix", "Cornell.Plus.Osmo", "SS8",
"BK25",
"BK25.plus", "BM7.Bark", "BM2.Germ", "BM2.NF.Wood",
"BM2.HP.NF.Wood"))
Arabidopsis_graph <- ggplot(data=data_decoded, aes(x= doe, y=area, group = Plant, color = soil.type))
Arabidopsis_graph <- Arabidopsis_graph + geom_line(alpha = 0.1)
Arabidopsis_graph <- Arabidopsis_graph + facet_wrap(~ soil.type)
Arabidopsis_graph <- Arabidopsis_graph + stat_summary(fun.data = mean_se, geom="ribbon", linetype=0, aes(group= soil.type), alpha=0.3)
Arabidopsis_graph <- Arabidopsis_graph + stat_summary(fun=mean, aes(group= soil.type), 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.Arabidopsis_graph <- Arabidopsis_graph + ylab("Rosette size (pixels)") + xlab("Days Of Experiment") + theme(legend.position = "none") + scale_color_jco()
Arabidopsis_graph
Growth Rate calculations
temp <- subset(data_decoded, data_decoded$Plant == unique(data_decoded$Plant)[1])
tempplot(temp$area ~ temp$doe) + abline(lm(temp$area ~ temp$doe))
## integer(0)model <- lm(temp$area ~ temp$doe)
model##
## Call:
## lm(formula = temp$area ~ temp$doe)
##
## Coefficients:
## (Intercept) temp$doe
## -1524 2188GR <- model$coefficients[2]
R2 <- summary(model)$r.squarednames <- c(text="ID", "soil.type", "GR", "R2")
growth_data <- data.frame()
for (k in names) growth_data[[k]] <- as.character()
growth_datagrowth_data[1,1] <- temp$Plant[1]
growth_data[1,2] <- as.character(temp$soil.type[1])
growth_data[1,3] <- GR
growth_data[1,4] <- R2
growth_dataall_plants <- unique(data_decoded$Plant)
for(i in 1:length(all_plants)){
temp <- subset(data_decoded, data_decoded$Plant == all_plants[i])
model <- lm(temp$area ~ temp$doe)
GR <- model$coefficients[2]
R2 <- summary(model)$r.squared
growth_data[i,1] <- temp$Plant[1]
growth_data[i,2] <- as.character(temp$soil.type[1])
growth_data[i,3] <- GR
growth_data[i,4] <- R2
}
growth_data library(ggpubr)
growth_data$GR <- as.numeric(as.character(growth_data$GR))
unique(growth_data$soil.type)## [1] "BM2.NF.Wood" "BM2.Germ" "BM2.HP.NF.Wood"
## [4] "BK25.plus" "BM7.Bark" "Cornell.Mix"
## [7] "Cornell.Plus.Osmo" "SS8" "BK25"growth_data$soil.type <- factor(growth_data$soil.type, levels = c("Cornell.Mix", "Cornell.Plus.Osmo", "SS8",
"BK25",
"BK25.plus", "BM7.Bark", "BM2.Germ", "BM2.NF.Wood",
"BM2.HP.NF.Wood"))
GR_overall <- ggerrorplot(growth_data, y = "GR", x = "soil.type", fill="soil.type", color="soil.type",
desc_stat = "mean_sd", add = "jitter",
add.params = list(color = "darkgray"),
xlab="Soil Type", ylab="Growth Rate (pix / day)")
#GR_overall <- GR_overall + rremove("legend") + stat_compare_means(method="t.test", ref.group = "Control",
# label = "p.signif", hide.ns = T)
GR_overall <- GR_overall + scale_color_jco() + theme(axis.text.x = element_text(angle=90, vjust=0.5)) + theme(legend.position = "none") + stat_compare_means(method = "aov", label.y = 7000)
GR_overall
library(stringr)
library(multcompView)
Output <- TukeyHSD(aov(GR ~ soil.type, data = growth_data))
Output## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = GR ~ soil.type, data = growth_data)
##
## $soil.type
## diff lwr upr p adj
## Cornell.Plus.Osmo-Cornell.Mix -793.967677 -2807.10899 1219.1736 0.9063523
## SS8-Cornell.Mix -703.311616 -2716.45293 1309.8297 0.9502405
## BK25-Cornell.Mix -1177.367677 -3329.50627 974.7709 0.6408807
## BK25.plus-Cornell.Mix -1335.234343 -3487.37293 816.9042 0.4863870
## BM7.Bark-Cornell.Mix -624.759596 -2776.89819 1527.3790 0.9830764
## BM2.Germ-Cornell.Mix 2065.150505 52.00919 4078.2918 0.0413819
## BM2.NF.Wood-Cornell.Mix -699.161616 -2712.30293 1313.9797 0.9518302
## BM2.HP.NF.Wood-Cornell.Mix -802.126768 -2815.26808 1211.0145 0.9015202
## SS8-Cornell.Plus.Osmo 90.656061 -1773.15063 1954.4628 1.0000000
## BK25-Cornell.Plus.Osmo -383.400000 -2396.54131 1629.7413 0.9989929
## BK25.plus-Cornell.Plus.Osmo -541.266667 -2554.40798 1471.8746 0.9895312
## BM7.Bark-Cornell.Plus.Osmo 169.208081 -1843.93323 2182.3494 0.9999979
## BM2.Germ-Cornell.Plus.Osmo 2859.118182 995.31149 4722.9249 0.0007362
## BM2.NF.Wood-Cornell.Plus.Osmo 94.806061 -1769.00063 1958.6128 1.0000000
## BM2.HP.NF.Wood-Cornell.Plus.Osmo -8.159091 -1871.96578 1855.6476 1.0000000
## BK25-SS8 -474.056061 -2487.19738 1539.0853 0.9956091
## BK25.plus-SS8 -631.922727 -2645.06404 1381.2186 0.9729156
## BM7.Bark-SS8 78.552020 -1934.58929 2091.6933 1.0000000
## BM2.Germ-SS8 2768.462121 904.65543 4632.2688 0.0010935
## BM2.NF.Wood-SS8 4.150000 -1859.65669 1867.9567 1.0000000
## BM2.HP.NF.Wood-SS8 -98.815152 -1962.62184 1764.9915 0.9999999
## BK25.plus-BK25 -157.866667 -2310.00526 1994.2719 0.9999993
## BM7.Bark-BK25 552.608081 -1599.53051 2704.7467 0.9922169
## BM2.Germ-BK25 3242.518182 1229.37687 5255.6595 0.0003951
## BM2.NF.Wood-BK25 478.206061 -1534.93525 2491.3474 0.9953447
## BM2.HP.NF.Wood-BK25 375.240909 -1637.90041 2388.3822 0.9991365
## BM7.Bark-BK25.plus 710.474747 -1441.66384 2862.6133 0.9637368
## BM2.Germ-BK25.plus 3400.384848 1387.24353 5413.5262 0.0002097
## BM2.NF.Wood-BK25.plus 636.072727 -1377.06859 2649.2140 0.9718546
## BM2.HP.NF.Wood-BK25.plus 533.107576 -1480.03374 2546.2489 0.9905006
## BM2.Germ-BM7.Bark 2689.910101 676.76879 4703.0514 0.0036699
## BM2.NF.Wood-BM7.Bark -74.402020 -2087.54333 1938.7393 1.0000000
## BM2.HP.NF.Wood-BM7.Bark -177.367172 -2190.50849 1835.7741 0.9999970
## BM2.NF.Wood-BM2.Germ -2764.312121 -4628.11881 -900.5054 0.0011135
## BM2.HP.NF.Wood-BM2.Germ -2867.277273 -4731.08397 -1003.4706 0.0007105
## BM2.HP.NF.Wood-BM2.NF.Wood -102.965152 -1966.77184 1760.8415 0.9999999P7 = Output$soil.type[,'p adj']
P7## Cornell.Plus.Osmo-Cornell.Mix SS8-Cornell.Mix
## 0.9063522881 0.9502405187
## BK25-Cornell.Mix BK25.plus-Cornell.Mix
## 0.6408806890 0.4863869775
## BM7.Bark-Cornell.Mix BM2.Germ-Cornell.Mix
## 0.9830763519 0.0413819327
## BM2.NF.Wood-Cornell.Mix BM2.HP.NF.Wood-Cornell.Mix
## 0.9518301773 0.9015201842
## SS8-Cornell.Plus.Osmo BK25-Cornell.Plus.Osmo
## 0.9999999723 0.9989929100
## BK25.plus-Cornell.Plus.Osmo BM7.Bark-Cornell.Plus.Osmo
## 0.9895311649 0.9999979250
## BM2.Germ-Cornell.Plus.Osmo BM2.NF.Wood-Cornell.Plus.Osmo
## 0.0007362188 0.9999999605
## BM2.HP.NF.Wood-Cornell.Plus.Osmo BK25-SS8
## 1.0000000000 0.9956090694
## BK25.plus-SS8 BM7.Bark-SS8
## 0.9729155625 0.9999999952
## BM2.Germ-SS8 BM2.NF.Wood-SS8
## 0.0010934844 1.0000000000
## BM2.HP.NF.Wood-SS8 BK25.plus-BK25
## 0.9999999451 0.9999992867
## BM7.Bark-BK25 BM2.Germ-BK25
## 0.9922169270 0.0003951410
## BM2.NF.Wood-BK25 BM2.HP.NF.Wood-BK25
## 0.9953446909 0.9991364734
## BM7.Bark-BK25.plus BM2.Germ-BK25.plus
## 0.9637368084 0.0002096896
## BM2.NF.Wood-BK25.plus BM2.HP.NF.Wood-BK25.plus
## 0.9718545580 0.9905006376
## BM2.Germ-BM7.Bark BM2.NF.Wood-BM7.Bark
## 0.0036699390 0.9999999969
## BM2.HP.NF.Wood-BM7.Bark BM2.NF.Wood-BM2.Germ
## 0.9999970021 0.0011134712
## BM2.HP.NF.Wood-BM2.Germ BM2.HP.NF.Wood-BM2.NF.Wood
## 0.0007104796 0.9999999240stat.test<- multcompLetters(P7)
stat.test## Cornell.Plus.Osmo SS8 BK25 BK25.plus
## "a" "a" "a" "a"
## BM7.Bark BM2.Germ BM2.NF.Wood BM2.HP.NF.Wood
## "a" "b" "a" "a"
## Cornell.Mix
## "a"test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"
testGR_overall <- ggerrorplot(growth_data, y = "GR", x = "soil.type", fill="soil.type", color="soil.type",
desc_stat = "mean_sd", add = "jitter",
add.params = list(color = "darkgray"),
xlab="Soil Type", ylab="Growth Rate (pix / day)")
#GR_overall <- GR_overall + rremove("legend") + stat_compare_means(method="t.test", ref.group = "Control",
# label = "p.signif", hide.ns = T)
GR_overall <- GR_overall + scale_color_jco() + theme(axis.text.x = element_text(angle=90, vjust=0.5)) + theme(legend.position = "none") + stat_compare_means(method = "aov", label.y = 7000)
GR_overall <- GR_overall + stat_pvalue_manual(test, label = "Tukey", y.position = 6000)
GR_overall
pdf("Soil.graph.Arabidopsis.pdf", height = 5, width = 10)
plot(GR_overall)
dev.off()## quartz_off_screen
## 2pdf("Soil.graph.Arabidopsis.size.progression.pdf")
plot(Arabidopsis_graph)
dev.off()## quartz_off_screen
## 2