\(if(toc)\)

\(toctitle\)

\(toc\)

\(endif\)

Loading Library

library(plyr)
library(ggplot2)
library(tidyr)
library(lme4)
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following object is masked from 'package:tidyr':
## 
##     expand
library(lmerTest)
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
library(MCMC.qpcr)
## Loading required package: MCMCglmm
## Loading required package: coda
## Loading required package: ape
library(Rmisc)
## Loading required package: lattice
library(nlme)
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:lme4':
## 
##     lmList
library(MASS)

Delta delta ct calculations

reading in data and parsing

x<-read.csv("../Data/20170427_final_dataset_2013_2014_years_HF_DF_del.csv",skip=10)
str(x)
## 'data.frame':    232 obs. of  20 variables:
##  $ Collection.Date: int  20140911 20140911 20140911 20130910 20140911 20130626 20130910 20130910 20130702 20130626 ...
##  $ Year_collect   : int  2014 2014 2014 2013 2014 2013 2013 2013 2013 2013 ...
##  $ Site           : Factor w/ 2 levels "DF","HF": 1 1 1 1 1 2 1 1 1 2 ...
##  $ Vial.me        : Factor w/ 218 levels "DF 1.1","DF 1.2",..: 2 30 13 44 40 139 17 29 67 142 ...
##  $ Cham           : int  1 5 13 1 8 11 14 4 15 2 ...
##  $ Sample         : int  2 2 2 2 2 3 2 2 3 3 ...
##  $ Window         : Factor w/ 9 levels "","1","2","3",..: 7 8 7 7 5 1 8 7 7 1 ...
##  $ BaitTemp1      : num  28.6 28.2 28.4 28.6 32.8 26.2 29.4 27.4 23.6 26.2 ...
##  $ BaitTemp2      : num  28.8 28 28.6 30.6 32.4 25.6 28.4 26.8 23.6 26.2 ...
##  $ BaitTemp3      : num  28.4 28.2 28.8 29.4 33 25.6 30.2 26.4 24 26.4 ...
##  $ BaitTemp4      : num  NA NA NA 30.4 NA 26.8 30.8 27.6 23.8 26.4 ...
##  $ R.conc.        : Factor w/ 190 levels "","10","10.4",..: 36 36 37 38 39 40 41 76 107 108 ...
##  $ Isolation.Date : int  20150624 20150710 20150824 20150729 20150817 20150805 20150723 20150723 20150811 20150810 ...
##  $ CT_18s         : num  17.6 16.3 15.5 18.2 16.9 ...
##  $ CT_40          : num  31.9 31.6 31.3 31.6 31.1 ...
##  $ CT_70          : num  27.5 26.8 26.9 25.1 27.5 ...
##  $ CT_83          : num  34.3 30 30.8 27.3 28.1 ...
##  $ CT_actin       : num  29 26.9 32 27.5 32.3 ...
##  $ CT_gapdh       : num  27.9 25.6 30.3 29.1 30.4 ...
##  $ RIN_Value      : num  7 6 5.5 1 2.9 1 1 1.1 2 1 ...
x$baittemp.ave<-apply(x[,8:11],1,mean,na.rm=TRUE)
x$Cham<-as.factor(as.character(x$Cham))
###Chamber data
chd<-read.csv("../Data/2017_chamber_dat.csv")
chd$Cham<-as.factor(chd$Cham)
chd$Cham2<-as.factor(chd$Cham2)
str(chd)
## 'data.frame':    30 obs. of  4 variables:
##  $ Cham : Factor w/ 15 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ Site : Factor w/ 2 levels "DF","HF": 2 2 2 2 2 2 2 2 2 2 ...
##  $ Delta: num  5.5 4 2 0 5 0 3 1.5 2.5 4.5 ...
##  $ Cham2: Factor w/ 30 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:MASS':
## 
##     select
## The following object is masked from 'package:nlme':
## 
##     collapse
## The following objects are masked from 'package:plyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
x2<-inner_join(x,chd,by=c("Cham","Site"))
## Warning in inner_join_impl(x, y, by$x, by$y, suffix$x, suffix$y): joining
## factors with different levels, coercing to character vector
#Calculating julian day
v2<-paste(substr(x2$Collection.Date,1,4),"-",substr(x2$Collection.Date,5,6),"-",substr(x2$Collection.Date,7,8))
x2$Date<-as.character(gsub(" ", "", v2, fixed = TRUE))
x2$JulianDay<-as.numeric(format(as.Date(x2$Date),"%j"))
x2$Jdaycont<-ifelse(x2$Year_collect==2013,x2$JulianDay,x2$JulianDay+365)
x2<-x2[!is.na(x2$RIN_Value),]
x3<-na.omit(x2[,-8:-11])

