This notebook describes data curation and analysis for the results collected for HapMap accessions of Arabidopsis grown on agar plates. The image data and experimental design is coming from Julkowska et al., 2017 paper. The image data was re-analyzed using a new tool for increase in root and shoot projected area and subsequently used in the Genome Wide Association Study (GWAS).

The below analysis was performed by Magdalena Julkowska, when she was still residing at King Abdullah University of Science and Technology as a PostDoc.

Data upload:

Let’s upload the original data - that has been only decoded for the genotype information:

rs <- read.csv("20160318_RS_all.csv", header=T)
head(rs)
rs<- na.omit(rs)
dim(rs)
## [1] 19830     9
unique(rs$Day)
## [1] 2 4 6 8
gens <- unique(rs$ABRC)
trs <- unique(rs$Condition)

vars <- colnames(rs)
cols <- c("experiment", "PlateNumber", "Condition", "Genotype", "ABRC", "Line", "Day", "Shoot", "Root")

Growth Modeling:

root1 <- data.frame(root = paste(rs$experiment, "-", rs$PlateNumber, "-", rs$Condition, "-", rs$ABRC, sep=""))
root1$ACCESSION <- factor(rs$ABRC)
root1$TREATMENT <- factor(rs$Condition)
# Those things CANNOT be a factor - otherwise it isnt going to work! 
root1$Day <-rs$Day
root1$Shoot <- rs$Shoot
root1$Root <- rs$Root
head(root1)
p <- unique(root1$root)

Let’s make data frame to put all the co-eficients for the linear function for each individual line:

how_fitting1a <- data.frame(root = unique(root1$root))

Let’s establish the models:

i=1
A <- root1$root == p[i]
temp <- root1[A,]
head(temp)
# for some models we will need to transform the data - so they will fit into a linear model:
temp$sqShoot <- sqrt(temp$Shoot)
temp$logShoot <- log(temp$Shoot)
temp$sqRoot <- sqrt(temp$Root)
temp$logRoot <- log(temp$Root)

head(temp)
# linear model for shoot
sl <- lm(temp$Shoot ~ temp$Day)
plot (temp$Shoot ~ temp$Day, main=p[i])

head(temp)
sl <- lm(temp$Shoot ~ temp$Day)
plot (temp$Shoot ~ temp$Day, main=p[i])
abline(lm(temp$Shoot ~ temp$Day))

sq <- lm(temp$sqShoot ~ temp$Day)
plot (temp$sqShoot ~ temp$Day, main=p[i])
abline(lm(temp$sqShoot ~ temp$Day))

se <- lm(temp$logShoot ~ temp$Day)
plot (temp$logShoot ~ temp$Day, main=p[i])
abline(lm(temp$logShoot ~ temp$Day))

now - let’s add all the values of the models into the table:

slSTART <- coef(sl)["(Intercept)"]
slGF <- coef(sl)["temp$Day"]
shoot_sqrSTART <- coef(sq)["(Intercept)"]
shoot_sqrGF <- coef(sq)["temp$Day"]
shoot_exSTART <- coef(se)["(Intercept)"]
shoot_exGF <- coef(se)["temp$Day"]

# transform the factors back to their "normal" version
sqSTART <- shoot_sqrSTART^2
sqGF <- shoot_sqrGF^2
seSTART <- log(shoot_exSTART)
seGF <- shoot_exGF

# get the r2 values for each model
sr2l <- summary(sl)$r.squared
sr2q <- summary(sq)$r.squared
sr2e <- summary(se)$r.squared

# make sure everything is numeric (should be but just in case) - actually not very sure how neccessary this part is
as.numeric(slGF)
## [1] 513.5577
as.numeric(slSTART)
## [1] -1123.692
as.numeric(sr2l)
## [1] 0.8712237
as.numeric(sqGF)
## [1] 52.77069
as.numeric(sqSTART)
## [1] 3.076293
as.numeric(sr2q)
## [1] 0.9547821
as.numeric(seGF)
## [1] 0.4676231
as.numeric(seSTART)
## [1] 1.509458
as.numeric(sr2e)
## [1] 0.9840748
# add calculated values from lm to the BIG list with all linear fitted functions
how_fitting1a[i,2] <- slSTART
how_fitting1a[i,3] <- slGF
how_fitting1a[i,4] <- sr2l
how_fitting1a[i,5] <- sqSTART
how_fitting1a[i,6] <- sqGF
how_fitting1a[i,7] <- sr2q
how_fitting1a[i,8] <- seSTART
how_fitting1a[i,9] <- seGF
how_fitting1a[i,10] <- sr2e

how_fitting1a

AND NOW let’s do the same for the root

rl <- lm(temp$Root ~ temp$Day)
plot (temp$Root ~ temp$Day, main=p[i])
abline(lm(temp$Root ~ temp$Day))

rq <- lm(temp$sqRoot ~ temp$Day)
plot (temp$sqRoot ~ temp$Day, main=p[i])
abline(lm(temp$sqRoot ~ temp$Day))

re <- lm(temp$logRoot ~ temp$Day)
plot (temp$logRoot ~ temp$Day, main=p[i])
abline(lm(temp$logRoot ~ temp$Day))

rlSTART <- coef(rl)["(Intercept)"]
rlGF <- coef(rl)["temp$Day"]
root_sqrSTART <- coef(rq)["(Intercept)"]
root_sqrGF <- coef(rq)["temp$Day"]
root_exSTART <- coef(re)["(Intercept)"]
root_exGF <- coef(re)["temp$Day"]

# transform the factors back to their "normal" version
rqSTART <- root_sqrSTART^2
rqGF <- root_sqrGF^2
reSTART <- log(root_exSTART)
reGF <- root_exGF

# get the r2 values for each model
rr2l <- summary(rl)$r.squared
rr2q <- summary(rq)$r.squared
rr2e <- summary(re)$r.squared

# make sure everything is numeric (should be but just in case) - actually not very sure how neccessary this part is
as.numeric(rlGF)
## [1] 631.1635
as.numeric(rlSTART)
## [1] -1373.462
as.numeric(rr2l)
## [1] 0.8535071
as.numeric(rqGF)
## [1] 62.66163
as.numeric(rqSTART)
## [1] 1.166595
as.numeric(rr2q)
## [1] 0.9521454
as.numeric(reGF)
## [1] 0.451193
as.numeric(reSTART)
## [1] 1.575757
as.numeric(rr2e)
## [1] 0.9875454
# add calculated values from lm to the BIG list with all linear fitted functions
how_fitting1a[i,11] <- rlSTART
how_fitting1a[i,12] <- rlGF
how_fitting1a[i,13] <- rr2l
how_fitting1a[i,14] <- rqSTART
how_fitting1a[i,15] <- rqGF
how_fitting1a[i,16] <- rr2q
how_fitting1a[i,17] <- reSTART
how_fitting1a[i,18] <- reGF
how_fitting1a[i,19] <- rr2e

