getwd()
## [1] "C:/Users/Julkowska Lab/Desktop/R codes by Maryam/20231215_Arabidopsis_soil_salt_IAA_ACC_spray"
list.files(pattern = ".csv")
## [1] "20231122FW_IAA_ACC_Arabidopsis_shoot_spray_soil_salt.csv"
## [2] "all-data-20231207.csv"                                   
## [3] "ICP-Arabidopsis-SOIL-SPRAY-R-Col-only-v2.csv"            
## [4] "ICP-Arabidopsis-SOIL-SPRAY-R-v2.csv"                     
## [5] "ICP-Arabidopsis-SOIL-SPRAY-R.csv"
ICP <- read.csv( "ICP-Arabidopsis-SOIL-SPRAY-R-Col-only-v2.csv")
ICP
colnames(ICP)
## [1] "Accession"               "Condition"              
## [3] "DW.g"                    "DW.mg"                  
## [5] "K.con..mg.mg.dry.weight" "Na.con.mg.mg.dry.weight"
## [7] "Na.K.ratio"
ICP$All.ID<-paste(ICP$Accession, ICP$Condition, sep="_")
ICP
library(ggplot2)
library(ggpubr)
library(multcompView)
## Warning: package 'multcompView' was built under R version 4.3.3
library(cowplot)
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggpubr':
## 
##     get_legend
aov(Na.con.mg.mg.dry.weight ~ All.ID, data = ICP)
## Call:
##    aov(formula = Na.con.mg.mg.dry.weight ~ All.ID, data = ICP)
## 
## Terms:
##                    All.ID Residuals
## Sum of Squares   73.02363 248.36408
## Deg. of Freedom         4        29
## 
## Residual standard error: 2.926479
## Estimated effects may be unbalanced
Output <- TukeyHSD(aov(Na.con.mg.mg.dry.weight ~ All.ID, data = ICP))
Output
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Na.con.mg.mg.dry.weight ~ All.ID, data = ICP)
## 
## $All.ID
##                               diff       lwr      upr     p adj
## col_accmock-col_acc      3.1825546 -1.550161 7.915270 0.3126892
## col_iaa-col_acc         -0.3656145 -4.912660 4.181431 0.9992978
## col_iaamock-col_acc      0.4327853 -4.114260 4.979831 0.9986372
## col_salt-col_acc         2.8104057 -1.736640 7.357451 0.3947218
## col_iaa-col_accmock     -3.5481691 -8.280884 1.184546 0.2158782
## col_iaamock-col_accmock -2.7497693 -7.482484 1.982946 0.4561159
## col_salt-col_accmock    -0.3721489 -5.104864 4.360566 0.9993570
## col_iaamock-col_iaa      0.7983998 -3.748646 5.345445 0.9856571
## col_salt-col_iaa         3.1760202 -1.371025 7.723066 0.2773851
## col_salt-col_iaamock     2.3776204 -2.169425 6.924666 0.5583605
P6 = Output$All.ID[,'p adj']
stat.test<- multcompLetters(P6)
stat.test
## $Letters
## col_accmock     col_iaa col_iaamock    col_salt     col_acc 
##         "a"         "a"         "a"         "a"         "a" 
## 
## $LetterMatrix
##                a
## col_accmock TRUE
## col_iaa     TRUE
## col_iaamock TRUE
## col_salt    TRUE
## col_acc     TRUE
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
#test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"
test
test$info <- strsplit(test$group1, "_")
test$info[[1]][2]
## [1] "accmock"
test$Accession <- "none"
test$Condition<- "none"
test
for(i in 1:nrow(test)){
  test$Accession[i] <- test$info[[i]][1]
     test$Condition[i] <- test$info[[i]][2]
  
}

test2 <- test[,c(4:5,1)]
test2$group1 <- test2$Condition
test2$group2 <- test2$Condition
ICP
ICP$Condition<- factor(ICP$Condition, levels=c("salt","iaamock", "iaa","accmock", "acc"))
ICP2 <- subset(ICP, ICP$Na.con.mg.mg.dry.weight < 15)

Na_content <- ggplot(data = ICP2, mapping = aes(x = Condition, y = Na.con.mg.mg.dry.weight, colour = Condition)) 
Na_content <- Na_content + geom_boxplot(alpha=0.2) + geom_jitter(width=0.1,alpha=0.2)