Actual delta delta CT calculations:

#hkg<-cbind(x3$CT_18s,x3$CT_actin,x3$CT_gapdh)
hkg<-cbind(x3$CT_actin,x3$CT_gapdh)
x3$geomean<-apply(hkg,1,function(x){exp(mean(log(x)))})

globalmean<-apply(x3[,10:15],2,mean)
deltaCTnum<-data.frame(t(apply(x3[,10:15],1,function(x){globalmean-x})))
deltaCTden<-mean(x3$geomean)-x3$geomean

##calculating delta delta ct
#gxp<-data.frame(log2(apply(deltaCTnum,2,function(x){2^(x)/2^(deltaCTden)})))
gxp1<-data.frame(log2(apply(deltaCTnum[,1:2],2,function(x){2^(x)/1.93^(deltaCTden)}))) # 18 eff = 2
gxp2<-data.frame(log2(apply(deltaCTnum[,2:4],2,function(x){1.97^(x)/1.93^(deltaCTden)}))) # hsp effic = 1.97
gxp3<-data.frame(log2(apply(deltaCTnum[,5:6],2,function(x){1.93^(x)/1.93^(deltaCTden)}))) # gapdh and actin eff = 1.93
gxp<-data.frame(CT_18s=gxp1[,1],gxp2,gxp3)
names(gxp)<-c("fc18s","fchsp40","fchsp70","fchsp83","fcactin","fcgapdh")

findat<-data.frame(x3,gxp)
str(findat)
## 'data.frame':    206 obs. of  29 variables:
##  $ Collection.Date: int  20140911 20140911 20140911 20130910 20140911 20130626 20130910 20130910 20130702 20130626 ...
##  $ Year_collect   : int  2014 2014 2014 2013 2014 2013 2013 2013 2013 2013 ...
##  $ Site           : Factor w/ 2 levels "DF","HF": 1 1 1 1 1 2 1 1 1 2 ...
##  $ Vial.me        : Factor w/ 218 levels "DF 1.1","DF 1.2",..: 2 30 13 44 40 139 17 29 67 142 ...
##  $ Cham           : chr  "1" "5" "13" "1" ...
##  $ Sample         : int  2 2 2 2 2 3 2 2 3 3 ...
##  $ Window         : Factor w/ 9 levels "","1","2","3",..: 7 8 7 7 5 1 8 7 7 1 ...
##  $ R.conc.        : Factor w/ 190 levels "","10","10.4",..: 36 36 37 38 39 40 41 76 107 108 ...
##  $ Isolation.Date : int  20150624 20150710 20150824 20150729 20150817 20150805 20150723 20150723 20150811 20150810 ...
##  $ CT_18s         : num  17.6 16.3 15.5 18.2 16.9 ...
##  $ CT_40          : num  31.9 31.6 31.3 31.6 31.1 ...
##  $ CT_70          : num  27.5 26.8 26.9 25.1 27.5 ...
##  $ CT_83          : num  34.3 30 30.8 27.3 28.1 ...
##  $ CT_actin       : num  29 26.9 32 27.5 32.3 ...
##  $ CT_gapdh       : num  27.9 25.6 30.3 29.1 30.4 ...
##  $ RIN_Value      : num  7 6 5.5 1 2.9 1 1 1.1 2 1 ...
##  $ baittemp.ave   : num  28.6 28.1 28.6 29.8 32.7 ...
##  $ Delta          : num  3.5 0 0 3.5 5 0 0 2 0 4 ...
##  $ Cham2          : Factor w/ 30 levels "1","2","3","4",..: 16 20 28 16 23 11 29 19 30 2 ...
##  $ Date           : chr  "2014-09-11" "2014-09-11" "2014-09-11" "2013-09-10" ...
##  $ JulianDay      : num  254 254 254 253 254 177 253 253 183 177 ...
##  $ Jdaycont       : num  619 619 619 253 619 177 253 253 183 177 ...
##  $ geomean        : num  28.5 26.2 31.1 28.3 31.3 ...
##  $ fc18s          : num  -3.356 -4.223 1.1571 -4.2239 0.0426 ...
##  $ fchsp40        : num  0.0887 -1.7144 3.1307 0.1767 3.5468 ...
##  $ fchsp70        : num  -0.196 -1.574 2.928 1.998 2.529 ...
##  $ fchsp83        : num  -5.109 -3.087 0.797 1.552 3.68 ...
##  $ fcactin        : num  -0.297 -0.365 -0.578 1.011 -0.673 ...
##  $ fcgapdh        : num  0.333 0.397 0.602 -0.988 0.693 ...
##convert to long format 
#findat.long<-gather(findat,gene,FC,fc18s:fcgapdh)
findat.long<-gather(findat,gene,FC,fchsp40:fchsp83) #taking out control genes

