Lovell — Jun 10, 2014, 12:31 PM
library(ggplot2)
library(GGally)
library(reshape)
rm(list=ls())
############
#Part I: Processing, vetting and organization of data
#1.1: Read in the data
setwd("~/Desktop/Adaptomics_Project/Drydown_Analysis")
rawD13<-read.csv("mr2014_d13Craw_edit1.csv", header=T, na.strings="#VALUE!")
spatD13<-rawD13[1:192,-7]
dd_info<-read.csv("mr2014_drydowninfo_jtledit.csv", header=T, na.strings="N/A") #weights at harvest
dd_wts<-read.csv("mr2014_allrydownweights.csv", header=T, na.strings="N/A") #daily weights
#names(dd_wts)[4:6]<-c("pot_wt","drysoil_wt","init_wt")
# Look at the data, make general plots of drydown conditions
summary(dd_wts)
Flat Position Pot.ID Experiment Treatment
Min. :1.00 Min. : 1.00 Min. : 9.0 Min. :1.0 Dry:128
1st Qu.:1.75 1st Qu.: 8.75 1st Qu.: 62.8 1st Qu.:1.0
Median :2.50 Median :16.50 Median :126.5 Median :1.5
Mean :2.50 Mean :16.50 Mean :125.8 Mean :1.5
3rd Qu.:3.25 3rd Qu.:24.25 3rd Qu.:176.8 3rd Qu.:2.0
Max. :4.00 Max. :32.00 Max. :300.0 Max. :2.0
NA's :8
Genotype Pot.cloth.wt..g. Dry.soil.wt..g.
Alv#2 :13 Min. : 8.94 Min. :66.0
Chicago#4 :12 1st Qu.: 9.80 1st Qu.:69.9
Alv#3 :11 Median : 9.89 Median :71.2
Chicago#2 : 8 Mean : 9.76 Mean :71.5
Chiquito#4: 8 3rd Qu.: 9.98 3rd Qu.:73.1
(Other) :68 Max. :10.14 Max. :77.4
NA's : 8 NA's :8 NA's :8
New.Weight..g..of.Transplant.Pots Day1 Day2
Min. :119 Min. :111 Min. :105
1st Qu.:134 1st Qu.:130 1st Qu.:122
Median :139 Median :134 Median :126
Mean :138 Mean :134 Mean :125
3rd Qu.:142 3rd Qu.:138 3rd Qu.:129
Max. :155 Max. :149 Max. :139
NA's :8 NA's :8 NA's :8
Day3 Day4 Day5 Day6
Min. :101 Min. : 95.1 Min. : 88.2 Min. : 79.8
1st Qu.:119 1st Qu.:112.4 1st Qu.:103.7 1st Qu.: 98.8
Median :123 Median :116.4 Median :109.7 Median :104.0
Mean :122 Mean :116.1 Mean :108.9 Mean :103.8
3rd Qu.:126 3rd Qu.:119.7 3rd Qu.:114.2 3rd Qu.:108.9
Max. :139 Max. :132.4 Max. :126.3 Max. :121.4
NA's :8 NA's :8 NA's :8 NA's :8
Day7 Day8 Day9 Day10
Min. : 76.0 Min. : 73.7 Min. : 84.2 Min. : 80.5
1st Qu.: 96.4 1st Qu.: 93.2 1st Qu.: 90.5 1st Qu.: 88.1
Median :100.8 Median : 97.3 Median : 93.2 Median : 90.6
Mean :100.5 Mean : 96.8 Mean : 93.3 Mean : 90.5
3rd Qu.:104.3 3rd Qu.:100.4 3rd Qu.: 95.9 3rd Qu.: 93.3
Max. :114.6 Max. :110.8 Max. :105.3 Max. :102.9
NA's :8 NA's :8 NA's :8 NA's :8
colnames(dd_wts)[7:9]<-c("pot_wt","dry_wt","wet_wt")
ddwts_long <- melt(dd_wts, id.vars = c("Flat","Position","Pot.ID","Experiment","Treatment","Genotype","pot_wt","dry_wt","wet_wt"))
colnames(ddwts_long)[10:11]<-c("day","weight")
#calculate grav soil moisture
#gsm=(wt-dry)/fullwet
ddwts_long<-ddwts_long[complete.cases(ddwts_long),]
ddwts_long$grav_moisture<-(ddwts_long$weight-ddwts_long$dry_wt)/(ddwts_long$wet_wt-ddwts_long$dry_wt)
ggplot(ddwts_long,aes(y=grav_moisture,x=day))+
geom_line(aes(group=Pot.ID,col=as.factor(Experiment)))+
geom_abline(intercept=0.4, slope=0)+
theme_bw()
#the drydown was a bit more gradual and in the end, less dry in the 2nd experiment
#genotypes are nested within experiment- no experiment * genotype interaction possible
#1.2: Combine relevant information
#spatifolia is in the 1st 192 rows
head(spatD13)
Tray_Name Well_Id Sample_ID d13C C_Amount_.ug. sample_wt
1 spatifolia_1 A1 147 -32.46 585.7 1760
2 spatifolia_1 A2 18 -32.99 762.4 2083
3 spatifolia_1 A3 148 -33.04 736.2 1932
4 spatifolia_1 A4 81 -32.29 501.6 1719
5 spatifolia_1 A5 115 -32.65 690.5 2020
6 spatifolia_1 A6 51 -33.03 703.5 1945
info<-dd_info[,c(5,2:4,10,1,6:7,9)]
names(info)[6:9]<-c("Sample_ID","pot_wt_atharvest","fresh_rosette_mass","dry_root_wt")
head(info)
Genotype Position Treatment Flat Experiment Sample_ID pot_wt_atharvest
1 Alvarado#2 8 Dry 1 2 167 89.59
2 Alvarado#2 12 Dry 1 2 165 95.31
3 Alvarado#2 21 Dry 1 2 168 88.85
4 Alvarado#2 22 Dry 1 2 224 93.21
5 Alvarado#2 26 Dry 1 2 193 93.29
6 Alvarado#2 28 Dry 1 2 80 89.52
fresh_rosette_mass dry_root_wt
1 29.2 6.4
2 57.7 10.1
3 94.3 19.3
4 99.3 14.5
5 36.9 7.5
6 45.0 8.1
alldata<-merge(info,spatD13,by="Sample_ID")
dd_lastday_GWC<-ddwts_long[ddwts_long$day=="Day10",c(1:6,12)]
names(dd_lastday_GWC)[3]<-"Sample_ID"
#1.3 Data quality check
#just for the drydown conditions:
dry_treatment_alldata<-merge(alldata,dd_lastday_GWC,by="Sample_ID", all.x=F)
dry_treatment_alldata<-dry_treatment_alldata[complete.cases(dry_treatment_alldata),]
pairs(dry_treatment_alldata[c(7:9,12:13,20)],pch = 21,
bg = c("red", "blue")[unclass(dry_treatment_alldata$Experiment.x)],
upper.panel=panel.smooth, main= "trait correlations in 1st (RED) and 2nd (BLUE) experiments")
#consistent with the second experiment being a bit wetter, there is higher
#gravimetric soil water content on average in the second experiment
tapply(dry_treatment_alldata$grav_moisture,dry_treatment_alldata$Experiment.x,mean)
1 2
0.2672 0.2996
#higher root mass
tapply(dry_treatment_alldata$dry_root_wt,dry_treatment_alldata$Experiment.x,mean)
1 2
3.283 8.062
#
#1.4 Process info, bring apomixis location etc. info into the fold
genotype_info<-alldata$Genotype
ids<-as.data.frame(do.call(rbind, strsplit(as.character(genotype_info),"#")))
colnames(ids)<-c("population","genotype")
apodata<-data.frame(c("Alvarado", "Alvarado","Chicago","Chicago","Chiquito" ,"Chiquito" ,"Cripple","Cripple" ,"Rosita","Rosita","Royal","Royal","Tiesiding","Tiesiding"),
c(2,3,2,4,4,7,6,7,3,4,1,2,2,7),c("apo","sex",'sex',"apo","apo","sex","sex","apo","sex","apo","apo","sex","sex","apo"))
names(apodata)<-c("population","genotype","mating_system")
alldata<-cbind(ids,alldata)
all<-merge(alldata,apodata,by=c("population","genotype"),all.x=T,all.y=T)
#1.4 Add in row/column data
all$position<-as.numeric(as.character(all$Position))
pos<-c(1:32)
row<-rep(1:8,4)
col<-rep(1:4,each=8)
all$row<-all$position
all$col<-all$position
for (i in 1:32){
all$row[all$row==pos[i]]<-row[i]
all$col[all$col==pos[i]]<-col[i]
}
head(all)
population genotype Sample_ID Genotype Position Treatment Flat
1 Alvarado 2 166 Alvarado#2 26 Wet 3
2 Alvarado 2 200 Alvarado#2 14 Wet 3
3 Alvarado 2 201 Alvarado#2 26 Wet 4
4 Alvarado 2 73 Alvarado#2 13 Dry 2
5 Alvarado 2 193 Alvarado#2 26 Dry 1
6 Alvarado 2 163 Alvarado#2 8 Wet 4
Experiment pot_wt_atharvest fresh_rosette_mass dry_root_wt Tray_Name
1 2 122.51 87.7 10.5 spatifolia_2
2 2 135.73 181.2 14.5 spatifolia_2
3 2 130.09 80.5 8.3 spatifolia_2
4 2 87.28 46.6 8.5 spatifolia_2
5 2 93.29 36.9 7.5 spatifolia_2
6 2 124.68 153.2 23.6 spatifolia_2
Well_Id d13C C_Amount_.ug. sample_wt mating_system position row col
1 E12 -32.34 792.9 1948 apo 26 2 4
2 E6 -33.18 880.3 2177 apo 14 6 2
3 G12 -31.21 749.3 1812 apo 26 2 4
4 C5 -31.65 763.8 1903 apo 13 5 2
5 B2 -30.23 763.8 1845 apo 26 2 4
6 F11 -32.66 742.1 1836 apo 8 8 1
#add in experiment information
head(info)
Genotype Position Treatment Flat Experiment Sample_ID pot_wt_atharvest
1 Alvarado#2 8 Dry 1 2 167 89.59
2 Alvarado#2 12 Dry 1 2 165 95.31
3 Alvarado#2 21 Dry 1 2 168 88.85
4 Alvarado#2 22 Dry 1 2 224 93.21
5 Alvarado#2 26 Dry 1 2 193 93.29
6 Alvarado#2 28 Dry 1 2 80 89.52
fresh_rosette_mass dry_root_wt
1 29.2 6.4
2 57.7 10.1
3 94.3 19.3
4 99.3 14.5
5 36.9 7.5
6 45.0 8.1
exp<-info[,c(6,5)]
all<-merge(all,exp)
#1.5 Analyze the conditions by row/column in the flats
wet<-all[all$Treatment=="Wet",]
dry<-all[all$Treatment=="Dry",]
ggplot(wet,aes(x=col,y=row)) +
geom_tile(aes(fill = pot_wt_atharvest),colour = "white") +
scale_fill_gradient(low = "white", high= "red")+
facet_wrap(~Flat+Experiment)+
ggtitle("wet treatment pot weight")+
theme_bw()
ggplot(dry,aes(x=col,y=row)) +
geom_tile(aes(fill = pot_wt_atharvest),colour = "white") +
scale_fill_gradient(low = "white", high= "red")+
facet_wrap(~Flat+Experiment)+
ggtitle("dry treatment pot weight")+
theme_bw()
########
# Part II: Basic analyses
#2.1: correlations between soil moisture and d13C
# across all treatments
lm_all<-lm(d13C~pot_wt_atharvest,data=all)
summary(lm_all)
Call:
lm(formula = d13C ~ pot_wt_atharvest, data = all)
Residuals:
Min 1Q Median 3Q Max
-1.6518 -0.5645 -0.0913 0.6149 2.2812
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -28.99474 0.33575 -86.4 < 2e-16 ***
pot_wt_atharvest -0.02662 0.00303 -8.8 1.1e-15 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.84 on 182 degrees of freedom
Multiple R-squared: 0.298, Adjusted R-squared: 0.294
F-statistic: 77.4 on 1 and 182 DF, p-value: 1.07e-15
ggplot(all, aes(y=d13C,x=pot_wt_atharvest))+
geom_point(aes(col=Treatment))+
stat_smooth(method=lm, aes(group=Treatment))+
ggtitle("correlations between soil GWC and d13C split by experiments")+
facet_wrap(~Experiment)+
theme_bw()
#strong effect of experiment on the results.
#esspecially the response in d13C
with(all,
by(all,Treatment,
function (x) summary(lm(d13C~pot_wt_atharvest,data=x))))
Treatment: Dry
Call:
lm(formula = d13C ~ pot_wt_atharvest, data = x)
Residuals:
Min 1Q Median 3Q Max
-1.8215 -0.3842 -0.0362 0.3253 1.6503
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -31.85387 1.64239 -19.39 <2e-16 ***
pot_wt_atharvest 0.00591 0.01852 0.32 0.75
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.657 on 88 degrees of freedom
Multiple R-squared: 0.00116, Adjusted R-squared: -0.0102
F-statistic: 0.102 on 1 and 88 DF, p-value: 0.75
--------------------------------------------------------
Treatment: Wet
Call:
lm(formula = d13C ~ pot_wt_atharvest, data = x)
Residuals:
Min 1Q Median 3Q Max
-1.627 -0.838 -0.257 0.720 2.344
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -31.45656 2.66454 -11.81 <2e-16 ***
pot_wt_atharvest -0.00765 0.02070 -0.37 0.71
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.979 on 92 degrees of freedom
Multiple R-squared: 0.00148, Adjusted R-squared: -0.00937
F-statistic: 0.136 on 1 and 92 DF, p-value: 0.713
#for the wet treatment
with(wet,
by(wet,Experiment,
function (x) summary(lm(d13C~pot_wt_atharvest,data=x))))
Experiment: 1
Call:
lm(formula = d13C ~ pot_wt_atharvest, data = x)
Residuals:
Min 1Q Median 3Q Max
-0.9818 -0.2683 0.0348 0.3122 1.1834
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -31.99777 1.71590 -18.65 <2e-16 ***
pot_wt_atharvest -0.00844 0.01336 -0.63 0.53
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.487 on 47 degrees of freedom
Multiple R-squared: 0.00841, Adjusted R-squared: -0.0127
F-statistic: 0.399 on 1 and 47 DF, p-value: 0.531
--------------------------------------------------------
Experiment: 2
Call:
lm(formula = d13C ~ pot_wt_atharvest, data = x)
Residuals:
Min 1Q Median 3Q Max
-1.9515 -0.5939 0.0651 0.6539 1.5897
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -27.9120 3.8523 -7.25 5.7e-09 ***
pot_wt_atharvest -0.0297 0.0299 -1.00 0.33
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.893 on 43 degrees of freedom
Multiple R-squared: 0.0225, Adjusted R-squared: -0.00022
F-statistic: 0.99 on 1 and 43 DF, p-value: 0.325
#for the dry treatment
with(dry,
by(dry,Experiment,
function (x) summary(lm(d13C~pot_wt_atharvest,data=x))))
Experiment: 1
Call:
lm(formula = d13C ~ pot_wt_atharvest, data = x)
Residuals:
Min 1Q Median 3Q Max
-1.5746 -0.2590 0.0551 0.2790 1.3638
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -30.5373 2.1434 -14.25 <2e-16 ***
pot_wt_atharvest -0.0107 0.0246 -0.44 0.67
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.516 on 40 degrees of freedom
Multiple R-squared: 0.00472, Adjusted R-squared: -0.0202
F-statistic: 0.19 on 1 and 40 DF, p-value: 0.665
--------------------------------------------------------
Experiment: 2
Call:
lm(formula = d13C ~ pot_wt_atharvest, data = x)
Residuals:
Min 1Q Median 3Q Max
-1.737 -0.547 0.000 0.541 1.540
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -30.31809 2.73554 -11.08 1.4e-14 ***
pot_wt_atharvest -0.00988 0.03035 -0.33 0.75
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.749 on 46 degrees of freedom
Multiple R-squared: 0.0023, Adjusted R-squared: -0.0194
F-statistic: 0.106 on 1 and 46 DF, p-value: 0.746
#even experiment #2 wet does not have a significant correlation between soil
#moisture and d13C
#but not significant within treatments
ggplot(wet, aes(y=d13C,x=pot_wt_atharvest))+
geom_point(aes(col=as.factor(Experiment)))+
stat_smooth(method=lm, aes(col=as.factor(Experiment),group=as.factor(Experiment)))+
ggtitle("WET treatment corr. between soil GWC and d13C split by experiments")+
theme_bw()
ggplot(dry, aes(y=d13C,x=pot_wt_atharvest))+
geom_point(aes(col=as.factor(Experiment)))+
stat_smooth(method=lm, aes(col=as.factor(Experiment), group=as.factor(Experiment)))+
ggtitle("DRY treatment corr. between soil GWC and d13C split by experiments")+
theme_bw()
#2.2 Basic fixed effect statistics
aov1<-aov(d13C~population*Treatment*mating_system, data=all)
summary(aov1)
Df Sum Sq Mean Sq F value Pr(>F)
population 6 43.7 7.3 22.17 < 2e-16 ***
Treatment 1 50.1 50.1 152.53 < 2e-16 ***
mating_system 1 0.1 0.1 0.28 0.60
population:Treatment 6 14.7 2.4 7.43 5.2e-07 ***
population:mating_system 6 21.5 3.6 10.92 3.9e-10 ***
Treatment:mating_system 1 0.0 0.0 0.09 0.77
population:Treatment:mating_system 6 1.6 0.3 0.79 0.58
Residuals 156 51.3 0.3
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(all,aes(x=Treatment,y=d13C,col=mating_system, group=mating_system))+
geom_point()+
theme_bw()+
facet_wrap(~population)
ggplot(all,aes(x=Treatment,y=d13C))+
geom_boxplot(aes(fill=mating_system))+
theme_bw()+
facet_wrap(~population)
#mixed effect ANOVA controlling for variation among experiments
#######