#Na_content <- Na_content + facet_grid(~Condition , scales = "free_y")

Na_content <- Na_content + stat_summary(fun=mean, geom="point", shape=95, size=6, color="black", fill="black")
Na_content <- Na_content + scale_color_manual(values = c("red","blue","blueviolet","cyan","brown"))
Na_content <- Na_content + ylab("Na content, mg/mg dry weight") + xlab("")+stat_pvalue_manual(test2, label = "Tukey", y.position = 15)
Na_content <- Na_content + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
Na_content <- Na_content + rremove("legend")
Na_content

library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
aov(Na.con.mg.mg.dry.weight ~ All.ID, data = ICP)
## Call:
##    aov(formula = Na.con.mg.mg.dry.weight ~ All.ID, data = ICP)
## 
## Terms:
##                    All.ID Residuals
## Sum of Squares   73.02363 248.36408
## Deg. of Freedom         4        29
## 
## Residual standard error: 2.926479
## Estimated effects may be unbalanced
Output <- TukeyHSD(aov(Na.con.mg.mg.dry.weight ~ All.ID, data = ICP))
Output
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Na.con.mg.mg.dry.weight ~ All.ID, data = ICP)
## 
## $All.ID
##                               diff       lwr      upr     p adj
## col_accmock-col_acc      3.1825546 -1.550161 7.915270 0.3126892
## col_iaa-col_acc         -0.3656145 -4.912660 4.181431 0.9992978
## col_iaamock-col_acc      0.4327853 -4.114260 4.979831 0.9986372
## col_salt-col_acc         2.8104057 -1.736640 7.357451 0.3947218
## col_iaa-col_accmock     -3.5481691 -8.280884 1.184546 0.2158782
## col_iaamock-col_accmock -2.7497693 -7.482484 1.982946 0.4561159
## col_salt-col_accmock    -0.3721489 -5.104864 4.360566 0.9993570
## col_iaamock-col_iaa      0.7983998 -3.748646 5.345445 0.9856571
## col_salt-col_iaa         3.1760202 -1.371025 7.723066 0.2773851
## col_salt-col_iaamock     2.3776204 -2.169425 6.924666 0.5583605
P6 = Output$All.ID[,'p adj']
stat.test<- multcompLetters(P6)
stat.test
## $Letters
## col_accmock     col_iaa col_iaamock    col_salt     col_acc 
##         "a"         "a"         "a"         "a"         "a" 
## 
## $LetterMatrix
##                a
## col_accmock TRUE
## col_iaa     TRUE
## col_iaamock TRUE
## col_salt    TRUE
## col_acc     TRUE
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
#test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"
test
test$info <- strsplit(test$group1, "_")
test$info[[1]][2]
## [1] "accmock"
test$Accession <- "none"
test$Condition<- "none"
test
for(i in 1:nrow(test)){
  test$Accession[i] <- test$info[[i]][1]
     test$Condition[i] <- test$info[[i]][2]
  
}

test2 <- test[,c(4:5,1)]
test2$group1 <- test2$Condition
test2$group2 <- test2$Condition
ICP
# Plot Ttest
ICP2 <- subset(ICP, ICP$Na.con.mg.mg.dry.weight < 15)

Na_content_col <- ggplot(ICP2, aes(x = Condition, y = Na.con.mg.mg.dry.weight, colour = Condition)) + 
  geom_boxplot(alpha = 0.2) + 
  geom_jitter(width = 0.1, alpha = 0.2) + 
  stat_summary(fun = mean, geom = "point", shape = 95, size = 6, color = "black", fill = "black") + 
  scale_color_manual(values = c("red","blue","blueviolet","cyan","brown")) + 
  scale_x_discrete(drop = TRUE) +  # 🔹 Forces ggplot to remove empty categories
  ylab("Na content, mg/mg dry weight") + 
  xlab("") + stat_compare_means(label = "p.signif", method = "t.test", ref.group = "accmock")
  theme(axis.text.x = element_text(angle = 90, hjust = 0.9, vjust = 0.5)) + 
  rremove("legend")