Statistics: Fitting linear mixed effects model

Full model analyzing all 3 Hsps together

  • Response variable: log2(RelativeExpression)
  • Fixed effects:
  1. Jdaycont - Julian day continuous
  2. Site - Harvard or Duke Forest
  3. baittemp.ave - average bait temperature in which ants were collected
  4. RIN_Value - RNA quality covariate
  5. gene - hsp83, hsp70, hsp40
  • Random effects: Chamber and vial nested within chamber
### fitting with lme
fullmod5<-lme(FC~RIN_Value+Jdaycont+gene*Site*baittemp.ave+gene*Site*Delta,random=~1|Cham2/Vial.me,data=findat.long)

# MOdel selection 
fullmod5<-lme(FC~RIN_Value+Jdaycont+gene*Site*baittemp.ave+gene*Site*Delta,random=~1|Cham2/Vial.me,data=findat.long,method="ML")
#anova(fullmod5)
#summary(stepAIC(fullmod5,direction="both"))
aicmod1<-anova(summary(stepAIC(fullmod5,direction="both")))
## Start:  AIC=2529.95
## FC ~ RIN_Value + Jdaycont + gene * Site * baittemp.ave + gene * 
##     Site * Delta
## 
##                          Df    AIC
## - gene:Site:Delta         2 2526.6
## - RIN_Value               1 2528.5
## - gene:Site:baittemp.ave  2 2528.7
## <none>                      2529.9
## - Jdaycont                1 2532.3
## 
## Step:  AIC=2526.57
## FC ~ RIN_Value + Jdaycont + gene + Site + baittemp.ave + Delta + 
##     gene:Site + gene:baittemp.ave + Site:baittemp.ave + gene:Delta + 
##     Site:Delta + gene:Site:baittemp.ave
## 
##                          Df    AIC
## - gene:Delta              2 2522.6
## - gene:Site:baittemp.ave  2 2524.7
## - RIN_Value               1 2525.1
## <none>                      2526.6
## - Site:Delta              1 2527.0
## - Jdaycont                1 2529.0
## + gene:Site:Delta         2 2529.9
## 
## Step:  AIC=2522.59
## FC ~ RIN_Value + Jdaycont + gene + Site + baittemp.ave + Delta + 
##     gene:Site + gene:baittemp.ave + Site:baittemp.ave + Site:Delta + 
##     gene:Site:baittemp.ave
## 
##                          Df    AIC
## - gene:Site:baittemp.ave  2 2520.8
## - RIN_Value               1 2521.1
## <none>                      2522.6
## - Site:Delta              1 2523.0
## - Jdaycont                1 2525.0
## + gene:Delta              2 2526.6
## 
## Step:  AIC=2520.77
## FC ~ RIN_Value + Jdaycont + gene + Site + baittemp.ave + Delta + 
##     gene:Site + gene:baittemp.ave + Site:baittemp.ave + Site:Delta
## 
##                          Df    AIC
## - RIN_Value               1 2519.3
## <none>                      2520.8
## - Site:Delta              1 2521.2
## - Site:baittemp.ave       1 2521.4
## - gene:Site               2 2522.0
## + gene:Site:baittemp.ave  2 2522.6
## - Jdaycont                1 2523.2
## + gene:Delta              2 2524.7
## - gene:baittemp.ave       2 2537.8
## 
## Step:  AIC=2519.29
## FC ~ Jdaycont + gene + Site + baittemp.ave + Delta + gene:Site + 
##     gene:baittemp.ave + Site:baittemp.ave + Site:Delta
## 
##                          Df    AIC
## <none>                      2519.3
## - Site:Delta              1 2520.0
## - gene:Site               2 2520.5
## + RIN_Value               1 2520.8
## - Site:baittemp.ave       1 2521.0
## + gene:Site:baittemp.ave  2 2521.1
## - Jdaycont                1 2522.1
## + gene:Delta              2 2523.2
## - gene:baittemp.ave       2 2536.3
aicmod1
##                   numDF denDF   F-value p-value
## (Intercept)           1   413  0.044025  0.8339
## Jdaycont              1   413  4.633790  0.0319
## gene                  2   413  0.000000  1.0000
## Site                  1    23 21.076865  0.0001
## baittemp.ave          1   413 29.836209  <.0001
## Delta                 1    23  0.169841  0.6841
## gene:Site             2   413  1.682588  0.1872
## gene:baittemp.ave     2   413 10.550810  <.0001
## Site:baittemp.ave     1   169  1.819052  0.1792
## Site:Delta            1    23  2.692626  0.1144