# I also need to have raw data somewhere to do correlations later on in the file with all the growth factor, what are they based on & co.
temp
B <- temp$Day ==2
Shoot_d2 <- mean(temp[B,5])
Root_d2 <- mean(temp[B,6])

C <- temp$Day ==4
Shoot_d4 <- mean(temp[C,5])
Root_d4 <- mean(temp[C,6])

D <- temp$Day ==6
Shoot_d6 <- mean(temp[D,5])
Root_d6 <- mean(temp[D,6])

E <- temp$Day ==8
Shoot_d8 <- mean(temp[E,5])
Root_d8 <- mean(temp[E,6])

how_fitting1a[i,20] <- Shoot_d2
how_fitting1a[i,21] <- Shoot_d4
how_fitting1a[i,22] <- Shoot_d6
how_fitting1a[i,23] <- Shoot_d8
how_fitting1a[i,24] <- Root_d2
how_fitting1a[i,25] <- Root_d4
how_fitting1a[i,26] <- Root_d6
how_fitting1a[i,27] <- Root_d8
how_fitting1a

Before we loop anything - let’s replace the column names as they should be:

colnames(how_fitting1a)[which(names(how_fitting1a) == "V2")] <- "slSTART"
colnames(how_fitting1a)[which(names(how_fitting1a) == "V3")] <- "slGF"
colnames(how_fitting1a)[which(names(how_fitting1a) == "V4")] <- "slr2"

colnames(how_fitting1a)[which(names(how_fitting1a) == "V5")] <- "sqSTART"
colnames(how_fitting1a)[which(names(how_fitting1a) == "V6")] <- "sqGF"
colnames(how_fitting1a)[which(names(how_fitting1a) == "V7")] <- "sqr2"

colnames(how_fitting1a)[which(names(how_fitting1a) == "V8")] <- "seSTART"
colnames(how_fitting1a)[which(names(how_fitting1a) == "V9")] <- "seGF"
colnames(how_fitting1a)[which(names(how_fitting1a) == "V10")] <- "ser2"

colnames(how_fitting1a)[which(names(how_fitting1a) == "V11")] <- "rlSTART"
colnames(how_fitting1a)[which(names(how_fitting1a) == "V12")] <- "rlGF"
colnames(how_fitting1a)[which(names(how_fitting1a) == "V13")] <- "rlr2"

colnames(how_fitting1a)[which(names(how_fitting1a) == "V14")] <- "rqSTART"
colnames(how_fitting1a)[which(names(how_fitting1a) == "V15")] <- "rqGF"
colnames(how_fitting1a)[which(names(how_fitting1a) == "V16")] <- "rqr2"

colnames(how_fitting1a)[which(names(how_fitting1a) == "V17")] <- "reSTART"
colnames(how_fitting1a)[which(names(how_fitting1a) == "V18")] <- "reGF"
colnames(how_fitting1a)[which(names(how_fitting1a) == "V19")] <- "rer2"

colnames(how_fitting1a)[which(names(how_fitting1a) == "V20")] <- "Shoot_d2"
colnames(how_fitting1a)[which(names(how_fitting1a) == "V21")] <- "Shoot_d4"
colnames(how_fitting1a)[which(names(how_fitting1a) == "V22")] <- "Shoot_d6"
colnames(how_fitting1a)[which(names(how_fitting1a) == "V23")] <- "Shoot_d8"

colnames(how_fitting1a)[which(names(how_fitting1a) == "V24")] <- "Root_d2"
colnames(how_fitting1a)[which(names(how_fitting1a) == "V25")] <- "Root_d4"
colnames(how_fitting1a)[which(names(how_fitting1a) == "V26")] <- "Root_d6"
colnames(how_fitting1a)[which(names(how_fitting1a) == "V27")] <- "Root_d8"
how_fitting1a

Awesome! Now - let’s loop it for all of the imaged plants:

for(i in 1:2497){
A <- root1$root == p[i]
temp <- root1[A,]

temp$sqShoot <- sqrt(temp$Shoot)
temp$logShoot <- log(temp$Shoot)
temp$sqRoot <- sqrt(temp$Root)
temp$logRoot <- log(temp$Root)

head(temp)
sl <- lm(temp$Shoot ~ temp$Day)
sq <- lm(temp$sqShoot ~ temp$Day)
se <- lm(temp$logShoot ~ temp$Day)

slSTART <- coef(sl)["(Intercept)"]
slGF <- coef(sl)["temp$Day"]
shoot_sqrSTART <- coef(sq)["(Intercept)"]
shoot_sqrGF <- coef(sq)["temp$Day"]
shoot_exSTART <- coef(se)["(Intercept)"]
shoot_exGF <- coef(se)["temp$Day"]

# transform the factors back to their "normal" version
sqSTART <- shoot_sqrSTART^2
sqGF <- shoot_sqrGF^2
seSTART <- log(shoot_exSTART)
seGF <- shoot_exGF

# get the r2 values for each model
sr2l <- summary(sl)$r.squared
sr2q <- summary(sq)$r.squared
sr2e <- summary(se)$r.squared


# add calculated values from lm to the BIG list with all linear fitted functions
how_fitting1a[i,2] <- slSTART
how_fitting1a[i,3] <- slGF
how_fitting1a[i,4] <- sr2l
how_fitting1a[i,5] <- sqSTART
how_fitting1a[i,6] <- sqGF
how_fitting1a[i,7] <- sr2q
how_fitting1a[i,8] <- seSTART
how_fitting1a[i,9] <- seGF
how_fitting1a[i,10] <- sr2e

# AND NOW let's do the same for the root
rl <- lm(temp$Root ~ temp$Day)
rq <- lm(temp$sqRoot ~ temp$Day)
re <- lm(temp$logRoot ~ temp$Day)

rlSTART <- coef(rl)["(Intercept)"]
rlGF <- coef(rl)["temp$Day"]
root_sqrSTART <- coef(rq)["(Intercept)"]
root_sqrGF <- coef(rq)["temp$Day"]
root_exSTART <- coef(re)["(Intercept)"]
root_exGF <- coef(re)["temp$Day"]

# transform the factors back to their "normal" version
rqSTART <- root_sqrSTART^2
rqGF <- root_sqrGF^2
reSTART <- log(root_exSTART)
reGF <- root_exGF

# get the r2 values for each model
rr2l <- summary(rl)$r.squared
rr2q <- summary(rq)$r.squared
rr2e <- summary(re)$r.squared

# add calculated values from lm to the BIG list with all linear fitted functions
how_fitting1a[i,11] <- rlSTART
how_fitting1a[i,12] <- rlGF
how_fitting1a[i,13] <- rr2l
how_fitting1a[i,14] <- rqSTART
how_fitting1a[i,15] <- rqGF
how_fitting1a[i,16] <- rr2q
how_fitting1a[i,17] <- reSTART
how_fitting1a[i,18] <- reGF
how_fitting1a[i,19] <- rr2e

# I also need to have raw data somewhere to do correlations later on in the file with all the growth factor, what are they based on & co.
temp
B <- temp$Day ==2
Shoot_d2 <- mean(temp[B,5])
Root_d2 <- mean(temp[B,6])

C <- temp$Day ==4
Shoot_d4 <- mean(temp[C,5])
Root_d4 <- mean(temp[C,6])

D <- temp$Day ==6
Shoot_d6 <- mean(temp[D,5])
Root_d6 <- mean(temp[D,6])

E <- temp$Day ==8
Shoot_d8 <- mean(temp[E,5])
Root_d8 <- mean(temp[E,6])

how_fitting1a[i,20] <- Shoot_d2
how_fitting1a[i,21] <- Shoot_d4
how_fitting1a[i,22] <- Shoot_d6
how_fitting1a[i,23] <- Shoot_d8
how_fitting1a[i,24] <- Root_d2
how_fitting1a[i,25] <- Root_d4
how_fitting1a[i,26] <- Root_d6
how_fitting1a[i,27] <- Root_d8
}
how_fitting1a