## List of 2
##  $ axis.text.x    :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : num 0.9
##   ..$ vjust        : num 0.5
##   ..$ angle        : num 90
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi FALSE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ legend.position: chr "none"
##  - attr(*, "class")= chr [1:2] "theme" "gg"
##  - attr(*, "complete")= logi FALSE
##  - attr(*, "validate")= logi TRUE
# Display plot
Na_content_col

aov(K.con..mg.mg.dry.weight ~ All.ID, data = ICP)
## Call:
##    aov(formula = K.con..mg.mg.dry.weight ~ All.ID, data = ICP)
## 
## Terms:
##                   All.ID Residuals
## Sum of Squares   20.6293  720.8996
## Deg. of Freedom        4        29
## 
## Residual standard error: 4.985841
## Estimated effects may be unbalanced
Output <- TukeyHSD(aov(K.con..mg.mg.dry.weight~ All.ID, data = ICP))
Output
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = K.con..mg.mg.dry.weight ~ All.ID, data = ICP)
## 
## $All.ID
##                                diff       lwr       upr     p adj
## col_accmock-col_acc     -1.29019399 -9.353319  6.772931 0.9898801
## col_iaa-col_acc         -1.24693876 -8.993738  6.499861 0.9896514
## col_iaamock-col_acc      0.13464211 -7.612157  7.881442 0.9999984
## col_salt-col_acc         0.66585882 -7.080941  8.412658 0.9990869
## col_iaa-col_accmock      0.04325523 -8.019869  8.106380 1.0000000
## col_iaamock-col_accmock  1.42483610 -6.638289  9.487961 0.9853117
## col_salt-col_accmock     1.95605281 -6.107072 10.019177 0.9536551
## col_iaamock-col_iaa      1.38158087 -6.365219  9.128380 0.9847998
## col_salt-col_iaa         1.91279759 -5.834002  9.659597 0.9507054
## col_salt-col_iaamock     0.53121671 -7.215583  8.278016 0.9996257
P5 = Output$All.ID[,'p adj']
stat.test<- multcompLetters(P5)
stat.test
## $Letters
## col_accmock     col_iaa col_iaamock    col_salt     col_acc 
##         "a"         "a"         "a"         "a"         "a" 
## 
## $LetterMatrix
##                a
## col_accmock TRUE
## col_iaa     TRUE
## col_iaamock TRUE
## col_salt    TRUE
## col_acc     TRUE
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
#test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"
test
test$info <- strsplit(test$group1, "_")
test$info[[1]][2]
## [1] "accmock"
test$Accession <- "none"
test$Condition<- "none"
test
for(i in 1:nrow(test)){
  test$Accession[i] <- test$info[[i]][1]
     test$Condition[i] <- test$info[[i]][2]
  
}

test2 <- test[,c(4:5,1)]
test2$group1 <- test2$Condition
test2$group2 <- test2$Condition
ICP
ICP2 <- subset(ICP, ICP$K.con..mg.mg.dry.weight < 24 & ICP$K.con..mg.mg.dry.weight > 16)

k_content <- ggplot(data = ICP2, mapping = aes(x = Condition, y = K.con..mg.mg.dry.weight, colour = Condition)) 
k_content <- k_content + geom_boxplot(alpha=0.2) + geom_jitter(width=0.1,alpha=0.2)

#Na_content <- Na_content + facet_grid(~Condition , scales = "free_y")

k_content <- k_content + stat_summary(fun=mean, geom="point", shape=95, size=6, color="black", fill="black")
k_content <- k_content + scale_color_manual(values = c("red","blue","blueviolet","cyan","brown"))
k_content <- k_content + ylab("K content, mg/mg dry weight") + xlab("")+stat_pvalue_manual(test2, label = "Tukey", y.position = 23)
k_content <- k_content + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
k_content <- k_content + rremove("legend")
k_content