Results: gene by bait temp effect, site effect

Separate linear mixed effects models for each Hsp (hsp70,hsp83,hsp40)

hsp83

########hsp83
hsp83<-subset(findat.long,findat.long$gene=="fchsp83")
fullmod7<-lme(FC~RIN_Value+Jdaycont+Site*baittemp.ave+Site*Delta,random=~1|Cham2/Vial.me,data=hsp83,method="ML")
#fullmod7<-lme(FC~RIN_Value+Jdaycont+Site*baittemp.ave+Site*Delta,random=~1|Cham2/Vial.me,data=hsp83)
#summary(fullmod7)
fullmod7sel<-summary(stepAIC(fullmod7,direction="both"))
## Start:  AIC=856.75
## FC ~ RIN_Value + Jdaycont + Site * baittemp.ave + Site * Delta
## 
##                     Df    AIC
## <none>                 856.75
## - Site:Delta         1 857.90
## - RIN_Value          1 858.14
## - Site:baittemp.ave  1 859.20
## - Jdaycont           1 861.90
fullmod7sel
## Linear mixed-effects model fit by maximum likelihood
##  Data: hsp83 
##        AIC      BIC    logLik
##   856.7462 893.3528 -417.3731
## 
## Random effects:
##  Formula: ~1 | Cham2
##         (Intercept)
## StdDev:   0.1946163
## 
##  Formula: ~1 | Vial.me %in% Cham2
##         (Intercept) Residual
## StdDev:    1.531984  1.03929
## 
## Fixed effects: FC ~ RIN_Value + Jdaycont + Site * baittemp.ave + Site * Delta 
##                         Value Std.Error  DF   t-value p-value
## (Intercept)         -9.915097  1.528096 169 -6.488532  0.0000
## RIN_Value            0.149878  0.082116   6  1.825203  0.1178
## Jdaycont             0.002731  0.001028   6  2.655774  0.0377
## SiteHF              -7.726510  4.135596  23 -1.868294  0.0745
## baittemp.ave         0.295245  0.058582   6  5.039819  0.0024
## Delta                0.008764  0.091198  23  0.096100  0.9243
## SiteHF:baittemp.ave  0.329418  0.157427 169  2.092515  0.0379
## SiteHF:Delta        -0.292487  0.166271  23 -1.759101  0.0919
##  Correlation: 
##                     (Intr) RIN_Vl Jdycnt SiteHF bttmp. Delta  StHF:.
## RIN_Value           -0.147                                          
## Jdaycont             0.089 -0.691                                   
## SiteHF              -0.404  0.291 -0.225                            
## baittemp.ave        -0.964  0.057 -0.164  0.377                     
## Delta                0.017 -0.070  0.040 -0.022 -0.129              
## SiteHF:baittemp.ave  0.401 -0.305  0.246 -0.993 -0.386  0.068       
## SiteHF:Delta        -0.019  0.111 -0.099  0.343  0.081 -0.553 -0.427
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.90917202 -0.24335576  0.07299595  0.29028657  1.92310971 
## 
## Number of Observations: 206
## Number of Groups: 
##              Cham2 Vial.me %in% Cham2 
##                 27                197

Results: Significant Site by Bait temp interaction

hsp70