Cool - and now - we need to dissect the “root” information, which contains all of the identifying meta-data for each plant.

how_fitting1a$experiment <- "none"
how_fitting1a$plate <- "none"
how_fitting1a$condition <- "none"
how_fitting1a$CS <- "none"

dim(how_fitting1a)
## [1] 2497   31
for(i in 1:2497){
  how_fitting1a[i,28] <- strsplit(how_fitting1a[i,1], "-")[[1]][1]
  how_fitting1a[i,29] <- strsplit(how_fitting1a[i,1], "-")[[1]][2]
  how_fitting1a[i,30] <- strsplit(how_fitting1a[i,1], "-")[[1]][3]
  how_fitting1a[i,31] <- strsplit(how_fitting1a[i,1], "-")[[1]][4]
}

head(how_fitting1a)
hf1 <- how_fitting1a[,c(28:31, 2:27)]
head(hf1)
dim(hf1)
## [1] 2497   30
hf1 <- na.omit(hf1)
dim(hf1)
## [1] 2493   30

Which model is fitting my data the best?

library(reshape2)
library(ggpubr)
## Loading required package: ggplot2
library(ggplot2)
library(ggsci)

colnames(hf1)
##  [1] "experiment" "plate"      "condition"  "CS"         "slSTART"   
##  [6] "slGF"       "slr2"       "sqSTART"    "sqGF"       "sqr2"      
## [11] "seSTART"    "seGF"       "ser2"       "rlSTART"    "rlGF"      
## [16] "rlr2"       "rqSTART"    "rqGF"       "rqr2"       "reSTART"   
## [21] "reGF"       "rer2"       "Shoot_d2"   "Shoot_d4"   "Shoot_d6"  
## [26] "Shoot_d8"   "Root_d2"    "Root_d4"    "Root_d6"    "Root_d8"
which_model <- hf1[,c(7,10,13,16,19,22)]
melt_model <- melt(which_model)
## No id variables; using all as measure variables
melt_model
model <- ggerrorplot(melt_model, y = "value", x = "variable", fill="variable", color="variable", 
                        desc_stat = "mean_sd", add = "jitter", 
                        add.params = list(color = "darkgray"),
                        xlab="model", ylab="Regression coefficient") 
model <- model + scale_color_jco()
model

as expected exponential model fits best - further follow on with only exponentially fitted data

how about we remove all of the other GR calculations - so it doesn’t confuse us anymore:

hf_clean <- hf1[,c(1:4, 23:30, 12:13, 21:22)]
hf_clean

also - let’s remove the samples for which regression coefficient was really low - like below 0.7:

dim(hf_clean)
## [1] 2493   16
hf_clean <- subset(hf_clean, hf_clean$rer2 > 0.8)
hf_clean <- subset(hf_clean, hf_clean$ser2 > 0.8)
dim(hf_clean)
## [1] 2045   16

Data cutation

DO I HAVE DIFFRERENCE BETWEEN EXPERIMENTAL BATCHES?

check if I have significant differences between my experiments first maybe

head(hf_clean)
hf_clean$condition <- factor(hf_clean$condition, levels = c(0, 75, 125))

Batch_rootGR <- ggerrorplot(hf_clean, y = "reGF", x = "experiment", fill="experiment", color="experiment", 
                        desc_stat = "mean_sd", add = "jitter", facet.by = "condition",
                        add.params = list(color = "darkgray"),
                        xlab="experiment", ylab="growth rate of the root (pix / day)") 
Batch_rootGR <- Batch_rootGR + rremove("legend") + stat_compare_means(method="aov")
Batch_rootGR <- Batch_rootGR + scale_color_aaas()
Batch_rootGR

Batch_shootGR <- ggerrorplot(hf_clean, y = "seGF", x = "experiment", fill="experiment", color="experiment", 
                        desc_stat = "mean_sd", add = "jitter",  facet.by = "condition",
                        add.params = list(color = "darkgray"),
                        xlab="experiment", ylab="growth rate of the shoot (pix / day)") 
Batch_shootGR <- Batch_shootGR + rremove("legend") + stat_compare_means(method="aov")
Batch_shootGR <- Batch_shootGR + scale_color_aaas()
Batch_shootGR

OK - so there doesnt seem to be a big effect of the experimental batch on the root and shoot growth rate. Great!

WHAT IS THE NATURAL VARIATION IN GROWTH RATES OF INDIVIDUAL ACCESSIONS?

Now - we would like to inspect if there are any accessions that show THREMENDOUS variation within their genotype:

Geno_rootGR <- ggerrorplot(hf_clean, y = "reGF", x = "CS", fill="CS", color="CS", 
                        desc_stat = "mean_sd", add = "jitter",  facet.by = "condition",
                        add.params = list(color = "darkgray"),
                        xlab="genotype", ylab="growth rate of the root (pix / day)") 
Geno_rootGR <- Geno_rootGR + rremove("legend") + stat_compare_means(method="aov")
Geno_rootGR <- Geno_rootGR + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank())
Geno_rootGR
## Warning: Removed 35 rows containing missing values (geom_segment).
## Warning: Removed 38 rows containing missing values (geom_segment).
## Warning: Removed 121 rows containing missing values (geom_segment).

