library(ggplot2) library(grid) library(gridExtra) library(ggpubr) library(sciplot) library(DHARMa) library(lme4) library(stats) ####################################################################### ####################################################################### ####################################################################### ####################### Graphs ############################### ####################################################################### ####################################################################### ####################################################################### ### all graphs require the data sheet in the excel file with the same figure label. ## R script for Fig. 2A needs data file Fig. 2A ## R script for Fig. 2B needs data file Fig. 2B ## R script for Fig. 3A needs data file Fig. 3A ## R script for Fig. 3B needs data file Fig. 3B ## R script for Fig. S1 needs data file Fig. S1 ## and so on ################################ Fig. 2A######################### ######################### Intracolony transmission. A. The proportion of cages in which experimentally infected donor bees (1st column of a treatment) and initially uninfected recipient bees (2nd column of a treatment) were infected with DWV by day 7. virus.transmission.part1_proportions<-read.table(file.choose(),header=T) cpPalette=c( "violetred3", "palegreen3", "steelblue") ggplot(virus.transmission.part1_proportions, aes(x=treatment, y=proportion,fill = actor)) + geom_bar(stat="identity",color="black",alpha=0.8,width=0.7, position=position_dodge())+ scale_fill_manual(values=cpPalette)+ scale_color_manual(values = cpPalette)+ scale_x_discrete(name="direction of transmission", labels=c("Apis to Apis","Apis to Bombus", "Bombus to Apis","Bombus to Bombus"))+ scale_y_continuous(name="proprtion of cages with infected individuals")+ theme(axis.line = element_line(colour = "black",size=.8), axis.ticks.x=element_line(colour="transparent",size=0.001), axis.ticks.y=element_line(colour="black",size=0.2), axis.ticks.length=unit(0.3,"cm"), axis.text.x=element_text(colour="black",size=10), axis.text.y=element_text(colour="black",size=15), axis.title.x=element_text(colour="black",size=17,margin=margin(15,0,0,0)), axis.title.y=element_text(colour="black",size=17,margin=margin(0,10,0,0)), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank(), legend.title=element_text(size=12), legend.text=element_text(size=12), legend.key.size=unit(0.8,"cm"), legend.key=element_blank(), legend.position="right") ################################ Fig. 2B######################### ################ Intracolony transmission. B. Viral (DWV) titres of donor and recipient bees at day 7. virus.transmission.part1<-read.table(file.choose(),header=T) ggplot(virus.transmission.part1, aes(x = treatment, y = virus35+1)) + geom_boxplot(outlier.shape=NA,aes(fill = status),lwd=0.5,alpha=0.6,width=0.8)+ geom_jitter(aes(bg = status),position=position_jitterdodge(dodge.width=0.8, jitter.width=0.3),pch = 21,col = "black", lwd = 3.0, alpha =6/10)+ scale_fill_manual(values=ctPalette)+ scale_color_manual(values = ctPalette)+ scale_x_discrete(name="transmission route",labels=c("Apis to Apis","Apis to Bombus", "Bombus to Apis","Bombus to Bombus"))+ scale_y_continuous( trans="log10",name="virus load per bee",breaks=c(1, 10^3, 10^6,10^9, 10^12, 10^15),limits=c(1,2000000000000000),labels=trans_format("log10",math_format(10^.x)))+ geom_segment(aes(x = 2.5, y = 0, xend = 2.5, yend = 200000000000000), col= "azure2", lwd = 0.5)+ geom_segment(aes(x = 3.5, y = 0, xend = 3.5, yend = 200000000000000), col= "azure2", lwd = 0.5)+ geom_segment(aes(x = 1.5, y = 0, xend = 1.5, yend = 2000000000000000), col= "azure2", lwd = 0.5)+ annotation_logticks(sides="l",size=0.2, short = unit(0.05, "cm"), mid = unit(0.1, "cm"), long = unit(0.2, "cm"))+ theme(axis.line = element_line(colour = "black",size=.8), axis.ticks.x=element_line(colour="transparent",size=0.001), axis.ticks.y=element_line(colour="black",size=0.2), axis.ticks.length=unit(0.3,"cm"), axis.text.x=element_text(colour="black",size=13), axis.text.y=element_text(colour="black",size=15), axis.title.x=element_text(colour="black",size=20,margin=margin(15,0,0,0)), axis.title.y=element_text(colour="black",size=20,margin=margin(0,10,0,0)), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank(), legend.title=element_text(size=12), legend.text=element_text(size=12), legend.key.size=unit(0.8,"cm"), legend.key=element_blank(), legend.position="right") ############################ Fig. 3A ############### ######################### A. The proportion of cages in which experimentally infected donor bees (1st column of a treatment) and initially uninfected recipient bees (2nd and 3rd columns of a treatment) were infected with DWV by day 7 (2nd column: green) and day 14 (3rd column: blue). virus.transmission.part2_proportions<-read.table(file.choose(),header=T) cpPalette=c( "violetred3", "palegreen3", "steelblue") ggplot(virus.transmission.part2_proportions, aes(x=treatment, y=proportion,fill = actor)) + geom_bar(stat="identity",color="black",alpha=0.8,width=0.7, position=position_dodge())+ scale_fill_manual(values=cpPalette)+ scale_color_manual(values = cpPalette)+ scale_x_discrete(name="direction of transmission", labels=c("Apis to Apis","Apis to Bombus", "Bombus to Apis","Bombus to Bombus"))+ scale_y_continuous(name="proprtion of cages with infected individuals")+ theme(axis.line = element_line(colour = "black",size=.8), axis.ticks.x=element_line(colour="transparent",size=0.001), axis.ticks.y=element_line(colour="black",size=0.2), axis.ticks.length=unit(0.3,"cm"), axis.text.x=element_text(colour="black",size=10), axis.text.y=element_text(colour="black",size=15), axis.title.x=element_text(colour="black",size=20,margin=margin(15,0,0,0)), axis.title.y=element_text(colour="black",size=20,margin=margin(0,10,0,0)), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank(), legend.title=element_text(size=12), legend.text=element_text(size=12), legend.key.size=unit(0.8,"cm"), legend.key=element_blank(), legend.position="right") ############################ Fig. 3B ############### ############################ B. Viral (DWV) titres of experimentally infected donor bees at day 7 and of recipient bees at days 7 (green) and 14 (blue). virus.transmission.part2<-read.table(file.choose(),header=T) ggplot(virus.transmission.part2, aes(x = treatment, y = virus35+1)) + geom_boxplot(outlier.shape=NA,aes(fill = statusday),lwd=0.5,alpha=0.6,width=0.8)+ geom_jitter(aes(bg = statusday),position=position_jitterdodge(dodge.width=0.8, jitter.width=0.3),pch = 21,col = "black", lwd = 3.0, alpha =6/10)+ scale_fill_manual(values=ctPalette)+ scale_color_manual(values = ctPalette)+ scale_x_discrete(name="transmission route",labels=c("Apis to Apis","Apis to Bombus","Bombus to Apis","Bombus to Bombus" ))+ scale_y_continuous( trans="log10",name="virus load per bee",breaks=c(1, 10^3, 10^6,10^9, 10^12, 10^15),limits=c(1,2000000000000000),labels=trans_format("log10",math_format(10^.x)))+ geom_segment(aes(x = 2.5, y = 0, xend = 2.5, yend = 200000000000000), col= "azure2", lwd = 0.5)+ geom_segment(aes(x = 1.5, y = 0, xend = 1.5, yend = 2000000000000000), col= "azure2", lwd = 0.5)+ geom_segment(aes(x = 3.5, y = 0, xend = 3.5, yend = 2000000000000000), col= "azure2", lwd = 0.5)+ annotation_logticks(sides="l",size=0.2, short = unit(0.05, "cm"), mid = unit(0.1, "cm"), long = unit(0.2, "cm"))+ theme(axis.line = element_line(colour = "black",size=.8), axis.ticks.x=element_line(colour="transparent",size=0.001), axis.ticks.y=element_line(colour="black",size=0.2), axis.ticks.length=unit(0.3,"cm"), axis.text.x=element_text(colour="black",size=13), axis.text.y=element_text(colour="black",size=15), axis.title.x=element_text(colour="black",size=20,margin=margin(15,0,0,0)), axis.title.y=element_text(colour="black",size=20,margin=margin(0,10,0,0)), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank(), legend.title=element_text(size=12), legend.text=element_text(size=12), legend.key.size=unit(0.8,"cm"), legend.key=element_blank(), legend.position="right") ############################ Fig. S1 ############### ###########Viral titre in recipient Bombus in the treatment Apis to Bombus of experiment 1: intracolony transmission in relation to the force of infection. Viral (DWV-A) load of bumble bees from the Apis to Bombus treatment (in which honey virus.transmission.viruspressure<-read.table(file.choose(),header=T) ggplot(virus.transmission.viruspressure, aes(x = after5dayswords, y = virus)) + geom_boxplot(outlier.shape=NA, aes(fill = "blue"), lwd=0.5,alpha=0.6,width=0.8)+ geom_jitter(aes(bg = "red"),position=position_jitterdodge(dodge.width=0.8, jitter.width=0.3),pch = 21,col = "black", lwd = 3.0, alpha =6/10)+ scale_fill_manual(values = c("paleturquoise3","deeppink3" ))+ scale_x_discrete(name="viral pressure at day 5",labels=c("high","low"))+ scale_y_continuous( trans="log10",name="virus load per bee",breaks=c(1, 10^2, 10^4,10^6, 10^8, 10^10),limits=c(1,20000000000),labels=trans_format("log10",math_format(10^.x)))+ annotation_logticks(sides="l",size=0.2, short = unit(0.05, "cm"), mid = unit(0.1, "cm"), long = unit(0.2, "cm"))+ theme(axis.line = element_line(colour = "black",size=.8), axis.ticks.x=element_line(colour="transparent",size=0.001), axis.ticks.y=element_line(colour="black",size=0.2), axis.ticks.length=unit(0.3,"cm"), axis.text.x=element_text(colour="black",size=13), axis.text.y=element_text(colour="black",size=15), axis.title.x=element_text(colour="black",size=20,margin=margin(15,0,0,0)), axis.title.y=element_text(colour="black",size=20,margin=margin(0,10,0,0)), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank(), legend.title=element_text(size=12), legend.text=element_text(size=12), legend.key.size=unit(0.8,"cm"), legend.key=element_blank(), legend.position="right") ############################ Fig. S2 A ############### ##################Control bee viral titres in experiments 1 (intracolony transmission). Viral (DWV-A) titre of untreated, control bees at day 7 virus.transmission.control1<-read.table(file.choose(),header=T) ccPalette=c("hotpink3", "cyan4") ggplot(virus.transmission.control1, aes(x = status, y = virus)) + geom_boxplot(outlier.shape=NA,aes(fill = weeks),lwd=0.5,alpha=0.8,width=0.8)+ geom_jitter(aes(bg = weeks),position=position_jitterdodge(dodge.width=0.8, jitter.width=0.3),pch = 21,col = "black", lwd = 4.0, alpha =8/10)+ scale_fill_manual(values=ccPalette)+ scale_color_manual(values = ccPalette)+ scale_y_continuous( trans="log10",name="virus load per bee",breaks=c(1, 10^3, 10^6,10^9),limits=c(1,2000000000),labels=trans_format("log10",math_format(10^.x)))+ annotation_logticks(sides="l",size=0.2, short = unit(0.05, "cm"), mid = unit(0.1, "cm"), long = unit(0.2, "cm"))+ theme(axis.line = element_line(colour = "black",size=.8), axis.ticks.x=element_line(colour="transparent",size=0.001), axis.ticks.y=element_line(colour="black",size=0.2), axis.ticks.length=unit(0.3,"cm"), axis.text.x=element_text(colour="black",size=13), axis.text.y=element_text(colour="black",size=15), axis.title.x=element_text(colour="black",size=20,margin=margin(15,0,0,0)), axis.title.y=element_text(colour="black",size=20,margin=margin(0,10,0,0)), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank(), legend.title=element_text(size=12), legend.text=element_text(size=12), legend.key.size=unit(0.8,"cm"), legend.key=element_blank(), legend.position="right") ############################ Fig. S2 B ############### ####################Figure S2B Control bee viral titres in experiments 2. Viral (DWV-A) titre of untreated, control bees after 7 days and 14 days virus.transmission.control2<-read.table(file.choose(),header=T) ccPalette=c("hotpink3", "cyan4") ggplot(virus.transmission.control2, aes(x = status, y = virus)) + geom_boxplot(outlier.shape=NA,aes(fill = weeks),lwd=0.5,alpha=0.8,width=0.8)+ geom_jitter(aes(bg = weeks),position=position_jitterdodge(dodge.width=0.8, jitter.width=0.3),pch = 21,col = "black", lwd = 4.0, alpha =8/10)+ scale_fill_manual(values=ccPalette)+ scale_color_manual(values = ccPalette)+ scale_y_continuous( trans="log10",name="virus load per bee",breaks=c(1, 10^3, 10^6,10^9),limits=c(1,2000000000),labels=trans_format("log10",math_format(10^.x)))+ annotation_logticks(sides="l",size=0.2, short = unit(0.05, "cm"), mid = unit(0.1, "cm"), long = unit(0.2, "cm"))+ theme(axis.line = element_line(colour = "black",size=.8), axis.ticks.x=element_line(colour="transparent",size=0.001), axis.ticks.y=element_line(colour="black",size=0.2), axis.ticks.length=unit(0.3,"cm"), axis.text.x=element_text(colour="black",size=13), axis.text.y=element_text(colour="black",size=15), axis.title.x=element_text(colour="black",size=20,margin=margin(15,0,0,0)), axis.title.y=element_text(colour="black",size=20,margin=margin(0,10,0,0)), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank(), legend.title=element_text(size=12), legend.text=element_text(size=12), legend.key.size=unit(0.8,"cm"), legend.key=element_blank(), legend.position="right") ##################### Fig. S3A ########################## #############################Viral prevalence and viral titres per cage for Experiment 2, mimicking floral transmission. A The proportion of recipient bees per cage which were infected with DWV by day 7 (2nd column: green) and day 14 (3rd column: blue) for all four treatments (AA: Apis to Apis; AB: Apis to Bombus; BA: Bombus to Apis; and BB: Bombus to Bombus). transmissionexperiment_part2_justrecipients_all_cages_all_treatments<-read.table(file.choose(),header=T) crPalette=c("palegreen3", "steelblue") ggplot(transmissionexperiment_part2_justrecipients_all_cages_all_treatments, aes(x=cagetreatment, y=infected,fill = days))+ geom_bar(stat="identity",color="black", alpha=0.8,width=0.7, position=position_dodge())+ scale_fill_manual(values=crPalette)+ scale_color_manual(values = crPalette)+ #scale_x_discrete(name="direction of transmission", labels=c("Apis to Apis","Apis to Bombus", "Bombus to Apis","Bombus to Bombus"))+ scale_y_continuous(name="proprtion of infected individuals per cage",limits=c(0,100))+ theme(axis.line = element_line(colour = "black",size=.8), axis.ticks.x=element_line(colour="transparent",size=0.001), axis.ticks.y=element_line(colour="black",size=0.2), axis.ticks.length=unit(0.3,"cm"), axis.text.x=element_text(colour="black",size=10), axis.text.y=element_text(colour="black",size=15), axis.title.x=element_text(colour="black",size=17,margin=margin(15,0,0,0)), axis.title.y=element_text(colour="black",size=17,margin=margin(0,10,0,0)), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank(), legend.title=element_text(size=12), legend.text=element_text(size=12), legend.key.size=unit(0.8,"cm"), legend.key=element_blank(), legend.position="right") ###################### Fig. S3B #################### ############### Viral prevalence and viral titres per cage for Experiment 2, mimicking floral transmission. B. Viral (DWV) titres of recipient bees per cage at days 7 (green) and 14 (blue) for all four treatments (AA, AB, BA, BB); donor bee icons are in pink and recipient bee icons are in green crPalette=c("palegreen3", "steelblue") virus.transmission.recipients2AA<-read.table(file.choose(),header=T) ggplot(virus.transmission.recipients2AA, aes(x = cage , y = virus)) + geom_boxplot(outlier.shape=NA,aes(fill = days),lwd=0.5,alpha=0.8,width=0.8)+ geom_jitter(aes(bg = days),position=position_jitterdodge(dodge.width=0.8, jitter.width=0.3),pch = 21,col = "black", lwd = 4.0, alpha =8/10)+ scale_fill_manual(values=crPalette)+ scale_color_manual(values = crPalette)+ scale_y_continuous( trans="log10",name="virus load per bee",breaks=c(1, 10^3, 10^6,10^9, 10^12),limits=c(1,20000000000000),labels=trans_format("log10",math_format(10^.x)))+ annotation_logticks(sides="l",size=0.2, short = unit(0.05, "cm"), mid = unit(0.1, "cm"), long = unit(0.2, "cm"))+ theme(axis.line = element_line(colour = "black",size=.8), axis.ticks.x=element_line(colour="transparent",size=0.001), axis.ticks.y=element_line(colour="black",size=0.2), axis.ticks.length=unit(0.3,"cm"), axis.text.x=element_text(colour="black",size=13), axis.text.y=element_text(colour="black",size=15), axis.title.x=element_text(colour="black",size=20,margin=margin(15,0,0,0)), axis.title.y=element_text(colour="black",size=20,margin=margin(0,10,0,0)), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank(), legend.title=element_text(size=12), legend.text=element_text(size=12), legend.key.size=unit(0.8,"cm"), legend.key=element_blank(), legend.position="right") virus.transmission.recipients2AB<-read.table(file.choose(),header=T) ggplot(virus.transmission.recipients2AB, aes(x = cage , y = virus)) + geom_boxplot(outlier.shape=NA,aes(fill = days),lwd=0.5,alpha=0.8,width=0.8)+ geom_jitter(aes(bg = days),position=position_jitterdodge(dodge.width=0.8, jitter.width=0.3),pch = 21,col = "black", lwd = 4.0, alpha =8/10)+ scale_fill_manual(values=crPalette)+ scale_color_manual(values = crPalette)+ scale_y_continuous( trans="log10",name="virus load per bee",breaks=c(1, 10^3, 10^6,10^9, 10^12),limits=c(1,20000000000000),labels=trans_format("log10",math_format(10^.x)))+ annotation_logticks(sides="l",size=0.2, short = unit(0.05, "cm"), mid = unit(0.1, "cm"), long = unit(0.2, "cm"))+ theme(axis.line = element_line(colour = "black",size=.8), axis.ticks.x=element_line(colour="transparent",size=0.001), axis.ticks.y=element_line(colour="black",size=0.2), axis.ticks.length=unit(0.3,"cm"), axis.text.x=element_text(colour="black",size=13), axis.text.y=element_text(colour="black",size=15), axis.title.x=element_text(colour="black",size=20,margin=margin(15,0,0,0)), axis.title.y=element_text(colour="black",size=20,margin=margin(0,10,0,0)), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank(), legend.title=element_text(size=12), legend.text=element_text(size=12), legend.key.size=unit(0.8,"cm"), legend.key=element_blank(), legend.position="right") virus.transmission.recipients2BA<-read.table(file.choose(),header=T) ggplot(virus.transmission.recipients2BA, aes(x = cage , y = virus)) + geom_boxplot(outlier.shape=NA,aes(fill = days),lwd=0.5,alpha=0.8,width=0.8)+ geom_jitter(aes(bg = days),position=position_jitterdodge(dodge.width=0.8, jitter.width=0.3),pch = 21,col = "black", lwd = 4.0, alpha =8/10)+ scale_fill_manual(values=crPalette)+ scale_color_manual(values = crPalette)+ scale_y_continuous( trans="log10",name="virus load per bee",breaks=c(1, 10^3, 10^6,10^9, 10^12),limits=c(1,20000000000000),labels=trans_format("log10",math_format(10^.x)))+ annotation_logticks(sides="l",size=0.2, short = unit(0.05, "cm"), mid = unit(0.1, "cm"), long = unit(0.2, "cm"))+ theme(axis.line = element_line(colour = "black",size=.8), axis.ticks.x=element_line(colour="transparent",size=0.001), axis.ticks.y=element_line(colour="black",size=0.2), axis.ticks.length=unit(0.3,"cm"), axis.text.x=element_text(colour="black",size=13), axis.text.y=element_text(colour="black",size=15), axis.title.x=element_text(colour="black",size=20,margin=margin(15,0,0,0)), axis.title.y=element_text(colour="black",size=20,margin=margin(0,10,0,0)), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank(), legend.title=element_text(size=12), legend.text=element_text(size=12), legend.key.size=unit(0.8,"cm"), legend.key=element_blank(), legend.position="right") virus.transmission.recipients2BB<-read.table(file.choose(),header=T) ggplot(virus.transmission.recipients2BB, aes(x = cage , y = virus)) + geom_boxplot(outlier.shape=NA,aes(fill = days),lwd=0.5,alpha=0.8,width=0.8)+ geom_jitter(aes(bg = days),position=position_jitterdodge(dodge.width=0.8, jitter.width=0.3),pch = 21,col = "black", lwd = 4.0, alpha =8/10)+ scale_fill_manual(values=crPalette)+ scale_color_manual(values = crPalette)+ scale_y_continuous( trans="log10",name="virus load per bee",breaks=c(1, 10^3, 10^6,10^9, 10^12),limits=c(1,20000000000000),labels=trans_format("log10",math_format(10^.x)))+ annotation_logticks(sides="l",size=0.2, short = unit(0.05, "cm"), mid = unit(0.1, "cm"), long = unit(0.2, "cm"))+ theme(axis.line = element_line(colour = "black",size=.8), axis.ticks.x=element_line(colour="transparent",size=0.001), axis.ticks.y=element_line(colour="black",size=0.2), axis.ticks.length=unit(0.3,"cm"), axis.text.x=element_text(colour="black",size=13), axis.text.y=element_text(colour="black",size=15), axis.title.x=element_text(colour="black",size=20,margin=margin(15,0,0,0)), axis.title.y=element_text(colour="black",size=20,margin=margin(0,10,0,0)), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank(), legend.title=element_text(size=12), legend.text=element_text(size=12), legend.key.size=unit(0.8,"cm"), legend.key=element_blank(), legend.position="right") ####################################################################### ####################################################################### ####################################################################### ####################### Statistics ############################### ####################################################################### ####################################################################### ####################################################################### ##################################### part1: first experiment ############################################################### ##########honeybees as donors virus.transmission.part1<-read.table(file.choose(),header=T) virus.transmission.part1.apis<-subset(virus.transmission.part1,infected_species=="Apis") virus.transmission.part1.apis<-drop.levels(virus.transmission.part1.apis) str(virus.transmission.part1.apis) ###calculation of summary statistics per treatment: ApisApis_donor=AAinf, ApisApis_receiver=AAclean, ApisBombus_donor=ABinf, ApisBombus_receiver=ABclean summmary.virus.transmission.part1.apis<-ddply(virus.transmission.part1.apis,"treatment2",summarize, N=length(virus), median=quantile(virus,0.5), lowerCI=quantile(virus,0.025), upperCI=quantile(virus,0.975), mean=mean(virus), sd=sd(virus), se=sd/sqrt(N)) summmary.virus.transmission.part1.apis # treatment2 N median lowerCI upperCI mean sd se #1 AAclean 9 1.100000e+11 3.937778e+08 1.715800e+12 3.634270e+11 6.492613e+11 2.164204e+11 #2 AAinf 6 2.180000e+12 1.781250e+12 1.143625e+13 4.000000e+12 4.247004e+12 1.733832e+12 #3 ABclean 9 9.294759e+07 1.024027e+07 1.572191e+09 3.441172e+08 5.864628e+08 1.954876e+08 #4 ABinf 6 1.122000e+13 7.330000e+12 4.856250e+13 1.974333e+13 1.772920e+13 7.237916e+12 ###calculation of second summary statistics per treatment: ApisApisandApisBombus_donor=ainf, ApisApisandApisBombus_receiver=clean summmary2.virus.transmission.part1.apis<-ddply(virus.transmission.part1.apis,"status",summarize, N=length(virus), median=quantile(virus,0.5), lowerCI=quantile(virus,0.025), upperCI=quantile(virus,0.975), mean=mean(virus), sd=sd(virus), se=sd/sqrt(N)) summmary2.virus.transmission.part1.apis # status N median lowerCI upperCI mean sd se #1 ainf 12 7.360000e+12 1.818750e+12 4.539750e+13 1.187167e+13 1.478748e+13 4.268778e+12 #2 clean 18 1.210063e+09 1.252469e+07 1.384825e+12 1.818856e+11 4.829783e+11 1.138391e+11 ###statistics #model m.virus.transmission.part1.apis<-lmer(viruslog~treatment2+(1|cage2),data=virus.transmission.part1.apis) #model validation via DHARMa res.m.virus.transmission.part1.apis<-simulateResiduals(fittedModel=m.virus.transmission.part1.apis) plot(res.m.virus.transmission.part1.apis)#ok #test of treatment2 effect m1.virus.transmission.part1.apis<-lmer(viruslog~1+(1|cage2),data=virus.transmission.part1.apis) anova(m.virus.transmission.part1.apis,m1.virus.transmission.part1.apis) #refitting model(s) with ML (instead of REML) #Data: virus.transmission.part1.apis #Models: #m1.virus.transmission.part1.apis: viruslog ~ 1 + (1 | cage2) #m.virus.transmission.part1.apis: viruslog ~ treatment2 + (1 | cage2) # npar AIC BIC logLik deviance Chisq Df Pr(>Chisq) #m1.virus.transmission.part1.apis 3 138.53 142.732 -66.264 132.53 #m.virus.transmission.part1.apis 6 86.27 94.677 -37.135 74.27 58.259 3 1.384e-12 *** #multiple comparisons between levels of treatment2 m.virus.transmission.part1.apis.glht<-glht(m.virus.transmission.part1.apis,linfct=mcp(treatment2="Tukey")) summary(m.virus.transmission.part1.apis.glht,test=adjusted("Westfall")) # Simultaneous Tests for General Linear Hypotheses # #Multiple Comparisons of Means: Tukey Contrasts # # #Fit: lmer(formula = viruslog ~ treatment2 + (1 | cage2), data = virus.transmission.part1.apis) # #Linear Hypotheses: # Estimate Std. Error z value Pr(>|z|) #AAinf - AAclean == 0 1.9181 0.4249 4.514 <0.001 *** #ABclean - AAclean == 0 -2.5819 0.3911 -6.602 <0.001 *** #ABinf - AAclean == 0 2.6229 0.4382 5.986 <0.001 *** #ABclean - AAinf == 0 -4.5000 0.4463 -10.083 <0.001 *** #ABinf - AAinf == 0 0.7048 0.4898 1.439 0.15 #ABinf - ABclean == 0 5.2048 0.4198 12.399 <0.001 *** ####################Bombus virus.transmission.part1.bombus<-subset(virus.transmission.part1,infected_species=="Bombus") virus.transmission.part1.bombus<-drop.levels(virus.transmission.part1.bombus) str(virus.transmission.part1.bombus) ###calculation of summary statistics per treatment: BombusApis_donor=BAinf, BombusApis_receiver=BAclean, BombusBombus_donor=BBinf, BombusBombus_receiver=BBclean summmary.virus.transmission.part1.bombus<-ddply(virus.transmission.part1.bombus,"treatment2",summarize, N=length(virus), median=quantile(virus,0.5), lowerCI=quantile(virus,0.025), upperCI=quantile(virus,0.975), mean=mean(virus), sd=sd(virus), se=sd/sqrt(N)) summmary.virus.transmission.part1.bombus # treatment2 N median lowerCI upperCI mean sd se #1 BAclean 9 0.1 0.1 1.000000e-01 1.000000e-01 0 0 #2 BAinf 6 3115969578.5 649760073.4 3.931947e+10 9.815621e+09 16883289064 6892573898 #3 BBclean 8 0.1 0.1 1.000000e-01 1.000000e-01 0 0 #4 BBinf 7 1070902160.0 173262849.5 4.682144e+10 1.361231e+10 21699037869 8201465413 ###calculation of second summary statistics per treatment: BombusApisandBombusBombus_donor=ainf, BombusApisandBombusBombus_receiver=clean summmary2.virus.transmission.part1.bombus<-ddply(virus.transmission.part1.bombus,"status",summarize, N=length(virus), median=quantile(virus,0.5), lowerCI=quantile(virus,0.025), upperCI=quantile(virus,0.975), mean=mean(virus), sd=sd(virus), se=sd/sqrt(N)) summmary2.virus.transmission.part1.bombus # status N median lowerCI upperCI mean sd se #1 ainf 13 2.51283e+09 184031491.0 4.644965e+10 1.185999e+10 18922843028 5248252370 #2 clean 17 1.00000e-01 0.1 1.000000e-01 1.000000e-01 0 0 ###statistics #model only with BombusApis_donor=BAinf and BombusBombus_donor=BBinf virus.transmission.part1.bombus.ainf<-subset(virus.transmission.part1,status=="ainf") virus.transmission.part1.bombus.ainf<-drop.levels(virus.transmission.part1.bombus.ainf) str(virus.transmission.part1.bombus.ainf) m.virus.transmission.part1.bombus.ainf<-lmer(viruslog~treatment2+(1|cage2),data=virus.transmission.part1.bombus.ainf) #model validation via DHARMa res.m.virus.transmission.part1.bombus.ainf<-simulateResiduals(fittedModel=m.virus.transmission.part1.bombus.ainf) plot(res.m.virus.transmission.part1.bombus.ainf)#ok #test of treatment2 effect m1.virus.transmission.part1.bombus.ainf<-lmer(viruslog~1+(1|cage2),data=virus.transmission.part1.bombus.ainf) anova(m.virus.transmission.part1.bombus.ainf,m1.virus.transmission.part1.bombus.ainf) #refitting model(s) with ML (instead of REML) #Data: virus.transmission.part1.bombus.ainf #Models: #m1.virus.transmission.part1.bombus.ainf: viruslog ~ 1 + (1 | cage2) #m.virus.transmission.part1.bombus.ainf: viruslog ~ treatment2 + (1 | cage2) # npar AIC BIC logLik deviance Chisq Df Pr(>Chisq) #m1.virus.transmission.part1.bombus.ainf 3 107.013 110.670 -50.507 101.013 #m.virus.transmission.part1.bombus.ainf 6 54.685 61.998 -21.343 42.685 58.328 3 1.338e-12 *** ##################################### part2: second experiment ############################################################### virus.transmission.part2<-read.table(file.choose(),header=T) #########Apis virus.transmission.part2.apis<-subset(virus.transmission.part2,infectedspecies=="Apis") virus.transmission.part2.apis<-drop.levels(virus.transmission.part2.apis) str(virus.transmission.part2.apis) ########## virus.transmission.part2.apis.summary<-ddply(virus.transmission.part2.apis,c("treatment","status_species"),summarize, N=length(virus35), median=quantile(virus35,0.5), lowerCI=quantile(virus35,0.025), upperCI=quantile(virus35,0.975), mean=mean(virus35), sd=sd(virus35), se=sd/sqrt(N)) virus.transmission.part2.apis.summary # treatment status_species N median lowerCI upperCI mean sd se #1 AA ainf1weekApis 6 2.120000e+13 9.4025e+12 4.011250e+13 2.252667e+13 1.126095e+13 4.597265e+12 #2 AA clean1weekApis 18 0.000000e+00 0.0000e+00 6.789500e+12 8.933333e+11 2.216760e+12 5.224953e+11 #3 AA clean2weeksApis 17 5.740803e+06 0.0000e+00 1.076617e+08 2.117996e+07 3.508640e+07 8.509703e+06 #4 AB ainf1weekApis 6 1.605000e+13 1.1375e+13 3.083750e+13 1.803333e+13 7.899283e+12 3.224869e+12 #5 AB clean1weekBombus 18 0.000000e+00 0.0000e+00 2.755900e+07 4.085784e+06 9.564236e+06 2.254312e+06 #6 AB clean2weeksBombus 18 0.000000e+00 0.0000e+00 1.154035e+08 1.127090e+07 4.578021e+07 1.079050e+07 ######### virus.transmission.part2.apis.wo0<-subset(virus.transmission.part2.apis,viruslog35!="0") virus.transmission.part2.apis.wo0<-drop.levels(virus.transmission.part2.apis.wo0) str(virus.transmission.part2.apis.wo0) virus.transmission.part2.apis.wo0.summary<-ddply(virus.transmission.part2.apis.wo0,c("treatment","status_species"),summarize, N=length(virus35), median=quantile(virus35,0.5), lowerCI=quantile(virus35,0.025), upperCI=quantile(virus35,0.975), mean=mean(virus35), sd=sd(virus35), se=sd/sqrt(N)) virus.transmission.part2.apis.wo0.summary # treatment status_species N median lowerCI upperCI mean sd se #1 AA ainf1weekApis 6 2.120000e+13 9.402500e+12 4.011250e+13 2.252667e+13 1.126095e+13 4.597265e+12 #2 AA clean1weekApis 3 5.720000e+12 2.927000e+12 7.487000e+12 5.360000e+12 2.420165e+12 1.397283e+12 #3 AA clean2weeksApis 13 8.102012e+06 1.453993e+06 1.106289e+08 2.769687e+07 3.802457e+07 1.054612e+07 #4 AB ainf1weekApis 6 1.605000e+13 1.137500e+13 3.083750e+13 1.803333e+13 7.899283e+12 3.224869e+12 #5 AB clean1weekBombus 3 2.666742e+07 1.905913e+07 2.814047e+07 2.451470e+07 5.130371e+06 2.962021e+06 #6 AB clean2weeksBombus 2 1.014381e+08 1.299036e+07 1.898858e+08 1.014381e+08 1.316674e+08 9.310288e+07 ##########summary of all recipients together virus.transmission.part2.apis.wo0.summary2<-ddply(virus.transmission.part2.apis.wo0,c("treatment","status"),summarize, N=length(virus35), median=quantile(virus35,0.5), lowerCI=quantile(virus35,0.025), upperCI=quantile(virus35,0.975), mean=mean(virus35), sd=sd(virus35), se=sd/sqrt(N)) virus.transmission.part2.apis.wo0.summary2 # treatment status N median lowerCI upperCI mean sd se #1 AA ainf 6 2.120000e+13 9.402500e+12 4.011250e+13 2.252667e+13 1.126095e+13 4.597265e+12 #2 AA clean 16 1.091684e+07 1.494347e+06 6.882500e+12 1.005023e+12 2.334410e+12 5.836026e+11 #3 AB ainf 6 1.605000e+13 1.137500e+13 3.083750e+13 1.803333e+13 7.899283e+12 3.224869e+12 #4 AB clean 5 2.666742e+07 9.367564e+06 1.779087e+08 5.528406e+07 7.824574e+07 3.499256e+07 ##########summary of all donors together virus.transmission.part2.apis.wo0.summary3<-ddply(virus.transmission.part2.apis.wo0,c("status"),summarize, N=length(virus35), median=quantile(virus35,0.5), lowerCI=quantile(virus35,0.025), upperCI=quantile(virus35,0.975), mean=mean(virus35), sd=sd(virus35), se=sd/sqrt(N)) virus.transmission.part2.apis.wo0.summary3 # status N median lowerCI upperCI mean sd se #1 ainf 12 1.95000e+13 9.168500e+12 3.94775e+13 2.028000e+13 9.56608e+12 2.761489e+12 #2 clean 21 1.86587e+07 1.561602e+06 6.65000e+12 7.657446e+11 2.06869e+12 4.514251e+11 ####### Fisher exact tests data.frame<-read.table(file.choose(),header=T) # create a dataframe, in this case to test effect of Apis as donor versus Bombus as donor, with cage as the unit of replication dfExp2 <- data.frame("infected" = c(8, 0), "noninfect" = c(4, 12), row.names = c("donorA", "donorB")) # display the dataframe dfExp2 # Fisher exact test for difference in proportions fisher.test(dfExp2)