###testing interaction for each hsp gene
hsp70<-subset(findat.long,findat.long$gene=="fchsp70")
fullmod6<-lme(FC~RIN_Value+Jdaycont+Site*baittemp.ave+Site*Delta,random=~1|Cham2/Vial.me,data=hsp70,method="ML")
fullmod6sel<-summary(stepAIC(fullmod6,direction="both"))
## Start:  AIC=821.91
## FC ~ RIN_Value + Jdaycont + Site * baittemp.ave + Site * Delta
## 
##                     Df    AIC
## - RIN_Value          1 820.04
## - Site:baittemp.ave  1 820.13
## - Site:Delta         1 820.62
## <none>                 821.91
## - Jdaycont           1 845.43
## 
## Step:  AIC=820.04
## FC ~ Jdaycont + Site + baittemp.ave + Delta + Site:baittemp.ave + 
##     Site:Delta
## 
##                     Df    AIC
## - Site:baittemp.ave  1 818.39
## - Site:Delta         1 818.83
## <none>                 820.04
## + RIN_Value          1 821.91
## - Jdaycont           1 861.91
## 
## Step:  AIC=818.39
## FC ~ Jdaycont + Site + baittemp.ave + Delta + Site:Delta
## 
##                     Df    AIC
## - Site:Delta         1 816.90
## <none>                 818.39
## + Site:baittemp.ave  1 820.04
## + RIN_Value          1 820.13
## - baittemp.ave       1 845.70
## - Jdaycont           1 860.64
## 
## Step:  AIC=816.9
## FC ~ Jdaycont + Site + baittemp.ave + Delta
## 
##                     Df    AIC
## - Delta              1 815.24
## <none>                 816.90
## + Site:Delta         1 818.39
## + RIN_Value          1 818.64
## + Site:baittemp.ave  1 818.83
## - Site               1 822.86
## - baittemp.ave       1 843.71
## - Jdaycont           1 859.22
## 
## Step:  AIC=815.24
## FC ~ Jdaycont + Site + baittemp.ave
## 
##                     Df    AIC
## <none>                 815.24
## + Delta              1 816.90
## + RIN_Value          1 817.02
## + Site:baittemp.ave  1 817.22
## - Site               1 820.90
## - baittemp.ave       1 842.26
## - Jdaycont           1 857.64
fullmod6sel
## Linear mixed-effects model fit by maximum likelihood
##  Data: hsp70 
##        AIC      BIC    logLik
##   815.2412 838.5363 -400.6206
## 
## Random effects:
##  Formula: ~1 | Cham2
##         (Intercept)
## StdDev:   0.2225052
## 
##  Formula: ~1 | Vial.me %in% Cham2
##          (Intercept) Residual
## StdDev: 0.0002263996 1.677996
## 
## Fixed effects: FC ~ Jdaycont + Site + baittemp.ave 
##                  Value Std.Error  DF   t-value p-value
## (Intercept)  -6.364371 1.3026365 170 -4.885761   0e+00
## Jdaycont     -0.004589 0.0006572   7 -6.981991   2e-04
## SiteHF        0.788088 0.2736238  25  2.880187   8e-03
## baittemp.ave  0.289517 0.0499854   7  5.792034   7e-04
##  Correlation: 
##              (Intr) Jdycnt SiteHF
## Jdaycont     -0.025              
## SiteHF        0.178  0.173       
## baittemp.ave -0.974 -0.168 -0.281
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.1274266 -0.6098562 -0.1545382  0.3930239  3.7842314 
## 
## Number of Observations: 206
## Number of Groups: 
##              Cham2 Vial.me %in% Cham2 
##                 27                197

Results: Main effect of Bait temp and Site

hsp40

#### hsp40

hsp40<-subset(findat.long,findat.long$gene=="fchsp40")
fullmod8<-lme(FC~RIN_Value+Jdaycont+Site*baittemp.ave+Site*Delta,random=~1|Cham2/Vial.me,data=hsp40,method="ML")