Geno_shootGR <- ggerrorplot(hf_clean, y = "seGF", x = "CS", fill="CS", color="CS", 
                        desc_stat = "mean_sd", add = "jitter", facet.by = "condition",
                        add.params = list(color = "darkgray"),
                        xlab="genotype", ylab="growth rate of the shoot (pix / day)") 
Geno_shootGR <- Geno_shootGR + rremove("legend") + stat_compare_means(method="aov")
Geno_shootGR <- Geno_shootGR + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank())
Geno_shootGR
## Warning: Removed 35 rows containing missing values (geom_segment).
## Warning: Removed 38 rows containing missing values (geom_segment).
## Warning: Removed 121 rows containing missing values (geom_segment).

Let’s remove above 0.5 values for 125 root growth and 0.4 shoot values for 125 mM

hf_clean2 <- subset(hf_clean, !(hf_clean$condition == "125" & hf_clean$reGF > 0.5))
hf_clean2 <- subset(hf_clean2, !(hf_clean2$condition == "125" & hf_clean2$seGF > 0.4))
hf_clean2 <- subset(hf_clean2, hf_clean2$Root_d8 < 20000)
Geno_rootGR2 <- ggerrorplot(hf_clean2, y = "reGF", x = "CS", fill="CS", color="CS", 
                        desc_stat = "mean_sd", add = "jitter",  facet.by = "condition",
                        add.params = list(color = "darkgray"),
                        xlab="genotype", ylab="growth rate of the root (pix / day)") 
Geno_rootGR2 <- Geno_rootGR2 + rremove("legend") + stat_compare_means(method="aov")
Geno_rootGR2 <- Geno_rootGR2 + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank())
Geno_rootGR2
## Warning: Removed 36 rows containing missing values (geom_segment).
## Warning: Removed 38 rows containing missing values (geom_segment).
## Warning: Removed 124 rows containing missing values (geom_segment).

Geno_shootGR2 <- ggerrorplot(hf_clean2, y = "seGF", x = "CS", fill="CS", color="CS", 
                        desc_stat = "mean_sd", add = "jitter", facet.by = "condition",
                        add.params = list(color = "darkgray"),
                        xlab="genotype", ylab="growth rate of the shoot (pix / day)") 
Geno_shootGR2 <- Geno_shootGR2 + rremove("legend") + stat_compare_means(method="aov")
Geno_shootGR2 <- Geno_shootGR2 + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank())
Geno_shootGR2
## Warning: Removed 36 rows containing missing values (geom_segment).
## Warning: Removed 38 rows containing missing values (geom_segment).
## Warning: Removed 124 rows containing missing values (geom_segment).

dim(hf_clean2)
## [1] 2041   16
head(hf_clean)
length(unique(hf_clean2$CS))
## [1] 348

great - this looks good. Moving on.

Data visualization throughout time:

Let’s explore the data and how it changes throughout time and with salt stress.

However, in order to do that - we would first need to calculate genotype specific mean for each accession and condition:

library(doBy)
colnames(hf_clean2)
##  [1] "experiment" "plate"      "condition"  "CS"         "Shoot_d2"  
##  [6] "Shoot_d4"   "Shoot_d6"   "Shoot_d8"   "Root_d2"    "Root_d4"   
## [11] "Root_d6"    "Root_d8"    "seGF"       "ser2"       "reGF"      
## [16] "rer2"
sum_data <- summaryBy(Shoot_d2 + Shoot_d4 + Shoot_d6 + Shoot_d8 + Root_d2  + Root_d4 + Root_d6  + Root_d8 + seGF + reGF ~ condition + CS, hf_clean2)
colnames(sum_data) <- gsub(".mean", "", colnames(sum_data)) 
head(sum_data)

To do that - I first would like to re-shape the data to plot it over time

time_data <- sum_data
time_data$plant.id <- paste(time_data$condition, time_data$CS, sep="_")
time_data
time_data2 <- time_data[,c(13,1,3:12)]
time_data2
time_melt <- melt(time_data2, id.vars = c("plant.id", "condition"))
time_melt

now - I need to split the time from tissue in the variable column:

time_melt$tissue <- "none"
time_melt$time <- "none"

for(i in 1:9690){
  time_melt[i,5] <- strsplit(as.character(time_melt[i,3]), "_")[[1]][1]
  time_melt[i,6] <- strsplit(as.character(time_melt[i,3]), "_")[[1]][2]
}
time_melt
unique(time_melt$tissue)
## [1] "Shoot" "Root"  "seGF"  "reGF"
unique(time_melt$time)
## [1] "d2" "d4" "d6" "d8" NA
time_melt$time <- gsub("d2", "2", time_melt$time)
time_melt$time <- gsub("d4", "4", time_melt$time)
time_melt$time <- gsub("d6", "6", time_melt$time)
time_melt$time <- gsub("d8", "8", time_melt$time)

time_melt$time <- as.numeric(as.character(time_melt$time))
unique(time_melt$time)
## [1]  2  4  6  8 NA
time_melt <- na.omit(time_melt)

OK - since we got too many comparisons to plot into one graph - we should be looking at root and shoot separately:

time_root <- subset(time_melt, time_melt$tissue == "Root")
time_shoot <- subset(time_melt, time_melt$tissue == "Shoot")

time_root_0 <- subset(time_root, time_root$condition == "0")
time_root_0$batch <- "75 mM NaCl"
time_root_75 <- subset(time_root, time_root$condition == "75")
time_root_75$batch <- "75 mM NaCl"
time_root2.1 <- rbind(time_root_0, time_root_75)
time_root_125 <- subset(time_root, time_root$condition == "125")
time_root_125$batch <- "125 mM NaCl"
time_root_0$batch <- "125 mM NaCl"
time_root2.2 <- rbind(time_root_0, time_root_125)
time_root2 <- rbind(time_root2.1, time_root2.2)

time_shoot_0 <- subset(time_shoot, time_shoot$condition == "0")
time_shoot_0$batch <- "75 mM NaCl"
time_shoot_75 <- subset(time_shoot, time_shoot$condition == "75")
time_shoot_75$batch <- "75 mM NaCl"
time_shoot2.1 <- rbind(time_shoot_0, time_shoot_75)
time_shoot_125 <- subset(time_shoot, time_shoot$condition == "125")
time_shoot_125$batch <- "125 mM NaCl"
time_shoot_0$batch <- "125 mM NaCl"
time_shoot2.2 <- rbind(time_shoot_0, time_shoot_125)
time_shoot2 <- rbind(time_shoot2.1, time_shoot2.2)

time_root2$batch <- factor(time_root2$batch, levels = c("75 mM NaCl", "125 mM NaCl"))
time_shoot2$batch <- factor(time_shoot2$batch, levels = c("75 mM NaCl", "125 mM NaCl"))

