\(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':
##
## expandlibrary(lmerTest)##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## steplibrary(MCMC.qpcr)## Loading required package: MCMCglmm## Loading required package: coda## Loading required package: apelibrary(Rmisc)## Loading required package: latticelibrary(nlme)##
## Attaching package: 'nlme'## The following object is masked from 'package:lme4':
##
## lmListlibrary(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, unionx2<-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 genesStatistics: Fitting linear mixed effects model
Full model analyzing all 3 Hsps together
- Response variable: log2(RelativeExpression)
- Fixed effects:
- Jdaycont - Julian day continuous
- Site - Harvard or Duke Forest
- baittemp.ave - average bait temperature in which ants were collected
- RIN_Value - RNA quality covariate
- 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.3aicmod1## 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.1144Results: 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.90fullmod7sel## 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 197Results: 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.64fullmod6sel## 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 197Results: 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.21fullmod8sel## 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 197Results: 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.0004433findat.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$doyhsp83
#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 197GSL$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.883915ggplot(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.4841771ggplot(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()
##tak ave
ggplot(GSL,aes(x=Jdaycont,y=hsp40pred,colour=Site))+geom_point()
##histrograms
ggplot(GSL,aes(x=hsp40pred,colour=Site,fill=Site))+geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(GSL,aes(x=hsp40pred,colour=Site,fill=Site))+geom_density(alpha=0.5)
###who is more stressed?
ddply(GSL,.(Site),summarize,netFC=sum(hsp40pred),mean=mean(hsp40pred))## Site netFC mean
## 1 DF -1650.8412 -0.28886111
## 2 HF 390.2097 0.04339521ggplot(GSL,aes(x=Site,y=hsp40pred,colour=Site))+geom_boxplot()
All figures
GSL.long<-gather(GSL,gene,PredictedFC,hsp83pred:hsp40pred)
##histograms/desnity
ggplot(GSL.long,aes(x=PredictedFC,colour=Site,fill=Site))+geom_density(alpha=0.5)+facet_grid(gene~.)
#scatter plot with predicted lines and the region where we constructed the model
ggplot(GSL.long,aes(x=baittemp.ave,y=PredictedFC,colour=Site))+geom_point(colour="gray")+stat_smooth(method="lm")+facet_grid(.~gene)+geom_vline(xintercept=range(findat.long$baittemp.ave))
###ID'ing where expression values for HF = 0
GSL.long.ave<-ddply(GSL.long,.(Site,gene))SessionInfo()
sessionInfo()## R version 3.3.3 (2017-03-06)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: macOS Sierra 10.12.3
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] dplyr_0.5.0 MASS_7.3-47 nlme_3.1-131 Rmisc_1.5
## [5] lattice_0.20-35 MCMC.qpcr_1.2.3 MCMCglmm_2.24 ape_4.1
## [9] coda_0.19-1 lmerTest_2.0-33 lme4_1.1-13 Matrix_1.2-10
## [13] tidyr_0.6.2 ggplot2_2.2.1 plyr_1.8.4
##
## loaded via a namespace (and not attached):
## [1] prettydoc_0.2.0 reshape2_1.4.2 splines_3.3.3
## [4] colorspace_1.3-2 htmltools_0.3.6 yaml_2.1.14
## [7] base64enc_0.1-3 survival_2.41-3 nloptr_1.0.4
## [10] foreign_0.8-68 DBI_0.6-1 RColorBrewer_1.1-2
## [13] stringr_1.2.0 munsell_0.4.3 gtable_0.2.0
## [16] htmlwidgets_0.8 evaluate_0.10 labeling_0.3
## [19] latticeExtra_0.6-28 knitr_1.15.1 parallel_3.3.3
## [22] htmlTable_1.9 Rcpp_0.12.10 acepack_1.4.1
## [25] corpcor_1.6.9 scales_0.4.1 backports_1.0.5
## [28] checkmate_1.8.2 Hmisc_4.0-3 gridExtra_2.2.1
## [31] tensorA_0.36 digest_0.6.12 stringi_1.1.5
## [34] grid_3.3.3 rprojroot_1.2 tools_3.3.3
## [37] magrittr_1.5 lazyeval_0.2.0 tibble_1.3.0
## [40] Formula_1.2-1 cluster_2.0.6 data.table_1.10.4
## [43] assertthat_0.2.0 minqa_1.2.4 rmarkdown_1.5
## [46] cubature_1.3-6 R6_2.2.0 rpart_4.1-11
## [49] nnet_7.3-12