fullmod8sel<-summary(stepAIC(fullmod8,direction="both"))
## Start:  AIC=872.57
## FC ~ RIN_Value + Jdaycont + Site * baittemp.ave + Site * Delta
## 
##                     Df    AIC
## - Site:Delta         1 871.61
## - RIN_Value          1 871.62
## - Site:baittemp.ave  1 871.70
## <none>                 872.57
## - Jdaycont           1 876.06
## 
## Step:  AIC=871.61
## FC ~ RIN_Value + Jdaycont + Site + baittemp.ave + Delta + Site:baittemp.ave
## 
##                     Df    AIC
## - Delta              1 869.92
## - Site:baittemp.ave  1 870.07
## - RIN_Value          1 870.44
## <none>                 871.61
## + Site:Delta         1 872.57
## - Jdaycont           1 875.67
## 
## Step:  AIC=869.92
## FC ~ RIN_Value + Jdaycont + Site + baittemp.ave + Site:baittemp.ave
## 
##                     Df    AIC
## - Site:baittemp.ave  1 868.26
## - RIN_Value          1 868.76
## <none>                 869.92
## + Delta              1 871.61
## - Jdaycont           1 873.98
## 
## Step:  AIC=868.26
## FC ~ RIN_Value + Jdaycont + Site + baittemp.ave
## 
##                     Df    AIC
## - RIN_Value          1 866.87
## <none>                 868.26
## + Site:baittemp.ave  1 869.92
## + Delta              1 870.07
## - Site               1 870.54
## - baittemp.ave       1 872.37
## - Jdaycont           1 873.33
## 
## Step:  AIC=866.87
## FC ~ Jdaycont + Site + baittemp.ave
## 
##                     Df    AIC
## <none>                 866.87
## + RIN_Value          1 868.26
## - Site               1 868.63
## + Delta              1 868.63
## + Site:baittemp.ave  1 868.76
## - baittemp.ave       1 870.44
## - Jdaycont           1 883.21
fullmod8sel
## Linear mixed-effects model fit by maximum likelihood
##  Data: hsp40 
##        AIC     BIC    logLik
##   866.8738 890.169 -426.4369
## 
## Random effects:
##  Formula: ~1 | Cham2
##          (Intercept)
## StdDev: 0.0007261298
## 
##  Formula: ~1 | Vial.me %in% Cham2
##          (Intercept) Residual
## StdDev: 0.0006572257 1.917732
## 
## Fixed effects: FC ~ Jdaycont + Site + baittemp.ave 
##                   Value Std.Error  DF   t-value p-value
## (Intercept)  -2.8557557 1.4687690 170 -1.944319  0.0535
## Jdaycont     -0.0032511 0.0007497   7 -4.336718  0.0034
## SiteHF        0.5809013 0.2956002  25  1.965159  0.0606
## baittemp.ave  0.1430383 0.0564215   7  2.535174  0.0389
##  Correlation: 
##              (Intr) Jdycnt SiteHF
## Jdaycont     -0.025              
## SiteHF        0.195  0.183       
## baittemp.ave -0.974 -0.170 -0.295
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.1019055 -0.6332219 -0.2321619  0.4740407  3.0532535 
## 
## Number of Observations: 206
## Number of Groups: 
##              Cham2 Vial.me %in% Cham2 
##                 27                197

Results: Main effect of Bait temp

Figures

FC vs bait temp

#linear fit
ggplot(findat.long,aes(x=baittemp.ave,y=FC,colour=Site))+geom_point()+stat_smooth(method="lm")+facet_grid(.~gene)

#ggplot(findat.long,aes(x=baittemp.ave,y=FC,colour=gene))+facet_grid(.~Site)+geom_point(colour="gray")+stat_smooth(method="lm",se=FALSE)
## trying to get residuals 
summary(lm(FC~RIN_Value+Jdaycont,data=findat.long))
## 
## Call:
## lm(formula = FC ~ RIN_Value + Jdaycont, data = findat.long)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.4612 -1.4788 -0.2408  1.1393  7.8652 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.0644959  0.2192484   0.294  0.76873    
## RIN_Value    0.1548627  0.0495557   3.125  0.00186 ** 
## Jdaycont    -0.0024262  0.0006234  -3.892  0.00011 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.125 on 615 degrees of freedom
## Multiple R-squared:  0.0248, Adjusted R-squared:  0.02163 
## F-statistic: 7.819 on 2 and 615 DF,  p-value: 0.0004433
findat.long$residuals<-residuals(lm(FC~RIN_Value+Jdaycont,data=findat.long))

ggplot(findat.long,aes(x=baittemp.ave,y=residuals,colour=Site))+geom_point()+stat_smooth(method="lm")+facet_grid(.~gene)

FC vs sites for each gene

findat.long$gene<-as.factor(findat.long$gene)
gg<-summarySE(findat.long,measurevar ="FC",groupvars=c("Site","gene"))
ggplot(gg,aes(x=Site,y=FC,colour=gene,group=gene))+geom_errorbar(aes(ymin=FC-se,ymax=FC+se),width=.4,position=position_dodge(.5))+geom_point(size=3,position=position_dodge(.5))+geom_line(position=position_dodge(.5))

Expression predictions in nature

Predicting expression levels given the temperatures ants experience in nature

growing season length data