unique(time_root2$condition)
## [1] 0   75  125
## Levels: 0 75 125
time_root2$condition <- as.character(time_root2$condition)
time_root2$condition <- gsub("0", "0 mM NaCl", time_root2$condition)
time_root2$condition <- gsub("75", "75 mM NaCl", time_root2$condition)
time_root2$condition <- gsub("125", "125 mM NaCl", time_root2$condition)

time_shoot2$condition <- as.character(time_shoot2$condition)
time_shoot2$condition <- gsub("0", "0 mM NaCl", time_shoot2$condition)
time_shoot2$condition <- gsub("75", "75 mM NaCl", time_shoot2$condition)
time_shoot2$condition <- gsub("125", "125 mM NaCl", time_shoot2$condition)

time_root2$condition <- factor(time_root2$condition, levels = c("0 mM NaCl", "75 mM NaCl", "125 mM NaCl"))
time_shoot2$condition <- factor(time_shoot2$condition, levels = c("0 mM NaCl", "75 mM NaCl", "125 mM NaCl"))

Now - let’s get plotting

Root_lgraph <- ggplot(data=time_root2, aes(x= time, y=value, group = plant.id, color = condition)) 
Root_lgraph <- Root_lgraph + geom_line(alpha = 0.1) + facet_grid(~ batch)
Root_lgraph <- Root_lgraph + stat_summary(fun.data = mean_se, geom="ribbon", linetype=0, aes(group= condition), alpha=0.3)
Root_lgraph <- Root_lgraph + stat_summary(fun=mean, aes(group= condition),  size=0.7, geom="line", linetype = "dashed")
Root_lgraph <- Root_lgraph + stat_compare_means(aes(group = condition), label = "p.signif", method = "aov", hide.ns = T)
Root_lgraph <- Root_lgraph + ylab("Root size (pixels)") + xlab("Days after stress") + scale_color_aaas() + theme_bw()
Root_lgraph <- Root_lgraph + guides(color=guide_legend(title="Treatment")) 
Root_lgraph

Shoot_lgraph <- ggplot(data=time_shoot2, aes(x= time, y=value, group = plant.id, color = condition)) 
Shoot_lgraph <- Shoot_lgraph + geom_line(alpha = 0.1) + facet_grid(~ batch)
Shoot_lgraph <- Shoot_lgraph + stat_summary(fun.data = mean_se, geom="ribbon", linetype=0, aes(group= condition), alpha=0.3)
Shoot_lgraph <- Shoot_lgraph + stat_summary(fun=mean, aes(group= condition),  size=0.7, geom="line", linetype = "dashed")
Shoot_lgraph <- Shoot_lgraph + stat_compare_means(aes(group = condition), label = "p.signif", method = "aov", hide.ns = T)
Shoot_lgraph <- Shoot_lgraph + ylab("Shoot size (pixels)") + xlab("Days after stress") + scale_color_aaas() + theme_bw()
Shoot_lgraph <- Shoot_lgraph + guides(color=guide_legend(title="Treatment"))
Shoot_lgraph

library(cowplot)
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggpubr':
## 
##     get_legend
plot_increase <- plot_grid(Shoot_lgraph, Root_lgraph, labels = "AUTO", ncol=1)
plot_increase

pdf("Fig.Supp.Increase.in.shoot.and.root.size.over.time.genotype.mean.pdf", width = 7, height = 8)
plot(plot_increase)
dev.off()
## quartz_off_screen 
##                 2

Growth rate data visualization

OK - so now its time to plot the histograms to visualize how the growth rate changes across the organs.

Before we go anywhere tho - I first need to isolate the growth rate data from the sum_data:

GR_data <- sum_data[,c(1:2, 11:12)]
GR_data
GR_data$RpS_GR <- GR_data$reGF / GR_data$seGF

shoot_GR_histo <- gghistogram(GR_data, x = "seGF", add = "mean", fill = "condition") + scale_color_aaas() + scale_fill_aaas() 
## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.
root_GR_histo <- gghistogram(GR_data, x = "reGF", add = "mean", fill = "condition") + scale_color_aaas() + scale_fill_aaas() 
## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.
rootpshoot_GR_histo <- gghistogram(GR_data, x = "RpS_GR", add = "mean", fill = "condition") + scale_color_aaas() + scale_fill_aaas() 
## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.
shoot_GR_histo

root_GR_histo

rootpshoot_GR_histo

shootGR2 <- ggerrorplot(GR_data, y = "seGF", , x = "condition", fill="condition", color="condition", 
                        desc_stat = "mean_sd", add = "jitter",
                        add.params = list(color = "darkgray"),
                        xlab="Condition", ylab="growth rate of the shoot (pix / day)") + scale_y_continuous(position = "right")
shootGR2 <- shootGR2 + rremove("legend") + stat_compare_means(method="aov", ref.group = "0") + coord_flip()
shootGR2 <- shootGR2 + scale_color_aaas() + scale_fill_aaas() 
shootGR2
## Warning: Computation failed in `stat_compare_means()`:
## Problem with `filter()` input `..1`.
## ℹ Input `..1` is `group1 == ref.group | group2 == ref.group`.
## x object 'group2' not found

rootGR2 <- ggerrorplot(GR_data, y = "reGF", , x = "condition", fill="condition", color="condition", 
                        desc_stat = "mean_sd", add = "jitter",
                        add.params = list(color = "darkgray"),
                        xlab="Condition", ylab="growth rate of the root (pix / day)") + scale_y_continuous(position = "right")
rootGR2 <- rootGR2 + rremove("legend") + stat_compare_means(method="aov", ref.group = "0") + coord_flip() 
rootGR2 <- rootGR2 + scale_color_aaas() + scale_fill_aaas() 
rootGR2
## Warning: Computation failed in `stat_compare_means()`:
## Problem with `filter()` input `..1`.
## ℹ Input `..1` is `group1 == ref.group | group2 == ref.group`.
## x object 'group2' not found

RpS.GR2 <- ggerrorplot(GR_data, y = "RpS_GR", , x = "condition", fill="condition", color="condition", 
                        desc_stat = "mean_sd", add = "jitter",
                        add.params = list(color = "darkgray"),
                        xlab="Condition", ylab="Root GR / Shoot GR (ratio)") + scale_y_continuous(position = "right")
RpS.GR2 <- RpS.GR2 + rremove("legend") + stat_compare_means(method="aov", ref.group = "0") + coord_flip() 
RpS.GR2 <- RpS.GR2 + scale_color_aaas() + scale_fill_aaas() 
RpS.GR2
## Warning: Computation failed in `stat_compare_means()`:
## Problem with `filter()` input `..1`.
## ℹ Input `..1` is `group1 == ref.group | group2 == ref.group`.
## x object 'group2' not found