aov(K.con..mg.mg.dry.weight ~ All.ID, data = ICP)
## Call:
##    aov(formula = K.con..mg.mg.dry.weight ~ All.ID, data = ICP)
## 
## Terms:
##                   All.ID Residuals
## Sum of Squares   20.6293  720.8996
## Deg. of Freedom        4        29
## 
## Residual standard error: 4.985841
## Estimated effects may be unbalanced
Output <- TukeyHSD(aov(K.con..mg.mg.dry.weight ~ All.ID, data = ICP))
Output
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = K.con..mg.mg.dry.weight ~ All.ID, data = ICP)
## 
## $All.ID
##                                diff       lwr       upr     p adj
## col_accmock-col_acc     -1.29019399 -9.353319  6.772931 0.9898801
## col_iaa-col_acc         -1.24693876 -8.993738  6.499861 0.9896514
## col_iaamock-col_acc      0.13464211 -7.612157  7.881442 0.9999984
## col_salt-col_acc         0.66585882 -7.080941  8.412658 0.9990869
## col_iaa-col_accmock      0.04325523 -8.019869  8.106380 1.0000000
## col_iaamock-col_accmock  1.42483610 -6.638289  9.487961 0.9853117
## col_salt-col_accmock     1.95605281 -6.107072 10.019177 0.9536551
## col_iaamock-col_iaa      1.38158087 -6.365219  9.128380 0.9847998
## col_salt-col_iaa         1.91279759 -5.834002  9.659597 0.9507054
## col_salt-col_iaamock     0.53121671 -7.215583  8.278016 0.9996257
P5 = Output$All.ID[,'p adj']
stat.test<- multcompLetters(P5)
stat.test
## $Letters
## col_accmock     col_iaa col_iaamock    col_salt     col_acc 
##         "a"         "a"         "a"         "a"         "a" 
## 
## $LetterMatrix
##                a
## col_accmock TRUE
## col_iaa     TRUE
## col_iaamock TRUE
## col_salt    TRUE
## col_acc     TRUE
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
#test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"
test
test$info <- strsplit(test$group1, "_")
test$info[[1]][2]
## [1] "accmock"
test$Accession <- "none"
test$Condition<- "none"
test
for(i in 1:nrow(test)){
  test$Accession[i] <- test$info[[i]][1]
     test$Condition[i] <- test$info[[i]][2]
  
}

test2 <- test[,c(4:5,1)]
test2$group1 <- test2$Condition
test2$group2 <- test2$Condition
ICP
# Plot Ttest
ICP2 <- subset(ICP, ICP$K.con..mg.mg.dry.weight < 24 & ICP$K.con..mg.mg.dry.weight > 16)

k_content_col <- ggplot(ICP2, aes(x = Condition, y = K.con..mg.mg.dry.weight, colour = Condition)) + 
  geom_boxplot(alpha = 0.2) + 
  geom_jitter(width = 0.1, alpha = 0.2) + 
  stat_summary(fun = mean, geom = "point", shape = 95, size = 6, color = "black", fill = "black") + 
  scale_color_manual(values = c("red","blue","blueviolet","cyan","brown")) + 
  scale_x_discrete(drop = TRUE) +  # 🔹 Forces ggplot to remove empty categories
  ylab("K content, mg/mg dry weight") + 
  xlab("") + stat_compare_means(label = "p.signif", method = "t.test", ref.group = "accmock")
  theme(axis.text.x = element_text(angle = 90, hjust = 0.9, vjust = 0.5)) + 
  rremove("legend")
## List of 2
##  $ axis.text.x    :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : num 0.9
##   ..$ vjust        : num 0.5
##   ..$ angle        : num 90
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi FALSE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ legend.position: chr "none"
##  - attr(*, "class")= chr [1:2] "theme" "gg"
##  - attr(*, "complete")= logi FALSE
##  - attr(*, "validate")= logi TRUE
# Display plot
k_content_col