GSL<-read.csv("gslave.csv")
str(GSL)
## 'data.frame':    14707 obs. of  7 variables:
##  $ X      : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Site   : Factor w/ 2 levels "DF","HF": 1 1 1 1 1 1 1 1 1 1 ...
##  $ year   : int  2011 2011 2011 2011 2011 2011 2011 2011 2012 2012 ...
##  $ doy    : int  1 1 1 2 2 2 2 2 7 7 ...
##  $ Cham   : int  12 8 9 12 3 4 8 9 3 9 ...
##  $ Delta  : num  4 5 5.5 4 4.5 2 5 5.5 4.5 5.5 ...
##  $ tempave: num  15.1 16.8 16.7 16.4 16.8 ...
names(GSL)[7]<-"baittemp.ave"
GSL$Jdaycont<-GSL$doy

hsp83

#GSL$RIN_Value<-sample(seq(0,10,1),size=length(GSL$X),replace=TRUE)

# the prediction model 
#predmod83<-lme(FC~Site*baittemp.ave,random=~1|Cham2/Vial.me,data=hsp83)

predmod83<-lme(FC~Site*baittemp.ave+Delta*Site+Jdaycont,random=~1|Cham2/Vial.me,data=hsp83)
summary(predmod83)
## Linear mixed-effects model fit by REML
##  Data: hsp83 
##        AIC      BIC    logLik
##   884.6457 917.5788 -432.3229
## 
## Random effects:
##  Formula: ~1 | Cham2
##         (Intercept)
## StdDev:   0.3567442
## 
##  Formula: ~1 | Vial.me %in% Cham2
##         (Intercept) Residual
## StdDev:    1.508096  1.11765
## 
## Fixed effects: FC ~ Site * baittemp.ave + Delta * Site + Jdaycont 
##                          Value Std.Error  DF   t-value p-value
## (Intercept)          -9.531546  1.539041 169 -6.193171  0.0000
## SiteHF              -10.091788  3.968694  23 -2.542849  0.0182
## baittemp.ave          0.289869  0.059404   7  4.879653  0.0018
## Delta                 0.020474  0.098398  23  0.208072  0.8370
## Jdaycont              0.004046  0.000741   7  5.458390  0.0009
## SiteHF:baittemp.ave   0.423721  0.150266 169  2.819811  0.0054
## SiteHF:Delta         -0.331660  0.175088  23 -1.894244  0.0708
##  Correlation: 
##                     (Intr) SiteHF bttmp. Delta  Jdycnt StHF:.
## SiteHF              -0.387                                   
## baittemp.ave        -0.966  0.382                            
## Delta               -0.011  0.005 -0.116                     
## Jdaycont            -0.016 -0.037 -0.170 -0.010              
## SiteHF:baittemp.ave  0.382 -0.992 -0.393  0.046  0.054       
## SiteHF:Delta         0.007  0.296  0.071 -0.562 -0.030 -0.391
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.94288531 -0.29149384  0.08241722  0.30198883  2.00057663 
## 
## Number of Observations: 206
## Number of Groups: 
##              Cham2 Vial.me %in% Cham2 
##                 27                197
GSL$hsp83pred<-as.vector(predict(predmod83,newdata=GSL,level=0))

#plot of pred vs bait
ggplot(GSL,aes(x=baittemp.ave,y=hsp83pred,colour=Site))+geom_point()

ggplot(GSL,aes(x=hsp83pred,colour=Site,fill=Site))+geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(GSL,aes(x=hsp83pred,colour=Site,fill=Site))+geom_density(alpha=0.5)

#color by delta
ggplot(GSL,aes(x=hsp83pred,colour=Site,fill=Site))+geom_histogram()+facet_grid(Delta~.)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(GSL,aes(x=hsp83pred,colour=Site,fill=Site))+geom_density(alpha=0.5)+facet_grid(Delta~.)

##ok let's plot delta on X axis
ggplot(GSL,aes(x=factor(Delta),y=hsp83pred,colour=Site))+geom_boxplot()

##potentially, who is more stressed? 
## summing across all values? 
ddply(GSL,.(Site),summarize,netFC=sum(hsp83pred),mean=mean(hsp83pred))
##   Site     netFC      mean
## 1   DF -12920.81 -2.260859
## 2   HF -43916.16 -4.883915
ggplot(GSL,aes(x=Site,y=hsp83pred,colour=Site))+geom_boxplot()

hsp70