Cool - now let’s fuse all of these into one figure:

shoot_part <- plot_grid(shoot_GR_histo, shootGR2, ncol = 1, rel_heights = c(1, 0.5))
## Warning: Computation failed in `stat_compare_means()`:
## Problem with `filter()` input `..1`.
## ℹ Input `..1` is `group1 == ref.group | group2 == ref.group`.
## x object 'group2' not found
root_part <- plot_grid(root_GR_histo, rootGR2, ncol = 1, rel_heights = c(1, 0.5))
## Warning: Computation failed in `stat_compare_means()`:
## Problem with `filter()` input `..1`.
## ℹ Input `..1` is `group1 == ref.group | group2 == ref.group`.
## x object 'group2' not found
RpS_part <- plot_grid(rootpshoot_GR_histo, RpS.GR2, ncol = 1, rel_heights = c(1, 0.5))
## Warning: Computation failed in `stat_compare_means()`:
## Problem with `filter()` input `..1`.
## ℹ Input `..1` is `group1 == ref.group | group2 == ref.group`.
## x object 'group2' not found
shoot_part

root_part

RpS_part

final_GR <- plot_grid(shoot_part, root_part, RpS_part, ncol=3, labels = "AUTO")
final_GR

pdf("Fig.Root_and_shoot_growth_rate_histo_and_scatter_plots.pdf", width = 10, height = 7)
plot(final_GR)
dev.off()
## quartz_off_screen 
##                 2

Stress tolerance index calculations for GR

Now - let’s calculate the Stress Tolerance Index (STI) for both root, shoot and root/shoot:

GR_data_0 <- subset(GR_data, GR_data$condition == "0")
GR_data_0
GR_data_0 <- GR_data_0[,c(2:5)]
colnames(GR_data_0)[2] <- "seGF.0"
colnames(GR_data_0)[3] <- "reGF.0"
colnames(GR_data_0)[4] <- "RpS.GR.0"

GR_data2 <- merge(GR_data, GR_data_0, by = "CS", all=TRUE)
dim(GR_data)
## [1] 969   5
dim(GR_data2)
## [1] 969   8
GR_data2$STI.sGR <- GR_data2$seGF / GR_data2$seGF.0
GR_data2$STI.rGR <- GR_data2$reGF / GR_data2$reGF.0
GR_data2$STI.RpS.GR <- GR_data2$RpS_GR / GR_data2$RpS.GR.0
GR_data2 <- subset(GR_data2, GR_data2$condition != "0")
GR_data2
sGR.STI <- ggerrorplot(GR_data2, y = "STI.sGR", , x = "condition", fill="condition", color="condition", 
                        desc_stat = "mean_sd", add = "jitter",
                        add.params = list(color = "darkgray"),
                        xlab="Condition", ylab="Shoot GR at salt / control (ratio)")
sGR.STI <- sGR.STI + rremove("legend") #+ stat_compare_means(method="aov")
sGR.STI <- sGR.STI + scale_color_futurama() + scale_fill_futurama() 
sGR.STI
## Warning: Removed 6 rows containing non-finite values (stat_summary).
## Warning: Removed 6 rows containing missing values (geom_point).

rGR.STI <- ggerrorplot(GR_data2, y = "STI.rGR", , x = "condition", fill="condition", color="condition", 
                        desc_stat = "mean_sd", add = "jitter",
                        add.params = list(color = "darkgray"),
                        xlab="Condition", ylab="Root GR at salt / control (ratio)")
rGR.STI <- rGR.STI + rremove("legend") #+ stat_compare_means(method="aov")
rGR.STI <- rGR.STI + scale_color_futurama() + scale_fill_futurama() 
rGR.STI
## Warning: Removed 6 rows containing non-finite values (stat_summary).
## Warning: Removed 6 rows containing missing values (geom_point).

rpsGR.STI <- ggerrorplot(GR_data2, y = "STI.RpS.GR", , x = "condition", fill="condition", color="condition", 
                        desc_stat = "mean_sd", add = "jitter",
                        add.params = list(color = "darkgray"),
                        xlab="Condition", ylab="Root/Shoot GR at salt / control (ratio)")
rpsGR.STI <- rpsGR.STI + rremove("legend") #+ stat_compare_means(method="aov")
rpsGR.STI <- rpsGR.STI + scale_color_futurama() + scale_fill_futurama() 
rpsGR.STI
## Warning: Removed 6 rows containing non-finite values (stat_summary).
## Warning: Removed 6 rows containing missing values (geom_point).

OK - let’s save all of the relative changes in root and shoot growth as a separate figure:

pdf("Fig.STI_Root.and.shoot.growth.pdf", width = 10, height = 5)
plot_grid(sGR.STI, rGR.STI, rpsGR.STI, ncol = 3)
## Warning: Removed 6 rows containing non-finite values (stat_summary).
## Warning: Removed 6 rows containing missing values (geom_point).
## Warning: Removed 6 rows containing non-finite values (stat_summary).
## Warning: Removed 6 rows containing missing values (geom_point).
## Warning: Removed 6 rows containing non-finite values (stat_summary).
## Warning: Removed 6 rows containing missing values (geom_point).
dev.off()
## quartz_off_screen 
##                 2

Correlations

OK - now we got all of the raw data ready - let’s see how individual traits are correlated.

First - lets make a correlation matrix. For this - all of the traits that we have extracted thus far should be in separate columns for each conditions:

head(sum_data)
colnames(sum_data)
##  [1] "condition" "CS"        "Shoot_d2"  "Shoot_d4"  "Shoot_d6"  "Shoot_d8" 
##  [7] "Root_d2"   "Root_d4"   "Root_d6"   "Root_d8"   "seGF"      "reGF"
sum_data$Total_d2 <- sum_data$Root_d2 + sum_data$Shoot_d2
sum_data$Total_d4 <- sum_data$Root_d4 + sum_data$Shoot_d4
sum_data$Total_d6 <- sum_data$Root_d6 + sum_data$Shoot_d6
sum_data$Total_d8 <- sum_data$Root_d8 + sum_data$Shoot_d8
sum_data$RpS_d2 <- sum_data$Root_d2 / sum_data$Shoot_d2
sum_data$RpS_d4 <- sum_data$Root_d4 / sum_data$Shoot_d4
sum_data$RpS_d6 <- sum_data$Root_d6 / sum_data$Shoot_d6
sum_data$RpS_d8 <- sum_data$Root_d8 / sum_data$Shoot_d8