aov(Na.K.ratio ~ All.ID, data = ICP)
## Call:
##    aov(formula = Na.K.ratio ~ All.ID, data = ICP)
## 
## Terms:
##                    All.ID Residuals
## Sum of Squares  0.2549755 1.3946201
## Deg. of Freedom         4        29
## 
## Residual standard error: 0.2192951
## Estimated effects may be unbalanced
Output <- TukeyHSD(aov(Na.K.ratio ~ All.ID, data = ICP))
Output
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Na.K.ratio ~ All.ID, data = ICP)
## 
## $All.ID
##                                 diff        lwr       upr     p adj
## col_accmock-col_acc      0.214890665 -0.1397544 0.5695357 0.4144046
## col_iaa-col_acc          0.002898551 -0.3378334 0.3436305 0.9999999
## col_iaamock-col_acc      0.011040288 -0.3296917 0.3517723 0.9999810
## col_salt-col_acc         0.144521158 -0.1962108 0.4852531 0.7326113
## col_iaa-col_accmock     -0.211992114 -0.5666372 0.1426530 0.4279007
## col_iaamock-col_accmock -0.203850376 -0.5584954 0.1507947 0.4667580
## col_salt-col_accmock    -0.070369507 -0.4250146 0.2842756 0.9774589
## col_iaamock-col_iaa      0.008141737 -0.3325902 0.3488737 0.9999944
## col_salt-col_iaa         0.141622607 -0.1991094 0.4823546 0.7467050
## col_salt-col_iaamock     0.133480870 -0.2072511 0.4742128 0.7848598
P7 = Output$All.ID[,'p adj']
stat.test<- multcompLetters(P7)
stat.test
## $Letters
## col_accmock     col_iaa col_iaamock    col_salt     col_acc 
##         "a"         "a"         "a"         "a"         "a" 
## 
## $LetterMatrix
##                a
## col_accmock TRUE
## col_iaa     TRUE
## col_iaamock TRUE
## col_salt    TRUE
## col_acc     TRUE
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
#test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"
test
test$info <- strsplit(test$group1, "_")
test$info[[1]][2]
## [1] "accmock"
test$Accession <- "none"
test$Condition<- "none"
test
for(i in 1:nrow(test)){
  test$Accession[i] <- test$info[[i]][1]
     test$Condition[i] <- test$info[[i]][2]
  
}

test2 <- test[,c(4:5,1)]
test2$group1 <- test2$Condition
test2$group2 <- test2$Condition
ICP
ICP2 <- subset(ICP, ICP$Na.K.ratio < 0.6 & ICP$Na.K.ratio > 0.25)

ratio <- ggplot(data = ICP2, mapping = aes(x = Condition, y = Na.K.ratio, colour = Condition)) 
ratio <- ratio + geom_boxplot(alpha=0.2) + geom_jitter(width=0.1,alpha=0.2)
ratio <- ratio + stat_summary(fun=mean, geom="point", shape=95, size=6, color="black", fill="black")
ratio <- ratio + scale_color_manual(values = c("red","blue","blueviolet","cyan","brown"))
ratio <- ratio + ylab("Na/K ratio") + xlab("")+stat_pvalue_manual(test2, label = "Tukey", y.position = 0.6)
ratio <- ratio + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
ratio <- ratio + rremove("legend")
ratio

# Plot Ttest
ICP2 <- subset(ICP, ICP$Na.K.ratio < 0.6 & ICP$Na.K.ratio > 0.25)

ratio_col <- ggplot(ICP2, aes(x = Condition, y = Na.K.ratio, colour = Condition)) + 
  geom_boxplot(alpha = 0.2) + 
  geom_jitter(width = 0.1, alpha = 0.2) + 
  stat_summary(fun = mean, geom = "point", shape = 95, size = 6, color = "black", fill = "black") + 
  scale_color_manual(values = c("red","blue","blueviolet","cyan","brown")) + 
  scale_x_discrete(drop = TRUE) +  # 🔹 Forces ggplot to remove empty categories
  ylab("Na/K ratio") + 
  xlab("") + stat_compare_means(label = "p.signif", method = "t.test", ref.group = "accmock")
  theme(axis.text.x = element_text(angle = 90, hjust = 0.9, vjust = 0.5)) + 
  rremove("legend")
## List of 2
##  $ axis.text.x    :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : num 0.9
##   ..$ vjust        : num 0.5
##   ..$ angle        : num 90
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi FALSE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ legend.position: chr "none"
##  - attr(*, "class")= chr [1:2] "theme" "gg"
##  - attr(*, "complete")= logi FALSE
##  - attr(*, "validate")= logi TRUE
# Display plot
ratio_col

library(cowplot)
pdf("Figure_all_ICP_AtLeafSpray.pdf", height = 3, width = 12)
plot_grid(Na_content_col, k_content_col, ratio_col, ncol=3,
          align = "hv", labels=c("AUTO"), 
          label_size = 24)
dev.off()
## png 
##   2