## prediction model
fullmod6sel
## Linear mixed-effects model fit by maximum likelihood
##  Data: hsp70 
##        AIC      BIC    logLik
##   815.2412 838.5363 -400.6206
## 
## Random effects:
##  Formula: ~1 | Cham2
##         (Intercept)
## StdDev:   0.2225052
## 
##  Formula: ~1 | Vial.me %in% Cham2
##          (Intercept) Residual
## StdDev: 0.0002263996 1.677996
## 
## Fixed effects: FC ~ Jdaycont + Site + baittemp.ave 
##                  Value Std.Error  DF   t-value p-value
## (Intercept)  -6.364371 1.3026365 170 -4.885761   0e+00
## Jdaycont     -0.004589 0.0006572   7 -6.981991   2e-04
## SiteHF        0.788088 0.2736238  25  2.880187   8e-03
## baittemp.ave  0.289517 0.0499854   7  5.792034   7e-04
##  Correlation: 
##              (Intr) Jdycnt SiteHF
## Jdaycont     -0.025              
## SiteHF        0.178  0.173       
## baittemp.ave -0.974 -0.168 -0.281
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.1274266 -0.6098562 -0.1545382  0.3930239  3.7842314 
## 
## Number of Observations: 206
## Number of Groups: 
##              Cham2 Vial.me %in% Cham2 
##                 27                197
#pretend doy is julian day
GSL$hsp70pred<-as.vector(predict(fullmod6sel,GSL,level=0))

#
ggplot(GSL,aes(x=baittemp.ave,y=hsp70pred,colour=Site))+geom_point()

ggplot(GSL,aes(x=hsp70pred,colour=Site,fill=Site))+geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(GSL,aes(x=hsp70pred,colour=Site,fill=Site))+geom_density(alpha=0.5)

###who is more stressed
ddply(GSL,.(Site),summarize,netFC=sum(hsp70pred),mean=mean(hsp70pred))
##   Site     netFC       mean
## 1   DF -4527.034 -0.7921319
## 2   HF -4353.721 -0.4841771
ggplot(GSL,aes(x=Site,y=hsp70pred,colour=Site))+geom_boxplot()

hsp40

#prediction model 
fullmod8sel
## Linear mixed-effects model fit by maximum likelihood
##  Data: hsp40 
##        AIC     BIC    logLik
##   866.8738 890.169 -426.4369
## 
## Random effects:
##  Formula: ~1 | Cham2
##          (Intercept)
## StdDev: 0.0007261298
## 
##  Formula: ~1 | Vial.me %in% Cham2
##          (Intercept) Residual
## StdDev: 0.0006572257 1.917732
## 
## Fixed effects: FC ~ Jdaycont + Site + baittemp.ave 
##                   Value Std.Error  DF   t-value p-value
## (Intercept)  -2.8557557 1.4687690 170 -1.944319  0.0535
## Jdaycont     -0.0032511 0.0007497   7 -4.336718  0.0034
## SiteHF        0.5809013 0.2956002  25  1.965159  0.0606
## baittemp.ave  0.1430383 0.0564215   7  2.535174  0.0389
##  Correlation: 
##              (Intr) Jdycnt SiteHF
## Jdaycont     -0.025              
## SiteHF        0.195  0.183       
## baittemp.ave -0.974 -0.170 -0.295
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.1019055 -0.6332219 -0.2321619  0.4740407  3.0532535 
## 
## Number of Observations: 206
## Number of Groups: 
##              Cham2 Vial.me %in% Cham2 
##                 27                197
#predic based on GSL
GSL$hsp40pred<-as.vector(predict(fullmod8sel,newdata=GSL,level=0))

##predict bait temp from 0 FC
#HF 
HFzero_hsp40<-subset(GSL,GSL$Site=="HF" & GSL$hsp40pred<0.01 & GSL$hsp40pred>-.01)
mean(HFzero_hsp40$baittemp.ave)
## [1] 20.28599
###DF 
DFzero_hsp40<-subset(GSL,GSL$Site=="DF" & GSL$hsp40pred<0.01 & GSL$hsp40pred>-.01)
mean(DFzero_hsp40$baittemp.ave)
## [1] 23.62782
#mean(DFzero_hsp40$hsp40pred)
### 23.62782-20.28599
#[1] 3.34183
#3.3 C diff!

##figures
ggplot(GSL,aes(x=baittemp.ave,y=hsp40pred,colour=Site))+geom_point()