sum_data$RpT_d2 <- sum_data$Root_d2 / sum_data$Total_d2
sum_data$RpT_d4 <- sum_data$Root_d4 / sum_data$Total_d4
sum_data$RpT_d6 <- sum_data$Root_d6 / sum_data$Total_d6
sum_data$RpT_d8 <- sum_data$Root_d8 / sum_data$Total_d8

colnames(sum_data)[11] <- "Shoot.GR"
colnames(sum_data)[12] <- "Root.GR"
sum_data$RpS.GR <- sum_data$Root.GR / sum_data$Shoot.GR

colnames(sum_data)
##  [1] "condition" "CS"        "Shoot_d2"  "Shoot_d4"  "Shoot_d6"  "Shoot_d8" 
##  [7] "Root_d2"   "Root_d4"   "Root_d6"   "Root_d8"   "Shoot.GR"  "Root.GR"  
## [13] "Total_d2"  "Total_d4"  "Total_d6"  "Total_d8"  "RpS_d2"    "RpS_d4"   
## [19] "RpS_d6"    "RpS_d8"    "RpT_d2"    "RpT_d4"    "RpT_d6"    "RpT_d8"   
## [25] "RpS.GR"
sum_data <- sum_data[,c(1:10,13:24,11:12,25)]
sum_data

OK - now I want to have all of the conditions in separate collumns too:

sum_data_0 <- subset(sum_data, sum_data$condition == "0")
sum_data_75 <- subset(sum_data, sum_data$condition == "75")
sum_data_125 <- subset(sum_data, sum_data$condition == "125")
head(sum_data_0)
sum_data_0 <- sum_data_0[,2:25]
sum_data_75 <- sum_data_75[,2:25]
sum_data_125 <- sum_data_125[,2:25]
colnames(sum_data_0)[2:24] <- paste(colnames(sum_data_0)[2:24], ".0mM.NaCl", sep="")
colnames(sum_data_75)[2:24] <- paste(colnames(sum_data_75)[2:24], ".75mM.NaCl", sep="")
colnames(sum_data_125)[2:24] <- paste(colnames(sum_data_125)[2:24], ".125mM.NaCl", sep="")
sum_data_0
sum_data_75
sum_data_125

ok - now let’s bind them all together into one dataset:

sum_data_all <- merge(sum_data_0, sum_data_75, by = "CS")
sum_data_all <- merge(sum_data_all, sum_data_125, by = "CS")
sum_data_all

Now - let’s calculate STI for Root and Shoot GR and for Total Area:

colnames(sum_data_all)
##  [1] "CS"                  "Shoot_d2.0mM.NaCl"   "Shoot_d4.0mM.NaCl"  
##  [4] "Shoot_d6.0mM.NaCl"   "Shoot_d8.0mM.NaCl"   "Root_d2.0mM.NaCl"   
##  [7] "Root_d4.0mM.NaCl"    "Root_d6.0mM.NaCl"    "Root_d8.0mM.NaCl"   
## [10] "Total_d2.0mM.NaCl"   "Total_d4.0mM.NaCl"   "Total_d6.0mM.NaCl"  
## [13] "Total_d8.0mM.NaCl"   "RpS_d2.0mM.NaCl"     "RpS_d4.0mM.NaCl"    
## [16] "RpS_d6.0mM.NaCl"     "RpS_d8.0mM.NaCl"     "RpT_d2.0mM.NaCl"    
## [19] "RpT_d4.0mM.NaCl"     "RpT_d6.0mM.NaCl"     "RpT_d8.0mM.NaCl"    
## [22] "Shoot.GR.0mM.NaCl"   "Root.GR.0mM.NaCl"    "RpS.GR.0mM.NaCl"    
## [25] "Shoot_d2.75mM.NaCl"  "Shoot_d4.75mM.NaCl"  "Shoot_d6.75mM.NaCl" 
## [28] "Shoot_d8.75mM.NaCl"  "Root_d2.75mM.NaCl"   "Root_d4.75mM.NaCl"  
## [31] "Root_d6.75mM.NaCl"   "Root_d8.75mM.NaCl"   "Total_d2.75mM.NaCl" 
## [34] "Total_d4.75mM.NaCl"  "Total_d6.75mM.NaCl"  "Total_d8.75mM.NaCl" 
## [37] "RpS_d2.75mM.NaCl"    "RpS_d4.75mM.NaCl"    "RpS_d6.75mM.NaCl"   
## [40] "RpS_d8.75mM.NaCl"    "RpT_d2.75mM.NaCl"    "RpT_d4.75mM.NaCl"   
## [43] "RpT_d6.75mM.NaCl"    "RpT_d8.75mM.NaCl"    "Shoot.GR.75mM.NaCl" 
## [46] "Root.GR.75mM.NaCl"   "RpS.GR.75mM.NaCl"    "Shoot_d2.125mM.NaCl"
## [49] "Shoot_d4.125mM.NaCl" "Shoot_d6.125mM.NaCl" "Shoot_d8.125mM.NaCl"
## [52] "Root_d2.125mM.NaCl"  "Root_d4.125mM.NaCl"  "Root_d6.125mM.NaCl" 
## [55] "Root_d8.125mM.NaCl"  "Total_d2.125mM.NaCl" "Total_d4.125mM.NaCl"
## [58] "Total_d6.125mM.NaCl" "Total_d8.125mM.NaCl" "RpS_d2.125mM.NaCl"  
## [61] "RpS_d4.125mM.NaCl"   "RpS_d6.125mM.NaCl"   "RpS_d8.125mM.NaCl"  
## [64] "RpT_d2.125mM.NaCl"   "RpT_d4.125mM.NaCl"   "RpT_d6.125mM.NaCl"  
## [67] "RpT_d8.125mM.NaCl"   "Shoot.GR.125mM.NaCl" "Root.GR.125mM.NaCl" 
## [70] "RpS.GR.125mM.NaCl"
sum_data_all$Total.STI.d2.75mMNaCl <- sum_data_all$Total_d2.75mM.NaCl / sum_data_all$Total_d2.0mM.NaCl
sum_data_all$Total.STI.d4.75mMNaCl <- sum_data_all$Total_d4.75mM.NaCl / sum_data_all$Total_d4.0mM.NaCl
sum_data_all$Total.STI.d6.75mMNaCl <- sum_data_all$Total_d6.75mM.NaCl / sum_data_all$Total_d6.0mM.NaCl
sum_data_all$Total.STI.d8.75mMNaCl <- sum_data_all$Total_d8.75mM.NaCl / sum_data_all$Total_d8.0mM.NaCl

sum_data_all$Total.STI.d2.125mMNaCl <- sum_data_all$Total_d2.125mM.NaCl / sum_data_all$Total_d2.0mM.NaCl
sum_data_all$Total.STI.d4.125mMNaCl <- sum_data_all$Total_d4.125mM.NaCl / sum_data_all$Total_d4.0mM.NaCl
sum_data_all$Total.STI.d6.125mMNaCl <- sum_data_all$Total_d6.125mM.NaCl / sum_data_all$Total_d6.0mM.NaCl
sum_data_all$Total.STI.d8.125mMNaCl <- sum_data_all$Total_d8.125mM.NaCl / sum_data_all$Total_d8.0mM.NaCl

sum_data_all$Root.GR.STI.75mM.NaCl <- sum_data_all$Root.GR.75mM.NaCl / sum_data_all$Root.GR.0mM.NaCl
sum_data_all$Root.GR.STI.125mM.NaCl <- sum_data_all$Root.GR.125mM.NaCl / sum_data_all$Root.GR.0mM.NaCl

sum_data_all$Shoot.GR.STI.75mM.NaCl <- sum_data_all$Shoot.GR.75mM.NaCl / sum_data_all$Shoot.GR.0mM.NaCl
sum_data_all$Shoot.GR.STI.125mM.NaCl <- sum_data_all$Shoot.GR.125mM.NaCl / sum_data_all$Shoot.GR.0mM.NaCl

sum_data_all
library("corrplot")
## corrplot 0.92 loaded
library(RColorBrewer)
library("lattice")
 
SDA <- na.omit(sum_data_all)
dim(SDA)
## [1] 276  82
colnames(SDA)
##  [1] "CS"                      "Shoot_d2.0mM.NaCl"      
##  [3] "Shoot_d4.0mM.NaCl"       "Shoot_d6.0mM.NaCl"      
##  [5] "Shoot_d8.0mM.NaCl"       "Root_d2.0mM.NaCl"       
##  [7] "Root_d4.0mM.NaCl"        "Root_d6.0mM.NaCl"       
##  [9] "Root_d8.0mM.NaCl"        "Total_d2.0mM.NaCl"      
## [11] "Total_d4.0mM.NaCl"       "Total_d6.0mM.NaCl"      
## [13] "Total_d8.0mM.NaCl"       "RpS_d2.0mM.NaCl"        
## [15] "RpS_d4.0mM.NaCl"         "RpS_d6.0mM.NaCl"        
## [17] "RpS_d8.0mM.NaCl"         "RpT_d2.0mM.NaCl"        
## [19] "RpT_d4.0mM.NaCl"         "RpT_d6.0mM.NaCl"        
## [21] "RpT_d8.0mM.NaCl"         "Shoot.GR.0mM.NaCl"      
## [23] "Root.GR.0mM.NaCl"        "RpS.GR.0mM.NaCl"        
## [25] "Shoot_d2.75mM.NaCl"      "Shoot_d4.75mM.NaCl"     
## [27] "Shoot_d6.75mM.NaCl"      "Shoot_d8.75mM.NaCl"     
## [29] "Root_d2.75mM.NaCl"       "Root_d4.75mM.NaCl"      
## [31] "Root_d6.75mM.NaCl"       "Root_d8.75mM.NaCl"      
## [33] "Total_d2.75mM.NaCl"      "Total_d4.75mM.NaCl"     
## [35] "Total_d6.75mM.NaCl"      "Total_d8.75mM.NaCl"     
## [37] "RpS_d2.75mM.NaCl"        "RpS_d4.75mM.NaCl"       
## [39] "RpS_d6.75mM.NaCl"        "RpS_d8.75mM.NaCl"       
## [41] "RpT_d2.75mM.NaCl"        "RpT_d4.75mM.NaCl"       
## [43] "RpT_d6.75mM.NaCl"        "RpT_d8.75mM.NaCl"       
## [45] "Shoot.GR.75mM.NaCl"      "Root.GR.75mM.NaCl"      
## [47] "RpS.GR.75mM.NaCl"        "Shoot_d2.125mM.NaCl"    
## [49] "Shoot_d4.125mM.NaCl"     "Shoot_d6.125mM.NaCl"    
## [51] "Shoot_d8.125mM.NaCl"     "Root_d2.125mM.NaCl"     
## [53] "Root_d4.125mM.NaCl"      "Root_d6.125mM.NaCl"     
## [55] "Root_d8.125mM.NaCl"      "Total_d2.125mM.NaCl"    
## [57] "Total_d4.125mM.NaCl"     "Total_d6.125mM.NaCl"    
## [59] "Total_d8.125mM.NaCl"     "RpS_d2.125mM.NaCl"      
## [61] "RpS_d4.125mM.NaCl"       "RpS_d6.125mM.NaCl"      
## [63] "RpS_d8.125mM.NaCl"       "RpT_d2.125mM.NaCl"      
## [65] "RpT_d4.125mM.NaCl"       "RpT_d6.125mM.NaCl"      
## [67] "RpT_d8.125mM.NaCl"       "Shoot.GR.125mM.NaCl"    
## [69] "Root.GR.125mM.NaCl"      "RpS.GR.125mM.NaCl"      
## [71] "Total.STI.d2.75mMNaCl"   "Total.STI.d4.75mMNaCl"  
## [73] "Total.STI.d6.75mMNaCl"   "Total.STI.d8.75mMNaCl"  
## [75] "Total.STI.d2.125mMNaCl"  "Total.STI.d4.125mMNaCl" 
## [77] "Total.STI.d6.125mMNaCl"  "Total.STI.d8.125mMNaCl" 
## [79] "Root.GR.STI.75mM.NaCl"   "Root.GR.STI.125mM.NaCl" 
## [81] "Shoot.GR.STI.75mM.NaCl"  "Shoot.GR.STI.125mM.NaCl"
SDA_m = SDA[,2:ncol(SDA)]
row.names(SDA_m) = SDA$CS
head(SDA_m)
SDA_cor = cor(SDA_m, method=c("pearson"))

rgb.palette <- colorRampPalette(c("blue", "white", "orange"), space = "rgb")

levelplot(SDA_cor, xlab="", ylab="", col.regions=rgb.palette(120), cuts=100, at=seq(-1,1,0.1))

pval_SDA <- cor.mtest(SDA_m, conf.level=0.95)


corrplot(SDA_cor, p.mat=pval_SDA$p, sig.level = c(0.001, 0.01, 0.05), method = "circle", type = "lower", insig = "blank", siag = FALSE)
## Warning in text.default(pos.xlabel[, 1], pos.xlabel[, 2], newcolnames, srt =
## tl.srt, : "siag" is not a graphical parameter
## Warning in text.default(pos.ylabel[, 1], pos.ylabel[, 2], newrownames, col =
## tl.col, : "siag" is not a graphical parameter
## Warning in title(title, ...): "siag" is not a graphical parameter

on